Dispersion forces between weakly disordered Van der Waals crystals

We describe a many-body theory for interlayer dispersion forces between weakly disordered atomically thin crystals and numerically investigate the role of disorder for different layer-separation distances and for different densities of induced electrons and holes. In contrast to the common wisdom that disorder tends to enhance the importance of Coulomb interactions in Fermi liquids, we find that short range disorder tends to {\it weaken} interlayer dispersion forces. This is in line with previous findings that suggest that transitioning from metallic to insulating propagation weakens interlayer dispersion forces. We demonstrate that disorder alters the scaling laws of dispersion forces and we comment on the role of the maximally crossed vertex-correction diagrams responsible for logarithmic divergences in the resistivity of two-dimensional metals.


I. INTRODUCTION
Even when two objects are each electrically neutral, forces between the two objects which are mediated by the electromagnetic field can still be present. These dispersion forces were named by London in his theoretical investigation of forces between molecules [1]. Although each molecule has zero total charge, quantum fluctuations in the charge density of each molecule lead to an effective dipole-dipole intermolecular force. This mechanism was later generalized by Lifshitz [2,3] to describe forces between solids, wherein he discovered a force which scales like 1/d 3 when the distance d between two thick slabs becomes large. Depending on the context, these forces also go under the name of van der Waals or Casimir forces, where the former (latter) often indicates that the force is mediated by the longitudinal (transverse) component of the electromagnetic gauge field [4].
Dispersion forces are relatively weak and short ranged compared to electrostatic forces, and are difficult to observe in experiments on solids. Recently however, advances in x-ray spectroscopy have allowed for atomiclevel precision measurements of interlayer strain in thin films and atomically thin crystals [5,6], and signatures consistent with interlayer dispersion forces among optically induced electrons and holes have been measured in transition-metal dichalcogenide multilayers [7]. This adds a new, experimentally measurable quantity to the class of phenomena which are sensitive to correlations among quasiparticles in neighboring layers of atomically thin crystals like transition-metal dichalcogenides, graphene, twisted bilayer graphene, and phosphorene. Coulomb drag [8] is a notable example of the type of phenomena which are sensitive to interlayer correlations. In * jvmilczewski@mpq.mpg.de these experiments a current is driven in one layer and as a result of interlayer Coulomb interactions an induced voltage drop appears in a second (otherwise passive) nearby layer. Drag experiments have led to a deeper understanding of the nature of the elementary excitations and ground state wave functions of complex phases of matter, from two-dimensional Fermi liquids to more exotic phases like exciton condensates [9,10] and Luttinger liquids [11]. Just like Coulomb drag, the interlayer dispersion force between atomically thin crystals offers an interesting test bed for the various many-body theories describing the complex behavior of solids.
In this paper we construct many-body approximations to explore the impact of weak disorder on the interlayer dispersion forces which act between layers of a bilayer heterostructure after a finite density of electrons and holes are induced in each layer as illustrated in Fig. 1. While ab initio methods for obtaining van der Waals contributions to the ground state en-ergy exist [12][13][14][15][16][17], our diagrammatic approach is sensitive to the exchange-correlation effects which densityfunctional theory usually deals with only on a meanfield level using variations of the local-density approximation; the approach discussed in this paper is complementary to these existing tools and allows for the treatment of systems with strongly correlated ground states or, as we investigate in detail below, random disorder. Quasiparticle-impurity interactions are known to be responsible for a number of fascinating properties of metals, from weak-localization corrections to the longitudinal conductivity [18] to anomalies in the tunneling conductivity [19,20], and we will make use of some of these welldeveloped many-body approximations in determining the role of weak disorder on interlayer dispersion forces.
Our paper is organized as follows. In Sec. II we describe a many-body theory for the interlayer dispersion force based on a linked-cluster expansion for the correlation energy of a bilayer in the absence of disorder as discussed previously [21,22]. In the limit of high quasiparticle density and large separation distance d, one recovers the well-known d −5/2 scaling behavior [22][23][24] in agreement with predictions from Quantum Monte Carlo methods [25]. In Sec. III we describe a leading-order-in-1/ε F τ theory for interlayer forces. We demonstrate that disorder qualitatively alters the scaling laws and demonstrate that disorder tends to reduce the magnitude of interlayer forces. In Sec. IV we discuss the impact on interlayer forces by a class of Feynman diagrams known to yield logarithmic divergences in the longitudinal resistivity of two-dimensional metals. Finally, in Sec. V we summarize our results and discuss interesting questions to be addressed in the future.

