Graphene Plasmonic Waveguides for Mid-Infrared Supercontinuum Generation on a Chip

: Using perturbation expansion of Maxwell equations with the nonlinear boundary condition, a generic propagation equation is derived to describe nonlinear effects, including spectral broadening of pulses, in graphene surface plasmon (GSP) waveguides. A considerable spectral broadening of an initial 100 fs pulse with 0 . 5 mW peak power in a 25 nm wide and 150 nm long waveguide is demonstrated. The generated supercontinuum covers the spectral range from 6 µ m to 13 µ m.


Introduction
Surface plasmon polaritons offer the unique platform to guide and manipulate light at a deeply sub-wavelength scale [1][2][3].Such tight confinement of light not only enables design of compact plasmonic circuits [3], but also boosts local field intensities and thus allows optical nonlinear processes, such as third harmonic generation [4] and nonlinear switching [5], to occur over short (typically sub-millimeter) propagation distances.However, high ohmic losses in metals impose serious limitations as to nonlinear functionality of plasmonic waveguides.In particular, generation of broad optical spectra (supercontinuum) in a compact plasmonic waveguide remains to be the challenging task.
To describe nonlinear effects, such as frequency mixing and spectral broadening, in graphene waveguides, one needs to take into proper account the purely two-dimensional nature of the optical response of graphene.Recently we developed a perturbation expansion procedure of Maxwell equations, in which graphene is treated as the nonlinear boundary condition [29].This method has been applied for analysis of self-focusing and switching of monochromatic GSP in single-and bi-layer graphene structures [29,30], and for the description of nonlinear pulse propagation in a graphene-clad dielectric fiber [31].In this work, the procedure is generalized to the case of a GSP waveguide.The generic pulse propagation equation and the expression for the effective nonlinear coefficient due to graphene are obtained, which enable to analyze nonlinear pulse dynamics in different GSP geometries.As an example, in this work a GSP waveguide formed by a modulation of the graphene Fermi level [17] is considered.We find that in such waveguides the characteristic nonlinear length can be several orders of magnitudes smaller than the propagation length dictated by damping.As the result, one can generate broadband spectra covering more than one frequency octave (6-13 µm) from spectral broadening of a 100 fs pulse with only a few hundred micro-Watt peak power and over propagation distances as short as few hundred nano-metres.

Derivation of Pulse Propagation Equation
Consider a GSP waveguide, in which graphene (finite width ribbon or a modulated layer) is located in plane x = 0, being embedded into a dielectric structure.The waveguide is homogeneous along the propagation direction z, see Figure 1.To describe nonlinear pulse propagation in the structure, it is convenient to use Fourier expansion of the total electric field: and similar expansions for other fields.Each Fourier component E solves Maxwell equations: Optical response of all bulk materials (dielectrics) is incorporated in the displacement vector D. Atom-thick graphene layer is described by means of the surface current J, the corresponding boundary condition is: n where n is the unit vector normal to the graphene layer and pointing from medium 1 to medium 2, which are on either side of the graphene layer.To simplify the derivations, in this work we focus on the cubic (Kerr) nonlinear response of graphene, leaving the polarization due to the dielectrics to be purely linear.It is a straightforward exercise to extend the described below procedure and include nonlinear dielectric polarization.At the same time, for GSPs in planar structures the relative contribution of dielectrics to the overall nonlinearity has been demonstrated to be minuscule [29], which justifies the simplification adopted here.We also neglect any effects associated with third harmonic generation, since they require specially engineered phase matching.In other words, we assume that the pulse spectrum is narrow enough to not accommodate the third harmonic.For the same reason, possible second harmonic generation and associated effects due to a non-zero second order nonlinear graphene conductivity, arising from spatial dispersion near the plasmon resonance [21], are also neglected.Under these assumptions, the Fourier components of the displacement vector and the graphene current are: In the above expression for J nl vertical dots stand for tensor product: a = Ô. ..bcd, For excitation frequencies below the single-photon absorption resonance in graphene, hω < 2E F , the structure of the third-order nonlinear conductivity tensor can be approximated as [19,32]: δ ij is the Kronecker's delta, the indexes i, p, j, s are each from the subset of in-plane coordinates (y, z).Also, the 2D symmetry of graphene and the assumption of zero transversal current J x = 0 still permit six additional non-zero tenzor components In absence of external magnetic fields, linear conductivity tensor has only two diagonal non-zero components: σ(1) ii = σ 1 , i = y, z.The relative dielectric permittivity = (x, y), together with the y-dependent tensors of linear and nonlinear graphene conductivity, define the transverse geometry of the waveguide, cf. Figure 1.

