Nonlinear oscillator metamaterial model: numerical and experimental verification

E. Poutrina; D. Huang; Y. Urzhumov; D. R. Smith

doi:10.1364/OE.19.008312

1. Introduction

Metamaterials are frequently composed of metal inclusions patterned on a dielectric substrate or embedded in a dielectric host medium. Being of subwavelength dimensions, the metamaterial inclusions can be thought of as quasi-static circuits that can be designed to exhibit a resonant response to applied electromagnetic fields over specific wavelength bands. The inhomogeneous local structure of metamaterial inclusions can in turn lead to a highly nonuniform field distribution within the inclusions; strong field enhancements, for example, can occur within the capacitive regions of an inclusion. A nonlinear material or element placed in the vicinity of such an intense field region would be expected to exhibit an enhanced nonlinear response [1], leading to the enhancement of the effective nonlinear response of the composite medium. Indeed, theoretical studies [2] and initial experiments have demonstrated a variety of nonlinear phenomena in such as tunability with bistable operation [3, 4], second-harmonic (SH) generation [5, 6], or parametric down-conversion [7].

In recent work by Poutrina et al [8], an analytical approach was developed allowing the description of a nonlinear metamaterial in terms of effective nonlinear susceptibilities, in analogy with the standard representation common in nonlinear optics [9]. (To accentuate this analogy, we refer to the artificial nonlinear structures discussed here as metacrystals.) The analytical model employs a perturbative solution to the nonlinear oscillator model that describes a meta-material inclusion in terms of an effective RLC circuit. As an application of the approach, analytic forms for the second- through fifth-order susceptibilities were derived for a varactor-loaded split ring resonator (VLSRR) medium. In initial experimental studies, the second harmonic generation (SHG) from a fabricated VLSRR medium was measured and, applying a retrieval procedure based on the nonlinear transfer matrix method (TMM) [11], the second order nonlinear susceptibility of the VLSRR structure was determined [10]. In a second experiment [12], the intensity-dependent shift of the resonance frequency of the VLSRR medium–an effect related to the third-order susceptibility term–was investigated at low power levels.

Though these initial experiments have been successful in illustrating the validity of a homogenized nonlinear description of a metamaterial, the validation approaches rely on approximations that might prove unreliable in certain situations. The nonlinear TMM, for example, makes use of the non-depleted pump approximation, limiting potentially its applicability for certain power levels; in addition, for the SHG retrieval process, only a second-order susceptibility was assumed in the model, whereas higher order even terms could also have contributed to observed second harmonic signal. Moreover, the previous testing considered only a single value of the applied power. Noting a noticeable amount of the parameters involved (e.g., accounting for losses before and after the media for the fundamental field, estimating losses for the second harmonic field, and accounting for the fraction of power contained between the plates of the transmission line where the plane wave approximation can be used), a question remained whether the good agreement obtained in the previous study was not an occasional result valid at the particular power level. To address such concerns, here we perform a sequence of experimental and numerical investigations of a VLSRR medium, determining the second order magnetization of the medium for several incident power levels, with the applied power being the only varying parameter. The approaches developed here, while providing further evidence of the applicability of the analytical model, also allow a path to characterization of nonlinear metacrystals beyond the approximations inherent in the prior methods.

2. Analytical approach and numerical simulations

Figures 1(a) and 1(b) show the orientation of the VLSRR-based unit cell with respect to the incident field and its equivalent representation as an effective RLC-circuit. The response of the SRR can be expressed in terms of the charge across the capacitive gap, which is described by the following nonlinear oscillator equation [8]

\ddot{\tilde{q}} + γ \dot{\tilde{q}} + ω_{0}^{2} V_{D} (\tilde{q}) = - ω_{0}^{2} A μ_{0} {\dot{\tilde{H}}}_{y}

where q̃(t) is the normalized charge, V_D(q̃) is the voltage across the effective capacitance, ω ₀ is the (linear) resonant frequency of the unit cell, gamma is the damping factor, A is the area encompassed by the circuit, μ ₀ is the permeability of vacuum, and H̃_y(t) is the incident magnetic field assumed to be composed of Λ plane waves with H-vector along the y-direction; the tilde sign signifies time-dependent variables while the dot denotes differentiation with respect to time. The normalization in Eq. (1) is chosen such that q̃(t) ≡ Q̃(t)/C ₀, where Q̃ is the charge and C ₀ is the linear value of the capacitance, i.e. the normalized charge has the units of Volts.

Fig. 1 (a) Orientation of the unit cell; (b) Equivalent effective circuit model; (c) Simulations geometry; (d) Second order effective nonlinear susceptibility χ ⁽²⁾(2ω) calculated from Eq. (3).

