††thanks: These authors contributed equally††thanks: These authors contributed
equally

# Geometric Stiffness in Interlayer Exciton Condensates

Nishchhal Verma  Department of Physics, Columbia University, New York, NY
10027, USA  Center for Computational Quantum Physics, Flatiron Institute, New
York, New York 10010, USA  Daniele Guerci  Center for Computational Quantum
Physics, Flatiron Institute, New York, New York 10010, USA  Raquel Queiroz
<EMAIL_ADDRESS>Department of Physics, Columbia University, New
York, NY 10027, USA  Center for Computational Quantum Physics, Flatiron
Institute, New York, New York 10010, USA

###### Abstract

Recent experiments have confirmed the presence of interlayer excitons in the
ground state of transition metal dichalcogenide (TMD) bilayers. The interlayer
excitons are expected to show remarkable transport properties when they
undergo Bose condensation. In this work, we demonstrate that quantum geometry
of Bloch wavefunctions plays an important role in the phase stiffness of the
Interlayer Exciton Condensate (IEC). Notably, we identify a geometric
contribution that amplifies the stiffness, leading to the formation of a
robust condensate with an increased BKT temperature. Our results have direct
implications for the ongoing experimental efforts on interlayer excitons in
materials that have non-trivial quantum geometry. We provide quantitative
estimates for the geometric contribution in TMD bilayers through a realistic
continuum model with gated Coulomb interaction, and find that the
substantially increased stiffness allows for an IEC to be realized at amenable
experimental conditions.

##### Introduction.—

Advancements in topological quantum matter have drawn attention to the crucial
role of Bloch wavefunctions in diverse condensed matter systems. While the
influence of Berry curvature on non-interacting electrons is well-understood
as an anomalous velocity [1], a closely related quantity, the quantum metric,
has recently gained significant attention particularly in the context of flat-
band superconductivity and related experiments in moiré heterostructures [2].
Derived from the geometric properties of Bloch wavefunctions, the quantum
metric has profound effects on various facets of superconductivity. Notably,
it modifies the mass of Cooper pairs [3, 4, 5], phase stiffness [6, 7, 8, 9,
10, 11, 12], spectral weight [13, 14, 15, 16], and potentially the critical
temperature [17, 18]. Other than superconductivity, quantum geometry is known
to appear in current noise spectrum [19], dielectric response [20], electron-
phonon coupling [21], plasmons [22] and nonlinear response [23, 24, 25, 26,
27]. In this work, we target Interlayer Exciton Condensates (IECs) and reveal
a significant geometric contribution to the phase stiffness that results in a
more robust condensate characterized by a higher Berezinskii-Kosterlitz-
Thouless (BKT) transition temperature.

Excitons, bound states of electrons and holes in semiconductors, are bosons
that have long been proposed to form a Bose-Einstein condensate (BEC) at low
temperatures [28, 29]. Unlike conventional BECs, exciton condensates conserve
total particle number and instead break a $U(1)_{e}\times U(1)_{h}$ symmetry
that corresponds to separate conservation of electrons and holes [30, 31].
This symmetry is experimentally realizable in a bilayer system with a spacer
that suppresses single-particle tunneling between the layers. If the electrons
reside on the top $(t)$ and holes in bottom $(b)$ layers, IEC is formed by the
spontaneous breaking of $U(1)_{t}\times U(1)_{b}$ symmetry [32, 33, 34, 35].
Its superfluid properties have been observed in quantum Hall systems [36, 37].

Although an exciton condensate arising intrinsically in a real material has
been a challenge, there has been progress in three-dimensional semimetal
1T-${\rm TiSe}_{2}$ [38],  monolayer ${\rm WTe}_{2}$ [39, 40], bilayer ${\rm
WSe}_{2}$ [41], and TMD bilayer ${\rm WSe}_{2}/{\rm MoSe}_{2}$ [42]. In
particular, ref. [43] has established the existence of interlayer excitons in
the ground state by capacitance measurements and characterized the exciton
Mott transition [44, 45, 46, 47, 48] as a function of the density of electron-
hole pairs. These interlayer excitons have finite dipole moment and interact
via dipole-dipole interaction. Thus, as a virtue of interacting bosons in 2D,
it is possible that there is condensation at low enough temperatures [49, 50,
51, 52, 53]. Confirming the existence of the condensate requires transport
experiments [31, 54], which have remained elusive  until two recent reports
[55, 56].

IECs display fascinating transport properties, most notably dissipationless
counterflow transport. When equal and opposite fields are applied to the two
layers, excitons flow without resistance in the condensate [36]. The
longitudinal counterflow conductivity, $\sigma_{\rm CF}(\omega)$, diverges in
the dc limit with $\sigma_{\rm CF}(\omega)=D_{s}\delta(\omega)+\cdots$, where
the weight of the delta function, $D_{s}$, represents the phase stiffness that
governs the free energy cost of phase fluctuations in the condensate [57]. In
this study, we demonstrate that $D_{s}$, in addition to a conventional
contribution from band dispersion, has a significant geometric contribution
that arises from the wavefunctions of the non-interacting electron and hole
bands.

The critical role of wavefunctions in exciton condensation is evident from
quantum Hall bilayers where excitons exhibit macroscopic coherence and
condense [58, 59] despite the non-interacting bands having flat dispersion and
infinite mass. The non-trivial wavefunctions of the Landau level endow
mobility to excitons even when the constituent electrons and holes are
immobile. Similar phenomenology is anticipated in materials without a magnetic
field, such as ${\rm WSe}_{2}/{\rm MoSe}_{2}$, as a virtue of their non-
trivial wavefunctions [60, 61, 62]. However, there are key distinctions due to
finite dispersion and asymmetric electron-hole bands, and it remains to be
seen if such effects preserve or diminish the geometric contribution. We
address these fundamental questions in this paper.

We begin by deriving the phase stiffness of a general exciton Hamiltonian,
separating it into conventional and geometric components. The latter arises
primarily from the non-interacting electron/hole wavefunctions.  To
demonstrate this effect, we consider a realistic continuum model endowed with
screened Coulomb interactions applicable to bilayer TMD devices and provide
estimates for the geometric contribution to stiffness in ${\rm MoTe}_{2}$
homobilayers [63] as well as ${\rm WSe}_{2}/{\rm MoSe}_{2}$ heterobilayers
[43]. Finally, we propose experimental setups that can serve as validation
platforms for our theory.

##### Phase Stiffness.—

Thermodynamic stability of a condensate depends on the free energy cost
associated with spatial fluctuations in the phase,
$\theta\rightarrow\theta({\bf r})$, which is quantified as
$\mathcal{F}=(1/2)D_{s}\int d{\bf r}|\boldsymbol{\nabla}\theta({\bf r})|^{2}$,
where $D_{s}$ represents the phase stiffness. For IECs, the phase stiffness is
calculated as a linear response coefficient when equal and opposite vector
potentials are applied to the two layers [37]. The antisymmetric vector
potential couples symmetrically to the exciton, since electron and holes have
opposite charges, inducing infinitesimal phase fluctuations in the condensate.

A general model Hamiltonian for the IEC is given by $\mathcal{H}_{\rm
ex}=\mathcal{H}_{0}+\mathcal{H}_{\rm int}$, where
$\mathcal{H}_{0}=\mathcal{H}_{t}+\mathcal{H}_{b}$ describes the non-
interacting properties with Bloch Hamiltonians $\\{\mathcal{H}_{\nu}\\}$ for
the two layers $\nu=\\{t,b\\}$, and $\mathcal{H}_{\rm int}$ is a density-
density interaction that gives rise to excitons. We write
$\mathcal{H}_{0}=\sum_{\bf k}\Psi^{\dagger}_{{\bf k}}H_{0}({\bf
k})\Psi^{\dagger}_{{\bf k}}$ where $H_{0}({\bf k})=[H_{t}({\bf k})\oplus
H_{b}({\bf k})]+V_{b}\tau_{z}$ and $\Psi^{\phantom{{\dagger}}}_{{\bf k}}$ is a
spinor that has internal labels $(\alpha,\beta)$ that are omitted for brevity
along with layer label $\nu$: $\Psi^{\phantom{{\dagger}}}_{{\bf
k}}=(c^{\phantom{\dagger}}_{{\bf
k},\alpha,t},\cdots,c^{\phantom{\dagger}}_{{\bf k},\beta,b},\cdots)^{T}$, and
${\bf k}\in{\rm BZ}$. The last term describes a bias voltage $V_{b}$ that
tunes the gap between the conduction and valence bands (see Fig. 1)  with
$\tau_{z}$ being the Pauli matrix in the layer subspace. With vanishing
single-particle tunneling between the layers, the non-interacting model has a
$U(1)_{t}\times U(1)_{b}$ symmetry that is spontaneously broken by the IEC.

Figure 1:  Device setup with the two layers, orange and blue, with electrons
in one and holes in the other. The hBN spacer (green) suppresses single-
particle tunneling between the layers. Bias Voltage, $V_{b}$, tunes the gap
between lowest unoccupied band and highest occupied band. It is important that
these are direct excitons. The bands are shifted horizontally to highlight the
layers.

Density-density interactions remain unaffected by the external vector
potential. Consequently, when equal and opposite ${\bf A}=A\hat{x}$ are
applied to the two layers, the current operators depend solely on the Bloch
Hamiltonian as $H_{0}({\bf k},A)=H_{t}({\bf k}-eA/\hbar)\oplus H_{b}({\bf
k}+eA/\hbar)$. We assume spatially isotropic systems to suppress the tensor
nature of currents and stiffness and leave the extension to anisotropic
systems to App. A6. The current operator can be expressed as
$j=-\delta\mathcal{H}_{0}/\delta A=j_{P}+Aj_{D}$, where
$j_{P}=(e/\hbar)\sum_{{\bf k}}\Psi^{\dagger}_{{\bf
k}}\left[\tau^{z}\partial_{k}\mathcal{H}_{0}({\bf
k})\right]\Psi^{\phantom{\dagger}}_{{\bf k}}$ and
$j_{D}=-(e/\hbar)^{2}\sum_{{\bf k}}\Psi^{\dagger}_{{\bf
k}}\partial_{k}^{2}\mathcal{H}_{0}({\bf k})\Psi^{\phantom{\dagger}}_{{\bf k}}$
are the paramagnetic and diamagnetic currents respectively. The stiffness is
determined by the Kubo formula, given by $D_{s}=-[\langle
j_{D}\rangle-\chi_{j_{P}j_{P}}({\bf q}_{\perp}\rightarrow 0,\omega=0)]/4$,
where $\chi_{j_{P}j_{P}}$ is the longitudinal current-current correlator [57,
64, 65]. Calculating the stiffness requires a complete enumeration of the
eigenstates of the full interacting Hamiltonian $\mathcal{H}_{\rm ex}$. To
make progress beyond the formal definition, we focus on mean-field theories
where the interaction term breaks down into fermionic bilinears. We can then
utilize the eigenspectrum $E_{m,{\bf k}}$ and eigenstates $|u_{m,{\bf
k}}\rangle$ of the mean-field Hamiltonian to express the stiffness as:

$D_{s}=\dfrac{\hbar^{2}}{4e^{2}A}\Big{[}\sum\limits_{{\bf k},m}f_{m,{\bf
k}}\left\langle u_{m,{\bf k}}\left|\partial^{2}_{k}H_{0}({\bf
k})\right|u_{m,{\bf k}}\right\rangle-\\\ \sum\limits_{{\bf
k},m,n}\dfrac{f_{m,{\bf k}}-f_{n,{\bf k}}}{E_{n,{\bf k}}-E_{m,{\bf
k}}}\left|\left\langle u_{m,{\bf k}}\left|\tau^{z}\partial_{k}H_{0}({\bf
k})\right|u_{n,{\bf k}}\right\rangle\right|^{2}\Big{]}$ (1)

where $f_{m,{\bf k}}\equiv f[E_{m,{\bf k}}]$ is the Fermi occupation factor
and $A$ is the volume normalization factor which is the area of the sample in
2D.  The stiffness measures coherence between transport in the two layers. In
the absence of inter-layer tunneling, the only way to get finite stiffness is
for the interaction to admit an off-diagonal exciton term
$\Psi^{\dagger}_{{\bf k}}[\tau^{i}\hat{\phi}({\bf
k})]\Psi^{\phantom{\dagger}}_{{\bf k}}$ with $i=x,y$ within mean-field. The
matrix function $\hat{\phi}({\bf k})$ is a mean-field ansatz that needs to be
evaluated self-consistently.

Although the dependence of the phase stiffness on wavefunctions is evident in
Eq. (1) from the matrix structure of the current operators, the geometric and
energetic terms are intertwined with no clear route to separation. To tackle
this challenge, we introduce a projected low-energy model [9], which is
applicable when excitons predominantly arise from the  lowest lying electron
band ($e$) and the highest hole band ($h$)

$\mathcal{H}_{\rm ex}=\sum\limits_{\bf k}\psi^{\dagger}_{{\bf
k}}\begin{pmatrix}\epsilon_{e}({\bf k})-\Sigma_{e}({\bf k})&\varphi({\bf
k})\\\ \varphi({\bf k})^{*}&\epsilon_{h}({\bf k})+\Sigma_{h}({\bf
k})\end{pmatrix}\psi^{\phantom{{\dagger}}}_{{\bf k}}$ (2)

where $\epsilon_{e/h}$ are the bare electron/hole dispersions,
$\psi^{\phantom{\dagger}}_{\bf k}=(c^{\phantom{\dagger}}_{{\bf
k},e},c^{\phantom{\dagger}}_{{\bf k},h})^{T}$ is the low-energy basis state,
$\Sigma_{e/h}({\bf k})$ are the self-energies (including $V_{b}$) and
$\varphi({\bf k})$ is a shorthand for

