---
abstract: 'We suggest a novel approach for generating second-harmonic radiation in subwavelength graphene waveguides. We demonstrate that quadratic phase matching between the plasmonic guided modes of different symmetries can be achieved in a planar double-layer geometry when conductivity of one of the layers becomes spatially modulated. We predict theoretically that, owing to graphene nonlocal conductivity, the second-order nonlinear processes can be actualized for interacting plasmonic modes with an effective grating coupler to allow external pumping of the structure and output of the radiation at the double frequency.'
author:
- Daria Smirnova
- 'Yuri S. Kivshar'
title: 'Second-harmonic generation in subwavelength graphene waveguides'
---

Introduction
============

Graphene plasmonics is the rapidly developing field of nanophotonics that attracted a lot of attention during the past two years [@RevGrigorenko; @RevBao; @JablanReview; @RevLuo; @Stauber_Review; @Low_Review; @Abajo_Review_ACSPhot]. One of the recent novel directions in this field is the study of nonlinear phenomena given the fact that strong confinement of surface p-plasmons in graphene can enhance light-matter interaction [@Abajo176] and facilitate nonlinear response [@PRL_SinglePhoton]. In particular, the recent theoretical studies included the analysis of spatial plasmon solitons in graphene layers [@LPR; @DissipSoliton_LPR], nonlinear difference frequency generation of THz surface plasmons [@PRL_Tokman_2014], and single-photon operation in a graphene cavity [@PRL_SinglePhoton].

Double-layer graphene waveguides are known to support symmetric and antisymmetric plasmonic guiding modes classified with respect to the in-plane electric field profile[@Hanson; @buslaev_JETP; @Coupler_PRB], so we may expect the coupling between different modes when the nonlinear response becomes important. In this paper, we suggest and study analytically a novel approach for the second-harmonic generation in double-layer graphene waveguides taking the advantage of the possibility to modulate spatially the conductivity of one of the graphene layers by doping or electrostatic gating [@Gating_doping; @Engheta; @Chemical_doping] thus assisting the coupling of freely propagating light to plasmons in graphene [@bludov_ScatMod_2012; @bludov_PhotCr_2012; @Nikitin_grating_2013; @bludov_primer_2013].

Phased matching between parametrically interacting waves is known to be a crucial requirement for the efficient second-order nonlinear effects [@Boyd], and in our geometry the phase matching can be achieved by tailoring the modal dispersion of the graphene waveguide enabling plasmon-to-plasmon frequency conversion. More specifically, we demonstrate below that the phase matching becomes possible between an antisymmetric fundamental frequency (FF) mode and a symmetric second-harmonic (SH) mode in a double-layer graphene waveguide.

![(Color online) Geometry of the problem. A free-standing waveguide created by two graphene layers is illuminated by an external wave being coupled to the phase-matched guided modes of the graphene waveguide. Red curve corresponds to an antisymmetric mode at the fundamental frequency, and a blue curve – to a second-harmonic symmetric mode (shown with the tangential electric field component). Conductivity of the upper layer is periodically modulated.[]{data-label="fig:fig1"}](Fig1_GW.eps){width="0.95\linewidth"}

We employ a theoretical analysis based on the perturbation theory and describe the second-order nonlinear process with graphene plasmons in terms of slowly varying modes of a waveguide structure, the approximation frequently used in the nonlinear optics [@book]. For definiteness, we consider a planar geometry shown schematically in Fig. \[fig:fig1\], where a graphene double-layer waveguide is placed into homogeneous surrounding medium with the dielectric permittivity $\varepsilon$ being illuminated by light from the upper half-space $x>d/2$. As a matter of fact, in our analysis we distinguish two parts of the problem and examine them sequentially,

-   [*Linear scattering.*]{} Radiation at the fundamental frequency normally incident onto the structure is scattered by the periodic conductivity grating, and then excites resonantly an antisymmetric mode of the graphene waveguide;

-   [*Nonlinear plasmon-to-plasmon conversion.*]{} Due to the second-order nonlocal nonlinearity of graphene [@Glazov2011_SHG; @Mikhailov2011_SHG; @Glazov2013], the induced current of the antisymmetric mode serves as a source for the phase-matched symmetric mode at the double frequency which eventually radiates into free space.

