Artificial spin-orbit coupling in ultra-cold Fermi superfluids

The effects of spin-orbit coupling in few body systems like the Earth-Moon complex in astrophysics or the electron spin and its orbital motion around the nucleus in isolated atoms of atomic physics are reasonably well understood due to the simplificity of these systems. However, in the setting of many identical particles, spin-orbit effects have revealed quite interesting surprises recently running from topological insulators in semiconductors kane-2005 to exotic superconductivity gorkov-2001 and non-equillibrium effects galitski-2007 depending on the precise form of the spin-orbit coupling. In atomic physics the coupling arises from the interaction of the magnetic moment of the electron and a magnetic field, present in the frame of electron, due to the electric field of the nucleus. Similarly in condensed matter physics, the coupling arises from the magnetic moment ${\bf m}$ of electrons, which move in the background of ions. In the electron’s reference frame, these ions are responsible for a magnetic field ${\bf B}$, which depends on the electron’s momentum ${\bf k}$ and couple to electron’s spin. The resulting spin-orbit coupling has the form ${H}_{SO}=-{\bf m}\cdot{\bf B}=-\sum_{j}h_{j}({\bf k})\sigma_{j},$ where $\sigma_{j}$ represents the Pauli matrices and $h_{j}({\bf k})$ describes the $j$-th component $(j=x,y,z)$ of the effective magnetic field vector ${\bf h}$. For some materials ${\bf h}$ can take the Dresselhaus dresselhaus-1955 form ${\bf h}_{D}({\bf k})=v_{D}(k_{y}{\hat{\bf x}}+k_{x}{\hat{\bf y}}),$ the Rashba rashba-1984 form ${\bf h}_{R}({\bf k})=v_{R}(-k_{y}{\hat{\bf x}}+k_{x}{\hat{\bf y}}),$ or more generally a linear combination of the two ${\bf h}_{\perp}({\bf k})={\bf h}_{D}({\bf k})+{\bf h}_{R}({\bf k}).$ In all these situations the type of spin-orbit coupling can not be changed arbitrarily and the magnitude can not be tuned from weak to strong, making the experimental control of spin-orbit effects very difficult.

Recently, however, it has been demonstrated experimentally that spin-orbit coupling can be engineered in a ultra-cold gas of bosonic atoms in their Bose-Einstein condensate phase spielman-2011 , when a pair of Raman lasers creates a coupling between two internal spin states of the atoms and its center-of-mass motion (momentum). Thus far, the type of spin-orbit field that has been created in the laboratory spielman-2011 has the equal-Rashba-Dresselhaus (ERD) form ${\bf h}_{\perp}({\bf k})={\bf h}_{ERD}({\bf k})=v{k_{x}}{\hat{\bf y}}$, where $v_{R}=v_{D}=v/2$. Other forms of spin-orbit fields require additional lasers and create further experimental difficulties dalibard-2010 . In ultra-cold bosons the momentum-dependent ERD coupling has been created in conjunction with uniform Zeeman terms, which are independent of momentum, along the z axis (controlled by the Raman coupling $\Omega_{R}$), and along the y-axis (controled by the detuning $\delta$). The simultaneous presence of $h_{z}$, $h_{y}$ and ${\bf h}_{ERD}({\bf k})$ leads to the Zeeman-spin-orbit (ZSO) Hamiltonian

for an atom with center-of-mass momentum ${\bf k}$ and spin basis $|\uparrow\rangle$, $|\downarrow\rangle$. The fields $h_{z}=-\Omega_{R}/2$, $h_{y}=-\delta/2$ and ${\bf h}_{ERD}=vk_{x}{\hat{y}}$ can be controlled independently, and thus can be used as tunable parameters to explore the available phase space and to investigate phase transitions, as achieved in the experiment involving a bosonic isotope of Rubidium (⁸⁷Rb). Although current experiments have focused on Bose atoms, there is no fundamental reason that impeeds the realization of a similar set up for Fermi atoms spielman-2011 ; sinova-2009 ; chapman-sademelo-2011 designed to study fermionic superfluidity chapman-sademelo-2011 . Considering possible experiments with fermionic atoms such as ⁶Li, ⁴⁰K, we discuss in this letter phase diagrams, topological phase transitions, spectroscopic and thermodynamic properties at zero and finite temperatures during the evolution from BCS to BEC superfluidity in the presence of controllable Zeeman and spin-orbit fields in three dimensions.

