Superfluidity enhanced by spin-flip tunnelling in the presence of a magnetic field

It is well-known that when the magnetic field is stronger than a critical value, the spin imbalance can break the Cooper pairs of electrons and hence hinder the superconductivity in a spin-singlet channel. In a bilayer system of ultra-cold Fermi gases, however, we demonstrate that the critical value of the magnetic field at zero temperature can be significantly increased by including a spin-flip tunnelling, which opens a gap in the spin-triplet channel near the Fermi surface and hence reduces the influence of the effective magnetic field on the superfluidity. The phase transition also changes from first order to second order when the tunnelling exceeds a critical value. Considering a realistic experiment, this mechanism can be implemented by applying an intralayer Raman coupling between the spin states with a phase difference between the two layers.

It is well-known that when the magnetic field is stronger than a critical value, the spin imbalance can break the Cooper pairs of electrons and hence hinder the superconductivity in a spin-singlet channel. In a bilayer system of ultra-cold Fermi gases, however, we demonstrate that the critical value of the magnetic field at zero temperature can be significantly increased by including a spin-flip tunnelling, which opens a gap in the spin-triplet channel near the Fermi surface and hence reduces the influence of the effective magnetic field on the superfluidity. The phase transition also changes from first order to second order when the tunnelling exceeds a critical value. Considering a realistic experiment, this mechanism can be implemented by applying an intralayer Raman coupling between the spin states with a phase difference between the two layers.
Magnetism is generally known to suppress superconductivity when the strength of the magnetic field exceeds a critical value. Survival of superfluidity in the presence of a strong magnetism has been a long-term interesting problem in the condensed matter physics. The central problem is that, in the Bardeen-Cooper-Schrieffer (BCS) theory of superconductivity, electrons form Cooper pairs in the spin singlet channel 1,2 . However, these pairs can be broken if the effective magnetic field is strong enough to flip the spin. This situation applies even if the Cooper pairs are mediated by magnetic fluctuations in some strongly correlated materials 3,4 . A possible exception is probably the theoretical prediction of a so called Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state [5][6][7] , where the Cooper pair has a finite center-of-mass momentum to form a spatially modulated order parameter [8][9][10][11] . Yet, the FFLO states have not yet been experimentally observed neither in condensed matter system 12,13 nor in the systems of ultracold atoms [14][15][16] . It is probably because the allowed parameter regime is in general too narrow to be observed. Another possible coexistence of magnetism and superconductivity arises in scenarios where the Cooper pairs become triplet states through the p-wave or f-wave interaction due to Pauli's exclusion principle [17][18][19][20][21][22][23] .
In this paper, we provide a new mechanism to greatly enhance superfluidity of ultracold Fermi gases in a bilayer system with a short range s-wave interaction within individual layers. The superfluidity then can survive in a much larger effective magnetic field even without going to the FFLO regime. This is possible by having a single particle spin-flip tunnelling between the layers. When the tunnelling amplitude exceeds a limiting value, the usual first order phase transition from the superfluid to normal state becomes second order, and the critical value of magnetic field increases almost proportionally to the tunnelling amplitude. Such a behavior can be understood from the fact that the spin-flip tunnelling couples atoms with two different spins in two different layers. This makes the Cooper pairs to include triplet contributions of spins in different layers to fulfil the Pauli exclusion principle. Similar results can be also observed in a multi-layer structure with a staggered effective magnetic field. Our results may be also relevant to the High T c superconducting material, where a strong anti-ferromagnetic correlation between nearest-neighboring CuO 2 planes is observed through the neutron-scattering experiment 24 .
In a realistic experiment, the outlined bi-layer scenario appears to be equivalent to a two-component (spinor) gas of ultracold atomic fermions loaded into a bi-layer trapping potential with a conventional tunnelling between the layers and a Zeeman magnetic field alternating in different layers, shown in Fig. 1. The alternating Zeeman field can be effectively generated by means of a Raman coupling [25][26][27] within individual layers with a properly chosen out of plane Raman recoil. The latter recoil provides the phase difference 2ϕ of the coupling amplitude in different layers needed for creating the alternating Zeeman field, as depicted in Fig. 1a. For ϕ = π/2 the scheme is mathematically equivalent to a setup involving a parallel Zeeman field and a spin-flip tunnelling (see Fig. 1b).