Similar to the assumptions employed earlier [@DissipSoliton_LPR; @Coupler_PRB; @Gorbach2013], namely low dissipation, weak nonlinearity, and small phase mismatch, here we solve Maxwell’s equations following the procedure of the asymptotic expansion, and demonstrate that the second-order nonlinear processes can be achieved with conductivity modulation to allow external pumping to couple to the guided modes of the structure generating output radiation at the doubled frequency.

Linear scattering
=================

We start our derivation by studying the linear resonant excitation of antisymmetric plasmons, and write the corresponding system of Maxwell’s equations in the form $$\begin{aligned} \label{eq:eqMaxw2}
&\nabla \times {\bf E} =  - \frac{1}{c} \frac{\partial {\bf H}}{\partial t}\:,\\
&\nabla \times {\bf H} = \frac{\varepsilon }{c} \frac{\partial {\bf E}}{\partial t}\!+\!\frac{4\pi }{c} \biggl[j^{(1)} \delta \left(x-\dfrac{d}{2}\right) + j^{(2)} \delta
\left(x+\dfrac{d}{2}\right) \biggr] {\bf z}_0\:,
\end{aligned}$$ where $j^{(1,2)}$ are the surface current densities induced in the graphene layers placed at $x= \pm d/2$, as indicated by Dirac’s delta functions $\delta$. Assuming the field to be p-polarized with the magnetic component ${\bf H} = H(x,z,t) {\bf y}_0$, from Eq. (\[eq:eqMaxw2\]) we find $$\label{eq:eqMaxwH}
\frac{\varepsilon}{c^2} \frac{\partial^2 H}{\partial t^2} - \Delta H \! = \!  - \frac{4\pi }{c}  \frac{\partial }{\partial x} \left[ j^{(1)} \delta \left(x-\dfrac{d}{2}\right) + j^{(2)} \delta \left(x+\dfrac{d}{2}\right) \right]\:,$$ where, if we assume the harmonic field dependencies $\sim \exp(-i \omega t)$, the currents are given by $$j^{(1,2)} (\omega, z) = \sigma^{(1,2)} (\omega, z) E_z (x = \pm {d}/{2}, z, \omega) \:.$$ Here $\sigma^{(2)} (\omega, z) = \sigma (\omega ) $, $\sigma^{(1)} (\omega, z) = \sigma (\omega ) [1 + f(z)] $ are surface conductivities, $\sigma (\omega ) \equiv \sigma^{(R)} (\omega )  + i \sigma^{(I)}(\omega ) $ is the linear frequency-dependent surface conductivity of graphene, while $ f(z) = f(z + a) = f_1 \cos \left(\dfrac {2 \pi } {a} z \right)+ f_2 \cos \left(\dfrac {4 \pi } {a} z \right)$ is assumed to be a periodic function of $z$, i.e. conductivity of one of the graphene layers is spatially modulated.

Assuming the dissipative losses and modulation amplitudes $f_{1,2}$ to be small, we formally introduce the smallness parameter $\mu$ $$\mu = \text{max}  \left\{   \left| \dfrac{\sigma^{(R)}}{\sigma^{(I)}}\right|, f_{1,2} \right\}\:,$$ and adopt the following asymptotic ansatz for the magnetic field $$\begin{aligned}  \label{eq:ansatz}
H(x,z,t)& = e^{-i \omega t} \biggl \{ \mu \bar{H}_0 (x, \mu z)  \\
&  + \left[  \mathcal{ A }_{1}(\mu z) h(\omega, x) + \mu \tilde {H}_1 (x, \mu z) \right] e^{i k_{\text{sp}}(\omega) z}  \\
& + \left[  \mathcal{ A }_{2}(\mu z) h(\omega, x) + \mu \tilde {H}_2 (x, \mu z) \right] e^{- i k_{\text{sp}}(\omega) z} \biggr \}\:,
\end{aligned}$$ where $h(x)$ is the transverse profile of the linear plasmonic mode of the guiding structure, $\mathcal{ A }_{1,2}$ are slowly varying mode amplitudes, where the subscripts $"1"$ and $"2"$ refer to the forward and backward propagation along the $z$ direction, respectively. Also $\tilde {H}_{1,2}$ are small corrections to the eigenmode profile, and the term $\bar{H}_0 $ does not contain a fast dependence on $z$.

