Stable spatial plasmon solitons in a dielectric-metal-dielectric geometry with gain and loss

Andrea Marini; Dmitry V. Skryabin; Boris Malomed

doi:10.1364/OE.19.006616

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ɛ″_m is the dielectric constant inside the metal (ɛ′_m < 0, ɛ″_m > 0). Inside the active layer at the top of the metal stripe we have $ɛ_{active} = (ɛ_{d} + \frac{1}{2} α) + i ɛ_{a}^{″} + χ_{3} {| \vec{E} |}^{2}$ . Here χ₃ = χ′₃ + iχ″₃, χ″₃ > 0 and ɛ″_a < 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]. ɛ in the passive dielectric layer is $ɛ_{passive} = ɛ_{d} - \frac{1}{2} α$ . We have expressed the real parts of the linear dielectric constants of the top layer and of the substrate as $ɛ_{d} \pm \frac{1}{2} α$ , 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.

Fig. 1 (a) Dielectric-metal-dielectric structure. (b) The coexistence domain of the stable zero and nonzero spatially homogeneous SPPs in the (g,Δβ)-plane. g_th is indicated with the full lines and g₀ with the dashed one. Other parameters are β = 1.43, l = 0.0026; fϒ = 3.5 × 10⁻³(1 + 0.1i); κ = 0.0028 corresponds to w = 98nm. (b) Subcritical dependence of the intensity |R|² of the spatially homogeneous SPP vs gain g, Δβ = −2.5 × 10⁻⁴.

Download Full Size | PDF

Assuming that ɛ″_m, α, ɛ″_a and |χ₃||E⃗|² 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:

{\vec{e}}_{a} = (\begin{matrix} e_{x}^{a} \\ e_{z}^{a} \end{matrix}) = θ (x + w / 2) (\begin{matrix} i \frac{β}{q_{m}} \\ 1 \end{matrix}) e^{- q_{m} (x + w / 2)} + θ (- x - w / 2) (\begin{matrix} - i \frac{β}{q_{d}} \\ 1 \end{matrix}) e^{q_{d} (x + w / 2)},

{\vec{e}}_{p} = (\begin{matrix} e_{x}^{p} \\ e_{z}^{p} \end{matrix}) = θ (x - w / 2) (\begin{matrix} i \frac{β}{q_{d}} \\ 1 \end{matrix}) e^{- q_{d} (x - w / 2)} + θ (- x + w / 2) (\begin{matrix} - i \frac{β}{q_{m}} \\ 1 \end{matrix}) e^{q_{m} (x - w / 2)},

where a and p mark the interfaces adjoint with the active and passive layers, respectively. θ is the Heaviside function,

q_{d} \equiv \sqrt{β^{2} - ɛ_{d}}

,

q_{m} \equiv \sqrt{β^{2} - ɛ_{m}^{'}}

are the decay rates of the plasmonic tails inside the dielectrics and metal respectively,

β = \sqrt{ɛ_{d} ɛ_{m}^{'} / (ɛ_{d} + ɛ_{m}^{'})}

is the dimensionless SPP propagation constant. A solution of the Maxwell’s equations is then sought in the form

\begin{array}{l} E_{x} = [ψ_{p} (y, z) e_{x}^{p} (x) + ψ_{a} (y, z) e_{x}^{a} (x) + δ E_{x} + \dots] e^{i β z}, \\ E_{y} = [ϕ_{p} (y, z) e_{y}^{p} (x) + ϕ_{a} (y, z) e_{y}^{a} (x) + \dots] e^{i β z}, \\ E_{z} = [ψ_{p} (y, z) e_{z}^{p} (x) + ψ_{a} (y, z) e_{z}^{a} (x) + δ E_{z} + \dots] e^{i β z}, \end{array}

where the amplitudes ψ_p,a and ϕ_p,a are slow functions of their arguments, such that ∂_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̂)(ψ_pe⃗_p + ψ_ae⃗_a) – b⃗, where

{\vec{e}}_{p, a} = {(e_{x}^{p, a}, e_{x}^{p, a})}^{T}

and

