2.1.

Solving PDEs via Circuit Models: Theory

In previous works, it has been shown how the solution of PDEs (such as the Poisson equation ${\nabla }^{2}g=0$, where $g$ is the function to be solved for) can be computed by exploiting a network of lumped circuit elements^{40,42,44,59,62} connected in a grid-like lattice. In those works,^40,59 the structure consisted of a unit cell made of periodic junctions of lumped elements (see also Supplementary Material). Each junction can be formed by four lumped elements having an impedance $Z_L$, which are, each of them, connected to adjacent junctions, in this way forming a grid. The governing equation for the voltage distribution of such a lumped element-based network can be found by considering the division of current at the junctions according to Kirchhoff’s current law, as follows (see Supplementary Material for further details)

Now, what would be an equivalent circuit model that could emulate Eq. (3)? To answer this, in this work we propose the circuit shown in Fig. 1(a). It consists of a square lattice of T-circuits connected in a series configuration to form a junction. As observed, and in line with other works using lumped elements,^40,59,62 the proposed circuit is periodic. However, the junctions are now changed from four interconnected lumped elements to four interconnected T-circuits. As will be discussed below, the impedances forming the T-circuits can be physically emulated using transmission lines (parallel plate waveguides in our case) loaded with dielectric slabs, as shown in Fig. 1(b). It is important to note that here we implement series junctions of T-circuits emulated by parallel plate waveguides as transmission lines. However, it is also possible to use $\mathrm{\Pi }$-circuits emulated by transmission lines connected in a parallel^64–66 configuration, as explained in the Supplementary Material. The selection of T-circuits is to enable the implementation of the metatronic elements used to emulate the electronic circuits, while maintaining a small junction cross section (compared to the operating wavelength of the structure), a requirement to ensure equal splitting of the incident signals at the junction according to Kirchhoff’s laws.^23,47,67 Now, as the T-circuits from Fig. 1(a) [and their equivalent waveguide-based model from Fig. 1(b)] are connected to a junction in a series configuration, the flow of current through each of the connected T-circuits will form a rotation around the center of the junction [see Fig. 1(a)]. Consider a junction in Fig. 1(a) where there is a rotating current $I_0$. By looking into one of the connected T-circuits at this junction, for example the top T-circuit, one can calculate the voltage $V_1$ across the circuit (at the input of the junction, see Fig. 1) by considering the difference between the rotating current at the central junction $I_0$ and the top junction $I_1$ as $V_1=Z_p(I_1-I_0)-Z_s{I}_0$, where $Z_p$ and $Z_s$ are the parallel and series impedance values of the T-circuit, as shown in Figs. 1(a) and 1(b) (full details of this calculation can be found in the Supplementary Material). The magnitude and direction (clockwise/counterclockwise) of the rotating current at the junction can be found by solving Kirchoff’s voltage law ($\sum _a{V}_a=0$) at the junction using the associated voltages from each T-circuit, with $a$ representing the T-circuit from which the voltage is obtained ($a=1,2,3,4$, each representing the top, right, bottom, and left T-circuit, respectively). After accounting for the division of current at the T-circuits, one can arrive at the following expression relating the rotating current at a junction to the rotating currents at each of the adjacent junctions in the network (see Supplementary Material for the full calculation):

2.2.

Emulating T-Circuit Lumped Elements via Metatronic Circuit Elements

In the previous section, it was described how periodic T-circuits formed by lumped circuit elements within a network may be used to calculate the solution of PDEs when arranged in a grid. To implement these lumped elements at high frequencies (microwaves in our case), here we opted for the alternative of exploiting metatronic elements: thin dielectric, or metallic inserts, that may be engineered to emulate the performance of lumped circuit elements at different frequency ranges (from microwaves up to visible).^49–52,68 To mimic the series and parallel impedances of the T-circuit from Fig. 1 (i.e., $Z_s$ and $Z_p$, respectively) we consider a metatronic circuit formed by three thin films (either dielectric or metallic) embedded within a host medium (vacuum in our case, ${\epsilon }_r=1$, ${\mu }_r=1$), as shown in Fig. 1(b). When considering a transverse EM signal (TEM, as it is in the case of the proposed parallel plate waveguides) impinging onto a single one of these dielectric or metallic slabs under normal incidence, the slab can be defined by a parallel impedance with a value⁵²

