Superconductivity of surfaces , edges and vertices in systems with imbalanced fermions

We describe boundary effects in superconducting systems with Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) superconducting instability, using Bogoliubov-de-Gennes and Ginzburg-Landau formalisms. Firstly, we show that in dimensions larger than one the standard Ginzburg-Landau (GL) functional formalism for FFLO superconductors is unbounded from below. This is demonstrated by finding solutions with zero Laplacian terms near boundaries. We generalize the GL formalism for these systems by retaining higher order terms. Next, we demonstrate that a cubic superconductor with imbalanced fermions, at a mean-field level has a sequence of the phase transitions. At low temperatures it forms Larkin-Ovchinnikov state in the bulk but has a different modulation pattern close to the boundaries. When temperature is increased the first phase transition occurs when the bulk of the material becomes normal while the surfaces remain superconducting. The second transition occurs at higher temperature where the system retains superconductivity on the edges. The third transition is associated with the loss of edge-superconductivity while retaining superconducting gap in the vertices. We obtain the same sequence of phase transition by numerically solving the Bogoliubov-de Gennes model.


I. INTRODUCTION
Fulde and Ferrell [1] and Larkin and Ovchinnikov [2] (FFLO) considered a superconducting state where, Cooper pair forms out of two electrons with different magnitude of momenta. The original considerations expected such a situation to arise in the presence of a strong magnetic field and thus Zeeman splitting of the Fermi surfaces for spin up and spin down electrons. Later it was shown that in other physical systems the fermionic imbalance occurs without any applied magnetic field. Indeed the FFLO state was discussed in cold-atom gases where one can create different imbalances of fermions on demand [3][4][5][6]. Similarly the difference in Fermi surfaces naturally occurs in dense quark matter. The resulting superconducting states of quarks are called color superconductivity which is suggested to realize FFLO state in the cores of neutron stars [7]. Even in electronic superconductors, finite momentum pairing may arise due to reasons other than application of an external magnetic field [8]. This state has for a long time been of great interest and was searched for in a number of superconducting materials [9][10][11][12][13][14][15][16][17][18][19].
In the recent work [20], using microscopically derived Ginzburg-Landau model, it was shown that systems that support FFLO superconductivity in the bulk do not undergo a direct superconductor-normal metal phase transition. Instead these systems have a different intermediate phase at elevated temperatures where superconductivity occurs only on the surface but the bulk of the systems is a normal metal.
In this work we answer three questions. (i) We demonstrate this effect in microscopic models for imbalanced fermions without relying on a Ginzburg-Landau expansion. (ii) For a Ginzburg-Landau approarch, we explain that, as alluded in [20], in dimensions larger than one, it is necessary to retain more terms in the Ginzburg-Landau expansion than what is done in the standard calculations [21] to describe systems with boundaries. (iii) We show that a three dimensional cubic superconductor with imbalanced fermions, undergoes a sequence of phase transitions where superconductivity survives at sub-domains of sequentially lower dimensions. That is, by increasing the temperature and fermionic population imbalance in a cubic system, at the meanfield level, the system will undergo the following sequence of the phase transitions: superconducting bulk → superconducting surfaces → superconducting edges → superconducting corners → normal state.
The plan of the paper is the following: In Section II we firstly recap canonical microscopic derivation of the Ginzburg-Landau model for non uniform FFLO superconductors. Then we demonstrate that the usual GL model is unbounded from below because it is unstable to formation of infinitely strong gradients near boundaries. In Section III we construct a Ginzburg-Landau functional which has energy bounded from below by retaining additional terms in the expansion. That allows us to in Section IV, obtain the surface, edges and corner states without divergent energies.
In Section V we solve numerically Bogoluibov-de-Gennes equation to show that the states exist in microscopic models that does not rely on any Ginzburg-Landau expansion.

