# Propagation and excitation properties of nonlinear surface plasmon
polaritons in a rectangular barrier

Xiangchun Tian, Yundong Zhang, Yu Duan, Yong Zhou  School of Physics and
Electronics, Shandong Normal University, Jinan 250014, China  Chaohua Tan
<EMAIL_ADDRESS>School of Physics and Electronics, Shandong Normal
University, Jinan 250014, China

###### Abstract

We propose a scheme to study the nonlinear propagation properties of nonlinear
surface plasmon polaritons (SPPs) in a three level $\Lambda$ type
electromagnetically induced transparency (EIT) system with modulation of a
rectangular barrier. Based on the multi scale method, the nonlinear
Schrödinger equation (NLSE) describing nonlinear propagation of SPPs is
derived, and the rectangular barrier affecting propagation of nonlinear SPPs
is provided by an off-resonance Stark field. For the single nonlinear SPPs
incident case, by adjusting the height and half width of the barrier, we can
realize transmission, trapping and reflection of the nonlinear SPPs. For two
nonlinear SPPs symmetrical incident case, we find that a periodic intensity
distribution in transverse direction mode can be excited in the rectangular
barrier, and we study the relationship between propagation properties of such
excited modes in the barrier with nonlinearity, half width of the barrier and
phase difference of the initial nonlinear SPPs. In addition, we design an
optical switch of nonlinear SPPs based on the above results. The results
obtained here not only provide a theoretical basis for the study of the
interaction between nonlinear SPPs and external potentials, but also have
broad application prospects in the field of optical information at micro/nano
scale.

##  I Introduction

Surface plasmon polaritons (SPPs) is a kind of electromagnetic eigenmode which
exist in the surface of metal and dielectric material [1]. It can confine the
electromagnetic fields to the spatial extent which scale is far smaller than
the wavelength, therefore, SPPs can overcome the diffraction limitation, and
produce strong local field enhancement effect. Recently, SPPs have been widely
researched and applied in the fields of light-matter nonlinear interaction
enhancement [2, 3, 4, 5], optical communication at micro/nano scale [6, 7, 8,
9], nano-focusing [10, 11, 12, 13], super resolution display [14, 15] and so
on.

With the further understanding of SPPs, how to manipulate SPPs effectively at
micro/nano scale has aroused wide attention. Meanwhile, several kinds of
methods have been proposed in order to better modulate the excitation and
propagation properties of SPPs:

Firstly, by changing or designing different structures of metals, one can
realize manipulation of SPPs [16, 17, 18, 19, 20, 21, 22, 24, 23], e.g.,
carving gratings-based structure on the metal surface is an effective way to
realize unidirectional control of SPPs [17, 16]; using the metal film covered
with a thin patterned dielectric layer can be used to restrict or redirect
SPPs [18, 19]; adopting metasurface can also manipulate the SPPs effectively
[23, 24].

Secondly, some new materials (such as graphene, sodium) have also been
proposed to make dynamic manipulation of SPPs possible [25, 26, 28, 27, 29,
30]. For example, using the negative-index metamaterial (NIMM) instead of
metal to excite lower loss surface polaritons (SPs) modes at the dielectric-
NIMM interface [30] can realize coherent control of SPs for a long distance.

At last, introducing external potential fields to change the electromagnetic
environment in the system can also be used to modulate properties of SPPs [31,
32, 33, 34, 35, 36, 37, 38, 39]. e.g., by using a beam of control light with a
specific polarization state to shine on the surface of aluminium, then the
refractive index of aluminum will change due to its absorption of control
light, thus, the SPPs can be manipulated by the incident light field [31]; the
phase information and trajectory of SPPs can be modulated by introducing an
external gradient magnetic field [32, 33, 34]. Besides, there are some other
kinds of external potential fields can be used to realize the modulation of
SPPs, such as electric field [35, 36, 37] and thermal field [38, 39].

From the above, one can see that the interaction between SPPs and external
potential fields has been widely studied experimentally, and some interesting
properties of SPPs has been found under the modulation of external potential
fields, e.g., the reflection, transmission and trapping of SPPs in the
external potentials or interfaces [40, 41]. However, theoretical exploration
of interaction mechanism between SPPs and external potentials is relative
deficiency, especially for some fundamental cases, such as SPPs with potential
wells and SPPs with potential barriers.

Thus, in this work, we propose a scheme to study interaction between SPPs and
a rectangular barrier, including nonlinear propagation and excitation of SPPs
in the barrier, in a NIMM-atomic gas interface waveguide system. We first
obtain the low loss stable propagation of nonlinear SPPs via
electromagnetically induced transparency (EIT). By introducing a Stark field
into the system, a rectangular barrier is produced, we investigate the
reflection, tranmission and trapping dynamics of single nonlinear SPPs in the
rectangular barrier. Meanwhile, we also obtain the stable propagation modes
existed in barrier by inputting two symmetrically incident nonlinear SPPs.
Then, we show the influence of several major factors on the propagation modes,
such as nonlinearity, width of the barrier and phase difference of initial
nonlinear SPPs. At last, we introduce its potential application in optical
switch.

The rest of the article is arranged as follows: In Sec. II, we introduce our
theoretical model and give the Maxwell-Bloch (MB) equations for the system. In
Sec. III, the nonlinear Schrödinger equation (NLSE) describing the envelope of
the probe field is given. In Sec. IV, we investigate the propagation
properties of SPPs under the manipulation of the rectangular barrier. In final
section, we summarize the main results obtained in this work.

##  II Theoretical Model

Figure 1: (Color online) (a) The schematic of SPPs is excited in the interface
between the NIMM and the hot atomic gas and propagates in the $z$ direction.
In this waveguide system, the control field (blue arrow) is incident along the
vertical direction. The purple arrow stands for the Stark field which
distributes as a rectangular function. (b) The distribution function of the
Stark field along the y direction. (c) The energy-level diagram of
$\Lambda$-type atomic system (the details are given in the text).

The waveguide system we proposed here consists of a layer of NIMM in the lower
half plane $x<0$ (with relative permittivity ${\varepsilon_{2}}$ and relative
permeability $\mu_{2}$) and a hot atomic gas filled above the NIMM, as shown
in Fig. 1. The hot atomic gas we using here has a $\Lambda$-type energy level
excited configuration, two ground states $\left|1\right\rangle$ and
$\left|2\right\rangle$ couple to the excited state $\left|3\right\rangle$ with
the weak probe field $\mathbf{E}_{p}$ with center angle frequency $\omega_{p}$
(which propagates along $z$ direction and is chosen from guide modes of SPPs
at the interface of NIMM) and the strong continuous control field
$\mathbf{E}_{c}$ with center angle frequency $\omega_{c}$ (which is incident
vertically to the waveguide interface), respectively. ${\Delta_{3}}$
(${\Delta_{2}}$) is the one (two) photon detuning. The atoms occupy on the
state $\left|3\right\rangle$ can spontaneously radiate to the ground states
with the rate $\Gamma_{j3}$ ($j=1,2$). The two ground states are chosen from
the hyperfine structure, thus, the transition between the two ground states is
forbidden.

The NIMM interface waveguide system can support both transverse electric (TE)
and transverse magnetic (TM) modes, the TE modes decay rapidly in the system,
however, the TM modes can propagate without loss under specific conditions
[30]. Thus, we choose the lowest order TM mode as the probe field of our
system, which can be expressed as
$\mathbf{E}_{p}(\mathbf{r},t)={{\mathcal{E}}_{p}}(y,z,t)\mathbf{u}_{p}(x){e^{i\left[{{k(\omega_{p})}z-{\omega_{p}}t}\right]}}+c.c.$,
where ${{\mathcal{E}}_{p}}$ is the envelope of the probe field,
$\mathbf{u}_{p}(x)$ is the mode function in the $x$ direction (which is given
in Appendix A),
${k}(\omega_{p})=\left({{\omega_{p}}/c}\right){\left[{{\varepsilon_{1}}{\varepsilon_{2}}\left({{\varepsilon_{1}}{\mu_{2}}-{\varepsilon_{2}}{\mu_{1}}}\right)/\left({\varepsilon_{1}^{2}-\varepsilon_{2}^{2}}\right)}\right]^{1/2}}$
is the propagate constant of the probe field without modulation of the atomic
gas. The control field reads
${\mathbf{E}_{c}}\left({\mathbf{r},t}\right)={{\mathcal{E}}_{c}}{\mathbf{e}_{c}}{e^{i\left({{k_{c}}x-{\omega_{c}}t}\right)}}+c.c.$,
in which ${{\mathcal{E}}_{c}}$ is strength, $\mathbf{e}_{c}$ is the unit
polarization vector and $k_{c}$ is the propagate constant of the control
field.

