## Abstract

Using a combination of numerical and analytical techniques we demonstrate that a metal stripe surrounded by the active and passive dielectrics supports propagation of stable spatial surface-plasmon solitons. Our analytical methods include the multiple scale reduction of the Maxwell’s equations to the coupled Ginzburg-Landau system, and the soliton perturbation theory developed in the framework of the latter.

© 2011 Optical Society of America

## 1. Introduction

Surface plasmon polaritons (SPPs) provide one of the favored approaches to realization of on-chip photonic devices and are a well established tool in sensing applications. While SPPs are exponentially localized in the direction perpendicular to the metal-dielectric interface by the natural boundary conditions, one should take a special care about suppression of their in-plane diffraction [1]. An interesting alternative to various geometrical methods providing lateral confinement of SPPs is to use the concept of spatial solitons, where diffraction is suppressed by the nonlinearity induced focusing, see, e.g., [2, 3]. Characteristic loss length of SPPs is typically order of few tens of microns. In order to talk about diffraction compensation over such distances one should work with the beams having diffraction length shorter than the loss length. Such strong diffraction requires very large powers to achieve the soliton effect. One can structure the metal surface and thus weaken the diffraction and hence power requirements, and work with the discrete SPP solitons [4, 5].

On the other hand, the spatial soliton concept can be extended by complementing the diffraction vs nonlinearity balance with the gain vs loss balance. Thus completely solving the problem of the soliton decay due to linear absorption. The paradigm model in this context is the complex Ginzburg-Landau (GL) equation [6, 7]. Cubic GL equation has been derived for the SPPs at the boundary with an active dielectric [7]. However, SPP solitons reported in the above work demonstrate substantial instabilities developing, both, in the soliton core and background [7]. The nature of these instabilities is the same as the one known in other physical contexts with the supercritical transition between the zero and non-zero spatially homogeneous solutions [6]. Making this transition subcritical is known to favor stabilization of solitons. This can be achieved, e.g., if the quintic nonlinearity is accounted for [6, 8, 9]. However, the higher order nonlinearities usually become important only for sufficiently high intensities. An alternative method to suppress the soliton instabilities is to drain the energy generated by the gain in the cubic GL equation using the coupling of the latter to the second purely lossy GL equation [10–12].

Below we consider a thin metal stripe, where SPPs at the two interfaces are weakly coupled. Dielectric at one of the interfaces is assumed active and at the other is passive. Using the method of multiple scales applied to the Maxwell equations we derive the system of two coupled GL equations for the SPP amplitudes and report stable spatial SPP solitons. Note, that the amplification of SPPs without considering transverse effects has been previously studied theoretically [13, 14] and experimentally using rhodamine 6G, quantum dots and erbium [15–17].

## 2. Reduction of the Maxwell’s equations to the coupled Ginzburg-Landau system

We start our analysis from the time independent Maxwell equations ∇⃗ × ∇⃗ × *E*⃗ = *ɛE*⃗. The spatial coordinates are normalized to *k* = *λ*/(2*π*), where *λ* is the vacuum wavelength. *x* is the coordinate perpendicular to the metal-dielectric interfaces, *z* is the propagation coordinate, while the in-plane diffraction of SPPs happens along *y* direction, see Fig. 1(a). The metal stripe has the width *w* along the *x*-axis and is infinite along *y* and *z*. *ɛ _{metal}* =

*ɛ*″

*+*

_{m}*iɛ*″

*is the dielectric constant inside the metal (*

_{m}*ɛ*′

*< 0,*

_{m}*ɛ*″

*> 0). Inside the active layer at the top of the metal stripe we have ${\varepsilon}_{\mathit{active}}=\left({\varepsilon}_{d}+\frac{1}{2}\alpha \right)+i{\varepsilon}_{a}^{\u2033}+{\chi}_{3}{\left|\overrightarrow{E}\right|}^{2}$. Here*

_{m}*χ*

_{3}=

*χ*′

_{3}+

*iχ*″

_{3},

*χ*″

_{3}> 0 and

*ɛ*″

*< 0 corresponding to the nonlinear absorption and linear gain, respectively. Our model of the active dielectric corresponds to the two-level atom medium in the regime far from saturation [7].*

