Analogue of double-$\Lambda$-type atomic medium and vector dromions in a plasmonic metamaterial

TL;DR

This paper demonstrates that a specially designed plasmonic metamaterial can mimic atomic systems with complex interactions, exhibit nonlinear effects like plasmon induced transparency, and support stable high-dimensional vector plasmonic dromions with low power requirements.

Contribution

It introduces a classical analogue of a double-Λ atomic medium in a plasmonic metamaterial and shows the emergence of vector dromions with low power and robustness.

Findings

Metamaterial exhibits plasmon induced transparency similar to atomic systems.

System supports high-dimensional vector plasmonic dromions with low power.

Giant second- and third-order Kerr nonlinearities are achieved.

Abstract

We consider an array of the meta-atom consisting of two cut-wires and a split-ring resonator interacting with an electromagnetic field with two polarization components. We prove that such metamaterial system can be taken as a classical analogue of an atomic medium with a double- $\Lambda$-type four-level configuration coupled with four laser fields, exhibits an effect of plasmon induced transparency (PIT), and displays a similar behavior of atomic four-wave mixing (FWM). We demonstrate that when nonlinear varactors are mounted onto the gaps of the split-ring resonators the system can acquire giant second- and third-order Kerr nonlinearities via the PIT and a longwave-shortwave interaction. We also demonstrate that the system supports high-dimensional vector plasmonic dromions [i.e. (2+1)-dimensional plasmonic solitons with two polarization components, each of which has a coupling between a longwave and a shortwave], which have very low generation power and are robust during propagation. Our work gives not only a plasmonic analogue of the FWM in coherent atomic systems but also provides the possibility for obtaining new type of nonlinear polaritons in plasmonic metamaterials.

Equations

$${\hat{H}}_{\mathrm{int}}=-\hslash \sum _{j=1}^4{\Delta }_j|j⟩⟨j|-\hslash \left[{\Omega }_{p1}|3⟩⟨1|+{\Omega }_{p2}|4⟩⟨1|+{\Omega }_{c1}|3⟩⟨2|+{\Omega }_{c2}|4⟩⟨2|+H.c.\right],$$

$${\Omega }_{p1}{\Omega }_{c2}-{\Omega }_{p2}{\Omega }_{c1}=0.$$

$$\left(i\frac{\partial }{\partial t}+d_{31}\right){\sigma }_{31}+{\Omega }_{c1}{\sigma }_{21}+{\Omega }_{p1}=0,$$

$$\left(i\frac{\partial }{\partial t}+d_{41}\right){\sigma }_{41}+{\Omega }_{c2}{\sigma }_{21}+{\Omega }_{p2}=0,$$

$$\left(i\frac{\partial }{\partial t}+d_{21}\right){\sigma }_{21}+{\Omega }_{c1}^{\ast }{\sigma }_{31}+{\Omega }_{c2}^{\ast }{\sigma }_{41}=0,$$

$$i\left(\frac{\partial }{\partial z}+\frac{1}{c}\frac{\partial }{\partial t}\right){\Omega }_{p1}+{\kappa }_{13}{\sigma }_{31}=0,$$

$$i\left(\frac{\partial }{\partial z}+\frac{1}{c}\frac{\partial }{\partial t}\right){\Omega }_{p2}+{\kappa }_{14}{\sigma }_{41}=0,$$

$$K_a^{\pm }\left(\omega \right)=\frac{\omega }{c}+\frac{-\left({\kappa }_{14}{D}_3+{\kappa }_{13}{D}_4\right)\pm \sqrt{{\left({\kappa }_{14}{D}_3-{\kappa }_{13}{D}_4\right)}^2+4{\kappa }_{13}{\kappa }_{14}|{\Omega }_{c1}{\Omega }_{c2}{|}^2}}{2\left[|{\Omega }_{c1}{|}^2\left(\omega +d_{41}\right)+|{\Omega }_{c2}{|}^2\left(\omega +d_{31}\right)-\left(\omega +d_{21}\right)\left(\omega +d_{31}\right)\left(\omega +d_{41}\right)\right]},$$

$$\frac{{\partial }^2{q}_1}{\partial t^2}+{\gamma }_1\frac{\partial q_1}{\partial t}+{\omega }_1^2{q}_1-{\kappa }_1{q}_3=g_1{E}_x,$$

$$\frac{{\partial }^2{q}_2}{\partial t^2}+{\gamma }_2\frac{\partial q_2}{\partial t}+{\omega }_2^2{q}_2-{\kappa }_2{q}_3=g_2{E}_y,$$

$$\frac{{\partial }^2{q}_3}{\partial t^2}+{\gamma }_3\frac{\partial q_3}{\partial t}+{\omega }_3^2{q}_3-{\kappa }_1{q}_1-{\kappa }_2{q}_2=0,$$

$$\kappa_{2}{\cal E}_{y0}-\kappa_{1}{\cal E}_{x0}=0.\$$

$${\nabla }^2{E}_{x(y)}-\frac{1}{c^2}\frac{{\partial }^2{E}_{x(y)}}{\partial t^2}=\frac{1}{{\epsilon }_0{c}^2}\frac{{\partial }^2{P}_{x(y)}}{\partial t^2},$$

$$\left(i\frac{\partial }{\partial t}+d_1\right){\tilde{q}}_1+\frac{{\kappa }_1}{2{\omega }_p}{\tilde{q}}_3+\frac{g_1}{2{\omega }_p}{\mathcal{E}}_x=0,$$

$$\left(i\frac{\partial }{\partial t}+d_2\right){\tilde{q}}_2+\frac{{\kappa }_2}{2{\omega }_p}{\tilde{q}}_3+\frac{g_2}{2{\omega }_p}{\mathcal{E}}_y=0,$$

$$\left(i\frac{\partial }{\partial t}+d_3\right){\tilde{q}}_3+\frac{{\kappa }_1}{2{\omega }_p}{\tilde{q}}_1+\frac{{\kappa }_2}{2{\omega }_p}{\tilde{q}}_2=0,$$

$$i\left(\frac{\partial }{\partial z}+\frac{n_{\mathrm{D}}}{c}\frac{\partial }{\partial t}\right){\mathcal{E}}_x+{\kappa }_0{\tilde{q}}_1=0,$$

$$i\left(\frac{\partial }{\partial z}+\frac{n_{\mathrm{D}}}{c}\frac{\partial }{\partial t}\right){\mathcal{E}}_y+{\kappa }_0{\tilde{q}}_2=0,$$

$$K_m^{\pm }\left(\omega \right)=\frac{n_{\mathrm{D}}}{c}\omega +{\kappa }_0\frac{-\left(R_1{g}_{f2}+R_2{g}_{f1}\right)\pm \sqrt{{\left(R_1{g}_{f2}-R_2{g}_{f1}\right)}^2+4{\kappa }_{f1}^2{\kappa }_{f2}^2{g}_{f1}{g}_{f2}}}{2\left[{\kappa }_{f1}^2\left(\omega +d_2\right)+{\kappa }_{f2}^2\left(\omega +d_1\right)-\left(\omega +d_3\right)\left(\omega +d_1\right)\left(\omega +d_2\right)\right]},$$

$${\mathcal{E}}_x(z,t)=\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }d\omega \left[F_0^{+}{e}^{i(K_m^{+}z-\omega t)}+F_0^{-}{e}^{i(K_m^{-}z-\omega t)}\right],$$

$${\mathcal{E}}_y(z,t)=\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }d\omega \left[G^{+}{F}_0^{+}{e}^{i(K_m^{+}z-\omega t)}+G^{-}{F}_0^{-}{e}^{i(K_m^{-}z-\omega t)}\right],$$

$${\mathcal{E}}_x(z,t)=\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }d\omega \frac{G^{+}{e}^{i(K_m^{-}z-\omega t)}-G^{-}{e}^{i(K_m^{+}z-\omega t)}}{G^{+}-G^{-}}{\tilde{\mathcal{E}}}_x(0,\omega ),$$

$${\mathcal{E}}_y(z,t)=\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }d\omega \frac{G^{+}{G}^{-}}{G^{+}-G^{-}}\left[e^{i(K_m^{-}z-\omega t)}-e^{i(K_m^{+}z-\omega t)}\right]{\tilde{\mathcal{E}}}_x(0,\omega ),$$

$${\mathcal{E}}_x(z,t)=\frac{G_0^{+}{\mathcal{E}}_x(0,{\tau }_{-})e^{i{K}_0^{-}z}-G_0^{-}{\mathcal{E}}_x(0,{\tau }_{+})e^{i{K}_0^{+}z}}{G_0^{+}-G_0^{-}},$$

$${\mathcal{E}}_y(z,t)=\frac{G_0^{+}{G}_0^{-}}{G_0^{+}-G_0^{-}}\left[{\mathcal{E}}_x(0,{\tau }_{-})e^{i{K}_0^{-}z}-{\mathcal{E}}_x(0,{\tau }_{+})e^{i{K}_0^{+}z}\right],$$

