Universal linear and nonlinear electrodynamics of the Dirac fluid

TL;DR

This paper derives a universal relation between linear and nonlinear ac conductivities in Dirac fluids within the hydrodynamic regime, applicable across various systems and conditions, with specific predictions for graphene.

Contribution

It introduces a universal theoretical framework linking linear and second-order nonlinear conductivities in Dirac systems, valid for diverse dispersions and conditions.

Findings

Derived a universal relation for ac conductivities in Dirac fluids.

Showed differences between hydrodynamic and kinetic regimes.

Predicted photon drag and second-harmonic generation in graphene.

Abstract

A general relation is derived between the linear and second-order nonlinear ac conductivities of an electron system in the hydrodynamic regime of frequencies below the interparticle scattering rate. The magnitude and tensorial structure of the hydrodynamic nonlinear conductivity are shown to differ from their counterparts in the more familiar kinetic regime of higher frequencies. Due to universality of the hydrodynamic equations, the obtained formulas are valid for systems with an arbitrary Dirac-like dispersion, ranging from solid-state electron gases to free-space plasmas, either massive or massless, at any temperature, chemical potential or space dimension. Predictions for photon drag and second-harmonic generation in graphene are presented as one application of this theory.

Tables

Table 1. Table 1: Summary for the general case. Notations: D h ( 2 ) = − 1 2 ​ n 3 ​ v 4 W 2 ​ ( 1 − C ise ) subscript superscript 𝐷 2 ℎ 1 2 superscript 𝑛 3 superscript 𝑣 4 superscript 𝑊 2 1 subscript 𝐶 ise D^{(2)}_{h}=-\dfrac{1}{2}\dfrac{n^{3}v^{4}}{W^{2}}(1-C_{\mathrm{ise}}) , D p ( 2 ) = − e 3 ​ n 2 ​ m 2 subscript superscript 𝐷 2 𝑝 superscript 𝑒 3 𝑛 2 superscript 𝑚 2 D^{(2)}_{p}=-\dfrac{e^{3}n}{2m^{2}} , and D 0 ( 2 ) = − g ​ e 3 ​ v 2 32 ​ π ​ ℏ 2 superscript subscript 𝐷 0 2 𝑔 superscript 𝑒 3 superscript 𝑣 2 32 𝜋 superscript Planck-constant-over-2-pi 2 D_{0}^{(2)}=-\dfrac{ge^{3}v^{2}}{32\pi\hbar^{2}} .

Regime/Dispersion	Parabolic band (2D)	graphene
Hydrodynamic
$G_1$	$\frac{D_p^{\left(2\right)}}{{\omega }_1^{+}{\omega }_2^{+}{\omega }_3^{+}}\left(-\frac{i{\Gamma }_d}{{\omega }_1}\right)$	$\frac{D_h^{\left(2\right)}}{{\omega }_1^{+}{\omega }_2^{+}{\omega }_3^{+}}\left(-2\frac{i{\Gamma }_d}{{\omega }_1}\right)$
$G_2$	$\frac{D_p^{\left(2\right)}}{{\omega }_1^{+}{\omega }_2^{+}{\omega }_3^{+}}\left(\frac{{\omega }_3^{+}}{{\omega }_1}\right)$	$\frac{D_h^{\left(2\right)}}{{\omega }_1^{+}{\omega }_2^{+}{\omega }_3^{+}}\left(\frac{{\omega }_3^{+}}{{\omega }_1}\right)$
$G_3$	$\frac{D_p^{\left(2\right)}}{{\omega }_1^{+}{\omega }_2^{+}{\omega }_3^{+}}\left(1+\frac{i{\Gamma }_d}{{\omega }_1}+\frac{\partial P}{\partial n_E}\frac{2i{\Gamma }_d}{\omega +i{\Gamma }_E}\right)$	$\frac{D_h^{\left(2\right)}}{{\omega }_1^{+}{\omega }_2^{+}{\omega }_3^{+}}\left[1+2\left(\frac{i{\Gamma }_d}{{\omega }_1}+\frac{\partial P}{\partial n_E}\frac{2i{\Gamma }_d}{\omega +i{\Gamma }_E}\right)\right]$
$T_1$	$-\frac{8D_p^{\left(2\right)}}{\omega \left({\omega }^2+{\Gamma }_d^2\right)}$	$-\frac{12D_h^{\left(2\right)}}{\omega \left({\omega }^2+{\Gamma }_d^2\right)}$
$T_2$	$0$	$\frac{4D_h^{\left(2\right)}}{\omega \left({\omega }^2+{\Gamma }_d^2\right)}$
${\tilde{T}}_1$	$0$	$0$
Kinetic nonconserving
$G_1$	$\frac{D_p^{\left(2\right)}}{{\omega }_1^{+}{\omega }_2^{+}{\omega }_3^{+}}\left(\frac{{\omega }_1^{+}+{\omega }_2^{+}}{{\omega }_3^{+}}-\frac{{\omega }_1^{+}}{{\omega }_1}\right)$	$\frac{2D_0^{\left(2\right)}}{{\omega }_1^{+}{\omega }_2^{+}{\omega }_3^{+}}\left(-1-\frac{{\omega }_2^{+}}{2{\omega }_1^{+}}+\frac{i\Gamma }{{\omega }_3^{+}}-\frac{2i\Gamma }{{\omega }_1}\right)$
$G_2$	$\frac{D_p^{\left(2\right)}}{{\omega }_1^{+}{\omega }_2^{+}{\omega }_3^{+}}\left(\frac{{\omega }_1^{+}+{\omega }_2^{+}}{{\omega }_3^{+}}+\frac{{\omega }_2^{+}}{{\omega }_1^{+}}\right)$	$\frac{2D_0^{\left(2\right)}}{{\omega }_1^{+}{\omega }_2^{+}{\omega }_3^{+}}\left(1+\frac{3}{2}\frac{{\omega }_2^{+}}{{\omega }_1^{+}}+\frac{i\Gamma }{{\omega }_3^{+}}\right)$
$G_3$	$\frac{D_p^{\left(2\right)}}{{\omega }_1^{+}{\omega }_2^{+}{\omega }_3^{+}}\frac{{\omega }_1^{+}}{{\omega }_1}$	$\frac{2D_0^{\left(2\right)}}{{\omega }_1^{+}{\omega }_2^{+}{\omega }_3^{+}}\left(1-\frac{{\omega }_2^{+}}{2{\omega }_1^{+}}-\frac{i\Gamma }{{\omega }_3^{+}}+\frac{2i\Gamma }{{\omega }_1}\right)$
$T_1$	$4D_p^{\left(2\right)}\frac{-3{\omega }^2-{\Gamma }^2}{\omega {\left({\omega }^2+{\Gamma }^2\right)}^2}$	$16D_0^{\left(2\right)}\frac{-{\omega }^2-{\Gamma }^2}{\omega {\left({\omega }^2+{\Gamma }^2\right)}^2}$
$T_2$	$4D_p^{\left(2\right)}\frac{-{\omega }^2+{\Gamma }^2}{\omega {\left({\omega }^2+{\Gamma }^2\right)}^2}$	$16D_0^{\left(2\right)}\frac{{\Gamma }^2}{\omega {\left({\omega }^2+{\Gamma }^2\right)}^2}$
${\tilde{T}}_1$	$8D_p^{\left(2\right)}\frac{\Gamma }{{\left({\omega }^2+{\Gamma }^2\right)}^2}$	$32D_0^{\left(2\right)}\frac{\Gamma }{{\left({\omega }^2+{\Gamma }^2\right)}^2}$
Kinetic
$G_1$		$\frac{2D_0^{\left(2\right)}}{{\omega }_1^{\left(1\right)}{\omega }_2^{\left(1\right)}{\omega }_3^{\left(1\right)}}\{-\frac{{\omega }_2^{\left(1\right)}}{2{\omega }_1^{\left(2\right)}}(1+{\epsilon }_F\frac{{\partial }_{\epsilon }i{\Gamma }_1}{{\omega }_3^{\left(1\right)}})-\frac{2{\omega }_1^{\left(1\right)}}{{\omega }_1}$
		$+\frac{{\omega }_1^{\left(1\right)}+{\omega }_2^{\left(1\right)}}{{\omega }_3^{\left(2\right)}}[1-\frac{1}{2}{\epsilon }_F(\frac{{\partial }_{\epsilon }i{\Gamma }_1}{{\omega }_3^{\left(1\right)}}+\frac{{\partial }_{\epsilon }i{\Gamma }_2}{{\omega }_3^{\left(2\right)}})]\}$
