Generalized additional boundary conditions for wire media

We generalize additional boundary conditions (ABCs) for wire media by including arbitrary wire junctions with impedance loading. Special attention is given to the conditions at the interface of two uniaxial wire media with metallic patches at the junction. The derived ABCs are validated against full-wave numerical simulations.


3
The purpose of this work is to further extend the theory of [12,13] and study the general case where metallic wires are connected to arbitrarily distributed or lumped loads or to a second WM with different parameters. We will show that it is possible to derive an ABC in a quasi-static approximation, in a manner similar to what has been done when developing a quasi-static model of spatial dispersion in wire media [19], and we will present numerical results that support our theoretical findings.

General form of additional boundary conditions (ABCs) in uniaxial wire media
In a recent paper by Maslovski and Silveirinha [19], it has been shown that spatial dispersion in dense wire media can be explained with simple quasi-static considerations. Namely, in this model a current I z , a charge q per unit length and an additional potential ϕ due to this charge are associated with each wire so that the field equations in a uniaxial WM with wires oriented along the z-axis can be written as ∂ I z ∂z ≡ −jω q = −jωC ϕ , where ε 0 and µ 0 are the permittivity and permeability of the host material in which the wires are immersed (e.g. vacuum in the simplest case); E and H are the averaged (macroscopic) electric and magnetic fields in the medium; the angular brackets . . . represent any suitable interpolating (averaging) operator that smoothens the microscopic quantities defined at discrete wires and makes them continuous through all volume; E z ≡ẑ · E; A cell is the area of the unit cell, which for a square lattice of wires is A cell = a 2 , where a λ is the lattice period; and L, C and Z w are the effective inductance, capacitance and loss impedance per unit length of a wire, respectively. The additional potential at each wire is defined as an integral of the radial component of the electric field in the vicinity of a wire. With a suitable choice of the coordinate system (see figure 1), this integral can be written as ϕ(z) = a/2 r 0 e x (x, z)dx, where r 0 is the radius of the wire, e x is the x-component of the microscopic electric field around the wire and the integration is done till the middle point in a pair of neighboring wires. With this definition of the additional potential, the effective capacitance is found as C = 2πε 0 /log[a 2 /4r 0 (a − r 0 )] and the effective inductance as L = (µ 0 /2π )log[a 2 /4r 0 (a − r 0 )] [19]. It is seen that in an unloaded uniaxial WM, LC = ε 0 µ 0 . More generally, the additional potential could be as well regarded as the average potential difference between a given wire and the boundary of the corresponding unit cell.
A very important property of the system (1)-(5)-which was not noticed in [19]-is that it corresponds to a local homogenization model. Indeed, the dynamics of the 'state variables' E, H, I z and ϕ are described by simple local relations: the pertinent spatial derivatives of the state variables at a given point only depend on the values of the state variables at that same point. This is seen from the fact that equations (1)- (5) can be written as a first-order differential equation for the state vector = [H, E, I z , ϕ ] T :L = −jωM , whereL =L(∂/∂ x, ∂/∂ y, ∂/∂z) is a first-order linear differential operator, and M = M(ε 0 , µ 0 , L , Z w , C, A cell ) is a material matrix that describes the response of the medium. Note that the material matrix does not include any integro-differential operators and, thus, describes a local response. This property may look surprising at first given that the WM is known to have a strongly non-local response [15], as discussed before and emphasized in many papers. How can we reconcile the non-local response of the WM with such a local model? It is simple: we managed to obtain a local model for the material only because we complemented the macroscopic Maxwell equations (1)-(2) with two other equations, (4) and (5), that describe the dynamics of I z and ϕ . These two additional state variables characterize internal degrees of freedom of the material and have a clear physical meaning as the current and potential associated with a given wire. Thus, as long as the effective medium is characterized by the four parameters E, H, I z and ϕ , it can be described using a local homogenization model, i.e. equations (1)- (5). On the other hand, as shown in [19], it is possible to write I z and ϕ in terms of the macroscopic electromagnetic fields; however, in such a case the homogenization model becomes spatially dispersive because then it is not possible to write J in terms of the macroscopic fields through a local relation. In other words, the effective medium can be described using a local model provided one introduces suitable additional state variables. The local homogenization model (1)-(5) provides the natural framework to study a problem involving interfaces. Indeed, in contrast to the traditional non-local model of [15], equations (1)- (5) are defined in the space domain and thus hold even relatively close to the interfaces. Moreover, the need for ABCs is quite obvious from the formulation (1)-(5): these are nothing more than the boundary conditions satisfied by the additional state variables I z and ϕ . In comparison to any local magnetodielectric or a vacuum, equations (1)-(5) add two more scalar functions that are the solutions of a pair of differential equations of first order. Therefore, in general, there will be a need for two ABCs imposed on these scalar quantities (in the case of a single-sided wire junction, a single ABC is sufficient, as discussed later). The boundary conditions on vector fields E and H can be obtained in the usual manner from (1) and (2), which leads to standard continuity conditions for the tangential components of the fields when there are no surface-bound currents at an interface.
Let us consider a junction between two uniaxial wire media with the same lattice period and with different wire radii, as shown in figure 1. For simplicity, first we will consider the case where the wires are electrically connected at the junction. To derive an ABC, we consider the circulation of the electric field over the path ABC D shown in the figure. The path starts at point A at the surface of the first wire, then goes down to the point B that is at the symmetry line in between the two wires, then follows along this symmetry line until point C and then goes up to point D that is at the surface of the second wire, and then goes back to point A along the wire surfaces. We assume that the segments along the z-axis are such that l 1,2 λ and l 1,2 a. Then, it can be seen that the integrals over the segments AB and DC are given by the additional potentials ϕ 1 (−l 1 ) = ϕ | z=−l 1 and ϕ 2 (l 2 ) = ϕ | z=l 2 .
Therefore, the circulation of the electric field over the whole path can be written as follows: where e is the microscopic electric field in the space surrounding the wires, e z (a/2, z) is the z-component of e along the line of symmetry and I z (z) is the microscopic (i.e. quickly varying close to the junction) current in the wires on both sides of the interface; away from the junction this current coincides with the macroscopic one: . To a first approximation, consistent with the analysis of [19], away from the junction, we can identify e z calculated at the line of the symmetry of the two wires with the macroscopic field E z : e z (a/2, −l 1 ) . Hence, because l 1,2 λ we can rewrite (6): where ϕ z = O(1), λ → ∞ is a quasi-static correction to the potential due to non-uniformity in the z-component of the microscopic electric field and the current in the vicinity of a junction. This correction may be neglected when both r 1,2 a. On the other hand, the circulation of the electric field equals the negative time derivative of the magnetic flux e · dl = −jω = jωL 1 l 1 I 1 | z=−l 1 + jωL 2 l 2 I 2 | z=l 2 + jω where = O(1), λ → ∞ is a quasi-static correction to the magnetic flux due to the nonuniformity of the microscopic magnetic field in the vicinity of the junction, which we may also neglect when r 1,2 a.
We note that the additional potential ϕ(z) is a function that varies smoothly with z away from the junction, i.e. it essentially behaves in the same way as its macroscopic average ϕ(z) . On the other hand, the boundary condition we are looking for must connect the values of ϕ (rather than the values of ϕ) exactly at the interface, i.e. at the point where the material properties change abruptly and the meaning of any macroscopic quantity is arguable. Therefore, we cannot simply let l 1,2 → 0 in equation (7), and instead we must keep l 1,2 ∼ a in the evaluation of the term that describes the fluctuations of the microscopic fields near the junction, i.e. ϕ z . Likewise, if we let l 1,2 → 0 in equation (8), the term must be evaluated for l 1,2 ∼ a. Taking  6 into account these results and defining ϕ 1 (0 − ) ≡ ϕ | z→0 − and ϕ 2 (0 + ) ≡ ϕ | z→0 + , we easily obtain the following boundary condition for the additional potential: When r 1,2 a the boundary condition (9) is simply the continuity of the additional potential across the interface. Note that if we had assumed l 1,2 → 0 already in (6), we would always obtain ϕ 2 (0 + ) = ϕ 1 (0 − ). Later in the text we will consider some cases where this condition does not hold and one must use (9). The boundary condition for the current can be obtained as follows. We integrate the divergence of the microscopic current density in the wires ∇ · j w = −jωρ w over the volume of the junction that is in between the planes z = −l 1 and z = l 2 (see figure 1; the integral is over the volume of a single wire): On the other hand, this integral equals where the quasi-static correction Q = O(1), λ → ∞ accounts for additional charge at the junction due to non-uniformity in the charge distribution. Combining (10) and (11) and using the same arguments as above for the potential, we obtain the following boundary condition for the current: Again, when r 1,2 a the correction term is small and can be neglected, in which case we obtain a simple continuity condition for the current: It is worth noting that the conditions (9) and (12) derived in this section are rather general and applicable also to the cases when there are bulk insertions or loads at the wire junction, because a voltage drop at such an insertion can be taken into account by the terms on the righthand side of (9). Indeed, inserting a load produces a local inhomogeneity in the distribution of the electric and magnetic fields close to the junction. Furthermore, if an insertion possesses a significant capacitance (like, for example, a metallic patch-see the next section) then it will also produce a local inhomogeneity in the charge distribution that will be taken into account by the term on the right-hand side of equation (12).
It is interesting to note that equations (9) and (12) yield two independent ABCs. At first sight, this may seem inconsistent with the results of [12,13], which considered a single ABC to model the interfaces of wire media with either a standard dielectric or a good conductor. However, the reason why we obtained an extra ABC is quite simple. In fact, here we considered the junction of two different wire media; that is, we have a spatially dispersive material on both sides of the interface (a double-sided wire junction). Quite differently, the configurations considered in [12,13] consist of single-sided wire junctions, because one of the semi-spaces separated by the interface is filled with a standard local material. Therefore, while a single ABC is sufficient to describe the electrodynamics of single-sided wire junctions, the general case of a double-sided junction requires two ABCs due to the increased number of degrees of freedom.
The very general form of the boundary conditions (9) and (12) poses the question of how the quantities ϕ z , and Q can be found in particular problems. Although giving a recipe that would be applicable to all imaginable kinds of wire junctions and terminations is clearly impossible, we would at least like to outline some guidelines. The key point is in understanding that conditions (9) and (12) are quasi-static conditions and therefore can be obtained by solving a pair of related electrostatic and magnetostatic problems. Let us illustrate it with an example junction shown in figure 2. This junction is different from what is shown in figure 1. In figure 2 we have added metallic patches at the points of wire junctions (we will consider such a configuration in more detail in the next section). Both the patch and the difference in the wire radii introduce irregularities in the charge and the current distributions close to the junction.
To calculate the additional charge Q (mostly situated on the patch) one considers a part of the junction at −l 1 z l 2 (l 1,2 must be large enough as discussed above) assuming that it is under a fixed potential ϕ 1 = ϕ 2 = V with respect to the remaining structure. This electrostatic problem is illustrated in figure 2, where the fixed potential is produced by the voltage gaps V . The gaps should be well outside of the region −l 1 z l 2 in order not to contribute to the non-uniformity of the field in the vicinity of the junction. The charge Q is then found as where Q is the total charge situated on the wire and the patch in the region −l 1 z l 2 (the region inside the dashed rectangle in figure 2). The effective capacitance associated with the junction is found as C 0 = Q/V . The term from condition (9) is found by solving a dual magnetostatic problem, in which the voltage generators V are replaced by a pair of equal current sources of the amplitude I 1 = I 2 = I . The non-uniformity contribution to the magnetic flux can be, therefore, found as = lim where is the total magnetic flux through the area −a/2 x 0, −l 1 z l 2 . Thus, the effective series inductance of the junction is L 0 = /I . Lastly, the term ϕ z is mostly due to the non-uniformity in the conductivity of the structure close to the junction and can be found from the quasi-static current-driven problem (the same as the magnetostatic problem above, but with ω > 0 and l 1,2 /λ 1) as where U is the total voltage drop on the segment −l 1 z l 2 that includes the junction. In the cases we consider in the following sections, ϕ z is due to lumped impedance insertions. For such insertions (lumped loads in technical terms), it is usually not necessary to solve for (15), as the form of ϕ z can be guessed directly from the topology of the loading circuits. Based on the above observations, the junctions we consider in this paper may be represented by an equivalent circuit depicted in figure 3. This circuit is analogous to the models of microstrip line junctions, or junctions of coaxial cables, etc [20]; therefore an alternative way of identifying the equivalent elements could be by adapting the known formulae of the microwave theory to the structures of our interest. However, this should be done with care as normally the microwave transmission lines are excited so that the two conductors forming the line are oppositely charged (the propagating mode is the odd mode), while in the structures we consider the neighboring wires are equally charged in the static limit (which corresponds to an even mode).
Yet another way of identifying the equivalent parameters is the so-called curve fitting, in which the elements of a given equivalent network are found by comparing the results of fullwave numerical simulations with the predictions of an analytical theory. We use this approach in section 7 to estimate the values of the parasitic inductance and capacitance of a bulk impedance insertion connected to a ground plane. However, such a purely numerical approach has an obvious drawback: it does not allow for analytical modeling of the effective parameters of the junction. In contrast, with the physical quasi-static approach outlined above, it is possible in many cases to obtain closed-form analytical approximations of the effective parameters.