Download Full Size | PDF

The voltage V_D can be expanded in a Taylor series in terms of the normalized charge according to V_D (q̃) = q̃ + aq̃ ² + bq̃ ³ where the Taylor coefficients a and b depend on the particular mechanism of nonlinearity [16]. In the case of a VLSRR- based medium, with varactors incorporated in the capacitive gaps of the SRRs as shown in Fig. 1(b), the Taylor coefficients can be found from the known nonlinear dependence of the capacitance of the varactor on the applied voltage [4, 8] as $a = \frac{K}{2 V_{p}}$ , $b = \frac{K (2 K - 1)}{6 V_{p}^{2}}$ , with the intrinsic potential V_p and the power coefficient K defined by the specifications of the varactor. The perturbation solution to Eq. (1) leads to the following expressions [8] for the linear and the second order effective susceptibilities characterizing the metacrystal:

χ_{y}^{(1)} (ω) = \frac{F ω^{2}}{D (ω)}

χ_{y y y}^{(2)} (ω_{r}; ω_{n}, ω_{m}) = - i a \frac{ω_{0}^{4} (ω_{n} + ω_{m}) ω_{n} ω_{m} μ_{0} A F}{D (ω_{n}) D (ω_{m}) D (ω_{n} + ω_{m})},

where ω_r ≡ ω_n + ω_m, each of the indices n and m taking values between ±Λ, Λ being the total number of distinct waves incident on the medium. The resonant denominator is defined as

D (ω) \equiv ω_{0}^{2} - ω^{2} - i γ ω

,

F \equiv ω_{0}^{2} N A^{2} C_{0} μ_{0}

is the amplitude factor in the expression for the linear susceptibility, and conventional notation is adopted for the arguments of the nonlinear susceptibility in which the first frequency term is the sum of the subsequent frequency arguments [9]. As seen from Eq. (2), the parameters F, ω ₀, and γ define the linear susceptibility, in qualitative agreement with the detailed studies [13]. These parameters can be obtained through standard retrieval methods from experimentally obtained or simulated transmission and reflection coefficients [14], thus not requiring the knowledge of the particular factors entering their equivalent expressions. We note that the parameters obtained in this manner account for all interactions between the inclusions. Once the linear parameters are known, the second and higher-order susceptibility terms can be easily and accurately determined since the same F, ω ₀, and γ enter these expressions; only the additional parameters strictly related to the nonlinear properties of the embedded material or device are required for the full description.

Equation (1) represents the nonlinear oscillator response of a metacrystal inclusion. For a dilute medium, the magnetization is approximately related to the magnetic moment as M̃(t) = Nm̃(t), where N is the volume density of moments, and the magnetic dipole moment of the effective circuit encompassing the effective area S is related to the oscillating charge as $\tilde{m} (t) = \frac{d \tilde{Q} (t)}{d t} A = A C_{0} \frac{d \tilde{q} (t)}{d t}$ . Using the above relations in Eq. (1), the latter can be rewritten as

{\ddot{\tilde{M}}}_{y} + γ {\dot{\tilde{M}}}_{y} + ω_{0}^{2} {\tilde{M}}_{y} = - F {\ddot{\tilde{H}}}_{y} - α {\tilde{M}}_{y} \int {\tilde{M}}_{y} d t

where

α \equiv \frac{2 ω_{0}^{2} a}{N A C_{0}}

and we include only the second order nonlinear response in the expansion of the nonlinear voltage V_D. The ”microscopic” equation of motion for a single inclusion has thus been converted into the “macroscopic” equation of motion for the effective medium polarization, which can then be used in conjunction with Maxwell’s equations to arrive at a complete description of the system.

We now consider propagation of a plane wave through a 1 cm sample of the VLSRR medium shown in Fig. 1, with unit cell size of 1 cm, ring diameter of 8 mm, and with the width of the metal strip and the gap size of 1 mm each. For the numerical simulations, we employ the finite-element time-domain (FETD) method available in the software package COMSOL. Our time domain simulations account simultaneously for all parametric interactions described by the nonlinear oscillator equation. By contrast, a frequency domain approach would usually require a separate equation for each participating frequency; obtaining a closed system of such coupled differential equations including all the processes is usually not feasible, and even consideration of more than a single nonlinear process complicates the system and intrinsically neglects many of the interactions. The geometry for the FETD simulations consists of a 1 cm thick and 6 cm wide layer of the effective nonlinear medium with the above parameters, with the input and the output ports located at 14.5 cm on each side of the medium, as shown in Fig. 1(c); the parameters correspond to the experimental transmission line used later in the experiment. Due to the field-invariance along the wavefront, it is sufficient to consider the two-dimensional case with the electric vector of the incident transverse electric (TE) polarized wave directed along the normal to the plane of the simulations.