$\varphi({\bf k})=\sum\limits_{\alpha\in t,\beta\in b}[U_{t}({\bf
k})]^{*}_{e,\alpha}\phi_{\alpha\beta}({\bf k})[U_{b}({\bf k})]_{\beta,h}.$ (3)

A crucial aspect to highlight is that $U_{t}({\bf k})$ and $U_{b}({\bf k})$
are independent unitary matrices that diagonalize $H_{t}({\bf k})$ and
$H_{b}({\bf k})$ respectively. The subscripts $e$ and $h$ indicate the
electron and hole bands involved in the exciton pairing. It is important to
note that, unlike in BCS theory where the redundancy in Nambu basis introduces
a particle-hole symmetry in the BdG Hamiltonian, $U_{t}({\bf k})$ and
$U_{b}({\bf k})$ do not have such constraints. Hence our framework is a
generalization of phase stiffness, that is designed for ${\rm WSe}_{2}/{\rm
MoSe}_{2}$ bilayers where the electron and hole bands originate from distinct
materials.

We next assume $\epsilon_{e}({\bf k})=-\epsilon_{h}({\bf k})=\epsilon({\bf
k})$, while allowing $U_{t}({\bf k})$ to differ from $U_{b}({\bf k})$. This
assumption is justified as the effect of asymmetric dispersion on exciton is a
well understood textbook problem [66] that only complicates our analysis
without offering additional insight (see also App. C). On the other hand, the
presence of $U_{t}({\bf k})$ and $U_{b}({\bf k})$ has a non-trivial
consequence for $\varphi({\bf k})$ in the presence of vector potential:

$\varphi({\bf k},A)\approx\varphi({\bf k})-\dfrac{eA}{\hbar}\mathcal{P}({\bf
k})-\dfrac{e^{2}A^{2}}{2\hbar^{2}}\mathcal{D}({\bf k}),$ (4)

where $\mathcal{P}({\bf k})$ and $\mathcal{D}({\bf k})$ are the paramagnetic
and diamagnetic terms that involve derivatives of $U_{\nu}({\bf k})$ (see App.
A2 for the full expression). These terms are in addition to the energetic
terms arising from $\partial_{k}\epsilon({\bf k})$ and
$\partial_{k}^{2}\epsilon({\bf k})$ in the current operator and are key to our
analysis. After performing a lengthy but straightforward calculation that is
detailed in App. A2, we find two main contributions to the stiffness
$D_{s}=D_{s}^{c}+D_{s}^{g}$ with

$\displaystyle D_{s}^{c}$ $\displaystyle=$
$\displaystyle\dfrac{1}{2A}\sum\limits_{{\bf
k}}\partial_{k_{i}}^{2}\epsilon({\bf k})v_{\bf k}^{2}$ (5) $\displaystyle
D_{s}^{g}$ $\displaystyle=$ $\displaystyle\dfrac{1}{4A}\sum\limits_{{\bf
k}}\dfrac{\mathcal{G}({\bf k})}{E({\bf k})}+\dfrac{1}{4A}\sum\limits_{{\bf
k}}\dfrac{{\rm Re}[\mathcal{P}({\bf k})\varphi({\bf k})^{*}]^{2}}{E({\bf
k})^{3}}$ (6)

where $v_{{\bf k}}^{2}=(1-\xi({\bf k})/E({\bf k}))/2$ is the momentum
occupation factor with $\xi({\bf k})=\epsilon({\bf k})-\Sigma({\bf k})$, and
$E({\bf k})=\sqrt{\xi({\bf k})^{2}+|\varphi({\bf k})|^{2}}$ is the mean-field
quasiparticle dispersion. The second part of the equation includes the
geometric quantity $\mathcal{G}({\bf k})={\rm Re}[\mathcal{D}({\bf
k})\varphi({\bf k})^{*}]-\left|\mathcal{P}({\bf k})\right|^{2}$. We note that
all terms in eq. (6) are individually gauge invariant (see App. A4) and
constitute the main result of our paper.

A few comments are in order. The geometric contribution in Eq. (6) arises
explicitly from the Taylor expansion of the off-diagonal term $\varphi({\bf
k},A)$. If the wavefunctions were trivial, say independent of ${\bf k}$, this
contribution would vanish since $\mathcal{P}({\bf k})=\mathcal{D}({\bf k})=0$.
As one might expect, there are some connections with superconductivity as
well. If the bands were particle-hole symmetric, eq. (2) would reduce to a BdG
matrix and the function $\mathcal{G}({\bf k})$ will become precisely the
quantum metric, $g({\bf k})$. This is not surprising since superfluid
stiffness is known to be enhanced by quantum metric in superconductors within
mean-field [2]. Since the particle-hole symmetry is not enforced in an exciton
mean-field theory, our results present a generalization of quantum geometric
phase stiffness beyond pairing between particle-hole symmetric states.  We
further note that orbital embedding modifies our calculation in a similar
fashion as in superconductivity [11]. It enforces the minimal condition where
the solution to gap equation is kept real in the presence of external vector
potential [67]. Our result demonstrate a generalization of excitonic phase
stiffness beyond particle-hole symmetry discussed in the context of quantum
Hall bilayers [58].

##### Hubbard Model with tunable quantum metric.—

To illustrate the phenomenology with a simple model, we consider an interlayer
on-site orbital-diagonal Hubbard interaction $\mathcal{H}_{\rm
int}=V\sum_{{\bf i},\alpha}\hat{n}_{{\bf i},\alpha,t}\hat{n}_{{\bf
i},\alpha,b}$. The labels ${\bf i},\alpha$ pertaining to unit cell index and
orbitals, are not crucial for our discussion as long as the interaction is
interlayer in character. Details about the mean-field theory can be found in
App. B1.

For the Bloch hamiltonian, we consider a tunable quantum metric model that has
two orbitals on a square lattice with the Bloch Hamiltonian

$H(\zeta,{\bf k})=2t(2-p_{\bf k})+t_{F}\cos(\zeta p_{\bf
k})\sigma^{x}+t_{F}\sin(\zeta p_{\bf k})\sigma^{y}$ (7)

where the Pauli matrices $\sigma^{i}$ act in orbital space, $p_{\bf k}$ is a
periodic function given by $p_{\bf k}=\cos k_{x}+\cos k_{y}$ and the parameter
$\zeta$ controls the quantum geometry by introducing long range hoppings [68].
More precisely, $\zeta$ is an overall scaling factor for the quantum geometric
tensor $\mathcal{Q}({\bf k},\zeta)_{\mu\nu}=(\zeta^{2}/4)\sin k_{\mu}\sin
k_{\nu}$, which is real because of inversion and time-reversal symmetries. The
Berry curvature $F({\bf k})=-{\rm Im}[\mathcal{Q}({\bf k})]_{xy}/2$ vanishes
identically while the quantum metric $g_{\mu\nu}({\bf k})={\rm
Re}[\mathcal{Q}_{{\bf k}}]_{\mu\nu}$ is finite. Another convenient aspect of
this model is that the band dispersions $\epsilon_{\pm,{\bf k}}=2t(2-p_{\bf
k})\pm t_{F}$ are independent of $\zeta$. This permits the use of $\zeta$ to
tune quantum geometry without affecting the band dispersion, capturing the
discussion around Eq. (6).

Figure 2: ${\bf a}$ Mean-field gap solution for the toy model with tunable
quantum metric at fixed interaction $U=6t$. ${\bf b}$. Mean-field Phase
stiffness. All energy scales, $\phi$, $D_{s}$, $V_{b}$ are plotted in units of
single-particle hopping $t$ (see Eq. (48)). Subscripts $c$ and $g$ denote
conventional and geometric contributions respectively.

Symmetric bilayers, where both layers are described by the same parameter
$\zeta$, give rise to excitons (see Fig.2a) where the stiffness has a
significant enhancement coming from the geometric term $D_{s}^{g}$. The extra
contribution is indeed proportional to the trace of the quantum metric in this
particular case, as Eq.(6) simplifies to $D_{s}^{g}=(\phi^{2}/4A)\sum_{\bf
k}g({\bf k})/E({\bf k})$ [69]. While the geometric contribution dominates at
higher bias voltage $|V_{b}|$, it is important to note that a high $|V_{b}|$
results in higher densities of electron-hole pairs and may ultimately drive
the system to the exciton Mott transition [45, 46, 47, 48]. The case of
asymmetric wavefunctions follows along similar lines and is outlined in App. C
and Fig. S3.

##### Continuum Model for TMDs.—

With the formalism in place, we now aim to quantify the role of quantum
geometric effects in TMD exciton bilayers. TMDs are semi-conductors with
valley optical selection rules [70] and hence their minimal model is that of a
gapped Dirac cone. We consider the continuum Hamiltonian introduced in ref.
[71] with spin-valley locking. The continuum $k.p$ Hamiltonian for a given
spin and valley, to first order in $k$, reads

$H_{\nu}({\bf k})=\Delta_{\nu}\sigma_{z}+v_{\nu}\mathbf{k}\cdot\sigma,$ (8)

where $\sigma^{i}$ is the Pauli matrix in the internal sublattice degrees of
freedom for the conduction and valence bands, and dispersion $\epsilon_{\bf
k\nu\pm}=\pm\sqrt{\Delta^{2}_{\nu}+v_{\nu}^{2}|{\bf k}|^{2}}$. The two-band
aspect of the model is justified for the valence band where the Ising spin-
orbit coupling (SOC) gap is of the order of $180$meV in MoX2 and $440$meV in
WX2 [72]. We also neglect the quadratic trigonal warping term and the
quadratic particle-hole mass imbalance since it does not play an important
role in our phenomenological discussion.

The geometric properties of the model are encoded in the form factors
$F^{\nu}_{{\bf k},{\bf q}}=U^{\dagger}_{\nu}({\bf k})U_{\nu}({\bf k}+{\bf q})$
which has non-trivial momentum dependence.  In order to isolate the geometric
contribution, we will compare our results with the gapped Dirac cone in eq.
(8) to that of a parabolic band $H_{\nu}(\bf k)=\pm\epsilon_{\bf
k\nu}\sigma_{z}$ with energies $\epsilon_{\bf
k\nu}\approx\Delta_{\nu}+k^{2}/2m_{\nu}$ with $m_{\nu}=v^{2}_{\nu}/m_{\nu}$
and trivial form factor $F^{\nu}_{{\bf k},{\bf q}}=\delta_{{\bf k},{\bf q}}$.

Figure 3: We show the interlayer exciton coherence $\phi$ and the phase
stiffness $D_{s}$ as a function of $V_{b}$. Green and Blue data corresponds to
a gapped Dirac cone and parabolic bands with mass $m=v^{2}/\Delta$,
respectively. ${\bf a}$: Homobilayer with velocity $v=2.16$eVÅ, energy gap
$\Delta=0.9$eV relevant for MoTe2. ${\bf b}$: Heterobilayer with velocity
$v_{t}=2.2$eVÅ and $v_{b}=2.6$eVÅ, energy gap $\Delta_{t}=1.21$eV and
$\Delta_{b}=1.16$eV representative values for WSe2/MoSe2. In both cases, we
find spontaneous formation of excitons around $V^{*}_{b}\approx 70/75$meV. The
calculations were done with gate distance $\xi=12$nm and spacing between the
two layers $d=1$nm.

We overlay this model with a gated Coulomb interaction

$\mathcal{H}_{\rm int}=\dfrac{1}{2A}\sum^{\Lambda}_{{\bf
q}}\sum_{\nu,\nu^{\prime}=t,b}V^{\nu\nu^{\prime}}_{\bf q}\hat{n}_{{\bf
q},\nu}\hat{n}_{-{\bf q},\nu^{\prime}}$ (9)

where $\hat{n}_{{\bf q},\nu}$ is the layer resolved density operator
$\hat{n}_{{\bf q},\nu}=\sum^{\Lambda}_{\bf
k}\sum_{a=K,K^{\prime}}\psi^{\dagger}_{{\bf k},a\nu}\psi_{{\bf k}+{\bf
q},a\nu}$ and

$V^{tt}_{\bf q}=V^{bb}_{\bf
q}=\frac{e^{2}}{\epsilon\epsilon_{0}}\frac{\tanh\frac{q\xi}{2}}{2q},\quad
V^{tb}_{\bf q}=V^{bt}_{\bf q}=e^{-dq}V^{tt}_{\bf q}.$ (10)

Here $d$ is the inter-layer distance, $\Lambda$ is the UV cutoff, and $\xi$ is
the screening length of the bilayer which is defined as the distance between
the bilayer and the metallic gates. As representative values [43, 41], we set
$\xi=12$nm, $\epsilon\approx 6$ and $d\approx 1$nm. We focus on the regime of
low-density of electron-hole pairs  $n_{\rm e-h}\leq 0.075$ to reduce
screening effects leading to the exciton Mott transition [44, 45, 46, 47, 48].
We now proceed with a Hartree-Fock calculation. We find that while the Hartree
term vanishes as a consequence of charge neutrality, the Fock term gives rise
to a self energy correction to the single-particle Hamiltonian. The self-
consistent equations for the Fock self-energy are solved employing an
iterative scheme, decribed in App. D.

