Giant Kerr nonlinearity and low-power gigahertz solitons via plasmon-induced transparency

We propose a method to enhance Kerr nonlinearity and realize low-power gigahertz solitons via plasmon-induced transparency (PIT) in a new type of metamaterial, which is constructed by an array of unit cell consisting of a cut-wire and a pair of varactor-loaded split-ring resonators. We show that the PIT in such metamaterial can not only mimic the electromagnetically induced transparency in coherent three-level atomic systems, but also exhibit a crossover from PIT to Autler-Townes splitting. We further show that the system suggested here also possess a giant third-order nonlinear susceptibility and may be used to create solitons with extremely low generation power. Our study raises the possibility for obtaining strong nonlinear effect of gigahertz radiation at very low intensity based on room-temperature metamaterials.

transparency window becomes wider and deeper as d is reduced. These phenomena are clear indications of PIT-ATS crossover and will be analyzed in detail below.
The dynamics of the bright and dark modes in the unit cell at the position r = (x, y, z) can be described by the coupled Lorentz oscillator model 4,7 γ ω κ where q 1 and q 2 are respectively the displacement from the equilibrium position of oscillators 26 (the dot over q j (j = 1, 2) denotes time derivative), with γ 1 and γ 2 respectively their damping rates; ω 0 = 2π × 32 GHz and ω 0 + Δ are respectively resonance frequencies of the bright and dark modes (γ 2 ≪ γ 1 ≪ ω 0 ); g is the parameter indicating the coupling strength of the bright mode with the incident radiation E; parameter κ denotes the coupling strength between the CW and SRR-pair. The last two terms on the left hand of Eq. (2) are provided by the hyperabrupt tuning varactors mounted onto gaps of the SRRs [18][19][20] . Thus the metamaterial structure suggested here is a coupled anharmonic oscillator system with the nonlinear coefficients α and β and driven by the electric field E. Note that when writing Eq. (1) and (2) we have taken a coordinate system in which the electric (magnetic) field E (H) is along y (x) direction and wavevector k f is along z direction. We assume the frequency ω f of the incident radiation is near ω 0 (the natural frequency of the oscillator q 1 in linear regime). Evolution of the electric field is governed by the Maxwell equation Here N the density of unit cells, e is the unit charge, and χ ( ) D 1 is the optical susceptibility of the background material (which is assumed to be linear). Note that the resonance of the two SRRs has nearly no contribution to the average field E since the asymmetric field distribution at the two SRRs cancels each other when PIT occurs 4 . In addition, since the wavelength of the incident field (8.5 mm) is much larger than the thickness of each unit cell (10 μm), the electric field E seen by each meta-atom (i.e. the unit cell) is nearly homogeneous. Thus we can take an "electric-dipole approximation", widely used in atomic physics and quantum optics 27 , to investigate the dynamics of the present system.