System and Methods
System Hamiltonian in original basis. We consider a spin-1/2 Fermi gas trapped in a bilayer potential.
In each layer the Raman beams induce spin-flip transitions with Rabi frequencies where the upper (lower) sign in ± corresponds to down (up) layer, with Ω x = Ω cos ϕ and Ω y = Ω sin ϕ. The phase difference for the Raman coupling in different layers 2ϕ ≡ |k R |d is achieved by taking a wave-vector of the Raman coupling k R perpendicular to the layers separated by a distance d (see Fig. 1a). The Pauli matrixes for the spin 1/2 atoms are denoted by σ x,y,z . On the other hand, it is convenient to treat the layer index as a pseudospin to be represented by the Pauli matrices τ x,y,z . As a result, the second quantized single particle Hamiltonian describing intralayer Raman transitions and interlayer tunneling can be written as x x x y z y k k k k 0 0 0 where σ 0 and τ 0 are identity matrixes, ε k = k 2 /2m − μ measures the kinetic energy with respect to the chemical potential μ, and t is the interlayer tunneling amplitude. The four component vector field operator , , , featured in Eq. (2) is a column matrix composed of operators ψ γ ∼ j k , annihilating an atom with a spin γ = ↑ , ↓ and a momentum k in a layer = j , , whereas Ψ ∼ † k is the corresponding raw matrix composed of the creation operators. For brevity in the following, we will omit the identity matrices τ 0 and σ 0 in tensor products like τ σ σ ⊗ ≡ The Hamiltonian (2) describes a quantum system of four combined layer-spin atomic states γ = =↑ ↓ j , ; , coupled in a cyclic way (see Fig. 1a): ↓ → ↑ → ↑ → ↓ → ↓ . The phase 2ϕ accumulated during such a cyclic transition allows to control the single particle spectrum 28 . The choice of the phase 2ϕ affects significantly also the many-body properties of the system, as we shall see later on.
We are considering a short range interaction between the atoms with opposite spins in the same layer. It is described by the following interaction Hamiltonian where g is the coupling strength. Note that H int has a symmetry group: 2  , where the two U(2)s describe the spin rotations in the first and second layer respectively, and Z 2 is the transpose transformation in pseudospin (layer) space, i.e., ψ ψ ↔ For ϕ = π/2, the Raman coupling can be represented by an effective Zeeman field antiparallel in each layer. This is mathematically equivalent to a parallel Zeeman field and a spin-flip tunnelling, as illustrated in a lower part of (b).
Scientific RepoRts | 6:33320 | DOI: 10.1038/srep33320 for the two layers. In order to have a better understanding of the following calculation results, it is convenient to represent the system in another basis. We first apply a unitary transformation rotating the spin σ around the z axis by the angle ϕ  for the up (down) layer. The resulting Zeeman field then becomes aligned along the x-axis in both layers. A subsequent spin rotation around the y axis by the angle − π/2 transforms σ x to σ z . After the two consecutive transformations the single particle Hamiltonian takes the form The transformed four component field operator is made of components ψ j↑,k and ψ j↓,k , which are superpositions of the original spin up and down field operators ψ ↑ ∼ j k , and ψ ↓ ∼ j k , belonging to the same layer = j , . Note that going to the new basis the spins are rotated differently in different layers.
The transformed single particle Hamiltonian (5) corresponds to a bilayer system subjected to a parallel Zeeman field along the z-axis for both layers, with the interlayer tunneling becoming spin-dependent for sin ϕ ≠ 0. For ϕ = π/2 the transformed Hamiltonian describes a completely spin-flip tunnelling, as illustrated in Fig. 1b. The interaction Hamiltonian given by Eq.
, which involves spin rotation within individual layers and thus does not change the form of H int .
Single-particle spectrum. The single-particle Hamiltonian H 0 given by Eq. (5) can be reduced to a diagonal form via a unitary transformation V ϕ for the field operator Ψ k , where c jγ,k is an annihilation operator for a normal mode characterized by the eigen-energy ξ γ , . In the following we shall concentrate on two specific cases of interest. (1) In the first case one has ϕ = 0, so that Ω y = 0 and Ω x = Ω. (2) In the second case the relative phase is ϕ = π/2, giving Ω x = 0 and Ω y = Ω. The first case corresponds to a spin-independent tunneling and non-staggered Zeeman field (in the transformed representation, Eq. (5)). The second case corresponds to the spin-flip tunneling and non-staggered Zeeman field along the z-axis, as shown in Fig. 1b. For these two cases the unitary transformation V ϕ diagonalizing the single particle Hamiltonian H 0 and the corresponding diagonal operator Ξ ϕ k of eigenenergies ξ γ )/ t . For ϕ = 0 the tunneling and Raman coupling are decoupled in the single particle Hamiltonian (5) or (2), so Ω and t are separable in single particle dispersion ξ γ j k , 0 . On the other hand, for ϕ = π/2 there is a term τ σ ⊗ y x in Eq. (5) which mixes the interlayer tunneling t and the Raman coupling Ω, so the single particle dispersion ξ γ π j k , /2 becomes non-separable. The latter case corresponds to a ring coupling scheme between four atomic states with an overall phase 2ϕ = π 28 . In such a situation the single particle eigenvalues ξ ε are twice degenerate with resect to the index = j , . This leads to significant differences in the BCS pairing for the two cases where ϕ = 0 and ϕ = π/2.
Without including the interaction effects, the chemical potential (Fermi energy) µ ϕ F satisfies where N is the total number of particles, A is the area of system and Θ is a unit step function. We use µ µ to represent the chemical potential for noninteracting particles at Ω = t = 0, with k F 0 being the corresponding Fermi momentum.
General framework in meanfield theory. In the present paper, we are interested in the effects due to attractive interaction between the atomic fermions (g < 0) in the bilayer system. As usual, a superfluid order parameter can be expected between fermions with opposite spins in the same layer, i.e.,  denotes the ground state expectation value. Without a loss of generality we can apply a U(1) transformation ψ jγ → e iα ψ jγ to make the order parameter complex conjugated in different layers ∆ = ∆ ≡ ∆ + ∆ ⁎ i R I . In general, there may be a phase difference between the order parameters in different layers described by Δ I . As it will be shown later, the imaginary part Δ I is zero for all cases to be considered.
Scientific RepoRts | 6:33320 | DOI: 10.1038/srep33320 Adopting the BCS mean-field approximation 29 , the interaction Hamiltonian (3) reduces to the following quadratic form of creation and annihilation field operators in the momentum space: Consequently we can express the total Hamiltonian, H = H 0 + H int , in the BCS form in terms of a set of normal single-particle operators C k and − † C k given by Eq. (7): y z y T R 0 I describes the mean-field atom-atom interaction responsible for the BCS pairing, and ϕ V T is a transposed diagonalization matrix (for a detailed derivation, see the Section I of the Supplementary Material). For ϕ = 0 and π/2, we have 1, 2, 3, 4) to be eigenvalues of the first term in Eq. (11), representing the Bogoliubov-DeGuinne (BdG) term, one arrives at the following total ground-state energy (see the Section II of the Supplementary Material), Finally, in a 2D Fermi gas, the fluctuation correction for the effective short-range interaction can be accounted by replacing , where b  is the two-body binding energy [30][31][32] . The superfluid gap equations are then determined by minimizing the total energy, i.e. ∂ ∂∆ = / 0 On the other hand, the equation , relates the atomic number N to the chemical potential μ. Our major aim is to study effects of the interlayer spin-flip tunneling (ϕ = π/2) on the superfluid properties for the bilayer system in the presence of magnetic field. We will not consider a possible FFLO phase that results from a mismatch in the Fermi energies for the two spins leading to the finite center-of-mass momentum for the Cooper pairs 5,6 . Including the spin mismatch should not substantially affect our major results, since the FFLO regions are in most cases too small to be observed 8,14,15,33,34 .
We note that the mean-field theory applied for the present 2D system is justified at zero temperature and in a weakly interacting regime considered here. Even at finite temperature of the 2D system, the critical temperature predicted by the mean-field theory is very close to the Kosterlist-Thouless transition temperature for a weakly interacting system 35 .

