Emission of plasmons by drifting Dirac electrons: where hydrodynamics matters

TL;DR

This paper demonstrates that drifting Dirac electrons in hydrodynamic regimes can emit plasmons via Cerenkov radiation, with the effect characterized through analytical conductivity models and instability growth rates in graphene-based setups.

Contribution

It reveals high-frequency hydrodynamic Cerenkov plasmon emission in Dirac materials and provides analytical tools to characterize this phenomenon across the hydrodynamic-to-ballistic transition.

Findings

Cerenkov emission occurs only in the hydrodynamic regime due to reduced plasmon velocity.

Analytical high-frequency non-local conductivity across the crossover is derived.

Hydrodynamic plasmon instabilities are shown to grow in specific graphene configurations.

Abstract

Direct current in clean semiconductors and metals was recently shown to obey the laws of hydrodynamics in a broad range of temperatures and sample dimensions. However, the determination of frequency window for hydrodynamic phenomena remains challenging. Here, we reveal a phenomenon being a hallmark of high-frequency hydrodynamic transport, the Cerenkov emission of plasmons by drifting Dirac electrons. The effect appears in hydrodynamic regime only due to reduction of plasmon velocity by electron-electron collisions below the velocity of carrier drift. To characterize the Cerenkov effect quantitatively, we analytically find the high-frequency non-local conductivity of drifting Dirac electrons across the hydrodynamic-to-ballistic crossover. We find the growth rates of hydrodynamic plasmon instabilities in two experimentally relevant setups: parallel graphene layers and graphene covered by subwavelength grating, further showing their absence in ballistic regime. We argue that the possibility of Cerenkov emission is linked to singular structure of non-local conductivity of Dirac materials and is independent on specific dielectric environment.

Equations

$$-i\omega \delta f_{\mathbf{p}}+i\mathbf{q}{\mathbf{v}}_{\mathbf{p}}\delta f_{\mathbf{p}}+i\mathbf{q}\delta {\phi }_{\mathbf{q}}{\partial }_{\mathbf{p}}{f}_{\mathbf{p}}^{(0)}={\mathcal{C}}_{ee}\{\delta f_{\mathbf{p}}\},$$

$${\mathcal{C}}_{ee}\{\delta f_{\mathbf{p}}\}=-{\gamma }_{ee}(\delta f_{\mathrm{p}}-\delta f_{\mathrm{hd}}),$$

$$\delta f_{\mathrm{hd}}=\delta \mu {\partial }_{\mu }{f}_{\mathbf{p}}^{(0)}+\delta \mathbf{u}{\partial }_{\mathbf{u}}{f}_{\mathbf{p}}^{(0)}+\delta T{\partial }_T{f}_{\mathbf{p}}^{(0)}.$$

$$\hat{M}=\left(\begin{array}{ccc}1-i{\tilde{\gamma }}_{ee}{I}_{02}& -i{\tilde{\gamma }}_{ee}{\partial }_{\beta }{I}_{02}& 0\\ 0& 1-\frac{2}{3}i{\tilde{\gamma }}_{ee}{\partial }_{\beta }{I}_{13}& {\beta }_0-\frac{2}{3}i{\tilde{\gamma }}_{ee}{I}_{13}\\ 0& \beta -i{\tilde{\gamma }}_{ee}{\partial }_{\beta }{I}_{03}& 1+\frac{{\beta }^2}{2}-i{\tilde{\gamma }}_{ee}{I}_{03}\end{array}\right)$$

$$\delta \mathbf{F}=-2\frac{e\delta \phi }{m{v}_0^2}\left(\begin{array}{c}{I}_{12}-{\beta }_0{I}_{02}\\ I_{23}-{\beta }_0{I}_{13}\\ \frac{3}{2}\left(I_{13}-{\beta }_0{I}_{03}\right)\end{array}\right)$$

$$I_{nm}\left(a,\beta \right)=\frac{{\left(1-{\beta }^2\right)}^{m-\frac{1}{2}}}{2\pi }\underset{0}{\overset{2\pi }{\int }}\frac{{\cos }^n\theta d\theta }{{\left(1-\beta \cos \theta \right)}^m\left(a-\cos \theta \right)}.$$

$$\Delta {\omega }_{\mathrm{bal}}^{\pm }=2q{u}_0(s^2-1)(\sqrt{s^2-1}-s{)}^2,$$

$$\delta f\propto \frac{\cos \theta \delta (\tilde{p}{v}_0-{\epsilon }_F)}{[1-(q{v}_0/\omega )\cos \theta ][1-\beta \cos \theta {]}^2},$$

$$\delta f\propto \frac{\cos \theta \delta (\tilde{p}{v}_0-{\epsilon }_F)}{[1-\beta \cos \theta {]}^3},$$

$${\beta }_{\mathrm{th}}^{\pm }=\frac{v_0}{\sqrt{2}}\sqrt{\frac{q^2{v}_0^2+2{\omega }_p^2(1\pm e^{-qd})}{q^2{v}_0^2+{\omega }_p^2(1\pm e^{-qd})}},$$