From Eq. (7), it can be seen how, due to the impedance transformation, materials with $\text{Re}({\epsilon }_r)>0$ and $\text{Re}({\epsilon }_r)<0$ are now expressed as series inductors and capacitors, respectively. Based on this, we can emulate the required T-circuit to solve the PDE described in the previous section [see Fig. 1(a) and Eq. (5)] by designing the metatronic elements using Eqs. (6) and (7). The proposed metatronic element-based structure that can emulate a lumped-element based T-circuit is shown in Fig. 1(b). The metatronic circuit consists of a parallel plate waveguide having two thin slabs (yellow slabs) placed in between two ${\lambda }_0/4$ regions filled with air, in so doing, the series impedances of the T-circuit are emulated. For the parallel impedance, a thin slab is also used (brown slab), as shown in Fig. 1(b). This structure is then arranged in a network configuration forming junctions as explained above [schematically shown in Fig. 1(a)]. Importantly, as the dielectric slabs are embedded within parallel plate waveguides [see Fig. 1(b)], Eqs. (6) and (7) for the calculation of metatronic elements are valid in our case, as we are also working with an incident TEM signal. Now, when selecting the dimensions of the waveguides in which the dielectric slabs are embedded, one must ensure that the separation between plates is small compared to the incident wavelength in order to limit the impact of fringing fields at the edges of the waveguide (small cross sections are also required when working with waveguides connected in a parallel configuration,^64,65,70 as is the case shown in the Supplementary Material). Finally, in this work, we consider that all waveguides have the same filling materials (vacuum ${\epsilon }_r=1$, ${\mu }_r=1$) and cross section (i.e., equal impedance).

With this configuration, the numerical results of the out-of-plane $H_z$-field distribution, calculated at the junctions for a $3\times 3$ waveguide network are shown in Fig. 1(c). A monochromatic (10 GHz, ${\lambda }_0=30\mathrm{mm}$) incident signal is excited from the left waveguide of the top-left junction, here labeled as junction 1 in Fig. 1(c). Using Eqs. (6) and (7), the metatronic elements were designed to emulate the parallel and series impedance values of $Z_p=2.5i{Z}_0$ and $Z_s=-0.9i{Z}_0$, respectively. These impedances correspond to $h$ and $k$ values of 0.4 and 3 (arbitrary units), respectively. It should be noted that the units of $h$ and $k$ are informed by the physics of the equation being solved. For example, in an EM problem, $h$ and $k$ have units of meters (m) and ${\mathrm{m}}^{-1}$, respectively. However, here we first consider both of them as unitless in order to calculate a general solution of a PDE. With these considerations, the metatronic elements have a thickness of 0.2 mm (${\lambda }_0/150$) each, with permittivity values of ${\epsilon }_r=9.54$ and ${\epsilon }_r=21.44$ for the slabs emulating parallel and series impedances, respectively. To corroborate the T-circuit design, a full ABCD matrix analysis³ of the structure presented in Fig. 1(b) (upper panel) was performed (see models in the Supplementary Material). Based on these calculations, the reflection and transmission coefficients varied slightly from the ideal design using the T-circuit model, as expected due to the nonzero thickness of the dielectric slabs for the structure in Fig. 1(b). To overcome this difference, one can simply optimize the geometrical or EM parameters of the structure shown in Fig. 1(b) by, for instance, slightly changing the distance ${\lambda }_0/4$ between the slabs, by changing $\mathit{w}$ of the slabs, or by modifying ${\epsilon }_{\mathrm{r}}$ of the slabs. However, in this first example, the structure from Fig. 1(b) is kept unchanged with the same dimensions as described above, but now we can retrieve the equivalent $Z_p$ and $Z_s$ values from the ideal T-circuit model that match the reflection and transmission coefficients from the ABCD method of Fig. 1(b) (details on this calculation can be found in the Supplementary Material). The calculated impedances are $Z_p=2.1522i{Z}_0$ and $Z_s=-0.9311i{Z}_0$, respectively (corresponding to $h=0.4646$, $k=2.8313$, i.e. a small variation of the PDE parameters). In this example, these new values of $Z_p$ and $Z_s$ are then used to calculate the theoretical values for the $H_z$-field based on the ideal circuit model and these results are compared to the numerically calculated results using simulations, as will be shown below.

