Interaction-disorder-driven characteristic momentum in graphene, approach of multi-body distribution functions.

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 (ζ-1) is related to the disorder strength (ni), the interaction parameter (rs) and the average density ([Formula: see text]) via the relation [Formula: see text] for which [Formula: see text] 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 qch 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.

Introduction

Since the discovery of graphene as a hybrid between metal and insulator, this material with relativistic energy spectrum of zero-gap Dirac Fermions[1,2] attracts attention due to its unique electronic properties and prospective applications in nanoelectronics[3-6]. 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[5]. There is a huge literature concerning the interplay of particle interactions and disorder in mono and multilayer graphene[4,7-13] which cause variety of phenomenons. Examples of the effect of disorder on the properties of graphene are its effect on: the compressibility[14], the electron-phonon interaction[15], the Fano factor and conductivity[5], the magnetoresistance of bilayer graphene[16], the graphene Hall bars[17], the optical properties of graphene quantum dots[18] and the effective electronic mass in bilayer graphene[19]. Disorder-based spintronics in graphene[20], and Levi-flight transport caused by anisotropically distributed on-site impurities[21] are other examples for which the disorder plays a dominant role. The disorder-mediated Kondo effect[22] 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.[23] and Adam et al.[24] and was also confirmed in experiments in the vicinity of the Dirac point[25-34]. 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[35]. In this case (for which the transport is governed by the complex network of small random puddles with semi-metal character[35]) 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[36] 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 [7,36]. 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[6].

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[37,38]. 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[6,39] and no collective modes are allowed within the PRA approximation. For finite μ however, the more general solution of Shung[40] is applicable which predicts that the potential is screened with the screening length (v ≡ Fermi velocity, k ≡ Fermi momentum, ≡ background dielectric constant). This predicts the Friedel oscillations in graphene to be of the form [40]. 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[6]. These results were obtained for disorder-free system and within the RPA approach which is questionable for the graphene as pointed out by Mishchenko[41]. 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[23,24,35,42]. 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 qch which is related to the disorder strength and inter-particle interactions via in which n is the disorder strength, r is the dimensionless interaction factor to be defined in the text and d 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. 2. Section 3 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. 4. In the conclusion section we close the paper by highlighting the main findings of the paper. Appendix 1 in the Supplementary Material contains the essential rules for calculating Feynman diagrams for the TFD theory. Appendix 2 and 3 help to facilitate some calculations over the paper.

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[41]. 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[6]. In the strong interaction side however, it is shown that the graphene has non-Fermi liquid behaviors with power-law quasiparticle dispersion[6,43], 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[6]. 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[7,36]. 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 T 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[44]:in which the inner is the quantum expectation value, whereas the outer one stands for the disorder averaging. is the solution of a quantum model, which is the Thomas-Fermi-Dirac theory in this paper. 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. In the disorder-free (and also quenched disordered) systems the contribution of the first part is apparently trivial, as the averaging over quenched disorder of disconnected diagrams does not produce any contributions to the physical observables, in contrast to connecting the loops by interaction lines (and in contrast to the averaging over annealed disorder). Here we decorate this term by some external (disorder) as well as Coulomb interaction lines according to TFD theory, whose effects are non-trivial. In the other words in our effective TFD diagrams we replace the first term by some diagrams, interacting with external sources and find its contribution self-consistently.

Consider the ground state, and 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. . Then, we define the first term of this equation as in which shows the averaging over the disorder. Figure 1 shows this equation diagrammatically, in which the second part is connected term which has extensively been analyzed in the literature[5,45]. In this figure the single-line circles show the carrier density at the point x, i.e. n(x) and the shaded area shows all interaction effects, such as Coulomb lines and external (disordered) interaction lines. In the present paper we mainly deal with which also carries the effects of disorder and Coulomb interactions. We show this function diagrammatically by a zigzag and two full (double) circles as shown in Fig. 2. The rules for the TFD diagrammatic expansion of this function have been developed in Appendix 1 (Supplementary Material). In this expansion the wavy lines between r and r′ are which is the electrostatic potential, the gray circles with a cross inside are the external potential (disorder), and the double circles with r are the full electronic density . The full-circles (double line circles) show the full density (see Appendix 1) which is obtained self-consistently according to the employed theory, which is the TFD equations here. We concentrate mainly on the contribution of this term and its effect on the response functions. We do not follow the diagrammatic expansion in Fig. 2, but instead we try to find this function self-consistently, and by using some stochastic phenomena techniques, and find its analytic form in some limits.

