L\'evy flights and hydrodynamic superdiffusion on the Dirac cone of Graphene

TL;DR

This paper demonstrates that electron collisions in graphene at neutrality point lead to superdiffusive behavior of collective excitations, modeled as Le9vy flights in momentum space, with broader implications for quantum-critical systems.

Contribution

It introduces a fractional Fokker-Planck equation to describe superdiffusion in graphene, linking electron collisions to Le9vy flights and highlighting their significance in quantum-critical dynamics.

Findings

Electron-electron collisions cause superdiffusive collective excitations.

Superdiffusion modeled as Le9vy flights in momentum space.

Implications for out-of-equilibrium quantum-critical systems.

Abstract

We show that hydrodynamic collision processes of graphene at the neutrality point can be described in terms of a Fokker-Planck equation with fractional derivative, corresponding to a L\'evy flight in momentum space. Thus, electron-electron collisions give rise to frequent small-angle scattering processes that are interrupted by rare large-angle events. The latter give rise to superdiffusive dynamics of collective excitations. We argue that such superdiffusive dynamics is of more general importance to the out-of-equilibrium dynamics of quantum-critical systems.

Figures

Equations

$$\frac{\partial \rho }{\partial t}+D_{\mu }{\left|△\right|}^{\frac{\mu }{2}}\rho =0,$$

$$\left({\partial }_t+{\mathbf{v}}_{\mathbf{k}\lambda }\cdot {\nabla }_{\mathbf{x}}-{\tau }_L^{-1}{\left(\frac{{\partial }^2}{\partial {\theta }^2}\right)}^{1/2}\right)f_{\mathbf{k}\lambda }=S_{\mathbf{k}\lambda },$$

$$\hslash {\tau }_L^{-1}\approx 11.66{\alpha }^2{k}_{B}T,$$

$$\left({\partial }_t+{\mathbf{v}}_{\mathbf{k}\lambda }\cdot {\nabla }_{\mathbf{x}}+\mathbf{F}\left(\mathbf{x},t\right)\cdot {\nabla }_{\mathbf{k}}+\mathcal{C}\right)f_{\mathbf{k}\lambda }\left(\mathbf{x},t\right)=0$$

$$\mathcal{C}{f}_1$$

$$f_{\mathbf{k}\lambda }\left(\mathbf{x},t\right)=f_{k\lambda }^{\left(0\right)}+f_k^{\left(0\right)}\left(1-f_k^{\left(0\right)}\right){\psi }_{\mathbf{k}\lambda }\left(\mathbf{x},t\right).$$

$$\mathcal{L}=\left({\partial }_t+{\mathbf{v}}_{\mathbf{k}\lambda }\cdot {\nabla }_{\mathbf{x}}\right)f_k^{\left(0\right)}\left(1-f_k^{\left(0\right)}\right)$$

$${\left(\mathcal{C}\psi \right)}_{\mathbf{k}\lambda }$$

$${\chi }_{\mathbf{k}\lambda }^{\left(m,s\right)}={\lambda }^m{e}^{im\theta }\left\{1,\lambda ,\lambda v_{0}k/\left(k_{B}T\right)\right\},$$

$$⟨m\left|\mathcal{C}\right|m^{\prime }⟩=\frac{\ln 2}{\pi }{\delta }_{m,m^{\prime }}{\tau }_m^{-1}.$$

$$\hslash {\tau }_1^{-1}=3.646{\alpha }^2{k}_{B}T$$

$${\tau }_m^{-1}={\tau }_1^{-1}\left(\kappa \left|m\right|-{\kappa }^{\prime }\right),$$

$$⟨\theta \left|\mathcal{C}\right|{\theta }^{\prime }⟩$$

$${\left(\frac{{\partial }^{2}g\left(\theta \right)}{\partial {\theta }^2}\right)}^{1/2}=\frac{\partial }{\partial \theta }{\int }_0^{2\pi }\frac{g\left({\theta }^{\prime }\right)}{\tan \left(\frac{\theta -{\theta }^{\prime }}{2}\right)}d{\theta }^{\prime },$$

$$S_{\mathbf{k}\lambda }\left(t\right)=\delta \left(t\right)f_k^{\left(0\right)}\left(1-f_k^{\left(0\right)}\right)\sum _m\delta h_{\lambda m}{e}^{im\theta }.$$

$${\psi }_{\lambda }\left(\theta ,t\right)=\delta h\Theta \left(t\right)\frac{\sinh \left(t/{\tau }_L\right)}{\cosh \left(t/{\tau }_L\right)-\lambda \cos \left(\theta \right)}.$$

$${\psi }_{+}\left(t,\theta =0\right)\approx \delta h\frac{{\tau }_L}{\pi t},$$

$${\psi }_{+}\left(t,\theta \ne 0\right)\approx \frac{\delta h}{4\pi {\sin }^2\left(\theta /2\right)}\frac{t}{{\tau }_L}.$$