To investigate artificial spin-orbit and Zeeman fields in ultra-cold Fermi superfluids, we start from the Hamiltonian density

where the single-particle term is simply

Here, ${\hat{K}}_{\alpha}=-\nabla^{2}/(2m)-\mu_{\alpha}$ is the kinetic energy in reference to the chemical potential $\mu_{\alpha}$, ${\hat{h}}_{j}({\bf r})$ is the combined effective field including Zeeman and spin-orbit components along the $j$-direction $(j=x,y,z)$, and $\psi^{\dagger}_{\alpha}({\bf r})$ are creation operators for fermions with spin $\alpha$ at position ${\bf r}$. Notice that we allow the chemical potential $\mu_{\uparrow}$ to be different from $\mu_{\downarrow}$, such that the number of fermions $N_{\uparrow}$ with spin $\uparrow$ may be different from the number of fermions with spin $\downarrow$. The interaction term is

where $g$ represents a contact interaction that can be expressed in terms of the scattering length via the Lippman-Schwinger relation $V/g=-Vm/(4\pi a_{s})+\sum_{\bf k}1/(2\epsilon_{\bf k}).$ The introduction of the average pairing field $\Delta({\bf r})\equiv g\langle\psi_{\downarrow}({\bf r})\psi_{\uparrow}({\bf r})\rangle\approx\Delta_{0}$ and its spatio-temporal fluctuation $\eta({\bf r},\tau)$ produce a complete theory for superfluidity in this system.

From now on, we focus on the experimental case where a) the Raman detuning is zero $(\delta=0)$ indicating that there is no component of the Zeeman field along the $y$ direction; b) the Raman coupling $\Omega_{R}$ is non-zero meaning that a Zeeman component along the $z$ direction exists, that is, $h_{z}=-\Omega_{R}/2$; and c) the spin-orbit field has components $h_{y}({\bf k})$ and $h_{x}({\bf k})$ along the $y$ and $x$ directions. To start our discussion, we neglect fluctuations, and transform $H_{0}({\bf r})$ into momentum space as $H_{0}({\bf k})$. Using the basis $\psi_{\uparrow}^{\dagger}({\bf k})|0\rangle\equiv|{\bf k}\uparrow\rangle,$ $\psi_{\downarrow}^{\dagger}({\bf k})|0\rangle\equiv|{\bf k}\downarrow\rangle,$ where $|0\rangle$ is the vacuum state, the Fourier-transformed Hamiltonian $H_{0}({\bf k})$ becomes the matrix