A. Ginzburg-Landau Model
The Ginzburg-Landau description of superconductors in the presence of fermionic population imbalance was derived from microscopic theory for superconductors in [21]. The free energy functional reads F [ψ] = Ω Fd d x arXiv:1905.03774v1 [cond-mat.supr-con] 9 May 2019 where the free energy density F is where ψ is a complex field referred to as the superconducting order parameter and c.c. denotes complex conjugation. The coefficients α, γ, and ν depend on the fermionic population imbalance H and temperature T accordingly [21] where N (0) is the electron density of states at the Fermi surface, T c is the critical temperature at zero H and we have defined the functions where z = 1/2 − iH/2πT and Ψ (n) is the polygamma function of order n. The remaining coefficients are given as β =βv 2 F γ, δ =δv 4 F ν, and µ =μv 2 F ν, where v F is the Fermi velocity andβ,δ,μ are positive constants that depend on the dimensionality d. In one dimension we havê β = 1,δ = 1/2, andμ = 4, in two dimensions we havê β = 1/2,δ = 3/16, andμ = 2, and in three dimensions we haveβ = 1/3,δ = 1/10 andμ = 4/3. The Ginzburg-Landau description is valid in the parameter regime in which the highest order terms are positive, that is where ν is positive. In the parameter regime in which β is negative, inhomogeneous order parameters may be energetically favorable. Typically considered structures of the order parameter are the so-called Fulde-Ferrell (FF) state ψ FF ∝ e ipx and the Larkin-Ovchinnikov (LO) state ψ LO ∝ cos px, with p 2 = −β/2δ and transition into the normal state at α = α bulk c = β 2 /4δ. Inhomogeneous states can appear when the term |∇ψ| 2 has a negative prefactor. In this parameter regime, it is necessary to include higher order terms, resulting in a free energy density expansion in Eq. (1). Where term |∇ψ| 2 favors creation of the gradients, while the positive term |∆ψ| 2 is added to bound gradients from above. However we will show that in some cases, the inclusion of the stabilizing term |∆ψ| 2 is not sufficient. In small systems and generically in finite systems in two dimensions or higher, there exist solutions that satisfy ∆ψ = 0. These states are often characterized by enhancement of ψ close to the surface. The associated energy of such state can potentially diverge and one needs to resort to a more general Ginzburg-Landau theory.

B. Small systems
Consider first a one dimensional domain Ω = [−L/2, L/2] of length L. The equation ∆ψ = 0 is satisfied by the real field ψ(x) = qx, for any parameter q. The modulus |ψ| increases linearly as the boundaries are approached, resulting in increasing potential energy density closer to the boundaries. As system size increases, the penalizing potential energy is non-negligible, making the linear solution less energetically beneficial. Therefore, the linear solution is not energetically preferable over conventional one dimensional inhomogeneous structures, such as ψ LO ∝ cos px or ψ FF ∝ e ipx if L 1/p, where p 2 = −β/2δ. However, for small system sizes, the potential term does become negligible in comparison to the beneficial gradient term β|∇ψ| 2 = βq 2 . To lowest order in L, minimizing with respect to q 2 gives with transition into the normal state at α = α c1 = −12β/L 2 . In the limit of infinitesimal system size L → 0, we find that α c1 → ∞, and the associated momentum q and energy F diverges. Note that the value of ψ are the surfaces, ψ(±L/2) = ±qL/2, does not vanish in this limit since q ∝ 1/L. Consequently the value of ψ at the boundary is independent of system size for sufficiently small systems.
The effect generalizes to systems in higher dimensions. Consider the d-dimensional cube with volume L d . The linear solution here generalizes to a multi-linear solution ψ n (x) = q n n j=1 x j , where n ≤ d. Analogously to the one-dimensional case we find with transition into the normal state at α cn = nα c1 . Consequently we expect that in higher dimensions, we will see a sequence of transitions from (n = 1) → (n = 2) → . . . → (n = d) before transitioning into the normal state, as seen in Figure 1 in three dimensions. Studying the energy F n , we see that the free energy would approach some negative constant in two dimensions, and zero from below in three dimensions, as the system size approaches zero. Similarly to the one-dimensional case, the multi-linear solution ψ n is constant and independent of system size at the vertices of the n-dimensional cube, for significantly small systems. x j in three dimensions, in small system sizes. As α is increased, we transition from the linear solution n = 1, to the bilinear solution n = 2 and finally the tri-linear solution n = 3 before entering the normal state.

C. Surface state instability in Ginzburg-Landau model
Consider a two-dimensional domain Ω in the presence of boundaries. We can represent the domain Ω as section of the complex plane, by defining a complex coordinate ζ = x + iy. The order parameter ψ is now a function of ζ. Clearly if ψ is an analytic complex valued function (i.e. it satisfies the Cauchy-Riemann equations), it follows automatically that ∆ψ = 0. The system becomes unstable towards states of the form ψ(ζ) = Ae kζ with |k| → ∞, associated with a negative divergent energy. Crucially, the formations of these states require the existence of a boundary. If the complete two-dimensional plane were considered, the solution ψ(ζ) = Ae kζ would diverge in some direction, which would be penalized by the potential terms. However, the divergence of |ψ| can be cut off using the boundary of the system, as seen in Figure 2.
For example, consider the half-infinite two-dimensional domain x ≥ 0. An analytic function which satisfies ∆ψ = 0 is ψ FF = Ae −kx+iky with k > 0. This state is characterized by phase-wave modulation (FF state) tangential to the boundary and exponential decay perpendicular to the boundary. We compute the energy density along the boundary (that is in the y-direction) by integrating over x and find We see that as long as the amplitude |A| 2 < −2β/µ, the energy diverges to −∞ as k increases. This shows that the model has an instability towards formation of surface states with rapid modulation along the boundary.
Analogous calculation for the state ψ LO = Ae −kx cos ky would yield the same conclusion, see Figure 2. Figure 2: Order parameter of two-dimensional surface state with density modulation along the boundary. In an infinite system, such a state would not be favorable due to divergent modulus in one direction. However, in the presence of boundaries, the segment of divergent modulus can be cut off, depicted on the image by the transparent surface.