$$\frac{\partial s\left(t\right)}{\partial t}=4k_B\sum _{\lambda }{\int }_k\log \left(\frac{1-f_{\mathbf{k}\lambda }}{f_{\mathbf{k}\lambda }}\right)\frac{\partial f_{\mathbf{k}\lambda }}{\partial t}.$$

$$\frac{\partial s\left(t\right)}{\partial t}=\frac{4\log 2}{9\zeta \left(3\right)}\frac{s_{\mathrm{eq}}}{{\tau }_L}\frac{{\left(\delta h\right)}^2}{{\sinh }^2\left(t/{\tau }_L\right)},$$

$$\delta q\left(t\right)\propto T{s}_{\mathrm{eq}}\delta h^2\frac{{\tau }_L}{t}.$$

$$\tilde{\psi }\left(\theta ,\Delta t\right)=\frac{\sinh \left(\Delta t/{\tau }_L\right)}{\cosh \left(\Delta t/{\tau }_L\right)-\cos \left(\theta \right)}.$$

$${\theta }_N=\sum _{i=0}^N\Delta {\theta }_i.$$

$${\tilde{\psi }}_{\text{normal}}\left(\theta ,\Delta t\right)=\frac{e^{\cos \left(\theta \right)/\Delta t}}{2\pi I_0\left(1/\Delta t\right)}.$$

$$H_0=v\hslash {\int }_k\sum _{\alpha \beta }{\psi }_{\alpha }^{†}\left(\mathbf{k}\right){\left(\mathbf{k}\cdot \sigma \right)}_{\alpha \beta }{\psi }_{\beta }\left(\mathbf{k}\right)$$

$$U_{\mathbf{k}}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{cc}\frac{k_{x}+ik_{y}}{k}&1\\ -\frac{k_{x}+ik_{y}}{k}&1\end{array}\right)$$

$$U_{\mathbf{k}}v\hbar\mathbf{k\cdot\sigma}U_{\mathbf{k}}^{-1}=\left(\begin{array}[]{cc}v\hbar k&0\\ 0&-v\hbar k\end{array}\right).$$

$$H_0=v\hslash {\int }_{\mathbf{k}}\sum _{\lambda =\pm }\lambda k{\gamma }_{\mathbf{k},\lambda }^{†}{\gamma }_{\mathbf{k},\lambda }.$$

$$H_{\mathrm{int}}=\frac{1}{2}{\int }_{k,k^{\prime },q}\sum _{\alpha \beta }V\left(\mathbf{q}\right){\psi }_{\alpha }^{†}\left(\mathbf{k}+\mathbf{q},t\right){\psi }_{\beta }^{†}\left({\mathbf{k}}^{\prime }-\mathbf{q},t\right){\psi }_{\beta }\left({\mathbf{k}}^{\prime },t\right){\psi }_{\alpha }\left(\mathbf{k},t\right)$$

$$H_{\mathrm{int}}=\frac{1}{2}{\int }_{k,k^{\prime },q}\sum _{\alpha \beta }{T}_{\lambda \mu {\mu }^{\prime }{\lambda }^{\prime }}\left({\mathbf{k},\mathbf{k}}^{\prime },\mathbf{q}\right){\gamma }_{{\lambda }^{\prime }}^{†}\left(\mathbf{k}+\mathbf{q},t\right){\gamma }_{\mu }^{†}\left({\mathbf{k}}^{\prime }-\mathbf{q},t\right){\gamma }_{{\mu }^{\prime }}\left({\mathbf{k}}^{\prime },t\right){\gamma }_{\lambda }\left(\mathbf{k},t\right)$$

Full text

The kinetics of large gravitational systems such as globular clusters in galaxies or of a classical charged plasma are governed by continuous collisions with small-angle scatterings. The origin for this behavior is the long-range character of the Newton or Coulomb force, respectively. Such small-angle collisions behave in velocity space like drag and diffusion events, where a Fokker-Planck equation offers an efficient descriptionChandrasekhar1943 ; Rosenbluth1957 ; Chernoff1990 . Collisions can thus be seen as a Gaussian random walk in phase space. The velocity of a plasma or gravitational dust particle undergoes ordinary Brownian motion.

Quantum many-body systems that are near a quantum-critical point are governed by soft modes that will also induce effective long-range interactionsSachdevbook . This begs the question whether such quantum-critical systems also allow for an effective Fokker-Planck description of the non-equilibrium kinetics; in the collision-dominated hydrodynamic regime and in the crossover regime from hydrodynamic to ballistic dynamics. Candidate systems are itinerant electrons near magnetic or nematic quantum phase transitionsHertz1976 ; Moriya1984 ; Millis1993 ; Laughlin2001 ; Millis2002 ; Abanov2003 ; Loehneysen2007 ; DellAnna2007 ; Metlitski2010 ; Schattner2016 , the superconductor-insulator phase transitionDamle1997 , or graphene near the Dirac pointSheehy2007 . Anomalous diffusion was even shown to be present in two-dimensional Fermi liquidsGurzhi1995 ; Gurzhi1996 ; Maslov2011 ; Pal2012 ; Maslov2017 ; Ledwith2019 .