Such matrix can be diagonalized in the helicity basis $\Phi_{\Uparrow}^{\dagger}({\bf k})|0\rangle\equiv|{\bf k}\Uparrow\rangle,$ $\Phi_{\Downarrow}^{\dagger}({\bf k})|0\rangle\equiv|{\bf k}\Downarrow\rangle,$ where the spins $\Uparrow$ and $\Downarrow$ are aligned or antialigned with respect to the effective magnetic field ${\bf h}_{\rm eff}({\bf k})={\bf h}_{\parallel}({\bf k})+{\bf h}_{\perp}({\bf k}).$ Here, $K_{+}({\bf k})=(K_{\uparrow}+K_{\downarrow})/2=\epsilon_{\bf k}-\mu_{+},$ is a measure of the average kinetic energy $\epsilon_{\bf k}=k^{2}/2m$ in relation to the average chemical potential $\mu_{+}=(\mu_{\uparrow}+\mu_{\downarrow})/2.$ While ${\bf h}_{\perp}({\bf k})=h_{x}({\bf k})\hat{\bf x}+h_{y}({\bf k})\hat{\bf y}$ is the spin-orbit field and ${\bf h}_{\parallel}({\bf k})=(h_{z}-K_{-})\hat{\bf z}$ is the effective Zeeman field, with $K_{-}=(K_{\uparrow}-K_{\downarrow})/2=-\mu_{-}$ where $\mu_{-}=(\mu_{\uparrow}-\mu_{\downarrow})/2$ is the internal Zeeman field due to initial population imbalance, and $h_{z}$ is the external Zeeman field. When there is no population imbalance the internal Zeeman field is $\mu_{-}=0$, and we have only $h_{z}$. In general, the eigenvalues of the Hamiltonian matrix ${\bf H}_{0}({\bf k})$ are $\xi_{\Uparrow}({\bf k})=K_{+}({\bf k})-|{\bf h}_{\rm eff}({\bf k})|$ and $\xi_{\Downarrow}({\bf k})=K_{+}({\bf k})+|{\bf h}_{\rm eff}({\bf k})|,$ where $|{\bf h}_{\rm eff}({\bf k})|=\sqrt{(\mu_{-}+h_{z})^{2}+|h_{\perp}({\bf k})|^{2}}$ is the magnitude of the effective magnetic field, with the transverse component being expressed in terms of the complex function $h_{\perp}({\bf k})=h_{x}({\bf k})+ih_{y}({\bf k}).$ In the limit where the internal $\mu_{-}$ and external $h_{z}$ Zeeman fields vanish and the spin-orbit field is null $(h_{\perp}=0)$, the energies of the helicity bands are identical $\xi_{\Uparrow}({\bf k})=\xi_{\Downarrow}({\bf k})$ producing no effect in the original energy dispersions shenoy-2011 .

When interactions are added to the problem, pairing can occur within the same helicity band (intra-helicity pairing) or between two different helicity bands (inter-helicity pairing). This leads to a tensor order parameter for superfluidity that has four components $\Delta_{\Uparrow\Uparrow}({\bf k})=-\Delta_{T}({\bf k})e^{-i\varphi},$ corresponding to the helicity projection $\lambda=+1$; $\Delta_{\Uparrow\Downarrow}({\bf k})=-\Delta_{S}({\bf k}),$ and $\Delta_{\Downarrow\Uparrow}({\bf k})=\Delta_{S}({\bf k}),$ corresponding to helicity projection $\lambda=0$; and $\Delta_{\Downarrow\Downarrow}({\bf k})=-\Delta_{T}({\bf k})e^{i\varphi},$ corresponding to helicity projection $\lambda=-1$. The phase $\varphi({\bf k})$ is defined from the amplitude-phase representation of the complex spin-orbit field $h_{\perp}({\bf k})=|h_{\perp}({\bf k})|e^{i\varphi({\bf k})},$ while the amplitude $\Delta_{T}({\bf k})=\Delta_{0}|h_{\perp}({\bf k})|/|{\bf h}_{\rm eff}({\bf k})|$ for helicities $\lambda=\pm 1$ are directly proportional to the scalar order parameter $\Delta_{0}$ and to the relative magnitude of the spin-orbit field $|h_{\perp}({\bf k})|$ with respect to the magnitude of the effective magnetic field $|{\bf h}_{\rm eff}({\bf k})|$. Additionally, $\Delta_{T}$ has the simple physical interpretation of being the triplet component of the order parameter in the helicity basis, which is induced by the presence of non-zero spin-orbit field $h_{\perp}$, but vanishes when $h_{\perp}=0$. Analogously the amplitude $\Delta_{S}({\bf k})=\Delta_{0}h_{\parallel}({\bf k})/|{\bf h}_{\rm eff}({\bf k})|$ for helicity $\lambda=0$ are directly proportional to the scalar order parameter $\Delta_{0}$ and to the relative magnitude of the total Zeeman field $h_{\parallel}({\bf k})=\mu_{-}+h_{z}$ with respect to the magnitude of the effective magnetic field $|{\bf h}_{\rm eff}({\bf k})|$. Additionally, $\Delta_{S}$ has the simple physical interpretation of being the singlet component of the order parameter in the helicity basis. It is interesting to note the relation $|\Delta_{T}({\bf k})|^{2}+|\Delta_{S}({\bf k})|^{2}=|\Delta_{0}|^{2},$ which, for fixed $|\Delta_{0}|$, shows that as $|\Delta_{S}({\bf k})|$ increases, $|\Delta_{T}({\bf k})|$ decreases and vice-versa. Such relation indicates that the singlet and triplet channels are not separable in the presence of spin-orbit coupling. Furthermore, the order parameter in the triplet sector $\Delta_{\Uparrow\Uparrow}({\bf k})$ and $\Delta_{\Downarrow\Downarrow}({\bf k})$ contains not only $p$-wave, but also $f$-wave and even higher odd angular momentum contributions, as long as the total Zeeman field $\mu_{-}+h_{z}$ is non-zero. Similarly, the order parameter in the singlet sector $\Delta_{\Uparrow\Downarrow}({\bf k})$ and $\Delta_{\Downarrow\Uparrow}({\bf k})$ contains not only only $s$-wave, but also $d$-wave and even higher even angular momentum contributions, as long as the total Zeeman field $\mu_{-}+h_{z}$ is non-zero. Higher angular momentum pairing in the helicity basis, occurs because the original local (zero-ranged) interaction in the original $(\uparrow,\downarrow)$ spin basis is transformed into a finite-ranged interaction in the helicity basis $(\Uparrow,\Downarrow)$. In the limiting case of zero total Zeeman field $\mu_{-}+h_{z}=0$, the singlet component vanishes $(\Delta_{S}({\bf k})=0)$, while the triplet component becomes independent of momentum $(\Delta_{T}({\bf k})=\Delta_{0})$, leading to order parameter $\Delta_{\Uparrow\Uparrow}({\bf k})=-h_{\perp}^{*}({\bf k})$, and $\Delta_{\Downarrow\Downarrow}({\bf k})=-h_{\perp}({\bf k})$ which contains only $p$-wave contributions chuanwei-2008 , since the components of $h_{\perp}({\bf k})$ depend linearly on momentum ${\bf k}$.

