Analysis of nonlinear electromagnetic metamaterials

We analyze the properties of a nonlinear metamaterial formed by integrating nonlinear components or materials into the capacitive regions of metamaterial elements. A straightforward homogenization procedure leads to general expressions for the nonlinear susceptibilities of the composite metamaterial medium. The expressions are convenient, as they enable an inhomogeneous system of scattering elements to be described as a continuous medium using the standard notation of nonlinear optics. We illustrate the validity and accuracy of our theoretical framework by performing measurements on a fabricated metamaterial sample composed of an array of split ring resonators (SRRs) with packaged varactors embedded in the capacitive gaps, in a manner similar to that of Wang et al (2008 Opt. Express 16 16058). Because the SRRs exhibit a predominantly magnetic response to electromagnetic fields, the varactor-loaded SRR composite can be described as a magnetic material with nonlinear terms in its effective magnetic susceptibility. Treating the composite as a nonlinear effective medium, we can quantitatively assess the performance of the medium to enhance and facilitate nonlinear processes, including second harmonic generation, three- and four-wave mixing, self-focusing and other well-known nonlinear phenomena. We illustrate the accuracy of our approach by predicting the intensity-dependent resonance frequency shift in the effective permeability of the varactor-loaded SRR medium and comparing with experimental measurements.


Introduction
Metamaterials consist of arrays of magnetically or electrically polarizable elements. In many common configurations, the functional metamaterial elements are planar, metal patterns, on which currents are induced that flow in response to incident electromagnetic fields. Since these effective circuits are typically much smaller than the wavelengths over which the metamaterial is expected to operate, the polarizability of each element can be determined by applying a relatively straightforward quasi-static circuit model in which standard circuit parameters such as inductance, capacitance and resistance are introduced [1]- [5]. From the circuit point-ofview, it is then easily appreciated that the much larger range of properties routinely observed in metamaterials-including negative permittivity, permeability and refractive index-relates directly to the underlying electrical 'LC' resonances that the metamaterial circuits support. The strong, resonant response of the metamaterial circuit when coupled to applied fields translates to a strong polarizability in the language of material science. Applying homogenization methods to a collection of metamaterial elements, we thus obtain conceptually a continuous medium whose effective constitutive parameters can be controlled by the design of the metamaterial geometry. With the expanded palette of response available in metamaterials, novel electromagnetic structures and devices can be designed and demonstrated, negative index materials (NIMs) [6]- [8] and transformation optical media [9,10] being striking examples.
Although circuit methods are usually associated with radio-frequency and microwave applications, the circuit description for metamaterials provides an accurate model for higher frequency and optical metamaterials as well. At frequencies significantly lower than the plasma frequency of a conductor, it is possible to associate inductance with a conducting path and capacitance with gaps in the conducting path. However, the mechanism of inductance changes towards visible wavelengths, since the inductance becomes dominated by the inertial response of the charge carriers [11]. 3 A distinct property of metamaterials is that the capacitive regions strongly confine and enhance the local electric field. The capacitive region of a metamaterial element thus provides a natural entry point at which exogenous frequency-dispersive, active, tunable or nonlinear materials can be introduced, providing a mechanism for coupling the geometric LC resonance to other fundamental material properties [1]. The hybridization of metamaterials with semiconductors [12], ferromagnetic materials [13] and other externally tunable materials [14] has already expanded the realm of metamaterials from the passive, linear media initially considered.
Nonlinear materials represent a potentially useful class of materials to consider for integration into metamaterial composites. Many crystals and polymers exhibit substantial nonlinearities that are exploited across the spectrum for devices [15]. In conventional materials, the linear and nonlinear susceptibilities and their dispersions are intrinsically set, for example, by the fundamental anharmonic resonances associated with molecular systems. In metamaterial hybrid composites, however, those susceptibilities can potentially be altered independently with considerable control, retaining and even enhancing the nonlinear response of the embedded material while engineering the propagation characteristics of the otherwise linear metamaterial structure. This independent control provides an alternative route to the optimization of nonlinear materials, which form the basis of such optical devices as mixers, frequency doublers and modulators [16].
The incorporation of inherently nonlinear materials into the capacitive gaps of split ring resonators was considered theoretically by Zharov et al [17], who predicted an intensitydependent resonance frequency of the composite medium with associated bistability. Nonlinear phenomena such as harmonic generation [18,19], parametric down-conversion [20] or tunability [21,22] have been demonstrated experimentally either in planar metamaterial structures or in a single unit cell element. Many current metamaterial studies consider wave propagation in NIMs, assuming a homogeneous NIM layer with presumed values of linear and nonlinear response [23]- [32]. A diverse spectrum of nonlinear phenomena has been analyzed for NIMs, including the enhancement of interface second-harmonic generation near the zero-n gap of a negative-index Bragg grating [24,29], the enhancement of the parametric processes in binary metamaterials [30], the possibility of satisfying the phase matching condition for the counter-propagating waves for SH generation in a NIM [31], optical bistability in a nonlinear optical coupler with a negative index channel [32] or gap solitons in a nonlinear negative-index cavity [33].
For practical applications, the design and optimization of nonlinear metamaterial-based devices requires a quantitative approach, one that relates the particular metamaterial geometry incorporating nonlinear elements to the nonlinear properties of the composite, effective medium.
Steps towards this goal have been taken in [30,34], where the possibility of representing the nonlinear response of a metamaterial with diode insertions in terms of the effective second-order susceptibility has been discussed. Our goal here is to extend the general analytical framework for a metamaterial based on nonlinear elements integrated into the underlying circuit of each unit cell. Our approach is to expand the effective medium polarizability in a power series in terms of the applied field amplitude, thus relating the response of the metamaterial to the strength of the applied field, in the same manner as is done for conventional nonlinear crystals [15]. To emphasize the analogy between the anisotropic metamaterial hybrid media and natural crystals, we refer to the metamaterial structures described herein as metacrystals.
For a system close to resonance, which is the case for a resonant metamaterial, care should be taken regarding convergence of the power series at high powers [15]. Regardless 4 of this limitation, the series expansion representation is particularly useful since it allows the interpretation and quantitative prediction of an immense variety of nonlinear phenomena known to originate from certain orders of nonlinear response [15,35]. As an example, we provide an estimate of the applicable power levels using the varactor-loaded split ring resonator (VLSRR)based medium. Within the range of applicability the approach allows the development of the model of a homogeneous analogue nonlinear medium that aligns with the same concept in nonlinear optics.
In the following section, we introduce the basic circuit model description of a nonlinear metamaterial possessing a resonant magnetic response, determine the linear constitutive parameters via an analytic effective medium theory and derive expressions for the secondand third-order magnetic susceptibilities characterizing the nonlinear medium formed by the metacrystal. While there are many materials that could serve as the nonlinear component of a hybrid metamaterial composite, at low frequencies (below a few gigahertz), the low cost and wide availability of varactor diodes make them a convenient source for integration. The use of varactor diodes as a means of implementing tunability [21] and nonlinear response [36] in metamaterial hybrid media has been demonstrated experimentally. We relate our theory to a particular case of VLSRR-based medium in section 3. In section 4, we discuss the range of applicability and the accuracy of the solution. To illustrate the accuracy of our analytical expressions for the nonlinear susceptibilities, we present the results of experimental measurements on a VLSRR structure and compare the observed nonlinear response to that predicted in section 5. Both the theoretical and the experimental results are summarized in section 6.