In the zero-order approximation in $\mu$, for the function $h(x)$ Eq. (\[eq:eqMaxwH\]) takes the form, $$\begin{aligned}
& \frac{d^2 h}{dx^2} -  (k_{\text{sp}}^2 - k_0^2 \varepsilon )h(x) = \\
& = \frac{4\pi }{c} i \sigma^{(I)} (\omega ) \frac{\partial }{\partial x}
\biggl\{ \left[{\delta} \left(x-\dfrac{d}{2}\right) + {\delta} \left(x+\dfrac{d}{2}\right) \right] e_z(x) \biggr\}
\:,
\end{aligned}$$ where $e_z(x) = \displaystyle{\frac{i}{k_0 \varepsilon} \frac{dh}{dx} }$. This equation yields the dispersion relation $$\label{eq:DispRel1}
\left[1 - \dfrac{2\pi \sigma^{(I)} \kappa^{{\text{(s,a)}}}}{\omega \varepsilon}(1 \pm e^{-\kappa^{\text{(s,a)}} d})\right] = 0 \:,$$ where $\kappa^{\text{(s,a)}} = \sqrt{k^{\text{(s,a)} 2}_{\text{sp}} - k_0^2 \varepsilon}$, $k_0 = \omega / c$, and the modal transverse profiles are $$h^{\text{(s,a)}}(\omega, x) = \frac{-ik_0 \varepsilon}{\kappa^{\text{(s,a)} 2}} \frac{de^{\text{(s,a)}}_z}{dx} \:,$$ where the continuous tangential components of the electric field are expressed as $$\label{eq:TrStr1}
 e^{\text{(s)}}_z(\omega, x) =
    \begin{cases}
        e^{-\kappa^{\text{(s)}} (x-d/2)},&\text{}\ x>  d/2\:,\\
         \displaystyle{\frac{\cosh(\kappa^{\text{(s)} }x)}{\cosh(\kappa^{\text{(s)}} d/2)}},&\text{}\ |x|<  d/2\:,\\
        e^{\kappa^{\text{(s)}} (x+d/2)},&\text{} x< - d/2\:
    \end{cases}$$ and $$\label{eq:TrStr1}
  e^{\text{(a)}}_z(\omega,x) =
    \begin{cases}
        e^{-\kappa^{\text{(a)}} (x-d/2)},&\text{}\ x>  d/2\:,\\
         \displaystyle{\frac{\sinh(\kappa^{\text{(a)}} x)}{\sinh(\kappa^{\text{(a)}} d/2)}},&\text{}\ |x|<  d/2\:,\\
        -e^{\kappa^{\text{(a)}} (x+d/2)},&\text{} x< - d/2\:
    \end{cases}$$ for both symmetric and antisymmetric eigenmodes guided by a double-layer waveguide [@buslaev_JETP; @MultiWg], as indicated by the superscripts.