II. INDUCED DISPERSION FORCES IN BILAYER SYSTEMS
We consider a system governed by the following Hamiltonian: which describes the kinetic energy of electrons and holes, the Coulomb interaction, and the interaction of electrons and holes with impurities, respectively. We assume, as is often the case experimentally, that the density of induced electrons and holes (quasiparticles) is such that the kinetic energy of electrons and holes can be described by an effective mass approximation where ε α (k) = 2 k 2 /2m α and α is a composite index which labels the spin, valley, and band (e.g., valence vs conduction band) quantum numbers. In the following, we will consider the limit in which interlayer hopping is weak compared to the exchange-correlation energy per electron. Thus, the single-particle wave functions have a which-layer quantum number denoted by I. Interlayer hybridization of the conduction and valence bands is notoriously weak in van der Waals crystals (as the name suggests) and is often further weakened by rotational misalignment of neighboring layers. The charged quasiparticles in the various layers of the system interact with each other via the Coulomb interaction where The material-specific parameter κ describes the dielectric contributions of the elementary excitations outside of our model [e.g., phonons and propagation of electric field outside of the two-dimensional (2D) material]. The strength of Coulomb interactions is traditionally [26] described by the value of a parameter r s which expresses the ratio of average interaction energy to average kinetic energy in a disorder-free two-dimensional electron gas (2DEG), r s ∝ H e-e / H 0 . The parameter depends on the total density of electrons (and holes) in each layer n I and is larger when the density is lower, r s = a * B √ πn I −1 .
Here, a * B = κa B /m ef f is the effective Bohr radius. When the system contains particle populations described by different effective masses it is useful to define a * B using the geometric mean of the masses m ef f → √ m e m h . Interactions of charged quasiparticles in different layers are ultimately responsible for the induced van der Waals forces we describe. In this paper we consider densities of induced quasiparticles which are large enough to form electron liquids and hole liquids rather than excitons, as was recently demonstrated at room temperature [27].
The interaction between impurities of the crystal and electrons as well as holes is obtained by assuming that each impurity creates a deviation in the perfectly periodic scalar potential created by the underlying lattice. This scalar potential couples linearly to the density of electrons and holes, where ρ I (Q) is the Fourier transform of the density of impurities in layer I, and u I (Q) is the Fourier transform of the scalar potential of each impurity. We assume that electrons and holes only scatter off the impurity potential in the same layer, and we assume that the scalar potential is short ranged so that u I (Q) is actually independent of wave vector. The quasiparticle-impurity scattering time τ k can be defined using the Born approximation for the self-energy [28] where Σ(k, ω) = −i /2τ k . In the presence of finite disorder, the scattering rate at the Fermi energy is used to define the small parameter of our perturbation theory 1/(τ ε F ) 1, where we here (and will continue to) drop the subscript on τ .
Our method for evaluating the force between two atomically thin crystals consists of first calculating the ground state energy per layer as a function of interlayer separation distance d, and then calculating the force by taking the first derivative The ground state energy can be evaluated by taking the zero-temperature limit of the thermodynamic free energy Ω. The latter has a well-known perturbative formulation in the linked-cluster expansion [28] where ρ 0 is the noninteracting density matrix, T τ is the (imaginary) time-ordering operator, andV (τ ) = H e-imp (τ ) + H e-e (τ ) is the sum of the two interactions in our model within the interaction picture of time evolution [28]. By applying Wick's theorem, all contributions at order can be expressed in terms of integrals over noninteracting Green's functions, the Coulomb interaction V , and the electron-impurity interaction u I . One can now make use of Feynman diagram techniques to efficiently calculate these contributions. We now have all the tools necessary to evaluate the interlayer force to any order in perturbation theory. Before we consider the effects of weak disorder on the interlayer forces, we reproduce the well-known d −5/2 scaling of energy [22][23][24][25]29] by examining the force between two two-dimensional electron gases within the randomphase approximation (RPA) [30][31][32][33] and taking the limit of large interlayer distance d. We thus ignore disorder and takeV (τ ) = H e-e (τ ) within Eq. (7). The RPA can be understood as an expansion of the ground state energy in powers of the small parameter r s , and therefore gives a criterion for selecting which subset of Feynman diagrams at each order in within Eq. (7) must be included in an approximation to a given order in r s . The four lowest-order diagrams which contribute to the correlation energy are shown in Fig. (2). The full RPA approximation consists of summing all diagrams of this type, which at each order in contain bubble subdiagrams. The degeneracy of the diagrams in Fig. (2) is such that the infinite series of these types of diagrams can be resummed into a logarithm of a simple function of the single bubble diagram. After taking the derivative of the RPA approximation for the correlation energy [22], we obtain Feynman diagrams for the correlation energy of a bilayer system whose quasiparticles interact via intralayer Coulomb interactions (single wavy lines) and interlayer Coulomb interactions (double wavy lines). Only the four lowest-order diagrams are shown here. Solid lines with arrows represent noninteracting Green's functions of quasiparticles.
the following integral expression for the force per layer between a bilayer system containing a finite density of electrons and holes in each layer Here, χ 0 is represented by the bubble subdiagrams found in the four diagrams in Fig. (2) and describes the noninteracting density-density response function of each layer. The zero-temperature limit of χ 0 can be evaluated for parabolic-band effective mass models, and in the presence of both valence and conduction bands, χ 0 = α χ α 0 , where χ α 0 is the Lindhard function [34] of the α−particle species. The integral over frequency in Eq. (8) is over the imaginary frequency axis, and the arguments of χ α 0 (q, iω) have been omitted for brevity.
The application of Eq. (8) assumes that thermal equilibrium has been reached among the electrons and holes, which is usually several orders of magnitude faster than the electron-hole recombination time, and does not limit experimental observations. For arbitrary electron/hole densities and interlayer separation distances, Eq. (8) must be evaluated numerically. Furthermore, it should be mentioned that Eq. (8) leads to a nonvanishing force even in the absence of holes.
In Fig. (3) we present the results of numerical calculations for the pressure (i.e., force per area) between two layers of atomically thin crystals with induced densities of electrons and holes parametrized by r s . We immediately notice that the force between layers is attractive and that the magnitude varies dramatically with interlayer separation distance as V 12 is dependent on d. This is a particular feature of the type of dispersion force that derives from the instantaneous Coulomb interaction (typically called van der Waals forces) instead of forces originating from the transverse and retarded parts of the elec- A plot of interlayer forces vs the density of induced quasiparticles in a disorder-free bilayer system. On the vertical axis is the force per area in units of effective Rydbergs per effective Bohr radius cubed. On the horizontal axis is the dimensionless parameter rs which is inversely proportional to the square root of the density of induced quasiparticles. Explicit definitions for rs and a * B can be found in the main text, while Ryd * = e 2 /κa * B .
tromagnetic field (typically called Casimir forces). While Casimir forces act at larger distances than van der Waals forces, they are significantly weaker and they are independent of the amount of impurities in the materials, and therefore are not addressed in this paper.
To demonstrate how the RPA theory obtains the known d −5/2 scaling for the energy [22][23][24][25]29], Eq. (8) is now evaluated in the limit of large interlayer separation. Specifically, we will find the leading-order contribution to the interlayer force in the small parameter F is the Fermi wave vector of the electron and hole Fermi seas which are present in each layer after excitation and thermalization. The presence of e −2qd in the numerator of Eq. (8) restricts the relevant range of q in the integral to q 1/d, which bears the physical interpretation that 2D in-plane charge perturbation waves at wavelengths which are short compared to the interlayer distance appear averaged out on the adjacent plate and thus will not contribute to forces. Long wavelengths, however, will not appear as averaged out and will therefore contribute to interlayer forces. In the limit k F d 1, the dominant contribution to interlayer forces will then come from long in-plane wavelengths and this thus restricts the relevant part of phase space to small values of q. In this region of phase space we are permitted to ap-proximate χ α 0 by its dynamic long-wavelength limit (i.e., ω > q, q → 0) which gives the leading-order contribution to the force. In the dynamic long-wavelength limit the noninteracting density-density response function of band α is given by where ρ α is the two-dimensional density of charged quasiparticles in band α. It is then straightforward to evaluate Eq. (8) analytically to obtain the leading order in 1/(k F d): which corresponds to the d −5/2 scaling for the energy.
Here, ξ 1 ≈ 0.315, ρ is the total two-dimensional quasiparticle density in each layer, and we have taken m h = m e = m, and κ = 1 for simplicity. Interestingly, in the case of infinitely many parallel plates (superlattice), the scaling of force per layer is identical to Eq. (10) up to redefinition of ξ 1 [7]. By randomly choosing two adjacent plates and identifying the gap between them as the gap between two semi-infinite thick slabs separated by a distance d, one can connect this result to Lifshitz' theory for thick, semi-infinite slabs. Introducing the threedimensional density ρ 3D = ρ/d in Eq. (10) to compare with Lifshitz' theory, we immediately see that we have reproduced the power law for the interlayer force in terms of interlayer separation and quasiparticle density (i.e., Despite the obvious utility of simple formulas like Eq. (10), the derivation demonstrates that only the longwavelength excitations (i.e., plasmons) are accounted for, while finite q excitations (e.g., noncoherent particle-hole excitations) are neglected. Indeed, Eq. (10) is only reasonable in the limit 1/(k F d) 1, and outside of this regime the interlayer forces are more accurately described by numerically evaluating Eq. (8). This is demonstrated in Fig. 4, where we show the ratio of pressure in the RPA approximation of Eq. (8) to the asymptotic form of Eq. (10). In the limit of k F d 1 the predictions coincide, while for smaller values of k F d the asymptotic form gives much higher interlayer attraction than the RPA form. In subsequent sections we will describe how these power law scalings are altered by the presence of impurities.

III. IMPACT OF DISORDER ON VDW FORCES: THE 'DIFFUSON'
In this section we lay out the basic elements of a manybody theory for the impact of weak disorder on the interlayer van der Waals (VDW) forces between atomically thin crystals. Specifically, we begin by introducing the small parameter (i.e., 1/ε F τ ) of the electron-impurity and hole-impurity interactions within the context of the first-order Born approximation (1BA) for the self-energy. We then identify the most relevant Feynman diagrams which contribute to interlayer dispersion forces within the regime of r s < 1/ε F τ . These diagrams contain an infinite series of ladder diagrams, and we discuss the solution of the Bethe-Salpeter equation for the vertex correction of the density-response function in the limit of short ranged impurity potentials. In contrast to the effect of disorder on other phenomena which arise due to interlayer interactions (e.g., Coulomb drag [35]), we find that disorder tends to weaken the magnitude of van der Waals forces.
The electron-impurity and hole-impurity scattering rates can be defined by the 1BA for the self-energy. In this approximation the self-energy is purely imaginary, Σ(k, ω) = −i /2τ k . For simplicity, we will take the hole's and electron's impurity scattering rates to be equal, although this condition is easily relaxed if required. The 1BA is given by the Feynman diagrams depicted in panel a) and b) of Figs. (5). Explicitly, the 1BA for the qindependent scattering rate at the Fermi energy is where ν α is the two-dimensional density of states at the Fermi surface of a single spin-and valley-resolved band, and ρ imp = lim Q→0 [ρ I (Q)]. In obtaining Eq. (11) we have made two standard approximations for treating quenched disorder in solids [28]. First, we assume that the impurity potential is short ranged, such that the Fourier transform of the potential which appears in Eq. (5), u I (Q), becomes independent of wave vector. Second, the impurity potential at any two different points is uncorrelated, such that the average over the probability distribution governing the impurity potential leads to ρ I (Q)ρ I (−Q) imp = N imp , where N imp is the number of impurities in layer I. Next, we consider how to incorporate quasiparticlequasiparticle interaction diagrams and quasiparticleimpurity interaction diagrams into an approximation for the dispersion force between atomically thin crystals. In the previous section we identified the leading-order-inr s contribution to interlayer forces as the derivative of the RPA diagrams for the ground state energy. In order to work with a well-controlled perturbation theory we will restrict our selection of diagrams to the case when r s 1/(τ ε F ). This allows us to obtain a well-controlled theory in both small parameters. The key is to not alter the order in r s of a diagram when adding any particular quasiparticle-impurity interaction line. We can accomplish this by adding to the RPA diagrams a nearly identical set of diagrams in which the noninteracting density-density response function bubble is dressed by quasiparticle-impurity interaction lines between the electron propagator and hole propagator which form each bubble. As long as these vertex-correction quasiparticle-impurity lines do not cross each other, they can be summed to infinite order and together they give the leading order in 1/(τ ε F ). The sum of all ladder Feynman diagrams for the density-density response function of each layer I is represented in panels c) and d) of Fig. (5). The latter is the diagrammatic representation of the Bethe-Salpeter equation and where G R/A (k, ω) = ω − −1 ξ kα ± i/2τ −1 and ξ kα = ε kα − ε F . The Bethe-Salpeter equation must usually be solved self-consistently for an arbitrary impurity potential. However, here it can be solved directly as a result of the bare-scattering amplitude being independent of momentum Γ 0 k,k = ρ imp |u I | 2 . In the regime where disorder gives significant contributions to the density-density response of a system, ω < 1/τ and q < 1/v F τ , straightforward calculations [38][39][40] yield Γ D (q, ω) = Γ 0 (q)/ −iωτ + τ Dq 2 where the diffusion constant is defined in two dimensions as D = v 2 F τ /2. The diffusion pole present in Γ D (q, ω) at ω = −iDq 2 is also present in the disordered density-density response function of layer I that is obtained by summing the diagrams in panel c) of Fig. (5) and yields where ν 0 is the total density of states at the Fermi energy in layer I. We can now evaluate the effect of weak disorder on the dispersion force between two atomically thin crystals by numerically evaluating Eq. (8) after replacing χ 0 (q, iω) by χ D (q, iω) in the region of phase space where ω < 1/τ and q < 1/v F τ . In Fig. (6) we plot the ratio of the interlayer force in the presence of disorder F dirty to the force in the absence of disorder F clean . We find that the interlayer attraction is reduced in magnitude by the presence of quasiparticle-impurity interactions, which we will analyze in more detail below. We also find that F dirty /F clean is reduced as d increases. This occurs due to the presence of e −2qd in Eq. (8) which originates from the form of the 2D in-plane Fourier transform of the interlayer Coulomb interaction. This factor restricts the density fluctuations which contribute to interlayer forces to wave vectors q 1/2d, and as d is increased, more of this region of phase space lies in the region governed by the disordered density-density response, q < 1/v F τ . We will now show that this phase space effect is also responsible for a change in the power-laws for the dispersion forces at large interlayer separation distances. In other words, in the presence of disorder, the asymptotic limit for forces  6. Ratio of the interlayer force in the presence of disorder F dirty to the interlayer force with no disorder F clean , plotted against the interaction parameter rs which is inversely proportional to the square root of the induced quasiparticle density in each layer of a bilayer. The three curves are for three different values of the interlayer separation distance d in units of the effective Bohr radius a * B . The degree of disorder is given by / F τ = 1/2. The values of both F dirty and F clean are calculated numerically using Eq. (8). In F dirty , the density response is given by the disordered limit χD(q, ω) for ω < 1/τ and q < 1/vFτ . between 2D planes presented in Eq. (10), F ∝ d −7/2 , is altered.
The numerical results presented in Fig. (6) show that disorder decreases the magnitude of interlayer forces. This is in contrast to the effect of disorder on other phenomena, like Coulomb drag [35], which also originates from interlayer quasiparticle-quasiparticle interactions. In the case of Coulomb drag, this conventional cartoon picture of the effect of disorder is that the change in the density-density response function from the noninteracting limit χ 0 (q, iω) to the disordered limit χ D (q, iω) represents a change from ballistic to diffusive motion of the quasiparticles. Indeed, the disordered density-density response function can be derived from semiclassical arguments using the diffusion equation [26], which is equivalent to the relaxation time approximation (RTA) [41] in the region ω < 1/τ, q < 1/v F τ in the dynamic limit. Since quasiparticles in neighboring layers which experience diffusive motion tend to spend longer periods of time near each other, they interact more strongly and this increases the Coulomb drag (i.e., disorder tends to enhance the transresistivity). However, since the interlayer forces are decreased in magnitude by the presence of disorder, we find that the cartoon picture of the ef-fect of disorder cannot be imported to understand our case of interest. The reason why disorder decreases interlayer forces while increasing the interlayer Coulomb drag is most simply identified by again examining the large-d limit of the two quantities. Specifically, while both Coulomb drag and the interlayer force depend on the density-density response function, the leading-orderin-1/(k F d) contribution to Coulomb drag comes from the static limit (ω < q, q → 0) of χ(q, iω) while the analogous contribution to the interlayer force comes from the dynamic limit (ω > q, q → 0) of χ(q, iω).
In the large-d limit our numerical results for the correlation energy per layer can be compared to previous investigations of disordered correlation energies within single-layer systems [42] where it was found that the introduction of disorder increases exchange energies in magnitude but decreases correlation energies in magnitude. By following similar steps as we took to derive the disorder-free expression presented in Eq. (10), we find the following leading-order expression: where ξ 2 ≈ 0.768 and ρ is the total two-dimensional density of quasiparticles in each layer and we have again taken m e = m h = m, κ = 1 for simplicity. Notice that the interlayer force now decays more quickly with distance than in the absence of disorder. This qualitative change is a direct result of the transition of electron and hole propagation from ballistic to diffusive. While it might be surprising that the effect of disorder on the interlayer forces is opposite to its effect on interlayer Coulomb drag, this behavior actually fits nicely into a trend observed in other systems [43,44]: the less metallic a system is, the faster its energy (and therefore its pressure) decreases with interlayer separation. Concretely, for a metallic sample, the energy scales as d −5/2 [22][23][24][25]29] while for a combination of a graphene and a metallic plate it scales as log(d)d −3 [43] and for two graphene plates is scales as d −3 [43]. Finally, for two insulator system it scales as d −4 [13]. The change of the scaling of the distance-dependent part of the correlation energy from d −5/2 to d −3 upon changing from ballistic to diffusive propagation thus confirms this picture.