Figure 1

The diagrammatic representation of Eq. 1 for which the first term shows the connected component of the polarization operator, whereas the second term shows the disconnected component. Note that the disconnected term is also effectively connected via the external legs corresponding to the disorder which is averaged[44].

Figure 2

The diagrammatic representation of according to TFD theory. The zigzags between two points r and r′ (with two double circles at the ends) are , the wavy lines between r and r′ are 1/|r − r′| which is the electrostatic potential, the gray circles with a cross inside are the external potential (disorder), and the double circles with r are the full electronic density n(r). For the details see Appendix 1. In this expansion the first term on the right hand (second line) is f3 and the second term (third line) is f1 in the appendix.

The non-trivial effect of this term on this expectation value can be seen from the following relation[44] (the full treatment is postponed to Appendix 1 in the Supplementary Material):from which we see that the change of energy due to this term is in which and G 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 x. Therefore the exchange-correlation term becomes a local term in the Hamiltonian, and the only non-local terms are the Hartree and disorder terms. Appendix 1 in the Supplementary Material 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. 2 of Appendix 1 in Supplementary Material. The coulomb potential is screened only via the mediator as depicted in Fig. 4 of Appendix 1 in Supplementary Material, 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 (Supplementary Material). 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 r′), the extra contribution comes from the quantities like (in which f is an arbitrary function). In such cases we define the pair (or two-body) distribution measure to write the mentioned correlation function as[46]:for which is a special case. is defined as the probability density of charge density to be n at r and n′ at r′. 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 t is its thickness and V is the gate voltage. The experimental data show a strong dependence on in which n is the impurity density. In ordinary densities, the conductivity is linear function of x and for very low x’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[5]:in which v 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[47]:in which k is the momentum cut-off and . The remote Coulomb disorder potential is calculated by the relation[42]:in which is the charged impurity density and is the distance between substrate and the graphene sheet. For the graphene on the SiO2 substrate, , so that , nm, where a0 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[23]: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 V. 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 V makes things difficult, since in which . Therefore the rescaled equation is[46]:in which . Therefore the first term survive marginally in the infra-red limit and the scale invariance is expected, even in the presence of V.

One-Body Probability Density

