Superfluid weight and Berezinskii–Kosterlitz–Thouless temperature of spin-imbalanced and spin–orbit-coupled Fulde–Ferrell phases in lattice systems

Fulde–Ferrell–Larkin–Ovchinnikov (FFLO) superfluid states, identified by finite center-of-mass Cooper pairing momenta [1, 2], have gained widespread interest since their existence was predicted in the 1960s [3]. Traditionally, FFLO states are considered in the context of spin-imbalanced degenerate Fermi gases where finite momenta of condensed Cooper pairs originate from the mismatch between the Fermi surfaces of two pairing Fermion species [4, 5]. In such spin-polarized systems magnetism and superfluidity, usually thought to be incompatible with each other, co-exist and the superfluid order parameter is spatially varying, in contrast to the conventional Bardeen–Cooper–Schrieffer (BCS) pairing states characterized by the uniform order parameter and the absence of magnetism.

Realizing such spin-polarized FFLO states is challenging due to the requirement for large imbalance which in turn yields small superconducting order parameters and low critical temperatures. In recent years, a very different physical mechanism for realizing FFLO phases, namely the introduction of spin–orbit coupling (SOC) and Zeeman fields, has been investigated in many theoretical studies [6–27], for a review see [5]. The advantage of these SOC-induced FFLO states is the absence of large spin polarizations as now finite Cooper pairing momenta originate from the deformation of the single-particle band dispersions and not from the mismatch of Fermi surfaces. As large polarizations are not needed, SOC-induced FFLO states might have higher critical temperatures than conventional imbalance-induced FFLO phases.

Despite many theoretical studies supporting the existence of FFLO phases, direct observation of such exotic superfluid states has been lacking [3, 28]. For studying the FFLO state experimentally, ultracold Fermi gas systems are promising as they provide exact control of system parameters such as the spatial dimensionality, interaction strengths between the particles, and the system geometry [29–32]. Ultracold gas experiments performed with quasi-one-dimensional population-imbalanced atomic gases have shown to be consistent with the existence of the FFLO state [33] but unambiguous proof is still missing.

In addition to conventional spin-imbalanced quantum gas experiments, recently also synthetic SOC and Zeeman fields have been realized in ultracold gas experiments [34–38] which makes it possible to investigate SOC-induced FFLO states as well. As SOC-induced FFLO states have been predicted to be stable in larger parameter regime than conventional spin-imbalanced FFLO phases [10], synthetic SOC could provide a way to realize FFLO experimentally in ultracold gas systems [15].

Low dimensionality has been predicted to favor FFLO-pairing [39, 40]. However, in two and lower dimensional systems thermal phase fluctuations of the Cooper pair wave functions prevent the formation of true superfluid long range order as stated by the Mermin–Wagner theorem [41]. Instead, only quasi-long range order is possible. In two-dimensions, the phase transition from a normal Fermi gas to a superfluid state of quasi-long range order is determined by the Berezinskii–Kosterlitz–Thouless (BKT) transition temperature T_BKT [42]. Below T_BKT the system is a superfluid and above T_BKT superfluidity is lost.

In recent years, SOC-induced FFLO phases in two-dimensional systems have gained considerable attention [7, 10, 13, 14, 20, 21, 25]. In these systems it has been argued that SOC accompanied with the in-plane Zeeman field would yield FFLO states. Furthermore, in [13, 14] it was predicted that in the presence of the out-of-plane Zeeman field, i.e. spin-imbalance, SOC-induced FFLO states could be topologically non-trivial and support Majorana fermions. Such topological FFLO states are conceptually new and exotic superconductive phases of matter. However, these studies were performed by applying mean-field theories which do not consider the stability of FFLO states against thermal phase fluctuations in terms of the BKT transition. Superfluidity and BKT transition temperatures of BCS phases in spin–orbit-coupled Fermi gases have been theoretically investigated previously in [43–46] but BKT transitions of FFLO states have remained largely unstudied. As an exception, T_BKT for FFLO states in case of a 2D continuum system was explicitly computed in [12, 47–49] where it was shown that SOC is required in order to have a non-zero T_BKT for FFLO states. However, in case of spin–orbit-coupled lattice systems, T_BKT of FFLO phases has not been studied before. Lattice systems are interesting since, due to Fermi surface nesting effects, the FFLO states are expected to be more stable and accessible than in continuum [5, 40, 50].

FFLO-pairing states can be classified to two main categories: Fulde–Ferrell (FF) and Larkin–Ovchinnikov (LO) phases. In case of FF, the Cooper pair wave function Δ(r) is a plane wave associated with a single pairing momentum so that it has a uniform amplitude but a spatially oscillating complex phase. The LO wave function, on the contrary, consists of two plane waves of opposite momenta and therefore has spatially varying amplitude. In spin-imbalanced systems without SOC, it has been shown, at the mean-field level, that in a square lattice the LO states should be slightly more energetically favorable than FF states [51], whereas in the presence of SOC both FF and LO states can exist as was shown in [10]. Moreover, in [20, 26, 27] the existence of topologically non-trivial FFLO phases in square and triangular lattices was predicted. However, studies presented in [10, 20, 26, 27] did not consider the stability of FFLO phases against thermal phase fluctuations.

In this work we investigate the stability of FF phases in lattice systems with and without SOC by calculating the BKT transition temperature T_BKT. For a superconducting system the BKT temperature depends on the superfluid weight D^s which is responsible for the dissipationless electric current and the Meissner effect—the fundamental properties of superconductors [52, 53]. In our study we develop a general theory for obtaining D^s in any kind of lattice geometry in the presence of SOC and Zeeman fields, and apply the theory to a square lattice. We show that FF states in a square lattice indeed have a finite T_BKT with and without SOC, which is of fundamental importance as well as a prerequisite for their experimental observation. Topological FF states created by the interplay of SOC and Zeeman fields are identified with the Chern numbers C = {±1, −2}, and we explain how different topological FF phases can be distinguished by investigating the momentum density profiles which are experimentally accessible quantities. Additionally, we compare the superfluid weight components in orthogonal spatial directions. We also compute the so-called geometric superfluid weight component which is just recently found new superfluid contribution that depends on the geometric properties of the single-particle Bloch functions [54, 55].

In our study we discard the existence of LO phases as the LO ansatzes break the translational invariance which is required for deriving the superfluid weight in a simple form. Ignoring LO states, however, is not an issue because we are interested in the stability and BKT transition temperatures of exotic superfluid states: if there exists more stable LO states than FF states that we find, it implies the BKT transition temperatures of these LO states being higher than the temperatures we obtain for FF states. Therefore, our results can be considered as conservative estimates. Furthermore, in [10, 26] LO states were argued to exist when the superfluid pairing occurs with in both helicity branches of a spin–orbit-coupled square lattice. Thus, by studying the pairing amplitude profiles, we can deduce in which parts of our parameter space LO states would be more stable than the FF states we study.

The rest of the article is structured as follows. In the next section we provide expressions for the superfluid weight and thus for T_BKT in the presence of SOC in case of an arbitrary lattice geometry. In section 3 we apply our equations for a spin–orbit-coupled square lattice and show T_BKT for various system parameters. We also discuss the topological properties of the system, and the different components of the superfluid weight. Lastly, in section 5 we present concluding remarks and an outlook for future research.

In this section we derive the expressions for the superfluid weight in the framework of BCS mean-field theory by applying linear response theory in a very similar way as was done in [55]. We consider the following two-dimensional Fermi–Hubbard Hamiltonian

$$\begin{eqnarray}&&H=\displaystyle \sum _{i,j,\alpha ,\beta ,\sigma ,\sigma ^{\prime} }{t}_{i\alpha \sigma ,j\beta \sigma ^{\prime} }{c}_{i\alpha \sigma }^{\dagger }{c}_{j\beta \sigma ^{\prime} }-\displaystyle \sum _{i\alpha \sigma }{\mu }_{\sigma }{c}_{i\alpha \sigma }^{\dagger }{c}_{i\alpha \sigma }+U\displaystyle \sum _{i\alpha }{c}_{i\alpha \uparrow }^{\dagger }{c}_{i\alpha \uparrow }{c}_{i\alpha \downarrow }^{\dagger }{c}_{i\alpha \downarrow },\end{eqnarray} \tag{ 1 }$$

where ${c}_{i\alpha \sigma }^{\dagger }$ creates a fermion in the α-orbital of the ith unit cell with spin $\sigma \in \{\uparrow ,\downarrow \}$. The first term describes the hopping processes which in addition to usual kinetic hopping terms ($\sigma =\sigma ^{\prime}$) can now also include spin-flipping terms ($\sigma \ne \sigma ^{\prime}$) required to take into account the SOC contribution. In the second term μ_σ is the spin-dependent chemical potential and the last term is the attractive on-site Hubbard interaction characterized by the coupling strength U < 0. The above Hamiltonian describes any two-dimensional lattice geometry with arbitrary hopping and spin-flip terms, including the Rashba spin–orbit-coupled two-component Fermi gases considered in this work.

We treat the interaction term by performing the standard mean-field approximation ${{Uc}}_{i\alpha \uparrow }^{\dagger }{c}_{i\alpha \uparrow }{c}_{i\alpha \downarrow }^{\dagger }{c}_{i\alpha \downarrow }\approx {{\rm{\Delta }}}_{i\alpha }{c}_{i\alpha \downarrow }{c}_{i\alpha \uparrow }+{{\rm{\Delta }}}_{i\alpha }^{\dagger }{c}_{i\alpha \uparrow }^{\dagger }{c}_{i\alpha \downarrow }^{\dagger }$ where ${{\rm{\Delta }}}_{i\alpha }=U\langle {c}_{i\alpha \downarrow }{c}_{i\alpha \uparrow }\rangle$ is the superfluid order parameter or in other words the wave function of the condensed Cooper pairs. To investigate the properties of the usual BCS and exotic inhomogeneous FF superfluid phases, we let the order parameter to have the form ${{\rm{\Delta }}}_{i\alpha }={{\rm{\Delta }}}_{\alpha }\exp [{\rm{i}}\tilde{{\bf{q}}}\cdot {{\bf{r}}}_{i}]$, where $\tilde{{\bf{q}}}$ is the Cooper pair momentum and r_i is the spatial coordinate of the ith unit cell. The momentum of the Cooper pairs in a FF phase is finite, in contrast to a normal BCS phase where the Cooper pairs do not carry momentum.