Meanwhile, we apply another off-resonance laser field (Stark field)
propagating along the interface of the NIMM in $z$ direction to generate a
rectangular barrier in $y$ direction, with center frequency ${\omega_{s}}$,
and reading

${\mathbf{E}_{s}}(y,t)={\mathbf{e}_{s}}\sqrt{2}{E_{s}}(y)\cos({\omega_{s}}t),$
(1)

where ${E_{s}}(y)$ is the amplitude of the Stark field and $\mathbf{e}_{s}$ is
the unit polarization vector of the Stark field. Such a Stark field can
produce a small energy level shift to state $\left|j\right\rangle$ in $y$
direction, as
$\Delta{E_{j}}=-{\alpha_{j}}{\left\langle{\mathbf{E}_{s}^{2}}\right\rangle_{t}}/2=-{\alpha_{j}}{\left|{{E_{s}}\left(y\right)}\right|^{2}}/2$,
in which $\alpha_{j}$ is the scalar polarizability of the state
$\left|j\right\rangle$,
${\left\langle{\mathbf{E}_{s}^{2}}\right\rangle_{t}}=\int_{0}^{T}{\mathbf{E}_{s}^{2}}dt/T$
is the time average in an oscillating cycle. After taking the energy level
shift caused by the Stark field into account, the optical detunings should be
corrected to
${\Delta_{j}}^{\prime}={\Delta_{j}}+{\alpha_{j1}}{\left|{{E_{s}}\left(y\right)}\right|^{2}}/2$,
with ${\alpha_{jl}}=\left({{\alpha_{j}}-{\alpha_{l}}}\right)/\hbar$.

In the interaction picture, under electric-dipole and rotating-wave
approximation, the Hamiltonian of the system is given by

${\hat{H}_{{\rm{int}}}}=-\hbar\left[{\sum\limits_{j=1}^{3}{{\Delta_{j}}^{\prime}\left|j\right\rangle\langle
j|}+{\Omega_{c}}\left|3\right\rangle\langle
2|+{\zeta}\left(x\right)e^{i\theta_{p}}{\Omega_{p}}\left|3\right\rangle\langle
1|+\rm H.\rm c.}\right],$ (2)

where
${\Omega_{p}}=\left|{{{\mathbf{p}}_{13}}}\right|{\mathcal{E}_{p}}/\hbar$,
${\Omega_{c}}=\left|{{{\mathbf{p}}_{23}}}\right|{\mathcal{E}_{c}}/\hbar$,
$\zeta\left(x\right)={\mathbf{e}_{13}}\cdot{\mathbf{u}_{p}}\left(x\right)$,
and $\mathbf{p}_{ij}$ is the electric-dipole matrix element associated with
the state $\left|i\right\rangle$ to state $\left|j\right\rangle$,
$\theta_{p}=[k(\omega_{p})+k_{1}-k_{3}]z$ is the phase mismatch caused by the
dispersion of SPPs, $k_{j}$ ($j$=1, 2, and 3) are the wave number of state
$\left|j\right\rangle$. The motion of atoms is governed by the optical Bloch
equation

$\frac{{\partial\sigma}}{{\partial
t}}=-\frac{i}{\hbar}\left[{{{\hat{H}}_{{\mathop{\rm
int}}}},\sigma}\right]-\Gamma\sigma,$ (3)

where $\sigma$ is a $3\times 3$ density matrix, $\Gamma$ is a $3\times 3$
relaxation matrix describing the spontaneous emission and dephasing effect of
the system. The explicit expressions of Eq. (3) are given in Appendix B.

The dynamics evolution of the probe field in this system is governed by the
Maxwell equations. Under the slowly varying envelope approximation, the
Maxwell equation can be reduced to

$i\left({\frac{\partial}{{\partial z}}+\frac{1}{c}\frac{1}{{{n_{\rm
eff}}}}\frac{\partial}{{\partial
t}}}\right){\Omega_{p}}{e^{i\theta_{p}}}+\frac{c}{{2{\omega_{p}}}}\frac{1}{{{n_{\rm
eff}}}}\frac{{{\partial^{2}}}}{{\partial{y^{2}}}}{\Omega_{p}}{e^{i\theta_{p}}}+{\kappa_{13}}\int_{-\infty}^{\infty}{d(kv)f\left(kv\right)\left\langle{{\sigma_{31}}}\right\rangle}=0,$
(4)

where ${n_{\rm eff}}=c{k(\omega_{p})}/{\omega_{p}}$ is the effective
refractive index,
${\kappa_{13}}={{{\mathcal{N}_{\alpha}}{{\left|{{\mathbf{p}_{13}}}\right|}^{2}}\omega_{p}^{2}}\mathord{\left/{\vphantom{{{\mathcal{N}_{\alpha}}{{\left|{{\mathbf{p}_{13}}}\right|}^{2}}\omega_{p}^{2}}{[2{\varepsilon_{0}}{c^{2}}\hbar\tilde{k}({\omega_{p}})]}}}\right.\kern-1.2pt}{[2{\varepsilon_{0}}{c^{2}}\hbar\tilde{k}({\omega_{p}})]}}$
is the coupling coefficient describing the interaction between the weak probe
field and the atoms, ${\mathcal{N}_{\alpha}}$ represents the number density of
the atoms and $\tilde{k}={\mathop{\rm Re}\nolimits}\left(k\right)$. The
expectation operator $\left\langle{...}\right\rangle$ is defined as
$\left\langle{\psi\left(x\right)}\right\rangle\equiv\int_{-\infty}^{\infty}{dx{\zeta^{*}}\left(x\right)\psi\left(x\right)}/\int_{-\infty}^{\infty}{dx{{\left|{\zeta\left(x\right)}\right|}^{2}}}$,
i.e., average over mode distribution function in $x$ direction. In addition,
the hot atomic gas is above the NIMM, thus, the Doppler effect must be taken
into account, which could cause an inhomogeneous broadening of the energy
level. In Eq. (4), we have averaged the inhomogeneous broadening due to
Doppler effect over the Maxwell velocity distribution $f(kv)$ [42], and
$f(kv)$ reads

$f\left(kv\right)=\frac{{\sqrt{\ln 2}}}{{\sqrt{\pi}{W_{D}}}}{e^{-\ln
2{{\left({\frac{{kv}}{{{W_{D}}}}}\right)}^{2}}}},$ (5)

where ${W_{D}}=k\sqrt{2\ln 2{k_{B}}T/m}$ is the full width at the half-maximum
(FWHM) of the inhomogeneous broadening, $k_{B}$ is the Boltzmann constant, $T$
and $m$ are the temperature and mass of the atom, respectively.

Equations (3) and (4) are known as MB equations which describe the interaction
dynamics of the system.

##  III Derivation of the Nonlinear Schrödinger Equation

In EIT system, the probe field is weaker than the control field, we can adopt
the standard multiple scales method [43] to solve the MB equations. The
asymptotic expansion are given as
$\sigma_{jl}-\sigma_{jl}^{(0)}=\epsilon\sigma_{jl}^{(1)}+\epsilon^{2}\sigma_{jl}^{(2)}+\epsilon^{3}\sigma_{jl}^{(3)}+...$,
$\Omega_{p}=\epsilon\Omega_{p}^{(1)}+\epsilon^{2}\Omega_{p}^{(2)}+...$, where,
$\epsilon$ is a dimensionless small parameter characterizing the typical
amplitude ratio of the probe and the control field, and the physical
quantities in the left of the expressions can be treated as the functions of
the multiple scales variables, with $z_{j}=\epsilon^{j}z$ $(j=0,2)$,
$y_{j}=\epsilon^{j}y$ $(j=1)$ and $t_{j}=\epsilon^{j}t$ $(j=0,2)$. In our
analysis, we assume that the amplitude of Stark field can be written as
$E_{s}=\epsilon E_{s}^{\left(1\right)}$. Then $d_{ij}$ can be expressed as
$d_{ij}=d_{ij}^{\left(0\right)}+{\epsilon^{2}}d_{ij}^{\left(2\right)}$, with
$d_{ij}^{\left(2\right)}={\alpha_{ij}}{\left|{E_{s}^{(1)}}\right|^{2}}/2$.
Substituting the above expansion into the MB equation, we can obtain a series
linear but inhomogeneous equations order by order.

