Identification of Compact Equivalent Circuit Model for Metamaterial Structures

In this work, a simple and efficient circuit modeling of metamaterial structures, providing a compact circuit able to describe accurately the device over a large bandwidth of operation, is proposed. The equivalent circuit model is obtained by the identification process of the load in terms of shunt branches constituted by reactive elements that can be both positive and negative. The circuit model is validated by analyzing the split-ring resonator (SRR) structure. The presence of negative elements in the non-Foster load is transformed into positive reactive elements by converting the load from shunt into series. Unlike the T or $\Pi $ circuit models, using this approach a circuit model can be constructed directly from the scattering parameters and valid for any circuit topologies.

Abstract-In this work, a simple and efficient circuit modeling of metamaterial structures, providing a compact circuit able to describe accurately the device over a large bandwidth of operation, is proposed. The equivalent circuit model is obtained by the identification process of the load in terms of shunt branches constituted by reactive elements that can be both positive and negative. The circuit model is validated by analyzing the split-ring resonator (SRR) structure. The presence of negative elements in the non-Foster load is transformed into positive reactive elements by converting the load from shunt into series. Unlike the T or circuit models, using this approach a circuit model can be constructed directly from the scattering parameters and valid for any circuit topologies.

I. INTRODUCTION
M ETAMATERIALS are artificially designed resonant structures which exhibit negative permittivity and/or negative permeability, realized with periodic metallic inclusions in a dielectric media [1]. These artificial materials are of great interest for their peculiar effective material properties and its usefulness in scalability and integration of microwave/mm-wave components. They became popular since the first demonstration of negative dielectric permittivity using a periodically arranged thin metallic wires as inclusions [2] and negative magnetic permeability using a resonating metallic ring structures [3], [4]. Accurate equivalent lumped element models are useful to predict the behavior of the metamaterial structures in a simple and computationally efficient way. Equivalent circuit model permits to establish straightforward relationships with physical parameters (i.e., geometrical dimensions, etc.) of metamaterials' structure and its frequencydependent reflection/transmission characteristics. The most studied resonant structure for the creation of metamaterial is the split-ring resonator (SRR) that is a ring-shaped subwavelength metallic structure with a split gap providing negative permeability near and above the resonance frequency. SRR and its variants are the widely studied metamaterial structures for the frequency range stretches from microwaves to terahertz (THz) optics [5], [6], [7], [8]. In a simplest form, SRRs can be modeled as resonant LC circuit [9], L accounting for the line inductance of the metallic ring and C coming from the dielectric split gap. At the resonance frequency, the magnetic permeability goes to negative and will create a stopband in the transmission spectra. Several works in the past were dedicated to the analytical/numerical modeling of SRRs to understand its resonance behavior [10], [11]. Researchers have already shown that SRRs exhibit not only magnetic resonance but also an electric resonance depending on the orientation of the SRRs to the incident field [12], [13], [14]. Some of the reported equivalent circuit models for SRRs use a quasi-static approach as it is acceptable given that the metamaterial structures are smaller than the free-space wavelength [15], [16], [17], but these methods suffer from their own limitations and require tedious calculations and evaluations before applying them into the design process. Some other works [18], [19] also reported similar quasi-static analysis to predict the resonance frequency of SRRs embedded in microstrip, or coplanar structures. The calculation of the capacitance per unit length in those analyses requires an accurate value of the effective dielectric permittivity of the substrate. For a multilayer high-frequency technology stack, it requires additional effort on calculating the effective permittivity considering all the layers of the substrate while using the simple formula that takes into account only one substrate layer will produce substantial error in predicting the resonance frequency. In addition, it should be noted that different kinds of excitations produce different spatial configurations of the electromagnetic fields impinging on the SRRs and, in this case, it may not be so simple to apply similar analysis to the prediction of the frequency of resonance. Moreover, some of those models are useful in predicting only the resonance frequency of the metamaterial structure but not the complete reflection/transmission behavior for the whole frequency band. Sometimes, due to the distributed character of the structures being modeled, they could yield physically unsuitable lumped elements with no physical correspondence to the modeled metamaterials. For practical design purposes, we need precision over a large bandwidth of operation other than just predicting only the resonance frequency.
In this work, we propose a simple and efficient circuit modeling of lossless metamaterial structures yielding a compact circuit to address high accuracy over a large bandwidth. The objective is mainly devoted to develop an equivalent network valid for a large bandwidth based on the "minimal" model defined by Marcuvitz [20] and Montgomery et al. [21] for microwave discontinuities that is constituted by a shunt load embedded in two transmission lines. This circuit is "minimal" because it contains only three electrical quantities, just as the number of the independent parameters of the scattering matrix S of the lossless device under study. An improvement of the model is obtained by the identification process of the load in terms of shunt branches, constituted by one branch with only a capacitance and a number of resonant branches made by all positive or all the negative inductances and capacitances. Although it is correct to use such negative elements from an identification point of view, we propose to substitute these negative reactive elements, producing "non-Foster" loads, with positive Foster reactive elements by means of a combination of transmission lines and a shunt-to-series transformation.
It should be stressed that the non-Foster elements obtained in the identification process proposed in this article are due to the "compression" of the electrical complexity of the lossless device in terms of an equivalent circuit made only with a "minimal" representation, based on three electrical parameters. Moreover, these non-Foster elements cannot be seen as standalone elements, as the previous one obtained by active devices, but they have to be considered as a constituting block of the overall equivalent circuit of the entire device that still has a Foster behavior at the input-output ports, being a passive lossless device.
This proposed approach is not limited only to SRR case but can be applied to any kind of lossless metamaterials' design to attain a compact circuit.

