Bilayer WSe₂ as a natural platform for interlayer exciton condensates in the strong coupling limit

The heterostructure is assembled using the van der Waals dry transfer technique, and pre-patterned Pt electrodes are used for electrical contacts to bilayer WSe₂. Hexagonal boron nitride are used as the dielectric, and graphite or metal as top and bottom gates for bilayer WSe₂. The top gate was lithographically shaped so that its overlap with bottom gate covers only bilayer WSe₂, and the overlap defines the device area. In order to achieve good electrical contact, we use an additional contact gate on top to heavily dope the contact area, which is isolated from the top gate by Al₂O₃ dielectric. Data from a different device are shown in SI.9. Penetration capacitance was measured with an FHX35X high electron mobility transistor serving as a low-temperature amplifier, in a similar setup as in Ref. Shi et al., 2019. Measurements were performed at $T=0.3$ K unless otherwise specified.

Here we discuss the suppression of interlayer tunneling in the context of conventional double quantum wells. In conventional GaAs double quantum wells, the interlayer tunneling is usually characterized by $\Delta_{\text{SAS}}$. The two single-particle eigenstates can be written as,

Here, $\phi_{t}$ and $\phi_{b}$ are the wavefunctions localized in the top and bottom quantum wells, respectively;

characterize the degree of wavefunction localization, where $V$ is the interlayer electrical potential difference.

For any finite tunneling gap $\Delta_{\text{SAS}}$, $\cos 2\alpha=0$ at zero bias and the single-particle wavefunctions are just the bounding and anti-bounding states that is an superposition of wavefunctions in both layers with equal weight ($\cos\alpha=\sin\alpha$). The single-particle ground state has a gap $\Delta_{\text{SAS}}$ to the excited state.

In contrast, the bilayer WSe₂ manifest a large degree of wavefunction localization in individual layers even at zero bias, and this is a result of the strong spin-orbit coupling. It was shown that the spin doublet single-particle eigenstates in the valence band can be written as Gong et al. (2013),

where

is the wavefunction localized in top and bottom layers and $\tau$ is the valley index (1 for K and -1 for K’),

$\lambda$ is the spin-orbit coupling strength and $t$ is the interlayer hopping amplitude.

For holes in bilayer WSe₂, $\cos 2\beta\approx 0.96$ according to ab initio calculations Gong et al. (2013). Therefore, the single-particle eigenstate wavefunctions are always highly localized in individual layers. The spin-orbital coupling strength $\lambda$ is effectively an intrinsic bias and only the low energy branches are accessible in our experiment. All four flavors described in equations above are degenerate in energy in the absence of external fields, and no single-particle gap is expected. Such energetically degenerate yet spatially separated single-particle wavefunctions set the basis for exciton condensate physics.

The above discussion suggests a high level of layer localization which persists to the LL regime. We therefore treat the LL structures in the two layers independently. In a monolayer WSe₂, the LL structure has two isospin flavors, separated by a large spin-splitting energy. The spin susceptibility, characterized by $E_{Z}/E_{C}$, where $E_{Z}$ is the spin-splitting energy and $E_{C}$ the cyclotron energy, increases as carrier density decreases and could be as large as 9 Gustafsson et al. (2018); Shi et al. (2019). In bilayer WSe₂, the same physics applies, but due to a different dielectric environment with larger dielectric constant, the spin susceptibility is slightly smaller. Further, the susceptibility in one layer mainly depends on the carrier density of this layer, and slightly on the density of the other layer. In our experiment, the spin of a certain LL is decided from similar analysis as in Shi et al. (2019), and assisted with transport measurements where different spin flavors exhibits very different conductivity Lin et al. (2019). At 15 T, five spin-up LLs are populated before the first spin-down LL is populated in each layer.