\begin{array}{l} \hat{L} = (\begin{matrix} β^{2} - \tilde{ɛ} & i β \partial_{x} \\ i β \partial_{x} & - \partial_{x x}^{2} - \tilde{ɛ} \end{matrix}), \hat{N} = (\begin{matrix} 2 i β \partial_{z} + \partial_{y y}^{2} + δ \tilde{ɛ} & - \partial_{x z}^{2} \\ - \partial_{x z}^{2} & \partial_{y y}^{2} + δ \tilde{ɛ} \end{matrix}), \\ \tilde{ɛ} = ɛ_{m}^{'} [θ (x + w / 2) - θ (x - w / 2)] + ɛ_{d} [θ (- x - w / 2) + θ (x - w / 2)], \\ δ ɛ = {i ɛ}_{m}^{″} [θ (x + w / 2) - θ (x - w / 2)] + (1 / 2) α [θ (x - w / 2) - θ (- x - w / 2)] \\ + [i ɛ_{a}^{″} + χ_{3} {| ψ_{a} |}^{2} {| {\vec{e}}_{a} |}^{2}] θ (x - w / 2), \vec{b} = (\begin{matrix} \partial_{y} ϕ_{p} \partial_{x} e_{y}^{p} + \partial_{y} ϕ_{a} \partial_{x} e_{y}^{a} \\ i β \partial_{y} ϕ_{p} e_{y}^{p} + i β \partial_{y} ϕ_{a} e_{y}^{a} \end{matrix}) . \end{array}

Note, that we have avoided (without any loss of rigour) the splitting of the operator 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 ψ_a can be derived by projecting the left and right hand side of the equation for δE⃗ onto e⃗_p,a and assuming that the interaction of the two SPPs through the overlap of their tails becomes important only in the 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_x and finite values of E_z at the interfaces [18]. The terms containing ϕ_p,a partly cancel out and partly can be expressed via ψ_p,a. The resulting system is the two coupled complex GL equations:

\begin{array}{l} i \partial_{z} ψ_{p} + (1 / 2 β) \partial_{y y}^{2} ψ_{p} + (i l - Δ β) ψ_{p} + κ ψ_{a} = 0, \\ i \partial_{z} ψ_{a} + (1 / 2 β) \partial_{y y}^{2} ψ_{a} + [i (l - g) + Δ β] ψ_{a} + f ϒ {| ψ_{a} |}^{2} ψ_{a} + κ ψ_{p} = 0, \end{array}

where l = β³ɛ″_m/[2(ɛ′_m)²],

g = - β^{3} ɛ_{a}^{″} / (2 ɛ_{d}^{2})

, and κ = 2β³e^−q_mw/(ɛ_d – ɛ′_m) are, respectively, the normalized loss, gain and coupling coefficients. The associated physical lengths are calculated as 1/(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 ψ_p,a, such that |ψ_p,a|² = 1 correspond to the linear power density of P Watts/m. The complex nonlinear SPP parameter ϒ = ϒ′ + iϒ″ is given by

ϒ = k^{2} n_{2} \frac{β^{3} (ɛ_{d} - ɛ_{m}^{'})}{2 ɛ_{d}^{2}},

and is measured in 1/W, where χ₃ = χ′₃ + iχ″₃ has been expressed in terms of the more conventional complex Kerr coefficient n₂ measured in m²/W. For the model of two level atoms χ₃ = 2cn₂ɛ₀ɛ_d, where c is the vacuum speed of light and ɛ₀ 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|² = (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 = ψ_p = 0 becomes unstable for g > g_th, where the latter solves the cubic equation

(g_{th} - l - κ^{2} / l) {(g_{th} - 2 l)}^{2} = 4 Δ β^{2} (l - g_{th}) .

The non-zero spatial homogeneous solution ψ_p = Seⁱ^μ^z, ψ_a = Re^iμz exists however not only for g > g_th, but also for g₀ < g < g_th. Thus, the transition in this case is subcritical. While an expression for g₀ 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|² 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 ɛ_d = 1.8, the metal is silver, and λ = 530nm, so that ɛ′_m = −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⁻¹. 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|²) has a realistic value ∼ 10⁻³.

Fig. 2 (a) Maximum of the soliton intensity, max|ψ_a|² vs gain g. The crosses correspond to the Newton-Raphson method, while the full line corresponds to the analytical results. The dotted line correspond to the spatially homogeneous SPPs. The other parameters are as in Fig. 1(c). (b) Dots show the numerically computed soliton profiles |ψ_p(y)| (blue) and |ψ_a(y)| (red). Full lines are the soliton profiles predicted by the perturbation theory. The solitons shown correspond to the large amplitude branch. g = 0.0042 and the other parameters as in Fig. 1(c).

Download Full Size | PDF

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} d y ({| ψ_{a} |}^{2} + {| ψ_{p} |}^{2})$ and using Eqs. (2) we find the following balance equation

\frac{d N}{d z} + 2 l \int_{- \infty}^{+ \infty} d y {| ψ_{p} |}^{2} - 2 (g - l) \int_{- \infty}^{+ \infty} d y {| ψ_{a} |}^{2} + 2 f ϒ^{″} \int_{- \infty}^{+ \infty} d y {| ψ_{a} |}^{4} = 0,

which transforms into the power conservation law in the absence of gain and linear an nonlinear losses. Taking a reasonable trial function for a soliton solution containing a single adiabatically evolving in 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| ≫ |ψ_p|. In this case

ψ_{a} = η \sqrt{1 / (β f ϒ^{'})} sech (η y) e^{i μ z} + O (κ^{2}),

ψ_{p} = \frac{2 κ}{η} \sqrt{\frac{β}{f ϒ^{'}}} {cosh (η y) ln [2 cosh (η y)] - η y sinh (η y)} e^{i μ z} + O (κ^{3}),

where μ = η²/2β and η > 0 parameterizes the soliton family. Using Eq. (4), we find the equation for η:

\frac{d η}{d z} + β^{2} κ^{2} l C η^{- 3} - 2 (g - l) η + \frac{4 ϒ^{″}}{3 β ϒ^{'}} η^{3} = 0,

where C ≃ 5.694.

Stationary (dη/dz = 0) values of η obey the cubic equation for η², 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 η² 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^iμz, approximate $\partial_{y y}^{2}$ with the 2nd order finite differences, and apply the Newton-Raphson iterations to the resulting set of the algebraic equations. Parameter μ 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ψ_a to zero at y = 0, which always can be achieved by an appropriate selection of the arbitrary common phase of ψ_p and ψ_a. 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μ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_th. 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 g < g_th, 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.

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]

Stable spatial plasmon solitons in a dielectric-metal-dielectric geometry with gain and loss

Abstract

1. Introduction

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

3. Subcritical transition to the spatially homogeneous SPPs

4. Spatial solitons

5. Summary

Acknowledgments

References and links

Cited By

Figures (2)

Equations (11)

Optics Express