In this paper we analyze the quantum kinetics of graphene near the Dirac point with electron-electron Coulomb interaction. We show that the kinetic theory at charge neutralityFritz2008 ; Kashuba2008 ; Mueller2009 ; Schuett2011 ; Kiselev2019 can be expressed in terms of a Fokker-Planck equation, yet with fractional derivative with respect to the momentum direction. The underlying random processes are Lévy flightsLevy1954 ; Mandelbrot1977 ; Feller1972 , non-Gaussian random walks whose step widths are distributed according to a powerlaw. The slowly decaying tail of the step-width distribution makes it impossible to define a diffusion constant or to use a conventional Fokker-Planck equation. However, a diffusion equation of the form

[TABLE]

with appropriately generalized fractional derivativeSamko1993 ; Hermann2014 can be used to describe such random walks. Lévy flights have been discussed to model the migration pattern of animals as they search for resourcesBartumeus2005 ; Reynolds2009 , the high-frequency index dynamics of the stock marketMantegna1997 , or to describe durations between consecutive earthquakesCorral2006 . In our system they correspond to random walks in momentum space with powerlaw weight for large momentum-transfer processes. We demonstrate that the collision operator due to electron-electron interactions in graphene takes the form of a fractional derivative. Then the Boltzmann equation becomes a fractional Fokker-Planck equation, similar to Eq.1 with exponent $\mu=1$:

[TABLE]

where $\theta$ determines the electron momentum direction: $\mathbf{k}=k\left(\cos\theta,\sin\theta\right)$. The precise definition of the fractional derivative is given below. This result implies that the out-of-equilibrium dynamics of graphene in the hydrodynamic regime is governed by a wrapped Cauchy flightLevy1939 ; Mardia1999 , a specific Lévy flight on the Dirac cone. In Fig. 1a we show a simulation of ordinary Brownian motion on a ring and of the wrapped Cauchy flight. Details of this simulation are summarized insupplementary material . The occurrence of rare large-angle jumps is clearly visible. The corresponding phase-space dynamics is sketched in Fig. 1b. While the direction of $\mathbf{k}$ undergoes anomalous diffusion, its magnitude $k\equiv\left|\mathbf{k}\right|$ is of the order of $k_{B}T/v_{0}$ with the graphene group velocity $v_{0}\approx 10^{8}{\rm cm}/{\rm s}$. The characteristic time of the process is $\tau_{L}$ with

[TABLE]

where the fine-structure constant of graphene is $\alpha=e^{2}/\left(\hbar\epsilon v_{0}\right)$. $\tau_{L}$ agrees up to a numerical coefficient with the collision time for the hydrodynamic transport behavior of graphene at the Dirac pointFritz2008 ; Kashuba2008 ; Mueller2009 . Below we discuss how $\tau_{L}$ is determined. Such a time scale was recently observed experimentally in THz spectroscopy of graphene at charge neutralityGallagher2019 .

Lévy flights in graphene have been discussed in Ref.Gattenloehner2016 , where an egineered distribution of adatoms was shown to result in a superdiffusive behavior of charge carriers, and in Ref.Briskot2014 in the context of highly photo-excited carriers that relax according to a cascade of processes

• a behavior with interesting implications for pump-probe experiments. This can be seen as a superdiffusion in energy space far from equilibrium. It affects the magnitude of the momentum. Here we focus on the low-energy hydrodynamic regime and find a very different behavior for the directional diffusion in momentum space. Nevertheless, these results strongly suggest that superdiffusive phase-space dynamics is a more common phenomenon in quantum-critical systems.

We start from the Boltzmann equation

[TABLE]

for the electron distribution function $f_{\mathbf{k}\lambda}\left(\mathbf{x},t\right)$ where $\mathbf{k}$ refers to the momentum and $\lambda=\pm 1$ labels the upper and lower cone of the Dirac spectrum $\varepsilon_{\mathbf{k}\lambda}=\lambda v_{0}\left|\mathbf{k}\right|$. $\mathbf{v}_{\mathbf{k}\lambda}=\partial\varepsilon_{\mathbf{k}\lambda}/\partial\mathbf{k}$ is the velocity vector and $\mathbf{F}\left(\mathbf{x},t\right)$ some external force, e.g. due to an external electric field. ${\cal C}$ is the Boltzmann collision operator due to electron-electron interactions and was derived to order $\alpha^{2}$ in Ref.Fritz2008 from a Keldysh-Schwinger approach; see also insupplementary material . It takes the usual form of a two-body interaction:

[TABLE]

The transition probability $W_{12,34}$ is due to the electron-electron Coulomb interaction $e^{2}/\epsilon$ of Dirac fermions that are confined to a two-dimensional system. $\epsilon$ is the dielectric constant determined by the substrate. For free standing graphene $\epsilon=1$ and the fine-structure constant $\alpha\approx 2.2$ is of order unity. A renormalization group analysis shows that $\alpha$ flows towards weak coupling, justifying our perturbative approachSheehy2007 .