$$f_0=[1+e^{({ϵ}_p-{\mathbf{pu}}_0-\mu )/T}{]}^{-1}$$

$${\partial }_{t}n+{\partial }_{\mathbf{r}}(n\mathbf{u})=0,$$

$${\partial }_t(\rho u_i)+{\partial }_{x_j}{\Pi }_{ij}=en{\partial }_{x_i}\phi ,$$

$${\partial }_t\epsilon +{\partial }_{\mathbf{r}}(\rho v_0^2\mathbf{u})-en\mathbf{u}{\partial }_{\mathbf{r}}\phi =0.$$

$$n=\frac{n_{\beta =0}}{[1-{\beta }^2{]}^{3/2}},\rho =\frac{{\rho }_{\beta =0}}{[1-{\beta }^2{]}^{5/2}},$$

$$\epsilon ={\epsilon }_{\beta =0}\frac{1+{\beta }^2/2}{[1-{\beta }^2{]}^{5/2}},$$

$${\rho }_{\beta =0}{v}_0^2=\frac{3}{2}{\epsilon }_{\beta =0}{\Pi }_{xx}=\frac{{\epsilon }_{\beta =0}}{2}\frac{1+2{\beta }^2}{[1-{\beta }^2{]}^{5/2}},{\Pi }_{yy}=\frac{{\epsilon }_{\beta =0}}{2}\frac{1-{\beta }^2}{[1-{\beta }^2{]}^{5/2}}$$

$${\partial }_t(\rho u_x)+\frac{1}{3}{\partial }_x[\rho (v_0^2+2u^2)]=en{\partial }_x\phi ,$$

$$\frac{2}{3}{\partial }_t(\rho [v_0^2+u^2/2])+{\partial }_{\mathbf{r}}(\rho v_0^2\mathbf{u})-en\mathbf{u}{\partial }_{\mathbf{r}}\phi =0.$$

$${\hat{M}}_{\mathrm{hd}}=\left(\begin{array}{ccc}-i(\omega -q{u}_0)& iq& 0\\ 0& -i(\omega -\frac{4}{3}q{u}_0)& -i({\beta }_0\omega -\frac{q}{3})\\ 0& -i({\beta }_0\omega -\frac{3}{2}q)& -i(\omega -\frac{3}{2}q{u}_0)\end{array}\right)$$

$${\delta \mathbf{F}}_{\mathrm{hd}}=\frac{e\delta \phi }{m{v}_0^2}\left(\begin{array}{c}0\\ iq{v}_0\\ \frac{3}{2}iq{u}_0\end{array}\right)$$

$$m=\frac{{\rho }_0}{n_0}\approx \frac{{\epsilon }_F}{v_0^2-u_0^2}.$$

$$\Pi =-\frac{n{q}^2(1-{\beta }_0^2)}{m\left[{\omega }^2(1-\frac{{\beta }_0^2}{2})-\frac{q^2{v}_0^2}{2}(1-2{\beta }^2)-2q{u}_0\omega \right]},$$

$$\sigma =-\frac{i{e}^2\omega }{q^2}\Pi .$$

$$\Pi =-\frac{n{q}^2/m}{{\omega }^2-v_s^2{q}^2},$$

$$\delta f=\frac{i{\nu }_{ee}\delta f_{\mathrm{hd}}-ie\delta \phi \mathbf{q}{\partial }_{\mathbf{p}}{f}_0}{\omega +i{\nu }_{ee}-\mathbf{q}{\mathbf{v}}_{\mathbf{p}}},$$

$$\delta f_{hd}=\delta \mu {\partial }_{\mu }{f}_0+\delta \mathbf{u}{\partial }_{\mathbf{u}}{f}_0+\delta T{\partial }_T{f}_0,$$

$$f_0={\left[1+\exp \left\{\frac{{ϵ}_p-{\mathbf{pu}}_0-\mu }{T}\right\}\right]}^{-1}.$$

$$\sum _{\mathbf{p}}[\delta f-\delta f_{\mathrm{hd}}]=0,\sum _{\mathbf{p}}\mathbf{p}[\delta f-\delta f_{\mathrm{hd}}]=0,\sum _{\mathbf{p}}{ϵ}_{\mathbf{p}}[\delta f-\delta f_{\mathrm{hd}}]=0.$$

$$\sum _{\mathbf{p}}\frac{p^{m-1}{\cos }^n{\theta }_{\mathbf{p}}{f}_0(\mathbf{p})}{a-\cos {\theta }_{\mathbf{p}}},a=\frac{\omega +i{\nu }_{ee}}{q{v}_0}.$$

$$\tilde{p}=p(1-u_0/v_0\cos {\theta }_{\mathbf{p}}),$$

$$\sum _{\mathbf{p}}\frac{p^{m-1}{\cos }^n{\theta }_{\mathbf{p}}{f}_0(\mathbf{p})}{a-\cos {\theta }_{\mathbf{p}}}=\frac{1}{(2\pi {)}^2}{\int }_0^{\infty }d\tilde{p}{\tilde{p}}^m{f}_F(\tilde{p}){\int }_0^{2\pi }\frac{{\cos }^n\theta d\theta }{(a-\cos \theta )(1-\beta \cos \theta {)}^m}.$$