_{a}*ɛ*in the passive dielectric layer is ${\varepsilon}_{\mathit{passive}}={\varepsilon}_{d}-\frac{1}{2}\alpha $. We have expressed the real parts of the linear dielectric constants of the top layer and of the substrate as ${\varepsilon}_{d}\pm \frac{1}{2}\alpha $, because it is more straightforward to assume that in the zero approximation of our perturbation theory the SPPs at both interfaces are identical, so that the

*α*terms appear only in the higher orders. This assumption helps to avoid unnecessary complexity of the analytical calculations.

Assuming that *ɛ*″* _{m}*,

*α*,

*ɛ*″

*and |*

_{a}*χ*

_{3}||

*E⃗*|

^{2}are of the order

*s*(where

*s*is the small dummy parameter) and that the width

*w*of the metal stripe is much greater than the skin depth [5], we have the following expressions for the SPPs at the two interfaces:

*a*and

*p*mark the interfaces adjoint with the active and passive layers, respectively.

*θ*is the Heaviside function, ${q}_{d}\equiv \sqrt{{\beta}^{2}-{\varepsilon}_{d}}$, ${q}_{m}\equiv \sqrt{{\beta}^{2}-{\varepsilon}_{m}^{\prime}}$ are the decay rates of the plasmonic tails inside the dielectrics and metal respectively, $\beta =\sqrt{{\varepsilon}_{d}{\varepsilon}_{m}^{\prime}/\left({\varepsilon}_{d}+{\varepsilon}_{m}^{\prime}\right)}$ is the dimensionless SPP propagation constant. A solution of the Maxwell’s equations is then sought in the form

*ψ*and

_{p,a}*ϕ*are slow functions of their arguments, such that

_{p,a}*∂*∼

_{z}*O*(

*s*),

*∂*∼

_{y}*O*(

*s*

^{1/2}). |

*ψ*| ∼

_{p,a}*O*(

*s*

^{1/2}) and |

*ϕ*| ∼

_{p,a}*O*(

*s*), so that the

*y*-component of the field existing due to the finite beam size is smaller than the other two [18]. Proceeding to the order

*s*

^{3/2}one obtains a system of equations for the correction

*δE⃗*= (

*δE*,

_{x}*δE*)

_{z}*(|*

^{T}*δE⃗*| ∼

*O*(

*s*

^{3/2})):

*L̂δE⃗*=(

*N̂*–

*L̂*)(

*ψ*+

_{p}e⃗_{p}*ψ*) –

_{a}e⃗_{a}*b⃗*, where ${\overrightarrow{e}}_{p,a}={\left({e}_{x}^{p,a},{e}_{x}^{p,a}\right)}^{T}$ and

*L̂*in the left hand side of the equation for

*δE⃗*into the two parts independently associated with the two interfaces.

The propagation equations for the amplitudes *ψ _{p}* and

*ψ*can be derived by projecting the left and right hand side of the equation for

_{a}*δE⃗*onto

*e⃗*and assuming that the interaction of the two SPPs through the overlap of their tails becomes important only in the

_{p,a}*s*

^{3/2}order. Note, that the calculation of the projection integrals and use of the boundary conditions for all the components of the electric and magnetic fields reveal that the left hand side yields a nontrivial contribution due to discontinuities of

*E*and finite values of

_{x}*E*at the interfaces [18]. The terms containing

_{z}*ϕ*partly cancel out and partly can be expressed via

_{p,a}*ψ*. The resulting system is the two coupled complex GL equations:

_{p,a}*l*=

*β*

^{3}

*ɛ*″

*/[2(*

_{m}*ɛ*′

*)*

_{m}^{2}], $g=-{\beta}^{3}{\varepsilon}_{a}^{\u2033}/\left(2{\varepsilon}_{d}^{2}\right)$, and

*κ*= 2

*β*

^{3}

*e*

^{−qmw}/(

*ɛ*–

_{d}*ɛ*′

*) are, respectively, the normalized loss, gain and coupling coefficients. The associated physical lengths are calculated as 1/(*

_{m}*lk*), 1/(

*gk*) and 1/(

*kκ*).