The zeroth order solution corresponds to the steady state of the system, i.e.,
the probe field is absent. The specific expression of the zeroth order
solution reads $\sigma^{(0)}_{11}=1$, and the zeroth order of all the other
elements of the density matrix are equal to zero, i.e., $\sigma_{ij}^{(0)}=0$,
which means that, when the system is prepared, due to the pumping effect of
the strong control field, all the atoms occupy on the state $|1\rangle$, and
there are no coherence between all the states.

In the first order solution, we can obtain the linear propagation property of
the probe field. When the probe field is incident, the population and
coherence in the system should be modified. However, due to $\Omega_{p}$ is
taken as a small quantity, thus, the population on state $|1\rangle$ is not
changed. We assume that the solution of first order of the probe field has the
following form $\Omega_{p}^{(1)}=Fe^{i\theta}$, where
$\theta=K(\omega)z_{0}-\omega t_{0}$, $F$ is the slowly varying envelop
function of multiscale variables, $K$ is the linear dispersion relation of the
probe field, the imaginary part of $K$ means linear absorption of the system
and the real part of $K$ means linear dispersion of the system. The expression
of the linear dispersion relation $K(\omega)$ reads

${K}\left(\omega\right)=\frac{\omega}{c}\frac{1}{{{n_{\rm
eff}}}}+{\kappa_{13}}\int_{-\infty}^{\infty}{d(kv)f\left(kv\right)\left\langle{\frac{{\left({\omega+d_{21}^{\left(0\right)}}\right)}}{{{D}}}\zeta\left(x\right)}\right\rangle},$
(6)

where
${{D}={{\left|{{\Omega_{c}}}\right|}^{2}}-\left({\omega+d_{21}^{\left(0\right)}}\right)\left({\omega+d_{31}^{\left(0\right)}}\right)}$,
$\omega$ is the frequency shift to the center frequency of the probe field
$\omega_{p}$, and $K(\omega)$ actually is also the wavenumber shift to
$k(\omega_{p})$.

The absorption spectrum Im$(K)$ of the probe field will open a transparency
window at the center frequency for a large control field, i.e., EIT effect.
Although the EIT effect can effectively inhibit the linear absorption of the
medium to the probe field and slow down the group velocity (which is defined
as ${V_{g}}={\left[{{{\partial{\rm
Re}({K}\left(\omega\right))}\mathord{\left/{\vphantom{{\partial{K}\left(\omega\right)}{\partial\omega}}}\right.\kern-1.2pt}{\partial\omega}}}\right]^{-1}}$)
of the probe field, the diffraction and nonlinear effects will lead to serious
deformation of the envelope of the probe field in the propagation process.
Therefore, we solve the MB equations up to third order to investigate the
nonlinear propagation and interaction of the probe field in the rectangular
potential produced by the Stark field.

In the third-order approximation equations, we obtain the solvable condition
of $F$, which reads

$i\left({\frac{\partial}{{\partial{z_{2}}}}+\frac{1}{{{V_{g}}}}\frac{\partial}{{\partial{t_{2}}}}}\right){F}+\frac{c}{{2{\omega_{p}}}}\frac{1}{{{n_{\rm
eff}}}}\frac{{{\partial^{2}}}}{{\partial
y_{1}^{2}}}{F}+{W_{1}}{\left|{{F}}\right|^{2}}{F}{e^{-2{{\tilde{\alpha}}}{z_{2}}}}+{W_{2}}{\left|{E_{s}^{\left(1\right)}\left(y\right)}\right|^{2}}{F}=0,$
(7)

with ${{\tilde{\alpha}}}={\epsilon^{-2}}\alpha={\epsilon^{-2}}{\mathop{\rm
Im}\nolimits}\left[{{K}\left(\omega\right)}{+k(\omega_{p})}\right]$, $W_{1}$
is the nonlinear coefficient describing the self-phase modulation of the probe
field and $W_{2}$ is the external potential modulation coefficient, specific
expressions are given in Appendix C. By defining $U=\epsilon Fe^{-\alpha z}$
and returning the variables in Eq. (7) to their original scale, one can obtain
the equation describing evolution of envelope $F$ of the probe field, which
reads

$i\left({\frac{\partial}{{\partial
z}}+\frac{1}{{{V_{g}}}}\frac{\partial}{{\partial
t}}}\right)U+\frac{c}{{2{\omega_{p}}}}\frac{1}{{{n_{\rm
eff}}}}\frac{{{\partial^{2}}}}{{\partial{y^{2}}}}U+{W_{1}}{\left|U\right|^{2}}U+{W_{2}}{\left|{{E_{s}}\left(y\right)}\right|^{2}}U=-i\alpha
U.$ (8)

Generally, the coefficients of Eq. (8) are complex, and $\alpha\neq 0$, which
means the linear absorption and nonlinear dissipation in our system will
hugely affect the propagation of the probe field. In order to study excitation
and propagation of nonlinear SPPs in the rectangular barrier, the probe field
need to propagate for a long distance stably. Thus, firstly, the linear and
nonlinear absorption of the system must be suppressed, secondly, the
diffraction effect of the probe field must be balanced by the nonlinearity
provided by the system.

However, when the system works under EIT condition, a realistic set of
physical parameters, which will be shown in the following, can be found to
make the imaginary part of the coefficients in Eq. (8) to be negligible small
compare to its corresponding real parts, together with $\alpha\approx 0$, and
a giant Kerr nonlinearity.

Based on the above condition, we define a set of dimensionless variables as
$\xi=\sqrt{{n_{\rm eff}}}y/{R_{y}}$, $s=z/{L_{D}}$, $\tau=t/{\tau_{0}}$,
$u=U/{u_{0}}$, $g={V_{g}}{\tau_{0}}/{L_{D}}$, ${w_{1}}={L_{D}}/{L_{N}}$,
${w_{2}}={L_{D}}{\rm Re}({W_{2}})E_{0}^{2}$,
$V\left(\xi\right)={\left({{E_{s}}/{E_{0}}}\right)^{2}}$, with
${L_{D}}={\omega_{p}}R_{y}^{2}/c$ being the typical diffraction length,
$R_{y}$ being the transverse width of the SPPs in the $y$ direction, and
$L_{N}=1/\left[{{\mathop{\rm
Re}\nolimits}\left({{W_{1}}}\right)u_{0}^{2}}\right]$ being the typical
nonlinear length, $E_{0}$ being a real parameter representing typical strength
of the Stark field. Then, Eq. (7) can be simplified as the form of the
dimensionless NLSE

$i\left({\frac{\partial}{{\partial
s}}+\frac{1}{g}\frac{\partial}{{\partial\tau}}}\right)u+\frac{1}{2}\frac{{{\partial^{2}}u}}{{\partial{\xi^{2}}}}+{w_{1}}{\left|u\right|^{2}}u+{w_{2}}V\left(\xi\right)u=0.$
(9)

To solve the above equation, we assume the solution as
$u\left({\tau,s,\xi}\right)=\lambda\left({\tau,s}\right)v\left({s,\xi}\right)$
based on separated variable method. For $\lambda\left({\tau,s}\right)$, we
take it as a Gaussian profile
$\lambda\left({\tau,s}\right)={\left[{1/({\rho_{0}}\sqrt{\pi})}\right]^{1/2}}\exp\left[{-\rho^{2}/\left({2\rho_{0}^{2}}\right)}\right]$,
with $\rho=\tau-s/g$, and ${\rho_{0}}$ is a free real parameter [44, 45].
Then, Eq. (9) can be further simplified as

$\left({i\frac{\partial}{{\partial
s}}+\frac{1}{2}\frac{{{\partial^{2}}}}{{\partial{\xi^{2}}}}}\right)v+{w_{1}}{\left|v\right|^{2}}v+{w_{2}}V\left(\xi\right)v=0.$
(10)

It has been proved that, Eq. (10) can be used to describe the propagation and
interaction of nonlinear optical pluse with the modulation of the external
potential.

##  IV Propagation properties of nonlinear SPPs interacting with the
rectangular barrier