$G_2$		$\frac{2D_0^{\left(2\right)}}{{\omega }_1^{\left(1\right)}{\omega }_2^{\left(1\right)}{\omega }_3^{\left(1\right)}}\{\frac{{\omega }_2^{\left(1\right)}}{{\omega }_1}(1-{\epsilon }_F\frac{{\partial }_{\epsilon }i{\Gamma }_1}{{\omega }_3^{\left(1\right)}})+$
		$\frac{{\omega }_2^{\left(1\right)}}{2{\omega }_1^{\left(2\right)}}\left(1+{\epsilon }_F\frac{{\partial }_{\epsilon }i{\Gamma }_1}{{\omega }_3^{\left(1\right)}}\right)+$
		$\frac{{\omega }_1^{\left(1\right)}+{\omega }_2^{\left(1\right)}}{{\omega }_3^{\left(2\right)}}[1-\frac{1}{2}{\epsilon }_F(\frac{{\partial }_{\epsilon }i{\Gamma }_1}{{\omega }_3^{\left(1\right)}}+\frac{{\partial }_{\epsilon }i{\Gamma }_2}{{\omega }_3^{\left(2\right)}})]\}$
$G_3$		$\frac{2D_0^{\left(2\right)}}{{\omega }_1^{\left(1\right)}{\omega }_2^{\left(1\right)}{\omega }_3^{\left(1\right)}}(2\frac{{\omega }_1^{\left(1\right)}}{{\omega }_1}-\frac{{\omega }_1^{\left(1\right)}+{\omega }_2^{\left(1\right)}}{{\omega }_3^{\left(0\right)}}{\epsilon }_F\frac{{\partial }_{\epsilon }i{\Gamma }_1}{{\omega }_3^{\left(1\right)}}$
		$-\frac{{\omega }_2^{\left(1\right)}}{2{\omega }_1^{\left(2\right)}}\left(1+{\epsilon }_F\frac{{\partial }_{\epsilon }i{\Gamma }_1}{{\omega }_3^{\left(1\right)}}\right)$
		$-\frac{{\omega }_1^{\left(1\right)}+{\omega }_2^{\left(1\right)}}{{\omega }_3^{\left(2\right)}}[1-\frac{1}{2}{\epsilon }_F(\frac{{\partial }_{\epsilon }i{\Gamma }_1}{{\omega }_3^{\left(1\right)}}+\frac{{\partial }_{\epsilon }i{\Gamma }_2}{{\omega }_3^{\left(2\right)}})]\}$
$T_1$		$8D_0^{\left(2\right)}\frac{\left({\omega }^2+{\Gamma }_2^2\right)\left(-3+{\epsilon }_F{\partial }_{\epsilon }\ln {\Gamma }_1\right)+{\omega }^2\left(1+{\Gamma }_2/{\Gamma }_1\right)\left(1+{\epsilon }_F{\partial }_{\epsilon }\ln {\Gamma }_1\right)}{\omega \left({\omega }^2+{\Gamma }_1^2\right)\left({\omega }^2+{\Gamma }_2^2\right)}$
$T_2$		$8D_0^{\left(2\right)}\frac{1+{\epsilon }_F{\partial }_{\epsilon }\ln {\Gamma }_1}{\omega \left({\omega }^2+{\Gamma }_1^2\right)}$
${\tilde{T}}_1$		$-4D_0^{\left(2\right)}\frac{\left(1+{\Gamma }_2/{\Gamma }_1\right){\Gamma }_2}{\left({\omega }^2+{\Gamma }_1^2\right)\left({\omega }^2+{\Gamma }_2^2\right)}\left(1+{\epsilon }_F{\partial }_{\epsilon }\ln {\Gamma }_1\right)$
Quantum	Ref. Stolz,1967	Refs. Cheng et al.,2017; Wang et al.,2016

Table 2. Table 2: Clean limit: Γ d = Γ 1 → 0 subscript Γ 𝑑 subscript Γ 1 → 0 \Gamma_{d}=\Gamma_{1}\to 0 , Γ 2 → 0 → subscript Γ 2 0 \Gamma_{2}\to 0 , Γ 2 / Γ 1 = 2 = const subscript Γ 2 subscript Γ 1 2 const \Gamma_{2}/\Gamma_{1}=2=\mathrm{const} .

Regime/Dispersion	General	Parabolic band (2D)	graphene
Hydrodynamic
$G_1$	$0$	$0$	$0$
$G_2$	$\frac{D_h^{\left(2\right)}}{{\omega }_1{\omega }_2{\omega }_3}\left(\frac{{\omega }_3}{{\omega }_1}\right)$	$\frac{D_p^{\left(2\right)}}{{\omega }_1{\omega }_2{\omega }_3}\left(\frac{{\omega }_3}{{\omega }_1}\right)$	$\frac{D_h^{\left(2\right)}}{{\omega }_1{\omega }_2{\omega }_3}\left(\frac{{\omega }_3}{{\omega }_1}\right)$
$G_3$	$\frac{D_h^{\left(2\right)}}{{\omega }_1{\omega }_2{\omega }_3}$	$\frac{D_p^{\left(2\right)}}{{\omega }_1{\omega }_2{\omega }_3}$	$\frac{D_h^{\left(2\right)}}{{\omega }_1{\omega }_2{\omega }_3}$
$T_1$	$4D_h^{\left(2\right)}\left(-\frac{1}{1-C_{\mathrm{ise}}}-1\right)\frac{1}{{\omega }^3}$	$-8D_p^{\left(2\right)}\frac{1}{{\omega }^3}$	$-12D_h^{\left(2\right)}\frac{1}{{\omega }^3}$
$T_2$	$4D_h^{\left(2\right)}\left(\frac{1}{1-C_{\mathrm{ise}}}-1\right)\frac{1}{{\omega }^3}$	$0$	$4D_h^{\left(2\right)}\frac{1}{{\omega }^3}$
${\tilde{T}}_1$	$0$	$0$	$0$
Kinetic
$G_1$		$0$	$\frac{D_0^{\left(2\right)}}{{\omega }_1{\omega }_2{\omega }_3}\left(-2-\frac{{\omega }_2}{{\omega }_1}\right)$
$G_2$		$\frac{D_p^{\left(2\right)}}{{\omega }_1{\omega }_2{\omega }_3}\frac{{\omega }_3}{{\omega }_1}$	$\frac{D_0^{\left(2\right)}}{{\omega }_1{\omega }_2{\omega }_3}\left(2+3\frac{{\omega }_2}{{\omega }_1}\right)$
$G_3$		$\frac{D_p^{\left(2\right)}}{{\omega }_1{\omega }_2{\omega }_3}$	$\frac{D_0^{\left(2\right)}}{{\omega }_1{\omega }_2{\omega }_3}\left(2-\frac{{\omega }_2}{{\omega }_1}\right)$
$T_1$		$-12D_p^{\left(2\right)}\frac{1}{{\omega }^3}$	$32D_0^{\left(2\right)}\frac{1}{{\omega }^3}$
$T_2$		$-4D_p^{\left(2\right)}\frac{1}{{\omega }^3}$	$16D_0^{\left(2\right)}\frac{1}{{\omega }^3}$
${\tilde{T}}_1$	$0$	$0$	$0$
Quantum			Refs. Cheng et al.,2017; Wang et al.,2016

Equations

$$\mathbf{E}(\mathbf{r},t)=\mathbf{E}({\mathbf{q}}_1,{\omega }_1)e^{i{\mathbf{q}}_1\mathbf{r}-i{\omega }_{1}t}+\mathbf{E}({\mathbf{q}}_2,{\omega }_2)e^{i{\mathbf{q}}_2\mathbf{r}-i{\omega }_{2}t}+c.c.$$

$${\sigma }_{ilm}^{(2)}=\frac{D_h^{(2)}}{{\omega }_1{\omega }_2{\omega }_3}\left(\frac{{\omega }_3}{{\omega }_1}{q}_{1l}{\delta }_{im}+q_{1i}{\delta }_{lm}\right)+\left(\begin{array}{ccc}1& \leftrightarrow & 2\\ l& \leftrightarrow & m\end{array}\right),$$

$$D_h^{(2)}=-\frac{1}{4{\pi }^2}{\left(\frac{\partial (D_h{)}^2}{\partial \rho }\right)}_{\mathrm{ise}}$$

$$D_h=\pi \frac{e^{2}n}{m^{\ast }},n\equiv \frac{\rho }{e}.$$

$$\sigma (q=0,\omega )=\frac{1}{\pi }\frac{D_h}{-i\omega +{\Gamma }_d}+{\sigma }_0,\omega \ll {\Gamma }_{ee}.$$

$$\begin{array}{rl}{F}_i^{(2)}(\mathbf{r})& =\frac{i}{2}{\partial }_i\left[\frac{\mathbf{E}(\mathbf{r},{\omega }_1)}{{\omega }_2}\frac{\partial \sigma ({\omega }_2)}{\partial \rho }\mathbf{E}(\mathbf{r},{\omega }_2)+(1\leftrightarrow 2)\right]\\ & =-\frac{i}{2}{\partial }_i\left[\frac{\mathbf{E}(\mathbf{r},{\omega }_1)}{{\omega }_1}\frac{\partial \sigma ({\omega }_2)}{\partial \rho }\mathbf{E}(\mathbf{r},{\omega }_2)+(1\leftrightarrow 2)\right].\end{array}$$