It is well known that for the valence band of WSe₂, the LL with orbital index $N=0$ only exists in one valley; in the other valley, the lowest LL has orbital index $N=1$ Rose et al. (2013). However, it is worth noting that the LL wavefunctions is not related to the orbital index in the simple way as in non-Dirac systems, but a mixture of the harmonic oscillator wavefunctions. Indeed, in the valley with LL orbital index $N=0$, the wavefunction is pure $|0\rangle$, while the next LL with $N=1$ is mainly $|1\rangle$ with a few percent of mixing of $|0\rangle$, where $|n\rangle$ stands for the harmonic oscillator wavefunction of index $n$. On the other hand, in the valley without a $N=0$ LL, its lowest LL with $N=1$ actually has the wavefunction dominated by $|0\rangle$ and only with a few percent of mixing from $|1\rangle$. Therefore, to the first order, the lowest LLs in the two valleys have similar wavefunctions despite of different LL orbital index. This is important to establish interlayer EC in the case of mismatched spins. In the main text, the orbital number $n$ stands for the index associated with the dominant component $|n\rangle$ of the single-particle wavefunction.

The strong correlation effects in TMDs result in a negative compressibility which is stronger at lower densities Pisoni et al. (2019). This effect manifests in several aspects of our data. In Fig.1 (d) and (f) of the main text we outline the polarization of the bilayer at zero field. Here the boundary of the fully polarized region is curved. In a simple electrostatic model,

where $\epsilon$ is the dielectric constant and $E$ is the unscreened interlayer electrical field. Without interaction effects considered, $E$ would not have any dependence on the total density, and the displacement field required to fully polarize the bilayer should increase linearly with the total density. Due to the negative compressibility, $E$ is of opposite sign to $D$, and its amplitude increases with decreasing $n$ - the bilayer electrons tend to build a larger density imbalance to“overscreen” the displacement field. Therefore, the displacement field required to fully polarize the bilayer increases superlinearly with the total density, giving rise to the curve shaped boundary in Fig.1 of main text.

Similarly, the negative compressibility effect also manifest in the data under a perpendicular magnetic field. Under a magnetic field, a single-particle picture would predict that the LL crossings happen at the same $D$ for a certain filling factor difference $\Delta\nu=\nu_{B}-\nu_{T}$, as shown in Fig.3(b) in the main text. In contrast, as shown in Fig.3(a) in the main text, for a fixed $\Delta\nu$, the LL crossings happen at a smaller $D$ for smaller $n$, and the polarizability segments at partial fillings also assume a slope.

The penetration capacitance is a sum of both “incompressibility” and “polarizability” contributions. First, we confirm that the polarizability contribution is only in the capacitance channel, and is not dissipation related. In Fig. 5 we plot the capacitive and dissipative signal in the measured penetration capacitance in (a) and (b), respectively. While most integer quantum Hall gaps show up in both, the features cutting the integer gap only appear in the capacitive channel, which are consistent with the polarizability contribution.

We further distinguish the two contributions according to their behavior at different filling factors. Fig. 6 (b) shows the schematic LL diagram corresponding to the data in Fig. 6 (a). In Fig. 6(c) we plot the linecut of $C_{P}$ along constant non-integer filling factor $\nu_{\text{tot}}=3.7$ [horizontal line in Fig. 6(a)]. Here, the Fermi level changes back and forth between two sets of LLs, as marked by the thick lines in Fig. 6(b). The bilayer system is always compressible and the peaks arises from the polarizability contributions . Indeed, whenever the Fermi level changes between the LLs localized in different layers, interlayer charge transfer happens, and gives rise to a peak in $C_{P}$. We observe a total of seven peaks in the whole $D$ range in Fig. 6(c), separating red and blue regions, which correspond to the cases that the Fermi level is in the bottom (red lines) and top (blue lines) layers, respectively. Further, a small peak corresponds to charge transfer amount of 0.3 of the Landau level degeneracy, while a large peak corresponds to 0.7, marked by the numbers in Fig. 6(c).

In Fig. 6(d) we plot the linecut of $C_{P}$ along constant integer filling factor $\nu_{\text{tot}}=4$. We observe that $C_{P}$ show four dips in the whole $D$ range, which corresponds to four LL crossings, separating regions of different filling factors marked in Fig. 6(d). At the level crossings, in principle there is “polarizability” contribution as well, but it is much smaller, as evidenced from different vertical scales for Fig. 6(c) and (d).

