1. Introduction

Plasmonics has developed into one of the most active research fields of modern optics. Novel phenomena have been observed at nanoscales in the structures incorporating metals, with many promising applications in lasing, sensing, and waveguiding at subwavelength scales [1]. Waveguiding in plasmonic nanoguides and circuits represents a great interest for telecommunication, computing, and information processing [2]. However, in order to increase efficiency of plasmonic effects in many of such applications, one should employ nonlinear effects which may enhance substantially capabilities of plasmonic nanoguides and circuits, also leading to novel effects such as plasmon self-focusing and frequency conversion.

To study the plasmon propagation in nonlinear media, we should analyze nonlinear Maxwell’s equations for the transverse magnetic (TM) waves in the presence of a metal-dielectric interface, and a relatively small number of studies has been devoted to the analysis of the corresponding nonlinear problems. In 1980-s, several groups suggested to employ simpler approximations for the analysis of Kerr nonlinearities in a single metal-dielectric interface geometry [3, 4, 5]. More specifically, the assumption that the dielectric tensor depends either on longitudinal [3] or transverse [4] field components simplifies the problem allowing its analytical description; the validity of both these approaches has been verified by numerical simulations [5]. Later studies considered more complicated waveguiding structures, including a metal film embedded into nonlinear Kerr media [6] and nonlinear slot waveguide [7, 8].

Most of the earlier studies considered nonlinear guided waves in metal-dielectric structures localized in the transverse dimension. The temporal effects and the formation of temporal surface-polariton solitons were first discussed by Boardman et al. [9], whereas recently Feigenbaum and Orenstein [10] considered spatial localization of plasmon waves in the planar metal-dielectric waveguide geometry with a slot width considerably smaller than the wavelength. However, the surface plasmon polaritons are usually introduced for a single metal-dielectric interface [1, 2], where the analysis of nonlinear localization is more complicated, and it should be carried out by taking into account the boundary conditions. More importantly, losses in the visible frequency range cannot be neglected in the study of the soliton generation and propagation in plasmonic structures.

In this paper, we study nonlinear propagation of surface plasmon polaritons along an interface between metal and nonlinear Kerr dielectric. We describe the plasmon self-focusing and generation of a spatial plasmon-soliton taking into account losses in metal, by direct numerical finite-difference time-domain (FDTD) simulations and also theoretically, by developing an analytical theory based on the effective nonlinear Schrödinger equation with losses.

2. Geometry and theoretical analysis

We study propagation of a surface plasmon-polariton beam along a metal-nonlinear dielectric interface, as shown in Fig. 1(a). We assume that dielectric is a Kerr-type nonlinear material with the intensity-dependent dielectric permittivity. The dielectric permittivity in the structure can be written in the form: ε (x)=ε _metal+iγ, for x<0, and ε(x)=ε _nln≡ε _lin+α|E|², for x>0, where ε _metal is the real part of the dielectric permittivity of a metal, γ is the imaginary part of dielectric permittivity that describes losses in metal, ε _lin is the linear part of the dielectric permittivity of the dielectric, and α is a nonlinear coefficient.

First, we study eigenmodes in the linear case, when α=0 and losses are neglected, γ=0. Then, only TM waves with magnetic field H in the plane of the interface can propagate, and the general wave equation can be written as [11]: ∇×∇×E-ε (x)E=0, where ∇ is a three-dimensional vector operator, E=(E_x,0,E_z) is the electric field of a plasmon, and the coordinates are normalized by 1/k ₀=λ/2π, where λ is the free-space wavelength.

 

Fig. 1. (a) Schematic of a metal-dielectric interface shown with a transverse structure of the surface plasmon. (b) Geometry used in the numerical FDTD calculations.

Download Full Size | PDF

The equation is invariant in the y-direction, so that we look for its solutions in the form

E=C(x E _x0+izE _z0)e ^-iβz,

where C is the amplitude, β is the plasmon wavenumber normalized to k ₀, which satisfies the dispersion relation, β ²=ε _metal ε _lin(ε _metal +ε _lin)^-1 [12]. The field components E _x0 and E _z0 describe the transverse structure, and they are found from the boundary conditions [12]:

Next, we take into account losses, nonlinearity and also include the y-dependence describing the beam self-focusing in the paraxial approximation. In this case, the general wave equation can not be solved analytically since all three field components interplay in the nonlinear wave propagation. Hence, we make several simplifications to the problem. First, we assume that the plasmon remains mainly TM polarized, thus the E_y component of the electric field is negligible, E_y≪E_x,E_z and we do not take it into account in our analysis. This approximation is valid for paraxial plasmonic beams. Thus, the general vector equation can be rewritten on the axis projections as follows:

We search for the solution of the equation in the following form:

where E ₀=(E _x0,0,E _z0) is the linear plasmon profile satisfying Eq. (1) with the plasmon wavenumber β, A and B being slowly varying amplitudes. An other assumption we make is that the plasmon transverse structure does not change significantly, so that ∂_xA⇍0 and ∂_xB⇍0. Substituting Eq. (3) into Eq. (2) and taking into account the linear plasmon profile of Eq. (1), we obtain,