$${\partial }_t\rho +\mathrm{div}1\mathbf{j}=0,$$

$${\partial }_t{n}_E+\mathrm{div}1\left({\gamma }^{2}W\mathbf{u}\right)=\mathbf{j}\cdot \mathbf{E},$$

$$({\gamma }^2{m}^{\ast }n){\mathcal{D}}_t\mathbf{u}=\rho {\mathbf{F}}_L-\frac{\mathbf{u}}{v^2}\left(\mathbf{j}\cdot \mathbf{E}\right)-\mathcal{D}P,$$

$$m^{\ast }=\frac{W}{n{v}^2},{\mathbf{F}}_L=\mathbf{E}+\frac{1}{c}\mathbf{u}\times \mathbf{B},\mathrm{curl}\mathbf{E}=-\frac{1}{c}{\partial }_t\mathbf{B}.$$

$$\sigma (q=0,\omega )$$

$$D_h^{(2)}=-\frac{1}{2}\frac{e^{3}n}{m^{\ast 2}}(1-C_{\mathrm{ise}}),$$

$$C_{\mathrm{ise}}=\frac{n}{W}{\left(\frac{\partial P}{\partial n}\right)}_{\mathrm{ise}}=\frac{1}{m^{\ast }{v}^2}{\left(\frac{\partial P}{\partial n}\right)}_{\mathrm{ise}}$$

$$D_0^{(2)}=-\frac{g}{32\pi }\frac{e^3{v}^2}{{\hslash }^2}.$$

$${\sigma }_{ilm}^{(2)}=\text{mbox}$$

$$\begin{array}{rl}{\Sigma }_{ilmn}=& -\left(1-3\frac{{\omega }_3}{{\omega }_1}\right){\delta }_{im}{\delta }_{nl}-\left(1+\frac{{\omega }_3}{{\omega }_1}\right){\delta }_{il}{\delta }_{nm}\\ & \text{mbox}+\left(3-\frac{{\omega }_3}{{\omega }_1}\right){\delta }_{in}{\delta }_{lm}.\end{array}$$

$$\frac{j_x}{C_j}$$

$$\frac{j_y}{C_j}$$

$$T_1=-3T_2,T_2=\frac{4D_h^{(2)}}{\omega \left({\omega }^2+{\Gamma }_d^2\right)},{\tilde{T}}_1=0.$$

$${\tilde{T}}_1=-\frac{48D_0^{(2)}{\Gamma }_d}{\left({\omega }^2+{\Gamma }_d^2\right)\left({\omega }^2+4{\Gamma }_d^2\right)}.$$

$$\sigma (0,\omega )=\frac{D_h/\pi }{-i\omega +{\Gamma }_d}+{\sigma }_0,$$

$$D_h=\frac{\pi e^{2}n}{m^{\ast }}.$$

$$D_h(T=0)=\frac{g}{4}\frac{e^2{k}_{F}v}{\hslash },$$

$${\sigma }_{ij}\left(\mathbf{q},\omega \right)={\sigma }_L(q,\omega )\frac{q_i{q}_j}{q^2}+{\sigma }_T(q,\omega )\left({\delta }_{ij}-\frac{q_i{q}_j}{q^2}\right).$$

$${\sigma }_L(q,\omega )=\frac{D_h/\pi }{-i\omega +{\Gamma }_d-\frac{v_d^2{q}^2}{\omega }}.$$

$$q=\frac{i\kappa \omega }{2\pi {\sigma }_L(q,\omega )}.$$

$$C_{\mathrm{ise}}=\frac{n}{W}\frac{\hslash k_F{v}_F}{d},$$

$$C_{\mathrm{ise}}=\frac{1}{d}\frac{\hslash k_F{v}_F}{{\epsilon }_F}.$$

$$C_{\mathrm{ise}}(T=0)=\frac{1}{d}\frac{v_F^2}{v^2}$$

$$C_{\mathrm{ise}}=\frac{1}{2}$$

$$\mathbf{j}=e\sum {\mathbf{v}}_{\mathbf{p}}{f}_{\mathbf{p}}(t)$$

$$f_{\mathbf{p}}-f_{\mathbf{p}}^{(0)}=a_v(t)g_v+a_p(t)g_p.$$

$$\left({\partial }_t+e\mathbf{E}{\partial }_{\mathbf{p}}\right)f_{\mathbf{p}}=-{\partial }_t{a}_v{g}_v+{\partial }_t{a}_p{g}_p+eE{g}_v=\hat{I}\left[f_{\mathbf{p}}\right].$$

Full text

There has been a renewed interest to hydrodynamic phenomena in electron systems with a Dirac-like energy-momentum dispersion $\varepsilon_{p}^{2}=(pv)^{2}+(mv^{2})^{2}$. This subject was revived by studies in quantum criticality and holographic field theory Hartnoll et al. (2007) and outspread in research on two-dimensional (2D) conductors, e.g., graphene where the massless dispersion $m=0$ is realized de Jong and Molenkamp (1995); Müller et al. (2008, 2009); Briskot et al. (2015); Narozhny et al. (2015); Principi and Vignale (2015a, b); Bandurin et al. (2016); Crossno et al. (2016); Moll et al. (2016); Sun et al. (2016); Lucas et al. (2016); Guo et al. (2017). An experimental observation of a viscous electron flow in graphene has been recently reported Bandurin et al. (2016). Although not uncommon in plasmas Tsytovich (1970), this type of transport is highly unusual in a solids. It may be possible only in a limited range of temperatures $T$ and chemical potentials $\mu$ in pristine samples where the combined rate of electron-impurity (ei) and electrons-phonon (ep) scattering $\Gamma_{d}=\Gamma_{ei}+\Gamma_{ep}$ is lower than the momentum-conserving electron-electron (ee) scattering rate $\Gamma_{ee}$ Gurzhi (1968); Andreev et al. (2011). The respective mean-free paths must obey the inequality $l_{d}\gg l_{ee}$. Under these conditions, the electron dynamics at frequencies $\omega\ll\Gamma_{ee}$ and momenta $q\ll l_{ee}^{-1}$ is governed by collective variables that obey hydrodynamic equations Landau and Lifshitz (1987). The frequency range $\Gamma_{d}<\omega<\Gamma_{ee}$ may be as wide as several THz in graphene (see below). Therefore exploring electrodynamics of Dirac fluids may be worthwhile.

In this Letter we focus on the second-order ac conductivity, which controls nonlinear optical phenomena such as sum (difference) frequency generation and also photon drag. In the Supplemental material SM , we also discuss the third-order conductivity important for the Kerr effect. Prior work Glazov (2011); Mikhailov (2011); Yao et al. (2014); Mikhailov (2016); Wang et al. (2016); Tokman et al. (2016); Cheng et al. (2017); Manzoni et al. (2015); Rostami et al. (2017) has indicated that in graphene such effects may or may not Khurgin (2014) be stronger than in typical metals and semiconductors. We find significant differences of our results from what one obtains at frequencies $\omega>\Gamma_{ee}$ where the dynamics is described by the Boltzmann kinetic equation. The still higher frequency quantum regime (Fig. 1) is beyond the scope of our investigation.

Recall that the second-order conductivity is a third-rank tensor $\sigma^{(2)}_{ilm}\left(\mathbf{q}_{1},\omega_{1};\mathbf{q}_{2},\omega_{2}\right)$ that describes the current of frequency $\omega_{3}=\omega_{1}+\omega_{2}$ and momentum $\mathbf{q}_{3}=\mathbf{q}_{1}+\mathbf{q}_{2}$ generated, to the order ${O}\left(\mathbf{E}^{2}\right)$, in response to an electric field

[TABLE]

By convention, $\sigma^{(2)}_{ilm}$ is symmetrized, i.e., invariant under the interchange $(1\leftrightarrow 2,l\leftrightarrow m)$. If the system preserves parity, which we assume to be the case, $\sigma^{(2)}_{ilm}$ must vanish if both $\mathbf{q}_{\nu}$, $\nu=1,2$ are zero. At small $\mathbf{q}_{\nu}$, relevant for optical/THz experiments, $\sigma^{(2)}_{ilm}$ should scale linearly with $\mathbf{q}_{\nu}$. In comparison, dissipative effects due to viscosity and heat conduction Andreev et al. (2011); Forcella et al. (2014), which scale as $|\mathbf{q}_{\nu}|^{2}$, are subleading: the fluid dynamics is approximately isentropic (ise) Landau and Lifshitz (1987). Below we show that in this regime the second-order conductivity has the universal form

[TABLE]

for an arbitrary mass $m$, equilibrium charge density $\rho$, temperature $T$, and space dimension $d$. All the material-specific parameters are contained in the second-order spectral weight $D_{h}^{(2)}$, which we find to be equal to the derivative

[TABLE]

of the squared linear-response (i.e., Drude) spectral weight

[TABLE]

