Orbital Fulde-Ferrell-Larkin-Ovchinnikov state in an Ising superconductor

The critical field behavior of a layered Ising superconductor with finite number of layers is studied. Under in-plane magnetic fields, the finite-momentum superconductivity dubbed as the orbital Fulde- Ferrell-Larkin-Ovchinnikov state is found in the regime of low field and high temperature. Our theory agrees well with the experimental results in Nature 619, 46 (2023).

Introduction-The dimensionality of a superconductor is usually derived from its critical field behavior, which reflects the spatial profile of the order parameter under external magnetic fields.When the superconductor is three-dimensional (3D), Abrikosov vortices [1] will be formed under fields above the lower critical field, where the magnitude of the order parameter is highly nonuniform and the phase has an integer winding around a vortex core.As a result, the upper critical field of a 3D superconductor is linear in temperature, and the critical exponent is independent of the field direction.
On the contrary, in a two-dimensional (2D) superconductor the order parameter is mostly uniform in magnitude and Abrikosov vortices cannot be formed under inplane fields.Since the characteristic size of the Abrikosov vortex core is the coherence length ξ, the criterion of 2D superconductivity is roughly d < ξ, where d is the thickness.Detailed calculations further reveal the critical thickness for 2D superconductivity d < d c ≈ 1.8ξ [2].In a 2D superconductor, near the zero-field critical temperature, the out-of-plane upper critical field is still linear in temperature, while the in-plane one is square-root in temperature [2,3], as verified in experiments [4][5][6][7].
The above arguments on dimensionality are based on continuum models of superconductors.Over the past several decades, the layered superconductors have been discovered [2,[7][8][9], where each layer is an atomically thin 2D superconductor, and different layers are weakly coupled by Josephson coupling [10][11][12][13][14].For a layered superconductor with infinite number of layers, a dimensional crossover can be realized by an in-plane magnetic field [2,9,[15][16][17].When the in-plane field is weak compared with the interlayer Josephson coupling, the layered superconductor can be treated as 3D with an upper critical field linear in temperature.As the field increases so that the interlayer coupling is relatively negligible, the layered superconductor behaves as if a 2D superconductor with an upturn in the upper critical field.Such a dimensional crossover has been experimentally observed in several layered superconductors [16][17][18][19].
Recently, it was experimentally found [20] that a layered Ising superconductor NbSe 2 with intermediate thickness can have unconventional critical field behavior.The upper critical field is square-root in temperature and hence 2D at low fields.As the field increases, an additional phase transition associated with a tricritical point is found instead of a dimensional crossover from 3D to 2D.These results are inconsistent with the critical field behavior of the layered superconductor with infinite layers [2,[15][16][17], but show similarities to bilayer superconductors [21,22] instead.In Refs.[21,22], it is found that in a bilayer superconductor linked by Josephson coupling, unconventional superconducting states similar to the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) states [23,24] can be realized under in-plane magnetic fields mainly through the orbital effect.
In this manuscript, we would like to study the layered superconductors with intermediate number of layers under in-plane fields, and try to generalize the FFLO-like states in Refs.[21,22] to multilayer cases.
It is well known that two depairing mechanisms are introduced by an external magnetic field, namely the orbital effect and the Zeeman effect.In a conventional superconductor, the Zeeman effect can be neglected at weak fields far below the Pauli limit B P ≡ ∆ 0 / √ 2µ [16,17,25], where ∆ 0 is the pairing gap at zero temperature and µ is the magnetic moment.In an Ising superconductor, the Zeeman effect is screened [5][6][7] for in-plane fields far below the Ising limit B so ≡ ∆ so /µ, where ∆ so is the Ising spin-orbit coupling (SOC) energy.As the materials we are interested in (e.g.Ising superconductors of transition metal dichalcogendies) are composed of bilayer unit cells, in the following we consider only even number of layers if not specified otherwise.
Model -We consider the N -layer Lawrence-Doniach (LD) model with the following free energy density per area [15,21,22] where m > 0 is the electron mass, e the electron charge, and J > 0 is the Josephson coupling between nearest neighbor layers.In the gradient terms, ∇ ∥ = (∂ x , ∂ y ) is the in-plane gradient operator, A l is the in-plane vector potential on layer l, and A z is the out-of-plane component of vector potential.The second order coefficient can be worked out as a function of temperature T and field (2) where N 0 is the density of states of each layer, the special function F (x) = Re Ψ( 12 ) − Ψ 1 2 (1 + ix) is defined in terms of the digamma function Ψ(x), T c1 is the critical temperature of the monolayer.The fourth order coefficient β > 0 can be treated as a positive constant independent of temperature and field.When B ≪ B so , the field-dependent term in Eq. ( 2) can be neglected, and we can write α = α 0 (T − T c1 ) when T ≲ T c1 , where α 0 = N 0 /T c1 .In other words, the Zeeman effect of inplane fields far below Ising limit can be neglected in a layered Ising superconductor.Now we consider finite number of layers.We employ the 2D ansatz ψ l = ∆ l e iq•r , whose magnitude is in-plane uniform |ψ l (r)| ≡ ∆ l , and the phase is characterized by Cooper pair momentum q.Correspondingly we choose the gauge Hence the free energy is f = f ({∆ l }, q, B).Then the upper critical field of the second order superconducting phase transition is given by the highest critical field among all q and order parameter profile {∆ l }, and the corresponding q is the optimal Cooper pair momentum.
In the following, we mainly focus on the upper critical field regime, hence B usually denotes the in-plane upper critical field if not specified otherwise.
Results-The upper critical field and optimal Cooper pair momentum are numerically calculated and plotted in Fig. 1a and b respectively, normalized by the three units of temperature, field and momentum respectively where Φ 0 = h/(2e) is the flux quantum.At zero field, superconductivity occurs at temperatures below the layer-dependent critical temperature T cN as plotted in Fig. 2a.Similar results can be found in Ref. [26], while the mechanism is interlayer Cooper pairs instead of interlayer Josephson coupling.The zero-field critical temperature data of few-layer NbSe 2 can be found in Ref. [6], which agree well with Eq. (B5) as shown in the inset of Fig. 2a.This result can be analytically obtained as shown in the next session.
As shown in Fig. 1b, at low fields, the Cooper pair momentum remains zero, and the layered superconductor behaves as a 2D superconductor with square-root temperature-dependence of the upper critical field as shown in Fig. 1a (dahsed lines).When the field further increases, the optimal Cooper pair momentum becomes finite and along ±B × ẑ directions, which applies to all even numbers of N ≥ 2. We define q ≡ q • ( B × ẑ), which is plotted in Fig. 1b as a function of field strength B, implying a phase transition with a tricritical point (T * N , B * N ).When B < B * N , the Cooper pair momentum is zero q = 0, When B > B * N , q ̸ = 0, and at the tricritical point (T * , B * ), two types of superconducting phases q = 0 and q ̸ = 0 coexist with the normal phase.The numerical tricritical field B * N is plotted for different N in Fig. 2b.In our finite-layer model, inversion symmetry is found to hold for the free energy, under which layer l with momentum q maps to layer N + 1 − l with momentum −q.Thus states with ±q are degenerate in calculating the upper critical field, and in Fig. 1b we only plot the non-positive branch of q.
The superconducting phase with finite-momentum Cooper pairs is not rare at least in theoretical studies.In 1964, Fulde and Ferrell [23], Larkin and Ovchinnikov [24] proposed that due to Zeeman effect, magnetic field can drive Cooper pairs into finite-momentum states.Since the energy saved by finite-momentum Cooper pairs per unit field is small, conventional FFLO states are expected at low temperatures and high fields, which turn out not easy to be detected experimentally.However, in our theory, the finite-momentum Cooper pairs are boosted by magnetic field via orbital effect, which can survive at relatively low field and high temperatures.We may dub such states as the orbital FFLO states.
Analysis-To figure out the origin and mechanism of the orbital FFLO states, we analyze the order parameter profiles near superconducting phase transitions.
We first review the inifinite layer limit N → ∞.In this limit, due to the translation symmetry l → l + 1, the order parameter at zero field is uniform, and the bulk critical temperature is T cN → T c ≡ T c1 + 2T 0 .At finite in-plane fields, Abrikosov vortices can be formed, and the bulk in-plane upper critical field is When N is finite, the translation symmetry is lost, and the middle layers become more favorable than other layers because of the finite size effect.Under an in-plane field, Abrikosov vortices become metastable, and Josephson vortices formed by middle layers are the stable state, which are the orbital FFLO states we found previously.
To be concrete, the order parameter at zero field reads which is nonuniform, and the middle two layers l = 1 2 N and l = 1 2 N + 1 have the highest weightings.From Eq. (B2), the critical temperatures at zero and weak fields are Detailed derivations of the order parameter Eq. (B2), the critical temperatures Eqs.(B5,B8), and the expression of the dimensionless coefficient c N can be found in the Appendix.Importantly, the quadratic dependence of the field in Eq. (B8) implies the 2D characteristic of the  10).Black dots are experimental data of upper critical field in Ref. [20], while orange and blue lines are analytical results.BCS state q = 0, orbital FFLO state q ̸ = 0 and normal phase coexist at the tricritical point (T * , B * ).
in-plane upper critical field in the weak field regime, as shown in the dashed lines of Fig. 1a and found experimentally in Ref. [20].
As shown in Fig. 3a, at zero field, the order parameter profile Eq. (B2) is peaked in the middle between layer l = 1 2 N and l = 1 2 N + 1, with a broad width.At weak field, the order parameter is still peaked at in the middle between layer l = 1 2 N and l = 1 2 N +1, but with a smaller width.As field increases, the width of order parameter keeps shrinking.When the field is strong enough, the width of order parameter becomes small enough, Cooper pairs tend to accumulate on layer l = 1 2 N or l = 1 2 N + 1.To be more precise, the order parameter profile at intermediate in-plane fields is a Gaussian profile [2] with field-dependent width l B = B 0 /B and center At zero field, the order parameter can be characterized by the typical width δl N ∼ √ N + 1/π (see Appendix).The order parameter becomes sufficiently localized when l B ≲ δl N , or equivalently when B ≳ B * N , with the tricritical field and tricritical temperature T * N ≈ T c − 1.6πT 0 /(N + 1).The critical behavior near the tricritical point (T * , B * ) is governed by an effective bilayer model of two localized modes ψ ± on layer 1  2 N and 1 2 N + 1 respectively.One then expands the free energy density as f = Ψ † HΨ + . . . on the basis Ψ = (ψ + , ψ − ) T , and up to q 2 and B 2 : where τ denote Pauli matrices acting on the Ψ basis, a = r(T −T a )+χB 2 , and r, χ, b, c, J , T a are phenomenonlogical parameters.For the stability of the free energy, we require r, χ, c to be positive.The q, B coupling term is due to the orbital effect, while the magnetic energy term χB 2 can be induced by both orbital and Zeeman effects of the in-plane magnetic field.At zero field, the critical temperature of the effective bilayer model Eq. ( 10) is T c = T a + J /r.Under an in-plane field, the optimal Cooper pair momentum is with the tricritical field The field-dependence of optimal Cooper pair momentum is plotted in Fig. 1b (dashed lines), which agrees well with numerical simulations (dots).Correspondingly, the upper critical field as a function of temperature is plotted in Fig. 3b (orange and blue lines), which shows a kink at the tricritical point, and agrees well with experimental data from Ref. [20] (black dots).Details of the fitting parameters can be found in Appendix.
Symmetry principles also allow us to write down higher order terms.In transition metal dichalcogenide such as NbSe 2 , the point group is D 3d , and the sixth order terms are H 1 (q, B) = λ 1 Re(q 2 + B 4 + ) + λ 2 Re(q 4 + B 2 + ) + λ 3 Im(q + B 2 + )τ y , where q ± = q x ±iq y , B ± = B x ±iB y , and the total Hessian matrix is H = H 0 + H 1 .These higher order terms leads to the anisotropy of the field-dependent critical temperature in the orbital FFLO phase with q ̸ = 0, while the field-dependent critical temperature in the conventional superconducting phase is isotropic since q = 0. To be concrete, by plugging Eq. (C7) into H 1 we can calculate the anisotropic part of the field-dependent critical temperature ∆T c (B, φ) = λ(B)θ(B − B * ) cos(6φ), where B = B(cos φ, sin φ), θ is the Heaviside theta function, and λ(B) is given in Appendix.
The leading anisotropy is sixfold due to the following symmetry reasons.Since there is no intrinsic timereversal symmetry breaking in our model, T c (−B) = T c (B).The threefold in-plane rotation symmetry C 3z of TMD implies T c (C 3z B) = T c (B), we derive that the leading anisotropy of the orbital FFLO states in TMD is sixfold T c (C 6z B) = T c (B) as shown in experiments [20].
The stable combination of two localized modes ψ ± will be determined by the quartic order free energy with coefficients β 1,2 , and β 1 > 0, β 2 > −2β 1 for the stability of the free energy (see Appendix).When β 2 > 2β 1 the stable state will break the inversion symmetry and choose one of the localized modes of ψ ± .When β 2 < 2β 1 the stable state will preserve the inversion symmetry to form a symmetric combination of ψ ± , while the in-plane translation symmetry along q ∝ B × ẑ direction is spontaneously broken.In the decoupling limit J → 0, the equilibrium state will preserve the inversion symmetry to compensate the orbital effect, hence we expect β 2 < 2β 1 generally holds when the Josephson coupling is not too strong.The orbital FFLO state would be a Josephson vortex formed by localized modes on middle layers.

Conclusion-
In this work, we generalize the orbital FFLO states theoretically derived in the bilayer superconductors [21,22] to multilayers.We find that due to the finite size effect, the middle layers become more favorable than other layers.Under an in-plane field, the multilayer Ising superconductor can be described by an effective bilayer model, with a field-driven phase transition from the conventional pairing state to the orbital FFLO state.Our theory can be applied in transition metal dichalcogenide layers such as NbSe 2 [6,20,27].

FIG. 1 :
FIG. 1: Orbital FFLO phases in N -layer superconductors.a) Upper critical field as a function of temperature.With the increasing zero-field critical temperature, even number N changes from 2 to 20, labeled by different colors.Solid lines are numerical results, dashed lines denote the analytical formula Eq. (B8).b) Upper critical field versus Cooper pair momentum.Colors denote N , the same as a).Dashed lines denote the analytical formula Eq. (C7).As N increases, upper critical field approaches the envelope B = B0(Tc − T )/T0 in a) and optimal Cooper pair momentum approaches the envelope q/q0 = − 1 2 B/B0 in b).Here B0, T0, q0 are defined in Eq. (A2).

FIG. 2
FIG. 2: a) Zero-field critical temperature TcN as a function of layer number N .Dots are from numerical simulations as in Fig. 1a and the black line is Eq.(B5).Inset: Fittings (lines) of experimental data (markers) in Ref. [6] by Eq. (B5).Different colors denote Tc values according to different criteria, and Rn is the normal resistance.b) Tricritical field B * N as a function of layer number N .Dots are from numerical simulations as in Fig. 1b and the black line is Eq.(9).

FIG. 3
FIG. 3: a) Order parameters at different in-plane fields with layer number N = 20.Open circles are numerical results of Eq. (A1) and solid lines are analytical formulae Eqs.(B2,7).b) Phase diagram of the effective model Eq.(10).Black dots are experimental data of upper critical field in Ref.[20], while orange and blue lines are analytical results.BCS state q = 0, orbital FFLO state q ̸ = 0 and normal phase coexist at the tricritical point (T * , B * ).