The numerical $H_z$-field values (amplitude and phase) at junctions 1 through 9 are presented in the middle panel of Fig. 1(c). These simulations were obtained using the full-wave simulation software CST Studio Suite^® (see the methods section for details of the simulation setup). As shown in Fig. 1(a) and discussed above, the solution of a PDE can be obtained by looking at the $H_z$-field (and hence the circulating current) at the center of a waveguide junction. However, as Eq. (2) requires four terms from neighboring junctions $[g(x,y+h)$, $g(x+h,y)$, $g(x,y-h)$, and $g(x-h,y)$ from the top, right, bottom, and left junctions, respectively], only a junction that is fully surrounded by four other junctions will be able to produce a solution to the PDE. By looking at the structure from Fig. 1(c) (left panel) only junction 5 of the $3\times 3$ network will satisfy Eq. (5). The phasor representation of the numerically calculated values for the $H_z$-field at junction 5 and its neighboring junctions are presented in the rightmost panel of Fig. 1(c) as “×” symbols, alongside the theoretical values for an ideal grid of metatronic elements calculated using the T-circuit model with the impedance values obtained via the ABCD matrix method (i.e., $Z_p=2.1522i{Z}_0$ and $Z_s=-0.9311i{Z}_0$), represented as “○” symbols (see Supplementary Material, as mentioned above). As observed, an excellent agreement between the numerical and theoretical results is obtained. Finally, the numerical results of the $H_z$-field for the adjacent junctions can also be used to calculate what would be the ideal value for the $H_z$-field at junction 5 that will satisfy Eq. (3). This can be calculated as $H_5=(H_2+H_4+H_6+H_8)/(4-h^2{k}^2)$, where the subscripts indicate the junction number [see Fig. 1(c)]. By doing this, the error between the numerical solution and the ideal solution of the Helmholtz equation can be calculated as the difference between the numerical results for $H_z$-field at junction 5 and the ideal value of $H_5$ that satisfies Eq. (3) (as described before). The difference between these results is $\sim 7.23\%$, demonstrating how a small error is achieved and the proposed structure is capable of calculating a solution to the Helmholtz equation at junction 5.

2.3.

Scaling of the Calculated PDE Solutions

As discussed in Fig. 1(c), the proposed structure can be used to calculate the solution of a PDE at the junction of waveguides as long as the junction is interconnected to an adjacent junction in all directions [for instance junction 5 in Fig. 1(c)]. Hence, if one wants to calculate a more detailed solution to the PDE (Helmholtz equation in this study, for instance), it is necessary to increase the size of the network by introducing more junctions that are interconnected by waveguides filled with slabs (to emulate T-circuits) at their top, bottom, left, and right directions. In other words, if a sampling point of the solution of the PDE is represented by a junction, it is possible to increase the resolution of the calculated solution by increasing the number of junctions. To address this, here we consider networks of waveguide-based metatronic T-circuits consisting of $N\times M$ junctions along the horizonal and vertical directions, respectively [see Fig. 2(a)]. A three-dimensional (3D) representation (top view) of a central section of a larger network of junctions is included in Fig. 2(b), showing how there are now more junctions fully interconnected to adjacent junctions in all directions compared to the scenario discussed in Fig. 1(c). As explained before in Eq. (5), the scaling and spatial sampling of the calculated solution are determined by the chosen $h$ value of the structure. This means that it is possible to obtain the solution of a PDE with an increased resolution by increasing the physical size of the network (number of junctions) while keeping the same simulation space. This will effectively represent a smaller separation between sampling points in the simulation space $h$. Alternatively, a higher resolution of the solution for a PDE could be achieved without altering the size of the network but instead reducing $h$. This will mean that the simulation space will become smaller, as will be discussed below.

Fig. 2

