Inhomogeneous phases in coupled electron-hole bilayer graphene sheets: Charge Density Waves and Coupled Wigner Crystals

Recently proposed accurate correlation energies are used to determine the phase diagram of strongly coupled electron-hole graphene bilayers. The control parameters of the phase diagram are the charge carrier density and the insulating barrier thickness separating the bilayers. In addition to the electron-hole superfluid phase we find two new inhomogeneous ground states, a one dimensional charge density wave phase and a coupled electron-hole Wigner crystal. The elementary crystal structure of bilayer graphene plays no role in generating these new quantum phases, which are completely determined by the electrons and holes interacting through the Coulomb interaction. The experimental parameters for the new phases lie within attainable ranges and therefore coupled electron-hole bilayer graphene presents itself as an experimental system where novel emergent many-body phases can be realized.

The interplay between superconducting and charge density wave (CDW) phases that is often observed in connection with High-Temperature superconductors, is attracting considerable attention. It has been argued that details of the very elaborate crystal structures typical of High-Temperature superconductors play a central role in determining the properties of the CDWs. A polarizable background is needed to drive the CDW phase, and the crystal lattice provides this background. However the complexity of the crystals in these materials, makes it challenging to identify the competing contributions [1][2][3][4][5] . We report in this manuscript on a, by far, simpler system exhibiting the same association of superfluid and CDW phases, but a system in which the polarizable background is uniform, so there is no intricate background structure that could determine the properties of the CDW.
There exist several studies pertinent to inhomogeneous phases in coupled two-dimensional-electron-gas (2DEG) layers [6][7][8] . Despite these early attempts, however, there are no studies of the CDW phase for small electron-hole layer separations where interlayer correlation effects will be strong. The earliest approaches could not treat strong interlayer coupling 7,8 , while for QMC calculations, it is impractical to study an inhomogeneous ground state like a CDW with oscillations of both unknown period and amplitude 6 .
Here we propose two strongly coupled two-dimensional (2D) bilayers of graphene (BLG), one bilayer containing electrons and the other holes, as an experimentally accessible system to observe the CDW phase. An electron and hole bilayer graphene has a quadratic dispersion around the Fermi level E F for densities ρ < 4 × 10 12 cm −2 9 . The electron and hole effective masses are matched, lying in the range . < < .   ⁎ m m 0 03 / 0 05 e 10,11 . A lower limit of experimental attainable densities is ρ > ∼ 10 10 cm −2 12-14 . The electron and hole densities are controlled by top and back metal gates. The electric field from the metal gates opens up an energy band gap between the parabolic conduction and valence bands 15 . We consider bilayer graphene embedded in a hexagonal Boron Nitride (h-BN) dielectric which is important to ensure high charge mobility and a large potential barrier between the graphene bilayers. The valley degeneracy is g v = 2 and the dielectric constant κ = 3 16 . The separation barrier between the bilayer sheets can be as little as d = 1 nm with no significant leakage from tunneling 17 . Such separations are much less than the effective Bohr radius, ∼ ⁎ a 11 B nm for bilayer graphene. At the densities of interest to us, the average spacing between the carriers in each bilayer is much greater than the lattice constant of graphene, making details of the graphene lattice structure unimportant. We will consider only equal electron and hole densities. Furthermore, high quality graphene lattices with extremely low disorder can be fabricated 18,19 . In this paper we find, besides the superfluid phase discussed in refs 20-23, that there are two new inhomogeneous ground-state phases: (i) a strongly-coupled one-dimensional charge density wave (1D-CDW), and (ii) a coupled electron-hole Wigner crystal (c-WC). Because of the simple structure of the system, the drivers for all these phases are very simple. First, there is the Coulomb attraction between electrons and holes from the two bilayers which is controlled by the layer separation d. Second, there is the Coulomb repulsion between pairs of electrons or pairs of holes within each bilayer which is controlled by the carrier density ρ in the bilayer.
We find that the 1D-CDW, with density modulations in one planar direction, is always more stable than the two-dimensional CDW. (We also find that in regions of phase space where the liquid phase is more stable than the 1D-CDW, that the liquid phase has a lower energy than the two-dimensional CDW). As would be expected, for the 1D-CDW phase, the planar density modulation in the two bilayers share the same periodicity and they are in phase (see Fig. 1(a)), since this ensures a maximum attractive potential energy gain from the electron-hole interactions. In contrast with the case of High-Temperature superconductors, the holes from one bilayer act as a perfectly symmetric polarizable background for the electrons in the other bilayer, and vice versa. This property makes our CDW phase uniquely different from CDW phases in other systems.
In the c-WC phase, an electron Wigner crystal in one bilayer couples to a hole Wigner crystal in the other bilayer. The electron and hole sites of the c-WC will align, again to maximize the attractive potential energy gain from the electron-hole interactions, so each lattice site in the hole Wigner crystal lies directly opposite a lattice site in the electron Wigner crystal (see Fig. 1(b)).
When the separation d between bilayers is not too large, the Coulomb attraction between the electrons and holes will generate strongly bound electron-hole pairs that can condense into a quantum coherent superfluid state 23, 24 . However, for densities above a cut-off density, strong screening is predicted to suppress the superfluidity 23, 25 .
As the layer separation d is increased, an interesting phase transition should occur from the superfluid phase to the 1D-CDW phase. While CDWs are conventionally treated as a classical phase, Bardeen 26 argued for extending the coherent quantum state interpretation from superconductors to include 1D-CDWs. Latyshev et al. 27 subsequently reported observation of Aharonov-Bohm-like oscillations in the original CDW material NbSe 3 28 .
At still lower densities, the transition with increasing d would, instead, be from the superfluid to the c-WC phase. This transition would be driven by the competition between the repulsive interactions between like-charge carriers in one bilayer, which favour the Wigner crystal, and the attractive interactions between electrons and holes, which favour pair formation and superfluidity. At very small densities ρ and for small d, the electron-hole pairs are compact and nearly neutral. The residual repulsive interactions within a bilayer are through weak dipolar interactions, so a superfluid phase would be favoured. However at larger d for the same ρ, the attractive electron-hole interactions between layers become weak compared with the repulsive interactions between like-carriers within each bilayer. The repulsive interactions favour the WC phase, and the system would make a transition from the superfluid to the quantum c-WC phase. At extremely low densities, the Wigner crystal would become classical, and at that point there would be an interesting transition from the quantum coherent superfluid condensate to a classical Wigner crystal phase.
At sufficiently large d, the bilayers will decouple into two independent bilayers, so for decreasing ρ, the Fermi liquid phase will make a transition directly to a Wigner crystal that is only very weakly coupled.