The Coulomb interaction induces intravalley excitons as well as charge
transfer between the two layers
$\phi_{z}=\langle{\tau^{z}\sigma^{0}}\rangle=2n_{\rm e-h}$. The operators
$\tau^{\pm}\sigma^{a}$ probe spontaneous interlayer coherence where
$\tau^{\pm}=(\tau^{x}\pm i\tau^{y})/2$ are raising/lowering operators in the
layer degree of freedom and $\sigma^{a}$ are Pauli matrices in the band space.
A finite expectation value $\langle\mathcal{\tau}^{\pm}\sigma^{a}\rangle\neq
0$ corresponds to an interlayer intravalley exciton condensate that breaks the
$U(1)_{e/h}$ symmetry. The first panel of Fig. 3 shows the evolution of
$\phi=\sqrt{\sum_{a=x,y,z}\langle\mathcal{\tau}^{x}\sigma^{a}\rangle^{2}}$ in
the gauge where interlayer symmetry is broken along $\tau^{x}$ as a function
of $V_{b}$ for the homobilayer (Fig. 3a) and heterobilayer (Fig. 3b) case. We
notice that for the continuum theory the energy $V_{b}$ is given by
$V_{b}=(\Delta_{t}+\Delta_{b}-2E_{z})/2$ with $E_{z}$ electric displacement
field applied to the bilayer. Above a critical value of the electrostatic
potential $V_{b}^{*}\approx 70/75$meV the system turns into a semiconductor
where the energy gap $E_{\text{gap}}$ grows linearly with the applied bias.
Finally, the phase stiffness is calculated using Eq. (1) where the expectation
value is taken with the Hartree-Fock wavefunctions. The results are shown in
the right panel of Fig. 3a-b blue data for parabolic bands and green ones for
a gapped Dirac cone. We emphasize that the non-trivial structure of the
wavefunctions appears in the geometric contribution of the superfluid
stiffness, while it does not change the exciton binding energy (energy gap)
and the size of the order parameter.

Remarkably, the geometric component increases $D_{s}$ and, correspondingly,
the BKT transition temperature $(T_{\rm BKT})$ related to the superfluid
stiffness by the Nelson-Kosterlitz relation  $k_{B}T_{\rm
BKT}=\pi/2(\lim_{T\rightarrow T_{\rm BKT}^{-}}D_{s}(T))$ [73]. Approximating
$D_{s}(T)$ with its $T=0$ value (from last panel in Fig. 3b), we infer that
geometric contribution increases $T_{\rm BKT}$ roughly from 12 to 36 Kelvin at
$V_{b}\approx 40$meV. This threefold increase should be observable in
experiments.

##### Conclusion.—

Phase stifness is significantly modified by the constituent electron and hole
wavefunctions of the exciton. We investigated two limiting scenarios, one with
local Hubbard interaction and the other with gated Coulomb interaction
projected onto the low-energy bands, and remarkably, we find geometric
contributions that play a vital role in either case. Our findings complement
existing research on the effects of quantum geometry in exciton spectrum [74,
75, 76], exciton wavefunction [77] and its possible realizations in moire
heterostructures [78, 69, 79]. Importantly, we find a novel gauge-invariant
quantum geometric quantity $\mathcal{G}({\bf k})$ whose properties require
further investigation.

Bilayer TMD devices offer an ideal platform to validate our predictions.
Quite recently, there have been two reports on exciton superfluidity in
WSe2/MoSe2 bilayer [55, 56]. The BKT temperature of $20$ K appears to be
consistent with the quantum geometric contribution, although a more
sophisticated calculation is needed for a more quantitative prediction.
Critical counterflow current from non-linear transport, similar to TBG
superconductor [80], may also provide clues for the said geometric
contribution.

##### Acknowledgment.—

Research on geometric properties of exciton condensates is supported as part
of Programmable Quantum Materials, an Energy Frontier Research Center funded
by the U.S. Department of Energy (DOE), Office of Science, Basic Energy
Sciences (BES), under award DE-SC0019443. We benefited from the 2023 Quantum
Geometry Working Group meeting that took place at the Flatiron institute. We
thank Päivi Törmä and Jonah Herzog-Arbeitman for their comments on orbital
embeddings and minimal quantum metric. We acknowledge discussions with Kin Fai
Mak, Pavel Volkov, Yongxin Zeng and Ilia Komissarov. D.G. and R.Q. acknowledge
support from the Flatiron Institute, a division of the Simons Foundation.

## References

  * Xiao  _et al._ [2010] D. Xiao, M.-C. Chang, and Q. Niu, Rev. Mod. Phys. 82, 1959 (2010).
  * Törmä  _et al._ [2022] P. Törmä, S. Peotta, and B. A. Bernevig, Nat. Rev. Phys. 4, 528 (2022).
  * Törmä  _et al._ [2018] P. Törmä, L. Liang, and S. Peotta, Phys. Rev. B 98, 220511 (2018).
  * Iskin [2018] M. Iskin, Phys. Rev. A 97, 033625 (2018).
  * Iskin [2021] M. Iskin, Phys. Rev. A 103, 053311 (2021).
  * Peotta and Törmä [2015] S. Peotta and P. Törmä, Nat. Commun. 6, 8944 (2015).
  * Julku  _et al._ [2016] A. Julku, S. Peotta, T. I. Vanhala, D.-H. Kim, and P. Törmä, Phys. Rev. Lett. 117, 045303 (2016).
  * Liang  _et al._ [2017] L. Liang, T. I. Vanhala, S. Peotta, T. Siro, A. Harju, and P. Törmä, Phys. Rev. B 95, 024515 (2017).
  * Xie  _et al._ [2020] F. Xie, Z. Song, B. Lian, and B. A. Bernevig, Phys. Rev. Lett. 124, 167002 (2020).
  * Herzog-Arbeitman  _et al._ [2022] J. Herzog-Arbeitman, V. Peri, F. Schindler, S. D. Huber, and B. A. Bernevig, Phys. Rev. Lett. 128, 087002 (2022).
  * Huhtinen  _et al._ [2022] K.-E. Huhtinen, J. Herzog-Arbeitman, A. Chew, B. A. Bernevig, and P. Törmä, Phys. Rev. B 106, 014518 (2022).
  * Mao and Chowdhury [2023a] D. Mao and D. Chowdhury, Proc. Natl. Acad. Sci. U.S.A. 120, e2217816120 (2023a).
  * Tovmasyan  _et al._ [2016] M. Tovmasyan, S. Peotta, P. Törmä, and S. D. Huber, Phys. Rev. B 94, 245149 (2016).
  * Verma  _et al._ [2021] N. Verma, T. Hazra, and M. Randeria, Proc. Natl. Acad. Sci. U.S.A. 118 (2021).
  * Ahn and Nagaosa [2021] J. Ahn and N. Nagaosa, Phys. Rev. B 104, L100501 (2021).
  * Mao and Chowdhury [2023b] D. Mao and D. Chowdhury, arXiv , 2304.07318 (2023b).
  * Hofmann  _et al._ [2020] J. S. Hofmann, E. Berg, and D. Chowdhury, Phys. Rev. B 102, 201112 (2020).
  * Peri  _et al._ [2021] V. Peri, Z.-D. Song, B. A. Bernevig, and S. D. Huber, Phys. Rev. Lett. 126, 027002 (2021).
  * Neupert  _et al._ [2013] T. Neupert, C. Chamon, and C. Mudry, Phys. Rev. B 87, 245103 (2013).
  * Komissarov  _et al._ [2023] I. Komissarov, T. Holder, and R. Queiroz, arXiv , 2306.08035 (2023).
  * Yu  _et al._ [2023] J. Yu, C. J. Ciccarino, R. Bianco, I. Errea, P. Narang, and B. A. Bernevig, arXiv , 2305.02340 (2023).
  * Arora  _et al._ [2022] A. Arora, M. S. Rudner, and J. C. W. Song, Nano Lett. 22, 9351 (2022).
  * Lapa and Hughes [2019] M. F. Lapa and T. L. Hughes, Phys. Rev. B 99, 121111 (2019).
  * Gao and Xiao [2019] Y. Gao and D. Xiao, Phys. Rev. Lett. 122, 227402 (2019).
  * Kozii  _et al._ [2021] V. Kozii, A. Avdoshkin, S. Zhong, and J. E. Moore, Phys. Rev. Lett. 126, 156602 (2021).
  * Ahn  _et al._ [2022] J. Ahn, G.-Y. Guo, N. Nagaosa, and A. Vishwanath, Nat. Phys. 18, 290 (2022).
  * Chaudhary  _et al._ [2022] S. Chaudhary, C. Lewandowski, and G. Refael, Phys. Rev. Res. 4, 013164 (2022).
  * Blatt  _et al._ [1962] J. M. Blatt, K. W. Böer, and W. Brandt, Phys. Rev. 126, 1691 (1962).
  * Halperin and Rice [1968] B. I. Halperin and T. M. Rice, Rev. Mod. Phys. 40, 755 (1968).
  * Lozovik and Yudson [1975] Y. E. Lozovik and V. Yudson, JETP Lett.(USSR)(Engl. Transl.) 22 (1975).
  * Lozovik and Yudson [1976] Y. E. Lozovik and V. Yudson, Zh. Eksp. Teor. Fiz 71, 738 (1976).
  * Yoshioka and MacDonald [1990] D. Yoshioka and A. H. MacDonald, Journal of the Physical Society of Japan 59, 4211 (1990).
  * Yang [2001] K. Yang, Phys. Rev. Lett. 87, 056802 (2001).
  * De Palo  _et al._ [2002] S. De Palo, F. Rapisarda, and G. Senatore, Phys. Rev. Lett. 88, 206401 (2002).
  * Joglekar  _et al._ [2006] Y. N. Joglekar, A. V. Balatsky, and S. Das Sarma, Phys. Rev. B 74, 233302 (2006).
  * Eisenstein and MacDonald [2004] J. P. Eisenstein and A. H. MacDonald, Nature 432, 691 (2004).
  * Eisenstein [2014] J. P. Eisenstein, Annu. Rev. Condens. Matter Phys. 5, 159 (2014), arXiv:1306.0584 .
  * Kogar  _et al._ [2017] A. Kogar, M. S. Rak, S. Vig, A. A. Husain, F. Flicker, Y. I. Joe, L. Venema, G. J. MacDougall, T. C. Chiang, E. Fradkin, J. Van Wezel, and P. Abbamonte, Science 358, 1314 (2017), arXiv:1611.04217 .
  * Jia  _et al._ [2022] Y. Jia, P. Wang, C.-L. Chiu, Z. Song, G. Yu, B. Jäck, S. Lei, S. Klemenz, F. A. Cevallos, M. Onyszczak, N. Fishchenko, X. Liu, G. Farahi, F. Xie, Y. Xu, K. Watanabe, T. Taniguchi, B. A. Bernevig, R. J. Cava, L. M. Schoop, A. Yazdani, and S. Wu, Nat. Phys. 18, 87 (2022).
  * Sun  _et al._ [2022] B. Sun, W. Zhao, T. Palomaki, Z. Fei, E. Runburg, P. Malinowski, X. Huang, J. Cenker, Y.-T. Cui, J.-H. Chu, X. Xu, S. S. Ataei, D. Varsano, M. Palummo, E. Molinari, M. Rontani, and D. H. Cobden, Nat. Phys. 18, 94 (2022).
  * Shi  _et al._ [2022] Q. Shi, E.-M. Shih, D. Rhodes, B. Kim, K. Barmak, K. Watanabe, T. Taniguchi, Z. Papić, D. A. Abanin, J. Hone, and C. R. Dean, Nat. Nanotechnol. 17, 577 (2022).
  * Wang  _et al._ [2019] Z. Wang, D. A. Rhodes, K. Watanabe, T. Taniguchi, J. C. Hone, J. Shan, and K. F. Mak, Nature 574, 76 (2019).
  * Ma  _et al._ [2021] L. Ma, P. X. Nguyen, Z. Wang, Y. Zeng, K. Watanabe, T. Taniguchi, A. H. MacDonald, K. F. Mak, and J. Shan, Nature 598, 585 (2021).
  * Mott [1949] N. F. Mott, Proc. Phys. Soc. Sec. A 62, 416 (1949).
  * Mott [1973] N. F. Mott, Contemp. Phys. 14, 401 (1973).
  * Brinkman and Rice [1973] W. F. Brinkman and T. M. Rice, Phys. Rev. B 7, 1508 (1973).
  * Smith and Wolfe [1986] L. M. Smith and J. P. Wolfe, Phys. Rev. Lett. 57, 2314 (1986).
  * Guerci  _et al._ [2019] D. Guerci, M. Capone, and M. Fabrizio, Phys. Rev. Mater. 3, 054605 (2019).
  * Nozieres and Schmitt-Rink [1985] P. Nozieres and S. Schmitt-Rink, J. Low Temp. Phys. 59, 195 (1985).
  * Fogler  _et al._ [2014] M. M. Fogler, L. V. Butov, and K. S. Novoselov, Nat. Commun. 5, 4555 (2014).
  * Wu  _et al._ [2015] F.-C. Wu, F. Xue, and A. H. MacDonald, Phys. Rev. B 92, 165121 (2015).
  * Berman and Kezerashvili [2016] O. L. Berman and R. Y. Kezerashvili, Phys. Rev. B 93, 245410 (2016).
  * Debnath  _et al._ [2017] B. Debnath, Y. Barlas, D. Wickramaratne, M. R. Neupane, and R. K. Lake, Phys. Rev. B 96, 174504 (2017).
  * High  _et al._ [2012] A. A. High, J. R. Leonard, A. T. Hammack, M. M. Fogler, L. V. Butov, A. V. Kavokin, K. L. Campman, and A. C. Gossard, Nature 483, 584 (2012).
  * Nguyen  _et al._ [2023] P. X. Nguyen, L. Ma, R. Chaturvedi, K. Watanabe, T. Taniguchi, J. Shan, and K. F. Mak, arXiv preprint arXiv:2309.14940 (2023).
  * Qi  _et al._ [2023] R. Qi, A. Y. Joe, Z. Zhang, J. Xie, Q. Feng, Z. Lu, Z. Wang, T. Taniguchi, K. Watanabe, S. Tongay, _et al._ , arXiv preprint arXiv:2309.15357 (2023).
  * Scalapino  _et al._ [1993] D. J. Scalapino, S. R. White, and S. Zhang, Phys. Rev. B 47, 7995 (1993).
  * Fertig [1989] H. A. Fertig, Phys. Rev. B 40, 1087 (1989).
  * Spielman  _et al._ [2001] I. B. Spielman, J. P. Eisenstein, L. N. Pfeiffer, and K. W. West, Phys. Rev. Lett. 87, 036803 (2001).
  * Xiao  _et al._ [2007] D. Xiao, W. Yao, and Q. Niu, Phys. Rev. Lett. 99, 236809 (2007).
  * Cao  _et al._ [2012] T. Cao, G. Wang, W. Han, H. Ye, C. Zhu, J. Shi, Q. Niu, P. Tan, E. Wang, B. Liu, and J. Feng, Nat. Commun. 3 (2012).
  * Wang  _et al._ [2018] G. Wang, A. Chernikov, M. M. Glazov, T. F. Heinz, X. Marie, T. Amand, and B. Urbaszek, Rev. Mod. Phys. 90, 021001 (2018).
  * Robert  _et al._ [2016] C. Robert, R. Picard, D. Lagarde, G. Wang, J. P. Echeverry, F. Cadiz, P. Renucci, A. Högele, T. Amand, X. Marie, I. C. Gerber, and B. Urbaszek, Phys. Rev. B 94, 155425 (2016).
  * Benfatto  _et al._ [2004] L. Benfatto, A. Toschi, and S. Caprara, Phys. Rev. B 69, 184510 (2004).
  * Cea  _et al._ [2015] T. Cea, C. Castellani, G. Seibold, and L. Benfatto, Phys. Rev. Lett. 115, 157002 (2015).
  * Fox [2001] M. A. M. Fox,  _Optical properties of solids_ (Oxford University Press, 2001) Chap. 4.
  * Peotta [2022] S. Peotta, New Journal of Physics 24, 113019 (2022).
  * Hofmann  _et al._ [2023] J. S. Hofmann, E. Berg, and D. Chowdhury, Phys. Rev. Lett. 130, 226001 (2023).
  * Hu  _et al._ [2022] X. Hu, T. Hyart, D. I. Pikulin, and E. Rossi, Phys. Rev. B 105, L140506 (2022).
  * Xiao  _et al._ [2012] D. Xiao, G.-B. Liu, W. Feng, X. Xu, and W. Yao, Phys. Rev. Lett. 108, 196802 (2012).
  * Zhu  _et al._ [1995] X. Zhu, P. B. Littlewood, M. S. Hybertsen, and T. M. Rice, Phys. Rev. Lett. 74, 1633 (1995).
  * Kormá nyos  _et al._ [2015] A. Kormá nyos, G. Burkard, M. Gmitra, J. Fabian, V. Zólyomi, N. D. Drummond, and V. Fal’ko, 2D Mater. 2, 022001 (2015).
  * Nelson and Kosterlitz [1977] D. R. Nelson and J. M. Kosterlitz, Phys. Rev. Lett. 39, 1201 (1977).
  * Zhou  _et al._ [2015] J. Zhou, W.-Y. Shan, W. Yao, and D. Xiao, Phys. Rev. Lett. 115, 166803 (2015).
  * Srivastava and Imamoğlu [2015] A. Srivastava and A. Imamoğlu, Phys. Rev. Lett. 115, 166802 (2015).
  * Yao and Niu [2008] W. Yao and Q. Niu, Phys. Rev. Lett. 101, 106401 (2008).
  * Sethi  _et al._ [2023] G. Sethi, M. Cuma, and F. Liu, Phys. Rev. Lett. 130, 186401 (2023).
  * Kwan  _et al._ [2021] Y. H. Kwan, Y. Hu, S. H. Simon, and S. A. Parameswaran, Phys. Rev. Lett. 126, 137601 (2021).
  * Hu  _et al._ [2023] X. Hu, E. Rossi, and Y. Barlas, arXiv preprint arXiv:2304.04825 (2023).
  * Tian  _et al._ [2023] H. Tian, X. Gao, Y. Zhang, S. Che, T. Xu, P. Cheung, K. Watanabe, T. Taniguchi, M. Randeria, F. Zhang, C. N. Lau, and M. W. Bockrath, Nature 614, 440 (2023).
  * Sun  _et al._ [2021] Z. Sun, T. Kaneko, D. Golež, and A. J. Millis, Phys. Rev. Lett. 127, 127702 (2021).
  * Simon and Rudner [2020] S. H. Simon and M. S. Rudner, Phys. Rev. B 102, 165148 (2020).
  * Chandrasekhar [1962] B. S. Chandrasekhar, Appl. Phys. Lett. 1, 7 (1962).
  * Clogston [1962] A. M. Clogston, Phys. Rev. Lett. 9, 266 (1962).
  * Cheiwchanchamnangij and Lambrecht [2012] T. Cheiwchanchamnangij and W. R. L. Lambrecht, Phys. Rev. B 85, 205302 (2012).
  * Liang  _et al._ [2013] Y. Liang, S. Huang, R. Soklaski, and L. Yang, Appl. Phys. Lett. 103, 042106 (2013).

