Coexistence of topological and nontopological Fermi-superfluid phases

The superfluidity of polarized, i.e., species-imbalanced Fermi gases underpins a wide range of topics focused on in current research [1, 2]. Pioneering works were motivated by the desire to understand how Zeeman spin splitting affects s-wave pairing of electrons in metals [3, 4, 5]. Since then, the possibility to tune the attractive two-particle interaction strength in ultracold atom gases using Feshbach resonances [6] has opened up the opportunity to explore the BCS–BEC crossover [7, 8, 9, 10] in polarized Fermi gases, revealing interesting features of their phase diagram [11, 12, 13, 14, 15, 16, 17, 18]. The recent advent of tunable spin-orbit couplings in both condensed-matter and ultracold-atom realizations of the spin-$1/2$ Fermi gas [19, 20] has galvanized interest in studying the combined effects of spin polarization and spin-orbit coupling on s-wave pairing [21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47], especially the emergence of topological superfluidity [48, 49, 50].

Our present study is focused on zero-temperature properties of the two-dimensional (2D) spin-$1/2$ Fermi gas with attractive interactions. The interaction strength can be parameterized in terms of the energy $E_{\mathrm{b}}$ of the two-particle bound state in vacuum without spin-orbit coupling and Zeeman splitting, which exists in a 2D system at any nonzero strength of s-wave interactions [51, 52]. In the absence of spin-orbit coupling, raising the Zeeman energy splitting $2h$ between opposite-spin states drives a first-order transition from the s-wave superfluid phase to the normal phase, regardless of the magnitude of $E_{\mathrm{b}}$ [16]. Introducing spin-orbit coupling drastically changes the $E_{\mathrm{b}}$-$h$ phase diagram [30, 31, 32, 40, 47]: the first-order transition is now between two superfluid phases and occurs only above a critical value $E_{\mathrm{b}}^{\mathrm{(c)}}$. In addition, homogeneous superfluid phases realize a topological superfluid (TSF) when $h$ is larger than a critical value $h_{\mathrm{c}}$ [22, 23, 24] determined by the chemical potential $\mu$ and the s-wave pair-potential magnitude $|\Delta|$ via

These features are illustrated in Fig. 1 for a 2D Fermi gas with fixed total particle density $n$, using the scales $E_{\mathrm{F}}=\pi\hbar^{2}n/m$ and $k_{\mathrm{F}}=\sqrt{2\pi n}$ as energy and wave-vector units, respectively ¹¹1Figure 1 shows results obtained for a system with relatively weak spin-orbit coupling. For comparison, see Fig. 2(a) from Ref. [16] depicting the case without spin-orbit coupling and Fig. 1 from Ref. [47] for situations with moderately large spin-orbit-coupling magnitude.. The shaded region shown in the figure corresponds to parameter combinations for which homogeneous ground-state phases cannot exist. Instead, the ground state of the system in this region will form domains of different coexisting homogeneous phases, similar to the familiar liquid-gas equilibrium in the thermodynamics of real gases [54]. Interestingly, the curve $h_{\mathrm{c}}(E_{\mathrm{b}})$ associated with the topological transition generically crosses into the phase-coexistence region at some binding energy $\geq E_{\mathrm{b}}^{\mathrm{(c)}}$. As a result, the first-order superfluid-superfluid transition becomes intertwined with the topological transition and coexistence now occurs between a nontopological superfluid (NSF) and a TSF [32, 47]. A conceptually different scenario for TSF–NSF coexistence was proposed for trapped Fermi gases [31, 35, 34, 36, 43] where the trap-potential-induced spatial variation of $\mu$ and $\Delta$ causes $h_{\mathrm{c}}$ to be position-dependent, thus creating the possibility for TSF and NSF regions to exist simultaneously within a trap at fixed Zeeman energy $h$.

Here we provide a detailed theoretical analysis of the first-order phase-coexistence region that has previously been mapped in mean-field phase diagrams [32, 47]. In particular, we obtain the chemical potential as well as the pair-potential magnitudes of coexisting superfluid phases as functions of the Zeeman energy $h$, the spin-orbit-coupling strength $\lambda$, and the interaction strength parameterized by $E_{\mathrm{b}}$. Systematic physical insight is gained from approximate analytical expressions that are derived assuming small spin-orbit coupling ($\lambda k_{\mathrm{F}}\ll h$) and validated by comparison with exact numerical results.