General expressions for the nonlinear metacrystal
We consider here a nonlinear metacrystal formed from resonant circuit elements that couple strongly to the magnetic field. Although we specifically consider the split ring resonator (SRR) medium [1]- [3], the general circuit model described here should apply to other effective magnetic metacrystals such as those formed from cut wire pairs [37,38] or the fishnet structure [39,40]. The principal geometry of an SRR structure, which is essentially a current loop with a capacitive gap that breaks the otherwise continuous current path, is illustrated in figure 1(a), with its equivalent inductively driven RLC circuit model shown in figure 1(b). The induced electromotive force, ε(t), drives the circuit, producing a current I that satisfies In equation (1), L, R and V D = 1 C t 0 I (τ ) dτ are, respectively, the distributed inductance, distributed resistance and the induced voltage across the effective capacitor C of the circuit, with the particular interpretation of these parameters determined by the geometry and the frequency range considered. The linear circuit parameters can be estimated from analytical approximations; however, as will be seen later in this section, the parameters can also be obtained by curve fitting the numerically simulated effective permeability. The latter method, based on well-established numerical retrieval procedures, has proven extremely accurate for linear materials and avoids any reliance on analytical approximations [41]. Generally, any or all of the effective circuit parameters (R, L or C) can possess a nonlinear response to the amplitude of the induced electromotive force (and, hence, to the power of the applied field), becoming functions of the amplitude of charge oscillations in the circuit [42]. Introducing the time-dependent charge Q(t) = t 0 I (τ ) dτ , equation (1) can be written in terms of the normalized charge q = Q/C 0 , C 0 being the value of the capacitance in the linear regime (at low powers), as where ω 2 0 = 1/L 0 C 0 denotes the resonant frequency of the circuit in the linear regime, γ = R 0 /L 0 the linear damping constant and R 0 and L 0 the linear values of resistance and inductance. Note that the normalized charge q is in units of volts. Equation (2) is that of a driven oscillator, with the term F(q,q) representing a normalized nonlinear restoring force. We note that, whereas the above formal expression for γ accounts for only the resistive loss in the effective circuit, in practice, all loss sources, including the losses in the substrate, will contribute to the value of γ . This value can be estimated numerically, for example, by the decay time of the free oscillations of the initial charge at the effective capacitor [43].
Although nonlinear response could be introduced into any of the circuit parameters, the field enhancement within the capacitive gaps suggests that the nonlinearities of materials or components will be best leveraged by integrating them into the high-field regions of the circuit. We thus assume that only the capacitance C(q) has a nonlinear response, so that equation (2) reduces to [22] We expect that at low power levels the nonlinear response will be generally weaker than the linear response, so that we are justified in expanding the voltage in a Taylor series, retaining a limited number of terms in q. In this section, we write the terms of the expansion to third order, with higher (up to fifth order) terms considered in the appendix in connection with discussing the limits of convergence. The expansion can be written as and equation (3) then takes the standard form of a nonlinear oscillator: where α ≡ aω 2 0 and β ≡ bω 2 0 .

