Nonlinear Surface Plasmons

# Nonlinear Surface Plasmons

Ryan G. Halabi Department of Mathematics, University of California at Davis  and  John K. Hunter Department of Mathematics, University of California at Davis
October 27, 2015
###### Abstract.

We derive an asymptotic equation for quasi-static, nonlinear surface plasmons propagating on a planar interface between isotropic media. The plasmons are nondispersive with a constant linearized frequency that is independent of their wavenumber. The spatial profile of a weakly nonlinear plasmon satisfies a nonlocal, cubically nonlinear evolution equation that couples its left-moving and right-moving Fourier components. We prove short-time existence of smooth solutions of the asymptotic equation and describe its Hamiltonian structure. We also prove global existence of weak solutions of a unidirectional reduction of the asymptotic equation. Numerical solutions show that nonlinear effects can lead to the strong spatial focusing of plasmons. Solutions of the unidirectional equation appear to remain smooth when they focus, but it is unclear whether or not focusing can lead to singularity formation in solutions of the bidirectional equation.

###### Key words and phrases:
Plasmon, surface wave, nonlinear optics, nonlocal evolution equation
###### 2010 Mathematics Subject Classification:
35Q60
The second author was partially supported by the NSF under grant number DMS-1312342.

## 1. Introduction

Surface plasmon polaritons (or SPPs) are electromagnetic surface waves that propagate on an interface between a dielectric and a conductor and decay exponentially away from the interface; for example, optical SPPs propagate on an interface between air and gold. Maradudin et. al. [9] give an overview of SPPs and their applications in plasmonics. Kauranen and Zayats [8] review nonlinear aspects of plasmonics.

We model SPPs by the classical macroscopic Maxwell equations, and consider the basic case of SPPs that propagate along a planar interface separating an isotropic dielectric in the upper half-space from an isotropic conductor in the lower half-space (see Figure 1).

Figure 2 shows a typical linearized dispersion relation for such SPPs. The phase speed of long-wavelength SPPs approaches a constant speed , and the frequency of short-wavelength SPPs approaches a constant frequency . This limiting frequency is determined by the condition that

 ^ϵ+(ω0)+^ϵ−(ω0)=0

where is the permittivity of the dielectric and is the permittivity of the conductor, expressed as functions of the frequency. In this short-wave limit, the electromagnetic field is approximately quasi-static, and we will refer to the corresponding oscillations on the interface as surface plasmons (SPs) for short.

Nonlinear SPs are interesting to study because their linearized frequency is independent of their wavenumber. As a result, they are nondispersive with zero group velocity, and weakly nonlinear SPs are subject to a large family of cubically nonlinear, four-wave resonant interactions between wavenumbers such that

 k1+k2=k3+k4.

The corresponding resonance condition for the frequencies, , is satisfied automatically.

The nonlinear self-interaction of an SP therefore leads to a wave with a narrow frequency spectrum and a wide wavenumber spectrum, meaning that nonlinearity distorts the spatial profile of the wave. This behavior is qualitatively different from NLS-type descriptions of dispersive waves in nonlinear optics, where the spatial profile consists of a slowly modulated, harmonic wavetrain [10]. Thus, the limit considered here is useful in understanding the spatial dynamics of short-wave optical pulses in the case of surface plasmons.

We remark that if the spatial dispersion of the optical media is negligible, as we assume here, then the behavior of the SP depends only on the response of the media at the frequency , so we do not need to make any specific assumptions about the frequency-dependence of their linear permittivity or nonlinear susceptibility.

A related example of constant-frequency surface waves on a vorticity discontinuity in an inviscid, incompressible fluid is analyzed in [1], which gives further discussion of constant-frequency waves.

Our asymptotic solution for the tangential electric field in the -direction of a weakly nonlinear SP on the interface has the form

 (1) E∥(x,0,t)=ε[A(x,ε2t)e−iω0t+A∗(x,ε2t)eiω0t]+O(ε3),