As usual, the kinetic distribution function $f_{\lambda,\bm{k}}$ is expanded around the local equilibrium distribution $f_{k\lambda}^{0}=\left(e^{\beta\left(\epsilon_{\lambda,\bm{k}}-\mu\right)}+1\right)^{-1}$ and parametrized as ($f_{k}^{\left(0\right)}=f_{k+}^{\left(0\right)}$):

[TABLE]

We linearize the Boltzmann equation with respect to $\psi_{\mathbf{k}\lambda}\left(\mathbf{x},t\right)$. With the Liouville operator

[TABLE]

we obtain a compact formulation of the Boltzmann equation: $\left({\cal L}+{\cal C}\right)\psi=S$. $S_{\mathbf{k}\lambda}\left(\mathbf{x},t\right)$ contains external perturbations, such as those due to a space and time dependent electric field or flow-velocity gradient. The operators ${\cal L}$ and ${\cal C}$ act on the momentum and band indices $\mathbf{k}$ and $\lambda$, respectively. Taking into account the kinematic constraints of the linear Dirac spectrum, the collision operator becomes:

[TABLE]

where the matrix elements $\gamma_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{q}}^{\left(1,2\right)}$ are given in Ref.supplementary material and $\int_{k}\cdots=\int\frac{d^{2}k}{\left(2\pi\right)^{2}}\cdots$. One easily finds the zero modes that correspond to the conservation lawsFritz2008 . Eq.4 was obtained by projecting the distribution function onto the helical eigenstates of the problem. The same projection was performed in the derivation of the collision operatorFritz2008 ; supplementary material .

The usual analysis of the Boltzmann equation proceeds as follows: One performs a Fourier transformation from $\left(\mathbf{x},t\right)$ to $\left(\mathbf{q},\omega\right)$ and introduces a complete set of states $\chi_{\mathbf{k}\lambda}^{\left(s\right)}$ to evaluate the matrix elements $\left\langle s\left|{\cal L}+{\cal C}\right|s^{\prime}\right\rangle$ with scalar product $\left\langle s|s^{\prime}\right\rangle=\sum_{\mathbf{k}\lambda}\chi_{\mathbf{k}\lambda}^{\left(s\right)*}\chi_{\mathbf{k}\lambda}^{\left(s^{\prime}\right)}$. The Liouville operator becomes ${\cal L}=\left(-i\omega+i\mathbf{v}_{\mathbf{k}\lambda}\cdot\mathbf{q}\right)f_{k}^{\left(0\right)}\left(1-f_{k}^{\left(0\right)}\right)$. The distribution function then follows as $\psi=\left({\cal L}+{\cal C}\right)^{-1}S$. For finite $\omega$ or $\mathbf{q}$ the operator ${\cal L}+{\cal{\cal C}}$ is nonsingular. This program is somewhat simplified for graphene at charge neutrality. As shown in Refs.Fritz2008 ; Kashuba2008 ; Mueller2009 ; Kiselev2019 , scattering processes where all momenta are collinear are enhanced by a factor $\log\left(1/\alpha\right)$. This can be used to identify the dominant modes, derived in the supplementary material:

[TABLE]

where $m\in\mathbb{Z}$ is the angular momentum quantum number while $s=1\cdots$3 labels the collinear modes for given $m$. We solve the kinetic equation by projecting it onto the dominant collinear modes $\chi_{\mathbf{k}\lambda}^{\left(m,s\right)}$, but checked that our key conclusions are unchanged if we chose a larger set of basis functions. Also, if we restrict our considerations to the transport of charge due to external electric fields, it suffices to consider the modes $\chi_{\mathbf{k}\lambda}^{\left(m,1\right)}=\lambda^{m}e^{im\theta}$ of Eq. (9). For simplicity we confine ourselves to electric-field source terms and only discuss this mode. The generalization to other modes is straightforward.

The low-energy Dirac Hamiltonian is rotationally invariant such that the collision operator becomes diagonal in the angular momentum representation

[TABLE]

The diagonal elements are, besides a convenient prefactor, the scattering rates of the corresponding angular momentum channel. $\tau_{0}^{-1}=0$ due to charge conservation, while the collision rate

[TABLE]

for $m=1$ was determined in Ref.Fritz2008 to yield the optical conductivity $\sigma\left(\omega\right)=\frac{e^{2}}{h}4\ln 2k_{B}T\left(-i\hbar\omega+\hbar\tau_{1}^{-1}\right)^{-1}.$ $\tau_{1}^{-1}$ was recently observed in Ref.Gallagher2019 using a waveguide setup; a demonstration of quantum-critical hydrodynamic transport. The dramatic violation of the Wiedemann-Franz law at charge neutrality is another important indication for electronic hydrodynamics at charge neutralityCrossno2016 .

We evaluated the matrix elements $\left\langle m\left|{\cal C}\right|m\right\rangle$ and obtain

[TABLE]