Perturbation Expansion of Nonlinear Maxwell Equations
It is convenient to decompose the surface current as J = J 0 + J p , so that the solution of Maxwell equations with J 0 gives the linear guided mode of the structure, while J p is treated as a perturbation and contains the nonlinear term J nl and damping due to non-zero real part of σ 1 (one can formally include the damping term into the leading order J 0 and consider quasi-guided modes with complex propagation constants, see e.g., Reference [33], however it is more convenient to work with real guided modes).
Developing perturbation expansion, we introduce a dummy small parameter s assuming J p ∼ s 3/2 .Each Fourier component of the electric field is expanded in the perturbation series as: and a similar expansion for the magnetic field is assumed.Here e ω is the mode of a given structure, β = β(ω) is the corresponding propagation constant, N ω is an optional normalization factor.We consider a single-mode situation, neglecting any nonlinear coupling to possible higher order modes of the waveguide, which is inefficient due to the high mismatch between propagation constants of different modes.Following substitution of the ansatz in Equation ( 8) into Maxwell equations, in the lowest order of the small parameter, O(s 1/2 ) the eigenvalue problem is obtained: where q 2 = β 2 − k 2 .The boundary conditions for the mode e ω are: where σ is the imaginary part of linear graphene conductivity, the operator ∆ is defined as: Solving the eigenvalue problem in Equation ( 9), we obtain the modal profile e(x, y) and the dispersion of the mode β(ω).
It is convenient to choose the normalization factor N ω as: where êz is the unit vector along z-axis.With this normalization, in the lowest order of the small parameter, O(s), the total energy carried by a pulse along the waveguide is given by: Here we assume that the pulse is localized in time, and the mode e ω (x, y) is localized in the transversal plane at all frequencies, so that the integrals in the above expression converge.
In the next order of the perturbation expansion of Maxwell equations, O(s 3/2 ), we obtain: together with the boundary conditions: where the perturbation current J p combines damping due to non-zero real part σ (R) 1 = Re(σ 1 ) and the nonlinear current: Next, we project Equation ( 17) onto the mode e ω to obtain: where a|b It is important to note that e ω and E 1 satisfy different boundary conditions, Equations (11)(12)(13) and (20)(21)(22), respectively.This removes the self-adjoint property of the operator L, so that e ω | L|E 1 = E 1 | L|e ω * .To proceed, we split integration in x in the l.h.s. of Equation (25) as and take the resulting integrals by parts.Following some algebra, it is possible to show that [29]: J py e * ω,y + J pz e * ω,z dy Computing integrals in the r.h.s. of Equation ( 25) we obtain: Taking into account that Le ω = 0, and combining the results in Equations ( 26) and ( 27), the propagation equation for the modal amplitudes A ω is obtained: where the damping and nonlinear coefficients are given by: In the above expressions the deformed scalar products are introduced: (ab) η = ηa x b x + a y b y + a z b z , (ab) η0 = a y b y + a z b z .The deformation factor η = σ 3 /σ 3 characterizes the relative impact of the orthogonal field component on the induced current in the graphene layer.
In the derivation of Equation ( 28) the following identity was used: which can be obtained directly from Equation ( 15) by using the relationship between magnetic h ω and electric e ω fields of the guided mode.