To distinguish the two contributions at LL crossings at integer filling factors, we make a linecut along the polarizability segments, cutting through integer filling lines, as shown in Fig. 6(e). Such a linecut represents the condition $\nu_{T}^{*}=\nu_{B}^{*}$, where $\nu_{T}^{*}$ ($\nu_{B}^{*}$) is the fractional part of the filling factor in the top (bottom) layer. While the polarizability contribution continuously changes from non-integer $\nu_{\text{tot}}$ to integer $\nu_{\text{tot}}$, an incompressible contribution would only be possible at integer $\nu_{\text{tot}}$, and shows up as a sharp peak at $\nu_{T}^{*}=\nu_{B}^{*}=0.5$. The two panels in Fig. 6(e) demonstrate the two situations, respectively. In the right panel, an EC gap exist as orbital numbers match in the two layers; in the left panel, an EC gap does not exist as the orbitals are mismatched.

In the presence of EC, the linecut along constant $\nu_{tot}$ is shown in Fig. 6(f). Here, while $C_{P}$ at other LL crossings demonstrate simply a dip at the closing of IQH gap, its behavior in the EC (grey shaded) region is different, exhibiting an “reentrant” increase of $C_{P}$, which corresponds to the emergence of a correlation gap.

In Fig. 7 we plot $C_{P}/C_{G}$ for LL crossings with the same orbital number but different spin. For $n=0$, as shown in Fig. 7 (a), no incompressible state is observed when both layers host partially filled LLs. In contrast, for $n=2$, as shown in Fig. 7 (b), incompressible state persists when both layers are partially filled, consistent with EC formation. The suppression of EC at $n=0$ is consistent with what we observe for the LL crossings between the same spins, at $D=0$, demonstrated in Fig.3 of the main text.

In Fig. 8 (a) and (b) we show the capacitance as a function of top and bottom gate voltages at 14.9 T and 30 T, focusing at $\nu_{tot}=1$. The charge gap at $\nu_{tot}=1$ disappears at the layer balanced point $D=0$ at 14.9 T, while it recovers at 30 T. The difference could be due to a larger interaction energy scale at 30 T, which enhances the EC gap that overcomes disorder effects.

In Fig. 9(a) we plot the same data as in Fig.3(c) in the main text, which shows $C_{P}$ at density balance as a function of the magnetic field, yet here versus the density and in a larger filling factor range. The transition at $\nu_{T}=\nu_{B}=5.5$ is clearly distinct from the behavior at other half-integer filling factors. This transition originates from a density dependent $g$-factor, which varies the ratio of spin-splitting energy $E_{z}$ and the cyclotron energy $E_{c}$ Gustafsson et al. (2018); Shi et al. (2019) and shuffles LLs of different orbital and spin index as density is varied.

For fixed excitation amplitude and frequency, the measured penetration signal $V_{P}$, to a good approximation, can be written as Dultz and Jiang (2000); Goodall et al. (1985),

where the imaginary and real part, denoted $V_{P}$ and $V_{P}^{disp}$ stands for capacitive and dissipative contribution, respectively. In the equation above, $C_{t}$ ($C_{b}$) is the geometry capacitance between WSe₂ and the top (bottom) gate, $C_{G}=C_{t}C_{b}/(C_{t}+C_{b})$ is the geometry capacitance between the two gates, $C_{q}=Ae^{2}(\partial\mu/\partial n)^{-1}$ is the quantum capacitance of the WSe₂, with $A$ being the area of the device, and $R_{s}$ is the resistance of the WSe₂ channel.

In the low frequency and low resistance limit ($\alpha\to 0)$, the real part - dissipative contribution is zero, and the penetration capacitance $C_{P}$ can be written as,

in the high resistance limit ($\alpha\to\infty)$ or when the WSe₂ is gapped, the penetration capacitance is equal to the geometry capacitance between the two gates,

and the dissipative contribution is also zero.

