# Interaction-disorder-driven characteristic momentum in graphene,

approach of multi-body distribution functions

###### Abstract

Multi-point probability measures along with the dielectric function of Dirac Fermions in mono-layer graphene containing particle-particle and white-noise (out-plane) disorder interactions on an equal footing in the Thomas-Fermi-Dirac approximation is investigated. By calculating the one-body carrier density probability measure of the graphene sheet, we show that the density fluctuation () is related to the disorder strength (), the interaction parameter () and the average density () via the relation for which leads to strong density inhomogeneities, i.e. electron-hole puddles (EHPs), in agreement with the previous works. The general equation governing the two-body distribution probability is obtained and analyzed. We present the analytical solution for some limits which is used for calculating density-density response function. We show that the resulting function shows power-law behaviors in terms of with fractional exponents which are reported. The disorder-averaged polarization operator is shown to be a decreasing function of momentum like ordinary 2D parabolic band systems. It is seen that a disorder-driven momentum emerges in the system which controls the behaviors of the screened potential. We show that in small densities an instability occurs in which imaginary part of the dielectric function becomes negative and the screened potential changes sign. Corresponding to this instability, some oscillations in charge density along with a screening-anti-screening transition are observed. These effects become dominant in very low densities, strong disorders and strong interactions, the state in which EHPs appear. The total charge probability measure is another quantity which has been investigated in this paper. The resulting equation is analytically solved for large carrier densities, which admits the calculation of arbitrary-point correlation function.

###### pacs:

72., 77.22.-d, 68., 64.## I Introduction

Since the discovery of graphene as a hybrid between metal and insulator, this material with relativistic energy spectrum of zero-gap Dirac fermions Wallace (1947); Novoselov et al. (2004) attracts attention due to its unique electronic properties and prospective applications in nanoelectronics Beenakker (2008); Neto et al. (2009); Sarma et al. (2011); Kotov et al. (2012). The understanding of the origin and the influence of disorder, as well as interactions in graphene seems to be essential in understanding of the experiments and also in designing graphene-based electronic devices. The effect of these quantities and the resulting screening are not as easygoing as ordinary simple metals Sarma et al. (2011). There is a huge literature concerning the interplay of particle interactions and disorder in mono and multilayer graphene Stauber et al. (2005); Ostrovsky et al. (2006); Aleiner and Efetov (2006); Fogler et al. (2008); Pereira et al. (2008); Mucciolo and Lewenkopf (2010); Sadhukhan et al. (2017); Neto et al. (2009) which cause variety of phenomenons. Examples of the effect of disorder on the properties of graphene are its effect on: the compressibility Asgari et al. (2008), the electron-phonon interaction Song et al. (2012), the Fano factor and conductivity Lewenkopf2008numerica; Sarma et al. (2011), the magnetoresistance of bilayer graphene Hatami et al. (2011), the graphene Hall bars Petrović and Peeters (2016), the optical properties of graphene quantum dots Altıntaş et al. (2017) and the effective electronic mass in bilayer graphene Li et al. (2016). Disorder-based spintronics in graphene Rocha et al. (2010), and Levi-flight transport caused by anisotropically distributed on-site impurities Gattenlöhner et al. (2016) are other examples for which the disorder plays a dominant role. The disorder-mediated Kondo effect Miranda et al. (2014) and the formation of electron-hole puddles (EHPs) are examples for which the disorder and the inter-particle interactions play vital roles simultaneously. EHPs are believed to be responsible for the observed minimum conductivity of graphene and was predicted theoretically by Hwang et al Hwang et al. (2007a) and Adam et al Adam et al. (2007) and was also confirmed in experiments in the vicinity of the Dirac point Martin et al. (2008); Rutter et al. (2007); Brar et al. (2007); Zhang et al. (2009); Deshpande et al. (2009a, b); Ishigami et al. (2007); Cho and Fuhrer (2008); Berezovsky et al. (2010); Berezovsky and Westervelt (2010). They are characterized by the state in which some strong carrier density inhomogeneities with density fluctuations much larger than the average density (for low densities) emerge Rossi and Sarma (2008). In this case (for which the transport is governed by the complex network of small random puddles with semi-metal character Rossi and Sarma (2008)) and the similar situations, a many body treatment is needed in which the effect of interaction and the disorder enter to the play with the same footing. Such an investigation can be found in Herbut et al. (2008) in which it has been shown that for the gauge field randomness, related to the formation of ripples characterized by a disorder strength , displaces the IR fixed point to Herbut et al. (2008); Stauber et al. (2005). This shows the vital role of disorder in electronic properties of graphene which even can change the role of interaction and cause the system to exhibit the non-Frmi liquid behaviors Kotov et al. (2012).