By performing the Fourier transform to the momentum space ${c}_{i\alpha \sigma }=(1/\sqrt{N}){\sum }_{{\bf{k}}}{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {{\bf{r}}}_{i}}{c}_{\sigma {\bf{k}}\alpha }$, where N is the number of unit cells, one can rewrite the Hamiltonian in the form (discarding the constant terms)

$$\begin{eqnarray}H & = & \displaystyle \sum _{{\bf{k}}}\left(\left[\begin{array}{cc}{c}_{\uparrow {\bf{k}}}^{\dagger } & {c}_{\downarrow {\bf{k}}}^{\dagger }\end{array}\right]\right.\left[\begin{array}{cc}{{ \mathcal H }}_{\uparrow }({\bf{k}})-{\mu }_{\uparrow } & {\rm{\Lambda }}({\bf{k}})\\ {{\rm{\Lambda }}}^{\dagger }({\bf{k}}) & {{ \mathcal H }}_{\downarrow }({\bf{k}})-{\mu }_{\downarrow }\end{array}\right]\left[\begin{array}{c}{c}_{\uparrow {\bf{k}}}\\ {c}_{\downarrow {\bf{k}}}\end{array}\right]\\ & & +{c}_{\uparrow {\bf{k}}}^{\dagger }{\rm{\Delta }}{c}_{\downarrow \tilde{{\bf{q}}}-{\bf{k}}}^{\dagger }+{c}_{\downarrow \tilde{{\bf{q}}}-{\bf{k}}}{{\rm{\Delta }}}^{\dagger }{c}_{\uparrow {\bf{k}}}),\end{eqnarray} \tag{ 2 }$$

where ${c}_{\sigma {\bf{k}}}^{\dagger }=[{c}_{\sigma {\bf{k}}1},{c}_{\sigma {\bf{k}}2},\,\ldots ,\,{c}_{\sigma {\bf{k}}M}]$ and ${\rm{\Delta }}={\rm{diag}}({{\rm{\Delta }}}_{1},{{\rm{\Delta }}}_{2},\,\ldots ,\,{{\rm{\Delta }}}_{M})$, M being the number of orbitals within a unit cell. Furthermore, ${{ \mathcal H }}_{\sigma }({\bf{k}})$ and ${\rm{\Lambda }}({\bf{k}})$ are the Fourier transforms of the kinetic hopping and the spin-flip terms, respectively.

To write our Hamiltonian in a more compact form, let us introduce a four-component spinor ${{\rm{\Psi }}}_{{\bf{k}}}$ and rewrite the Hamiltonian as follows:

$$\begin{eqnarray}&&H=\displaystyle \frac{1}{2}\displaystyle \sum _{{\bf{k}}}{{\rm{\Psi }}}_{{\bf{k}}}^{\dagger }{{ \mathcal H }}_{{\bf{k}}}{{\rm{\Psi }}}_{{\bf{k}}},\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{{\bf{k}}}=\left[\begin{array}{c}{c}_{\uparrow {\bf{k}}}\\ {c}_{\downarrow {\bf{k}}}\\ {c}_{\downarrow \tilde{{\bf{q}}}-{\bf{k}}}^{\dagger }\\ -{c}_{\uparrow \tilde{{\bf{q}}}-{\bf{k}}}^{\dagger }\end{array}\right]\equiv \left[\begin{array}{c}{\psi }_{{\bf{k}}}\\ {\rm{i}}{\tau }^{y}{({\psi }_{\tilde{{\bf{q}}}-{\bf{k}}}^{\dagger })}^{T}\end{array}\right]\equiv \left[\begin{array}{c}{\psi }_{{\bf{k}}}\\ {\psi }_{2,{\bf{k}}}\end{array}\right],\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}{{ \mathcal H }}_{{\bf{k}}}=\left[\begin{array}{cc}{{ \mathcal H }}_{p}({\bf{k}})-\tilde{\mu } & \tilde{{\rm{\Delta }}}\\ {\tilde{{\rm{\Delta }}}}^{\dagger } & -{{ \mathcal H }}_{h}({\bf{k}}-\tilde{{\bf{q}}})+{\tau }^{y}\tilde{\mu }{\tau }^{y}\end{array}\right],\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}{{ \mathcal H }}_{p}({\bf{k}})=\left[\begin{array}{cc}{{ \mathcal H }}_{\uparrow }({\bf{k}}) & {\rm{\Lambda }}({\bf{k}})\\ {{\rm{\Lambda }}}^{\dagger }({\bf{k}}) & {{ \mathcal H }}_{\downarrow }({\bf{k}})\end{array}\right],\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{{ \mathcal H }}_{h}({\bf{k}})=-{\rm{i}}{\tau }^{y}{{ \mathcal H }}_{p}^{* }(-{\bf{k}}){\rm{i}}{\tau }^{y},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}\tilde{{\rm{\Delta }}}=\left[\begin{array}{cc}{\rm{\Delta }} & 0\\ 0 & {\rm{\Delta }}\end{array}\right],\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}\tilde{\mu }=\left[\begin{array}{cc}{\mu }_{\uparrow }{I}_{M} & 0\\ 0 & {\mu }_{\downarrow }{I}_{M}\end{array}\right].\end{eqnarray} \tag{ 9 }$$

Here ${\tau }^{y}={\hat{\sigma }}_{y}\otimes {I}_{M}$, where I_M is a M × M identity matrix and $\hat{\sigma }=[{\hat{\sigma }}_{x},{\hat{\sigma }}_{y},{\hat{\sigma }}_{z}]$ are the Pauli matrices. One should note that now the single-particle Hamiltonian is not anymore simply ${{ \mathcal H }}_{\sigma }$ but ${{ \mathcal H }}_{p}$ in which the two spin components are coupled via ${\rm{\Lambda }}({\bf{k}})$.

In two-dimensions the total superfluid weight D^s is a 2 × 2 tensor which reads

$$\begin{eqnarray}{D}^{s}=\left[\begin{array}{cc}{D}_{{xx}}^{s} & {D}_{{xy}}^{s}\\ {D}_{{yx}}^{s} & {D}_{{yy}}^{s}\end{array}\right],\end{eqnarray} \tag{ 10 }$$

where x and y are the spatial dimensions. To compute the superfluid weight tensor elements ${D}_{\mu \nu }^{s}$, we exploit the fact that at the mean-field level ${D}_{\mu \nu }^{s}$ is the long-wavelength, zero-frequency limit of the current–current response function ${K}_{\mu \nu }$ [53], that is

$$\begin{eqnarray}{D}_{\mu \nu }^{s} & = & \mathop{\mathrm{lim}}\limits_{{\bf{q}}\to 0}\mathop{\mathrm{lim}}\limits_{\omega \to 0}{K}_{\mu \nu }({\bf{q}},\omega )\\ & = & \mathop{\mathrm{lim}}\limits_{{\bf{q}}\to 0}\mathop{\mathrm{lim}}\limits_{\omega \to 0}\left[\left\langle {T}_{\mu \nu }\rangle -{\rm{i}}{\displaystyle \int }_{0}^{\infty }{{\rm{d}}{t}{\rm{e}}}^{{\rm{i}}\omega t}\langle [{j}_{\mu }^{p}({\bf{q}},t),{j}_{\nu }^{p}(-{\bf{q}},0)]\right\rangle \right],\end{eqnarray} \tag{ 11 }$$

where ${j}^{p}({\bf{q}})$ and T are the paramagnetic and diamagnetic current operators, respectively. The current operators can be derived by applying the Peierls substitution to the single-particle Hamiltonian ${{ \mathcal H }}_{p}$ such that the hopping elements, both kinetic and spin-flipping terms, are modified by a phase factor of $\exp [-{\rm{i}}{\bf{A}}\cdot ({{\bf{r}}}_{j}-{{\bf{r}}}_{i})]$ where ${\bf{A}}$ is the vector potential. By assuming the phase factor to be spatially slowly varying, we can expand the Hamiltonian up to second order in A to obtain $H={j}_{\mu }^{p}{A}_{\mu }+{T}_{\mu \nu }{A}_{\mu }{A}_{\nu }/2$. In our case the μ-component of the paramagnetic and diamagnetic current operators can be cast in the form

$$\begin{eqnarray}{j}_{\mu }^{p}({\bf{q}}) & = & \displaystyle \sum _{{\bf{k}}}{\psi }_{{\bf{k}}+{\bf{q}}}^{\dagger }{\partial }_{\mu }{{ \mathcal H }}_{p}({\bf{k}}+{\bf{q}}/2){\psi }_{{\bf{k}}}\\ & = & \displaystyle \sum _{{\bf{k}}}{{\rm{\Psi }}}_{{\bf{k}}+{\bf{q}}}^{\dagger }{\partial }_{\mu }{ \mathcal H }({\bf{k}}+{\bf{q}}/2){P}_{+}{{\rm{\Psi }}}_{{\bf{k}}}\end{eqnarray} \tag{ 12 }$$

and

$$\begin{eqnarray}{T}_{\mu \nu }({\bf{q}}) & = & \displaystyle \sum _{{\bf{k}}}{\psi }_{{\bf{k}}}^{\dagger }{\partial }_{\mu }{\partial }_{\nu }{{ \mathcal H }}_{p}({\bf{k}}){\psi }_{{\bf{k}}}\\ & = & \displaystyle \sum _{{\bf{k}}}{{\rm{\Psi }}}_{{\bf{k}}}^{\dagger }{\partial }_{\mu }{\partial }_{\nu }{ \mathcal H }({\bf{k}}){P}_{+}{{\rm{\Psi }}}_{{\bf{k}}},\end{eqnarray} \tag{ 13 }$$