where the two numerical constants are given as $\kappa\approx 3.199$ and $\kappa^{\prime}\approx 4.296$, see also Fig. 2. This behavior is asymptotically exact at large $m$ but valid with good accuracy already for $m>2$. The most important aspect of this result is that the dependence of the scattering rate on the angular momentum $m$ is non-analytic. To simplify the analysis we assume in the following that $\tau_{m}^{-1}=\tau_{L}^{-1}\left|m\right|$, where $\tau_{L}^{-1}=\kappa\tau_{1}^{-1}$ is the characteristic time of the of the Lévy flight process, given in Eq. (3).

The implication of the $\left|m\right|$-dependence of $\tau_{m}^{-1}$ becomes evident if we consider the scattering between two distinct momentum directions. Fourier transformation of $\tau_{m}^{-1}$ yields:

[TABLE]

Thus, we obtain a slowly-decaying powerlaw $\sim\left(\theta-\theta^{\prime}\right)^{-2}$ for scattering processes away from forward scattering. Using this result for $\left\langle\theta\left|{\cal C}\right|\theta^{\prime}\right\rangle$ we can rewrite the Boltzmann equation in the form Eq. (2) with characteristic time $\tau_{L}$ of Eq. (3) for the Lévy flight. To arrive at Eq. (2) we used that the convolution of the distribution function with $\left\langle\theta\left|{\cal C}\right|\theta^{\prime}\right\rangle$ can be expressed as a fractional derivative

[TABLE]

a special case of the Riesz-Feller derivative $\triangle^{\mu/2}$Samko1993 ; Hermann2014 .

There are some profound implications that this fractional Fokker-Planck formulation immediately reveals. For example, we consider a scenario where we inject a highly directed excitationLedwith2019 . To this end we consider a source term in the Boltzmann equation that causes this excitation:

[TABLE]

We assumed that we will only inject excitations in a window $\pm k_{B}T$ near the Dirac point, hence the factor $f_{k}^{\left(0\right)}\left(1-f_{k}^{\left(0\right)}\right)$. In addition we decomposed the source term into its angular momentum modes. The linearized Boltzmann equation is applicable if $\left|\delta h_{\lambda m}\right|\ll 1$. To describe an excitation that is peaked along an axis given by a certain momentum direction, we use $\delta h_{\lambda,m}=\delta h\lambda^{m}$, which has a $\lambda$-dependence of the $s=1$ mode of Eq. (9). The solution of the fractional Fokker-Planck equation for a homogeneous case $\mathbf{q}=0$ is then given as

[TABLE]

This function is known as wrapped Cauchy distribution with circular variance $1-e^{-t/\tau_{L}}$Levy1939 ; Mardia1999 . $\Theta\left(t\right)$ is the step function. $\psi_{+}\left(\theta,t\right)$ is shown in the upper panel of Fig. 3.

For $t=0$, $\psi_{\lambda}\left(\theta,t\right)$ corresponds to two delta functions due to particle and hole flows in opposite directions. Let us concentrate on the particle channel $\lambda=+1$. For short times $t\ll\tau_{L}$, the peak in the initial current direction decays as

[TABLE]

while the distribution function grows linearly for all non-zero angles:

[TABLE]

The same behavior occurs for $\lambda=-1$ if we shift $\theta\rightarrow\theta+\pi$. This behavior in contrast to the one that follows from usual Fokker-Planck diffusion. The latter we obtain for example from collision rates $\tau_{m}^{-1}\sim m^{2}.$ Then the usual spreading of a Gaussian wave package occurs with $\psi_{+}\left(t,\theta=0\right)\propto t^{-1/2}$ and $\psi_{+}\left(t,\theta\neq 0\right)\propto t^{2}$ (lower panel of Fig. 3). While the forward direction of a Levy flight decays more slowly than in usual diffusion, the growth at larger angles is much faster, hence the name superdiffusion.

A tangible implication of this superdiffusive charge motion is the heating of the system after the injection. To this end we determine the time dependence of the entropy density

[TABLE]

The heat density caused by the injection is given by $\delta q\left(t\right)=T\left(s_{{\rm eq}}-s\left(t\right)\right).$ Inserting the distribution function of Eq. (16) we obtain

[TABLE]

where $s_{{\rm eq}}$ is the equilibrium entropy density. In order to stay within the regime of linear response, we are confined to $t>\delta h\tau_{L}$. For $t\rightarrow\infty$ one finds $s\rightarrow s_{{\rm eq}},$ and we obtain $s\left(t\right)=s_{{\rm eq}}\left(1-\left(\delta h\right)^{2}\frac{4\log(2)}{9\zeta\left(3\right)}\left(\coth\left(\frac{t}{\tau_{L}}\right)-1\right)\right)$. Thus, initial heating occurs according to

[TABLE]

This result is a direct consequence of the superdiffusive behavior, in particular of the slow decay along the forward direction. In case of ordinary diffusion follows $\delta q\left(t\right)\propto t^{-1/2}$ which is much faster (see Fig. 3). The $m$-dependence of $\tau_{m}^{-1}$ that is responsible for the Lévy flight behavior can also be seen in non-local transport coefficients since the conductivity at finite momentum $\mathbf{q}$ couples the different harmonics of the distribution function. As an example we show in the supplementary material the transverse optical conductivity at finite $\mathbf{q}$. Nevertheless, experiments with directed electron beams Wang2019 , which in the past have been used to investigate electron-electron scattering effects Predel2000 , seem to offer a more direct way of testing the short time behavior of Eq. (21).