We now assume that the period of the conductivity grating is chosen to almost compensate the momentum mismatch with the antisymmetric plasmon, $\displaystyle{|k_{\text{sp}}^{\text{(a)}}(\omega) - {2 \pi}/{a} | \equiv |\Delta k| \ll {2 \pi}/{a} }$ so that $|\Delta k /k_{\text{sp}}^{\text{(a)}}(\omega)|  \lesssim \mu $. Substituting the expression of the total field (\[eq:ansatz\]) into Eq. (\[eq:eqMaxwH\]) and keeping only resonant terms, in the first order in $\mu$ for the field perturbations $\tilde{H}_{1,2}$, we obtain $$\label{eq:H2}
\begin{aligned}
& \dfrac{d^2 \tilde{H}_{1,2}}{dx^2} -  (k_{\text{sp}}^{{\text{(a)}}2} (\omega)- k_0^2 \varepsilon ) \tilde{H}_{1,2} - \dfrac{4 \pi}{c} i \sigma^{(I)} (\omega )  \\
& \times \frac{\partial }{\partial x}
\biggl\{ \left[{\delta}\left(x - \dfrac{d}{2}\right) + {\delta}\left(x + \dfrac{d}{2}\right) \right]  \tilde{E}_{1,2z} \biggr\} = F_{1,2}\:,\\
& F_{1,2} (x,z) = \mp 2ik_{\text{sp}}^{\text{(a)}}(\omega) \dfrac{\partial \mathcal A_{1,2}}{\partial z} h^{\text{(a)}}(\omega, x) + \frac{4 \pi}{c}  \sigma^{(R)} (\omega ) \mathcal A_{1,2} \\
&   \times  \frac{\partial }{\partial x}  \left\{ \left[ {\delta}\left(x - \dfrac{d}{2}\right) + {\delta} \left(x+\dfrac{d}{2}\right) \right]  e^{\text{(a)}}_z(\omega, x) \right\}  \\
&   + \frac{4 \pi}{c}  i \sigma^{(I)} (\omega )  \frac{\partial }{\partial x} \biggl\{  {\delta}\left(x - \dfrac{d}{2}\right) \left[ \frac{f_1} {2}\bar{E}_{0z} e^{\mp i\Delta k z} \right.  \biggr.   \\
& + \frac{f_2} {2}\mathcal A_{2,1}
\left. \biggl. e^{\mp 2i\Delta k z} e^{\text{(a)}}_z(\omega, x)  \right] \biggr\} \:,
\end{aligned}$$ whereas $\bar{H}_{0}$ satisfies the equation, $$\label{eq:H0}
\begin{aligned}
\dfrac{d^2 \bar{H}_{0}}{dx^2} &\!+\! k_0^2 \varepsilon \bar{H}_{0} \! = \! \dfrac{4 \pi}{c} i \sigma^{(I)} (\omega ) \frac{\partial }{\partial x}
\biggl \{ \left [{\delta}\left(x + \dfrac{d}{2}\right) \!+ \!{\delta}\left(x - \dfrac{d}{2}\right) \right]\\
&  \times  {E}_{0z} (\omega,x)
\biggl. +\frac{f_1} {2} \left(\mathcal A_{1} e^{ i\Delta k z}  + \mathcal A_{2} e^{-i\Delta k z} \right) {\delta} \left(x-\dfrac{d}{2}\right)  \biggr \} \:.
\end{aligned}$$ Equation (\[eq:H0\]) includes reflection of a plane wave incident from the semi-space $x>d/2$, written as $H_0{\rm exp} [-ik_{0}\sqrt{\varepsilon}(x-d/2)]$, where we assume $H_0 \sim \mu$, $E_0 = H_0/\sqrt{\varepsilon}$, ${E}_{0z} (\omega,x)$ is the tangential electric field distribution at $f_{1,2}=0$. Taking into account the boundary conditions at the graphene interface, we find the amplitude of the electric field at $x=d/2$ as follows $$\label{eq:E0}
\begin{aligned}
 & \bar{E}_{0z} (\omega, x=d/2) = E_0(1+r) \\
 & - \frac{2 \pi}{c \sqrt{\varepsilon}} i  \sigma^{(I)}(\omega )  \frac{f_1} {2} \left(\mathcal A_{1} e^{ i\Delta k z}  + \mathcal A_{2} e^{-i\Delta k z} \right) \:.
\end{aligned}$$ where $r = - \displaystyle{{ \frac{4\pi} {c \sqrt{\varepsilon}} i \sigma^{(I)}(\omega ) }\left({1+\frac{4\pi} {c \sqrt{\varepsilon}}  i \sigma^{(I)} (\omega ) } \right)^{-1} } \approx - { \frac{4\pi} {c \sqrt{\varepsilon}} i \sigma^{(I)}(\omega ) }$ stands for the reflection coefficient from two closely spaced graphene layers, $k_0 \sqrt{\varepsilon} d  \ll 1$. Substituting the expression (\[eq:E0\]) into Eq. (\[eq:H2\]), in order for the corrections $\tilde{H}_{1,2}$ to be nondiverging, in accord with the Fredholm alternative [@Korn], we have to satisfy the orthogonality condition for the right-hand side of Eq. (\[eq:H2\]) with the solution of its homogeneous part, or equivalently, with the plasmonic mode itself. This overintegrating is mathematically written as $$\int\limits_{-\infty}^{+\infty}F_{1,2}(x,z){h^{\text{(a)}}}^{*}(\omega, x)dx = 0,$$ and it leads to the nonlinear equations for the slowly varying amplitudes $\mathcal A_{1,2}$ $$\begin{gathered}
 \label{eq:EqA}