Pulse Propagation Equation
The propagation Equation (28) takes into full account material and geometrical dispersion of linear and nonlinear coefficients.However, numerical propagation within this model is a computationally challenging task.For certain applications, including propagation of a relatively narrow-band pulse (with ∆ω/ω 0 1, where ω 0 is the central frequency of the pulse), it is possible to neglect the dispersion of nonlinear coefficients and replace γ ωω 1 ω 2 ω 3 in Equation ( 28) with a constant coefficient evaluated at the reference frequency (A more accurate procedure, which still allows the preservation of the dispersion of nonlinearity, is based on the factorization of the nonlinear coefficients and is described in detai in Reference [34]): Furthermore, assuming that the spectrum of the pulse is localized in a vicinity of ω 0 , we can consider the analytical extension of the spectrum A ω onto the negative frequency half-axis, and the corresponding complex time-domain function defined as: Let us introduce the pulse envelope function Ψ(z, t): where A δ = A ω−ω 0 is the shifted spectrum of the pulse, β 0 = β(ω 0 ) and v −1 g = ∂β/∂ω(ω 0 ) are the propagation constant and the inverse group velocity at the reference frequency, respectively.It is easy to see, that the relationship between Ψ(z, t) and A(z, t) is: Extending integration in the r.h.s. of Equation ( 28) onto negative frequency domain, we thus can re-write the propagation equation in terms of the spectral amplitudes Ψ δ of the envelope function Ψ(z, t) as follows: where β δ = β − β 0 − δ/v g and F {f (t)} δ is the Fourier amplitude of the function f (t) at frequency δ.It is easy to see that: The derived Equation ( 36) is the spectral representation of the renowned generalized Nonlinear Schrödinger Equation, widely used for the description of pulse propagation in nonlinear waveguides and fibres [35].Unlike the Equation ( 28), the reduced model in Equation ( 36) is easy to integrate numerically with the help of the well-developed split-step method, whereby the nonlinear part is evaluated in time domain [35].Note that the definition of the pulse envelope function Ψ in Equation (34) does not involve any additional renormalizations.In particular, the pulse energy W , Equation (16), is computed as: The derived propagation Equation ( 28) and its simplified version in Equation ( 36), together with the set of coefficients defined in Equations ( 29) and (30) represent a comprehensive set of tools to analyze nonlinear pulse dynamics and frequency conversion processes in any GSP waveguide geometry.

Spectral Broadening in a GSP Waveguide
In this section we consider a particular example of a GSP waveguide created by a local modulation of the graphene Fermi level [17], see Figure 1.The graphene is embedded into a dielectric matrix, for which the constant permittivity = 10 was assumed.
The linear conductivity of graphene is given by [36]: For frequencies below the optical phonon resonance hω < 0.2eV (λ > 6 µm) and below the absorption threshold hω < 2E F the relaxation time due to impurities is typically below 1 ps [6,7].In the simulations we used τ = 0.5 ps and T = 300 K.
The profile of Fermi energy, shown in the top panel of Figure 1, was modeled as: Taking E F 1 = 0.25 eV and E F 2 = 0.2 eV, this modulation corresponds to the variation of conductivity from σ (I) 1,hi = 1.291 × 10 −4 S to σ (I) 1,lo = 0.973 × 10 −4 S at λ 0 =9 µm.While the propagation constant of GSP on a homogeneous graphene is inversely proportional to σ 1 [6,7]: such a local decrease of the conductivity effectively forms a waveguide for plasmons [17].
In Figure 2 linear and nonlinear parameters for waveguides of different widths w are summarized.The effective index of the mode n ef f = β/k 0 , Figure 2a, and the mode profile e ω were computed with the help of Comsol finite-element (FEM) Maxwell solver, where graphene was introduced as the surface current.This data was used to compute the attenuation length L p = 1/α ω and the nonlinear coefficient γ according to Equations ( 29) and (30), see Figure 2b,c.For the nonlinear conductivity of graphene the expression derived by Mikhailov and Ziegler [20] was used: and the orthogonal nonlinear factor η was set to zero.
To check the accuracy of the perturbation expansion procedure, the computed attenuation lengths were compared against numerical data from Comsol (where the non-zero real part of σ 1 was introduced), see circles in Figure 2b.The agreement is nearly perfect, a slight discrepancy towards large wavelengths is due to the finite thickness of dielectric layers used in Comsol, which starts to play a role as the mode becomes less localized.The attenuation length of a localized GSP in a waveguide is generally larger than that of a GSP on a homogeneous graphene sheet, and it grows as the waveguide width w decreases [17].This tendency is due to GSP de-localization in narrow waveguides.However, typical propagation distances dictated by the damping still remain relatively short, well below 1 µm.To observe any nonlinear effects over such short distances, one needs to have the characteristic nonlinear length L nl = 1/(γP ) to be considerably shorter than L p .While the nonlinear length is inversely proportional to input power P , the important characteristic of a waveguide is the critical power, for which local field intensities are still below the damage threshold of graphene and dielectrics involved.Setting the threshold intensity to I d = 10 14 W/m 2 [37], the minimal nonlinear distance is found to be several orders of magnitude below the attenuation length over the broad wavelength range in mid-infrared, see Figure 2d.In particular, for the waveguide of width w = 25 nm the optimal wavelength is identified at around λ = 10 µm, where max(L nl /L P ) ≈ 9000 and the maximal efficiency of nonlinear processes is expected.To illustrate the nonlinear functionality of GSP waveguides, here we consider the spectral broadening process of a short pulse propagating in a 25 nm wide waveguide.We set the initial pulse to be T 0 = 100 fs long (FWHM, sech pulse), centered at λ 0 = 9 µm and with a variable peak power P 0 .The nonlinear coefficient at this wavelength is γ ≈ 1.3 × 10 13 W −1 m −1 , and the critical power is P cr ≈ 4 mW, cf. also bottom panel if Figure 1.The attenuation constant is α ≈ 5.6 µm −1 , which gives the attenuation length L p ≈ 180 nm.The waveguide length was set to L = 150 nm.The dispersion at λ 0 = 9 µm is computed as β 2 = ∂ 2 β/∂ω 2 = 9.2 × 10 −21 s 2 /m.For the chosen pulse duration, this gives the dispersion length L d = T 2 0 /β 2 ≈ 1.1 µm [35].Hence the dispersive spreading of the pulse is practically negligible over the full waveguide length, which helps to sustain high efficiency of the nonlinear spectral broadening.In Figure 3 output spectra for different input peak powers are plotted.For the peak power as low as 0.5 mW a considerable spectral broadening can be observed.The generated mid-infrared supercontinuum covers more than one optical octave and spans from λ = 6 µm to λ > 13 µm.While this peak power is well below the estimated critical power P cr , a similar spectral broadening could be observed in samples with shorter attenuation lengths, e.g., due to reduced relaxation times τ in graphene, by increasing P 0 further.

