Nonlinear plasmonics in a two-dimensional plasma layer

The nonlinear electron dynamics in a two-dimensional (2D) plasma layer are investigated theoretically and numerically. In contrast to the Langmuir oscillations in a three-dimensional (3D) plasma, a well-known feature of the 2D system is the square root dependence of the frequency on the wavenumber, which leads to unique dispersive properties of 2D plasmons. It is found that for large amplitude plasmonic waves there is a nonlinear frequency upshift similar to that of periodic gravity waves (Stokes waves). The periodic wave train is subject to a modulational instability, leading to sidebands growing exponentially in time. Numerical simulations show the breakup of a 2D wave train into localized wave packets and later into wave turbulence with immersed large amplitude solitary spikes. The results are applied to systems involving massless Dirac fermions in graphene as well as to sheets of electrons on liquid helium.


Introduction
Two-dimensional (2D) conducting layers have unique plasmonic properties, where the wave frequency of lowenergy plasmons (in the long wavelength limit) typically depends on the square root of the wavenumber, and observed 2D plasmons have wave-frequencies most commonly in the THz range or below. Examples of experimental observations include 2D plasmons in a sheet of electrons on liquid helium [1], in silicon inversion layers [2], monolayer graphite [3], single layer graphene on SiC substrates [4,5], epitaxial graphene [6], etc. Graphene has unique electronic and material [7][8][9][10] properties, and plasmons in graphene is now a field of intense research with important applications in optics, electronics, metamaterials, light harvesting, energy storage, THz technology, and so on [11,12]. Experiments have shown that Dirac plasmons in graphene on SiC are hybridized with optical phonons of the SiC substrate [13], while scattering of plasmons off capillary waves (riplons) were observed in an electron sheet on liquid helium [1]. Nonlinear photonics have important applications in graphene with respect to harmonic generation [14], optical and plasmonic bistability [15,16], broadband optical limiting application [17], four-wave mixing frequency conversion [18], and other nonlinear photonic applications [19]. On the other hand, the self-interaction among large amplitude 2D plasmons and their sidebands could be an important source of instability and localization of wave energy. For example, the instability of large amplitude plasma waves have been suggested to lead to 'rogue waves' in semiconductor plasmas [20], and optical rogue waves have been observed in semiconductor lasers [21] and other optical systems [22,23]. It has also been proposed that graphene supports optical solitons in a graphene monolayer [24], dissipative plasmon solitons in multilayer graphene [25] and graphene nanodisk arrays [26].
The aim of this paper is to theoretically and numerically investigate the nonlinear behavior of largeamplitude plasmons in a 2D plasma layer. The basic fluid equations are presented and motivated in section 2. As discussed in section 3, the unique dispersive properties of 2D plasmons lead to a similar behavior to surface gravity waves [27][28][29], such as a nonlinear upshift of the frequency due to harmonic generation, and a modulational instability resulting in the growth of sidebands of the periodic wavetrain. In section 4, the theoretical results are compared with a direct numerical simulation of the dynamic system, which also demonstrates the breakup of the periodic wavetrain into localized envelope solitary waves and later into spiky large amplitude structures immersed in wave turbulence. Numerical examples are given for systems involving massless Dirac fermions in graphene and for sheets of electrons on liquid helium. Finally, some conclusions are drawn in section 4.

Nonlinear fluid model of a two-dimensional plasma layer
Low-energy plasmons, i.e. in the long wavelength and low-frequency limit, may be studied using a hydrodynamic model [11,30], in which the conservation of electrons (charges) is governed by a continuity equation and their fluid velocity (or alternatively the current density) by a momentum equation. Due to the dynamics of the electrons in a 2D sheet of plasma, the electrostatic potential, which extends also to the out-ofplane dimension, has a particular form which leads to a unique dispersive behavior of the system. The potential is obtained from Poisson's equation with the electron density concentrated to a 2D plane, and is given by where the integral is taken over the 2D plane of the plasma layer. Here n e is the areal electron number density, n 0 represents the neutralizing ion background density, e is the magnitude of the electron charge,  0 is the electric vacuum permittivity, and ò is the mean dielectric constant of the surrounding medium [11]. We choose a coordinate system such that the plasma layer is in the x-y plane. For a plane wave, where the electron density given by = + n n n e e 0 complex conjugate, and the potential is the magnitude of the wave vector. The factor ( ) k 1 2 in equation (2) leads to significantly different behavior of 2D plasmons compared to the three-dimensional (3D) case, where the corresponding factor k 1 2 gives rise to Langmuir oscillations. The electron dynamics is governed by the continuity and momentum equation x y e e x e y is the electron fluid velocity, x  and y  are unit vectors in the x and y-directions, D is the Drude weight [11,31], n 0 is the mean electron density, and the gradient operator , while for massless Dirac fermions, the Drude weight is 1 is the Fermi speed, and p = k n F 0 is the Fermi wavenumber. For analytic convenience in the following section, we assume irrotational potential flow, y = v e e , where y e is the flow potential, which inserted into equations (3) and (4) gives the two scalar equations, and y y f , 6 e e 2 0 0 for the unknowns n e and y e . Equations (1), (5) and (6) are coupled systems for f, n e and y e . In equation (6), B(t) is an arbitrary function of time, which we can choose to set to zero, since the fluid velocity v e is the spatial gradient of ψ and will not depend on B(t), and hence, B(t) will not influence the dynamics of n e . The nonlinearity in equation (6) can also give rise to a mean component of y e depending only on time, which likewise will not contribute to the dynamics of the system.