where ${P}_{+}=({I}_{4M}+{\hat{\sigma }}^{z}\otimes {I}_{2M})/2$ and more generally ${P}_{\pm }=({I}_{4M}\pm {\hat{\sigma }}^{z}\otimes {I}_{2M})/2$.

We are interested in computing the current–current response function ${K}_{\mu \nu }({\bf{q}},\omega )$ which at the limit of ${\bf{q}}\to 0$, ω = 0 yields the superfluid weight ${D}_{\mu \nu }^{s}$. To this end, we first define a Green's function $G(\tau ,{\bf{k}})=-\langle T{{\rm{\Psi }}}_{{\bf{k}}}(\tau ){{\rm{\Psi }}}_{{\bf{k}}}^{\dagger }(0)\rangle$. In the Matsubara frequency space this reads $G({\rm{i}}{\omega }_{n},{\bf{k}})=1/({\rm{i}}{\omega }_{n}-{ \mathcal H }({\bf{k}}))$ which follows from the quadratic form of the Hamiltonian (3). Now, the current operators (12), (13), the Green's function and the Hamiltonian all have the same structure as those for conventional BCS theory developed in [55]. Thus one can compute, by applying the Matsubara formalism and analytic continuation, the current–current response function in a similar fashion as done in [55]. One starts from (11), inserts the expressions (12), (13) for the current operators, deploys the Matsubara formalism, applies the diagrammatic expansion up to first order diagrams and obtains

$$\begin{eqnarray}{K}_{\mu \nu }({\bf{q}},{\rm{i}}{\omega }_{n}) & = & \displaystyle \frac{1}{\beta }\displaystyle \sum _{{\bf{k}}}\displaystyle \sum _{{{\rm{\Omega }}}_{m}}{\rm{Tr}}[{\partial }_{\mu }{\partial }_{\nu }{ \mathcal H }({\bf{k}}){P}_{+}G({\rm{i}}{{\rm{\Omega }}}_{m},{\bf{k}})\\ & & +{\partial }_{\mu }{ \mathcal H }({\bf{k}}+{\bf{q}}/2){P}_{+}G({\rm{i}}{\omega }_{n}+{\rm{i}}{{\rm{\Omega }}}_{m},{\bf{k}}+{\bf{q}})\\ & & \times {\partial }_{\nu }{ \mathcal H }({\bf{k}}+{\bf{q}}/2){\hat{\gamma }}_{z}G({\rm{i}}{{\rm{\Omega }}}_{m},{\bf{k}})],\end{eqnarray} \tag{ 14 }$$

where $\beta =1/{k}_{B}T$, ${\hat{\gamma }}_{z}={\hat{\sigma }}_{z}\otimes {I}_{2M}$, and ω_n (Ω_m) are bosonic (fermionic) Matsubara frequencies. From (14) one eventually obtains (see appendix A):

$$\begin{eqnarray}{D}_{\mu \nu }^{s} & = & {K}_{\mu \nu }({\bf{q}}\to 0,0)\\ & = & 2\displaystyle \sum _{{\bf{k}},i,j}\displaystyle \frac{n({E}_{j,{\bf{k}}})-n({E}_{i,{\bf{k}}})}{{E}_{i,{\bf{k}}}-{E}_{j,{\bf{k}}}}(\langle {\phi }_{i}({\bf{k}})| {\partial }_{\mu }{ \mathcal H }({\bf{k}}){P}_{+}| {\phi }_{j}({\bf{k}})\rangle \\ & & \times \langle {\phi }_{j}({\bf{k}})| {P}_{-}{\partial }_{\nu }{ \mathcal H }({\bf{k}})| {\phi }_{i}({\bf{k}})\rangle ),\end{eqnarray} \tag{ 15 }$$

where n(E) is the Fermi–Dirac distribution and $| {\phi }_{i}({\bf{k}})\rangle$ are the eigenvectors of ${ \mathcal H }({\bf{k}})$ with the eigenvalues ${E}_{i,{\bf{k}}}$. For i = j, the prefactor should be understood as $-{\partial }_{{E}_{i}}n({E}_{i})$, which vanishes at zero temperature if the quasi-particle spectrum is gapped. For gapless excitations, $-{\partial }_{{E}_{i}}n({E}_{i})$ gives finite contribution even at zero temperature. We have benchmarked our superfluid weight relation (15) to earlier studies as discussed in appendix C.

The BKT transition temperature T_BKT can be obtained from the superfluid weight tensor by using the generalized KT-Nelson criterion [56] for the anisotropic superfluid [12, 49]:

$$\begin{eqnarray}&&{T}_{{\rm{BKT}}}=\displaystyle \frac{\pi }{8}\sqrt{\det [{D}^{s}({T}_{{\rm{BKT}}})]}.\end{eqnarray} \tag{ 16 }$$

In the computations presented in this work D^s is at low temperatures nearly a constant and therefore we can safely use the following approximation

$$\begin{eqnarray}&&{T}_{{\rm{BKT}}}\approx \displaystyle \frac{\pi }{8}\sqrt{\det [{D}^{s}(T=0)]}.\end{eqnarray} \tag{ 17 }$$

In [54, 55] it was shown that in case of conventional BCS states the superfluid weight can be divided to two parts: the so-called conventional and geometric contributions, ${D}_{\mu \nu }^{s}={D}_{{\rm{conv}},\mu \nu }^{s}+{D}_{{\rm{geom}},\mu \nu }^{s}$. The conventional superfluid term ${D}_{{\rm{conv}},\mu \nu }^{s}$ depends only on the single-particle energy dispersion relations, whereas the geometric part ${D}_{{\rm{geom}},\mu \nu }^{s}$ comprises the geometric properties of the Bloch functions. In a similar fashion than in [55], also in our case the superfluid weight can be split to conventional and geometric parts so that ${D}_{{\rm{conv}},\mu \nu }^{s}$ is a function of the single-particle dispersions of ${{ \mathcal H }}_{p}$ and ${{ \mathcal H }}_{h}$, and correspondingly ${D}_{{\rm{geom}},\mu \nu }^{s}$ depends on the Bloch functions of ${{ \mathcal H }}_{p}$ and ${{ \mathcal H }}_{h}$. The separation of D^s to ${D}_{{\rm{geom}}}^{s}$ and ${D}_{{\rm{conv}}}^{s}$ terms is shown in appendix B.

The above expression (15) of the superfluid weight holds for an arbitrary multiband lattice system. Here we focus on the simplest possible case, namely the square lattice geometry where the so-called Rashba SOC is applied to induce FF phases. By computing the superfluid weight and thus the BKT transition temperature, one can investigate the stability of SOC-induced FF phases versus the conventional FF phases induced by the spin-imbalance. We start by writing the Hamiltonian in the form

$$\begin{eqnarray}H & = & -t\displaystyle \sum _{\langle i,j\rangle ,\sigma }{c}_{i\sigma }^{\dagger }{c}_{j\sigma }-\mu \displaystyle \sum _{i\sigma }{c}_{i\sigma }^{\dagger }{c}_{i\sigma }+U\displaystyle \sum _{i}{c}_{i\uparrow }^{\dagger }{c}_{i\uparrow }{c}_{i\downarrow }^{\dagger }{c}_{i\downarrow }\\ & & +{H}_{z,{\rm{in}}}+{H}_{z,{\rm{out}}}+{H}_{{\rm{SOC}}},\end{eqnarray} \tag{ 18 }$$

where the first term is the usual nearest-neighbor hopping term (we discard the orbital indices as in a square lattice there is only one lattice site per unit cell). The last three terms are the in-plane Zeeman field, out-of-plane Zeeman field and the Rashba-coupling, respectively. They are

$$\begin{eqnarray}&&{H}_{z,{\rm{in}}}={h}_{x}\displaystyle \sum _{i}{c}_{i}^{\dagger }{\hat{\sigma }}_{x}{c}_{i},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{H}_{z,{\rm{out}}}={h}_{z}\displaystyle \sum _{i}{c}_{i}^{\dagger }{\hat{\sigma }}_{z}{c}_{i},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{H}_{{\rm{SOC}}}={\rm{i}}\lambda \displaystyle \sum _{\langle i,j\rangle }{c}_{i}^{\dagger }{({{\bf{d}}}_{{ij}}\times \hat{\sigma })}_{z}{c}_{j}.\end{eqnarray} \tag{ 21 }$$

Here ${{\bf{d}}}_{{ij}}$ is the unit vector connecting the nearest-neighbor sites i and j, $\hat{\sigma }={[{\hat{\sigma }}_{x},{\hat{\sigma }}_{y},{\hat{\sigma }}_{z}]}^{T}$ are the Pauli matrices and ${c}_{i}={[{c}_{i\uparrow },{c}_{i\downarrow }]}^{T}$. The out-of-plane Zeeman fields can be included to the spin-dependent chemical potentials by writing ${\mu }_{\uparrow }=\mu +{h}_{z}$ and ${\mu }_{\downarrow }=\mu -{h}_{z}$. Furthermore, due to the in-plane Zeeman field and the Rashba spin-flipping terms, ${\rm{\Lambda }}({\bf{k}})$ in (2) has the form ${\rm{\Lambda }}({\bf{k}})={h}_{x}-2\lambda (\sin {k}_{y}+{\rm{i}}\sin {k}_{x})$. We determine the order parameter amplitude Δ and the Cooper pair momentum $\tilde{{\bf{q}}}$ self-consistently by minimizing the grand canonical thermodynamic potential ${\rm{\Omega }}({\rm{\Delta }},\tilde{{\bf{q}}})=-{k}_{B}T\mathrm{log}[{\rm{Tr}}({{\rm{e}}}^{-\beta H})]$ which in the mean-field framework at T = 0 reads as

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{M.F.}}}=-\displaystyle \frac{{{\rm{\Delta }}}^{2}}{U}+\displaystyle \frac{1}{2}\displaystyle \sum _{{\bf{k}},\nu ,\eta }{E}_{{\bf{k}},\nu }^{\eta }{\rm{\Theta }}(-{E}_{{\bf{k}},\nu }^{\eta }),\end{eqnarray} \tag{ 22 }$$