The polarization function of the graphene as an important quantity which is sensitive to both the inter-particle interaction and the disorder has been obtained and analyzed in the cone approximation in the bubble expansion Barlas et al. (2007); Hwang and Sarma (2007). For the undoped graphene in one loop approximation of the polarization function (for which the intra-band excitations are forbidden due to the Pauli principle), it is shown that the static dielectric function is a constant González et al. (1999); Kotov et al. (2012) and no collective modes are allowed within the PRA approximation. For finite however, the more general solution of Shung Shung (1986) is applicable which predicts that the potential is screened with the screening length ( Fermi velocity, Fermi momentum, background dielectric constant). This predicts the Friedel oscillations in graphene to be of the form Shung (1986). Based on these results, in low-densities (in intrinsic or weakly disordered graphene) due to the lack of screening one should take into account the long-range part of the potential which leads to the issue of logarithmic phase shifts for Coulomb scatterers Kotov et al. (2012). These results were obtained for disorder-free system and within the RPA approach which is questionable for the graphene as pointed out by Mishchenko Mishchenko (2007). Despite of this huge literature, there is a limited understanding concerning the properties of the polarization function in the general form and also the behavior of multi-body charge density probability measures, especially in the low-density case in which EHPs appear.

A more reliable approach should involve both interaction and disorder simultaneously which is the aim of the present paper. For this, the Thomas-Fermi-Dirac (TFD) is a proper candidate which has the capacity to bring the effect of all LDA energies in the calculations as well as the tunable disorder in the same time. We employ this approach involving exchange-correlation energies and white noise out-plane disorder. For the EHPs, TFD techniques have proved to be very useful and are widely used Adam et al. (2007); Hwang et al. (2007a, b); Rossi and Sarma (2008). To obtain disorder corrections to the dielectric function we develop a diagramatic technique for the TFD theory. We obtain the two-body correlation function using some stochastic analysis and find the equation governing it. We show that, in addition to the screening length found by others, there is a characteristic length scale for the screening, namely which is related to the disorder strength and inter-particle interactions via in which is the disorder strength, is the dimensionless interaction factor to be defined in the text and is the substrate distance. For small densities we show that there occurs an instability in which the imaginary part of the dielectric function becomes negative and some charge density oscillations with disorder-dependent wave length emerge. At these densities screening-anti-screening transition occurs along with changing sign of the real screened potential. The disorder significantly changes the Fourier component of the screened potential for large wave lengths showing the fact that the large-scale behaviors of the system is mostly affected by the disorder. We characterize the system in the vicinity of this transition in detail. In addition to the one- and the two-body (static) correlation functions, the equation governing the total probability measure of the system is introduced and analyzed and a closed expression is proposed for it. The proposed expression is valid for large densities.

The paper has been organized as follows: In the next section we present the general construction of the problem and introduce the TFD theory for determining the ground state of graphene. The one-body probability measure is obtained in SEC. III. Section IV is devoted to finding the two-body probability measure and the dielectric function. The screened potential as well as the charge density oscillations are argued in this section. The total probability measure of the system is analyzed in SEC. V. In the conclusion section we close the paper by highlighting the main findings of the paper. Appendix1 contains the essential rules for calculating Feynman diagrams for the TFD theory. Appendix2 and Appendix3 help to facilitate some calculations over the paper.

## Ii General Construction; Ground state of Graphene

The effect of interaction and disorder in graphene has a long story in the literature. It is known that the random-phase approximation (RPA) fails to describe the interaction effects in graphene which was firstly pointed out by Mishchenko Mishchenko (2007). The renormalization group analysis in the weak coupling limit in the first order approximation shows that the Coulomb interactions are marginally irrelevant due to the logarithmically divergent velocity Kotov et al. (2012). In the strong interaction side however, it is shown that the graphene has non-Fermi liquid behaviors with power-law quasiparticle dispersion Son (2007); Kotov et al. (2012), the limit which cannot be reached for finite densities away from the Dirac point. In the weak coupling side, it is known that in contrast to ordinary metals in which the screening makes the interactions short-ranged, in low-density (intrinsic or weakly disordered) graphene due to the lack of screening one should take into account the long-range part of the potential which leads to the issue of logarithmic phase shifts for Coulomb scatterers Kotov et al. (2012). The effect of disorder is however of special importance since, as stated in the introduction, it may completely change the behavior of the Dirac Fermions Herbut et al. (2008); Stauber et al. (2005). TFD theory, as an approach which has the potential to bring the interaction and disorder in the calculations in the same footing is introduced and analyzed in this section. We mention some points concerning the TFD theory as a coarse grained method and the perturbation diagrams which are used in our analysis. The multi-point probability measures and their relations to the many body disordered response functions are of special importance in this paper which is described in the next sections.

