Searching for Unconventional Superfluid in Excitons of Monolayer Semiconductors

In two-dimensional (2D) homogeneous systems, it is well known that continuous symmetry cannot be broken spontaneously according to Mermin-Wagner theorem [1, 2], and there is no true long-range order. As a special example, Bose-Einstein condensation (BEC) critical temperature in 2D is zero. Nevertheless, superfluid is still possible at finite temperature, and the transition from superfluid to normal phase is described by the Berezinskii-Kosterlitz-Thouless (BKT) theory [3, 4], where the underlying mechanism is the vortex-antivortex pair unbinding at high temperature. This generic paradigm is successful in the description of a variety of 2D superfluids, including liquid helium films [5], superconductors [6], cold atomic gases [7, 8], exciton-polariton condensates [9], and dipolar excitons [10].

Generically, for a 2D bosonic system, the nature of low temperature phases depends crucially on the density of states. Deviation from homogeneous space or parabolic dispersion may lead to superfluids beyond the conventional paradigm. For example, when ultracold bosonic atoms are confined in a harmonic trap, the BEC critical temperature is nonzero [11]. Similar examples include bosons confined on the surface of a sphere, where a finite BEC critical temperature also exists due to finite size effect [12]. On the other hand, for homogeneous bosons with quartic dispersion realized at the transition point of spin-orbit coupled gases, BKT transition temperature vanishes and the low temperature phase is characterized by an algebraic order [13]. This is an interesting example of enhanced low energy fluctuations brought by increased density of states. The contrary case, i.e., interacting bosons with single-particle dispersion $\epsilon(k)\sim k^{\nu}$ ($\nu<2$) in realistic experimental systems are still lacking. This type of system is of fundamental importance and interesting, as one of the three exhaustive cases of an isotropic 2D bosonic system in homogeneous space, e.g., $\nu<2$, $\nu=2$, and $\nu>2$. In a word, 2D homogeneous bosons beyond the conventional paradigm are quite rare and highly interesting, and will surely enrich our understanding of such fundamental concepts as BEC and superfluidity.

In recent years, there has been growing interest in the realization of exciton condensation in bilayer 2D materials [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25], such as transition metal dichalcogenide (TMD). Besides the long-sought interlayer exciton condensation in heterobilayers, the possibility of intralayer exciton condensation in monolayer TMD can not be excluded either [26], despite of its relatively short lifetime. In particular, significant experimental progress has been achieved in the exciton fluid at thermal equilibrium in monolayer TMD [27], which paves the way towards the realization of exciton superfluid in this system.

Here we point out that in the system of monolayer TMD, intralayer exciton with linear dispersion [28, 29] provides another example beyond conventional paradigm of superfluids. Note that system anisotropy can lead to interesting anisotropic superfluids in 2D materials [30, 31, 32], also beyond the conventional paradigm of isotropic superfluids, but here we only consider the isotropic superfluids. The center-of-mass motion of intralayer exciton in monolayer TMD features intrinsic valley-orbit coupling [28, 29], originating from electron-hole exchange interaction. With additional valley Zeeman energy introduced by either magnetic field or valley-selective optical Stark effect, the dispersion of intralayer exciton at lower branch becomes linear (see Fig. 1(b)). We will reveal that this special dispersion endows intralayer excitons with a variety of novel features in low energy physics. Our main findings include the following: (i) With a valley Zeeman energy, the BEC critical temperature of this 2D homogeneous system is nonzero and exhibits rapid increase with the field strength. (ii) The Bogoliubov excitation spectrum of this exciton condensate shows $\varepsilon(\boldsymbol{q})\sim\sqrt{q}$ behaviour at long wavelength in the presence of valley Zeeman energy, which is an usual form of gapless Goldstone mode, in contrast to conventional phonon excitation. (iii) The vortex energy deviates from logarithmic form with respect to system size, resulting in a non-BKT-type phase transition for large system size. There exists a characteristic system size, beyond which the superfluid to normal phase transition evolves from BKT type to 3D-like without discontinuous jump in superfluid density. This crossover can also be observed in a single system by tuning the magnitude of valley Zeeman energy. (iv) With the increase of temperature, the system undergoes a two-step phase transition, first from a BEC with long-range order to a superfluid with quasi-long-range order, and then to a normal phase. These novel phases can be experimentally detected by measuring the spatial and temporal coherence of emitted photons by excitons.