The eigenvalues $E_{j}({\bf k})$ of the Hamiltonian including the order parameter contribution emerge from the diagonalization of a $4\times 4$ matrix (see supplementary material). The two eigenvalues for quasiparticles are

corresponding to the highest-energy quasiparticle band, and

corresponding to the lowest-energy quasiparticle band, while the eigenvalues for quasiholes are $E_{3}({\bf k})=-E_{2}({\bf k})$ for highest-energy quasihole band and $E_{4}({\bf k})=-E_{1}({\bf k})$ for the lowest-energy quasihole band. The energy $\xi_{h_{-}}=\left[\xi_{\Uparrow}({\bf k})-\xi_{\Downarrow}({\bf k})\right]/2$ is momentum-dependent, corresponds to the average energy difference between the helicity bands and can be written as $\xi_{h_{-}}=-|{\bf h}_{\rm eff}({\bf k})|,$ while the energy $\xi_{h_{+}}=\left[\xi_{\Uparrow}({\bf k})+\xi_{\Downarrow}({\bf k})\right]/2$ is also momentum dependent, corresponds to the averaged energy sum of the helicity bands and can be written as $\xi_{h_{+}}=K_{+}({\bf k})=\epsilon_{\bf k}-\mu_{+}.$

There are a few important points to notice about the excitation spectrum of this system. First, notice that $E_{1}({\bf k})>E_{2}({\bf k})\geq 0$. Second, that the eigenergies are symmetric about zero, such that we can regard quasiholes (negative energy solutions) as anti-quasiparticles. Third, that only $E_{2}({\bf k})$ can have zeros (nodal regions) corresponding to the locus in momentum space satisfying the following conditions: a) $\xi_{h_{-}}=-\sqrt{\xi_{h_{+}}^{2}+|\Delta_{S}({\bf k})|^{2}},$ which corresponds physically to the equality between the effective magnetic field energy $|{\bf h}_{\rm eff}({\bf k})|$ and the excitation energy for the singlet component $\sqrt{\xi_{h_{+}}^{2}+|\Delta_{S}({\bf k})|^{2}};$ and b) $|\Delta_{T}({\bf k})|=0,$ corresponding to zeros of the triplet component of the order parameter in momentum space.