6
It can be shown [42] that equation (2) can equally be brought to the form of equation (5) if the effective inductance or both the inductance and the capacitance were to possess a nonlinear response, with the only difference occurring in the particular expressions for the nonlinear coefficients a and b. The form of the equation would also be similar in the case of a nonlinear resistive response, with the Taylor expansion in terms of the derivative of the normalized charge.
The nonlinear oscillator form of equation (5) can be solved by standard methods. Following Boyd [15], we assume that the response of the oscillator will be in proportion to the strength of the applied field, motivating a perturbative expansion of the form Here, λ is a dimensionless perturbation parameter (λ 1) indicating the field strength λε(t).
The tilde sign signifies variables that are time dependent. Substituting equation (6) into (5) and equating terms of the same order in λ, one can obtain the following system of coupled equations: Equations (7a)-(7c) describe, respectively, the first-(linear), second-and third-order responses of the normalized charge. Note that the first term in equation (6) is the same whether we assume a linear or nonlinear oscillator and that all subsequent higher-order terms depend on the amplitude of the lower terms. Because the SRR exhibits a predominantly magnetic response, we assume that the incident driving field corresponds to the magnetic component of the linear-polarized incident field or where the subscript y indicates the y-component of the magnetic field B and B y (ω n ) = B ry 2 e ikz with the factor of 1/2 arising from the reality of the field. Note that in the latter expression, B y is a complex amplitude dependent on the spatial coordinate, z, while B r y is a real (constant) number representing the amplitude of the real field at frequency ω n [15,35,44]. The axis of the circuit is taken to lie along the y-direction. We thus assume that the incident field can be expressed as a sum of components, each of which oscillates with an angular frequency ω n . According to Faraday's law, the driving electromotive force ε(t) is related to the amplitude of the incident magnetic flux according toε(t) = iA n=− ω n B y (ω n ) e −iω n t , where A represents the effective area enclosed by the circuit. In the above expression, the summation is taken over both positive and negative frequencies, where we have assumed B y (−ω n ) = B * y (ω n ) to arrive at a more compact notation. For compactness, henceforward we use the notation ω −n ≡ −ω n and omit the subscript y.
We look for a steady-state solution to equations (7) of the form where ω r ≡ ω n + ω m , ω s ≡ ω n + ω m + ω p and the summations are taken over both positive and negative frequencies, with n, m and p each taking values between ± . Note that, according to this notation, the summation over index r assumes the summation over all possible values of frequency ω r ≡ ω n + ω m ; for example, ω r includes the values of ±2ω 1 and 0 if = 1; similar considerations hold for the summation over s in equation (9c). Note also that, according to equation (8), the term with n = 0 is absent from the summation over . We point out that the term with the zero subscript, even if present in the summation, is not to be confused with the resonant frequency ω 0 of the circuit. We also do not explicitly include the dc term in the incident electromagnetic field (while it is not excluded that any of the frequencies ω n in the sum can take a 0 value, i.e. represent a dc field). The dc magnetic field would be expected to not produce a magnetic dipole response in the medium. The latter statement is confirmed later by the solution toequation(5),ifoneofthefrequenciesω n is taken to have a zero value. The dc term can, however, arise as a result of nonlinear frequency mixing, e.g. the fields with frequencies ω n and ω −n can formally produce a response at ω n + ω −n = 0. Using equation (9a) in (7a), one can write the following expression for each ω n : from which we obtain where the resonant denominator is defined as D(ω n ) ≡ ω 2 0 − ω 2 n − iγ ω n . Using equations (9a) and (11), equation (7b) can be written as where the indices m and n vary independently between ± . The right-hand side of equation (12) contains many combinatorial frequencies because of the summation over m and n. Assuming the form of the solution given by equation (9b), equation (12) can be transformed into a set of independent equations for each frequency ω r = ω n + ω m with the second-order response q (2) (ω r ) at each frequency ω r satisfying the equation In equation (13), the parentheses under the sum indicate that the sum ω n + ω m remains fixed, while the indices n and m vary. From equation (13), the resulting expression for the secondorder response at each ω r is The complete steady-state solution to equation (7b), describing the second-order response, is given by equation (9b), summing equation (14) over all possible values of ω r .
Following the same steps as above and using equations (11) and (14) in combination with equations (9) in equation (7c), we obtain for the third-order response at a frequency The expression for q (3) (ω s ) then follows directly from equation (15). However, the factor D(ω n + ω m ) in the first term on the right-hand side of equation (15) results from the secondorder response q (2) (ω r ) at the frequency ω r = ω n + ω m with n and m taking any values between ± that produce in combination the frequency ω r . To account for the various contributions, this factor must be permuted over all possible combinations of indices n and m. We arrive then at the following expression for the third-order response at the frequency ω s accounting for the second-and the third-order nonlinearity: The complete steady-state solution to equation (7c) is given by equation (9c) with q (3) (ω s ) determined from equation (16). Equations (11), (14) and (16), in combination with equations (9), provide general expressions for the first-, second-and third-order responses, respectively, accounting for the second-and the third-order nonlinearities. These expressions can be used to derive the first-, second-and third-order magnetic susceptibilities that will characterize the metacrystal.
For a dilute medium, the magnetization,M, is approximately related to the magnetic momentm according toM(t) = Nm(t), where N is the volume density of moments. The magnetic dipole moment of the effective circuit encompassing the effective area A is given approximately dt . Following the perturbative expansion for q(t) given by equations (6) and (9), the magnetization can be written as where, as before, ω r ≡ ω n + ω m , ω s ≡ ω n + ω m + ω p and each of indices n, m and p takes values between ± .
We assume that the magnetization can be expressed as a power series in terms of the strength of the applied field, according tõ where and where we have adopted the conventional notation for the arguments of the nonlinear susceptibility in which the first frequency term is the sum of the subsequent frequency arguments [15]. The subscript y in this notation reflects the Cartesian coordinates of each participating field component and of the resulting magnetization, which are all polarized along the y-axis in the present case. Equating the terms with equal powers of the exponents in equations (17) and (18) and employing equations (11), (14) and (16) for q (1) (ω n ), q (2) (ω r ) and q (3) (ω s ), we obtain, respectively, the linear, second-and third-order magnetic susceptibilities as follows: where each frequency can take both positive and negative values with the indices n, m and p each varying between ± ; µ 0 is the permeability of vacuum. From an inspection of equation (20), we see that the linear susceptibility of the SRR can be expressed in terms of its geometrical and electrical parameters as where F ≡ N ω 2 0 µ 0 A 2 C 0 , in qualitative agreement with more detailed analytical studies [3]. We do not seek more accurate analytical expressions for the effective circuit parameters, since we can easily find the linear properties through numerical retrievals, performing a fitting to determine the unknown coefficients F and γ (entering through the factor D(ω n )) in equation (23). The higher-order terms do not add to this parameter set, so in practice it suffices to have the initial expression in equations (20) or (23). We note also that the presented model considers a single resonance in the SRR, which is a known limitation of the nonlinear oscillator presentation, i.e. the existence of higher-order resonances is not accounted for. Higher-order resonances can strongly couple to the higher harmonic field components; such interactions would serve as the topic of a useful future study.