Scaling of PDE results. (a) Nodal representation of an arbitrarily sized network of interconnected metatronic elements. (b) (Top view) $3\times 3$ sections of a larger waveguide network extending in all directions. The perfect electric conductor (PEC) regions for the waveguides are represented as gray blocks. The red and yellow slabs represent the elements that enable us to emulate the parallel and series metatronic elements, respectively. (c)–(e) Full-wave numerical results of the $H_z$-field distribution of a $25\times 25$ waveguide network. As in Fig. 1(c), a 10 GHz monochromatic input signal is excited from the left waveguide of the top-left junction. These results have been normalized so that the out-of-plane $H_z$-field seen at the top left junction is unity. Here, ${\lambda }_1$ is the wavelength of the PDE in the simulation space, while ${\lambda }_0$ is the wavelength of the incident signal in free space. (c) Numerical results for the case where there are no dielectric slabs present within the connecting waveguides; see inset. (d), (e) Analytical (left), theoretical (middle), and numerical (right) results of the same setup from (c) (same color scale applies here) but now when the waveguides are loaded with the dielectric slabs emulating metatronic elements with effective impedances $Z_p=2.498i{Z}_0$, $Z_s=-0.9003i{Z}_0$ and $Z_p=5.001i{Z}_0$, $Z_s=-0.4501i{Z}_0$, respectively. These values correspond to PDEs with parameters $h=0.4003$, $k=2.999$ $({\lambda }_1=2.095)$, and $h\approx 0.2$, $k=3.001$ $({\lambda }_1=2.094)$, respectively. The EM and geometrical parameters after optimization are $L_1=L_4=7.394$ mm ($\sim 0.2465{\lambda }_0$), $L_2=L_3=7.307\mathrm{mm}$ ($\sim 0.2436{\lambda }_0$), $w_s=0.2111\mathrm{mm}$ ($\sim 7.0366\times {10}^{-3}{\lambda }_0$), $w_p=0.1741$ mm ($\sim 5.8033\times {10}^{-3}{\lambda }_0$), ${\epsilon }_s=21.5$, and ${\epsilon }_p=12$ for the results presented in panel (d), where $L_1$, $L_2$, $L_3$, and $L_4$ are the lengths of the impedance transforming waveguides from left to right, as shown in Fig. 1(b). $w_s$, $w_p$, ${\epsilon }_s$, and ${\epsilon }_p$ are the widths and permittivity values of the slabs representing the series and parallel elements, respectively. These parameters are $L_1=L_4=7.390\mathrm{mm}$ ($\sim 0.2463{\lambda }_0$), $L_2=L_3=7.294\mathrm{mm}$ ($\sim 0.2431{\lambda }_0$), $w_s=0.2201\mathrm{mm}$ ($\sim 7.3367\times {10}^{-3}{\lambda }_0$), $w_p=0.1911\mathrm{mm}$ ($\sim 6.37\times {10}^{-3}{\lambda }_0$), ${\epsilon }_s=10.80$, and ${\epsilon }_p=6.000$ for the results presented in panel (e). The line plots in the rightmost panels of (d), (e) show the numerical (green triangles), analytical (red squares), and theoretical (gray circles) results of the magnitude of the $H_z$-field taken along a straight line from the top-left to bottom-right corners of the simulation space, respectively.

Consider for example, an $N=M=25$ network of junctions where 96 out of the total 625 junctions are placed at the edges and used as boundary junctions; i.e., as they are not surrounded by junctions in all directions, they cannot be used as part of the solution of the PDE. The analytical (calculated using the Huygens–Fresnel principle), theoretical (calculated by solving the governing equations of an ideal T-circuit-based network), and numerical (CST Studio Suite) results of the out-of-plane $H_z$-field distribution at the junction of waveguides emulating the T-circuits are shown in Fig. 2(d). Here we consider the incident signal as in Fig. 1(c) (10 GHz monochromatic incident signal excited at the left waveguide of the top left junction). The metatronic elements were designed using Eqs. (6) and (7) to emulate the parallel and series impedance values $Z_p=2.5i{Z}_0$ and $Z_s=-0.9i{Z}_0$ (impedances corresponding to $h=0.4$, $k=3$). After optimizing the geometric (distance between slabs, $w_s$ and $w_p$) and EM parameters (permittivity values of the slabs ${\epsilon }_r$) of the design, the emulated impedance values calculated from the ABCD matrix method where $Z_p=2.498i{Z}_0$ and $Z_s=-0.9003i{Z}_0$ ($h=0.4003$, $k=2.999$), respectively (see dimensions and EM parameters in the caption of Fig. 2). As is shown in Fig. 2(d), there is an excellent agreement between all the results with the solution of the PDE resembling the radiation of a dipole. Note that this is an expected result due to a monochromatic signal being applied only from the top-left junction. As a result, the analytical calculations using the Huygens–Fresnel principle will only have a single point as the radiating source, as observed in Fig. 2(d). The results presented in Fig. 2(d) are the calculated $H_y$-field values once the structure has been allowed to settle into a steady state. A study of the time taken to reach this state ($\approx 30\mathrm{ns}$) is shown in the Supplementary Material.