Methods
It is impractical to investigate the CDW phase in this system using quantum Monte Carlo (QMC) techniques, since the amplitude, periodicity and dimensionality, are not a priori known. In a recent work 29 , a new and fast interpolation scheme was introduced to obtain accurate correlation energies for the one-layer 2DEG. The method gives excellent agreement with QMC correlation energies for both one-valley and two-valley two-dimensional-electron-gas systems. In this approach the correlation energy E c [ρ] is determined by an interpolation of W α [ρ], the potential energy functional of ρ without the Hartree contribution, of a fictitious system that interacts with a Coulomb interaction scaled by a coupling constant factor α. (strongly-interacting) 30 , Here, we extend the approach from ref. 29 (see Sec. S1 of Supplementary for a brief discussion of the approach) to coupled electron-hole bilayer graphene sheets. For the two bilayers, the exact Random Phase Approximation for the above expressions in the weakly-interacting limit, α → 0, are given by (3) is determined from the contribution of the second-order correction to the correlation energy E c (2) . It is given by Hartree, independent of the layer separation.
is the non-interacting density-density response function for layer  and is given by, , ,  f k is the Fermi-Dirac function for the wave vector k, g s and g v are the spin and valley degeneracies, η > 0 is an infinitesimal number, and is the wavefunction overlap factor in BLG with φ the angle between k and k′. ε ′ k s s ( ) are the electron (s = +1) and hole (s = −1) energy bands in BLG defined as, which are the eigenenergies of the 2 × 2 Hamiltonian, This can properly describe BLG at low energies 9,31 . γ ≈ 400 meV is the dominant interlayer coupling between the sublattices A and B′ from the upper and lower layers, and v F = 10 6 m/s is the Fermi velocity. U 1 and U 2 are the electrostatic potentials applied to the upper and lower layers which are required to induce electrons (or holes) in each BLG. We note that the transition to the inhomogeneous phases occurs at very low densities, corresponding to low energies, e.g. for < . ×  n 1 5 10 11 cm −2 , <  E 5 F meV. We show in the Supplementary Information that for such low energies, whenever there is a significant gap in the BLG spectrum, , χ 0 can be evaluated by considering only the conduction band E + . For this reason we neglect the influence of the hole band in the calculation of Eq. (5), considering only the conduction band in the BLG with quadratic dispersion.
In the limit of strong interactions, α → ∞, the ground state is the classical WC. We can treat the classical crystal as a collection of neutral unit cells, each cell with an electron or hole at its centre and surrounded by a charged disk of uniform neutralizing background of radius πρ = r 1/ 0 . Then W ∞ is obtained from a straightforward electrostatic calculation,  where ρ  r ( ) is the charge density distribution in layer . We use Eq. (10) to obtain the ground state energy per particle [ ]/ 2 0 for the liquid phase with uniform density ρ(r) = ρ 0 , and for the non-uniform density distribution ρ(r) of the c-WC and the CDW phases.
For the c-WC phase, we take the variational form for the density distribution, For the 2D-CDW we use the same variational form for ρ(r), Eq. (12), but now we take a as an additional variational parameter.
For the 1D-CDW phase, we take the variational form with the amplitude β CDW and the periodicity γ CDW as the two variational parameters. Figure 2 shows the phase diagram as a function of the inter-particle spacing r s within each bilayer sheet and the separation d between the bilayer sheets. We see that for large , the transition to the c-WC phase initially moves to lower density. This is caused by the increasing strength of the interlayer attractions when <  d r 0 . As a result the positions of the electrons and holes become more tightly correlated, and so there is increasing cancellation between their charges. This has the effect of reducing the electron-electron and hole-hole repulsion within each bilayer, which is what drives the transition to the WC. The long-range Coulomb repulsion between like-charges in the same bilayer is screened and tends towards the much weaker dipolar repulsion. The weaker repulsion makes it more difficult to form the c-WC, and the transition moves to lower density.

