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

We study the superfluid weight $D^s$ and Berezinskii-Kosterlitz-Thouless (BKT) transition temperatures $T_{BKT}$ in case of exotic Fulde-Ferrell (FF) superfluid states in lattice systems. We consider spin-imbalanced systems with and without spin-orbit coupling (SOC) accompanied with in-plane Zeeman field. By applying mean-field theory, we derive general equations for $D^s$ and $T_{BKT}$ in the presence of SOC and the Zeeman fields for 2D Fermi-Hubbard lattice models, and apply our results to a 2D square lattice. We show that conventional spin-imbalanced FF states without SOC can be observed at finite temperatures and that FF phases are further stabilized against thermal fluctuations by introducing SOC. We also propose how topologically non-trivial SOC-induced FF phases could be identified experimentally by studying the total density profiles. Furthermore, the relative behavior of transverse and longitudinal superfluid weight components and the role of the geometric superfluid contribution are discussed.


Introduction
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 , 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][30][31][32]. Ultracold gas experiments 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.

Derivation of the superfluid weight in the presence of SOC for an arbitrary lattice geometry
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- , M being the number of orbitals within a unit cell. Furthermore,  s ( ) k and L( ) 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 Y k and rewrite the Hamiltonian as follows: Here t s = Ä I y y M , where I M is a M×M identity matrix and s s s s = [ˆˆˆ] , , x y z are the Pauli matrices. One should note that now the single-particle Hamiltonian is not anymore simply  s but  p in which the two spin components are coupled via L( ) k . In two-dimensions the total superfluid weight D s is a 2×2 tensor which reads where x and y are the spatial dimensions. To compute the superfluid weight tensor elements mn D s , we exploit the fact that at the mean-field level mn D s is the long-wavelength, zero-frequency limit of the current-current response function mn K [53], that is . We are interested in computing the current-current response function w mn ( ) K q, which at the limit of  q 0, ω=0 yields the superfluid weight mn D s . To this end, we first define a Green's function In the Matsubara frequency space this reads 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 currentcurrent response function in a similar fashion as done in [55]. One starts from (11), inserts the expressions (12), (13)  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]: 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 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, = + 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 mn D s conv, is a function of the single-particle dispersions of  p and  h , and correspondingly

Rashba spin-orbit-coupled fermions in a square lattice
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. 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 Here d ij is the unit vector connecting the nearest-neighbor sites i and j, s s s s = [ˆˆˆ] , , x . We determine the order parameter amplitude Δ and the Cooper pair momentumq self-consistently by minimizing the grand canonical thermodynamic potential log Tr e B H which in the mean-field framework at T=0 reads as where Q( ) x is the Heaviside step function and n h E k, are the eigenvalues of  k . Here h = + -{ } , labels the quasiparticle and quasi-hole branches, respectively and ν={1, 2} the helicity branches split by the SOC. The quasiparticle branches are taken to be the two highest eigenvalues of  k . In (22) we have discarded the constant term which is not needed when one minimizes W M.F. . Consistent with previous lattice studies [10,26,27], the Cooper pair momentum is in the y-direction, i.e. =˜q q e y 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 Δ andq 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 G n h ( ) k associated with the quasi-hole branches η=− over the first Brillouin zone as follows: The explicit form for the Berry curvature can be expressed with the eigenvalues n h E k, of  k and the corresponding eigenvectors

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 m 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 letq y to have only discrete values in the first Brillouin zone such that , 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 ofq y are presented as a function of h x and h z at m = 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 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 t 0.25 BKT and roughly » T t 0.1 BKT 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.
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 −1 ) and C=−2 (tFF −2 ). By applying large h z one is able to reach topological BCS and FF phases, tBCS −2 and tFF −2 , 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 momentumq y : for intermediate Zeeman field strengths the FF state is characterized by rather smallq 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-q y region T BKT is around t 0.3 and becomes only smaller for large-q y region where T BKT at largest is roughly » T t 0.17 BKT . 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 andq 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  annihilation operator for the nth Bloch function of the single-particle Hamiltonian  ( ) k p . In case of a square lattice,  ( ) k p 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), . 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 nonsymmetric 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˜q e y 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 workq has been numerically checked to be always in the y-direction.
With figures 2(e), (f) one can also understand the fundamental differences between conventional spinimbalanced-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  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.
The situation is drastically different when the system enters to the large-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-q y region than for the small-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 +q y and some with -q y which results in an LO phase. Thus, the small-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 forq 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.
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.

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 n + E k, and quasi-holes n -E k, 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. 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.
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 n n c c c c k k k k k k k for six different values of μ, corresponding to six yellow dots depicted in figure 5(c). The topological transition corresponding to the gap closing at k 3 is studied in figures 5(d), (e), and correspondingly closings at k 2 and 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 k 3 (white line in figure 5(c)), the momentum distribution changes qualitatively in the vicinity of k 3 . This is further shown in figure 5(f) where n 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 k 2 (red line in figure 5(c)) is identified as an emergence of a prominent density peak around k 2 as clearly illustrated in figure 5(i). A similar peak can be also observed for the topological transition corresponding to 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.

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. In all three cases, D s diff more or less vanishes in large parts of the phase diagrams. However, especially when entering the large-q y FF region from the small-q y region, D s diff 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 s diff is positive and that the tFF −2 -phase in figure 6(c) near half-filling is clearly distinguishable from the neighboring phases. Therefore, by measuring D s diff 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].
In addition to D s diff , 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][13][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 quasiparticle energy at the same momentum k such that both are also the eigenvalues of  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 momentaq 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 s diff . 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 s diff than at smaller values of μ. In addition to different spatial components, one can also investigate the role of the geometric superfluid weight contribution D s geom which is presented for 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 s geom 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 all three cases the geometric contribution is smaller than the total superfluid weight and more or less vanishes when the system enters the large-q y FF regime.

Conclusions and outlook
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 spinflipping 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 atT t 0.1 . 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 t 0. 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 weakcoupling 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.
In the diamagnetic term there exists a double derivative  ¶ ¶ m n ( ) k which can be transformed to a single derivative via integrating by parts:

Appendix E. Direction of the Cooper pair momentum
In our computations the Cooper pair momentumq is in the y-direction, i.e. q 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 thatq is taken to be in the ydirection, in (b)q is along the diagonal direction ( =q q x y ) and in (c)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 + ˆq e e x y and q e y is plotted and correspondingly in figure E1(e) the thermodynamic potential difference of cases q e x and q e y is depicted. White lines show the phase boundaries between the BCS, FF and normal phases in case of q 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 q e y . Thus the states shown in figure E1(a) with q e y are energetically more stable than the states with the Cooper pair momentum in the diagonal or x-direction.
In figure E1 we have only presented three different options for the direction ofq 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 whenq 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˜||q e y . Therefore we are confident that our statements and results are correct within the mean-field theory framework.