As stated above, these formulas hold for either massless or massive electrons. Conventional metals and semiconductors have a parabolic dispersion. This case is exemplified by the nonrelativistic limit $|\mu|,T\ll mv^{2}$ of our equations, yielding $m^{\ast}=m$. This result can also be understood as the consequence of Galilean invariance, which demands that the ee interactions affect the linear and nonlinear conductivities only in higher orders in $\mathbf{q}_{\nu}$. This is why the effective mass $m^{\ast}$ in Eq. (S2) is equal to the bare mass $m$ and the leading $\mathbf{q}_{\nu}$-linear terms of $\sigma^{(2)}_{ilm}$ [Eq. (2)] are the same in hydrodynamic Tsytovich (1970), kinetic Aliev et al. (1992), and quantum Stolz (1967) domains. The equality of $m^{*}$ and $m$ does not hold if either $|\mu|$ or $T$ are comparable or larger than the energy gap $2mv^{2}$, e.g., in the case of graphene. (For linear conductivity at $T=0$ this has been discussed at length Kotov et al. (2012); Basov et al. (2014); Link et al. (2016).)

In the hydrodynamic regime of graphene, frequent collisions force electrons and holes to move together, causing cancellation of their partial currents. This enhances $m^{\ast}$ and reduces $D_{h}$ below its kinetic counterpart $D_{k}$ at all $T>0$, see Fig. 2(a). Similarly, $D_{h}^{(2)}$ decreases with $T$ at fixed $\rho$ much faster in the hydrodynamic regime than in the previously studied kinetic one, see Fig. 2(b).

Let us now present a qualitative argument for Eq. (2). Consider the expansion of a given Fourier harmonic of the electric current $j_{i}(\mathbf{q},\omega)=j_{i}^{(1)}(\mathbf{q},\omega)+j_{i}^{(2)}(\mathbf{q},\omega)+\ldots$ in power series of the driving electric field $\mathbf{E}$. The first term is given by $j_{i}^{(1)}(\mathbf{q},\omega)=\sigma_{ij}(\mathbf{q},\omega)F_{i}^{(1)}(\mathbf{q},\omega)$ where $F_{i}^{(1)}(\mathbf{q},\omega)=E_{i}(\mathbf{q},\omega)+\mathcal{O}(\mathbf{q})$ is the driving force per unit charge and $\sigma_{ij}$ is the linear-response conductivity tensor. It suffices to consider the $\mathbf{q}\to 0$ limit in which $F_{i}^{(1)}\to E_{i}$, $\sigma_{ij}\to\delta_{ij}\sigma$. The scalar $\sigma$ can be in general separated into the Drude pole and a nonsingular correction $\sigma_{0}$ (to be discussed below):

[TABLE]

Next, to the second order we expect $j_{i}^{(2)}=\sigma^{(1)}F_{i}^{(1)}+\sigma F_{i}^{(2)}$. Here $\sigma^{(1)}=(\partial\sigma/\partial\rho)\rho^{(1)}$ and $\rho^{(1)}$ are the perturbations of the conductivity and charge density. The latter perturbation can be found from the continuity equation (S23a), which gives $\rho^{(1)}(\mathbf{q},\omega)=\mathbf{q}\cdot\mathbf{j}^{(1)}(\mathbf{q},\omega)/\omega$. Calculation of the second-order driving force $F_{i}^{(2)}$ is the difficult part of the problem. We glean the answer from the $\omega_{1}\simeq-\omega_{2}\gg\omega_{3}$ case where it is equal to the sum of the pondermotive and Abraham forces Landau and Lifshitz (1984). The former is of order $(\omega_{3})^{0}$, the latter is the leading $(\omega_{3})^{1}$ correction. Following Landau and Lifshitz (1984), Sec. 81, we find the real-space representation of the pondermotive force to be

[TABLE]

The replacement of $\omega_{2}$ by $-\omega_{1}$ in the second line cannot be strictly justified if $\omega_{3}\neq 0$. However, it is a natural way to ensure the triangular permutation symmetry of $\sigma^{(2)}_{ilm}$, which follows from the energy conservation Il’inskii and Keldysh (1994) in the dissipationless limit $\Gamma_{d},\sigma_{0}\to 0$. Assembling all the terms of $j_{i}^{(2)}$, we can read off $\sigma^{(2)}_{ilm}$ and see it coincides with Eq. (2). One can verify that for a nonrelativistic electron gas our formulas agree with those in literature Tsytovich (1970); Aliev et al. (1992).

The case of a Lorentz-invariant Dirac fluid can be studied rigorously. Proposed solid-state examples of such fluids Hartnoll et al. (2007) actually lack true Lorentz invariance. Their matter and field components have different limiting velocities, $v$ and $c$. However, if Coulomb interactions are weak, the approximate Lorentz invariance with velocity $v$ holds. In graphene this is so if the dielectric constant $\kappa$ of the environment is large, so that the interaction constant $e^{2}/(\hbar\kappa v)$ is small. We will use relativistic hydrodynamics to derive $D_{h}$ and $D_{h}^{(2)}$ for this model and verify our key result (2).

Let us introduce two additional quanitites. One is the flow velocity $\mathbf{u}$ that defines the electric current $\mathbf{j}=\rho\mathbf{u}$. The other is the energy density $n_{E}=\gamma^{2}W-P$ related to the pressure $P=P(\mu,T)$ and enthalpy density $W=W(\mu,T)$ at thermal equilibrium, $u=0$ Landau and Lifshitz (1987). Here $\gamma\equiv{1}/\sqrt{1-(u/v)^{2}}$ and $n_{E}$ is referenced to the $\mu=T=u=0$ state. Relativistic hydrodynamic equations admit many equivalent formulations Landau and Lifshitz (1987); Müller et al. (2008, 2009); Kovtun (2012); Briskot et al. (2015), e.g.,

[TABLE]

The first pair is the charge continuity equation and the energy conservation equation sans the subleading viscous and thermal conductivity terms. Equation (7d) for the Lorentz force $\mathbf{F}_{L}$ includes the force from the ac magnetic field $\mathbf{B}$ induced by $\mathbf{E}$. (We assume that no static magnetic field is present.) This term is important if $\mathbf{E}$-field has a transverse component. Equation (S23c) is the relativistic Euler equation written in “covariant derivatives” $\mathcal{D}_{i}=\partial_{i}+({u}_{i}/{v^{2}})\partial_{t}$, $\mathcal{D}_{t}=\partial_{t}+\Gamma_{d}+{u}_{i}\partial_{i}$ with the scattering rate $\Gamma_{d}$ accounting for momentum dissipation. We solve these equations for $\mathbf{j}$ perturbatively in $\mathbf{E}$ to get the desired conductivities.

The linear response has already been treated at length Hartnoll et al. (2007); Müller et al. (2008, 2009); Kovtun (2012); Briskot et al. (2015); Narozhny et al. (2015); Sun et al. (2016); Lucas et al. (2016). For massless particles, $W=(d+1)P=\frac{d+1}{d}n_{E}\propto T^{d+1}$. The hydrodynamic Drude weight $D_{h}=\pi\rho^{2}v^{2}/W$ [cf. Eqs. (S2) and (7d)] decreases as $T^{-d-1}$ with $T$, i.e., as $T^{-3}$ in graphene. The usual, kinetic Drude weight $D_{k}(\mu,T)=(g/2)({e}/{\hbar})^{2}T\ln\left[2\cosh({\mu}/{2T})\right]$ where $g=4$ is the number of Dirac cones Basov et al. (2014) behaves differently. After some initial drop, $D_{k}$ increases with $T$ because of thermal excitation of carriers, see Fig. 2(a). The question how the opposite trends of $D_{h}$ and $D_{k}$ could be reconciled has not been given proper attention in prior literature. As a tentative answer, we suggest the interpolation formula:

[TABLE]

This formula can be derived from the Boltzmann kinetic equation with the ee scattering rate $\Gamma_{ee}$ added to the collision integral SM . Matching it with Eq. (5) at $\omega\ll\Gamma_{ee}$, we deduce the parameter $\sigma_{0}=(D_{k}-D_{h})/(\pi\Gamma_{ee})$ therein (we assume $\Gamma_{d}\ll\Gamma_{ee}$) Briskot et al. (2015). According to Eq. (S17), the effective Drude weight $\pi\omega\operatorname{Im}\sigma$ as a function of $\omega$ exhibits two plateaus, see Fig. 3, and $\sigma$ as a function of $T$ at fixed $\omega$ may look like as sketched in Fig. 2(a). A quantitative theory of these crossover behaviors is a challenge for future work. Meanwhile, Fig. 3 indicates the existence of two separate frequency intervals where $\operatorname{Im}\sigma(\omega)\gg\operatorname{Re}\sigma(\omega)$. In these intervals weakly damped collective modes are possible: sound waves Kovtun (2012) (or energy waves Phan et al. or “demons” Sun et al. (2016) ) in the hydrodynamic regime and plasmons in the kinetic one, see also Fig. 1.

