Pair-density-wave superconductivity of faces, edges and vertices in systems with imbalanced fermions

TL;DR

This paper investigates boundary effects in FFLO superconductors with imbalanced fermions, revealing a sequence of phase transitions from bulk to edges and vertices, using both analytical and numerical methods.

Contribution

It extends the Ginzburg-Landau formalism for FFLO systems by including higher order terms and characterizes boundary-induced phase transitions in imbalanced fermionic superconductors.

Findings

Bulk FFLO state with boundary modulation patterns

Sequential phase transitions: bulk normal, edge superconducting, vertex superconducting

Numerical confirmation of phase transition sequence

Abstract

We describe boundary effects in superconducting systems with Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) superconducting instability, using Bogoliubov-de-Gennes and Ginzburg-Landau (GL) formalisms. First, we show that in dimensions larger than one the standard 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 cuboid sample of a 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 faces 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.

Figures

Equations

$$\mathcal{F}=$$

$$\displaystyle\mu|\psi|^{2}|\nabla\psi|^{2}+\frac{\mu}{8}\big{[}(\psi^{*}\nabla\psi)^{2}+\text{c.c.}\big{]}+\nu|\psi|^{6},$$

$$\alpha$$

$$\gamma$$

$$\nu$$

$$K_n(H,T)=\frac{2T}{(2\pi T{)}^n}\frac{(-1{)}^n}{(n-1)!}\text{Re}\left[{\Psi }^{(n-1)}(z)\right],$$

$$q^2=-\frac{2(\alpha L^2+12\beta )}{5\mu L^2},F=-\frac{(\alpha L^2+12\beta {)}^2}{60\mu L},$$

$$q_n^2$$

$$F_n$$

$$\psi =A{e}^{-kx+\mathrm{i}ky}$$

$${\mathcal{F}}_{\parallel }=\frac{6\alpha |A{|}^2+3\gamma |A{|}^4+2\nu |A{|}^6}{2k}+\frac{2\beta |A{|}^2+\mu |A{|}^4}{2}k.$$

$$\psi =A{e}^{-kx}\cos ky$$

$$\kappa |\nabla \psi {|}^4.$$

$$\kappa \simeq \left(-\frac{29}{32}\right)\frac{\pi N(0)v_{\text{F}}^4{\Omega }_d}{2}{K}_7(H,T),$$

$$\kappa |A{|}^4{k}^3,$$

$$\tilde{\kappa }=-\frac{1}{v_{\text{F}}^4}\frac{{\hat{\beta }}^2}{2{\hat{\delta }}^2}\frac{\kappa \gamma }{{\nu }^2}.$$

$${\int }_{\Omega }|\Delta \psi {|}^2\text{d}V={\int }_{\Omega }\sum _{i,j=1}^3|{\partial }_i{\partial }_j\psi {|}^2\text{d}V+{\oint }_{\partial \Omega }f\text{d}S,$$

$$\displaystyle f=\sum_{i,j=1}^{3}\Big{(}\big{(}\partial_{i}^{2}\psi^{*}\big{)}\partial_{j}\psi-\big{(}\partial_{i}\partial_{j}\psi^{*}\big{)}\partial_{i}\psi\Big{)}n_{j},$$

$$\psi (x,y)=\eta (x)+\eta (y),\eta (x)=e^{-kx}\cos (qx+\varphi ).$$

$$S={\int }_0^{\frac{1}{T}}d\tau \sum _{i,j=1,\sigma =\downarrow }^{i,j=N,\sigma =\uparrow }{\psi }_{i\sigma }^{†}({\partial }_{\tau }+h_{ij\sigma }){\psi }_{j\sigma }$$

$$-V\sum _{i=1}^N{\psi }_{i\uparrow }^{†}{\psi }_{i\downarrow }^{†}{\psi }_{i\downarrow }{\psi }_{i\uparrow }$$

$$Z=\int D[{\psi }^{†},\psi ]e^{-S}.$$

$$\exp \left[V\int d\tau \sum _i{\psi }_{i\uparrow }^{†}{\psi }_{i\downarrow }^{†}{\psi }_{i\downarrow }{\psi }_{i\uparrow }\right]=\int D[{\Delta }^{†},\Delta ]$$

$$\exp \left[-\int d\tau \sum _i\left(\frac{{\Delta }_i^{†}{\Delta }_i}{V}+{\Delta }_i^{†}{\psi }_{i\downarrow }{\psi }_{i\uparrow }+{\Delta }_i{\psi }_{i\uparrow }^{†}{\psi }_{i\downarrow }^{†}\right)\right].$$

