Interplay between orbital-quantization effects and the Fulde-Ferrell-Larkin-Ovchinnikov instability in multiple-band layered superconductors
Kok Wee Song Alexei E. Koshelev


	

	Materials Science Division,
	Argonne National Laboratory,
	9700 South Cass Avenue, Lemont, Illinois 60639, USA



	Materials Science Division,
	Argonne National Laboratory,
	9700 South Cass Avenue, Lemont, Illinois 60639, USA




	We explore superconducting instability for a clean two-band layered superconductor with deep and shallow bands in the magnetic field applied perpendicular to the layers. In the shallow band, the quasiclassical approximation is not applicable, and Landau quantization has to be accounted for exactly. 
	The electronic spectrum of this band in the magnetic field is composed of the one-dimensional Landau-level minibands. With increasing magnetic field the system experiences a series of Lifshitz transitions when the chemical potential enters and exits the minibands. These transitions profoundly influence the shape of the upper critical field at low temperatures.
	In addition, the Zeeman spin splitting may cause the nonuniform state with interlayer modulation of the superconducting order parameter [Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state]. Typically, the quantization effects in the shallow band strongly promote the formation of this state. The uniform state remains favorable only in the exceptional resonance cases when the spin-splitting energy exactly matches the Landau-level spacing. Furthermore, for specific relations between electronic spectrum parameters, the alternating FFLO state may realize, in which the order parameter changes sign between the neighboring layers. For all above cases, the reentrant high-field superconducting states may emerge at low temperatures if the shallow band has significant contribution to the Cooper pairing. 


Interplay between orbital-quantization effects and the Fulde-Ferrell-Larkin-Ovchinnikov instability in multiple-band layered superconductors
    Alexei E. Koshelev
    December 30, 2023
============================================================================================================================================





§ INTRODUCTION


The nature of superconducting instability in the magnetic field is a long-standing fundamental problem. In clean type-II superconductors, the upper critical field H_C2 and its temperature dependence are sensitive to the electronic band properties as well as to the gap structure which contains crucial information on how the correlated Cooper pairs are formed. The magnitude of H_C2 is mostly determined by the suppression of superconductivity due to the quasiparticles orbital motion and the Zeeman spin splitting of the Fermi surface by the magnetic field.
The problem of orbital H_C2 for a single-band system was solved in the seminal papers 
<cit.> within the quasiclassical approximation
which neglects the Landau quantization of the orbital motion. This approach works with high accuracy in common metals with large Fermi energies and describes very well most of known conventional superconductors.
Nevertheless, the exact Landau quantization calculation was discussed in Refs. Gunther:SSComm4.1966,*Gruenberg:PRev176.1968,Rajagopal:PLett23.1966 shortly after the quasiclassical work and the topic was further elaborated in great details later<cit.>.

Landau-quantization effects are most relevant when the Fermi energy is comparable to the cyclotron frequency so that only a few lowest Landau levels are occupied. Such extreme quantum limit has not been viable in most of known superconducting materials. The situation changed with the recent discovery of the iron-based superconductors (FeSC). It is very common in FeSCs that some of the bands have very small Fermi energies and they can be driven through the Lifshitz transition by chemical dopings<cit.> and, importantly, these materials usually remain superconducting after one of the bands is completely depleted.
As the Fermi energy of the shallow band may be smaller than the pairing energy scale, the proper treatment of this band requires revision of the conventional BCS approach. Such situation was first considered in the context of shape resonances in confined superconductors (thin films and nanowires) <cit.>
The influence of the shallow band on the superconducting transition temperature  in multiple-band materials near Lifshitz transition has been considered in Ref. <cit.> and, more recently, 
physics of FeSCs motivated detailed investigations of this problem   
<cit.>. 
It was demonstrated that, in spite of low quasiparticle density, shallow bands may strongly influence the Cooper pairing. 
Another key property facilitating the extreme quantum limit near superconducting instability of FeSCs is that they are characterized by the very high upper critical fields, up to 100 T.

As the orbital upper critical field is inversely proportional to the Fermi energy, it is strongly affected by the shallow bands. Moreover, in the vicinity of the Lifshitz transition, the extreme quantum limit can be reached for these bands near H_C2 meaning that the quasiclassical consideration cannot be used. 
It has been indeed demonstrated that the Landau quantization has a dramatic effect on superconducting instability in the magnetic field for two-dimensional materials near the Lifshitz transition<cit.>. The most spectacular prediction is the emergence of the pronounced reentrant states at the magnetic field corresponding to the matching of the Landau levels with the chemical potential.

In addition to the orbital effect, shallow bands also promote the Zeeman spin-splitting suppression of superconductivity. The relative role of the spin and orbital mechanisms is usually characterized by the Maki parameter α_M defined as  √(2)H_C2^O/H_C2^P, where H_C2^O and H_C2^P are the orbital and spin upper critical fields, respectively.  In a clean isotropic single-band material with the free-electron g factor, the Maki parameter is proportional to the ratio of the superconducting gap Δ and the Fermi energy ϵ_F, α_M=π^2Δ/4ϵ_F.
In this case, unless the band is very shallow, α_M is very small, meaning that the orbital effect strongly dominates. The Maki parameter is greatly enhanced in special cases of weak orbital effect such as quasi-one-dimensional materials and layered superconductors for magnetic field directed along the layers.

An important consequence of the strong Zeeman effect (i.e., large α_M) is the emergence of the nonuniform Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state in very clean materials <cit.>. In this state, the order parameter is periodically modulated which allows for a gain in the Zeeman energy which may exceed the kinetic-energy loss due to the nonzero center-of-mass momentum pairing. Rich physics of the FFLO state has been extensively investigated in many subsequent theoretical studies within quasiclassical approximation <cit.>. 
The FFLO instability may only realize in very pure materials with weak scattering of quasiparticles.
That is why,
even though the FFLO state was predicted more than half century ago, the experimental data consistent with this state have been reported only relatively recently in organic superconductors <cit.> and, less convincingly, in the heavy-fermion compounds <cit.>.
The conditions for the FFLO instability with modulation along the magnetic field in isotropic material were studied by Gruenberg and Gunther <cit.> within the quasiclassical framework. They found that the nonuniform state appears only for very large Maki parameters, α_M> 1.8, which is unlikely to realize in any single-band isotropic material. That is why most FFLO studies have been focused on quasi-low-dimensional materials with strongly reduced orbital effects. 

The FFLO state may realize in the iron-based superconductors due to the presence of shallow bands and huge upper critical fields which are likely to be limited by the Zeeman effect. 
This motivated recent investigations of the conditions for the emergence of this state in multiple-band materials within the quasiclassical approach in different situations <cit.>. However, in the presence of very shallow bands, this approach may be insufficient and quantization effects have to be accounted for.

In this paper, we investigate the superconducting instability for a two-band layered material in the magnetic field applied perpendicular to the layers. We consider the case of the material near the Lifshitz transition when one of the bands is very shallow. In this case, the Landau-quantization effects strongly influence the formation of the superconducting state. Even though our consideration is motivated by physics of iron-based superconductors, it is very general, and our goal is not to describe  any particular compound but, instead, to develop a general understanding of how the orbital-quantization effects influence the onset of superconductivity in clean two-band layered materials.  The major new feature in comparison with a pure two-dimensional case <cit.> is that, due to the large Zeeman effect in the shallow band, this system is prone to the formation of the FFLO state with interlayer modulation of the order parameter. The quantum effects have a profound influence on this FFLO instability.

The interlayer tunneling lifts all the degenerate Landau levels to dispersive minibands along the out-of-plane momentum direction. As a consequence, the system experiences series of Lifshitz transitions with increasing the magnetic field corresponding to crossing of the chemical potential with the miniband edges. Every Landau-level miniband has two van Hove singularities at the reduced z-axis momentums k_z=0 and π, at which the density of states (DOS) is enhanced. In special situations, when two such singular points for spin-up and spin-down bands simultaneously match the Fermi level, the pairing strongly enhanced. Two distinct cases of such resonance matching are possible. The first well-known case is realized when spin-splitting energy is equal to the Landau-level spacing <cit.>. In this case, the van Hove points of the same kind (either k_z=0 or π) may match leading to the standard uniform superconducting state. The second case corresponds to matching of the opposite van Hove points, e. g., spin-down/k_z=0 and spin-up/k_z=π points, which may occur only for a certain relation between the electronic band parameter. In this situation, the alternating FFLO state may emerge, in which the superconducting order parameter changes sign between the neighboring layers.  These matching effects may generate the high-field reentrant superconducting states which are somewhat less pronounced than for two-dimensional case <cit.> due to the DOS spreading by the interlayer tunneling.

On the other hand, the Landau-level spreading somewhat mitigates the Zeeman pair-breaking effect, 
since the dispersive spin-up and -down minibands can cross the Fermi level simultaneously within a finite energy range for arbitrary spin splitting. In such generic situation, the shallow band favors the formation of the FFLO state, in which the optimal modulation wave vector equals to the difference between the spin-up and -down Fermi momenta. Such FFLO instability leads to different kind of reentrant states with the field-dependent modulation wave vectors. In contrast to the reentrant states caused by the matching of the van Hove singularities, the latter states can extend over a broad magnetic-field range.

This paper is organized as follows. In Sec. <ref>, we introduce the microscopic Hamiltonian describing a two-band layered superconductor the in magnetic field and derive the corresponding linearized gap equations. In Sec. <ref>, we briefly discuss these equations for the transition temperature in zero magnetic field providing the reference for the further investigation of instability in finite magnetic field. 
In Sec. <ref>, we derive the equations for the upper critical field H_C2 and evaluate the pairing kernels in these equations. 
In Sec. <ref>, we discuss the dependences of the quantum pairing kernel on relevant parameters.
In Sec. <ref>, we present the typical magnetic field versus temperature phase diagrams. We conclude the paper in Sec. <ref>.



§ THE MODEL OF A TWO-BAND LAYERED SUPERCONDUCTOR


We will investigate the shallow-band effects in layered superconductors using the simple tight-binding Hamiltonian with only the nearest-neighbor interlayer hopping term,

    ℋ=∑_ jα∫d^2𝐫[c^†_α j s(𝐫)(ε^α(𝐤̂) σ^0_ss'-μ_zHσ^z_ss')c_α js'(𝐫)
       -t_z^αexp[ie/c∫_l_jdzA_z]c^†_α j s(𝐫)c_α,j+1,s(𝐫)+H.c.
       -∑_βU_αβc^†_α j↓(𝐫)c^†_α j↑(𝐫)
    				c_β j↓(𝐫)c_β j↑(𝐫)],

where 𝐫=(x, y) is the in-plane coordinate, j is the layer index, s represents spin, and α=e (h) represents the e-band (h-band). Furthermore, t^α_z is the (nearest) interlayer hopping energy, ε^e(𝐤̂)=𝐤̂^2/(2m_e) and ε^h(𝐤̂)=-𝐤̂^2/(2m_h)+ε_0 are intralayer energy dispersions with the band masses m_h and m_e and the momentum operator 𝐤̂=-i∇_𝐫-e𝐀/c. In the model, we also consider the Zeeman spin splitting and assume, for simplicity, that the band electron's magnetic moment, μ_z, is the same in each band. 
In the second line of Eq. (<ref>), ∫_l_j=∫^(j+1)a_z_j a_z is the line integration along the out-of-plane direction and a_z is the interlayer spacing (set to unity in the later calculation) and A_z is the z-component of the vector potential. In zero magnetic field, the three-dimensional energy dispersions are ϵ^α(𝐤,k_z)=ε^α(𝐤)-2t^α_zcos k_z, see Fig. <ref>(a).


In this paper, we limit ourselves to the case of the magnetic field applied in the z direction. In this case, A_z=0 with standard gauge choices and the magnetic phase factor in Eq. (<ref>) tight-binding terms drops out. For definiteness, we consider the case when the hole band is shallow meaning that the chemical potential μ is located near its edge μ_h=ε_0-μ≪ε_0. We note that whether the shallow band is hole-like or electron-like 
does not have any influence on the results of this paper.
The two Lifshitz transitions in this model occur at μ_h=± 2t^h_z [see Fig. <ref>(b)] <cit.>. At μ_h= 2t^h_z the neck near k_z=0 interrupts and μ_h= -2t^h_z the hole band got completely depleted. 

To study the superconducting instabilities in the model Eq. (<ref>), we follow the standard approach and write the linearized gap equation as

    Δ^α_𝐫, j=T∑_ω_n∑_β j' U_αβ∫_𝐫'K^β_ω_n(𝐫 j,𝐫' j')Δ^β_𝐫',j',