IV. QUANTUM INTERFERENCE EFFECTS ON VDW FORCES: THE 'COOPERON'
In the previous section we developed a theory for interlayer dispersion forces between the layers of a bilayer system of atomically thin crystals which have uncorrelated and short ranged disorder. We summed an infinite set of Feynman diagrams by solving the Bethe-Salpeter equation and thus obtained the diffuson vertex correction of the density-density response function to leading-order in 1/( F τ ). In this section we will sum the class of diagrams which corresponds to the subleading-order terms for the interlayer dispersion force in powers of 1/ε F τ . These diagrams are familiar from the theory of weaklocalization and together they constitute the cooperon vertex-correction. Despite being of lower order in the small parameter governing the impurity-quasiparticle interaction, they are known to be responsible for a logarithmic divergence in the longitudinal resistivity of twodimensional conductors [18], which motivates us to consider them here as well.
The cooperon contributions to the density-density response function are obtained by summing the 'maximally crossed' vertex-correction; this infinite set of diagrams is illustrated in panel c) of Fig. (7). These diagrams represent the quantum interference of a wave packet of charge density which interferes with itself while traversing along the time-reversed path. This requires the system to have a time-reversal symmetry present in order for phase coherence to be maintained in-between collisions of the wave packet with different impurities. As previously mentioned, these diagrams give a logarithmic divergence in the resistivity (which is proportional to the current-density response function), and indeed a similar phenomenon happens in our case of interest. Specifically, the subleading-order contribution to the density-density response function yields a logarithmic divergence in the diffusion constant. When both the diffuson and cooperon contributions to the density-density response function are included [45], the functional form of χ D (q, ω) remains the same as presented in the last section except that D gets The fractional change in the interlayer dispersion force when the maximally crossed (i.e., weak-localization) diagrams are included. Notably, the logarithmic divergence which appears in the longitudinal resistivity of twodimensional conductors is not present here. Instead, the cooperon diagrams have a similar, but weaker, effect as the diffuson diagrams, where both tend to reduce the magnitude of interlayer attractive forces. The fractional changes are shown at / F τ = 1/2, kF d = 10 as a function of rs for different values of τ /τ0. an additional contribution which depends on frequency where ν 0 is the total two-dimensional density of states of all quasiparticles in layer I. Just as in the case of the cooperon contribution to the longitudinal resistivity, the logarithmic divergence we obtain is cutoff at long distances, or small momenta, by the inelastic scattering time of the quasiparticles, τ 0 . This time-scale is determined, for example, by the quasiparticle-quasiparticle scattering rate, and is responsible for destroying the phase coherence of the propagating (and time-reversed propagating) wave packet on very long time scales τ 0 > τ . This form of the disordered response function is only a reasonable approximation in the range where ω < 1/τ and q < 1/v F τ . We numerically evaluate the interlayer dispersion forces using the disordered density-density response function including the renormalized diffusion constant D → D + δD(ω). The results are shown in Fig. (8). They demonstrate that the maximally crossed diagrams tend to further reduce the magnitude of interlayer forces. More surprisingly, perhaps, there is no logarithmic divergence in the interlayer force, in contrast to what happens when using the analogous approximation for the longi-tudinal conductivity. This is surprising in light of the well-known relationship σ dc = lim q→0 (q 2 /ω 2 )χ D (q, ω), which follows from the presence of global gauge symmetry. However, while the conductivity is the response of the system to an external electric field whose frequency we can always fix to zero, in contrast, the interlayer dispersion force is an integral over all frequencies of density fluctuations in both layers (it is the Coulomb interaction between these density fluctuations which yields the dispersion force). And when the logarithmic divergence in χ D (q, ω) is integrated over frequency, it results simply in a finite reduction (on the order of 10 percent) of the interlayer force's magnitude.