II. THEORY
The evaluation of an equivalent circuit for the SRR can be performed with various approaches, as discussed in the introduction, but they mainly refer to the lossless case or the quasi-static one. The hypothesis of lossless materials may not be correct in principle, because any realized structure is made by lossy dielectric and lossy metal but considering at first a lossless structure permits to obtain the general properties of the equivalent circuit related to the SRR geometric parameters and to the SRR electromagnetic field distribution.
The unit cell of the analyzed structure is shown in Fig. 1(a) and (b). It is constituted by an SRR, placed on the x y plane, with N rings (two in the figure) embedded in the Si 3 N 4 layer (height h Si 3 N 4 = 4 µm) and fabricated on a Si substrate (height h Si = 650 µm) covered with a SiO 2 (height h SiO 2 = 300nm) layer. The external ring has a width w, length a, and thickness t = 1 µm; spacing between the concentric rings is s and the width of the gap in each ring is g. A multilayer stack is taken into account here for a practical consideration, because the results discussed in this article are the preliminary step of a more global project that  would joint together the presence of metamaterials, antennas, and active device at different planes in the same technological stack [31], [32].
Another important aspect that must be considered in deriving an equivalent circuit is the kind of excitation of the device. In this work, we suppose that a plane wave with E x and H y component is traveling in the z-direction, impinging from the air on a periodic distribution of SRR's along the x-direction, with periodicity P x . Similarly, a periodic distribution of the SRR's in a vertical stack with periodicity P y = h Si 3 N 4 + h SiO 2 + h Si = 654.3 µm is considered. Under this hypothesis, we can define a unit cell, with dimensions P x × P y and use the Floquet modes' approach (with no phase shift between the walls of the unit cell) that reduces the simulations effort. The electromagnetic software CST Microwave Studio has been used to evaluate the S matrix of the unit cell. Such an arrangement permits to analyze the stopband properties of the SRR in the presence of an exciting electric field, parallel to the plane of the SRR. This arrangement could be used, for example, to suppress the surface mode arising in a stack planar structure.
The extraction of an equivalent circuit could also be performed in a quasi-static approach, but it is very difficult to obtain a reasonable and simple equivalent circuit that could be used to predict all the resonances at high frequencies with an acceptable precision. Moreover, these models cannot be considered reliable for parametric analysis of the devices, except for particular cases.
In the hypothesis of lossless material, the simplest way to extract an equivalent circuit is the method presented in the Marcuvitz [20] or Montgomery et al. [21] handbooks that, starting from the knowledge of the scattering parameters, allows to obtain the equivalent circuit shown in Fig. 2. The proposed equivalent circuit takes into account the discontinuities at the air-dielectric interfaces at z = 0 and z = L r that can be represented by the connection of two different transmission lines, being the first relative to the wave propagation in air and the second to the propagation in a dielectric medium with effective dielectric constant ε e r and wave impedance η r = (η 0 / √ ε e r ), where η 0 is the air wave impedance. This circuit is obtained in the hypothesis that the evaluation of ε e r of the multilayer structure shown in Fig. 1(b) is possible by means of the effective dielectric constant method [33] or by the evaluation of the modes properties with CST or other full-wave electromagnetic simulation software. In fact, in this case, the scattering matrix S air obtained by CST with the Floquet ports placed at z = −L p , z = L r + L p (black dashed box in Fig. 2) can be de-embedded in phase by adding two negative lines of length θ 0 = −k 0 L p (this evaluation can also be done directly by CST) obtaining the matrix S 0 (blue dashed box in Fig. 2), which can be used to obtain the scattering matrix "seen" at the dielectric ports S r (red dashed box in Fig. 2), by denormalizing S 0 with respect to η 0 and normalizing with respect to η r .
It should be noted that the circuit representing S r (red box) is minimal, being composed by only three unknown parameters (i.e., the electrical lengths of the two external transmission lines, θ r 1 , θ r 2 , and the shunt load B) which correspond to the number of independent scattering parameters, for example, S r 11 , ϕ r 11 , and ϕ r 22 . The circuit of Fig. 2 is efficient, because the reflection/transmission properties of the device are directly and only related to the susceptance B that can be simply evaluated by [34] B = ± 2 S r 11 η r S r