Nonlinear wave-wave interactions
Two major nonlinear effects of large amplitude plasmons are here investigated theoretically: a nonlinear frequency upshift due to harmonic generation, and a modulational instability of the wavetrain leading to sidebands growing exponentially in time during the linear phase. These effects have been investigated in the past with respect to general nonlinear dispersive media [32,33], plasma waves [34][35][36][37][38][39], and nonlinear water waves [28]. A full investigation of the nonlinear system would involve the following steps: i) to find the profile of the nonlinear periodic wavetrain governed by equations (5)- (6) and the corresponding nonlinear frequency (or phase velocity) depending on the wavelength and amplitude. This wavetrain would form a 'lattice' on which small-amplitude, unstable sidebands in the form of Bloch waves q q ( ) ( ) F i exp s 0 , can develop, where q 0 and q s are carrier and perturbation phases (see definitions below), and q ( ) F 0 is periodic with respect to its argument. ii) To derive a linearized system of differential equations for the Bloch functions q ( ) F 0 and corresponding perturbation frequencies Ω. The set of Bloch functions for each of the perturbation wave vector would then work as eigenfunctions with corresponding eigenvalues Ω, whose imaginary part gives the growth rate. While such a direct approach could be implemented numerically, we here follow a somewhat simpler route where the carrier wave and the Bloch functions are expanded up to the second harmonics, and where the results are more easy to interpret. We mention that in a previous study of the parametric instabilities of electron oscillations due to intense laser light [40], higher order Fourier expansions of the sidebands were employed using Floquet theory (analogous to Bloch theory).

Stokes waves
Stokes' investigation of water waves (see e.g. [27]) was the starting point of nonlinear theory of dispersive waves, where he among other results found that the dispersion relation involves the wave amplitude. For large amplitude periodic waves, the nonlinearities lead to harmonic generation, where the second harmonics are generated by products of the first harmonics as is the phase, k 0 is the wave vector, and ω is the frequency which may depend on the wave amplitude. On the other hand, the beating of the second harmonics by the first harmonics gives . The crucial point is that second harmonic components multiplied by first harmonic components give again components proportional to q ( ) sin 0 and q ( ) cos 0 . Such 'secular terms,' which seemed to drive the linear system resonantly, led Stokes to the conclusion that the frequency must change due to nonlinearity. Using a Stokes expansion including first and second harmonics . Third and higher harmonics were neglected in obtaining equations (7)-(10). Equations (7)-(10) constitute a non-linear eigenvalue problem forn e01 ,n e02 , ŷ e01 and ŷ e02 , with non-trivial solutions only if the nonlinear dispersion relation is fulfilled, where the nonlinearity is given by The linear dispersion relation w = g k e 0 0 is of the same form as for deep water gravity waves if g e is replaced by the gravitational constant g [27][28][29]. In deriving equation (11), we have used that equations (8)-(10) give the relations which inserted into equation (7) gives equation (11). Equation (11) relates the frequency ω to the wave amplitude ŷ e01 . To the lowest order, the frequency depends on the wave amplitude as which shows a nonlinear upshift of the frequency due to the finite wave amplitude.

