Coexistence of topological and nontopological Fermi-superfluid phases

The two-dimensional spin-imbalanced Fermi gas subject to s-wave pairing and spin-orbit coupling is considered a promising platform for realizing a topological chiral-p-wave superfluid. In the BCS limit of s-wave pairing, i.e., when Cooper pairs are only weakly bound, the system enters the topological phase via a second-order transition driven by increasing the Zeeman spin-splitting energy. Stronger attractive two-particle interactions cause the system to undergo the BCS-BEC crossover, in the course of which the topological transition becomes first-order. As a result, topological and nontopological superfluids coexist in spatially separated domains in an extended region of phase space spanned by the strength of s-wave interactions and the Zeeman energy. Here we investigate this phase-coexistence region theoretically using a zero-temperature mean-field approach. Exact numerical results are presented to illustrate basic physical characteristics of the coexisting phases and to validate an approximate analytical description derived for weak spin-orbit coupling. Besides extending our current understanding of spin-imbalanced superfluid Fermi systems, the present approach also provides a platform for future studies of unconventional Majorana excitations that, according to topology, should be present at the internal interface between coexisting topological and nontopological superfluid parts of the system.

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 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 * uli.zuelicke@vuw.ac.nz states drives a first-order transition from the s-wave superfluid phase to the normal phase, regardless of the magnitude of E b [16]. Introducing spin-orbit coupling drastically changes the E b -h phase diagram [30-32, 40, 47]: the first-order transition is now between two superfluid phases and occurs only above a critical value E (c) b . In addition, homogeneous superfluid phases realize a topological superfluid (TSF) when h is larger than a critical value h c [22][23][24] determined by the chemical potential µ and the s-wave pair-potential magnitude |∆| via h c = µ 2 + |∆| 2 . (1) These features are illustrated in Fig. 1 for a 2D Fermi gas with fixed total particle density n, using the scales E F = π 2 n/m and k F = √ 2πn as energy and wavevector units, respectively [53]. 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 c (E b ) associated with the topological transition generically crosses into the phase-coexistence region at some binding energy ≥ E (c) b . 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,[34][35][36]43] where the trap-potential-induced Mean-field phase diagram for a spin-orbit-coupled 2D Fermi gas at zero temperature and fixed density n = mEF/(π 2 ) ≡ k 2 F /(2π) in the parameter space of attractiveinteraction strength (quantified in terms of the two-particle binding energy E b in vacuum without spin-orbit coupling and Zeeman spin splitting) and Zeeman energy h. The blue and green curves correspond to the critical Zeeman energies h< and h> that delimit the first-order phase-transition region h< < h < h>, indicated by light-blue shading. In this region, homogeneous ground-state phases are not possible and coexistence between domains of different superfluid phases is expected. The second-order transition between nontopologicalsuperfluid (NSF) and topological-superfluid (TSF) phases occurs at the red curve where h = hc. Orange dashed lines illustrate the parameter range for which further data are shown in Fig. 3. Results presented in the figure have been calculated for the spin-orbit-coupling strength λ = 0.35 EF/kF. spatial variation of µ and ∆ causes h c to be positiondependent, 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-orbitcoupling strength λ, and the interaction strength parameterized by E b . Systematic physical insight is gained from approximate analytical expressions that are derived assuming small spin-orbit coupling (λk F 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 [50,66] and can thus be expected to emerge also at the interface between the TSF coexisting with a NSF in the firstorder 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 phasecoexistence 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.

II. FIRST-ORDER SUPERFLUID-SUPERFLUID TRANSITION IN A SPIN-IMBALANCED 2D FERMI GAS WITH SPIN-ORBIT COUPLING
We utilize the spin-resolved version of Bogoliubovde 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 ≡ (u ↑ , u ↓ , v ↑ , v ↓ ) T in the representation of 2D position vector r ≡ (x, y). Here and in the following, the superscripts ' * ', 'T ', and ' †' are used to denote complex conjugation, transposition in spinor space, and Hermitian conjugation, respectively. The single-particle HamiltonianĤ 0 (k) and pair potential ∆(k) are 2 × 2 matrices in spin-1/2 space, which for our system of interest are given bŷ Here we utilize the standard notation where the σ j with j = x, y, z denote Pauli matrices in spin-1/2 space, and σ 0 is the identity operator in that subspace. Furthermore, the combinations σ ± = (σ x ± iσ y )/2 are ladder operators for the eigenstates of σ z ,k ≡ (k x ,k y ) = (−i∂ x , −i∂ y ) is the operator for the 2D wave vector, h denotes the Zeeman energy, and ∆ is the, in general, complex-number-valued s-wave pair potential. To be specific, we assume a quadratic single-particle dispersion k = 2 (k 2 x +k 2 y )/(2m) and Rashba-type [83] spin-orbit coupling λk = λ(ik x +k y ) in the following, but our results apply more generally. For a homogeneous system, the Nambu eigenspinors are plane waves in 2D coordinate (r) space, i.e., with k indicating 2D-wave-vector eigenvalues. The associated quasiparticle eigenenergies are grouped into four bands with dispersions E kα,η with α ∈ {+, −} and η ∈ {<, >}, 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 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 selfconsistent determination of the chemical potential µ and pair-potential magnitude |∆| such that the conditions are fulfilled, ensuring also that E (MF) gs (|∆|, µ) taken as a function of |∆| at fixed µ 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 val- (|∆|, µ) at the self-consistently determined value of |∆| becomes degenerate with a local minimum that is present for h h < (h h > ) and, when h < < h < h > , the self-consistent |∆| 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 µ 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 b exceeds a critical value E (|∆|, µ), taken as a function of the pairpotential magnitude |∆| for fixed chemical potential µ, when the Zeeman energy h is increased beyond h<. For h h< [panel (a)], the global minimum (marked by the filled black circle) is at the value of |∆| that satisfies the fixed-density condition Eq. (7b) with the given value of µ, thus representing a valid single-phase equilibrium ground state of the system. As h is increased, the energy difference between the global minimum and another (local) minimum (marked with a cross) reduces and, for h = h<, both minima are degenerate [panel (b)]. For h< < h < h>, the minimum associated with the self-consistent value of |∆| [indicated by the empty black circle in panel (c)] then ceases to be the global minimum. It therefore cannot any longer represent a viable equilibrium ground state but remains a possible metastable state. Adjusting µ to maintain degeneracy between the two minima of E strength λ. This explains the 'disappearance' of the firstorder phase transition noted in earlier works [29,31,40] when λ was increased while keeping E b fixed.
We proceed to describe details of the theoretical approach and present relevant results. In the phasecoexistence region, the value of the chemical potential is determined via the requirement [85] that the two minima, appearing at pair-potential magnitudes |∆ s | and |∆ w |, respectively, in the ground-state energy density taken as a function of |∆| at fixed µ are degenerate, i.e., The condition (8) (8) is satisfied. To be specific, we assume |∆ s | > |∆ w |, thus associating |∆ s | and |∆ 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 ≤ x ≤ 1 of the weak superfluid in the coexistence state from the modified fixed-density condition [30,85] 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 < µ 2 + |∆ s | 2 for h < < h < h > . In contrast, the weak superfluid is nontopological only for small-enough E b and not-too-large λ -specifically λ 0.7 E F /k F [47]. The left column of panels in Fig. 3 [Figs. 3(a,d,g)] depicts such a situation. Conversely, when E b /E 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 λ ≥ 0.7 E F /k F [as for Fig. 3(c,f,i)], h > µ 2 + |∆ 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 w of the weak superfluid phase to be always smaller than the density n 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 µ(h) curve in the coexistence region exhibits features similar to the familiar Maxwell construction [54]. This observation emerges from a comparison with the µ 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 |∆|-dependence of E (MF) gs (|∆|, µ) at fixed µ, 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 µ(h < ) and µ(h > ). In contrast, the magnitudes of the pair potentials |∆ s | and |∆ w | vary much more slowly as h is tuned across the phase-coexistence region. Finiteness of |∆ w | is a direct consequence of finite spin-orbit coupling, and its magnitude is enhanced monotonically as λ 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 λk F min{E F , h} holds for the relevant range of Zeeman energies within the phase-coexistence region, it is possible to derive approximate analytical expres- , with EF ≡ 2 k 2 F /(2m) = π 2 n/m. Empty circles (squares, triangles) indicate values for µ and |∆| that represent self-consistent homogeneous nontopological (topological, critical) superfluid states at the given h but are not associated with a global minimum of the ground-state energy, such as in the situation depicted in Fig. 2(c). These states are metastable and, hence, will not be realized in equilibrium.
sions that generalize those obtained previously [16] when λ = 0. To leading order in explicit small-λ corrections, we find (see the Appendix for details of the derivation) .
For compactness, we introduced the effective spin-orbitcoupling-renormalized binding energy that is relevant for the strong-superfluid phase. The form (12) captures the general tendency of spin-orbit coupling to strengthen pairing [90,91], but the influence of that on the magnitude of |∆ s | is counteracted by a concomitant decrease in µ as per Eq. (11a). Figure 4 illustrates the applicability of the approximate analytical expressions from (11a-c) by comparing them with the exact numerical values for µ, |∆ s | and |∆ w |.
Our results for µ [Eq. (11a)] and the pair-potential magnitude |∆ 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 λ = 0 situation. In particular, the dependence of |∆ s | on λ 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 µ. Thus the expression (11b) together with (11a) captures the trend of |∆ s | to increase approximately linearly with h while varying only insignificantly with λ as long as λ √ E F h/k F . As the existence of the weak-superfluid phase is predicated on spin-orbit coupling being finite, |∆ w | depends materially on λ. 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 c > µ, assuming also spin-orbit coupling to be small. According to this formula, which constitutes one of our main results, |∆ w | vanishes in a singular fashion for λ → 0, as surmised in a previous numerical analysis [31]. The comparatively weak h dependence of |∆ w | is also reproduced by the intricate functional form of Eq. (11c).
The formula (11c) relies on the condition 2h > E b being satisfied in the TSF-NSF phase-coexistence region. This is indeed the case because 2h < > E b holds, as can be seen from both the exact numerical results and the approximate analytical expressions for critical fields, (1). Phase coexistence is between a TSF and an NSF when h> > hc. This is the regime for which the analytical formulae were derived and where the agreement is excellent. The inset shows results for a reduced parameter range together with the λ = 0 critical fields plotted as the short-dashed blue and green curves.
The formulae (13) are generalizations of the previously obtained [16] λ = 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 µ > 0 for the BCS regime.

III. CONCLUSIONS
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 F ≡ √ 2πn and E F ≡ 2 k 2 F /(2m) as wave-vector and energy scales, respectively. Thus, relevant parameters controlling the system properties are the dimensionless spin-orbit-coupling strength λk F /E F , Zeeman energy h/E F , and the interaction strength parameterized in terms of E b /E F , where E 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 F , keeping all other parameters fixed. We focus on this scenario because increasing h/E F above a critical value h 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 b /E F , the topological transition gets intertwined with the firstorder 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,[34][35][36]43] where both the chemical potential and the pair potential have an r-dependent profile, causing a spatially varying critical Zeeman energy h c (r) ≡ µ(r) 2 + |∆(r)| 2 for the second-order topological transition. When the Zeeman energy h is adjusted to satisfy h c (0) > h > h c (r F ), where r F is the trap's Thomas-Fermi radius, the system consists of a central region where h < h c (r), which is thus an NSF, surrounded by an outer region where h > h 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].