$$J_{nm}=\frac{1}{2\pi }\underset{0}{\overset{2\pi }{\int }}\frac{{\cos }^n\theta d\theta }{(a-\cos \theta )(1-\beta \cos \theta {)}^m}.$$

$$\delta n=\frac{i{\nu }_{ee}}{q{v}_0}\delta [n_{\beta =0}{J}_{02}]+i\frac{e\delta \phi }{q{v}_0}\frac{\partial n_{\beta =0}}{\partial \mu }[q{v}_0{J}_{12}-{\mathbf{qu}}_0{J}_{02}],$$

$$\delta \mathbf{P}{v}_0=\frac{i{\nu }_{ee}}{q{v}_0}\delta [{\epsilon }_{\beta =0}{J}_{13}]+i\frac{e\delta \phi }{q{v}_0}2n_{\beta =0}[q{v}_0{J}_{23}-{\mathbf{qu}}_0{J}_{13}],$$

Full text

The realm of hydrodynamic transport spans at length scales exceeding the particle free path Lifshits and Pitaevskii (1981). Experimental demarcation of hydrodynamics and ballistics is conveniently performed by measuring the flow through a pipe between two reservoirs. The flow through wide pipes is limited by viscosity (Poiseuille flow), while in narrower pipes it is limited by particle injection (Knudsen flow). Recently, similar experiments were performed in ultra-clean solid-state systems, including thin metal wires Moll et al. (2016), Weyl semimetals Gooth et al. (2018), GaAs-based quantum wells Levin et al. (2018), and graphene Bandurin et al. (2016); Kumar et al. (2017); Sulpizio et al. (2019). They have revealed quite a broad window of temperatures and sample dimensions where electrons obey the laws of hydrodynamics Ho et al. (2018) but not ballistics, as thought previously Mayorov et al. (2011).

While the place for dc hydrodynamic (HD) phenomena on temperature and length scales is established Bandurin et al. (2018); Shytov et al. (2018), the bounds for hydrodynamics on frequency scale are less probed Gallagher et al. (2019). Generally, electron-electron (e-e) collisions being the prerequisite of HD transport affect neither dc nor ac conductivity in uniform fields, though they may affect the properties of waves in solids – plasmons. Still, the spectra of plasmons in HD and ballistic regimes are almost identical as they are dictated by long-range Coulomb forces insensitive to microscopic details of e-e interactions Bhatnagar et al. (1954); Abrikosov and Khalatnikov (1959). The character of damping due to e-e scattering in ballistic and HD regimes is different Conti and Vignale (1999); Atwal and Ashcroft (2002), still, it is often masked by extrinsic damping.

In this Letter, we theoretically reveal a plasmonic phenomenon serving as a hallmark of hydrodynamic transport, and is fully prohibited in collisionless ballistic regime. The effect is emission of plasmons by drifting Dirac electrons or, in other words, Cerenkov plasmon instability of electron drift. Our emphasis on Dirac electron systems, especially graphene, is motivated by numerous observations of hydrodynamic phenomena therein Bandurin et al. (2016); Kumar et al. (2017); Sulpizio et al. (2019); Crossno et al. (2016); though some fingerprints of effect can be found in systems of massive electrons.

The possibility of Cerenkov instability in the HD regime is not merely due to reduction of viscous dissipation. It appears due to softening of plasmon velocity by e-e collisions down to the value sufficient to provide phase synchronism between drifting carriers and waves. More precisely, the lower bound on ballistic plasmon velocity in materials with Dirac spectra $\epsilon_{p}=\pm pv_{0}$ is exactly the carrier velocity $v_{0}$ Ryzhii (2006); Lundeberg et al. (2017). The velocity of drift $u_{0}<v_{0}$ thus never satisfies the Cerenkov criterion. In the HD regime, the lower bound on plasmon velocity is only $v_{0}/\sqrt{D}$, where $D$ is the dimension of space Svintsov et al. (2012); Svintsov (2018). The carriers accelerated to drift velocity $u_{0}>v_{0}/\sqrt{D}$ are thus capable of plasmon emission.

Previous studies of current-driven plasmon emission in graphene were a field of delusions, with spurious instabilities in ballistic regime predicted Gumbs et al. (2015); Morgado and Silveirinha (2017). The inadequacy of these predictions stems from breakdown of Galilean invariance in Dirac systems Levitov et al. (2013) which makes the Doppler transform inapplicable for plasmon frequencies in a moving reference frame Svintsov and Ryzhii (2019). More accurate studies Duppen et al. (2016); Wenger et al. (2018); Sabbaghi et al. (2015) revealed no Cerenkov-type instabilities, but were limited to the collisionless case 111The instability obtained in Sabbaghi et al. (2015) is related to current-induced interband population inversion which can be obtained near the charge neutrality point. Here, we consider only the intra-band plasmon emission..