For completeness, and in order to demonstrate the impact of the chosen $h$ value as a scaling parameter, the metatronic elements were also modified so that the calculated solution of the PDE would resemble a zoom-in image of the top-left quadrant of the solution calculated in Fig 2(d). To do this, the metatronic elements were designed to emulate the impedance values $Z_p=5i{Z}_0$ and $Z_s=-0.45i{Z}_0$ ($h=0.2$, $k=3$). As observed, now $Z_p$ and $Z_s$ have been doubled and halved [compared to the values for Fig. 2(d)], respectively. This in turn can enable a solution to the same PDE as Fig. 2(d) (same $k$ value) but now with the separation between sampling points $h$ being halved. As a result, the simulation space is reduced to a quarter, effectively zooming into the top-left quadrant of Fig. 2(d) (dashed white square). The emulated impedance values after optimization (see dimensions, EM parameters in the caption of Fig. 2) are $Z_p=5.001i{Z}_0$ and $Z_s=-0.4501i{Z}_0$ ($h\approx 0.2$, $k=3.001$); the results are shown in Fig. 2(e). Observing the number of oscillations present in Fig. 2(e) and comparing it to the top-left quadrant of Fig. 2(d), it can be qualitatively confirmed that the results presented in Fig. 2(e) indeed represent a zoomed-in solution to the same PDE as Fig. 2(d). Furthermore, as the same number of sampling points (junctions) are used to display a smaller area in the simulation space [one-quarter of Fig. 2(d)], the resolution of the calculated solution is quadrupled. For completeness, line plots showing the analytical, theoretical, and numerical magnitude of the $H_z$-field along a line from the top-left to bottom-right of each simulation space are shown in the last column of Figs. 2(d) and 2(e). In these figures, all three plots show the characteristic 1/distance decay associated with a radiating dipole in 2D. However, in the numerical and theoretical plots, there are oscillations. As will be discussed in more detail later, this is due to the impact of reflections produced at the boundary junctions. Finally, to demonstrate that the waveguides along with the loaded dielectric slabs are indeed emulating the T-circuits from Fig. 1(a) and are responsible for the PDE-solving performance of the proposed network, the case when the waveguides are not loaded with the dielectric slabs is also shown in Fig. 2(c). This is done by replacing the dielectric slabs with air so that the physical length of the connection between junctions is the same as in Figs. 2(d) and 2(e). As observed, the solution does not resemble the solutions to the Helmholtz equation shown in Figs. 2(d) and 2(e), indicating that it is indeed the dielectric slabs and the waveguides (as a whole) that are responsible for the PDE-solving behavior of the structure.

2.4.

Calculating Solutions to Dirichlet Boundary Value Problems

The examples presented in the previous sections have considered only a single input signal applied from a waveguide connected at one boundary junction (top-left junction from Figs. 1 and 2). Here, it is shown how the proposed structure for PDE solving can also be used to calculate the solutions to Dirichlet boundary value problems.⁷¹ This is done by implementing simultaneous excitations from the different outer waveguides connected to the boundary junctions (from now on called boundary waveguides) around the network. For this, the left waveguide from the top-left junction is used as a reference. In this way, Dirichlet boundary conditions, such as $g=1$ or $g=0$, can be considered at each of the boundary junctions. These boundary conditions are physically implemented in the structure from Fig. 1 by designing the inputs from the boundary waveguides such that the $H_z$-field at the center of each of the boundary junctions is 1 or 0, respectively. To implement a specific boundary condition at a boundary junction, it is important to carefully engineer the amplitude ratios and phase differences of the incident signals applied from the different boundary waveguides. This means that the boundary value at a given boundary junction (calculated by the $H_z$-field at the junction center) is determined by the superposition of the input signal applied from the connected boundary waveguide and those signals coming from the other boundary junctions from the network.