Examples for particular combinatorial frequencies
As an example of the use of equations (20)-(22), we write in a more explicit form the susceptibilities for some particular combinatorial frequencies for second-and third-order processes. For example, for the process of second-harmonic generation, we have ω m = ω n and, from equation (21), Similarly, for difference frequency generation, the frequency ω m is negative and, using explicitly a negative sign in equation (21), we obtain where n and m can take any values between ± . The susceptibilities χ (2) yyy (−2ω n ) and χ (2) yyy (−(ω n − ω m )) are straightforwardly seen to be the conjugates of equations (24a) and (24b). Note that, according to equation (21), the nonlinear susceptibility that is responsible for the effect of optical rectification, χ (2) (0; ω n , −ω n ), vanishes, despite the fact that the response q (2) (0) is nonzero according to equation (14). This is an expected result since, for the type of metacrystals we consider, the magnetic susceptibility depends on the ac response of the medium, while the response q (2) (0) does not vary in time. The zero-frequency second-order response, in combination with the linear response or a higher-order response, does contribute, however, to the higher-order susceptibilities. This contribution is addressed in the appendix, where we provide a more intuitive derivation of the higher-order response for some combinatorial frequencies.
As can be seen from equation (22), many more combinatorial frequencies arise as the result of third-order nonlinear processes. We provide, as an example, the explicit expressions for the third-order susceptibilities χ (3) yyyy (3ω n ) responsible for third-harmonic generation, and χ (3) yyyy (ω n ), relating the nonlinear response at the fundamental frequency ω n (equivalently described as a power-dependent refractive index).
For the susceptibility at the third harmonic, one has ω n = ω m = ω p in equation (22), which gives The response at the fundamental frequency ω n can result both from the interaction of a field at ω n with itself, leading to the effect of self-phase modulation, and from the interaction of two distinct fields, leading to a cross-phase modulation process [15,35]. For the first case, we set ω m = ω n and ω p = −ω n in equation (22), resulting in while for a similar response leading to cross-phase modulation, using ω p = −ω m in equation (22), we obtain Note that equation (25c) approaches (25b) as ω m approaches ω n , as expected for the susceptibilities χ (3) yyyy (ω n ; ω n , ω n , −ω n ) and χ (3) yyyy (ω n ; ω n , ω m , −ω m ) [15]. Using equations (25b) and (25c) in (19c) and performing the summation to account for the possible number of permutations leading to each response, we can also see that, as expected, the nonlinear magnetization produced by two distinct fields is twice as large as the one produced by a single field: The effective permeability characterizing the response at other combinatorial frequencies can be obtained from equations (21) and (22) in a similar manner. We emphasize again that no new geometrical or linear circuit parameters are introduced into the higher-order terms; aside from factors relating to the nonlinear component (a and b in this case), all unknown coefficients can be found by analyzing the linear response.
According to equations (20), (21) and (22), we can also alternatively express χ (2) (ω r ) and χ (3) (ω s ) as products of linear magnetic susceptibilities at the contributing frequencies, allowing for a more explicit relation between the linear and nonlinear susceptibilities: and In particular, for the susceptibilities at the second harmonic and at the fundamental frequency, 12 and χ (3) yyyy (ω n ; ω n , ω n , −ω n ) =

Intensity-dependent permeability
Besides the effects of self-and cross-phase modulation, the third-order susceptibility at the fundamental frequency also leads to an intensity-dependent refractive index [15] (in our case, due to the intensity-dependent permeability) and, as discussed later, to a shift in the resonance frequency. The expression for the intensity-dependent permeability can be obtained by considering the total nonlinear magnetization at the fundamental frequency expanded up to the third order. Assuming a single incident field at ω n we have according to equations (18) and (19) M tot so that the intensity-dependent permeability is given by We can expect from equation (32) a shift in the spectral position of the resonance with increased intensity compared to the linear case. Numerical examples of this shift will be presented in the following section, in which we apply the present theory to a particular case of a VLSRR-based medium.

Application to the varactor-loaded split ring resonator (VLSRR) medium
A convenient means of achieving a nonlinear medium is to connect a varactor diode across the gap of an SRR. Depending on the number and orientation of the diodes integrated into the medium, different nonlinear properties will result. We investigate two specific configurations, one in which a single diode is inserted into the gap of an SRR and the other in which two diodes are inserted into two gaps in the SRR, as shown in figures 2(a) and (b). In the latter case, the varactors are oriented in a back-to-back configuration to prevent dissipative currents from flowing at large power levels, as occurs in the single-gap VLSRR [22].