ACKNOWLEDGMENTS
This work was supported by the Marsden Fund Council (contract no. VUW1713), from New Zealand government funding managed by the Royal Society Te Apārangi. The authors thank Philip Brydon for useful comments.
Appendix: Derivation of approximate analytical expressions for chemical potential and pair-potential magnitudes in the TSF-NSF-coexistence regime Our analytical approach is based on approximating the true mean-field ground-state energy E (MF) gs (|∆|, µ) in the phase-coexistence region by expressions E (s) gs (|∆|, µ) and E (w) gs (|∆|, µ) that accurately capture the functional dependence on the pair-potential magnitude around its degenerate minima |∆ s | and |∆ w |, respectively. We also consider only the weak-spin-orbit-coupling limit, which means practically that nonleading-order corrections in λk F /h and λk F /E F are neglected.
We start by determining E (s) gs (|∆|, µ), which describes the strong-superfluid part of the system, i.e., the one with the larger pair-potential magnitude |∆ s |. As phase coexistence requires a sufficiently large E b ≥ E (c) b [47], it turns out that the condition h < |∆ s | holds in situations where the first-order transition also constitutes the topological transition. Using this condition alongside the one for small λ in calculating the expression (6) for the meanfield ground-state energy, we obtain to leading order in both small quantities (λk F /E F as well as h/|∆|) (A.1) In the limit λ → 0, Eq. (A.1) reduces to the familiar h-independent ground-state energy for the 2D Fermi gas with finite Zeeman spin splitting when it is in its unpolarized superfluid phase [16,51]. For h → 0, (A.1) also recovers the previously obtained [90] leading-order correction to the ground-state energy due to spin-orbit coupling. For the purposes of our following calculations, the term ∝ h in E (s) gs (|∆|, µ) will be neglected, as it constitutes a further small correction to a contribution that is already small in the parameter λk F /E F .
We now proceed to obtain E (w) gs (|∆|, µ), which has its minimum at the pair-potential magnitude |∆ w | of the weaksuperfluid part of the phase-coexistence state. In the absence of spin-orbit coupling, this phase would be a normal (nonsuperfluid) Fermi gas [16]. This implies that |∆ w | vanishes for λ → 0 and, thus, will generally be small in the weak-spin-orbit-coupling limit. Focusing on situations where the weak superfluid phase is topological, we assume h > µ 2 + |∆| 2 together with λk F h when evaluating the r.h.s. of Eq. (6), which yields the result The result (A.2) generalises the formula for the groundstate energy of a fully polarised (h > µ) 2D Fermi gas to the situation of having a finite |∆| due to λ = 0. The expressions for |∆ s,w | are determined from the minimum-ground-state-energy conditions ∂ E can also be utilized to find the chemical potentials µ s and µ w in the homogeneous nontopological (strong) and topological (weak) superfluid phases that adjoin the phase-separation region. These can then be used to find analytical expressions for the critical Zeeman energies h < and h > , by equating µ as given in Eq. (11a) with µ s and µ w , respectively. For implementing this procedure [16], we apply the fixeddensity conditions ∂ E gs in (A.6), inserting Eq. (11b) for the pair-potential magnitude, and collecting leadingorder terms in the spin-orbit-coupling strength yields the h-independent result [40,46] that coincides with the expression derived previously for the spin-orbit-coupled 2D superfluid at h = 0 [90]. Turning to determining µ w , we use Eq. (A.2) for E (w) gs in (A.6) and insert Eq. (11c) for |∆|. Keeping only leadingorder λ-dependent corrections (which implies neglecting any terms with exponentials that vanish faster than any power law as λ → 0), we find consistent with previous results given for the chemical potential of the topological 2D Fermi superfluid [46]. Equating µ s from Eq. (A.7) and µ w from Eq. (A.8), respectively, with µ from Eq. (11a) yields Eqs. (13) for the critical Zeeman energies. Figure 5 illustrates the validity of the analytical description within its region of applicability, i.e., when the weak superfluid phase is a TSF.