Let us start by analyzing the zero-temperature polarization operator which is defined by in which is the time ordering operator, is the space-time, is the ground state density operator, and are the electron annihilation and creation operators and is the ground state of the system. Defining the density fluctuation operator as (setting ) one finds:

(1) |

in which the inner is the quantum expectation value, whereas the outer one stands for the disorder averaging. Note that although due to the presence of disorder this function (and other correlation functions) is not translational invariant, its disorder-averaged form is invariant. The first term is the disconnected part which has been shown in the first term in Fig. 1, whereas the second part is connected term which has extensively been analyzed in the literature Sarma et al. (2011); Kim and Neto (2008). In this figure the single-line circles show the carrier density at the point , i.e. and shows the averaging over the disorder. In the same diagrams in Appendix1 the full-circles (double line circles) show the full density which is obtained self-consistently according to the employed theory, which is the TFD equations here. In that disorder-free systems the contribution of the first part is apparently trivial, since it generates an additional irrelevant constant. For the disordered systems however its effect is non-trivial as is extensively analyzed throughout this paper. In this paper we concentrate mainly on the contribution of this term and its effect on the response functions. As an example consider the ground state expectation value of the potential energy of the system. Let us define which is the ground state density for a particular configuration of disorder and is not time-dependent (since has been chosen to be the exact ground state of the system), i.e. . Also let us define the second part as . The non-trivial effect of this term on this expectation value can be seen from the following relation Fetter and Walecka (2012) (the full treatment is postponed to Appendix1):

(2) |

from which we see that the change of energy due to this term is in which and are Furrier components of bare coulomb potential and respectively.

To investigate the effect of this term we turn to the TFD theory which contains coarse-graining of the space. In this case the contribution of the connected (second) term of the Fig. 1 becomes local, i.e. all the contributions are localized in a region in the close vicinity of the original spatial point, namely . Therefore the exchange-correlation term becomes a local term in the Hamiltonian, and the only non-local terms are the Hartree and disorder terms. Appendix1 contains some diagrammatic analysis of the problem and the screened coulomb interaction as well as the spatial charge screening of the coulomb impurities have been analyzed. The Feynman diagrams need an especial care in this case since in the diagrammatic expansion of the energy (TFD energy), the only non-local terms are the Hartree and the external disorder interaction terms, as shown in Fig. 5. The coulomb potential is screened only via the mediator as depicted in Fig. 6 in which the full (double) lines involve simultaneously the disorder and coulomb interaction lines. The overall result is that considering such a mediator results in a change in the dielectric function in which is proportional to the Fourier component of , is proportional to (wave vector) and is a characteristic wave vector which has been defined in this appendix. All of these quantities will be defined and analyzed in the following sections.

The determination of need some information about the two-body distribution function. The multi-body distribution functions play an important role in the condensed matter physics, since the transport parameters are expressed in terms of these functions via the Kubo formulas. The other example in which we need such functions is the mean field theories in which the quantum degrees of freedom are reduced (the quantum fluctuations are killed) and the disorder averaging becomes the most important and challenging problem. The behavior of multi-point functions depend on the type and strength of the disorder.

For the more general case in which we are interested in calculating multi-point correlation functions with two insertion points (r and ), the extra contribution comes from the quantities like (in which is an arbitrary function). In such cases we define the pair (or two-body) distribution measure to write the mentioned correlation function as:

(3) |

for which is a special case. is defined as the probability density of charge density to be at r and at . It is notable that carries simultaneously the effect of the quantum system (Hamiltonian of the system) as well as the disorder statistics. Therefore finding the multi-point correlation functions reduces to finding multi-body probability densities. In the following we describe the TFD theory of the graphene monolayer to be used in the next sections as a quantum mechanical model to calculate these multi-body probability measures.