Summary
We introduce a comprehensive set of analytical and numerical tools for the description of pulse propagation and frequency conversion processes in GSP waveguides of arbitrary geometry.The generic propagation Equation (28) and its simplified version (36), together with the nonlinear and damping coefficients due to graphene conductivity are derived following perturbation expansion of Maxwell equations with graphene being treated as the nonlinear boundary condition.Unlike the conventional and widely used models, originally derived for bulk materials [35], our approach does not require introduction of an artificial graphene thickness, and thus does not suffer from any related errors.
As an example, we considered GSP waveguides formed by a local modification of Fermi energy.Remarkably, such waveguides are found to be an outstanding platform for nonlinear functional and extremely compact plasmonic circuits in the mid-infrared to terahertz range.In particular, a supercontinuum generation covering more than one optical octave from a 100 fs input pulse with the peak power as low as 0.5 mW is demonstrated in a compact 25 nm wide and 150 nm long waveguide.The compactness of the waveguide and the low powers of operation could play the decisive role for future development of on-chip integrated light sources in the spectral range, which is of high importance for many applications in spectroscopy, chemical and bio-sensing, medical imaging.

Figure 1 .
Figure 1.Graphene surface plasmon (GSP) waveguide formed by a modulation of the graphene Fermi energy E F (see top panel), graphene is embedded into a dielectric (Al 2 O 3 ) with the relative permittivity = 10.The bottom panel illustrates field intensity distribution in the fundamental guided mode at λ = 9µm for the case of E F 1 = 0.25 eV, E F 2 = 0.2 eV and w = 25 nm, the total power of the mode is P = 1 mW.Cyan arrows indicate the in-plane direction and magnitude of electric field.

Figure 2 .
Figure 2. Linear and nonlinear parameters for GSP waveguides of different widths: effective index (a); attenuation distance (b); nonlinear coefficient (c); and the maximal ratio between nonlinear and attenuation distances (d).Dashed line in panels (a,b) indicate characteristics of the GSP on a homogeneous graphene sheet with σ 1 = σ 1,lo , Equation (41).In panel (b) circles indicate attenuation length for w = 25 nm as computed directly from Comsol, see text for details.

Figure 3 .
Figure 3. Numerically computed output spectra for propagation of a 100 fs sech pulse with the peak power P 0 in the 25 nm wide and 150 nm long GSP wavegide.