Accurate projective two-band description of topological superfluidity in spin-orbit-coupled Fermi gases

The interplay of spin-orbit coupling and Zeeman splitting in ultracold Fermi gases gives rise to a topological superfluid phase in two spatial dimensions that can host exotic Majorana excitations. Theoretical models have so far been based on a four-band Bogoliubov-de Gennes formalism for the combined spin-1/2 and particle-hole degrees of freedom. Here we present a simpler, yet accurate, two-band description based on a well-controlled projection technique that provides a new platform for exploring analogies with chiral p-wave superfluidity and detailed future studies of spatially non-uniform situations.


Introduction
Topological superfluids and superconductors [1,2] are the focus of great current interest because of their ability to host unconventional Majorana excitations [3]. An attractive route towards realization that was suggested early on [4][5][6][7][8] utilizes two-dimensional (2D) s-wave superfluids with spin-orbit coupling. The transition from the non-topological superfluid phase to the topological superfluid phase in these systems is driven by increasing the Zeeman energy splitting between spin-↑ and spin-↓ single-particle states beyond the critical value where the excitation gap for spin-↓ particles closes. The resulting effectively spinless superfluid state is expected to show all the hallmarks associated with chiral p-wave pairing [9], including exotic Majorana states in vortex cores [10][11][12]. Promising experimental efforts are currently undertaken in condensed-matter systems [13][14][15] and ultracold-atom gases [16,17], which have the potential to provide complementary insight and crucial ingredients for topological quantum information devices [18]. One of the important technical differences between these two platforms is the way how the Zeeman spin splitting is introduced. For superconductors, the required magnetic-field strengths are typically deleterious to superconductivity, motivating a search for alternative approaches [19][20][21]. Being unencumbered by this drawback, ultracold-atom realizations may offer a more direct avenue towards implementation of topological superfluidity. Our present study is intended to provide a useful tool for investigating topological effects in superfluid spin-orbit-coupled Fermi gases and ultimately enable the design and optimization of proof-of-concept devices.
Previous theoretical studies [22][23][24][25][26][27][28][29][30][31][32][33][34] of s-wave pairing in Fermi gases with spin-orbit coupling and Zeeman splitting have examined the physical properties of these paradigmatic systems using a mean-field treatment in a four-dimensional Nambu space. While the breaking of spin-rotational invariance generally requires such a more complicated [35], and in general only numerically accessible, treatment, the subspace of the spin-↑ degrees of freedom that are relevant for topological properties is only two-dimensional. See Fig. 1(a) for an illustration. Thus it would be desirable to have an effective description based on projecting into the spin-↑ subspace. However, to discuss manipulations of the system by controllable physical parameters and to make predictions for experimentally accessible observables, the influence of the spin-↓ degrees of freedom cannot be ignored. Here we present a fully self-consistent effective theory that is based on an application of the Feshbach-partitioning technique [36,37] to the spin-resolved Bogoliubov-de Gennes (BdG) Hamiltonian [35,38] under the assumption that the s-wave pair potential |∆| is small compared to the Zeeman energy shift h that favors spin-↑ over spin-↓ configurations. Our effective theory is designed to reproduce salient features of the Bogoliubov-quasiparticle spectrum [ Fig. 1(b)] and all relevant parametric dependences associated with the topological phase transition. This formalism is therefore ideally suited to be a platform for further studies of topological effects, including those associated with non-uniform superfluid phases [39][40][41][42][43]. In order to be specific, and also because it is the physically most interesting case, we develop the theory for a 2D superfluid, but the formalism lends itself to be easily applied to 1D or 3D situations as well.
One of the main benefits associated with having a spin-↑-projected effective theory is that it facilitates the numerical treatment of inhomogeneous and time-dependent situations. For example, in the time-dependent study of soliton or vortex dynamics (e.g., similar to recent work reported in Ref. [40]), the reduction of numerical complexity obtained by moving from the original four-spinor formalism to the projected two-spinor approach is significant. But already for the homogeneous superfluid, the projected-theory results are very useful because, e.g., they provide simpler expressions for the wave functions of the Bogoliubov quasi-particle excitations and thus facilitate convenient analytical approxima- x + k 2 y in the uniform topological-superfluid phase of a 2D Fermi gas with spin-orbit coupling λ k ≡ λ(k x − i k y ) and Zeeman splitting h. In panel (a), the dashed blue (red) curves are the spin-↑ (spin-↓) dispersions in the absence of s-wave pairing (∆ = 0) and spin-orbit coupling (λ = 0) but with large Zeeman splitting h = 2E 0 . A finite ∆ couples spin-↑ and spin-↓ states, resulting in gaps opening at E = ±h and k = 2mµ/ 2 , where µ denotes the chemical potential. In the situation depicted here, µ = E 0 with E 0 ≡ 2 k 2 0 /(2m) being an arbitrary energy scale. When λ is also finite, a third gap opens at E = 0, and the system is a topological superfluid for h > µ 2 + |∆| 2 . The black solid curves are the exact dispersions for 2mλ/( 2 k 0 ) = 0.4 and |∆|/E 0 = 0.8. Panel (b) again depicts these exact dispersions as black solid curves, together with the approximate dispersions obtained by us using a Feshbach projection onto the spin-↑ (spin-↓) subspace shown as the dashed blue (red) curves. tions. As an example, we derive a simple analytic expression for the chemical potential of the spin-orbit-coupled two-dimensional superfluid [see Eq. (20)].
The results presented below can be compared with, and also extend, those of previous studies of superfluidity in Fermi gases with spin-orbit coupling and Zeeman splitting where self-consistency implied fixing the total particle density [23][24][25][26][30][31][32][33][34]. (In contrast, Refs. [27][28][29] consider the situation with fixed chemical potential. See also early work [22] that focused on a lattice realization.) Most relevant benchmarking for our present context is provided by previous works pertaining to uniform 2D systems [32][33][34], but there are also useful connections to be made with known results for trapped 2D systems [30] and 3D systems with 2D Rashba-type spin-orbit coupling [23][24][25][26]31]. In a slight variation on our situation of interest, Ref. [32] considers the population imbalance between spin-↑ and spin-↓ particles as a control parameter rather than the Zeeman energy.
This article is organised as follows. In the following Section 2, we introduce the microscopic model for superfluid 2D Fermi gases with spin-orbit coupling and Zeeman splitting and apply Feshbach partitioning to derive effective theories describing the spin-↑ and spin-↓ sectors separately. The obtained formalism is applied in Section 3 to devise a fully self-contained procedure for finding the chemical potential µ and s-wave pair potential ∆ for uniform systems at fixed total particle density n ≡ n ↑ + n ↓ entirely within the 2 × 2-projected theory for the spin-↑ sector. The efficacy of this approach is demonstrated in Section 4 by presenting a comparison of predictions for phase boundaries and thermodynamic quantities obtained within the effective two-band and exact four-band theories. Following the usual convention, we measure relevant parameters in units of the densitydefined magnitude of the 2D Fermi wave vector k F = √ 2πn and associated Fermi energy E F = 2 k 2 F /(2m) ≡ π 2 n/m, with m denoting the single-particle mass. Our conclusions and an outlook toward future work are presented in the final Section 5.