Results and Discussion
Single layer limit. To better understand results for our bilayer system, it is instructive to consider first a familiar single layer limit [30][31][32] corresponding to zero interlayer tunneling (t = 0) in the present model. This will allow one to see how the superfluidity is affected by the effective magnetic field provided by the Raman coupling Ω.
For t = 0 the positive eigenvalues of the BdG operator ε = +∆ ± Ω α ϕ E k k , 2 2 exhibit a two-fold degeneracy corresponding to different layers and are independent of ϕ as expected. The ground state energy Eq. (13) can then be calculated analytically to be (see the Section IIA of the Supplementary Material): is the ground-state energy for Ω = 0 30-32 . , however, results solely from the finite effective magnetic field, Ω, and comes into play only when Ω > Δ . Therefore the superfluidity can not be affected by a relatively small effective magnetic field field acting on the singlet Cooper pairs. In Fig. 2, we show how the ground state energy changes as a function of the order parameter, Δ , for various values of Ω. We take µ = .
0 01 b 0  in this and the subsequent calculations.
One can determine several important regimes for Ω where the superfluid order parameter, Δ , can be analytically determined by looking for the global minimum of   = ∆ ( ): Regime I corresponds to a limit of small Raman coupling (small effective magnetic field)  µ < Ω < 0 / 2 b . In this limit, we have . In other words, the last term of Eq. (14) is effectively zero and therefore the superfluid properties are identical to those for a usual 2D BCS state.
Regime II appears for µ µ < Ω < . In other words, the superfluid state starts to compete in energy with the normal state as the effective magnetic field is increased. Note that  µ In that case the last term of Eq. (14) is relevant. The obtained ground state corresponds to the normal state with Δ = 0. The superfluid state becomes then a meta-stable state with a finite stiffness.
Regime IV is reached for  µ Ω > 2 b . In that case the meta-stable superfluid state disappears, and therefore the system transforms to the completely normal state.
We note that the true first order phase transition occurs at the border between the Regimes II and III for  µ Ω = b . Yet the appearance of the meta-stable state in the Regimes II and III effectively broadens the phase transition making it not easily measurable. As we will see later, the inter-layer tunneling can completely change the situation. Zero Raman coupling limit. Next let us suppose there is a non-zero inter-layer tunneling t and no Raman coupling Ω = 0) [36][37][38] , so there is no effective magnetic field. In that case one arrives at spin-degenerate eigenvalues of the , which goes to the known single layer result, µ ∆ = 2 b  30 , in the zero tunneling limit (t → 0). Note that the single particle spectrum implies that both of the two bands are occupied. In the other limit, , only the states from the lower band could be occupied at zero temperature, and we have Figure 3 shows a behavior of order parameter as a function of the tunneling strength t for Ω = 0. Obviously, the superfluidity decreases with an increase of the tunneling strength, because the inter-layer tunneling plays a role of an effective Zeeman field in the pseudo-spin (layer) space. Yet now the order parameter decays in a power law in the limit of larger t, whereas in the previous case it goes abruptly to zero with increasing the Zeeman field Ω .
Raman coupling with ϕ = 0. Now let us consider a more general case with a finite Raman coupling and a finite interlayer tunneling for ϕ = 0. In such a situation eigenvalues of the BdG operator have no degeneracy,  Scientific RepoRts | 6:33320 | DOI: 10.1038/srep33320 of the momentum k, and goes to + ∞ in the limit of large k. If for some k the function becomes negative, it must cross the zero continuously. In such a situation the BdG spectrum will not open the gap at the Fermi surface.) This implies that |Δ R | should exceed Ω to have the superfluid phase. Since in that case ∑ α α ϕ= E k , 0 is independent of Ω, the superfuid ground energy would be the same as in the limit of zero Raman coupling. As in the previous cases, from the gap equations we have Δ I = 0 and thus Δ = |Δ R | for all Ω. To evaluate a possibility of a metastable state and a realistic border of phase transitions, we will consider analytical results for the effect of the Raman coupling Ω in two regimes.
0 Similar to the single layer limit (t = 0), one goes through four regimes with increasing the Raman coupling, Ω: , the ground state is superfluid with the order parameter being   , the ground state is still a superfluid phase with µ ∆ = 2 b t  , but the normal state becomes metastable. (III) When µ µ < Ω < 2 b t b t   , the ground state becomes a normal state with a metastable superfluid order parameter:  µ ∆ = 2 b t . (IV) When  µ Ω > 2 b t , the superfluid order disappears completely. The four regimes are shown in Fig. 4.  In the strong tunneling regime, t ≥ μ + Ω, the ground energy can be expressed to be (see the Section IIC of the Supplementary Material) 0   Similar to previous discussion, the four regimes as a function of Raman coupling, Ω can be also obtained analytically. Since this does not provide essentially new results, there is no need to present such analytic expression here. However, as one can see in the numerical phase diagram shown in Fig. 4, the regimes II and III are shrinking in the large tunneling limit, because the superfluid order parameter is also decreasing. In other words, for ϕ = 0, the ground state phase diagram is qualitatively similar to the single layer case (t = 0), because the inter-layer tunneling couples the two layers in the same way for both spin states (without a phase difference).
Raman coupling with ϕ = π/2. Now we consider the case where ϕ = π/2, with a finite Raman coupling Ω and interlayer tunneling t. In such a situation the tunneling involves a spin-flip (in the rotated basis). Eigenvalues of the BdG operator now are given by However, when the tunneling amplitude becomes larger, the range of superfluid phase increases significantly. This is very different from the phase diagram for ϕ = 0 shown in Fig. 4. Therefore a much stronger Raman coupling (effective magnetic field) is now required to destroy the superfluid phase (Regime I) which now goes directly to the normal phase (Regime IV) without passing the metastable phases (Regimes II and III). Such a phase transition is of the second order, a feature absent in the previously considered cases where the phase transition is of the first order. The first order phase transition now occurs only for small tunneling (t < 0.04 μ 0 ) where the superfluid state (Regime I) first goes to the metastable states (Regimes II and III) before reaching the normal state (Regime IV).
Note that the nature of the phase transition between Regime I (superfluid) and Regime IV (normal) is determined by the meta stable solutions in-between them, i.e., Regimes II and III, in the small t limit. When the interlayer tunneling is stronger than 0.04 μ 0 , the intermediate regime disappears because opposite spins in different layers are mixed due to the spin-flip inter-layer tunneling ( Fig. 1(b)), making the superfluid state in the s-wave pairing channel hardly to form in Regimes II and III. As a result, the phase transition for large t becomes fully determined by the curvature of free energy  ∂ ∂∆ ( / ) 2 2 at Δ = 0, the same as the condition determining the boundary between Regimes I and II.
This result is very unusual. Normally the effective magnetic field Ω and tunneling t reduce the superfluid properties by breaking the Cooper pairs through the Zeeman effects. As we can see from the single particle shows an even larger effective Zeeman field + Ω t 2 2 when including both effects, the Raman coupling and interlayer tunneling. The magnetism can be defined as = , Finally, we explore a situation where 0 < ϕ < π/2, so that both Ω x = Ω cos ϕ and Ω y = Ω sin ϕ are non-zero. In Fig. 7, we show the phase diagram for three different finite values of Ω x , i.e., for different magnitudes of the parallel Zeeman field. Although the increase in the superfluid pairing is still significant, the extent of the superfluid regime reduces for larger Ω x . In the rotated basis Ψ k , an increase of Ω x enhances the importance of the conventional tunneling with respect to the spin-flip tunneling. The conventional tunnelling determined by Ω x has a tendency to destroy the degenerate structure in the spectrum, so that it prevents formation of triplet pairing and reduces the superfluidity, unlike the spin-flip tunneling which is determined by Ω y . Furthermore, the phase boundary due to the co-existence of a meta-stable state becomes much broader, compared to the case with zero Ω x . In Fig. 8, we show the order parameter with respect to t for a fixed Ω y = 0.2 μ 0 . When t = 0 and Ω x = 0, the Cooper pair is a complete spin singlet, so the finite Raman coupling of Ω y prevents formation of Cooper pairs. Increasing the tunneling amplitude (t) enhances the triplet pairing. Thus for a sufficient large t, the system undergoes a transition to the superfluid from the normal state. For finite Ω x and Ω y = 0.2 μ 0 , there is a conventional tunneling in addition to the spin-flip tunneling in the rotated basis. In that case the superfluid formes in a narrow range of tunneling values t. Scientific RepoRts | 6:33320 | DOI: 10.1038/srep33320

Conclusions
We have explored a new mechanism to greatly enhance superfluidity of ultracold Fermi gases in a large range of the effective magnetic field. The mechanism can be implemented for a bilayer atomic system subjected to an interlayer tunneling. Additionally a Raman coupling induces intralayer spin-flip transitions with a phase difference between the two layers. Such a Raman coupling serves as a magnetic field staggered in different layers. After introducing a proper gauge transformation, one arrives at a non-staggered magnetic field and a spin-flip tunnelling between the layers. In such a situation the Cooper pairs were shown to acquire a component due to the triplet pairing. This supports a co-existence of the superfluidity for a much stronger effective magnetism. Our findings are helpful for understanding and controlling the superconductivity in the presence of the magnetic fields.