Since $E_{2}({\bf k})<E_{1}({\bf k})$, and only $E_{2}({\bf k})$ can have zeros, the low energy physics is dominated by this eingenvalue. In the case of equal Rashba-Dresselhaus (ERD) where $h_{\perp}({\bf k})=v|k_{x}|$, zeros of $E_{2}({\bf k})$ can occur when $k_{x}=0$, leading to the following cases: (a) two possible lines (rings) of nodes at $(k_{y}^{2}+k_{z}^{2})/(2m)=\mu_{+}+\sqrt{(\mu_{-}+h_{z})^{2}-|\Delta_{0}|^{2}}$ for the outer ring, and $(k_{y}^{2}+k_{z}^{2})/(2m)=\mu_{+}-\sqrt{(\mu_{-}+h_{z})^{2}-|\Delta_{0}|^{2}}$ for the inner ring, when $(\mu_{-}+h_{z})^{2}-|\Delta_{0}|^{2}>0$; (b) doubly-degenerate line of nodes $(k_{y}^{2}+k_{z}^{2})/(2m)=\mu_{+}$ for $\mu_{+}>0$, doubly-degenerate point nodes for $\mu_{+}=0$, or no-line of nodes for $\mu_{+}<0$, when $(\mu_{-}+h_{z})^{2}-|\Delta_{0}|^{2}=0$; (c) no line of nodes when $(\mu_{-}+h_{z})^{2}-|\Delta_{0}|^{2}<0$. In addition, case (a) can be refined into cases (a2), (a1) and (a0). In case (a2), two rings indeed exist provided that $\mu_{+}>\sqrt{(\mu_{-}+h_{z})^{2}-|\Delta_{0}|^{2}}$. However, the inner ring disappears when $\mu_{+}=\sqrt{(\mu_{-}+h_{z})^{2}-|\Delta_{0}|^{2}}$. In case (a1), there is only one ring when $|\mu_{+}|<\sqrt{(\mu_{-}+h_{z})^{2}-|\Delta_{0}|^{2}},$ In case (a0), the outer ring disappears at $\mu_{+}=-\sqrt{(\mu_{-}+h_{z})^{2}-|\Delta_{0}|^{2}}$, and for $\mu_{+}<-\sqrt{(\mu_{-}+h_{z})^{2}-|\Delta_{0}|^{2}}$ no rings exist.

We choose our momentum, energy and velocity scales through the Fermi momentum $k_{F+}$ defined from the total density of fermions $n_{+}=n_{\uparrow}+n_{\downarrow}=k_{F_{+}}^{3}/(3\pi^{2}).$ This choice leads to the Fermi energy $\epsilon_{F_{+}}=k_{F_{+}}^{2}/2m$ and to the Fermi velocity $v_{F_{+}}=k_{F_{+}}/m$, as energy and velocity scales respectively. In Fig. 1, we show the phase diagram of Zeeman field $h_{z}/\epsilon_{F_{+}}$ versus chemical potential $\mu_{+}/\epsilon_{F_{+}}$ describing possible superfluid phases according to their quasiparticle excitation spectrum. We label the uniform superfluid phases with zero, one or two rings of nodes as US-0, US-1, and US-2, respectively. Non-uniform (NU) phases also emerge in regions where uniform phases are thermodynamically unstable. The US-2/US-1 phase boundary is determined by the condition $\mu_{+}=\sqrt{(\mu_{-}+h_{z})^{2}-|\Delta_{0}|^{2}}$, when $|\mu_{-}+h_{z}|>|\Delta_{0}|$; the US-0/US-2 boundary is determined by the Clogston-like condition $|(\mu_{-}+h_{z})|=|\Delta_{0}|$ when $\mu_{+}>0$, where the gapped US-0 phase disappears leading to the gapless US-2 phase; and the US-0/US-1 phase boundary is determined by $\mu_{+}=-\sqrt{(\mu_{-}+h_{z})^{2}-|\Delta_{0}|^{2}}$, when $|\mu_{-}+h_{z}|>|\Delta_{0}|$. Furthermore, with the US-0 boundaries, a crossover line between an indirectly gapped and a directly gapped US-0 phase occurs at $\mu_{+}=0$. Lastly, some important multi-critical points arise at the intersections of phase boundaries. First the point $\mu_{+}=0$ and $|(\mu_{-}+h_{z})|=|\Delta_{0}|$ corresponds to a tri-critical point for phases US-0, US-1, and US-2. Second, the point $|\Delta_{0}|=0$ and $\mu_{+}=|(\mu_{-}+h_{z})|$ corresponds to a tri-critical point for phases N, US-1 and US-2. In the limit where both $\mu_{-}$ and $h_{z}$ vanish no phase transitions take place and the problem is reduced to a crossover zhai-2011 ; hu-2011 ; han-2011 .