In graphene the carrier density is controlled by the gate voltage in which is the substrate dielectric constant and is its thickness and is the gate voltage. The experimental data show a strong dependence on in which is the impurity density. In ordinary densities, the conductivity is linear function of and for very low ’s, it reaches a minimum of order which is linked with the formation of EHP’s. Using local density approximation one can prove that the total energy of the graphene for a disorder configuration and a density profile is Sarma et al. (2011):

(4) |

in which is the Fermi velocity, is the dimensionless interaction coupling constant, is the chemical potential, is the total spin and valley degeneracy. The exchange-correlation potential is calculated to be Polini et al. (2008):

(5) |

in which is the momentum cut-off and . The remote Coulomb disorder potential is calculated by the relation:

(6) |

in which is the charged impurity density and is the distance between substrate and the graphene sheet. For the graphene on the SiO substrate, , so that , nm, where is the graphene lattice constant nm corresponding to energy cut-off eV. One can easily minimize the energy with respect to and obtain the following equation:

(7) |

which should be solved self-consistently. In this paper we consider the disorder to be white noise with Gaussian distribution and .

Let us now concentrate on the scaling properties of this equation excluding . By zooming out of the system, i.e. the transformation , we see that the equation remains unchanged if we transform as expected from the spatial dimension of . This is because of the fact that . This symmetry is very important, since it causes the system to be self-affine. This scale-invariance in two dimension may lead to power-law behaviors and some exponents. It may also lead to conformal invariance of the system, and if independent of type of disorder, brings the graphene surface into a class in the minimal conformal series. The existence of makes things difficult, since in which . Therefore the rescaled equation is:

(8) |

in which . Therefore the first term survive marginally in the infra-red limit and the scale invariance is expected, even in the presence of .

## Iii One-Body probability density

In this section we concentrate on calculating one-body probability density by focusing on Eq. 7. From now on stands for the disorder averaging. Using the identities and we reach to the following identity:

(9) |

which is essential in the following sections. In this equation and and . An important point is superficial contradiction that we have used which is ignored in the TFD theory, since it is valid only when ( is the Fermi wave number). We turn back to this point at the end of this section.

The differential form of the charge profile is:

(10) |

Now let us perform some Ito calculations to obtain the probability measure of . We consider a one-body function and let be the probability measure of charge density to be at the position r, so that . The differential of is defined as , according to which ():

(11) |

Let us consider the above expression term by term. Firstly we note that which arises from the fact that . The remaining two terms, can easily be transformed to the following form by using the integration by parts:

(12) |

If we note that this equation should hold for all arbitrary s, we can demand that the equality holds for the integral kernels and obtain:

(13) |

Now let us concentrate on and . The calculation of these quantities is done to appendix VI from which we see that , and also in which . If we substitute the above relations in the Eq. 12, i.e. replacing and with their averages, we obtain:

(14) |

Therefore for the homogeneous system, which is independent of the observation point r, the left hand side of the above equation equals to zero which results to:

(15) |

By equating the expression inside the bracket to a constant (which is zero for the symmetry considerations) we obtain the equation governing the distribution of as follows:

(16) |

in which and . It is worth noting that in obtaining the above differential equation we have divided both sides of the equation by a task which is true only for non-zero finite s. In fact is everywhere non-zero and well-defined and finite except at and for which becomes zero and negative infinity respectively. Therefore the equation 16 is everywhere reliable except at for which, as we will see, becomes divergent and for which the solution is not reliable. This divergence at the Dirac point is due to the non-analytical behavior of in Eq. 9, i.e. the singular behavior of carrier density at the Dirac point. Such a behavior is missing in the other ordinary two-dimensional electronic systems and is specific to the Dirac electrons systems, since their kinetic energy densities are proportional to whose derivative is not well-behaved in .

With the above limitation in mind, let us first turn to the asymptotic behaviors of Eq. 16. One may be interested in the answer for weak coupling limit , or the weak disorder limit , i.e. large limit. In this limit, and considering to be nearly constant, we have in which . The solution is therefore ():

(17) |

in which is the area of the sample and is the normalization constant. This relation may seem to be unsuited, since it grows unboundedly for negative values of . Actually there is no contradiction, due to the presence of whose amount grows negatively for negative values, which returns the above equation to the expected form. In fact the original charge equation has electron-hole symmetry for the case which should result to an electron-hole symmetric form of . Our approximation (considering as a constant) violated this symmetry. Re-considering this quantity as a dynamical variable retains the mentioned symmetry. It is also notable that the second term in the exponent () has been inserted due to some symmetry considerations and the above equation satisfies the original equation of .