Experimental realizations of a TSF are vigorously pursued [49, 55, 56, 57, 58] because such systems host unconventional zero-energy excitations [59, 60, 61, 62, 63] that mimic properties of Majorana fermions [64] and could be used as building blocks for fault-tolerant quantum bits [65]. Such Majorana excitations occur generically at a TSF’s boundaries with vacuum or other nontopological matter [66, 50] and can thus be expected to emerge also at the interface between the TSF coexisting with a NSF in the first-order transition region explored in our present work. The detailed understanding of the system in the coexistence regime developed in this Article can inform a realistic theoretical description [67] of Andreev bound states [68] localized at the TSF–NSF interface. Treating such an unconventional version of a superfluid-superfluid (SS^′) hybrid structure where pair-potential magnitudes are different on opposite sides requires generalization of models applied previously to study Josephson junctions [69, 70, 71, 72, 73], solitons [74, 75, 76] and vortices [77].

We present a detailed study of the first-order phase-coexistence region in 2D Fermi superfluids with Zeeman spin splitting and spin-orbit coupling in Sec. II. Exact numerical and approximate analytical results are presented and compared with each other, followed by a detailed discussion of their physical meaning and implications. Section III contains our conclusions and an outlook on their application. Mathematical details of the derivations yielding our approximate analytical results are provided in the Appendix.

We utilize the spin-resolved version of Bogoliubov–de Gennes (BdG) theory [78, 79] applicable to the description of unconventional [80] or noncentrosymmetric [81, 82] superfluids. It is based on solving the general time-independent BdG equation

to obtain Bogoliubov-quasiparticle energies $E$ and associated Nambu eigenspinors $(u,v)^{T}\equiv(u_{\uparrow},u_{\downarrow},v_{\uparrow},v_{\downarrow})^{T}$ in the representation of 2D position vector ${\bm{\mathrm{r}}}\equiv(x,y)$. Here and in the following, the superscripts ’$\ast$’, ’$T$’, and ’$\dagger$’ are used to denote complex conjugation, transposition in spinor space, and Hermitian conjugation, respectively. The single-particle Hamiltonian $\hat{H}_{0}(\hat{{\bm{\mathrm{k}}}})$ and pair potential $\hat{\Delta}(\hat{{\bm{\mathrm{k}}}})$ are $2\times 2$ matrices in spin-$1/2$ space, which for our system of interest are given by

Here we utilize the standard notation where the $\sigma_{j}$ with $j=x,y,z$ denote Pauli matrices in spin-1/2 space, and $\sigma_{0}$ is the identity operator in that subspace. Furthermore, the combinations $\sigma_{\pm}=(\sigma_{x}\pm i\sigma_{y})/2$ are ladder operators for the eigenstates of $\sigma_{z}$, $\hat{{\bm{\mathrm{k}}}}\equiv(\hat{k}_{x},\hat{k}_{y})=(-i\partial_{x},-i\partial_{y})$ is the operator for the 2D wave vector, $h$ denotes the Zeeman energy, and $\Delta$ is the, in general, complex-number-valued s-wave pair potential. To be specific, we assume a quadratic single-particle dispersion $\epsilon_{\hat{{\bm{\mathrm{k}}}}}=\hbar^{2}(\hat{k}_{x}^{2}+\hat{k}_{y}^{2})/(2m)$ and Rashba-type [83] spin-orbit coupling $\lambda_{\hat{{\bm{\mathrm{k}}}}}=\lambda(i\hat{k}_{x}+\hat{k}_{y})$ in the following, but our results apply more generally.

For a homogeneous system, the Nambu eigenspinors are plane waves in 2D coordinate (${\bm{\mathrm{r}}}$) space, i.e.,

with ${\bm{\mathrm{k}}}$ indicating 2D-wave-vector eigenvalues. The associated quasiparticle eigenenergies are grouped into four bands with dispersions $E_{{\bm{\mathrm{k}}}\alpha,\eta}$ with $\alpha\in\{+,-\}$ and $\eta\in\{<,>\}$, given explicitly by [30, 31]