Below we construct the theory of plasmon instabilities in Dirac systems that can handle the subtle issues of Galilean invariance breakdown. Moreover, it is capable to trace the evolution of instabilities across the hydrodynamic-to-ballistic crossover analytically. It is based on solution of kinetic equation with model e-e collision integral satisfying the conservation laws. The obtained conductivity $\sigma(q,\omega)$ of drifting Dirac electrons has a number of unexpected features, including the absence of dissipation at special frequencies and wave vectors satisfying $\omega u_{0}=qv_{0}^{2}$. It is subsequently used as a building block for analysis of current-driven plasmon instabilities in experimentally relevant setups, including parallel graphene layers and graphene covered by sub-wavelength gratings, shown in Fig. 1.

The conductivity $\sigma(q,{\bf\omega})$ of drifting Dirac electron fluid is obtained by solving the kinetic equation for distribution function $f_{\bf p}=f^{(0)}_{\bf p}+\delta f_{\bf p}$ within linear response to the external electric field $\delta{\rm E}_{\bf q}=-i{\bf q}\delta\varphi_{\bf q}e^{i({\bf{qr}}-\omega t)}$:

[TABLE]

here ${\bf v_{p}}=\partial_{\bf p}\epsilon_{\bf p}$ is the quasiparticle velocity. The carrier drift is encoded in zero-order distribution function, which we take in the local-equilibrium (hydrodynamic) form with velocity ${\bf u}_{0}$, Fermi energy $\varepsilon_{F}$ and temperature $T$, $f_{0}=[1+e^{(\epsilon_{p}-{\bf pu}_{0}-\varepsilon_{F})/T}]^{-1}$. We restrict ourselves to collinear propagation of waves and carrier drift.

The crucial step of the solution is the approximation of e-e collision integral $\mathcal{C}_{ee}\{\delta f_{\bf p}\}$ which does not enable any analytical treatment in its original form. We adopt $\mathcal{C}_{ee}\{\delta f_{\bf p}\}$ that pulls all perturbations of distribution function toward a local equilibrium $\delta f_{\rm hd}$ (but not to zero! Bhatnagar et al. (1954); Atwal and Ashcroft (2002)) with a characteristic rate $\gamma_{ee}=\tau_{ee}^{-1}$:

[TABLE]

The main properties of true e-e scattering are encoded into the model, as the distribution modes corresponding to shift of particle number, momentum, and energy are not relaxed. The weights of these modes $\delta\mu$, $\delta{\bf u}$, and $\delta T$ are obtained from respective conservation laws for e-e collisions. These requirements lead us to a linear system of generalized HD equations which can be written symbolically as ${\hat{M}}{\delta\bf x}={\delta\bf F}$. The vector $\delta{\bf x}$ contains unknown hydrodynamic parameters that can be arbitrary linear combinations $\delta\mu$, $\delta{u}$, and $\delta T$; $\delta\bf F$ is the vector of generalized force densities, and $\hat{M}$ is the dynamic matrix. The simplest form is achieved when relative perturbations of particle density $\delta n/n_{0}$, ’relativistic’ velocity $\delta\beta=\delta u/v_{0}$ and mass density $\delta\rho/\rho_{0}$ are treated as unknowns. In this representation, the HD matrix and force vector take the form (Supporting Section II)

[TABLE]

where we have introduced ’relativistic mass’ $m\approx\varepsilon_{F}/v_{0}^{2}$, the inverse Knudsen number $\tilde{\gamma}_{\rm ee}=(qv_{0}\tau_{\rm ee})^{-1}$, and dimensionless functions $I_{nm}(a,\beta)$ of scaled frequency $a=(\omega+i\gamma_{ee})/qv_{0}$ and drift velocity $\beta$:

[TABLE]

The system (4-5) is the central result of this Letter. It provides an explicit expression for high-frequency non-local graphene conductivity $\sigma({\bf q},\omega)$ in the presence of carrier drift across the hydrodynamic-to-ballistic crossover (Supporting section II). It encloses numerous previous studies of graphene ac conductivity as limiting cases Svintsov (2018); Kukhtaruk et al. (2016); Briskot et al. (2015). Classical Navier-Stokes equations along with microscopic expression for viscosity are restored in the HD limit $\tilde{\gamma}_{ee}\gg 1$ by expanding $I_{nm}$. In the opposite ballistic limit, $I_{nm}$ diverge at the boundary of single-particle excitations $\omega\rightarrow qv_{0}$. This divergence translates into singular conductivity at $\omega=qv_{0}$ and absence of plasmon modes below the singular line Lundeberg et al. (2017).

Drift-induced Doppler shift and plasmon undamping. Several non-trivial plasmonic effects appear already in isolated graphene layer in the presence of drift due to the breakdown of Galilean invariance. The latter is readily seen from the generalized hydrodynamic system (4-5) as the wave frequency $\omega$ and drift velocity $u_{0}$ appear therein not only in combination $\omega-qu_{0}$, as it should be for massive electrons.

The first such effect is anomalous Doppler splitting betwеen frequencies of up- and downstream plasmons $\Delta\omega^{\pm}$. It is always below the conventional value of $2qu_{0}$; in the hydrodynamic limit it is exactly one half of it. In the ballistic limit

[TABLE]