\pm 2ik^{\text{(a)}}_{\text{sp}} (\omega) \left(\dfrac{\partial \mathcal A_{1,2}}{\partial z} \pm (\gamma^{\text{(a)}}  + \gamma^{\text{(a)}} _{r} ) \mathcal A_{1,2} \right)  {} \\
 {}\!=\!(\!  - 2ik^{\text{(a)}}_{\text{sp}} (\omega) \gamma^{\text{(a)}}_{r}\!+Q^{\text{(a)}} _2 ) \mathcal A_{2, 1} e^{ \mp 2 i\Delta k z} + Q^{\text{(a)}} _1 E_0 (1+r) e^{ \mp i\Delta k z} \:,\end{gathered}$$ where linear damping due to losses in graphene $\gamma^{\text{(a)}}$, grating-induced damping due to radiation $\gamma^{\text{(a)}}_r$, and coupling coefficients $Q^{\text{(a)}} _{1,2} $ are derived, respectively, in the form $$\begin{gathered}
\gamma^{\text{(a)}}  =  \dfrac{4\pi}{c}  \sigma^{(R)} (\omega ) \dfrac{k_0 \varepsilon}{k^{\text{(a)}}_{\text{sp}}(\omega) q^{\text{(a)}} } \:,\nonumber \\
\gamma^{\text{(a)}}_r =  \left(\dfrac{\pi \sigma^{(I)}(\omega ) f_{1} }{ c} \right)^2 \dfrac{k_0 \varepsilon^{1/2}}{k^{\text{(a)}}_{\text{sp}} (\omega) q^{\text{(a)}} } \:, \nonumber \\
Q^{\text{(a)}}_{1,2} = \dfrac{2\pi}{c} \sigma^{(I)} (\omega )  f_{1,2} \dfrac{k_0 \varepsilon}{ q^{\text{(a)}}} \:,\end{gathered}$$ where in all the expressions we employ the following definition, $$q^{\text{(a)}} = {k_0^2} \dfrac{\varepsilon^2}{\kappa^{\text{(a)}3} } \left[ 1 + \dfrac{\sinh{(\kappa^{\text{(a)}} d)} + \kappa^{\text{(a)}} d }{2 \sinh^2 {(\kappa^{\text{(a)}} d /2)}} \right].$$

When we assume a simple homogeneous case, $\dfrac{\partial }{\partial z} =0$, and also the exact momentum matching, $\Delta k=0$, from Eq. (\[eq:EqA\]) we obtain the slowly varying amplitudes, $$\mathcal A_{1} = \mathcal A_{2} \equiv \mathcal A = \dfrac{Q^{\text{(a)}}_{1}  E_0 (1+r)}{2ik^{\text{(a)}}_{\text{sp}}  (\omega) [\gamma^{\text{(a)}} (\omega) + 2 \gamma^{\text{(a)}}_r] - Q^{\text{(a)}}_{2}}$$ and substitute them into Eq. (\[eq:E0\]). Thus, the excitation of the antisymmetric mode results in changing the reflection coefficient compared to the case without a conductivity grating, $r_{\mathcal A} \approx r(1 + f_1 \mathcal A /2 E_0)$.

Nonlinear plasmon-to-plasmon conversion
=======================================

First we notice that, for the case of TM-polarized waves with the tangential electric field of the form $\mathbf{E}_{\tau \omega} = \mathbf{z}_0 E e^{iqz}$ (the monochromatic time dependence $\sim \exp (-i \omega t)$ is omitted), the induced second-harmonic current in graphene depends on the in-plane momentum $q$ manifesting its nonlocal nature. It is expressed as [@Mikhailov2011_SHG; @SHG_Wrapped] $${\bf j}_{2 \omega} =  - \mathbf{z}_0 \dfrac{3} {8} \dfrac{e^3 v_{F}^2} {\pi \hbar^2 \omega^3} q E^2 e^{i2qz} \:,$$ where $v_{F}\approx c/300$ is the Fermi velocity. This result is derived in the semiclassical limit from the Boltzmann kinetic equation written for Dirac electrons with a linear energy dispersion under the approximations $\hbar \omega \leq \mathcal{E}_{F}$, $ k_B T \ll \mathcal{E}_{F}$ and $q v_{F}/ \omega \ll 1$, where $\mathcal{E}_{F}$ is the Fermi energy, $k_B$ is the Boltzmann constant, $T$ is the temperature. Thereby, a nonlinear source associated with the antisymmetric mode will generate a synchronous symmetric mode. Importantly, in general, when discussing the parametric interaction that involves the modes of different symmetries, the phase matching is not a sufficient condition, and the overlaps of the modes with the respective nonlinear sources should be examined. Since a quadratic nonlinear source localized at the graphene layers is symmetric, the conversion between the FF antisymmetric mode and SH symmetric mode becomes possible.