12
(1) The presence of the external transmission lines θ r 1 , θ r 2 ensures the correct values of the phase of the global scattering parameters of the equivalent circuit with respect to S r . Therefore, the equivalent circuit shown in Fig. 2 is very effective being able to mix the effect of the reflection at the air-dielectric discontinuities and the effect of the rings in the dielectric layers. It is noteworthy to mention that this equivalent circuit can always be defined for any reciprocal microwave structure because it is directly related to the scattering parameters by (1)-(3) that can be obtained numerically or experimentally for any device. On the contrary, T or equivalent circuit coming from the Z or Y matrix is not defined for some circuit topologies, as, for example, it occurs for a transformer that cannot be expressed in terms of Z or Y but only with T and S matrix. Moreover, the scattering properties of the device in those representations are related to the series and shunt reactances, and a quick evaluation of the reflection/transmission properties of the equivalent circuit is not possible due to the combination of three reactances.
It is evident that the reflection/transmission properties of the device, the SRR in our case, are included only in the shunt load and it can be frequency-dependent. For very simple structures, such as thin inductive or capacitive windows in waveguides, the frequency dependence of the susceptance B corresponds exactly to a pure inductive or capacitive susceptance [20], [21]. For more complex shapes of the aperture between the two sides of the waveguide, the shunt load can be represented by resonant LC circuit. Hence, the global susceptance B can be represented by some combinations of inductance and capacitance matching with the frequency behavior of B. This is the problem known as identification, i.e., how to identify the best representation of B in terms of lumped inductances or capacitances in a proper network combination.
The same B identification holds also for the SRR but, in this case, the complexity of the structure that depends on many geometrical parameters does not allow a simple evaluation of inductances and capacitances in a particular combination to be used to represent the frequency behavior of B. To this aim, the identification of B can be done by means of the classical approach based on the pole-zero identification starting from the frequency response [35]. The typical implementation of such technique is based on an optimization process (least squares, maximum likelihood estimation, vector fitting, or others. . . ) which try to evaluate the coefficients of a transfer function made by the ratio of two polynomials, which approximates the frequency response of B. The classical evaluation imposes some constraints on the polynomial coefficients to obtain a positive rational function (PRF) and it can be synthesized by positive inductances and capacitances. In this case, the frequency behavior of the obtained network satisfies the condition (d B( f )/d f ) > 0. This kind of network representing the shunt load B( f ) is defined as a Foster network. Removing the constraints of the classical approach, non-Foster networks can be included in the identification process for the lossless case under study. The main difference between a non-Foster and a Foster circuit is related to the slope of (d B( f )/d f ) that has to be positive for Foster and negative for non-Foster susceptance (or reactance). Consequently, the identification of the corresponding equivalent circuit is based on the presence of inductances and capacitances both positive or both negative for the two cases, respectively. Moreover, a non-Foster circuit can also be identified with positive inductance/capacitance but with positive and negative resistors [36], [37], [38].
Hence, the proposed approach to perform the identification of the shunt load B in the equivalent network is to include both the Foster and non-Foster networks with the only requirements that the identification with positive and negative resistors and positive inductance/capacitance for the non-Foster part is not taken into account for the lossless SRR case because, to the authors' opinion, the presence of positive and negative resistors for a lossless device can be considered as questionable.
The presence of non-Foster reactances/susceptances in a lossless circuit could be surprising. Actually, there are some cases where the identification process yields to the presence of elements with the non-Foster behavior, even if the actual device shows a Foster behavior at its input and output ports. For example, the representation of a simple transmission line with a two-port circuit requires the presence of such elements. In fact, we know that the scattering matrix of a transmission line is simply S 11 = 0, S 12 = e − jθ , where θ = 2π f √ LC L is the electrical length, f is the frequency, L, C are the inductance and capacitance per unit length, respectively, L is the physical length, and Z 0 = √ L/C is its characteristic impedance. From the S matrix, we can obtain the Z matrix that can be represented in terms of a T equivalent circuit, made by a series reactance X line s = −Z 0 cot(θ) and a shunt reactance X line p = −Z 0 csc(θ). If we plot X line s and X line p versus frequency f , we can observe that X line s corresponds to a classic Foster reactance, characterized by (d X line has a behavior corresponding to a Foster case for 0 ≤ θ ≤ π/2 and 3π/2 ≤ θ ≤ 2π and to a non-Foster case for π/2 ≤ θ ≤ 3π/2. X line s can be identified with a network containing Foster elements while X line p with a network that must contain non-Foster elements (or positive RLC elements with additional negative R).
We can explain this different behavior for X line p in terms of the circuit representation we are considering. In fact, the actual circuit representation of a transmission lines should be a ladder network with infinite elements L Z and C Z , with Z , length of a unit cell, approaching to zero. To obtain a more compact and easy-to-use equivalent circuit of the transmission line, the infinite ladder network, made by positive inductances/capacitances (with the Foster behavior), can be replaced by the T network that "compresses" the infinite ladder network in a simple two-port representation made by X line s and X line p . This "compression" yields to the presence of a combination of the Foster and non-Foster elements in some frequency band, while the overall behavior at the input-output ports is always of Foster type. On the other hand, as stated by Marcuvitz [20] in his handbook, we can define, for the same device under study, a number of different equivalent circuits. In fact, ". . . an equivalent circuit of a device can represent a simple frequency dependence or the minimization of the electrical parameters, or the effects of evanescent modes. . . . Hence, any device can be represented with various equivalent circuits, and each one represents a particular characteristic and is 'correct'." This is just the case for the transmission line for which we can find at least two different equivalent circuits, each of one with its own characteristics.
1) The T network representation implies "compactness" and the possibility to simply manipulate the cascade with other discontinuities, but with the presence of the Foster and non-Foster elements, which do not change the behavior at the input-output port of the overall transmission line that still has a Foster behavior.
2) The ladder network representation implies the presence of positive inductances/capacitances (Foster behavior) but with infinite elements in the ladder network and not a simple circuit manipulation when connected to other circuit.
It should be stressed that the non-Foster behavior of X line p in some θ range does not imply that we can actually realize a non-Foster reactance by the T equivalent circuit of the transmission line, because the only way to realize a non-Foster reactance is based on the use of active elements [22], [23], [24], [25], [26], [27], [28], [29], [30]. In fact, the T equivalent circuit must be taken as a "whole," where any elements contained in it must be seen as "internal" elements strictly connected to the other elements, i.e., the series reactance X line s for this case. Furthermore, for the same passive lossless device we could find many circuit representations containing only Foster "internal" elements (just as the ladder network with an infinite number of positive L and C) and many equivalent circuits containing also non-Foster "internal" elements, but the latter can greatly reduce the topological complexity of the obtained equivalent circuit. In this sense, the presence of "internal" elements with the non-Foster behavior (just as X line p ) can be a key to define simple equivalent circuits based on "minimal" representation.
The same "compression" effect can be found in the equivalent circuit shown in the blue dashed box in Fig. 2, represented by (1)-(3). As previously discussed, this is a "minimal" circuit representation for a complex device under study, which is represented by its two-port scattering matrix, which contains in its numerical values the complexity of the device, just as a "black" box. Hence, the global electromagnetic behavior of the device is "compressed" in a two-port representation and we have no other information about the complexity of the "black" box and which could be the "best" representation in terms of an equivalent circuit. As an example, let us think to a passband Chebichev filter realized in lossless waveguide with many cavities for which we have to define an equivalent circuit from the knowledge of its S matrix evaluated at the input and output ports. By the transformation of S in the Z or Y matrix, we can define a T or equivalent circuit which contains the global properties of the filter in the reactances/susceptances of the series or shunt branches that could also have the non-Foster behavior, as it occurs in the previously discussed case of the simple transmission line. Similarly, from the S matrix we can obtain the "minimal" circuit representation (1)-(3), where the reflection and transmission properties are "translated" only in the frequency behavior of the shunt susceptance, being the two external line lossless and useful only to the match phase condition on the global scattering parameters. Being the filter characterized by zeros of reflection and limited ripples in the passband, we can expect that the shunt susceptance B is characterized by an oscillating frequency behavior crossing zero at the zeros of reflection, yielding to the Foster and non-Foster behaviors between two consecutive zeros of reflection. Hence, whichever the equivalent circuit we choose, we can expect the presence of the Foster and non-Foster "internal" elements in the circuit representation due to the "compression" of the twoport representation. The "internal" elements, with their Foster, non-Foster, or dispersive behavior, must be seen as bricks of the whole equivalent circuit because they must represent the passive lossless device with a Foster behavior at the input-output ports and they cannot be "extracted" to realize a real non-Foster element. This consideration is applicable to all the cases analyzed in this article.
As previously discussed, once the "minimal" representation shown in Fig. 2 has been obtained by the global lossless S matrix, obtaining the frequency behavior of B from (1) to (3), the next step is to identify the shunt load B by means of a rational function, obtained by the combination of a number of branches containing one positive capacitance/inductance and resonant LC series, made by both positive (Foster type) or both negative (non-Foster type) inductance and capacitance, to reproduce the non-Foster behavior in some frequency bands due to the "compression" of the proposed circuit representation. The rational function representing the reactive lossless part including positive inductances and capacitances satisfies the conditions established by Brune for the identification of an admittance Y(s), s = σ + jω [35].
1) No poles lie in the right half-plane.
2) No zeros lie within the right half-plane.
3) Poles on the boundary are simple and have positive real residues. 4) Zeros on the boundary are simple and at them To include the presence of the negative inductances/ capacitances (non-Foster type) in the rational function that identifies B, we have to add two other conditions. 6) Poles on the boundary are simple and have negative real residues. 7) Zeros on the boundary are simple and at them (d B(ω)/dω) < 0.
It should be noted that the network proposed to identify B, based on resonant LC branches containing only positive or only negative L and C, satisfies conditions 1 and 2 as the poles are always purely imaginary, and condition 5 because the network is purely reactive. Hence, no problem arises in terms of stability of the network. The identification process can be made by minimizing the differences (least squares, maximum likelihood estimation, vector fitting, or others. . . ) between the obtained susceptance B of the equivalent network shown in Fig. 2 and the fitting rational function B fit obtained by the previous Brune's modified conditions. The obtained fitting rational function B fit will be transformed in a complex circuit by the usual technique of extraction of poles and zeros and their identification in terms of shunt or series branches with positive and negative inductances or capacitances that compose the global network of the prescribed B fit .
The last step is to try to redraw a shunt branch of inductances/capacitances eventually identified in the global network in a series branch, using the properties of the transmission lines. This transformation can be obtained by means of simple network arrangements.
To this aim, let us evaluate the scattering matrix for a generic shunt load, represented by its normalized susceptance b = Bη r The transformation from the shunt normalized susceptance b shown in Fig. 3 into the series normalized reactance x has been done by introducing two transmission lines with characteristic impedance η r and electrical length −θ bx 1 , −θ bx 2 that transform ensuring at the same time that S x is the scattering matrix of a series reactance. This condition is obtained in terms of the normalized currents i x 1 , i x 2 , as defined in Fig. 3. In fact, let us define the normalized currents and the wave amplitudes as a x The correct values of θ bx 1 , θ bx 2 that transform S b into S x , by guaranteeing that S x is the scattering matrix of a series reactance, are obtained by imposing that for any excitation coming from the two ports. Recalling that for excitation coming from ports 1 and 2, respectively. The previous equations can be cast in the following form: (13) and (14) have the following solutions: The first solution, (15), gives that is the scattering matrix of a series normalized impedance Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
while the second solution, (16), gives that is the scattering matrix of a series normalized impedance Hence, the two solutions produce a transformation of the shunt load into a series load with the difference that by the first solution (15) Similarly, if the shunt load is a non-Foster type, with (db( f )/d f ) < 0, we can obtain a series non-Foster-or Fostertype reactance, by choosing (15) or (16).
It should be noted that the analyzed transformation can also be applied in the reverse case to transform a series load into a shunt one. This can be done by adding two negative lines −θ bx 1 , −θ bx 2 to the left and right sides of each circuit (shunt or series) in Fig. 3. In so doing, a series reactance is transformed in a shunt load with similar properties as in (18)-(20) with two lines with electrical length equal to −θ bx 1 , −θ bx 2 as in (15) and (16).
Therefore, the use of transmission lines of proper electrical lengths permits to overcome the presence of negative inductances/capacitances in the global network, obtaining a more comprehensive network able to represent the property of the device under study over a large bandwidth up to very high frequencies. The proposed identification technique can be summarized as follows: given a scattering matrix S 0 representing the SRR under test excited by a plane wave, the first step is the identification of the equivalent network shown in Fig. 2 by means of the denormalization of S 0 obtaining the scattering matrix seen at the dielectric ports, S r . The second step is the identification of the electrical lengths θ r 1 , θ r 2 and the shunt susceptance B by means of (1)-(3). The third step is the identification of B by means of a rational function B fit ( f ) satisfying the modified Brune's conditions that can be found with the minimization of a functional based on the least-square differences between B and the fitting rational function B fit ( f ). The choice of the best method to minimize the differences is not analyzed in this article where we use the least-squares method. The fourth step is the transformation of B fit ( f ) in terms of series or shunt branches composed by positive and negative inductances by the classical extraction of poles and zeros. If only positive components have to be defined in the global network, the last step is the transformation of the shunt negative components in terms of two equal transmission lines embedding series positive components by means of (16) and (20). Similar equations can be used to transform negative series components into positive shunt ones between two transmission lines.