Single-gap VLSRR
The nonlinear properties of individual VLSRRs have been analyzed theoretically and experimentally by Wang et al [22], who demonstrated power-dependent resonance frequency shifting and bistable behavior. The nonlinearity associated with the diode elements occurs through the voltage-dependent capacitance of the varactors. For the particular varactors studied in [22] and used here (Skyworks SMV1231), the junction capacitance has the form [45] C where C 0 is the zero bias capacitance, K is the gradient coefficient and V P is the intrinsic potential. For small field amplitudes (and hence voltages), the time-dependent charge that accumulates as a function of the voltage can be found as Inverting equation (34) The diode thus provides a nonlinear charge/voltage relationship expressed as V D (q) that can be Taylor expanded for small arguments. Assuming a small argument expansion in equation (35), the expressions for the nonlinear coefficients that appear in equation (4) are the following: For the particular case of the Skyworks SMV1231 varactor, according to the specifications, K = 0.8, V P = 1.5 V [45], and hence a = −0.2667 s −1 V −1 and b = 0.0356 s −2 V −2 . If only one gap is present in the VLSRR, then the expressions derived in section 2 for the linear and higher-order susceptibilities can be used in their existing form, with the coefficients given by equation (36). The orientation of the unit cell relative to the incident field is such that the vector of magnetic inductance is normal to the plane of the ring and the electrical field is parallel to that plane, as shown in figure 2(c).

Double-gap VLSRR
The single-gap VLSRR is not optimal, since the current flowing in the SRR increases substantially with increasing power, leading to a bias voltage appearing across the diode large enough to provide considerable forward bias, increasing the resistive loss. To avoid the excitation of large currents, the diodes can be integrated into the SRR in an opposing configuration, such that one of the diodes is always reverse biased. The orientation of the backto-back diodes is illustrated in figure 2(b). Equation (1) in this case becomes where V D1 and V D2 are the voltages at each of the effective capacitors. Due to the opposite configuration of the varactors, the same incident field will induce a forward bias on one of the varactors and the reverse bias on the other. Thus, the back-to-back configuration leads to the opposite sign of the induced voltage at the varactors and also to the opposite sign of the accumulated charge at the effective capacitors representing the varactors. Assuming that both varactors are otherwise identical leads to the condition V D2 (q) = −V D1 (−q) in equation (37). Accounting for this condition, equation (3) can be written in terms of q as follows: As a result, the sum V D1 + V D2 becomes an odd function of the normalized charge q, as illustrated in figure 3(b) and the even-order terms in the Taylor expansion of V D1 + V D2 in equation (37) are canceled. Consequently, there is no contribution to χ (3) yyyy from the secondorder susceptibility. Equation (5) in this case is where ω d = √ 2ω 0 is the resonant frequency of the circuit in the double-gap configuration and ω 0 is the resonant frequency from the single-gap case. Following the same procedure as that used before, we arrive at the following expression for the third-order susceptibility: where D d (ω) ≡ ω 2 d − ω 2 − iγ ω. In particular, equation (25b) for the third-order susceptibility at the fundamental frequency becomes χ (3) yyyy (ω n ; ω n , ω n , −ω n ) = −2b   Figure 4 shows examples of the second-and the third-order nonlinear magnetic susceptibilities and the effective permeability obtained according to equations (24a) for the single-gap case and according to equations (25a), (25b) and (29) for both single-and double-gap configurations, with the nonlinear coefficients a and b given by equation (36). In order to estimate the values of γ and F for the construction of figure 4, we numerically simulate the linear propagation of a plane wave through a single layer of the metacrystal with the following parameters: the unit cell size 10 mm, the diameter of the ring and the width of the metal strip 8 and 0.5 mm, respectively, and the gap size 1 mm, as shown in figure 2. The substrate of thickness 0.2 mm is assumed to be made of FR4 circuit board material with = 4.4(1 + 0.02i). The simulations are performed using the finite integration-based frequency domain solver, contained in the commercial package CST Microwave Studio. From the computed transmission and reflection coefficients of a single layer of metacrystal, we employ a standard retrieval procedure [41] to find the frequencydependent constitutive parameters. We calculate F and γ , as well as the resonant frequency, by fitting the resonance curve of the effective permeability retrieved via this procedure.