as with , where is complex-valued amplitude. This solution consists of an oscillation of frequency whose spatial profile evolves slowly in time.

To write the evolution equation for in a convenient form, we introduce the projections , onto positive and negative wave numbers, and define

 (2) u(x,τ)=P[A](x,τ)=∫∞0^A(k,τ)eikxdk,v(x,τ)=Q[A](x,τ)=∫0−∞^A(k,τ)eikxdk,

where

 A(x,τ)=∫∞−∞^A(k,τ)eikxdk.

For spatially periodic solutions, the Fourier integrals are replaced by Fourier series. Then

 E∥=ε[ue−iω0t+u∗eiω0t]+ε[ve−iω0t+v∗eiω0t]+O(ε3),

where is the amplitude of the right-moving, positive-frequency waves and is the amplitude of the left-moving, positive-frequency waves.

In the absence of damping and dispersion, the complex-valued amplitudes , satisfy the following cubically nonlinear, nonlocal asymptotic equations

 (3) uτ=i∂xP[αu∗∂−1x(u2)+βv∂−1x(uv∗)],P[u]=u,vτ=i∂xQ[αv∗∂−1x(v2)+βu∂−1x(u∗v)],Q[v]=v,

where the coefficients are proportional to sums of the nonlinear susceptibilities of the media on either side of the interface, and the inverse derivative is defined spectrally by

 ∂−1x(eikx)=1ikeikx.

Here, we assume that the Fourier transforms of , vanish to sufficient order at in the spatial case , or that , have zero mean in the spatially periodic case . As we show in Section 5, the spatially periodic form of (3) is an ODE in the -Sobolev space for .

The inclusion of weak short-wave damping and dispersion in the asymptotic expansion leads to a more general equation (21), which consists of (3) with additional lower-order, linear terms. The spectral form of this equation for is given in (17).

A basic feature of (3) is that the actions of the left and right moving waves, given in (28), are separately conserved. In quantum-mechanical terms, this means that the numbers of left and right moving plasmons are conserved. In billiard-ball terms, the collision of two right-moving plasmons produces two right-moving plasmons; recoil into a low-momentum, left-moving plasmon and a high-momentum, right-moving plasmon is not allowed.

This property is explained by the fact that surface plasmons have helicity [2], even though they are plane-polarized, and the helicities of left and right moving surface plasmons have opposite signs. As a result, conservation of angular momentum implies that nonlinear interactions cannot convert right-moving plasmons into left-moving plasmons, or visa-versa.

In particular, setting in (3) and normalizing , we get an equation for unidirectional surface plasmons

 (4) uτ=i∂xP[u∗∂−1x(u2)],P[u]=u.

As we discuss further in Section 6, this equation is related to the completely integrable Szegö equation,

 (5) iuτ=P[|u|2u],P[u]=u,

which was introduced by Gérard and Grellier [7] as a model for totally nondispersive evolution equations. In our context, the projection operator arises naturally from the condition that nonlinear interactions between positive-wavenumber components do not generate negative-wavenumber components.

In the rest of this paper, we describe the linearized equations for SPPs and derive the asymptotic equations for nonlinear SPs. We prove a short-time existence result for the asymptotic equation (3), describe its Hamiltonian structure, and prove the global existence of weak solutions for the unidirectional equation (4). We also present some numerical solutions, which show that nonlinear effects can lead to strong spatial focusing of the SPs. Solutions of the unidirectional equation (4) appear to remain smooth through the focusing, but it is unclear whether or not singularities form in solutions of the system (3).

## 2. Maxwell’s equations

The macroscopic Maxwell equations for electric fields , and magnetic fields , are

 ∇⋅B=0,∇×E+Bt=0, ∇⋅D=0,∇×H−Dt=0,

where we assume that there are no external charges or currents.