Monolayer exciton dispersion. There are two inequivalent valleys in the Brillouin zone corner of monolayer TMD, denoted as $\pm K$ valley, respectively. Valley excitons in monolayer TMD behave as a pseudospin-$1/2$ bosonic system, whose valley pseudospin is coupled with center-of-mass momentum of excitons [28, 29], giving rise to the valley-orbit coupling. The effective Hamiltonian describing the center-of-mass motion of intralayer excitons with valley-orbit coupling reads [28, 33, 29, 34, 35, 36]

	$\displaystyle\hat{H}_{0}=$	$\displaystyle\frac{\hbar^{2}Q^{2}}{2m}+AQ+AQ\cos(2\theta_{\boldsymbol{Q}})\sigma_{x}+AQ\sin(2\theta_{\boldsymbol{Q}})\sigma_{y}$
		$\displaystyle+\delta\sigma_{z}.$		(1)

Here $\boldsymbol{Q}$ is the center-of-mass momentum of exciton with magnitude $Q=|\boldsymbol{Q}|$ and angle $\theta_{\boldsymbol{Q}}$, $m\approx 1.1m_{0}$ is exciton effective mass, $\sigma_{i}$ $(i=x,y,z)$ are the Pauli matrices of valley pseudospin, and $A\approx 0.9$ eV$\cdot\mathrm{\AA}$ is the valley-orbit coupling strength [29], related with electron-hole exchange interaction. With the exchange interaction modeled by the screened Keldysh potential [33, 35], which gives a good description of the exciton Rydberg series in monolayer TMD [37, 48], the $AQ$ term is an approximation of the more accurate $AQ/(1+r_{0}Q)$ term, where both $A$ and the screening length $r_{0}$ depend on the effective dielectric constant $\epsilon_{d}$, determined by surrounding environment. Unless otherwise specified, the value of $A$ adopted here corresponds to a free-standing TMD monolayer with $r_{0}$ neglected, realizable with TMD placed on top of circular holes prepatterned on SiO₂ substrate, using the experimental setup in Ref. [38]. For TMD placed on substrate materials, the screening will reduce the value of $A$ and hence modifies the quantitative results such as the critical temperature calculated below (see Fig. 2 and Ref. [39]). We also introduce a $\delta\sigma_{z}$ term, known as the valley Zeeman energy, with magnitude $\delta$ tunable by applying a magnetic field [42, 43] or utilizing the valley-selective optical Stark effect [44, 45]. The non-interacting dispersion has two branches, given by $\xi_{\pm}(\boldsymbol{Q})=\hbar^{2}Q^{2}/2m+AQ\pm\sqrt{A^{2}Q^{2}+\delta^{2}}$, as shown in Figs. 1(a) and (b). For finite $\delta$, the lower branch of dispersion at small momentum features a linear spectrum, distinct from usual parabolic dispersion and brings dramatic change to low energy properties of this system. The range of linear dispersion defines a characteristic momentum $Q_{c}$, at which the parabolic and linear terms are comparable in weight, and corresponds to a characteristic system size $L_{c}=1/Q_{c}$, both shown in Fig. 3(b).

Non-interacting BEC critical temperature. In the absence of valley Zeeman energy, the lower branch of dispersion is parabolic, and obviously the BEC critical temperature vanishes, as in conventional 2D bosonic systems. With valley Zeeman energy ($\delta\neq 0$), the lower branch of dispersion is linear at small momentum, which modifies the low energy density of states and renders a BEC at finite temperature possible. We first neglect the exciton-exciton interaction and calculate the non-interacting condensation temperature $T_{\mathrm{BEC}}$ through the relation,