Modulational instability
The stability of the system (5)-(6) is investigated by assuming an expansion of the carrier wave and its sidebands as illustrated in figure 1. The sidebands at w W   ( ) K k , 0 are almost resonant, i.e. they are located almost on the dispersion surface of linear waves. The coupling between these almost resonant sidebands by beating of the non-resonant sidebands at (W K , ) and 2 0 with the large amplitude wave train at w ( ) k , 0 and its second harmonic w ( ) k 2 , 2 0 forms a feedback loop which gives rise to a modulational instability for certain values of K. The couplings to second harmonics are crucial to accurately model the modulational instability, and it should be noted that Benjamin and Feir [28] in a similar manner as in figure 1 included the second harmonics in their treatment of the instability of nonlinear water waves. The nonlinear frequency shift is also important, where the carrier frequency ω is given by the nonlinear dispersion relation (11). The dynamical variables in equations (5) and (6) The respective Fourier component of the potential is obtained by using equation (2). It is assumed that the sideband components  n e , y e ,  n e2 , y e2 ,n s , and ŷ s have much smaller amplitudes than n 0 , n e01 , ŷ e01 ,n e02 , and ŷ e02 , and thus products among sideband components are neglected. The nonlinear wave 2 0 , and a non-resonant driven mode at W ( ) K , .
Since the dynamic variables are real-valued, the carrier frequency and its second harmonic also occur at w -- are also produced but are numerically found to give only minor contributions to the instability and are therefore neglected here. Separating different Fourier modes gives equations (A1)-(A10) in appendix. The coefficient matrix is on the form of a linear eigenvalue problem for the unknowns ( + n e , y +e , -n e , y -e ,n es , ŷ es , + n e2 , y +e2 , -n e2 , y -e2 ) where the complex-valued perturbation frequency W = W + G i R takes the role of an eigenvalue and where its imaginary part Γ is the growth-rate of the instability.
In the numerical solution of the eigenvalue problem, it is assumed (with no loss of generality) that ŷ e01 is real-valued, which leads to that alson e01 ,n e02 and ŷ e02 are also real-valued and are related to ŷ e01 according to equations (14)- (16); any complex phase factor can be eliminated by a simple change of sideband variables. The wave frequency ω is obtained by solving equation (11) iteratively. A coordinate system is chosen such that the carrier wave propagates in the x-direction, = k k x 0 0 . Amongst the roots given by the system in appendix, we here concentrate on the instability near resonance, where the sidebands (w  ,  k ) approximately fulfill the linear dispersion relation, so that w w W » -» - . The obtained growth-rate for the instability is shown in figure 2 1 8 x y x y The resonance curve is indicated with dashed lines in figure 2. Similar conditions as equation (18) have also been obtained in the theory of nonlinear water waves [41][42][43]; see for example equation (28) of [43], which in the limit of infinite depth  ¥ h reduces to our equation (18).

Nonlinear evolution of the system
The nonlinear evolution of the system is investigated by solving the system is resolved on a numerical grid with 250 intervals in each direction, and with periodic boundary conditions. A pseudospectral method is employed to accurately calculate spatial derivatives, and to obtain the potential by using the formula (2) in Fourier space. The standard 4th-order Runge-Kutta method is used to advance the solution in time with the time-step w D = t 0.2 0 1 . Density waves propagating from left to right are excited resonantly to an amplitude about 10% of the background density, seen in figure 3(a). Small amplitude random numbers of the order 10 −5 of the background density are added to the electron density to seed the instability. After a linear growth-phase, the wavetrain is amplitude modulated and breaks up into a train of obliquely propagating solitary waves, seen in figure 3(b), and later into more chaotic behavior (figure 3(c)), and into wave turbulence on different lengthscales ( figure 3(d)).  . The theoretical and numeric growth-rates show excellent agreement for small wavenumbers. Due to beating with the harmonics of the carrier wave, (  · e k r i 0 ,  · e k r 2i 0 , etc), the sidebands appear in a quasiperiodic pattern in figure 4(b). Figure 5 shows a time-series of the electron density n e at = ( ) ( ) x y , 0,0. The main features are that the initially periodic wavetrain (see figures 5(a), (b)) becomes unstable and breaks up into a train of localized modulated pulses ( figure 5(c)) and density spikes immersed in wave turbulence ( figure 5(d)). These spikes have amplitudes several times larger than the mean amplitude of the surrounding waves, which is similar to rogue  waves on the ocean [29]. Figure 6 shows frequency spectra in time, averaged over all points in the x-and ydirections. The spectrum for earlier times t = 0-w -5000 0 1 (see figure 6(a)) shows signatures of harmonic generation and smaller amplitude sidebands of the carrier wave. A small upshift of the carrier frequency compared to the linear carrier frequency w w = 0 is also visible in the inset of figure 6(a), which is consistent with the predicted frequency shift due to second harmonic generation given by equation (17) as indicated with a vertical dashed line in the inset. In the fully nonlinear regime shown in figure 6(b) for t = 5000-w -10000 0 1 , the frequency spectrum has evolved into a continuous spectrum with significant components at both higher and lower frequencies compared to w w = 0 .
In figure 7, the spatial wave spectrum as a function of (k x , k y ) is shown at different times. Due to the harmonic generation at multiples of =  k k 1 x 0 and k y = 0, the wave spectrum initially ( figure 7(a)) attains a quasiperiodic pattern along the propagation direction. As shown in figures 7(b), (c), the sidebands later grow to large amplitudes and the spectrum widens to both higher and lower wavenumbers. Finally, in figure 7(d), the spectrum is continuous and broad in particular at large wavenumbers.  The regions of instability in figure 2 are to lowest order similar to the ones of the Benjamin-Feir instability of large amplitude water waves modeled by a nonlinear Schrödinger equation; see e.g. figure 12 of [29] and models taking into account higher order nonlinearity and dispersion [43,44]. Based on the growth-rate of the instability as shown in figure 2, the dynamics may in the limit of long wavelengths/small wavenumbers be phenomenologically modeled by a nonlinear Schrödinger equation of the form . The nonlinear Schrödinger equation (19) takes only into account the lowest order dispersive effects, and an assumed cubic nonlinearity amplitude with a coefficient of nonlinearity α. For a monochromatic wave with amplitude = | | A A 0 , the nonlinear Schrödinger equation (19) supports a modulational instability with the growth-rate x y x y 0 0 2 0 2 0 2 2 2 2 2 2 , which has a maximum of aw x y 2 2 0 2 0 2 . By studying instability diagrams such as figure 2 and recording the maximum growth-rate as a function of amplitude, it is found that the coefficient of nonlinearity α is weakly dependent on nonlinearity and decreases with A 0 as a » + ( ) A 0.53 1 6.5 0 for , giving a » 0.5 for A 0 = 0.01 and a » 0.32 for A 0 = 0.1. It should be mentioned that formal methods of the derivation of the nonlinear Schrödinger equation for general nonlinear systems of equations can be found in [32,45,46], for nonlinear water waves in [43,44], and for plasma waves in [33,[35][36][37][38][39]47]; we postpone such a formal derivation of the nonlinear Schrödinger equation (19) to a future work.
We next discuss specific examples of experimental parameters in systems of nonlinear 2D plasmons. As a first example, we use parameters for massless Dirac Fermions in graphene [11] with an areal number density of electrons =   12 . Similar electric field amplitudes (~-10 V cm 3 1 ) have been found theoretically for dissipative plasmon solitons in multilayer graphene [25] and graphene nanodisk arrays [26]. As a second example, we consider a plasma layer consisting of an electron sheet on liquid helium [1] with an electron density = n 10 cm . Typical experimental plasmas involving sheets of electrons on liquid helium have sizes of a few centimeters, which is probably too small to observe the modulational instability in a single pass. However, the modulational instability could potentially develop in the standing wave generated by repeated reflections of the wave against the edges of the plasma.