In this section we concentrate on calculating one-body probability density P by focusing on Eq. 7. From now on stands for the disorder averaging. Using the identities and we reach to the following identity[46]:which is essential in the following sections. In this equation and and , n. An important point is superficial contradiction that we have used which is ignored in the TFD theory, since it is valid only when (k 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:

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

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:

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

Now let us concentrate on and . The calculation of these quantities is done to Appendix 1 (Supplementary Material) 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:

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:

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 n(r) as follows: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 Fs. In fact F is everywhere non-zero and well-defined and finite except at n = 0 and for which becomes zero and negative infinity respectively. Therefore Eq. 16 is everywhere reliable except at n = 0 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 sgn(n) |n|3/2 whose derivative is not well-behaved in n = 0.

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 ζ0 limit. In this limit, and considering to be nearly constant, we have in which . The solution is therefore ():in which S is the area of the sample and P0 is the normalization constant. This relation may seem to be unsuited, since it grows unboundedly for negative values of n. Actually there is no contradiction, due to the presence of G whose amount grows negatively for negative n 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 P. Our approximation (considering G 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 P.

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 G favors large charge fluctuations. To see how G controls the charge fluctuations, we note that the limit (zero-gated graphene) has a direct effect on G. In fact μ controls 〈n〉 which directly affects which is G. In the limit , we expect that G 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:

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. 3a we have shown a plot which compares rF 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. 3b we have shown and P for positive ns for . Both of these functions have decreasing behavior in terms of the carrier density. In the case , P approaches continuously to which admits larger carrier densities that is consistent with larger mean densities 〈n〉. Before closing this section, let us turn to calculating the compressibility (κ). Knowing that the compressibility is proportional to 〈n2〉 − 〈n〉2 we can easily do the integrals and show that, for the un-gated intrinsic graphene , . This shows that in the limit 〈n〉 → 0, (since ) for which a transition is possible. Based on these observations, one may conclude that this transition is the formation of EHP landscapes.

Figure 3

(a) A comparison between rF and to see the validity region of the solution. (b) Log-log plot of the solutions and P in a valid range.

Two-body probability density

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

Just like the calculations of the previous section we write (setting and defining , , ):from which, after iteration by parts and symmetrization and using also the definition , we reach at the master equation for P:in which:

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

Let us consider the case or , for which and G ≈ 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:in which and . We try the solution . Using this relation and the fact that + we find that P0 satisfies the following equation:in which , , and are the dimensionless parameters. To solve this equation we consider P0 to be composed of two parts, i.e. in which and . Then we obtain two linear differential equations as follows: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 R. 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 P0 is hypergeometric and s:in which 1F1 is the Kummer confluent hypergeometric function and is the Meijer G function[48], and A, B, and are some constants to be determined by the boundary conditions. For λ = 0, and , and for , 1F1. These functions for λ = 0, 1 and 2 have been shown in Fig. 4a. 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 λ = 0 and we see that for . The same solution can be done for the case with interchanging the roles of x and x′. This shows that P has the following important behavior:in which A is some constant to be determined by normalization of P and and are the solutions of Eq. 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 Eq. 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 P tends to zero except for the case as expected.

Figure 4

(a) The plot of Kummer confluent hypergeometric and the Meijer G functions for various λ’s. Higher λ functions have some nodes which is forbidden. (b) The plot of in terms of R for various rates of ζ. Inset: The same plot in the semi-logarithmic form in which the − is evident. (c) The log-log plot of δ in terms of ζ. Inset: The log-log plots of G(Rmax) and Rmax in terms of ζ. The power-law behaviors are evident in all of these graphs.

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 10−5 in each direction (one for the U-direction and another for n-direction). The result has been shown in Fig. 4b in which has been sketched for various rates of ζ (Rmax is the amount of R in which G(R) shows a peak). It can be shown that as stated in the end of the previous section. This function tends to zero as as expected. G(R) has some interesting features. Firstly the tail of this function (for large R’s) behaves in a logarithmic fashion. To show this we have shown / in terms of re-scaled R/Rmean in a semi-logarithmic graph in the inset of Fig. 4b for which a − slope is evident. Except for small ζs (ζ = 1 in this graph), the graphs fit properly to each other with a nearly the same slope. Clearly this behavior can’t remain for very large Rs, since G(R) should be positive in all R range. The logarithmic behavior in graphene has been observed in many aspects, e.g. the logarithmic enhancement of Fermi velocity[39], or exchange-driven singularity in the Fermi velocity in intrinsic graphene[49], but up to the author’s knowledge such a behavior in two-point correlation functions has not been reported before. The logarithmic behavior is the characteristics of free Boson systems that our effective model for intermediate spatial scales corresponds to. The other feature of G(R) is the narrowing of its peak as ζ increases, and that Rmax and G(Rmax) are decreasing functions of ζ. The change of these quantities has been shown in the Fig. 4c. It is interestingly seen that these quantities show some clean power-law behaviors in terms of ζ. The power-law behavior arises from the scale-invariance of the TFD equation which was analyzed at the end of SEC. 1. To characterize the behavior of G(R) let us analyze the width of the G(R) defined by in which , and also G(Rmax) and Rmax in terms of ζ. The behaviors are , and in which , and , which arises to the hyper-scaling relation:

In the limit in which the interaction is strong enough and (or) the fraction 〈x〉 is small enough, i.e. ζ is small, Rmax becomes large and therefore δ increases unboundedly which leads to large density inhomogeneity, i.e. formation of EHPs.

The other important quantity is the Fourier transform of G(R), i.e. G(q) which is shown to be:in which J(x) is the Bessel function of the first kind, and . Note that the imaginary part of G(q) vanishes. The presence of disorder induces a characteristic wave vector (so that ), which implies that . In this relation we have used and . has been sketched in Fig. 5a for various rates of ζ in a semi-logarithmic scale. For large ζ values this function is nearly constant for small values of and starts to falling off rapidly at some q value which we name q0. For q values in the order of (and larger than) qch, we enter the phase in which is vanishingly small. q0 is ζ-dependent and increases as ζ increases. It is interestingly seen that and q0 show power-law behaviors in terms of ζ with exponents shown in the figure, i.e. and in which and .

Figure 5

(a) The semi-log plot of . In this graph q0 is defined as the point above which the graph falls off rapidly. (b) The power-law behavior of in terms of ζ with the exponent . Inset: The log-log plot of in terms of ζ with the exponent . (c) The real-space screened potential which is deformed significantly for smaller ζ’s. For the potential becomes negative for some R intervals which is the finger print of some instability. Inset: The semi-log plot of the Fourier component of the screened potential. The potential if deformed most significantly for small q’s and is magnified for smaller ζ’s, showing that this deformation is disorder-driven. (d) The real-space screening of charge impurity and its oscillatory behavior for small ζ’s. The right inset shows the screening at R = 0. It is seen that for changes sign which signals the screening-anti-screening transition. The average wave length of the oscillatory behavior has been shown in the left inset.

To understand where these power-law behaviors come from, we should look at the main equation of n(r), i.e. Eq. 7. This behavior has its roots in the scale invariance of this equation (to see this symmetry let and which yields the same equation with a renormalized logarithmic-enhanced coefficient) which cause a scale symmetry in Eq. 28. In fact if we do the transformations , and , then we see that .

The effect of G(q) on the screened potential is very interesting. The Fourier component of the screened potential (vsc(q)) in long wave lengths (small q with respect to qch) is crucially changed, so that , whereas the form of vsc is not changed with respect to the bare interaction potential for large and intermediate wave numbers. It has been shown in the inset of Fig. 5c from which it is seen that all curves are fitted to 1/q for large q’s. The cross-over between these two behaviors occur in some q interval. It is seen that there is a characteristic q there (the peak), which results to some oscillations in the real space. The potential in the real space v(R) has some interesting features. For large R’s and large ζ’s this quantity falls off just like the bare potential, i.e. 1/R with some oscillations which is evident in Fig. 5c (which results from the peak in the Fourier component). In the small R limit the potential tends to a constant value, which cause to be finite. We see that as ζ increases, also increases, as expected for the free system. Inversely, when ζ decreases () the behaviors change crucially as is seen in the figure and becomes negative in some R intervals, signaling the instability of the Dirac gas. We think that this signals the EHP phase in which the density fluctuations become larger than the average density. This change of behavior can also be seen in the response function which has been shown in Fig. 5d. In the right inset − has been shown in terms of ζ, from which a screening-anti-screening transition is observed for . This corresponds to the change of potential sign, as explained above. An oscillatory behavior appears with the wave length λaverage shown in the left inst of this figure for small ζ′s. The corresponding wave numbers are consistent with the ones for which the potential vsc(q) shows a peak. The emergent oscillatory behavior is the effect of characteristic wave number for which the screened potential shows peak.

Before closing this section, it is worth mentioning some points concerning the form of G(q). It is well-known that for zero temperature case the polarization given by the RPA for mono-layer garphene is constant for and increases nearly linearly for . This is in sharp contrast to the ordinary two-dimensional electron gas (2DEG) in which for , it is a decreasing function of q with a discontinuity in the derivative at [5]. G(q) (as the contribution of the second term in Fig. 1) shows actually an opposite behavior with the other reference point qch which is more similar to the behavior of ordinary 2DEG. Therefore one may be interested in the case in which the second term dominates the first term in Fig. 1, which cause the graphene to behave like ordinary 2DEG. This may be possible for small enough ζs (strong interactions, strong disorder and small densities) for which is larger and electron-hole puddles are present. The total information concerning the two-body charge density distribution function has been gathered in Table 1.

Table 1

The critical exponents of graphene related to the two-body charge density distribution function with their definitions.

Exponent	Definition	Closest fractional amount
a G	G_x(R)∝a_Glog(R)	13
τ G	G_x(R_max)∝ζ^−τ_G	72
τ R	R_max∝ζ^−τ_R	2
τ δ	δ∝ζ^−τ_δ	τG/2
τ_G˜	α^2G˜(q=0)∝ζ^−τ_G˜	92
τ q	α^−1q_0∝ζ^τ_q	12
γ δ, R	δ∝Rmaxγ_δ,R	12

The charge probability measure

In this section we seek for the equation governing P({n}) which is the probability measure for the charge configuration {n}. This function actually depends not only to n configuration, but also to its gradient: . But as we will see the dependence to the local charges is enough when we are looking at the large carrier densities. Let us suppose that f is an arbitrary local or non-local smooth function of the density n, i.e. in which r0 is some reference point. Without loose of generality we set to facilitate the calculations, having in mind that there is some reference point from which the positions are calculated. Using the the equations of SEC 1, it can therefore be expanded in terms n (defining and ):where in the second line we have used the fact that is equivalent to due to homogeneity of the system (of course after averaging over disorder). We can calculate the probability measure of the density, noting that the average value of f should not depend on r0 due to homogeneity of the system. To do averages let us define . The change of the average, due to changing the origin is:

Noting again that , and representing the averages as the integrals over probability measures, and doing integration by parts, we reach at:Let us return to the Eq. 16 and for convenience define again the used functions as follows:in which , and . Now let us try the solution:in which depends on r′ and r″ via the densities, i.e. . We obtain:in which we have used the symmetry . The right hand side of the above equation should be equal to the right hand side of Eq. 34, i.e.Therefore defining , , , , , , , and we obtain:in which the second equation is obtained by symmetry. Let us first consider the case for which we obtain and the same equation for . By summing two terms we obtain:in which and . To continue some relations are essential, namely ():which is used to re-write Eq. 39 in the following form:in which and and . The general solution of the above equation is as follows:in which is an arbitrary function having the property to restore the symmetry of H and c is some constant. It is simply seen that there is no choice but .

Now let us to turn to the Eq. 38. To satisfy these equations, one should have an asymmetric term to generate the delta function. Satisfying this discontinuity (which is proportional to ) is a difficult task and the inclusion of n is needed. Fortunately for the limit or we can neglect this term and the term corresponding to the discontinuity becomes very small. Therefore the above analysis is applicable and the final form of H for small densities is as follows:With this solution, the exponent of the Eq. 35, i.e. is well-defined and smooth in all range of n and n′.

Conclusion

In this paper we have investigated the properties of the disorder-averaged dielectric function by analyzing the one- and two- body charge density distribution functions of mono-layer graphene. For calculating them, we have used the Thomas-Fermi-Dirac (TFD) approach, taking into account the tunable disorder and inter particle interactions on an equal footing. To use the TFD theory for calculating the linear screening, we developed a diagrammatic technique and in the first order approximation in disorder strength, the polarization operator and charge modulation due to external potential have been obtained. The one- and the two- body charge density distribution functions were obtained using some stochastic analysis which carries both the effect of white-noise out-plane disorder and the interaction. The closed form of the functions were obtained for high density and low inter-particle interaction limits. By analyzing the one-body distribution function we found that the charge density fluctuations is proportional to in which n is the disorder strength, r is the interaction strength and is the average density of the system, therefore the electron-hole puddles can appear in the low density limit . By analyzing the two-body charge density distribution function we obtained that:in which is the one-body distribution function in the large ζ limit. Using this function and some other analytical investigations we showed that some power-law behaviors in terms of ζ raise, specially for Rmax which is the length scale at which the static density-density correlation function G(R) shows a peak. It was found that G(R) falls off logarithmically for large R’s, i.e. . The exponents of the scaling relations have been reported in the text. An interesting effect of disorder and interaction is changing the screened potential both in the Fourier and direct spaces. It was that in addition to k a disorder-driven characteristic momentum (qch) emerges in the system that controls the behaviors of the screened potential. The Fourier component of screened potential changes significantly for small wave numbers. for small enough ζ’s (small densities where EHPs can appear) we have observed an instability in which the imaginary part of the dielectric function becomes negative and some oscillations are observed for the charge density δn(r). The oscillations are absent for larger ζ values (larger densities or smaller inter-particle interactions and disorder strengths). Along with this change, a screening-anti-screening transition occurs in which δn(0) changes the sign. The observed instability is certainly affected by the higher order contribution, but we think that it should show itself for small enough densities, since the effect increases in a power-law fashion with ζ.

In the last part of the paper we have calculated the total distribution measure of the mono-layer graphene sheets. The resulting equation has been analytically solved for large ζ limit and was shown to be of the quadratic form. For smaller ζ’s we think that derivatives of charge densities should come to calculations which is the subject of our future research.

Interaction-disorder-driven characteristic momentum in graphene, approach of multi-body distribution functions