A. Closed Ring Resonator (CRR)
The proposed identification process has been applied to the analysis of CRR unit cell containing only one closed ring (without the gap) with g = 0, w = 2, P x = L r = 174, L x = 172 (hereafter, all the geometric dimensions are expressed in µm). The effective dielectric constant extracted by CST is ε e r ≈ 8. The scattering parameters' amplitudes seen from the air (S 0 ) and in the dielectric (S r ) are shown in Fig. 4(a). From S r , and (1)-(3), we can evaluate the electrical lengths θ r 1 = θ r 2 and the normalized susceptance b = Bη r of the equivalent circuit as shown in Fig. 4(b). The closed ring CRR totally reflects at about 99.5 GHz and this is due to some resonant part contained in the identified susceptance. The behavior of b [black line in Fig. 4(b)] shows a resonance that can be related to a Foster network, for which (db( f )/d f ) > 0, and synthesized by a resonant series LC, being L and C both positive. Moreover, the behavior of the lower and higher parts of the band suggests the presence of a capacitance also in shunt with the resonant part of the circuit. In fact, by applying the identification process previously discussed to b with a realizable shunt generic load made by two branches, we have obtained that b can be realized with the global shunt load shown in Fig. 5(a) where C = 7.95 fF, L F = 0.57 nH, C F = 4.49 fF. The comparison between the exact susceptance b (black line) and the fitting b fit (red line) that is approximated by the circuit with two branches is shown in Fig. 4(b) and it matches well. Some explanations are needed for the two branches' circuit whose contributions are shown in Fig. 5(b). The first branch is constituted by a capacitance (blue line) that is mainly related to the capacitor formed between the sides of the left and right unit cell conductors [39], as shown in Fig. 5(c). This can be verified by the following analysis. The proposed identification approach has been applied to the CRR with a rectangular shape by fixing the length L x = 172 and varying the conductor length L z [see Fig. 1(a)]. The identified capacitance C is almost linear with L z , Fig. 6(a), while C F shows a log behavior versus L z , Fig. 6(b), as in the following approximated models: Hence, the length L z of the side conductors affects mainly C and slightly C F . While the capacitance in the first branch is mainly related to the length of the side conductors as previously discussed, the capacitance C F in the LC branch is mainly related to the change in direction that the electric field undergoes when it crosses the dielectric and impact on the ring conductor. In fact, Fig. 6. (a) Identified capacitance C of the first branch versus the lateral conductor length L z or the front conductor width L x and its approximated model as in (21) and (23). (b) Identified capacitance C F of the second branch versus L z or L x and its approximated model as in (22) and (24). the electric field of the plane wave impinging from the air on the CRR is directed horizontally along the x-direction and it must satisfy the boundary conditions on the CRR conductors. Due to that, it must change its direction to be orthogonal to the conductors. The change in direction of the electric field is clearly shown in the simplified and not in scale sketch of the electric field and charge distribution shown in Fig. 7, obtained from the CST simulations. The reference plane for the sketch is in the plane AA' indicated in Fig. 1(b). The field intensity is proportional to the arrow strength while the legend for the charge intensity is shown inside the figure. The wave is impinging from the air (bottom) with the horizontal electric field and the presence of the horizontal conductor constrains the electric field to change its direction to satisfy the boundary conditions as shown by the field lines that are orthogonal to the conductors even if it is not evident in the sketch. It is clear that the end parts of the horizontal conductors contribute to the capacitance C of the first branch, whereas the inner part of the same conductors contributes to the capacitance C F of the second branch. To verify this hypothesis, the proposed identification approach has been applied to the CRR with rectangular shape, fixing the length L z = 172 and varying the conductor width L x [see Fig. 1(a)]. The identified capacitance of the first branch C shows a log behavior, Fig. 6(a), while the capacitance C F of the second branch is almost linear with L x , apart from a weak quadratic dependence, Fig. 6(b), as in the following approximated models: Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
Hence, the width L x of the front conductor affects mainly C F and slightly C. It is not simple to evaluate the correct relationships for C, L F , C F with respect to the geometrical parameters other than L x , L z for the analyzed plane wave excitation. Many researchers in the past have proposed some approximated formulas [40], [41], which refer to the case of oblique incidence on the plane of the CRR. In this article, we are analyzing the effect of a plane wave impinging at 90 • with respect to the CRR plane; so, those formulas are not valid to the case under study. Moreover, our interest is not to provide an exact evaluation of such relationships but the evaluation of an equivalent network suitable to describe the effects of the CRR over a large bandwidth taking into account also the electromagnetic field behaviors.
Concerning the inductance L F = 0.57 nH, it should be useful if it can be obtained by an equation similar to the one valid for the inductance of a square loop with rectangular conductors of section wt [42]. Unfortunately, this is not the case because equation was obtained in the hypothesis of constant current along the square sides while in the CRR under study, the current distribution is completely different as the magnetic field. Hence, the inductance of the loop cannot be evaluated as in [42].
This aspect can be clarified by the sketch of the surface currents on the conductors obtained by the CST simulations at the resonance frequency, f = 99.5 GHz, as shown in Fig. 8. The surface current density distribution in the top and bottom surfaces of the ring conductor is shown as blue lines with arrows (currents are very similar). The surface current density distributions in the external and internal surfaces of the ring conductor are shown as red lines with arrows. Their intensity is proportional to the strength of the lines/arrows. The blue and red circles represent positions where currents vanish. The current distribution follows a not constant behavior and it goes in opposite directions in the upper and lower parts of the ring, vanishing in the middle of the structure. The magnetic field is also shown in the same sketch and it is evident that it has a behavior related to the currents and it vanishes in the middle of the CRR. Therefore, the current is far to be constant along the ring and the related magnetic field is completely different from Fig. 9. Inductance L F versus the side conductor length L z or width L x and its approximated model, as in (25) and (26).
that of a loop with constant current, making [42] inapplicable. Moreover, the value of the inductance for a radiating circular loop is related to the frequency as discussed in [43], due to that a prediction of the value of the inductance L F is not easy to perform.
Anyway, starting from the behavior of L F versus L z or L x , as shown in Fig. 9, an approximating equation for L F can be obtained for the variation with respect to L z for L x = 172 (25) and for the variation with respect to L x for L z = 172 Hence, we can conclude that the second branch with the LC resonator is due to the loop inductance and mainly to the front capacitance of the CRR.