In the US-1 and US-2 phases near the zeros of $E_{2}({\bf k})$, quasiparticles have linear dispersion and behave as Dirac fermions. Such change in nodal structures is associated with bulk topological phase transitions of the Lifshitz class as noted for $p$-wave volovik-1992 and $d$-wave duncan-2000 ; botelho-2005 superfluids. Such Lifshitz topological phase transitions are possible here because the spin-orbit coupling field induces the triplet component of the order parameter $\Delta_{T}({\bf k})$. The loss of nodal regions correspond to annihilation of Dirac quasiparticles with opposite momenta, which lead to the disappearance of rings. The transition from phase US-2 to indirect gapped US-0 occurs through the merger of the two-rings at the phase boundary followed by the immediate opening of the indirect gap at finite momentum. However, the transition from phase US-2 to US-1 corresponds to the disappearance of the inner ring through the origin of momenta, similarly the transition from US-1 to the directly gappped US-0 corresponds to the disappearance of the last ring also through the origin of momenta. In the case of Rashba-only coupling rings of nodes are absent and it is possible to have at most nodal points  chuanwei-2011 ; iskin-2011 . The last two phase transitions are special because the zero-momentum quasiparticles at these phase boundaries correspond to true Majorana zero energy modes if the phase $\varphi({\bf k})$ of the spin-orbit field $h_{\perp}({\bf k})=|h_{\perp}({\bf k})|e^{i\varphi({\bf k})}$ and the phase $\theta({\bf k})$ of the order parameter $\Delta_{0}=|\Delta_{0}|e^{i\theta({\bf k})}$ have opposite phases at zero momentum: $\varphi({\bf 0})=-\theta({\bf 0})$ $[{\rm mod}(2\pi)]$. This can be seen from an analysis of the quasiparticle eigenfunction

corresponding to the eigenvalue $E_{2}({\bf k})$. The emergence of zero-energy Majorana fermions requires the quasiparticle to be its own anti-quasiparticle: $\Phi_{2}^{\dagger}({\bf k})=\Phi_{2}({\bf k})$. This can only happen at zero momentum ${\bf k}={\bf 0}$, where the amplitudes $u_{1}({\bf 0})=u_{3}^{*}({\bf 0})$ and $u_{2}({\bf 0})=u_{4}^{*}({\bf 0})$. Such requirement leads to the conditions $\mu_{+}^{2}=(\mu_{-}+h_{z})^{2}+|\Delta_{0}|^{2}$, and $\varphi({\bf 0})=-\theta({\bf 0})$ $[{\rm mod}(2\pi)]$, showing that Majorana fermions can exist only at the US-0/US-1 and US-2/US-1 phase boundaries. It is important to emphasize that the Majorana fermions found here exist in the bulk, and thus their emergence or disappeareance affect bulk thermodynamic properties, unlike Majorana fermions found at the edge (surfaces) of topological insulators and some topological superfluids. The common ground between bulk and surface Majorana fermions is that both exist at boundaries: the bulk Majorana zero-energy modes may exist at the phase boundaries between two topologically distinct superfluid phases, while surface Majorana zero-energy modes may exist at the spatial boundaries of a topologically non-trivial superfluid.