In terms of these, the grand-canonical ground-state energy density can be obtained via [47]

where $A$ denotes the area occupied by the 2D Fermi gas, and $E_{\mathrm{b}}>0$ is the magnitude of the two-particle bound-state energy in vacuum, for zero spin-orbit coupling and also with no Zeeman spin splitting, introduced already in Sec. I. A complete description of the homogeneous Fermi superfluid at fixed density $n$ requires self-consistent determination of the chemical potential $\mu$ and pair-potential magnitude $|\Delta|$ such that the conditions

are fulfilled, ensuring also that $\mathcal{E}_{\mathrm{gs}}^{\mathrm{(MF)}}(|\Delta|,\mu)$ taken as a function of $|\Delta|$ at fixed $\mu$ has a global minimum [84].

As illustrated in the phase diagram shown in Fig. 1, a homogeneous superfluid ground state is found for the 2D Fermi gas with spin-orbit coupling except for values of $h$ within a region $h_{<}(E_{\mathrm{b}})<h<h_{>}(E_{\mathrm{b}})$. This is because, for $h=h_{<}$ ($h_{>}$), the global minimum of $\mathcal{E}_{\mathrm{gs}}^{\mathrm{(MF)}}(|\Delta|,\mu)$ at the self-consistently determined value of $|\Delta|$ becomes degenerate with a local minimum that is present for $h\lesssim h_{<}$ ($h\gtrsim h_{>}$) and, when $h_{<}<h<h_{>}$, the self-consistent $|\Delta|$ is no longer associated with the global minimum of the ground-state energy as required by thermodynamic stability but rather corresponds to a local minimum or even a maximum. Panels (a)–(c) of Fig. 2 illustrate this behavior of the homogeneous, i.e., single-phase ground-state energy as the Zeeman energy is raised beyond $h_{<}$. The impossibility to realize a homogeneous ground state signals the breakdown of the single-phase description for the system in this parameter range. Instead, coexistence of two phases having the same value of $\mu$ but different densities and pair-potential magnitudes has to be assumed [85], in the spirit of the familiar treatment of the liquid-gas phase transition [54]. Such a first-order phase transition and associated coexistence region occur generically in our system of interest when the two-particle s-wave binding energy $E_{\mathrm{b}}$ exceeds a critical value $E_{\mathrm{b}}^{\mathrm{(c)}}$. As our previous systematic study has shown [47] [see Fig. 8(a) in that article], $E_{\mathrm{b}}^{\mathrm{(c)}}$ increases monotonically with increasing spin-orbit-coupling strength $\lambda$. This explains the ’disappearance’ of the first-order phase transition noted in earlier works [31, 29, 40] when $\lambda$ was increased while keeping $E_{\mathrm{b}}$ fixed.

We proceed to describe details of the theoretical approach and present relevant results. In the phase-coexistence region, the value of the chemical potential is determined via the requirement [85] that the two minima, appearing at pair-potential magnitudes $|\Delta_{\mathrm{s}}|$ and $|\Delta_{\mathrm{w}}|$, respectively, in the ground-state energy density taken as a function of $|\Delta|$ at fixed $\mu$ are degenerate, i.e.,

The condition (8) replaces the one formulated in Eq. (7b) for the homogeneous single-phase ground state. Figure 2(d) depicts the situation that emerges from the one shown in panel (c) after $\mu$ is adjusted to ensure (8) is satisfied. To be specific, we assume $|\Delta_{\mathrm{s}}|>|\Delta_{\mathrm{w}}|$, thus associating $|\Delta_{\mathrm{s}}|$ and $|\Delta_{\mathrm{w}}|$ with the strong and weak superfluids in the coexistence state, respectively. Using these values as input for the homogeneous-single-phase relation (7b) yields the particle densities in the two coexisting phases;

enabling us to determine the proportion $0\leq x\leq 1$ of the weak superfluid in the coexistence state from the modified fixed-density condition [85, 30]

Results obtained from this procedure for three representative parameter combinations are shown in Fig. 3.