In our experiment, the resistance at EC states are large, giving rise to finite dissipative contribution. To obtain the gap value, we obtain the quantum capacitance $C_{q}$ and resistance $R_{s}$ by numerically solving Eq. (13) and Eq. (14). An example of the obtained values $C_{q}$ and $R_{s}$ are shown in Fig. 10. The chemical potential jump due to the formation of incompressible states was then found by integrating $\partial\mu/\partial n$ over the peaks in the EC region. Such gap values are plotted in Fig.4(a) of the main text. In addition, we fit the temperature dependence of $R_{s}$ at the EC state by an activation behavior $R_{s}\sim\exp(-\Delta/2T)$. Such obtained gaps are plotted in Fig. 11 as circle markers. The quantitative differences in the gap values may be due to the contact resistance, which are not accounted for in our analysis. Nevertheless, they demonstrate the same trend, with a maximum at $n=2$.

EC with similar behavior as discussed in the main text have been observed in several devices. Here we show the data in a different device. In Fig. 12(a) we display $C_{P}/C_{G}$ versus the top and bottom gate voltages. In this device, the top gate also covers the contact region, so $V_{T}$ is kept at a large negative value in order to highly dope the contact region and maintain good electrical contact. As a result, the top layer always host a large filling factor. For $V_{B}\gtrsim-1.5$ V, the bottom layer is empty, and only one set of LLs gaps from the top layer are observed; for $V_{B}\lesssim-1.5$ V, both layers are populated, and the polarizability features appear, cutting the incompressible features along integer $\nu$ to discontinuous segments. It can be seen that the layer transitions are sharp when the bottom layer filling $\nu_{B}$ is low; that is, the carriers in the topmost LL transfers from the top to bottom layer within a small voltage change. The sharp segments smears out to a diamond shape as $V_{B}$ and the filling factor difference $\Delta\nu=\nu_{T}-\nu_{B}$ decreases. The white circles mark where the LLs in the two layers host the same orbital index, which corresponds to n = 2, 3, 4 and 5 respectively. Within these “diamonds”, we observe an incompressible feature when $\nu$ is integer while $\nu_{B}$ and $\nu_{T}$ are both fractions. Fig. 12(b) shows a highlight for $n=4$, with the LL structure sketched in Fig. 12(c). In contrast, the incompressible features are not present in the cases where the LL orbital indices are mismatched across the two layers.

In addition, in Fig. 13 we plot $C_{P}$ as a function of the bottom gate voltage $V_{B}$ and $B$ with the top gate voltage $V_{T}$ fixed. Three fans of high penetration signal are revealed. First, a main fan tracks constant integer fillings of $\nu=\nu_{T}+\nu_{B}$ and emerges from the band edge. For $V_{B}<1.5$ V, as both layers are populated, this fan is cut by layer transition features and becomes discontinuous. Second, a side fan emerges at the population onset of bottom layer, $V_{B}=1.5$ V, indicated by dotted lines. Third, a fan tracks constant $\Delta\nu=\nu_{T}-\nu_{B}$, marked by dashed lines. While the first two fans are both indications incompressible states at integer LL filling, the third fan is an indication of charge transfer between the layers. The black circles mark the position where the top most LLs in the two layers host the same LL orbital index, and an incompressible feature appear which signals the emergence of EC.

Our theoretical analysis of bilayer $\mathrm{WSe_{2}}$ is based on a simplified model that comprises two copies of a two-dimensional electron gas (2DEG), separated by distance $d$ in the perpendicular $z$-direction. The two layers are threaded by the magnetic field $B\hat{z}$ and they can have arbitrary Landau level (LL) indices, $n_{1}$ and $n_{2}$ (for simplicity, we neglect the spin degree of freedom). The case of $n_{1}=n_{2}=0$ and $d/\ell_{B}\gtrsim 1$ has been studied extensively in previous work on GaAs electron bilayers Chakraborty and Pietiläinen (1987); Moon et al. (1995); Scarola and Jain (2001); Schliemann et al. (2001); Sheng et al. (2003); Park (2004); Shibata and Yoshioka (2006); Möller et al. (2009); Milovanović et al. (2015); Zhu et al. (2017); Lian and Zhang (2018). Although the single-particle orbitals in $\mathrm{WSe_{2}}$ are spinors, the admixture of different LL orbitals is strongly suppressed compared to other materials such as graphene Rose et al. (2013).