The occurrence of Lévy flights to describe scattering processes in momentum space is a more general phenomenon and not restricted to graphene at the neutrality point. In two-dimensional Fermi liquids with characteristic rate $\hbar\tau_{FL}^{-1}\sim k_{B}T^{2}/T_{F}$, it holds for $\left|m\right|<M\sim\sqrt{T_{F}/T}$ that \tau_{m}$${}^{-1}\sim\tau_{FL}^{-1}\frac{m^{p}}{M^{p}}\log\left|m\right| with $p=2\left(1+\left(-1\right)^{m}\right)$, while $\tau_{m}^{-1}\sim\tau_{FL}^{-1}$ for $\left|m\right|>M$Gurzhi1996 ; Ledwith2019 . $T_{F}$ is the Fermi temperature. This yields superdiffusive behavior in a wide time window. Another system that also shows $\tau_{m}^{-1}\propto\left|m\right|$ for arbitrarily large $m$ consists of electrons in a random magnetic field, important for the description of composite fermions in the fractional quantum Hall regimeMirlin1997 . Our analysis implies that this system should also undergo a wrapped Cauchy flight in momentum space. Large classes of quantum-critical systems, discussed e.g. in Refs.Hertz1976 ; Moriya1984 ; Millis1993 ; Laughlin2001 ; Millis2002 ; Abanov2003 ; Loehneysen2007 ; DellAnna2007 ; Metlitski2010 ; Schattner2016 ; Damle1997 are governed by long-ranged soft-mode interactions. An analysis of collision processes along the lines discussed here may reveal a non-analytic dependence of the scattering rates on angular momentum quantum number according to $\tau_{m}^{-1}\propto\left|m\right|^{\mu/2}$. This would give rise to a more general class of wrapped Lévy flights, a consequence of the power-law behavior $\left\langle\theta\left|{\cal C}\right|\theta^{\prime}\right\rangle\propto\left|\theta-\theta^{\prime}\right|^{-1-\frac{\mu}{2}}$ near forward scattering. This could occur on the Fermi surface for itinerant quantum critical systems or near a soft momentum in critical bosonic systems. If a fractional Fokker-Planck formulation, along the lines of our Eq. (2), can be derived, it will be significantly easier to draw conclusions about the out-of-equilibrium dynamics of the system such as a focussed injection of collective excitations. Finally we mention that the formulation of the Boltzmann equation presented here can also be used to study the non-local electric and thermal conductivities and viscosities, allowing insight into the diffusive and sound excitations in the hydrodynamic regimeKiselev2019b .

Acknowledgements: We are grateful to Andrey V. Chubukov, Leonid S. Levitov, Dimitrii L. Maslov, and Alexander D. Mirlin for stimulating discussions and to the European Commission’s Horizon 2020 RISE program Hydrotronics for support.

Supplementary material

I Summary of the Simulations shown in Fig.1

Superdiffusion on the Dirac cone can be seen as a random walk of particles, where the step sizes are distributed according to a wrapped Cauchy distribution. This distribution solves the fractional Fokker-Planck equation (2) of the main text (see Chechkin2004 ; Metzler2012 ). The anglular distance on the Dirac cone travelled by an electron during a time interval $\Delta t$ is therefore distributed according to

[TABLE]

To generate Figs. 1 a) and b) of the main text we created a sequence of random angles $\Delta\theta_{i}$ using the distribution (22). The position of the particle after $N$ steps, i.e. after a time intervall $N\Delta t$, then is

[TABLE]

Fig. 4 depicts a wrapped Cauchy random walk with $N=500$ steps.

In the case of ordinary diffusion, the step size distribution of Eq. (22) must be replaced by a wrapped normal distribution, which is written in terms of the Jacobi theta function, but can be closely approximated by the van Mises distribution (see e.g. Ref. Kurz2014 ):

[TABLE]

Using the described procedure we obtain the wrapped Gaussian random flight shown in Fig. 5.

II Collision operator due to electron-electron Coulomb interaction

We briefly summarize the main steps in deriving the collision operator of the Boltzmann equation used in this paper. The collision operator is determined from the larger and smaller self energies on the Keldysh contour. For further details, see Ref. Fritz2008-1 .

The non-interacting part of the Hamiltonian is

[TABLE]

which is diagonalized by the unitary transformation

[TABLE]

with

[TABLE]

The eigenvalues of the Hamiltonian are $\pm v\hbar k$. Thus we obtain quasiparticle states for the two bands: $\gamma_{\mathbf{k}}=U_{\mathbf{k}}\psi_{\mathbf{k}}$ with

[TABLE]

The electron-electron Coulomb interaction is

[TABLE]

with $V\left(q\right)=\frac{e^{2}}{2\pi\varepsilon\left|\mathbf{q}\right|}$. Transforming the interaction into the band, or helical representation, which takes into account the locking between momentum and pseudo-spin that originates from the two sub-lattice structure of graphene. It follows