III. GENERALISATION OF THE GINZBURG-LANDAU MODEL
In contrast to the example of linear solutions in small systems, see Section II B, the two-dimensional solutions with modulation along the boundary and exponential decay from the boundary exist in any system in the presence of boundaries, regardless of the system size. The divergence in energy and Cooper pair momentum is indicative towards the existence of superconducting surface states. However, in order to find a suitable theory which is bounded from below, we have to include additional terms in the Ginzburg-Landau expansion.
We would like to find the lowest order term in both momentum and amplitude of the order parameter, which is sufficient to include in order to stabilize the theory. The origin of the inhomogeneous state is the negative term |∇ψ| 2 , which is of second order in both momentum and amplitude. Conventionally it is stabilized by the term |∆ψ| 2 , which is of fourth order in momentum and second order in amplitude. Since there exists states that satisfy ∆ψ = 0, we have to resort to higher order in either momentum or amplitude. Terms such as |∇∆ψ| 2 , which is of sixth order in momentum and second order in amplitude, would not be helpful since it is automatically zero if ∆ψ = 0. Therefore we resort to introducing the term which is of fourth order in momentum and fourth order in amplitude, which reads Since this term is of higher order in amplitude than the lowest order beneficial gradient term, there will exist inhomogeneous superconducting states for any α, as long as the gradient coefficients are negative. This term was used to find two-dimensional solutions in [20]. Here we derive from the microscopic theory the following estimate for the coefficient κ in Eq. (9) κ − 29 32 where factor Ω d depends on the dimension d, where Ω 1 = 1, Ω 2 = 3/8 and Ω 3 = 1/5. Studying the functions K 5 (H, T ) and K 7 (H, T ), we see that there exist an overlapping region in which both functions are negative, which implies that the additional term is positive simultaneously as the previously higher order terms. In order to study a minimal model, we proceed with only retaining the term in Eq. (9) in the regularized free energy expansion, even though additional terms also proportional to K 7 (H, T ) could be included. For example, in principle we could have included a potential term proportional to |ψ| 8 . However, since the sixth order potential term remains non-zero and positive, the inclusion of this term is not so important.
Let us now consider the previous example of the surface state ψ FF = Ae −kx+iky , which previously was associated with negative infinite energy and infinite momentum k.
With the inclusion of the term in Eq. (9), the associated free energy in Eq. (8) obtains the additional term which makes the free energy bounded from below and the optimal value of k is now finite.