where $s=\omega/qv_{0}$ is the ratio of wave phase velocity and Fermi velocity. The ballistic Doppler shift approaches zero as the wave velocity approaches $v_{0}$; it stems from singular ballistic conductivity at $\omega=qv_{0}$.

Much more surprising is the wave damping due to e-e collisions, which is shown on the bottom panel of Fig. 2. The damping of upstream wave continuously increases with the drift speed. The damping of downstream wave for finite $u_{0}$ approaches zero at some peculiar frequencies satisfying $\omega u_{0}=qv_{0}^{2}$, and then continues to grow. Far away from ’undamping point’, the imaginary part of frequency is proportional to $q^{2}$, as it should be for the viscous damping.

The origin of undamping points can be traced back to the excitation of distribution modes $\delta f$ that are insensitive to e-e collisions. This is readily seen in ’boosted coordinates’ $(\tilde{p},\theta)$, where $\tilde{p}=p(1-\beta\cos\theta)$ Moors et al. (2019). In the absence of collisions and at $T/\varepsilon_{F}\ll 1$, electric field excites the distributions

[TABLE]

which do not generally coincide with zero modes of $\mathcal{C}_{ee}$ and are therefore relaxed. But at special points $qv_{0}/\omega=\beta$, the excited distribution coincides with the hydrodynamic momentum mode of $\mathcal{C}_{ee}$:

[TABLE]

It implies that collisions do not have any effect on these modes, and relaxation is absent. We note here that undamping occurs not only for plasmons, but the whole conductivity becomes dissipationless ($\sigma^{\prime}(q,\omega)=0$) at these special frequencies and wave vectors.

Instabilities of counter-streaming flows. Cerenkov-type plasmon instability is most simply achieved in a double layer setup where electron velocities in two layers point in opposite direction. The theory of such instabilities developed for massive electrons Krasheninnikov and Chaplik (1980) was recently erroneously applied to graphene in the ballistic regime Morgado and Silveirinha (2017); Gumbs et al. (2015). Here, we find that the very presence of such instabilities depends strongly on transport regime in the two layers.

When both layers are ballistic, $\omega\tau_{ee}\gg 1$, instabilities are absent despite negative real part of conductivity in the Cerenkov domain, $\omega\lesssim qu_{0}$ Svintsov and Ryzhii (2019). The symmetric (optical) plasmon mode is unaffected by drift, while the frequency of asymmetric one is pulled by current toward lower frequencies. However, it cannot be decreased below the boundary of single-particle excitations $\omega=qv_{0}$ due to singular non-local response of graphene, as given by Eq. 7. As a result, the gain region encloses no plasmon eigenmode.

The acoustic mode, however, readily reaches the region of Cerenkov gain in the hydrodynamic regime, as shown in Fig. 3 A. The non-local dielectric response of graphene in this regime is no more singular, and the wave frequency unimpededly passes through $\omega=qv_{0}$ border. Above the critical velocity, the mode aperiodically growing, i.e. ${\rm Im}\omega>0$, ${\rm Re}\omega=0$. The range of velocities for observation of instabilities in HD regime is located between $\beta^{-}_{\rm th}$ and $\beta^{+}_{\rm th}$,

[TABLE]

where $\omega_{p}=(2\pi ne^{2}|q|/m)^{1/2}$ is the plasma frequency in an isolated graphene layer. $\beta^{-}_{\rm th}$ has a natural lower bound $v_{0}/\sqrt{2}$ that coincides with sound velocity. Interestingly, the threshold velocity weakly depends on carrier density as far as layers are closely bound ($qd\ll 1$).

The growth rate of unstable modes is going down as the e-e collision frequency is reduced, shown in Fig. 3 B. In the weak HD regime, $ql_{\rm ee}\lesssim 1$, the reduced growth rate can be attributed to increased viscous damping. However, outside of the hydrodynamic domain $ql_{\rm ee}\gtrsim 1$, the instabilities do not re-appear as the velocity of acoustic modes is forced to lie above the velocity of carrier drift.

The full stability diagram of counter-streaming double-layer system is calculated in Fig.3 C: the values of drift velocity and e-e collision frequency above the threshold lines correspond to unstable modes. Remarkably, the e-e collision frequency in the dispersion relation appears scaled to $qv_{0}$. It implies that the only parameter governing the transition between HD and ballistic regimes is the Knudsen number $qv_{0}/\gamma_{\rm ee}=ql_{\rm ee}$, where $l_{\rm ee}$ is electron free path with respect to e-e collisions. As a result, instabilities can always be observed in clean systems of sufficiently large length.

Distributed-feedback plasmon lasing. A highly resonant instability leading to electromagnetic emission can be observed in graphene with conducting substrate covered by a metal grating. Such setup is commonly used for spectroscopy of plasmon resonance in 2d electron systems Allen et al. (1977). The reflectance spectrum of such setup is calculated using the formalism of Mikhailov (1998); Matov et al. (1993) with graphene conductivity found from Eqs. (4-5) as a building block.