$${\mathcal{E}}_x(z,t)=\frac{G_0^{-}}{G_0^{-}-G_0^{+}}{\mathcal{E}}_x(0,{\tau }_{+})e^{i{K}_0^{+}z},$$

$${\mathcal{E}}_y(z,t)=\frac{G_0^{-}{G}_0^{+}}{G_0^{-}-G_0^{+}}{\mathcal{E}}_x(0,{\tau }_{+})e^{i{K}_0^{+}z}.$$

$$\eta (L)=\frac{|G_0^{+}{G}_0^{-}{|}^2}{|G_0^{+}-G_0^{-}{|}^2}|\exp (i{K}_0^{+}L){|}^2.$$

$$\frac{{\partial }^2{q}_3}{\partial t^2}+{\gamma }_3\frac{\partial q_3}{\partial t}+{\omega }_3^2{q}_3-{\kappa }_2{q}_1-{\kappa }_2{q}_2+\alpha q_3^2+\beta q_3^3=0,$$

$$i\frac{\partial F_{+}}{\partial z_2}-\frac{K_2^{+}}{2}\frac{{\partial }^2{F}_{+}}{\partial t_1^2}+\frac{c}{2{\omega }_p{n}_{\mathrm{D}}}\left(\frac{{\partial }^2}{\partial x_1^2}+\frac{{\partial }^2}{\partial y_1^2}\right)F_{+}+\frac{{\omega }_p{R}_0}{2c{n}_{\mathrm{D}}}{\chi }_{+}^{(2)}{Q}_{+}{F}_{+}+\frac{{\omega }_p}{2c{n}_{\mathrm{D}}}{\chi }_{+}^{(3)}{\left|F_{+}\right|}^2{F}_{+}=0,$$

$${\chi }_{+}^{(2)}=$$

$${\chi }_{+}^{(3)}=$$

$$\times \left[\left(\frac{4{\alpha }^2{\omega }_1^2{\omega }_2^2}{{\omega }_1^2{\omega }_2^2{\omega }_3^2-{\omega }_2^2{\kappa }_1^2-{\omega }_1^2{\kappa }_2^2}+\frac{2{\alpha }^2{H}_1{H}_2}{H_1{H}_2{H}_3-H_2{\kappa }_1^2-H_1{\kappa }_2^2}\right)-3\beta \right].$$

$$\left(\frac{{\partial }^2}{\partial x_1^2}+\frac{{\partial }^2}{\partial y_1^2}\right)Q_{+}+\left[{\left(\frac{1}{V_g^{+}}\right)}^2-{\left(\frac{1}{V_p^{+}}\right)}^2\right]\frac{{\partial }^2{Q}_{+}}{\partial t_1^2}-\frac{{\chi }_{+}^{(2)}}{c^2}\frac{{\partial }^2|F_{+}{|}^2}{\partial t_1^2}=0,$$

$$\frac{1}{V_p^{+}}=\frac{n_{\mathrm{D}}}{c}+\frac{N_{m}e}{2{\epsilon }_{0}c{n}_{\mathrm{D}}}\frac{(g_1{X}_2+g_2{X}_1)+\sqrt{(g_1{X}_2-g_2{X}_1{)}^2+4{\kappa }_1^2{\kappa }_2^2{g}_1{g}_2}}{{\omega }_1^2{\omega }_2^2{\omega }_3^2-{\omega }_1^2{\kappa }_2^2-{\omega }_2^2{\kappa }_1^2},$$

$${\chi }_{\mathrm{eff}}^{(3)}={\chi }_{+}^{(3)}+\frac{R_0{\left[{\chi }_{+}^{(2)}\right]}^2}{c^2\left[\left(\frac{1}{V_p^{+}}\right)-{\left(\frac{1}{V_g^{+}}\right)}^2\right]},$$

$$i\frac{\partial u}{\partial s}+\left(\frac{{\partial }^2}{\partial {\xi }^2}+g_{d0}\frac{{\partial }^2}{\partial {\eta }^2}+g_{d1}\frac{{\partial }^2}{\partial {\tau }^2}\right)u+2g_1|u{|}^{2}u+g_{2}vu=0,$$

Full text

I Introduction

Electromagnetically induced transparency (EIT), a very intriguing phenomenon occurring in atomic gases, has been intensively investigated due to its interesting physical properties and promising practical applications. The basic mechanism of EIT is the existence of a destructive quantum interference effect between two pathways of atomic transitions induced by a control laser field, through which the absorption of a probe laser field can be largely cancelled Fleischhauer .

In recent years, there are tremendous efforts seeking for the classical analogue of EIT in solid systems, including coupled resonators Alzar ; Spring ; Harden , electric circuits Alzar ; Harden ; Souza , optomechanical devices Weis ; Kronwald , whispering-gallery-mode microresonators Peng , and various metamaterials Zhangs ; Papa ; Tassin ; Gu ; han ; Chen0 ; Dong ; NLiu ; Liu ; TNakanishi ; Sun ; Chen ; Mats ; NakanishiPRApp ; bai ; Bai2 ; Bai3 . Especially, the plasmonic analogue of EIT, called plasmon-induced transparency (PIT) Zhangs ; Papa ; Tassin , has become a very important platform for exploring EIT-like physical properties of plasmonic polaritons and for designing new types of metematerials.

Similar to EIT, PIT is resulted from a destructive interference effect between wideband bright and narrowband dark modes in artificial atoms (called meta-atoms). A typical character of PIT is the opening of a deep transparency window within broadband absorption spectrum, together with a steep dispersion and greatly reduced group velocity of plasmonic polaritons. PIT metamaterials can work in different frequency regions (including micro Papa and terahertz Tassin ; Gu ; Dong waves, infrared and visible radiations Zhangs ; NLiu ; han ), and may be used to design novel, chip-scale plasmonic devices (including highly sensitive sensors Dong ; Chen0 , optical buffers Gu ; Liu , and ultrafast optical switches Gu , etc.) in which the radiation damping can be significantly eliminated, very intriguing for practical applications.

However, the PIT in plasmonic metamaterials reported up to now Zhangs ; Papa ; Tassin ; Gu ; han ; Chen0 ; Dong ; NLiu ; Liu ; TNakanishi ; Sun ; Chen ; Mats ; NakanishiPRApp ; bai ; Bai2 ; Bai3 is only for the classical analogue of the simplest atomic EIT, i.e. the one occurring in a coherent three-level atom gas resonantly interacting with two laser fields. We know that atoms possess many (energy) levels, quantum interference effect may occur in atomic systems with level number larger than three and the number of laser fields larger than two Fleischhauer . In fact, in the past two decades the EIT has been extended into the atomic systems with various multi-level configurations, such as four-level systems of double $\Lambda$-type Lukin ; Korsunsky1 ; Korsunsky2 ; Merriam ; DengPRA1 ; Kang ; Ying ; DengPRA3 ; YingSoliton ; Kang1 ; IteYu , tripod-type Petrosyan ; Rebic ; Beck , Y-type Joshi ; GaoJY ; Kha ; LiuYM , five-level systems of M-type Ott ; Mat ; Rebic1 ; HangC , and six-level systems of double-tripod-type Rus ; Lee , etc. Thus it is natural to ask the question: Is it possible to get a classical analogue of the atomic EIT with the level number more than three in a metamaterial?

In this article, we give a positive answer for the above question. The metamaterial we consider is assumed to be an array of meta-atoms [see Fig. 2(a)], i.e. the unit cells consisting of two cut-wires (CWs) and a split-ring resonator (SRR) [see Fig. 2(b)], interacting with an electromagnetic (EM) field with two polarization components. We show that such plasmonic metamaterial system may be taken as a classical analogue of an atomic medium with a double-$\Lambda$-type four-level configuration coupled with four (two probe and two control) laser fields [see Fig. 1(a) ], exhibits an effect of PIT and displays a similar behavior of atomic four-wave mixing (FWM).

Based on this classical analogue, we further show that, if nonlinear varactors are mounted onto the gaps of the SRRs, the system can acquire giant second- and third-order Kerr nonlinearities many orders of magnitude larger than conventional nonlinear optical media. Using a method of multiple scales, we derive coupled envelope equations, which include dispersion, diffraction, and the Kerr nonlinearities and govern the evolution of the two polarization components of the EM field. We demonstrate that the system supports a new type of nonlinear plasmonic polaritons, i.e. high-dimensional vector plasmonic dromions, which are (2+1)-dimensional plasmonic solitons with two polarization components. Each polarization component has a coupling between a longwave and a shortwave, which have very low generation power and are robust during propagation. The results presented here not only gives a close metamaterial analogue of the EIT and FWM in multi-level atomic systems, useful to illustrate and find novel interference and nonlinear properties in solid systems, but also provides a way to obtain new type of plasmonic polaritons via suitable design of plasmonic metamaterials.