Let us move on to the second-order conductivity, ignoring the momentum dissipation for now, $\Gamma_{d}\to 0$. In the hydrodynamic regime we have two ways to derive $\sigma^{(2)}_{ilm}$. The quick one is via Eq. (2). The only unknown parameter is $D_{h}^{(2)}$, which we can calculate from Eq. (3) applied to $D_{h}=\pi\rho^{2}v^{2}/W$. This yields

[TABLE]

where

[TABLE]

is the dimensionless isentropic bulk modulus. Note that for massless electrons $C_{\mathrm{ise}}=1/d$. The second derivation we can do is from hydrodynamic Eqs. (S23), which is more tedious SM but gives the same result. This verifies the validity of our universal formula (2) for Dirac fluids.

Let us examine the $T$-dependence of the spectral weight $D_{h}^{(2)}$. As one can anticipate, $D_{h}^{(2)}$ rapidly decreases at high $T$, e.g., $D_{h}^{(2)}\propto T^{-2d-2}=T^{-6}$ for graphene. At $T\to 0$, Eq. (9) predicts $D_{h}^{(2)}\to 2D_{0}^{(2)}\mathrm{sign}\,n$, where

[TABLE]

It may seem unusual that $D_{h}^{(2)}$ becomes doping-independent in this limit (except for the overall sign) but this can be rationalized by the dimensional analysis. Of course, at $T=0$ the system must be in the kinetic not hydrodynamic regime. Surprisingly, in the kinetic regime of graphene, $\sigma^{(2)}_{ilm}$ has a different tensorial structure:

[TABLE]

This result can be obtained from either the Boltzmann kinetic equation Manzoni et al. (2015) or the semiclassical limit $q\ll k_{F}$, $\omega\ll\varepsilon_{F}$ of the quantum random-phase approximation Cheng et al. (2017); Wang et al. (2016); Rostami et al. (2017). Here $\varepsilon_{F}$ and $k_{F}$ are the Fermi energy and momentum. Note that some of the related formulas in prior literature, e.g., Eq. (A.8) of Manzoni et al. (2015) and Eq. (42) of Rostami et al. (2017) are valid only for response to a longitudinal $\mathbf{E}$-field. If $\operatorname{curl}\mathbf{E}\neq 0$, the correct result is obtained only if the induced $\mathbf{B}$-field is included Cheng et al. (2017); Wang et al. (2016); SM . When extended further SM , such calculations show that $D^{(2)}_{k}=D_{0}^{(2)}$ at $T=0$ and $D^{(2)}_{k}\propto T^{-2}$ at high $T$. Hence, $D_{h}^{(2)}$ is twice larger than $D_{k}^{(2)}$ at $T=0$ but becomes smaller at high $T$, see Fig. 2(b).

A direct experimental probe of the second-order spectral weight is the second harmonic generation (SHG), which corresponds to $\omega_{2}=\omega_{1}=\omega$, $\omega_{3}=2\omega$, see Fig. S2(a). As explained above, the hydrodynamics predicts the SHG signal that is twice larger at low $T$ and much smaller at high $T$ compared to the standard kinetic theory Mikhailov (2011, 2016), see Fig. S2(b). The crossover from the kinetic regime to the hydrodynamic one would occur at temperature $T^{\ast}$ such that $\Gamma_{ee}\left(T^{\ast}\right)=\omega$. The measured SHG signal may look like as sketched by the dashed curve in Fig. S2(b).

Another effect controlled by $\sigma^{(2)}_{ilm}$ is the photon drag (PD), the generation of a dc current in response to a monochromatic beam of frequency $\omega$, see Fig. S2(a). (A recent work Tomadin and Polini (2013) studied a similar phenomenon for a surface plasmon playing the role of the incident beam.) To the second order in the in-plane field $\mathbf{E}(\mathbf{q},\omega)=(E_{x},E_{y})$ the PD is described by $\sigma^{(2)}_{ilm}$ evaluated at $\omega_{2}=-\omega_{1}=\omega$, and $\mathbf{q}_{1}=-\mathbf{q}_{2}=\mathbf{q}$. The PD in graphene has been previously studied in the kinetic regime Glazov and Ganichev (2014); Jiang et al. (2011); Karch et al. (2010). It was shown that the dc current can be parametrized by three constants $T_{1}$, $T_{2}$ and $\tilde{T}_{1}$, which multiply the three Stokes parameters of the incident beam. Coefficients $T_{1}$ and $T_{2}$ quantify the linear PD, $T_{2}$ and $\tilde{T}_{1}$ characterize the circular PD. Instead of the Stokes parameters, we find it convenient to use the incident angle $\theta$ and the $E_{y}$–$E_{x}$ phase delay $\psi$, so that $E_{x}=E\cos\alpha\cos\theta$, $E_{y}=E\sin\alpha e^{i\psi}$. Note that $\alpha=0$ means p-polarization and $\alpha=\pi/2$ means s-polarization. For a beam with the in-plane momentum $\mathbf{q}=(q_{x},0)$, the longitudinal and transverse current components are:

[TABLE]

where $C_{j}=\frac{1}{2}q_{x}|E|^{2}$, cf. Eq. (10) of Glazov and Ganichev (2014). To compute $T_{1}$, $T_{2}$ and $\tilde{T}_{1}$ for graphene in the hydrodynamic regime, we use the dissipative version of Eq. (2), which corresponds to retaining $\Gamma_{d}$ in the Euler equation (S23c). The resultant expression for $\sigma^{(2)}_{ilm}$ at arbitrary $\omega$ is ponderous SM . We present only the formulas for the drag coefficients:

[TABLE]

They are quite unlike those in the kinetic regime in which $\tilde{T}_{1}$ is nonzero, e.g.,

[TABLE]

This expression, which is a particular case of a general formula given in SM ; Glazov and Ganichev (2014), assumes that the scattering rate $\Gamma_{d}$ is due to short-range scatterers. The difference between the two regimes is illustrated in Fig. S2(c,d).

The following estimates suggest that the hydrodynamic regime $\Gamma_{d}<\omega<\Gamma_{ee}$ could be fairly wide in ultra clean graphene where electrons are scattered primarily by acoustic phonons, $\Gamma_{d}\approx\Gamma_{ep}$. The electron-phonon scattering rate $\Gamma_{ep}(T_{l},T,n)$ Principi et al. (2014) is a function of the lattice temperature $T_{l}$, electron temperature $T$, and doping $n$. From Ni et al. (2016, ) we estimate $\Gamma_{ep}(150\,\mathrm{K},150\,\mathrm{K},2\times 10^{12}\,\mathrm{cm}^{-2})\sim 0.3\,\mathrm{THz}$. On the other hand, $\Gamma_{ee}(T,n)$ is a function of $T$ and $n$. (In the kinetic regime $\omega\gg\Gamma_{ee}$, it may also depend on frequency.) Recent dc transport experiments Bandurin et al. (2016) indicate $\Gamma_{ee}(150\,\mathrm{K},10^{12}\,\mathrm{cm}^{-2})\sim 0.5\,\mathrm{THz}$, so the hydrodynamic region is narrow. There are two possible schemes to diminish $\Gamma_{ep}$ or enhance $\Gamma_{ee}$. The first one is to reduce $n$ to make electron gas non-degenerate, which should bring $\Gamma_{ee}$ to the theoretical maximum Schütt et al. (2011) of $4(e^{2}/\hbar\kappa v)^{2}T\sim 10\,\mathrm{THz}$. The other route is ultrafast pump-probe experiments Ni et al. (2016) that can keep the lattice cold, perhaps, at $T_{l}\sim 30\,\mathrm{K}$ but heat electrons to $T\sim 3000\,\mathrm{K}$.

The universal relation (3) between linear and nonlinear ac conductivities is the most important result of this Letter. Although we have used graphene as the example, this and our other formulas Eqs. (2), (9), etc., should apply as well to ultrapure metals and semiconductors Moll et al. (2016); de Jong and Molenkamp (1995), to surface states of topological insulators and Dirac/Weyl semimetals, provided they are in the hydrodynamic regime.

This work is supported by the DOE under Grant DE-SC0012592, by the ONR under Grant N00014-15-1-2671, by the NSF under Grant ECCS-1640173, and by the SRC. D. N. B. is an investigator in Quantum Materials funded by the Gordon and Betty Moore Foundation’s EPiQS Initiative through Grant No. GBMF4533. We thank G. Falkovich, M. Glazov, and G. Ni for discussions.

I Linear ac conductivity

I.1 Drude weight and demons in the hydrodynamic regime

As shown in literature Hartnoll et al. (2007); Müller et al. (2008, 2009); Kovtun (2012); Briskot et al. (2015); Narozhny et al. (2015); Sun et al. (2016); Lucas et al. (2016), the linear-response ac conductivity of a Dirac fluid at $q=0$ is given by

[TABLE]

which is Eq. (5) of the main text. The hydrodynamic Drude weight that enters Eq. (S1) is