The electromagnetic response of grating-coupled graphene differs for hydrodynamic and ballistic regimes already in the absence of dc current. Namely, the frequencies of plasmonic dips are reduced by e-e collisions (Fig. 4). When passing direct current in 2DES, the absorption peak is split by Doppler effect Matov et al. (1993); Mikhailov (1998); Wenger et al. (2018). With increasing current, the distinctions between HD and ballistic regimes become more drastic.

The Doppler shift in the ballistic regime is so weak and the singular non-local response at $\omega=qv_{0}$ is so strong that the plasmon frequency is almost unaffected by current (shown in Fig. 4 B). On the contrary, the plasmon frequency in the hydrodynamic regime passes to zero frequency unimpededly (Fig. 4 A). At higher current, the resonant frequency grows again, but the reflection coefficient exceeds unity. The negative absorption is associated with generation of evanescent waves by the grating that fall into the negative conductivity domain of 2DES; the effect is enhanced if frequency corresponds to plasmon resonance.

Both absorption and amplification effects are enhanced if the distance $D$ between 2DES and conducting substrate satisfies the anti-reflection condition, $D=\lambda_{0}/4\sqrt{\kappa}$, where $\kappa$ is the background dielectric constant. The enhancement occurs once the eigen-frequency of Fabry-Perot cavity formed in the vertical direction coincides with the frequency of plasmon. The interaction of plasmon at high current with Fabry-Perot mode leads to divergent reflection coefficient. The divergence implies that such a mode can grow without external stimulus until it is stabilized by nonlinear effects, as it occurs in the distributed feedback lasers.

Discussion and conclusions. We now argue that Cerenkov-type plasmon instability in Dirac materials in HD regime and its absence in ballistic regime are linked to singular structure of conductivity $\sigma(q,\omega)$ and independent of particular dielectric environment. These considerations should be applicable both to 2D and 3D Dirac materials, where conductivity has square-root and log-singularities Lv and Zhang (2013), respectively, at $\omega=qv_{0}$.

The TM plasmon modes of arbitrary structure exist in the domain of positive imaginary part of conductivity ${\rm Im}\sigma(q,\omega)>0$ which, in the ballistic regime, lies above the singularity, $\omega>qv_{0}$. The Cerenkov domain $\omega<qu_{0}$ lies below the singularity. The position of singularity is insensitive to carrier drift. Therefore, one cannot thread the plasmon modes through the singularity by a continuous change of parameter $u_{0}$. This, however, becomes possible in the hydrodynamic regime, where the singularity is removed from real frequency axis due strong to e-e collisions. These arguments do not apply to combinations of Dirac and parabolic-band materials (e.g. graphene parallel to bulk collisionless plasma), where joint plasmon modes can exist at $\omega<qu_{0}$.

The predicted effects can be readily tested experimentally. The anomalous Doppler shifts of plasmons in graphene can be measured with Raman spectroscopy, as it was done for III-V based 2DES Ó Súilleabháin et al. (1994). The plasmon instabilities, both in double layer Gribnikov et al. (2000) and grating-gate setup can result in oscillatory current regimes and emission of terahertz radiation. Such emission can be distinguished from hot-plasmonic emission Tsui et al. (1980) by the presence of well-defined threshold current El Fatimy et al. (2010). Reflectance spectroscopy of grating-gated 2DES is another convenient tool to study Doppler shift and wave amplification Boubanga-Tombet et al. (2018).

Acknowledgement. This work was supported by the grants 18-37-20058/18 and 16-29-03402/18 of the Russian Foundation for Basic Research. The author thanks Denis Fateev for helpful discussions and Mikhail K. Maslov for assistance at the early stage of the work.

I Supporting information

I.1 Hydrodynamic transport regime

The hydrodynamic equations for electron drift velocity $\bf u$, chemical potential $\mu$ and temperature $T$ can be obtained by supplying the local-equilibrium distribution function

[TABLE]

into kinetic equation and integrating it by phase space multiplied by $1$, $\bf p$ and $\epsilon_{\bf p}$. Though a more accurate derivation with explicit account of collisions is possible (see section II), such simplistic derivation may be convenient for analysis of nonlinear effects. The resulting equations can be presented as

[TABLE]

Above, $n$ is the density of electrons, $\rho$ is the equivalent of mass density, $\varepsilon$ is the internal energy density, and $\Pi_{ij}$ is the stress tensor. At given value of chemical potential $\mu$, all these quantities depend on drift velocity $\beta=u/v_{0}$. At the same time, they can be expressed via their values in the absence of drift, i.e. at $\beta=0$:

[TABLE]

To avoid dealing with ’relativistic factors’ of the type $[1-\beta^{2}]^{\alpha}$, it is convenient to consider density $n$, mass density $\rho$ and velocity $\beta$ as hydrodynamic variables. In these variables, the Euler and heat balance equations take on the closed form (assuming one-dimensional motion along $x$-axis)

[TABLE]

The equations can now be linearized to find the conductivity of drifting Dirac electrons in the hydrodynamic regime, $n=n_{0}+\delta ne^{i(kx-\omega t)}$, $u=u_{0}+\delta ue^{i(kx-\omega t)}$, $\rho=\rho_{0}+\delta\rho e^{i(kx-\omega t)}$. This results in the system of equations $\hat{M}_{\rm hd}\delta{\bf x}=\delta{\bf F}_{\rm hd}$, with the matrix and right-hand side given by