[TABLE]

where

[TABLE]

Within second order perturbation theory it holds for the self energies for occupied and unoccupied states, respectively.

[TABLE]

$N$ combines the valley and spin degrees of freedom and takes the value $N=4$. Next we use the fact that within a quasiparticle description the upper and lower propagators are expressed in terms of the distribution functions $f_{\lambda}\left(\mathbf{k},\mathbf{r},t\right)$ as

[TABLE]

As usual, $\mathbf{k}$ and $\omega$ stand for the Fourier-transformed variables of the relative coordinates and times while $\mathbf{r}$ and $t$ stand for the center of gravity or mean time.

The collision operator can now we determined from the self energies $\Sigma^{<}$ and $\Sigma^{>}$:

[TABLE]

For simplicity we only keep the momentum $\mathbf{k}$ and band index $\lambda$. Inserting $G^{>}$ and $G^{<}$ into the self energies yields with the linearization Eq.(6) of the main paper the result for the collision operator given in Eq.(8) of the main paper. The matrix elements $\gamma_{\mathbf{k},\mathbf{k}^{\prime},\mathbf{q}}^{\left(1,2\right)}$ of that equation are given as:

[TABLE]

with

[TABLE]

and

[TABLE]

Upper-case letters like $K=k_{x}+ik_{y}$ etc. combine the two components of the momentum onto a complex variable.

From the same unitary transformation also follows that

[TABLE]

This can be used to analyze the current

[TABLE]

of Dirac particles which consists of intra- and inter-band contributions:

[TABLE]

The two terms are given as

[TABLE]

Thus, the velocity used in our Eq. (4) of the main paper is precisely the expression $\mathbf{v}_{\mathbf{k}\lambda}=v\frac{\lambda\mathbf{k}}{k}$ of the intraband current $\mathbf{j}_{\mathrm{intra}}$. Spin-momentum locking is included naturally, if one goes to the helical states of the upper and lower Dirac cone. The hydrodynamic response is governed by strong collisions of intraband excitations.

III Identification of the collinear modes at finite angular momentum

In this section we determine the zero modes of the collision operator of Eq.(8) if we confine ourselves to collinear collision processes. To this end we need to find under what conditions the two expressions

[TABLE]

that occur in Eq.(8), vanish. Here we have to include the additional constrain

[TABLE]

that follows from energy conservation.

By collinear modes we mean that all involved momenta are either parallel or antiparallel. As discussed in Ref.Fritz2008-1 we consider such zero modes of collinear processes because all other processes are suppressed by $1/\log\left(1/\alpha\right)$ where $\alpha$ is the fine-structure constant. Of course, the analysis allows for scattering processes that are not collinear; the issue is merely that distribution functions that become zero modes for collinear scattering are enhanced relative to those that are no such zero modes. Finally we comment that the main conclusion of our paper, namely that the scattering rate depends on angular momentum like $\tau_{m}^{-1}\propto k_{B}T\left|m\right|$, is unchanged if we go beyond the collinear mode regime.

One immediately finds that $A^{\left(1\right)}=B^{\left(1\right)}=0$ subject to Eq.41 is obeyed by $\psi_{\mathbf{k},\lambda}=1$, $\psi_{\mathbf{k},\lambda}=\mathbf{k}$, and $\psi_{\mathbf{k},\lambda}=\lambda\left|\mathbf{k}\right|$, regardless of whether we confine ourselves to collinear modes. These modes correspond to the conservation of charge, momentum, and energy, respectively. In addition to these modes one also finds $\psi_{\mathbf{k},\lambda}=\lambda$ is a zero mode. It corresponds to the fact that second order perturbation theory does not relax a charge imbalance between the upper and lower Dirac cone.

Next we consider distribution functions

[TABLE]

where $k=\left|\mathbf{k}\right|$ is the magnitude of the momentum and $\theta_{\mathbf{k}}$ its polar angle: $\mathbf{k}=k\left(\cos\theta_{\mathbf{k}},\sin\theta_{\mathbf{k}}\right)$. Collinear scattering corresponds to

[TABLE]

where even or odd $s_{i}$ correspond to parallel and antiparallel momenta relative to $\mathbf{k}$. We first show that all $s_{i}$ are even due to energy conservation. To this end we assume without restriction that $\mathbf{k}=\left(k,0\right)$ with $k>0$. Then $\mathbf{k}^{\prime}=u\left(k,0\right)$ and $\mathbf{q}=w\left(k,0\right)$, where we do not assume that $u$ and $w$ are positive. Energy conservation now implies

[TABLE]

We need to fulfill this condition for an extended set of variables, not just for an isolated point in momentum space. This implies that $1+w>0$ so we can cancel the “$1$” on both sides. Then, to be able to cancel $w$ it must hold that $u>w$, which in turn implies $u>0$ to cancel $u$. Thus, we find that the momenta $\mathbf{k}^{\prime}$, $\mathbf{k}+\mathbf{q}$, and $\mathbf{k}^{\prime}-\mathbf{q}$ point in the same direction as $\mathbf{k}$ even though $\mathbf{q}$ is allowed to point in the opposite direction. It follows that we can assume without restriction that