[TABLE]

At zero temperature the hydrodynamic mass $m^{\ast}$ is no different from the Fermi-liquid effective mass $m^{\ast}=\hbar k_{F}/v_{F}$, where $v_{F}$ is the Fermi velocity. Hence, $D_{h}$ is equal to the conventional (kinetic) Drude weight $D_{k}$. For example, for parabolic band, $m^{\ast}$ is simply the band mass $m$. For graphene with weak ee interactions,

[TABLE]

where $g=4$ is the total spin-valley degeneracy Basov et al. (2014).

As usual, at finite $q$, the conductivity becomes a tensor

[TABLE]

Neglecting two subleading dissipative effects [$\sigma_{0}$ in Eq. (S1) and viscous damping], the longitudinal conductivity $\sigma_{L}$ is given by Kovtun (2012); Briskot et al. (2015)

[TABLE]

The longitudinal conductivity enters the equation for the dispersion of longitudinal collective modes. In 2D case, this equation reads Basov et al. (2014)

[TABLE]

The longitudinal mode in the hydrodynamic regime has been variously referred to as the sound Kovtun (2012), the energy wave Phan et al. ; Briskot et al. (2015), and finally, the demon Sun et al. (2016), which is our preference here. Equations (S5) and (S6) imply that away from charge neutrality, $n\neq 0$, the dispersion of the demon varies from $\omega\propto\sqrt{q}$ at low $q$ to $\omega\simeq v_{d}q$ at large $q$, see Fig. 1 of the main text. For neutral fluid, $D_{h}\to 0$, the demon dispersion is acoustic starting from $q=0$. The (asymptotic) speed of the demon is given by $v_{d}=v\sqrt{C_{\mathrm{ise}}}$ where $C_{\mathrm{ise}}$ is the dimensionless isentropic bulk modulus [Eq. (10) of the main text or Eq. (S24) below]. For a degenerate Fermi gas

[TABLE]

where $W\simeq n\varepsilon_{F}$ is the enthalpy density and $d$ is the space dimension; therefore,

[TABLE]

For Dirac dispersion $\varepsilon_{p}^{2}=(pv)^{2}+(mv^{2})^{2}$, the relation $\hbar k_{F}=\varepsilon_{F}v_{F}/v^{2}$ holds; thus,

[TABLE]

and $v_{d}=v_{F}/\sqrt{d}$, same as the speed of the first sound in a neutral Fermi liquid. For graphene,

[TABLE]

at any $T$ (see below). Therefore, $v_{d}=v/\sqrt{2}$ Kovtun (2012); Phan et al. ; Briskot et al. (2015); Sun et al. (2016).

I.2 Interpolation formula for the ac conductivity

In this Section we derive Eq. (8) of the main text, which smoothly connects the hydrodynamic and kinetic regimes of the linear-response theory. We start with the formula for the current

[TABLE]

in terms of the quasiparticle distribution function $f_{\mathbf{p}}$ and velocity $\mathbf{v}_{\mathbf{p}}=\partial\varepsilon_{p}/\partial\mathbf{p}$ as a function of momentum $\mathbf{p}$. For simplicity of notations, all the other quantum numbers such as spin, valley, and band index are omitted. We use $\sum\ldots$ to denote the summation over these quantum numbers combined with the integration $\int d^{2}p/(2\pi)^{2}\ldots$ over momentum.

Let us assume that the electric field in the system is position-independent and directed along $x$, i.e., $\mathbf{E}=\hat{\mathbf{x}}E(t)$. We want to compute the current $\mathbf{j}$ to the first order in $E(t)$. The result is different in the two regimes because the deviation $f_{\mathbf{p}}-f_{\mathbf{p}}^{(0)}$ of the distribution function from the equilibrium value $f_{\mathbf{p}}^{(0)}=\left[e^{\left(\varepsilon_{p}-\mu\right)/T}+1\right]^{-1}$ has different forms. It is proportional to $g_{p}=p_{x}\partial_{\varepsilon}f_{\mathbf{p}}^{(0)}$ in the hydrodynamic limit but to $g_{v}=v_{x}\partial_{\varepsilon}f_{\mathbf{p}}^{(0)}$ in the kinetic limit. To obtain the desired interpolation, we postulate that in general, $f_{\mathbf{p}}-f_{\mathbf{p}}^{(0)}$ is a certain linear combination

[TABLE]

To find the coefficients $a_{v}$ and $a_{p}$ we consider the Boltzmann kinetic equation

[TABLE]

We further assume that the linearized collision operator $\hat{I}$ acts within the space of functions given by Eq. (S12) and is characterized by two parameters: $\Gamma_{d}$, the scattering rate due to disorder and phonons, and $\Gamma_{ee}$, the electron-electron (ee) scattering rate. Mode $g_{v}$ is damped by both types of scattering but $g_{p}$ is immune to the ee one, which implies

[TABLE]

The condition that ee scattering conserves momentum fixes the coefficient $a_{pv}=D_{h}/(\pi e^{2}n)=1/m^{\ast}$ [Eq. (S2)], leading us to

[TABLE]

For $E(t)\propto e^{-i\omega t}$, the solution is

[TABLE]

Combining Eqs. (S11), (S12), and (S16), we get the linear conductivity

[TABLE]

which is Eq. (8) of the main text. Note that the obtained $\sigma(\omega)$ can be recast in the form of an extended Drude model Basov and Timusk (2005):

[TABLE]

The complex memory function $M(\omega)$ appearing in this equation is illustrated by Fig. S1. Both the effective scattering rate $1/\tau=\operatorname{Re}M(\omega)$ and the mass renormalization factor $\lambda=-\operatorname{Im}M(\omega)/\omega$ show step-like crossovers at the boundary $\omega\sim\Gamma_{ee}$ of the hydrodynamic and kinetic regimes.

II Second-order conductivity: general

The second-order nonlinear conductivity $\sigma_{ilm}^{(2)}$ determines the second-order current

[TABLE]

in response to the total electric field $\mathbf{E}$ in the system. By convention, $\sigma_{ilm}^{(2)}\left(\mathbf{q_{1}},\omega_{1};\mathbf{q_{2}},\omega_{2}\right)$ is chosen to be symmetrized, i.e., invariant under the interchange $(1\leftrightarrow 2,l\leftrightarrow m)$. Expanded to the linear order in momenta, the second-order conductivity must have the form

[TABLE]

where $\Sigma_{ilmn}$ is some isotropic rank-$4$ tensor. Any such tensor is a linear combination of the following three:

[TABLE]

In other words, $\sigma^{(2)}_{ilm}$ is fully characterized by three functions $G_{1}$, $G_{2}$, and $G_{3}$ such that

[TABLE]

Below we derive $\sigma^{(2)}_{ilm}$ and show it has a different form in the hydrodynamic and the kinetic regimes.

III Second-order conductivity in the hydrodynamic regime

To derive $\sigma_{ilm}^{(2)}$ in the hydrodynamic regime we solve the equations

[TABLE]

These equations are the same as Eqs. (7) of the main text, except we added phenomenological energy dissipation rate $\Gamma_{E}$ in Eq. (S23b) and chose the units $e=v=1$ to lighten the notations. Hence, the Lorentz factor in Eq. (S23c) is now $\gamma=1/\sqrt{1-u^{2}}$. The derivation of Eqs. (S23) can be found in literature Landau and Lifshitz (1984); Müller et al. (2008, 2009); Kovtun (2012); Briskot et al. (2015). The definitions of pressure $P$, energy density $n_{E}$, and enthalpy density $W$ deserve a comment. Whereas the current $\mathbf{j}$ is proportional to the actual charge density $n$, the pressure $P=P(n_{0},n_{E0})$ is the equilibrium thermodynamic parameter, which is a function of the proper density $n_{0}=n/\gamma$ and the proper energy density $n_{E0}$. The actual energy density is $n_{E}$ [Eq. (S23b)] and the enthalpy density is $W=n_{E0}+P$. Another key thermodynamic parameter is the dimensionless isentropic (ise) bulk modulus $C_{\mathrm{ise}}$. It is defined by Eq. (10) of the main text:

[TABLE]

The second equation in Eq. (S24) follows from the thermodynamic relation

[TABLE]

for the quantity $s_{n}=s/n_{0}$, with $s=s(n_{0},n_{E0})$ being the entropy density.

Suppose $\mathbf{E}(\mathbf{r},t)\propto e^{i\mathbf{q}\mathbf{r}-i\omega t}$ and define $\omega^{+}=\omega+i\Gamma_{d}$. To the first order in $\mathbf{E}$ we obtain, for $\Gamma_{E}=0$:

[TABLE]

Now let us assume that the electric field consists of two plane waves:

[TABLE]

To the second order in $\mathbf{E}$, various quantities of interest develop Fourier amplitudes of frequency and momenta $(\mathbf{q}_{3},\,\omega_{3})=(\mathbf{q}_{1}+\mathbf{q}_{2},\,\omega_{1}+\omega_{2})$. These amplitudes are given by