The main body of the article is arranged as follows. In Sec II, we give a simple introduction of the four-level atomic model allowing EIT and FWM, describe the metamaterial model, and show the similarity between the two models. The propagation of linear plasmonic polaritons in the metamaterial is discussed in detail. In Sec. III, we derive the coupled nonlinear envelope equations and present the vector plasmonic dromion solutions when the nonlinear varactors are mounted onto the gaps of the SRRs. Lastly, in Sec IV we give a discussion and a summary of our work. Details of some calculating results are given in five appendixes.

II EIT-based atomic FWM and its metamaterial analogue

II.1 EIT-based FWM in a double-$\Lambda$-type four-level atomic system

For a detailed comparison with the metamaterial model presented in the next subsection, we first give a brief introduction on a lifetime-broadened atomic gas with a double-$\Lambda$-type four-level configuration, shown in Fig. 1(a).

In this system, two weak probe laser fields with central angular frequencies $\omega_{p1}$ and $\omega_{p2}$ and wavevectors ${\bf k}_{p1}$ and ${\bf k}_{p2}$ drive respectively the transitions $|1\rangle\leftrightarrow|3\rangle$ and $|1\rangle\leftrightarrow|4\rangle$, and two strong control laser fields with central angular frequencies $\omega_{c1}$ and $\omega_{c2}$ and wavevectors ${\bf k}_{c1}$ and ${\bf k}_{c2}$ drive respectively the transitions $|2\rangle\leftrightarrow|3\rangle$ and $|2\rangle\leftrightarrow|4\rangle$. The total electric fields in this system is given by ${\bf E}={\bf e}_{p1}{\cal E}_{p1}\exp[i(k_{p1}z-\omega_{p1}t)]+{\bf e}_{p2}{\cal E}_{p2}\exp[i(k_{p2}z-\omega_{p2}t)]+{\bf e}_{c1}{\cal E}_{c1}\exp[i(k_{c1}z-\omega_{c1}t)]+{\bf e}_{c2}{\cal E}_{c2}\exp[i(k_{c2}z-\omega_{c2}t)]+{\rm c.c.}$, where ${\bf e}_{jn}$ and ${\cal E}_{jn}$ ($j=p,c;n=1,2$) are respectively the unit vector denoting the polarization direction and the envelope of the corresponding laser field. Note that for simplicity all the laser fields are assumed to be injected in the same (i.e. $z$) direction (which is also useful to suppress Doppler effect). Under electric-dipole approximation and rotating-wave approximation (RWA), the Hamiltonian of the system in interaction picture reads

[TABLE]

where $\Delta_{1}=0$; $\Delta_{3}=\omega_{p1}-(E_{3}-E_{1})/\hbar$, $\Delta_{2}=\omega_{p1}-\omega_{c1}-(E_{2}-E_{1})/\hbar$, and $\Delta_{4}=(\omega_{p1}-\omega_{c1}+\omega_{c2})-(E_{4}-E_{1})/\hbar$ are respectively one-, two-, and three-photon detunings, with $E_{l}$ the eigenenergy of the atomic state $|l\rangle$ ($l=1,2,3,4$); $\Omega_{p1}=({\bf e}_{p1}\cdot{\bf p}_{31}){\cal E}_{p1}/\hbar$, $\Omega_{p2}=({\bf e}_{p2}\cdot{\bf p}_{41}){\cal E}_{p2}/\hbar$, $\Omega_{c1}=({\bf e}_{c1}\cdot{\bf p}_{32}){\cal E}_{c1}/\hbar$, and $\Omega_{c2}=({\bf e}_{c2}\cdot{\bf p}_{42}){\cal E}_{c2}/\hbar$ are respectively the half Rabi frequencies of the probe and the control laser fields, with ${\bf p}_{jl}$ the electric dipole moment related to the transition $|j\rangle\leftrightarrow|l\rangle$. The Hamiltonian (1) allows three bright states and one dark state note00 . The dark state reads $|\psi_{\rm dark}\rangle=(\Omega_{c1}|1\rangle-\Omega_{p1}|2\rangle)/\sqrt{|\Omega_{p1}|^{2}+|\Omega_{c1}|^{2}}$, which is a superposition of only the two lower states $|1\rangle$ and $|2\rangle$ and has a zero eigenvalue. The condition yielding the dark state is Korsunsky2

[TABLE]

The dynamics of the atoms is governed by the optical Bloch equation $i\hbar\left(\partial/\partial t+\Gamma\right)\sigma=[{\hat{H}}_{\rm int},\sigma]$, where $\sigma$ is a $4\times 4$ density matrix, $\Gamma$ is a $4\times 4$ decoherence (relaxation) matrix describing spontaneous emission and dephasing. The explicit expression of the Bloch equation is given in Appendix A. We assume that initially the probe fields are absent, thus for substantially strong control fields the atoms are populated in the ground state $|1\rangle$. The solution of the Bloch equation reads $\sigma_{11}=1$ and all other $\sigma_{jl}$ are zero.

When the two weak probe fields are applied, the ground state $|1\rangle$ is not depleted much. In this case, the Bloch equation reduces to

[TABLE]

with $d_{j1}=\Delta_{j}+i\gamma_{j1}$ with $\gamma_{j1}=\Gamma_{1j}/2$ ($j=2,3,4$). Equations (3a)-(3c) describe the dynamics of three coupled harmonic oscillators note01 , where $\sigma_{31}$ and $\sigma_{41}$ are bright oscillators due to their direct coupling to the probe fields $\Omega_{p1}$ and $\Omega_{p2}$, but $\sigma_{21}$ is a dark oscillator because it has no direct coupling to any of the two probe fields.

The dynamics of the probe fields is governed by the Maxwell equation $\nabla^{2}{\bf E}-(1/c^{2})\partial^{2}{\bf E}/{\partial t^{2}}=1/(\varepsilon_{0}c^{2})\partial^{2}{\bf P}/{\partial t^{2}}$. Here the polarization intensity is given by ${\bf P}=N_{0}[\sigma_{31}e^{i(k_{p1}{z}-\omega_{p1}t)}+\sigma_{41}e^{i(k_{p2}{z}-\omega_{p2}t)}+{\rm c.c.}]$, with $N_{0}$ the atomic density. Under a slowly-varying envelope approximation (SVEA), the Maxwell equation reduces to

[TABLE]

with $\kappa_{13}=N_{0}|{\bf e}_{p1}\cdot{\bf p}_{13}|^{2}\omega_{p1}/(2\hbar\varepsilon_{0}c)$ and $\kappa_{14}=N_{0}|{\bf e}_{p2}\cdot{\bf p}_{14}|^{2}\omega_{p2}/(2\hbar\varepsilon_{0}c)$. For simplicity, we assume the two control fields are strong enough and thus have no depletion during the evolution of the probe fields; additionally, the diffraction effect is negligible, which is valid for the probe fields having large transverse size.

It is easy to understand the basic feature of the propagation of the probe fields through solving the Maxwell-Bloch (MB) equations (3) and (4) with $\sigma_{l1}$ ($l=1,2,3$) and $\Omega_{pj}$ ($j=1,2$) proportional to the form $\exp[i(K_{a}z-\omega t)]$ note02 . We obtain

[TABLE]

with $D_{3}=|\Omega_{c1}|^{2}-(\omega+d_{21})(\omega+d_{31})$ and $D_{4}=|\Omega_{c2}|^{2}-(\omega+d_{21})(\omega+d_{41})$. We see that the MB equations allow two normal modes, with the linear dispersion relations given by $K_{a}^{+}$ and $K_{a}^{-}$, respectively.

Fig. 1(b) shows ${\rm Im}(K_{a}^{+})$ [i.e. the imaginary part of $K_{a}^{+}$] as a function of $\omega$ for $\Omega_{c1}=\Omega_{c2}=20\,{\rm MHz}$ (red dashed line) and $\Omega_{c1}=\Omega_{c2}=60\,{\rm MHz}$ (green dashed-dot line). When plotting the figure, $\Delta_{j}$ ($j=1,2,3$) are set to be zero, and realistic parameters from ${}^{87}{\rm Rb}$ atoms are taken, given by $\Gamma_{13}=\Gamma_{23}=\Gamma_{14}=\Gamma_{24}=16\,{\rm MHz}$, $\kappa_{13}=\kappa_{14}=1\times 10^{10}\,{\rm cm}^{-3}$ Rb87 . We see that a transparency window is opened in the profile of ${\rm Im}(K_{a}^{+})$ near $\omega=0$; the transparency window becomes larger when the control fields are increased. The opening of the transparency window (called EIT transparency window) is due to the EIT effect contributed by the control fields. The blue solid cure in the figure is ${\rm Im}(K_{a}^{-})$ as a function of $\omega$, which however has always a large absorption peak near $\omega=0$ irrespective of the value of the control fields. Below, for convenience we shall call the normal mode with the linear dispersion relation $K_{a}^{+}$ ($K_{a}^{-}$) as EIT-mode (non-EIT-mode).