*k*Δ

*β*is the small difference between the SPP propagation constants at the two interfaces. Choosing the scaling factor

*f*=

*P/k*, we have dimensionless

*ψ*, such that |

_{p,a}*ψ*|

_{p,a}^{2}= 1 correspond to the linear power density of

*P*Watts/m. The complex nonlinear SPP parameter ϒ = ϒ′ +

*i*ϒ″ is given by

*χ*

_{3}=

*χ*′

_{3}+

*iχ*″

_{3}has been expressed in terms of the more conventional complex Kerr coefficient

*n*

_{2}measured in m

^{2}/W. For the model of two level atoms

*χ*

_{3}= 2

*cn*

_{2}

*ɛ*

_{0}

*ɛ*, where

_{d}*c*is the vacuum speed of light and

*ɛ*

_{0}is the vacuum permittivity.

## 3. Subcritical transition to the spatially homogeneous SPPs

If coupling between the two interfaces is neglected, *κ* = 0, then the trivial solution *ψ _{a}* = 0 loses its linear stability at

*g*=

*l*, in favor of |

*ψ*|

_{a}^{2}= (

*g*–

*l*)/(

*f*ϒ″). The latter obviously exists only for

*g*>

*l*, so that the transition between the trivial and nontrivial states is supercritical. Thus any SPP soliton existing on the zero background for

*g*>

*l*suffers from the instability of its tails [7].

The transition between the zero and non-zero SPPs is however more involved for *κ* ≠ 0. The trivial solution *ψ _{a}* =

*ψ*= 0 becomes unstable for

_{p}*g*>

*g*, where the latter solves the cubic equation

_{th}*ψ*=

_{p}*Se*,

^{i}^{μ}^{z}*ψ*=

_{a}*Re*exists however not only for

^{iμz}*g*>

*g*, but also for

_{th}*g*

_{0}<

*g*<

*g*. Thus, the transition in this case is subcritical. While an expression for

_{th}*g*

_{0}can not be found in a simple analytical form, the Fig. 1(b) shows the range of

*g*and Δ

*β*values, where the non-zero solution coexist with the zero one. One can see, that the range of the suitable gain is maximal for Δ

*β*= 0. An example of the subcritical dependence of |

*R*|

^{2}vs

*g*is shown in Fig. 1(c). From this analysis we conjecture that in the subcritical case the soliton solutions on the zero background, if they exist, will have stable tails. In our numerical examples we assume

*ɛ*= 1.8, the metal is silver, and

_{d}*λ*= 530nm, so that

*ɛ*′

*= −15,*

_{m}*ɛ*″

*= 0.4. Under these condition the dimensionless gain values used in Figs. 1 and 2 are scaled back to the physical units by multiplying them with ∼ 12*

_{m}*μ*m

^{−1}. Thus to see the subcritical transition the gain should reach values boosting the intensity by the factor

*e*after ∼ 20

*μ*m of the propagation distance. For comparison, the gain length demonstrated with dyes [15] was close to this value, while in the quantum dot case [16] it was one order of magnitude larger. The intensity plot shown in Fig. 1(c) corresponds to the nonlinear response of the material and pump powers such that the nonlinearity induced correction to the refractive index (given by

*f*ϒ|

*R*|

^{2}) has a realistic value ∼ 10

^{−3}.

## 4. Spatial solitons

To analyse the existence of the soliton solutions we develop an approximate analytical theory and use a numerical approach based on the Newton-Raphson iterations. Differentiating the power integral $N={\int}_{-\infty}^{+\infty}dy\left({\left|{\psi}_{a}\right|}^{2}+{\left|{\psi}_{p}\right|}^{2}\right)$ and using Eqs. (2) we find the following balance equation

*z*parameter, one can use Eq. (4) to derive an evolution equation for this parameter and determine its stationary values. For the single cubic-quintic GL, the soliton of the ideal nonlinear Schrodinger equation has been successfully used in this approach [8, 9]. Neglecting all the dissipative terms

*l*=

*g*= ϒ″ = 0 in Eqs. (2), and assuming Δ