Similar to the analysis of the fundamental frequency field, we employ here the perturbational method and derive the equations for the envelopes of the second-harmonic fields $\mathcal B_{1, 2}$, $$\begin{gathered}
 \label{eq:EqB}
\pm 2ik^{\text{(s)}}_{\text{sp}} (2\omega)  \left(\dfrac{\partial \mathcal B_{1,2}}{\partial z} \pm (\gamma^{\text{(s)}}  + \gamma^{\text{(s)}} _{r} ) \mathcal B_{1,2} \right)  {} \\
 {} =  - 2ik^{\text{(s)}}_{\text{sp}} (2\omega) \gamma^{\text{(s)}}_{r}  \mathcal B_{2, 1} e^{ \mp 2 i (\widetilde{\mathtt{\Delta} k } + 2 \Delta k) z} + g \mathcal A_{1, 2}^2 e^{ \mp  i \widetilde{\mathtt{\Delta} k } z}    \:,\end{gathered}$$ where $\widetilde{ \mathtt{\Delta} {k} }= k^{\text{(s)}}_{\text{sp}}(2\omega) - 2 k^{\text{(a)}}_{\text{sp}}(\omega)$ is a small detuning, and other parameters are $$\begin{gathered}
q^{\text{(s)}} = 4 {k_0^2} \dfrac{\varepsilon^2}{\kappa^{\text{(s)}3}} \left[ 1 + \dfrac{ \sinh{(\kappa^{\text{(s)}} d)} - \kappa^{\text{(s)}} d }{2 \cosh^2 {(\kappa^{\text{(s)}} d /2)}} \right]\:, \nonumber   \\
\gamma^{\text{(s)}}  =  \dfrac{8\pi}{c}  \sigma^{(R)} (2 \omega ) \dfrac{k_0 \varepsilon}{k^{\text{(s)}}_{\text{sp}} (2\omega) q^{\text{(s)}} }\:,\nonumber \\
\gamma^{\text{(s)}}_r =  \left(\dfrac{\pi \sigma^{(I)}(2 \omega ) f_{2} }{ c} \right)^2 \dfrac{2 k_0 \varepsilon^{1/2}}{k^{\text{(s)}}_{\text{sp}} (2\omega) q^{\text{(s)}} }\:,\nonumber \\
g = -i \dfrac{16\pi}{c} \sigma_{2}(\omega) \dfrac{k_0 \varepsilon}{ q^{\text{(s)}} }\:,\end{gathered}$$ with $\sigma_{2} (\omega)$ denoting the quadratic conductivity of graphene defined as $$\sigma_{2} (\omega) = -  \dfrac{3} {8} \dfrac{e^3 v_{F}^2} {\pi \hbar^2 \omega^3} k^{\text{(a)}}_{\text{sp}} (\omega).$$ Within the undepleted pump approximation, when the amplitude at fundamental frequency is not affected by nonlinearity, for the exact phase matching, $\widetilde{\mathtt{\Delta} k }= \Delta k =0$, and homogeneous case discussed above, we have $$\mathcal B_{1} = \mathcal B_{2} \equiv \mathcal B = \dfrac{g \mathcal A^2}{2ik^{\text{(s)}}_{\text{sp}} (2\omega) (\gamma^{\text{(s)}} + 2 \gamma^{\text{(s)}}_r) }.$$

Since the radiation associated with SH field is emitted equally to both directions, the normalized conversion efficiency can be introduced as a ratio of the doubled SH energy flux density in the upper half-space to the energy flux of the incident fundamental wave. Equivalently, this definition takes the following form $$\eta = \dfrac{2|\bar{E}_{0z} (2\omega, x=d/2)|^2}{|{E}_{0}|^2},$$ where $\bar{E}_{0z} (2\omega, x=d/2) = - \dfrac{2 \pi}{c \sqrt{\varepsilon}} i  \sigma^{(I)}(2\omega )  f_2 \mathcal B $.