It is evident that the transition between different superfluid phases occurs without a change in symmetry in the order parameter $\Delta_{0}$, and thus violates the symmetry-based Landau classification of phase transitions. In the present case, the simultaneous existence of spin-orbit and Zeeman fields (internal or external) couple the singlet $\Delta_{S}({\bf k})$ and triplet $\Delta_{T}({\bf k})$ channels and all the superfluid phases US-0, US-1 and US-2 just have different weights from each order parameter component. However a finer classification based on topological charges can be made via the construction of topological invariants. Since the superfluid phases US-0, US-1, US-2 are characterized by different excitation spectra corresponding to the eigenvalues of the Hamiltonian matrix including interactions ${\bf H}({\bf k})$, we can use the resolvent matrix ${\bf R}(\omega,{\bf k})=\left[-\omega{\bf 1}+{\bf H}({\bf k})\right]^{-1}$ and the methods of algebraic topology nakahara-1990 to construct the topological invariant

where ${\bf\Lambda}_{k_{\mu}}={\bf R}\partial_{k_{\mu}}{\bf R}^{-1}.$ The topological invariant is $\ell=0$ in the gapped US-0 phase, is $\ell=1$ in the gapless US-1 phase and $\ell=2$ in the gapless US-2 phase, showing that, for ERD spin-orbit coupling, $\ell$ counts the number of rings of zero-energy excitations in each superfluid phase. The integral above has a hyper-surface measure $dS_{\gamma}$ and a domain ${\cal D}$ that encloses the region of zeros of $\omega=E_{j}({\bf k})=0$. Here $\mu,\nu,\lambda,\gamma$ run from 0 to 3, and $k_{\mu}$ has components $k_{0}=\omega$, $k_{1}=k_{x}$, $k_{2}=k_{y}$, and $k_{3}=k_{z}$. The topological invariant measures the flux of the four-dimensional vector $F^{\gamma}=\epsilon^{\mu\nu\lambda\gamma}{\rm Tr}\left[{\bf\Lambda}_{k_{\mu}}{\bf\Lambda}_{k_{\nu}}{\bf\Lambda}_{k_{\lambda}}\right]/{24\pi^{2}},$ through a hypercube including the singular region of the resolvent matrix ${\bf R}(\omega,{\bf k})$, much in the same way that the flux of the electric field ${\bf E}$ in Gauss’ law of classical electromagnetism measures the electric charge $q$ enclosed by a Gaussian surface: $\oint d{\bf S}\cdot(\epsilon_{0}{\bf E})=q$. Thus, the topological invariant defined above defines the topological charge of fermionic excitations, in the same sense as Gauss’ law for the electric flux defines the electric charge.

A full phase diagram can be constructed only upon verification of thermodynamic stability of all the proposed phases. For this purpose it becomes imperative to investigate the maximum entropy condition (see supplementary material). Independent of any microscopic approximations, the necessary and sufficient conditions for thermodynamic stability of a given phase are: positive isovolumetric heat capacity $C_{V}=T\left(\partial S/\partial T\right)_{V,\{N_{\alpha}\}}\geq 0;$ positive chemical susceptibility matrix $\xi_{\alpha\beta}=\left(\partial\mu_{\alpha}/\partial N_{\beta}\right)_{T,V},$ i.e, eigenvalues of the matrix $\left[\xi\right]$ are both positive; and positive bulk modulus $B=1/\kappa_{T}$ or isothermal compressibility $\kappa_{T}=-V^{-1}\left(\partial V/\partial P\right)_{T,\{N_{\alpha}\}}.$ Using these conditions, we construct the full phase diagram described in Fig. 1 for equal Rasha-Dresselhaus (ERD) spin-orbit coupling. The regions, where the uniform superfluid phases are unstable are labeled by the abbreviation NU to indicate that non-uniform phases such as phase separation or modulated superfluid (supersolid) may emerge. In Fig. 2, we show the phase diagram of Zeeman field $h_{z}/\epsilon_{F_{+}}$ versus interaction parameter $1/(k_{F+}a_{s})$, for population balanced fermions, where the number of spin-up fermions $N_{\uparrow}$ is equal to the number of spin-down fermions $N_{\downarrow}$.