Next, we choose a realistic system to carry out our study. The NIMM is
selected as a sliver based material, the relative permittivity and
permeability of the NIMM can be described by the Drude model in optical
region, i.e.,
${\varepsilon_{2}}\left(\omega\right)={\varepsilon_{\infty}}-\omega_{e}^{2}/\left[{\omega\left({\omega+i{\gamma_{e}}}\right)}\right]$,
${\mu_{2}}\left(\omega\right)={\mu_{\infty}}-\omega_{m}^{2}/\left[{\omega\left({\omega+i{\gamma_{m}}}\right)}\right]$,
where ${\omega_{e}}$ ($\omega_{m}$) is the electric (magnetic) plasma
frequency, ${\gamma_{e}}$ ($\gamma_{m}$) is the electric (magnetic) decay
rate, ${\varepsilon_{\infty}}$ and ${\mu_{\infty}}$ are the background
constant. The related parameters are ${\varepsilon_{1}}=1,{\mu_{1}}=1$,
${\varepsilon_{\infty}}={\mu_{\infty}}=1$, $\omega_{e}/2\pi=2.18\times
10^{15}~{}\rm{s^{-1}}$, $\gamma_{e}/2\pi=4.35\times 10^{12}~{}\rm{s^{-1}}$,
$\omega_{m}=\omega_{e}/5$, $\gamma_{m}=\gamma_{e}/1000$ (as for Ag) [30]. The
$\Lambda$ type energy level system is chosen from D2-line of ${}^{87}\rm{Rb}$
atomic gas at room temperature, with the levels
$\left|1\right\rangle=\left|{5{}^{2}{S_{1/2}},~{}F=1,~{}{m_{F}}=0}\right\rangle$,
$\left|2\right\rangle=\left|{5{}^{2}{S_{1/2}},~{}F=2,~{}{m_{F}}=0}\right\rangle$
and
$\left|3\right\rangle=\left|{5{}^{2}{P_{3/2}},~{}F=2,~{}{m_{F}}=1}\right\rangle$
[46]. Thus, ${\lambda_{c}}\simeq{\lambda_{p}}=780~{}\rm{nm}$,
$\left|{{\mathbf{p}_{13}}}\right|\simeq\left|{{\mathbf{p}_{23}}}\right|\simeq
1.27\times{10^{-27}}~{}\rm{C\cdot cm}$, ${\Gamma_{3}}\simeq 2\pi\times
6~{}\rm{MHz}$, $\Gamma_{13}=\Gamma_{23}=0.5\Gamma_{3}$, and the Doppler width
is $2{W_{D}}=3.21~{}\rm{GHz}$. We assume that the atomic density is
${\mathcal{N}_{\alpha}}=7.18\times{10^{13}}~{}\rm{c{m^{-3}}}$ and obtain
$\kappa_{13}=5\times{10^{11}}~{}\rm{c{m^{-1}}{s^{-1}}}$.

In this case, the wavenumber of the probe field without modulation of the
atomic gas $k(\omega_{p})=(8.05\times 10^{4}+0.053i)$ cm-1. We choose the
control field ${\Omega_{c}}=2\pi\times 0.2~{}\rm{GHz}$, the detunings
${\Delta_{2}}=-2\pi\times 40.25~{}\rm{MHz}$, ${\Delta_{3}}=-2\pi\times
31.11~{}\rm{MHz}$. Then, the coefficients of Eq. (8) are
${W_{1}}=\left({7.53-0.89i}\right)\times{10^{-17}}~{}\rm{c{m^{-1}}{s^{2}}}$,
${W_{2}}=\left({1.26-0.02i}\right)\times{10^{-9}}~{}\rm{c{m^{-1}}{V^{2}}}$ and
the group velocity is ${V_{g}}=7.69\times{10^{-5}}c$. The total linear
absorption of the probe field is
$\alpha={\rm{Im}}[K(\omega)+k(\omega_{p})]=0.37$ cm-1, thus, the typical
absorption length $L_{A}=1/\alpha=2.69$ cm. In addition, by taking
$\tau_{0}=0.25\times{10^{-7}}~{}\rm{s}$ and $R_{y}=8.46~{}\rm{\mu m}$, one can
obtain the typical diffraction length and the nonlinear length
${L_{D}}\simeq{L_{N}}=0.058~{}\rm{cm}$, and $w_{1}=w_{2}=1$. It is easy to see
that the imaginary parts of these coefficients in Eq. (8) are negligible small
compared to their corresponding real parts, thus we can ignore the effect of
these imaginary parts, together with linear absorption of the system. The
following discussion of SPPs propagation and excitation in a rectangular
barrier is based on the above set of parameters.

As we mentioned above, the influence of the Stark field is reflected as the
external potential term $V\left(\xi\right)$ in Eq. (10), which can produce a
spatial modulation to the propagation of nonlinear SPPs. In the following, we
will investigate the propagation properties of nonlinear SPPs interacting with
the external potential based on the Eq. (10). A intensity uniform Stark field
with a transversely finite width can produce a $V(\xi)$ with the form as a
rectangular function (a finite high barrier), which reads