PIT-ATS crossover.
To treat them analytically, we assume where q dj , q fj , q sj , q tj (j = 1, 2) are respectively amplitudes of the direct-current, fundamental, second harmonic, third harmonic waves of the jth oscillator, with k 0 (ω 0 ) the wavenumber (frequency) of the fundamental wave; E d , E f , E s , E t are respectively amplitudes of the direct-current, fundamental, second harmonic, third harmonic waves of the electric field, with k f ≈ k 0 (ω f ), k s (ω s ), and k t (ω t ) respectively the wavenumber (frequency) of the fundamental, second harmonic, and third harmonic waves. Under rotating-wave and slowly varying envelope approximations, from Eqs (1) (2) and (3) we can obtain a series of equations for q μj and E μ (μ = d, f, s, t), which are presented in Methods.
We assume that the transverse spatial distribution of the radiation is large enough so that the diffraction effect (i.e. its dependence on x and y) of the system can be neglected. We use the standard method of multiple scales 28,29 to derive the envelope equation of the radiation field based on the asymptotic expansion of q μj and E μ (see Methods). The first-order solution of the asymptotic expansion reads Here F is a yet to be determined envelope function depending on the "slow" variables z 1 , z 2 and t 1 , δ = ω f − ω 0 is frequency detuning, and K is the linear dispersion relation given by Shown in panels (d), (e), (f) of Fig. 2 is the absorption spectrum Im(K) (the imaginary part of K) as a function of frequency. When plotting the figure we used the damping rates γ 1 ≈ 60 GHz and γ 2 ≈ 10 GHz, which are nearly independent of d, whereas κ decreases from 152.5 GHz at d = 0.02 mm to 69 GHz at d = 0.38 mm. We see that the analytical result (the lower part of Fig. 2) fits well with the numerical one (the upper part of Fig. 2).
To obtain a clear physical explanation on the variation of the PIT transparency window for different d shown in Fig. 2, we make a detailed analysis on the linear dispersion relation. Near the resonance point with A = κ 0 g/(2ω 0 ) and B = κ 2 /(2ω 0 ). The linear dispersion relation (7) is similar to that obtained in Λ -type three-level atomic systems 1 . The PIT condition of the system reads B 2 ≫ γ 1 γ 2 /4. To illustrate the interference effect between the bright and dark modes apparently, we make a spectrum decomposition on Im(K) using the technique developed recently in refs 30-32 for different coupling constant κ.
1. Weak coupling region (κ 2 < (γ 1 − γ 2 )ω 0 /2)): In this region, the absorption spectrum can be expressed as 2 . We see that Im(K) consists of two Lorentzians terms. Figure 3(a) shows Im(K) as a function of δ for κ = 69 GHz. The dashed (dotted-dashed) line is the result of the first (second) Lorentzian term in Eq. (8), which is negative (positive). Because the two Lorentzians terms have the same center position but opposite sign, their superposition gives a destructive interference between CW and the SRR-pair. As a result, a small dip in Im(K) curve (i.e. the solid line) appears. Such phenomenon belongs PIT in nature.