IV. SEQUENTIAL PHASE TRANSITIONS IN THREE DIMENSIONS
Having obtained the Ginzburg-Landau free energy expansion without the spurious divergences enables us to study higher dimensional systems. We will numerically minimize the free energy functional using the nonlinear conjugate gradient method, parallelized on the CUDAenabled NVIDIA graphical processing unit (GPU). It is convenient to introduce rescaled coordinates and parameters defined accordingly:ψ = ψ/|ψ U |,α = α/α U , x = px, where |ψ U | 2 = −γ/2ν, α U = γ 2 /4ν and p 2 = −β/2δ. Expressed in these quantities, the free energy reads where the rescaled free energy densityF is identical to Eq. (1), but the coefficients have been replaced accordingly: α →α, β →β, and so on, whereγ = −2ν = −2, β = −2δ = −2β 2 /δ, andμ =βμ/δ. Among these coefficients, all are constant exceptα, which parametrizes both temperature T and fermionic population imbalance H. With the inclusion of the additional term in Eq. (9), in rescaled coordinates, the coefficientκ reads Specifically in three dimensions we have thatκ = 145K 3 K 7 /18K 2 5 , which is not constant in the rescaled units. However, studying the functions K n in Eq. (5), we see that K n ∝ T /T n · f n (H/T ), where f n is some elementary function. Therefore, if we study a line in the T H-plane where H/T is constant, the rescaled parameterκ will also be constant. Particularly, we will study the line where H/T = 2π/3, along whichκ 0.5752. On a technical note, since we are working in rescaled units, where we measure length in units of p −1 ∝ K 5 /K 3 , we also have to alter the rescaled length accordingly in order to appropriately describe the sample of fixed size.
The obtained solutions are shown in Figure 3. We see that as temperature is increased, there are multiple phase transitions. The dimensionality of the transition sequentially decreases. At the lowest temperature, the system is in a superconducting state with a non-uniform order parameter in the bulk. Note that deep in the bulk the solution is of Larkin-Ovchinnikov type. On the surface, due to discussed above reasons we expect larger gradients. This is indeed clearly seen on the upper panel of Fig.  3. Note, that the order parameter configurations on the surfaces are very different from the Larkin-Ovchinnikov solutions in the bulk.
As temperature is increased, we first transition into a superconducting state which allows for modulation on the surface of the cube, while the bulk now has transition into a normal state. Second, when the temperature is elevated further, the surfaces become normal too, but superconductivity survives on the edges of the cube. Thirdly, the edges become normal and only the vertices remain superconducting, before finally the fully normal state is entered at an even higher temperature. Figure 3: Pair-density-wave states in a three-dimensional cube obtained from Ginzburg-Landau theory, for various temperatures T and fermionic population imbalances H. As temperature is increased, we observe a sequence of phase transitions. For low temperatures, superconductivity survives in the whole system. In the interior the solution is of Larkin-Ovchinnikov type but has higher gradients and different pattern on the surface due to reasons explained in the text. As temperature is increased, superconductivity vanishes first in the bulk, secondly on the surface and thirdly on the edges, such that eventually superconducting gap only exists in the vertices of the cube. The fermionic population imbalance H = 2πT /3 for each of the illustrations (a-e). The minimal and maximal valuesψ min andψ max are different for each of the different illustrations. In a −ψ min =ψ max = 0.7. In b −ψ min =ψ max = 0.4. In c and d −ψ min =ψ max = 0.2.