In order to couple Eq. (4) to Maxwell’s equations, we use the standard TE wave propagation solver provided by COMSOL RF Module, coupled with a user-defined differential equation describing temporal evolution of the local magnetic response (Eq. (4)). We note that the Finite elements time domain method has been implemented in the present approach, rather than the regularly used finite difference time domain procedure [18]. The additional degrees of freedom describing the magnetization vector are introduced only within the domain representing the slab. Perfect magnetic conductor boundary conditions are used for the top and bottom boundaries, scattering boundary conditions — for the incoming and outgoing ports in the TE wave equation, while the homogeneous Neumann boundary condition is employed on all boundaries of the slab domain for the magnetization equation (Eq. (4)). After solving the problem for the fields H⃗ (r⃗,t), M⃗(r⃗,t), we extract the amplitude of the H_y component oscillating at the second harmonic frequency from the output port by Fourier transforming the steady-state portion of the time-dependent solution calculated over 128 wave periods with 64-point per period resolution, and compare with the measured second harmonic generated from an equivalent slab of the VLSRR metamaterial.

For the experiments, we use a transmission line apparatus shown in the inset in Fig. 2(a) that is similar to that in Refs. [8, 12] but having the geometry optimized for improved transmission at harmonic frequencies. We first perform a measurement of linear S-parameters by using Agilent 8364B network analyzer and applying a low (−15 dBm) power on the sample. Fitting the spectral shape of the linear experimental S-parameters by simulating the linear response of the unit cell and applying the standard retrieval procedure [14], we obtain values of F = 0.215, γ = 1.57 × 10⁸ s⁻¹, and ω ₀ = 2π × 0.816 GHz. We employ the Skyworks SMV1231 varactors having K = 0.8, V_p = 1.5 V according to the specifications [16,17], thus leading to the value of $α = - 4.14 \times 10^{9} \frac{m}{s A}$ . In subsequent SH measurements, the Agilent E8267D signal generator and Agilent E440A spectrum analyzer were used with the same transmission line set-up. The plane-wave assumption used in the simulations is assumed to approximate the transverse electromagnetic (TEM) mode propagating in the transmission line: the particular waveguide geometry here confines the incident and scattered waves, such that for a sample that fills the cross section, the geometry equates to that of a plane-wave scattering from an infinite sheet of medium [19]. An example of field distribution inside the transmission line before and after the material, obtained by the full-wave CST Microwave studio simulations is shown in Fig. 2(b), thus justifying the retrieval procedure that relies on the plane wave assumption.

Fig. 2 (a) Field distribution inside the transmission line at the crossection before the sample; the inset shows the transmission line structure with the material inside. (b) Field distribution after the sample inside the transmission line.

Download Full Size | PDF

Comparison between the experimental and analytically determined SHG for several values of the input power are presented in Fig. 3(a). The indicated powers correspond to the output values from the signal generator. In estimating the actual power incident on the metacrystal, we account for the fact that, according to the simulations, only about 49% of the power entering the port is concentrated between the metal plates where the TEM approximation is valid, and we assumed the 4.5 dB loss coming from the cables and the connector before the medium as well as 1 dB loss for the first harmonic and 1.3 dB loss for the second harmonic after the medium, according to the measurements. The corrected incident power P is related to the real-valued amplitude of the corresponding magnetic field as $H_{y} = \sqrt{2 P / Z_{0} S}$ [12], where Z ₀ = 377 Ω is the impedance of vacuum and S = 18 cm² is the crossectional area between the plates of the transmission line. We use the latter field value as the excitation field in the simulations. The figure reveals very good agreement between analytic theory and measurements for incident powers up to 7 dBm. The slightly higher numerical value of the SHG at 7 dBm occurs as a result of dissipative currents that increase as the diode bias increases. The presence of this power-dependent loss mechanism is not accounted for in the present nonlinear oscillator model. At larger power levels, higher losses lead to a more noticeable increase in the spectral width of the resonance and, with the present model, require a separate fitting procedure to determine a corrected value of γ. While this can be overcome by including a nonlinear resistance in the oscillator model, the appearance of this current at high powers is specific to using the single-gap VLSRR-based medium, and is not expected to be present with most other mechanisms of introducing a nonlinear response, e.g. by using a nonlinear dielectric inclusion. Very good agreement for all the powers below the values where a significant dissipative current takes place demonstrates the validity of the present nonlinear oscillator model. We note that in the present study the system became limited by the appearance of the dissipative current before it could reach the depleted pump regime. The developed approach, however, can be used at arbitrary power levels, thus paving the way for the analysis of nonlinear metamaterials in the depleted pump regime