Feshbach partitioning of the BdG Hamiltonian
Our starting point is the Bogoliubov-de Gennes (BdG) equation describing quasiparticle excitations in a superfluid Fermi gas without spin-rotational invariance. It reads [35,38] with complex spinor entries u σ (v σ ) denoting quantum amplitudes of spin-σ particle (hole) states in a Bogoliubov excitation of the superfluid. (Here and in the following, σ ∈ {↑, ↓} is used as a compact label for the spin degree of freedom.) The Hamiltonian matrix in four-dimensional particle-hole (Nambu) space is where ' * ' indicates complex conjugation, k ≡ (k x , k y ) is the 2D wave vector, and k↑(↓) = k − (+) h with k = 2 (k 2 x + k 2 y )/2m. For spatially inhomogeneous configurations, k j ≡ −i∂ j is to be treated as an operator while, for a homogeneous superfluid, it can be replaced by its wave-number eigenvalue. The Zeeman energy splitting is denoted by h, and λ k is the spin-orbit coupling. Examples of typically considered k-linear spin-orbit couplings are the 2D-Dirac [44], 2D-Rashba [45,46], and 2D-Dresselhaus [46,47] types that correspond to different functional forms λ k = λ(k x − ik y ), λ i(k x − ik y ), and λ(k x + ik y ), respectively, but are all unitarily equivalent. In particular, the eigenvalue spectrum of H for the homogeneous superfluid depends on the spin-orbit coupling only via the quantity |λ k | 2 ≡ λ 2 (k 2 x + k 2 y ) and therefore has the same functional form for all three of the abovementioned spin-orbit-coupling types. The eigenvalue spectrum E kα,η is characterised by four bands of dispersion relations as shown in Fig. 1(a) for a particular set of parameters. In order to be able to refer to a specific band, we introduce the following naming convention (for the fully gapped case), where α = +(−) indicates states that have energy E ≥ 0 (E ≤ 0), whereas the index η = > (<) labels the higher(lower)-energy pair of excitation branches; i.e., |E kα,> | ≥ |E kα,< |. For what follows, it will be useful to know the asymptotic behavior of the dispersions in the limit of large 2D-wave-vector magnitude k ≡ k 2 x + k 2 y . We find The matrix equation (1a) can be reorganized by forming 2 × 2 sub-blocks on the diagonal that are associated with subspaces for fixed spin degree of freedom, with the definitions and ' †' denoting Hermitian conjugation. Simple algebra yields 2 × 2-matrix equations that formally decouple the individual spin sectors, where 1 is the 2 × 2 identity matrix andσ denotes the opposite of σ; i.e.,σ = ↓ (↑) if σ = ↑ (↓). While formally exact and a 2 × 2 BdG-like equation in spin-σ space, Eq. (4b) is not really useful without approximations, since the unknown energy eigenvalue E appears on both sides of the equation. Assuming ∆ to be small compared to other relevant energy scales, we replace E in the denominator by the exact ∆ = 0 solution from Eq. (2) to obtain the approximation Here 2h k = h + h 2 + |λ k | 2 can be thought of as a k-dependent effective Zeeman splitting that is modified by the presence of the spin-orbit coupling. The approximation becomes good if 2h k is large compared to the neglected term, i.e. when |∆| 2 k h k , which is always asymptotically true for large k. Substituting the approximations (5) into (4b) yields truly decoupled BdG equations for the individual spin sectors, where We adopted an anticommutator notation {A , B} = (AB + BA)/2 to be able to incorporate situations with spatially inhomogeneous s-wave pair potential 1 ∆. Notice that the off-diagonal matrix element in (7a) that is responsible for opening the gap in the quasiparticle spectrum is given by λ k ∆/h k , which is proportional to both the spin-orbitcoupling strength λ and ∆. In the case of k-linear spin-orbit coupling and uniform ∆, its leading-order dependence on k resembles the pair potential for a chiral-p-wave superfluid. The emergence of p-wave pairing in the present context was inferred in earlier works [5,7,8,31] through a transformation into the so-called 'helicity' basis of the singleparticle Hamiltonian in the presence of spin-orbit coupling. Although instructive at the time, this transformation does not provide a useful basis for in-depth quantitative studies. In contrast, our present work enables a precise derivation of the effective pairing potential of spin-σ states, including systematic corrections to the chiral-p-wave form. Assuming a spatially uniform pair potential ∆, straightforward diagonalization of (7a) yields approximate energy dispersions and corresponding eigenspinors for the spin-σ subspace as The N σ kα are normalization factors that have to be found from the four-spinor normalization condition 1 = σ w σ † w σ . Using the substitution (4a), this condition translates into separate normalization conditions for the individual spin sectors, Further application of the approximations (5) in (9b) yields Figure 1(b) shows a comparison between the exact dispersions E kα,η , obtained by diagonalizing the original 4 × 4 BdG Hamiltonian (1b), and the approximate results E σ kα from (8a), calculated using the 2×2-subspace projections. While the projected theory does not reproduce the s-wave pairing gaps around E = ±h and k = 2mµ/ 2 , it describes very well the region around the topological gap at E = 0. Most crucially, as it turns out, the dispersions for large k are correctly given by the effective 2 × 2-projected approach. As shown in the following, this enables a faithful description of relevant thermodynamic properties. Further comparisons between the exact 4 × 4-theory dispersions and the 2 × 2projection approximations are explored in Fig. 2. The low-energy gap structure turns out to be well-described even in the non-topological regime, as illustrated in Fig. 2(a). Overall good agreement is achieved in situations when the chemical is negative 2 due to the absence of any crossing points between opposite-spin dispersions [ Fig. 2(b)]. The topological gap ceases to be well-described for quite large spin-orbit-coupling strengths [ Fig. 2 Our Feshbach-projection approach embodied in Eqs. (4) and (5) differs from common perturbative methods such as the Schrieffer-Wolf transformation 3 in two crucial aspects.
The case depicted in Panel (b) is for a topological superfluid having a negative chemical potential, which typically occurs for large two-particle binding energies and/or large spin-orbit-coupling strengths Firstly, the dispersions (8a) derived from our 2×2-projected Bogoliubov-de Gennes Hamiltonians (7a) are well-behaved at all k, whereas those obtained, e.g., from the Schrieffer-Wolf transformation become singular at the s-wave pairing gap because of a degeneracy between the eigenvalues of H ↑↑ and H ↓↓ at this point. Secondly, by careful choice of the approximation (5), we are able to reproduce the large-|k| asymptotics of the exact energy dispersions within the 2×2-projected approach, which cannot be achieved by perturbation theory because it treats the spin-sector couplings given in Eq. (3d) as small quantities.
By construction, the projected-theory results for quasiparticle dispersions and wave functions become strictly exact in the limit of vanishing s-wave pair potential ∆, i.e., when no avoided crossings occur. Thus, for finite ∆, we may expect all quantities that do not explicitly depend on the avoided crossing between the spin-↑ and spin-↓ dispersions to be described correctly as long as |∆| h.

Self-consistency for uniform systems with fixed density
Knowledge of the Bogoliubov-quasiparticle excitations permits calculation of all physical quantities of interest [35,38]. For simplicity, we focus on the zero-temperature limit in the following. Generalization to the case of finite temperatures is straightforward [35,38] but does not add any crucial insights for our present purpose. The pair potential ∆ can be expressed in terms of the eigenspinors of the BdG equation (1a) and the strength g of attractive interactions in the s-wave channel as where Ω denotes the system volume (here: area). As the quasiparticle excitation energies and associated spinor amplitudes are themselves functions of ∆, Eq. (11) constitutes a self-consistency condition [35,38]. However, the expression (11) is formally divergent and needs to be regularized using the relation [48] where E b > 0 is the absolute value of the binding energy of the two-particle bound state in 2D [52][53][54] in the absence of spin-orbit coupling, i.e., for λ = 0. (Modifications of the two-body bound state in a quasi-2D Fermi gas due to spin-orbit coupling are discussed in Refs. [50,55,56], but these are not relevant for the pairing-gap regularisation procedure [57].) The binding energy is related to the 2D scattering length via the expression . is the Euler constant 4 . Recent experimental realizations of low-temperature 2D Fermi gases [61][62][63] have been able to access a wide parameter range −7 ln(k F a 2D ) 4. Combination of Eqs. (11) and (12) yields the practically relevant s-wave pair-potential self-consistency condition The densities n σ of spin-σ particles are also implicit functions of system parameters via the expressions For a uniform Fermi gas with fixed total particle number density n ≡ σ n σ , we thus have a second self-consistency condition given by Explicit analytical expressions for the momentum-space density distributions n kσ defined in (14b) and the quantity entering the r.h.s. of Eq. (13) have been derived within the exact 4 × 4 BdG theory [27,33]. We now discuss in some detail the corresponding results provided by the approximate 2×2projected approach developed here.

Momentum-space density distributions and chemical potential
Using results for the spinor amplitudes given in Eq. (8b), we obtain the momentum-space distribution n kσ of the spin-σ particle density as a sum of contributions from the two projected 2 × 2 sectors, n kσ = n σ kσ + nσ kσ , where 4 Generally, E b ∼ 2 /(ma 2 2D ) for shallow dimers, but values given in the literature for the prefactor on the r.h.s. of that relation vary. This is due to different conventions being used when defining the two-dimensional scattering length a2D [58][59][60]. Our choice follows related previous work [33,34]. This artefact of our approximations is remedied by neglecting n ↓ k↑ (i.e., setting n k↑ ≡ n ↑ k↑ within the 2 × 2-projected theory) from now on. In contrast, n ↑ k↓ is well-behaved [as the denominator of Eq. (17b) for σ = ↓ is always finite] and contributes importantly to n k↓ . Figure 3 shows a comparison between the density distributions n kσ thus obtained within the 2 × 2-projected theory with those calculated within the exact 4 × 4 formalism 5 . There is excellent agreement for the spin-resolved density distributions from both approaches as long as spin-orbit coupling is not too strong. For fixed spin-orbit-coupling strength, deviations are greater for smaller values of the Zeeman splitting h, i.e., these tend to be more pronounced in the non-topological regime.
Interestingly, the concept of separate spin-↓ and spin-↑ Fermi spheres with radius k ↓ and k ↑ , respectively, turns out to be useful even at significant levels of spin-orbit coupling, since the exact density distributions satisfy very accurately the approximate relation as is illustrated in Figs. 3(c) and 3(f). Here Θ(·) denotes the Heaviside step function, k ↑ > k ↓ generically, and k ↓ ≡ 0 in the topological regime. Motivated by observing the apparent broad validity of relation (18), we insert it into the number-density equation (15) and straightforwardly derive the result as an equivalent self-consistency condition. Furthermore, Eq. (2) together with the fact that E k+,< ≈ 0 for |k| = k ↑ (generally valid to leading order in small |∆|) implies For the case where k ↓ = 0, we can extrapolate Eq. (2) to the point E k+,> ≈ 0 when |k| = k ↓ and find In the topological regime (realized for h > h c , where h c denotes the value for the Zeeman energy at the transition), k ↓ = 0 so that Eq. (19) implies k ↑ = √ 2 k F and (20a) yields the approximate relation between chemical potential and particle density. In the non-topological regime (realized for h < h c ), k ↓ = 0 and we need to simultaneously solve Eqs. (20a), (20b) and (19). Adding (20b) to (20a), using (19), and expanding to first sub-leading order in large h, we find Furthermore, subtracting (20b) from (20a) and expanding again to first sub-leading order in large h yields where we have again also used Eq. (19). Combining the results from Eqs. (22a) and (22b), we find the relation between chemical potential and particle number that is valid in the non-topological regime, to leading order in large h. Figure 4 shows a detailed comparison between the self-consistent chemical potential obtained from the exact 4×4 approach, from the effective 2×2-projected theory developed here, and the approximate analytical expressions (21) and (23). The situation depicted in panel (a) is the same as in Fig. 1(b) of Ref. [33]. Capitalizing on the weak |∆| dependence, we used fixed values for the s-wave gap in our calculation, corresponding to the selfconsistent results at h = h c for E b /E F = 0.0462. The curve for the 2 × 2-projected theory is seen to agree very well with the exact result for large-enough h/E F , which includes not only the topological regime but also part of the non-topological regime near the transition. Deviations between the 2×2-projected theory and the exact 4×4 results become significant in the limit of small h, where the Feshbach-projection approach is indeed expected to fail. As illustrated by the situation shown in panel (b), the agreement between the effective-2 × 2 and exact-4 × 4 results becomes excellent for smaller values of λ, reaching also deeper into the non-topological regime. The observation that the exact 4 × 4-theory results for µ/E F and the approximate analytical expressions given in Eqs. (21) and (23) are practically indistinguishable at small-enough magnitude of spin-orbit coupling suggests the possibility to utilize these analytical formulae, for better efficiency and greater insight, as input into the self-consistent calculation of the s-wave pair potential. The analytical expression (23) could also be useful to more accurately represent the chemical potential in the low-h limit where the 2 × 2-projected results deviate significantly from those obtained from the exact 4 × 4 approach.

Self-consistency of s-wave pair potential: Spin-↑-projected theory
The approximate description based on the 2 × 2-subspace projections gave Eq. (8b) for the Nambu-spinor amplitudes. Inserting these expressions into (16)  With k 0 = √ 2 k F , these values coincide with those used for/obtained by numerical calculations whose results are shown in Fig. 1 of Ref. [33]. The energy and momentum scales E 0 and k 0 are related via E 0 ≡ 2 k 2 0 /(2m).
Using our approximation (10) for the normalization factors and, for consistency, replacing Hence, we find that the self-consistency condition for the s-wave pair potential can be formulated entirely in terms of quantities relating to the projected spin-↑ degrees of freedom. Figure 5 shows the k dependence of the quantityΥ k = Υ k − 2/(2 k + E b ) that is the summand in the self-consistency equation (13) for the s-wave pair potential. Parameters are chosen to coincide with those from a recent numerical study [33]. The system is deep in the topological-superfluid regime for panel (a), still topological but close to the transition in panel (b), and a non-topological superfluid close to the transition in panel (c). Perhaps not surprisingly, the agreement between the projected 2 × 2 theory and the exact 4 × 4 formalism is best deep in the topological-superfluid regime, as the fidelity of the projected theory should improve for increasing h. Generally, the small-k and large-k behaviors of Υ k are captured almost perfectly within the projected 2 × 2 theory, with deviations at smaller h occurring chiefly at intermediate values of k. However, as the self-consistency condition (13) involves a sum over all k, the overall effect of such deviations cannot be easily ascertained without explicitly finding the self-consistent s-wave pair potentials within both the 4 × 4 and 2 × 2 approaches. Such a detailed comparison is one of the foci of the next Section.
procedures that are based on a fully explicit knowledge of the exact four-band spectrum E kα,η and the associated eigenspinors. Although such a procedure has the advantage of yielding exact results, its complicated formal structure obscures possibilities for gaining a deeper intuitive understanding of relevant physical effects. In contrast, the formalism developed in Sec. 3 offers the attractive alternative to be able to describe the system entirely in terms of a conceptually simpler theory based on the spin-↑-projected (two-band) spectrum. We now investigate in greater detail the physical picture provided by the effective two-band approach, where the self-consistency condition (13) is solved using Υ k ≡ Υ ↑ k with (25a) and approximating the chemical potential by Eqs. (21) and (23) as appropriate.

Boundary between topological and non-topological phases
We start by considering relevant thermodynamic quantities at the transition between the non-topological and topological superfluid regimes. This transition occurs at the value h ≡ h c of the Zeeman energy that satisfies the condition For a given system with fixed uniform total particle density n and s-wave interaction strength measured in terms of the two-body bound-state energy E b , both µ and |∆| are implicit functions of h and n via the self-consistency conditions and, thus, their values µ c ≡ µ(h c ) and ∆ c ≡ |∆(h c )| are also fixed. In Table 1, we summarize these critical values obtained using the exact 4 × 4 theory and compare with those calculated within the effective 2 × 2 approach using two different methods. To obtain µ 2×2 c and ∆ 2×2 c , we simultaneously solve the self-consistency conditions for the number density and pair potential, Eqs. (15) and (13), assuming also n k↑ ≡ n ↑ k↑ , n k↓ ≡ n ↓ k↓ + n ↑ k↓ , and Υ k ≡ Υ ↑ k with relevant expressions given in Eqs. (17a), (17b), and (25a). In contrast,∆ 2×2 c is the result of a simpler routine where only the pair-potential self-consistency condition (13) is solved, setting Υ k ≡ Υ ↑ k with Eq. (25a) and approximating µ c /E F by Eq. (23). Inspection of Table 1 shows that the values obtained for ∆ 2×2 c and∆ 2×2 c are generally very close, even in the regime where the approximation (23) for µ c /E F is not accurate. [For easy reference, values for µ c /E F and µ 2×2 c /E F that agree to within 5% with the analytical approximation Eq. (23) are indicated in green.] Thus, at least to determine critical values within the effective 2 × 2-projected theory, using the simpler routine yielding∆ 2×2 c is a viable approach. Interestingly, the agreement between values for ∆ c and∆ 2×2 c turns out to be generally better for larger λ. (Values for∆ 2×2 c indicated in magenta are close to within 25% to the exact 4 × 4 results.) More specifically, even though the assumption (23) made for µ when determining∆ c is more broadly valid across the range of accessible E b at smaller λ, the projected 2 × 2-theory's self-consistency equations appear to fail for small |∆|. As a rule of thumb, the condition 2mλ/( 2 k F ) 1 is needed for 2 × 2-theory results to be in reasonable agreement with the exact 4 × 4 values ∆ c . Surprisingly, at larger λ, the rather good agreement between ∆ c and∆ 2×2 c extends even to situations where µ c differs significantly from the approximation (23).
In the regime of small |∆|, for which Fig. 4(b) is an illustration, the approximations Eqs. (21) and (23) are accurate over the entire range of Zeeman energies h, including the critical value h c where both expression yield coinciding values. Thus, from the condition that the right-hand sides of (21) and (23) are equal, we can obtain an approximate expression for the phase boundary in h-λ space, The result (27) is consistent with the expectation that h c ≈ |µ c | for |∆| µ, which follows straightforwardly from (26), in conjunction with the validity of the approximation (23). Figure 6 shows the phase boundary calculated within the projected 2 × 2 theory by solving the self-consistency condition (13) by setting Υ k ≡ Υ ↑ k with (25a) and approximating µ/E F by (23) while also enforcing the relation (26). For comparison, the approximation (27) and exact results obtained from the 4 × 4 formalism are also included in these plots. [The phase boundary found within the 2 × 2-projected theory by simultaneously self-consistent determination of ∆ 2×2 c and µ 2×2 c differs only imperceptibly from the more easily obtained 2×2-theory curve shown in Fig. 6 where µ c /E F is approximated by (23) in the calculation of the critical Zeeman energy.] We restrict ourselves to showing the phase boundary only for intermediate values of 2mλ/( 2 k F ) where superfluidity is not expected to be destabilized by phase separation [27,32].
Interestingly, the approximated 2 × 2 approach turns out to predict the phase bound-   Figure 6: Phase boundary between topological and non-topological superfluid states that occur for h > h c and h < h c , respectively. Results labeled 2 × 2 (4 × 4) were calculated by finding the Zeeman energy satisfying (26) from solution of the self-consistency equation (13) for the s-wave pair potential using Υ k ≡ Υ ↑ k with (25a) and approximating the chemical potential µ by (23) (by simultaneous solution of the exact 4 × 4-theory selfconsistency equations for ∆ and µ), using the indicated values for E b /E F . The latter correspond to ln(k F a 2D ) = 1.00, 2.00, 3.00, respectively (cf. Table 1). Dashed curves show the approximate expression (27) that is expected to be valid in the regime where |∆| µ.
ary between topological and non-topological phases correctly over a broader range of spin-orbit-coupling strengths than naïvely expected when considering the deviations between ∆ c and∆ 2×2 c given in Table 1. This is the result of h c being generally dominated either by the value of µ c or that of ∆ c . Although the 2 × 2-projected theory significantly underestimates ∆ c for small λ, h c is dominated by the chemical potential in this parameter range, in which the expression (23) for µ c is highly accurate. On the other hand, for large-enough λ when ∆ c starts to become more important than µ c for determining h c , the 2 × 2 approach yields quite accurate values for ∆. As a result, the h c (λ) dependence obtained within the projected 2×2 theory faithfully reproduces known qualitative features such as the minimum at 2mλ/( 2 k F ) 1 [32].

Parametric dependences of the self-consistent s-wave pair potential
The magnitude |∆| of the s-wave pair potential depends intricately on the tunable system parameters n, h, λ, and E b through the self-consistency conditions (13) and (15). As it turns out, the dependence on the particle density n is most conveniently absorbed by using the Fermi wave vector k F and Fermi energy E F as units for all other quantities to be measured in. Figure 7 illustrates the λ and h dependence of |∆| and provides a comparison between results obtained within the approximate 2×2 approach and the exact 4×4 theory.
[Both the simplified 2 × 2-theory self-consistency routine where µ c /E F is approximated by (23) and the simultaneously self-consistent determination of µ and |∆| within the 2 × 2projected approach yield practically indistinguishable results for the parameters chosen in the Figure.] We show numbers pertaining to fixed E b /E F = 0.0462, corresponding to ln(k F a 2D ) = 2, to enable direct comparison also with previous works [33,34] that give numerical results for |∆| vs. h 6 .
From the derivation of the main decoupling approximation (5) of the projected approach, we may expect good agreement between the 2 × 2 and 4 × 4 results when |∆| h,  Figure 7: Magnitude |∆| of the s-wave pair potential obtained self-consistently as a function of spin-orbit coupling strength λ and Zeeman energy h. Data points labeled 2 × 2 (4×4) were calculated by solving the self-consistency equation (13) for the s-wave pair potential using Υ k ≡ Υ ↑ k with (25a) and approximating the chemical potential µ by Eqs. (21) and (23)  which is generally supported by the results reported in Fig. 7. It should be noted, though, that obtaining a small |∆| from the self-consistent 2 × 2 theory is not sufficient to guarantee this situation, as can be seen from Fig. 7(a), where |∆|/h > 1 according to the 4 × 4 equations but the accidental compensation of positive and negative parts of the summand (as shown in Fig. 5) results in small values for |∆| within the approximate 2 × 2 theory. In such situations, the projected 2 × 2 approach typically tends to underestimate the value of the self-consistent |∆|. Figure 7(a) [7(b)] shows the λ dependence of |∆| for a situation where the system is in the non-topological [topological] superfluid phase. The same qualitative behavior of |∆| increasing for increased λ is exhibited in both panels (a) and (b), for both the 2 × 2 and 4 × 4 data points. However, the much weaker |∆|-vs.-λ dependence in the non-topological phase is not reproduced correctly by the approximate 2 × 2 formalism, whereas there is quite good agreement with the exact results in the topological phase. This is expected, as the projected 2 × 2 theory should be more accurate at larger h. The situation shown in our Fig. 7(a) corresponds reasonably closely to the case for which |∆| vs. λ is plotted in Figs. 4(a) and 4(b) in Ref. [32] (they have a larger E b /E F and smaller h/E F ), and there is excellent qualitative agreement between their results and ours.
The exact results for the |∆|-vs.-h dependence given in Fig. 7(c) agree with those available from Refs. [33,34]. For small h, deviations between the values calculated within the projected 2 × 2 formalism and the exact 4 × 4 theory are significant, and even the qualitative behavior exhibited by the respective |∆|-vs.-h dependences is seen to be quite different. However, the agreement becomes quite good in the topological regime realized for h > h c = 0.507 E F . In contrast, the projected 2×2 theory is seen to become overall very accurate, even in the non-topological phase, for the larger value of λ for which results are given in Fig. 7(d). Thus, as already indicated by the numbers in Table 1, the effective 2×2 theory yields quantitatively satisfactory results for sufficiently large values of 2mλ/( 2 k F ).

Conclusions and outlook
We have derived an accurate effective theoretical description for superfluidity in 2D Fermi gases with broken spin-rotational invariance due to the presence of spin-orbit coupling and Zeeman spin splitting. Starting from the usually applied self-consistent Bogoliubovde Gennes (BdG) mean-field theory for s-wave pairing in four-dimensional Nambu space [Eq. (1a) with (1b) and (13)], we performed a Feshbach projection onto subspaces associated with fixed spin-σ degrees of freedom [Eqs. (4)]. Using also the approximations given in Eqs. Our subsequent investigations focusing on uniform systems at fixed total particle density have demonstrated that the effective two-band descriptions for individual spin subspaces provide a useful theoretical framework for studying the unusual physical properties of this superfluid, including topological effects. In particular, we found that the effective theory faithfully reproduces the relevant physical aspects of the Bogoliubov-quasiparticle dispersion with the chiral-p-wave-like gap [see Figs. 1(b) and 2] and the occupation-number distribution in reciprocal space (see Fig. 3). For both the dispersions and reciprocal-space density distributions, the projected-2 × 2-theory's accuracy is excellent in the topological regime but generally very good even within a finite range on the non-topological side of the transition. As the Zeeman spin-splitting energy h decreases, so does the accuracy of the projected 2 × 2 approach. This is most apparent in the comparison of the chemical potentials self-consistently obtained within the exact 4 × 4 and approximate 2 × 2 approaches, respectively, shown in Fig. 4. Based on the observation of Fermi-surface-like features in the reciprocal-space occupation-number distribution [Figs. 3(c) and 3(f)], we derived analytical formulae for the chemical potential [Eqs. (21) and (23)] that agree very well with the exact results (see Fig. 4).
The self-consistency condition for the s-wave pair potential within the 2×2 theory turns out to be given entirely in terms of quantities relating to the spin-↑ states [Eq. (13) where Υ k as defined in (16) is replaced by Υ ↑ k from (25a)]. We devise two routines for achieving full self-consistency within the 2 × 2-projected theory. One is based on the simultaneous solution of the self-consistency conditions (13) and (15) using 2 × 2-theory results as input: n k↑ ≡ n ↑ k↑ , n k↓ ≡ n ↓ k↓ + n ↑ k↓ , and Υ k ≡ Υ ↑ k with relevant expressions given in Eqs. (17a), (17b), and (25a). The other, simpler routine solves the self-consistency condition (13) using Υ k ≡ Υ ↑ k as given in Eq. (16) and with the chemical potential approximated by the analytical expressions from Eqs. (21) and (23). For the parameter ranges explored in this work, both routines yield practically indistinguishable results for |∆|, thus making it possible to adopt the simpler routine for further exploration of the physical ramifications of the 2×2-projected theory. Overall, the combination of the projected two-band description for spin-↑ states with the analytical formulae for the chemical potential is seen to provide a reliable theoretical description of the system, with impressive quantitative agreement achieved in the limit of sufficiently large, but entirely realistic, values of the Zeeman splitting and spin-orbit coupling (see Figs. 6, 7 and Table 1).
The ability to utilize an effective two-band (2 × 2) theory for describing superfluidity in the 2D spin-split Fermi gas opens up the opportunity to explore in greater detail suggested analogies with chiral-p-wave pairing [9]. In particular, based on the demonstrated accuracy of the projected 2 × 2 approach for the case of uniform systems, we expect this formalism to also be effective for describing situations with non-uniform order-parameter configurations [39][40][41] or in the presence of disorder [42,43]. These scenarios are interesting because they offer possibilities to create and manipulate exotic Majorana excitations spatially [12,30] or temporally [40]. Future work will address in detail the question of applicability of the 2 × 2 approach in such instances and, as appropriate, apply it to inform the design of basic building blocks for fault-tolerant quantum information processing [3,18].