The required incident signals from the boundary waveguides can be directly calculated from the scattering matrix of the structure $\mathit{A}$.^3,64,66,70 As an example, consider a vector of incident signals defined by their $E$-field (the same process could be done considering voltages) $\mathit{x}=[x_1,x_2\dots ,x_{2(M+N)}{]}^T$, where $T$ indicates the transpose operation. These signals are used to excite the boundary waveguides connected to the boundary junctions. Due to the interaction of these incident signals with the PDE solving structure, a vector of output signals defined by the $E$-fields $\mathit{y}=\mathit{Ax}$^64,66,70 is created (observed at the same boundary waveguides), with $\mathit{y}=[y_1,y_2,\dots ,y_{2(N+M)}{]}^T$. The vector containing the complex values of the instantaneous $H_z$-field (which relates to the current rotating around the junction, as explained above) at each of these boundary junctions can be written as³ ${\mathit{H}}_b=(\mathit{x}-\mathit{y})/Z_0$, where $Z_0$ is the characteristic impedance of the boundary waveguides (as before, all the waveguides in the network have the same dimensions and filling materials, i.e., they have the same characteristic impedance) and $b=[\mathrm{1,2},\dots 2(N+M)]$ represents each of the boundary waveguides connected to the boundary junctions. Combining this expression for ${\mathit{H}}_b$ with the scattering equation for $\mathit{y}$ allows us to define a relation between the desired boundary values at the boundary junctions and the vector of incident signals from the boundary waveguides required to produce them:

With this configuration, two examples of Dirichlet boundary value problems that can be solved by our proposed structure using EM waves are presented in Fig. 3. For both cases, we exploit the same $25\times 25$ network of waveguide-based metatronic structures as in Fig. 2(e) (right panel). In the first example [Figs. 3(a) and 3(b)], the incident signals from the boundary waveguides have been selected using Eq. (8) to produce a Dirichlet boundary value of $g=1$ at the left boundary, while a value of $g=0$ is implemented for the top, right, and bottom boundaries junctions, represented by the red and gray-scale lines in Fig. 3(a), respectively. These boundary values mean that the calculated solution to the PDE will resemble a standing wave produced when a wave propagates from left to right and is reflected by the right boundary. For the analytical calculations in Fig. 3, the solution of the PDEs is calculated using the FEM from the build PDE toolbox of a commercial software^72,73 (MATLAB in our case; see the Appendix for more details about this calculation). To further compare the results shown in Fig. 3(b), the $H_z$-field values were extracted along the horizontal and vertical lines (at 13 junctions down from the top boundary and 12 junctions along from the left boundary, respectively) as shown in the top panels from Fig. 3(b); the results are shown in the bottom-left and bottom-right panels respectively, demonstrating an agreement among the analytical, theoretical, and numerical results. The second example of a Dirichlet boundary value problem is presented in Fig. 3(c), in which the incident signals are selected such that the magnitude of the $H_z$-field at each of the boundary junctions is the same. However, the phase is varied such that, along a closed path around the entire boundary, the $H_z$-field completes a full $2\pi$ phase loop. The analytical, theoretical, and numerical results of the $H_z$-field are shown in the top panels from Fig. 3(d). For completeness, the magnitude of the $H_z$ field was recorded along the diagonal lines from these panels; the results are shown on the bottom panels from Fig. 3(d), again, showing an agreement between the results and demonstrating the potential of the proposed structure for PDE solving.

Fig. 3