Fig. 3 (a) Experimental amplitude of the H_y field component at the second harmonic, versus the simulation results; (b) Magnetization at the second harmonic obtained by the numerical solution to NL oscillator model (black lines) and analytically (red lines), for several values of the driving field (c) Same as in (b), with the numerical solution obtained from Comsol simulations, for the lower three values of the excitation field shown in (b). The grey gradient color areas indicate the range of the second-order magnetization amplitude that could be produced with a variation of the fundamental driving magnetic field shown on the inset.

Download Full Size | PDF

The nonlinear oscillator equation accounts simultaneously for all the included nonlinear parametric interactions, while more approximations are introduced by the perturbation approach solution leading to Eq. (3) [8]. We next verify the accuracy of this expression by comparing the value of the nonlinear magnetization $M_{y}^{(2)} (2 ω)$ extracted from the simulations with the one obtained using Eq. (3) according to $M_{y}^{(2)} (2 ω) = χ_{y y y}^{(2)} (2 ω, ω, ω) H_{y}^{2} (ω)$ . While the reflections of the input field from the boundaries of the metacrystal due to the impedance mismatch were included in the previous simulations emulating the experimental configuration, they are not accounted for in the perturbation approach solution. Therefore, to compare with the theoretical predictions, we extract the actual value of the fundamental magnetic field H_y inside the slab from the simulations, average it over the length of the slab for each excitation frequency, and use the resulting values in the above analytical expression for the magnetization as an approximation for the driving magnetic field. We note that due to the resonant response of the medium, the reflected power and hence the field inside the medium strongly depend on the frequency of the applied field, hence the driving field varies with frequency in this case, even with a constant input from the port. For completeness, we also solve Eq. (4) using ODE45 solver provided by MATLAB, assuming a constant driving magnetic field for all frequencies taken to be equal to the averaged H_y amplitude obtained from the simulations at the frequency f = 0.818 GHz. The averaged H_y amplitude is determined for each of the excitation powers used in the FETD simulations.

The comparison between the numerical magnetization obtained by solving the nonlinear oscillator equation in MATLAB and the analytical one is shown in Fig. 3(b). The three lower curves correspond the excitation powers of 5, 6, and 7 dBm used in FETD simulations in Figs. 3(a) and 3(c). We see a very good agreement for all excitation field values. As power increases, the shift in the resonant frequency towards lower values becomes pronounced in the numerical solution resulting from the contribution of second-order nonlinearity to the third-order nonlinear response [8, 9], which, in turn, leads to the power-tuning of the resonant frequency of the magnetization at the second harmonic. These effects are not accounted for in the expression for the second order susceptibility given by Eq. (3), hence the analytical curves have the same resonant frequency for all excitation fields. Also, the amplitude difference between the numerical and the analytical curves becomes pronounced at higher field values, staying however below 5% for the all the excitation field values considered.

Similar conclusions follow from Fig. 3(c) presenting the numerical results of magnetization distribution within the slab versus the analytical predictions. A small decrease in the fundamental magnetic field amplitude along the slab due to absorption loss is shown in the inset in Fig. 3(c). Since even a small field variation leads to a rapid change in the amplitude of a produced magnetization near the resonance, we indicated the possible variation in the magnetization amplitude by the grey gradient areas on the graph. Also, since the driving magnetic field is different in this case for every frequency used in the simulations, the resonance for each input power is wider than in Fig. 3(b), which, in turn, leads to a larger difference in the resonance positions between the numerical and analytical curves, as well as to a slightly misplaced resonance of the analytical magnetization from the assumed value of 0.818 GHz. The overall prediction of the amplitude of the magnetization at the SH is still very good within 10% for all power values considered, suggesting that the theory can be used with good accuracy for the quantitative prediction of the nonlinear properties of the effective medium. We also note that, according to Eq. (3), the value of χ ⁽²⁾(2ω) is many orders of magnitude larger than in most optical materials (χ ⁽²⁾(2ω) ∼ 10 as seen in Fig. 1(d)), which is the consequence of both a strong nonlinearity of varactor elements and of a resonant response of the medium. It nevertheless agrees well with the experimental data, thus demonstrating quantitatively the large nonlinear enhancement in the present metamaterial geometry.

Within the range of validity, the presented results thus suggest the applicability of the model for the prediction, analysis, and possible engineering of the enhanced nonlinear response of a metacrystal.

Acknowledgments