Results and Discussions
When d decreases below = ⁎ d a / 12 B , a new CDW phase interposes itself between the Fermi liquid and c-WC phases. With increasing r s for fixed d, there are then two transitions: (i) from the Fermi liquid to the CDW phase, and (ii) from the CDW phase to the c-WC phase. As d is further decreased, both these transitions move to higher densities. The reason is that an increase in the strength of the interlayer electron-hole attraction favours formation of the inhomogeneous phases 7,8 .
At sufficiently small d, the electron-hole attraction becomes strong enough for electron-hole bound pairs to form in significant numbers, and these should condense into a coherent superfluid state 6,23 . It has been predicted that there is no BCS superfluid regime because for small < .  r 2 5 s there is very strong screening and this suppresses superfluidity 23 . As r s increases above 2.5, the superfluid phase would first form directly in the BCS-BEC cross-over regime, and for larger r s it would then evolve into the Bose-Einstein regime.
For . < <   r 2 5 1 0 s , the transition temperature for superfluidity at small d~1 nm can be large, i.e. >10 K 23 . However when > ∼ r 20 s the superfluid phase would occur only at extreme low temperatures, since the Kosterlitz-Thouless temperature T KT , which establishes an upper limit on the transition temperature for a 2D superfluid of only T KT ≤ 40 mK by r s = 20 even for the smallest d. We recall that T KT decreases with r s as − r s 2 . In practice, for such low transition temperatures, residual disorder in the graphene lattice is likely to kill any superfluidity. Furthermore, for r s = 20, the predicted superfluid condensate fraction is already extremely small for all values of d. Even at = ⁎ d a / 1 B , the condensate fraction is only c = 0.02 6 , and c decreases even further for larger d. We conclude that for > ∼ r 20 s , even if any superfluid should survive, it would most likely coexist in a phase separated state with the CDW or c-WC phase. In Fig. 2, a possible superfluid phase (not calculated here) has been schematically represented by the transparent area. Figure 3(a-d) show the CDW and c-WC electron bilayer density distributions at r s = 60 for two values of the layer separation d, chosen close to the two CDW-c-WC phase boundaries. It is interesting to note how the periodicity of the CDW changes with d. The periodicity is longer for smaller d, and it evolves towards a WC as d is increased. However even at the upper boundary of the CDW phase region at = ⁎ d a / 10 B , the value is still ∼ a a 2 CDW W C , (Fig. 3(d)). In Fig. 3(a,c), notice that the carriers in the WC become more localized on the lattice sites as d is decreased, due to the fact that for r s = 60 the WC lies closer to the liquid phase boundary when = ⁎ d a / 10 Fig. 2). The differences in the total energies of the Fermi liquid, CDW, and c-WC phases provide us with an upper limit estimate of the transition temperatures. Figure 4(a,b) show the energy differences between the Fermi liquid and c-WC, and the Fermi liquid and CDW, respectively. Both sets of energy differences are consistently of the order of 50-100 K. This remains the case even for the largest layer spacings, where the c-WC approaches the uncoupled Wigner crystal. We thus expect that the transition temperatures for the CDW and c-WC phases to be comparable with the transition temperature for the uncoupled Wigner crystal.

