{'fraction_non_alphanumeric': 0.043304497394647645, 'fraction_numerical': 0.01492386300399656, 'mean_word_length': 4.326866073834546, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Relaxation of charge in monolayer graphene: fast non-linear diffusion vs Coulomb effects

Pristine monolayer graphene exhibits very poor screening because the density of states vanishes at the Dirac point. As a result, charge relaxation is controlled by the effects of zero-point motion (rather than by the Coulomb interaction) over a wide range of parameters. Combined with the fact that graphene possesses finite intrinsic conductivity, this leads to a regime of relaxation described by a non-linear diffusion equation with a diffusion coefficient that diverges at zero charge density. Some consequences of this fast diffusion are self-similar superdiffusive regimes of relaxation, the development of a charge depleted region at the interface between electron- and hole-rich regions, and finite extinction times for periodic charge profiles.


I. INTRODUCTION
Understanding how excess charge carriers approach equilibrium is important to both condensed matter and device physics. In an isotropic three-dimensional medium of electrical conductivity σ and dielectric constant κ the driving force of relaxation is the Coulomb interaction; relaxation is exponentially fast and characterized by the Maxwell relaxation rate 4πσ/κ. In earlier pioneering studies Dyakonov and Furman, and Govorov and Chaplik (DFGC) 1 have recognized that the Maxwell relaxation law does not apply to low-dimensional systems. In the two-dimensional case relevant to semiconductor heterostructures 2 , the two-dimensionally confined charges interact according to the three-dimensional Coulomb law; then the counterpart of the Maxwell relaxation rate (2πσ/κ) is dimensionally a velocity, suggesting that the charged region expands with constant velocity proportional to σ/κ 1 . Non-Maxwellian relaxation in twodimensional systems has been experimentally detected 3 .
Here we will discuss the relaxation of charge in monolayer graphene 4 . This problem was posed by Efros 5 , who conjectured that Coulomb forces continue to dominate the physics, so that the DFGC theory 1 would apply. However, we will find that in the low doping limit the situation is somewhat delicate. Graphene is distinguished from other two-dimensional electron systems in that reduction in doping causes a decrease in the density of states (DOS). This has its origin in the pseudo-relativistic dispersion law of graphene's elementary excitations 4 , and translates 2 into a decreasing screening response as the Dirac point is approached. This can be quantified as an increase in the Debye screening length with decrease in doping. Screening is an important consideration here because the Coulomb interaction dominates relaxation only on the scale exceeding the screening length. When the screening length is large compared to other characteristic scales of the problem, quantum-mechanical effects can take over. This situation is most pronounced at the Dirac point, where the DOS vanishes and the screening length is infinite. Here in a relevant range of parameters we identify a novel regime of relaxation described by a non-linear diffusion equation of the type first obtained in theories of thermal waves and filtration 6 and later encountered in studies of dynamics of crystal interfaces below the roughening phase transition 7 , soliton dynamics in one-dimensional conductors 8 , diffusion in a plasma 9 , spreading of liquid drops in the presence of Van der Waals interactions 10 , and self-organized criticality 11 . Relaxation in graphene has more in common with the last three examples as they all are characterized by having a divergent diffusion coefficient in the limit that the density of the diffusing quantity vanishes. In graphene the effect of the "fast" non-linear diffusion manifests itself in the self-similar superdiffusive regimes of relaxation, the development of charge depleted regions at the interface between electron-rich and hole-rich regions, and the finite extinction times of periodic charge profiles. Recently demonstrated high-resolution non-invasive imaging of charge currents in graphene structures 12 can be employed to observe these effects.
Below we focus on the interesting case of monolayer graphene at the Dirac point. Our continuum theory is limited to zero temperature and applicable in the longwavelength, low-frequency limit. We additionally neglect the effects of retardation as all the relevant velocity scales that we will encounter are significantly smaller than the speed of light. Whenever it cannot cause confusion, we also suppress the time dependence of the dynamical variables.

II. STATEMENT OF THE PROBLEM
The total potential ϕ(r) felt by a carrier of charge e at a position r is due to the potential of the other carriers of number density n(r) where κ is due to graphene's own electrons and surrounding environment. In thermodynamic equilibrium the electrochemical potential is zero across the system. Here ζ(n) = v F √ πn is the chemical potential for charge carriers in the absence of a perturbing potential ϕ(r), v F is the Fermi velocity, and the density n is assumed to be small enough that the low-energy Dirac dispersion law, ε = v F k, holds 4 . The equilibrium state is carrier-free, n = 0, and ϕ = 0. Out of equilibrium the driving force for evolution is the negative of the gradient of the electrochemical potential: where the first term, inversely proportional to the DOS ∂n/∂ζ, is the consequence of Fermi statistics (increasing density causes an increase in the total kinetic energy of the charge carriers), while the second term is the Coulomb force. The DOS determines the Debye screening length of the two-dimensional electron gas 2 which diverges as n → 0. Here α is the fine structure constant for graphene 4 ; throughout this work we assume that α ≪ 1 or equivalently κ ≫ 1. Suspended graphene, characterized by α ≈ 1/2 13 , is at the verge of applicability of our theory. According to Ohm's law the force (3) induces an electric current density where the conductivity of graphene σ 0 is approximated by its intrinsic value σ(n → 0) = 4e 2 /π 14 . Conservation of charge within a region implies a continuity equation

III. RELAXATION DOMINATED BY COULOMB EFFECTS
The quantum-mechanical effects accumulated in the first term of the expression for the driving force of relaxation (3) can be neglected when the condition holds. This requires that the characteristic length scale of the problem L be significantly larger than the local Debye screening length q −1 s (n) (4). Then the theory is linear and can be solved by a Fourier transform; it is equivalent to the DFGC theory 1 , whose hallmark is the linear evolution equation where the velocity scale v 0 owes its existence to the intrinsic conductivity of graphene 14 . We now quote two solutions 1 of Eq.(8) supplemented with their ranges of applicability according to Eq. (7). For an initially localized one-dimensional distribution of number density N = ∞ −∞ n(x, 0)dx (a stripe of excess charge), the largetime asymptotic solution to Eq.(8) has the form while for an initially localized two-dimensional distribution of strength Z = n(r, 0)d 2 r, one similarly finds These results show that most of the charge in the expanding clouds is localized within a distance L ≈ v 0 t from the origin; the clouds spread with the constant velocity v 0 . Eq.(9) applies at times exceeding the cross-over time scale t x which can be very large for N small. In this case the intermediate asymptotic behavior for t ≪ t x is dominated by quantum-mechanical effects. Similarly, Eq.(10) fails to describe relaxation of clouds of net dimensionless charge Z smaller than the crossover charge Z x . For example, for α = 1/10 one finds Z x ≃ 100. Excess charges smaller than a 100 are certainly easier to realize in practice than charges in excess of it. Existence of a wide range of parameters where quantum-mechanical effects dominate relaxation of charge in graphene is a direct consequence of the vanishing DOS, ∂n/∂ζ ∝ √ n → 0 as n → 0.

IV. RELAXATION DOMINATED BY QUANTUM EFFECTS
Having established the range of applicability of the DFGC theory 1 to relaxation in graphene, we now turn to the regime when physics is dominated by quantummechanical effects, and assume that the condition opposite to Eq.(7) holds.
Neglecting the Coulomb term in Eq.(3), the force exerted on a carrier causes an electric current (5) where we have introduced the diffusion coefficient The first representation is the Einstein relation, which makes it clear that the n −1/2 divergence of D(n) is the consequence of the finite intrinsic conductivity of graphene 14 and of the vanishing DOS at the Dirac point 4 . Curiously, Planck's constant which is present in both the conductivity (σ 0 = 4e 2 /π 14 ) and in the inverse DOS, ∂ζ/∂n ∝ , cancels out. Thus the outcomes (11) and (12) are independent of .
A. Steady state in the presence of fixed current The simplest non-equilibrium effect to consider is the steady state in the presence of a fixed current j flowing through the system. This can be due to the hole flow along j, the electron flow opposite to j or both flows present at once and meeting at a line annihilation front that is perpendicular to j. Focusing on the last case, choosing the direction of j to coincide with the negative x axis the charge density distribution ρ(x) = en(x) can then be found by integrating Eq. (11): (13) where the holes (e > 0) present in the x > 0 region flow toward x = 0 where they annihilate with the electrons present at x < 0 and flowing toward x = 0. Since graphene is a gapless material, there is no energy released during the event of the electon-hole annihilation, and thus there is no additional force associated with the process. The x 2 behavior of the density is a necessary outcome of having a finite current density at the annihilation line x = 0, made possible by the interplay between the diffusion coefficient (12) (which diverges as 1/x) and the gradient of the density in Eq.(11) (which vanishes as x). The rapid x 2 drop in density near the annihilation front is an example of a charge depletion region or a sink brought about by fast singular diffusion. The result (13) is valid within a range |x| L satisfying the condition αL 2 j/v F e ≪ 1.

B. Fast diffusion equation
Substituting the expression for the current density (11) into the continuity equation (6) and employing Eq. (12) we arrive at our central result This is a non-linear diffusion equation with a diffusion coefficient that diverges as n → 0. Not only does it not contain Planck's constant, but it also does not depend on the electron charge or the dielectric constant κ, because the equation describes the regime where the Coulomb interaction is negligible. The only material parameter entering Eq. (14) is the Fermi velocity v F . The n −1/2 singularity of the diffusion coefficient is the same as that found in the problem of diffusion of plasma particles in the presence of magnetic field of a toroidal multipole 9 . The remarkable property of the non-linear diffusion equation with a power-law dependence D(n) is the existence (for a narrow class of problems) of exact selfsimilar solutions capturing the asymptotic large-time behavior of more general problems 6 . The cases of "slow", D(n → 0) → 0, and "fast", D(n → 0) → ∞, diffusion are qualitatively different, and graphene belongs to the latter category. Below we analyze a few particular cases of late stage relaxation in graphene which illustrate the salient features of fast diffusion. These can be understood in simple terms via a combination of scaling arguments and conservation laws 6 , and are likely to be realized in practice.

One-dimensional localized charge distribution
When the initial and/or boundary conditions depend only on one coordinate, Eq. (14) becomes An initially localized distribution (such as ∞ −∞ n(x, 0)dx = N ) remains localized through the evolution.
Then dimensional considerations determine L, the width of the region where most of the charge is localized at a time t, to have the form Combining this with the estimate nL ≃ N following from conservation of linear charge density (N = ∞ −∞ n(x, t)dx), we find where the displayed numerical factors hereafter correspond to the asymptotically exact propagation speed of the charge wave dL/dt inferred from the exact solutions given below. We see that the charged region expands with time faster than t 1/2 (linear diffusion) because the diffusion coefficient (12) grows in the direction of spreading. It is straightforward to verify that at the crossover time scale t ≃ t x , Eq.(9), the size of the cloud (16) has the same order of magnitude as that supplied by the DFGC theory, v 0 t x . In the bulk of the distribution (|x| ≪ L(t)), the density decreases with time as n(0, t) ≃ N/L(t) ≃ N 4/3 (v F t) −2/3 ; well outside it (|x| ≫ L(t)), the density can be directly estimated from Eq.(15) as n/t ≃ v F √ n/x 2 or n(x, t) ≃ (v F t) 2 /x 4 . The combination of these observations implies that the solution to Eq.(15) has the self-similar form where the scaling function f (ξ) satisfies the condition of conservation of charge density ( ∞ −∞ f (ξ)dξ = 1) and behaves as f ′ (ξ → 0) = 0 and f (ξ → ∞) ≃ 1/ξ 4 . The exact result f (ξ) = 2/π(ξ 2 + 1) 2 due to Zel'dovich, Kompaneets, and Barenblatt 6 exhibits these features.

Two-dimensional localized charge distribution
For an initially localized two-dimensional distribution (Z = n(r, 0)d 2 r), similar arguments hold except that we seek a solution to the radially-symmetric version of Eq.(14) and conservation of charge (Z = n(r, t)d 2 r) now implies the estimate nL 2 ≃ Z. As a result, for the size of the charge cloud at large time we find We see that the charged region spreads out with the constant velocity 8v F / √ Z significantly exceeding that of the DFGC theory, v 0 = 8αv F , Eq.(8); both velocities have the same order of magnitude at the cross-over charge Z x . Moreover the expansion is faster than in one dimension (16) because the diffusion coefficient grows in all directions away from the cloud center. Ballistic spreading with a velocity proportional to v F can be anticipated from the outset because Z is dimensionless while v F is the only parameter of the problem having dimensionality (of velocity). For a localized two-dimensional distribution the solution to Eq.(18) has the self-similar form where the scaling function g(ξ) satisfies the condition of conservation of charge (2π ∞ 0 g(ξ)ξdξ = 1) and behaves as g ′ (ξ → 0) = 0 and g(ξ → ∞) ≃ 1/ξ 4 . The exact solution g(ξ) = 1/π(ξ 2 + 1) 2 , also due to Zel'dovich, Kompaneets, and Barenblatt 6 , exhibits these features.

One-dimensional localized dipole charge distribution
Relaxation of one-dimensional charge profiles can be further accelerated compared to the t 2/3 result (16) by introducing a sink into the problem. Indeed, let us consider a charge distribution composed of electrons (x < 0) and holes (x > 0) with localized number density, so that the total charge is zero while the first (dipole) moment ∞ 0 xn(x, 0)dx = P is finite. Due to annihilation events at x = 0, the total number of electrons or holes is no longer conserved; however, the dipole moment P = ∞ 0 xn(x, t)dx is conserved by the equation of motion (15) 6 . This conservation law implies the estimate P ≃ nL 2 . When combined with the relationship L 2 ≃ v F t/ √ n, we find We see that similar to the case of a localized twodimensional distribution (Eq. (19)), the region occupied by carriers of each sign expands ballistically with a velocity of the order v F / √ P . The total number of carriers (per unit length) of a given kind decreases with time as nL ≃ P 3 /v F t. The analytic solution in this case has a self-similar form resembling Eq.(20), where the scaling function h(ξ 0) satisfies the condition of conservation of dipole moment( ∞ 0 ξh(ξ)dξ = 1) and behaves as h(ξ → ∞) ≃ 1/ξ 4 and h(ξ → 0) ≃ ξ 2 . The latter property follows from the observation that the sink at x = 0 only differs from its steady-state counterpart (13) by time-dependent current j(0, t). The exact solution in this case, h(ξ) = 9 √ 3ξ 2 /2π(1+ξ 3 ) 2 (a relative of the "dipole" solution due to Barenblatt and Zel'dovich for the slow diffusion case 6 ) exhibits these properties.

One-dimensional periodic charge distribution
One-dimensional relaxation can be drastically accelerated in the presence of a series of sinks. A relevant example is a periodic charge profile with both electrons and holes present in equal amounts. The spatial periodicity of the distribution is preserved by time evolution which means that only a half-period occupied by carriers of one sign needs to be considered. We further assume that the evolving density profile is symmetric about its maximum at x = L, where the current density (11) vanishes. It is then sufficient to solve Eq.(15) in the 0 x L quarter-period range subject to the boundary conditions n(0, t) = 0, ∂n(x, t)/∂x| x=L = 0, and then employ symmetry to extend the results beyond it. If the initial density profile has amplitude n 0 , position is measured in units of L and time in units of the characteristic diffusion time √ πn 0 L 2 /2v F , we seek a solution to the problem in the form n(x, t) = n 0 ν(x/L, 2v F t/L 2 √ πn 0 ) where the dimensionless density ν(ξ, τ ) evolves according tȯ Assuming a separable form for the solution (ν(ξ, τ ) = a(τ )s(ξ)), the amplitude a(τ ) and shape s(ξ) functions satisfy the equationṡ a a 1/2 = 2 where λ 0 is a separation constant; without loss of generality we can set s(1) = 1 9 . The solution for the amplitude is then a(τ τ ext ) = λ 2 (τ ext − τ ) 2 /4, and a(τ > τ ext ) = 0, thus implying that the charge profile goes extinct in a finite time τ ext ≃ 1. This corresponds to the diffusion time scale √ n 0 L 2 /v F in physical units, and represents another manifestation of fast diffusion.
The equation for u = s 1/2 can be integrated; it is similar to the Newtonian motion of a particle of unit mass that has position u at time ξ in a potential V = λu 3 /6: where B(x, y) is Euler's beta function, and the value of λ ≈ 5.898 is fixed by evaluation of the integral at ξ = 1. Collecting these observations, the solution to the problem (23) is ν(ξ, τ ) ≈ 8.698(τ ext − τ ) 2 s(ξ) for τ τ ext , and ν = 0 for τ > τ ext . It has been proven that an arbitrary initial profile ν(ξ, 0) on an interval satisfying Dirichlet boundary conditions evolves into the separable form (26) and rigorous bounds on τ ext in terms of ν(ξ, 0) can be given 9 . Its implication for relaxation of periodic charge profiles in graphene is that after an initial transient they all "settle" into the fixed shape s(ξ), Eq.(25), with amplitude falling off as (τ ext − τ ) 2 ; relaxation ends in a finite time. It is straightforward to infer from Eq.(25) that near the ξ = 0 sink the profile exhibits the already familiar ξ 2 density drop.

V. INTERPLAY OF COULOMB AND QUANTUM EFFECTS AT FINITE DOPING
Quantum-mechanical effects have their strongest influence on relaxation at zero doping, but they continue to play a role at finite doping n 0 . When the deviation δn(x, t) = n(x, t) − n 0 is small, a linear theory captures the physics of relaxation. It is then straightforward to show that small-wavevector q ≪ q s (n 0 ) modes of δn evolve according to the DFGC equation (8) while largewavevector q ≫ q s (n 0 ) perturbations relax according to the linear diffusion equation with the diffusion constant D(n 0 ), Eq. (12). The divergence of the screening length q −1 s (n 0 ) (4) as n 0 → 0 is a clear indicator that quantummechanical effects control relaxation of a wide range of modes.

VI. ACKNOWLEDGEMENT
We are grateful to Eva Andrei for a valuable discussion of the topics covered in this work.