Conclusions and discussion
In conclusion, we have investigated the nonlinear behavior of plasmons in a 2D plasma layer. For a periodic large-amplitude wave, the harmonic generation leads to a nonlinear upshift of the frequency. The large amplitude wave-train is subject to a modulational instability that leads to the growth of sidebands and modulation of the periodic wavetrain. The linear dispersion of the 2D plasmons, as well as their nonlinear frequency upshift and modulational instability, have an interesting analogy in gravity waves such as ocean waves. Direct numerical simulations of the dynamical system confirm both the predicted nonlinear upshift of the frequency and the growth-rate of the modulational instability. The wavetrain eventually breaks up into a series of modulated pulses and more complicated wave turbulence. The spatial and temporal wave spectra show a dual cascade of wave energy to both higher and lower frequencies associated with components at both smaller and larger length-scales.
It has been suggested in the past [34] that large amplitude electron plasma waves in 3D cold electron gaseous plasmas with fixed ions also can undergo a parametric instability with a redistribution of energy from large scales to small-scale sidebands. We carried out a simulation (not shown) similar to that in figure 3 for this case but could not see any sign of instability. Hence, it seems that the 2D geometry is crucial for the modulational instability of plasmons involving only electrons. In a 3D electron-ion plasma, due to the very low group velocity of the electron oscillations, the coupling to the ion dynamics instead leads to Langmuir wave collapse and strong turbulence [48,49] via the ponderomotive pressure acting on the electrons on a slow time-scale.
Our results show that a modulational instability for plasmons involving massless Dirac fermions in graphene with 10% density variations compared to the background density could give rise to a modulational instability on a typical typical time-scale of the order of a nanosecond and a length-scale of the order 10 −3 cm. The instability would give rise to large amplitude spikes emersed in wave turbulence, similar as has been observed in optical systems [22,23]. With the advances in graphene fabrication [50], such experiments could potentially be realized in near future.

Acknowledgments
BE acknowledges the hospitality of University of Macau, Macau, PR China, where this work was carried out. This work was partially supported by the EPSRC (UK), grant EP/M009386/1. Discussions with V K Tripathi are gratefully appreciated. The simulation data associated with this paper are available at: http://dx.doi.org/ 10.15129/fb2ab071-2ed2-4718-aafd-c56315764b97.