Discussions
===========

![(Color online) Plasmonic guide index vs. frequency at the fixed separation between the graphene layers, $d=123$ nm, calculated for the doping level of $\mathcal{E}_{F}=0.3$ eV, $N = 5$, $\varepsilon = 1$. Blue and red solid curves correspond to the symmetric and antisymmetric modes, respectively. Dashed gray and red lines show dispersion of an isolated graphene sheet and half-frequency antisymmetric mode, for reference. Color dots in the horizontal axis mark the fundamental and second-harmonic frequencies.[]{data-label="fig:PhMatch"}](Fig2_GrWg_SHG){width="0.95\linewidth"}

![(Color online) Operational fundamental frequency as a function of the separation between the layers for the exact phase-matching at $\mathcal{E}_{F}=0.3$ eV, $N = 5$, $\varepsilon = 1$. The point corresponds to the parameters used in Fig. \[fig:PhMatch\]. []{data-label="fig:dDep"}](Fig3_GrWg_SHG_v4){width="0.95\linewidth"}

In the framework of the nonlinear amplitude equations, we can now analyze the parametric effects and estimate the frequency conversion efficiency. To make a period of the conductivity modulation realistic, in our calculations we employ multilayer graphene that effectively increases the equivalent surface conductivity leading to a reduction of the wavenumber of the p-polarized plasmons supported by multilayer graphene structures [@DissipSoliton_LPR; @MultiWg]. We take ${\sigma(\omega)}= N {\sigma_s}(\omega)$ for a randomly stacked multilayer graphene film consisting of $N$ layers, with each layer characterized by a surface conductivity [@sptdep] $$\sigma_{s}(\omega) = \displaystyle{\frac{ie^2}{\pi \hbar}
\left\{\frac{\mathcal{E}_{F}}{\hbar\left(\omega + i\tau_{\text{intra}}^{-1}\right)} + \frac{1}{4} \text{ln}\left|\frac{2\mathcal{E}_{F} -\hbar \omega}{2\mathcal{E}_{F} + \hbar \omega}\right|\right\}}\:,  \label{eq:sigma}$$ where for doped graphene we assume $\hbar\omega < 1.67\mathcal{E}_{F} $ and $k_B T \ll \mathcal{E}_{F}$, and $\displaystyle{\tau_{\text{intra}} ^{-1}}$ is the relaxation rate.

To reveal a possibility of plasmonic frequency conversion in the graphene waveguide, in Fig. \[fig:PhMatch\] we plot the frequency-dependence of the plasmonic guide indices for the symmetric and antisymmetric modes. We observe that the first harmonic antisymmetric and second-harmonic symmetric modes of the waveguide can be matched at $\hbar \omega = 0.12$ eV ($\lambda=2\pi/k_0\approx 10.6$ $\mu$m, corresponding to one of the radiation lines of high power $\text{CO}_2$ lasers), $ k^{\text{(a)}}_{\text{sp}}  \approx 10.5 k_0$.

Figure \[fig:dDep\] shows the dependence of the fundamental frequency on the separation between the layers. Note while discussing the dispersion relations and the possibility of phase-matching conditions we take $\sigma = i  \sigma^{(I)} $ and do not take losses into account, however these losses are included into our considerations perturbatively. Calculated for $f_1=f_2=0.1$, $a \approx 1 $ $\mu$m, $\tau_{\text{intra}} = 0.3$ ps, the incident wave intensity of $1$ MW/cm$^2$ and the other parameters of Fig. \[fig:PhMatch\], efficiency is $\eta \sim 8.7 \times 10^{-7}$. Overall, for this structure the smallness parameter is found to be $\mu \sim 0.1$, so that our asymptotic approach is well justified.

Conclusions
===========

By employing the approximation of slowly-varying field amplitudes, we have derived nonlinear equations for describing the second-harmonic generation from a double-layer graphene structure with modulated conductivity. We have revealed that the quadratic phase matching between the plasmonic modes of different symmetries becomes possible in a planar waveguide geometry with conductivity modulation playing a role of an effective grating that couples the external pumping radiation to the waveguiding modes with the subsequent parametric conversion into radiation at the double frequency.