The junction of uniaxial wire media with patches at the interface
One case of practical interest, which illustrates the application of the boundary conditions derived in the previous section, is the junction of two uniaxial wire media with metallic patches at the interface. We assume the wires to be thin so that r 1,2 a, and the host material to be the same on both sides of the interface. At the interface the wires are attached to metallic patches of width w < a such that the gap between the neighboring patches is small: d = a − w a. In [19] the following expression for the capacitance of a patch in a regular array of patches has 9 been derived: C 0 = 2π ε 0 w/log[sec(π d/2a)]. Because this capacitance is large, C 0 aC 1,2 , the additional charge Q associated with the junction is situated mostly on the patch. In other words, the microscopic current I z (z) in the wires practically coincides with the macroscopic currents I 1,2 (z) on both sides of the interface.
Next, because the wires are thin and the patch grid is symmetric with respect to the plane z = 0 it can be shown that = 0 and also ϕ z = 0, so that ϕ 1 (0 − ) = ϕ 2 (0 + ) = [ϕ 1 (0 − ) + ϕ 2 (0 + )]/2. Expressing Q through the known capacitance and the potential, we obtain the following pair of ABCs: It can be seen that if d → 0 (the gap between the patches closes), then C 0 → ∞, and the above conditions reduce to ϕ 1 (0 − ) = ϕ 2 (0 + ) = 0. It can also be seen that in this limit the wire currents on both sides of the interface become independent. The boundary condition in the limit d → 0 can be reformulated in terms of the charge density q 1,2 (z) = C 1,2 ϕ 1,2 (z), in which case it requires the (macroscopic) charge on wires to vanish at the patch surface, consistent with the theory of [13]. However, when d > 0, the macroscopic charge q(z) does not completely vanish at the plane z = 0, in contrast to the microscopic surface charge density, which vanishes exactly at the points where the wires meet the patch.
The obtained conditions are in perfect agreement with the quasi-static model of a uniaxial WM periodically loaded with patches [19]. Indeed, this model holds for the case when the period of the loading (the separation of patches on the same wire) h λ. In this case, the current distribution in wire segments that connect the patches is approximately linear and consequently the charge distribution (as well as the additional potential) on the same segments is approximately uniform. Considering a single period −h/2 z h/2 with a patch located at z = 0, we can write, with the help of (4) and (17), where I z (z) is a piecewise linear function representing the current in the wire segments and ϕ(0) = ϕ(0 − ) = ϕ(0 + ) is the potential at the patch. From (18)- (20), we obtain which after introducing a smoothly varying macroscopic current will result in equation (5) with the wire effective capacitance per unit length C eff = C + C 0 / h. This is exactly the same result as that obtained using our quasi-static model [19]. It is worth noting, however, that using the ABCs derived in this section, it is possible to model wire media with a loading period (h) of the order of a wavelength or even larger.