Solving Dirichlet boundary value problems. (a), (c) Schematic representations of a $25\times 25$ network of junctions of waveguide-based metatronic circuits excited with different Dirichlet boundary conditions. The metatronic elements are chosen considering $Z_p=5.001i{Z}_0$ and $Z_s=-0.4501i{Z}_0$ (see Fig. 2 caption for EM and geometrical parameters), corresponding to $h\approx 0.2$ and $k=3.001$ (${\lambda }_1=2.094$). (a) The left-hand boundary is set to $g=1$ while the top, right, and bottom boundaries are $g=0$. (c) The boundary conditions are such that the magnitude at each boundary junction is 1 but the phase along the boundary spatially varies from 0 to $2\pi$, counterclockwise. (b), (d) Analytical, theoretical, and numerical results of the scenarios from (a) and (c), respectively. The top panels from (b) and (d) represent the normalized instantaneous $H_z$-field values calculated at each of the junctions inside the network. The bottom panels from (b) and (d) represent the magnitude of the calculated analytical (red squares), theoretical (gray circles), and numerical (green triangles) $H_z$-field values along the dashed lines from the top panels.

2.5.

Open Boundary Value Problems

As a final study, here we show that it is also possible to use the proposed technique to solve open boundary value problems, such as the two examples presented in Fig. 4. As could be expected, the signals will need to be absorbed by the boundary waveguides mimicking an open boundary, similar to the well-known perfectly matching layers, PMLs.⁷⁴ In attempting to do this using the proposed metatronic-based network, using waveguide ports connected to the boundary waveguides, unwanted large reflections can be obtained in the calculated solution. This is due to the impedance mismatch between the boundary junctions and the rest of the junctions given that they are connected to three and four adjacent junctions, respectively. To tackle this, a solution is to extend the simulation space by including more junctions into the waveguide network and then extracting the PDE solution from a smaller portion of the simulation space, rather than the entire area. This allows for the signals within the network to decay as they propagate through the extended network, reducing the impact of reflections on the overall solution. With this configuration, we provide in Fig. 4 examples of open boundary value problems. Here we use a $150\times 100$ metatronic network (with the same $Z_p$ and $Z_s$ values as the examples presented in Fig. 2(e) and Fig. 3, $Z_p=5.001i{Z}_0$, $Z_s=-0.4501i{Z}_0$, $h=0.2$ and $k=3$). The simulation space used to produce the results in Fig. 4 is a $50\times 50$ subnetwork of the total $150\times 100$ network. The position of the subnetwork is selected such that the network extends for 50 junctions from the top, right, and bottom boundaries of the subnetwork, respectively, effectively emulating the behavior of a PML at these boundaries. No extra junctions are connected to the left-hand boundary of the subnetwork, allowing for incident signals to be excited at the boundary waveguides connected to these junctions.

Fig. 4

Lensing and particle scattering. (a), (d) Schematic representations of a $50\times 50$ subnetwork of waveguide-based metatronic circuits used to solve open boundary value problems. The $50\times 50$ grid is constructed with $Z_p=5.001i{Z}_0$ and $Z_s=-0.4501i{Z}_0$, corresponding to $h\approx 0.2$ and $k=3.001$ (${\lambda }_1=2.094$, where ${\lambda }_1$ is the wavelength seen inside the simulation space; see Fig. 2 caption for EM and geometrical parameters). (a) A 10 GHz monochromatic wave is excited at each of the boundary waveguides connecting to the left-hand boundary junctions of the subnetwork. The amplitude and phase of these signals is selected such that the $H_z$ values at the boundary junctions resemble the output signal from a lens designed to produce a focus at $x=1.432{\lambda }_1$ and $y=2.387{\lambda }_1$ (15 and 25 junctions in the $x$ and $y$ directions, respectively), represented by a gray spot. (d) As in panel (a), a 10 GHz monochromatic signal is excited at the left-hand boundary waveguides, now with amplitude and phase selected such that the $H_z$ values at the boundary junctions resemble a planewave. A $0.9549\times 0.9549{\lambda }_1$ ($10\times 10$ junctions) $g=0$ insert is placed at the center of the simulation space, represented by a white block. (b), (e) Theoretical (left) and analytical (right), power distribution for the scenarios presented in panels (a) and (b), respectively. In panels (b) and (e), results are normalized with respect to the power distribution at the focus and the maximum standing wave, respectively. (c), (f) Theoretical (gray circles) and analytical (red squares) values of power distribution along the vertical (top) and horizontal (bottom) lines drawn in panels (b) and (e), respectively.