We consider nonmagnetic, isotropic media without spatial dispersion, in which case we have the constitutive relations

 (6) D(x,t)=∫t−∞ϵ(t−t′)E(x,t′)dt′+∫t−∞χ(t−t1,t−t2,t−t3)[E(x,t1)⋅E(x,t2)]E(x,t3)dt1dt2dt3+O(|E|5),B(x,t)=μH(x,t),

where is the permittivity, is a third-order nonlinear susceptibility, and the magnetic permeability is a constant.

We suppose that the interface between the dielectric and the conductor is located at , with

 (7) ϵ,χ={ϵ+,χ+if z>0,ϵ−,χ−if z<0,

and assume that both media have the same magnetic permeability . The jump conditions across are

 (8) [n⋅B]=0,[n×E]=0,[n⋅D]=0,[n×H]=0,

where is the unit normal, denotes the jump in across the interface, and we assume that there are no surface charges or currents. We further require that the fields approach zero as .

We consider transverse-magnetic SPPs that propagate in the -direction and decay in the -direction, with in the -plane and in the -direction. Solution of the linearized equations for Fourier-Laplace modes gives an electric field of the form [9]

 E(x,z,t)={^Φ+eikx−β+|k|z−iωtin z>0,^Φ−eikx+β−|k|z−iωtin z<0,

where are positive constants and are constant vectors in the -plane. The frequency satisfies the dispersion relation

 (9) k2=μω2[^ϵ+(ω)^ϵ−(ω)^ϵ+(ω)+^ϵ−(ω)],

where denotes the Fourier transform of .

As a typical example, the permitivities for an interface between a vacuum and a loss-less Drude metal are given by

 ^ϵ+=ϵ0,^ϵ−(ω)=ϵ0(1−ω2pω2),

where is the plasma frequency of the metal. In that case, the dispersion relation (9) becomes

 (10) ω2=ω20+c20k2−√ω40+c40k4,

where the limiting frequency and speed are given by

 ω0=ωp√2,c0=1√ϵ0μ.

This dispersion relation is plotted in Figure 2.

For example, one estimate of the plasma frequency for the Drude model of gold, cited in [11], is . The corresponding limiting frequency of an SP on a vacuum-gold interface is , which lies in the near ultraviolet. The wavenumber is given by , and the frequency of an SP is close to when , corresponding to wavelengths nm. Ultraviolet frequencies are outside the range of usual applications of SPs, but, with the inclusion of weak dispersive effects, our asymptotic solution applies to SPs with somewhat smaller frequencies. Moreover, different parameter values are relevant for other plasmonic materials than gold.

If damping effects are included in the free-electron Drude model for the metal, then one gets a complex permittivity

 ^ϵ(ω)=ϵ0(1−ω2pω2+iγω)=ϵ0(1−ω2pω2+iγω2pω3+O(γ2)).

Typical values of are of the order , corresponding to relatively weak damping of an SP over the time-scale of its period.

## 3. Asymptotic expansion

In order to carry out an asymptotic expansion for quasi-static SPs, we first non-dimensionalize Maxwell’s equations. Let , , , and be a characteristic frequency, wave speed, permittivity, and permeability for the SP, with , and let denote a characteristic electric field strength at which the response of the media is nonlinear. We define corresponding wavenumber and magnetic field parameters by

 k0=ω0c0,B0=E0c0.

The quasi-static limit applies to wavenumbers . If is a characteristic length-scale for variations of the electric field, then we assume that

 ε=k0λ≪1

is a small parameter. If is a typical magnitude of the electric field strength in the SP, then we also assume that

 δ=E∗E0≪1.

In the expansion below, we take , which leads to a balance between weak nonlinearity and weak short-wave dispersion. The nondispersive equations (3) apply in the limit .

We define dimensionless variables, written with a tilde, by

 ~x=xλ,~t=ω0t,~E=EE0,~D=Dϵ0E0,~B=BB0,~H=μ0HB0.

The corresponding dimensionless permittivity, nonlinear susceptibility, and permeability are

 ~ϵ=ϵϵ0ω0,~χ=E20χϵ0ω30,~μ=μμ0.

After dropping the tildes, we get the nondimensionalized Maxwell equations

 (11) ∇⋅B=0,∇×E+εBt=0,∇⋅D=0,∇×H−εDt=0,

with the constitutive relations (6), where and are the nondimensionalized permitivities and nonlinear susceptibilities in and , and is the nondimensionalized permeability.

The simplest form of the equations, in which we neglect dispersion and include only nonlinearity, arises when we take in (11). In that case, and , satisfy the purely electrostatic equations

 ∇×E=0,∇⋅D=0,

with the time-dependent constitutive relation in (6).

We assume that, in the frequency range of interest, the permittivity has an expansion

 (12) ϵ=ϵr+iε2ϵi,

where is the leading-order real part of and is a small imaginary part that describes weak damping of the SP. We also assume that the susceptibility is real-valued to leading order in .

We denote the Fourier transforms of , , and by , , and , where

 ^ϵr(ω)=∫ϵr(t)eiωtdt,^ϵi(ω)=∫ϵi(t)eiωtdt, ^χ(ω1,ω2,ω3)=∫χ(t1,t2,t3)eiω1t1+iω2t2+iω3t3dt1dt2dt3.

These transforms have the symmetry properties

 (13) ^ϵr(−ω)=^ϵr(ω),^ϵi(−ω)=−^ϵi(ω),^χ(−ω1,−ω2,−ω3)=^χ(ω1,ω2,ω3),^χ(ω2,ω1,ω3)=^χ(ω1,ω2,ω3).

Using the method of multiple scales, we look for a weakly nonlinear asymptotic solution for SPs of the form

 (14) E=εE1(x,z,t,τ)+ε3E3(x,z,t,τ)+O(ε5),D=εD1(x,z,t,τ)+ε3D3(x,z,t,τ)+O(ε5),B=ε2B2(x,z,t,τ)+O(ε4),H=ε2H2(x,z,t,τ)+O(ε4),

where the ‘slow’ time variable is evaluated at .

The details of the calculation are given in Section 9. We summarize the result here. At the order , we find that the leading-order electric field is given by the linearized solution

 (15) E1(x,z,t,τ)=Φ1(x,z,τ)e−iω0t+Φ∗1(x,z,τ)eiω0t,Φ1(x,z,τ)=∫⎡⎢⎣[]c10isgn(kz)⎤⎥⎦^A(k,τ)eikx−|kz|dk,

where satisfies

 ^ϵ+r(ω0)+^ϵ−r(ω0)=0,

and is a complex-valued amplitude function. In particular, the tangential -component of the electric field on the interface is given by (1) with

 (16) A(x,τ)=∫^A(k,τ)eikxdk.

Writing the solution in real form, we have

 E∥(x,0,t,τ)=R(x,τ)cos(ω0t+ϕ(x,τ)),A=12Re−iϕ.

This electric field consists of oscillations of frequency with a spatially-dependent amplitude and phase that vary slowly in time.

The imposition of a solvability condition at the order to remove secular terms from the expansion gives the following spectral equation for :

 (17) ^Aτ(k,τ)=i|k|∫T(k,k2,k3,k4)^A∗(k2,τ)^A(k3,τ)^A(k4,τ)δ(k+k2−k3−k4)dk2dk3dk4−γ^A(k,τ)+iνk2^A(k,τ).

The nonlinear term in (17) describes four-wave interactions of wavenumbers , , into , as indicated by the Dirac -function in the integral. The interaction coefficient is given by

 (18) T(k1,k2,k3,k4)=2a(k1k3+|k1k3|)(k2k4+|k2k4|)k1k2k3k4(|k1|+|k2|+|k3|+|k4|)+b(k1k2−|k1k2|)(k3k4−|k3k4|)k1k2k3k4(|k1|+|k2|+|k3|+|k4|),

where

 (19) a=^χ+(ω0,−ω0,ω0)+^χ−(ω0,−ω0,ω0)^ϵ′+r(ω0)+^ϵ′−r(ω0),b=^χ+(ω0,ω0,−ω0)+^χ−(ω0,ω0,−ω0)^ϵ′+r(ω0)+^ϵ′−r(ω0).

Here, the prime denotes a derivative with respect to . The coefficients of damping and dispersion are given by

 (20) γ=^ϵ+i(ω0)+^ϵ−i(ω0)^ϵ′+r(ω0)+^ϵ′−r(ω0).ν=12μω20[^ϵ+r(ω0)2+^ϵ−r(ω0)2^ϵ′+r(ω0)+^ϵ′−r(ω0)],

These coefficients agree with the expansion of the full linearized dispersion relation (9) as , which is

 ω=ω0−νk2−iε2γ+….

The relative strength of the two terms in is proportional to . The coefficient appearing in is associated with cubically nonlinear interactions that produce waves of the same polarization, while the coefficient appearing in is associated with interactions that produce waves of opposite polarization [3]. A typical value of the ratio

 r=^χ(ω0,ω0,−ω0)^χ(ω0,−ω0,ω0)

for nonlinearity due to nonresonant electronic response is . If are the values of this ratio in and , then

 ba=r−+sr+1+s,s=^χ+(ω0,−ω0,ω0)^χ−(ω0,−ω0,ω0).

For example, if the dielectric in is a vacuum, then and .

It is convenient to write the interaction coefficient in (17) as an asymmetric function. We could instead use the symmetrized interaction coefficient

 Ts(k1,k2,k3,k4)=12[T(k1,k2,k3,k4)+T(k1,k2,k4,k3)],

which satisfies

 Ts(k1,k2,k3,k4)=Ts(k2,k1,k3,k4)=Ts(k1,k2,k4,k3)=Ts(k3,k4,k1,k2).

The above solution generalizes in a straightforward way to two-dimensional SPs that depend on both tangential space variables , with corresponding tangential wavenumber vector . In that case, we write the equations in terms of a potential variable instead of a field variable , where . One finds that

 E1(x,z,t,τ) =Φ1(x,z,τ)e−iω0t+Φ∗1(x,z,τ)eiω0t, Φ1(x,z,τ) =∫[ik−|k|sgn(z)]^a(k,τ)eik⋅x−|kz|dk,

where satisfies

 ^aτ(k,τ) =i|k|∫S(k,k2,k3,k4)^a∗(k2,τ)^a(k3,τ)^a(k4,τ) δ(k+k2−k3−k4)dk2dk3dk4−γ^a(k,τ)+iν|k|2^a(k,τ),

with

 S(k1,k2,k3,k4)

In this paper, we restrict our attention to one-dimensional SPs.

## 4. Spatial form of the equation

In this section, we write the spectral equation (17)–(18) in spatial form for given in (16). To simplify the notation, we do not show the time-dependence of .

First, consider a nonlinear term proportional to the one in with coefficient :

 ^F(k1) =∫(k1k4+|k1k4|)(k2k3+|k2k3|)2k1k2k3k4(|k1|+|k2|+|k3|+|k4|) δ12−34^A∗(k2)^A(k3)^A(k4)dk2dk3dk4,

where we write

 δ(k1+k2−k3−k4)=δ12−34.

The interaction coefficient in this integral is nonzero only if , have the same sign and , have the same sign. Furthermore, in that case, we have

 (k1k4+|k1k4|)(k2k3+|k2k3|))2k1k2k3k4(|k1|+|k2|+|k3|+|k4|)=2|k1|+|k2|+|k3|+|k4| =⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩1/(k3+k4)if k1,k4>0 and k2,k3>0,1/(k2−k4)if k1,k4<0 and k2,k3>0,−1/(k2−k4)if k1,k4>0 and k2,k3<0,−1/(k3+k4)if k1,k4<0 and k2,k3<0,

on . We decompose into its positive and negative wavenumber components,

 ^A(k) =^u(k)+^v(k), ^u(k) =^P[A](k)={^A(k)if k>0,0if k<0, ^v(k) =^Q[A](k)={0if k>0,^A(k)if k<0.

The corresponding spatial decomposition is , where , are given by (2).

It then follows that

 ^F(k1) =^P∫δ12−34^u∗(k2)^u(k3)^u(k4)k3+k4dk2dk3dk4 +^Q∫δ12−34^u∗(k2)^u(k3)^v(k4)k2−k4dk2dk3dk4 −^P∫δ12−34^v∗(k2)^v(k3)^u(k4)k2−k4dk2dk3dk4 −^Q∫δ12−34^v∗(k2)^v(k3)^v(k4)k3+k4dk2dk3dk4.

The integral in the first term,

 ^I(k1)=∫δ12−34^u∗(k2)^u(k3)^u(k4)k3+k4dk2dk3dk4,

has inverse Fourier transform

 I(x) =∫^I(k1)eik1xdk1 =∫^u∗(k2)e−ik2x^u(k3)^u(k4)ei(k3+k4)xk3+k4dk2dk3dk4 =iu∗∂−1x(u2).

The other terms in are treated similarly, and we find that

 F(x)=iP[u∗∂−1x(u2)+v∂−1x(uv∗)]−iQ[u∂−1x(u∗v)+v∗∂−1x(v2)].

The nonlinear term proportional to ,

 ^G(k1) =∫(k1k2−|k1k2|)(k3k4−|k3k4|)2k1k2k3k4(|k1|+|k2|+|k3|+|k4|) δ12−34^A∗(k2)^A(k3)^A(k4)dk2dk3dk4,

is zero unless , and , have opposite signs, in which case

 (k1k2−|k1k2|)(k3k4−|k3k4|)2k1k2k3k4(|k1|+|k2|+|k3|+|k4|) =⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩1/(k2−k3)if k1<0, k2>0, k3<0, k4>0,1/(k2−k4)if k1<0, k2>0, k3>0, k4<0,−1/(k2−k4)if k1>0, k2<0, k3<0, k4>0,−1/(k2−k3)if k1>0, k2<0, k3>0, k4<0,

on . It then follows that

 G(x)=2iP[v∂−1x(uv∗)]−2iQ[u∂−1x(u∗v)].

Projecting (17) onto positive and negative wavenumbers and using the previous equations to take the inverse Fourier transform, we find that

 (21) uτ=i∂xP[αu∗∂−1x(u2)+βv∂−1x(uv∗)]−γu−iν∂−2xu,P[u]=u,vτ=i∂xQ[αv∗∂−1x(v2)+βu∂−1x(u∗v)]−γv−iν∂−2xv,Q[v]=v,

where , are given by (20) and , are given by

 (22) α=4a,β=4(a+b),

where , are defined in (19).

In the presence of damping, SPs require external forcing to maintain their energy. Direct forcing by electromagnetic radiation is not feasible, since the wavenumber, or momentum, of an electromagnetic wave with same frequency as the SP is smaller than that of the SP. Instead, SPs are typically forced by an evanescent wave (in the Otto or Kretschmann configuration). At least heuristically, this forcing can be modeled by the inclusion of nonhomogeneous terms in (21). For example, a steady pattern of forcing at a slightly detuned frequency leads to equations of the form

 uτ=i∂xP[αu∗∂−1x(u2)+βv∂−1x(uv∗)]−γu−iν∂−2xu+f(x)e−iΩτ, vτ=i∂xQ[αv∗∂−1x(v2)+βu∂−1x(u∗v)]−γv−iν∂−2xv+g(x)e−iΩτ,

where , . We do not study the behavior of the resulting damped, forced system in this paper.

## 5. Short-time existence of smooth solutions

For definiteness, we consider spatially periodic solutions of the asymptotic equations (21) with zero mean. For , let denote the Sobolev space of -periodic, zero-mean functions (or distributions)

 f(x)=∑n∈Z∗^f(n)einx,

where , such that

 ∥f∥s=∥|