$${\Psi }_i^{†}=\left({\psi }_{i\uparrow }^{†},{\psi }_{i\downarrow }\right),{\Psi }_i=\left(\begin{array}{c}{\psi }_{i\uparrow }\\ {\psi }_{i\downarrow }^{†}\end{array}\right)$$

$$S=\int d\tau \sum _i\frac{{\Delta }_i^{†}{\Delta }_i}{V}+\sum _{ij}{\Psi }_i^{†}\left({\delta }_{ij}{\partial }_{\tau }+H_{ij}\right){\Psi }_j$$

$$H_{ij}=\left(\begin{array}{cc}{h}_{ij\uparrow }& {\delta }_{ij}{\Delta }_i\\ {\delta }_{ij}{\Delta }_i^{†}& -h_{ji\downarrow }\end{array}\right).$$

$$Z=\int D[{\Delta }^{†},\Delta ]\exp \left[\ln \det \left({\partial }_{\tau }+H\right)-\int d\tau \sum _i\frac{{\Delta }_i^{†}{\Delta }_i}{V}\right].$$

$$\ln \det \left({\partial }_{\tau }+H\right)=\text{Tr}\sum _{n=-\infty }^{\infty }\ln \left(\mathrm{i}{\omega }_n+H\right)=-\text{Tr}\ln f(H)$$

$$S=\sum _i\frac{|{\Delta }_i{|}^2}{TV}+\text{Tr}\ln f(H).$$

$${\Delta }_i=-V{\mathbf{e}}_i^{T}f(H){\mathbf{h}}_i=-V⟨c_{i\uparrow }{c}_{i\downarrow }⟩$$

$${\Delta }_i=-V\sum _{n=1}^{2N}{u}_i^n{v}_i^{n\ast }f(E_n)$$

$$\Omega =\text{mbox}tr\left[\mathcal{H}{\rho }_{MF}+T{\rho }_{MF}\ln {\rho }_{MF}\right]$$

$$U_{ij}=\left(\begin{array}{cc}{\delta }_{ij}{U}_{i\uparrow }& 0\\ 0& -{\delta }_{ij}{U}_{i\downarrow }\end{array}\right)$$

Full text

I Introduction