The conditions at the interface between a uniaxial wire medium (WM) and air
The case when the wires end abruptly at the plane z = 0 (a single-sided wire junction) is essentially the limiting case of the junction shown in figure 1 with r 2 → 0. In this limit, We are interested in the case when r 1 ≡ r 0 a. It is known that the behavior of the surface charge density close to a wire tip strongly depends on the geometry of the tip. For instance, if a wire end can be approximated with a very elongated ellipsoid, it can be shown that the linear charge density q(z) on such a wire does not have a singularity at the wire tip and behaves smoothly in the limit λ → ∞. Indeed, such a wire (approximated by a conducting ellipsoid with the main half-axes a e = b e = r 0 , c e r 0 along the x-, yand z-axis, respectively) is characterized by the surface charge density [21] σ ∝ (1/r 0 )/ (x 2 + y 2 )/r 4 0 + z 2 /c 4 e ∼ 1/ 1 − z 2 /c 2 e , when c e r 0 , when kept under a uniform potential. Because the linear charge density is q = σ dS/dz, where the surface element dS (a thin ring on the surface of the ellipsoid) is dS ∼ 2πr 0 dz 1 − z 2 /c 2 e , c e r 0 , one can see that indeed q(z) does not have a singularity (varies smoothly) at the ends of the wire at z = ±c e . Therefore, in the case of thin wires modeled by ellipsoids, there is no additional contribution due to strong non-uniformity in the charge distribution close to the junction, Q = 0, and the macroscopic current and charge on the wires coincide with the microscopic ones up to the point z = 0. Thus, from (12) the boundary condition for the current in this case is simply because I z (0) = 0. Equation (22) is a well-known condition that was used successfully by different authors [12]- [14], [22,23]. On the other hand, it is known that abruptly cut cylindrical wires have a singularity in the charge distribution close to the cut. This singularity results in an additional charge sitting at the wire tip and thus Q = 0 in this case. For instance, in [24] the surface charge density on a long, charged conducting tube was approximated by where l is the half-length of the tube. Interpreting this result in our terms, we identify the first term c 0 as the average surface charge density on a wire and the second term proportional to c 1 as the strongly non-uniform component of the same density that results in Q = 0. Nevertheless, it can be shown [24] that (c 1 /c 0 ) ∼ 1/log 2 (16l/r 0 ) → 0 when r 0 → 0. Therefore, for sufficiently thin wires the correction due to a charge at the end of wires may be neglected, and the ABC (22) may be used with an accuracy acceptable for most practical applications [12]. As already discussed in section 2, the electrodynamics of a single-sided wire junction is described by a single ABC and thus equation (22) is sufficient to fully characterize the case where the thin wires abruptly end at the junction, consistent with the results of [12]. For wires of a significant thickness, the ABC (22) has to be modified according to (12). Specifically, taking into account the non-uniform component of the charge situated close to the 11 wire tip, we may write I 1 (0 − ) = jω Q = jωC tip ϕ(0 − ), where C tip is the capacitance associated with the wire tip. On the other hand, ϕ(0 − ) = −(jωC) −1 (∂ I 1 /∂z)| z=0 − (equation (5)). Thus, for abruptly cut thick wires, this ABC can be written in the form where χ = C tip /C. The tip capacitance C tip can be calculated in the quasi-static approximation provided a certain model for the wire tips is available (like, for example, a conducting tube model mentioned above). For thin wires, χ → 0, and one obtains (22). Condition (24) obviously applies also to the case of a single-sided wire junction terminated with patches. In this case, the tip capacitance is replaced by the capacitance of the patch: C tip = C 0 .