We thank Alec Rose for help with testing the numerical code. This work was supported by the Air Force Office of Scientific Research(Contract No. FA9550-09-1-0562).

References and links

1. J. B. Pendry, A. J. Holden, D. J. Robbins, and W. J. Stewart, “Magnetism from conductors and enhanced nonlinear phenomena,” IEEE Trans. Microwave Theory Tech. 47, 2075–2084 (1999). [CrossRef]

2. M. Lapine, M. Gorkunov, and K. H. Ringhofer, “Nonlinearity of a metamaterial arising from diode insertions into resonant conductive elements,” Phys. Rev. E 67, 065601 (2003). [CrossRef]

3. I. V. Shadrivov, S. K. Morrison, and Y. S. Kivshar, “Tunable split-ring resonators for nonlinear negative-index metamaterials,” Opt. Express 14, 9344–9349 (2006). [CrossRef] [PubMed]

4. B. Wang, J. Zhou, T. Koschny, and C. M. Soukoulis, “Nonlinear properties of split-ring resonators,” Opt. Express 16, 16058–16063 (2008). [CrossRef] [PubMed]

5. I. V. Shadrivov, A. B. Kozyrev, D. W. van der Weide, and Y. S. Kivshar, “Tunable transmission and harmonic generation in nonlinear metamaterials,” Appl. Phys. Lett. 93, 161903 (2008). [CrossRef]

6. M. W. Klein, M. Wegener, N. Feth, and S. Linden, “Experiments on second- and third- harmonic generation from magnetic metamaterials,” Opt. Express 15, 5238–5247 (2007). [CrossRef] [PubMed]

7. A. B. Kozyrev and D. W. van derWeide, “Nonlinear left-handed transmission line metamaterials,” J. Phys. D: Appl. Phys. 41, 173001 (10pp) (2008).

8. E. Poutrina, D. Huang, and D. R. Smith, “Analysis of nonlinear electromagnetic metamaterials,” N. J. Phys. 12, 093010 (2010). [CrossRef]

9. R. W. Boyd, Nonlinear Optics, 3rd ed. (Academic Press, 2008).

10. S. Larouche, A. Rose, E. Poutrina, D. Huang, and D. R. Smith, “Experimental determination of the quadratic nonlinear magnetic susceptibility of a varactor-loaded split ring resonator metamaterial,” Appl. Phys. Lett. 97, 011109 (2010). [CrossRef]

11. S. Larouche and D. R. Smith, “A retrieval method for nonlinear metamaterials,” Opt. Commun. 283, 1621–1627 (2010). [CrossRef]

12. D. Huang, E. Poutrina, and D. R. Smith, “Analysis of the power dependent tuning of a varactor-loaded metamaterial at microwave frequencies,” Appl. Phys. Lett. 96, 104104 (3pp) (2010).

13. J. Garcia-Garcia, F. Martin, J. D. Baena, R. Marques, and L. Jelinek , “On the resonances and polarizabilities of split ring resonators,” J. Appl. Phys. 98, 033103 (2005). [CrossRef]

14. D. R. Smith, S. Schultz, P. Markos, and C. M. Soukoulis, “Determination of effective permittivity and permeability of metamaterials from reflection and transmission coefficients,” Phys. Rev. B 65, 195104 (2002).

15. E. Poutrina, S. Larouche, and D. R. Smith, “Parameteric oscillator based on a single-layer resonant metamaterial,” Opt. Commun. 283, 1640–1646 (2010). [CrossRef]

16. V.V. Migulin, V.I. Medvedev, E. R. Mustel, and V.N. Parygin, Theory of Oscillations. Mir Publishers, Moscow, 1983.

17. Skyworks, “Skyworks smv123x series hyperabrupt junctiontuning varactors,” Technical report, Skyworks, September 2009. Data sheet. http://pdf1.alldatasheet.com/datasheet-pdf/view/155290/SKYWORKS/SMV1231-079.html).

18. I. Shadrivov, N. Zharova, A. Zharov, and Y. Kivshar, “Nonlinear transmission and spatiotemporal solitons in metamaterials with negative refraction,” Opt. Express 13, 1291–1298 (2005). [CrossRef] [PubMed]

19. B. Popa and S. Cummer, “Compact dielectric particles as a building block for low-loss magnetic metamaterials,” Phys. Rev. Lett. 100, 207401 (2008). [CrossRef] [PubMed]

Nonlinear oscillator metamaterial model: numerical and experimental verification

Abstract

1. Introduction

2. Analytical approach and numerical simulations

Acknowledgments

References and links

Cited By

Figures (3)

Equations (4)

Optics Express