where ${\rm{\Theta }}(x)$ is the Heaviside step function and ${E}_{{\bf{k}},\nu }^{\eta }$ are the eigenvalues of ${{ \mathcal H }}_{{\bf{k}}}$. Here $\eta =\{+,\,-\}$ labels the quasi-particle and quasi-hole branches, respectively and ν = {1, 2} the helicity branches split by the SOC. The quasi-particle branches are taken to be the two highest eigenvalues of ${{ \mathcal H }}_{{\bf{k}}}$. In (22) we have discarded the constant term ${\sum }_{{\bf{k}}}{\rm{Tr}}[{{ \mathcal H }}_{h}({\bf{k}}-\tilde{{\bf{q}}})-{\tau }^{y}\tilde{\mu }{\tau }^{y}]$ which is not needed when one minimizes ${{\rm{\Omega }}}_{{\rm{M.F.}}}$. Consistent with previous lattice studies [10, 26, 27], the Cooper pair momentum is in the y-direction, i.e. $\tilde{{\bf{q}}}={\tilde{q}}_{y}{\hat{{\bf{e}}}}_{y}$ as the in-plane Zeeman field in the x-direction deforms the single-particle dispersions in the y-direction. We have numerically checked that the solutions with the Cooper pair momentum in the y-direction minimize the thermodynamic potential, as discussed in appendix E. When the correct values for Δ and ${\tilde{q}}_{y}$ are found, the superfluid weight can be computed with (15).

We investigate the topological properties by computing the Chern number C for our interacting system by integrating the Berry curvature ${{\rm{\Gamma }}}_{\nu }^{\eta }({\bf{k}})$ associated with the quasi-hole branches η = − over the first Brillouin zone as follows:

$$\begin{eqnarray}&&C=\displaystyle \frac{1}{2\pi }\displaystyle \sum _{\nu =1}^{2}{\int }_{-\pi }^{\pi }{\int }_{-\pi }^{\pi }{{\rm{d}}{k}}_{x}{{\rm{d}}{k}}_{y}{{\rm{\Gamma }}}_{\nu }^{-}({\bf{k}}).\end{eqnarray} \tag{ 23 }$$

The explicit form for the Berry curvature can be expressed with the eigenvalues ${E}_{{\bf{k}},\nu }^{\eta }$ of ${{ \mathcal H }}_{{\bf{k}}}$ and the corresponding eigenvectors $| n({\bf{k}})\rangle$, where $n=(\eta ,\nu )$, in the form

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\nu }^{\eta }({\bf{k}})={\rm{i}}\displaystyle \sum _{n\ne n^{\prime} }\displaystyle \frac{\langle n| {\partial }_{{k}_{x}}{{ \mathcal H }}_{{\bf{k}}}| n^{\prime} \rangle \langle n^{\prime} | {\partial }_{{k}_{y}}{{ \mathcal H }}_{{\bf{k}}}| n\rangle -({k}_{x}\leftrightarrow {k}_{y})}{{({E}_{{\bf{k}},\nu }^{\eta }-{E}_{{\bf{k}},\nu ^{\prime} }^{\eta ^{\prime} })}^{2}}.\end{eqnarray} \tag{ 24 }$$

4.1. Phase diagrams and the BKT temperature

By deploying our mean-field formalism we determine the phase diagrams and T_BKT as functions of the Zeeman fields and the average chemical potential $\mu =({\mu }_{\uparrow }+{\mu }_{\downarrow })/2$. We fix the temperature to T = 0 as, according to (17), the zero temperature superfluid weight gives a good estimate for T_BKT. In all the computations we choose t = 1 and U = −4. Furthermore, we let ${\tilde{q}}_{y}$ to have only discrete values in the first Brillouin zone such that ${\tilde{q}}_{y}\in \{\tfrac{\pi n}{L},n=1,2,\,\ldots L\}$, where L is the length of the lattice in one direction, i.e. the total number of lattice sites is N = L × L. In all of our computations we choose L = 104 and deploy periodic boundary conditions.

In figures 1(a), (b) the superfluid phase diagrams in terms of the magnitude of ${\tilde{q}}_{y}$ are presented as a function of h_x and h_z at $\mu =0.95$ for λ = 0 and λ = 0.75, respectively, and the corresponding BKT transition temperatures T_BKT are shown in figures 1(c), (d). From figure 1(a) we see that in the absence of SOC the phase diagram is symmetric with respect to the Zeeman field orientation. This is due to the SO(2) symmetry as under the rotation ${ \mathcal U }{[{c}_{i\uparrow },{c}_{i\downarrow }]}^{T}{{ \mathcal U }}^{-1}=\tfrac{1}{\sqrt{2}}{[{c}_{i\uparrow }+{c}_{i\downarrow },{c}_{i\uparrow }-{c}_{i\downarrow }]}^{T}\equiv {[{d}_{i\uparrow },{d}_{i\downarrow }]}^{T}$ the Hamiltonian remains invariant except ${h}_{x}\to {h}_{z}$ and ${h}_{z}\to {h}_{x}$. For small Zeeman fields, the BCS phase is the ground state and becomes only unstable against the FF phase for larger Zeeman field strengths. One can see from figure 1(c) that the BKT temperature for the BCS phase is ${T}_{{\rm{BKT}}}\approx 0.25t$ and roughly ${T}_{{\rm{BKT}}}\approx 0.1t$ for the FF phase. This implies that conventional imbalance-induced FF phases without SOC could be observed in lattice systems, in contrast to continuum systems where it is shown that T_BKT = 0 [47]. This is the first time that the stability against the thermal phase fluctuations of spin-imbalanced FF states in a lattice system is confirmed.

Zoom In Zoom Out Reset image size

Unlike in the case of without SOC, the phase diagram shown in figure 1(b) for λ = 0.75 depends on the direction of the total Zeeman field, as SOC together with the in-plane Zeeman field breaks the SO(2) symmetry. The interplay of the SOC and the Zeeman fields stabilize inhomogeneous superfluidity in larger parameter regions than in case of conventional spin-imbalanced FF states. Furthermore, by introducing SOC one is able realize topologically distinct BCS and FF phases. As with λ = 0, at small Zeeman fields there exist topologically trivial BCS states. When h_x is increased, the system enters non-topological FF phase and eventually for large enough h_x topological FF states of C = −1 (tFF₋₁) and C = −2 (tFF₋₂). By applying large h_z one is able to reach topological BCS and FF phases, tBCS₋₂ and tFF₋₂, characterized by C = −2. For large enough Zeeman fields the superfluidity is lost and the system enters normal (N) state.

From figure 1(b) we see that in addition to topological classification, FF phases can be further distinguished by the magnitude of the Cooper pair momentum ${\tilde{q}}_{y}$: for intermediate Zeeman field strengths the FF state is characterized by rather small ${\tilde{q}}_{y}$, in contrast to region of large Zeeman fields where the pairing momenta are comparable to those of FF states of λ = 0. The same behavior can be seen by observing T_BKT presented in figure 1(d). We see that for small-${\tilde{q}}_{y}$ region T_BKT is around $0.3t$ and becomes only smaller for large-${\tilde{q}}_{y}$ region where T_BKT at largest is roughly ${T}_{{\rm{BKT}}}\approx 0.17t$. Therefore, by deploying SOC, one is able to stabilize FF phases considerably against thermal phase fluctuations and increase T_BKT. This is similar to continuum studies [12, 48, 49] where it was proposed that FF states could be observed with the aid of SOC. The difference of λ = 0 and λ = 0.75 is further demonstrated in figure 1(e), where T_BKT and ${\tilde{q}}_{y}$ for both the cases are plotted as a function of h_x at h_z = 0. We see that the phase diagram becomes richer and T_BKT is increased when SOC is deployed.

To understand why in the presence of SOC there exist distinct FF regions of considerably different BKT temperatures, we investigate the inter- and intraband pairing functions $\langle {c}_{{\bf{k}},n}{c}_{\tilde{{\bf{q}}}-{\bf{k}},n^{\prime} }\rangle$, where ${c}_{{\bf{k}},n}$ is the annihilation operator for the nth Bloch function of the single-particle Hamiltonian ${{ \mathcal H }}_{p}({\bf{k}})$. In case of a square lattice, ${{ \mathcal H }}_{p}({\bf{k}})$ is a 2 × 2 matrix so we have two energy bands, called also helicity branches. As an example, in figure 2 the single-particle energy dispersion bands have been plotted at h_z = 0 for λ = 0, h_x = 0 (figures 2(a), (b)), $\lambda \ne 0$, h_x = 0 (figures 2(c), (d)) and $\lambda \ne 0$, ${h}_{x}\ne 0$ (figures 2(e), (f)). Without SOC, the single-particle dispersions for spin-up and down components are degenerate (figures 2(a), (b)). By turning on the SOC, this degeneracy is lifted (figures 2(c), (d)) and when also h_x is applied, the dispersion becomes deformed in a non-symmetric way with respect to k_y = 0 (figures 2(e), (f)). This deformation of the dispersions results in the intraband pairing of finite momentum in the y-direction when h_x is large enough as there exists a momentum mismatch of ${\tilde{q}}_{y}{\hat{{\bf{e}}}}_{y}$ between the pairing fermions. If in addition the interband pairing occurs, the momentum mismatch can exist also in the x-direction and consequently the Cooper pair momentum is not necessarily in the y-direction. However, in the computations presented in this work $\tilde{{\bf{q}}}$ has been numerically checked to be always in the y-direction.

Zoom In Zoom Out Reset image size

With figures 2(e), (f) one can also understand the fundamental differences between conventional spin-imbalanced-induced and SOC-induced FF states in terms of spontaneously broken symmetries. Both cases break the time-reversal symmetry spontaneously and in case of spin-imbalanced FF also the rotational symmetry within the lattice plane is spontaneously broken. In other words, for imbalance-induced FF states, it is energetically equally favorable for the Cooper pair momentum to be in the x- or y-direction. However, SOC and the in-plane Zeeman field break the rotational symmetry explicitly, and therefore the Cooper pair wavevector is forced to be in the perpendicular direction with respect to the in-plane Zeeman field as the dispersions are deformed in that direction (figures 2(e), (f)).