In momentum space, our model is described by the Hamiltonian

where $\sigma=\uparrow,\downarrow$ labels the two layers, $N_{\phi}$ is the number of magnetic flux quanta through the system, $\bar{\rho}_{\sigma}$ is the projected electron density in a layer $\sigma$. $V_{\sigma\sigma^{\prime}}$ is the Fourier transform of the Coulomb interaction which is layer-dependent:

where $d$ is the distance between the layers which “softens” the repulsion between $\uparrow$ and $\downarrow$ electrons since $V_{\uparrow\downarrow}$ corresponds to interaction $1/\sqrt{r^{2}+d^{2}}$ in real space. Thus, $\uparrow$ and $\downarrow$ are sometimes referred to as “pseudospin” degrees of freedom Moon et al. (1995) to distinguish them from real electron spin. Furthermore, the effective interaction is also dependent on which LLs the layers are in, which modifies the projected density operators,

Here $\mathbf{R}_{j\sigma}$ stands for the guiding center coordinate of $j$th electron with the layer index $\sigma$ Prange and Girvin (1987), and $L_{n_{\sigma}}$ is the $n$th Laguerre polynomial for the layer $\sigma$ Macdonald (1994).

Below we work in units $e^{2}/\epsilon\ell_{B}=1$. We focus on the total electron filling factor $\nu=1$, corresponding to the total number of electrons being equal to the number of flux quanta, $N=N_{\phi}$. If a layer has LL index $n>0$, we assume that the LLs $0,1,\ldots,n-1$ are fully occupied and inert, while the $n$th LL is partially occupied. We study the model in Eq. (17) using exact diagonalization of the Hamiltonian matrix appropriate for systems of finite size, either on the surface of a sphere Haldane (1983) or torus Yoshioka et al. (1983). Various symmetries, either spatial (e.g., translation, inversion) or internal (e.g., pseudospin U(1) symmetry in the absence of interlayer tunneling), are used to block-diagonalize the Hamiltonian and classify its eigenstates. We also interpret our results against Hartree-Fock theory Côté et al. (1992), which at $\nu=1$ captures the leading effects as quantum fluctuations are fairly weak given the small $d/\ell_{B}$.

Experimental results in the main text demonstrate the quantum Hall effect at integer filling factors when the two layers have matched LL indices, $n_{1}=n_{2}$. Here we provide theoretical evidence that identifies this phase as exciton condensate, where each exciton consists of an electron and a hole in the opposite layer (see the recent review Eisenstein (2014)). While the exciton phase exhibits quantized Hall effect, its distinctive signature is the Goldstone mode which governs the gapless excitations of the system in the neutral channel, probed by driving oppositely directed electrical currents through the two layers Spielman et al. (2000); Tutuc et al. (2004).

In the numerics we can access the Goldstone mode directly by computing the energy dispersion, $E(k)$, as a function of momentum $k$ in the torus geometry. The torus spectra for a fixed interlayer distance $d/\ell_{B}=0.1$ and both layers having the same LL index $n=0,1,2,3$ are shown in Fig. 14(a)-(d). We compare exact results against the prediction of Hartree-Fock theory Fertig (1989), shown by the blue line in Fig. 14. We observe, overall, an excellent agreement between the exact results and the Hartree-Fock prediction, in particular in view of the fact that Hartree-Fock results do not contain any adjustable parameters. The agreement could be anticipated given that we are looking at small $d/\ell_{B}$ and the Hartree-Fock theory captures the ground state exactly for $d/\ell_{B}=0$. In comparison, much bigger discrepancy between exact results and Hartree-Fock theory are seen already at $d/\ell_{B}=0.5$, even in the long-wavelength limit (data not shown).