$$n=\frac{1}{(2\pi)^{2}}\int\sum_{\tau=\pm}\frac{d^{2}\boldsymbol{Q}}{e^{\beta[\xi_{\tau}(\boldsymbol{Q})-\mu(T)]}-1},$$

(2)

where $\mu$ is the exciton chemical potential and $\beta=1/k_{B}T$. When $T\rightarrow T_{\mathrm{BEC}}$, $\mu\rightarrow-\delta$. The calculated relation between $T_{\mathrm{BEC}}$ and $\delta$ is shown in Fig. 2 for different exciton densities in the two cases illustrated by the insets. Clearly, $T_{\mathrm{BEC}}$ monotonically increases with $\delta$. For typical exciton density $n=10^{10}\sim 10^{12}$ cm^-2 below the Mott limit ($\sim 10^{13}$ cm^-2) [27], a moderate $\delta$ can lead to a relatively high $T_{\mathrm{BEC}}$ on the order of $10\sim 100$ K in the free-standing case. This is in stark contrast to the conventional bosonic system with parabolic dispersion, where $T_{\mathrm{BEC}}=0$ and only a quasi-condensate exists at finite temperature. Here one has a true 2D condensate with long-range order, in homogeneous space. As will be shown below, we find that the inclusion of exciton-exciton interaction can further enhance the BEC critical temperature.

Bogoliubov excitation spectrum. The linear dispersion of non-interacting exciton at small momentum also implies unusual low energy excitations of exciton condensate. Assuming that excitons condense at zero momentum state, the low energy excitation can be calculated by the standard Bogoliubov theory (see for example, [51]). The mean-field energy functional reads

	$\displaystyle E[\Psi_{\sigma}]$	$\displaystyle=\int d^{2}\boldsymbol{r}\Bigg{\{}\left(\Psi_{\uparrow}^{*},\Psi_{\downarrow}^{*}\right)\Big{[}\left(\frac{\hbar^{2}{Q}^{2}}{2m}+AQ\right)+\delta\sigma_{z}$
		$\displaystyle+AQ\cos(2\theta_{\boldsymbol{Q}})\sigma_{x}+AQ\sin(2\theta_{\boldsymbol{Q}})\sigma_{y}\Big{]}\left(\begin{array}[]{r}\Psi_{\uparrow}\\ \Psi_{\downarrow}\end{array}\right)$		(5)
		$\displaystyle+\frac{c_{1}}{2}\left(\left|\Psi_{\uparrow}\right|^{4}+\left|\Psi_{\downarrow}\right|^{4}\right)+{{c}_{2}}{\left|{{\Psi_{\uparrow}}}\right|^{2}}{\left|{{\Psi_{\downarrow}}}\right|^{2}}\Bigg{\}},$		(6)

where $\Psi_{\sigma}$ is the condensate wave function for pseudospin $\sigma=\uparrow,\downarrow$, corresponding to $\pm K$ valley. At the qualitative level, the exciton-exciton interaction is modelled as a contact interaction, which is valid when the exciton density is low [52, 17], with $c_{1}$ and $c_{2}$ being exciton-exciton interaction strengths between the same and different pseudospin states, respectively.

Within the same framework, the condensate dynamics can be described by the spinor Gross-Pitaevskii (GP) equation,

$$i\hbar\frac{\partial}{{\partial t}}\left({\begin{array}[]{*{20}{c}}{{\Psi_{\uparrow}}}\\ {{\Psi_{\downarrow}}}\\ \end{array}}\right)=\left({\begin{array}[]{*{20}{c}}{H_{\uparrow\uparrow}}&{{H_{A}}}\\ {H_{A}^{*}}&{H_{\downarrow\downarrow}}\\ \end{array}}\right)\left({\begin{array}[]{*{20}{c}}{{\Psi_{\uparrow}}}\\ {{\Psi_{\downarrow}}}\\ \end{array}}\right),$$

(7)

