The spatio-temporal dependencies of the plasmon fields in the groove waveguide are assumed in the form \(\,e^{{ – i\omega t + ik_{x} x + ik_{y} y}} \cdot \left( {A^{ + } e^{{ik_{z} z}} + A^{ – } e^{{ik_{z} z}} } \right)\) above graphene (\(A^{ + }\) and \(A^{ – }\) are the amplitudes of electric or magnetic fields of the forward and the counter-propagating waves, respectively, along z axis) and \(e^{{ – i\omega t + ik_{x} x}} \cdot \left[ {\left( {A^{ + + } e^{{ik_{z} z}} + A^{ + – } e^{{ik_{z} z}} } \right) \cdot e^{{ik_{y} y}} + \left( {A^{ – + } e^{{ik_{z} z}} + A^{ – – } e^{{ik_{z} z}} } \right) \cdot e^{{ – ik_{y} y}} } \right]\) in the substrate (\(A^{ + + } ,\,\)\(A^{ + – } ,\,\)\(A^{ – + } ,\,\) and \(A^{ – – }\) are the amplitudes of fields of the forward and the counter-propagating waves along the y and z axes, where the first superscript indicates the direction along the y axis, the second index does for the z axis). We use the following boundary conditions. The components of the electric field tangential to the metal surfaces are set to zero at the walls and bottom of the groove waveguide: \(E_{x,z} \left( {y = – d} \right) = 0,\) \(E_{x,y} \left( {z = 0} \right) = 0,\) \(E_{x,y} \left( {z = L} \right) = 0.\) The tangential components of the electric field are continuous across graphene: \(E_{{x{\text{a}}}} \left( {y = 0} \right) = E_{{x{\text{s}}}} \left( {y = 0} \right),E_{{z{\text{a}}}} \left( {y = 0} \right) = E_{{z{\text{s}}}} \left( {y = 0} \right).\) The discontinuity of the tangential components of the magnetic field across graphene are equal to the current density components in graphene \(H_{{z{\text{a}}}} \left( {y = 0} \right) – H_{{z{\text{s}}}} \left( {y = 0} \right) = \sigma_{{{\text{gr}}}} \left( \omega \right) \cdot E_{x} \left( {y = 0} \right);\)\(H_{{x{\text{a}}}} \left( {y = 0} \right) – H_{{x{\text{s}}}} \left( {y = 0} \right) = – \sigma_{{{\text{gr}}}} \left( \omega \right) \cdot E_{z} \left( {y = 0} \right).\) Finally, we derive the following dispersion equation for the LSM modes
$$\frac{{\varepsilon_{{\text{s}}} }}{{k_{{y\,{\text{s}}}} }}{\text{coth}}\;\left( {k_{{y\,{\text{s}}}} d} \right) + \frac{{\varepsilon_{{\text{a}}} }}{{k_{{y\,{\text{a}}}} }} = – \frac{{\sigma_{{{\text{gr}}}} \left( \omega \right)}}{{\omega \varepsilon_{0} }},$$
(2)
where \(k_{{y\,{\text{a,s}}}} = \sqrt {\varepsilon_{{\text{a,s}}} \left( {{\omega \mathord{\left/ {\vphantom {\omega c}} \right. \kern-\nulldelimiterspace} c}} \right)^{2} – k_{x}^{2} – k_{z}^{2} }\) are the wavevector components normal to the graphene plane, the sign before the radical for the normal to graphene wavevector component above graphene is chosen to meet the condition of the surface wave, \(k_{z} = {\pi \mathord{\left/ {\vphantom {\pi L}} \right. \kern-\nulldelimiterspace} L}\) is the transverse (along the z axis) wavevector component (we study the fundamental mode along z axis in this paper, only), ε0 is the electric constant, \(\sigma_{{{\text{gr}}}} \left( \omega \right)\) is the dynamic conductivity of active graphene given by6:
$$\sigma_{{{\text{gr}}}} \left( \omega \right) = \frac{{e^{2} }}{4\hbar }\left\{ {\frac{{8ik_{{\text{B}}} T}}{{\pi \hbar {\kern 1pt} \left( {\omega + {i \mathord{\left/ {\vphantom {i \tau }} \right. \kern-\nulldelimiterspace} \tau }} \right)}}\ln } \right.\left( {1 + \exp \;\left( {\frac{{E_{F} }}{{k_{{\text{B}}} T}}} \right)} \right)\left. { + \tanh \;\left( {\frac{{\hbar {\kern 1pt} \omega – 2E_{F} }}{{4k_{{\text{B}}} T}}} \right) + \frac{4i\hbar \omega }{{\pi }}\int\limits_{0}^{\infty } {\frac{{G\left( {\varepsilon ,E_{F} } \right) – G\left( {{{\hbar {\kern 1pt} \omega } \mathord{\left/ {\vphantom {{\hbar {\kern 1pt} \omega } 2}} \right. \kern-\nulldelimiterspace} 2},E_{F} } \right)}}{{\left( {\hbar \omega } \right)^{2} – 4\varepsilon^{2} }}d\varepsilon } } \right\},$$
(3)
where e is the elementary charge, \(\hbar\) is the reduced Planck constant, \(k_{{\text{B}}}\) is the Boltzmann constant τ, and T are the mean free time and temperature of the charge carriers in graphene, respectively, \(E_{{\text{F}}} > 0\) is the quasi-Fermi energy in graphene determining the inversion of the charge carriers (+ EF and –EF for electrons and holes, respectively), \(G\left( {\varepsilon ,\varepsilon ^{\prime}} \right) = {{\sinh \;\left( {{\varepsilon \mathord{\left/ {\vphantom {\varepsilon {k_{{\text{B}}} T}}} \right. \kern-\nulldelimiterspace} {k_{{\text{B}}} T}}} \right)} \mathord{\left/ {\vphantom {{\sinh \;\left( {{\varepsilon \mathord{\left/ {\vphantom {\varepsilon {k_{{\text{B}}} T}}} \right. \kern-\nulldelimiterspace} {k_{{\text{B}}} T}}} \right)} {\left[ {\cosh \;\left( {{\varepsilon \mathord{\left/ {\vphantom {\varepsilon {k_{{\text{B}}} T}}} \right. \kern-\nulldelimiterspace} {k_{{\text{B}}} T}}} \right) + \cosh \;\left( {{{\varepsilon ^{\prime}} \mathord{\left/ {\vphantom {{\varepsilon ^{\prime}} {k_{{\text{B}}} T}}} \right. \kern-\nulldelimiterspace} {k_{{\text{B}}} T}}} \right)} \right]}}} \right. \kern-\nulldelimiterspace} {\left[ {\cosh \;\left( {{\varepsilon \mathord{\left/ {\vphantom {\varepsilon {k_{{\text{B}}} T}}} \right. \kern-\nulldelimiterspace} {k_{{\text{B}}} T}}} \right) + \cosh \;\left( {{{\varepsilon ^{\prime}} \mathord{\left/ {\vphantom {{\varepsilon ^{\prime}} {k_{{\text{B}}} T}}} \right. \kern-\nulldelimiterspace} {k_{{\text{B}}} T}}} \right)} \right]}}.\) The integral in the right-hand side of dynamic conductivity formula has to be interpreted as the principal value integral. The first term in the curly braces in Eq. (3) describes a Drude-model response for the intraband processes involving the phenomenological electron and hole scattering time τ. The remaining terms in Eq. (3) arise from the interband transitions. The second term in the curly braces becomes negative for \(\hbar {\kern 1pt} \omega < 2E_{F}\) which corresponds to the population inversion in graphene. At sufficiently large \(E_{F}\) (strong pumping), the interband stimulated emission of photons (plasmons) can prevail over the intraband (Drude) absorption. In this case, the real part of the dynamic conductivity of graphene can be negative in the THz range which leads the plasmon amplification.
Note that the dispersion expression for the TM mode in a layered graphene structure reads the same as Eq. (2) given there is only one wavevector component transverse to the plasmon propagation direction: \(k_{{y\,{\text{a,s}}}} = \sqrt {\varepsilon_{{\text{a,s}}} \left( {{\omega \mathord{\left/ {\vphantom {\omega c}} \right. \kern-\nulldelimiterspace} c}} \right)^{2} – k_{x}^{2} .}\)
Figures 7a and 7b show color raster maps with plots of Eq. (3) for the real and imaginary parts of graphene conductivity as the functions of frequency and Fermi energy for τ = 1 ps27. It is seen from Fig. 7a that the real part of the graphene conductivity can be negative at THz frequencies for \(E_{{\text{F}}} > 20\) meV. The imaginary part of the graphene conductivity is positive (inductive) in the entire THz frequency range which is the necessary, as well as sufficient, condition for existing plasmons in graphene28.
Numerical solution of Eq. (2) was performed for the structure with the following.parameters: \(\varepsilon_{{\text{s}}} = 5\) (corresponds to hBN), T = 300 K, τ = 1 ps, \(E_{{\text{F}}} = 50\) meV, d = 10 μm unless different value is specified. These parameter values correspond to a negative real part of the graphene conductivity at THz frequencies yielding the possibility of THz plasmon amplification.