Results in Fig. 14 confirm the existence of a gapless mode with a dispersion in quantitative agreement with the Hartree-Fock theory based on easy-plane XY ferromagnetism Girvin and Macdonald (1995). Unfortunately, accessing the long-wavelength of this mode, $k\to 0$, is challenging for exact numerics where the smallest available (non-zero) $k$ scales as $\propto 1/\sqrt{N_{\phi}}$, thus we cannot resolve finer features of the mode, such as whether the dispersion is quadratic (as expected at $d=0$ when the system is SU(2)-invariant) or linear (for $d>0$). Interestingly, the exact spectra in Fig. 14 suggest that the Goldstone mode is increasingly better separated from the continuum of the excitation spectrum as LL index is increased, at all wavelengths – see Fig. 14(d). While this could be due to stronger finite-size effects at larger $n$, the fully exposed collective mode in bilayer $\mathrm{WSe_{2}}$ may lead to a different optical response compared to GaAs and graphene materials, where typical $d/\ell_{B}$ are several times larger and the collective mode is only exposed below the continuum in the range of momenta $k\ell_{B}\lesssim 1$.

Further evidence for the exciton condensate follows from the study of the system’s ground state properties. At $d=0$, the exact ground state of the system is given by the first-quantized wave function known as the 111 Halperin state Halperin (1983):

where $z_{i}$, $w_{i}$ are complex 2D coordinates of electrons in top and bottom layers, respectively. The Jastrow factor, $\prod_{i,j}(z_{i}-w_{j})$, builds in the exciton correlation between an electron in one layer ($z_{i}$) and the corresponding hole at the same position $w_{i}$ in the opposite layer. Although $\psi_{111}$ is only an exact ground state when strictly $d=0$, it remains a good representative wave function for the entire exciton phase extending to $d\sim\ell_{B}$ Girvin and Macdonald (1995). Thus, by evaluating the inner product between the exact ground state of a bilayer and $\psi_{111}$ we can test the existence of the exciton phase in bilayers with $n>0$. These overlaps are shown in Fig. 15(a) as a function of $d/\ell_{B}$ in the sphere geometry. We observe that for any value of $n$, as long as $d/\ell_{B}\lesssim 0.3$, the overlap with 111 state is very close to unity. At larger values of $d$, the overlap drops rapidly, and the decay is faster for larger values of $n$, indicating the sensitivity of the exciton phase in higher LLs to the increase in layer separation.

Further, we do not need to rely on a microscopic candidate wave function to detect the exciton phase. In Fig. 15(b) we have computed the exciton order parameter Shibata and Yoshioka (2006),

expressed in terms of an operator that creates an exciton at a given point $\mathbf{r}$,

where $\hat{\Psi}_{\sigma}$ is the standard electron field operator. The normalization of $\hat{\Psi}_{\mathrm{exc}}$ is chosen such that the expectation value is 1 when $r\gg\ell_{B}$ in the ideal exciton case (i.e., for the exact 111 Halperin state). In contrast, a value close to 0 indicates no exciton correlations. Specifically, on a sphere geometry, we can pick the origin to be the north pole and the point far away is the south pole. The value of the correlator in Eq. (22) is shown in Fig. 15(b); when $d/\ell_{B}$ is small ($\lesssim 0.4$), the exciton order parameter is small in any of the LLs ($n=0,1,2,3$). For $d/\ell_{B}=0.6$ the correlator decays sharply with $n$, and the exciton phase is only likely to exist in $n=0$ LL.