Even if the in-plane Zeeman field causes the single-particle dispersion to be non-centrosymmetric, it is still not a sufficient condition to reach the FF state as can be seen in figure 1(b) where the ground state is BCS for small enough values of h_x. Homogeneous BCS states can be still more favorable than FF states if for example the chemical potential is such that the shapes and the density of states of the Fermi surfaces prefer the Cooper pairing with zero momentum. However, when the in-plane Zeeman field becomes strong enough, the deformation of the dispersion results in the FF pairing.

In figures 3(a)–(i) we present $| \langle {c}_{{\bf{k}},1}{c}_{\tilde{{\bf{q}}}-{\bf{k}},1}\rangle |$, $| \langle {c}_{{\bf{k}},1}{c}_{\tilde{{\bf{q}}}-{\bf{k}},2}\rangle |$ and $| \langle {c}_{{\bf{k}},2}{c}_{\tilde{{\bf{q}}}-{\bf{k}},2}\rangle |$ for h_x = 0 (a)–(c), h_x = 0.8 (d)–(f) and h_x = 0.9 (g)–(i) in case of λ = 0.75, μ = 0.95 and h_z = 0. These three cases correspond to three red squares of figure 1(e). For clarity, also the non-interacting Fermi surfaces are depicted as red (blue) contours for the upper (lower) branch. The case h_x = 0 shown in figures 3(a)–(c) corresponds to conventional BCS phase for which intraband pairing takes place within both bands and interband pairing is vanishingly small. When h_x is finite, the system enters first to the small-${\tilde{q}}_{y}$ region (figures 3(d)–(f)) where both intraband pairing contributions are still prominent and the interband pairing is finite but small. Due to the contribution of both bands, T_BKT is more or less the same as for h_x = 0, see figure 1(e). The only qualitative difference is the asymmetric pairing profiles of h_x = 0.8 which causes the finite momentum pairing to be more stable than the zero momentum BCS pairing.

Zoom In Zoom Out Reset image size

The situation is drastically different when the system enters to the large-${\tilde{q}}_{y}$ region at h_x = 0.9 (figures 3(g)–(i)). In contrast to cases with smaller h_x, the prominent intraband pairing contribution comes now from the upper band alone. As the pairing occurs only in one of the bands instead of both bands, T_BKT is significantly lower for the large-${\tilde{q}}_{y}$ region than for the small-${\tilde{q}}_{y}$ phase, as seen in figure 1(e).

It should be reminded that we consider FF states only and ignore LO states. In recent real-space mean-field studies [10, 26], it was pointed out that LO states are associated with finite pairing amplitudes occurring within both bands and correspondingly FF phases are a consequence of the pairing occurring within a single helicity band only. This is easy to understand as the in-plane Zeeman field shifts the other helicity band to $+{k}_{y}$ and the other to $-{k}_{y}$ direction. Therefore, when the pairing occurs within both bands, some pairing occurs with Cooper pair momentum $+{\tilde{q}}_{y}$ and some with $-{\tilde{q}}_{y}$ which results in an LO phase. Thus, the small-${\tilde{q}}_{y}$ region we find is likely the one where LO states are more stable than FF states and hence T_BKT is considerably higher for LO states than for FF states. Unfortunately, accessing LO states directly is not possible with our momentum space study as LO phases break the translational invariance which is utilized in the derivation of the superfluid weight as shown in section 2. For computing the superfluid weight also in case of LO ansatzes, one should derive the expressions for the superfluid weight by using real-space quantities only.

For completeness, in figure 4 we provide the phase diagrams for ${\tilde{q}}_{y}$ and T_BKT as functions of μ and h_z (figures 4(a), (b)) and of μ and h_x (figures 4(c), (d)) at λ = 0.75. In case of the (μ, h_z)-phase diagram the in-plane Zeeman field is fixed to h_x = 0.658 and in case of the (μ, h_x)-diagram the out-of-plane Zeeman field is h_z = 0.8. As in figure 1 with (h_x, h_z)-diagram, also here we find various topologically non-trivial FF and BCS phases identified with the Chern numbers C = −1 and C = −2 near the half-filling. However, for higher chemical potential values we also find topological FF and BCS phases characterized by C = 1. Furthermore, we can once again identify FF phases with high T_BKT but considerably small Cooper pair momenta existing near the half-filling with moderately low Zeeman field values. From figures 4(b) and (d) we see that for a non-topological FF phase T_BKT is 0.1–0.3t at relatively large parameter regime. For topological FF states T_BKT is somewhat lower, the maximum transition temperature being T_BKT ∼ 0.15t.

Zoom In Zoom Out Reset image size

In previous FFLO studies [5, 40, 51] it has been shown that Van Hove singularities associated with the divergent behavior of the density of states near the Fermi surface can enlarge the parameter regime of FFLO states. In our spin–orbit-coupled square lattice system there are six different Van Hove singularities for fixed μ. In figures 4(a), (b) two of these singularities are depicted with red dashed–dotted lines, the other four occurring near the depicted two. One can see that in the vicinity of the Van Hove singularities the FF phases can exist at higher values of h_z than away from the singularities. However, in (μ, h_x)-diagrams depicted in figures 4(c), (d) the Van Hove singularities are not playing a role and therefore they are not shown.

4.2. Topological phase transitions

Topological phase diagrams presented here and in [26] for a square lattice are relatively rich compared to the topological phase diagrams of Rashba-coupled 2D continuum where they are characterized by C = 1 only. This can be explained by considering possible topological phase transitions which occur when the bulk energy gap E_g between the quasi-particle eigenvalues ${E}_{{\bf{k}},\nu }^{+}$ and quasi-holes ${E}_{{\bf{k}},\nu }^{-}$ closes and reopens. Because of the intrinsic particle–hole symmetry present in our system, topological phase transitions can occur when the gap closes and reopens in particle–hole symmetric points [57]. In continuum there exists only one particle–hole symmetric point, i.e. ${\bf{k}}=({k}_{x},{k}_{y})=(0,{\tilde{q}}_{y}/2)$. However, in a square lattice there are four different particle–hole symmetric points, namely ${{\bf{k}}}_{1}=(0,{\tilde{q}}_{y}/2)$, ${{\bf{k}}}_{2}=(0,-\pi +{\tilde{q}}_{y}/2)$, ${{\bf{k}}}_{3}=(\pi ,{\tilde{q}}_{y}/2)$ and ${{\bf{k}}}_{4}=(\pi ,-\pi +{\tilde{q}}_{y}/2)$ which yields four different gap closing equations instead of only one. Therefore, it is reasonable to find more distinct topological phases in a lattice system than in continuum. For similar reasons, topological phase diagrams studied in [27] in case of triangular lattices possessed many distinct topological states characterized by different Chern numbers. Analytical gap closing equations for the square lattice geometry are provided in appendix D.

In figures 5(a)–(c) we plot the minimum energy gap E_g for (h_x, h_z), (μ, h_z) and (μ, h_x)-phase diagrams shown previously in figures 1(b), 4(a) and (c). One can see that E_g goes to zero at the topological phase boundaries as expected. In figures 5(a)–(c) we also depict the fulfilled analytical gap closing conditions which match with numerically computed topological boundaries. Analytical gap closing conditions can be thus used to identify distinct topological transitions in terms of the gap closing locations in the momentum space.

Zoom In Zoom Out Reset image size

From figures 5(a)–(c) we see that the Chern invariant changes by one when the gap closes in one of the particle–hole symmetric momenta. However, when the system enters from the trivial C = 0 phase to C = −2 phase, the gap closes simultaneously in two different momenta. This is consistent with the theory presented in [57] considering the connection between the Chern number and gap closings at particle–hole symmetric points: if the Chern number changes by an even (odd) number at a topological phase transition, then the number of gap closing particle–hole symmetric momenta is even (odd).

We further investigate the topological phase transitions in figures 5(d)–(l), where we present the momentum density distributions ${n}_{{\bf{k}}}={n}_{\uparrow {\bf{k}}}+{n}_{\downarrow {\bf{k}}}=\langle {c}_{\uparrow {\bf{k}}}^{\dagger }{c}_{\uparrow {\bf{k}}}\rangle +\langle {c}_{\downarrow {\bf{k}}}^{\dagger }{c}_{\downarrow {\bf{k}}}\rangle$ for six different values of μ, corresponding to six yellow dots depicted in figure 5(c). The topological transition corresponding to the gap closing at ${{\bf{k}}}_{3}$ is studied in figures 5(d), (e), and correspondingly closings at ${{\bf{k}}}_{2}$ and ${{\bf{k}}}_{4}$ are investigated in figures 5(g)–(i) and (j)–(l), respectively.

By comparing the momentum distributions in figures 5(d), (e) shown for μ = 0.792 and μ = 0.912, we observe that once the system goes through the topological transition identified by the gap closing and reopening at ${{\bf{k}}}_{3}$ (white line in figure 5(c)), the momentum distribution changes qualitatively in the vicinity of ${{\bf{k}}}_{3}$. This is further shown in figure 5(f) where ${n}_{{\bf{k}}}$ for both cases is plotted at k_y = 0 along the blue dashed–dotted line depicted in figures 5(d), (e). In a similar fashion, one sees from figures 5(g)–(i) that the topological transition corresponding to the gap closing at ${{\bf{k}}}_{2}$ (red line in figure 5(c)) is identified as an emergence of a prominent density peak around ${{\bf{k}}}_{2}$ as clearly illustrated in figure 5(i). A similar peak can be also observed for the topological transition corresponding to ${{\bf{k}}}_{4}$ though less pronounced as shown in figures 5(j)–(l).

Drastic qualitative changes in the momentum distributions at the topological phase boundaries imply that one could experimentally measure and distinguish different topological phases and phase transitions in ultracold gas systems by investigating the total density distributions with the time-of-flight measurements. A similar idea to measure topological phase transitions were proposed in [14] in case of a simpler continuum system. Our findings show that density measurements could be applied also in lattice systems to resolve different topological phases.