Fulde and Ferrell Fulde and Ferrell (1964) and Larkin and Ovchinnikov Larkin and Ovchinnikov (1964) (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 Zwierlein (2006); Radzihovsky (2011); Kinnunen et al. (2018); Samoilenka et al. (2020). 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 Alford et al. (2001). Even in electronic superconductors, finite momentum pairing may arise due to reasons other than application of an external magnetic field Barkman et al. (2019a). This state has for a long time been of great interest and was searched for in a number of superconducting materials Bianchi et al. (2003); Cho et al. (2017); Uji et al. (2018); Mayaffre et al. (2014); Coniglio et al. (2011); Norman (1993); Cho et al. (2011); Matsuda and Shimahara (2007); Lortz et al. (2007); Singleton et al. (2000); Martin et al. (2005).

In the recent work Barkman et al. (2019b), using microscopically derived Ginzburg-Landau (GL) 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 boundary but the bulk of the systems is a normal metal.

In this work we answer two questions. (i) We demonstrate this effect in microscopic models for imbalanced fermions without relying on a Ginzburg-Landau expansion. We show that in dimensions larger than one, the superconductor undergoes a sequence of phase transition where superconductivity survives in sub-domains of sequentially lower dimension. That is, by increasing the temperature and fermionic population imbalance in a square system, at the mean-field level, the system will undergo the following sequence of the phase transitions: $\text{superconducting bulk}\to\text{superconducting edges}\to\text{superconducting corners}\to\text{normal state}$. (ii) For a Ginzburg-Landau approach, we explain that, as alluded in Barkman et al. (2019b), in dimensions larger than one, it is necessary to modify the Ginzburg-Landau expansion compared to what is done in the standard calculations Buzdin and Kachkachi (1997) to describe systems with boundaries. By retaining higher order derivative terms, we demonstrate the same sequence of phase transition obtained in the microscopic model. That is, a three dimensional cuboid superconductor with imbalance fermions undergoes the sequence of phase transitions: $\text{superconducting bulk}\to\text{superconducting faces}\to\text{superconducting edges}\to\text{superconducting vertices}\to\text{normal state}$, as temperature and fermionic population imbalance is increased. We also demonstrate an alternative, simpler Ginzburg-Landau expansion in the presence of boundaries, which also exhibits the boundary pair-density-wave state, but does not capture the same sequence of phase transitions.

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 higher order derivative terms in the expansion. That allows us to in Section IV, obtain the face, edge and vertex states without divergent energies. In Section V we construct an alternative form for the Ginzburg-Landau expansion which does not suffer from the divergence near the boundary, but does not capture the full set of different types of boundary states.

In Section VI we solve numerically Bogoluibov-de-Gennes (BdG) equation to show that the states exist in microscopic models that does not rely on any Ginzburg-Landau expansion. In the Appendix we confirm the existence of the boundary pair-density-wave by solving a more general BdG including Hartree terms.

II Breakdown of the standard Ginzburg-Landau approach in finite systems

II.1 Ginzburg-Landau Model

The Ginzburg-Landau description of superconductors in the presence of fermionic population imbalance was derived from microscopic theory for superconductors in Buzdin and Kachkachi (1997). The free energy functional reads $F[\psi]=\int_{\Omega}\mathcal{F}\text{d}^{d}x$ where the free energy density $\mathcal{F}$ is

[TABLE]

where $\psi$ is a complex field referred to as the superconducting order parameter and c.c. denotes complex conjugation. The coefficients $\alpha,\gamma$, and $\nu$ depend on the fermionic population imbalance $H$ and temperature $T$ accordingly Buzdin and Kachkachi (1997)

[TABLE]

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

[TABLE]

where $z=\frac{1}{2}-\mathrm{i}\frac{H}{2\pi T}$ and $\Psi^{(n)}$ is the polygamma function of order $n$. The remaining coefficients are given as $\beta=\hat{\beta}v_{\text{F}}^{2}\gamma$, $\delta=\hat{\delta}v_{\text{F}}^{4}\nu$, and $\mu=\hat{\mu}v_{\text{F}}^{2}\nu$, where $v_{\text{F}}$ is the Fermi velocity and $\hat{\beta},\hat{\delta},\hat{\mu}$ are positive constants that depend on the dimensionality $d$. In one dimension we have $\hat{\beta}=1$, $\hat{\delta}=1/2$, and $\hat{\mu}=4$, in two dimensions we have $\hat{\beta}=1/2$, $\hat{\delta}=3/16$, and $\hat{\mu}=2$, and in three dimensions we have $\hat{\beta}=1/3$, $\hat{\delta}=1/10$ and $\hat{\mu}=4/3$. The Ginzburg-Landau description is valid in the parameter regime in which the highest order terms are positive, that is where $\nu$ is positive. In the parameter regime in which $\beta$ is negative, inhomogeneous order parameters may be energetically favorable. The first considered structures of the order parameter are the so-called Fulde-Ferrell (FF) state $\psi_{\text{FF}}\propto\mathrm{e}^{\mathrm{i}px}$ and the Larkin-Ovchinnikov (LO) state $\psi_{\text{LO}}\propto\cos px$, with $p^{2}=-\frac{\beta}{2\delta}$ and transition into the normal state at $\alpha=\alpha_{c}^{\text{bulk}}=\frac{\beta^{2}}{4\delta}$.

Inhomogeneous states can appear when the term $|\nabla\psi|^{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 $|\nabla\psi|^{2}$ favors creation of the gradients, while the positive term $|\Delta\psi|^{2}$ is added to bound gradients from above. However we will show that in some cases, the inclusion of the stabilizing term $|\Delta\psi|^{2}$ is not sufficient. In small systems and generically in finite systems in two dimensions or higher, there exist solutions that satisfy $\Delta\psi=0$. These states are often characterized by enhancement of $\psi$ close to the boundary. The associated energy of such state can potentially diverge and one needs to resort to a more general Ginzburg-Landau theory.

II.2 Small systems

Consider first a one dimensional domain $\Omega=[-L/2,L/2]$ of length $L$. The equation $\Delta\psi=0$ is satisfied by the real field $\psi(x)=qx$, for any parameter $q$. The modulus $|\psi|$ 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 $\psi_{\text{LO}}\propto\cos px$ or $\psi_{\text{FF}}\propto e^{\mathrm{i}px}$ if $L\gtrsim 1/p$, where $p^{2}=-\frac{\beta}{2\delta}$. However, for small system sizes, the potential term does become negligible in comparison to the beneficial gradient term $\beta|\nabla\psi|^{2}=\beta q^{2}$. To lowest order in $L$, minimizing with respect to $q^{2}$ gives

[TABLE]

with transition into the normal state at $\alpha=\alpha_{c1}=-12\beta/L^{2}$. In the limit of infinitesimal system size $L\to 0$, we find that $\alpha_{c1}\to\infty$, and the associated momentum $q$ and energy $F$ diverges. Note that the value of $\psi$ at the boundary, $\psi(\pm L/2)=\pm qL/2$, does not vanish in this limit since $q\propto 1/L$. Consequently the value of $\psi$ 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 $\psi_{n}(\bm{x})=q_{n}\prod_{j=1}^{n}x_{j}$, where $n\leq d$. Analogously to the one-dimensional case we find

[TABLE]

with transition into the normal state at $\alpha_{cn}=n\alpha_{c1}$. Consequently we expect that in higher dimensions, we will see a sequence of transitions from $(n=1)\to(n=2)\to\ldots\to(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 $\psi_{n}$ is constant and independent of system size at the vertices of the $n$-dimensional cube, for significantly small systems.

II.3 Boundary state instability in Ginzburg-Landau model

Consider a two-dimensional domain $\Omega$ in the presence of boundaries. We can represent the domain $\Omega$ as section of the complex plane, by defining a complex coordinate $\zeta=x+\mathrm{i}y$. The order parameter $\psi$ is now a function of $\zeta$. Clearly if $\psi$ is an analytic complex valued function (i.e. it satisfies the Cauchy-Riemann equations), it follows automatically that $\Delta\psi=0$. The system becomes unstable towards states of the form $\psi(\zeta)=Ae^{k\zeta}$ with $|k|\to\infty$, 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 $\psi(\zeta)=Ae^{k\zeta}$ would diverge in some direction, which would be penalized by the potential terms. However, the divergence of $|\psi|$ 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\geq 0$. An analytic function which satisfies $\Delta\psi=0$ is

[TABLE]

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

[TABLE]

We see that as long as the amplitude $|A|^{2}<-2\beta/\mu$, the energy diverges to $-\infty$ as $k$ increases. This shows that the model has an instability towards formation of boundary states with rapid modulation along the boundary. Analogous calculation for the state

[TABLE]

would yield the same conclusion, see Figure 2.

III Generalisation of the Ginzburg-Landau model

In contrast to the example of linear solutions in small systems, see Section II.2, 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 boundary states. Note, however that GL solutions can only represent the long-wave-length physics of how order parameter decays away from the $1/k_{F}$ length scales from the boundary. Hence, 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 $|\nabla\psi|^{2}$, which is of second order in both momentum and amplitude. Conventionally it is stabilized by the term $|\Delta\psi|^{2}$, which is of fourth order in momentum and second order in amplitude. Since there exists states that satisfy $\Delta\psi=0$, we have to resort to higher order in either momentum or amplitude. Terms such as $|\nabla\Delta\psi|^{2}$, which is of sixth order in momentum and second order in amplitude, would not be helpful since it is automatically zero if $\Delta\psi=0$. Therefore we resort to retaining the term which is of fourth order in momentum and fourth order in amplitude, which reads

[TABLE]

Since this term is of higher order in amplitude than the lowest order beneficial gradient term, there will exist inhomogeneous superconducting states for any $\alpha$, as long as the gradient coefficients are negative. This term was used to find two-dimensional solutions in Barkman et al. (2019b). Here we derive from the microscopic theory the following estimate for the coefficient $\kappa$ in Eq. (11)

[TABLE]

where factor $\Omega_{d}$ depends on the dimension $d$, where $\Omega_{1}=1$, $\Omega_{2}=3/8$ and $\Omega_{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. (11) 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 $|\psi|^{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 boundary state $\psi=Ae^{-kx+\mathrm{i}ky}$, which previously was associated with negative infinite energy and infinite momentum $k$. With the inclusion of the term in Eq. (11), the associated free energy in Eq. (9) obtains the additional term

[TABLE]

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 CUDA-enabled NVIDIA graphical processing unit (GPU). It is convenient to introduce rescaled coordinates and parameters defined accordingly: $\tilde{\psi}=\psi/|\psi_{\text{U}}|$, $\tilde{\alpha}=\alpha/\alpha_{\text{U}}$, $\tilde{x}=px$, where $|\psi_{\text{U}}|^{2}=-\frac{\gamma}{2\nu}$, $\alpha_{\text{U}}=\frac{\gamma^{2}}{4\nu}$ and $p^{2}=-\frac{\beta}{2\delta}$. Expressed in these quantities, the free energy reads $F[\psi]=\frac{\alpha_{\text{U}}|\psi_{\text{U}}|^{2}}{p^{d}}\tilde{F}[\tilde{\psi}]$, where $\tilde{F}[\tilde{\psi}]=\int_{\Omega}\tilde{\mathcal{F}}\text{d}^{d}\tilde{x}$, where the rescaled free energy density $\tilde{\mathcal{F}}$ is identical to Eq. (1), but the coefficients have been replaced accordingly: $\alpha\mapsto\tilde{\alpha}$, $\beta\mapsto\tilde{\beta}$, and so on, where $\tilde{\gamma}=-2\tilde{\nu}=-2$, $\tilde{\beta}=-2\tilde{\delta}=-2\hat{\beta}^{2}/\hat{\delta}$, and $\tilde{\mu}=\hat{\beta}\hat{\mu}/\hat{\delta}$. Among these coefficients, all are constant except $\tilde{\alpha}$, which parametrizes both temperature $T$ and fermionic population imbalance $H$. With the inclusion of the additional term in Eq. (11), in rescaled coordinates, the coefficient $\tilde{\kappa}$ reads

[TABLE]

Specifically in three dimensions we have that $\tilde{\kappa}=\frac{145K_{3}K_{7}}{18K_{5}^{2}}$, which is not constant in the rescaled units. However, studying the functions $K_{n}$ in Eq. (5), we see that $K_{n}\propto\frac{T}{T^{n}}\cdot f_{n}(H/T)$, where $f_{n}$ is some elementary function. Therefore, if we study a line in the $TH$-plane where $H/T$ is constant, the rescaled parameter $\tilde{\kappa}$ will also be constant. Particularly, we will study the line where $H/T=2\pi/3$, along which $\tilde{\kappa}\simeq 0.5752$. On a technical note, since we are working in rescaled units, where we measure length in units of $p^{-1}\propto\sqrt{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 boundary, 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 faces 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 faces of the cube, while the bulk now has transition into a normal state. Second, when the temperature is elevated further, the faces 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.

V Resolving the boundary instability in an alternative Ginzburg-Landau expansion

Let us consider how to obtain a model which does not suffer from divergences near a boundary without introduction of the higher order terms, as in Section III.

Firstly, observe that when one uses the usual microscopic derivation of the model Eq. (1) one does not obtain fully unambiguously the form of the gradient terms. The ambiguity is due to freedom in integration by parts. For example, the term $-\psi\Delta\psi^{*}+c.c.$ is not discriminated from $|\nabla\psi|^{2}$. Other way to put it is that the standard microscopic derivation fixes Euler-Lagrange equations but does not fully fix the boundary conditions.

The GL free energy in Eq. (1) is unbounded from below due to the existence of boundary states satisfying $\Delta\psi=0$. We can rewrite the Laplacian term $|\Delta\psi|^{2}$ using integration by parts as

[TABLE]

where $\bm{n}$ is the normal to the boundary $\partial\Omega$. Since the standard microscopic derivation does not unambiguously determine the boundary contribution, we consider now the free energy expansion where $|\Delta\psi|^{2}$ is replaced by $\sum_{ij}|\partial_{i}\partial_{j}\psi|^{2}$ and the boundary term $f$ is discarded. The boundary state in Eq. (8) does not nullify $\sum_{ij}|\partial_{i}\partial_{j}\psi|^{2}$, and only linear function $\psi=a+\sum_{i=1}^{3}b_{i}x_{i}$ with constants $a,\ b_{i}$ does. Note that in this alternative form, the free energy does not suffer from boundary divergences and the inclusion of higher order terms is unnecessary.

We simulated numerically the GL model given by Eq. (1) with Laplacian term replaced by $\sum_{ij}|\partial_{i}\partial_{j}\psi|^{2}$. The typical configuration of the order parameter is presented in Figure 4. The main difference from the GL model with usual Laplacian term and higher order terms, see Figure 3, is that we do not recover the full four-step sequence (bulk-faces-edges-vertices) of the phase transitions in the cuboid geometry. Instead in this approximation it retains a two-step transition only from bulk to boundary superconductivity and then to normal state everywhere.

Now we can construct an analytical guess for the edge state, which is the same as in Barkman et al. (2019b). Namely, let us consider superconductor positioned at $x,\ y>0$ and for all $z$. Then configuration of $\psi$ in Figure 4b can be approximated by

[TABLE]

Even though this configuration has enhanced $\psi$ on the edges, due to the overlap of $\eta(x)$ and $\eta(y)$, it will, in this approximation, have the same critical temperature, as the critical temperature of the face superconductivity in the limit of infinitely large sample.

VI Microscopic Bogoluibov-de Gennes (BdG) Approach

VI.1 Derivation of BdG from path integral

We start from Fermi-Hubbard model in path integral formulation defined by action (first part of the derivation follows Altland and Simons (2010))

[TABLE]

with Grassman fields $\psi_{i\sigma}(\tau),\psi^{\dagger}_{i\sigma}(\tau)$ corresponding to fermionic creation and annihilation operators defined on some lattice indexed by $i,j\in[1,N]$. Kinetic part of action is given by $h_{ij\sigma}=-(\mu+\sigma h)\delta_{i,j}-t\delta_{|i-j|,1}$, with chemical potential $\mu$, spin-population imbalance $h$ and nearest-neighbor hopping parameter $t$. Interaction is set by on-site attractive potential $V>0$, which leads to $s$-wave pairing. Partition function is given by path integral

[TABLE]

We perform Hubbard-Stratonovich transformation, which introduces auxiliary bosonic field $\Delta_{i}(\tau)$. Up to a constant the interaction term becomes:

[TABLE]

By introducing Nambu spinors

[TABLE]

the action becomes

[TABLE]

With $2N\times 2N$ matrix $H$ defined by:

[TABLE]

Integrating out fermionic degrees of freedom we obtain

[TABLE]

Note, that up till this point this model exactly corresponds to Fermi-Hubbard model Eq. (17). To simplify it and proceed with mean field approximation we make two assumptions: 1) $\Delta_{i}$ is classical, i.e. it does not depend on $\tau$ 2) $\Delta_{i}$ does not fluctuate thermally, i.e. it extremizes the action $S$. Note that this mean-field nature of BdG approach will give us a sequence of second order transitions for the superconductivity of the bulk, face, edge, and vertex states. Beyond mean-field approximation the phase fluctuations make the face phase transitions of Berezinskii-Kosterlitz-Thouless type.

Next, one can simplify the action by reducing the term

[TABLE]

where last equality holds up to a constant, summation is performed over Matsubara frequencies $\omega_{n}=2\pi T\left(n+1/2\right)$ and Fermi factor $f(x)=1/(e^{x/T}+1)$. Note, that trace sums over $2N$ diagonal elements of function of $2N\times 2N$ dimensional matrix $H$ and hence equals to just sum over a function evaluated at the eigenvalues of $H$. Finally, the action is given by

[TABLE]

Taking variation of Eq. (25) with respect to $\Delta^{*}_{i}$ we obtain the self consistency equation111Note, that this formula is equivalent to the usual one with hyperbolic tangent since $-\bm{e_{i}}^{T}f(H)\bm{h_{i}}=\frac{1}{2}\bm{e_{i}}^{T}\tanh\frac{H}{2T}\bm{h_{i}}$.

[TABLE]

where vectors, following Nagai et al. (2012), are denoted as $(\bm{e}_{i})_{j}=\delta_{i,j}$ and $(\bm{h}_{i})_{j}=\delta_{i+N,j}$. Note, that in this approach, based on path integral formulation, we obtain equation only for $\Delta_{i}$ Eq. (26). If one instead employs variational method self consistency equations are also obtained for Hartree terms, which shift the chemical potentials. In the appendix we show that the boundary states, discussed here, exist in the model with Hartree term as well.

We solve the problem numerically by taking a recursive approach. That is, we start with some initial guess for the pairing field $\Delta_{i}$. Secondly, we diagonalize the Hamiltonian in Eq. (22), i.e. we calculate its $2N$ eigenvalues $E_{n}$ and eigenvectors $w_{i}^{n}=\begin{pmatrix}u_{i}^{n}\\ v_{i}^{n}\end{pmatrix}$. Then $\Delta_{i}$ is updated through the self consistency equation Eq. (26), which explicitly can be written as

[TABLE]

The process is repeated until convergence.

VI.2 Numerical results

In this study we will investigate the surface properties, which requires significantly large 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 5. The remaining parameters are set to $t=1$, $V=2$, $\mu=0.5$ and $h=0.3$. Note that this parameter set does not correspond to the microscopic physics where Ginzburg-Landau expansion was obtained, so we should not expect quantitative similarities. Rather the main question we ask if the boundary states exist beyond the Ginzburg-Landau regime. To emphasize that the results do not require finetuning in the Appendix we give an example of a boundary state in a different Bogoliubov-de-Gennes model with a Hartree term retained. We can notice, as temperature is increased, a transition from the uniform superconducting state with upshoot on the boundaries, into the non-uniform LO state and eventually the boundary state, referred to as the surface pair-density wave state, before finally transitioning into the normal state. The solutions are very similar to those obtain from the Ginzburg-Landau approach in Barkman et al. (2019b), shown on the same figure. Note however that the microscopical considerations using BdG are not performed in the weak-coupling limit where Ginzburg-Landau theory applies. Therefore, a quantitative agreement between these two models are neither necessary nor trivially expected. We emphasize that the similarity between the two approaches concern the order parameter structure and not quantitative values in the phase diagram.

Note, that in the tight-binding model, non-zero values of the pairing field $\Delta_{i}$ on the boundary, as in Figure 5, does not mean that it is discontinuous. Instead $\Delta_{i}$ decays to zero outside the superconductor over length scale corresponding to the localization of the Wannier function’s tail or lattice spacing.

At a microscopic level, the mechanism of formation of boundary states has many complex aspect beyond the simple energy argument given in Barkman et al. (2019b). As pointed out in Barkman et al. (2019b), 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. The energetic argument clearly indicates that Caroli-de Gennes-Matricon (CdGM) boundary conditions that conjecture zero normal derivative of the gap modulus de Gennes (1964); Caroli et al. (1962, 1963); De Gennes and Guyon (1963) should not be expected to hold for pair-density-wave case because the system has clear energetic preference for having high gradients of the gap modulus near the boundary. In fact that the Caroli-de Gennes-Matricon (CdGM) boundary conditions do not hold even for ordinary superconductors, at least for clean surfaces, as was recently shown microscopically in Samoilenka and Babaev (2020). 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 Samoilenka and Babaev (2020). 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 Samoilenka and Babaev (2020). By varying both the temperature and the spin-population imbalance, we obtained the complete phase diagram shown in Figure 6. The phase diagram shows both the PDW boundary state regime and the non-PDW boundary state regime.

Next we show the existence of the boundary states in two dimensions by solving the BdG equation in Eq. (27). Obtained pairing fields are shown in Figure 7, Figure 8 and Figure 9. Consideration of two-dimensional samples are computationally much more demanding than one-dimensional. Hence we do not perform full phase diagram exploration, but demonstrate boundary states by fixing the spin-population imbalance $h=0.4$ and gradually increase the temperature, see Figure 7. First, we study a square system which is periodic in all directions – in this case no boundary states can form and we obtain bulk critical temperature $T_{c1}$ by identifying the temperature at which the maximal pairing field goes to zero. Next, for a sample with periodic boundary conditions only in one direction, edge states are formed, which gives edge critical temperature $T_{c2}$, see Figure 7. Finally, for open boundaries corner states survive even for higher temperatures, until the corner critical temperature $T_{c3}$ is reached.

For selected values of $T$ we plot configurations of order parameter in Figure 8. The state in Figure 8a is the usual LO modulation in one direction. States in Figure 8b, c are 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 8d. Finally, at even higher temperatures, the gap field vanishes completely and the system transitions to the normal state. These transitions would be the two-dimensional 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 GL theory, which is not surprising given that they correspond to different physical systems. Indeed, BdG is considered in tight binding limit and hence corresponds to localized electrons with much stronger interactions than free electron approximation conventionally used to obtain GL model. The microscopic BdG results in Figure 8 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.

On the other hand, when replacing the Laplacian term by $\sum_{ij}|\partial_{i}\partial_{j}\psi|^{2}$ we observe a strong resemblance between the GL and BdG simulations - compare for example we obtain configuration of the order parameter in GL, which is very close to that of BdG: compare Figure 4b and Figure 8c. In this case modulation of $\psi$ is mostly orthogonal to the boundary. However, the modification of GL model with the Laplacian term replaced by $\sum_{ij}|\partial_{i}\partial_{j}\psi|^{2}$, in contrast to the BdG solutions, does not have the phases with gap formation on the edges and vertices in the absence of the gaps on the faces.

To demonstrate the presence of edge-modulated states on the BdG level we consider a sample with edges forming an angle of $\pi/4$ with the crystalline 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 9 shows this result, characterized by the order parameter suppression in the bulk and enhancement in the corners.

In context of cold atoms a phase separation between FFLO-like bulk and normal or BCS-like boundary was predicted in a one-dimensional harmonic trap Orso (2007); Hu et al. (2007); Liu et al. (2008). Note, that these states appear as different phase due to spatially variating local density approximation chemical potential $\mu(z)=\mu_{0}-V(z)$ inside the trap making the system form different states at different distances from the trap center. Where the potential is deep, i.e. in the center of a trap, FFLO state appears while near the boundary of the trap, normal or BCS states are formed. By contrast, the states that we find not only are represented by a very different solution but appear under different conditions: induced by hard wall potential, when particles are not permitted to hop outside the sample, which is equivalent to infinite chemical potential at a boundary. It does not have a direct counterpart in the phase separation picture in traps. Another principal difference is that our states rely on the density of state structure near boundaries, as can be most clearly seen in a BCS superconductor with zero population imbalance Samoilenka and Babaev (2020).

VII Conclusions

In conclusion, we considered boundary 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 with higher order derivative terms were retained to guarantee finiteness of the free energy density. 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 the superconducting state. Namely as temperature is increased the system first looses bulk superconductivity but retains two-dimensional superconductivity on its faces. The next phase transition is associated with the loss of two-dimensional superconductivity on the faces but retaining one-dimensional superconductivity on the edges. When the temperature is increased further the superconducting gaps survives only at the vertices. The transition from the bulk to surface pair-density-wave superconductivity is also captured by an alternative and simpler form of the Ginzburg-Landau functional described in Section V.

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 boundary superconductivity.

The results provide a new route to identify and study properties of pair-density-wave superconductors. The results predict a discrepancy between the bulk-dominated specific heat identification of critical temperature and surface-dominated current and diamagnetic response probes. The predicted possibility of a different modulation of the order parameter on the surface also has direct implication for scanning tunneling microscopy and Josephson experiments.

The current work also opens up the question of the physics of boundary states beyond mean-field approximation. Beyond mean-field approximation, quasi-one-dimensional superconductors are always phase disordered in thermodynamic limit at any finite temperature. However, again as in any superconducting wire that statement applies only to thermodynamic limit and it’s expected to see the mean-field sequence of the phase transitions in local probes such as scanning tunneling microscopy. The surface phase transition beyond mean-field approximation should be of the Berezinskii-Kosterlitz-Thouless type. However multiple transitions are also possible beyond mean-field in cases where the modulation is along the surface, associated with melting of stripes and loss of crystalline modulation order. This will be a counterpart of multiple transitions beyond mean-field approximation that attract currently substantial interest in the conventional two-dimensional Larkin-Ovchinnikov cases Agterberg and Tsunetsugu (2008); Berg et al. (2009) that yields additional paired and charge-$4e$ phases. Our results indicate that beyond mean-field approximation, the charge charge-$4e$ superconductivity via Berezinskii-Kosterlitz-Thouless type transition can occur also at the surfaces of bulk superconductors (for a detailed general discussion of the fluctuation-induced paired phases in low-dimensional multicomponent superconductors see Babaev ; Svistunov et al. (2015)).

Acknowledgements.

We thank Robert Vedin for his valuable assistance. The work was supported by the Swedish Research Council Grants No. 642-2013-7837, 2016-06122, 2018-03659 and Göran Gustafsson Foundation for Research in Natural Sciences and Medicine. A part of this work was performed at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611 We gratefully acknowledge the support of NVIDIA Corporation with the donation of the Quadro P6000 GPU used for this research. Published with the support from the Längman Culture Foundation.

Appendix

It is common to use variational approach to derive BdG equations De Gennes (2018). Namely, the self consistency equations are obtained by setting the variation of the free energy $\Omega$ with respect to mean fields to zero.

[TABLE]

where $\mathcal{H}$ is Hamiltonian corresponding to exact Fermi-Hubbard model Eq. (17) and $\rho_{MF}=\frac{e^{-H_{MF}/T}}{\mbox{tr}e^{-H_{MF}/T}}$ with mean field Hamiltonian $H_{MF}$.

The main difference from our derivation in Section VI is that one has to include additional mean fields in $H_{MF}$ via Hartree term $\sum_{i,\sigma}U_{i\sigma}\psi^{\dagger}_{i\sigma}\psi_{i\sigma}$. It amounts to replacing matrix $H$ in self consistency Eq. (26) for $\Delta_{i}$ with $M=H+U$, where

[TABLE]

In addition, the fields $U_{i\sigma}$ have to be solved for self consistently using

[TABLE]

In Figure 10 we present results in one dimension which incorporate self-consistent solution for both the pairing field and the Hartree fields, for various temperatures, keeping the remaining parameters fixed. We show result using both periodic and open boundary conditions, in order to investigate the existence of boundary states. We find that boundary states also exist in this model but with a smaller difference between the critical temperatures. Note, that modulation in spin population imbalance $\langle n_{\uparrow}-n_{\downarrow}\rangle$ is much smaller and is twice as frequent as that of $\Delta$. The later is easy to explain by the fact that each of the populations $\langle n_{\sigma}\rangle$ in Eq. (30) is even function of $\Delta$. Hence they depend on $|\Delta|$, which has twice shorter period than $\Delta$. Note that the spin-population imbalance modulation is enhanced close to the boundary, but the decay length significantly shorter than for the pairing potential. This enhancement remains even for higher temperatures, when superconductivity is completely lost and is related to Friedel oscillations.