The double-$\Lambda$-type four-level system can be used to describe a resonant FWM process in atomic systems Fleischhauer ; Merriam ; DengPRA1 ; Kang ; Ying . The first laser field (i.e. the control field tuned to the $|2\rangle\leftrightarrow|3\rangle$ transition with the half Rabi frequency $\Omega_{c1}$) and the second laser field (i.e. the probe field tuned to the $|1\rangle\leftrightarrow|3\rangle$ transition with the half Rabi frequency $\Omega_{p1}$) can adiabatically establish a large atomic coherence of the Raman transition, described by the off-diagonal density matrix element $\sigma_{21}$. The third laser field, i.e. the control field tuned to the $|2\rangle\leftrightarrow|4\rangle$ transition with the half Rabi frequency $\Omega_{c2}$, can mix with the coherence $\sigma_{21}$ to generate a fourth field with the half Rabi frequency $\Omega_{p2}$ resonant with the $|1\rangle\leftrightarrow|4\rangle$ transition. For details, see Refs.Fleischhauer ; Merriam ; DengPRA1 ; Kang ; Ying and references therein.

II.2 Metamaterial analogue of the double-$\Lambda$-type four-level atomic system

We now seek for a possible classical analogue of the above four-level atomic model by using a metamaterial, which is assumed to be an array [Fig. 2(a)] of unit cells (i.e. meta-atoms) [Fig. 2(b)] consisting two CWs (indicated by “A” and “B”) and a SRR.

The CW A and CW B are, respectively, positioned along $x$ and $y$ direction, while the SRR is formed by a square ring with a gap at the center of each side. 20-${\rm\mu m}$-thick copper forming the CWs and the SRR is etched on a substrate with a height of $H=1.5\,{\rm mm}$. Geometrical parameters of the meta-atom are $L_{x}=L_{y}=8\,{\rm mm}$, $w_{s}=1.2\,{\rm mm}$, $w_{g}=0.6\,{\rm mm}$, and $w_{b}=1.2\,{\rm mm}$ note1000 .

We assume that an incident gigahertz radiation ${\bf E}={\bf e}_{x}E_{x}+{\bf e}_{y}E_{y}$ [with $E_{j}={\cal E}_{j0}e^{-i\omega_{p}t}+{\rm c.c.}$ ($j=x,y$)] is collimated on the array of the meta-atoms, with polarization component $E_{x}$ ($E_{y}$) parallel to the CW A (CW B). In order to understand the EM property of the system, a numerical simulation is carried out by using the commercial finite difference time domain software package (CST Microwave Studio) and probing the imaginary part of the radiative field amplitude at the center of the end facet of the CW $A$ [red arrow in Fig. 2(b)] Zhangs . The blue dashed lines in Fig. 2(c) are normalized absorption spectrum as a function of the incident wave frequency by taking the excitation condition ${\cal E}_{y0}=-{\cal E}_{x0}$ for $d_{x}=d_{y}=3.4\,{\rm mm}$ (first panel), and ${\cal E}_{y0}=-{\cal E}_{x0}$ for $d_{x}=d_{y}=4.0\,{\rm mm}$ (second panel). Here and below, ${\cal E}_{x0}$ and ${\cal E}_{y0}$ are taken to be real for simplicity. Fig. 2(d) shows the normalized absorption spectrum under the excitation condition ${\cal E}_{y0}={\cal E}_{x0}$ for $d_{x}=d_{y}=4.0\,{\rm mm}$. The red solid lines in the figure are analytical results obtained from ${\rm Im}(q_{10})$ based on Eq. (32) in Appendix B. We see that the absorption spectrum profile depends on excitation condition, which is quite different from the PIT absorption spectrum considered before.

The dependence on the excitation condition for the absorption spectrum can be briefly explained as follows. A sole CW in the meta-atoms is function as an optical dipole antenna and thus serves as a bright (or radiative) oscillator, which can be directly excited by the incident radiation. The surface current in an excited SRR can be clockwise or anticlockwise direction, indicating that there is no direct electric dipole coupling with the incident radiation and hence the SRR serves as a dark or trapped oscillator with long dephasing time note03 . For the excitation condition ${\cal E}_{y0}=-{\cal E}_{x0}$ [Fig. 2(c)], the surface current is cooperatively induced through the near-field coupling between SRR and CWs, resulting in a maximum enhancement of the dark-oscillator resonance and thus the substantial suppression of the absorption of the incident radiation, acting like a typical PIT metamaterial. However, for the excitation condition ${\cal E}_{y0}={\cal E}_{x0}$ [Fig. 2(d)], the surface current is suppressed due to an opposite excitation direction, leading to a complete suppression of the dark-oscillator resonance. As a result, the radiation absorption is significant (acting like a sole CW) and hence no PIT behavior occurs. For convenience, in the following we called the excitation mode under the ${\cal E}_{y0}=-{\cal E}_{x0}$ [Fig. 2(c)] as the PIT-mode, and the excitation mode under the ${\cal E}_{y0}={\cal E}_{x0}$ [Fig. 2(d)] as the non-PIT-mode (or absorption mode).

The dynamics of the bright oscillators (i.e. CW A and CW B) and dark oscillator (i.e. SRR) in the meta-atoms can be described by the coupled Lorentz equations Zhangs ; Gu ; Chen ; bai ; Bai2 ; Bai3

[TABLE]

where $q_{l}$ are displacements from the equilibrium position of the bright oscillators ($l=1,2$) and the dark oscillator ($l=3$), with $\gamma_{l}$ and $\omega_{l}$ note03 respectively the damping rate and the natural frequency of $l$th oscillator; $g_{1}$ ($g_{2}$) is the parameter describing the coupling between the CW A (CW B) and the $x$-polarization ($y$-polarization) component of the EM wave, and $\kappa_{1}$ ($\kappa_{2})$ is the parameter describing the coupling between CW A (CW B) and SRR. The numerical values of these coefficients can be obtained from the numerical result presented in Fig. 2, by using the method described in Appendix B.

Based on the solution given by Eq. (32), we deduce that, in the case of $\omega_{3}=\omega_{p}$, $\gamma_{3}=0$ and $g_{1}=g_{2}$, Eq. (6) allows a “dark state” (i.e. the state where both the bright oscillators are not excited, i.e. $q_{10}=q_{20}=0$) exists, if

[TABLE]

This “dark state” condition is equivalent to the one obtained in the four-level double-$\Lambda$-type atomic system, given by Eq. (2). Obviously, the PIT-mode shown in Fig. 2(c) corresponds to the case $\kappa_{2}=-\kappa_{1}$, where the minus symbol can be understood as a $\pi$-phase difference resulting in a cooperative coupling effect, which is assumed in all the numerical calculations carried out below.

The equation of motion of the EM wave is governed by the Maxwell equation

[TABLE]

where $P_{x(y)}=\varepsilon_{0}\chi_{\rm D}^{(1)}E_{x(y)}+N_{m}eq_{1(2)}$, with $N_{m}$ the unit-cell density, $e$ the unit charge, and $\chi_{\rm D}^{(1)}$ the optical susceptibility of the hosting material. We assume the distance between the meta-atoms is large so that the interaction between them can be neglected.

Assuming the central frequency of the incident radiation $\omega_{p}$ is near the natural frequencies of the Lorentz oscillators described by Eq. (6) note03 , a resonant interaction occurs between the incident radiation and these oscillators. To deal with the propagation problem of the plasmonic polaritons in the system analytically, we assume $E_{j}({\bf r},t)={\cal E}_{j}(z,t)e^{i\left(k_{p}z-\omega_{p}t\right)}+{\rm c.c.}$ and $q_{l}({\bf r},t)={\tilde{q}}_{l}(z,t)\exp[i(k_{p}z-\omega_{l}t-\Delta_{l}t)]+{\rm c.c.}$, where ${\cal E}_{j}(z,t)$ and ${\tilde{q}}_{l}(z,t)$ are slowly-varying envelopes and $\Delta_{l}=\omega_{p}-\omega_{l}$ is a small detuning. With this ansatz and under RWA, Eq. (6) is simplified into the reduced Lorentz equation

[TABLE]