[TABLE]

From Eq. (S23c) and the Faraday law $\epsilon_{ikl}\partial_{k}E_{l}=-c^{-1}\partial_{t}B_{i}$ we obtain

[TABLE]

Keeping only the terms linear in $\mathbf{q}_{\nu}$, we find

[TABLE]

The current to the second order in field is

[TABLE]

Hence, the second-order conductivity is

[TABLE]

Here we added the neglected earlier energy dissipation rate $\Gamma_{E}$ in one of the terms in Eq. (S33). In principle, $\Gamma_{E}$ should appear in more than one place. However, we assume that $\Gamma_{E}$ is very small and its sole role is to resolve the indeterminacy of the ratio $q_{3}/\omega_{3}$ in the context of the photon drag problem where $q_{3},\omega_{3}\to 0$. The second-order spectral weight appearing in Eq. (S33) is

[TABLE]

where we restored physical units and replaced $n_{0}$ by $n$, $n_{E0}$ by $n_{E}$ to simplify notations. Let us discuss the value of $D^{(2)}_{h}$ in representative cases, assuming ee interaction corrections to pressure and enthalpy density are negligible. The result for particles with a parabolic dispersion can be obtained taking the nonrelativistic limit, in which $P\simeq({2}/{d})(n_{E}-nmv^{2})$ and $W\simeq nmv^{2}$. This gives

[TABLE]

In the massless case, one finds $P=n_{E}/d$ and $W=(1+1/d)n_{E}$, so that

[TABLE]

Taking $d=2$ for graphene, we get

[TABLE]

Functions $G_{1}$, $G_{2}$, and $G_{3}$ [Eq. (S22)] corresponding to Eq. (S33) are

[TABLE]

In the collisionless limit $\Gamma_{d}\to 0$, these formulas simplify to

[TABLE]

while Eq. (S33) reduces to

[TABLE]

which is equivalent to Eq. (2) of the main text.

IV Second-order conductivity in the kinetic regime

The kinetic regime corresponds to the frequency range $\Gamma_{ee}\ll\omega\ll\varepsilon_{F}$. The linear and nonlinear conductivities in this regime can be computed by solving the Boltzmann kinetic equation

[TABLE]

In this section, we again set $e=1$ and suppress the subscripts ${\mathbf{p}}$ in $\mathbf{v}_{\mathbf{p}}$, $f_{\mathbf{p}}$. One should not confuse the quasiparticle velocity $\mathbf{v}$ at finite $\mathbf{p}$, a vector, with $v$, the limiting velocity at $p=\infty$, a scalar.

The magnetic-field term in Eq. (S42) can be expressed with the help of the kernel

[TABLE]

similar to Eq. (S29), leading to

[TABLE]

The second-order current we want to compute is

[TABLE]

To do this we need to specify the collision integral $\hat{I}\left[f\right]$.

IV.1 Nonconserving relaxation-time approximation

It is useful to consider first the approximation $\hat{I}\left[f\right]=-\Gamma\left(f-f^{(0)}\right)$, with $\Gamma$ being an energy-independent relaxation rate. This is probably the simplest model one can study. However, one should keep in mind that this approximation is flawed because it may not conserve the particle number.

Assuming the electric field is composed of two plane waves [Eq. (S26)], we expand $f$ to the first and second order in field:

[TABLE]

Hence, the second-order conductivity is

[TABLE]

The evaluation of this expression for Dirac electrons in graphene is tedious but straightforward. The final result is

[TABLE]

Its collisionless limit $\Gamma\to 0$ is

[TABLE]

which implies [cf. Eq. (S22)]

[TABLE]

The second-order spectral weight in these equations is

[TABLE]

where [Eq. (S38)]

[TABLE]

In the limit of $T\to 0$, we have $D_{k}^{(2)}=D_{0}^{(2)}\,\mathrm{sign}(\mu)$. At $T\gg\mu$, the asymptotic behavior of the chemical potential is $\mu\propto n/T$, see, e.g., Supplemental material of Ref. Sun et al., 2016. Therefore, $D_{k}^{(2)}\propto n/T^{2}$ as mentioned in the main text. However, at high $T$, an interband contribution to $\sigma^{(2)}_{ilm}$, not included in our semiclassical approach, may become important.

IV.2 Multiple relaxation-time approximation

Let us assume now that the collision operator $\hat{I}\left[f\right]$ is linear and diagonal in the angular momentum basis, so that the Boltzmann equation can be written as

[TABLE]

where $\hat{L}_{0}$ is the operator

[TABLE]

and $\Gamma_{l}$ is the scattering rate for the angular momentum $l$. This rate may depend on the quasiparticle energy $\varepsilon$. The model conserves the number of particles if $\Gamma_{0}=0$. The action of $\hat{L}_{0}$ can be written in terms of the complex frequencies

[TABLE]

Instead of Eq. (S47) we now get a more complicated expression:

[TABLE]

To calculate $\sigma_{ilm}^{(2)}(\mathbf{q}_{1},\omega_{1};\mathbf{q}_{2},\omega_{2})$ we need the $(\mathbf{q}_{3},\omega_{3})$ Fourier harmonic of $f^{(2)}$:

[TABLE]

[The argument $(\mathbf{q}_{3},\omega_{3})$ of $f^{(2)}$ is omitted.] For our purpose of computing the terms linear in gradients the expansion

[TABLE]

suffices. It yields

[TABLE]

To do the summation over the angular directions, we expand all the variables in the angular momentum basis. To this end, we do a set of unitary transformations. For example, the velocity goes from $\mathbf{v}=\left(v_{x},v_{y}\right)^{T}=\left|\mathbf{v}\right|\left(\cos\phi,\,\sin\phi\right)^{T}$ to $\left(\tilde{v}_{-},\tilde{v}_{+}\right)^{T}$:

[TABLE]

We do the same transformation for the momentum-space derivatives:

[TABLE]

To the electric fields and spatial momenta we apply a conjugate transformation, $\tilde{E}_{m}=E_{j}U^{\dagger}_{jm}$, $\tilde{q}_{m}=q_{j}U^{\dagger}_{jm}$, in order to leave the scalar products $E_{i}v_{i}$, $E_{i}\partial_{p_{i}}$, $q_{i}v_{i}$ invariant. The net effect on Eq. (S60) is simply to add tildes for every variable. The second-order current becomes

[TABLE]

Only the terms of zero net angular momentum, i.e., $i+n+l+m=0$ survive after the summation. Since $i$, $n$, $l$, $m$ have values $\pm 1$, they have to appear in opposite-sign pairs. This constraint can be implemented with the help of the transformed rank-2 and rank-4 isotropic tensors

[TABLE]

The subsequent calculations are done for $T=0$ where $\partial_{\varepsilon}f^{(0)}=-\delta\left(\varepsilon-\varepsilon_{F}\right)$. We obtain

[TABLE]

To convert back to the $(x,\,y)$ coordinates, one simply needs to drop the tildes everywhere. Therefore, the second-order nonlinear optical conductivity is

[TABLE]

Completing the symmetrization step $\left(\begin{smallmatrix}1&\leftrightarrow&2\\ l&\leftrightarrow&m\end{smallmatrix}\right)$, we get the following:

[TABLE]

The corresponding functions $G_{1}$, $G_{2}$, and $G_{3}$ [Eq. (S22)] are

[TABLE]

Setting the particle number relaxation rate $\Gamma_{0}$ to zero, which is the physical case, we get

[TABLE]

In the collisionless limit, $\Gamma_{l}\to 0$, these formulas reduce to Eq. (S51).

V Third-order conductivity in the hydrodynamic regime

The third-order ac conductivity $\sigma_{ilmn}^{(3)}\left(\mathbf{q}_{1},\omega_{1};\mathbf{q}_{2},\omega_{2};\mathbf{q}_{3},\omega_{3}\right)$ is defined as

[TABLE]

Unlike the second-order conductivity, $\sigma^{(3)}_{ilmn}$ can approach a nonzero value at $q=0$ in inversion-symmetric systems. We will compute this value and disregard $\mathcal{O}(q^{2})$ nonlocal corrections. The calculation is simplified by the observation that Eqs. (S25) and (S27) yield $n^{(1)}=n^{(2)}=n_{E}^{(1)}=W^{(1)}=P^{(1)}=\mathcal{O}(q)\to 0$ in this approximation. An alternative way to get the same result is to neglect spatial gradients in Eqs. (S23c), (S23b) and (S23a), after which the hydrodynamic equations reduce to

[TABLE]

The last equation entails $n=n^{(0)}$, and so $j_{i}^{(3)}=n^{(0)}u^{(3)}_{i}$. The third-order velocity can be found from

[TABLE]

Since $\gamma^{2}W=n_{E}+P$, we have

[TABLE]