where $H_{\uparrow\uparrow}=\hbar^{2}Q^{2}/2m+A{Q}+\delta+{c_{1}}{{\left|{{\Psi_{\uparrow}}}\right|}^{2}}+{c_{2}}{{\left|{{\Psi_{\downarrow}}}\right|}^{2}}$, $H_{\downarrow\downarrow}=\hbar^{2}Q^{2}/2m+A{Q}-\delta+{c_{1}}{{\left|{{\Psi_{\downarrow}}}\right|}^{2}}+{c_{2}}{{\left|{{\Psi_{\uparrow}}}\right|}^{2}}$, and $H_{A}=AQ\exp(-2i\theta_{\boldsymbol{Q}})$. In the presence of valley Zeeman energy ($\delta>0$), minimization of the mean-field energy functional (Eq. 6) gives the wave function of an exciton condensate at ground state. The ground state is found to be $(|\Psi_{\uparrow}|^{2},|\Psi_{\downarrow}|^{2})=n_{0}(1/2-\delta/(n_{0}c_{1}-n_{0}c_{2}),1/2+\delta/(n_{0}c_{1}-n_{0}c_{2}))$ for $0\leq\delta\leq(n_{0}c_{1}-n_{0}c_{2})/2$, while $(|\Psi_{\uparrow}|^{2},|\Psi_{\downarrow}|^{2})=n_{0}(0,1)$ for $\delta>(n_{0}c_{1}-n_{0}c_{2})/2$, with $n_{0}$ being the condensate density. For intralayer exciton, $c_{1}>c_{2}$ ensures that when $\delta=0$, exciton densities in two valleys are equal at ground state. With the estimated values $c_{1}\sim 1.0$ eV$\cdot$nm², $c_{2}\sim 0.6$ eV$\cdot$nm² [46, 47], and typical exciton density $n_{0}=1.0*10^{11}$ cm^-2, the critical $\delta_{c}\equiv(n_{0}c_{1}-n_{0}c_{2})/2\sim 0.2$ meV, beyond which the excitons are completely valley polarized. In the following, we consider the case $\delta>\delta_{c}$, which is readily accessible in experiment with moderate magnetic field and greatly simplifies the calculations. Results within the small range $0<\delta<\delta_{c}$ are also calculated with the wave function $n_{0}(0,1)$ as an approximation.

By expanding the wave function around the stationary state, $\Psi^{\prime}_{\sigma}=\Psi_{\sigma}+\delta\Psi_{\sigma}$, and assuming the form of perturbation $\delta\Psi_{\sigma}=\Psi_{\sigma}e^{-i\mu/\hbar}[u_{\sigma}(\boldsymbol{q})e^{-i\varepsilon(\boldsymbol{q})t}-v_{\sigma}^{\ast}(\boldsymbol{q})e^{i\varepsilon(\boldsymbol{q})t}]$, one arrives at the Bogoliubov equation for the low energy excitation,

$\displaystyle\mathcal{M}\left(\begin{array}[]{c}u_{\uparrow}\\ u_{\downarrow}\\ v_{\uparrow}\\ v_{\downarrow}\end{array}\right)=\varepsilon\left(\begin{array}[]{c}u_{\uparrow}\\ u_{\downarrow}\\ v_{\uparrow}\\ v_{\downarrow}\end{array}\right),$

(16)

with

$\displaystyle\mathcal{M}=\begin{pmatrix}H_{1}^{+}&B^{+}_{1}&0&0\\ B^{+}_{2}&H_{2}^{+}&0&-n_{0}c_{1}\\ 0&0&H_{1}^{-}&B^{-}_{1}\\ 0&n_{0}c_{1}&B^{-}_{2}&H_{2}^{-}\\ \end{pmatrix}.$

(17)