with $d_{j}=\Delta_{j}+\gamma_{j}/2$. We see that the reduced Lorentz equation (9) describing the unit cell has the same form as the optical Bloch equation (3) describing the four-level double-$\Lambda$ atom. Consequently, each unit cell in the metamaterial is analogous to a four-level double-$\Lambda$-type atom in the atomic gas presented in the last subsection. That is to say, the unit cell is indeed a meta-atom, where the bright-oscillator excitation in the CW A (CW B) driven by ${\cal E}_{x}$ (${\cal E}_{y}$) is equivalent to the dipole-allowed transition $|1\rangle\leftrightarrow|3\rangle$ ($|1\rangle\leftrightarrow|4\rangle$) driven by the probe field $\Omega_{p1}$ ($\Omega_{p2}$), and the dark-oscillator excitation in the SRR is equivalent to the dipole-forbidden transition $|1\rangle\rightarrow|2\rangle$ in the four-level double-$\Lambda$-type atom. We also see that the coupling between the CW A (CW B) and the SRR, described by $\kappa_{1}$ ($\kappa_{2}$), is equivalent to the control field $\Omega_{c1}$ ($\Omega_{c2}$) driven the atomic transition $|2\rangle\leftrightarrow|3\rangle$ ($|2\rangle\leftrightarrow|4\rangle$).

Under SVEA, the Maxwell equation in the metamaterial reads

[TABLE]

with $\kappa_{0}=N_{m}e\omega_{p}/(2\varepsilon_{0}cn_{\rm D})$, $n_{\rm D}=\sqrt{1+\chi_{\rm D}^{(1)}}$. Obviously, Eq. (10) has similar structure as Eq. (4). Thus, a complete correspondence between the MB equations (3) and (4) described the four-level double-$\Lambda$-type atomic gas and the Maxwell-Lorentz (ML) equations (9) and (10) described the plasmonic metamaterial is established.

The propagation feature of a plasmonic polariton in the metamaterial can be obtained by assuming all quantities in the ML equations (9) and (10) proportional to $\exp[i(K_{m}z-\omega t)]$. It is easy to get the linear dispersion relation

[TABLE]

where $R_{j}=\kappa_{fj}^{2}-(\omega+d_{j})(\omega+d_{3})$, with $\kappa_{fj}=\kappa_{j}/(2\omega_{p})$ and $g_{fj}=g_{j}/(2\omega_{p})$ ($j=1,2$). As expected, the metamaterial system allows two normal modes with the linear dispersion relation respectively given by $K_{m}^{+}$ and $K_{m}^{-}$. In fact, $K_{m}^{+}$ ($K_{m}^{-}$) is a PIT-mode (non-PIT-mode) of the system, as explained below.

The character of the above two normal modes can be clearly illustrated by plotting $K_{m}^{+}$ and $K_{m}^{-}$ as functions of $\omega$. Shown in Fig. 3(a)

are Im($K_{m}^{+}$) (blue dashed line) and Re($K_{m}^{+}$) (red solid line) for $\kappa_{2}=-\kappa_{1}=50\,{\rm GHz}^{2}$ (first panel; corresponding to $d_{x}=d_{y}=4.0\,{\rm mm}$) and $\kappa_{2}=-\kappa_{1}=250\,{\rm GHz}^{2}$ (second panel; corresponding to $d_{x}=d_{y}=3.4\,{\rm mm}$), respectively. When plotting the figure, the system parameters are taken from Appendix B, and additional parameters are chosen by $\kappa_{0}=10^{10}\,{\rm kg/(cm\cdot s^{2}\cdot C)}$ Mikhail and $\Delta_{j}=0\,(j=1,2,3)$. We see that ${\rm Im}(K_{m}^{+})$ displays a transparency window (called PIT transparency window) near $\omega=0$, analogous to the EIT transparency window in ${\rm Im}(K_{a}^{+})$ of the four-level double-$\Lambda$-type atomic system [red dashed line and green dashed-dot line in Fig. 1(b)]. The steep slope of ${\rm Re}(K_{m}^{+})$ indicates a normal dispersion and a slow group velocity of the plasmonic polariton. As the coupling strength between the CWs and the SRR gets larger (i.e. the separations $d_{x}$ and $d_{y}$ is reduced), the PIT transparency window becomes wider and deeper, and the slope of ${\rm Re}(K_{m}^{+})$ gets flatter. The opening of the PIT transparency window is attributed to the destructive interference between the two bright oscillators and the dark oscillator through cooperative near-field coupling. Shown in Fig. 3(b) is the imaginary (red solid line) and the real (blue dashed line) of $K_{m}^{-}$, which is nearly independent on the coupling constant $\kappa_{1}$ ($\kappa_{2}=-\kappa_{1}$). We see that ${\rm Im}(K_{m}^{-})$ has a single, large absorption peak and ${\rm Re}(K_{m}^{-})$ has an abnormal dispersion near $\omega=0$, analogous to ${\rm Im}(K_{a}^{-})$ of the double-$\Lambda$-type atomic system [blue solid line in Fig. 1(b)].

II.3 Propagation of linear plasmonic polaritons via an analogous FWM process of atomic system

As indicated above, the meta-atoms in the present metamaterial system are analogous to the four-level atoms with the double-$\Lambda$-type configuration, and hence an analogous resonant FWM phenomenon for the plasmonic polaritons is possible. That is to say, if initially only one polarization-component of the EM wave (e.g. $x$-component) is injected into the metamaterail, a new polarization-component (e.g. $y$-component) will be generated through two equivalent control fields (i.e. the couplings between the SRR and CWs, described by $\kappa_{1}$ and $\kappa_{2}$). To illustrate this, we present the solution of the ML equations (9) and (10)

[TABLE]

which can be obtained by using Fourier transform Ying ; LiHJ . Here $G^{\pm}=[-A\pm(A^{2}+4\kappa_{f1}^{2}\kappa_{f2}^{2}g_{f1}g_{f2})^{1/2}]/(2\kappa_{f1}\kappa_{f2}g_{f2})$ with $A=R_{1}g_{f2}-R_{2}g_{f1}$, and $F_{0}^{\pm}$ is the initial amplitude of the normal mode $K_{m}^{\pm}$ determined by given excitation condition. We assume initially only the $x$-component of the EM field in input to the system, i.e. the initial condition for the EM field is given by ${\cal E}_{x}(0,t)\neq 0$, ${\cal E}_{y}(0,t)=0$. By Eq. (12) we have

[TABLE]

where ${\cal\tilde{E}}_{x}(0,\omega)=\int^{+\infty}_{-\infty}dt{\cal E}_{x}(0,t)e^{i\omega t}$. For simplicity, we consider the adiabatic regime where the power series of $K_{m}^{\pm}$ and $G^{\pm}$ on $\omega$ converge rapidly. By taking $K_{m}^{\pm}=K_{0}^{\pm}+\omega/V_{g}^{\pm}+O(\omega^{2})$ and $G^{\pm}=G_{0}^{\pm}+O(\omega)$, we readily obtain

[TABLE]

where $\tau_{\pm}=t-z/V_{g}^{\pm}$, with $V_{g}^{\pm}\equiv(\partial K_{m}^{\pm}/\partial\omega)^{-1}|_{\omega=0}$ being the group-velocity of the normal mode $K_{m}^{\pm}$.

The conversion efficiency of the FWM is given by $\eta(L)\equiv\int_{-\infty}^{+\infty}dt|{\cal E}_{y}(L,t)/{\cal E}_{x}(0,t)|^{2}$, where $L$ is the medium length. For the case $\kappa_{2}=-\kappa_{1}$, one has ${\rm Im}(K_{0}^{-})\gg{\rm Im}(K_{0}^{+})$, which means that the $K_{m}^{-}$ mode decays away rapidly during propagation and hence can be safely neglected. Then Eq. (14) is simplified as

[TABLE]

We see that the $x$-and $y$-polarization components of the EM wave have matched group velocity $V_{g}^{+}$. The expression of the FWM conversion efficiency reduces into

[TABLE]

Shown in Fig. 4

is the FWM conversion efficiency $\eta$ as a function of the dimensionless optical depth $(\kappa_{0}g_{f1}/\gamma_{1})L$ for $\Delta_{1}=\Delta_{2}=0$ (blue dashed line) and for $\Delta_{1}=\Delta_{2}=5\gamma_{1}$ (red solid line). When plotting this figure, we have set $\Delta_{3}=0$ and $\gamma_{3}\approx 0$ in order for a better analogue to the atomic system. The influence of $\gamma_{3}$ can be effectively reduced by introducing a gain element into the gaps of the SRRs (see the discussion in Sec. IV). From the figure, we see that for the case of exact resonance (i.e. $\Delta_{1}=\Delta_{2}=0$), the FWM efficiency $\eta$ increases and rapidly saturates to $25\%$ when the dimensionless optical depth $(\kappa_{0}g_{f1}/\gamma_{1})L\approx 5$ (i.e. $L\approx 0.9\,{\rm cm}$), indicating a unidirectional energy transmission from ${\cal E}_{x}$ to ${\cal E}_{y}$. For the case of far-off resonance (i.e. $\Delta_{1}=\Delta_{2}=5\gamma_{1}$), the FWM efficiency displays a damped oscillation in the interval $0<(\kappa_{0}g_{f1}/\gamma_{1})L<250$, indicating a back-and-forth energy exchange between ${\cal E}_{x}$ and ${\cal E}_{y}$; eventually the efficiency reach to the steady-state value $25\%$ when $(\kappa_{0}g_{f1}/\gamma_{1})L\geq 300$ (see the inset). Interestingly, the value of the FWM conversion efficiency may reach to $\eta\approx 76\%$ at $(\kappa_{0}g_{f1}/\gamma_{1})L\approx 15$ (i.e. $L\approx 3\,{\rm cm}$).