$E_{x0}\left[A_{\mathrm{yy}}^{″}+A_{\mathrm{zz}}^{″}-2i\beta {A\text{'}}_z\right]-i{B}_z^{\prime }\frac{\partial E_{z0}}{\partial x}+i\gamma A{E}_{x0}+\alpha A\left({\mid A\mid }^2\mid E_{x0}^2+{\mid B\mid }^2\mid E_{z0}^2\right)E_{x0}=0.$ $E_{z0}{B}_{\mathrm{yy}}^{″}+i{A}_z^{\prime }\frac{\partial E_{x0}}{\partial x}+i\gamma B{E}_{z0}+\alpha B\left({\mid A\mid }^2\mid E_{x0}^2+{\mid B\mid }^2\mid E_{z0}^2\right)E_{z0}=0.$.

 

Fig. 2. Magnetic field distribution at the metal-dielectric interface for: (a) Linear diffraction (low power), R ₀=0.35µm. (b) Self-focusing and soliton generation for |E|²⇍3×10⁶ V ²/µm ² (Media 1). (c) Same as in (b) but for low losses (Media 2). Animations corresponding to (b) and (c) show that soliton formation with gradual increase of the incident power.

Download Full Size | PDF

Assuming that A⇍B and the field amplitude A varies slowly along the propagation direction, i.e |A_zz|≪|2iβ A_z|, we neglect the term A_zz. Multiplying the first equation by E _x0 and the second equation by E _z0, thus we include the effect of the losses and nonlinearity on the evolution of the corresponding mode components. Taking the sum and integrating over the transverse dimension x, we derive the effective nonlinear Schrödinger equation with losses,

where D=∫E ² _x0 dx/∫|E ₀|² dx, Γ=∫γ(x)|E ₀|² dx/∫|E ₀|² dx and I=1/2∫α(x)|E ₀|⁴ dx/∫|E ₀|² dx, γ(x) is the imaginary part of dielectric constant of metal, α(x) - Kerr nonlinear coefficient. Equation (4) describes the propagation of the plasmon polariton beam localized at the interface and expanding or focusing in the transverse direction y. The similar equation but without losses was derived in Ref. [10] for a nonlinear slot waveguide. We have verified that the effective coefficient Γ describes the attenuation of surface plasmons with a good accuracy (as compared to the imaginary part of the plasmon propagation constant) when the absolute value of metal dielectric permittivity is much larger then the permittivity of dielectric layer, i.e. in our case for the wavelengths above 600nm.

In the linear regime (I=0), Eq. (4) describes plasmon diffraction with the diffraction coefficient D _diff=1/2βDR ² ₀, where R ₀ is the beam width. In the lossless (Γ=0) nonlinear (I≠0) regime, Eq. (4) has a stationary soliton solution. This solution can be presented as A=ϕ(y)e ^{-iΔβ z}, where ϕ(y) satisfies the equation, ϕ_yy-2βDΔβϕ+2Iϕ ³=0, Δβ being a nonlinear correction to the plasmon wavenumber. Finally, the spatial plasmon-polariton soliton profile is ϕ(y)=qI ^-1/2sech(qy), where $q=\sqrt{2\beta D\Delta \beta }$, and the soliton width is Δy≈1/q, which is much larger than the plasmon transverse width, Δx⇍1/κ ₂. In the lossless case, the soliton peak intensity is found as |E|²=q ²/I|E ₀|².

In the presence of losses (Γ≠0), Eq. (4) has no stationary solutions, and the soliton can be treated in a generalized sense as a nonlinear localized wave with varying parameters [13].

3. Numerical FDTD simulations

To evaluate the applicability of our analytical results, we employ numerical FDTD simulations and study the propagation of plasmonic beams along a metal-dielectric interface for different powers. For the FDTD simulations, we use the commercial software RSoft with Fullwave package implemented at the Australian National Supercomputer Facility.

The simulated structure consists of a metal layer and nonlinear dielectric having a common interface [see Figs. 1(a,b)]. We choose the nonlinear dielectric with ε _lin=4.84 and χ=1.4×10^-19 m ²/V ². Such material parameters can be found in chalcogenite glasses [14]. The metal is chosen to be silver, since it is widely used for plasmonic applications [1, 2]. The dispersion of metal is modeled by the Drude-Lorentz formula [1], ε _metal=ε_inf+∑δ_k/(akω ²+ib_kω-ω ² _0k), k=1..6 with the coefficients chosen in such a way that the metal permittivity fits the experimental data [15].

The whole simulation domain is surrounded by perfectly matching layers to minimize reflection from the edges of the simulation domain. Both the metal and dielectric layers are 2µm thick, which is much larger than the plasmon decay into either medium for the wavelength of interest. The transverse width of the layers is chosen 8µm so that the diffraction and soliton formation processes can be clearly seen without the effect of the boundaries used in the FDTD simulations, see Fig. 1(b).

To excite plasmon-polaritons at the metal-dielectric interface we use a setup shown in Fig. 1(b). We illuminate the structure with a continuous TM polarized light (H⇈y) with the free space wavelength λ=800nm. The source is placed in a close proximity to the structure, and its transverse dimensions are 0.5×0.3µm ² in order to create a relatively narrow plasmonic beam with R ₀⇍0.35µm. Plasmons are excited by the light scattered on the metal corner. To avoid exciting plane waves in the dielectric, we place a thin metal screen parallel to the front dielectric facet, 250nm above the metal surface, see Fig. 1(b).

To observe the plasmon beam propagation, we plot the field in the interface plane. Distribution of the magnetic field in this plane is shown in Figs. 2(a–c). Figure 1(a) shows the transverse plasmon structure in the linear regime, where the plasmon diffracts in the propagation direction. Due to losses in metal, the plasmon beam decays while propagating along the surface, the estimated propagation distance is about 5µm. The plasmon is highly localized near the metal-dielectric interface, the maximum decay length into the dielectric is about 15 nm. The plasmon wavelength obtained from the simulations is λ_sp=0.325µm, which agrees well with the guided index β=2.46 (note that at λ=800nm, we have ε _metal⇍-25+i1.74).

At low powers, the plasmon beam diffracts, the diffraction length is about 2.5µm. For larger powers, e.g. |E|²⇍2×10⁶ V ²/µm ², we observe the beam narrowing and self-focusing. Further increase of the beam intensity leads to the soliton generation, see Fig. 2(b); this corresponds to the maximum change of the nonlinear index α|E|² ε ^-1/2=ε ^-1/23/4χ|E|²≈0.05, being consistent with the previous nonlinear plasmonic results [10, 16]. This change of the refractive index cannot be archived in conventional dielectric materials, however semiconductor multi-quantum wells can potentially provide such a high-index change [17].

In our case, self-focusing is clearly observed at 1µm. The generated spatial soliton propagates about 700–800nm and then starts loosing its power due to dissipation, propagating for about 3µm. Since the process of the soliton formation is overshadowed by the losses, solely for the presentation purposes we demonstrate the same dynamics at lower value of losses, taking b ^k _new=b_k/100. The results of the corresponding FDTD simulations are presented in Fig. 2(c), where the generated soliton propagates practically without loosing its energy, and thus preserving its shape for longer distances, similar to spatial solitons in optics [13]. Although the spatial soliton exists for longer distances, eventually it decays due to losses, which cannot be compensated by any value of nonlinearity.

To analyze the effect of losses on the plasmon self-focusing, we trace the variation of the effective beam width, R, and present the results in Fig. 3(a), similar to the physics of magnetic solitons [18]. In the linear regime (curve 1), the plasmon beam diffracts, this is indicated by a monotonic growth of the beam width with the propagation distance. In the nonlinear regime, the beam width decreases reaching a minimum (curve 2), where self-trapping takes place and the spatial plasmon soliton is formed. For longer propagation distances, the soliton broadens due to dissipation decreasing its amplitude, but this broadening is reduced dramatically for low losses (curve 3). We also model the dissipation processes in the system by solving the Schrödinger equation (4) with the help of beam propagation method [19]. The theoretical results (dashed lines) are in a good agreement with numerical simulations. However, our theoretical model works well only for highly localized plasmons, i.e. with relatively small free-space wavelengths. For low-frequency plasmons the effective dynamics is not described well by the simplified model and the assumption A⇍B is not valid.

 

Fig. 3. (a) Evolution of the beam width for: 1 - linear regime, 2 - nonlinear regime with realistic losses, 3 - soliton formation in nonlinear regime with low losses. Dashed lines correspond to the results obtained from Eq. (4). (b) Soliton profiles: analytical (solid) and FDTD simulations (dots). (c) Soliton intensity vs. the effective index: solid - analytical results; crosses - FDTD data.

Download Full Size | PDF

To demonstrate the accuracy of our theoretical analysis we compare analytical and numerical results by calculating the soliton profiles in the stationary lossless regime. From our FDTD simulations, we find the normalized intensity χ|E|² of the generated soliton at z⇍1µm and, using the analytical expression for the soliton amplitude, we extract the nonlinear correction to the wavenumber, Δβ, and plot the corresponding soliton profile in Fig. 3(b). It is clearly seen that the theoretical (solid) and numerical (dots) results are in a good agreement. In addition, we calculate the dependence of the soliton intensity on the nonlinear effective index β _NL=β+Δβ and compare it to the values obtained in numerical simulations. For different soliton intensities we find the distance over which the phase incursion is equal to 2π and find the plasmonic wavelength and nonlinear effective index β _NL. The normalized intensity vs. effective index is presented in Fig. 3(c) demonstrating a good agreement with the FDTD simulations.

4. Conclusions

We have studied nonlinear self-focusing of surface plasmon polaritons propagating along an interface between metal and nonlinear dielectric. While for low powers, we have observed linear diffraction of plasmons, the beam propagation becomes different for larger powers with self-focusing of plasmon beams, even in the presence of strong losses in metal. These feature of the plasmon propagation in a nonlinear medium have been shown to agree with an analytical theory based on the nonlinear Schrödinger equation with losses.

The authors acknowledge a support of the Australian Research Council.