Discussion and Conclusions
The CDW phase could be identified using scanning tunneling microscopy (STM). The formation of the CDW along with the associated periodic lattice distortion opens a gap in the Fermi surface, modifying the local density-of-states, and this could be detected by tunneling [34][35][36][37] . A CDW phase can be also detected by transport measurements. 1D stripes may be pinned by disorder, in which case the CDW could be identified either with standard threshold voltage conductivity measurements or with frequency threshold measurements of ac conductivity 38 . A 1D-CDW breaks the rotational symmetry of the 2D plane, so it could manifest itself as a highly anisotropic transport response. Anisotropic resistivity as evidence of a CDW state has been experimentally studied, for example, in single crystals of R 2 CoGa 8 (R = Gd, Tb, Dy, Ho, Er, Tm, Y, and Lu) 39 and RTe 3 (R = Y, La, Ce, Figure 2. Zero temperature phase diagram for the ground state of the coupled electron-hole bilayer graphene system as a function of layer density parameter r s and the separation between the layers d. The top x-axis gives the corresponding electron (hole) density. The transparent region at small d, schematically represents the superfluid phase (not calculated here) that would interplay with the 1D charge-density-wave (1D-CDW) and coupled Wigner crystal phases.
Pr, Nd, Sm, Gd, Tb, Dy, Ho, Er, and Tm) 40 . However, there would not necessarily be any observable anisotropy, because the lack of an underlying anisotropic crystal structure means that no unique orientation of the 1D-CDW would be favored. Thus formation of randomly oriented domains of 1D-CDWs could occur. Further experimental evidence for the formation of a CDW state would be the observation of negative electronic compressibility in the bilayers. This can be measured via the difference between the actual differential capacitance and the classical geometric capacitance. The negative compressibility has been previously studied in double graphene monolayers in the presence of a perpendicular magnetic field 41 .
STM could also be used to observe a WC density profile. The c-WC phase could be also observed with transport measurements. Wigner crystallization should be accompanied by a transition to an insulating state, caused by pinning of the Wigner lattices by residual disorder 42,43 . The c-WC phase can be distinguished from the CDW phase by observing differences in their low-lying excitation spectrum, since unlike the plasmon of the Fermi liquid phase, the low-energy collective mode for the WC is centred on a momentum transfer q equal to the reciprocal vector of the WC lattice 44 . This characteristic dispersion property of the collective mode should be readily observable in Raman spectra. Experimental investigation of optical phonon modes is an alternative approach to detect a c-WC phase as has been suggested for coupled 2D layers in semiconductor structures 45 .
Next we comment on the puzzling results reported from recent Coulomb drag experiments in coupled electron-hole graphene bilayers 46,47 . A negative Coulomb drag was observed in two different temperature regimes, i.e. at low temperatures down to T = 1.5 K (ref. 46), and at high temperatures up to T = 170 K (ref. 47). Because the drag reported in ref. 47 was observed to be symmetric across the electron-hole and electron-electron systems, our inhomogeneous phases resulting from electron-hole correlations are unlikely to be the origin of the anomalous drag observed in coupled graphene bilayers. In summary, we have proposed a system of two strongly coupled electron-hole bilayer graphene sheets as a promising candidate to observe new inhomogeneous c-WC and 1D-CDW phases which would interplay with the previously predicted electron-hole superfluid 23, 24 . We find in the strong interlayer coupling regime, that a 1D-CDW occurs at significantly higher densities than the c-WC phase.
There has been a long-standing issue in coupled electron-hole systems of whether a CDW phase would be two-dimensional or one-dimensional. A two-dimensional CDW phase would most likely retain the hexagonal structure of the 2D WC phase, but with a longer periodicity. We find that a hexagonal two-dimensional CDW phase has always a larger energy than both the Fermi liquid and 1D-CDW phases. We conclude from this that the hexagonal two-dimensional CDW phase would not be found as the ground state in coupled electron-hole layers.