III Vector plasmonic dromions in the PIT metamaterial

Note that when deriving Eq. (10), the diffraction effect has been neglected, which is invalid for the plasmonic polaritons with small transverse size or long propagation distance; furthermore, because of the highly resonant (and hence dispersive) character inherent in the PIT metamaterial, the linear plasmonic polaritons obtained above inevitably undergo significant distortion during propagation. Hence it is necessary to seek the possibility to obtain a robust propagation of the plasmonic polaritons in the PIT metamaterial. One way to solve this problem is to make the PIT system work in a nonlinear propagation regime.

In recent years, nonlinear metamaterials have attracted much attention due to their potential applications (see Ref. Lapine and references therein). One suitable way to design a nonlinear PIT metamaterial in microwave and lower THz ranges is to use nonlinear insertions onto the meta-atoms Shadrivov ; Wang . Here, as suggested in Refs. Shadrivov ; Wang ; Bai2 , we assume the nonlinear insertion in the PIT metamaterial are varactor diodes, which are mounted onto the gaps of the SRRs Wang [see Fig. 2(b)].

III.1 Nonlinear envelope equations

Since the introduction of the nonlinear element onto the SRRs, Eq. (6c) should be replaced by Wang

[TABLE]

where $\alpha$ and $\beta$ are nonlinearity coefficients, described in Appendix E.

Due to the quadratic and cubic nonlinearities in Eq. (17), the input EM field (with only a fundamental wave) will generate longwave (rectification), and second harmonic components, i.e. $E_{l}({\bf r},t)={\cal E}_{dl}({\bf r},t)+[{\cal E}_{fl}({\bf r},t)e^{i(k_{p}z-\omega_{p}t)}+{\rm c.c.}]+[{\cal E}_{sl}({\bf r},t)e^{i\theta_{p}}+{\rm c.c.}]$ ($l=x,y$), with $\theta_{p}=(2k_{p}+\Delta k)z-2\omega_{p}t$ and $\Delta k$ a detuning in wavenumber. The oscillations of the Lorentz oscillators in the meta-atoms have the form $q_{j}({\bf r},t)=q_{dj}({\bf r},t)+[q_{fj}({\bf r},t)e^{i\theta_{j}}+{\rm c.c.}]+[q_{sj}({\bf r},t)e^{2i\theta_{j}}+{\rm c.c.}]$ ($j=1,2,3$), with $\theta_{j}=k_{j}z-\omega_{j}t-\Delta_{j}t$. Substituting these expressions into Eqs. (6a), (6b), (8), and (17), and adopting RWA and SVEA, we obtain a series of equations for the motion of $q_{\mu j}$ and ${\cal E}_{\mu l}$, listed in Appendix C.

We solve the equations for $q_{\mu j}$ and $E_{\mu l}$ by using the standard method of multiple scales Jeffery . Take the asymptotic expansion $q_{fj}=\epsilon q_{fj}^{(1)}+\epsilon^{2}q_{fj}^{(2)}+\cdots$, $q_{dj}=\epsilon^{2}q_{dj}^{(2)}+\cdots$, $q_{sj}=\epsilon^{2}q_{sj}^{(2)}+\cdots$, ${\cal E}_{fl}=\epsilon{\cal E}_{fl}^{(1)}+\epsilon^{2}{\cal E}_{fl}^{(2)}+\cdots$, and ${\cal E}_{dl}=\epsilon^{2}{\cal E}_{dl}^{(2)}+\cdots$ (here $\epsilon$ is a dimensionless small parameter characterizing the amplitude of the incident EM field), and assume all quantities on the right sides of the asymptotic expansion as functions of the multiscale variables Jeffery $z_{l}=\epsilon^{l}z\,(l=0,1,2)$ and $t_{l}=\epsilon^{l}t\,(l=0,1)$. Substituting the expansion into the equations for $q_{\mu j}$ and ${\cal E}_{\mu j}$ and comparing powers of $\epsilon$, we obtain a chain of linear but inhomogeneous equations which can be solved order by order.

The leading order [i.e. $O(\epsilon)$] solution reads ${\cal E}_{fx}^{(1)}=F_{+}e^{i\theta_{+}}$ and ${\cal E}_{fy}^{(1)}=G^{+}F_{+}e^{i\theta_{+}}$, where $\theta_{+}=K_{m}^{+}z_{0}-\omega t_{0}$ and $F_{+}$ is a slowly-varying envelope function to be determined in higher-order approximations. The expression of $G^{+}$ is given in Sec. II.3. Here we consider only the PIT (i.e. $K_{m}^{+}$) mode because the non-PIT (i.e. $K_{m}^{-}$) mode decays rapidly during propagation, as indicated in the last section. The solution for $q_{\mu j}^{(1)}$ is presented in Appendix D.

At the second order [i.e. $O(\epsilon^{2})$], a solvability condition yields $i[\partial/\partial z_{1}+(1/V_{g}^{+})\partial/\partial t_{1}]F_{+}=0$, where $V_{g}^{+}$ is the group-velocity of the fundamental wave. Solution for the longwave is ${\cal E}_{dx}^{(2)}=Q_{+}$ and ${\cal E}_{dy}^{(2)}=G^{+}Q_{+}$, with $Q_{+}$ the slowly-varying envelope to be determined yet. Explicit expressions of the solutions for other quantities at this order are listed in Appendix D.

At the third order [i.e. $O(\epsilon^{3})$], a solvability condition results in the equation

[TABLE]

where $K_{2}^{+}\equiv[\partial^{2}K_{m}^{+}/\partial\omega^{2}]|_{\omega=0}$ describes the group-velocity dispersion of the fundamental wave; $R_{0}$ is a coefficient characterizing the coupling between the fundamental and long waves, with the expression given in Appendix D; $\chi_{+}^{(2)}$ and $\chi_{+}^{(3)}$ are, respectively, the second-order and third-order nonlinear susceptibilities of the $K_{m}^{+}$ mode, with the form

[TABLE]

At the fourth order [i.e. $O(\epsilon^{4})$], a solvability condition results in the equation for the longwave $Q_{+}$

[TABLE]

here $V_{p}^{+}\equiv[n_{m}^{+}(0)/c]^{-1}$ is the phase-velocity of the longwave, defined by $n_{m}^{+}(0)=n_{m}^{+}|_{\omega_{p}=0,\omega=0}$ with $n_{m}^{+}(\omega;\omega_{p})=c[k_{p}(\omega_{p})+K_{m}^{+}(\omega;\omega_{p})]/(\omega_{p}+\omega)$. It is easy to obtain

[TABLE]

with $X_{j}=\omega_{j}\omega_{3}-\kappa_{j}^{2}$.

Under the PIT condition [i.e. $(\kappa_{j}/2\omega_{p})^{2}\gg\gamma_{j}\gamma_{3}/4$, $j=1,2$], the real parts of the nonlinearity susceptibilities $\chi_{+}^{(2)}$ and $\chi_{+}^{(3)}$ are greatly enhanced and their imaginary parts are much smaller than their real parts. Interestingly, $\chi_{+}^{(3)}$ can be further enhanced via a strong coupling between the longwave and the fundamental wave. This can be seen from the expression of the effective third-order nonlinearity susceptibility of the $K_{m}^{+}$ mode Newell

[TABLE]

obtained by neglecting the diffraction term in Eq. (20) and then plugging the derived expression for $Q_{+}$ into Eq. (18). The second term in Eq. (22) is due to the longwave-shortwave interaction, where one can clearly see that when the group-velocity of the shortwave is nearly equal to the phase-velocity of the longwaves (i.e. $V_{g}^{+}\approx V_{p}^{+}$), the effective third-order nonlinear susceptibility $\chi_{\rm eff}^{(3)}$ can be further enhanced, which can be realized by adjusting the coupling parameters $\kappa_{1}$ and $\kappa_{2}$. Compared with the conventional PIT-based metamaterials bai ; Bai2 with one bright and one dark oscillators, the present FWM-based metamaterial consisting double bright and one dark oscillators allows for a stronger bright-dark coupling, resulting in an enhanced Kerr nonlinearity. As a result, the coupling strength for the longwave-shortwave interaction in the present FWM-based metamaterial is $\sqrt{2}$ time larger than that in the conventional PIT-based metamaterials.