V. MICROSCOPIC BOGOLUIBOV-DE GENNES APPROACH
The existence of solutions with rapid oscillations in Ginzburg-Landau theory motivates investigation in microscopic models. We consider the mean-field Hamiltonian with spin-population imbalance with spin σ =↑, ↓, cite index i, fermionic annihilation and creation operators c iσ and c † iσ , chemical potential µ, spin-population imbalance h and nearest-neighbor hopping parameter t. We consider only on-site interaction potential V , leading to s-wave pairing, where the pairing field∆ i is defined as We diagonalize the Hamiltonian by performing the Bogoliubov canonical transformation where γ n and γ † n are annihilation and creation operators of quasi-particle excitations with energy E n . The pairing field satisfies the Bogoliubov-de-Gennes (BdG) self-consistency equation where denotes summation restricted over positive energies [22]. We will solve the problem numerically by taking a recursive approach. That is, we start with some initial guess for the pairing field, diagonalize the Hamiltonian in Eq. (13) and update the pairing field using Eq. (16), and repeat the process until convergence. In this study we will investigate the surface properties, which requires significantly larger system sizes. We consider both one-and two-dimensional systems, using open boundary conditions. In two dimensions we use a GPU-based approach, while in one dimension we use CPUs.
The obtained order parameter for a one-dimensional system, while varying the temperature, is shown in Figure 4. The remaining parameters are set to t = 1, V = 2, µ = 0.5 and h = 0.3. We can notice, as temperature is increased, a transition from the uniform superconducting state with upshoot on the boundaries, into the nonuniform LO state and eventually the surface-pair-density wave state, before finally transitioning into the normal T , for some fixed non-zero spin-population imbalance h = 0.3. As temperature is increased, we observe a sequence of phase transitions. Firstly, the system transitions from a uniform superconducting state with boundary enhancement, to the pair-density-wave FFLO state. Secondly, the pairing field vanishes in the bulk, while the pair-density-wave state remains non-zero close to the boundary, before finally transitioning fully into the normal conducting state. The behavior is of accordance with the prediction from Ginzburg-Landau theory in [20], from which the right panels have been taken. Remaining parameters are set to t = 1, V = 2 and µ = 0.5.
state. The solutions are very similar to those obtain from the Ginzburg-Landau approach in [20], shown on the same figure.
Note however that, at a microscopic level, the mechanism of formation of surface states is substantially more complex than the simple energy argument given in [20]. As pointed out in [20], at the level of Ginzburg-Landau theory, the state has oscillatory energy density and boundaries are accompanied with beneficial energy segments. However the ability of the system to start with such a segment depends on microscopic boundary conditions applied at a single-electron level. Moreover in order for surface states to appear, the Caroli-de Gennes-Matricon (CdGM) boundary conditions [23][24][25][26] should be violated. The fact that these boundary conditions in general do not hold was recently shown microscopically in [27]. It leads to the fact that for the simplest non-FFLO superconductors, there is an enhancement of the superconducting gap at the boundaries, and superconductivity of clean surfaces [27]. Therefore by varying the fermionic imbalance in the BdG formalism, one can investigate how the surface-pair-density-wave states are connected with the boundary states in non-FFLO regime discussed in [27]. By varying both the temperature and the spin-population imbalance, we obtained the complete phase diagram shown in Figure 5. The phase diagram shows there is a phase transition from the the PDW to non-PDW surface state as a function of temperature and fermionic imbalance.   (1) the surface-pair-density-wave state discussed in [20], and the surface superconducting state which can appear in conventional superconductors [27]. The crosses mark the points for which the pairing field is shown in Figure 4. Remaining parameters are set to t = 1, V = 2 and µ = 0.5.
Next we show the existence of the boundary states in two dimensions by solving the BdG equation in Eq. (16). The obtained pairing fields are shown in Figure 6, Figure  7 and Figure 8. In Figure 6, we show two states exhibit- ing bulk-modulated pairing fields for increasing temperature. For fixed spin-population imbalance, we see at low temperatures the appearance of an egg-crate pattern, shown in Figure 6a). As temperature is increased, the state transitions into the typical one-dimensional Larkin-Ovchinnikov solution, shown in Figure 6b). By fixing the spin-population imbalance and gradually increasing the temperature, we obtain the surface states shown in Figure 7. The state in Figure 7a) is identified as the edge-pair-density-wave, in which the pairing field is enhanced along the sample edges and decays as an exponentially dampened oscillation into the bulk. As temperature is increased, the pairing field vanishes on the edges but remain non-zero in the corners, as shown in Figure 7b). Finally, at even higher temperatures, the gap field vanishes completely and the system transitions to the normal state. These transitions would be the twodimensional analogue of the multiple phase transitions that occur in the three-dimensional cube in Figure 3. The BdG simulations confirm the existence of the boundary states.
Note however, that the precise structures of the edge and corner states obtained from the BdG calculations differ from those predicted by Ginzburg-Landau theory, which is not surprising given that Ginzburg-Landau model was derived from a different microscopic model than BdG model. The microscopic BdG results in Figure 7 exhibit parallel modulation and exponential decay, while the Ginzburg-Landau results in Figure 3 indicate modulation along the edges. The apparent conclusion obtained from both simulations is the existence of the edge and corner states and sequence of phase transitions. To demonstrate the presence of edge-modulated states on the BdG level we consider a sample with edges forming an angle of π/4 with the cristalline axis. The effect of this choice influences the interactions of the edge sites. In this case there appears modulation of the gap field along the boundaries. Figure 8 shows this result, characterized by the order parameter suppression in the bulk and enhancement in the corners.

VI. CONCLUSIONS
In conclusion, we considered surface effects in a superconductor where Cooper pairing involves electrons with finite center-of-mass momentum (Fulde-Ferrel-Larkin-Ovchinnikov instability). Firstly, we demonstrated that the standard microscopically-derived Ginzburg-Landau model cannot be used to describe such superconductors in the presence of boundaries in dimension larger than one since it gives a spurious divergence of the free energy near the boundaries. To describe such states we constructed a generalized functional where higher order derivative terms were retained to guarantee finiteness of the free energy density. We found that when the system forms FFLO states in the bulk, on its surface the modulation is different and has higher gradients and different patterns. The reservation that should be made in connection with the obtained patterns, is that calculations were made only for smooth and clean surfaces. We showed that at the mean-field-level the system undergoes a sequence of the phase transitions where each transition is associated with decreased dimensionality of superconducting state. Namely as temperature is increased the system first looses bulk superconductivity but retains two-dimensional superconductivity on its surfaces. The next phase transition is associated with the loss of two-dimensional superconductivity on the surfaces but retaining one-dimensional superconductivity on the edges. When the temperature is increased further the superconducting gaps survives only at the vertices. In order to demonstrate the existence of multiple phase transitions in a model that does not rely on Ginzburg-Landau expansion we solved numerically Bogoliubov-de-Gennes model in two dimensions. Clearly solution demonstrate the phase transition from the bulk to surface superconductivity.