The strong superfluid phase turns out to be always nontopological, i.e., $h<\sqrt{\mu^{2}+|\Delta_{\mathrm{s}}|^{2}}$ for $h_{<}<h<h_{>}$. In contrast, the weak superfluid is nontopological only for small-enough $E_{\mathrm{b}}$ and not-too-large $\lambda$ — specifically $\lambda\lesssim 0.7\,E_{\mathrm{F}}/k_{\mathrm{F}}$ [47]. The left column of panels in Fig. 3 [Figs. 3(a,d,g)] depicts such a situation. Conversely, when $E_{\mathrm{b}}/E_{\mathrm{F}}$ is above a critical value [as is the case for the parameter combination used to calculate results presented in the middle column, i.e., panels (b), (e) and (h) of Fig. 3], or when $\lambda\geq 0.7\,E_{\mathrm{F}}/k_{\mathrm{F}}$ [as for Fig. 3(c,f,i)], $h>\sqrt{\mu^{2}+|\Delta_{\mathrm{w}}|^{2}}$ throughout the phase-coexistence region and the weak-superfluid phase is topological.

Across the first order phase transition, we find the density $n_{\mathrm{w}}$ of the weak superfluid phase to be always smaller than the density $n_{\mathrm{s}}$ of the strong-superfluid phase, regardless of whether the weak superfluid phase is a TSF or an NSF. It is thus possible to predict the spatial distribution of coexisting phases in physically realistic situations when the system is subject to finite potential gradients, e.g., because it is in a trap or subject to gravity. According to the local density approximation, an inhomogeneous external potential gives rise to a spatially varying effective chemical potential of local equilibrium, which is maximised at the minimum of the potential. Since the chemical potential is a monotonously growing function of the density (a stability condition of the homogeneous gas), it follows that the high-density strong superfluid phase (always NSF) will accumulate in the central region of a trapping potential under coexistence conditions, while the low-density weak superfluid phase (TSF or NSF) will occupy the shallower region of the potential. This is consistent with the results obtained from the local density approximation in Ref. [31]. The expected layering of phases is indeed the configuration found in numerical mean-field studies of spin-orbit-coupled Fermi superfluids in 2D [34] and 1D [36, 43] traps. In the case that the weak superfluid is topological, the internal boundary between the two coexisting phases constitutes a TSF–NSF interface. More generally, a boundary between the two coexisting phases can be expected to form under any circumstances and, to minimize the associated energy cost [86, 87, 88, 89], will typically have a simply connected shape and shortest-possible length as allowed by sample geometry.

Empirical observation suggests that the $\mu(h)$ curve in the coexistence region exhibits features similar to the familiar Maxwell construction [54]. This observation emerges from a comparison with the $\mu$ values of metastable single-phase states [shown as empty symbols in Figs. 3(a,b,c)] that correspond to only local, not global minima in the $|\Delta|$-dependence of $\mathcal{E}_{\mathrm{gs}}^{\mathrm{(MF)}}(|\Delta|,\mu)$ at fixed $\mu$, such as the one indicated by the empty circle in Fig. 2(c). We explored metastable-state characteristics in our earlier work [47], including the possibility to have more than a single one of these at a given value of $h$ in certain parameter regimes, as is the case for the situation depicted in Figs. 3(b,e). Thus the system lends itself to further experimental and theoretical study of processes akin to superheating and supercooling in gases [54], dynamic bistability, and relaxation into phase coexistence.

The $h$-dependence of the chemical potential is observed to be quite strong in the phase-coexistence regime, interpolating approximately linearly between the typically very different single-phase values $\mu(h_{<})$ and $\mu(h_{>})$. In contrast, the magnitudes of the pair potentials $|\Delta_{\mathrm{s}}|$ and $|\Delta_{\mathrm{w}}|$ vary much more slowly as $h$ is tuned across the phase-coexistence region. Finiteness of $|\Delta_{\mathrm{w}}|$ is a direct consequence of finite spin-orbit coupling, and its magnitude is enhanced monotonically as $\lambda$ increases. Such features and further trends observed in the numerically obtained self-consistent values for the chemical potential and pair-potential magnitudes can be discussed more systematically and quantitatively based on our approximate analytical results. Here we focus on the situation where the first-order phase transition coincides with the topological transition, i.e., when the weak superfluid is a TSF. The middle and right columns of panels in Fig. 3 pertain to examples for such a scenario. For the case of weak spin-orbit coupling that is realized when the condition $\lambda k_{\mathrm{F}}\ll\mathrm{min}\{E_{\mathrm{F}},h\}$ holds for the relevant range of Zeeman energies within the phase-coexistence region, it is possible to derive approximate analytical expressions that generalize those obtained previously [16] when $\lambda=0$. To leading order in explicit small-$\lambda$ corrections, we find (see the Appendix for details of the derivation)

