{'fraction_non_alphanumeric': 0.04010583034878844, 'fraction_numerical': 0.01581002161778466, 'mean_word_length': 4.258568035290126, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Evidence from quantum Monte Carlo of large gap superfluidity and BCS-BEC crossover in double electron-hole layers

We report quantum Monte Carlo evidence of the existence of large gap superfluidity in electron-hole double layers over wide density ranges. The superfluid parameters evolve from normal state to BEC with decreasing density, with the BCS state restricted to a tiny range of densities due to the strong screening of Coulomb interactions, which causes the gap to rapidly become large near the onset of superfluidity. The superfluid properties exhibit similarities to ultracold fermions and iron-based superconductors, suggesting an underlying universal behavior of BCS-BEC crossovers in pairing systems.

We report quantum Monte Carlo evidence of the existence of large gap superfluidity in electron-hole double layers over wide density ranges. The superfluid parameters evolve from normal state to BEC with decreasing density, with the BCS state restricted to a tiny range of densities due to the strong screening of Coulomb interactions, which causes the gap to rapidly become large near the onset of superfluidity. The superfluid properties exhibit similarities to ultracold fermions and iron-based superconductors, suggesting an underlying universal behavior of BCS-BEC crossovers in pairing systems. There are intense ongoing experimental efforts to observe superfluidity in electron-hole double layer systems, including double quantum wells in GaAs-AlGaAs heterostructures, double graphene monolayers, double graphene bilayers, and hybrid graphene-GaAs structures [1][2][3][4][5]. A recent very significant experimental advance by multiple groups has been the fabrication of closely spaced electron-hole double graphene bilayers with carrier densities tunable by metal gates [6][7][8][9]. Insertion of a few layers of hexagonal boron nitride between the bilayers creates an insulating barrier. This can be as thin as 1 nm while still blocking tunneling of carriers between the bilayers, permitting the electrons and holes to interact with a very strong Coulomb attraction.
Condensation of electron-hole pairs into a BEC superfluid state in double graphene bilayers has been experimentally demonstrated in the quantum Hall regime with a magnetic field [8,9], opening the way to the generation of quantum coherent macroscopic states in spatially separated twodimensional sheets.
Room temperature superfluidity had in fact been predicted earlier in electron-hole double graphene monolayers [10], but it was later established that since a graphene monolayer remains blocked in the weakly-coupled regime, strong screening always suppresses the superfluidity [11]. Double electronhole graphene bilayers were subsequently proposed to overcome this problem [12], combining ideas for the realization of high-T c superfluidity with the ability to move across the BCS-BEC crossover by changing the carrier densities using metal gates. This ability to tune the system into the strong-coupling regime is key to obtaining a superfluid in a solid-state device at experimentally accessible conditions [12], offering a fascinating alternative to ultracold fermionic atoms for studying superfluid physics across the BCS-BEC crossover.
With superfluidity in cold atoms, quantum Monte Carlo (QMC) results [13][14][15] were from the outset integrated closely with experiments and theory to understand and control the phenomenon. QMC simulations are extremely useful in any strongly correlated condensed matter system where there is no small parameter that can be used in perturbative expansions or controllable diagrammatic approximations. For double layer electron-hole systems, QMC simulations [16,17] are of key importance since the superfluid phase arises from a complicated competition between the long-ranged Coulomb interlayer attraction and intralayer repulsion. Our work provides a comprehensive characterization of the superfluid properties of electron-hole double layers at zero temperature, and establishes the density range in which the condensate exists. Our predictions that too high a density kills superfluidity and of the existence of a BCS-BEC crossover regime at accessible intermediate densities will be extremely useful in guiding experiments. We use state-of-the-art QMC methods to determine the dependence of BCS-like effective parameters on the equal carrier densities, with the distance separating the layers fixed at the value that maximizes exciton formation [17][18][19].
The boundary between the crossover and BEC regimes is of great interest because the pseudogap transforms into a real gap associated with the molecular phase near this boundary [20]. The BCS-BEC crossover has been recently observed in the shallow Fermi surface pockets of iron-based superconductors [21,22].
The zero temperature BCS-BEC crossover [23,24] can be traced by following the evolution of the condensate fraction c, the average radius of the superfluid pairs r ex [25], or the superfluid gap ∆. Here we use explicit many-body wave-functionbased methods to track the evolution of the superfluid with carrier density in symmetric electron-hole layers. We directly obtain total energies, condensate fractions, and pair radii, and we extract superfluid gaps and chemical potentials from total energies using BCS-like relations with effective parameters. We employ methods similar to those applied in Ref. 26 for coupled layers of dipolar fermions.
The condensate fraction c measures the fraction of carriers in the condensate [27][28][29][30]. In the BCS regime c < 0.2, with only a small fraction of the carriers near the Fermi surface forming the condensate, while in the BEC regime c > 0.8, so almost all the carriers form local molecular bosonic pairs, and condense. Using the pair radius r ex , the BCS regime is characterized by k F r ex ≫ 1, where k F is the Fermi wave vector, while in the BEC regime the compact pairs correspond to k F r ex < 1. The parameter k F r ex was the first to be studied in investigations of the BCS-BEC crossover in the high-T c cuprate superconductors [31]. r ex determines the correlation length of the pairs, which enters the expressions for the properties of the vortex state of the superfluid and all other quantities affected by the spatial structure of the superfluid wave-function, such as Josephson and Andreev reflection effects [32].
In contrast to c and r ex , the superfluid gap is experimentally straightforward to measure using angle-resolved photoemission spectroscopy (ARPES), scanning tunnelling microscopy, or measurements of the specific heat. Knowledge of the evolution of the superfluid gap with the external parameters is of fundamental relevance in designing experiments to detect and characterize electron-hole superfluidity. It is possible to link the entry into the BCS-BEC crossover regime as determined by c and k F r ex to ∆/E F ∼ 1 at zero temperature [33]. In our calculations we also monitor the evolution of the pseudo-Luttinger wave vector k min at which the k-dependent excitation energy passes through its minimum. k min can be traced by ARPES, because it directly affects the shape of the remnant Fermi surface in the broken symmetry phase at zero temperature.
With our results for ∆, µ, c, and r ex as functions of density, we are then in a position to follow the evolution of the system through the weak-coupling regime, the superfluid BCS-BEC crossover regime, and the BEC regime, enabling comparisons with predictions from various microscopic theories.
In our calculations we use excitonic Hartree units, = |e| = m * e = κ4πǫ 0 = 1, where m * e is the effective electron mass and κ is the relative permittivity of the system, and we obtain energies in units of Ha * = (m * e /m e )κ −2 Ha and distances in units of a * 0 = κ(m e /m * e ) a 0 . For reference, the relative permittivity for bilayer graphene (BLG) encapsulated in few-layer hexagonal boron nitride (hBN) is κ = 2 [34], and the effective mass is m * e = 0.04m e [35]. We simulate a finite version of the paramagnetic, equalmass electron-hole double layer system with parabolic singleparticle energy dispersion using square simulation cells of area A subject to periodic boundary conditions, with N particles in each layer. The in-layer particle density is defined via the density parameter r s = A/(πN ). We run calculations for systems containing N = 58 electron-hole pairs; tests with systems of N = 114 electron-hole pairs show that finite-size errors are small in our results [36]. We use a fixed interlayer separation of d/a * 0 = 0.4, slightly greater than the largest d at which biexciton formation is favorable [17][18][19], and vary the density between r s /a * 0 = 1.75 and 15. At all densities considered we evaluate total energies, condensate fractions, and pair-correlation functions (PCFs) using the variational quantum Monte Carlo (VMC) method [37,38]. Wave function parameters can be optimized within VMC [39,40], and the accuracy of VMC expectation val-ues depends on the quality of the resulting wave function. The more computationally costly diffusion quantum Monte Carlo (DMC) method [38,41] employs stochastic projection to extract the lowest-energy state compatible with the nodal surface of a VMC-optimized trial wave function. Once the time-step and population-control biases are eliminated [42][43][44], the accuracy of the DMC method depends only on the accuracy of the nodes of the trial wave function. We have performed DMC calculations at selected densities representative of the weak-coupling, crossover, and strong-coupling regimes. Thus, our DMC calculations serve as quantitative corrections to our VMC results throughout the density range considered.
We use trial wave functions of the form Ψ = e J D e↑h↓ D e↓h↑ , where e J is a Jastrow correlation factor [45,46], imposing the Kato cusp conditions [47], and D e↑h↓ and D e↓h↑ are pairing determinants [16,17]. In systems with an additional (up-spin) electron, we complete the corresponding determinant with a plane-wave orbital of wave vector k. Details of our trial wave functions are given in the Supplemental Material [36]. We use the CASINO code for our calculations [48].
The main properties of the BCS-BEC crossover at low temperature are captured by BCS theory, as demonstrated with ultracold fermions. When the gap equation is coupled to the density equation, the excitation energy ε(k) corresponding to the addition of an electron of wave vector k follows the BCS dispersion relation [26], where µ, m * , and ∆ are the chemical potential, effective mass, and superfluid gap of the electron quasiparticle, respectively. ε(k) therefore contains the parameters that characterize the superfluid state. This excitation energy can be obtained from ab initio total energies as where E A (N ) is the energy of a system of area A containing N electron-hole pairs, E A (N + 1/2; k) is the energy of a system of area A containing N electron-hole pairs and one additional electron of associated wave vector k, and µ QMC (N ) is the chemical potential of the system. Note that µ QMC differs from the mean-field µ due to many-body effects. We simulate systems with N , N + 1, and N − 1 electronhole pairs, and systems with N electron-hole pairs and an unpaired electron at several wave vectors k. We then compute ε(k) at each k using Eq. 2 and fit the resulting values to Eq. 1 with m * , µ, and ∆ as fitting parameters. Following tests [36], we use wave vectors such that 0 ≤ k < k cut ; we set k cut /k F = 1.5 at high densities and use larger cut-off values at low densities. Although the superfluid gap is expected to be k-dependent [11,12], our calculations do not yield any significant variation of ∆ with k in the ranges of k we have considered [36]. Figure 1 shows a plot of the DMC values of ε(k) for densities of r s /a * 0 = 2, 5, and 10, along with the resulting fits to Eq. 1. The fits follow the DMC data remarkably well, indicating that the BCS dispersion relation provides a robust description of the ab initio results. In Fig. 2 the zero-temperature superfluid gap is reported as a function of r s , both in excitonic Hartree units and relative to E F . At high densities there is no superfluidity because the electron-hole pairing interaction is strongly screened [12]. Near the onset density, the electron-hole condensate is already close to the BCS-BEC crossover boundary. As r s /a * increases above 2 there is a very steep increase in the gap, which exceeds ∆/E F ∼ 1 by r s /a * 0 ∼ 3. ∆/E F > 1 signals entry into the BCS-BEC crossover regime, so this occurs practically immediately after the onset of superfluidity. The steep rise in ∆ is associated with strong screening at high densities [12]. Consequently the weakly coupled BCS superfluidity regime, for which ∆/E F ≪ 1, exists at most in a tiny range of densities. Figure 3 (a) shows µ as a function of r s . µ becomes negative by r s /a * 0 ∼ 6, which signals entry into the BEC regime. In Figs. 3 (b) and 3 (c) we plot the location of the minimum of ε(k), k min = arg min k ε(k), and m * as functions of r s . The value of k min /k F tracks the collapse of the Fermi surface, going from unity in the weak-coupling regime to zero in the BEC regime. The DMC results suggest that the Fermi surface fully collapses at a somewhat lower density than predicted by VMC. In Fig. 3 (c) the quasiparticle mass m * has a minimum of less than the effective electron mass m * e near where the superfluid gap is maximal, r s /a * 0 ∼ 5. This is indicative of the interplay between the intralayer repulsion and the interlayer attraction, leading to quasiparticles in the superfluid state with masses m * < m * e for intermediate r s . This behavior of m * differs from theoretical and experimental findings in ultracold fermions in two dimensions, where the interaction is purely attractive. There, m * m * e always, and it varies monotonically with r s . The quasiparticle mass of the two-dimensional electron gas is extensively discussed in Ref. 49. Experimental measurements indicate a regime of small r s /a * 0 < 3 in which the quasiparticle mass is smaller than the effective electron mass, see Fig. 4 of Ref. 49. Thus, competition between intralayer repulsion and interlayer attraction in the electron-hole double layer can lead to small m * < m * e , as we find. The boundaries between the BCS, BCS-BEC crossover, and BEC regimes can be determined from the condensate fraction c, which is defined as where ρ (2) eh (r) is the translational-rotational average of the two-body density matrix for electron-hole pairs [36]. We have evaluated c using the estimator of Ref. 50 which removes onebody contributions to ease extrapolation to the r → ∞ limit. The results for c shown in Fig. 4 (a) are consistent with the conclusions drawn from the behavior of µ. The condensate fraction is negligible for r s /a * 0 1.5 [17]. As r s /a * 0 increases, c grows rapidly to ∼ 0.2 by r s /a * 0 = 2, signaling entry into the BCS-BEC crossover regime. As r s /a * 0 is further increased, c increases substantially and by r s /a * 0 = 8 it exceeds c = 0.8, thus entering the BEC regime.
We also compute the translational-rotational average of the electron-hole PCF g eh (r) [36], which allows us to evaluate the exciton radius r ex as where r 1 is the radius of a circle centred on a hole containing on average one electron. Figure 4 (b) shows the pair radius r ex for the condensate. As r s increases, r ex converges to the isolated-exciton limit, which for d/a * 0 = 0.4 is r ex /a * 0 = 2.002. The values of k F r ex are always less than unity, indicating pair sizes of the order of, or smaller than, the interparticle distance, confirming that the superfluidity is always in the strongly coupled crossover or BEC regimes.
In contrast with the short-range interactions typical of ultracold fermions and superconductivity in general, electron-hole superfluidity should be affected by screening because of the long-range nature of the Coulomb pairing attraction. The nature and effectiveness of this screening has been a source of controversy in the past, with predictions from mean-field calculations ranging from negligible screening resulting in roomtemperature superfluidity [10,11,51,52], to full screening by the normal state which would essentially completely suppress superfluidity [53]. Reference 54 compared the dramatically different mean-field predictions for the density dependence of the condensate fraction with the QMC values of the condensate fraction for the same system [17]. The conclusion was that the best mean-field approximation for screening was selfconsistent screening by the superfluid state introduced by Lozovik [11] and applied in Ref. 12 to double bilayer graphene. Using a similar argument, Ref. 54 concluded that the vertex corrections and intralayer correlations in the superfluid state make relatively small contributions. Our present results for the superfluid gap, Fig. 2 (a), are in good agreement with The trends of our QMC results are consistent with the complete suppression of superfluidity at high densities predicted in Ref. [12]. After the onset of superfluidity, r s /a * 0 = 1.5, the system is in the weak-coupled BCS regime, but the condensate fraction rises rapidly and by r s /a * 0 = 2 it reaches c = 0.2. By r s /a * 0 = 2 the superfluid gap has reached ∼ E F . Thus the BCS regime, for which ∆/E F ≪ 1, is restricted to the very small density range 1.5 < r s /a * 0 < 2. This confirms the effects of the highly non-trivial competition between Coulomb screening, which tends to suppress electron-hole pairing induced by the interlayer Coulomb attraction, and the opening of a large superfluid gap, which suppresses the particle-hole processes near the Fermi surface responsible for screening, thus severely weakening the screening.
When r s is further increased, the superfluid gap first increases and then reaches a flat maximum around r s /a * 0 ∼ 5 -6 with a very large value ∆/Ha * ∼ 0.2. The large gap indicates that the superfluidity is robust with a high transition temperature. The chemical potential is still positive at this density, µ/E F ∼ 0.5, the condensate fraction is ∼ 0.7, and k min /k F ∼ 0.5. Thus the system retains its fermionic properties with a Fermi surface intact but smeared out by the large gap, ∆/E F ∼ 7.
The BEC superfluid regime is reached at larger r s /a * 0 ∼ 8, where the condensate fraction acquires values c > 0.8. In this regime µ/E F is large and negative, the ratio ∆/E F > 10 is very large, the Fermi surface has completely collapsed, and the average pair size approaches the radius of an isolated exciton. The electron-hole superfluid can then be regarded as an ensemble of well-formed electron-hole dipoles, which are indirect excitons. The excitons will behave as a twodimensional bosonic gas with a repulsive interaction, with a Kosterlitz-Thouless transition [55] governing the critical temperature for superfluidity. Thus in the BEC regime the critical temperature should diminish with decreasing density.
We can thus conclude that the optimal density for experimental realization of the electron-hole condensate is around r s /a * 0 ∼ 5, which is deep inside the BCS-BEC crossover regime with large values of ∆/E F . For BLG encapsulated in hBN, r s /a * 0 ∼ 5 corresponds to a density of 2 × 10 11 cm −2 , the maximum gap corresponds to ∆ ∼ 54 meV ∼ 630 K, and the interlayer distance d/a * 0 = 0.4 corresponds to 1 nm. Our QMC results are consistent with a universal behavior of materials in the BCS, BCS-BEC, and BEC regimes not depending on the details of the microscopic interactions, and thus they point to a very general physics. The ground state properties and their evolution with coupling strength appear to be universal for (i) long-range Coulomb interactions, (ii) contact interactions in fermions [20], and (iii) spin fluctuations and phonons in iron-based superconductors [21,22,33].
Our results confirm that, unlike for fermionic superfluids with short-range pairing interactions, the BCS regime in Coulomb systems with their long-range interactions and screening is restricted to a very small range of densities. This is due to competition between screening and the superfluid gap [11,12], with strong screening suppressing the small-gap BCS regime in Coulomb systems. At high densities, the onset of superfluidity is delayed by screening, so that when the onset density is eventually reached, the pairs are relatively compact, and the superfluid gap, which rapidly becomes large both in absolute value and relative to E F , will strongly suppress screening. Thus the system almost immediately enters the strong-coupling BCS-BEC crossover regime. For this reason the superfluidity is likely to be robust against potential detrimental effects like disorder, density imbalance, and low dimensional fluctuations, and we expect the largest gaps and highest critical temperatures not to be far from the onset density.
e i and h j are the in-plane position vectors of the ith electron and the jth hole and d is the distance between the layers. Note that, in terms of the in-layer density parameter, r s , the magnitude of the Fermi wave vector is k F = √ 2/r s and the Fermi energy is E F = 1/r 2 s . We simulate systems consisting of a finite number of electrons and holes subject to periodic boundary conditions, thus the Coulomb interaction of Eq. S1 must be replaced with the two-dimensional version of Ewald summations [S2].
Each of the pairing determinants in our trial wave function is of the form where e σ i and h σ i are the position vectors of the ith σ-spin electron and hole, respectively, and φ(r) is a pairing orbital containing optimizable parameters of the form used in Ref. S1, p l cos(k l · r) where k l is the lth shortest reciprocal lattice vector [excluding one of each (k, −k) pair due to symmetry] and {p l }, {c m }, and L are optimizable parameters, with p l = p l ′ when |k l | = |k l ′ |. The pairing orbital is constrained to be cuspless at r = 0. There are 22 optimizable parameters in the pairing orbitals (one of the coefficients is fixed due to normalization, and another is determined by the cusplessness condition). At r s /a * 0 ≥ 8 we find that the cosine expansion does not provide useful variational freedom to the wave function and we set p l = 0, in which case the pairing orbitals effectively contain 8 optimizable parameters.
We also run calculations for systems with an additional (upspin) electron for which we replace the first pairing determi-nant with where k is a wave vector commensurate with the reciprocal lattice of the simulation cell associated with the additional electron.
We use an isotropic two-body Jastrow factor [S3] of the Drummond-Towler-Needs parametrization [S4], where r ij is the distance between particles i and j, Θ is the Heaviside step function, and {α m } and L u are optimizable parameters. Different parameter values are used for samespin same-particle pairs, opposite-spin same-particle pairs and electron-hole pairs, and the two-dimensional Kato cusp conditions are applied for same-particle pairs [S5]. The electronhole u function is constrained to be cuspless, and we do not use the "quasi-cusp" electron-hole term decribed in Ref. S1. There are 27 optimizable parameters in the Jastrow factor (the cusp/cusplessness conditions determine one of the coefficients for each particle-pair type).

EVALUATION OF EXPECTATION VALUES
The translational-rotational average of the one-body density matrix for an electron (or hole) is Similarly, the translational-rotational average of the two-body density matrix for electron-hole pairs is (S7) The evaluation of ρ (1) e (r) and ρ (2) eh (r) in QMC is performed by computing the corresponding wave function ratio at random values of r and accumulating the resulting values in bins. We use these density matrices to compute the function [S6] from which the condensate fraction can be obtained as c = lim r→∞ c(r). In Fig. S1 we present VMC plots of c(r) at the densities considered, where the plateau at large r that determines c can be seen. The slight increase of the value of c(r) at large r for r s /a * 0 = 1.75 is typical of the weak-pairing regime. The translational-rotational average of the electron-hole pair-correlation function (PCF) is (S9) We use the PCF to evaluate the average size of the exciton pair. In Fig. S2 we present VMC plots of the radial charge density of electrons around a hole, g eh (r)r −2 s 2πr, at the densities considered, as well as the radial charge density of an isolated exciton.

DETAILS OF THE DIFFUSION MONTE CARLO CALCULATIONS
We have performed DMC calculations at selected densities of r s /a * 0 = 2, 5, and 10. Each DMC total energy is obtained from a DMC calculation consisting of M 1 steps with a time step of τ 1 and a target walker population of P 1 , and a second DMC calculation consisting of M 2 = M 1 /2 steps with a time step of τ 2 = 4τ 1 and a target walker population of P 2 = P 1 /4. These ratios give the most efficient extrapolation of the results to the zero-time-step, infinite-population limit for a fixed computational cost [S7, S8].
The DMC time steps must be set according to the smallest length scale of the system λ. At high densities λ is of the order of r s , while at low densities λ is of the order of the exciton radius. Thus, we approximate λ = min(r s , 3a * 0 ) and set τ 1 = 0.01λ 2 . We use P 1 = 2048 walkers and adjust M 1 to obtain the desired statistical accuracy.
We compare the VMC and DMC excitation energies for r s /a * 0 = 2 as a function of k in Fig. S3. The absolute differences between the VMC and DMC excitation energies are of the order of 0.02 Ha * , with the shape and basic features of both curves being largely identical.

FITTING PROCEDURE
We evaluate the excitation energy using Eq. 2 in the manuscript at values of k commensurate with the reciprocal lattice of the simulation cell in the range 0 ≤ k < k cut , and fit the results to the BCS dispersion relation, Eq. 1 in the manuscript. The standard errors in the fit parameters are determined by a stochastic process where ε(k) at each k is replaced by a random number drawn from a normal distribution centred at the mean value of ε(k) of variance its standard error, and the fit is carried out. We obtain the standard error in each fit parameter as the square root of the variance of the values of the parameter in 10,000 realizations of this process.
While the quality of a least-squares fit increases with the number of data points available, we find that the fitting pa-rameters depend strongly on k cut at small r s . In particular, parameters µ and m * are intrinsic properties of the electron quasiparticle and are expected to be k-independent, but are found to vary significantly with k cut . This dependence indicates that Eq. 1 models the ab initio results inconsistently over different k ranges, and it is therefore critical to restrict k cut so that the values of µ and m * are constant over the whole range 0 ≤ k < k cut .
We determine k cut at each value of r s by analyzing subsets of VMC data in different k ranges. From the full set of k values at which we have obtained ε(k), k 1 < k 2 < . . . < k n k , we construct a subset of n consecutive data points from k i = k L to k i+n−1 = k R , and perform a fit to Eq. 1. We then compare the resulting fit parameters for different k ranges, and choose k cut to be the largest k R that yields parameters µ and m * consistent with those from ranges with smaller k R .
In Fig. S4 we demonstrate this procedure for r s /a * 0 = 2 using n = 11 points in each subset. Also shown are fits using smaller windows containing the first 6, 7, 8, 9, and 10 data points to visualize the low-k behavior of the fit parameters. For all densities in the range 1.75 ≤ r s /a * 0 ≤ 6 we find that µ and m * remain constant up to k cut /k F = 1.5. For r s /a * 0 = 8 we set k cut /k F = 2, and for larger r s we do not find any significant variation of the fit parameters with k cut , and we use n k = 40 and 60 data points for r s /a * 0 = 10 and 15, respectively, corresponding to k cut /k F = 3.11 and 3.97. We apply the values of k cut obtained from this VMC assessment to our DMC calculations, with the only difference that in our DMC runs for r s /a * 0 = 10 we skip values of k so as to require 20 data points instead of 40, due to computational cost considerations.
The superfluid gap ∆ is predicted to be a function of k peaked at k 0 of half-width w and a 1/k tail at large k [S9, S10]. The centre is expected to shift from k 0 ∼ k F in the BCS regime to k 0 = 0 in the BEC regime, and the halfwidth is expected to increase with r s . However, due to the small k region over which we are able to fit our data to Eq. 1, we do not detect any significant variation of ∆ with k, and we can only infer that its half-width is of at least the order of k F at high densities. We note that the introduction of a k-dependent ∆ of the aforementioned properties does not account for the variation of µ and m * with k beyond k cut .
In Tables S1 and S2 we give the parameters obtained using this fitting procedure with VMC and DMC data, respectively, as plotted in the manuscript. The value of the effective mass parameter appears to tend asymptotically to m * e as r s increases. At r s /a * 0 = 15, since we obtain a VMC effective mass within uncertainty of m * /m * e = 1, we perform a second fit with m * /m * e fixed at unity to improve the estimation of µ and ∆, and we report the parameters obtained from this fit as our final results for this density.  TABLE S2. Superfluid parameters ∆, µ, and m * obtained by fitting DMC values of ε(k) to Eq. 1, along with DMC estimates of the condensate fraction c and pair size rex.

FINITE-SIZE TESTS
The results reported in the main manuscript have been obtained by simulating a 58-exciton system. In this section we quantify the bias incurred by this choice of system size. We expect this bias to be small, since the superfluid parameters are obtained via ε(k), which is a difference of energy differences, ε(k) = E A (N + 1/2; k) − E A (N ) − 1 4 [E A (N + 1) − E A (N − 1)], potentially enabling substantial error cancellations in the final results.
We have performed VMC tests using a simulation cell containing N = 114 electron-hole pairs at selected densities of r s /a * 0 = 2, 5, and 10 to quantify the finite-size error in our results at N = 58. To illustrate these results we plot the N = 114 VMC excitation energies for r s /a * 0 = 2 as a function of k in Fig. S5. The VMC values of ε(k) at both system sizes appear to lie on the same curve as a function of k/k F . The fit parameters obtained from the N = 114 VMC data are tabulated in Table S3 and plotted in Fig. S6. The superfluid parameters obtained for N = 114 are nearly identical to those at N = 58, and we conclude that the finite-size errors of our results are negligible.