B. Lateral Gap SRR
The second CST simulation refers to a one-ring SRR with a gap g = 2 in the ring. As in the previous case, the scattering parameters seen from the air and in the dielectric are shown in Fig. 10(a). In this case, there are two resonances at around 98.48 and 77.06 GHz. The introduction of the gap produces a little change in the higher resonance about 1 GHz compared with the closed ring case, less than 1%, and a new resonance appears at 77.06 GHz. The highest one corresponds to the resonance of the closed ring without the gap. This is due to the fact that the gap has been placed in the position where the current vanishes at the resonance of the CRR, as shown in Fig. 8. The placement of the gap in that position permits to maintain almost the resonance of the full closed ring case. The corresponding surface current behavior for the gap case at 98.48 GHz is shown in Fig. 11(b) that is very similar to the closed ring case shown in Fig. 8. A little asymmetry appears for currents and magnetic fields and it is due to the gap that is placed in one side only.
The new resonance is mainly due to the effect of the gap between the conductors and this is confirmed by the surface current behavior at 77 GHz, Fig. 11(a), that is totally different from its behavior at 98.48 GHz, Fig. 11(b). In fact, at the lower frequency resonance, all the components of the surface current  It is interesting to compare the equivalent network parameters, θ r 1 = θ r 2 = θ r and b obtained by (1)-(3), for the CRR [closed ring, Fig. 4(b)] and the SRR [see Fig. 10(b)]. The first remark concerns the different behavior of the electrical lengths θ r for the two cases: while θ r for the closed ring case is linear versus frequency [see Fig. 4(b)], the gap case shows a sharp change in the slope at 77.06 GHz (see Fig. 10(b), black line), while it remains quite linear at 98.48 GHz. This behavior is strictly related to the evaluation of b (see Fig. 10(b), red line): the gap case shows two infinite values at the resonances but the main effect is that (db( f )/d f ) is negative around 77 GHz, while is positive around 98.48 GHz. Hence, b can be identified with a non-Foster circuit around 77.06 GHz and a Foster circuit around 98.48 GHz and the corresponding identification will be based both on the negative inductance/capacitance, for the non-Foster part, and positive inductance/capacitance for the Foster part.
Starting from these considerations, the proposed identification procedure has been applied to b obtaining a fitting b fit that can be represented by the schematic circuit shown in Fig. 12(a), where the first branch is related to the capacitance between the lateral side of the rings of two adjacent unit cells (as previously discussed for the closed ring case), and the second and third branches to the non-Foster and Foster parts, respectively. The corresponding values obtained by the identification procedure are: C = 8.11 fF, L n F = −1.99 nH, C n F = −2.14 fF, L F = 0.52 nH, C F = 4.98 fF. The values of C, L F , C F are very similar to the case of closed ring, because the presence of the gap has a little effect on them, as previously discussed.
The frequency behavior of each branch of the identified b fit is shown in Fig. 12(b). The blue line in Fig. 12(b) corresponds to the pure capacitance, while the green and black lines refer to the non-Foster and Foster branches, respectively. As previously discussed, the non-Foster part shows a behavior with negative slope at any frequency with very low effect in the lower and upper parts of the band, while the Foster part, with positive slope everywhere, has an effect on the whole frequency band. Also in this case, the presence of "internal" non-Foster elements is strictly related to the "compression" effects of the proposed equivalent circuit. Hence, the non-Foster elements must be considered only as a part and cannot be separated from the other components of the overall equivalent circuit of the passive lossless device. For example, near the non-Foster resonance the strongly dispersive behavior of the transmission line length θ r shown in Fig. 10(b) together with the behavior of the non-Foster branch "compensate" each other producing a global behavior at the input-output ports of the overall equivalent circuit that still has a Foster behavior being related to the passive lossless SRR. It should be noted that the transmission line with electrical length θ r shown in Fig. 10(b) could be replaced by a more complex network containing as the Foster and non-Foster elements (or positive RLC and negative R), "hidden" by the strong frequency dispersive behavior near the non-Foster resonance. Such a replacement will produce a very complex global network, making it very challenging to "read" the overall equivalent circuit. On the other hand, recalling that the transmission lines with electrical length θ r 1 = θ r 2 = θ r change only the phase of the scattering parameters to obtain the correct relationships between them at the input-output ports, we prefer to maintain the actual representation with the transmission lines, even if they show a frequency dispersive behavior, to obtain a model that is simple and easy to handle.
It should be helpful to understand whether the representation with three shunt branches could be improved. To this intent, the electric field behavior in the neighborhood of the conductor gap is shown in Fig. 13(a)-(e) at five different frequency values, 10, 74.96, 77.06, 98.48, and 150 GHz, corresponding to the lowest (10) and highest (150) frequencies of the analyzed band, the total transmission frequency (74.96), the non-Foster (77.06), and Foster (98.48) total reflection frequencies. The plots are evaluated in the reference plane A-A' placed in between the ring conductors as shown in Fig. 1(b). The electric field behavior shows that near the total transmission and non-Foster resonance, Fig. 13(b) and (c), there is a change in polarity of the charges in the conductors of the gap, while, at the other frequencies, the polarity is always the same in the two gap conductors. Similar behavior occurs for the charges of the straight conductor of the left unit cell that have the opposite polarity with respect to that in the gap conductors of the right unit cell. At these frequencies, the null of charge at about the middle of the left straight conductor of the right unit cell in Fig. 13(b) and (c) induces to think that there should be a "voltage drop" across the gap as highlighted in the sketch (not in scale) of the electric field shown in Fig. 14(a), obtained by the CST simulations.
Due to this fact, the circuit representation of this behavior with a shunt non-Foster branch will be suitable to take into account this "voltage drop," and to overcome this problem, the shunt non-Foster branch is transformed into a series reactance as in Fig. 3 and previously discussed. Two equal transmission lines of electrical lengths θ bx evaluated by (16) embed a series impedance Z s = j X p F = − j B n F with resulting in the Foster series reactance with behavior shown in Fig. 15(b) and identified by an LC shunt with L F p = −C n F = 2.14 fH and C F p = −L n F = 1.99 nF. The complete identification of the shunt load B of the equivalent network in Fig. 12(a) becomes as shown in Fig. 15(a). It is interesting to analyze the behavior of θ bx shown in Fig. 15(b): the effects of the transmission lines vanish at the edges of the bandwidth, where they tend to π or 0 and also the series reactance X p F → 0, just as the non-Foster shunt load shown in Fig. 12(b) (green dotted line). In this representation, where the series branch has positive L F p , C F p , the non-Foster behavior of L n F , C n F has been converted in the anomalous dispersion of the transmission line θ bx , which could be replaced by a more complex network containing the Foster and non-Foster elements (or positive RLC and negative R), "hidden" by the anomalous frequency dispersive behavior.  Summing up, the presence of a non-Foster load in the equivalent network is due to the proposed identification process of B( f ) that is very easy to perform and produces both the Foster and non-Foster branches. While the Foster branches are directly identified with shunt capacitances or shunt LC series blocks with all the positive values, the non-Foster branches are identified by shunt LC series blocks with negative values for L and C. This last identification is due to a "compression" of the circuit representation in a shunt non-Foster load that is not able to correctly represent the "voltage drop" across the gap. Therefore, the transformation from the non-Foster shunt into the Foster series load permits to model more correctly the behavior of the SRR versus frequency with positive L and C for any branch.
As an example of the proposed identification process applied to SRR with gap, the variation in L F p and C F p with respect to L z or L x is shown in Fig. 16(a) where the comparison with the following approximations is also reported: C F p,app (L x ) ≈ 0.627 + 5.68 · 10 −3 L x + 12.8 · 10 −6 L 2 x nF (29) L F p,app (L z ) ≈ −0.276 + 5.24 · 10 −3 L z + 51.7 · 10 −6 L 2 z fH (30) An interesting remark can be done if we plot the product LC for shunt and series Foster cases versus L x or L z as in Fig. 16(b): the product LC is linear with L z or L x in the range 100 ÷ 250 and this can be verified by a series expansion around 175 of the LC product, obtained from (22), (24), and (25)-(26) or (28)- (31). For example This linear behavior of LC permits to predict the resonance frequency, f c = (1/2π √ LC), for any rectangular SRR by means of only two numerical simulations of two different SRRs with different length L z , for L x = 172 (or the opposite). In fact, by means of the identification process of the two resonant branches (L F , C F , L F p , C F p ) for only two values of L z , the parameters m and n of the straight line L F C F = m F + n F L z and L F p C F p = m F p + n F p L z can be easily found and the prediction of the resonant frequency for each branch can be done for any other values of L z . This can be a very fast and efficient approach if the SRR has to be developed at a certain resonance frequency.