For compactness, we introduced the effective spin-orbit-coupling-renormalized binding energy

that is relevant for the strong-superfluid phase. The form (12) of $E_{\mathrm{b}}^{(\lambda)}$ captures the general tendency of spin-orbit coupling to strengthen pairing [90, 91], but the influence of that on the magnitude of $|\Delta_{\mathrm{s}}|$ is counteracted by a concomitant decrease in $\mu$ as per Eq. (11). Figure 4 illustrates the applicability of the approximate analytical expressions from (11a–c) by comparing them with the exact numerical values for $\mu$, $|\Delta_{\mathrm{s}}|$ and $|\Delta_{\mathrm{w}}|$.

Our results for $\mu$ [Eq. (11)] and the pair-potential magnitude $|\Delta_{\mathrm{s}}|$ of the strong-superfluid phase [Eq. (11b)] exhibit only small corrections due to finite spin-orbit coupling. This is expected because these quantities should smoothly recover known results [16] for the $\lambda=0$ situation. In particular, the dependence of $|\Delta_{\mathrm{s}}|$ on $\lambda$ is extremely weak (as was previously also observed in the $h=0$ limit [90, 91]), and it becomes a function of $h$ only implicitly through its dependence on $\mu$. Thus the expression (11b) together with (11) captures the trend of $|\Delta_{\mathrm{s}}|$ to increase approximately linearly with $h$ while varying only insignificantly with $\lambda$ as long as $\lambda\lesssim\sqrt{E_{\mathrm{F}}h}/k_{\mathrm{F}}$.

As the existence of the weak-superfluid phase is predicated on spin-orbit coupling being finite, $|\Delta_{\mathrm{w}}|$ depends materially on $\lambda$. This is embodied by the analytical expression given in Eq. (11c), which was derived for the situation when the weak superfluid is topological; i.e., for $h>h_{\mathrm{c}}>\mu$, assuming also spin-orbit coupling to be small. According to this formula, which constitutes one of our main results, $|\Delta_{\mathrm{w}}|$ vanishes in a singular fashion for $\lambda\to 0$, as surmised in a previous numerical analysis [31]. The comparatively weak $h$ dependence of $|\Delta_{\mathrm{w}}|$ is also reproduced by the intricate functional form of Eq. (11c).

The formula (11c) relies on the condition $2h>E_{\mathrm{b}}$ being satisfied in the TSF–NSF phase-coexistence region. This is indeed the case because $2h_{<}>E_{\mathrm{b}}$ holds, as can be seen from both the exact numerical results and the approximate analytical expressions for critical fields,

The formulae (13) are generalizations of the previously obtained [16] $\lambda=0$ expressions for the critical fields. Details of our derivation are given in the Appendix, and Fig. 5 illustrates the applicability of Eqs. (13) by a comparison with our exact numerical results. Although the leading corrections to $h_{<}$ and $h_{>}$ due to finite spin-orbit coupling are of the same order, their actual magnitudes turn out to be quite different within the shown parameter range — relevant for $h_{<}$ but insignificant for $h_{>}$.

Our theoretical approach is based on BdG mean-field theory whose applicability to low-dimensional systems is known to be limited [92, 93] to zero temperature and weak-enough interactions. Our choice of parameters is designed to reflect these limitations, ensuring in particular the condition $\mu>0$ for the BCS regime.