$V\left(\xi\right)=\left\\{{\begin{array}[]{*{20}{c}}{{V_{0}}}&{{\rm{for}}}&{\left|\xi\right|\leq
l}\\\ 0&{{\rm{for}}}&{\left|\xi\right|>l}.\end{array}}\right.$ (11)

where $V_{0}$ and $l$ are the height and half width of the potential barrier
$V(\xi)$, respectively.

In the case of $w_{1}>0$ and $V(\xi)=0$, Eq. (10) supports bright soliton
solutions with the form as

$v={\varsigma_{0}}{\mathop{\rm
sech}\nolimits}\left[{{\varsigma_{0}}\left({\xi-{\eta_{0}}s-{\xi_{0}}}\right)}\right]{e^{i\left[{{\eta_{0}}\xi-\left({\eta_{0}^{2}-\varsigma_{0}^{2}}\right)s/2-{\varphi_{0}}}\right]}},$
(12)

where $\varsigma_{0},\eta_{0},{\xi_{0}}$ and $\varphi_{0}$ are the free
parameters [47]. Thus, we take this solution as the initial condition of the
probe field to investigate the propagation and excitation properties of SPPs
interacting with the rectangular barrier.

###  IV.1 The transmission, trapping and reflection of a nonlinear SPPs

Shown in Fig. 2 are simulation results of inputting the nonlinear SPPs with
the form of single bright soliton. In Fig. 2(a), the specific initial
condition of the nonlinear SPPs is $v={\rm sech}(\xi-10)e^{-i0.4\xi}$, we can
see that, when the rectangular barrier is absent, i.e., $V(\xi)=0$, as red
line in inset of Fig. 2(a), the nonlinear SPPs can propagate stably for a long
distance due to the balance between the diffraction effect and the
nonlinearity of the SPPs.

Shown in Fig. 2(b) is the transmission of the nonlinear SPPs through the
rectangular barrier. The specific initial condition of the nonlinear SPPs is
the same as Fig. 2(a). The height and half width of the rectangular barrier
are $V_{0}=0.2$ and $l=3$, respectively, the profile of $V(\xi)$ is as the red
line in the top of Fig. 2(b), and the white dashed lines denote the boundary
of the barrier. The nonlinear SPPs first contacts the right boundary, and
partial energy of the SPPs is reflected out of the system, but we can see that
the energy loss is negligible, the nonlinear SPPs can keep its shape and pass
through the barrier with a phase shift.

However, when enhancing the strength of the rectangular barrier and decrease
the incident angle of the nonlinear SPPs, as shown in Fig. 2(c), $V_{0}=0.4$
and $v={\rm sech}(\xi-10)e^{-i0.25\xi}$. The nonlinear SPPs can also go though
the right boundary of the barrier and its most energy can go into the barrier,
then the SPPs keep going, and on the left boundary, due to decay of amplitude
of the SPPs, it is mainly reflected, and the pulse coming out from the left
boundary part quickly broadens due to the diffraction effect. Inside the
barrier, the SPPs is reflected back and forth on the left and right boundaries
with continuous attenuation of amplitude, thus, we can realize the trapping of
SPPs, but the effective propagation length of SPPs is short due to the energy
loss on the boundary. In Fig. 2(d), the height and half width of the
rectangular barrier are $V_{0}=2$ and $l=1.2$, respectively. For such a high
and narrow barrier, we choose a specific initial condition which is the same
as Fig. 2(b), we can see that the nonlinear SPPs is mainly reflected on the
right boundary. Therefore, by adjusting the height and half width of the
rectangular barrier, we can realize the transmission, trapping and reflection
of the nonlinear SPPs.

Figure 2: The propagation properties of single nonlinear SPPs in the
rectangular barrier. (a) The stable evolution process of a nonlinear SPPs
without $V(\xi)$. (b) The tranmission process of a nonlinear SPPs through the
rectangular barrier. (c) The trapping process of a nonlinear SPPs in the
rectangular barrier. (d) The reflection process of a nonlinear SPPs at the
right boundary of the rectangular barrier.

###  IV.2 Excitation of stable propagation modes of nonlinear SPPs in the
rectangular barrier

Figure 3: (a) The periodic propagation modes which are excited by two
symmetrically incident SPPs in the rectangular barrier. (b) The profile of
intensity $|v|^{2}$ as a function of $\xi$ at $s=30$.

In the previous subsection, we can see that even the SPPs can be trapped by
the barrier, it can not propagate stably in the barrier, i.e., the stable
propagation modes of the barrier can not be excited by inputting a single
nonlinear SPPs on one side of the barrier. However, when inputting two
nonlinear SPPs symmetrically on the two boundaries, as shown in Fig. 3(a), the
initial condition of the nonlinear SPPs are $v={\rm
sech}(\xi-12)e^{-i0.5\xi}+{\rm sech}(\xi+12)e^{i0.5\xi}$, with $V_{0}=5$ and
$l=4$, we can see that, after the two nonlinear SPPs incident on the two
boundaries of the barrier, partial energy reflects out of the system, but the
parts transit into the barrier can produce a periodic intensity distribution
in transverse direction $\xi$, i.e., a ’standing wave’ shape mode can be
excited inside the barrier, and it can propagate stably for a long distance.
The blue solid line shown in Fig. 3(b) is the profile of intensity $|v|^{2}$
as a function of $\xi$ at $s=30$, the red solid line on top of the figure is
the corresponding barrier, and the black dashed lines denote the boundary of
the barrier. Such a ’standing wave’-like SPPs can be used to design phase or
intensity SPPs gratings, which have great application potentials in SPPs based
precision measurement at micro/nano scale.

The physical reason for this stable propagation and periodic distributed mode
is that, inside the barrier, the two SPPs are reflected back and forth,
respectively, thus, they can meet periodically, and lead to quantum
interference, such interference can generate a periodic intensity
distribution.

However, our in-depth study find that, not all the periodic modes excited in
the barrier can propagate stably, it depends on three major factors, which are
nonlinearity, half width of the barrier and phase difference of the initial
nonlinear SPPs, respectively.

Figure 4: (a) and (b) express the detail propagation circumstance of the
propagation modes in the potential barrier for the case with the nonlinearity
modulation or without it, respectively. (c) The variation tendency about the
average intensity of the propagation modes with the growth of propagation
distance $s$ under the circumstance of nonlinearity modulation or without it.

Shown in Fig. 4(a) and (b) are intensity $|v|^{2}$ as functions of $s$ and
$\xi$ with and without nonlinearity, respectively. The initial conditions are
the same as that used in Fig. 3 but with $V_{0}=5$ and $l=3$. We can see that,
when the two nonlinear SPPs enter the barrier and interact around the position
at $s=20$, the periodic modes can both be excited, for the case with
nonlinearity (Fig. 4(a)), intensity of the mode out of the barrier at $s=50$
is much larger than that in the case without nonlinearity (Fig. 4(b)). To be
more intuitive, we define the average intensity of the mode inside the barrier
in $\xi$ direction as

$\bar{|v|^{2}}(s)=\frac{\int_{-l}^{l}|v(s,\xi)|^{2}d\xi}{\int_{-l}^{l}d\xi}.$
(13)

In Fig. 4(c), the blue solid (red dashed) line is $\bar{|v|^{2}}$ as a
function of $s$ for the case with nonlinearity (without nonlinearity). The
first peaks from left denote the position where the modes be excited for the
two cases, but after the excitation process, we can see that the blue solid
line first drops, then oscillates and flattens as increase of $s$ due to the
modulation of the nonlinearity, which means the periodic mode be excited can
propagate stably. However, when the nonlinearity of the system is absence, the
average intensity $\bar{|v|^{2}}$ keeps decaying to a quite small value with
the increase of $s$, shown as red dashed line after $s=20$, which means the
periodic mode even can be excited in the system without nonlinearity, it can
not propagate stably in the barrier, which can be treated as decay modes, and
the parts of the nonlinear SPPs out of the barrier are broadening due to the
diffraction effect.

Figure 5: (a) The changing relations between the average intensity
$\bar{|v|^{2}}$ and the half width of the barrier potential $l$ for the case
with different phase difference. (b) - (i) show the corresponding profile of
$|v|^{2}$ with $\Delta\varphi=0$ (blue solid line) and $\Delta\varphi=\pi$
(red solid line), respectively.

In addition to the effect of nonlinearity, we study the influence of width of
the barrier to the propagation modes in the next. Shown in Fig. 5(a) are the
average intensity $\bar{|v|^{2}}$ at $s=20$ as functions of half width of the
barrier $l$, in which, we have selected two typical phase difference
$\Delta\varphi=0$ (blue solid line, attractive interaction of nonlinear SPPs
inside the barrier) and $\Delta\varphi=\pi$ (red dashed line, repulsive
interaction of nonlinear SPPs inside the barrier), and the initial condition
of the two nonlinear SPPs are the same as that given in Fig. 3(a), the height
of the barrier is $V_{0}=5$ for both cases. From blue solid line in Fig. 5(a),
we can find that the average intensity $\bar{|v|^{2}}$ has quasi periodic
peaks with the increase of $l$, and the value of the peaks continuously
decrease with the increase of $l$, which means with the increase of $l$, quasi
periodic coherence enhancement occurs in the barrier. The peaks mean that
propagation modes can be excited in the barrier with corresponding half
widths, but the dips in the blue solid line mean decay modes are excited,
i.e., even modes can be excited, they can not propagate for a long distance.
In addition, we can obtain the similar results for $\Delta\varphi=\pi$, shown
as red dashed line in Fig. 5(a).

Shown in Fig. 5(b)-(e) are the wave shapes $|v|^{2}$ as functions of $\xi$ at
$s=20$, the corresponding half width $l$ are selected from positions of the
first four peaks of the blue solid line in Fig. 5(a), the purple dash lines
denote the boundaries of the barriers. We can see that there are three
complete peaks inside the barrier in Fig. 5(b), and the number of complete
peaks are 5, 7, 9 in Fig. 5(c), (d) (e), respectively, we define those modes
as odd modes. For $\Delta\varphi=0$, we can find that the number of complete
peaks are odd numbers. In contrary, the number of complete peaks are even
numbers for $\Delta\varphi=\pi$, as shown in Fig. 5(f)-(i). Thus, by adjusting
the half width $l$ and phase difference $\Delta\varphi$, we can modulate
structure of the propagation modes in the barrier.

In addition, it is worth noting that we have selected two typical phase
differences in the simulation of last figure, and peaks for $\Delta\varphi=0$
approximately correspond to dips for $\Delta\varphi=\pi$ and vice versa. Such
a character can be used to design an optical switch or XNOR logic gate based
on phase difference. Next, we will study the average intensity $\bar{|v|^{2}}$
when the phase difference $\Delta\varphi$ changes continuously, and give the
theoretical scheme to realize the optical switch.

Figure 6: (a) The average intensity $\bar{|v|^{2}}$ as a function of the phase
difference $\Delta\varphi$ at $s=20$. (b) and (c) are the intensity $|v|^{2}$
as functions of $s$ and $\xi$, which describe the excitation of nonlinear SPPs
in the rectangular barrier.

Shown in Fig. 6(a) is the average intensity $\bar{|v|^{2}}$ as a function of
the phase difference $\Delta\varphi$ at $s=20$. Initial condition of the
nonlinear SPPs is $v={\rm sech}(\xi-12)e^{-i0.5\xi}+{\rm
sech}(\xi+12)e^{i0.5\xi+\Delta\varphi}$, and the height and half width are
$V_{0}=5$ and $l=4$, respectively. With the phase difference increase, we can
see that the average intensity $\bar{|v|^{2}}$ is large at first, and changes
slowly, then it drops rapidly to a small value (close to 0) in the region
$[\pi/6,\pi/2]$, and then keeps changing slowly, rises sharply in the region
$[3\pi/2,11\pi/6]$. This step change curve can be used to design an optical
switch.

Shown in Fig. 6(b) and (c) are intensity $|v|^{2}$ as functions of $s$ and
$\xi$, parameters are the same as that used in Fig. 6(a). In Fig. 6(b), we can
see that a large partial energy of the nonlinear SPPs couples into the barrier
when the phase difference $\Delta\varphi=0$, and the output strength in the
end of the barrier is quite large, and we denote this state as output ’1’.
However, when the phase difference $\Delta\varphi=\pi$, as shown in Fig. 6(c),
most energy of the nonlinear SPPs reflects out of the system, and the output
strength in the end of the barrier is approach to 0, and we denote this state
as output ’0’. Thus, a phase depended optical switch is obtained.

##  V Conclusion

In summary, we have proposed a scheme to study the propagation and excitation
properties of nonlinear SPPs under modulation of the rectangular potential
barrier. Based on multi scale method, we have derived the NLSE in the system,
and the Stark field applied into the system produced an equivalent rectangular
barrier to affect the propagation properties of the nonlinear SPPs. When
incident a single nonlinear SPPs, we have realized transmission, trapping and
reflection of the nonlinear SPPs. When two nonlinear SPPs symmetrically
inputting the system, we found that a periodic mode can be excited in the
rectangular barrier, and excitation of such modes are depending on
nonlinearity in the system, half width of the barrier and phase difference of
the initial nonlinear SPPs. Based on the related properties, we have designed
an optical switch based on phase difference of the nonlinear SPPs. These
results not only provide a theoretical basis for the study of the interaction
between nonlinear SPPs and external potentials, but also have broad
application prospects in the field of optical information at micro/nano scale.

## Acknowledgments

This work was supported by Natural Science Foundation of Shandong Province
(Grant No. ZR2021MA035). Xiangchun Tian and Yundong Zhang contribute to this
work equally.

## Appendix

##  Appendix A THE TM MODE IN THE WAVEGUIDE SYSTEM

The quantized probe field can be expressed as [30, 48]

${{\bf{E}}}\left({{\bf{r}},t}\right)=\left\\{{\begin{array}[]{*{20}{l}}{\left({k{{\bf{e}}_{x}}-i{k_{1}}{{\bf{e}}_{z}}}\right)\frac{c}{{{\varepsilon_{1}}{\omega_{l}}}}\sqrt{\frac{{\hbar{\omega_{l}}}}{{{\varepsilon_{0}}{L_{x}}{L_{y}}{L_{z}}}}}\hat{a}\left({{\omega_{p}}}\right){e^{-{k_{1}}x+i\left({kz-{\omega_{l}}t}\right)}}+c.c.}&{z>0}\\\
{\left({k{{\bf{e}}_{x}}+i{k_{2}}{{\bf{e}}_{z}}}\right)\frac{c}{{{\varepsilon_{2}}{\omega_{l}}}}\sqrt{\frac{{\hbar{\omega_{l}}}}{{{\varepsilon_{0}}{L_{x}}{L_{y}}{L_{z}}}}}\hat{a}\left({{\omega_{p}}}\right){e^{{k_{2}}x+i\left({kz-{\omega_{l}}t}\right)}}+c.c.}&{z<0.}\end{array}}\right.$
(14)

Here the expression of the wave vector is
$k_{j}^{2}=k^{2}-\omega_{l}^{2}{\varepsilon_{j}}{\mu_{j}}/{c^{2}}$( $j=1$ for
the hot atomic gas and $j=2$ for the NIMM) satisfies the boundary condition
${\varepsilon_{2}}{k_{1}}=-{\varepsilon_{1}}{k_{2}}$ at $x=0$. The propagation
constant of probe field can be expressed as
${k}(\omega_{l})=\left({{\omega_{l}}/c}\right){\left[{{\varepsilon_{1}}{\varepsilon_{2}}\left({{\varepsilon_{1}}{\mu_{2}}-{\varepsilon_{2}}{\mu_{1}}}\right)/\left({\varepsilon_{1}^{2}-\varepsilon_{2}^{2}}\right)}\right]^{1/2}}$.
$\hat{a}\left({{\omega_{p}}}\right)$ is the creation operator of TM photons.
$L_{y}$ and $L_{z}$ are the lengths of the NIMM-hot atomic gas interface in
the $y$ and $z$ directions, respectively. $L_{x}$ is the effective mode length
characterizing the electromagnetic field confinement in the $x$ direction
which is defined as

${L_{x}}\equiv\sum\limits_{j=1,2}{\frac{1}{{2{\mathop{\rm
Re}\nolimits}\left({{k_{j}}}\right)}}}\left[{\frac{{{{\tilde{\varepsilon}}_{j}}}}{{{{\left|{{\varepsilon_{j}}}\right|}^{2}}}}\frac{{{c^{2}}}}{{\omega_{l}^{2}}}\left({{{\left|k\right|}^{2}}+{{\left|{{k_{j}}}\right|}^{2}}}\right)+{{\tilde{\mu}}_{j}}}\right],$
(15)

with $\tilde{\varepsilon}\equiv{\mathop{\rm
Re}\nolimits}\left[{\partial\left({{\omega_{l}}\varepsilon}\right)/\partial{\omega_{l}}}\right]$,
$\tilde{\mu}\equiv{\mathop{\rm
Re}\nolimits}\left[{\partial\left({{\omega_{l}}\mu}\right)/\partial{\omega_{l}}}\right]$.

##  Appendix B THE BLOCH EQUATIONS IN THE INTERACTION PICTURE

The explicit expressions of the Bloch equations in the interaction picture can
be expressed as follows:

$\displaystyle{i\frac{\partial}{{\partial
t}}{\sigma_{11}}-i{\Gamma_{13}}{\sigma_{33}}+{e^{-i\theta_{p}^{*}}}{\zeta^{*}}\left(x\right)\Omega_{p}^{*}{\sigma_{31}}-{e^{i\theta_{p}}}\zeta\left(x\right){\Omega_{p}}\sigma_{31}^{*}=0},$
(16a) $\displaystyle{i\frac{\partial}{{\partial
t}}{\sigma_{22}}-i{\Gamma_{23}}{\sigma_{33}}+\Omega_{c}^{*}{\sigma_{32}}-{\Omega_{c}}\sigma_{32}^{*}=0},$
(16b) $\displaystyle{i\frac{\partial}{{\partial
t}}{\sigma_{33}}+i{{\Gamma_{3}}}{\sigma_{33}}-{e^{-i\theta_{p}^{*}}}{\zeta^{*}}\left(x\right)\Omega_{p}^{*}{\sigma_{31}}+{e^{i\theta_{p}}}\zeta\left(x\right){\Omega_{p}}\sigma_{31}^{*}-\Omega_{c}^{*}{\sigma_{32}}+{\Omega_{c}}\sigma_{32}^{*}=0},$
(16c) $\displaystyle{\left({i\frac{\partial}{{\partial
t}}+{d_{21}}}\right){\sigma_{21}}+\Omega_{c}^{*}{\sigma_{31}}-{e^{i\theta_{p}}}\zeta\left(x\right){\Omega_{p}}\sigma_{32}^{*}=0},$
(16d) $\displaystyle{\left({i\frac{\partial}{{\partial
t}}+{d_{31}}}\right){\sigma_{31}}+{e^{i\theta_{p}}}\zeta\left(x\right){\Omega_{p}}\left({{\sigma_{11}}-{\sigma_{33}}}\right)+{\Omega_{c}}{\sigma_{21}}=0},$
(16e) $\displaystyle{\left({i\frac{\partial}{{\partial
t}}+{d_{32}}}\right){\sigma_{32}}+{\Omega_{c}}\left({{\sigma_{22}}-{\sigma_{33}}}\right)+{e^{i\theta_{p}}}\zeta\left(x\right){\Omega_{p}}\sigma_{21}^{*}=0},$
(16f)

where
${d_{21}}=-\left({{\mathbf{k}_{p}}-{\mathbf{k}_{c}}}\right)\cdot\mathbf{v}+{{\Delta}_{2}}-{{\Delta}_{1}}+i{\gamma_{21}}+\frac{1}{2}{\alpha_{21}}{\left|{{E_{s}}}\right|^{2}},~{}{d_{31}}=-{\mathbf{k}_{p}}\cdot\mathbf{v}+{{\Delta}_{3}}-{{\Delta}_{1}}+i{\gamma_{31}}+\frac{1}{2}{\alpha_{31}}{\left|{{E_{s}}}\right|^{2}},~{}{d_{32}}=-{\mathbf{k}_{c}}\cdot\mathbf{v}+{{\Delta}_{3}}-{{\Delta}_{2}}+i{\gamma_{32}}+\frac{1}{2}{\alpha_{32}}{\left|{{E_{s}}}\right|^{2}}$,
with
${\Gamma_{j}}=\sum\nolimits_{E_{l}<E_{j}}{{\Gamma_{lj}}},~{}{\gamma_{jl}}=\left({{\Gamma_{j}}+{\Gamma_{l}}}\right)/2$.
Here ${{\Gamma_{j}}}$ is the total spontaneous emission decay rate and
${\gamma_{jl}}$ is the atomic decay rate.

##  Appendix C EXPRESSIONS OF EACH ORDER SOLUTIONS

###  C.1 First-order solutions

The first-order solutions reads

$\displaystyle\Omega_{p}^{\left(1\right)}={F}{e^{i{\theta}}},$ (17a)
$\displaystyle\sigma_{21}^{\left(1\right)}=-\frac{{\Omega_{c}^{*}}}{{{D}}}{F}{{e^{i{(\theta+\theta_{p})}}}}\zeta\left(x\right),$
(17b)
$\displaystyle\sigma_{31}^{\left(1\right)}=\frac{{\left({\omega+d_{21}^{\left(0\right)}}\right)}}{{{D}}}{F}{{e^{i{(\theta+\theta_{p})}}}}\zeta\left(x\right),$
(17c)

where
${\theta}={K}\left(\omega\right){z_{0}}-\omega{t_{0}},{{D}={{\left|{{\Omega_{c}}}\right|}^{2}}-\left({\omega+d_{21}^{\left(0\right)}}\right)\left({\omega+d_{31}^{\left(0\right)}}\right)}$,
with all other $\sigma_{ij}^{(1)}$ being zero. Here $F$ is the envelope
function of the slow variables $y_{1},z_{2}$ and $t_{2}$. The expression of
the linear dispersion relation $K(\omega)$ is

${K}\left(\omega\right)=\frac{\omega}{c}\frac{1}{{{n_{\rm
eff}}}}+{\kappa_{13}}\int_{-\infty}^{\infty}{d(kv)f\left(kv\right)\left\langle{\frac{{\left({\omega+d_{21}^{\left(0\right)}}\right)}}{{{D}}}\zeta\left(x\right)}\right\rangle}.$
(18)

###  C.2 Second-order solutions

The expressions of the second-order solutions are given by

$\displaystyle{\sigma_{11}^{\left(2\right)}=a_{11}^{\left(2\right)}{{\left|{\zeta\left(x\right)}\right|}^{2}}{{\left|{{F}}\right|}^{2}}{e^{-2{{\tilde{\alpha}}}{z_{2}}}}},$
(19a)
$\displaystyle{\sigma_{22}^{\left(2\right)}=a_{22}^{\left(2\right)}{{\left|{\zeta\left(x\right)}\right|}^{2}}{{\left|{{F}}\right|}^{2}}{e^{-2{{\tilde{\alpha}}}{z_{2}}}}},$
(19b)
$\displaystyle{\sigma_{32}^{\left(2\right)}=a_{32}^{\left(2\right)}{{\left|{\zeta\left(x\right)}\right|}^{2}}{{\left|{{F}}\right|}^{2}}{e^{-2{{\tilde{\alpha}}}{z_{2}}}}},$
(19c)
$\displaystyle{\sigma_{33}^{\left(2\right)}=a_{33}^{\left(2\right)}{{\left|{\zeta\left(x\right)}\right|}^{2}}{{\left|{{F}}\right|}^{2}}{e^{-2{{\tilde{\alpha}}}{z_{2}}}}},$
(19d)

where ${{\tilde{\alpha}}}={\epsilon^{-2}}{\mathop{\rm
Im}\nolimits}\left[{{{K}\left(\omega\right)+k(\omega_{p})}}\right]$ and

$\displaystyle
a_{11}^{\left(2\right)}=-a_{22}^{\left({\rm{2}}\right)}-a_{33}^{\left(2\right)},$
(20a) $\displaystyle
a_{22}^{\left(2\right)}=\frac{{\frac{{i{\Gamma_{23}}}}{{{{\left|{{\Omega_{c}}}\right|}^{2}}}}a_{33}^{\left(2\right)}+\frac{1}{{\omega+d_{32}^{*\left(0\right)}}}\left[{a_{33}^{\left(2\right)}+\frac{1}{{{D}}}}\right]-\frac{1}{{\omega+d_{32}^{\left(0\right)}}}\left[{a_{33}^{\left(2\right)}+\frac{1}{{D^{*}}}}\right]}}{{\frac{1}{{\omega+d_{32}^{*\left(0\right)}}}-\frac{1}{{\omega+d_{32}^{\left(0\right)}}}}},$
(20b) $\displaystyle
a_{32}^{\left(2\right)}=\frac{{{\Omega_{c}}}}{{\omega+d_{32}^{\left(0\right)}}}\left[{a_{33}^{\left(2\right)}-a_{22}^{\left(2\right)}+\frac{1}{{D^{*}}}}\right],$
(20c) $\displaystyle
a_{33}^{\left(2\right)}=\frac{1}{{i{\Gamma_{13}}}}\left[{\frac{{\omega{\rm{+}}d_{21}^{\left(0\right)}}}{{{D}}}-\frac{{\omega+d_{21}^{*\left(0\right)}}}{{D^{*}}}}\right].$
(20d)

###  C.3 Third-order solutions

In the third-order approximation, we obtain the explicit expression of $W_{1}$
and $W_{2}$ in Eq. (8), which read

$\displaystyle{{W_{1}}={\kappa_{13}}\int_{-\infty}^{\infty}{d(kv)f\left(kv\right)\left\langle{\frac{{\left({\omega+d_{21}^{\left(0\right)}}\right)\left({a_{11}^{\left(2\right)}-a_{33}^{\left(2\right)}}\right)+{\Omega_{c}}a_{32}^{*\left(2\right)}}}{{{D}}}\zeta\left(x\right){{\left|{\zeta\left(x\right)}\right|}^{2}}}\right\rangle}},$
(21a)
$\displaystyle{{W_{2}}={\kappa_{13}}\int_{-\infty}^{\infty}{d(kv)f\left(kv\right)\left\langle{\frac{{{{\left({\omega+d_{21}^{\left(0\right)}}\right)}^{2}}{\alpha_{31}}+{{\left|{{\Omega_{c}}}\right|}^{2}}{\alpha_{21}}}}{{2D^{2}}}\zeta\left(x\right)}\right\rangle}}.$
(21b)

## References

  * [1] P. Dombi, Z. Pápa, J. Vogelsang, S. V. Yalunin, M. Sivis, G. Herink, S. Schäfer, P.Groß, C. Ropers, and C. Lienau, Strong-field nano-optics, Rev. Mod. Phys. 92, 025003 (2020).
  * [2] E. Verhagen, L. Kuipers, and A. Polman, Enhanced nonlinear optical effects with a tapered plasmonic waveguide, Nano Letters 7, 334 (2007).
  * [3] J. A. Schuller, E. S. Barnard, W. Cai, Y. C. Jun, J. S. White, and M. L. Brongersma, Plasmonics for extreme light concentration and manipulation, Nature Mat. 9, 193 (2010).
  * [4] Q. Ren, J. W. You, and N. C. Panoiu, Large enhancement of the effective second-order nonlinearity in graphene metasurfaces, Phys. Rev. B 99, 205404 (2019).
  * [5] V. Kravtsov, S. AlMutairi, R. Ulbricht, A. R. Kutayiah, A. Belyanin, and M. B. Raschke, Enhanced third-order optical nonlinearity driven by surface-plasmon field gradients, Phys. Rev. Lett. 120, 203903 (2018)
  * [6] K. C. Y. Huang, M.-K. Seo, T. Sarmiento, Y. Huo, J. S. Harris, and M. L. Brongersma, Electrically driven subwavelength optical nanocircuits, Nature Photon. 8, 244 (2014).
  * [7] Y. Luo, E. D. Ahmadi, K. Shayan, Y. Ma, K. S. Mistry, C. Zhang, J. Hone, J. L. Blackburn, and S. Strauf, Purcell-enhanced quantum yield from carbon nanotube excitons coupled to plasmonic nanocavities, Nature Commun. 8, 1413 (2017).
  * [8] W. Choi, Y. Jo, J. Ahn, E. Seo, Q.-H. Park, Y. M. Jhon, and W. Choi, Control of randomly scattered surface plasmon polaritons for multiple-input and multiple-output plasmonic switching devices, Nature Commun. 8, 14636 (2017).
  * [9] H. C. Zhang, L. P. Zhang, P. H. He, J. Xu, C. Qian, F. J. Garcia-Vidal, and T. J. Cui, A plasmonic route for the integrated wireless communication of subdiffraction-limited signals, Light: Sci. & Appl. 9, 1 (2020).
  * [10] C. Ropers, C. C. Neacsu, T. Elsaesser, M. Albrecht, M. B. Raschke, and C. Lienau, Grating-coupling of surface plasmons onto metallic tips: a nanoconfined light source, Nano Letters 7, 2784 (2007).
  * [11] S. F. Becker, M. Esmann, K. Yoo, P. Gross, R. Vogelgesang, N. Park, and C. Lienau, Gap-plasmon-enhanced nanofocusing near-field microscopy, ACS Photonics 3, 223 (2016).
  * [12] L. Zhu, Y. Yin, L. Dai, Y. Hu, J. Nan, Z. Zheng, C. Cai, W. Zhao, and M. Ding, Round-tower plasmonic optical microfiber tip for nanofocusing with a high field enhancement, Optics Communications 453, 124358 (2019).
  * [13] C. Tan, N. Li, D. Xu, and Z. Chen, Spatial focusing of surface polaritons based on cross-phase modulation, Results in Physics 27, 104531 (2021).
  * [14] N. Fang, H. Lee, C. Sun, and X. Zhang, Sub-diffraction-limited optical imaging with a silver superlens, Science 308, 534 (2005).
  * [15] T. Umakoshi, M. Tanaka, Y. Saito, and P. Verma, White nanolight source for optical nanoimaging, Science Advances 6, 4179 (2020).
  * [16] B. Bai, X. Meng, J. Laukkanen, T. Sfez, L. Yu, W. Nakagawa, H. P. Herzig, L. Li, and J. Turunen, Asymmetrical excitation of surface plasmon polaritons on blazed gratings at normal incidence, Phys. Rev. B 80, 035407 (2009).
  * [17] H. Caglayan and E. Ozbay, Surface wave splitter based on metallic gratings with sub-wavelength aperture, Opt. Express 16, 19091 (2008).
  * [18] P. Berini and I. De Leon, Surface plasmon–polariton amplifiers and lasers, Nature Photon. 6, 16 (2012).
  * [19] A. Hohenau, J. R. Krenn, A. L. Stepanov, A. Drezet, H. Ditlbacher, B. Steinberger, A. Leitner, and F. R. Aussenegg, Dielectric optical elements for surface plasmons, Opt. Lett. 30, 893 (2005).
  * [20] M. A. Schmidt, L. N. Prill Sempere, H. K. Tyagi, C. G. Poulton, and P. S. J. Russell, Waveguiding and plasmon resonances in two-dimensional photonic lattices of gold and silver nanowires, Phys. Rev. B 77, 033417 (2008).
  * [21] J. Chen, C. Sun, H. Li, and Q. Gong, Ultra-broadband unidirectional launching of surface plasmon polaritons by a double-slit structure beyond the diffraction limit, Nanoscale 6, 13487 (2014).
  * [22] J. Lin, J. P. B. Mueller, Q. Wang, G. Yuan, N. Antoniou, X.-C. Yuan, and F. Capasso, Polarization-controlled tunable directional coupling of surface plasmon polaritons, Science 340, 331 (2013).
  * [23] I. Epstein and A. Arie, Arbitrary bending plasmonic light waves, Phys. Rev. Lett. 112, 023903 (2014).
  * [24] S. Dong, Y. Zhang, H. Guo, J. Duan, F. Guan, Q. He, H. Zhao, L. Zhou, and S. Sun, Highly efficient wave-front reshaping of surface waves with dielectric metawalls, Phys. Rev. Applied 9, 014032 (2018).
  * [25] B. Wang, X. Zhang, F. J. García-Vidal, X. Yuan, and J. Teng, Strong coupling of surface plasmon polaritons in monolayer graphene sheet arrays, Phys. Rev. Lett. 109, 073901 (2012).
  * [26] J. A. Gerber, S. Berweger, B. T. O’Callahan, and M. B. Raschke, Phase-resolved surface plasmon interferometry of graphene, Phys. Rev. Lett. 113, 055502 (2014).
  * [27] S. S. Sunku, G. X. Ni, B. Y. Jiang, H. Yoo, A. Sternbach, A. S. McLeod, T. Stauber, L. Xiong, T. Taniguchi, K. Watanabe, P. Kim, M. M. Fogler, and D. N. Basov, Photonic crystals for nano-light in moiré graphene superlattices, Science 362, 1153 (2018).
  * [28] L. Xiong, C. Forsythe, M. Jung, A. McLeod, S. Sunku, Y. Shao, G. Ni, A. Sternbach, S. Liu, J. Edgar, E. Mele, M. Fogler, G. Shvets, C. Dean, and D. Basov, Photonic crystal for graphene plasmons, Nature Commun. 10, 4780 (2019).
  * [29] Y. Wang, J. Yu, Y.-F. Mao, J. Chen, S. Wang, H.-Z. Chen, Y. Zhang, S.-Y. Wang, X. Chen, T. Li, L. Zhou, R.-M. Ma, S. Zhu, W. Cai, and J. Zhu, Stable, high-performance sodium-based plasmonic devices in the near infrared, Nature 581, 401 (2020).
  * [30] A. Kamli, S. A. Moiseev, and B. C. Sanders, Coherent control of low loss surface polaritons, Phys. Rev. Lett. 101, 263601 (2008).
  * [31] K. F. MacDonald, Z. L. Sámson, M. I. Stockman, and N. I. Zheludev, Ultrafast active plasmonics, Nature Photon. 3, 55 (2009).
  * [32] V. V. Temnov, G. Armelles, U. Woggon, D. Guzatov, A. Cebollada, A. Garcia-Martin, J.-M. Garcia-Martin, T. Thomay, A. Leitenstorfer, and R. Bratschitsch, Active magneto-plasmonics in hybrid metal–ferromagnet structures, Nature Photon. 4, 107 (2010).
  * [33] Q. Liu, N. Li, and C. Tan, All-optical logic gate based on manipulation of surface polaritons solitons via external gradient magnetic fields, Phys. Rev. A 101, 023818 (2020).
  * [34] Y. Zhou, Q. Liu, C.-k. Wang, and C. Tan, Trapping effect and trajectory control of surface plasmon polaritons in a metal-dielectric-metal waveguide, Phys. Rev. A 102, 063516 (2020).
  * [35] A. Melikyan, L. Alloatti, A. Muslija, D. Hillerkuss, P. C. Schindler, J. Li, R. Palmer, D. Korn, S. Muehlbrandt, D. Van Thourhout, B. Chen, R. Dinu, M. Sommer, C. Koos, M. Kohl, W. Freude, and J. Leuthold, High-speed plasmonic phase modulators, Nature Photon. 8, 229 (2014).
  * [36] C. Haffner, W. Heni, Y. Fedoryshyn, J. Niegemann, A. Melikyan, D. L. Elder, B. Baeuerle, Y. Salamin, A. Josten, U. Koch, C. Hoessbacher, F. Ducry, L. Juchli, A. Emboras, D. Hillerkuss, M. Kohl, L. R. Dalton, C. Hafner, and J. Leuthold, All-plasmonic Mach–Zehnder modulator enabling optical high-speed communication at the microscale, Nature Photon. 9, 525 (2015).
  * [37] I. V. Dzedolik, and V. Pereskokov, Control of plasmon-polariton vortices on the surface of a metal layer, J. Opt. Soc. Am. A 35, 1420 (2018).
  * [38] J.-C. Weeber, T. Bernardin, M. G. Nielsen, K. Hassan, S. Kaya, J. Fatome, C. Finot, A. Dereux, and N. Pleros, Nanosecond thermo-optical dynamics of polymer loaded plasmonic waveguides, Opt. Express 21, 27291 (2013).
  * [39] Q. Li, L. Chen, H. Xu, Z. Liu, and H. Wei, Photothermal modulation of propagating surface plasmons on silver nanowires, ACS Photonics 6, 2133 (2019).
  * [40] Y. Kominis and K. Hizanidis, Power-dependent reflection, transmission, and trapping dynamics of lattice solitons at interfaces, Phys. Rev. Lett. 102, 133903 (2009).
  * [41] Z. Bai, Q. Zhang, and G. Huang, Quantum reflections of nonlocal optical solitons in a cold rydberg atomic gas, Phys. Rev. A 101, 053845 (2020).
  * [42] E. Figueroa, F. Vewinger, J. Appel, and A. I. Lvovsky, Decoherence of electromagnetically induced transparency in atomic vapor, Opt. Lett. 31, 2625 (2006).
  * [43] G. Huang, L. Deng, and M. G. Payne, Dynamics of ultraslow optical solitons in a cold three-state atomic system, Phys. Rev. E 72, 016617 (2005).
  * [44] Y. Guo, L. Zhou, L.-M. Kuang, and C. P. Sun, Magneto-optical stern-gerlach effect in an atomic ensemble, Phys. Rev. A 78, 013833 (2008).
  * [45] C. Hang and G. Huang, Stern-gerlach effect of weak-light ultraslow vector solitons, Phys. Rev. A 86, 043809 (2012).
  * [46] D. A. Steck, Rubidium 87 d line data, (2001).
  * [47] Z. Chen, H. Xie, Q. Li, and G. Huang, Stern-gerlach deflection of optical thirring solitons in a coherent atomic system, Phys. Rev. A 100, 013827 (2019).
  * [48] C. Tan and G. Huang, Surface polaritons in a negative-index metamaterial with active raman gain, Phys. Rev. A 91, 023803 (2015).