C. Front Gap SRR
As a further study, a one-ring SRR with a gap aligned at front has been simulated, as shown in the inset of Fig. 17(a). The obtained scattering parameters seen from the air and in the dielectric are shown in Fig. 17(a) and only one resonance appears at f = 61.5 GHz. The identified equivalent circuit is no more symmetric with respect to the direction of propagation of the incident wave, and two different electrical lengths for the transmission lines are found as shown in Fig. 17(b). The identified normalized susceptance b is shown in Fig. 17(b) and only one resonance appears, due to the gap at front. In fact, the resonance of the closed ring completely disappeared in the presence of the front gap that drastically changes the current behavior around the ring, imposing a null of current exactly in the gap position. Moreover, the susceptance b is identified with a capacitor in the first branch and a Foster series LC resonator in the second branch as in Fig. 5(a), because the gap is in the front side of the conductor and it is correctly seen as a pure shunt load with positive L and C values. In fact, in this case there is no "voltage drop" along the direction of propagation of the incident wave as in the lateral gap case and, for the front gap case, the identified values are: C = 3.36 fF, L F = 0.97 nH, andC F = 6.88 fF.

D. Circular SRR
As a very brief example of the possibility to apply the proposed approach to other kind of geometries, we have analyzed the case of one circular ring without and with gap for the same excitation as the previous lateral gap SRR case. The unit cell contains only one circular ring with an inner radius of 128 µm, w = 2 µm, P x = 258 µm, and g = 2 µm (with gap case). The approach is similar and in this case we have obtained a frequency of resonance around 99 GHz for the closed ring (CRR). The equivalent circuit identified for B is constituted by two branches as in Fig. 5(a) with C = 7.78 fF, L F = 0.84 nH, andC F = 2.41 fF. The circular SRR case with lateral gap exhibits a resonance frequency around 99 GHz, similar to the CRR, and introduces also a new resonance around 62 GHz, corresponding to a non-Foster case. The corresponding identified circuit is similar as shown in Fig. 12(a) with C = 7.48 fF, L n F = −3.1 nH, C n F = −1.54 fF, L F = 0.78 nH, andC F = 2.59 fF.