In the above equation, the effects of disorder and interaction and have been actually coded in . Large amounts of results to a very large charge fluctuations, showing that the interaction and the disorder favor small deviations from mean density, whereas favors large charge fluctuations. To see how controls the charge fluctuations, we note that the limit (zero-gated graphene) has a direct effect on . In fact controls which directly affects which is . In the limit , we expect that becomes vanishingly small, so that implying large scale density fluctuations, which in turn leads to EHPs.

The total solution of Eq 16 is , in which :

(18) |

In the above equation is the solution 17. As stated above, this solution is reliable for intermediate values of carrier density. To see the limitations of this solution, let us mention that the TFD solutions are valid only when . Approximating , we find that . In Fig. 1(a) we have shown a plot which compares and for , and . We see that the validity of the solution 18 is limited. Therefore for very high densities we do not expect that the one-body probability density be of the form 18. By the similar argument, we see that this solution is not the exact solution for limit. In the Fig. 1(b) we have shown and for positive s for . Both of these functions have decreasing behavior in terms of the carrier density. In the case , approaches continuously to which admits larger carrier densities that is consistent with larger mean densities .

Before closing this section, let us turn to calculating the compressibility (). Knowing that the compressibility is proportional to we can easily do the integrals and show that, for the un-gated intrinsic graphene , . This shows that in the limit , (since ) for which a transition is possible. Based on these observations, one may conclude that this transition is the formation of EHP landscapes.

## Iv Two-body probability density

In this section we consider a two-body functional. i. e. . As stated in previous sections, . The infinitesimal change of caused by the infinitesimal growth of the charge density is . According to the previous section we have:

(19) |

Just like the calculations of the previous section we write (setting and defining , , ):

(20) |

from which, after iteration by parts and symmetrization and using also the definition , we reach at the master equation for :

(21) |

in which:

(22) |

In the previous section we have calculated and to be respectively and , leading to the relation ( and ):

(23) |

Let us consider the case or , for which and constant. To facilitate the procedure we restrict the calculations to positive densities (so that ), having in mind that the same calculations should be done for negative densities (for which a minus sign is necessary). In this case become very large and we have:

(24) |

in which and . We try the solution . Using this relation and the fact that we find that satisfies the following equation:

(25) |

in which , , and are the dimensionless parameters. To solve this equation we consider to be composed of two parts, i.e. in which and . Then we obtain two linear differential equations as follows:

(26) |

in which and so on, and is an arbitrary real number required in the separation of variables method. This form is not completely a factorized form since, as is seen the equations are not independent due to the presence of the common factor . The second equation has been written in a form which is more suitable for our analysis. The initial form has been , and we have used the reflection symmetry (note that has been which changes sign under the mentioned operation). This means that we have two types of solution which is important in our analysis. For a while we suppose that and seek for the solution of the above equation. Let us consider the limit of very small distances, namely for which one can safely ignore in the first equation and in the second equation with respect to the other terms. The general solution of the equation and the same for is hypergeometric and s:

(27) |

in which is the Kummer confluent hypergeometric function and is the Meijer G function Andrews (1985), and , , and are some constants to be determined by the boundary conditions. For and , and for . These functions for and have been shown in Fig. 2(a). To decide about the solutions we note that as () the statistics of two disjoint parts become independent and we expect . In the opposite limit, for () , so that for , in this limit, i.e. . These properties are only satisfied for , showing that all for all ’s and for . Noting that for and we see that for . The same solution can be done for the case with interchanging the roles of and . This shows that has the following important behavior:

(28) |

in which is some constant to be determined by normalization of and and are the solutions of equation 17, i.e. one-body distribution functions. It is interestingly seen that for the limit (which is beyond the approximation which was assumed above) this solution reduces to a factorized form of equation 17 which is expected since in this limit the particles behave as disjoint independent particles and the two-body distribution function should be multiplication of two one-body distribution function. The other important limit is , for which the factor guarantees that in this limit tends to zero except for the case as expected.

A very important quantity is the density-density correlation function which, as stated in the previous sections has very important information about the density response of the system to an external potential. It can easily be seen that the density-density correlation function is obtained to be , in which:

in which and we have used the new variable . We have evaluated numerically this integral using the iterative integration method with summation steps and precision ratio in each direction (one for the -direction and another for -direction). The result has been shown in Fig. 2(b) in which