The theoretical and analytical results of the normalized power distribution calculated at the center of the junctions for the two scenarios from Figs. 4(a) and 4(d) are shown in Figs. 4(b) and 4(e), respectively (numerical simulations are not shown due to the large size of the whole structure). The analytical results from Figs. 4(b) and 4(e) are calculated using the Huygens–Fresnel principle, considering each boundary junction as a radiating dipole and using the FEM from the PDE toolbox of commercial software (MATLAB,^72,73 as in Fig. 3), respectively. Let us first focus on the results from Figs. 4(a)–4(c), which show an example of focusing/lensing. In this case, the boundary waveguides connected to the boundary junctions on the left-hand side of the network are excited such that the $H_z$-field distribution along the boundary junctions resembles an output signal traveling away from a lens (designed to produce a focus at a position $x=1.432{\lambda }_1$, $y=2.387{\lambda }_1$, where ${\lambda }_1=2.094$ is the wavelength of the PDE. In the simulation space, the focus would appear 15 and 25 junctions along the $x$ and $y$ axes of the network, respectively; see Figs. 4(a) and 4(b). As shown in Fig. 4(b), a clear focus is produced at the positions $x=1.623{\lambda }_1$, $y=2.387{\lambda }_1$ (17 and 25 junctions along the $x$ and $y$, respectively), and $x=1.432{\lambda }_1$, $y=2.387{\lambda }_1$ for the theoretical and analytical calculations, respectively, demonstrating a good agreement. Note that from the theoretical results, there are some slight spatial ripples. This can be attributed to the reflections produced at the boundary junctions of the entire network, as mentioned above. As one would expect, the signal inside the network would asymptotically decay to zero if the network was infinitely long. However, as the network is finite (extended by 50 junctions from the top, right and bottom of the $50\times 50$ subnetwork), reflections, while small, are still present. For completeness, the theoretical and analytical results extracted along horizontal and vertical lines from Fig. 4(b) are shown in Fig. 4(c). As is shown, despite the presence of small reflections in the PDE solution, the focus is still well defined, with both results in good agreement.

A second and final example of an open value problem in shown in Figs. 4(d)–4(f), where the simulation space has been divided into two regions. Here, an insert is considered to be placed within a background medium [with the latter being the same as the one used for Figs. 4(a)–4(c)]. The insert is modeled by removing the waveguide-based metatronic T-circuits (i.e., effectively removing any junctions within that region). By doing this, any incident signal upon this region will be reflected, enabling the junctions at the boundary to that region to act as boundary junctions with $g=0$. Physically, this could be implemented by using a perfect magnetic conductor such as doped ENZ materials,⁷⁵ or by replacing the metatronic T-circuits connecting to this region with PEC-ended stubs of length ${\lambda }_0/4$. The insert, shown in Figs. 4(d) and 4(e), is a $0.9549{\lambda }_1\times 0.9549{\lambda }_1$ square centered at the middle of the simulation space, i.e., a $10\times 10$ grid in the center of the $50\times 50$ subnetwork. The theoretical and analytical results of the power distribution for this scenario are shown in Fig. 4(e) when considering a planewave excitation (boundary waveguides are excited such that the boundary junctions on the left-hand side of the network fulfill a boundary condition of $g=1$). As observed, there is a clear agreement between the results, demonstrating how the incident wave is reflected and scattered by the insert. For completeness, the power distribution along the vertical and horizontal lines from Fig. 4(e) are shown in Fig. 4(f). The small differences between the analytical and theoretical results are due to the finer mesh used in the FEM method for the analytical results compared to the number of junctions to discretize the simulation region using the network of waveguides. This is an expected result, which may be reduced by increasing the resolution of the calculated PDE solution in the metatronic network using the methods discussed in the previous sections. However, the results from Figs. 4(e) and 4(f) can be considered as a good approximation of the analytical solution. The proposed technique can be translated to different frequency ranges such as the optical regime via the emulation of metatronic elements within waveguides with the ability to split the incident signal in all directions,⁴⁸ allowing its potential experimental implementation in integrated photonic circuits.