Further generalizations
The ABCs derived in the previous sections can be further generalized. Let us, for example, consider the case when thin wires at a junction are connected through a resistive, or in general, any impedance insertion Z load . Such an impedance will produce an additional voltage drop at the junction. Therefore, when calculating the circulation of the electric field over the path depicted in figure 1, one will obtain an additional term I z (0)Z load on the right-hand side of (6). We have also seen above that the current at the junction of two thin wires must be continuous: Hence, the pertinent boundary conditions are As can be seen, these ABCs are dual with respect to equations (16) and (17). This is because a metallic patch at a junction acts in a similar way to a parallel load in a transmission line analogy, while a series insertion is a series load in the same analogy. When adding bulk loads to the wires of finite thickness, one may expect the current and the charge distributions in the vicinity of a load to be affected by such an insertion. This introduces non-uniformities in the current and charge distributions in the same manner as in a junction of two wires of different radii. These effects can be taken into account by correction terms ϕ z and . In many cases, however, these correction terms can be included directly into the load impedance Z load . For example, if inserting a bulk load into a wire results in a parasitic capacitance C par at the point of the insertion, then this capacitance should be considered connected in parallel with the load. Likewise, the correction can be described by a parasitic inductance, L par , connected in series with the load. Thus, the impedance Z load that appears in the boundary condition (25) must be taken to be equal to Z load,eff = jωL par + 1 jωC par + (1/Z load ) .
An even more general case is when the wires on the both sides of the interface are terminated with distinct impedance loads Z load1,2 that are attached to the same metallic patch.
Denoting the potential of the patch as V 0 and taking into account the charges collected on the patch and the voltage drops across the loads, we may write Eliminating V 0 from the above relations, we obtain the following generalized ABCs: It is readily seen that these conditions reduce to (16) and (17) when Z load1,2 = 0 and to (25) and (26) when C 0 → 0 and Z load = Z load1 + Z load2 . They also reduce to simple continuity conditions for the current and potential when C 0 → 0 and Z load1,2 = 0. Finally, the generalized ABC for a single-sided junction with wires in the half-space z > 0 (z < 0) can be obtained from (31) and (32) by adding (31) to (32) (subtracting (31) from (32)) and letting I 1 (0 − ) = 0 (I 2 (0 + ) = 0) in the resulting equation. For example, for the case of wires in the half-space z > 0, the ABC reads When the gap between patches closes and C 0 → ∞, we obtain an ABC for the wires connected to a perfectly conducting ground plane through a generic lumped load Z load : The same ABC applies to wires connected to a non-ideally conducting ground plane or to an impedance surface. In this case Z load represents the contact impedance of a wire of a given cross-section connected to a non-ideally conducting ground. We note that because the considered interface is associated with a single-sided junction, a single ABC of the form (33) or (34) is sufficient to characterize the electrodynamics of the problem, similar to the case considered in section 4.