Here $H_{1}^{\pm}=\pm\left(\hbar^{2}q^{2}/{2m}+Aq+{n_{0}c_{2}}-\mu+\delta\right)$, $H_{2}^{\pm}=\pm\left(\hbar^{2}q^{2}/{2m}+Aq+{2}{n_{0}c_{1}}-\mu-\delta\right)$, $B^{\pm}_{1}=\pm A\left(q^{2}_{x}-q_{y}^{2}-2iq_{x}q_{y}\right)/q$, $B^{\pm}_{2}=\pm A\left(q^{2}_{x}-q_{y}^{2}+2iq_{x}q_{y}\right)/q$, and $\mu=-\delta+n_{0}c_{1}$. There are four groups of eigenvalues and only the two whose corresponding eigenvectors satisfy $|u_{\sigma}|^{2}-|v_{\sigma}|^{2}=1$ are physical. The two branches of excitations are

$$\varepsilon_{\pm}(\boldsymbol{q})=\frac{1}{2}\sqrt{\mathcal{A}+\mathcal{B}\pm 2\sqrt{\mathcal{C}^{2}-\mathcal{D}}},$$

(18)

where the expressions of $\mathcal{A}$, $\mathcal{B}$, $\mathcal{C}$ and $\mathcal{D}$ are lengthy and listed in the Supplementary Material [39]. As shown in Figs. 1(c) and 1(d), for $\delta=0$, $\varepsilon_{-}(\boldsymbol{q})\sim q$ at small $q$, while $\varepsilon_{-}(\boldsymbol{q})\sim\sqrt{q}$ for $\delta>0$ [39], which is a new form of gapless Goldstone mode unreported before.

Interacting BEC critical temperature. With exciton-exciton interaction, the dispersion of low energy excitation is modified, and the resulting change of density of states also affects the BEC critical temperature. We quantitatively calculate the interacting critical temperature using the standard Hartree-Fock-Bogoliubov-Popov theory [53, 54, 55], by solving self-consistently the total exciton density [39]

$$n=n_{0}+\sum_{{\tau=\pm}\atop{\sigma=\uparrow,\downarrow}}\int\frac{d^{2}\boldsymbol{q}}{(2\pi)^{2}}\left\{|v_{\sigma}^{\tau}(\boldsymbol{q})|^{2}+\frac{|u_{\sigma}^{\tau}(\boldsymbol{q})|^{2}+|v_{\sigma}^{\tau}(\boldsymbol{q})|^{2}}{e^{\varepsilon_{\tau}(\boldsymbol{q})/k_{B}T}-1}\right\},$$

(19)

where the condensate density $n_{0}$ also enters the excitation spectrum and quasiparticle amplitudes $u_{\sigma}(\boldsymbol{q})$ and $v_{\sigma}(\boldsymbol{q})$. By calculating the dependence of $n_{0}$ on $T$, one can extrapolate to $n_{0}\rightarrow 0$ and find the critical temperature $T_{\mathrm{BEC}}$ with interaction. As shown in Fig. 2, the BEC critical temperature is enhanced by exciton-exciton interaction, which can be qualitatively understood by considering that the density of states is reduced at low energy and thereby the condensate fraction at given temperature is increased compared with non-interacting case.