Examples
From the numerical simulations applied to the unit cell in figure 2, we obtain the values of γ = 0.1202 × 10 9 s −1 , F = 0.1426 and ω 0 = 5.1459 × 10 9 rad s −1 for the singlegap metacrystal. Note that, assuming that the area of the ring is given by the parameters in figure 2 (A ≈ 5 × 10 −6 m 2 ), the estimated value of F leads to a slightly different zero bias capacitance value of C 0 = 1.6 pF than in the varactor specifications (2.4 pF [45]). This discrepancy is not unreasonable since the actual value of C 0 is modified by the packaging capacitance, the parasitic capacitance of the varactor and the capacitance of the SRR gap itself.
In general, additional considerations would be necessary to account for the impact of spatial dispersion in the theoretical approach. Since, according to equations (27) and (28), the nonlinear susceptibility can be expressed as a product of linear susceptibilities at the participating frequencies, the first-order correction to the nonlinear susceptibility in terms of spatial dispersion effects could be made by accounting for spatial dispersion using the analytic correction factor in each of the linear susceptibilities present in the expression [46]. Additional corrections might also be necessary to the forms of equations of (27) and (28). Fitting the retrieved effective parameters obtained from the simulations thus partially allows the effects of spatial dispersion to be included if such are present [46], by accounting for the spatial dispersion in the linear response. In the present example, however, the size of the unit cell was small enough compared to the wavelength such that spatial dispersion effects were unpronounced. We verified that accounting for the spatial dispersion factor in the linear effective parameters retrieval made almost no difference compared to using the Lorenz model in fitting the linear effective parameters retrieved from the S-parameters.
Note also that, although we do not take into account any interactions between the metamaterial elements in the above theory (i.e. no Lorentz factor), the simulations intrinsically take into account the entire periodic system (at least in the plane perpendicular to propagation). Thus, the numerical solution combined with the fitting procedure accounts for the effective medium properties of the metacrystal, which mainly leads to a renormalization of the various parameters of the resonant form (i.e. oscillator strength, damping factor, resonance frequency, etc). Such a renormalization is certainly not a rigorous theoretical treatment but rather a phenomenological procedure justified by the agreement between the final theoretical results and the experimental data, as discussed in the next section. In the future, a more detailed theoretical study accounting for the interactions between the unit cells would be beneficial.
Due to the resonant nature of the nonlinear magnetization of the metacrystal, the effective susceptibilities are complex-valued functions exhibiting a strong dispersion in the frequency range near the fundamental resonance, as shown in figure 4. The presented frequency range corresponds to the frequency of the excitation field, i.e. the resonant nature of the susceptibilities is due to the system resonance at the fundamental frequency. Hence, for example, there is no linear absorption for the field at the harmonic frequencies. The imaginary part of the susceptibility of the second or the third harmonic does not contribute to the time-averaged absorbed power or to the stored energy of a component of a polarization at the harmonic frequency (2ω or 3ω) in a field at the fundamental frequency, and just relates the phase of the harmonic polarization to the phase of light at the fundamental frequency [44]. The imaginary part of the third-order susceptibility χ (3) (ω) gives a nonlinear addition to the linear absorption at the fundamental frequency. Note that it differs both in sign and in phase from the analogous nonlinear susceptibility due to electrical polarization [15]-an expected result due to the delayed nature of magnetic response. The detailed comparison of these two cases would be an interesting analysis to pursue but is beyond the scope of the present discussion.
As seen from figure 4, the values of χ (2) and χ (3) near the resonance are many orders of magnitude larger than in a standard nonlinear medium, which is the consequence of both the proximity to the resonance and the strong nonlinearity the varactor diodes possess. Noting that the metacrystals are inherently resonant materials and hence the resonance region is important to consider, it is necessary to estimate the range of the field strength for which the truncation of the series expansion of the normalized voltage q or the magnetization M at the third-order of the field strength is valid. We provide such an estimate in the following section. The maximum applicable field amplitude ensuring the convergence of the series expansion in each of the singlegap or the double-gap configurations was used in constructing figure 4.
As seen from figure 4(e) and the inset in figure 4(f) and as discussed earlier in connection with equation (32), the position of the effective permeability resonance shifts with increased power. According to equation (32), the theory predicts the opposite direction of the resonant frequency shift for the two types of media considered, noting the opposite signs of χ (3) (5) and (39) and with the experimental results, as discussed later in sections 4 and 5.
Another observation from figure 4 is that the third-order susceptibility for third-harmonic generation is an order of magnitude smaller than the χ (3) (ω), as seen in figures 4(b) and (c). This indicates that we can expect a much weaker effect of third-harmonic generation than the resonant frequency shift. A weaker response at the third harmonic clearly follows from equations (25a) and (25b), since all the factors D(ω n ) are resonant in the denominator in the susceptibility χ (3) (ω), while this is not the case for the third-harmonic susceptibility χ (3) (3ω), in agreement with the discussion in [30] for the second-order response. The latter observation is general for the response at the fundamental frequency compared to any other combinatorial frequency for a medium possessing a single magnetic resonance: while the nonlinear response of a metacrystal is resonant at any frequency combination, a generally larger response can always be expected at the fundamental resonance.