where “$\mathrm{perm}.$” stands for permutations among subscripts $1$, $2$, and $3$, corresponding to frequencies $\omega_{1}$, $\omega_{2}$, and $\omega_{3}$, respectively. The equation for the Fourier amplitude $u^{(3)}(\omega)$ of the combined frequency $\omega=\omega_{1}+\omega_{2}+\omega_{3}$ becomes

[TABLE]

Therefore,

[TABLE]

(For brevity, we omitted “$\mbox{}+\mathrm{perm}.$” in the above equations.) In the dissipationless limit Eq. (S79) simplifies to [cf. Eq. (S24)]

[TABLE]

Therefore,

[TABLE]

where $e$ and $v$ were restored.

We can compare our formula for the third-order ac conductivity in the hydrodynamic regime with other results in the literature for the case $\omega_{1}=\omega_{2}=\omega_{3}$, which corresponds to the third harmonic generation. This effect is controlled by the conductivity $\sigma^{(3)}_{ilmn}(0,\omega_{1};0,\omega_{1};0,\omega_{1})$. Applied to graphene at $T=0$, our result for $D_{h}^{(3)}$ is twice larger than the third-order spectral weight from the collisionless Boltzmann transport theory Mikhailov (2016). Compared to the linear response, the third-order current is suppressed by the small parameter $\xi=\left(\frac{-eE/\omega}{m^{\ast}v}\right)^{2}$. At zero temperature, neglecting exchange-correlation corrections, $m^{\ast}v$ is just the Fermi momentum $p_{F}$, so that $\xi=({\delta p}/{p_{F}})^{2}$. The quantity $\delta p=-eE/\omega$ is equal by the order of magnitude to the change in electron momentum caused by the electric field during one half cycle of the sum-frequency oscillations, $\delta t\sim\pi/\omega$. The ratio of $\xi$ factors for a nonrelativistic and ultrarelativistic Dirac fluids is $\sim(v_{F}/v)^{2}\ll 1$. This factor vanishes for a system with a parabolic dispersion corresponding to $v\to\infty$. Indeed, for such a system all nonlinearities at zero $q$ should be absent because of the Galilean invariance. On the other hand, the linear and second-order conductivities, $\sigma$ and $\sigma^{(2)}_{ilm}$, do not show this contrasting behavior because they do not contain $v$ explicitly.

VI Applications and summary

VI.1 Photon drag

The photon drag effect is the generation of dc current by a light incident on the sample. Unlike optical rectification and photogalvanic effect, the photon drag current is the result of the transfer of the linear momentum of photons $\mathbf{q}$ to free carriers Glazov and Ganichev (2014). This is why photon drag can appear only if $\mathbf{q}=\hat{\mathbf{x}}q_{x}+\hat{\mathbf{z}}q_{z}$ is not strictly normal to the $x$–$y$ plane of the sample, see Fig. S2(a). An alternative classical picture of the photon drag is the carrier drift in the crossed electric and magnetic fields of the electromagnetic wave, and so the photon drag is also sometimes referred to as the dynamical Hall effect. The drag current can have both longitudinal $j_{x}$ and transverse $j_{y}$ components. Let the in-plane component of the electric field be $\mathbf{E}=\mathbf{E}_{0}e^{i(\mathbf{q}\mathbf{r}-\omega t)}+\mathrm{c.c.}$ where $\mathbf{E}_{0}=\hat{\mathbf{x}}E_{x}+\hat{\mathbf{y}}E_{y}$. The polarization of the incident wave in the $x$–$y$ plane is important. This polarization can be specified in terms of the Stokes parameters $s_{0}=|E_{x}|^{2}+|E_{y}|^{2}$, $s_{1}=|E_{x}|^{2}-|E_{y}|^{2}$, $s_{2}=E_{x}E_{y}^{\ast}+E_{y}E_{x}^{\ast}$, and $s_{3}=-i(E_{x}E_{y}^{\ast}-E_{y}E_{x}^{\ast})$. From Eq. (S33), we can calculate the induced dc current components as Jiang et al. (2011)

[TABLE]

The coefficients $T_{1}$ and $T_{2}$ are as follows:

[TABLE]

where

[TABLE]

and $G_{a}$ are the functions introduced in Eq. (S22). For $\tilde{T}_{1}$, we get

[TABLE]

For the hydrodynamic regime, we take $G_{a}(\omega,-\omega)$ from Eq. (S39) and obtain

[TABLE]

For the case of graphene, $C_{\mathrm{ise}}=1/2$, these equations give

[TABLE]

which is Eq. (15) of the main text.

In the kinetic regime, the photon drag coefficients are more complicated. Equations (S71)–(S73) for $G_{a}$ can be used to compute them for graphene at zero temperature. We get the following:

[TABLE]

in agreement with Refs. Glazov and Ganichev, 2014; Karch et al., 2010; Jiang et al., 2011. If the dominant electron scattering in graphene is due to short-range impurities, then the scattering rates $\Gamma_{1}$ and $\Gamma_{2}$ for the $p$- and $d$-wave angular deformations of the Fermi surface obey the relations

[TABLE]

When substituted into the general formulas above, followed by the notation change $\Gamma_{1}\to\Gamma_{d}$, these relations lead to Eq. (16) of the main text.

Instead of the Stokes parameters, we can use two angles $\psi$ and $\alpha$ such that $E_{x}=E\cos\alpha\cos\theta$, $E_{y}=E\cos\alpha e^{i\psi}$. Note that $\alpha=0$ means p-polarization and $\alpha=\pi/2$ means s-polarization, see Fig. S2(a). The formulas for $j_{x}$ and $j_{y}$ become

[TABLE]

In the hydrodynamic regime where $\tilde{T}_{1}=0$, the transverse current $j_{y}$ has no component proportional to $\sin\psi$. However, $j_{y}$ does have such a component in the kinetic regime, as illustrated by Fig. S2(d). This distinction may be used to identify the two regimes in experiments.

VI.2 Second harmonic generation

The second-harmonic generation (SHG) signal is proportional to the $\left(2\mathbf{q},2\omega\right)$ Fourier harmonic of the second-order ac current, which is given by Glazov (2011)

[TABLE]

where

[TABLE]

Neglecting the damping, for graphene in the kinetic regime, we get

[TABLE]

In the hydrodynamic regime, we find

[TABLE]

This implies that SHG signal has the same polarization dependence in the hydrodynamic and kinetic regimes but the magnitude of the response is different because it is controlled by either $D^{(2)}_{h}$ or $D^{(2)}_{k}$. At zero temperature the ratio $D^{(2)}_{h}/D^{(2)}_{k}$ is equal to $2$, but at high temperature it rapidly decreases, see Fig. S2(b). Experimentally, this difference may be observed as the electron temperature is increased and the system crosses over from the kinetic to the hydrodynamic regime at some $T^{\ast}=T^{\ast}(\omega)$. This crossover temperature is the solution of the equation $\Gamma_{ee}(T^{\ast})=\omega$. As $T^{\ast}$ is approached from below, the SHG signal may increase, by up to a factor of two from its $T=0$ value, as sketched by the dashed line in Fig. S2(b). When the temperature is raised beyond $T^{*}$, the system enters the hydrodynamic regime where the SHG signal should drop due to decreasing $D^{(2)}_{h}$. The transient states of high electron temperatures can be realized with intense photoexcitation.

VI.3 Summary

tables for the second-order conductivity

As we pointed out earlier, $\sigma^{(2)}_{ilm}$ is fully characterized by three functions $G_{1}$, $G_{2}$, and $G_{3}$. Shown in Table 1 are $G_{1}$, $G_{2}$, $G_{3}$ and $T_{1}$, $T_{2}$, $\tilde{T}_{1}$ in different regimes and for different band dispersions. In addition, the same formulas in the clean limit are summarized in Table 2.

Electron systems with parabolic dispersion is an interesting case. In such systems $\sigma^{(2)}_{ilm}$ has the same form in the kinetic and hydrodynamic regimes but only in the absence of momentum dissipation. As mentioned in the main text, in this system the random-phase approximation (RPA) also gives the same $\sigma^{(2)}_{ilm}$ to the linear order in $q$ in the absence of dissipation. In the diagrammatic derivation Stolz (1967) of this RPA result only the “diamagnetic” terms contribute to $\sigma^{(2)}_{ilm}$. Those diamagnetic terms are all determined by the linear-response Drude weight. The “paramagnetic” term, that is, a single-loop diagram with three current vertices vanishes to the first order in $q$. This is superficially similar yet apparently unrelated to Furry’s theorem in quantum electrodynamics, which says that fermion loops with odd number of photon vertices vanish because of electron-positron symmetry. For Dirac electrons in graphene one may invoke Furry’s theorem to explain vanishing of the spectral weight at $\mu=0$ [see Eq. (S52)]. However, in the case of interest, $\mu\neq 0$, the three-point current correlation function Cheng et al. (2017); Wang et al. (2016); Rostami et al. (2017) is finite. In fact, it is the diamagnetic contribution that vanishes, so that $\sigma^{(2)}_{ilm}$ is determined solely by this paramagnetic term.