By generalising the Hartree-Fock theory of Ref. Fertig (1989), we can conveniently estimate the critical interlayer distance $d_{c}/\ell_{B}$ when the exciton condensate is destroyed by a density wave that breaks translation symmetry. In any LL, we find that the exciton dispersion develops a “roton”-like minimum as $d/\ell_{B}$ becomes sufficiently large – see an example of $n=2$ in Fig. 16(a). The minimum becomes more pronounced by increasing $d/\ell_{B}$ and eventually goes soft, signalling a quantum phase transition where the ground state changes its momentum from 0 to a finite value $k^{*}\sim 1/\sqrt{2n+1}$ Brey and Fertig (2000). We numerically perform the golden section search and determine the critical value of $d_{c}/\ell_{B}$, which is shown in Fig. 16(b) as a function of $n$. The critical distance decays sharply for low values of $n$, reaching $d_{c}/\ell_{B}\sim 0.3$ for $n\lesssim 10$. Note that this estimate does not include the effect of quantum fluctuations, which become significant near the transition and potentially further decrease $d_{c}/\ell_{B}$. Nevertheless, this estimate is consistent with the possibility of the exciton phase in bilayer $\mathrm{WSe_{2}}$ where we estimate $d/\ell_{B}\approx 0.1$. We note that in recent works Shi et al. (2008); Zhu et al. (2019) the breakdown of the exciton phase was studied at GaAs total filling $\nu=5$, which corresponds to $n=1$ in our case. Refs. Shi et al. (2008); Zhu et al. (2019) estimated $d_{c}/\ell_{B}\approx 0.9-1.2$, which is somewhat larger than our estimate in Fig. 16(b), possibly due to finite size effects.

We have also studied the response of the exciton phase to adding or removing charge, which is characterized by the chemical potential discontinuity Yoshioka et al. (1983); Moon et al. (1995):

where $E_{N,N_{\phi}}$ is the ground state energy for the given values of $N$ and $N_{\phi}$. This quantity plays the role of the activation gap of the system: the non-zero value of $\Delta\mu$, in the limit $N\to\infty$, signals an incompressible ground state that would give rise to the quantized Hall effect. Alternatively, one could keep $N$ fixed and vary flux $N_{\phi}\to N_{\phi}\pm 1$. At $\nu=1$, it is more convenient to keep the geometry of the system unchanged, i.e., fix $N_{\phi}$ and vary $N$.

As we are adding or removing charge, it is important to include the relevant electrostatic corrections when evaluating $\Delta\mu$ in finite-size systems. For example, the Coulomb interaction gives rise to electrostatic charging energy when particles are shuffled between the layers Moon et al. (1995). Formally, this term originates from the finite $\mathbf{q}=0$ contribution in the pseudospin Coulomb interaction, $(V_{\uparrow\uparrow}(0)-V_{\uparrow\downarrow}(0))/2=\pi d$, and it has a form $\epsilon_{c}\propto dS_{z}^{2}$, where the total pseudospin $S_{z}\equiv(N_{\uparrow}-N_{\downarrow})/2$. On a finite sphere, the charging energy can be shown to be given by

where the radius of the sphere is $R=\sqrt{N_{\phi}/2}$ (in units $\ell_{B}=1$).

In Fig. 17 we have computed $\Delta\mu$ as a function of density imbalance, $-D(N_{\uparrow}-N_{\downarrow})$. Panel (a) shows the charge gap for $n_{1}=n_{2}=0$, while panel (b) illustrates the case of a higher LL with $n_{1}=n_{2}=2$. We find the behaviour of the gap of any $n>0$ to be qualitatively similar to each other and distinct from $n=0$ LL. Although we have previously established the ground state for any $n_{1}=n_{2}=n$ is the same exciton condensate, approximately described by the 111 state in Eq. (21), we now see that the charged excitations are very different in $n=0$ compared to $n>0$.

In $n=0$ case, the gap rises sharply upon density imbalance, which is reminiscent of a skyrmion excitation Sondhi et al. (1993); Schmeller et al. (1995) that involves a large number of pseudospin flips. In the absence of Zeeman field ($D=0$), the number of flipped pseudospins is extensive in system size; for a finite value of $D$, the number of flips is optimized to minimize the total energy due to pseudospin stiffness and the competing Zeeman term. When $D$ is sufficiently large, all pseudospins orient in the $z$-direction [shaded regions in Fig. 17(a)], and the gap then trivially grows linearly with $D$. Adding a small amount of interlayer tunneling, $-\Delta_{\mathrm{SAS}}S_{x}$, does not qualitatively change the behaviour of $\Delta\mu$. We note that the increase of the gap in $n=0$ GaAs bilayers has been observed experimentally in Refs. Spielman et al. (2004) and Tutuc et al. (2003). In these cases, however, the gap has a much slower parabolic rise with $D$ due to the larger value of $d/\ell_{B}$.