Limits of validity
The diode system we consider experimentally exhibits strong nonlinearity, and we are working close to the resonance of the effective medium, so that it is a legitimate concern as to whether the series expansions used in equations (4)-(6) and in (17) and (18) are convergent or appropriate. Specifically, since we seek to describe a medium by its second-and third-order nonlinear susceptibilities, it is important to determine whether the higher-order terms can be neglected. To address this concern, we perform a rough estimate by finding the value of the field strength that ensures that the values of the fourth and fifth terms in the series expansion are small relative to the second-and third-order terms at the corresponding frequencies (the necessary truncation criteria for the third-order approximation to be valid). The estimated value of the field strength can then be verified by comparing the solution obtained by the perturbation approach with the exact numerical solution of equation (5). We proceed in this way, deriving the expressions for the fourth and fifth orders of the magnetic susceptibility of the SRR medium in a similar way to that above. The derivation and expressions for the fourth-order response at 2ω and the fifth-order response at ω are provided in the appendix.
Comparing the contributions to the amplitude of the magnetization response at ω coming from the terms 3χ (3) (ω)H 3 and 10χ (5) (ω)H 5 (the factors before these terms coming from the number of distinct permutations over n, m, p leading to the response at ω), we find that the maximum applicable field satisfying the necessary truncation criteria for the series truncated at the third order of the field strength is about 27 mA m −1 . The same field strength ensures the necessary truncation condition for the response at 2ω truncated at the second order of the field strength. From an examination of the higher-order terms in the series expansion for the doublegap VLSRR metacrystal, the maximum applicable H for up to the third-order expansion in field magnitude is estimated to be about 100 mA m −1 . The larger validity range is expected in this case since the higher-order effective permeability terms decrease more rapidly in the absence of contributions from the even orders.
To investigate further the accuracy of the solution, we solve equation (3) numerically using the exact expression for the voltage V D given by equation (35) and using the expansion of V D up to the third order in the power of q, according to and compare both results with the perturbation solution for the normalized voltage q. In equation (42), the driving voltage E dr is given by E dr = AB r ω = 2µ 0 H Aω. We solved differential equation (42) using the ODE15 subroutine in Matlab (Mathworks), scanning the frequency ω of the applied field over the resonant region. For each ω, we determine the amplitude of the response at the fundamental frequency ω and at each of the harmonics to obtain the resonant curve. We then repeat the procedure for various power levels. We also solve an equation similar to equation (42) for the double-gap case, assuming that V D (q) is the sum of voltages at the two capacitors, using a = 0 in its Taylor expansion, and replacing ω 0 with ω d in the linear and the cubic terms, as discussed above. Figure 5 shows a comparison of the numerical and analytical solutions for the fundamental frequency for several values of the applied field amplitude. We present the results for the field strength values corresponding to the range of power levels used in the experiment as discussed later in section 5. For the single-gap case ( figure 5(a)), the agreement is excellent for low magnetic field amplitudes of up to about 27 mA m −1 , in agreement with the value obtained earlier by analyzing the convergence of the series expansion. We also see good agreement for amplitudes up to about 22 mA m −1 for the double-gap VLSRR metacrystal. For larger field strength, the perturbative solution starts to deviate considerably from the numerical one near the resonance region, although it still exhibits very good agreement off-resonance. Accounting for the fifth-order term in the perturbative solution, as indicated by the green dash-dot lines in figure 5, improves slightly the resulting amplitude and the position of the resonance, while leading to a stronger discrepancy in the resonance shape. The numerical solutions obtained using the exact expression for V D and using its expansion in the power series just up to the third order, shown in gray dashed lines, agree very well with each other for all power values (the curves overlap in figure 5(b)), indicating that the inaccuracy in analytical solution at high powers is due to the approximations involved with using the perturbation approach rather than with the Taylor expansion of the V D itself.
As discussed before and as seen from both figure 5(a) and (b), the third-order response at the fundamental frequency produces a shift in the resonant frequency of the metacrystal. The resonance shifts in opposite directions for the single-gap and double-gap media, in agreement with the results predicted by the expression for intensity-dependent permeability (equation (32)) and discussed in the previous section. The opposite direction of the resonance shift obtained from the numerical solution for the two types of media can also be directly predicted from equations (5) and (39), noting the opposite signs of the nonlinear coefficients in these equations (since a < 0, b > 0 and a b, the resulting nonlinear term is negative in the single-gap case but is positive, due to a = 0, in the double-gap case).
The resonance in the solution for the normalized charge q and, consequently, in the effective nonlinear permeability characterizing the medium manifests itself as the minimum in the transmission spectrum of a plane wave propagating through the metacrystal. Therefore, the accuracy of the obtained solutions can be verified experimentally by observing the shift in transmission minimum as the intensity of the plane wave propagating through the varactorloaded metacrystal is increased and comparing the amount of resonance shift with that predicted by the theory and the numerical solution. We provide such a comparison with an experiment in the following section.