where the gap function is defined as Δ^α_𝐫, j=∑_β U_αβ⟨ c_β, j↓(𝐫)c_β, j↑(𝐫)⟩ and we used the notation ∫_𝐫=∫d𝐫, ω_n=2π T(n+1/2) are the Matsubara frequencies, and the kernel is 

    K^α_ω_n(𝐫 j,𝐫' j')
    =G^α_ω_n,+(𝐫 j,𝐫' j')G^α_-ω_n,-(𝐫' j',𝐫 j),

where G^α_ω_n,±(𝐫 j,𝐫' j') is the one-particle Green's function in the normal phase and the subscripts + or - describe spin orientation. Without magnetic field, the normal Green's functions can be expanded into the plane-wave basis 

    G^α_ω_n,±(𝐫 j,𝐫' j')=∑_𝐤k_ze^-i[𝐤·ρ+
    k_z( j- j')]/iω_n-ξ^α_𝐤+2t_z^αcos k_z,

where ξ^α_𝐤=ε^α(𝐤)-μ and ρ=𝐫-𝐫'.
In the presence of out-of-plane magnetic fields,  in the symmetric gauge 𝐀=H/2(-y,x,0) 
the one-particle Green's function can be represented as

    G^α_ω_n,±(𝐫 j,𝐫' j')=   exp(i[𝐫×𝐫']_z/2l^2)
    ×   ∑_k_ze^-ik_z( j- j')g^α_ω_n,±(ρ,k_z)

where l=√(c/(eH)) is the magnetic length and  ρ=|ρ|. We use the quasiclassical approximation for the Green's function of the deep e band 

    g^e_ω_n,±≈∑_𝐤e^-i𝐤·ρ/iω_n∓μ_zH-ξ^e_𝐤+2t^e_zcos k_z,

and expand the Green's function of the shallow h band over the exact Landau-level basis <cit.>

    g^h_ω_n,±=1/2π l^2∑_ℓ=0L_ ℓ(ρ^2/2l^2)exp(-ρ^2/4l^2)/iω_n∓μ_zH+ℰ_ℓ-μ_h+2t^h_zcos k_z,

where μ_h=ε_0-μ, ℰ_ℓ=ω_c( ℓ+1/2), ω_c=eH/m_hc is the cyclotron frequency, and L_ℓ(x) are the Laguerre polynomials.

We note that in the single-band case the chemical potential of the shallow band may depend on the magnetic field <cit.>.  However, in our case the deep band serves as the particle reservoir that stabilizes the chemical potential. This allows us to solve the gap function assuming fixed chemical potential (see also Ref. <cit.>).



§ TRANSITION TEMPERATURE IN ZERO MAGNETIC FIELD (T_C)


For zero magnetic field, the superconducting order parameter is homogeneous, Δ^α_𝐫 j=Δ^α_0.
This gives the linearized gap equation, which we represent as 

    Δ^α_0=∑_βΛ_αβΛ^-1_0,βΔ^β_0,

where we have introduced the notations for the coupling matrix
Λ_αβ=U_αβN_β with N_α=m_α/(2π) is the densities of states per layer (for shallow band, N_h is the true density of states only in the regime of open Fermi surface, μ_h>2t_z^h), and

    Λ^-1_0,α≡ N^-1_α T∑^∞_ω_n=-∞∑_j'∫_𝐫'K^α_ω_n(𝐫 j,𝐫' j').
 
The kernels are determined by the zero-field Green's function in Eq. (<ref>) giving

    Λ^-1_0,α=N^-1_α T∑_ω_n∑_𝐤k_z[ω^2_n+(ξ^α_𝐤-2t_zcos k_z)^2 ]^-1,

(see Appendix <ref>).
We further integrate out the wave vectors in Λ^-1_0,α by using ∑_𝐤=N_e∫^Ω_-Ωdξ^e for e band (∑_𝐤=N_h∫^μ_h_-Ωdξ^h for h band) and ∑_k_z=∫^π_-πdk_z/2π. Here, Ω is the energy cutoff and we have assumed that Ω≫ t_z, T_C,μ_h. We therefore obtain the following gap equation near T_C:

    [Λ̂^-1-Λ̂^-1_0][ Δ^h_0; Δ^e_0 ]=0,

where Λ̂^-1_0=diag[Λ^-1_0,h,Λ^-1_0,e] with


    Λ^-1_0,e   ≈∑^∞_ω_n>02 T_C/ω_ntan^-1Ω/ω_n≈ln𝐴Ω/ T_C,
    Λ^-1_0,h   ≈1/2ln𝐴Ω/T_C+Υ_C,


𝐴=2e^γ_E/π≈ 1.134, and γ_E≈ 0.5772 is the Euler-Mascheroni constant.
Here, the parameter Υ_C=Υ_T_C is the value of the temperature-dependent function

    Υ_T=∑^∞_ω_n>0∫^π_-πdk_z/2π2 T/ω_ntan^-1[μ_h-2t_z^hcos k_z/ω_n]

at T=T_C. This parameter appears due to the cut off at the band edge for
the shallow band. We can carry out k_z integration and the Matsubara-frequency
sum in Υ_T using the relation
tan^-11/a=-∫^∞_0dx/xe^-axsin x. This
gives

    Υ_T
    =-∫^∞_0ds/π slntanh(π T s)sin(2μ_h s)J_0(4t_z^hs),

where J_0(x) is the Bessel's function. 
In the limit of low temperatures, T≪μ_h, this function diverges logarithmically, 
Υ_T≈12ln[𝐴(μ_h+√(μ_h^2-(2t_z^h)^2))/(2 T)].

The transition temperature  T_C is determined by the condition of degeneracy of the matrix Ŵ≡Λ̂^-1-Λ̂^-1_0, i. e., [Ŵ]=0. This equation determines Λ_0,e as

    Λ _0,e^-1   =Λ _ee+Λ _hh/2/𝒟_Λ-Υ _C
    
    +   δ _Λ√(( Λ _ee-Λ _hh/2/𝒟_Λ-Υ _C) ^2+2Λ _ehΛ _he/𝒟_Λ^2)

with δ_Λ=-sign[ (1-Υ_C
Λ_hh) 𝒟_Λ] and 𝒟_Λ=Λ _eeΛ _hh-Λ _ehΛ _he. This parameter
directly determines  T_C by Eq. (<ref>). In the following
consideration, we will use the zero-field equation, Eq.  (<ref>), to
eliminate the logarithmic divergences from the gap equation at finite magnetic
field.



§ SUPERCONDUCTING INSTABILITY IN MAGNETIC FIELD


In the presence of magnetic field, the superconducting order parameter is nonuniform. Near the onset of the superconducting instabilities in the magnetic field, the solutions for the linear problem in Eq. (<ref>) are given by the Landau-level eigenfunctions. Typically, the lowest Landau-level eigenfunction yields the optimal H_C2, where the superconducting instability develops first
[For strong Zeeman spin splitting, the largest H_C2 may be realized for higher Landau-level eigenfunction <cit.>. 
]. 
This solution has the following shape<cit.>

    Δ^α_𝐫, j=Δ^α_0exp(-r^2/2l^2+iQ_z j).

We have assumed the possibility of the layer-to-layer modulation in the order parameters with the wave vector Q_z. 
[To be accurate, the complex order parameter in Eq. (<ref>) was suggested by Fulde and Ferrel <cit.>. As demonstrated by Larkin and Ovchinnikov<cit.>, for purely paramagnetic case the ground state below the transition is actually given by the order parameter with the amplitude modulation ∝cos(Q_z j).  
In the presence of orbital effects, both states may realize <cit.>.
We only investigate the instability location here, which is the same for both states.] 
Such modulation is the realization of the nonuniform FFLO state <cit.>.
Superconducting instability is determined by the condition that Eq. (<ref>) provides a solution of the gap equation at least for one value of Q_z and then this optimal modulation wave vector realizes in the emerging superconducting state. 
The ansatz in Eq. (<ref>) is the kernel eigenfunction in the magnetic field,

    ∑_ j'∫_𝐫'K^α_ω_n(𝐫 j,𝐫' j')Δ^α_𝐫' j'=π N_αλ^α_ω_n,Q_zΔ^α_𝐫, j

with

    λ^α_ω_n,Q_z=2/N_α∫^∞_0ρdρexp(-ρ^2/2l^2) 
       ×⟨ g^α_ω_n,+(ρ,k_z-Q_z2)
    g^α_-ω_n,-(ρ,k_z+Q_z2)⟩_z.

Here and below we use notation ⟨ F(k_z)⟩_z≡∫^π_-πdk_z/2π F(k_z) for the averaging with respect to k_z. We omit the dependence of λ^α_ω_n,Q_z on T, H, and the electronic parameters.
As follows from the definitions of the Green's functions, Eqs. (<ref>) and (<ref>), (λ^α_ω_n,Q_z)^∗=λ^α_-ω_n,Q_z. 
Next, we first discuss the kernel eigenvalue for the deep e-band and then for the shallow h-band. 

For the deep e-band, the Landau quantization does not play a role. Using Eqs. (<ref>) and (<ref>), 
linearizing the band dispersion at Fermi level [ε^e_𝐤_F+𝐪≈ε^e_𝐤_F+v_e·𝐪, where v_e=∇_𝐤ε^e_𝐤|_𝐤^F=𝐤^F/m_e with the Fermi wave vector 𝐤^F=(k^F_x, k^F_y)], 
we obtain the following quasiclassical result for the kernel eigenvalue

    λ_ω_n,Q_z^e=2∫_0^∞ds ⟨exp[-s^2(μ+2t_z^ecos k_zcosQ_z/2)/m_el^2
       -2sζ_ω(ω_n+iμ_zH-2it_z^esin k_zsinQ_z/2)]⟩_z,

where ζ_ω≡sign(ω_n). The derivation details are described in the Appendix <ref>.

For the shallow h band, substituting  the Green's function, Eq. (<ref>), into the general presentation in Eq. (<ref>), and 
using relation ∫^∞_0 d x L_ℓ(x)L_ℓ'(x)e^-2x=(ℓ+ℓ')!/2^ℓ+ℓ'+1ℓ!ℓ'!,
we derive

    λ^h_ω_n,Q_z= -1/2πω_c
    ×   ∑_ℓℓ'⟨(ℓ+ℓ')!/(2^ℓ+ℓ'ℓ!ℓ'!)/(iω̅_n+ℓ+1/2-γ̃_z+μ̃_h)
    	(iω̅_n-ℓ'-1/2-γ̃_z+μ̃_h)⟩_z,

where we introduced the following notations


    μ̃_h(k_z,Q_z)   =μ̅_h-2t̅_z^hcos k_zcos (Q_z/2),
    γ̃_z(k_z,Q_z)   =γ_z-2t̅_z^hsin k_zsin (Q_z/2).


Here all “barred” normalized quantities are defined as a̅≡ a/ω_c (with a=ω_n,μ_h,t_z). Furthermore,
γ_z=μ_zm_hc/e=gm_h/4m_0 is the reduced spin-splitting parameter, where m_0 is the free-electron mass and g is the band-electrons g-factor. 

Therefore, the gap equation (<ref>), in the presence of the magnetic
field, becomes

    Λ̂^-1[ Δ^h_0; Δ^e_0 ]
    =2π T Re∑^Ω_ω_n>0[ λ^h_ω_n,Q_zΔ^h_0; λ^e_ω_n,Q_zΔ^e_0 ].

The Matsubara-frequency sums, 𝒮_α(H,T,Q_z)≡2π T∑^Ω_ω_n>0λ^α_ω_n,Q_z, are logarithmically divergent and have to be cut at ω_n∼Ω. Similarly to the two-dimensional case <cit.>, we can regularize the gap equation using its zero-field counterpart, Eq. (<ref>). Namely, we decompose 𝒮_α(H,T,Q_z) as 𝒮_α(H,T,Q_z)=[𝒮_α(H,T,Q_z)-𝒮_α(0,T,0)]+[𝒮_α(0,T,0)-𝒮_α(0,T_C,0)]
+𝒮_α(0,T_C,0) with 𝒮_α(0,T_C,0)≡Λ_0,α^-1 being the only log-diverging term. 
This leads to

    Ŵ[ Δ^h_0; Δ^e_0 ]
    +[ 𝒜_1Δ^h_0; 𝒜_2Δ^e_0 ]
    =[ 𝒥_1Δ^h_0; 𝒥_2Δ^e_0 ].

Here, 𝒜_α(T)≡𝒮_α(0,T_C,0)-𝒮_α(0,T,0)
are the temperature-dependent parts of the pairing eigenvalues,
𝒜_2(T)=ln(T/T_C) and  𝒜_1(T)=12ln(T/T_C)-Υ_T+Υ_C,
where the function Υ_T(T/μ_h, t_z^h/μ_h) is defined in Eq. (<ref>). 
The field-dependent parts, 
𝒥_α (H,T,Q_z)≡𝒮_α(H,T,Q_z)-𝒮_α(0,T,0), are


    𝒥_1(H,T,Q_z)=2π T∑^∞_ω_n>0Re(λ^h_ω_n,Q_z-1/2ω_n)-Υ_T,
       𝒥_2(H,T,Q_z)=2π T∑^∞_ω_n>0Re(λ^e_ω_n,Q_z-1/ω_n).


Note that, by definition, 𝒜_α(T_C)=0 and 𝒥_α(0,T,0)=0.
Therefore, the UV cutoffs are explicitly removed and the Matsubara-frequency sums in the right-hand side in Eqs. (<ref>) and (<ref>) converge now in the limit of Ω→∞. 
We can represent the functions 𝒥_α in this limit as  


	
    𝒥_1(H,T,Q_z)    =1/4∑_m=0^∞∑_ℓ=0^mm!/2^m(m-ℓ)!ℓ!⟨𝒯(ℓ+γ̃_z-μ̃_h)+𝒯(m-ℓ-γ̃_z-μ̃_h)-2𝒯(m/2-μ̃_h)/m+1-2μ̃_h⟩ _z
        -1/2⟨∫ _0^1/2dx𝒯(x-1/2-μ̃_h0)/x-2μ̃_h0+∑_m=0^∞∫ _-1/2^1/2dx[𝒯(m+x/2-μ̃_h0)/m+1+x-2μ̃_h0-𝒯(m/2-μ̃_h)/m+1-2μ̃_h]⟩ _z,
    𝒥_2(H,T,Q_z)=     2∫_0^∞dslntanh(π T/ω_c^es)⟨exp(-μ̃_es^2)[μ̃_escos(2γ̃_z^es)+γ̃_z^esin(2γ̃_z^es)]⟩ _z,



where we introduced notations 𝒯(x)≡tanh[ω_c(x+12)/2T], μ̃_h0(k_z) =μ̅_h-2t̅_z^hcos k_z. Furthermore,   ω_c^e=eH/cm_e is the e-band cyclotron frequency, μ̃_e≡μ̅+2t̅_z^ecos k_zcosQ_z/2,  γ̃_z^e≡γ_z-2t̅_z^esin k_zsinQ_z/2, μ̅=μ/ω_c^e, and t̅_z^e=t_z^e/ω_c^e.
We remind that the parameters μ̃_h and γ̃_z also depend on k_z [see Eqs. (<ref>) and (<ref>)]. We describe derivation of 𝒥_1(H,T,Q_z) in Appendix <ref>.  The double sum in Eq. (<ref>) collects the contributions to pairing coming from the quasiparticles located at the Landau levels ℓ and m-ℓ. The quantum kernel eigenvalue 𝒥_1 depends on five independent dimensionless parameters: the reduced magnetic field ω_c/μ_h, the reduced temperature T/μ_h, the modulation wave vector Q_z, the ratio t_z^h/μ_h, and the spin-splitting factor γ_z.

The H_C2 problem is reduced to the solution of equation

    [ W_11+𝒜_1-𝒥_1         W_12;         W_21 W_22+𝒜_2-𝒥_2 ]=0

for a given T and Q_z, and we need to find the optimized Q_z in H_C2(T,Q_z) 
for which the instability develops first.
As the matrix Ŵ is degenerate, this equation can be rewritten as

    ∏_α=1,2( 1+𝒜_α(T)-𝒥_α(H,T,Q_z)/W_αα) =1.

This is our main equation for determination of the upper critical field in two-band layered superconductors. 
All information about the coupling matrix is contained in the two parameters, W_11 and W_22.
The analytical expressions for these parameters can be derived from Eq. (<ref>) (see also Ref. <cit.>),


    W_11   =Λ_ee-1/2Λ_hh/2𝒟_Λ-Υ_C/2+δ_WR/2
    
    W_22   =-Λ_ee-1/2Λ_hh/𝒟_Λ-Υ_C+δ_WR


with 𝒟_Λ=Λ_eeΛ_hh-Λ_heΛ_eh, δ_W=sign[𝒟_Λ(1-Υ_CΛ_hh)], and

    R=√((Λ_ee-1/2Λ_hh/𝒟_Λ-Υ_C)^2+2Λ_heΛ_eh/𝒟_Λ^2).

As follows from Eq. (<ref>) the weights with which the bands contribute to the pairing near the upper critical field scale as 1/|W_αα|. 

The behavior of the upper critical field is mostly determined by the shape of the kernel
eigenvalues 𝒥_α(H,T,Q_z). In general, larger values of 𝒥_α correspond to stronger pairing strength. 
In the next section, we will explore in detail the quantum kernel eigenvalue
𝒥_1. 



§ BEHAVIOR OF  THE QUANTUM PAIRING KERNEL EIGENVALUE
 



The shape of magnetic field–temperature phase diagrams is determined by the behavior of the pairing kernel eigenvalues 𝒥_α. While the quasiclassical kernel eigenvalue 𝒥_2 is well studied and has monotonic dependence on the magnetic field for all temperatures, the quantum kernel eigenvalue 𝒥_1 has rather complicated nonmonotonic dependence on the magnetic field which is very sensitive to the electronic-spectrum parameters μ_h and t_z^h as well as the spin-splitting factor γ_z. 
The electronic spectrum of the shallow band in the magnetic field is composed of the Landau-level minibands with width 4t_z^h [see Fig. <ref>(a)]. The system has series of Lifshitz transitions with increasing magnetic field when the chemical potential enters or exits a particular miniband. 
At low temperatures, the magnetic-field dependences of 𝒥_1 have features at these transitions whenever the chemical potential crosses the van Hove points at the miniband edges. There are two such points for every miniband corresponding to two values of k_z, 0 and π. In addition, there are two minibands for every Landau level for two spin orientations.
This gives four miniband-edge magnetic fields per Landau level, H_ℓ, σ,δ_t=(cm_h/e)ω_ℓ, σ,δ_t, corresponding to cyclotron frequencies

    ω_ℓ, σ,δ_t
    =μ_h+2δ_t t_h^h/ℓ+12+σγ_z,
   
where σ=± 1 (↑/↓) describes spin orientation and δ_t=1 and -1 corresponds to k_z=π and 0, respectively. 
In addition, we have to consider the behavior of 𝒥_1  for different modulation wave vectors Q_z. In general, as the larger 𝒥_1 corresponds to stronger pairing strength, the shallow band favors the states which maximize 𝒥_1.

Before proceeding to the investigation of the magnetic field–temperature phase diagrams for two-band systems, it is very instructive to study the analytical properties of the function 𝒥_1(H,T,Q_z), particularly, in the low-temperature limit. We first focus on its singularities when the magnetic field crosses the typical values in Eq. (<ref>) and on
identifying the possible divergences for T→0, since these features have important implications on the superconducting instabilities. Next, we consider the magnetic field dependences of 𝒥_1 for representative cases. Furthermore, by studying the Q_z dependences of 𝒥_1 for different Zeeman spin-splitting parameters, we identify possible parameter ranges for FFLO instabilities. 



 §.§ Low-temperature limit and its divergences



In this subsection, we investigate the leading divergences of the quantum kernel eigenvalues 𝒥_α as T→ 0 for different cases.
We review first the behavior of the quasiclassical kernel eigenvalue 𝒥_2 [Eq. (<ref>)]. For all magnetic fields, 𝒥_2 has the same logarithmic divergence as 𝒜_2 so that 𝒥_2-𝒜_2 approaches a finite value. In the typical situation of a moderate spin-splitting factor, γ_z≪√(μ/ω_c^e), it can be treated perturbatively. For the uniform case, Q_z=0, we derive from Eq. (<ref>) the result for the 
zero-temperature limit of the full quasiclassical kernel eigenvalue

    𝒥_2(H,T→0,0 ) -𝒜_2(T→0 )
    =-1/2ln r_C^(0),
       r_C^(0)≈e^γ_Eω_c^e/π^2T_C^2μ+√(μ^2-(2t_z^e)^2)/2(1+2γ_z^2ω_c^e/√(μ^2-(2t_z^e)^2)).

The parameter r_C^(0) is just the ratio H/H_c2^e(0), where H_c2^e(0) is the upper critical field of the deep band at zero temperature.
In the case of finite Q_z, the parameter r_C^(0) has to be replaced by the function r_C(Q_z). The closed analytical result for  r_C(Q_z) is not available even for γ_z=0. One can only derive an approximate result in the limits,  γ_z≪√(μ/ω_c^e) and t_z^e/μ, t_z^esin(Q_z/2)/√(ω_c^eμ)≪ 1:

    r_C(Q_z)≈   e^γ_Eω_c^eμ/π^2T_C^2[1+2ω_c^eγ_z^2/μ+4(t_z^e)^2/μω_c^esin^2Q_z/2
       -(t_z^e)^2/μ^2cos^2Q_z/2].

Generally, the quasiclassical kernel eigenvalue 𝒥_2-𝒜_2 is a monotonically decreasing function of H at all temperatures and spin-splitting parameters  meaning that the magnetic field always suppresses  superconductivity. As expected, it has maximum at Q_z=0 in the limit  γ_z≪√(μ/ω_c^e) meaning that the deep band favors the uniform state. 

The quantum kernel eigenvalue 𝒥_1 typically behaves similarly to 𝒥_2, i. e., it  has the same logarithmic divergence as 𝒜_1, 𝒥_1 ∝ln (T/T_C)  so that the total kernel eigenvalue 𝒥_1-𝒜_1 approaches a finite value in the zero-temperature limit.
This zero-temperature value, however, has singular contributions when the magnetic field crosses the typical values given by Eq. (<ref>), corresponding to Lifshitz transitions for the Landau-level minibands. We discuss these singularities in the next subsection.  
In several exceptional resonant cases, when two miniband-edge fields with opposite spin orientations are identical, the Landau quantization leads to faster divergence 𝒥_1∝ 1/√(T).  In the case when these two fields originate from  the van Hove singular points of the same type (either k_z=0 or π), the divergence occurs in the uniform channel Q_z=0. On the other hand, if the two fields correspond to the opposite  van Hove points, the divergence takes place in the alternating channel Q_z=π. We discuss both these cases below.
Another divergence appears within the field ranges where the chemical potential simultaneously crosses two minibands with opposite spin orientations. In this case, the total kernel eigenvalue 𝒥_1-𝒜_1 diverges logarithmically for T→ 0 at the optimal wave vector Q_z=Q_op connecting  the minibands' Fermi momenta [see Fig. <ref>(b)].
All these low-T divergences may lead to the high-field superconducting states.



  §.§.§ Square-root singularity of the pairing kernel at the miniband-edge fields 


As discussed above, the system has series of the miniband Lifshitz transitions at the magnetic fields given by Eq. (<ref>). In this subsection we discuss singularity of the pairing kernel eigenvalue 𝒥_1(H,T→0,Q_z) at these transitions, when the cyclotron frequency ω_c crosses the miniband-edge value ω_ℓ_0,σ,δ_t. We consider here only a general nondegenerate situation when the corresponding transition magnetic field, H_ℓ_0,σ,δ_t, is separated from other typical fields. The derivation in Appendix <ref> gives the result for the singular contribution at zero temperature

    𝒥_1(ω_c)-𝒥_1(ω_ℓ_0,σ,δ_t)
       ≈-  δ_t/2π√(μ_h+2δ_tt_z^h/t_z^h)√(|1-ω_c/ω_ℓ_0,σ,δ_t|) θ[δ_t(1-ω_c/ω_ℓ_0,σ,δ_t)]
       ×𝒢_ℓ_0[μ_h+2δ_tt_z^hcos^2Q_z2/ω_c]

with θ(x) being the step function and

    𝒢_ℓ(x)≡∑_m=ℓ^∞m!/2^m(m-ℓ)!ℓ!1/m+1-2x.

We see that the square-root singularity 𝒥_1(H)∝√(|H-H_ℓ_0,σ,δ_t|) appears near the transition point when the chemical potential is inside the Landau-level miniband, i.e., H≳ H_ℓ_0,σ,- or H≲ H_ℓ_0,σ,+. It reflects the pairing enhancement caused by the square-root divergence of the density of state at the edge of one-dimensional miniband. Finite temperature smears this singularity.



  §.§.§ Resonant cases for the uniform state (Q_z=0)


In the uniform state, Q_z=0, the resonant condition corresponds to the matching of the Zeeman
spin-splitting energy, 2γ_zω_c, and the Landau-level energy spacing, j_zω_c,
giving 2γ_z=j_z with integer j_z≥ 0. 
The divergence occurs when the Fermi level matches the van Hove singular point at k_z=0
or π corresponding to the magnetic field

    ω_c=μ_h∓2t_z/ℓ_0+(j_z+1)/2.

Figure <ref>(a) illustrates the simplest case with j_z=0.
We derive in Appendix <ref> the following asymptotic behavior of 𝒥_1(H,T,0) for T≪ω_c, t_z^h and arbitrary ℓ_0 and j_z:

    𝒥_1(H,T,0)∼𝐶(2ℓ_0+j_z)!/2^2ℓ_0+j_z(ℓ_0+j_z)!ℓ_0!ω_c/√(t_z^hT),

where 𝐶=(2√(2)-1)ζ(32)/(8π^3/2)≈ 0.1072 and ζ(x) is the Riemann zeta function [ζ(3/2)≈ 2.6124].
We can see
that the 𝒥_1 diverges as ∼1/√(T) manifesting the enhancement of pairing. The
coefficient decreases with increasing interlayer hoping t_z^h and with increasing of the
resonance order described by the integers ℓ_0 and j_z. We note that for any small deviations from the condition (ℓ_0+1/2j_z+1/2)ω_c=μ_h±2t_z^h, the total eigenvalue  𝒥_1-𝒜_1 approaches a finite limit at T→ 0. 
We can also note that this divergence is somewhat weaker than the 1/T divergence for the 2D case <cit.>, due to the smearing of δ-function singularity in the density of states at the Landau levels by the interlayer hopping. Such 1/√(T) divergence in the uniform state for zero or resonant spin-splitting factor has also been reported in Refs. Gruenberg:PRev176.1968,Maniv:RMP73.2001. 

For the strongest resonance at ℓ_0=0, j_z=0, and ω_c=2(μ_h∓ 2t_z^h) for k_z=0/π, we derived in Appendix <ref> the more accurate asymptotic for the total kernel eigenvalue

    𝒥_1-𝒜_1   ≈𝐶ω_c/√(t_z^hT)+ℛ_∓(t̅_z^h)-1/2lnω_c/π T_C-Υ_C,
    	
    ℛ_∓(t̅_z^h)     =1/2∑_m=1^∞(1/√(m(m∓8t̅_z^h))-1/m).

For k_z=0, the function ℛ_-(t̅_z^h) diverges for ω_c=8t_z^h. As ω_c=2(μ_h-2t_z^h), meaning that this result is only valid for μ_h>6t_z^h.

We mention that the spin splitting is a fixed material's parameter and therefore the resonance cases, 2γ_z=j_z, are exceptional. To some extent, the effective spin splitting can be tuned by tilting the magnetic field <cit.>. We note, however, that in contrast to the two-dimensional case <cit.>, in layered superconductors the in-plane magnetic field also influences the orbital motion of quasiparticles, meaning that the problem of the upper critical field for this case requires separate consideration.



  §.§.§ Resonant cases for alternating FFLO state (Q_z=π)



The resonance enhancement of pairing may also take place in the alternating FFLO state because in this case the modulation wave vector Q_z=π couples the van Hove singular points at k_z=0 and π. Such enhancement is a unique property of a layered superconductor with open Fermi surface and it only exists for specific relations between the band parameters μ_h and t_z^h. 
Writing the full quasiparticle energy in the magnetic field as  
ξ_±^h(ℓ,k_z)=μ_h-ω_c(ℓ+1/2±γ_z)-2t_z^hcos k_z, 
we can identify that the resonance conditions are realized when the spin-down (spin-up) energy at k_z=0 at the Landau level ℓ_0 simultaneously matches with the spin-up (spin-down) energy at k_z=π at the Landau level ℓ_π and with the chemical potential, i. e.,
ξ_∓^h(ℓ_0,0) =ξ_±^h(ℓ_π,π)=0 [see the first-case example with ℓ_0=ℓ_π=0 in Fig. <ref> (c)].
This gives the conditions


    ω_c     =4t_z^h/ℓ_π-ℓ_0±2γ_z=2μ_h/ℓ_0+ℓ_π+1,
    μ_h    =ω_c/2(ℓ_0+ℓ_π+1)
     =2t_z^h(ℓ_0+ℓ_π+1)/ℓ_π-ℓ_0± 2γ_z.



To derive the low-temperature behavior of 𝒥_1,  we substitute the conditions in Eqs. (<ref>) and (<ref>) into Eq. (<ref>) with  Q_z=π. The derivation described in Appendix <ref>
yields the low-temperature asymptotic for generic case of noninteger 2γ_z:

    𝒥_1(H,T,π)≃𝐶/2(ℓ_0+ℓ_π)!/2^ℓ_0+ℓ_πℓ_0!ℓ_π!ω_c/√(Tt_z^h).

where the numerical constant 𝐶 is defined above, after Eq. (<ref>).
We note that this additional divergence is purely a consequence of the interplay between the Landau quantization and interlayer tunneling. In particular, such FFLO instability is absent in the 3D case with k^2_z dispersion. 

For the strongest resonance, ℓ_0=ℓ_π=0,  μ_h=t_z^h/γ_z (γ_z<0.5), and ω_c=2μ_h, we derived the more accurate asymptotic

    𝒥_1-𝒜_1≃𝐶/2ω_c/√(Tt_z^h)-1/2ln2ω_c/π T_c-Υ_C.

In the real systems such instability can be obtained by tuning the chemical potential μ_h  by doping or pressure. 



  §.§.§ Logarithmic divergence at the optimal modulation wave vector


In the case when the chemical potential crosses two Landau-level branches with opposite spin orientations, the kernel eigenvalue also has logarithmic divergence at the wave vector connecting the crossing points [see Fig. <ref>(b)]. This divergence originates from the one-dimensional character of the electronic spectrum. The similar result was also found in the study of a three-dimensional model with parabolic Landau minibands <cit.>. For illustration, we consider such divergence for such branches belonging to the same Landau level with index ℓ_0. In this case, the diverging term in Eq. (<ref>) is the one with m=2ℓ_0 and ℓ=ℓ_0:

    𝒥_1,ℓ_0    =ω_c/4(2ℓ_0)!/2^2ℓ_0(ℓ_0!)^2
    ×   ∫ _-π^πdk_z/2π∑_δ=±1tanh[ω_c(ℓ_0+1/2+δγ_z)-μ_h+2t_z^hcos(k_z+δQ_z2)/2T]/ω_c(2ℓ_0+1)-2(μ_h-2t_z^hcos k_zcosQ_z2).

For the optimal modulation wave vector Q_z=Q_op shown
in Fig. 2(b), we have the relations 

    ω_c(ℓ_0+12+δγ_z)-μ_h+2t_z^hcos(k_z0+δQ_op2)=0,

for δ=±1, where k_z0± Q_op/2 are the Fermi wave vectors for the two considered branches. These relations determine 
Q_op as

    Q_op=   ∑_δ=± 1δarccos[μ_h-ω_c(ℓ_0+12+δγ_z)/2t_z^h].

 

In the limit T≪ t_z^h the dominating contribution to 𝒥_1,0
is coming from the regions near the Fermi momenta meaning that we
can use expansion k_z=k_z0+k̃ with k̃≪ k_z0 and approximate

    𝒥_1,ℓ_0   ≈(2ℓ_0)!/32π2^2ℓ_0(ℓ_0!)^2ω_c/t_z^hsin k_z0cosQ_op2
    ×   ∑_δ=±1∫ _-π^πdk̃/k̃tanh[t_z^h/Tsin(k_z0+δQ_op2)k̃].

Evaluating the integral, we obtain 

    𝒥_1,ℓ_0   ≈(2ℓ_0)!/8π2^2ℓ_0(ℓ_0!)^2ω_c/t_z^hsin k_z0cosQ_op2
    ×   ln[
    t_z^h/T√(|sin(k_z0-Q_op2)sin(k_z0+Q_op2)|)].

This logarithmic divergence for T→ 0 is similar to the well-known Peierls divergence of the
electronic susceptibility for one-dimensional metals at the nesting wave vector which, in
particular, leads to the charge-density-wave transition. In the simplest scenario when the chemical potential crosses only two minibands,  the kernel eigenvalue diverges only at one wave vector. In more complicated cases with many minibands at the Fermi level, the divergence takes place at several wave vectors and one can expect competition between multiple ground states, similar to the situation considered in Ref. <cit.>.




 §.§ Magnetic-field dependences of quantum kernel eigenvalue 𝒥_1 for different cases





Now, we discuss the general behavior of the kernel eigenvalue 𝒥_1. In contrast to the monotonic behavior of the quasiclassical result, the low-temperature divergences discussed in the Sec. <ref> may lead to the emergence of strong peaks in 𝒥_1 at the specific magnetic fields and modulating wave vectors Q_z which indicate the enhancement of pairing. The magnitude and location of these peaks are sensitive to the details of electronic band properties. 



  §.§.§ Uniform state: Q_z=0


Figure <ref> illustrates the effects of interlayer tunneling, temperature, and Zeeman spin-splitting in the uniform state (Q_z=0). Without spin splitting, the interlayer tunneling broadens the Landau levels and splits a single divergent peak into two smaller peaks at low temperatures [see Fig. <ref>(a)]. The smaller peaks correspond to the matching between the chemical potential and the van Hove singularities in the spectrum at k_z=0 and π [Eq. (<ref>)]. Their values are estimated by Eq. (<ref>). The shaded regions mark the field ranges for which the chemical potential is located within the broadened Landau levels, and their boundaries correspond to the crossing of the van Hove singularities.

Figure <ref>(b) shows that the peaks are rapidly smeared by the thermal effects and disappear at T∼ω_c/2. In this figure, we compare the exact kernel eigenvalue with the quasiclassical results (dotted-dashed lines), which are calculated by using

    𝒥^qc_1(H,T,Q_z)   =2∫^∞_0ds⟨lntanh(π T/ω_cs)exp(-μ̃_hs^2)
    ×   [μ̃_hscos(2γ̃_zs)+γ̃_zsin(2γ̃_zs)]⟩_z.

This quasiclassical approximation of 𝒥_1 can be obtained from the quantum kernel in Eq. (<ref>) by assuming small-field limit<cit.>, ω_c≪μ_h with the temperature range, ω_c≪ T≪μ_h. One can see that this approximation yields a monotonic behavior for the kernel eigenvalue in the magnetic field and gives a good approximation at low fields. 

We can see in Fig. <ref>(c) that the divergent peaks are also rapidly suppressed by finite spin splitting γ_z meaning that the uniform state is highly susceptible to the Zeeman effect. The smaller peaks reappear if the resonant conditions are met (the solid red line for γ_z=0.5). 
We also observe a noticeable enhancement of 𝒥_1 in the field ranges where the chemical potential crosses both spin up and spin down minibands for the zero Landau level (e. g., in the range 1.5<ω_c/μ_h<2.4 for γ_z=0.1). We point, however, that in this range the maximal  𝒥_1 is not at Q_z=0. We discuss the Q_z dependences in the next subsection.



  §.§.§ FFLO modulated state: Q_z≠0


We now consider the behavior of the quantum kernel eigenvalue at finite modulation wave vectors Q_z.
Figure <ref> shows the dependences of 𝒥_1 on Q_z for fixed electronic-spectrum parameters, fixed magnetic field, and for three spin-splitting factors representing different cases.
In the resonant cases with half-integer spin-splitting factors the kernel eigenvalue always has the maximum at Q_z=0, as illustrated by the curve with γ_z=0. The finite-Q_z modulated states are not favorable in these cases, due to the absence of strong Zeeman pair breaking at the Landau levels [see Fig. <ref>(a)], and the uniform state always dominates. Away from the uniform-state resonances, 2γ_z≠ j_z, the Zeeman pair breaking favors developing of the finite-Q_z states. The second curve with  γ_z=0.2 illustrates the situation when the chemical potential crosses both spin-up and -down zero Landau-level minibands. In this case, the maximum is realized at the modulation wave vector connecting the minibands Fermi momenta [see Fig.  <ref>(b)]. This modulation vector varies with the magnetic field. The third curve with γ_z=0.3 illustrates a special situation when the chemical potential matches two opposite van Hove points at k_z=0 and π [see Fig.  <ref>(c)]. Such matching is realized when the parameters satisfy the relation given by Eqs.  (<ref>) and (<ref>) (the plot is made for ℓ_0=ℓ_π=0). In this case, the maximum is realized at Q_z=π, favoring the alternating FFLO state.


Figure <ref> illustrates the magnetic field dependences of  𝒥_1 for different Q_z in two situations with noninteger 2γ_z. We also show the maximal 𝒥_1 with respect to Q_z. The light blue and pink regions mark the crossing of the chemical potential with the mismatching spin-up and -down Landau-level minibands. This regions are limited by the miniband-edge magnetic fields given by Eq. (<ref>) with δ_t=± 1 while ℓ and σ are fixed. The simultaneous crossing of both such bands at the chemical potential only occurs in the overlapping regions in which the pairing favors finite-Q_z states as illustrated in Fig. <ref>(b). Figure <ref>(a) corresponds to a general situation without any resonances. In this case, the states with different Q_z become favorable in different fields. In particular, for Q_z=π/2 the maximum at ω_c/μ_h≈ 2.1 corresponds to the condition for optimum modulation vector Q_op(ω_c)=Q_z [Eq. (<ref>)].  We can also see that the pairing strength is noticeably enhanced in the region 1.3<ω_c/μ_h<2.3 when the chemical potential crosses zero-Landau level minibands for both spin orientations. The choice of parameters in Fig. <ref>(b) allows for Q_z=π singularity [see Eq. (<ref>) and the illustration in Fig. <ref>(c)]. At ω_c/μ_h=2, the chemical potential matches simultaneously two van Hove singular points: the spin-up energy at k_z=0 and spin-down energy at k_z=π. This matching leads to the additional divergent peak at this point strongly favoring Q_z=π instability and may lead the formation of the superconducting state with the alternating sign of the order parameter between the layers. In both cases, the optimal kernel eigenvalue (blue circles) significantly exceeds the one for the uniform state (black line) almost in the whole field range. 


To obtain a better presentation for the range of parameters where the shallow band favors nonuniform states,
we expand 𝒥_1 with respect to small Q_z:

    𝒥_1(H,T,Q_z)≈𝒥_1(H,T,0)+1/2𝒥”_1(H,T,0)Q_z^2.

If 𝒥”_1(H,T,0) is positive, then Q_z=0 does not maximize 𝒥_1 and the finite non-zero Q_z state is favorable. For better visual impression of the parameter range where the shallow band favors nonuniform states, we present in Fig. <ref> the regions of positive  (shaded) and negative (unshaded) 𝒥”_1(H,T,0) in the magnetic field–spin splitting plane for representative parameter t_z^h/μ_h=0.2 and two values  of temperature T=0.02μ_h(a) and 0.2μ_h(b).  For comparison, we also show in Fig. <ref>(c) the same diagram for the quasiclassical kernel eigenvalue, Eq. (<ref>). In the latter case an analytical analysis <cit.> suggests that the boundary should behave as ω_c/μ_h∝ 1/γ_z^2 while the numerical fitting gives ω_c/μ_h≈ 1/(0.14+1.07γ_z)^2. If the upper critical field with Q_z=0 state falls into the shaded region, this means that the shallow band favors the nonuniform state, and the order parameter has a tendency to develop nonzero Q_z modulation in real space. We see that at low temperatures the shaded regions  cover substantial part of the phase diagram except locations around the integer 2γ_z. Especially surprising is that they extend to very low fields, down to ω_c/μ_h∼ 0.1. These low-field regions are eliminated at higher temperatures. 

To avoid misunderstanding, we note that for two-band systems the nonuniform states do not automatically appear in the shaded regions, because the ground state is determined by both bands and the deep band typically favors the uniform state, 𝒥_2^''(H,T,0)<0. The FFLO modulation may appear below certain temperature T_FFLO, at which the second derivative ∂^2 H_C2/∂ Q_z^2 changes sign. 
From the general equation for the upper critical field [Eq. (<ref>)] we can derive the equation for  this temperature, T=T_FFLO,

    𝒥_1^''(H,T,0)/W_11+𝒜_1(T)-𝒥_1,0(H,T)
    +𝒥_2^''(H,T,0)/W_22+𝒜_2(T)-𝒥_2,0(H,T)=0

with 𝒥_α,0(H,T)≡𝒥_α(H,T,0), 
in which we have to substitute H→ H_C2(T) at Q_z=0. This equation may have solution only if (i) 𝒥_1^''(H,T,0)>0 at low temperatures [i.e., the parameters are in the shaded region in Fig. <ref>(a)] and (ii) the ratio W_11/W_22  characterizing the relative weight of the shallow band is not too small.



§ SUPERCONDUCTING INSTABILITIES AND PHASE DIAGRAMS


The exact shape of the H-T phase diagram of a layered multiple-band superconductor depends on many parameters: electronic spectra properties of the bands, spin-splitting factors, and structure of the coupling matrix. The most crucial factor is the quantum pairing kernel eigenvalue for the shallow band, 𝒥_1, whose shape at low temperatures is extremely sensitive to the electronic parameters of the shallow band, μ_h, t_z^h, and γ_z, as we discussed in the previous section. With the finite interlayer tunneling, the k_z dispersion in the Landau-quantization spectrum can lead to complex behavior of the superconducting state in the magnetic field.

In addition, the Zeeman spin splitting in the shallow band favors the formation of the FFLO states (see Figs. <ref> and <ref>). As the Landau quantization strongly influences the Cooper pairing in the bands with relatively low Fermi energies, this brings a natural general question about the role of quantum effects in the FFLO instability, even in the common case of a single-band superconductor in the quasiclassical limit. Surprisingly, this important topic has been never addressed before and  we leave it to a separate paper. Here, we focus on the general behavior of the H-T phase diagrams for the two-band system where the shallow band is close to the extreme quantum limit. We illustrate this behavior for several representative cases.

The most spectacular consequence of the Landau quantization is the emergence of the high-field superconducting states at the magnetic fields corresponding to matching of the chemical potential with the Landau levels. These states are most pronounced in the two-dimensional case and for the lowest Landau level <cit.>.  The interlayer tunneling smears the delta-function singularity in the density of states, leading to suppression of the high-field superconductivity. 
To quantify these effects, we first calculate the transition temperatures in magnetic field, T_C2, in the resonance cases for which the kernel eigenvalue diverges as ∼ 1/√(T), see Sec. <ref>. For this we find solution of Eq. (<ref>) assuming that T_C2≪ T_C, t_z. 


First, we evaluate the transition temperature for the strongest resonance in the uniform state without Zeeman spin-splitting (γ_z=0), T^(0)_C2, which is realized when the chemical potential matches the van Hove singular points at k_z=0,π at the magnetic fields ω_c=2(μ_h∓ 2t^h_z). 
In this case, substituting the low-T asymptotics of 𝒥_α-𝒜_α given by Eqs. (<ref>) and (<ref>) into Eq. (<ref>), we derive

    T^(0)_C2=𝐶^2ω_c^2/t^h_z[W_11ln r_C^(0)/2W_22+ln r_C^(0)-ℛ_∓(t̅_z^h)
    +1/2lnω_c/π T_C+Υ_C]^-2.

This result allows us to understand better the typical sizes of the high-field reentrant regions.
A similar result can be derived for the strongest alternating-state resonance with Q_z=π realized when the parameters satisfy the relation  μ_h=t_z^h/γ_z.  In this case, the low-temperature asymptotics of 𝒥_1-𝒜_1 is given by Eq. (<ref>), leading to the π-state transition temperature at ω_c=2μ_h:

    T^(π)_C2=𝐶^2ω_c^2/4t^h_z[
    W_11ln r_C^(π)/2W_22+ln r_C^(π)+1/2ln2ω_c/π T_C+Υ_C]^-2

with r_C^(π)≡ r_C(Q_z=π). 
Even though the overall scale of the above transition temperatures is given by ω_c^2/t^h_z with the small numerical factor 𝐶^2≈ 0.0115, their absolute values are very sensitive to the structure of the coupling matrix. The latter dependence is given by the first term in the square brackets in Eqs. (<ref>) and (<ref>). In general, these transition temperatures are not vanishingly small only if the shallow band gives substantial contribution to pairing. 

In particular, in the case when the deep band dominates pairing, Λ_ee>Λ_hh,|Λ_eh|,|Λ_he|, the constants
W_αα can be estimated as W_11≈Λ_ee/𝒟_Λ
and  W_22≈Λ_ehΛ_he/(Λ_ee𝒟_Λ), meaning that
|W_22|≪ |W_11|. In addition, inequalities ln r_C ≪ |W_22| and W_11ln r_C/W_22≈Λ_ee^2ln r_C/Λ_ehΛ_he≫Υ_C,ln(μ_h/T_c) are typically satisfied, leading to the estimate 
T^(0)_C2≈(𝐶^2ω_c^2/t^h_z) [2Λ_ehΛ_he/(Λ_ee^2ln r_C^(0)) ]^2. In this case, the transition temperature has the additional small factor which rapidly decreases with decreasing of the interband couplings constants, Λ_eh and Λ_he, and with increasing the deep-band coupling constant, Λ_ee. On the other hand, we see that T^(0)_C2 increases when the resonance field approaches the orbital field of the shallow band, r_C^(0)→1. In the interband-coupling scenario, |Λ_eh|,|Λ_he|>Λ_ee,Λ_hh, the ratio W_11/W_22≈ 1/2, meaning that there is no additional smallness caused by the coupling factor. 

We illustrate the typical behavior of T_C2 in Fig. <ref>, by plotting the transition temperatures for uniform and alternating states as function of the ratio Λ_h,e/Λ_0,e which controls the shallow-band weight. The plots are made for fixed effective coupling constant Λ_0,e=0.2 which fixes the zero-field transition temperature T_C. As expected, in all cases the transition temperatures rapidly decrease with decreasing of the shallow-band contribution to pairing. For comparison, we present the T_C2-dependences for two values of the interband hopping energy t_z^h=0.45T_C and  0.15T_C. We see that T_C2's are much higher for smaller  t_z^h.

We also note that, in spite of the kernel low-temperature divergences in the resonance cases, the
two-band system does not always become superconducting at the corresponding magnetic fields. Indeed, the
above solutions for T^(0)_C2 and T^(π)_C2 are only valid if the arguments in the
square brackets of Eqs. (<ref>) and (<ref>) are positive. When the interband
coupling dominates, 𝒟_Λ<0, the parameters W_αα are
negative. In this case, T^(α)_C2 may vanish if ln r_C^(α)=2|W_22|. 
As |W_22|≫ 1, this may only realize if the resonance magnetic field is much larger than the
deep-band orbital field,  r_C^(α)≫ 1, corresponding to the very large Fermi energy, μ >
100 T_C. 



We now present numerically computed magnetic field–temperature phase diagrams for different typical situations. We first illustrate the influence of the interlayer tunneling on the reentrant high-field states in the resonant cases with integer 2γ_z. Figure <ref> shows the H-T diagrams for  γ_z=0.5, μ_h=3T_C, and different hopping energies t_z^h. The diagrams are made for the case of dominating interband coupling scenario (Λ_hh=Λ_ee=0, Λ_he=Λ_eh=0.3), for which the discussed reentrant effects are very pronounced. 
The vertical scale of the magnetic field cm_hT_C/e in the plot can be rewritten as 7.44 T(m_h/m_0)(T_C/10K). Two reentrant states are realized for selected parameters, around the zeroth and first Landau levels. We can see that the interlayer tunneling smears these states. It splits the single peak in T_C2(H) around the field ω_c=μ_h(ℓ+1/2+1/2j_z)^-1 into the two peaks with maximums approximately located at ω_c=(μ_h±2t_z^h)(ℓ+1/2+1/2j_z)^-1 (in the plot ℓ=0,1 and j_z=1). The latter fields correspond to the crossing of the Fermi level and the van Hove singular points of the shallow h-band given by Eq. (<ref>), as illustrated by the insets. Since the increase of t_z^h reduces the low-T divergent peaks in 𝒥_1(H) dependence [Eq. (<ref>)] the reentrant T_C2 are always smaller than in the two-dimensional case. We emphasize again that for the resonant cases, 2γ_z=j_z, the FFLO state is not favorable and Q_z=0 always gives the largest T_C2.

We now discuss the formation of the interlayer FFLO states due to spin splitting in the shallow band for a more general case of  noninteger  2γ_z. Figure <ref> shows the representative H-T phase diagram for the spin-splitting factor γ_z=0.3, the interlayer hopping energy t_z=0.45T_C, and three values of the Fermi energy, μ_h=T_C (a), 1.25T_C (b), and 1.5T_C (c).  We again consider the interband-coupling scenario with the same coupling constants as in the previous figure. We plot the H_C2 lines with Q_z at the optimal values (Q_op) and the lines are color coded by Q_op. For comparison, we also show the Q_z=0 transition lines. 

In the case  μ_h=T_C [Fig. <ref>(a)], the Zeeman spin splitting in the Landau levels of the shallow h band leads to the pairings that favor the FFLO modulation along the out-of-plane direction for T<0.18T_C. The formation of this state leads to substantial enhancement of the upper critical field at low temperatures. The modulation wave vector at the transition  Q_op rapidly increases with decreasing temperature, reaching 2.17 at T→ 0. At somewhat higher Fermi energy μ_h=1.25T_C [Fig. <ref>(b)] the pairing at overlapping spin-up and spin-down zero-Landau-level minibands leads to the formation of the separated reentrant high-field state. This state also has z-axis modulation with the wave vector Q_op shown in the inset and its value is very close to the distance between spin-up and -down Fermi momenta given by Eq. (<ref>). Note that the reentrant FFLO state occupies a quite extended field range corresponding to 1.5<ω_c/T_C<2.75, where the spin-up and -down zero-Landau-level minibands overlap.


As demonstrated in Sec. <ref>, a layered superconductor with shallow band may have special resonance situations leading to appearance of the alternating FFLO state. Namely, if μ_h satisfies the condition in Eq. (<ref>), 𝒥_1 diverges for T→ 0 at the magnetic field that is given by Eq. (<ref>). This divergence typically gives rise to the FFLO modulation with Q_z=π. Figure <ref>(c) shows an example of the H-T diagram for such situation with ℓ_0=ℓ_π=0 in Eqs. (<ref>) and (<ref>) corresponding to the strongest resonance. We can see that the small reentrant region indeed appears at the high magnetic field. However, only part of this region is occupied by the alternating state. In contrast to the uniform state, T_C2(H) has only single maximum, as in the two-dimensional case. We see that, depending on parameters, the reentrant states can be either well separated from or very close to the main superconducting region.


For better understanding of general trends, we plot in  Fig. <ref> the dependences of the low-temperature upper critical field on the shallow-band Fermi energy μ_h. The upper curve is made for the same parameters as in Fig. <ref>. We can see that this dependence has several salient features. The H_C2 curve has a clear tendency to follow one of the miniband-edge field lines and it sharply turns away from this line close to crossing points with the second miniband-edge field line. As a consequence, the regions of multiple H_C2 values corresponding to reentrant behavior appear below some of these crossing points. Near the crossing of  0↓-t_z and 0↑+t_z lines corresponding to the matching of the opposite van Hove singularities for the zero Landau level, the region of the alternating FFLO state is realized. 

So far, we only considered the purely interband-coupling scenario for which the quantum effects from the shallow band are quite pronounced. The behavior, however, is sensitive to the structure of the coupling-constants matrix reflecting the pairing mechanism. 
In addition to the interband-coupling case, Fig. <ref> also shows the H_C2 - μ_h dependences at low temperatures for two coupling-constant sets with finite Λ_ee. We can see that with increasing deep-band pairing weight, the upper critical field decreases and moves to the region of miniband-edge field lines for higher Landau levels, especially for '+t_z' lines corresponding to larger-area cross section of the Fermi surface at k_z=π. Nevertheless, the main trends remain: the H_C2 curves still tend to follow the miniband-edge field lines and the reentrant regions appear near the crossing points (at least, for the case Λ_ee=0.2 and Λ_eh=0.3). Note that the reentrance completely vanishes when the deep-band coupling is too strong (the lowest curve in Fig. <ref>). The ground state at low temperature is modulated for all studied cases, but the optimal modulation wave vector Q_op  progressively decreases with increasing the deep-band weight. 
The middle curve includes the crossing of 0↓-t_z and 1↑+t_z lines,  where the shallow band strongly favors the alternating state. This state, however, is not formed due to the large weight of the deep band. The modulation wave vector is sharply enhanced when μ_h approaches the crossing-point value 1.125T_C but it only reaches ∼ 1.1 at the maximum.


The temperature range of reentrant states rapidly reduces with increasing deep-band coupling.
We illustrate this trend in Fig. <ref> in which we present the evolution of the reentrant region  of Fig. <ref>(b) with increasing deep-band coupling constant Λ_ee. We can see that the size of this region rapidly shrinks with increasing Λ_ee. 
This counterintuitive behavior is again caused by the reduction of the shallow-band weight in the pairing. In the regime of dominating deep-band coupling, Λ_ee>Λ_eh, the reentrant region becomes practically invisible for these electronic-spectrum parameters.
 
We discuss in Appendix <ref> the scenario in which superconductivity is dominated by the deep-band coupling with Λ_ee>Λ_he. The resulting H-T diagrams are not qualitatively different from the representative case in Figs. <ref> and <ref>, but there are substantial quantitative differences. In the resonant case, the sizes of reentrant regions are much smaller than for the interband-coupling case due to the strong reduction of the shallow-band contribution to the pairing (see also Fig. <ref>). In the nonresonant case, the reentrant behavior only appears either for smaller interlayer tunneling or smaller Fermi energy of the deep band in comparison with the parameters considered in this section. Nevertheless, the FFLO instability in the main superconducting region always appears at low temperatures.



§ SUMMARY AND DISCUSSION


In summary, we have investigated the pairing instabilities of a clean two-band layered superconductor in the magnetic field oriented perpendicular to the layers.  
In this system, the interlayer tunneling lifts the degeneracies in all Landau levels and transforms them into the dispersive minibands along the out-of-plane momentum (k_z) with two van Hove singularities at k_z=0 and π. We explored the possible pairing instabilities as the Landau minibands cross the chemical potential in the vicinity of the first Lifshitz transition (μ_h≳ 2t_z^h in Fig. <ref>).

Solving the linearized gap equation [Eq. (<ref>)], we found that the Landau quantization leads to strong 1/√(T) divergences in the shallow band's pairing kernel if either the same [Eq. (<ref>)] or opposite [Eqs. (<ref> and <ref>)] van Hove points match at the chemical potential. The former matching condition gives the uniform superconducting state and the latter one favors the alternating FFLO state with Q_z=π modulation. For general Zeeman spin-splitting energy, the pairing kernel has ln T divergence in the pairing channel with the optimal wave vector Q_op [Eq. (<ref>)] corresponding to the difference between the spin-up and spin-down Fermi momenta.
These logarithmic divergences are similar to the FFLO instabilities in quasi-one-dimensional superconductors <cit.> and they strongly promote the formation of the FFLO states with Q_z=Q_op. 

Furthermore, we studied the magnetic field versus temperature phase diagrams and found that the shallow-band divergences yield a complex reentrant behavior in the high-magnetic-field superconductivity. The properties of these reentrant states are highly sensitive to the shallow-band parameters, μ_h, t_z^h, γ_z, as well as to the coupling-matrix structure. In general, the interlayer tunneling smears the Landau-level densities of states which reduces the quantum effects in the resonant cases. Therefore, the reentrant transition temperatures are smaller than those in the two-dimensional monolayer<cit.>. On the other hand, the interlayer tunneling enables the formation of the FFLO modulations which mitigate the Zeeman suppression in the reentrant states. As a result, the reentrant region is stretched to a wider magnetic-field range [see Figs. <ref>(a) and (b)]. Furthermore, we also found the reentrant alternating-FFLO states [see Fig. <ref>(c)] due to the interplay between the Landau quantization and interlayer tunneling. 
The similar alternating-FFLO state has been also suggested for the layered system in the magnetic field applied parallel to the layers<cit.>. The formation of these alternating-FFLO states is the consequence of the compensation between the Zeeman energy and the interlayer hopping energy.


In this study, we only performed the linear stability analysis and, strictly speaking, the computed phase lines describe instability of the normal state. These lines would correspond to true second-order phase transitions only if the coefficient for the quartic order-parameter term in the free energy is positive. Otherwise, superconducting state will emerge via a first-order transition.
Such first-order scenario for the emergence of the FFLO state indeed realizes in quasi-two-dimensional materials with very strong Zeeman effect for the magnetic field perpendicular to conducting layers <cit.>. This analysis, however, has been performed only for the quasiclassical regime and it does not take into account the effect we stress here, promotion of the FFLO instability by the Landau quantization. A full nonlinear consideration with proper accounting for the quantum effects remains to be done. 

The described behavior can only be observed in very clean materials in which the scattering time τ satisfies the condition τω_c ≫ 1. This condition implies that superconducting transition should take place in the region of noticeable quantum magnetic oscillations of magnetization and conductivity (de Haas–van Alphen and Shubnikov–de Haas effects). This is obviously the most stringent requirement for observation of the described anomalous behavior of the upper critical field.
We expect that impurity scattering suppresses such behavior similar to the thermal noise. Indeed, in the interpretation of the quantum-oscillation phenomena scattering is frequently accounted for by introducing a fictitious temperature proportional to the scattering rate, so-called Dingle temperature.  

We only considered the simplest nearest-neighbor interlayer tunneling process. Additional strong pairing instability may arise if new van Hove singularities appear in the Landau band due to the complicated interlayer tunneling processes. Furthermore, for sufficiently strong Zeeman splitting and/or small Fermi surfaces for all bands, the higher-Landau-level gap eigenfunction may be relevant <cit.>. Influence of the quantization effects on this gap solution and modification of the reentrant behavior are other interesting topics for future study.

The quantum effects considered in this paper may be relevant to the high-magnetic-field behavior of several iron-based superconductors, in particular, for FeSe and LiFeAs single crystals. The simplest compound FeSe has the transition temperature ∼ 8K<cit.> and the low-temperature upper critical field ∼17T<cit.>. The material can be made clean allowing for observation of quantum oscillations at fields > 20T <cit.>. Its band structure is composed of hole and electron bands with rather small Fermi surfaces. In particular, Shubnikov–de Haas oscillations show that the smallest Fermi energy for the electron and hole bands are only 3.9 and 5.4 meV, respectively <cit.> (the first energy corresponds to the ratio ϵ_F/ω_c∼ 4 at
H_c2). Similar Fermi energies are also observed by ARPES <cit.>. Moreover, the additional field-induced transition inside the superconducting state has been observed in this compound near 13 T by thermal conductivity measurement <cit.>. The plausible intriguing interpretation of this transition is the onset of the FFLO modulation along the magnetic field <cit.>. Such explanation is also consistent with the low-temperature upturn of the upper critical field. If this interpretation is correct, then this FFLO transition is likely influenced by quantum effects. Furthermore, the transition temperature in FeSe can be increased by pressure up to 20 K at 25 kbar and the upper critical field increases above 35 T at 14 kbar. As a consequence, at high pressures, the superconducting transition in the magnetic field takes place in the region of pronounced quantum oscillations <cit.>. As for the most shallow band ϵ_F/ω_c∼ 2 at H_C2 in this region, it is clear that the quantum effects strongly influence the superconducting instability.

The 111 compound LiFeAs is one of the few stoichiometric iron-based superconductors with transition temperature ∼18 K<cit.> and zero-temperature c-axis H_C2 ∼24T<cit.>.
The quantum oscillations also have been reported for this material in Refs. <cit.>. 
It has the hole-like shallow band located at the zone center with 
closed Fermi surface and tiny Fermi energy μ_h∼3 meV <cit.>. It is important to note that this band actually has the largest superconducting gap <cit.>, meaning strong participation in the formation of the superconducting state.
In addition, this band can be completely depleted by small Co doping <cit.>. The band with so small Fermi energy definitely should cause noticeable quantum effects near the upper critical field if the material can be made sufficiently clean. 
The quantitative predictions for a particular material require a detailed knowledge of the electronic spectrum, spin-splitting factors for all bands, and coupling matrix.

The high-field reentrant superconductivity has been reliably found at least in one material, Eu-doped Chevrel phases, Eu_xSn_1-xMo_6S_8 <cit.>, where an isolated semi-elliptical superconducting region has been observed in the H-T diagram for T<1 K and a very wide magnetic field range, 4T<H<22.5T. The interpretation of this strong reentrance was based on the Jaccarino-Peter effect<cit.>, the compensation of the Zeeman spin splitting due to the interaction with the local magnetic moments of the Eu dopants. However, such interpretation also required an assumption of the extremely weak orbital effect of the magnetic field (a huge Maki parameter α_M≈ 4.8 was used in the theoretical fits). The reason for this assumption is not very clear,  as the material is isotropic. We cannot, therefore, exclude that the quantization effects may also play a role in the formation of the reentrant region.

The newly discovered Dirac/Weyl semimetals <cit.> (see also review <cit.>) with very small and tunable Fermi energies are another interesting systems to search for the high magnetic field induced superconductivity. Of particular interests to this paper are the superconducting states that are found in some compounds under high pressure<cit.>. It has been demonstrated in Ref. Rosenstein:PRB96.2017 that the orbital quantization in superconducting materials with such Dirac-type spectrum can also lead to the reentrant superconductivity in high magnetic fields. Furthermore, the similar reentrant behavior has also been found in the theoretical study on the Dirac-type surface states of a topological insulator <cit.>. In contrast to the conventional quadratic electronic dispersion, the linearly-dispersive band yields nonuniform Landau-level spacing and it has been argued that this intriguing feature makes the quantum limit much easier to attain <cit.>.

In addition to this feature, due to the small cyclotron effective mass near the Dirac point, the Landau quantization is more robust against disorders <cit.>. 
So far, the above studies only considered the Landau-quantization effects in the uniform states. The possible enhancement of the FFLO instabilities in the three-dimensional Dirac/Weyl semimetal due to quantum effects is an interesting problem to be explored in the future. 

Finally, we note that the study in this paper only concentrated on the correlation effects in the BCS pairing channels. For real systems, other correlation effects in spin and charge channels are also likely to play an important role. Particularly, the Landau-level miniband is effectively a one-dimensional system, and the spin and charge fluctuations may have significant influence. Indeed, the strong competition between the pairing and density-wave instabilities may have dramatic effects and give rise to a new class of strongly correlated quantum states<cit.>. The rich strongly correlated phenomena due to the interplay between these fluctuations in quantum limit are certainly important problems for future consideration.


The authors would like to thank A. Buzdin, K. Matveev, T. Shibauchi, L. Balicas, Y. Kopelevich, and A. Gurevich for useful discussions. The work was supported by the U.S. Department of Energy, Office of Science, Materials Sciences and Engineering Division. K.W.S. was supported by the Center for Emergent Superconductivity, an Energy Frontier Research Center funded by the U.S. DOE, Office of Science, under Award No. DEAC0298CH1088.






§ CALCULATION OF THE KERNEL EIGENVALUES

In this appendix, we give details of the calculations of the kernels in zero and finite magnetic field.



 §.§ Zero magnetic field



The kernels, Eq. (<ref>) are composed by the zero-field Green's function, Eq. (<ref>). The effective coupling constant Λ_0,α [Eq. (<ref>)] is determined by




    Λ_0,α^-1=T/N_α∑_ω_n,j'∑_𝐤𝐤'
    k_zk_z'∫_𝐫'e^-i[(𝐤-𝐤')·(r-r')+(k_z-k'_z)(j-j')]
       ×1/(iω_n-ξ^α_𝐤+2t_zcos k_z)(iω_n+ξ^α_𝐤'-2t_zcos k'_z).

Since the system is translationally invariant, it is straightforward to integrate out 𝐫' and sum all j'. This yields





    Λ_0,α^-1=N_α^-1T∑_ω_n∑_𝐤k_z1/ω^2_n+(ξ^α_𝐤-2t_zcos k_z)^2 .

We further integrate out the (𝐤,k_z) by using ∑_𝐤=N_e∫^Ω_-Ωdξ^e for e band (∑_𝐤=N_h∫^μ_h_-Ωdξ^h for h band) with N_e=m_e/2π (N_h=m_h/2π) and ∑_k_z=∫^π_-πdk_z/2π. We therefore obtain






    Λ_0,e^-1=∑_ω_n∫^π_-πdk_z/2π[T/ω_ntan^-1x/ω_n]^x=Ω_x=-Ω,
       Λ_0,h^-1=∑_ω_n∫^π_-πdk_z/2π[T/ω_ntan^-1x/ω_n]^x=μ_h-2t_zcos k_z_x=-Ω.

Note that we assumed Ω≫ t_z. Further simplifying the above expressions, we obtain Eqs. (<ref>) and (<ref>) in the main text.



 §.§ Finite magnetic field





  §.§.§ Deep band



For the deep e band in which the Landau quantization does not play a role, we have the following quasiclassical result for the kernel eigenvalue:

    λ^e_ω_n,Q_z=   2/N_e∑_𝐤𝐤'k_z∫_0^∞ρdρexp[-i(𝐤-𝐤')·ρ-ρ^2/2l^2]/[iω_n-μ_zH-ξ^e_𝐤+2t_z^ecos (k_z-1/2Q_z)]
       ×1/[iω_n-μ_zH+ξ^e_𝐤'-2t_z^ecos (k_z+1/2Q_z)].

Setting 𝐤→𝐤-1/2𝐪 and 𝐤'→𝐤+1/2𝐪, and keeping only the linear order in 𝐪, we can approximate the dispersion near Fermi level as
ξ^e_𝐤→ξ^e_𝐤-12v^e·𝐪  and ξ^e_𝐤'→ξ^e_𝐤+12v^e·𝐪.
We now integrate out 𝐤 in Eq. (<ref>) with the above approximation by using ∑_𝐤=N_e∫^∞_-∞dξ^e. Extending the energy integration to a closed contour in the complex plane, this yields

    λ^e_ω_n,Q_z≈-πi∑_k_z𝐪∫_0^∞ρdρexp[-i𝐪·ρ-ρ^2/2l^2]{iω_n-μ_zH
       
    +12v^e·𝐪+t_z^e[cos (k_z-1/2Q_z)-cos (k_z+1/2Q_z)]}^-1.

The absolute value of the momentum 𝐤 in the Fermi velocity v^e(𝐤) is approximately determined by ξ^e_𝐤=t_z^e[cos (k_z-1/2Q_z)+cos (k_z+1/2Q_z)] giving 

    v^2_e≈2/m_e( μ+2t_z^ecos k_z cosQ_z/2).

Exponentiating the denominator as follows,

    λ^e_ω_n,Q_z≈2π∑_k_z𝐪∫_0^∞ρdρexp[-i𝐪·ρ-ρ^2/2l^2]
    ∫_0^∞dsexp{2ζ_ω s
       ×[ω_n+i[μ_zH-12v^e·𝐪-2t_z^esin (k_z)sin (12Q_z)]]},

with ζ_ω=sign(ω_n), we integrate out 𝐪, which gives

    λ^e_ω_n,Q_z   ≈2π∑_k_z∫_0^∞ρdρ∫_0^∞ds
     δ(ρ+ζ_ω sv^e)exp[-ρ^2/2l^2]
    
       ×exp{2ζ_ω s[ω_n+i[μ_zH-2t_z^esin (k_z)sin (12Q_z)]]}.

Finally, taking the ρ integral, we obtain the eigenvalue of the kernel 

    λ^e_ω_n,Q_z=2∫^∞_0ds⟨exp[-2ζ_ω s(ω_n+iγ̃_z^e)-v^2_es^2/2l^2]⟩_z,

where γ̃_z^e=γ_z-2t_z^esin k_zsinQ_z/2.
Substituting the Fermi velocity from Eq. (<ref>), we obtain Eq. (<ref>).



  §.§.§ Shallow band



In the case of shallow band, the quantum kernel eigenvalue 𝒥_1 is defined
by Eqs. (<ref>) and (<ref>). The summation over the Matsubara frequencies in Eq.
(<ref>) can be performed using the relation
T∑_ω_n1/(iω_n+z)(iω_n-z')
=-1/2(tanhz/2T+tanhz'/2T)/(z+z'), 
which gives



    π T∑_ω_n=-∞^∞λ_ω_n^h   =1/4∑_m=0^M_Ω∑_ℓ=0^mm!/2^m(m-ℓ)!ℓ!
       ×⟨𝒯(ℓ+γ̃_z-μ̃_h)+𝒯(m-ℓ-γ̃_z-μ̃_h)/m+1-2μ̃_h⟩ _z,

where  and the
functions μ̃_h(k_z,Q_z) and γ̃_z(k_z,Q_z) are defined by
Eqs. (<ref>) and (<ref>) correspondingly. Here, we introduced a new summation index
m=ℓ+ℓ' and cut off the diverging sum over m at M_Ω=2Ω/ω_c.
The converging function 𝒥_1 is obtained by subtracting the zero-field limit
of this sum at Q_z=0 which yields


    𝒥_1=1/4∑_m=0^M_Ω∑_ℓ=0^mm!/2^m(m-ℓ)!ℓ!⟨𝒯(ℓ+γ̃_z-μ̃_h)+𝒯(m-ℓ-γ̃_z-μ̃_h)/m+1-2μ̃_h⟩ _z
    -1/2∫ _0^M_Ωdx⟨𝒯(x-1/2-μ̃_h0)/x-2μ̃_h0⟩ _z


 with μ̃_h0(k_z)≡μ̃_h(k_z,0)=μ̅_h-2t̅_z^hcos k_z.
This quantity remains finite in the limit Ω→∞. Note that the subtracted term can also be represented as

    1/2∫ _0^M_Ωdx⟨𝒯(x-1/2-μ̃_h0)/x-2μ̃_h0⟩ _z =π T∑_0<ω_n<Ω1/ω_n+Υ_T
       ≈1/2ln𝐴Ω/T+Υ_T.

To derive a presentation better suited for numerical evaluation, we split the integration
over x as ∫ _0^∞dxF(x)=∫
_0^1/2dxF(x)+∑_m=0^∞∫_-1/2^1/2dxF(m+1+x) and subtract the
term
⟨𝒯(m/2-μ̃_h)/m+1-2μ̃_h⟩ _z from both m sums making them converging independently. This results in Eq. (<ref>) of the main text.



§ LOW-TEMPERATURE ASYMPTOTICS OF 𝒥_1(H,T,Q_Z)



In this appendix, we analyze the low-temperature behavior of the function 𝒥_1(H,T,Q_z). In most cases it behaves ∝ln(T) which gives the finite zero-temperature limit of the total kernel eigenvalue 𝒥_1-𝒜_1. This zero-temperature limiting value has square-root singularities when magnetic field crosses the miniband-edge fields [Eq. (<ref>)]. In addition,
in several resonance cases 𝒥_1 diverges as 1/√(T). Such behavior may be realized only in the uniform state with Q_z=0 or in the alternating FFLO state with Q_z=π and for special values of the spin-splitting factor γ_z.



 §.§ Derivation of the pairing kernel near the miniband-edge fields 




In this appendix, we derive the square-root contribution to the pairing kernel eigenvalue 𝒥_1, Eq. (<ref>), which appears when the chemical potential enters the Landau-level miniband at H=H_ℓ_0,σ,δ_t [see Eq. (<ref>)]. The singular behavior is coming from the terms with ℓ=ℓ_0 for σ=1 (m-ℓ=ℓ_0 for σ=-1) in the Landau-level sum, Eq. (<ref>). The corresponding singular term can be written as

    r_m,ℓ_0,σ(ω_c)=⟨tanh[ω_c(ℓ_0-σγ_z+1/2)-μ_h+2t_z^hcos(k_z-12Q_z)/2T]/m+1-2(μ_h-2t_z^hcos k_zcosQ_z/2)/ω_c⟩ _z.

At T→0 we can replace tanh(A/T) with the sign function sign(A).
For ω_c near ω_ℓ_0,σ,δ_t, we
use the presentation 

    ω_c(ℓ_0-σγ_z+1/2)-μ_h+2t_z^hcos(k_z-Q_z2)
       =(μ_h+2δ_tt_z^h)(ω_c/ω_ℓ_0,σ,δ_t-1)+2t_z^h(δ_t+cos(k_z-Q_z2)),

which indicates that the main contribution to the difference r_m,ℓ_0,σ(ω_c)-r_m,ℓ_0,σ(ω_ℓ_0,σ,δ_t)
comes from the miniband-edge region (i.e., near k_z=π+Q_z2
for δ_t=1 and k_z=Q_z2 for δ_t=-1),
where the quadratic expansion of the cosine can be used. This allows us
to evaluate

    r_m,ℓ_0,σ(ω_c)-r_m,ℓ_0,σ(ω_ℓ_0,σ,δ_t)
       ≈-2δ_t/π/m+1-2(μ_h+2δ_tt_z^hcos^2Q_z2)/ω_c
       ×√(μ_h+2δ_tt_z^h/t_z^h)√(|1-ω_c/ω_ℓ_0,σ,δ_t|)θ[δ_t(1-ω_c/ω_ℓ_0,σ,δ_t)],

where θ(x) is step function. We see that r_m,ℓ_0,σ(ω_c)
has a singular square-root behavior when ω_c approaches
ω_ℓ_0,σ,δ_t from the side, at which the
chemical potential is inside the miniband. Collecting the singular
terms, we obtain the corresponding result for the kernel eigenvalue

    𝒥_1(ω_c)-𝒥_1(ω_ℓ_0,σ,δ_t)
       ≈1/4∑_m=ℓ_0^∞m!/2^m(m-ℓ_0)!ℓ_0![r_m,ℓ_0,σ(ω_c)-r_m,ℓ_0,σ(ω_ℓ_0,σ,δ_t)].

Substituting the difference r_m,ℓ_0,σ(ω_c)-r_m,ℓ_0,σ(ω_ℓ_0,σ,δ_t) from the previous equation, we arrive at Eq. (<ref>) of the main text.




 §.§ Uniform state (Q_z=0) and  2γ_z=j_z



We analyze the low-temperature divergence of the function 𝒥_1(H,T,Q_z) in the
uniform case when the spin-splitting energy matches the Landau-level separation, i. e.,
γ_z=j_z/2. In this case, a singular behavior takes place when the chemical potential
matches the Landau level energy at k_z=0 or π corresponding to the condition
(ℓ_0+1/2j_z+1/2)ω_c=μ_h∓2t_z^h. To extract the leading
low-T divergence, we will use Eqs. (<ref>) and (<ref>). First, near T=0, the
second term of Eq. (<ref>) diverges as ∼ln T. The dominating  low-T divergent contribution
is coming from the first term with m=2ℓ_0+j_z and ℓ=ℓ_0,

    𝒥_1  ≃   ω_c/8t_z^h(2ℓ_0+j_z)!/2^2ℓ_0+j_z(ℓ_0+j_z)!ℓ_0!∫^π_0dk_z/πtanht_z^h(cos k_z∓1)/T/cos k_z∓1.

Making the substitution 
x=t_z^h/T(cos k_z∓1), we obtain


    𝒥_1≃ω_c/8π√(Tt_z^h)(2ℓ_0+j_z)!/2^2ℓ_0+j_z(ℓ_0+j_z)!ℓ_0!
    ×   ∫^2t_z^h/T_0dx
    tanh x
    	/x^3/2√(2-Tx/t_z^h).

In the limit T≪ t_z^h, using
∫_0^∞dxx^-3/2tanh x=√(2)(2√(2)-1)ζ(3/2)/√(π), 
we obtain Eq. (<ref>) in the main text.

We derive a more accurate asymptotics for the strongest resonance, ℓ_0=0 and j_z=0. 
Substituting ω_c=2(μ_h∓ 2t_z) for k_z=0/π and Q_z=0 into the Eq. (<ref>), we obtain 

    𝒥_1    =1/2∑_m=0^2Ω/ω_c∑_ℓ=0^mm!/2^m(m-ℓ)!ℓ!⟨tanhℓ-2t̅_z^h(±1-cos k_z)/2T/ω_c/m-4t̅_z^h(±1-cos k_z)⟩ _z
    
    -   1/2ln𝐴Ω/T-Υ_T.

At low temperatures we can replace tanh→ 1 in all terms except m=0. 

This gives,







    𝒥_1≃   1/2⟨tanht_z^h(1∓cos k_z)/T/4t̅_z^h(1∓cos k_z)+∑_m=1^∞[1/m+4t̅_z^h(cos k_z∓1)-1/m]⟩ _z
       + 1/2(∑_m=1^2Ω/ω_c1/m-ln𝐴Ω/ T)-Υ_T.


Using  lim_N→∞(∑_m=1^N1/m-ln N)=γ_E and performing averaging with respect to k_z, we finally obtain

    𝒥_1   ≈𝐶ω_c/√(t_z^hT)+ℛ_∓(t̅_z^h)-1/2lnω_c/π T-Υ_T,
    ℛ_∓(t̅_z^h)     =1/2∑_m=1^∞(1/√(m(m∓8t̅_z^h))-1/m),

where 𝐶=(2√(2)-1)ζ(32)/(8π^3/2)≈ 0.1072 and ζ(x) is the Riemann zeta function. Subtracting 𝒜_1 gives  Eq. (<ref>) of the main text.



 §.§ Alternating state (Q_z=π)



For the alternating case, Q_z=π, the resonance conditions are given
by Eqs. (<ref>) and (<ref>) corresponding
to t̅_z^h=1/4(ℓ_π-ℓ_0±2γ_z)
and μ̅_h=1/2(ℓ_0+ℓ_π+1).
Substituting these relations into the first term of Eq. (<ref>)
and using μ̃_h(k_z,π)=μ̅_h, γ̃_z(k_z,π)=γ_z-2t̅_z^hsin k_z,
we obtain 


    𝒥_1(H,T,π)≃1/4∑_m=0^∞∑_ℓ=0^mm!/2^m(m-ℓ)!ℓ!(m-ℓ_0+ℓ_π)
       ×⟨tanh2ℓ-ℓ_0-ℓ_π±2γ_z+(ℓ_π-ℓ_0±2γ_z)sin k_z/4T/ω_c+tanh2(m-ℓ)-ℓ_0-ℓ_π∓2γ_z-(ℓ_π-ℓ_0±2γ_z)sin k_z/4T/ω_c⟩ _z.


The singular term m=ℓ_0+ℓ_π requires resolution of
a “zero over zero” uncertainty which leads to 

    𝒥_1(H,T,π)≃ω_c/8T∑_ℓ=0^ℓ_0+ℓ_π(ℓ_0+ℓ_π)!/2^ℓ_0+ℓ_π(ℓ_0+ℓ_π-ℓ)!ℓ!
       ×⟨sech^2ℓ_0+ℓ_π-2ℓ∓2γ_z-(ℓ_π-ℓ_0±2γ_z)sin k_z/4T/ω_c⟩ _z.

The terms for which |ℓ_0+ℓ_π-2ℓ∓2γ_z|=|ℓ_π-ℓ_0±2γ_z|
are divergent. This condition is always satisfied for ℓ=ℓ_π
and this term is

    𝒥_1(H,T,π)≃ω_c/16π T(ℓ_0+ℓ_π)!/2^ℓ_0+ℓ_πℓ_0!ℓ_π!∫_-π^πdk_zsech^2t_z^h(1+sin k_z)/T.

At low temperatures, T≪ t_z^h, the dominating contribution
to the integral comes from the region near the inflection point k_z=-π/2
and the integration can be approximately evaluated using the substitution
k_z=-π/2+√(2Tx/t_z^h), giving 

    ∫_-π^πdk_zsech^2t_z^h(1+sin k_z)/T   ≈√(2T/t_z^h)∫_0^∞dx/√(x)cosh^2x
       =√(2T/t_z^h)(2√(2)-1)ζ(3/2)/√(2π).

Therefore, we obtain

    𝒥_1(H,T,π)≃𝐶/2(ℓ_0+ℓ_π)!/2^ℓ_0+ℓ_πℓ_0!ℓ_π!ω_c/√(Tt_z^h),

where the constant 𝐶 is defined after Eq. (<ref>). We also mention that in the exceptional cases when 2γ_z is integer,  the term with ℓ=ℓ_0-2γ_z also diverges and this gives additional diverging contribution

    𝒥_1(H,T,π)    ≃𝐶/2(ℓ_0+ℓ_π)!/2^ℓ_0+ℓ_π(ℓ_π+2γ_z)!(ℓ_0-2γ_z)!ω_c/√(Tt_z^h).


We again derive a more accurate asymptotic for the strongest resonance of this kind with t_z^h/μ_h=γ_z and ω_c=2μ_h. Near the first Lifshitz transition, μ_h>2t_z^h which implies that γ_z<0.5. For these relations, the function 𝒥_1 becomes

    𝒥_1  =1/4∑_m=0^M_Ω∑_ℓ=0^mm!/2^m(m-ℓ)!ℓ!m⟨tanhℓ+γ_z(1+sin k_z)/2T/ω_c.
       .+tanhm-ℓ-γ_z(1+sin k_z)/2T/ω_c⟩ _z -1/2ln𝐴Ω/T-Υ_T.

Note that the m=0 term again has “zero over zero” uncertainty. 
All terms in the sum remain finite for T→0
except the m=ℓ=0 term. Therefore, we separate this diverging
term from the sum and take the limit T→0 for the rest of the
terms. This yields 





    𝒥_1≃   ω_c/8T⟨sech^2γ_z(1+sin k_z)/2T/ω_c⟩ _z
       -1/2∑_m=1^∞1/2^mm+1/2∑_m=1^M_Ω1/m-1/2ln𝐴Ω/ T-Υ_T
    ≈   𝐶/2ω_c/√(Tt_z^h)-1/2ln2ω_c/π T-Υ_T.

This corresponds to Eq. (<ref>) in the main text.




§ THE PHASE DIAGRAMS FOR THE DEEP-BAND DOMINATING SCENARIO




In the main text, we mostly presented results for the interband-coupling scenario (Λ_ee=Λ_hh=0), for which the shallow band is essential for the formation of the superconducting state and the quantization effects caused by this band are very pronounced. In this appendix, we consider the phase diagrams for the alternative scenario when the deep band itself has strong pairing strength and the interband coupling induces superconductivity into the shallow band, i. e., Λ_ee>Λ_he,Λ_eh,Λ_hh. In such scenario, the Lifshitz transition weakly affects superconductivity as is indeed realized in several iron-based superconductors. On the other hand, the influence of the shallow band is only noticeable if it has substantial contribution to pairing, i. e., the interband coupling constants  Λ_he and Λ_eh should be comparable with Λ_ee. Note that we focus in this paper on the scenario when there is no intraband pairing in the shallow band, Λ_hh=0. Such pairing would further enhance shallow-band quantization effects. For illustration, we use the coupling constants Λ_ee=0.3 and Λ_he=Λ_eh=0.2 in this appendix. This corresponds to the ratio Λ_he/Λ_0,e∼ 0.54 in Fig. <ref>.

We consider the resonant case first. Figure <ref> shows the representative H-T diagrams for γ_z=0.5. The diagrams are presented for three values of the interlayer hopping energies t_z^h shown in the plot. We can see that the overall behavior is similar to the interband-coupling scenario shown in Fig. <ref>. The quantitative difference is that the size of the reentrant region is rather small even in the 2D case and it further diminishes with increasing the interlayer tunneling. Already for small interlayer hopping energy, t_z^h=0.1T_C, the maximum T_C2 drops below 0.02 T_C.

In a more generic nonresonant case (noninteger 2γ_z), two features of the H-T diagrams have been emphasized in the main text: (i)the interlayer FFLO instability accompanied by the enhancement of the upper critical field at low temperatures and (ii)the emergence of reentrant states for the magnetic fields slightly below the crossings of two miniband-edge field lines.
We find that the FFLO instability in the main superconducting region is robust. For example, Fig. <ref>(a) shows the magnetic field – temperature phase diagram for the same electronic parameters as in Fig. <ref>(a) but for large pairing strength in the deep band 
(the same coupling constants Λ_ee=0.3, Λ_he=Λ_eh=0.2, and Λ_hh=0).
For these parameters, the shallow band induces the FFLO instability at T≈ 0.1T_C accompanied by the upturn of the upper critical field at lower temperatures. This upturn is much smaller than the one for the interband-coupling scenario presented in Fig. <ref>(a). It is enhanced, however, when the miniband-edge field approaches the zero-temperature H_C2. This can be seen from the inset which also shows the low-temperature behavior of the H_C2 line for the smaller Fermi energy μ_h=0.9T_C, where
the H_C2 line approaches the miniband-edge field “2↓+t_z” at T→ 0.
Such enhancement is only observed for few special values of the Fermi energy, while the behavior illustrated in the main figure is typical. Another difference from the interband-coupling scenario is that the modulation wave vector remains rather small, Q_op<1, much smaller than the optimal wave vector favored by the shallow band. 

We find that in the dominating deep-band regime the reentrant behavior is not realized for the electronic parameters used in Fig. <ref>. Nevertheless, this behavior does appear for smaller interlayer hopping energy and/or smaller Fermi energy of the deep band. For example, Fig. <ref>(b) presents of the phase diagram for t_z^h=0.3T_C and μ_h=2.05T_C. In this case, the reentrant region connected with the main domain exists. As one can see from the inset, it separates at slightly larger μ_h. This region is caused by the closely located miniband-edge fields for ℓ=3/spin down/k_z=π and ℓ=1/spin up/k_z=0. Such reentrance only exists within narrow range of the Fermi energies.


For better representation of the overall behavior, we show in Fig. <ref> the dependences of the low-temperature upper critical field on the Fermi energy μ_h. 
The lower curve in the lower plot is made for parameters t_z^h=0.45T_C and ε_0=12T_C (first set). The same plot is also shown in Fig. <ref>.  There are no reentrant regions for this parameter set. To illustrate the emergence of such regions and other general trends, we also present H_C2(μ_h) dependences for (i) smaller energy separation between the band edges  ε_0=8T_C (the upper curve in the lower plot, second set) and (ii) smaller interlayer hopping energy  t_z^h=0.3T_C (upper plot, third set).
Note that the parameter ε_0 controls the Fermi energy of the deep band μ=ε_0-μ_h.
In all cases, the H_C2(μ_h) curves have pronounced oscillations and have a clear tendency to follow the miniband-edge field lines. As we already mentioned, for the first set H_C2(μ_h) is a single-valued function reconfirming the absence of the reentrant behavior. We see, however, that the H_C2(μ_h) line has the pronounced kinklike features near the crossing of the miniband-edge lines. For the Fermi energies at the maxima, the H_C2(T) line has a pronounced low-temperature upturn, as illustrated in the inset of Fig. <ref> for μ_h=0.9T_C. For the second and third sets multiple solutions for H_C2(μ_h) appear below some crossing points corresponding to existence of the reentrant behavior in these regions. For example, for the t_z^h=0.3T_C (upper plot), such regions are located near μ_h=1.28 and 2.1T_C.  We also see that for these sets a kinklike feature also exists below the first Lifshitz transition at μ_h=0.9T_C, in the region of closed Fermi surface, where the H_C2 curve separates from the “1↓+t_z” line. It is clear that the reentrant regions will become more pronounced with further decreasing of either t_z^h or ε_0. 