ABCs written in terms of electric and magnetic fields
The ABCs derived in the previous sections are written in terms of averaged wire currents and potentials. We believe that this form is the most natural one and allows for easy modification of the ABC when the physical conditions at the junction change. Nevertheless, it is also relevant to express the ABCs directly in terms of E the and H fields, so that a scattering problem can be solved in a standard way by matching the modal fields (associated with plane waves) on both sides of the interface. In this section, we give formulae that relate the (macroscopic) current and the (macroscopic) additional potential with the vector fields. Using these formulae one can easily write the ABCs derived in previous sections in terms of electromagnetic fields.
We are interested in the boundary conditions at a planar interface z = 0. After a Fourier transform in the x y-plane, the operator ∇ in equations (1) and (2) becomes ∇ = −jk t +ẑ(∂/∂z), where k t is the wave-vector component in the x y-plane. Then, using the Fourier-transformed 13 Maxwell equations, the z-component of the current density and its partial derivative with respect to z can be expressed as where k 2 0 = ω 2 ε 0 µ 0 and E t and H t denote the vector field components tangential to the interface. As can be seen, the current is non-zero only for waves that have an electric field component along k t and the magnetic field component orthogonal to k t , i.e. for P-polarized (TM) waves. In particular, it can be seen that the ABCs are only relevant when k t is nontrivial. Finally, the wire current is expressed as I z = A cell J z and the additional potential as ϕ = −(A cell /jωC)(∂ J z /∂z) (equations (1)-(5)).
As an example, we give a complete set of boundary conditions written in terms of the vector fields for an interface of two identical uniaxial wire media connected to metallic patches at z = 0. For this case conditions (16) and (17) can be rewritten in terms of the current densities Thus, a complete set of boundary conditions written in terms of the tangential components of the electric and magnetic fields reads in this case (assuming that the host material is uniform across the junction) where Y g = jω(2ε 0 w/π )log[csc(π d/2a)] is the effective admittance of the patch grid [19,25], The relations (39) and (40) are the standard boundary conditions for the tangential components of the electric and magnetic fields on two sides of a thin grid of admittance Y g . On the other hand, equations (41) and (42) are the ABCs derived from (37) and (38) with the help of (35) and (36).
Note that the expressions on the right-hand side of (40) and (42) are continuous across the interface. When there are no patches (C 0 → 0, |Y g /C 0 | < ∞) one can see that the ABCs simply require the continuity of the first and second normal derivatives of the tangential electric field (for P-polarized waves).

