Quantum corrections to transport in graphene: a trajectory-based semiclassical analysis

We review a calculation of the quantum corrections to electrical transport in graphene using the trajectory-based semiclassical method. Compared to conventional metals, for graphene the semiclassical propagator contains an additional pseudospin structure, which influences the results for weak localization, and interaction-induced effects, such as the Altshuler-Aronov correction and dephasing. Our results apply to a sample of graphene that is doped away from the Dirac point and subject to a smooth disorder potential, such that electrons follow classical trajectories. In such system, the Ehrenfest time enters as an additional timescale.


Introduction
The discovery of a method to isolate single layers of graphene [1] has created a field of intense research over past years (see Refs. [2][3][4][5][6] for reviews). The remarkable electronic properties of graphene stem from a quasirelativistic dispersion with a fourfold degeneracy due to spin and valley degrees of freedom, that is characteristic for the underlying two-dimensional honeycomb lattice [7]. For pristine graphene, the Fermi level lies at the Dirac point, where the density of states vanishes. Electric transport is then dominated by evanescent modes, giving rise to a finite minimum conductivity 4e 2 /πh [8][9][10][11][12][13][14]. Remarkably, when smooth disorder is added, that does not couple the valleys, no Anderson localization takes place [15][16][17][18], but instead the conductivity rises with disorder strength [15,16,19].
When graphene is doped away from the Dirac point, its physical properties resemble those of a metal, but certain intriguing features from the Dirac spectrum remain. The reason for this lies in the "pseudospin" degree of freedom, which is connected to the two-atom basis of the two-dimensional honeycomb lattice. Sufficiently close to the Dirac point, where the electronic dispersion is still linear, the direction of the pseudospin is aligned with the momentum. This helicity of charge carriers strongly influences the electronic properties. Two important consequences are the absence of backscattering at a potential barrier (Klein tunneling) [20,21] and the half-integral quantum Hall effect [22,23].
The type of disorder [24,25] and the dielectric environment [26][27][28] play an important role for the transport properties of graphene. When charged impurities in the substrate are the main source of disorder, the disorder potential is smooth on the scale of the lattice constant. For graphene doped away from the Dirac point, one may eventually reach a situation, where the spatial scale ξ on which the disorder changes is much larger than the Fermi wavelength λ F . (Strictly speaking, this condition is met only if the graphene sheet is embedded in an insulating medium with a high dielectric constant, such as HfO 2 [29][30][31].) In this limit, it is justifiable to work with the semiclassical approximation, and to utilize classical trajectories for the calculation of physical quantities. In the past years, trajectory-based semiclassical methods have been successfully applied to the calculation of quantum corrections to the transport of normal metallic systems (see, e.g., [32][33][34][35][36][37]). The present article discusses the application of such methods to quantum transport in graphene.
In the limit that the Fermi wavelength λ F is much smaller than the transport mean free path l tr , the main contribution to the conductivity is given by the Drude conductivity, for which the quantum phase coherence of the electrons plays no role [28,38]. Corrections to the Drude conductivity are called "quantum corrections". They are the result of quantum interference of electrons moving along different trajectories. The two quantum corrections to the conductivity are the weak localization correction, which results from interference of time-reversed trajectories [39][40][41][42][43], and the Altshuler-Aronov correction, which is interaction-induced and originates from elastic scattering of electrons on Friedel oscillations of the electron density [44][45][46]. Inelastic interaction processes cause a loss of phase coherence, and set an upper limit on the timescale at which weak localization occurs [47][48][49]. Although the quantum corrections are a large factor ∼ l tr /λ F smaller than the Drude conductivity, they can be experimentally detected by their characteristic dependence on magnetic field and temperature. The measurement of such quantities provides important information about the type of disorder and interactions for a specific sample.
Soon after the initial experiments, the theory of quantum corrections in disordered metals was extended to graphene [25,[50][51][52][53][54][55]. Of particular interest was the potential valley-mixing effect of short range disorder (correlation length ξ λ F ), which leads to a transition between weak localization and weak antilocalization [25,[50][51][52] and strongly affects the magnitude of the Altshuler-Aronov correction [54][55][56]. Such transitions were also observed experimentally [54][55][56][57][58][59][60][61][62][63]. In the present article, we consider the quantum corrections to the conductivity of graphene in the presence of a long-range impurity potential, which is smooth on the scale of the Fermi wavelength λ F . Such a smooth random potential does not mix the valleys, but instead leads to a number of modifications of the quantum corrections because of its smoothness. Such modifications are known from the theory of conventional electron gases [32-34, 37, 64]. It is the goal of the present article to extend and collect those results for the case of graphene.
In order to understand why a smooth random potential modifies quantum corrections, one first recalls that electrons follow well-defined classical trajectories if the potential is smooth. The randomness of the potential then ensures that the classical dynamics is chaotic: Two nearby trajectories separate exponentially with time, the divergence rate being described by a Lyapunov coefficient λ. For quantum interference corrections this exponential divergence then leads to the notion of the "Ehrenfest time" τ E [32]: The Ehrenfest time is defined as the time after which two classical trajectories, initially a quantum distance λ F apart, get separated by a reference distance L c characteristic of the classical dynamics The significance of the Ehrenfest time is that it puts a short-time cutoff for the occurrence of quantum corrections, as it sets the minimum time for interference effects to take place. In the case of a non-smooth potential with variations on the scale of the Fermi wavelength, no such short-time cutoff appears, and the results of the semiclassical theory agree with those of standard diagrammatic methods. The trajectory-based semiclassical calculation of quantum corrections to the conductivity of graphene differs from the same calculation for conventional metals by the additional pseudospin structure. The problem of extending the semiclassical formalism to system with a spinor degree of freedom, such as metals with spin-orbit coupling or Dirac Hamiltonians, has received considerable attention in the literature [65][66][67][68][69][70][71][72]. The application of the formalism to the case of graphene by Carmier and Ullmo [73] will serve as a starting point for our calculation. A key element in the trajectory-based approaches is that the pseudospin can be reconstructed along the trajectories, where it remains aligned with the momentum. Associated with the transport of the pseudospin along the trajectory is an additional phase in the semiclassical propagator, which can be identified as the Berry phase. One example where this phase plays an important role is the semiclassical calculation of the Landau levels, where the electrons acquire a Berry phase of π during the cyclotron motion, ultimately leading to the half-integral quantum Hall effect.
The semiclassical theory presented here is specifically aimed at the leading order quantum corrections (weak localization and Altshuler-Aronov correction, as well as the effect of dephasing on weak localization) for graphene in a smooth random potential. Typical systems to which the trajectory-based semiclassical method has been applied in the literature are quantum billiards, ultraballistic systems, where particles scatter only at the boundary of the sample, or antidot arrays, high mobility two-dimensional electron gases with artificially superimposed antidots, that act as classical scatterers [64,74,75]. While for standard semiconductor structures, quantum billiards can be shaped by means of gate potentials, such procedure is problematic for graphene, as it is a gapless material. Quantum billiards in graphene can be realized in etched structures, where the edges are atomically sharp [76][77][78][79]. The scattering on such edges then depends on the precise atomic configuration, and deserves a careful consideration [80,81]. The same applies to antidot arrays in graphene [82], but also here, the boundaries of the antidots are so sharp, that they lead to scattering between the valleys. Such atomically sharp boundaries invalidate a description that is solely based on classical trajectories, although for sufficiently well-defined boundaries a theoretical description involving coupled valleys is possible [80,81]. Such a limitation does not exist for the effect of an impurity potential in a graphene sheet on a substrate with high dielectric constant, which is the scenario we consider here. In this case, the high dielectric constant ensures that the screening length is larger than λ F at sufficiently high doping. Hence, the potential has a correlation length ξ λ F and classical paths are well-defined objects. For the sake of readability, we have tried to make this article self contained, only referring to the literature for the finer technical details of the semiclassical approach. Starting point of our discussion is the semiclassical Green function, that will be introduced in Sec. 2. To set the stage, we then calculate the Drude conductance in Sec. 3. We then turn to the quantum corrections, where the weak localization is calculated in Sec. 4, and the interaction-induced corrections are treated in Sec. 5 (Altshuler-Aronov) and 6 (dephasing). We conclude in Sec. 7.

Semiclassical Green function
In this article, we consider graphene subject to a smooth disorder potential V (r) that does not couple the valleys. In the vicinity of the Dirac point, electrons are described by the Hamiltonian where v F is the Fermi velocity, µ is the chemical potential, and σ = (σ x , σ y ) are Pauli matrices for the pseudospin degree of freedom. The eigenvalues of the kinetic energy term v F p · σ are K ± = ±v F |p|. We will be interested in the case of electron-doped graphene for which the chemical potential µ is larger than the potential V (r). In this case we may restrict our attention to the conduction band and set K = v F |p|. The corresponding eigenspinor of the kinetic energy is where the angle φ p denotes the direction of the momentum p.
Starting point for our semiclassical analysis of transport is the semiclassical expression for the retarded Green function G R (r, r ; ε) at energy ε derived by Carmier and Ullmo [73], where the summation is over all classical trajectories α that connect the points r and r , S α is the classical action of the trajectory, A α the stability amplitude, γ α an additional phase shift to be defined below, and p α and p α are the initial and final momenta of the trajectory α, respectively. Equation (4) generalizes the corresponding expression for a system without spin or pseudospin degrees of freedom [83]. (In that case the projection factor |χ(p α ) χ(p α )| and the phase shift γ α are absent.) The classical trajectories are determined by the classical Hamilton function and the classical action of a trajectory satisfies the equations where τ α is the duration of the trajectory. The stability amplitude A α is found as Finally, the phase shift γ α contains the Berry phase where in the second equality we chose the initial and final angles to be 0 ≤ φ p α , φ pα ≤ 2π and add 2πn, n being integer, to account for the phase winding along the path α. We have not included the phase corresponding to the Maslov index [83], which can be disregarded for the calculation of transport.
In the next Sections we also need the advanced Green function, which follows from the relation

Drude conductance
We now turn to the calculation of the conductivity. Hereto, we consider a rectangular sample of graphene of dimensions L × W , calculate its conductance G, and obtain the conductivity σ from the relation G = σW/L. The conductance G is calculated from the Kubo formula where f (ε) = 1/(e ε/T + 1) is the Fermi function and d g = 4 denotes the degeneracy due to spin and valley. Further, the velocity operator for graphene readŝ and the trace indicates a summation over pseudospin indices. For a semiclassical calculation of the conductance, we insert the semiclassical Green function (4) into the Kubo formula, so that G is expressed as a double sum over trajectories α and β. Restricting the summation to diagonal terms α = β, the so-called diagonal approximation, then gives the Drude conductance. For α = β the semiclassical approximation Eq. (10) contains matrix elements of the form with v x = ∂H cl /∂p x . Apart from the factor d g , the resulting expression is the same as in the case of a standard two-dimensional electron gas, with the initial (final) classical velocity v x (v x ). The remaining summation over trajectories α can be transformed to an integral over initial and final momentum and the duration of the trajectories [84,85], Here, f is an arbitrary function, and the integration over momenta is restricted to the energy shell dp ε = dpδ(ε − H cl (r, p)). The trajectory density ρ ε (X → X; t) selects only those phase space points X = (r, p) that are connected via a classical trajectory with the initial point X . Since the initial phase space point together with the classical Hamilton function uniquely determines the classical trajectory, ρ is expressed as a δ-function where X(X , t) = (r(r , p ; t), p(r , p ; t)) is the phase space point that a trajectory starting out of X = (r , p ) reaches after a time t. The calculation proceeds by performing a statistical average of the conductance, where we average over small displacements in the disorder potential V (r). Such procedure replaces the exact trajectory density ρ ε (X → X; t) by the classical propagator P (X, X ; t), which is a smooth function of initial and final phase space coordinates as well as time. The classical propagator has only a weak dependence on energy which we neglect in the following. We then obtain for the Drude conductance For the further analysis of the Drude conductance, one needs to specify the classical propagator P (X, X ; t). Following [15,86], we consider a random Gaussian potential with the correlation function where ξ is the correlation length, and K 0 is the dimensionless strength of the potential. On spatial scales much longer than the correlation length, the electronic motion becomes diffusive, with a diffusion coefficient that we will calculate in the following. The random potential of Eq. (17) has been used to describe the impurity potential in experiments [87,88]. We start by considering a particle that moves in the x-direction at time t = 0 and consider how the direction of motion changes under the influence of the potential V (r). Following the classical dynamics, the angle φ(t), at which the electron propagates at time t is given by where k F is the Fermi wavenumber. This equation is valid for times t short enough, such that |φ(t)| 1, so that the motion of the electron is mainly along the x-direction, r(t) = v F te x . We then find for the mean quadratic deflection Using the correlation function Eq. (17), this gives provided t is much longer than the "correlation time" t ξ = ξ/v F . Our derivation required the time t to be short enough such that the deflection is small. Such time interval exists, as long as K 0 (k F ξ) 2 , which is a condition that can be met if the disorder is smooth on the scale of the Fermi wavelength (k F ξ 1). Equation (20) describes a linear-in-time increase of the quadratic deflection, a characteristic property of the diffusive motion for the angle φ. Continuing the diffusive process beyond small angles, the validity of this equation can be extended to all times longer than the correlation time t ξ . Further, extending the result to arbitrary starting times and arbitrary directions in the beginning of the propagation, we conclude that the angle difference φ(t) − φ(t ) has a Gaussian distribution with zero mean and with variance Now we can calculate the electron's mean square displacement. Since the mean square displacement is given by The average can be performed using the Gaussian distribution of φ(t ) − φ(t ) and one finds, in the long-time limit, where the diffusion constant D is given by The diffusion constant of Eq. (25) corresponds to a transport mean free path which is parametrically larger than the correlation length ξ in the limit k F ξ 1 of a smooth potential.
We now turn back to the evaluation of Eq. (16) in the diffusive limit. In this case, the direction of the momentum does not influence the propagation, which amounts to the replacements and where ν is the density of states per spin and valley. The relevant classical propagator P (r, r ; t) is then a function of position only, and satisfies a diffusion equation where D is the diffusion constant, determined above. The solution to the diffusion equation Eq. (29) in two dimensions, with perfect leads at x = 0 and x = L, and insulating boundaries at y = 0 and y = W is given by where we introduced the function The summation over q extends over q x = nπ/L with n = 1, 2, ... and q y = mπ/W with m = 0, 1, .... In our calculation of the Drude conductance, as well as in the forthcoming calculation of the quantum corrections, we also encounter expressions, where the diffusion propagator that connects to the leads is multiplied by the velocity v x . Such expressions are handled with the help of the diffusive flux at position r and time t, for a particle starting from r at time t = 0. We find The quantity P LR appears in the Drude conductance, while the quantities P L and P R , which represent the probability that a particle at position r has entered via the left lead or will exit via the right lead, respectively, are introduced for later use. For the Drude conductivity σ 0 = G 0 L/W we then obtain the standard expression where the factor d g = 4 accounts for the degeneracy for spin and valley. Importantly, pseudospin does not enter as an additional degeneracy, since it is locked to the momentum. Taking the expression for the diffusion coefficient for the Gaussian random potential, Eq. (25), as well as the density of states at graphene, ν = k F /2π v F , one obtains The same result was obtained in a quantum-mechanical calculation using the Boltzmann equation in [86].

Weak antilocalization
Deviations from the Drude conductance are termed quantum corrections. Without interactions and for conventional metals, the leading correction to the classical conductance results in a small reduction of the conductance, and is called "weak localization", since it describes the onset of Anderson localization. In graphene, the Berry phase is responsible for a different sign of this quantum correction, which gives rise to an enhanced conductance (when effects of intervalley scattering and trigonal warping are neglected), and is therefore called weak antilocalization [50-52, 58, 60]. In the following, we will show how the weak antilocalization correction is derived in the semiclassical formalism, and discuss the effect of a finite Ehrenfest time. Again, we will give explicit results for the case of a smooth random Gaussian-correlated potential.
In the semiclassical framework, weak (anti)localization results from configurations of retarded and advanced trajectories α and β as shown in Fig. 1. The trajectories can be divided into four segments: The entrance and exit segments, where the trajectories α and β are correlated or "paired" -i.e., the difference between the two trajectories is sufficiently small, that the chaotic classical dynamics can be linearized on that scale -, the loop segment, where the trajectory α is paired with the time-reversed of trajectory β, and the encounter region (or Lyapunov region), where trajectories α and β as well as their time-reversed are correlated. At the beginning of their first passage through the encounter region, the trajectories α and β are located within a Fermi wavelength λ F . Due to the chaotic motion, this phase-space distance increases exponentially along the encounter region as d(t) = λ F e λt , where λ is the Lyapunov coefficient characteristic of the chaotic motion. For the random Gaussian potential (17), the Lyapunov exponent is see [32] and Appendix Appendix A. At the end of the encounter region, the distance has reached a classical size L c , beyond which classical motion is considered uncorrelated -i.e., the classical dynamics can no longer be linearized. For the smooth random potential (17), we may identify L c ξ. The duration of the encounter is set by the Ehrenfest time Our final results can be expressed in terms of Ehrenfest time only, which depends logarithmically on L c , so that a more precise definition of the cutoff L c is not needed. For the calculation of the weak localization, one starts from the Kubo formula, Eq. (10), inserts the semiclassical expressions for the Green function, and then restricts the summation to configurations of trajectories as explained in the previous paragraph. As long as the duration of the encounter region is τ E or larger, the trajectories α and β acquire an action difference ∆S [32,35]. For graphene, we also have to keep track of the influence of the pseudospin, which has two effects: First, the spinor structure of the semiclassical Green function changes the velocity operator to the classical velocity, in the same way as before for the Drude conductance, see Eq. (12). Second, since the trajectories are no longer equal, we have to pay attention to the Berry phase collected along the trajectories α and β. At this stage, we can write where the summation is restricted to the configurations of trajectories shown in Fig. 1, The difference of the Berry phase ∆γ = γ α − γ β , is collected in the loop segment only. (In the encounter region, the trajectories α and β differ on a sub-macroscopic scale only, which adds a negligible contribution to the Berry phase difference.) Since the momenta of trajectory α are opposite in the beginning and the end of the loop segment, we have from Eq. (8) with integer n depending on the total winding of the momentum along the trajectory.
Since the Berry phase is expressed as integral along the trajectory, the Berry phase collected by trajectory β along the loop is just γ β,Loop = −γ α,Loop . Hence, we find for all configurations of trajectories contributing to the quantum correction. This minus sign is responsible for the change from weak localization (in conventional twodimensional electron gases without spin-orbit coupling) to weak antilocalization. The remaining calculation then proceeds as in the standard case, and we find [37] where P (r , r ; t) is the diffusion propagator, see Eq. (30), and P L (r) and P R (r) are defined in Eq. (33). For the further evaluation of Eq. (41), we write and perform two partial integrations on r. Making use of the explicit form of P L (r) and P R (r), we arrive at In two dimensions, the time integral in this equation is divergent for small times, and the appropriate cutoff is set by the transport time τ tr , below which the diffusive approximation breaks down. In the limit of large aspect ratio W/L, and small τ tr /τ D , where τ D = L 2 /Dπ 2 is the dwell time, we then find (see Appendix Appendix B) where the function h(x) is defined as It has the asymptotic behavior At zero Ehrenfest time, h(0) = 1 and we arrive at the well-known result for weak antilocalization of a symplectic metal [41,42]. At finite Ehrenfest time, our calculation results in a suppression of the weak antilocalization by the additional factor h(τ E /τ D ), shown in Fig. 2. (The same multiplicative factor h(τ E /τ D ) describes the suppression of weak localization or weak antilocalization in a conventional two-dimensional electron gas. We are not aware of a calculation of the function h in this context.)

Altshuler-Aronov correction
We now turn to the effects of interactions on the conductivity. Interactions modify the conductivity in two physically distinct ways. First, interactions cause the so-called Altshuler-Aronov correction [45,46], which has its origin in the interference of elastic scattering off impurities and off Friedel oscillations of the electron density around an impurity. Second, inelastic electron-electron scattering is responsible for a loss of phase coherence or "dephasing", which sets an upper limit on the time at which weak (anti)localization can occur.

Lowest-order interaction correction
The semiclassical treatment of interaction corrections proceeds via two steps. First, one considers a specific random potential and includes interactions to first order diagrammatic perturbation theory. Such procedure is rather standard, and results in expressions in terms of Green functions G(r, r ; ε) for the given disorder realization. The second step is to take the disorder average, where we employ the semiclassical framework. Hereto, we insert the semiclassical expressions for the Green functions and identify the relevant configurations of trajectories that contribute to the interaction corrections. Our results take into account the effects of a finite Ehrenfest time.
We now outline the calculation in more detail. For the inclusion of interactions using diagrammatic perturbation theory, we can adapt the results of [48], where the authors divide their expressions into contributions to the Altshuler-Aronov correction and to dephasing effects. For the calculation of the conductance, as done in this work, we may further simplify the expressions of [48] by keeping only terms, where one retarded and one advanced Green function is attached to the current vertex (in contrast to the calculation of the conductivity, where also diagrams with two Green functions of the same kind connected to the current vertex need to be kept, see discussion in [34]).
Explicit expressions for the Altshuler-Aronov conductance correction δG AA are obtained from Eq. (10) upon replacing the retarded Green function G R (r, r ; ε) by G R (r, r ; ε) + δG R,F (r, r ; ε) + δG R,H (r, r ; ε), and a similar replacement for the advanced Green function G A (r, r ; ε), keeping terms to first order in the interaction only. The functions δG R,F (r, r ; ε) and δG R,H (r, r ; ε) are Fock and Hartree corrections to the singleparticle Green function, respectively, In these expressions we wrote the pseudospin indices explicitly. We further allow for a frequency dependence of the interaction propagator U R (r 1 , r 2 ; ω) to include the effect of dynamical screening. To first order in interaction we also obtain additional corrections to the conductance which are contributing to dephasing. These will be discussed in the next section. Insertion of the semiclassical Green functions leads to a sum over four trajectories. Systematic contributions to δG AA are obtained only if trajectories originating from retarded and advanced Green functions are paired up or if they undergo a small-angle encounter [34]. Configurations of trajectories that are in line with these requirements for the Fock contribution to δG AA are shown in Fig. 3. Here, configuration (a) originates from a term with three advanced Green functions and one retarded Green function. In this case, the three "advanced" trajectories must join to a single trajectory that can be paired up with the "retarded" trajectory. Configurations (b)-(e) correspond to a term with two retarded and two advanced trajectories. In this situation, due to the specific requirements on start and end point of the Green functions, the trajectories cannot be paired one by one, instead the four trajectories undergo a small-angle encounter. The subdivision into configurations (b)-(e) reflects the possibilities to have none, one or both interaction points within the encounter region. For each one of the configurations shown . For graphene, the pseudospin structure contributes additional factors associated with the interaction vertex. In our notation, the retarded Green function G R (r 2 , r 1 ) is associated with trajectories running from r 1 to r 2 , while the advanced Green function G A (r 2 , r 1 ) is associated with trajectories running from r 2 to r 1 . This amounts to the following possibilities: If the interaction vertex is associated with two Green functions of the same kind, then one trajectory is pointing towards the vertex, while the other one is pointing away. For an interaction vertex associated with one retarded and one advanced Green function, both trajectories either point towards or away from the vertex. In the figure, we show four possibilities that result in a factor χ(p)|χ(p ) , with the associated labelling of momenta.
in Fig. 3, there is a counterpart for which the role of retarded and advanced trajectories is interchanged.
In close analogy to the calculation for conventional electron gases [34], we find, for a random potential with Gaussian correlations as in Eq. (17), with the kernel K(r 1 , r 2 ; ω) given by where is the Fourier transform of the diffusion propagator P (r, r ; t), and where the brackets . . . p 1 ,p 2 denote an average over the directions of the momenta p 1 and p 2 . Equation (50) demonstrates that the Ehrenfest time is acting as a short-time cutoff for the Altshuler-Aronov correction. The spinor structure of the semiclassical Green function contributes the factor Σ F (p 1 , p 2 ), which is not present in the calculation for the conventional two-dimensional electron gas. As explained in Fig. 4, it depends on the overlap of pseudospinors of the two trajectories at the interaction points. For the diffusive motion, only the quantity averaged over momenta p 1 and p 2 with |p 1 | = |p 2 | = p F is relevant. For the Fock Figure 5. Configurations of trajectories that contribute to the Hartree contribution to the Altshuler-Aronov correction. There is a one-to-one correspondence between the configurations for the Hartree and the Fock contribution, Fig. 3. contribution, momentum does not change at the interaction points, and the spinors at the interaction points combine to a factor Σ F (p 1 , p 2 ) = χ(p 1 )|χ(p 1 ) χ(p 2 )|χ(p 2 ) = 1.
Thus, up to the degeneracy d g , the result for the Fock contribution remains unchanged as compared to the conventional metal. The spinor structure will however influence the Hartree contribution δG H AA , as we now show. The relevant trajectory configurations for the Hartree correction are shown in Fig. 5. There is a one-to-one correspondence between the trajectory configurations for the Fock and Hartree contributions. However, unlike for the Fock contribution, the configurations for the Hartree contribution involve a finite-angle crossing at momenta p 1 and p 2 , which has two important consequences: First, it leads to an additional difference in the classical actions of the trajectories, resulting in the fast-oscillating factor e i(p 1 −p 2 )(r 1 −r 2 )/ . Such factor also enforces the interaction points r 1 and r 2 to remain close together on a scale of the Fermi wavelength. Since the function K(r 1 , r 2 ; ω) is built from classical propagators that are smooth on the scale of Fermi wavelength, we can identify r 1 = r 2 for this function. Second, the spinor structure from the interaction vertices now results in the nontrivial factor This result indeed reflects the chiral nature of the charge carriers, leading to a suppression of backward scattering processes. Combining everything, the Hartree contribution can be obtained from the Fock contribution by the replacement (54) where the factor d g comes from the existence of a closed trajectory-loop in the configurations of Fig. 5, and with the Fourier-transformed interaction For a short-range potential U (r 1 − r 2 ) ∝ δ(r 1 − r 2 ), we then find δG H AA = −(d g /2)δG F AA : The spin and valley degeneracies enhance the Hartree contribution by an extra factor d g = 4 compared to the Fock contribution, while chirality reduces it by a factor two [54,55].
Substituting the explicit expressions for the diffusion propagators, the final result for the interaction correction to the conductivity reads with the effective interaction kernel

Coulomb interaction
The Coulomb interaction, U C (r 1 , r 2 ) = e 2 /|r 1 − r 2 | is long-ranged, and screening effects need to be incorporated into the results of the previous section. The treatment of screening requires us to slightly reformulate the results of the previous Section. This discussion closely follows that of [34,46,[89][90][91], where the case of Coulomb interactions in a normal metal was considered. Since all relevant trajectories are paired in the semiclassical analysis, it is sufficient to consider interaction processes in which the transferred momentum q is small (eventually after an exchange of the particles). To be specific, we consider the scattering of two particles with initial momenta p 1 and p 2 , and final momenta p 1 + q and p 2 − q, respectively, and assign combined spin/valley indices α and γ of the incoming particles, and the indices β and δ to the two outgoing particles, respectively, see Fig.  6. Then the unscreened Coulomb interaction has the matrix elements U αβγδ (q) = U C (q)δ α,β δ γ,δ − U C (p 1 − p 2 )Σ H (p 1 , p 2 ) p 1 ,p 2 δ α,δ δ β,γ , where we neglected the small momentum q in the argument of U C for the second term and performed the average over the directions of p 1 and p 2 , consistent with the diffusive motion. It is the trace of this averaged matrix element that appears in the effective interaction kernel (57) of the perturbative analysis. Figure 6. Scattering processes of two particles with initial momenta p 1 (p 2 ) and spin+valley α (γ), and final momenta p 1 + q (p 2 − q) and spin+valley β (δ) to lowest order interaction.
The theoretical treatment of screening requires a different structure of the valley and spin indices. The interaction matrix elements are decomposed as where matrices Γ (i) act in spin and valley space, and read {s 0 ⊗τ 0 , s 0 ⊗τ k , s j ⊗τ 0 , s j ⊗τ k }, where s j (τ k ) are Pauli matrices acting in spin (valley) space, with j, k ∈ {x, y, z}, and s 0 (τ 0 ) are 2 × 2-unit matrices. Using the equality i Γ (i) αβ Γ (i) γδ = 4δ αδ δ βγ one then finds that the bare Coulomb interaction has Screening renormalizes each of the interaction channels U (i) in Eq. (57) separately. For weak interactions (gas parameter r s 1) one may use the random phase approximation, which gives where Π(q; ω) is the polarization operator Π(q; ω) = d g ν 1 + iω P (q; ω) , P (q; ω) = (Dq 2 − iω/ ) −1 being the diffusion propagator. The extension of the discussion to stronger interactions requires the use of Fermi liquid theory. Following [46,[89][90][91], the bare interaction acquires a short-range contribution set by Fermi-liquid constants, where we assume the same Fermi-Liquid constant F σ 0 in all non-singlet channels. The screened interaction is then still given by Eq. (59).
Applying this procedure to the Altshuler-Aronov correction, we find the same result as the first-order correction (56), but a different expression for the interaction kernel U R (q; ω), The singlet-channel Fermi liquid constant F ρ 0 does not enter into the correction because of the divergence of the Coulomb interaction at small momenta.
We plot the Altshuler-Aronov correction for various values of the interaction constant F σ 0 in Fig. 7. For small Ehrenfest time, one finds with c = d 2 g − 1. Such an expression is well-known from diagrammatic perturbation theory [92], where in our case, the Ehrenfest time takes over the role of the elastic scattering time as a short-time cutoff. The graphene-specific physics enters the result in two ways [54,55]: First, the constant c = d 2 g − 1 is 15 for graphene with only smooth disorder, in contrast to c = 3 for conventional metals without valley degeneracy. Second, chirality affects the interaction constant F σ 0 , as will be explained in more detail below. For small values of F σ 0 , the singlet contribution is dominant in Eq. (63), giving rise to a negative correction to the conductance. On the other hand, for graphene, for F σ 0 −0.12, the non-singlet channels render the interaction correction positive. For large Ehrenfest time, we find an exponential suppression where the prefactor of the exponential is determined by the singlet channel only to leading order in /T τ E , hence at large Ehrenfest times, δG AA is negative and has a universal behavior.
A striking consequence of this asymptotics is a sign-change of the interaction correction as a function of Ehrenfest time, provided the Fermi-liquid-type interactions in the non-singlet channels are strong enough. For graphene (c = 15) this sign change already takes place at F σ 0 −0.12, in contrast to a conventional metal (c = 3) where the sign change is observed for F σ 0 −0.45. On the other hand, the values for F σ 0 are typically somewhat smaller in graphene, as can be seen using Thomas-Fermi approximation [54,55]. For conventional metals, one has where e is the charge screened by the substrate [54,55], and κ = 2×2πνe 2 is the inverse screening length resulting from the metal electrons (a factor 2 accounts for spin). We then find where θ is the angle between the directions of momenta p 1 and p 2 . We further used k F = 2πν v F , as well as the "effective fine structure constant" α = e 2 / v F . The gas parameter r s is related to α as r s = √ 2α. For a value r s ≈ 1 we obtain F σ 0 ≈ −0.28 as a typical size for the Fermi liquid parameter.
For graphene, this calculation needs to be modified in two respects. First, the inverse screening length is twice larger, due to the valley degree of freedom. Second, chirality contributes the additional factor Σ H (p 1 , p 2 ) = cos 2 (θ/2). Both effects reduce the interaction in the non-singlet channel, so that now for r s ≈ 1 we find F σ 0 ≈ −0.1, which is close to the transition point for a sign-change as function of Ehrenfest time. The measurements of [54][55][56] report F σ 0 in a range between −0.05 and −0.15. However, we note that our theory requires graphene with a smooth disorder potential, but no other perturbations, as a necessary condition for the value c = 15 in Eq. (63), since there are 16 diffusion channels. Trigonal warping or ripples, while not invalidating the semiclassical analysis, reduce the number of diffusive channels to 8, resulting in a prefactor c = 7. In case of strong intervalley scattering, only four diffusion modes are present resulting in a prefactor c = 3 (see Ref. [93,94])although in that case the conditions for the semiclassical analysis are no longer valid. The aforementioned experiments on interaction corrections report to be in a regime where c = 3 or c = 7.

Dephasing
We now turn to the second type of interaction correction, responsible for dephasing. The way of calculating dephasing in this section follows that of Ref. [47]. An alternative discussion based on perturbation theory, similar to that of the previous Section, is given in the appendix.
Because of the interactions, the electrons are subject to a time-dependent potential V (r, t). This potential affects the phase that electrons accumulate while propagating through the sample. These phase fluctuations can be included into the classical action S α of the trajectory α as it appears in the semiclassical expression (4) for the Green function by the substitution where the correction δS α (t) depends on the time t at which the electron exits the sample. The shift reads [33,47] δS where τ α is the duration of the trajectory α. Such a shift of the classical actions does not affect the Drude conductivity, because the actions from the "retarded" and "advanced" trajectories cancel. It does, however, affect the weak antilocalization correction. Equation (38) acquires an additional factor e i(δSα(t)−δS β (t)) which, when averaged over the time t, reduces the contribution from the trajectory pair α, β by a factor The time average can be calculated using the quantum fluctuation-dissipation theorem, (Note that τ α = τ β for the trajectory pairs that contribute to weak antilocalization.) One immediately concludes that for the trajectories α and β that contribute to the weak antilocalization correction δG WAL only points r α or r β in the loop or encounter segments of Fig. 1 contribute to ((δS α (t) − δS β (t)) 2 t . To find an explicit expression for the dephasing correction in the limit of weak dephasing, we expand the correction factor (70) to lowest order in the interaction U R (q, ω) and calculate the leading interaction correction δG deph to the weak antilocalization correction δG WAL . We consider contributions from positions r α,β (t 1 ) and r α,β (t 2 ) in the loop and encounter regions separately.
The calculation for the dephasing in the loop segment is very similar to the one carried out in standard diagrammatic perturbation theory. The discussion below closely follows that of [49]. With both positions r α,β (t 1,2 ) in the loop region, see Fig. 8, we find that dephasing in the loop segment leads to the replacement P (r , r ; t) → P (r , r ; t) + δP (r , r ; t) for the loop propagator in Eq. (43), with where × P (r , r 2 , τ 3 )P (r 2 , r 1 , τ 2 )P (r 1 , r , τ 1 ).
Inserting the diffusion propagator P (q, τ ) = e −Dq 2 τ , one finds with t = τ 1 + τ 2 + τ 3 . After insertion of the interaction (62) and evaluation of the integrals over time, frequency, and momentum in Eq. (72) (see Appendix Appendix D), one obtains δP (r , r ; t) P (r , r ; t) with the dimensionless conductance g 0 = 2π d g νD, and the constant This result signifies that at large times the loop propagator gets suppressed by interactions. Here we calculated the leading order correction, describing the onset of an exponential suppression. (In fact, in two dimensions the decay is not purely exponential [49], but contains an additional logarithm e −at ln t , as can be seen from Eq. (75)). We estimate the dephasing rate as the time when the leading correction becomes unity, i.e., 1 The leading logarithmic dependence in this expression agrees with that obtained in [95] for a standard two-dimensional electron gas (d g = 2). We now turn to the encounter region, where dephasing leads to an additional suppression of weak antilocalization, if the typical time for the encounter passage, the Ehrenfest time τ E , is sufficiently long. As discussed before, dephasing is ineffective, as long as the trajectories coincide. Within the encounter region, the trajectories α and β are separated by a small distance, which does not exceed the classical correlation scale L c . Dephasing then only plays a role for interaction that transfers a momentum larger than inverse mean free path, and therefore can resolve such small distance [33,96,97]. On the other hand, for low temperatures T τ tr 1 one has ωτ tr 1. In this limit, the imaginary part of the screened interaction reads [46] Figure 8. Dephasing in the loop segment: for the calculation of the lowest-order interaction correction the propagation around the loop is split into three segments r → r 1 (duration τ 1 ), r 1 → r 2 (duration τ 2 ), and r 2 → r (duration τ 3 ). For our figure, the trajectories α (solid) and β (dashed) are travelled in clockwise and counterclockwise direction, respectively. The left diagram represents the contribution from the term proportional to e iq·rα(t1)−iq·rα(t2) in Eq. (71); the right diagram represents the contribution from the term proportional e iq·r β (t1)−iq·rα(t2) . The remaining two contributions, from terms proportional to e iq·rα(t1)−iq·r β (t2) and e iq·r β (t1)−iq·r β (t2) are not shown. The figures have been drawn for the case that 0 < t 2 < min(t 1 , t − t 1 ), which is the domain of integration in Eq. (72).
where ν = k F /2π v F is the density of states and we abbreviated Note that Im U R (q; ω) is proportional to q −1 , which is different from the dependence Im U R (q; ω) ∝ q −2 of the diffusive limit. This difference will result a different Tdependence of the dephasing rate in comparison to the loop contribution [33]. We proceed by the integration over ω, which can be done explicitly using the known ω dependence of ImU R (q; ω) [49], dω 2π Here the function w(x) = (x coth x − 1)/ sinh 2 x is peaked around x = 0, normalized ∞ −∞ dxw(x) = 1, and w(0) = 1/3. Hence, the times t 1 and t 2 need to be close together on the scale of inverse temperature. On the other hand, dephasing sets in on times much larger than T −1 , as we will show below. In the following, we therefore may assume that |t 1 − t 2 | τ E , when we consider encounters that are long enough to be affected by dephasing. (In fact, the calculation below will show that the main contribution stems from time differences |t 1 − t 2 | much smaller than the elastic mean free time.) In particular this allows us to consider the effect of interaction during the first and second passage through the encounter separately, since they are separated by a loop of long duration. The same observation also allows us to neglect contributions where t 1 is in the encounter, whereas t 2 is in the loop or vice versa. We therefore focus on the first passage through the encounter region, where the trajectories α and β are separated by a distance d(t) = r β (t) − r α (t), with the magnitude where t is varying form 0 to τ E . We can use this to rewrite the last two factors of Eq. (71) as wherer(t) = [r α (t) + r β (t)]/2 represents a trajectory intermediate between α and β. After performing the average over disorder configurations, we find that inclusion of the leading-order dephasing correction amounts to the replacement P (r , r; τ E ) → P (r , r; τ E )+2δP (r , r; τ E ) in Eq. (43), where the factor two accounts for the two passages through the encounter region, with and see Fig. 9. Because of the smallness of |t 1 − t 2 |, the propagator P (r 2 , r 1 , t 1 − t 2 ) is the ballistic propagator, whereas the propagators P (r , r 2 , τ E − t 1 ) and P (r 1 , r, t 2 ) can be taken in the diffusion approximation. We change the integration variables to the mean timet and the difference time t = t 1 − t 2 . Again using the smallness of |t 1 − t 2 |, we replace d(t 1 ) and d(t 2 ) by d(t). We neglect correlations between d(t) and the direction of the velocity at timet. Using P (r 1 , r, t 2 ) P (r 2 , r, t 1 ), again because of the smallness of |t 1 − t 2 |, we find where we inserted the Fourier transform of the ballistic propagator. Since the integration over t converges for |t| ∼ 1/v F q, the argument of the function w may be set to zero in Eq. (83). Finally, the angular average over the direction of d(t) gives a factor 1 − J 0 (qd(t)), so that we find We cut off the logarithmic divergence of the q integration at large q at λ −1 F , which gives The remaining time-integration is easily evaluated with the help of Eqs. (81), (37), and we obtain Figure 9. Dephasing in the encounter segment: For the calculation, the encounter is split into three segments r → r 1 (duration τ 1 = t 2 ), r 1 → r 2 (duration τ 2 = t 1 − t 2 ), and r 2 → r (duration τ 3 = τ E − t 2 ). Only configurations in which |t 1 − t 2 | τ E contribute to the interaction correction δG deph . In the middle segment, the distance between the trajectories α and β is d(t) = λ F e λt , wheret = (t 1 + t 2 )/2. The figure has been drawn for the case 0 < t 2 < t 1 , which is the domain of integration in Eq. (83) Hence, our the final result reads One may identify the right-hand side of Eq. (89) with τ E /τ enc ϕ , where τ enc ϕ is an effective dephasing time for the encounter region. With this identification, Eq. (89) describes the onset of an exponential suppression of the weak localization ∝ e −2τ E /τ enc ϕ at large Ehrenfest times. Note that the time τ enc ϕ is twice the dephasing time τ ball ϕ that one finds from dephasing in the loop region in the ballistic regime [95], consistent with the theory of [33]. [To compare with [95], take the energy-dependent dephasing time τ ϕ (ε) from Eqs. (18) and (19a) of Ref. [95] in the limit T τ and calculate (τ ball ϕ ) −1 = dε (−∂f /∂ε) τ ϕ (ε) −1 with, for a conventional metal, d g = 2. The lowmomentum cut-off in [95] is the inverse mean free path, whereas it is the classical correlation length L c in our case. The two lengths need not be equal, see Eq. (26).]

Conclusion
In this article we have presented a trajectory-based semiclassical theory of the quantum corrections to transport in graphene in the presence of a random potential that is smooth on the scale of the Fermi wavelength. A prominent role is played by the Ehrenfest time, which serves as a short-time threshold for the appearance of quantum interference effects. The Ehrenfest time also plays an important role for electrons in a conventional two-dimensional electron gas (with quadratic dispersion) if they are subject to a smooth random potential.
Compared to the conventional case, charge carriers in graphene have an additional pseudospin degree of freedom and they have an additional valley degeneracy, which leads to a few subtle modifications of the quantum corrections with respect to the conventional case. The pseudospin vector always points along the direction of motion, reflecting the chiral nature of the charge carriers in graphene. The evolution of the pseudospin along the trajectory is associated with a Berry phase of the spin transport, that additionally enters the semiclassical Green function. This Berry phase is responsible for a sign change in the weak localization correction, giving antilocalization behavior. The presence of a finite Ehrenfest time reduces the magnitude of this correction, but with a multiplicative factor that is the same for weak localization and weak antilocalization. We also considered the suppression of weak (anti)localization from dephasing at finite temperatures, and identified there, too, the role of the Ehrenfest time.
For the interaction correction there are two important differences with the case of the conventional two-dimensional electron gas: The Hartree-type processes (or, more precisely, interaction non-singlet channels) contain an additional angular dependence, as a result of chirality. Moreover, importance of screening is changed, because of the presence of the valley degeneracy. A finite Ehrenfest time suppresses the Altshuler-Aronov correction, in a similar way as for conventional metals, but unlike for weak (anti)localization the suppression is not simply a multiplicative factor. Interestingly, the interaction correction may undergo a sign change as a function of Ehrenfest time, for sufficiently strong interaction in the non-singlet channels. For graphene, the interaction strength at which this sign change takes place is smaller than in conventional electron gases, which may place it within experimental reach, as discussed in Sec. 5.
where V (t) is shorthand notation for V (r(t)). Upon integrating the evolution equations for an infinitesimal time interval δt the solution may be cast in the form of a transfer matrix equation, which we write as where z 2 = √ K 0 /k F ξ 1 and the transfer matrix M(t, t + δt) reads a stochastic function that contains all information on the random potential. The function f has zero mean, and its fluctuations in a time interval ∆t long in comparison to the correlation time (The condition ∆t t ξ is consistent with the smallness of the parameter z.) Sofar we have calculated the transfer matrix for an infinitesimal time interval δt. The result can be easily extended to calculate the transfer matrix for time intervals of arbitrary duration, via successive multiplication of transfer matrices valid for the infinitesimal segments. This results in a stochastic evolution of the transfer matrix, which can be analyzed using an explicit parameterization of the transfer matrix, where σ 2 and σ 3 are the Pauli matrices. The exponential divergence of the trajectories follows from the radial parameter l, For the calculation of l it is sufficient to consider the matrix product M T M, which has eigenvalues e ±2l and no longer depends on the angular variable ϕ. The time-evolution of the remaining parameters l and φ is given by a Langevin-type process which, for large l, reads where terms of higher order than δt are neglected.
It is helpful to introduce the variable y via y = (2π) 1/6 z −1/3 (2/3) 1/3 cot φ. (A.10) After averaging over fluctuations of f , we find that mean and variance of the change δy in an infinitesimal time interval δt read (A.11) The parameter y acquires a stationary probability distribution P s (y), which satisfies the equation [  Keeping leading terms in the parameter z only, we find that the average Since the angular variable y evolves statistically independent from the radial variable l for large l, we may average y with the help of the stationary distribution (A.14), for which we find From the definition (A.8) we then obtain the Lyapunov coefficient where, in the second equality, we inserted the transport mean free time and the mean free path of Eq. (26). This result agrees with the Lyapunov exponent calculated by Aleiner and Larkin [32]. (One has to identify the short-length cut-off a of Ref. [32] with ξ/ √ 3, see the text below Eq. (A3) of [32].
where τ D = L 2 /Dπ 2 is the dwell time, and r = W/L is the aspect ratio. For large aspect ratios r, the summation over k can be replaced by an integration, We are interested in leading terms in the small parameter τ tr /τ D only, for which one finds In the limit of small τ tr /τ D , the behavior of the summand for large l is relevant. We then may simplify the summation over n as follows: The main contributions arise for n ≈ 0 and n ≈ 2l, where the summand has poles. If l is large, the poles are well separated, and the dominant contribution comes from the pole at n ≈ 0, ∞ n=1, n odd which results in the expression For small l, this expression is not accurate, hence this summation has a lower cutoff, which is not relevant for small τ tr /τ D , however, where the l-summation results in ln τ D /τ tr . Hence, we find Eq. (44) from the main text.
The calculation proceeds by inserting the semiclassical expressions for the Green functions and identifying the relevant configurations of trajectories. Only configurations where "advanced" and "retarded" trajectories are paired up (where we also allow for small angle encounters or pairing of time-reversed trajectories) contribute systematically to the conductance. For the first term inside the trace in Eq. (C.1), this is only possible, if the three "retarded" trajectories join together to a single trajectory connecting the points r with r, that can be paired up with the advanced trajectory. In the semiclassical approximation, we then evaluate the integration over r 1 and r 2 within stationary phase approximation, where we keep only stationary configurations that join to a single trajectory. The result of such a calculation is dr 1 dr 2 G R (r, r 1 ; ε)G R (r 1 , r 2 ; ε − ω)G R (r 2 , r ; ε)Im U R (r 1 , r 2 ; ω) = − 1 2 2π (2πi ) 3/2 α:r →r;ε A α e iSα/ |χ(p α ) χ(p α )|e iγα × τα 0 dt t 0 dt Im U R (r α (t), r α (t ); ω) e −iω(t−t )/ , (C.2) We here restrict ourselves to explain, how the structure of this result can be understood, and refer to [34] for the detailed calculation. The first step is to identify points r 1 and r 2 , which make the total phase of the integrand stationary. Such configurations are obtained, whenever there exists a single classical trajectory α the connects the points r and r via r 2 and r 1 . Since the position of the intermediate points can be anywhere along the trajectory α, the summation over stationary configurations of the intermediate points is expressed a summation over trajectories α as well as two time integrations along the trajectory α. The Green function connecting the intermediate points is taken at a different energy, resulting in the additional factor e −iω(t−t )/ , as follows from Eq. (6). Furthermore, the actions of the three subpaths sum up to the action S α of the joined path. Similarly, the individual Berry phases for the subpaths combine to the Berry phase of the joined path γ α , as the Berry phase is expressed as integral along the trajectory. Since the momenta are smoothly connected at the intermediate points, the intermediate spinors match together, and only the spinors at the endpoints remain in the final expression. The second step in the evaluation of the integral is to integrate over quadratic fluctuations around the stationary configurations. This, in turns, renders the proper stability amplitude A α and the prefactor, see [34].
Similar considerations apply to the second and third term of the trace of Eq. (C.1). Therefore, we can write Eq. (C.1) as a sum over one "retarded" trajectory α and one "advanced" trajectory β, connecting source and drain. Since only paired trajectories are of relevance, we may simplify A α = A β and τ α = τ β . Since the only dependence on the propagation energy ε is in the factor between square brackets on the first line of Eq. (C.1), we may perform the integration over ε and find Inserting the Fourier transformed interaction we then obtain δG deph = e 2 d g (2π ) 2 dy dy

(D.6)
Here we introduced the dimensionless conductance g 0 = 2π d g νD. Our expression now reads where F(x) = x/ sinh 2 x. For small ω/T , we may expand the function F(x) ≈ 1/x, which gives a logarithmic divergent ω-integral. This integral should be cut at high energies by temperature, and at small energies by /t, where Eq. (D.5) ceases to be valid. So we find The result for the singlet channel can be obtained from the latter equation by formally sending F σ 0 → ∞.