4.3. Components of the superfluid weight

As the single-particle energy dispersions are deformed in the y-direction but not in the x-direction, the rotational symmetry of the lattice is broken. This manifests itself as different superfluid weight components in the x- and y-directions, i.e. ${D}_{{xx}}^{s}\ne {D}_{{yy}}^{s}$. As the Cooper pair momentum is in the y-direction, we call D^s_yy as the longitudinal and D^s_xx as the transverse component. Because ${D}_{{xx}}^{s}\ne {D}_{{yy}}^{s}$, the system has different current response in these directions when exposed to an external magnetic field. Therefore, it is meaningful to investigate the difference of the longitudinal and transverse components, ${D}_{{\rm{diff}}}^{s}\equiv {D}_{{yy}}^{s}-{D}_{{xx}}^{s}$, to see how it behaves as a function of our system parameters. We focus only on the diagonal elements of D^s as the off-diagonal elements in our case are always zero, i.e. ${D}_{{xy}}^{s}={D}_{{yx}}^{s}=0$.

In figures 6(a)–(c) we present ${D}_{{\rm{diff}}}^{s}$ for $({h}_{x},{h}_{z})$, (μ, h_z) and (μ, h_x)-phase diagrams, respectively, shown above in figures 1(b), 4(a) and (c). In all three cases, ${D}_{{\rm{diff}}}^{s}$ more or less vanishes in large parts of the phase diagrams. However, especially when entering the large-${\tilde{q}}_{y}$ FF region from the small-${\tilde{q}}_{y}$ region, ${D}_{{\rm{diff}}}^{s}$ reaches local minima and becomes negative. On the other hand, from figures 6(b), (c) we see that there also exists a parameter region where ${D}_{{\rm{diff}}}^{s}$ is positive and that the tFF₋₂-phase in figure 6(c) near half-filling is clearly distinguishable from the neighboring phases. Therefore, by measuring ${D}_{{\rm{diff}}}^{s}$ one could in principle distinguish some of the phase transitions existing in the system. It is interesting to note that, in the presence of SOC, the transverse component can be larger than the longitudinal component, in contrast to 2D continuum where the absence of SOC results in the vanishing transverse component and thus the vanishing BKT temperature T_BKT = 0 [47].

Zoom In Zoom Out Reset image size

In addition to ${D}_{{\rm{diff}}}^{s}$, in figures 6(a)–(c) we also plot with solid white lines the boundaries of gapped and gapless superfluid phases. Consistent with previous literature [12–14, 48], we call the system gapless (or nodal) if one or more of the Bogoliubov quasi-hole branches reach the zero-energy in some part of the momentum space, i.e. the quasi-particle excitation energy vanishes for some momenta. Note that this does not (necessarily) mean that the topological energy gap E_g closes as E_g is the difference of the highest quasi-hole and the lowest quasi-particle energy at the same momentum ${\bf{k}}$ such that both are also the eigenvalues of ${{ \mathcal H }}_{{\bf{k}}}$, whereas the highest quasi-hole and the lowest quasi-particle energy are not necessarily at the same momentum.

From figures 6(a) and (c) we see that the system stays gapped at low in-plane Zeeman field strengths which is consistent with continuum results [12, 48]. For larger h_x the system becomes eventually gapless and one can observe topologically trivial and non-trivial nodal FF phases. By comparing figures 1(b), 4(a) and (c) to 6(a)–(c) we can make a remark that FF states with small momenta ${\tilde{q}}_{y}$ are gapped. Furthermore, we observe from figures 6(a)–(c) that the transitions between the gapped and gapless states at moderate Zeeman fields and chemical potentials coincide with the prominent minima of ${D}_{{\rm{diff}}}^{s}$. This is consistent with the findings of [48] where it was shown that the longitudinal component exhibits a clear minimum when the system becomes gapless. However, in figures 6(b), (c) we see the system reaching a gapped region again at large enough μ without such a drastic change of ${D}_{{\rm{diff}}}^{s}$ than at smaller values of μ.

In addition to different spatial components, one can also investigate the role of the geometric superfluid weight contribution ${D}_{\mathrm{geom}}^{s}$ which is presented for (h_x, h_z), (μ, h_z) and (μ, h_x)-phase diagrams in figures 6(d)–(f). We see that for BCS states and gapped FF states of small Cooper pair momenta, the geometric contribution is notable but is otherwise vanishingly small. In all the cases the geometric contribution is relatively small compared to the total superfluid weight D^s which is, as an example, illustrated in the inset of figure 6(f) where ${D}_{\mathrm{geom}}^{s}$ and D^s are both plotted for h_x = 0. At largest, the geometric contribution is responsible up to 18% of the total superfluid weight which is fairly similar to what was reported in [58], where the geometric part was found to contribute up to a quarter of the total superfluid weight in case of a spin–orbit-coupled 2D BCS continuum model. In more complicated multiband lattices, such as honeycomb lattice or Lieb lattice (which also possesses a flat band), the geometric contribution in the presence of SOC might be more important than in our simple square lattice example as the geometric contribution is intrinsically a multiband effect [54].

In this work we have investigated the stability of exotic FF superfluid states in a lattice system by computing the superfluid weight and BKT transition temperatures systematically for various system parameters. The derivation of the superfluid weight is based on the linear response theory and is an extension of the previous studies of [54, 55] where only BCS ansatzes without spin-flipping terms were considered. Our method applies to BCS and FF states in the presence of arbitrary spin-flipping processes and lattice geometries. We find that, as previously in case of conventional BCS theory without the spin-flipping contribution, also in case of FF phases and with spin-flipping terms one can divide the total superfluid weight to conventional and geometric superfluid contributions.

We have focused on a square lattice geometry in the presence of the Rashba-coupling. One of the main findings of this article is that conventional spin-imbalance-induced FF states, in the absence of SOC, indeed have finite BKT transition temperatures in a lattice geometry. For our parameters they could be observed at $T\sim 0.1t$. In earlier theoretical studies it has been predicted that FF states could exist in two-dimensional lattice systems [5, 50, 51, 59] but the stability in terms of the BKT transition has never been investigated in lattice systems. By computing T_BKT we show that two-dimensional FFLO superfluids should be realizable in finite temperatures. By applying SOC, we show that FF states in a lattice can be further stabilized and for our parameter regime BKT temperatures as high as $T\sim 0.17\mbox{--}0.3t$ can be reached. SOC also enables the existence of topological nodal and gapped FF states, for which we show the BKT transitions to occur at highest around ${T}_{{\rm{BKT}}}\sim 0.15t$.

For literature comparison, we estimated that ${T}_{{\rm{BKT}}}\approx 0.25t$ at $U=-4t$ for usual spin-balanced BCS state at half-filling without SOC, see figure 1(c), whereas in [60] the corresponding estimate obtained by Monte Carlo simulations was ${T}_{{\rm{BKT}}}\sim 0.10\mbox{--}0.13t$. Thus, our mean-field approach probably overestimates T_BKT in case of a simple square lattice. However, in [55, 61] the superfluid weights of BCS states, derived in the framework of mean-field theory, were shown to agree reasonably well with more sophisticated theoretical methods in case of multiband systems. Thus, it is expected that our mean-field superfluid equations are in better agreement with beyond-mean-field methods when considering multiband lattice models.

We have also shown that different topological FF phases and phase transitions could be observed by investigating the total momentum density profiles. When the system goes through a topological phase transition, the momentum distribution develops peaks or dips in the vicinity of momenta in which the energy gap closes and reopens. In addition to density distributions, also the relative behavior of the longitudinal and transverse superfluid weight components yields implications about the phase transitions, especially near the boundaries of gapless and gapped superfluid phases. Therefore, our work paves the way for stabilizing and identifying exotic topological FF phases in lattice systems.

In future studies it would be interesting to see how stable FF states are in multiband models. This could be investigated straightforwardly with our superfluid weight equations as they hold for an arbitrary multiband system. Especially intriguing could be systems which possess both dispersive and flat bands such as kagome or Lieb lattices. In these systems the conventional spin-imbalanced FF states were recently shown to exhibit exotic deformation of Fermi surfaces due to the presence of a flat band [62]. In multiband systems one could also expect the geometric superfluid contribution to play a role, in contrast to our square lattice system where the geometric contribution was only non-zero for BCS and gapped FF phases. Furthermore, in flat band systems mean-field theory is shown to be in good agreement with more advanced beyond-mean-field approaches [55, 61, 63]. Flat band systems are tempting also because it is expected that their superfluid transition temperatures in the weak-coupling region are higher than in dispersive systems [54, 55, 61, 64, 65] and thus they could provide a way to realize exotic FFLO phases at high temperatures.

This work was supported by the Academy of Finland through its Centres of Excellence Programme (2012–2017) and under project NOs. 284621, 303351 and 307419, and by the European Research Council (ERC-2013-AdG-340748-CODE). AJ acknowledges support from the Vilho, Yrjö and Kalle Väisälä foundation.

Here we briefly go through how one obtains the final form for the superfluid weight D^s shown in (15) from the intermediate result (14). As one can see from (14), there exists two terms in K_μν, the first being the diamagnetic and the second one the paramagnetic contribution, ${K}_{\mu \nu ,{\rm{dia}}}$, ${K}_{\mu \nu ,{\rm{para}}}$, respectively. We focus on the diamagnetic term and after that just give the result for the paramagnetic term as the derivation for both terms is essentially the same.

In the diamagnetic term there exists a double derivative ${\partial }_{\mu }{\partial }_{\nu }{ \mathcal H }({\bf{k}})$ which can be transformed to a single derivative via integrating by parts:

$$\begin{eqnarray}{K}_{\mu \nu ,{\rm{dia}}} & = & \displaystyle \frac{1}{\beta }\displaystyle \sum _{{\bf{k}},{{\rm{\Omega }}}_{m}}{\rm{Tr}}[{\partial }_{\mu }{\partial }_{\nu }{ \mathcal H }({\bf{k}}){P}_{+}G({\rm{i}}{{\rm{\Omega }}}_{m},{\bf{k}})]\\ & = & -\displaystyle \frac{1}{\beta }\displaystyle \sum _{{\bf{k}},{{\rm{\Omega }}}_{m}}{\rm{Tr}}[{\partial }_{\mu }{ \mathcal H }({\bf{k}}){P}_{+}{\partial }_{\nu }G({\rm{i}}{{\rm{\Omega }}}_{m},{\bf{k}})].\end{eqnarray} \tag{ A1 }$$

Because $G({\rm{i}}{{\rm{\Omega }}}_{m},{\bf{k}})=1/({\rm{i}}{{\rm{\Omega }}}_{m}-{ \mathcal H }({\bf{k}}))$, we have ${\partial }_{\nu }{G}^{-1}=-{\partial }_{\nu }{ \mathcal H }$ and because ${\partial }_{\nu }({{GG}}^{-1})=0$ we also have ${\partial }_{\nu }G=-G{\partial }_{\nu }{G}^{-1}G$ so that (A1) can be written as

$$\begin{eqnarray}{K}_{\mu \nu ,{\rm{dia}}} & = & -\displaystyle \frac{1}{\beta }\displaystyle \sum _{{\bf{k}},{{\rm{\Omega }}}_{m}}{\rm{Tr}}[{\partial }_{\mu }{ \mathcal H }({\bf{k}}){P}_{+}G({\rm{i}}{{\rm{\Omega }}}_{m},{\bf{k}}){\partial }_{\nu }{ \mathcal H }({\bf{k}})G({\rm{i}}{{\rm{\Omega }}}_{m},{\bf{k}})]\\ & = & -\displaystyle \frac{1}{\beta }\displaystyle \sum _{{\bf{k}},{{\rm{\Omega }}}_{m}}\displaystyle \sum _{i}^{4M}\langle {\phi }_{i}({\bf{k}})| {\partial }_{\mu }{ \mathcal H }({\bf{k}}){P}_{+}G({\rm{i}}{{\rm{\Omega }}}_{m},{\bf{k}}){\partial }_{\nu }{ \mathcal H }({\bf{k}})G({\rm{i}}{{\rm{\Omega }}}_{m},{\bf{k}})| {\phi }_{i}({\bf{k}})\rangle ,\end{eqnarray} \tag{ A2 }$$

where $| {\phi }_{i}({\bf{k}})\rangle$ are the eigenvectors of ${{ \mathcal H }}_{{\bf{k}}}$. By using the completeness relation ${\sum }_{j}| {\phi }_{j}({\bf{k}})\rangle \langle {\phi }_{j}({\bf{k}})| =1$ and the alternative form for $G({\rm{i}}{{\rm{\Omega }}}_{m},{\bf{k}})$

$$\begin{eqnarray}&&G({\rm{i}}{{\rm{\Omega }}}_{m},{\bf{k}})=\displaystyle \sum _{l=1}^{4M}\displaystyle \frac{| {\phi }_{l}({\bf{k}})\rangle \langle {\phi }_{l}({\bf{k}})| }{{\rm{i}}{{\rm{\Omega }}}_{m}-{E}_{l,{\bf{k}}}}\end{eqnarray} \tag{ A3 }$$

we obtain

$$\begin{eqnarray}{K}_{\mu \nu ,{\rm{dia}}} & = & -\displaystyle \frac{1}{\beta }\displaystyle \sum _{{\bf{k}},{{\rm{\Omega }}}_{m}}\displaystyle \sum _{i,j}^{4M}\langle {\phi }_{i}({\bf{k}})| {\partial }_{\mu }{ \mathcal H }({\bf{k}}){P}_{+}| {\phi }_{j}({\bf{k}})\rangle \\ & & \times \langle {\phi }_{j}({\bf{k}})| {\partial }_{\nu }{ \mathcal H }({\bf{k}})| {\phi }_{i}({\bf{k}})\rangle \displaystyle \frac{1}{({\rm{i}}{{\rm{\Omega }}}_{m}-{E}_{j,{\bf{k}}})({\rm{i}}{{\rm{\Omega }}}_{m}-{E}_{i,{\bf{k}}})}.\end{eqnarray} \tag{ A4 }$$

The summation over the Matsubara frequencies Ω_m can be carried out analytically yielding

$$\begin{eqnarray}{K}_{\mu \nu ,{\rm{dia}}} & = & \displaystyle \sum _{{\bf{k}},{ij}}\langle {\phi }_{i}({\bf{k}})| {\partial }_{\mu }{ \mathcal H }({\bf{k}}){P}_{+}| {\phi }_{j}({\bf{k}})\rangle \\ & & \times \langle {\phi }_{j}({\bf{k}})| {\partial }_{\nu }{ \mathcal H }({\bf{k}})| {\phi }_{i}({\bf{k}})\rangle \displaystyle \frac{n({E}_{j,{\bf{k}}})-n({E}_{i,{\bf{k}}})}{{E}_{i,{\bf{k}}}-{E}_{j,{\bf{k}}}}.\end{eqnarray} \tag{ A5 }$$

In a similar fashion one derives the following result for the paramagnetic term:

$$\begin{eqnarray}{K}_{\mu \nu ,{\rm{para}}}({\bf{q}}\to 0,0) & = & -\displaystyle \sum _{{\bf{k}},{ij}}\langle {\phi }_{i}({\bf{k}})| {\partial }_{\mu }{ \mathcal H }({\bf{k}}){P}_{+}| {\phi }_{j}({\bf{k}})\rangle \\ & & \times \langle {\phi }_{j}({\bf{k}})| {\partial }_{\nu }{ \mathcal H }({\bf{k}}){\hat{\gamma }}_{z}| {\phi }_{i}({\bf{k}})\rangle \displaystyle \frac{n({E}_{j,{\bf{k}}})-n({E}_{i,{\bf{k}}})}{{E}_{i,{\bf{k}}}-{E}_{j,{\bf{k}}}}.\end{eqnarray} \tag{ A6 }$$

As ${D}_{\mu \nu }^{s}={K}_{\mu \nu }({\bf{q}}\to 0,0)={K}_{\mu \nu ,{\rm{dia}}}+{K}_{\mu \nu ,{\rm{para}}}({\bf{q}}\to 0,0)$ and ${P}_{-}=({I}_{4M}-{\hat{\gamma }}_{z})/2$, one readily obtains the final result presented in (15).

In this appendix we show how the total superfluid weight D^s presented in (15) can be split to the so-called conventional and geometric contributions, ${D}_{{\rm{conv}}}^{s}$ and ${D}_{{\rm{geom}}}^{s}$. We start by expressing the eigenvectors $| {\phi }_{i}({\bf{k}})\rangle$ of ${ \mathcal H }({\bf{k}})$ in terms of the eigenvectors of ${{ \mathcal H }}_{p}({\bf{k}})$ and ${{ \mathcal H }}_{h}({\bf{k}})$ as follows

$$\begin{eqnarray}&&| {\phi }_{i}({\bf{k}})\rangle =\displaystyle \sum _{m=1}^{2M}({w}_{p,{im}}| +\rangle \otimes | m{\rangle }^{p}+{w}_{h,{im}}| -\rangle \otimes | m{\rangle }^{h}),\end{eqnarray} \tag{ B1 }$$

where $| m{\rangle }^{p}$ ($| m{\rangle }^{h}$) are the eigenvectors of ${{ \mathcal H }}_{p}$ (${{ \mathcal H }}_{h}$) and $| \pm \rangle$ are the eigenvectors of $\hat{{\sigma }_{z}}\otimes {I}_{2M}$ with the eigenvalues±1. By noting that

$$\begin{eqnarray}{\partial }_{\mu }{ \mathcal H }({\bf{k}})=\left[\begin{array}{cc}{\partial }_{\mu }{{ \mathcal H }}_{p}({\bf{k}}) & 0\\ 0 & -{\partial }_{\mu }{{ \mathcal H }}_{h}({\bf{k}}-\tilde{{\bf{q}}})\end{array}\right]\end{eqnarray} \tag{ B2 }$$

we can rewrite (15) as

$$\begin{eqnarray}{D}_{\mu \nu }^{s} & = & \displaystyle \sum _{{\bf{k}},{ij}}\displaystyle \frac{n({E}_{j})-n({E}_{i})}{{E}_{i}-{E}_{j}}\\ & & \times \displaystyle \sum _{{m}_{1},{m}_{2}}^{2M}[{w}_{p,{{im}}_{1}}^{* }{w}_{p,{{jm}}_{2}}{}^{p}\langle {m}_{1}| {\partial }_{\mu }{{ \mathcal H }}_{p}({\bf{k}})| {m}_{2}{\rangle }^{p}]\\ & & \times \displaystyle \sum _{{m}_{3},{m}_{4}}^{2M}[{w}_{h,{{jm}}_{3}}^{* }{w}_{h,{{im}}_{4}}{}^{h}\langle {m}_{3}| -{\partial }_{\nu }{{ \mathcal H }}_{h}({\bf{k}}-\tilde{{\bf{q}}})| {m}_{4}{\rangle }^{h}]\\ & = & \displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{\genfrac{}{}{0em}{}{{\bf{k}}}{{m}_{1},{m}_{2},}}{{m}_{3},{m}_{4}}}{W}_{{m}_{1}{m}_{2}}^{{m}_{3}{m}_{4}}({}^{p}\langle {m}_{1}| {\partial }_{\mu }{{ \mathcal H }}_{p}| {m}_{2}{\rangle }^{{ph}}\langle {m}_{3}| -{\partial }_{\nu }{{ \mathcal H }}_{h}| {m}_{4}{\rangle }^{h}),\end{eqnarray} \tag{ B3 }$$

where

$$\begin{eqnarray}&&{W}_{{m}_{1}{m}_{2}}^{{m}_{3}{m}_{4}}=\displaystyle \sum _{{ij}}\displaystyle \frac{n({E}_{j})-n({E}_{i})}{{E}_{i}-{E}_{j}}{w}_{p,{{im}}_{1}}^{* }{w}_{p,{{jm}}_{2}}{w}_{h,{{jm}}_{3}}^{* }{w}_{h,{{im}}_{4}}\end{eqnarray} \tag{ B4 }$$