Vortex energy. In 2D bosonic systems, besides the non-singular excitations calculated above, there are also singular topological excitations, i.e., vortices, which play a decisive role in conventional BKT theory. In a single-component condensate, the appearance of a free vortex causes an energy increase $E_{\mathrm{v}}\sim mn_{s}/2\int d\boldsymbol{r}v^{2}(\boldsymbol{r})\simeq\pi\hbar^{2}n_{s}/m\ln(L/\zeta)$, where $v(\boldsymbol{r})$ is the magnitude of velocity, $L$ is the system size, $\zeta=\sqrt{\hbar^{2}/2mn_{s}c_{1}}$ is the condensate healing length, and $n_{s}$ is the superfluid density. This energy function has the same form with increase of entropy $S\simeq 2k_{B}\ln(L/\zeta)$ associated with a free vortex, thereby leading to the free energy change $\Delta\mathcal{F}=E_{\mathrm{v}}-TS\simeq(\pi\hbar^{2}n_{s}/m-2k_{B}T)\ln(L/\zeta)$, whose turning point gives $T_{\mathrm{SF}}$. This is the case for parabolic particle dispersion, and since intralayer exciton exhibits different dispersion here, the vortex energy is modified, as well as the nature of superfluid phase transition. In this pseudospin-1/2 system, vortices generally have two components with respective vorticity or circulation. Here we are interested in the case of $\delta>\delta_{c}$, where the condensate at ground state is pseudospin polarized, and thus we make the approximation by considering a single-component vortex with modified dispersion, i.e., $\xi_{-}(\boldsymbol{Q})=\hbar^{2}Q^{2}/2m+AQ-\sqrt{A^{2}Q^{2}+\delta^{2}}$. With similar argument, the vortex energy should be modified as $E_{\mathrm{v}}\sim mn_{s}/2\int d\boldsymbol{r}\left\{v^{2}(\boldsymbol{r})+(mA/\hbar)v(\boldsymbol{r})\right\}\sim(\pi\hbar^{2}n_{s}/m)\left[\alpha\ln(L/\zeta)+\beta(L/\zeta)\right]$, where $\alpha$ and $\beta$ are two constants dependent on $A$ and $\delta$. In particular, $\beta$ should vanish at $\delta=0$ and increase with $\delta$. Similar form of vortex energy is found in superconductors or two-component BECs with Josephson coupling [56, 57, 58, 59, 60], and the additional linear term $\beta(L/\zeta)$ implies the breakdown of BKT theory for large system size.

To be on a firmer ground, we adopt the trial wave function of a vortex as in the single-component case, i.e., $\psi(\boldsymbol{r})\approx\sqrt{n_{s}}re^{i\phi_{\boldsymbol{r}}}/{\sqrt{r^{2}+2}}$, where $r$ is in unit of $\zeta$, being a good fit to the numerical solution by Gross-Pitaevskii equation [61]. We numerically calculate the kinetic energy increase associated with the vortex within the new dispersion, through $E_{\mathrm{v}}\simeq\int d\boldsymbol{r}\psi^{*}(\boldsymbol{r})\hat{H}_{\mathrm{kin}}\psi(\boldsymbol{r})$, with $\hat{H}_{\mathrm{kin}}$ being the kinetic energy operator corresponding to dispersion $\xi_{-}(\boldsymbol{Q})$. The operation is carried out in momentum space, bypassing the difficulty in expressing the valley-orbit coupling term in coordinate space [39]. The obtained dependence of $E_{\mathrm{v}}$ on $L/\zeta$ can be fitted well with the relation above, which gives the specific value of $\alpha$ and $\beta$. Within the experimentally feasible range of $\delta\lesssim 7.5$ meV, we find $\alpha\approx 1$ and the change of $\beta$ with $\delta$ is shown in Fig. 3(a). The condition of equal contribution from these two terms $\alpha\ln(L/\zeta)+\beta(L/\zeta)$ also determines a characteristic system size, in quantitative agreement with $L_{c}$ defined above.

Superfluid to normal phase transition. The linear term of vortex energy will dominate for large system size, and since this term grows faster than entropy, proliferation of free vortices will always be suppressed even at high temperature. This means the main contribution to the depletion of superfluid density will come from non-singular excitations, including the Bogoliubov excitation with $\sqrt{q}$ dispersion, and vortex-antivortex bound pairs. These types of excitations will lead to continuous decrease of superfluid density with temperature, without discontinuous jump in superfluid density at the phase transition point, similar to 3D case. So we conclude that, for system size beyond a characteristic scale $L_{c}$, the superfluid phase transition is no longer BKT type. On the other hand, for small system size with large lower bound of momentum, the effect of linear dispersion will be minor, and in this case we still expect a BKT-type finite-size crossover. Note also that the characteristic system size $L_{c}$ depends on $\delta$, and thus the two limiting case above can also be observed in a single system by tuning $\delta$. For small $\delta$, $L_{c}$ is on the order of $\mu$m, fully within the range of sample size in current experiments.