[TABLE]

If we use that $e^{im\theta_{-\mathbf{p}}}=\left(-1\right)^{m}e^{im\theta_{\mathbf{p}}}$we obtain

[TABLE]

It is now easy so find that there are the following solutions that yield $A^{\left(1\right)}=B^{\left(1\right)}=0$ subject to Eq.41:

[TABLE]

In addition one finds $a_{\lambda,m}\left(k\right)=\lambda\left|\mathbf{k}\right|$ if $m$ is even and $a_{\lambda,m}\left(k\right)=\left|\mathbf{k}\right|$ if $m$ is odd.

Thus, we can write that the following modes are zero modes in the collinear scattering limit

[TABLE]

which is Eq.(9) of the main paper.

IV Superdiffusion and non-local transport

The Fokker-Planck equation (2) of the main text

[TABLE]

can be used to calculate the response of graphene electrons to an external perturbation, such as for example an electric field. In this case the force term is given by

[TABLE]

where $\mathbf{E\left(\mathbf{q},\omega\right)}=E_{0}e^{i\left(\mathbf{q}\cdot\mathbf{x}-\omega t\right)}\mathbf{e}_{x}$ is the electric field. To first order in the electric field it is

[TABLE]

We perform a Fourier transform $t\rightarrow\omega$, $\mathbf{x}\rightarrow\mathbf{q}$ and project the equation (49) onto the collinear zero modes (48) using the scalar product $\Braket{\phi}{\chi}=\sum_{\lambda}\int\frac{d^{2}k}{\left(2\pi\right)^{2}}\phi_{\bm{k},\lambda}\chi_{\bm{k},\lambda}$. The result is a simplified version of the Boltzmann equation:

[TABLE]

where $m$ labels the angular harmonic of the collinear zero mode and $\vartheta_{\bm{q}}$ is the angle of the wave vector $\mathbf{q}$ with respect to the $x$-axis. This equation is exact in the limit of a small coupling constant $\alpha$, where collinear events dominate the electron-electron scattering Fritz2008-1 . Notice, that the electric field only couples to angular harmonics with $\left|m\right|=1$, however for a spatially inhomogeneous field with $q\neq 0$, the second right hand side term of Eq. (50) couples all angular harmonics. Therefore, information on the $m$-dependence of the scattering times can be extracted from the non-local (i.e. $\mathbf{q}$-dependent) electric conductivity $\sigma_{\alpha\beta}\left(\mathbf{q},\omega\right)$, which is defined through

[TABLE]

The conductivity tensor $\sigma_{\alpha\beta}\left(\mathbf{q},\omega\right)$ can be decomposed into longitudinal and transverse parts according to

[TABLE]

where $\sigma_{\parallel/\bot}\left(\omega,q\right)$ only depend on the magnitude of $\mathbf{q}$. Fig 6 shows the influence of the $m$-dependence of the scattering time $\tau_{m}$ on the real part of the transverse conductivity $\sigma_{\bot}\left(\omega,q\right)$. We conclude that the non-local conductivity, playing an important role in experiments on surface acoustic waves Wixforth1986 ; Simon1996 , provides a possibility to detect the peculiar dependence of the scattering times $\tau_{m}$ on $m$, and to confirm the Lévy flight behavior predicted in the main text.

For completeness, we mention that the Boltzmann equation (50) can be solved exactly, and the expressions for the non-local conductivities can be written down in closed form EK2019 :

[TABLE]

Here, $\sigma_{0}=N\frac{e^{2}\log\left(2\right)\tau_{\sigma,1}}{2\pi\beta\hbar^{2}}$ is the quantum critical conductivity calculated in Ref. Fritz2008-1 and $M\left(q,\omega\right)$ is the memory function summerizing the effects of higher angular harmonics:

[TABLE]

where $I_{\nu}\left(z\right)$ is the modified Bessel function.

References

• (1) A. V. Chechkin, V. Y. Gonchar, J. Klafter, R. Metzler, L. V. Tanatarov, J. Stat. Phys. 115, 1505 (2004).
• (2) R. Metzler, A. V. Chechkin, J. Klafter, Levy Statistics and Anomalous Transport, Computational Complexity, Springer,

• (3) G. Kurz, I. Gilitschenski and U. D. Hanebeck, "Efficient evaluation of the probability density function of a wrapped normal distribution" in 2014 Sensor Data Fusion: Trends, Solutions, Applications (SDF), IEEE (2014).
• (4) L. Fritz, J. Schmalian, M. Müller and S. Sachdev, Quantum critical transport in clean graphene, Phys. Rev. B 78, 085416 (2008).
• (5) A. Wixforth, J. P. Kotthaus, and G. Weimann Phys. Rev. Lett. 56, 2104 (1986).
• (6) H. S. Simon, Phys. Rev. B 54, 13878 (1996).
• (7) E. I. Kiselev, J. Schmalian, unpublished.