2.
Intermediate coupling region (κ 2 > (γ 1 − γ 2 )ω 0 /2)): In this case, we obtain and h = (γ 2 − γ 1 )/(8δ 0 ). We see that the absorption spectrum is made of four terms. The first two terms are Lorentzians with the same width but different center position, Scientific RepoRts | 5:13780 | DOi: 10.1038/srep13780 which are two resonance peaks belonging respectively to the CW and the SRR-pair. The dip between the two Lorentzians can be interpreted as a gap between two resonances, which is a typical character of ATS. The next two terms are interference terms. Because they lowers the dip formed by the first two terms, a destructive interference (i.e. PIT) occurs. Since both PIT and ATS occur simultaneously in the system, we assign such phenomenon as PIT-ATS crossover. Plotted in Fig. 3(b) is the result of various terms in Eq. (9) and the total absorption spectrum Im(K) for κ = 82 GHz.
3. Strong coupling region (κ 2 ≫ (γ 1 − γ 2 )ω 0 /2)): The spectrum decomposition of Im(K) is still given by Eq. (9), but the destructive interference effect contributed by the last two terms plays a negligible role. Figure 3(c) shows the result for κ = 152.5 GHz. We see that in the strong coupling region the absorption spectrum of the electric field displays mainly the character of ATS.
In Fig. 3(d) we show the "phase diagram" of the system, where three regions (i.e. PIT, PIT-ATS crossover, and ATS) are indicated clearly. We see that a transition from PIT to ATS indeed occurs when κ changes from small to large values, which provides a satisfactory explanation for the phenomenon observed in our numerical simulation presented in the upper part of Fig. 2. Giant Kerr nonlinearity and low-power solitons. In order to explore the property of the nonlinear propagation of the radiation field, we go to high-order approximations. The solvability condition at the second order of the asymptotic expansion is i[∂F/∂z 1 + (∂K/∂ω)∂F/∂t 1 ] = 0, which means that the probe-pulse envelope F travels with the group velocity V g = (∂K/∂ω) −1 . The second-order solution is given in Methods.
With the above result we proceed to third order. The divergence-free condition in this order yields the nonlinear equation for the radiation field envelope F where α ε α = −2 1 , α 1 ≡ Im(K) is the coefficient describing linear absorption, K 2 = ∂ 2 K/∂δ 2 is the coefficient describing group-velocity dispersion, and W is the coefficient describing Kerr effect. We have (ii) Im(χ (3) ), which contributes a nonlinear absorption to the radiation field, is much less than Re(χ (3) ) when the system works in PIT transparency window (i.e., δ takes the values from − 20 GHz to 20 GHz). Such suppression of the nonlinear absorption is also due to the PIT effect. Although χ (3) is originated from the hyperabrupt tuning varactors mounted onto the slits of the SRRs (i.e. α ≠ 0, β ≠ 0) which can be large at GHz frequency, the giant enhancement of χ (3) to the order of 10 −7 m 2 V −2 obtained above is due to the contribution of the PIT effect. The physical reason is that the system we consider is a resonant one and works under the PIT condition, which can not only largely suppress the linear and nonlinearity absorptions of the radiation field, but also greatly strengthen the Kerr effect of the system, similar to the giant enhancement of the Kerr nonlinearity occurred in resonant atomic systems working under EIT condition 1 . Note that since Eqs (1) and (2) have quadratic nonlinearity resulting from the property of hyperabrupt tuning varactors, one is able to get a non-zero second-order nonlinear susceptibility χ (2) , which can also be enhanced via PIT effect. However, in the present work we do not discuss χ (2) process (e.g. second harmonic generation) for which a phase matching condition is required, an interesting topic that will be considered elsewhere.
Based on the above results and returning to the original variables, we obtain the NLS equation where τ = t − z/V g and U = ε F exp(− α 1 z). Generally, Eq. (12) has complex coefficients and hence it is a Ginzburg-Landau equation. However, due to the PIT effect the imaginary part of the complex coefficients can be made much smaller than their real part. Thus Eq. (12) can be converted into the dimensionless form are typical pulse length, dispersion length, absorption length, and amplitude of the pulse, respectively (the tilde above the quantity means taking real part). Notice that for forming solitons L D has been set to equal typical nonlinearity length If d 0 is small (this is the case we have here), the term d 0 u can be taken as a perturbation. One can use the perturbation method for solitons 33 to solve the Eq. (13) to obtain a single-soliton solution under the perturbation. The result is given by where η, ς, σ 0 and φ 0 are real free parameters determining the amplitude (as well as the width), propagating velocity, initial position, and initial phase of the soliton, respectively. When taking η = 1/2, ς = σ 0 = φ 0 = 0 and noting that s = − z/L D , we obtain the electric field corresponding the above single-soliton solution , which describes a damped bright soliton traveling with velocity ∼ V g . For δ = 15 GHz and with other parameters given above, we obtain numerical values of the coefficients in Eq. (12), given by α 1 = 0.0278 cm −1 , K 2 = (5.02 + 0.57i) × 10 −23 cm −1 s 2 , W = (4.62 + 0.72i) × 10 −2 C 2 kg −2 cm −3 s 4 . We see that, as expected, the imaginary part of these coefficients are indeed much smaller than their real part. By taking τ 0 = 1.5 × 10 −11 s, we obtain U 0 = 2.2 V cm −1 , L D = L NL = 4.48 cm. In order to have such dispersion and nonlinearity lengthes, 20 layers of the unit cell array are needed to form the gigahertz soliton. With the above system parameters one has L A = 36 cm, and hence d 0 = 0.12, i.e. the dissipation in the system can be indeed taken as a perturbation. The physical reason of the long absorption length for the soliton propagation in the present system is also due to the PIT effect.
Note that the diffraction length of the system is given by with R ┴ the transverse spatial width of the radiation pulse. With the system parameters given above and if taking R ┴ = 1.4 cm, we have L diff = 43 cm, which is much larger than the dispersion length and nonlinearity length (4.48 cm). Thus in this case the diffraction effect of the system is indeed negligible.
The power associated with the soliton is given by Poynting's vector integrated over the cross-section of the radiation in the transverse directions, i.e. P = ∫(E × H) · e z dS, where e z is the unit vector in the propagation direction 29 . We obtain the peak power of the soliton ε τ = / ( ) ∼ ∼ P cn S K W 2 p max 0 0 2 0 2 , here ω = + / ∼ n c K 1 p 0 is the refractive index and S 0 is the area of the intensity distribution of the radiation field in the transverse directions. Using the above parameters and taking S 0 = 6 × 10 −4 m 2 , the peak power for generating the soliton is found to be = P 568 mW max , which means that to generate the soliton in present system a very low input power is needed. This is a drastically contrast to the case in conventional media such as optical fibers, where pico-or femto-second laser pulses are needed to reach a very high peak power (usually at the order of several hundred kW) to stimulate enough nonlinearity for the formation of solitons.
The stability of the gigahertz soliton is tested by using numerical simulation. Figure 5(a) shows the result of the radiation intensity |E/U 0 | 2 of the soliton as a function of t/τ 0 and z/L D . The solution is obtained by numerically solving Eq. (12) with the complex coefficients taken into account, with the initial condition given by E(0, t)/U 0 = sech(t/τ 0 ). We see that the shape of the soliton undergoes no apparent deformation during propagation.
The collision between two solitons is also studied numerically, with the result shown in Fig. 5(b). The initial condition is given by E(0, t)/U 0 = sech[(t − 3.0)/τ 0 ] + 1.2sech[1.2(t + 3.0)/τ 0 ]. We see that the both solitons can resume their original shapes after the collision, indicating that solitons in the PIT-metamaterial are robust during interaction.