Supplementary material for “”

Nishchhal Verma,1,∗ Daniele Guerci,2,∗ and Raquel Queiroz1,2

1Department of Physics, Columbia University, New York, NY 10027, USA

2Center for Computational Quantum Physics, Flatiron Institute, New York, New
York 10010, USA

##  Appendix S-I Mean-field Phase Stiffness

###  S-I.1 Kubo formula

Exciton condensates exhibit the remarkable phenomenon of dissipationless
counterflow transport when equal and opposite fields are applied to the two
layers. The asymmetric field drives excitons in a coherent fashion. The
response kernel to layer anti-symmetric vector potential is defined as [37]

$j_{{\rm CF}}=-\dfrac{4e^{2}}{\hbar^{2}}D_{s}A$ (1)

where $D_{s}$ is the phase stiffness and ${\bf A}=A\hat{x}$ is the external
vector potential that is applied along $+\hat{x}$ in top layer and $-\hat{x}$
in bottom layer. As a matter of fact, this can be an in-plane magnetic field
as ${\bf A}=zA\hat{x}$ gives $B_{y}=\partial_{z}A_{x}$ where $z$ is the
direction normal to the interface [81]. The stiffness is the weight of the
delta function in counterflow conductivity ${\rm Re}[\sigma_{{\rm
CF}}(\omega)]=D_{s}\delta(\omega)+\cdots$. By utilizing the Kramers-Kronig
relation, $D_{s}$ can be extracted from the $1/\omega$ tail of ${\rm
Im}[\sigma_{CF}(\omega)]$, which follows ${\rm Im}[\sigma_{CF}(\omega)]\propto
D_{s}/\omega$. Stiffness controls the cost of phase fluctuations of the
condensate and is formally defined by the Kubo formula [57]

$D_{s}=-\dfrac{1}{4}\Big{(}\langle j_{D}\rangle-\chi_{j_{P}j_{P}}({\bf
q}_{\perp}=0,\omega=0)\Big{)}.$ (2)

The diamagnetic $(j_{D})$ and paramagnetic $(j_{P})$ currents are given by a
Taylor expansion of the Hamiltonian

$\displaystyle\mathcal{H}_{\rm ex}(A)$ $\displaystyle=$
$\displaystyle\sum\limits_{{\bf k}}\Psi^{\dagger}_{{\bf
k}}\left[\begin{pmatrix}H_{t}({\bf k}-eA/\hbar)&0\\\ 0&H_{b}({\bf
k}+eA/\hbar)\end{pmatrix}+V_{b}\tau^{z}\right]\Psi^{\dagger}_{{\bf
k}}+\mathcal{H}_{\rm int}$ (3) $\displaystyle=$ $\displaystyle\mathcal{H}_{\rm
ex}(0)-\dfrac{eA}{\hbar}\sum\limits_{{\bf k}}\Psi^{\dagger}_{{\bf
k}}\begin{pmatrix}\partial_{k}H_{t}({\bf k})&0\\\ 0&-\partial_{k}H_{b}({\bf
k})\end{pmatrix}\Psi^{\phantom{\dagger}}_{{\bf
k}}+\dfrac{e^{2}A^{2}}{2\hbar^{2}}\sum\limits_{{\bf k}}\Psi^{\dagger}_{{\bf
k}}\begin{pmatrix}\partial_{k}^{2}H_{t}({\bf k})&0\\\
0&\partial_{k}^{2}H_{b}({\bf k})\end{pmatrix}\Psi^{\phantom{\dagger}}_{{\bf
k}}$ (4)

which gives

$\displaystyle j_{P}=\dfrac{e}{\hbar}\sum\limits_{{\bf k}}\Psi^{\dagger}_{{\bf
k}}\begin{pmatrix}\partial_{k}H_{t}({\bf k})&0\\\ 0&-\partial_{k}H_{b}({\bf
k})\end{pmatrix}\Psi^{\phantom{\dagger}}_{{\bf k}},\quad
j_{D}=-\dfrac{e^{2}}{\hbar^{2}}\sum\limits_{{\bf k}}\Psi^{\dagger}_{{\bf
k}}\begin{pmatrix}\partial_{k}^{2}H_{t}({\bf k})&0\\\
0&\partial_{k}^{2}H_{b}({\bf k})\end{pmatrix}\psi^{\phantom{\dagger}}_{{\bf
k}}.$ (5)

Next, in order to calculate the expectation values in Eq. (2), we assume a
mean-field Hamiltonian with eigenspectrum $E_{m,{\bf k}}$ and eigenstates
$|u_{m,{\bf k}}\rangle\\}$ to find

$\langle j_{D}\rangle=\dfrac{1}{A}\sum\limits_{{\bf k},m}f[E_{m,{\bf
k}}]\left\langle u_{m,{\bf k}}\left|\partial^{2}_{k}H_{0}({\bf
k})\right|u_{m,{\bf
k}}\right\rangle,\quad\chi_{j_{P}j_{P}}=\dfrac{1}{A}\sum\limits_{{\bf k},m\neq
n}\dfrac{f[E_{m,{\bf k}}]-f[E_{n,{\bf k}}]}{E_{n,{\bf k}}-E_{m,{\bf
k}}}\left|\left\langle u_{m,{\bf k}}\left|\tau^{z}\partial_{k}H_{0}({\bf
k})\right|u_{n,{\bf k}}\right\rangle\right|^{2}$ (6)

which gives Eq. (1) of the main text. The expression not only captures the
conventional contribution coming from energies but also includes all geometric
contributions that are hidden inside the matrix elements. These geometric
factors arise from the underlying electron and hole wavefunctions that are not
directly visible in this equation. We will now describe a projected mean-field
model which systematically uncovers these geometric contributions.

###  S-I.2 Projected Phase Stiffness

If the electron and hole bands are sufficiently isolated from other bands, we
can focus on the low-energy Hamiltonian and couple to the external vector
potential

$\mathcal{H}_{\rm ex}(A)=\sum\limits_{{\bf k}}\psi^{\dagger}_{\bf
k}\begin{pmatrix}\epsilon({\bf k}-eA/\hbar)-\Sigma({\bf k})&\varphi({\bf
k},A)\\\ \phi\;\varphi({\bf k},A)^{*}&-\Big{(}\epsilon({\bf
k}+eA/\hbar)-\Sigma({\bf k})\Big{)}\end{pmatrix}\psi^{\phantom{\dagger}}_{{\bf
k}}=\sum\limits_{{\bf k}}\psi^{\dagger}_{\bf k}H_{\rm
ex}(A)\psi^{\phantom{\dagger}}_{{\bf k}}.$ (7)

The vector potential dependence on the off-diagonal component is key to the
geometric contribution

$\displaystyle\varphi({\bf k},A)$ $\displaystyle=$
$\displaystyle\sum_{\alpha\beta}[U_{t}({\bf
k}-eA/\hbar)]^{*}_{e,\alpha}\phi_{\alpha\beta}({\bf k})[U_{b}({\bf
k}+eA/\hbar)]_{\beta,h}=\varphi({\bf k})-\dfrac{eA}{\hbar}\mathcal{P}({\bf
k})-\dfrac{e^{2}A^{2}}{2\hbar^{2}}\mathcal{D}({\bf k})$ (8)

where $\mathcal{P}({\bf k})$ and $\mathcal{D}({\bf k})$ are defined as

$\displaystyle\mathcal{P}({\bf k})=\sum_{\alpha\beta}\partial_{k}[U_{t}({\bf
k})]^{*}_{e,\alpha}\phi_{\alpha\beta}({\bf k})[U_{b}({\bf
k})]_{\beta,h}-[U_{t}({\bf k})]^{*}_{e,\alpha}\phi_{\alpha\beta}({\bf
k})\partial_{k}[U_{b}({\bf k})]_{\beta,h}\hskip 155.06743pt$ (9)
$\displaystyle\mathcal{D}({\bf k})=\sum_{\alpha\beta}2\partial_{k}[U_{t}({\bf
k})]^{*}_{e,\alpha}\phi_{\alpha\beta}({\bf k})\partial_{k}[U_{b}({\bf
k})]_{\beta,h}-\partial_{k}^{2}[U_{t}({\bf
k})]^{*}_{e,\alpha}\phi_{\alpha\beta}({\bf k})[U_{b}({\bf
k})]_{\beta,h}-[U_{t}({\bf k})]^{*}_{e,\alpha}\phi_{\alpha\beta}({\bf
k})\partial_{k}^{2}[U_{b}({\bf k})]_{\beta,h}.$ (10)

Moving onto the stiffness calculation, we will use the thermodynamic
definition instead of using the Kubo formula. Focusing on $T=0$ and assuming
that the mean-field band is completely gapped (even for $H_{\rm ex}({\bf
k},A)$), the stifness is given by

$D_{s}=\dfrac{\hbar^{2}}{4e^{2}}\left(\dfrac{\partial^{2}}{\partial
A^{2}}\mathcal{F}(A)\right)=\dfrac{\hbar^{2}}{4e^{2}}\sum\limits_{{\bf
k}}\dfrac{\partial^{2}}{\partial A^{2}}\mathcal{E}_{-}({\bf k},(A))$ (11)