WM periodically loaded with patches
In this section, we calculate the band structure of a uniaxial WM with thin wires (r 0 a) periodically loaded with patches, using the novel boundary conditions (39)-(42). We have considered this problem in a quasi-static approximation earlier [19], but that previous study was limited to the case of a small loading period: h λ. Here we will present a solution that is free from such a restriction.
The geometry of the structure is shown in figure 4. We are interested only in the propagation of the P-polarized (TM) waves; for simplicity we also assume that the loss impedance of the wires is negligible: Z w = 0. In the chosen coordinate system, E t (z) = E x (z)x, H t (z) = H y (z)ŷ and k t = k xx . Then, the fields in a single period of the structure located at 0 < z < h can be decomposed as [19] H y (z) where η 0 = √ µ 0 /ε 0 , γ TM is the complex propagation factor of the TM mode in a uniaxial WM, γ TM = k 2 p + k 2 x − k 2 0 , with k p = a −1 √ µ 0 /L = a −1 2π/ log[a 2 /4r 0 (a − r 0 )] (the plasma wavenumber), γ TEM = jk 0 is the complex propagation factor of the TEM mode, and A ± TM and B ± TEM are unknown complex amplitudes of the waves propagating in a single period in two opposite directions with respect to the z-axis.
The eigenwaves must be Bloch-periodic and thus where γ z is the (unknown) complex propagation factor of an eigenwave. Using (45) together with the boundary conditions (39)-(42), we relate the fields in two neighboring periods and obtain a system of four equations: Once the expressions (43) and (44) are substituted into the above system, one obtains a homogeneous system of linear equations for the unknowns A ± TM and B ± TEM . The dispersion branches are then found by equating the determinant of this system to zero and solving for The results of this calculation are presented in figure 5. In this figure, in panels (a) and (b), we compare the dispersion curves obtained with the model of this section and the results obtained with full-wave numerical simulations calculated with a commercial solver [26] for two cases of different patch widths: (i) w = 0.5a and (ii) w = 0.9a. The wire radius in both cases is r 0 = 0.05a and the loading period has been set equal to h = a. We consider the Ppolarized waves propagating in the xoz plane along the directions of θ = 0 • , 30 • and 60 • with respect to the wires (tan θ = k x /k z ). The horizontal axis of figure 5 represents ka = k 2 x + k 2 z a. The agreement in all cases is rather good. As discussed in our previous work [19], the quasistatic model overestimates the plasma frequency by about 3% when r 0 = 0.05a, which explains the small difference in the dispersion of the plasmon mode at low ka. When compared to the simple quasi-static model of [19], the model based on the ABCs developed in this paper correctly predicts the appearance of a bandgap when the patches are large enough (w = 0.9a). For the waves propagating along the wires (θ = 0 • ) the present model also agrees much better with the full-wave simulations (particularly, in the dispersion of the plasmon mode) than the homogenization model based on the ABC of [13] (see also appendix B of [19]).
It should also be noted that the parameters used in this example of (figures 5(a) and (b)) are rather extreme for a quasi-static model. For instance, the wire period a may be comparable with λ/2, and the loading period is set to h = a, while one might expect the model to be valid only for a λ and h a. Nevertheless, even with these extreme parameters we observe good agreement with the results of the full-wave simulations. This indicates that the quasi-static model of [19], as well as the ABCs developed in this paper, describes very well the real physical processes taking place in wire media.
To further illustrate this, figure 5(c) represents the results of calculations for thinner wires with radius r 0 = 0.01a loaded with patches with moderate width: w = 0.8a. As expected, the agreement between the full-wave simulations and the quasi-static ABC-based model is better for wire media with lower plasma frequencies.
The numerical examples show that the accuracy of the developed approximate model is within few per cent when k p a 1.5, ka 1.5.