[TABLE]

the hydrodynamic mass $m$ si given by

[TABLE]

As a final step, it is possible to solve the system (20-21) for density $\delta n$ and velocity $\delta u$ to find the polarizability and conductivity

[TABLE]

At zero drift velocity, the polarzation of graphene acquires a simple form common for bulk solids:

[TABLE]

where $v_{s}=v_{0}/\sqrt{2}$ is the sound velocity.

I.2 Polarizability and conductivity for drifting electrons at the hydrodynamic-to-ballistic crossover

The solution of kinetic equation with electron-electron collisions reads

[TABLE]

where the perturbed hydrodynamic distribution reads

[TABLE]

and the unperturbed drifting distribution is

[TABLE]

The local-equilibrium Fermi energy $\mu$, drift velocity ${\bf u}_{0}$ and temperature $T$ should be determined from the solution of dc transport equations from the known drain and gate voltages. We leave this solution for further work. The variations of local-equilibrium parameters due to ac field $\delta\mu$, $\delta{\bf u}$ and $\delta T$ should be obtained by requiring particle, momentum and energy conservation upon e-e collisions:

[TABLE]

Substituting the distribution function (26) into conservation laws (29), we can obtain a set of generalized hydrodynamic equations for determination of Fermi energy, drift velocity and temperature. In the course of evaluation, one encounters the following integrals

[TABLE]

Anisotropy of $f_{0}({\bf p})$ introduces complications that can be handled via the change of momentum

[TABLE]

after which the integrals over momentum modulus and angle are decoupled

[TABLE]

Above $f_{F}$ is already an isotropic Fermi function with the same values of $\mu$, $T$ and zero drift velocity. We introduce the notation

[TABLE]

The integral is evaluated by setting $z=e^{i\theta}$ and computing the residues in the poles inside the unit circle $|z|=1$. The remaining integrals over momentum modulus can be expressed via electron density $n$ and energy density $\varepsilon$ at zero drift velocity $\beta\equiv u_{0}/v_{0}$. This results in generalized hydrodynamic system:

[TABLE]

Finally, to justigy the form (4-5) of the main text, we express the quantities $n_{\beta=0}$ and $\varepsilon_{\beta=0}$ through $n$ and $\varepsilon$ at finite drift velocity using the equations of state (15). We also introduce another set of dimensionless integrals

[TABLE]

and inverse Knudsen number $\tilde{\gamma}_{ee}=(qv_{0}\tau_{ee})^{-1}$. Generally, the right-hand side of generalized hydrodynamic equations has the form

[TABLE]

where the ’kinetic’ and ’hydrodynamic’ masses have the form

[TABLE]

In the main text, we neglect the difference between these masses which is justified in the degenerate limit $T/\epsilon_{F}\ll 1$. A small difference between these masses may result in extra plasmon damping due to relaxation of velocity modes by e-e collisions in non-parabolic bands Svintsov (2018).

The formal solution for polarizability is readily derived from generalized hydrodynamic system (4-5). Denoting the elements of hydrodynamic matrix as $M_{ij}$ and components of generalized force vector as $\delta F_{i}$, we find

[TABLE]

I.3 Evaluation of auxiliary spatially dispersive integrals

Upon obtaining the generalized hydrodynamic equations, on encounters the following integrals

[TABLE]

At any given integers $n$ and $m$, they are evaluated by passing to the complex variable $z=e^{i\theta}$. This results in

[TABLE]

The integrand has poles at the points:

[TABLE]

Among these points, $z^{(\beta)}_{-}$ and $z_{0}$ lie inside the unit circle $z=1$, $z^{(a)}_{-}$ lies inside the unit circle for $\rm{Re}a>0$ and $z^{(a)}_{+}$ – for $\rm{Re}a<0$. These statements are independent of sign of imaginary part of $a$.

Evaluation is completed by computation of integrand residues at these poles. As a result, we arrive at the following expressions:

[TABLE]

The remaining necessary integrals can be obtained with recurrence relations

[TABLE]

As apparent from the forms (44), there is a singularity at $a=\rightarrow 1$. As the parameter $a$ has finite imaginary part ($a=(\omega+i\gamma_{ee})/qv_{0}$), the singularity is present only in the ballistic regime (i.e. at $\omega\gg\gamma_{ee}$). Finite strength of e-e collisions softens the singularity.

The integrals also diverge at $\beta\rightarrow 1$, a situation close to that in special relativity. However, this divergence can be re-absorbed into definitions of particle and mass density (15), so that the resulting generalized hydrodynamic equations are free of divergences.

There is a spurious singularity at $a\beta\rightarrow 1$, however, a closer inspection reveals that it is compensated by the zero value of the numerator. An only special property of the computed dielectric response at $a\beta\rightarrow 1$ is the absence of dissipation, as discussed in the main text.

I.4 Analysis of instabilities in the double-layer system

The dispersion law for plasmons in a double-layer structure separated by distance $d$ reads

[TABLE]