E. Two-Ring SRR
The last studied SRR is the two-ring case shown in Fig. 1(a) with the same geometric parameters used for the one-ring case. The obtained equivalent circuit parameters corresponding to the scattering parameters of Fig. 18(a) are shown in Fig. 18(b). In this case, the presence of two rings moves the first two resonances to lower values (44.4 and 96.8 GHz) and this is due to the mutual capacitive coupling between the rings. A new third non-Foster resonance appears at about 125.8 GHz and this is due to the gap of the inner ring that acts as the one ring case as previously discussed. The b fit contains four shunt branches with one capacitance, two non-Foster, and one Foster LC series in the branches. Both the non-Foster resonances are related to the lateral gap and they are due to a "voltage drop" as in the one-ring case. Hence, it is possible to transform both the non-Foster branches in two series load yielding to the identified circuit for the total susceptance as shown in Fig. 18(d), with θ bx 1 , θ bx 2 as in Fig. 18(e). The identified values for the elements in the circuit of Fig. 18(d) are: C = 9.2 fF, L F p1 = 1.47 fH, C F p1 = 9.01 nF, L F = 0.443 nH, C F = 6.09 fF, L F p2 = 0.087 fH, and C F p2 = 18.06 nF. Also in this last case, it is interesting to compare the electric field maps and current distributions to better understand the relationships between the capacitances/inductances and electromagnetic fields.
The Foster resonance shows a behavior of the electric field similar to the one-ring cases as can be seen in Fig. 19(b). The two non-Foster resonances, Fig. 19(a) and (c), have a similar behavior of the electric field with two different kinds of capacitive coupling, the first between the straight conductors and the second between the arms of the inner conductors separated by the gap, as in the one-ring case. The difference between the two non-Foster fields is mainly in the polarity of the charge on the inner conductors.   A simplified sketch (not in scale) of the current density distributions on the conductors' surfaces of the external and internal rings and the magnetic field is shown in Fig. 20(a)-(c). Also in this case, the amplitude and direction of the current density and the related magnetic field are different for the Foster and non-Foster cases, drastically changing with frequency. Hence, the different behaviors of the electric and magnetic fields at the three frequencies justify the different values for the capacitances and inductances identified in the equivalent circuit of Fig. 18(d).

IV. CONCLUSION
A simple equivalent circuit for SRR has been proposed based on the improved approach of a shunt load embedded in two transmission lines. The reflection and transmission properties are contained in the shunt susceptance contained in the circuit. An identification process of a fitting susceptance based on Brune's synthesis has been applied, taking into account the presence of the Foster and non-Foster reactive loads. The non-Foster branches are at first identified with negative inductance and capacitance; then, they are related to a "voltage drop" across the ring gap and transformed into a series Foster loads, with positive inductance and capacitance, to make more effective the circuit representation of the actual electromagnetic fields' behavior in the SRR. The proposed approach is then applied to the two-ring case, and an effective equivalent circuit containing all the positive reactive elements is obtained, valid over a large bandwidth. The proposed method can be applied to the identification of equivalent circuits for any type of metamaterials with different geometries and electromagnetic nature as well as different frequency bands, from microwaves to sub-THz.