In contrast, in higher $n>0$ LLs, the gap is largely insensitive to imbalance until full pseudospin polarization is reached, see Fig. 17(b). This is consistent with the behavior found in Hartree-Fock approximation which does not include a possibility of pseudospin-texture excitations Jungwirth et al. (1998); Jungwirth and MacDonald (2000). In the Hartree-Fock picture, the ground state pseudospin smoothly rotates as bias potential is varied. The full pseudospin polarization is reached when the bias potential reaches the value (in the notation of Ref. Jungwirth et al. (1998))

where $V_{\sigma\sigma}$ denotes the interaction in the pseudospin channel. It follows that $U_{\sigma\sigma}^{(n)}\propto\sqrt{2n+1}$, thus it is harder to polarise the pseudospin in higher LLs, which is consistent with the numerics as well as the experimental data presented in the main text. In Fig. 17(b) tiny wiggles can be observed near $D=0$ – those correspond to the ground state pseudospin changing in integer steps as a single electron is transferred between the two layers.

The difference in types of charged excitations between $n=0$ and $n>0$ cases is reminiscent of integer quantum Hall ferromagnets, where it has been established that skyrmions are energetically competitive with ordinary electron-hole excitations only in $n=0$ LL Wu and Jain (1994); Wu and Sondhi (1995). In Fig. 18(a) we analytically compare the energy gap for creating a spin texture with that for creating an ordinary particle-hole excitation (i.e., with a single pseudospin flipped). Indeed, in $n=0$ we find the skyrmion excitation to be significantly lower in energy Sondhi et al. (1993), while in $n=1$ the ordinary particle-hole excitation is already more stable and continues to be energetically favored in all $n\geq 1$.

Experimental data in the main text shows the gap in bilayer $\mathrm{WSe_{2}}$ varies non-monotonically with $n$ and reaches a maximum in $n=2$. This suggests that excitations in both $n=0$ (at sufficiently high magnetic fields) and in $n=1$ LL are pseudospin textures. This discrepancy with respect to theory could be explained by disorder, which affects the two types of excitations very differently due to their different spatial profiles, but it could be also be an interaction effect, e.g., due to screening of the Coulomb interaction. In Fig. 18(b) we have explored the difference in the gap for creating the particle-hole excitation vs. a pseudospin texture in $n=1$ LL, by varying the interlayer distance $d$ and the amount of screening of the interaction. The latter is phenomenologically modelled by finite thickness ansatz, i.e., we model the intralayer Coulomb repulsion as $1/\sqrt{r^{2}+w^{2}}$ , while the interlayer potential is $1/\sqrt{r^{2}+(w+d)^{2}}$. We observe that either the increase in $d$ or some screening $w\gtrsim 0.1$ is sufficient to stabilize the skyrmion excitation, in line with previous results for the integer quantum Hall ferromagnets Cooper (1997); Wójs and Quinn (2002). Thus, it is likely that these interactions effects, perhaps in combination with disorder, stabilize the pseudospin texture in $n=1$ LL of $\mathrm{WSe_{2}}$ over the particle-hole excitation.

Finally, we comment on the case of unmatched LLs, $n_{1}\neq n_{2}$. In this case, experimental data in the main text suggest an absence of the exciton phase. In Hartree-Fock theory Jungwirth et al. (1998); Jungwirth and MacDonald (2000), the system in this case behaves an easy-axis ferromagnet. When bias potential is introduced, the gap varies linearly with $D$ until a critical value when the pseudospin flips, $S_{z}=N/2\to-N/2$. The excitation spectrum on the torus is found to be very different from the Goldstone collective mode seen in Fig. 14, possibly due to the existence of domain wall type excitations Rezayi et al. (2003). Similarly, the overlap with the 111 state as well as the exciton order parameter are found to be small and rapidly diminishing with system size.