V. SUMMARY AND DISCUSSION
We developed a many-body theory for the dispersion forces between atomically thin crystals with weak disorder. Such systems can be realized within van der Waals crystals [46] (e.g., graphene, transition-metal dichalcogenides, etc.) which form multilayer systems with very weak interlayer hybridization, a property which has allowed for optically induced interlayer strain, originating from dispersion forces, to be observed recently [7]. In these systems dispersion forces arise due to Coulomb interactions between fluctuations in the charge density of neighboring layers. The linked-cluster expansion method was used to approximate the correlation energy of a bilayer system and the force between the layers of the bilayer system was obtained by taking a derivative of the correlation energy with respect to interlayer separation distance.
In the high-density limit, the random-phase approximation bubble diagrams give the leading-order contribution to the disorder-free interlayer dispersion force. To account for disorder, we have summed an infinite series of ladder diagrams by solving the Bethe-Salpeter equation. These ladder diagrams form the diffuson contribution to the vertex correction of the density-density response function (i.e., the bubble), and yield the leading-orderin-1/(ε F τ ) theory. Numerical evaluation of the interlayer dispersion force shows that interlayer forces are weakened by disorder. On one hand, this is in contrast to the more conventional case [47] in which Coulomb interactions become more important when electron motion becomes diffusive rather than ballistic. On the other hand, this behavior is in accordance with previously observed changes in scaling laws as one transitions from metallic to insulating electron propagation.
We explain this behavior by considering the analytic structure of the density-response function in the small frequency and wave-vector limit. We find that the diffusive motion of electrons and holes leads to a qualitative change in the scaling laws for the interlayer dispersion force as a function of quasiparticle density and interlayer separation distance. Subsequently, the impact of the higher-order vertex-correction diagrams was investi-gated. Specifically, maximally crossed diagrams which are known to produce logarithmic divergences in the longitudinal resistivity of two-dimensional metals (i.e., weak localization diagrams) are found to be much less important for interlayer dispersion forces.
All the calculations shown in this paper were carried out within a bilayer system consisting of two parallel plates. It should be mentioned, however, that the effects of the theories developed in this paper were all at the level of "same-layer"-density-density response functions. As the theory of the bilayer system can easily be generalized to the theory of a superlattice system [7], the results of this paper can easily be transferred to the superlattice system with similar effects (e.g., same power laws and qualitative effects).
Optical control of electron and hole populations yields a convenient control knob for manipulating the interlayer separation distance of van der Waals crystals. In future calculations one may investigate the possibility of inducing interlayer dispersion forces by doping heterostructures with electrostatic gates. While these systems include interlayer electrostatic forces which compete with dispersion forces, the latter are not reliant on equal populations of electrons and holes and can hopefully still be observed. Through electrostatic gating the role of the excitonic spectrum in the formation of strains could be differentiated from the induced strains presented in our work. In order to complement this investigation of the role of excitons, it would furthermore be interesting to study the qualitative changes in interlayer dispersion forces which are present in multilayer systems with more exotic ground state wave functions, such as are present in bilayer exciton condensates.