We have performed a zero-temperature mean-field study of the first-order superfluid-superfluid transition in a 2D spin-$1/2$ Fermi gas with attractive two-particle interactions in the s-wave channel and also subject to Zeeman spin splitting and spin-orbit coupling. Both numerical and analytical results are presented in Sec. II to characterize the regime where two superfluid phases with different magnitudes of the s-wave pair potential coexist. The total fermion density $n$ is assumed to be fixed, and we absorb density dependences of physical quantities by using $k_{\mathrm{F}}\equiv\sqrt{2\pi n}$ and $E_{\mathrm{F}}\equiv\hbar^{2}k_{\mathrm{F}}^{2}/(2m)$ as wave-vector and energy scales, respectively. Thus, relevant parameters controlling the system properties are the dimensionless spin-orbit-coupling strength $\lambda k_{\mathrm{F}}/E_{\mathrm{F}}$, Zeeman energy $h/E_{\mathrm{F}}$, and the interaction strength parameterized in terms of $E_{\mathrm{b}}/E_{\mathrm{F}}$, where $E_{\mathrm{b}}$ is the two-particle binding energy in vacuum in the absence of spin-orbit coupling and Zeeman spin splitting. Figure 3 shows results for situations where the first-order transition is driven by changing $h/E_{\mathrm{F}}$, keeping all other parameters fixed. We focus on this scenario because increasing $h/E_{\mathrm{F}}$ above a critical value $h_{\mathrm{c}}$ [see Eq. (1)] is associated with the establishment of a topological-superfluid (TSF) phase in our system of interest. Above a critical value for $E_{\mathrm{b}}/E_{\mathrm{F}}$, the topological transition gets intertwined with the first-order phase transition (see Fig. 1), and the TSF coexists with a nontopological superfluid (NSF) over an extended range $h_{<}<h<h_{>}$ of Zeeman energies. The boundary between spatial domains occupied by different phases will generally be shaped by the requirement of energy minimization and/or the form of external potentials.

Coexistence of TSF and NSF phases arises in our system of interest as a consequence of the first-order phase transition driven by raising the Zeeman energy $h$. It is directly linked to the coexistence between superfluid and nonsuperfluid phases in a polarized Fermi gas without spin-orbit coupling. In a different context, TSF–NSF coexistence was proposed to occur in harmonically trapped Fermi gases [31, 35, 34, 36, 43] where both the chemical potential and the pair potential have an $r$-dependent profile, causing a spatially varying critical Zeeman energy $h_{\mathrm{c}}(r)\equiv\sqrt{\mu(r)^{2}+|\Delta(r)|^{2}}$ for the second-order topological transition. When the Zeeman energy $h$ is adjusted to satisfy $h_{\mathrm{c}}(0)>h>h_{\mathrm{c}}(r_{\mathrm{F}})$, where $r_{\mathrm{F}}$ is the trap’s Thomas-Fermi radius, the system consists of a central region where $h<h_{\mathrm{c}}(r)$, which is thus an NSF, surrounded by an outer region where $h>h_{\mathrm{c}}(r)$, thus realizing a TSF. Numerical BdG mean field studies of 2D [34] and 1D [36, 43] trap geometries have elucidated general features of the regime where both TSF and NSF phases are present. There is overall qualitative agreement between our results and those presented in Refs. [34, 36, 43], even though the physical origin of TSF–NSF phase coexistence is ostensibly different in both. However, looking closely at the numerically obtained pair-potential profiles (Fig. 2 in Ref. [36] and Fig. S1 in the Supplemental Material for Ref. [43]), one recognizes a quite sharp drop at the inner edge (the TSF–NSF interface), which is flanked by regions of basically constant pair-potential magnitudes. Such behavior is more evocative of the first-order phase coexistence considered in our work than a system divided into TSF and NSF regions by virtue of a trap-induced continuous variation of both chemical and pair potentials. This observation leads us to surmise that the numerical studies from Refs. [34, 36, 43] encountered the first-order phase-coexistence scenario in a finite-size system.

One of the intriguing characteristics of TSFs is the presence of unconventional quasiparticle excitations at boundaries. The TSF–NSF interface in the coexistence regime discussed in the present work constitutes such a boundary. Our results can feed directly into a theoretical description of the interface modeled as an SS^′ junction where the pair potential changes abruptly between its configuration in the two coexisting homogeneous phases [67]. This approach can provide useful guidance for the experimental investigation and potential manipulation of Andreev bound states at the TSF–NSF interface. Pursuing such avenues for in-depth study within well-controlled ultracold-atom setups can be expected to yield crucial insights about the, at this point, equally intriguing and elusive, Majorana quasiparticles [94, 95, 96, 97].

*