We finally have enough information on the finite temperature phase diagram of this system in $\delta-T$ plane, as shown in Fig. 4. At $\delta=0$, the BEC critical temperature is zero, while the BKT transition temperature, given by the well-known Nelson-Kosterlitz relation [62], is nonzero. With the increase of $\delta$, the BEC critical temperature increases and the superfluid critical temperature should also increases, since both free vortices and non-singular excitations are asymptotically suppressed. For large $\delta$, where the free vortices are completely suppressed for a given system size, the superfluid density should continuously drop to zero, similar to 3D case, where the superfluid critical temperature coincides with BEC critical temperature. So both $T_{\mathrm{BEC}}$ and $T_{\mathrm{SF}}$ increase with $\delta$, and asymptotically approach each other. One immediately realizes that the phase diagram consists of three different phases, i.e., BEC phase with long-range order (also a superfluid), superfluid phase with quasi-long-range order, and a trivial normal phase. In other words, at finite $\delta$, with the increase of temperature, the system undergoes a two-step phase transition, first from BEC to a non-BEC superfluid, and then to a normal phase. Note that the superfluid critical temperature is only qualitatively demonstrated in Fig. 4, by interpolation between two known limiting cases ($\delta=0$ and large $\delta$). A quantitative treatment of superfluid phase transition, taking full account of vortex-antivortex pair excitation induced screening, using methods such as Monte Carlo simulation [63, 64] will be postponed to a future work.

Experimental observation. We now comment on the possible experimental realization of exciton superfluid in monolayer TMD. The major obstacle is the short exciton lifetime, whose magnitude relative to exciton thermalization time determines whether thermal equilibrium can be reached. Recently, it is reported experimentally that an exciton fluid at thermal equilibrium in monolayer MoS₂ was observed [27]. This important finding points to the exciting possibilities for studying the exciton superfluid in monolayer system and also demonstrates the experimental relevance of our theoretical model. Using an experimental setup similar to Ref. [38], our predictions can be readily verified based on existing techniques. The three different phases predicted in this work can be distinguished by measuring the spatial and temporal coherence of emitted photons, from which we can infer the decay behaviour of exciton correlation function. Similar techniques have been used in previous experimental demonstration of BKT transition in exciton-polariton condensate [9] and dipolar excitons [10]. Note that in the presence of electron/hole doping, the screened Coulomb interaction is better treated by the Thomas-Fermi screening, which will turn the single exciton dispersion from linear to the conventional parabolic type [39], and hence should be avoided in experiments to observe the phenomena predicted here.

In summary, we have studied the unusual properties of exciton superfluid in monolayer TMD, possibly realizable based on recent experimental progress [27], which is an unique bosonic system beyond conventional paradigm. With the valley-orbit coupling and valley Zeeman energy, one can manipulate the dispersion of exciton center-of-mass motion. The attainable linear dispersion in lower branch brings a variety of exotic properties of excitons associated with BEC and superfluidity. The BEC critical temperature is nonzero, thereby realizing a true condensate with long-range order in 2D homogeneous space. The form of vortex energy has an additional linear term, giving rise to superfluid transition different with BKT type. There is a two-step phase transition with decrease of temperature from a normal exciton phase: at $T_{\mathrm{c1}}=T_{\mathrm{SF}}$ the system first enters a superfluid phase with quasi-long range order, and then at $T_{\mathrm{c2}}=T_{\mathrm{BEC}}$ it enters the BEC phase with both long-range order and superfluid properties. These interesting features can be possibly verified with the recent experimental progress on the exciton fluid in this system, by measuring spatial and temporal coherence of emitted photons.

Acknowledgements. We are grateful to Yi-Cai Zhang for valuable discussions. W.C. and Q.Z. are supported by NKRDPC (Grant No. 2022YFA1405304), NSFC (Grant No. 12004118), and Guangdong Basic and Applied Basic Research Foundation (Grants No. 2020A1515110228 and No. 2021A1515010212). C.-J.H. gratefully acknowledges research support from Gang Chen by the Research Grants Council of Hong Kong with General Research Fund Grant No. 17306520.