Experimental comparisons
To provide an experimental verification of the expressions derived above, a VLSRR medium was constructed, one unit cell in thickness (along the propagation direction) and 3 × 15 elements along the lateral directions. VLSRR metacrystals were made with the same unit cell parameters as those used in the above analysis, for both the single-and double-gap VLSRR cases. The metacrystal samples of the two kinds are shown in figure 6(a). A transmission line supporting TEM wave propagation below 2 GHz was used to measure the transmittance through the metacrystal samples, as shown in figure 6(b). An Agilent vector network analyzer (PNA N5230A) was used to launch microwaves into the transmission line and detect the transmitted fields. The frequency-dependent transmission properties associated with both the waveguide and connecting cables were removed using a standard calibration method. The excitation power P wg from the network analyzer ranged from −10 dBm to +15 dBm. The total loss from the transmission line structure, including the connection cables and adapters, was measured to be between 4 and 6 dB, so that the actual power exciting the sample ranged from −14 to 11 dBm, which spans the field range considered valid in the analytical calculations above. The  incident power can be related to the applied magnetic field as P = 2S 0 Z 0 |H | 2 [47], where Z 0 is the wave impedance (377 ) and S 0 is the effective mode area. In the calculations presented here, we assume that the mode area is approximately equal to the cross-sectional area of the waveguide, S 0 ≈ 4.5 × 10 −3 m 2 . We note, however, that the actual power inside the waveguide and the corresponding magnetic field are accurate up to about 2 dBm due to the imperfect loss measurement and the approximate mode area used in the calculation. The resulting transmission spectrum for the above range of powers is depicted in figure 7. In agreement with the theoretical predictions, we see that, when the power is increased, the transmission minimum shifts toward lower frequencies for the single-gap metacrystal and toward higher frequencies for the double-gap medium. In the case of a single-gap metacrystal, we also see a decrease in the amplitude of the response with increased power. This is a consequence of the dissipative current running through the varactor [22], which is not accounted for in the present analytical model.
In figure 8, we present a comparison of the experimentally obtained resonant frequency dependence of the single-and double-gap VLSRR metacrystals on the incident power with the exact numerical solution of equation (42) and with the perturbative solution expanded up to the third and fifth orders in field strength. The agreement between the exact numerical solution and the experimentally obtained resonant frequency is very good for the single-gap medium, validating the use of the effective circuit model leading to the nonlinear oscillator equation for the analytical description of the metacrystals. The slightly higher rate of decrease of the resonant frequency with power seen in the experimental data in figure 8(a) compared to the numerical results can be attributed to the contribution from the nonlinear resistance to the resonant frequency shift as mentioned above. For the double-gap metacrystal, a slightly higher rate of shift of the resonant frequency is observed in experiment than predicted by the numerical solution.
The agreement with the solution obtained by the perturbative approach holds for the range of powers discussed in connection with figures 4 and 5, and accounting for the fifth-order term in the series expansion leads to a small increase in the range of validity of the approach, as seen from the green dashed-dotted lines in figure 8. According to figure 8(b), the agreement between the experimental and the analytical resonant shift for the double-gap metacrystal holds up to about 2 dBm, corresponding to a field value of 21 mA m −1 , while the analytically predicted value is 100 mA m −1 . This difference can be attributed to the fact that the actual double-gap unit cells are not ideally symmetric, leading to a non-ideal cancellation of the even-order terms in the voltage expansion. The contribution from these terms becomes pronounced once the field strength increases, producing a stronger resonance shift.
One can verify that a similar discrepancy would hold near the resonance when comparing the perturbative and exact solutions of a Drude-Lorentz oscillator at high power levels, using the expressions for the electrical nonlinear susceptibility in the case of natural materials at optical frequencies. For example, assuming the perturbative solution for a non-centrosymmetric material with the regular for optical case nonlinear parameters in the Lorentz model [16], the noticeable discrepancy near the resonance would start at intensity levels of about 1 TW m −2 . Normally, the Lorentz atom model is not the approach of choice for describing the resonant effects in natural materials at optical frequencies, and usually other techniques such as the twolevel atom model [15] are used for the analysis of the resonant interaction of an optical field with a nonlinear medium. However, in the case of the VLSRR metacrystals considered here, as seen from the comparison of the experimental transmission resonance shift and the exact numerical 22 solution shown in figures 8, the nonlinear oscillator equation itself (equations (5) and (39)) provides a good model describing the interaction of the electromagnetic wave with a metacrystal even at large power levels, both off-and near the resonance. Noting that the metacrystals considered here are inherently resonant materials, a modification to the perturbation approach allowing better agreement with the numerical solution at high powers could be justified, but is beyond the scope of the present paper. For low power levels as discussed in this and previous sections, the perturbative approach solution works well for characterizing nonlinear metacrystals in terms of second-and third-order susceptibilities.

Conclusions
We have employed a perturbative solution to the nonlinear oscillator model of the effective RLC circuit of the unit cell to characterize the nonlinear properties of a metacrystal formed by unit cells that couple resonantly to the magnetic field. The nonlinear response of the metacrystals is characterized in terms of the series expansion of the nonlinear magnetization. We have provided general expressions for the effective magnetic nonlinear susceptibilities up to fifth order and discussed the valid power ranges for which the series expansion can be truncated at the third order of the field strength for the case of a varactor-loaded nonlinear metacrystal.
The expressions are convenient for the prediction, analysis and possible enhancement of the nonlinear response of a metacrystal. In particular, while, according to the theory, for a medium exhibiting a single magnetic resonance the nonlinear response is resonant at all combinatorial frequencies, a generally stronger (an order of magnitude in the case of a VLSRR-based medium) nonlinear response is expected at the fundamental than at any other frequency, in agreement with the discussion in [30] for the case of the second-order susceptibility. The absence of optical rectification in the medium with a magnetic nonlinear response also follows directly from the presented theory.
We analyzed in greater detail the third-and the fifth-order responses at the fundamental resonant frequency, which leads to the effects of self-phase modulation, intensity-dependent effective permeability and resonant frequency shift. We compared the results for the response obtained by perturbation solution with that found by numerical solution of the exact nonlinear oscillator equation, as well as with experimental transmission measurements of the shift of the resonant frequency of the metacrystal with increased incident power. The comparisons showed good agreement between the experimental frequency shift and the one obtained by numerical solution of the exact nonlinear oscillator equation for all power ranges considered (up to 10 dBm). The comparison also enabled an estimation of the range of validity of the series expansion representation of the nonlinear magnetization for the type of medium considered, and demonstrated that, near the resonance, the varactor-loaded nonlinear metacrystals could be described by the effective second-and third-order nonlinear susceptibilities for values of the excitation magnetic field amplitude up to 27 and 21 mA m −1 for the single-and doublegap metacrystals, respectively. The latter values correspond to incident powers of about 4 and 2 dBm. There is also very good off-resonance agreement between the perturbation approach and the numerical solutions for all the power ranges considered.
While the perturbative solution deviates considerably near resonance both from the numerical solution and from the experimental measurements at higher powers, the good agreement between the numerical solution of the same nonlinear oscillator equation and the experimental transmission results indicates that, in general, the nonlinear oscillator model provides a very close description of the interaction of the electromagnetic wave with a metacrystal at microwave frequencies even at high powers, both off-and near the resonance. In the range of validity of the obtained solution, there is very good agreement between the analytical, numerical and experimental results for the metacrystal response, indicating that the analytical method can be potentially used to investigate other nonlinear phenomena in metacrystals.