Here $\Pi_{+}$ and $\Pi_{-}$ are the polarizabilities of individual top and bottom layers, $V_{0}=2\pi e^{2}/\kappa|q|$ is the Fourier transform of Coulomb interaction and $\kappa$ is the background dielectric constant. Substitution of hydrodynamic polarizability (23) results in biquadratic equation with two eigenmodes

[TABLE]

here we have introduced the dimesionless ’plasmon phase velocity’

[TABLE]

The signs $+$ and $-$ in the absove dispersion can be traced back to optical and acoustic modes of the double-layer structure in the absence of drift, respectively. Indeed, at $\beta=0$ one obtains

[TABLE]

The instability emerges as the frequency of acoustic mode in (49) passes through zero with increasing the drift velocity. This occurs at

[TABLE]

The solution of dispersion relation (48) with polarizability including e-e collisions shows that instability sets on once the acoustic mode frequency crosses zero, i.e. at $\rm Re\omega_{-}=0$. A direct verification of this fact is challenging, but the experience of numerical solutions tells that it is the case. The stability diagram in Fig. 3 (C) was therefore obtained by numerical solution of $\epsilon_{2l}(q,0,\beta)=0$.

We note that the pattern of instabilities can be much richer if the carrier densities in the two layers are non-equal. In this case, the instability does not necessarily set on once the mode frequency crosses zero. A detailed analysis of these cases will be presented elsewhere.

I.5 Diffraction on grating-gated graphene

We consider the diffraction of an electromagnetic wave normally incident on graphene covered by a metal grating. The electric field is polarized along the $x$-axis, i.e. perpendicular to the gratings. The metal grating is assumed to be infinitely thin, its surface conductivity $\sigma_{\rm m}$ exceeds the velocity of light, $\sigma_{\rm m}\gg c$, and is set to infinity in the following calculation. The grating-to-graphene distance $d$ is well below the grating period $a$ and with $W=fa$, where $f$ is the filling factor. This ensures efficient coupling of evanescent waves generated by the grating to the surface plasmons. We also assume that the structure is globally gated with a highly conducting substrate, the distance to back gate $z_{0}$ is the largest length scale in the problem, $z_{0}\gg W\gg d$. We tune $z_{0}$ to accommodate nearly quarter of wavelength in the substrate material, $z_{0}\sim\lambda_{0}/4n_{\rm sub}$, where $\lambda_{0}$ is the free-space wavelength.

Due to uniformity of 2DES is the $x$-direction, the diffraction problem can be formulated on surface current in the grating $j_{s}(x)$, $x\in[0;W]$, while all the information about 2DES is accommodated in the Green’s function of electromagnetic problem. This results in the following integral equation

[TABLE]

where $\mathcal{E}_{0}$ is the field in the grating plane $z=0$ in the absence of grating. It can be presented as

[TABLE]

where $r_{2d}$ is the reflection coefficient of bottom-gated 2DES without grating, and $\mathcal{E}_{\rm inc}$ is the electric field in the incident wave.

The impedance kernel $Z(x-x^{\prime})$ is obtained as follows. First, one finds electric field induced at $z=0$ by $G$-th spatial Fourier harmonic of surface current passing in the grating plane:

[TABLE]

the function $Z_{G}$ is easily obtained by plane-wave matching or transfer-matrix methods. Then $Z(x-x^{\prime})$ is the inverse Fourier transform of $Z_{G}$ with the wave vectors running across the reciprocal wave vectors of the grating, $G_{n}=2\pi n/a$:

[TABLE]

One should distinguish the cases $n=0$ (corresponding to the normally incident propagating wave) and $|n|\geq 1$ (corresponding to evanescent waves generated by grating):

[TABLE]

here $k=\omega/c$ and $k_{1}=\omega\sqrt{\varepsilon_{\rm sub}}/c$ are the wave vectors in vacuum and the substrate material.

It is important to note that impedance kernel (56) with Fourier components (57-58) diverges at large $G$. This is associated with singularities of electric field near the keen edges of thin metal stripe carrying uniform current. The divergence can be cured in two ways. As a first possibility, one can transform Eq. 53 into a second-order differential with respect to $x$ and impose zero boundary conditions on current at the edges

[TABLE]

The Fourier components of modified kernel would have extra two powers of $G$ in the denominator, while the coordinate representation will be non-divergent.

As a second possibility, one can expand the unknown current over the orthogonal basis functions $\phi_{n}(x)$ that already satisfy the zero boundary condition:

[TABLE]

The resulting matrix equation

[TABLE]

would have matrix elements $Z_{nm}$ that quickly converge at large $G$.

We have chosen the basis functions

[TABLE]

that are orthogonal with respect to the inner product:

[TABLE]

In this basis, the elements of the impedance matrix are evaluated analytically at each $G$:

[TABLE]

Expansion of current density over Chebyshev polynomials $U_{n}(x)$ multiplied by their respective weight functions $w(x)$ is advantageous in predicting the character of the field at the edges. However, it comes at the cost of numerical approximation to matrix elements $Z_{nm}$.

In actual calculations, we have truncated the linear system at $n=m=10$, and evaluated the sums up to $G_{\max}=20(2\pi/a)$.