ACKNOWLEDGMENTS {#acknowledgments .unnumbered}
===============

The authors are grateful to I.V. Shadrivov and A.I. Smirnov for in-depth comments, and also thank S.A. Mikhailov for useful discussions at an early stage of this project. This work was supported by the Australian National University.

[100]{} A. N. Grigorenko, M. Polini, and K. S. Novoselov, Nature Photon. **6**, 749 (2012).

Q. Q. Bao and K. P. Loh, ACS Nano **6**, 3677 (2012).

M. Jablan, M. Soljacic, and H. Buljan, Proc. IEEE **101**, 1689 (2013).

X. Luo, T. Qiu, W. Lu, and Z. Ni, Mater. Sci. Eng., R **74**, 351 (2013).

T. Stauber, J. Phys.: Cond. Mat. **26**, 123201 (2014). T. Low and P. Avouris, ACS Nano **8**, 1086 (2014).

F. J. Garcia de Abajo, ACS Photonics **1**, 135 (2014).

F. H. L. Koppens, D. E. Chang, and F. J. Garcia de Abajo, Nano Lett. **11**, 3370 (2011).

M. Gullans, D. E. Chang, F. H. L. Koppens, F. J. Garcia de Abajo, and M. D. Lukin, Phys. Rev. Lett. **111**, 247401 (2013).

M. L. Nesterov, J. Bravo-Abad, A. Yu. Nikitin, F. J. Garcia-Vidal, and L. Martin-Moreno, Laser Phot. Rev. [**7**]{}, L7 (2013).

D.A. Smirnova, I.V. Shadrivov, A.I. Smirnov, and Y.S. Kivshar, Laser Phot. Rev. **8**, 291 (2014).

X. Yao, M. Tokman, and A. Belyanin, Phys. Rev. Lett. **112**, 055501 (2014).

G.W. Hanson, J. Appl. Phys. **104**, 084314 (2008).

P.I. Buslaev, I.V. Iorsh, I.V. Shadrivov, P.A. Belov, and Yu.S. Kivshar, JETP Lett. **97**, 535 (2013).

D.A. Smirnova, A.V. Gorbach, I.V. Iorsh, I.V. Shadrivov, and Y.S. Kivshar, Phys. Rev. B **88**, 045443 (2013).

F. Wang, Y. Zhang, C. Tian, C. Girit, A. Zettl, M. Crommie, and Y. R. Shen, Science **320**, 206 (2008). A. Vakil and N. Engheta, Science **332**, 1291 (2011).

H. Liu, Y. Liu, and D. Zhu, J. Mater. Chem. **21**, 3335–3345 (2011). Y. V. Bludov, N. M. R. Peres, and M. I. Vasilevskiy, Phys. Rev. B **85**, 245409 (2012).

N. M. R. Peres, A. Ferreira, Y. V. Bludov, and M. I. Vasilevskiy, J. Phys.: Cond. Mat. **24**, 245303 (2012).

T. M. Slipchenko, M. L. Nesterov, L. Martin-Moreno, and A. Yu. Nikitin, J. Opt. **15**, 114008 (2013).

Yu. V. Bludov, A. Ferreira, N. M. R. Peres, and M. I. Vasilevskiy, Int. J. Mod. Phys. B **27**, 1341001 (2013).

R. W. Boyd, [*Nonlinear Optics*]{} (Elsevier Science, USA, 2003).

Yu. S. Kivshar and G. A. Agrawal, [*Optical Solitons: From Fibers to Photonic Crystals*]{} (Academic, San Diego, 2003).

M.M. Glazov, JETP Lett. **93**, 366 (2011).

S.A. Mikhailov, Phys. Rev. B **84**, 045432 (2011).

M.M. Glazov and S.D. Ganichev, Phys. Rep. [**535**]{}, 101 (2014).

A.V. Gorbach, Phys. Rev. A **87**, 013830 (2013).

D.A. Smirnova, I.V. Iorsh, I.V. Shadrivov, and Y.S. Kivshar, JETP Lett. **99**, 456 (2014).

G.A. Korn and T.M. Korn, [Mathematical Handbook for Scientists and Engineers]{} (McGraw-Hill, New York, 1961).

D.A. Smirnova, I.V. Shadrivov, A.E. Miroshnichenko, A.I. Smirnov, and Y.S. Kivshar, Phys. Rev. B **90**, 035412 (2014).

L. A. Falkovsky and A. A. Varlamov, , **56**, 281 (2007).