*β*= 0, still leaves us with the system, where the solitons can not be found in a closed analytical form. Therefore we need to make further simplifying assumption that

*κ*is a small parameter and |

*ψ*| ≫ |

_{a}*ψ*|. In this case

_{p}*μ*=

*η*

^{2}/2

*β*and

*η*> 0 parameterizes the soliton family. Using Eq. (4), we find the equation for

*η*:

*C*≃ 5.694.

Stationary (*dη*/*dz* = 0) values of *η* obey the cubic equation for *η*^{2}, which roots predict the soliton families. Within the range of the subcritical behavior of the homogeneous solution we have found that there can be up to two positive solutions for *η*^{2} corresponding to the two different solitons. Maximum of |*ψ _{a}*| corresponding to these two branches is plotted as a function of the gain parameter

*g*, in the Fig. 2(a) with the full lines. The soliton profile corresponding to the large amplitude branch is shown in Fig. 2(b) with the full line.

Although the soliton profiles are predicted by the perturbation theory quite accurately, they are not exact. Therefore, the use of numerical methods is necessary. We seek solitons in the form *ψ _{a,p}* =

*f*(

_{a,p}*y*)

*e*, approximate ${\partial}_{yy}^{2}$ with the 2nd order finite differences, and apply the Newton-Raphson iterations to the resulting set of the algebraic equations. Parameter

^{iμz}*μ*is treated as one of the unknowns, while the perturbation theory provides an initial guess for the iterative procedure. Note that, since

*μ*is an unknown, the system formally contains

*N*equations and

*N*+ 1 unknowns. However, we reduce the number of unknowns by one, fixing Im

*ψ*to zero at

_{a}*y*= 0, which always can be achieved by an appropriate selection of the arbitrary common phase of

*ψ*and

_{p}*ψ*. Numerically found branches of the solitons and an example of the soliton profile are shown in Fig. 2(a) and 2(b), respectively, alongside with the ones found using our analytical approach. Scaling back into physical units one finds that the soliton width is order of 1

_{a}*μ*m.

Very good quantitative agreement between our numerical and analytical approaches is evident. Note, however, that the perturbation theory predicts that the branch of the small amplitude solitons exists not only in the subcritical range (*g* < *g*_{th}), but also for *g* > *g*_{th}. While the former prediction is well confirmed by the numerical results, the latter one is not. The large amplitude solitons exist on both sides from *g _{th}*, while their tails remain stable only for

*g*<

*g*. We have also solved Eqs. (2) directly using a split-step method and have found that the large amplitude branch of the solitons remains stable over and well beyond practical propagation distances providing

_{th}*g*<

*g*, while the small amplitude ones are unstable. The instability of the latter is not unexpected, since the small amplitude solitons separate the attraction basins of the stable zero solution and of the stable large-amplitude soliton.

_{th}## 5. Summary

In this work we demonstrate that the ubiquitous dielectric-metal-dielectric geometry can support propagation of *stable* spatial plasmon solitons. The stability is achieved by means of the coupling of the actively pumped SPP at one interface to the passive (lossy) one at the other. The lossy SPP absorbs the perturbations that would otherwise destabilize the zero background of the soliton. We have used the multiple scale approach to demonstrate the reduction of the Maxwell’s equations to the analytically and numerically tractable system of the coupled Ginzburg-Landau equations. We have found the gain interval of the coexistence of the stable zero amplitude solution with the stable spatially homogeneous SPP. Applying perturbative and numerical methods we have proved that within this interval there exist stable spatial plasmon solitons, such that the in-plane diffraction and absorption are both balanced through the interplay of the linear gain and complex nonlinearity, which was assumed coming from the dielectric. From our analysis it follows that, unlike the dielectric-metal-dielectric case, the metal-dielectric-metal geometry is not favorable for the existence of this type of solitons, while future studies in the multilayered structures can lead to interesting results.

## Acknowledgments

We acknowledge useful comments from A. V. Gorbach and W. J. Firth.

## References and links

**1. **D. K. Gramotnev and S. I. Bozhevolnyi, “Plasmonics beyond the diffraction limit,” Nat. Photonics **4**, 83–91 (2010). [CrossRef]

**2. **E. Feigenbaum and M. Orenstein, “Plasmon-solitons,” Opt. Lett. **32**, 674–676 (2007). [CrossRef] [PubMed]

**3. **A. R. Davoyan, I. V. Shadrivov, and Y. S. Kivshar, “Self-focusing and spatial plasmon-polariton solitons,” Opt. Express **17**, 21732–21737 (2009). [CrossRef] [PubMed]

**4. **Y. Liu, G. Bartal, D. A. Genov, and X. Zhang, “Subwavelength Discrete Solitons in Nonlinear Metamaterials,” Phys. Rev. Lett. **99**, 153901 (2007). [CrossRef] [PubMed]

**5. **A. Marini, A. V. Gorbach, and D. V. Skryabin, “Coupled-mode approach to surface plasmon polaritons in nonlinear periodic structures,” Opt. Lett. **35**, 3532–3534 (2010). [CrossRef] [PubMed]

**6. **N. N. Ahmediev and A. Ankiewicz, *Solitons: Nonlinear Pulses and Beams* (Chapman and Hall, 2007).

**7. **A. Marini and D. V. Skryabin, “Ginzburg-landau equation bound to the metal-dielectric interface and transverse nonlinear optics with amplified plasmon polaritons,” Phys. Rev. A **81**, 033850 (2010). [CrossRef]

**8. **B. A. Malomed, “Evolution of nonsoliton and quasiclassical wavetrains in nonlinear Schrdinger and Korteweg -de Vries equations with dissipative perturbations,” Physica D **29**, 155–172 (1987). [CrossRef]

**9. **S. Fauve and O. Thual, “Solitary waves generated by subcritical instabilities in dissipative systems,” Phys. Rev. Lett. **64**, 282–284 (1990). [CrossRef] [PubMed]

**10. **B. A. Malomed and H. G Winful. “Stable solitons in two-component active systems,” Phys. Rev. E **53**, 5365 (1996). [CrossRef]

**11. **J. Atai and B. A. Malomed, “Stability and interactions of solitons in two-component active systems,” Phys. Rev. E **54**, 4371 (1996). [CrossRef]

**12. **W. J. Firth and P. V. Paulau, “Soliton lasers stabilized by coupling to a resonant linear system,” Eur. Phys. J. D **59**, 13–21 (2010). [CrossRef]

**13. **D. J. Bergman and M. I. Stockman“Surface plasmon amplification by stimulated emission of radiation: Quantum generation of coherent surface plasmons in nanosystems,” Phys. Rev. Lett. **90**, 027402 (2003). [CrossRef] [PubMed]

**14. **M. P. Nezhad, K. Tetz, and Y. Fainman, “Gain assisted propagation of surface plasmon polaritons on planar metallic waveguides,” Opt. Express **12**, 4072–4079 (2004). [CrossRef] [PubMed]

**15. **M. A. Noginov, V. A. Podolskiy, G. Zhu, M. Mayy, M. Bahoura, J. A. Adegoke, B. A. Ritzo, and K. Reynolds, “Compensation of loss in propagating surface plasmon polariton by gain in adjacent dielectric medium,” Opt. Express **16**, 1385 (2008). [CrossRef] [PubMed]

**16. **P. M. Bolger, W. Dickson, A. V. Krasavin, L. Liebscher, S. G. Hickey, D. V. Skryabin, and A. V. Zayats, “Amplified spontaneous emission of surface plasmon polaritons and limitations on the increase of their propagation length,” Opt. Lett. **35**, 1197–1199 (2010). [CrossRef] [PubMed]

**17. **M. Ambati, S. H. Nam, E. Ulin-Avila, D. A. Genov, G. Bartal, and X. Zhang, “Observation of stimulated emission of surface plasmon polaritons,” Nano Lett. **8**, 3998–4001 (2008). [CrossRef] [PubMed]

**18. **D. V. Skryabin, A. Gorbach, and A. Marini, “Surface induced nonlinearity enhancement of TM-modes in planar subwavelength waveguides,” J. Opt. Soc. Am. B **28**, 109–114 (2011). [CrossRef]