III.2 Vector plasmonic dromions

Equations (18) and (20) show that the self-interaction of the shortwave $F_{+}$ can stimulate the generation of the longwave $Q_{+}$ [Eq. (20)], and at the same time the longwave $Q_{+}$ has a back-action to the shortwave $F_{+}$ [Eq. (18)]. Such equations admit solutions describing the excitation of high-dimensional, two-component (vector) nonlinear plasmonic polaritons in the PIT metamaterial. To demonstrate this, we convert them into the dimensionless form

[TABLE]

where $u=\epsilon F_{+}/U_{0}$, $v=\epsilon^{2}Q_{+}/V_{0}$, $s=z/(2L_{\rm diff})$, $\tau=(t-z/V_{g}^{+})/\tau_{0}$, $\xi=x/R_{x}$, $\eta=y/R_{y}$, and $g_{d0}=(R_{x}/R_{y})^{2}$, $g_{d1}=L_{\rm diff}/L_{\rm disp}$, $g_{d2}=[R_{x}^{2}/\tau_{0}^{2}][(1/V_{p}^{+})^{2}-(1/V_{g}^{+})^{2}]$, $g_{1}=L_{\rm diff}/L_{\rm nln}$, $g_{2}=L_{\rm diff}[\omega_{p}/(n_{\rm D}c)]R_{0}V_{0}\chi_{+}^{(2)}$, $g_{3}=R_{x}^{2}U_{0}^{2}\chi_{+}^{(2)}/[c^{2}\tau_{0}^{2}V_{0}]$. Here $R_{x}\,(R_{y})$ is the typical radius of the incident EM field in the $x\,(y)$ direction; $\tau_{0}$ is the typical pulse duration of the probe field; $U_{0}\,(V_{0})$ is the typical amplitude of the longwave (shortwave) envelope; $L_{\rm disp}=-\tau_{0}^{2}/{\rm Re}(K_{2}^{+})$, $L_{\rm diff}\equiv n_{\rm D}\omega_{p}R_{x}^{2}/c$ and $L_{\rm nln}=(2n_{\rm D}c)/[\omega_{p}U_{0}^{2}{\rm Re}(\chi_{+}^{(3)})]$ are, respectively, the typical dispersion length, the typical diffraction length, and the typical nonlinearity length. Note that in obtaining Eq. (23), we have neglected the small imaginary parts of $\chi_{+}^{(3)}$ and $K_{2}^{+}$, which is reasonable under the PIT condition as discussed above.

In favor of the formation of plasmonic dromions, we take the following two assumptions. First, we shall assume $R_{x}\ll R_{y}$ and thus $g_{d0}\ll 1$ so that the original (3+1)-dimensional nonlinear problem can be reduced into a (2+1)-dimensional one. Second, we assume the contribution of the dispersion, diffraction and the nonlinearity effects are of the same level, which can be achieved by taking $L_{\rm diff}=L_{\rm disp}=L_{\rm nln}$ and thus we obtain $\tau_{0}=R_{x}\sqrt{-\omega_{p}n_{\rm D}{\rm Re}(K_{2}^{+})/c}$ and $U_{0}=[c/(\omega_{p}R_{x})]\sqrt{2/{\rm Re}[\chi_{+}^{(3)}]}$. In fact, these two assumptions can be realized by taking the realistic set of parameters, i.e. $\Delta_{1}=\Delta_{2}=-5\gamma_{1}$, $\kappa_{2}=-\kappa_{1}=4260\,{\rm GHz}^{2}$, $R_{x}=1.8\,{\rm cm}$, $R_{y}=10.2\,{\rm cm}$, $U_{0}=8.25\,{\rm V/cm}$, $V_{0}=1.3\,{\rm V/cm}$, $\tau_{0}=4.1\times 10^{-11}\,{\rm s}$, and hence we have $K_{2}^{+}=(-1.30+i0.10)\times 10^{-19}\,{\rm cm^{-2}s}$, $\chi_{+}^{(3)}=1.20\times 10^{-3}+i2.26\times 10^{-7}\,{\rm cm^{2}/V^{2}}$, $L_{\rm diff}=L_{\rm disp}=L_{\rm nln}=13.4\,{\rm cm}$, $g_{1}=g_{2}=g_{d1}=g_{d2}=1$, $g_{3}=4$, and $g_{d0}\ll 1$. For such case, Eq. (23) can be simplified into standard Davey-Stewartson-I (DS-I) equation

[TABLE]

with $v_{1}=v+2|u|^{2}$, where, for convenience, we have performed a 45-degree rotation of coordinates $\xi_{1}=(\xi+\tau)/\sqrt{2}$ and $\tau_{1}=(\xi-\tau)/\sqrt{2}$. The DS-I equation (24) can be exactly solved via the Hirota’s bilinear method Hirota and various dromion solutions can be obtained. A single-dromion solution reads Dromion $u=G/F$ and $v_{1}=V_{11}+V_{12}$ with $V_{11}=2\partial^{2}{\rm ln}(F)/\partial\xi_{1}^{2}$ and $V_{12}=2\partial^{2}{\rm ln}(F)/\partial\tau_{1}^{2}$, and

[TABLE]

where $\eta_{j}=(k_{jr}+ik_{ji})(r_{j}-r_{j0})+(\Omega_{jr}+i\Omega_{ji})s\,(j=1,2)$, with $(r_{1(0)},r_{2(0)})=(\xi_{1(0)},\tau_{1(0)})$, $\Omega_{jr}=-2k_{jr}k_{ji}$, $\Omega_{1i}+\Omega_{2i}=k_{1r}^{2}+k_{2r}^{2}-k_{1i}^{2}-k_{2i}^{2}$, and $\rho=2\sqrt{2k_{1r}k_{2r}(\gamma-1)}$. Here $k_{jr}$, $k_{ji}$, $r_{j0}$ and $\gamma$ are free real parameters.

Shown in Fig. 5(a) and Fig. 5(b) are, respectively, intensity distributions of the shortwave profile $|u|^{2}$ and the longwave profile $|v_{1}|^{2}$ as functions of $\xi_{1}$ and $\tau_{1}$ at $s=0$, when taking $\gamma=9$, $\xi_{10}=\tau_{10}=0$, $k_{1r}=k_{2r}=1$, and $k_{1i}=k_{2i}=0$. Obviously, the shortwave profile $|u|^{2}$ denotes a localized envelope function, which decays exponentially in all spatial directions [Fig. 5(a)]; the longwave profile $|v_{1}|^{2}$ denotes two interacting plane solitons (kinks), which decay in their respective traveling directions [Fig. 5(b)]. Such high-dimensional nonlinear excitation is called dromion Hirota .

We proceed with the investigation on the interaction between two plasmonic dromions. In doing so, we first integrate Eq. (24b), yielding

[TABLE]

where the nonzero boundary conditions $V_{11}|_{\tau_{1}\rightarrow-\infty}$ and $V_{12}|_{\xi_{1}\rightarrow-\infty}$ are given by Dromion

[TABLE]

In the numerical simulation, the integral Eq. (26) is conducted by fourth-order Runge-Kutta method with the boundary condition (27), and then Eq. (24a) with the obtained result $v_{1}$ is solved by using split-step Fourier method. Shown in Fig. 5(c1)-Fig. 5(c4) [Fig. 5(d1)-Fig. 5(d4)] are numerical results for the evolution of the shortwave $|u|^{2}$ (longwave $|v_{1}|^{2}$) during the collision between two dromions at $s\equiv z/(2L_{\rm diff})=0,1,2,3$, respectively. When doing the simulation we have taken a superposition of two dromion solutions as an initial input [i.e. Fig. 5(c1) and Fig. 5(d1)], and initial speeds and positions of the two dromions are setting to be $k_{1i}=-k_{2i}=1.8$ and $r_{10}=-r_{20}=3.2$. From the figure we see that two initial dromions become four dromions after the collision, which are, respectively, located around the four intersections of the longwave $v_{1}$, and these four dromions gradually separate and propagate almost stably, indicating that the collision between dromions is inelastic. The reason is that the four intersections of $v_{1}$ are the most attractive points in the entire region, which attract the EM wave intensities of the main peaks of the shortwaves $u$ during the collision, resulting in the appearance of four pulses for the shortwave located around the four intersections after the collision Dromion .

Within the forth-order approximation, the explicit expression for the EM field in the metamaterial takes the form

[TABLE]

When $u$ and $v$ are taken as the dromion solution given above, we obtain a vector plasmonic dromion since the EM field (28) has two polarization components, with each component a plasmonic dromion. Note that, different from the result in the scalar model considered before Bai2 , the polarization of the EM field obtained here can be actively selected by adjusting the separation between the CWs and SRR [i.e. $d_{x}$ and $d_{y}$ in Fig. 2(b) and hence the coupling constants $\kappa_{1}$ and $\kappa_{2}$], which can be served as a polarization selector for practical applications Pol_Sel1 ; Pol_Sel2 .