WM connected through lumped loads to a ground plane
In the second numerical example, we consider a WM slab defined by the region −h < z < 0. The region z > 0 is filled with air, and the WM slab is backed by a (perfectly conducting) ground plane placed at z = −h (see the inset of figure 6). It is assumed that the metallic wires are connected to the ground plane through a lumped load Z load . We are interested in studying the reflection of P-polarized waves by the grounded slab. It is assumed that the incoming plane wave propagates in the xoz plane and illuminates the slab along the direction θ inc , measured with  . (a, b) The wire radius is r 0 = 0.05a; the loading period is h = a; the patch width is w = 0.5a (a) and w = 0.9a (b). (c) r 0 = 0.01a, w = 0.8a, ka < 1.7. The lines represent the result of the homogenization model based on the ABCs developed in this work, and the symbols correspond to the full-wave simulations obtained with CST Microwave Studio [26]. respect to the normal direction. Thus, the tangential fields in the region z > 0 can be written as where ρ is the reflection coefficient, γ 0 = k 2 x − ω 2 /c 2 , k x = (ω/c)sin θ inc , and c = 1/ √ ε a µ 0 and η a = √ µ 0 /ε a are the speed of light in vacuum and the intrinsic impedance, respectively (the permittivity of the WM host, ε 0 , is not necessarily the same as the permittivity of the air region ε a ). As in the previous subsection, inside the WM slab, −h < z < 0, the electromagnetic fields are calculated using formulae (43) and (44). The reflection coefficient can be obtained by matching the fields at the air interface, exactly in the same manner as discussed in [13], and in addition by considering the following boundary conditions at the ground plane (z = −h): where J z and ∂ J z /∂z are defined as in (35) In figure 6 we plot the phase of the reflection coefficient for different loads as a function of the normalized frequency for the cases θ inc = 60 • , ε 0 /ε a = 3.0, h = a, r 0 = 0.025a and a = 2 mm. We considered both inductive loads (L = 0.2 and 0.6 nH) and capacitive loads (C = 0.1 and 0.5 pF). We have also considered the limit cases of a short circuit (SC) and an open circuit (OC). It is assumed that the load is connected to the ground plane through a gap of 0.2 mm. By comparing the results of the analytical model with the results of full-wave simulations performed with CST Microwave Studio [26], we estimated that such a gap is characterized by parasitic inductance L par ≈ 0.11 nH and parasitic capacitance C par ≈ 0.003 pF. Figure 6 reveals good agreement between the analytical model (solid lines) and the numerical results (dashed lines) over the considered frequency range. The results show that the reflection characteristic depends strongly on the value of the load. The points where the phase crosses 0 • and −360 • correspond to the resonances where the metamaterial slab effectively behaves as a high-impedance surface, mimicking in part the response of a perfect magnetic conductor [27]. The amplitude of the reflection coefficient is equal to unity (not shown) because, for simplicity, the materials were assumed lossless.
Typically, as the reactance of the load becomes more positive, the frequency where the phase crosses zero decreases. Notwithstanding this property, it can be seen in figure 6 that the curve associated with the SC is shifted to lower frequencies with respect to the curve associated with C = 0.5 pF. The reason for this curious property is the effect of the parasitic inductance of the gap, which is relatively large for the case of a lumped load but vanishes for the case of an SC (supposing that the SC is attained by extending the pins to the ground plane, as considered in our model).