Discussion
We have suggested a method for enhancing Kerr nonlinearity and realizing low-power gigahertz solitons via PIT in a new type of metamaterial, which can be constructed by an array of unit cell consisting of a cut-wire and a pair of varactor-loaded SRRs. We have demonstrated that the PIT in such system can not only mimic the electromagnetically induced transparency in coherent three-level atomic systems, but also exhibit a crossover from PIT to Autler-Townes splitting. We have further demonstrated that our Scientific RepoRts | 5:13780 | DOi: 10.1038/srep13780 system may possess a giant third-order nonlinear susceptibility and can be used to produce solitons with extremely low generation power.
The research results reported here raises the possibility for obtaining strong nonlinear effect of gigahertz radiation at very low intensity based on room-temperature metamaterials. The method presented in this work can be used to investigate other types of PIT systems, such as that proposed in ref. 6; furthermore, they can also be generalized to terahertz and even optical frequency regimes if related nonlinear elements are available. Asymptotic expansion. We make the asymptotic expansion where ε is a dimensionless small parameter characterizing the amplitude of the incident electric field. All quantities on the right-hand side of the expansion are assumed as functions of the multi-scale variables z l = ε l z (l = 0, 1, 2) and t l = ε l t (l = 0, 1). Note that with the expansion based on the multiple-scale variables given here, a short-wavelength approximation is implied, i.e. the expansion is valid only for a propagation of a weak nonlinear pulse with finite K and finite δ (see Eq. (30) below). One can also consider other possibilities (e.g. long-wavelength approximation) that will result in different envelope equations, which will be discussed elsewhere. Substituting the expansion (28) and (29) into the Eqs (16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27) and comparing the coefficients of ε l (l = 1, 2, 3, …), we obtain a chain of linear but inhomogeneous equations which can be solved order by order.
At the third order (l = 3), a solvability condition results in the equation for F, i.e. Eq. (10).