Since these superfluid phases exhibit major changes in momentum-frequency space as evidenced by their single particle excitation spectrum, it is important to explore additional spectroscopic quantitities to characterize further the nature of these phases and the phase transitions between them. An important quantity is the $4\times 4$ resolvent matrix

from where the spectral density ${\cal A}_{\alpha}(\omega,{\bf k})=-(1/\pi){\rm Im}G_{\alpha\alpha}(i\omega=\omega+i\delta,{\bf k})$ for spin $\alpha=\uparrow,\downarrow$ can be extracted. The spectral function ${\cal A}_{\alpha}(\omega,{\bf k})$ in the plane of momenta $k_{y}$-$k_{z}$ with $k_{x}=0$ and frequency $\omega=0$ reveals the existence of rings of zero-energy excitations in the US-1 and US-2 phases. The density of states ${\cal D}_{\alpha}(\omega)=\sum_{\bf k}{\cal A}_{\alpha}(\omega,{\bf k})$ for spin $\alpha$ as a function of frequency $\omega$ is also an important spectroscopic quantity which is shown in Fig. 3 along with excitation spectra $E_{j}({\bf k})$ for phases US-1 and US-2 at fixed ERD spin-orbit coupling $v/v_{F_{+}}=0.28$. The parameters used for phase US-1 are $h_{z}/\epsilon_{F_{+}}=0.5$ and $1/(k_{F_{+}}a_{s})=-0.4$, while for phase US-2 they are $h_{z}/\epsilon_{F_{+}}=2.0$ and $1/(k_{F_{+}}a_{s})=1.0$. Notice that, even though the excitation spectrum $E_{j}({\bf k})$ is symmetric, the coherence factors appearing in the matrix ${\bf G}$ are not, such that the density of states ${\cal D}_{\alpha}(\omega)$ is not an even function of $\omega$, and thus it is not particle-hole symmetric. The main feature of ${\cal D}_{\alpha}(\omega)$ at low frequencies is the linear behavior due to the existence of Dirac quasiparticles and quasiholes in the US-1 and US-2 phases, which are absent in the direct-gap and the indirect-gap US-0 phases. The peaks and structures in ${\cal D}_{\alpha}(\omega)$ mostly emerge due to the maxima and minima of $E_{j}({\bf k})$. Notice that for finite Zeeman field $h_{z}$, the density of states ${\cal D}_{\uparrow}(\omega)\neq{\cal D}_{\downarrow}(\omega)$ because the induced population imbalance $P=(N_{\uparrow}-N_{\downarrow})/(N_{\uparrow}+N_{\downarrow})$ is non-zero. For the US-2 case shown in Fig. 3b, the induced population imbalance $P\ll 1$ since $h_{z}/\epsilon_{F_{+}}$ is small, while for the US-1 case shown in Fig. 3e, $P\approx 1$ as the spins are almost fully polarized since $h_{z}/\epsilon_{F_{+}}$ is large.

In summary, we have discussed the effects of spin-orbit and Zeeman fields in ultra-cold Fermi superfluids, obtained the phase diagrams of Zeeman field versus interaction parameter or versus chemical potential, and identified several bulk topological phase transitions between gapped and gapless superfluids as well as a variety of multi-critical points. We have shown that the presence of simultaneous Zeeman and spin-orbit fields induces higher angular momentum pairing, as manifested in the emergence of momentum dependence of the singlet and triplet components of order parameter expressed in the helicity basis. Finally, we have characterized topological phases and phase transitions between them through their excitation spectra (existence of Dirac quasiparticles or Majorana zero-energy modes), topological charges, and spectroscopic and thermodynamic properties, such as density of states and isothermal compressibility.