and

$$\begin{eqnarray}&&{}^{p}\langle {m}_{1}| {\partial }_{\mu }{{ \mathcal H }}_{p}| {m}_{2}{\rangle }^{p}={\delta }_{{m}_{1},{m}_{2}}{\epsilon }_{{m}_{1}}+({\epsilon }_{{m}_{1}}-{\epsilon }_{{m}_{2}})\,{}^{p}\langle {\partial }_{\mu }{m}_{1}| {m}_{2}{\rangle }^{p}.\end{eqnarray} \tag{ B5 }$$

Here ${\epsilon }_{{m}_{i}}$ are the eigenvalues for ${{ \mathcal H }}_{p}$. Similar expression holds also for the ${}^{h}\langle {m}_{3}| -{\partial }_{\nu }{{ \mathcal H }}_{h}| {m}_{4}{\rangle }^{h}$ elements. From (B3)–(B5) we note that there exists two superfluid weight components. The component which is called the conventional contribution ${D}_{{\rm{conv}}}^{s}$ consists of matrix elements with m₁ = m₂ and m₃ = m₄. As can be seen from (B5), the conventional contribution depends only on the single-particle dispersions ${\epsilon }_{{m}_{i}}$. The remaining part is the geometric contribution ${D}_{{\rm{geom}}}^{s}$ and it depends on the geometric properties of the Bloch functions, $| {m}_{i}{\rangle }^{p}$ and $| {m}_{i}{\rangle }^{h}$.

As our equations for the superfluid weight hold for arbitrary geometries in the presence and absence of SOC, we can make direct comparisons to previous studies. As the first benchmark, we reproduced the superfluid weight results of [61] where BCS states in the Lieb lattice geometry without the SOC are studied by applying mean-field theory and exact diagonalization (ED) methods. One should emphasize that mean-field equations used in [61] to compute the superfluid weight are derived by not using the linear response theory as in our study but by using an alternative approach based on the definition given in [54]. Our method yields exactly the same results as the alternative mean-field and ED approaches of [61]. Furthermore, we have checked that in the continuum limit our expression for the superfluid weight reduces to the expressions presented in [58] where BCS states in spin–orbit-coupled 2D continuum were considered.

We also benchmarked our equations by computing T_BKT in case of BCS phases for a 2D square lattice geometry with the same parameters that were used in [66] where topological BCS states in the presence of the SOC were studied. With our equations we find the same functional behavior for T_BKT as a function of U but our results are exactly a factor of two larger than those presented in [66]. The reason for this difference is because in [66], the phase fluctuations of the order parameter are rescaled by a factor of $1/\sqrt{2}$ (see equation (33) in [66]). With this rescaling, the periodicity of the ϕ field in (38) becomes $2\sqrt{2}\pi$ and therefore the expression for the BKT transition temperature (equation (39)) should be multiplied by a factor of 2.

In this appendix we show the analytical equations that were used to depict the topological phase transitions in figure 5. The energy gap E_g between the quasi-particle eigenvalues ${E}_{{\bf{k}},\nu }^{+}$ and quasi-holes ${E}_{{\bf{k}},\nu }^{-}$ can only close and reopen at particle–hole-symmetric points which in our case are ${{\bf{k}}}_{1}=(0,{\tilde{q}}_{y}/2)$, ${{\bf{k}}}_{2}=(0,-\pi +{\tilde{q}}_{y}/2)$, ${{\bf{k}}}_{3}=(\pi ,{\tilde{q}}_{y}/2)$ and ${{\bf{k}}}_{4}=(\pi ,-\pi +{\tilde{q}}_{y}/2)$. The single-particle Hamiltonian ${{ \mathcal H }}_{p}$ in these four points can be diagonalized analytically which yields four eigenvalues, namely ${E}_{{\bf{k}},1}^{-}\leqslant {E}_{{\bf{k}},2}^{-}\leqslant {E}_{{\bf{k}},2}^{+}\leqslant {E}_{{\bf{k}},1}^{+}$. By demanding ${E}_{{\bf{k}},2}^{-}={E}_{{\bf{k}},2}^{+}$ at each of the four particle–hole symmetric momenta, one obtains the four gap closing equations which read

$$\begin{eqnarray}{h}_{z}^{2} & = & 6+{{\rm{\Delta }}}^{2}+4\mu +{\mu }^{2}+4(2+\mu )\cos ({\tilde{q}}_{y}/2)+2\cos ({\tilde{q}}_{y})-{h}_{x}^{2}+2{\lambda }^{2}[\cos ({\tilde{q}}_{y})-1]\\ & & +4{h}_{x}\lambda \sin ({\tilde{q}}_{y}/2),\end{eqnarray} \tag{ D1 }$$

$$\begin{eqnarray}{h}_{z}^{2} & = & 6+{{\rm{\Delta }}}^{2}\,+\,4\mu +{\mu }^{2}-4(2+\mu )\cos ({\tilde{q}}_{y}/2)+2\cos ({\tilde{q}}_{y})-{h}_{x}^{2}+2{\lambda }^{2}[\cos ({\tilde{q}}_{y})-1]\\ & & -4{h}_{x}\lambda \sin ({\tilde{q}}_{y}/2),\end{eqnarray} \tag{ D2 }$$

$$\begin{eqnarray}{h}_{z}^{2} & = & 6+{{\rm{\Delta }}}^{2}-4\mu +{\mu }^{2}+4(2+\mu )\cos ({\tilde{q}}_{y}/2)+2\cos ({\tilde{q}}_{y})-{h}_{x}^{2}+2{\lambda }^{2}[\cos ({\tilde{q}}_{y})-1]\\ & & +4{h}_{x}\lambda \sin ({\tilde{q}}_{y}/2),\end{eqnarray} \tag{ D3 }$$

$$\begin{eqnarray}{h}_{z}^{2} & = & 6+{{\rm{\Delta }}}^{2}-4\mu +{\mu }^{2}-4(2+\mu )\cos ({\tilde{q}}_{y}/2)+2\cos ({\tilde{q}}_{y})-{h}_{x}^{2}+2{\lambda }^{2}[\cos ({\tilde{q}}_{y})-1]\\ & & -4{h}_{x}\lambda \sin ({\tilde{q}}_{y}/2).\end{eqnarray} \tag{ D4 }$$

By solving these equations for different values of h_x, h_z and μ, one obtains the topological boundaries shown in figures 5(a)–(c).

In our computations the Cooper pair momentum $\tilde{{\bf{q}}}$ is in the y-direction, i.e. $\tilde{{\bf{q}}}\parallel {\hat{{\bf{e}}}}_{y}$, consistent with earlier studies concerning lattice systems [10, 26, 27]. We have extensively tested numerically that indeed the wavevector in the y-direction minimizes the thermodynamic potential with and without SOC for all the used input parameters. As an example, we have demonstrated this in figure E1. In figures E1(a)–(c) we plot the (μ, h_x)-phase diagram for three different cases: in (a) the thermodynamic potential Ω is minimized so that $\tilde{{\bf{q}}}$ is taken to be in the y-direction, in (b) $\tilde{{\bf{q}}}$ is along the diagonal direction (${\tilde{q}}_{x}={\tilde{q}}_{y}$) and in (c) $\tilde{{\bf{q}}}$ is in the x-direction. The out-of-plane Zeeman field is chosen to be h_z = 0.8, the spin–orbit-coupling is λ = 0.75 and the interaction strength is U = −4 so the phase diagram in figure E1(a) is the same as in figure 4(c) in the main text. We see how gradually the FF region becomes smaller when the wavevector is forced to deviate from the y-direction. In figures E1(d), (e) we compare the thermodynamic potentials Ω of these three different cases. In figure E1(d) the thermodynamic potential difference of cases $\tilde{{\bf{q}}}\parallel {\hat{{\bf{e}}}}_{x}+{\hat{{\bf{e}}}}_{y}$ and $\tilde{{\bf{q}}}\parallel {\hat{{\bf{e}}}}_{y}$ is plotted and correspondingly in figure E1(e) the thermodynamic potential difference of cases $\tilde{{\bf{q}}}\parallel {\hat{{\bf{e}}}}_{x}$ and $\tilde{{\bf{q}}}\parallel {\hat{{\bf{e}}}}_{y}$ is depicted. White lines show the phase boundaries between the BCS, FF and normal phases in case of $\tilde{{\bf{q}}}\parallel {\hat{{\bf{e}}}}_{y}$. We see that within the BCS phase the thermodynamic potential is the same regardless of the direction of the wavevector as in the BCS phase the Cooper pair momentum is zero. When entering the FF phase, it is clear that phase diagrams shown in figures E1(b), (c) do not depict the true ground states as their thermodynamic potentials are higher than in case of $\tilde{{\bf{q}}}\parallel {\hat{{\bf{e}}}}_{y}$. Thus the states shown in figure E1(a) with $\tilde{{\bf{q}}}\parallel {\hat{{\bf{e}}}}_{y}$ are energetically more stable than the states with the Cooper pair momentum in the diagonal or x-direction.

Zoom In Zoom Out Reset image size

In figure E1 we have only presented three different options for the direction of $\tilde{{\bf{q}}}$ and only (μ, h_x)-phase diagram. However, they represent the general trend of all the computations of our work: the thermodynamic potential reaches its minimum when $\tilde{{\bf{q}}}$ is in the y-direction. We have confirmed this by choosing 20 other directions between the x and y-axes. Alternatively, we also minimized the thermodynamic potential by letting q_x and q_y be independent parameters. As the thermodynamic potential can have many local minima as a function of q_x and q_y, this procedure is not the most trustworthy for finding the global minimum. However, we did not find a single local minimum lying outside the y-axis that would have lower energy than the solutions we find by assuming $\tilde{{\bf{q}}}| | {\hat{{\bf{e}}}}_{y}$. Therefore we are confident that our statements and results are correct within the mean-field theory framework.