The threshold of the power density of the vector plasmonic dromion given above can be estimated by using Poynting vector. Based on the above system parameters, the average power of the vector plasmonic dromion is estimated as

[TABLE]

We see that due to the resonant character of the PIT effect in the system, extremely low generation power is required for generating the vector plasmonic dromion.

IV Discussion and summary

It should be mentioned that in writing the dark-state condition (7), the damping coefficient $\gamma_{3}$ of the dark oscillator in the meta-atoms is assumed to be small. However, due to the Ohmic loss inherent in the metal that construct the metamaterial, by our numerical calculation the numerical value of $\gamma_{3}$ is about $0.18\,{\rm GHz}$, which, though smaller than damping coefficients of the CWs ($\gamma_{1}=\gamma_{2}=2.1\,{\rm GHz}$), is still large and has inevitably detrimental impact on the PIT quality. In order to improve the performance of the PIT, one can suppress $\gamma_{3}$ by introducing a gain element into the SRR of the meta-atoms. One possible way is the use of tunneling diodes that have negative resistance and hence may provide gain to the PIT-based metamaterial Jiang ; Ye . Such method has been recognized to be useful for suppressing and even cancelling $\gamma_{3}$, particularly in microwave and THz regimes.

In conclusion, in this article we have considered a plasmonic metamaterial interacting with an EM field with two polarization components. We have proved that such metamaterial can be taken as a classical analogue of an atomic gas with a double-$\Lambda$-type four-level configuration coupled with four laser fields, displays an PIT effect and an equivalent process of atomic FWM. We have shown that, when the nonlinear varactors are mounted onto the gaps of the SRRs, the metamaterial system can acquire giant second- and third-order Kerr nonlinearities via the PIT and the longwave-shortwave interaction. We have also shown that the system supports high-dimensional vector plasmonic dromions, which have very low generation power and are robust during propagation. Our work not only contributes a plasmonic analogue of atomic EIT and FWM but also provides a way for generating novel plasmonic polaritons, and hence opens a new avenue on the exploration of PIT effect in metamaterials.

Acknowledgments

The authors thank Shuang Zhang for helpful discussions. This work was supported by the NSF-China under Grant No. 11474099.

Appendix A Explicit Expressions of the atomic Bloch equation

Explicit expressions of the Bloch equation for the density-matrix elements $\sigma_{jl}$ of the four-level double-$\Lambda$-type atoms are given by

[TABLE]

for diagonal elements, and

[TABLE]

for non-diagonal elements, where $\Gamma_{j}=\sum_{l<j}\Gamma_{jl}$, with $\Gamma_{jl}$ the spontaneous emission decay rate from the state $|l\rangle$ to the state $|j\rangle$; $d_{jl}=\Delta_{j}-\Delta_{l}+i\gamma_{jl}$, $\gamma_{jl}=(\Gamma_{j}+\Gamma_{l})/2+\gamma_{jl}^{\rm dep}$, with $\gamma_{jl}^{\rm dep}$ the dephasing rate between the state $|j\rangle$ and the state $|l\rangle$.

Appendix B Determination of the system parameters in the Lorentz equation

The coefficients in the coupled Lorentz equation Eq. (6) is determined by fitting the numerical result in Fig. 2(c) and Fig. 2(d) obtained by using the finite difference time domain software package (CST Microwave Studio) stated in the main text, and the analytical result of Eq. (6). Assuming the solution of Eq. (6) has the form $q_{l}(z,t)=q_{l0}\exp(ik_{p}z-i\omega_{p}t)+{\rm c.c.}\,(l=1,2,3)$ and $E_{j}(z,t)={\cal E}_{j0}\exp(ik_{p}z-i\omega_{p}t)+{\rm c.c.}\,(j=x,y)$, we have

[TABLE]

with $D_{j}=\omega_{j}^{2}-\omega_{p}^{2}-i\gamma_{j}\omega_{p}$.

The red solid lines in Fig. 2(c) and Fig. 2(d) (where ${\cal E}_{x0}$ and ${\cal E}_{y0}$ have been taken to be real) show the analytical result based on Eq. (32a) in the cases ${\cal E}_{y0}=-{\cal E}_{x0}$ (the excitation condition of the PIT-mode) and ${\cal E}_{y0}={\cal E}_{x0}$ (the excitation condition of the non-PIT-mode), respectively. A better fitting yields $\omega_{1}=\omega_{2}=2\pi\times 13.13\,{\rm GHz}$, $\omega_{3}=2\pi\times 13.08\,{\rm GHz}$ $\gamma_{1}=\gamma_{2}=2.1\,{\rm GHz}$, $\gamma_{3}=0.15\,{\rm GHz}$, $\kappa_{2}=-\kappa_{1}=50\,{\rm GHz}^{2}$ for $d_{x}=d_{y}=4.0\,{\rm mm}$ and $\kappa_{2}=-\kappa_{1}=250\,{\rm GHz}^{2}$ for $d_{x}=d_{y}=3.4\,{\rm mm}$, and $g_{1}=g_{2}=1.79\times 10^{11}\,{\rm C/kg}$.

Appendix C Equations for $q_{\mu j}$ and ${\cal E}_{\mu j}$

The equations of motion for $q_{\mu j}$ read

[TABLE]

with $d_{fj}=\Delta_{j}+i\gamma_{j}/2$ and $d_{sj}=2\Delta_{j}+i\gamma_{j}/2$, and equations for ${\cal E}_{\mu j}$ are given by

[TABLE]

Appendix D Solutions of $q_{\mu j}^{(m)}$ at each order approximations

The first-order [i.e. $O(\epsilon)$] solution has only the fundamental wave, which reads

[TABLE]

with $D_{j}=-2\omega_{j}(\omega+d_{fj})\,(j=1,2,3)$.

At the second order [i.e. $O(\epsilon^{2})$], solution for the fundamental wave is given by

[TABLE]

The expression of $q_{f3}^{(2)}$ is long and omitted for saving space. Solution for the longwave (rectification) reads

[TABLE]

where $n_{m}^{+}(0)=n_{m}^{+}|_{\omega_{p}=0,\omega=0}$, with $n_{m}^{+}(\omega;\omega_{p})=c[k_{p}(\omega_{p})+K_{m}^{+}(\omega;\omega_{p})]/(\omega_{p}+\omega)$ the linear refractive index of the metamaterial. Solutions of the second-harmonic wave is

[TABLE]

with $H_{j}=-4\omega_{j}(2\omega+d_{sj}+3\omega_{j}/4)$, where $\delta_{ji}$ is the Kronecker symbol. Expressions for $q_{s1}^{(2)}$ and $q_{s2}^{(2)}$ are omitted here for saving space.

At the third order [i.e. $O(\epsilon^{3})$], solution of the fundamental wave reads

[TABLE]

The expression of the coefficient $R_{0}$ in Eq. (18) is given by

[TABLE]

Appendix E Nonlinear coefficients in Eq. (17)

The nonlinear property of the SRRs has been theoretically analyzed and experimentally measured in Ref. Wang . The value of $q$ in Ref. Wang represents the renormalized voltage, which has the unit of volt (V), while the value of $q_{3}$ in our work represents the amplitudes of the dark modes, which has the unit of centimeters (cm). To make a comparison we switch the unit of Eq. (6), reading

[TABLE]

where $q_{j}=Q_{0}u_{j}\,(j=1,\,2,\,3)$. $u_{j}$ has the unit of V and $Q_{0}$ has the unit of $\rm cm\cdot V^{-1}$. Thus, nonlinear coefficients $\alpha$ and $\beta$ in our model (6) can be calculated through the dimension transformation $\alpha=Q_{0}^{-1}\alpha_{0}$ and $\beta=Q_{0}^{-2}\beta_{0}$, where $\alpha_{0}=-M\omega_{3}^{2}/(2V_{p})$ and $\beta_{0}=M(2M-1)\omega_{3}^{2}/(6V_{p}^{2})$ can be found in Ref. Wang , readily given by $\alpha_{0}=-2.3503\times 10^{3}\,{\rm V^{-1}GHz^{2}}$ and $\beta_{0}=4.8211\times 10^{2}\,{\rm V^{-2}GHz^{2}}$ with $\omega_{3}=2\pi\times 13.08\,{\rm GHz}$ given in context. Taking a typical value $Q_{0}=1.0\times 10^{-9}\,{\rm V^{-1}cm}$ bai ; Bai2 for dimension transformation, we finally obtain $\alpha=-2.3503\times 10^{12}\,{\rm cm^{-1}GHz^{2}}$ and $\beta=4.8211\times 10^{20}\,{\rm cm^{-2}GHz^{2}}$.