Conclusions
In this paper, we have developed a rather simple but powerful approach to the problem of additional boundary conditions in wire media. It was discussed that, despite the inherent nonlocality of wire media, it is possible to formulate a local effective medium theory by introducing two additional state variables. The approach is based on a quasi-static model of wire media that introduces two additional parameters: the wire current and the wire potential [19]. The conditions at single-and double-sided junctions of wire media with or without metallic patches or loads at an interface are then formulated in terms of these quantities. The developed approach is in a sense similar to the transmission line models used, for instance, in microwave theory. Because of this similarity, it seems rather easy to adjust the general form of the ABCs presented in this paper to a wide range of configurations of great practical interest. Some of these configurations have already been considered in this paper.
We have also shown that the obtained conditions could be reformulated in terms of the electric and magnetic field components tangential to an interface. Using such ABCs, we have calculated the dispersion of the low-frequency P-polarized eigenmodes in a uniaxial WM periodically loaded with metallic patches. The solution presented in this paper is not limited by the assumption that the loading period is small compared to the wavelength. This solution has been compared with a full-wave simulation and found to be in good agreement. In addition, we have also applied the developed ABCs to study the scattering of a plane wave by a WM slab connected to a ground plane through reactive lumped loads, demonstrating that the reflection characteristic is strongly dependent on the loads.