where $\mathcal{E}_{-}({\bf k},A)$ is the lower band corresponding to $H_{\rm
ex}({\bf k},A)$, given by

$\mathcal{E}_{\pm}({\bf k},A)=\xi_{-}({\bf k},A)\pm\sqrt{\xi_{+}({\bf
k},A)^{2}+|\varphi({\bf k},A)|^{2}}$ (12)

where

$\xi_{\pm}({\bf k},A)=\dfrac{1}{2}\Big{[}(\epsilon({\bf
k}-eA/\hbar)-\Sigma({\bf k}))\pm(\epsilon({\bf k}+eA/\hbar)-\Sigma({\bf
k}))\Big{]}$ (13)

that has convenient properties like

$\displaystyle\xi_{-}({\bf k},0)=0,$ $\displaystyle\xi_{+}({\bf
k},0)=\epsilon({\bf k})-\Sigma({\bf k})\equiv\xi({\bf k})$ (14)
$\displaystyle\partial_{A}\xi_{-}({\bf
k},A)|_{A=0}=-\dfrac{e}{\hbar}\partial_{k}\epsilon({\bf k}),$
$\displaystyle\partial_{A}\xi_{+}({\bf k},A)|_{A=0}=0$ (15)
$\displaystyle\partial_{A}^{2}\xi_{-}({\bf k},A)|_{A=0}=0,$
$\displaystyle\partial_{A}^{2}\xi_{+}({\bf
k},A)|_{A=0}=\dfrac{e^{2}}{\hbar^{2}}\partial_{k}^{2}\epsilon({\bf k}).$ (16)

We are now ready to evaluate $D_{s}$:

$\displaystyle D_{s}$ $\displaystyle=$
$\displaystyle\dfrac{\hbar^{2}}{4e^{2}}\dfrac{\partial^{2}}{\partial
A^{2}}\left(\sum\limits_{\bf k}\mathcal{E}_{-}({{\bf
k}},A)\right)_{A=0}=\dfrac{\hbar^{2}}{4e^{2}}\dfrac{\partial^{2}}{\partial
A^{2}}\left(\sum\limits_{\bf k}(\xi_{-}({\bf k},A)-\sqrt{\xi_{+}({\bf
k},A)^{2}+|\varphi({\bf k},A)|^{2}}\right)_{A=0}$ (17) $\displaystyle=$
$\displaystyle\dfrac{\hbar^{2}}{4e^{2}}\sum\limits_{\bf
k}\partial_{A}^{2}\xi_{-}({\bf k},A)|_{A=0}-\dfrac{\partial^{2}}{\partial
A^{2}}\left(\sqrt{\xi_{+}({\bf k},A)^{2}+|\varphi({\bf
k},A)|^{2}}\right)_{A=0}$ (18) $\displaystyle=$
$\displaystyle-\dfrac{1}{4}\sum\limits_{\bf k}\partial_{k}^{2}\epsilon({\bf
k})\dfrac{\xi({\bf k})}{E({\bf k})}+\dfrac{\hbar^{2}}{8e^{2}}\sum\limits_{\bf
k}\dfrac{(-\partial_{A}^{2}|\varphi({\bf k},A)|^{2})}{E({\bf
k})}+\dfrac{\hbar^{2}}{4e^{2}}\sum\limits_{\bf
k}\dfrac{(\partial_{A}|\varphi({\bf k},A)|^{2}/2)^{2}}{E({\bf k})^{3}}.$ (19)

Using the fact that $\sum_{\bf k}\partial_{k}^{2}\epsilon({\bf k})=0$, we add
and subtract this term to write

$D_{s}=\dfrac{1}{4}\sum\limits_{\bf k}\partial_{k}^{2}\epsilon({\bf
k})\left[1-\dfrac{\xi({\bf k})}{E({\bf
k})}\right]+\dfrac{\hbar^{2}}{8e^{2}}\sum\limits_{\bf
k}\dfrac{(-\partial_{A}^{2}|\varphi({\bf k},A)|^{2})}{E({\bf
k})}+\dfrac{\hbar^{2}}{4e^{2}}\sum\limits_{\bf
k}\dfrac{(\partial_{A}|\varphi({\bf k},A)|^{2}/2)^{2}}{E({\bf k})^{3}}$ (20)

where we can further use

$|\varphi({\bf k},A)|^{2}=|\varphi({\bf k})|^{2}-2\dfrac{eA}{\hbar}{\rm
Re}[[\mathcal{P}({\bf k})\varphi({\bf
k})^{*}]-\dfrac{e^{2}A^{2}}{\hbar^{2}}\Big{(}{\rm Re}[[\mathcal{D}({\bf
k})\varphi({\bf k})^{*}]-|\mathcal{P}({\bf k})|^{2}\Big{)}$ (21)

to arrive at the result quoted in the main text

$D_{s}=\dfrac{1}{2}\sum\limits_{\bf k}\partial_{k}^{2}\epsilon({\bf k})v_{\bf
k}^{2}+\dfrac{1}{4}\sum\limits_{\bf k}\dfrac{{\rm Re}[\mathcal{D}({\bf
k})\varphi({\bf k})^{*}]-|\mathcal{P}({\bf k})|^{2}}{E({\bf
k})}+\dfrac{1}{4}\sum\limits_{\bf k}\dfrac{{\rm Re}[\mathcal{P}({\bf
k})\varphi({\bf k})^{*}]^{2}}{E({\bf k})^{3}}.$ (22)

###  S-I.3 Reduction to Quantum Metric

We can validate our calculations by making two assumptions:
first,$\hat{\phi}({\bf k})=\phi\mathbb{1},\;\phi\in\mathbb{R}$ and second, the
two layers are identical, giving $[U_{t}({\bf k})]_{e,\alpha}=[U_{b}({\bf
k})]_{h,\alpha}\;\forall\;\alpha,{\bf k}$. We can express the matrix elements
in a more concise form (dropping $e$, $h$ labels):

$\displaystyle\varphi({\bf k})=\phi,\quad\mathcal{P}({\bf
k})=2\phi\sum_{\alpha}\partial_{k}U^{*}_{\alpha}({\bf k})U_{\alpha}({\bf
k})\equiv 2\phi\langle\partial_{k}u_{\bf k}|u_{\bf
k}\rangle,\quad\mathcal{D}({\bf
k})=4\phi\sum_{\alpha}\partial_{k}U^{*}_{\alpha}({\bf
k})\partial_{k}U_{\alpha}({\bf k})\equiv 4\phi\langle\partial_{k}u_{\bf
k}|\partial_{k}u_{\bf k}\rangle$ (23)

and the resulting geometric contribution to stiffness greatly simplifies

$\displaystyle D_{s}^{g}$ $\displaystyle=$
$\displaystyle\dfrac{1}{4A}\sum\limits_{\bf k}\dfrac{1}{E({\bf k})}\left({\rm
Re}[\mathcal{D}({\bf k})\varphi({\bf k})^{*}]-\left|{\rm Im}[\mathcal{P}({\bf
k})]\right|^{2}\right)=\dfrac{1}{A}\sum\limits_{\bf k}\dfrac{\phi^{2}}{E({\bf
k})}\left({\rm Re}\left[\langle\partial_{k}u_{\bf k}|\partial_{k}u_{\bf
k}\rangle\right]-\left|\langle\partial_{k}u_{\bf k}|u_{\bf
k}\rangle\right|^{2}\right)$ (24) $\displaystyle=$
$\displaystyle\dfrac{1}{A}\sum\limits_{\bf k}\dfrac{\phi^{2}}{E({\bf k})}{\rm
Re}[\mathcal{Q}({\bf k})].$ (25)

Here ${\rm Re}[\mathcal{Q}({\bf k})]=g({\bf k})$ is the quantum metric of the
electron (or hole) band.

###  S-I.4 Gauge Invariance

There is a gauge redundancy in non-interacting electron/hole wavefunctions
since Bloch wavefunctions are only defined upto a phase. Physical quantities
must be invariant under such gauge transformations. Specifically, introducing
phases in the wavefunctions should not modify any observable

$U_{t}({\bf k})_{\alpha,e}\rightarrow e^{i\phi_{t}({\bf k})}U_{t}({\bf
k})_{\alpha,e},\quad U_{b}({\bf k})_{\beta,h}\rightarrow e^{i\phi_{b}({\bf
k})}U_{b}({\bf k})_{\beta,h}$ (26)

where $\phi_{\nu}({\bf k})$ are real but arbitrary functions in the BZ. It
turns out that $\varphi({\bf k}),\mathcal{P}({\bf k})$ and $\mathcal{D}({\bf
k})$ are sensitive to the gauge choice

$\displaystyle\varphi({\bf k})$ $\displaystyle\rightarrow$ $\displaystyle
e^{-i(\phi_{t}({\bf k})-\phi_{b}({\bf k}))}\varphi({\bf k})$ (27)
$\displaystyle\mathcal{P}({\bf k})$ $\displaystyle\rightarrow$ $\displaystyle
e^{-i(\phi_{t}({\bf k})-\phi_{b}({\bf k}))}\Big{[}\mathcal{P}({\bf
k})-i\varphi({\bf k})[\partial_{k}\phi_{t}({\bf k})+\partial_{k}\phi_{b}({\bf
k})]\Big{]}$ (28) $\displaystyle\mathcal{D}({\bf k})$
$\displaystyle\rightarrow$ $\displaystyle e^{-i(\phi_{t}({\bf
k})-\phi_{b}({\bf k}))}\Big{[}\mathcal{D}({\bf k})-\varphi({\bf
k})[\partial_{k}\phi_{t}({\bf k})+\partial_{k}\phi_{b}({\bf
k})]^{2}+2i[\partial_{k}\phi_{t}({\bf k})+\partial_{k}\phi_{b}({\bf
k})]\mathcal{P}({\bf k})\Big{]}.$ (29)

However, the combinations that enter geometric contributions to phase
stiffness are gauge invariant. Eq. (21) offers a simple proof. The left-hand
side is gauge invariant and hence the right-hand side has to be gauge
invariant for arbitrary expansions in $A$.

It can also be checked by an explicit calculation. The first term, ${\rm
Re}[\mathcal{P}({\bf k})\varphi({\bf k})^{*}]$, gets an additional term which
is purely imaginary

$\displaystyle\mathcal{P}({\bf k})\varphi({\bf
k})^{*}\rightarrow\mathcal{P}({\bf k})\varphi({\bf k})^{*}-i|\varphi({\bf
k})|^{2}[\partial_{k}\phi_{t}({\bf k})-\partial_{k}\phi_{b}({\bf k})]$ (30)

and hence ${\rm Re}[\mathcal{P}({\bf k})\varphi({\bf k})^{*}]$ is gauge
invariant. The second term follows along similar lines

$\displaystyle\mathcal{D}({\bf k})\varphi({\bf k})^{*}$
$\displaystyle\rightarrow$ $\displaystyle\mathcal{D}({\bf k})\varphi({\bf
k})^{*}-|\varphi({\bf k})|^{2}[\partial_{k}\phi_{t}({\bf
k})+\partial_{k}\phi_{b}({\bf k})]^{2}+2i\mathcal{P}({\bf k})\varphi({\bf
k})^{*}[\partial_{k}\phi_{t}({\bf k})+\partial_{k}\phi_{b}({\bf k})]$ (31)
$\displaystyle|\mathcal{P}({\bf k})|^{2}$ $\displaystyle\rightarrow$
$\displaystyle|\mathcal{P}({\bf k})|^{2}+|\varphi({\bf
k})|^{2}[\partial_{k}\phi_{t}({\bf k})+\partial_{k}\phi_{b}({\bf
k})]^{2}-2{\rm Im}[\mathcal{P}({\bf k})\varphi({\bf
k})^{*}][\partial_{k}\phi_{t}({\bf k})+\partial_{k}\phi_{b}({\bf k})]$ (32)

where the last two terms exactly cancel and ${\rm Re}[\mathcal{D}({\bf
k})\varphi({\bf k})^{*}]-|\mathcal{P}({\bf k})|^{2}$ is gauge invariant.

###  S-I.5 Orbital Embedding

In addition to the gauge freedom mentioned above, there is an another
transformation that changes the wavefunctions, but does not affect the
energies

$U_{t}({\bf k})_{\alpha,e}\rightarrow e^{i\phi_{t,\alpha}({\bf k})}U_{t}({\bf
k})_{\alpha,e},\quad U_{b}({\bf k})_{\beta,h}\rightarrow
e^{i\phi_{b,\beta}({\bf k})}U_{b}({\bf k})_{\beta,h}$ (33)

where the phases depend on orbital index $\alpha$ (and $\beta$). The
transformation corresponds to physically moving the orbitals within the unit
cell while preserving all hopping amplitudes. It changes Berry curvature [82],
quantum metric [11] and the geometric contribution to stiffness. It has been
argued that orbital embedding should not change the phase stiffness [67] and
the correct prescription is discussed in detail in ref. [11]. The gauge
physically corresponds to orbital locations that keep order parameter
solutions real in the presence of vector potential.

###  S-I.6 Spatially Anisotropic Systems

Underneath all our discussion so far, we have assumed isotropic systems where
the vector potential is applied along $\hat{x}$ with ${\bf A}=A\hat{x}$ and
the longitudinal response is calculated along $\hat{x}$ while suppressing all
spatial indices. In this section, we will outline the extension to spatially
anisotropic systems. Introducing indices $i,j=\\{x,y\\}$, we find that the
phase stiffness in Eq. (A6) becomes a tensor

$\displaystyle[D_{s}]_{ij}$ $\displaystyle=$
$\displaystyle-\dfrac{1}{4A}\Big{(}\langle[j_{D}]_{ij}\rangle-\chi_{[j_{P}]_{i}[j_{P}]_{j}}({\bf
q}_{\perp}=0,\omega=0)\Big{)}$ $\displaystyle=$
$\displaystyle\dfrac{\hbar^{2}}{4e^{2}A}\Big{[}\sum\limits_{{\bf
k},m}f[E_{m,{\bf k}}]\left\langle u_{m,{\bf k}}\left|[\hat{j}_{D}({\bf
k})]_{ij}\right|u_{m,{\bf k}}\right\rangle$
$\displaystyle\quad\quad\quad\quad-\sum\limits_{{\bf k},m\neq
n}\dfrac{f[E_{m,{\bf k}}]-f[E_{n,{\bf k}}]}{E_{n,{\bf k}}-E_{m,{\bf
k}}}\left\langle u_{m,{\bf k}}\left|[\hat{j}_{P}({\bf k})]_{i}\right|u_{n,{\bf
k}}\right\rangle\left\langle u_{n,{\bf k}}\left|[\hat{j}_{P}({\bf
k})]_{j}\right|u_{m,{\bf k}}\right\rangle\Big{]}$ (35)

where the current operators are

$\displaystyle[\hat{j}_{P}({\bf
k})]_{i}=\dfrac{e}{\hbar}\begin{pmatrix}\partial_{k_{i}}H_{t}({\bf k})&0\\\
0&-\partial_{k_{i}}H_{b}({\bf k})\end{pmatrix},\quad[\hat{j}_{D}({\bf
k})]_{ij}=-\dfrac{e^{2}}{\hbar^{2}}\begin{pmatrix}\partial_{k_{i}}\partial_{k_{j}}H_{t}({\bf
k})&0\\\ 0&\partial_{k_{i}}\partial_{k_{j}}H_{b}({\bf k})\end{pmatrix}.$ (36)

The projected stiffness in Eq. (22) follows along similar lines

$\displaystyle[D_{s}]_{ij}=\dfrac{1}{4A}\sum\limits_{\bf
k}\left[\partial_{k_{i}}\partial_{k_{j}}\epsilon({\bf
k})\left(1-\dfrac{\xi({\bf k})}{E({\bf k})}\right)+\sum\limits_{\bf
k}\dfrac{[\mathcal{G}({\bf k})]_{ij}}{E({\bf k})}+\sum\limits_{\bf
k}\dfrac{{\rm Re}[[\mathcal{P}({\bf k})]_{i}\varphi({\bf k})^{*}]{\rm
Re}[[\mathcal{P}({\bf k})]_{j}\varphi({\bf k})^{*}]}{E({\bf k})^{2}}\right]$
(37)

where $[\mathcal{G}({\bf k})]_{ij}={\rm Re}[\mathcal{D}({\bf
k})]_{ij}\varphi({\bf k})^{*}]-[\mathcal{P}({\bf k})]_{i}[\mathcal{P}({\bf
k})]_{j}^{*}$ and the functions $[\mathcal{P}({\bf k})]_{i}$,
$[\mathcal{D}({\bf k})]_{ij}$ are

$\displaystyle[\mathcal{P}({\bf k})]_{i}$ $\displaystyle=$
$\displaystyle\sum_{\alpha\beta}\Big{[}\partial_{k_{i}}[U_{t}({\bf
k})]^{*}_{e,\alpha}\phi_{\alpha\beta}({\bf k})[U_{b}({\bf
k})]_{\beta,h}-[U_{t}({\bf k})]^{*}_{e,\alpha}\phi_{\alpha\beta}({\bf
k})\partial_{k_{i}}[U_{b}({\bf k})]_{\beta,h}\Big{]}$ (38)
$\displaystyle[\mathcal{D}({\bf k})]_{ij}$ $\displaystyle=$
$\displaystyle\sum_{\alpha\beta}\Big{[}\partial_{k_{i}}[U_{t}({\bf
k})]^{*}_{e,\alpha}\phi_{\alpha\beta}({\bf k})\partial_{k_{j}}[U_{b}({\bf
k})]_{\beta,h}+\partial_{k_{j}}[U_{t}({\bf
k})]^{*}_{e,\alpha}\phi_{\alpha\beta}({\bf k})\partial_{k_{i}}[U_{b}({\bf
k})]_{\beta,h}$ (39)
$\displaystyle-\partial_{k_{i}}\partial_{k_{j}}[U_{t}({\bf
k})]^{*}_{e,\alpha}\phi_{\alpha\beta}({\bf k})[U_{b}({\bf
k})]_{\beta,h}-[U_{t}({\bf k})]^{*}_{e,\alpha}\phi_{\alpha\beta}({\bf
k})\partial_{k_{i}}\partial_{k_{j}}[U_{b}({\bf k})]_{\beta,h}\Big{]}.$

Lastly, the BKT temperature is determined by the determinant of the stiffness
tensor $k_{B}T_{c}=(\pi/2)\sqrt{{\rm det}[D_{s}(T_{c}^{-})]_{ij}}$.

##  Appendix S-II Mean-Field theory for the Hubbard interaction

We consider an interlayer on-site orbital-diagonal Hubbard interaction
$\mathcal{H}_{\rm int}=V\sum_{{\bf i},\alpha}\hat{n}_{{\bf
i},\alpha,t}\hat{n}_{{\bf i},\alpha,b}$. The labels ${\bf i},\alpha$
pertaining to unit cell index and orbitals, are not crucial for our discussion
as long as the interaction is interlayer in character. We decouple the four-
fermion interaction with a spatially uniform ansatz

$\displaystyle V_{\alpha\beta}\;c^{\dagger}_{{\bf
i},\alpha,t}c^{\phantom{\dagger}}_{{\bf i},\alpha,t}\;c^{\dagger}_{{\bf
i},\beta,b}c^{\phantom{\dagger}}_{{\bf i},\beta,b}$ $\displaystyle\approx$
$\displaystyle V_{\alpha\beta}\;n_{\alpha,t}c^{\dagger}_{{\bf
i},\beta,b}c^{\phantom{\dagger}}_{{\bf
i},\beta,b}+V_{\alpha\beta}c^{\dagger}_{{\bf
i},\alpha,t}c^{\phantom{\dagger}}_{{\bf
i},\alpha,t}n_{\beta,b}+\phi_{\alpha\beta}c^{\dagger}_{{\bf
i},\beta,b}c^{\phantom{\dagger}}_{{\bf
i},\alpha,t}+\phi_{\alpha\beta}^{\dagger}c^{\dagger}_{{\bf
i},\alpha,t}c^{\phantom{\dagger}}_{{\bf i},\beta,b}$ (40)

where $\\{n_{\alpha,\nu}\\}$ and $\\{\phi_{\alpha\beta}\\}$ are given by

$n_{\alpha,\nu}=\dfrac{1}{A}\sum\limits_{{\bf i}}\left\langle
c^{\dagger}_{{\bf i},\alpha,\nu}c^{\phantom{\dagger}}_{{\bf
i},\alpha,\nu}\right\rangle=\dfrac{1}{N}\sum\limits_{{\bf k}}\left\langle
c^{\dagger}_{{\bf k},\alpha,\nu}c^{\phantom{\dagger}}_{{\bf
k},\alpha,\nu}\right\rangle,\quad\phi_{\alpha\beta}=-\dfrac{V_{\alpha\beta}}{A}\sum\limits_{{\bf
i}}\left\langle c^{{\dagger}}_{{\bf i},\alpha,t}c^{\phantom{{\dagger}}}_{{\bf
i},\beta,b}\right\rangle=-\dfrac{V_{\alpha\beta}}{N}\sum\limits_{{\bf
k}}\left\langle c^{{\dagger}}_{{\bf k},\alpha,t}c^{\phantom{{\dagger}}}_{{\bf
k},\beta,b}\right\rangle.$ (41)

Here $A=N\Omega$ is the area of the system, $N$ is the number of unit cells
and $\Omega$ is the area of unit cell. This pre-factor serves as a
normalization factor. When we substitute this result back into the full
Hamiltonian, we obtain the Hartree shifts

$\displaystyle[\widetilde{H}(\zeta_{t},{\bf
k})]_{\alpha\alpha^{\prime}}=[H(\zeta_{t},{\bf
k})]_{\alpha\alpha^{\prime}}+\delta_{\alpha,\alpha^{\prime}}V_{\alpha\beta}n_{\beta,b},\quad[\widetilde{H}(\zeta_{b},{\bf
k})]_{\beta\beta^{\prime}}=[H(\zeta_{b},{\bf
k})]_{\beta\beta^{\prime}}+\delta_{\beta,\beta^{\prime}}V_{\alpha\beta}n_{\alpha,t}$
(42)

in addition to the interlayer exciton term
$\propto\phi_{\alpha\beta}c^{\dagger}_{{\bf
i},\beta,b}c^{\phantom{\dagger}}_{{\bf i},\alpha,t}$. Combining the two, we
get the mean-field Hamiltonian $\mathcal{H}_{\rm MF}=\sum_{{\bf
k}}\Psi^{\dagger}_{{\bf k}}H_{\rm MF}({\bf k})\Psi^{\phantom{{\dagger}}}_{{\bf
k}}$ where

$H_{\rm MF}({\bf k})=[\widetilde{H}(\zeta_{t},{\bf
k})\oplus(-\widetilde{H}(\zeta_{b},{\bf k}))]+\tau^{z}V_{b}+\tau^{x}\phi.$
(43)

The last term breaks the $U(1)_{t}\times U(1)_{b}$ symmetry as discussed
before.

To establish the notation, let us review the simplest model consisting of two
identical bands with mass $m$: one particle-like and one hole-like, that
belong to two different layers. We consider the system at charge neutrality,
which means that the sum of the particle density $n_{t}$ and the hole density
$n_{b}$ equals 1, and introduce an inter-layer bias $V_{b}$. The interaction
term is $U\hat{n}_{t}(\hat{n}_{b}-1)$, specifically chosen to vanish when the
lower band is fully occupied ($n_{b}=1$). We then decouple the interaction
into two components: a real exciton field denoted by $\phi$, which can be
expressed as $\phi=-U(\langle
c^{\phantom{\dagger}}_{b}c^{\phantom{\dagger}}_{t}\rangle+\langle
c^{\dagger}_{t}c^{\dagger}_{b}\rangle)/2$, and an inter-layer polarization
field denoted by $\phi_{z}$ which is $\phi_{z}=U\langle
c^{\dagger}_{t}c^{\phantom{\dagger}}_{t}\rangle$. Combining these, we obtain
the mean-field Hamiltonian

$\mathcal{H}_{MF}=\sum\limits_{{\bf k}}\psi^{\dagger}_{{\bf
k}}\Big{[}\Big{(}\dfrac{{\bf
k}^{2}}{2m}+V_{b}-\phi_{z}\Big{)}\tau_{z}+\phi\tau_{x}\Big{]}\psi^{\phantom{\dagger}}_{{\bf
k}}.$ (44)

where $\psi^{\phantom{\dagger}}_{\bf k}=(c^{\phantom{\dagger}}_{e,{\bf
k}},c^{\phantom{\dagger}}_{h,{\bf k}})^{T}$ (since electron/hole labels are
identical to the layer labels). Next, we define $\xi_{{\bf k}}={\bf
k}^{2}/2m+h-\phi_{z}$ and $E_{{\bf k}}=\sqrt{\xi_{{\bf k}}^{2}+\phi^{2}}$ to
write the self-consistent equations (at $T=0$) for $\phi$ and $\phi_{z}$

$\dfrac{1}{V}=\dfrac{1}{N}\sum\limits_{{\bf k}}\dfrac{1}{2E_{{\bf
k}}},\quad\phi_{z}=\dfrac{V}{2N}\sum\limits_{{\bf k}}\left(1-\dfrac{\xi_{{\bf
k}}}{E_{{\bf k}}}\right).$ (45)

The solutions can be obtained numerically with a lattice regularization ${\bf
k}^{2}/2m\rightarrow\epsilon_{{\bf k}}=2t(2-\cos k_{x}-\cos k_{y})$ where $t$
can be thought of as a nearest-neighbor hopping amplitude on a square lattice.

The multi-band generalization follows along similar lines. First we compute
the eigenvalues and eigenvectors $\\{E_{m}({\bf k}),|u_{m,{\bf k}}\rangle\\}$
of the mean-field Hamiltonian. The mean-field states follow the thermal
occupation

$\left\langle\gamma^{\dagger}_{{\bf k},m}\gamma^{\phantom{\dagger}}_{{\bf
k},m^{\prime}}\right\rangle=f[E_{m}({\bf k})]\delta_{m,m^{\prime}}$ (46)

where $f$ is the Fermi function and the operators $\\{\gamma_{{\bf k}}\\}$ are
related to $\\{c_{\bf k}\\}$ by a unitary transformation. We use the new
operators $\\{\gamma_{{\bf k}}\\}$ to rewrite Eq. (41) which turns out to have
a concise matrix structure

$\phi_{\alpha\beta}=-\dfrac{V_{\alpha\beta}}{N}\sum\limits_{{\bf
k}}f[E_{m}({\bf k})]\left\langle u_{m,{\bf
k}}\left|\dfrac{\partial\mathcal{H}_{MF}}{\partial\phi_{\alpha\beta}^{*}}\right|u_{m,{\bf
k}}\right\rangle,\quad
n_{\alpha,t}=\dfrac{V_{\alpha\beta}}{N}\sum\limits_{{\bf k}}f[E_{m}({\bf
k})]\left\langle u_{m,{\bf k}}\left|\dfrac{\partial\mathcal{H}_{MF}}{\partial
n_{\beta,b}}\right|u_{m,{\bf k}}\right\rangle.$ (47)

The equation is dependent on $\\{n_{\alpha,t},\phi_{\alpha\beta}\\}$ on both
sides, thus requiring an iterative approach. We start with a guess and use
FindRoot function in Mathematica to find the solution. The solutions are shown
in Fig. S1.

###  S-II.1 Tunable Quantum Metric Model

We reprint the Bloch hamiltonian for the tunable quantum metric model for
convenience

$H(\zeta,{\bf k})=2t(2-p_{\bf k})+t_{F}\cos(\zeta p_{\bf
k})\sigma^{x}+t_{F}\sin(\zeta p_{\bf k})\sigma^{y}.$ (48)

Here the Pauli matrices $\sigma^{i}$ act in orbital space, $p_{\bf k}$ is a
periodic function given by $p_{\bf k}=\cos k_{x}+\cos k_{y}$ and $\zeta$ is an
overall scaling factor for the quantum geometric tensor $\mathcal{Q}({\bf
k},\zeta)_{\mu\nu}=(\zeta^{2}/4)\sin k_{\mu}\sin k_{\nu}$.

We find a significant simplification with $\phi_{\alpha\beta}({\bf
k})=\phi\delta_{\alpha\beta}$. This relatively simple form can be attributed
to three key factors. The on-site nature makes momentum independent
$\hat{\phi}$ a natural ansatz. The orbital-diagonal aspect
$\hat{n}_{\alpha,t}\hat{n}_{\alpha,b}$ removes the possibility of off-diagonal
elements in $\hat{\phi}$. Lastly, the two orbitals in the toy model are
related by an in-plane inversion which ultimately reduces the exciton field to
be proportional to identity with a single scalar $\hat{\phi}=\phi\mathbb{1}$.

When both layers are described by the same parameter $\zeta$ (see Fig.2a), the
stiffness has a significant enhancement coming from the geometric term
$D_{s}^{g}$. The extra contribution is indeed proportional to the trace of the
quantum metric in this particular case $D_{s}^{g}=(\phi^{2}/4A)\sum_{\bf
k}g({\bf k})/E({\bf k})$ [69]. While the geometric contribution dominates at
higher bias voltage $|V_{b}|$, it is important to note that a high $|V_{b}|$
results in higher densities of electron-hole pairs and may ultimately drive
the system to the exciton Mott transition [45, 46, 47, 48].

Figure S1: Mean-field theory for symmetric excitons in the minimal model. a
Both the particle and hole bands have the same dispersion and wavefunction.
The non-interacting bands (dashed line) have a gap $V_{b}$ that increases
further with the formation of excitons (solid line). b Self-consistent
solution of $\phi$ together with $D_{s}$. Both require a finite interaction
strength $U/t>0$ to activate. c Geometric contribution to stiffness scales
with $\zeta^{2}$ as expected from trace of the quantum metric. d and e show
the temperature dependence of $\phi$ and $D_{s}$ respectively. The mean-field
temperature (at which $\phi$ vanishes) is clearly larger than the BKT
temperature (intersection of $D_{s}(T)$ and $2T/\pi$). The condensate is
limited by phase fluctuations.

##  Appendix S-III Excitons with Asymmetric Bands

In this section, we discuss various aspects of particle-hole asymmetry in the
geometric contribution to stiffness.

Let us begin by ignoring the effects of wavefunctions and recall the theory of
excitons with asymmetric electron and hole masses. This is a common scenario
in most semiconductors where the conduction and valence bands exhibit
different dispersions, characterized by $m_{e}\neq m_{h}$. It leads to a
slightly generalized form of Eq. (44) with the mean-field Hamiltonian:

$\mathcal{H}_{\rm MF}=\sum\limits_{{\bf k}}\psi^{\dagger}_{{\bf
k}}\Big{[}\dfrac{{\bf k}^{2}}{2M}\mathbb{1}+\Big{(}\dfrac{{\bf
k}^{2}}{2m}+V_{b}-\phi_{z}\Big{)}\tau_{z}+\phi\tau_{x}\Big{]}\psi^{\phantom{\dagger}}_{{\bf
k}}.$ (49)

Here, $M^{-1}=(m_{e}-m_{h})/2m_{e}m_{h}$ and
$m^{-1}=(m_{e}+m_{h})/2m_{e}m_{h}$ are effective masses. The inclusion of the
identity term modifies the mean-field dispersion relation $E_{\pm}({\bf
k})=\xi_{M}({\bf k})\pm E({\bf k})$, where $E({\bf k})=\sqrt{\xi({\bf
k})^{2}+\phi^{2}}$ represents the previous quasiparticle dispersion and
$\xi_{M}({\bf k})={\bf k}^{2}/2M$. The self-consistent order parameter
equation changes as well, mainly via the thermal factors, which were
suppressed in Eq. (45),

$\dfrac{1}{V}=\dfrac{1}{N}\sum\limits_{{\bf k}}\dfrac{f[E_{-}({\bf
k})]-f[E_{+}({\bf k})]}{2E({\bf k})}.$ (50)

It is easy to check that the equation reduces to Eq. (45) if $E_{+}({\bf
k})=-E_{-}({\bf k})=E({\bf k})$. The effects stemming from the asymmetry can
be visualized as a momentum-dependent Zeeman term in analogy to the one
encountered in BCS theory.  It is well known that large anisotropy can act
against exciton binding [83, 84] and destroy the condensate. Phase stiffness
is affected as well since it is a convolution of $\phi$ with the counterflow
current operator $\hat{j}_{\rm CF}$. If $\phi$ decreases, $D_{s}$ will
decrease as well.

Next, we consider the more subtle scenario of symmetric energies but
asymmetric wavefunctions. We again utilize the minimal model presented in the
previous section, eq. (48), with $\zeta_{t}=1$ for the top layer and
$\zeta_{b}=\zeta$ in the bottom layer. This configuration results in a
particle-hole symmetric spectrum but with distinct electron and hole
wavefunctions. The overlap $\theta({\bf k})=\cos[(\zeta-1)p_{\bf k}/2]$
exhibits significant variation across the Brillouin zone (see Fig. S2a). This
variation plays a crucial role in the quasiparticle energy $E({\bf
k})=\sqrt{\xi({\bf k})^{2}+\phi^{2}|\theta({\bf k})|^{2}}$, which, in turn,
affects the gap equation $1/V=(1/N)\sum_{\bf k}1/[E({\bf k})]$. The reduction
in $|\theta({\bf k})|$ from unity leads to a decrease in $\phi$ compared to
the particle-hole symmetric case where $\theta({\bf k})=1$ (see Fig. S2b and
Fig. S2d). The stiffness is also affected, as shown in Fig. S2c. Two factors
contribute to this alteration: pre-factor $\phi$ and $\mathcal{G}({\bf k})$
which explicitly depends on $\zeta$ through $\theta({\bf k})$,
$\mathcal{P}({\bf k})$, and $\mathcal{D}({\bf k})$. Although it may be model
specific, we find that the first contribution is symmetric around $\zeta=1$
while the geometric contibution is skewed towards higher $\zeta$ (see Fig.
S2c). Additionally, it is important to highlight that the geometric
contribution is not necessarily positive definite, see Fig. S2e. This stands
in contrast to the geometric contributions reported in superfluid stiffness
[2] and electron-phonon coupling [21].

Figure S2: Asymmetric excitons in the minimal model with tunable quantum
metric. a The overlap function between two bands defined by $\zeta_{t}=1$,
$\zeta_{b}=\zeta=2$. The gap equation (b) and stiffness (c) are both modified
as the asymmetry in wavefunction is tuned by $\zeta$ (with $V_{b}=0.1t$).
Notice that the conventional contribution is symmetric around
$\zeta_{t}=\zeta_{b}=1$. d Overlap function averaged over the BZ. Departure
from unity results in a smaller mean-field solution for $\phi$. The isotropic
case $\zeta_{t}=\zeta_{b}$ results in the strongest exciton. e Sign of the
geometric contribution in the minimal model. There are parameter regions where
the contribution is negative. Intuitively, quantum geometry can be detrimental
to stiffness if the bands are very different.

Motivated by these two simple scenarios, we can deduce what happens in a
realistic situation with both asymmetric dispersions and wavefunctions. The
dispersion asymmetry results in redefinition of effected masses $m$ and $M$.
On top of that, the wavefunction asymmetry modifies both $\phi$ and
$\mathcal{G}({\bf k})$. Both these effects are at play simultaneously, the
details of which depend on the model and we are left limited to numerical
results.

##  Appendix S-IV Continuum model for exciton condensation in AA-stacked TMDs

To start with we consider spinless electrons described by the effective
Hamiltonian:

$H=\sum^{\Lambda}_{{\bf k},a}{\psi}^{\dagger}_{{\bf
k},a}\begin{bmatrix}H_{t,a}({\bf k})-\mu+E_{z}&0\\\ 0&H_{b,a}({\bf k})-\mu-
E_{z}\end{bmatrix}{\psi}_{{\bf k},a}+\frac{1}{2A}\sum^{\Lambda}_{{\bf
q}}\sum_{\nu,\nu^{\prime}=t,b}V^{\nu\nu^{\prime}}_{\bf q}n_{{\bf
q},\nu}n_{-{\bf q},\nu^{\prime}},$ (51)

where $\psi_{k,a}$ is a 4-dim spinor in the sublattice and layer indeces,
$\psi_{k,a}=[\psi_{k,a+t},\psi_{k,a-t},\psi_{k,a+b},\psi_{k,a-b}]$, $K$ and
$K^{\prime}$ refers to the two different valleys, $\Lambda$ is the UV cutoff
and $E_{z}$ the electric displacement field controlling the energy offset
between the two layers. The interaction is obtained by taking the Fourier
transform of $V^{tt}_{r}=V^{bb}_{r}=e^{2}/(4\pi\epsilon_{0}\epsilon r)$ and
$V^{tb}_{r}=V^{bt}_{r}=e^{2}/(4\pi\epsilon_{0}\epsilon\sqrt{r^{2}+d^{2}})$
with $d$ interlayer distance. From

$\int\frac{d^{2}{\bf q}}{2\pi}\frac{e^{-dq+i{\bf q}\cdot{\bf
r}}}{q}=\oint_{|z|=1}\frac{2dz}{2\pi
r}\frac{1}{z^{2}+1-2dz/(ir)}=\frac{1}{\sqrt{r^{2}+d^{2}}},$ (52)

we readily find:

$V^{tt}_{\bf q}=V^{bb}_{\bf q}=V_{\bf q}=\frac{e^{2}}{2\epsilon_{0}\epsilon
q},\quad V^{tb}_{\bf q}=e^{-qd}V_{\bf q}.$ (53)

Accounting for the presence of the metallic gates at distance $\xi$ from the
sample we replace $1/q\to\tanh(q\xi/2)/q$, which regularizes the singularity
at $q=0$. We use $\sigma$ for $d_{z^{2}}/d_{{x^{2}-y^{2}}}\pm id_{{xy}}$ at
valley $K/K^{\prime}$ while $\tau$ for layer. The chemical potential $\mu$
enforces the charge neutrality condition:

$\displaystyle\frac{1}{N}\sum^{\Lambda}_{\bf k}\langle\psi^{\dagger}_{{\bf
k},a}\psi_{{\bf k},a}\rangle-n_{\rm CNP}=0.$ (54)

####  S-IV.0.1 Model for TMDs

We assume that due to the strong Ising spin orbit coupling the conduction (CB)
and valence band (VB) of each layer around $K$ and $K^{\prime}$ points are
spin polarized. This assumption is realized under the assumption of large
Ising spin-orbit coupling $\tau\Delta_{vb}s^{z}$ and $\tau\Delta_{cb}s^{z}$
with $\tau=\pm$ for the two valleys. For MX2 the SOC gap is
$2\Delta_{vb}=148,186,219$meV with
X=$\\{\text{S}_{2},\text{Se}_{2},\text{Te}_{2}\\}$. For WX2 the gap is larger
due to the larger SOC $2\Delta_{vb}=429,466,484$meV for VB and
$2\Delta_{cb}=-32,-37,-53$meV for CB with
X=$\\{\text{S}_{2},\text{Se}_{2},\text{Te}_{2}\\}$. In the conduction band the
SOC is small and corrections from the nearby conduction band might play an
important role. Keeping terms up to the second order in momentum the $k.p$
Hamiltonian reads:

$H_{\nu,K}({\bf
k})=\begin{bmatrix}\Delta_{\nu}+\alpha_{\nu}k^{2}&v_{\nu}k_{-}\\\
v_{\nu}k_{+}&-m_{\nu}-\gamma_{\nu}k^{2}\end{bmatrix},\quad
H_{t/b,K^{\prime}}({\bf k})=H^{*}_{t/b,K}(-{\bf k}),$ (55)

where the latter relation enforces the time reversal symmetry. The typical
values of the parameters of the model are given in Table 1. Different
configurations are characterized by a typical value of the ratio $v/m\approx
2\text{\AA}$. Finally, in the main text we compare results obtained with
gapped Dirac cones and parabolic bands. The latter is obtained expanding the
Dirac dispersion around $K$ and $K^{\prime}$,
$\pm\sqrt{\Delta^{2}+v^{2}k^{2}}\simeq\pm\Delta\pm k^{2}/2m$ with
$m=\Delta/v^{2}$ and gap $2\Delta$.

| MoS2 | MoSe2 | MoTe2 | WS2 | WSe2 | WTe2
---|---|---|---|---|---|---
$2m$[eV] | 2.76 | 2.33 | 1.82 | 2.88 | 2.42 | 1.82
$v$[eV Å] | 2.22 | 2.20 | 2.16 | 2.59 | 2.60 | 2.79
$\alpha$[eV Å2] | 0.52 | -0.54 | -1.66 | 2.03 | 1.08 | 0.17
$\gamma$[eV Å2] | -6.65 | -6.20 | -5.78 | -7.96 | -7.58 | -7.31
Table 1: Relevant parameters of the $k.p$ Hamiltonian for different TMDs [72]
computed performing GW calculations [85, 86].

####  S-IV.0.2 Symmetries

The model is invariant under the transformation:

$\psi_{{\bf k},at}\to e^{i\phi_{t}}\psi_{{\bf k},at},\quad\psi_{{\bf k},at}\to
e^{i\phi_{b}}\psi_{{\bf k},at},$ (56)

where $\phi_{b}+\phi_{t}$ is the total phase associate to U(1)charge while the
relative phase $\phi_{t}-\phi_{b}$ is related to the U(1)layer symmetry. The
latter is spontaneously broken in the excitonic state which is characterized
by interlayer coherence. In addition, we observe that the model in Eq. 55
breaks the particle-hole symmetry for $\alpha\neq\gamma$. The low-energy model
is also invariant under the continuous space rotational rotational symmetry
${\bf k}\to R_{\theta}{\bf k}$ and
$U_{\theta}=\text{diag}\left(e^{-i\theta/2},e^{i\theta/2}\right)$ which acts
on the sublattice degree of freedom:

$U^{\dagger}_{\theta}H_{\nu,K/K^{\prime}}(R_{\theta}{\bf
k})U_{\theta}=H_{\nu,K/K^{\prime}}({\bf k}).$ (57)

Including higher order terms in the $k.p$ expansion lower the aforementioned
symmetry to $C_{3z}$. Finally, $\mathcal{M}$ acting as $(x,y)\to(-x,y)$ and
$\uparrow\leftrightarrow\downarrow$ is a symmetry of the model.

####  S-IV.0.3 Diagonalizing the single-particle Hamiltonian

The eigenvalues of the non-interacting Hamiltonian are given by
$E_{\pm,t/b}(k)=\pm\epsilon_{t/b}(k)$,
$\epsilon_{t/b}(k)=\sqrt{v^{2}_{t/b}k^{2}+\Delta_{t/b}^{2}}$ while the matrix
of eigenstates is:

$U({\bf
k})=\frac{1}{\sqrt{2}}\begin{bmatrix}\cos\frac{\varphi_{k}}{2}+\sin\frac{\varphi_{k}}{2}&\cos\frac{\varphi_{k}}{2}-\sin\frac{\varphi_{k}}{2}\\\
e^{i\phi_{k}}\left(\cos\frac{\varphi_{k}}{2}-\sin\frac{\varphi_{k}}{2}\right)&-e^{i\phi_{k}}\left(\cos\frac{\varphi_{k}}{2}+\sin\frac{\varphi_{k}}{2}\right)\end{bmatrix},$
(58)

where $\tan\varphi_{k}=\Delta/(vk)$ and we dropped the layer and valley
labels. Applying the unitary transformation:

$U({\bf k})=\begin{bmatrix}U_{t}({\bf k})&0\\\ 0&U_{b}({\bf k})\end{bmatrix},$
(59)

the kinetic part becomes:

$H_{\rm kin}=\sum_{a=K,K^{\prime}}\sum^{\Lambda}_{\bf k}\Psi^{\dagger}_{{\bf
k},a}\begin{bmatrix}\varepsilon_{t}(k)\sigma^{z}+(V_{b}-\mu)\sigma^{0}&0\\\
0&\varepsilon_{b}(k)\sigma^{z}+(-V_{b}-\mu)\sigma^{0}\end{bmatrix}\Psi_{{\bf
k},a},$ (60)

where $\Psi_{{\bf k},a\nu}=U^{\dagger}_{a\nu}(k)\psi_{{\bf k},a\nu}$ with
$a=K/K^{\prime}$ and $\nu=t,b$. We now look at the interaction projected in
the basis of the eigenstates of the kinetic term. The density operator for the
layer $\nu=t/b$ reads:

$\rho^{\nu}_{\bf q}=\sum_{\bf k}\psi^{\dagger}_{{\bf k}\nu}\psi_{{\bf k}+{\bf
q}\nu}=\sum_{\bf k}\sum_{\alpha\alpha^{\prime}}\left[F^{\nu}_{{\bf k},{\bf
q}}\right]_{\alpha\alpha^{\prime}}\Psi^{\dagger}_{{\bf k}\nu\alpha}\Psi_{{\bf
k}+{\bf q}\nu\alpha^{\prime}},$ (61)

where we drop for the moment the valley index, and the form factor is obtained
projecting the density operator in the Bloch basis

$F^{\nu}_{{\bf k},{\bf k}}=U^{\dagger}_{\nu}({\bf k})U_{\nu}({\bf k}+{\bf
q}),\quad\left[F^{\nu}_{{\bf k},{\bf
q}}\right]_{\alpha\alpha^{\prime}}=\innerproduct{u^{\nu}_{{\bf
k},\alpha}}{u^{\nu}_{{\bf k}+{\bf q},\alpha^{\prime}}}.$ (62)

The orthogonality condition implies that $\left[F^{\nu}_{{\bf
k},0}\right]_{\alpha\alpha^{\prime}}=\delta_{\alpha\alpha^{\prime}}$. As a
result the interaction takes the form:

$H_{\rm int}=\frac{1}{2A}\sum^{\Lambda}_{{\bf q}{\bf k}{\bf
k}^{\prime}}\sum_{\nu\nu^{\prime}}\sum_{aa^{\prime}}\left[F^{\nu,a}_{{\bf
k},{\bf q}}\right]_{\alpha\alpha^{\prime}}V^{\nu\nu^{\prime}}_{\bf
q}\left[F^{\nu^{\prime},a^{\prime}\dagger}_{{\bf k}^{\prime},{\bf
q}}\right]_{\gamma\gamma^{\prime}}\Psi^{\dagger}_{{\bf
k},a\nu\alpha}\Psi_{{\bf k}+{\bf q},a\nu\alpha^{\prime}}\Psi^{\dagger}_{{\bf
k}^{\prime}+{\bf q},a^{\prime}\nu^{\prime}\gamma}\Psi_{{\bf
k}^{\prime},a^{\prime}\nu^{\prime}\gamma^{\prime}},$ (63)

At mean-field level we decouple the interaction in the Hartree and Fock
(exchange) terms that are derived in the following sections.

###  S-IV.1 Hartree term

The Hartree term corresponds to:

$H_{H}=\sum^{\Lambda}_{\mathbf{k}}\sum_{a=K,K^{\prime}}\sum_{\nu=t,b}\delta\mu_{\nu}\,\Psi^{\dagger}_{{\mathbf{k}},a\nu}\Psi_{{\mathbf{k}},a\nu}-\frac{1}{2A}\sum_{\nu\nu^{\prime}}V^{\nu\nu^{\prime}}_{0}\sum^{\Lambda}_{{\mathbf{k}}{\mathbf{k}}^{\prime}}\sum_{aa^{\prime}=K,K^{\prime}}\delta\rho^{a\nu}({\mathbf{k}})\delta\rho^{a^{\prime}\nu^{\prime}}({\mathbf{k}}^{\prime}),$
(64)

where we have introduced the charge density measured with respect to charge
neutrality:

$\delta\rho^{a\nu}({\mathbf{k}})=\sum_{\alpha=\pm}\langle\Psi^{\dagger}_{{\mathbf{k}},a\nu\alpha}\Psi_{{\mathbf{k}},a\nu\alpha}\rangle-N\delta_{{\mathbf{k}},0}.$
(65)

the contribution $\delta\mu_{\nu}$ is given by:

$\delta\mu_{\nu}=\frac{1}{N}\sum^{\Lambda}_{k}\sum_{a=K,K^{\prime}}\sum_{\alpha=\pm}\sum_{\nu^{\prime}}\delta\rho^{a\nu^{\prime}}_{\alpha\alpha}({\mathbf{k}})V^{\nu^{\prime}\nu}_{0}.$
(66)

Before moving on we observe that

$\frac{1}{N}\sum^{\Lambda}_{k}\sum_{a=K,K^{\prime}}\sum_{\alpha=\pm}\delta\rho^{a}_{\alpha\alpha}({\mathbf{k}})=n_{\nu}-1=\text{sign}(\nu)\phi_{z},$
(67)

where $\phi_{z}$ measure the charge density transferred between the layers
$\phi_{z}=(n_{t}-n_{b})/2=n_{t}-1=1-n_{b}$ since at charge neutrality
$n_{t}+n_{b}=2$ with this definition $\text{sign}(t)=1$ and
$\text{sign}(b)=-1$. Thus, we have:

$\delta\mu_{t}=\phi_{z}(V^{tt}_{0}-V^{bt}_{0}),\quad\delta\mu_{b}=\phi_{z}(V^{tb}_{0}-V^{bb}_{0}).$
(68)

From Eq. (53) we find:

$\delta\mu_{t}=\delta\mu,\quad\delta\mu_{b}=-\delta\mu,\quad\delta\mu=\phi_{z}(V^{tt}_{0}-V^{bt}_{0}),\quad\phi_{z}=(n_{t}-n_{b})/2,$
(69)

which corresponds to electrostatic repulsion. We notice that
$V^{tt}_{0}=V^{tb}_{0}$ and $\delta\mu=0$.

###  S-IV.2 Fock term

The exchange term is given by:

$H_{X}=-\sum^{\Lambda}_{k}\sum_{a=K,K^{\prime}}\sum_{\nu\nu^{\prime}}\Psi^{\dagger}_{{\mathbf{k}},a\nu\alpha}\Sigma^{\nu\nu^{\prime},a}_{F,\alpha\beta}({\mathbf{k}})\Psi_{{\mathbf{k}},a\nu^{\prime}\beta},$
(70)

where we have introduced the quantity:

$\Sigma^{\nu\nu^{\prime},a}_{X,\alpha\beta}({\mathbf{k}})=\frac{1}{A}\sum_{{\mathbf{k}}^{\prime}}V^{\nu\nu^{\prime}}_{{\mathbf{k}}^{\prime}-{\mathbf{k}}}\sum_{\alpha^{\prime}\beta^{\prime}}\left[F^{\nu^{\prime},a\dagger}_{{\mathbf{k}},{\mathbf{k}}^{\prime}-{\mathbf{k}}}\right]_{\beta^{\prime}\beta}\left[F^{\nu,a}_{{\mathbf{k}},{\mathbf{k}}^{\prime}-{\mathbf{k}}}\right]_{\alpha\alpha^{\prime}}\delta\rho^{\nu^{\prime}\nu,a}_{\beta^{\prime}\alpha^{\prime}}({\mathbf{k}}^{\prime}),$
(71)

where we have introduced the quantity:

$\delta\rho^{\nu^{\prime}\nu,a}_{\beta^{\prime}\alpha^{\prime}}({\mathbf{k}})=\langle\Psi^{\dagger}_{{\mathbf{k}},a\nu^{\prime}\beta^{\prime}}\Psi_{{\mathbf{k}},a\nu\alpha^{\prime}}\rangle-\delta_{\nu,\nu^{\prime}}\delta_{\alpha^{\prime}\beta^{\prime}}\delta_{\alpha^{\prime}-}.$
(72)

$\Sigma({\mathbf{k}})$ is a 4-dim matrix in the layer and band basis. It
includes an intralayer terms which renormalize the electronic band structure
and an interlayer term giving rise to interlayer phase coherence and exciton
condensation. The Hartree-Fock Hamiltonian takes the form:

$H_{\rm HF}=\sum^{\Lambda}_{k}\Psi^{\dagger}_{{\mathbf{k}}}\left[\hat{H}_{\rm
kin}({\mathbf{k}})-\hat{\Sigma}({\mathbf{k}})\right]\Psi_{{\mathbf{k}}}.$ (73)

The self-consistency equation 71 for the self-energy is solved at every
${\mathbf{k}}$-point in momentum space employing an iterative approach. Fig.
S3 shows the typical band structure and the interlayer gap function, defined
as $\sqrt{\Tr\Delta^{\dagger}({{\mathbf{k}}})\Delta({{\mathbf{k}}})}$ with
$\hat{\Delta}({\mathbf{k}})$ interlayer component of the self-energy
$\hat{\Delta}({\mathbf{k}})=\text{Tr}\left[\tau^{+}\hat{\Sigma}({\mathbf{k}})\right]$.
Finally, we remark that the excitonic order parameter $\phi$ has been computed
as:

$\phi_{x}=\frac{T}{A}\sum_{{\mathbf{k}}}\sum_{i\omega}\text{Tr}\left[\tau^{x}\otimes\sigma^{a}\hat{G}({\mathbf{k}},i\omega)\right],$
(74)

in the gauge where the interlayer coherence is broken along $\tau^{x}$ and the
Green’s function $\hat{G}$ is obtained including all the bands of the model
from Eq. (73).

Figure S3:  a) Non-interacting band structure (dashed) and interacting
excitonic energy (solid). The color denotes the layer polarization showing the
mixing between the different bands in the excitonic state. The inset shows the
intralayer and interlayer Coulomb interaction as a function of $q/q_{\Lambda}$
in unit of $V_{\bf q=0}$. Notice that $V^{tb}_{\bf q}$ decays exponentially
with $q$. b) Gap function $\sqrt{\Tr\Delta^{\dagger}({\bf k})\Delta({\bf k})}$
denoting interlayer phase coherence. Calculations are performed for $\Delta
E=(\Delta_{t}+\Delta_{b}-2V_{b})/2=55$meV where the single particle spectrum
is gapped, $v_{t}=v_{b}=2.19$eVÅ, $\Delta_{t}=\Delta_{b}=0.9$eV, we set the
screening length $\xi=12$nm and $d=$1nm. We fix the UV cutoff
$q_{\Lambda}a=2\pi/3$ corresponding to half-distance between the two corners
$K$ and $K^{\prime}$ of the TMD Brillouin zone. These parameters are the same
ad those employed for the results on homobilayer reported in the maintext.