Superconductivity model for a spin-vortex checkerboard

TL;DR

This paper proposes a microscopic model for fermionic excitations in a spin-vortex checkerboard, aiming to explain spin and charge modulations in cuprates, and finds a temperature-dependent energy gap consistent with experiments.

Contribution

It introduces a generalized fermionic model for spin-vortex textures and analyzes its mean-field solution, connecting theoretical predictions with experimental observations in cuprate superconductors.

Findings

Model captures fermionic correlations in spin-vortex structures.

Energy gap temperature dependence aligns with experimental data.

Superconductivity emerges only with added small hopping terms.

Abstract

We introduce a microscopic model aimed at describing the behavior of fermionic excitations in the background of a magnetic texture called "spin-vortex checkerboard". This texture was proposed previously as a possible alternative to stripes to interpret the experimental phenomenology of spin and charge modulations in 1/8-doped lanthanum cuprates. The model involves two kinds of interacting fermionic excitations residing in spin-rich and spin-poor regions of the modulated structure. It is a generalization of another model developed earlier for the so-called "grid checkerboard". The principal terms of our model describe the decay of fermionic pairs belonging to spin-poor regions into single fermions occupying spin-rich regions and vice versa. These terms induce intricate fermionic correlations throughout the system but fall short of inducing superconductivity unless arbitrarily small hopping terms are added to the model Hamiltonian. We present the mean-field solution of the model, including, in particular, the temperature dependence of the energy gap. The latter is found to be in a good overall agreement with available experimental data for high-$T_c$ cuprate superconductors.

Figures

Equations

$$\begin{gathered} H=\sum _{i,\eta }{ϵ}_a{a}_{i\eta }^{†}{a}_{i\eta }+\sum _{i,\zeta }{ϵ}_c{c}_{i\zeta }^{†}{c}_{i\zeta }+\sum _{i,\alpha ,\eta }{ϵ}_b{b}_{i\alpha \eta }^{†}{b}_{i\alpha \eta } \\ +g\sum _{i,\alpha }[\left(b_{i\alpha \uparrow }^{†}{b}_{i\alpha \downarrow }^{†}{a}_{j{\text{e}}_{[i,\alpha ]}}{a}_{k{\text{o}}_{[i,\alpha ]}}+h.c\right) \\ +\left(b_{i\alpha \uparrow }^{†}{b}_{i\alpha \downarrow }^{†}{c}_{m{\text{e}}_{[i,\alpha ]}}{c}_{n{\text{o}}_{[i,\alpha ]}}+h.c\right)], \end{gathered}$$

$$\begin{array}{rl}{b}_{i\alpha \uparrow }& =s{B}_{i\alpha \uparrow }+w{e}^{i{\phi }_{\alpha }}{B}_{i\alpha \downarrow }^{†},\\ b_{i\alpha \downarrow }& =s{B}_{i\alpha \downarrow }-w{e}^{i{\phi }_{\alpha }}{B}_{i\alpha \uparrow }^{†},\end{array}$$

$$s^2+w^2=1,$$

$$H_{a}=8\epsilon_{b}N\left[s^{2}n_{B}+w^{2}(1-n_{B})\right]\\ +\epsilon_{a}\sum_{i,\eta}a^{\dagger}_{i\eta}a_{i\eta}+\epsilon_{c}\sum_{i,\zeta}c^{\dagger}_{i\zeta}c_{i\zeta}\\ +gsw(1-2n_{B}){{\sum}}_{\alpha}\bigg{[}(e^{-i\varphi_{\alpha}}\sum_{ij}a_{j\text{e}_{[i,\alpha]}}a_{k\text{o}_{[i,\alpha]}}+h.c)\\ +(e^{-i\varphi_{\alpha}}\sum_{ij}c_{m\text{e}_{[i,\alpha]}}c_{n\text{o}_{[i,\alpha]}}+h.c)\bigg{]},$$

$$n_B=\frac{1}{\exp \frac{{ϵ}_B}{T}+1},$$

$$\mathbf{L}=(0,4),$$

$$\begin{array}{rlr}{\mathbf{R}}_1& =(4,4),& \\ {\mathbf{R}}_2& =(-4,4),& \\ {\mathbf{R}}_3& =(-4,-4),& \\ {\mathbf{R}}_4& =(4,-4),& \end{array}$$

$$H_{a}=8\epsilon_{b}N\left[s^{2}n_{B}+w^{2}(1-n_{B})\right]\\ +\sum_{\mathbf{r}_{\text{e}}}\left\{\epsilon_{a}a^{\dagger}({\bf r}_{\text{e}})a({\bf r}_{\text{e}})+\epsilon_{a}a^{\dagger}({\bf r}_{\text{e}}+{\bf R}_{1})a({\bf r}_{\text{e}}+{\bf R}_{1})\right.\\ \left.+\epsilon_{c}c^{\dagger}({\bf r}_{\text{e}}+{\bf L})c({\bf r}_{\text{e}}+{\bf L})\right.\\ \left.+\epsilon_{c}c^{\dagger}({\bf r}_{\text{e}}+{\bf L}+{\bf R}_{1})c({\bf r}_{\text{e}}+{\bf L}+{\bf R}_{1})\right\}\\ +gsw(1-2n_{B})\sum_{\alpha}\bigg{[}\left(e^{-i\varphi_{\alpha}}\sum_{{\bf r}_{\text{e}}}a({\bf r}_{\text{e}})a({\bf r}_{\text{e}}+{\bf R}_{\alpha})+h.c\right)\\ +\left(e^{-i\varphi_{\alpha}}\sum_{{\bf r}_{\text{e}}}c({\bf r}_{\text{e}}+{\bf L})c({\bf r}_{\text{e}}+{\bf L}+{\bf R}_{5-\alpha})+h.c\right)\bigg{]},$$

$$a_{\text{e}}(\mathbf{k})$$

$$a_{\text{o}}(\mathbf{k})$$

$$c_{\text{e}}(\mathbf{k})$$

$$c_{\text{o}}(\mathbf{k})$$

$$H=8\epsilon_{b}N\left[s^{2}n_{B}+w^{2}(1-n_{B})\right]+\epsilon_{a}\sum_{\mathbf{k}}a^{\dagger}_{\text{e}}(\mathbf{k})a_{\text{e}}(\mathbf{k})+\\ \epsilon_{a}\sum_{\mathbf{k}}a^{\dagger}_{\text{o}}(\mathbf{k})a_{\text{o}}(\mathbf{k})+\epsilon_{c}\sum_{\mathbf{k}}c^{\dagger}_{\text{e}}(\mathbf{k})c_{\text{e}}(\mathbf{k})+\epsilon_{c}\sum_{\mathbf{k}}c^{\dagger}_{\text{o}}(\mathbf{k})c_{\text{o}}(\mathbf{k})\\ +gsw(1-2n_{B})\sum_{\mathbf{k}}\bigg{[}\left(a_{\text{e}}(\mathbf{k})a_{\text{o}}(\mathbf{-k})V(\mathbf{k})+h.c\right)\\ +(c_{\text{e}}(\mathbf{k})c_{\text{o}}(\mathbf{-k})\tilde{V}(\mathbf{k})+h.c)\bigg{]},$$

$$V(\mathbf{k})$$

$$\tilde{V}(\mathbf{k})$$

$$a_{\text{e}}(\mathbf{k})$$

$$a_{\text{o}}(-\mathbf{k})$$

$$c_{\text{e}}(\mathbf{k})$$

$$c_{\text{o}}(-\mathbf{k})$$

$$u(k{)}^2+v(k{)}^2$$

$$p(k{)}^2+q(k{)}^2$$

$$\cos [{\varphi }_V(\mathbf{k})+{\varphi }_a(\mathbf{k})]$$

$$\cos [{\varphi }_{\tilde{V}}(\mathbf{k})+{\varphi }_c(\mathbf{k})]$$

$$E=8\epsilon_{b}N\left[s^{2}n_{B}+w^{2}(1-n_{B})\right]\\ +2\epsilon_{a}\sum_{\mathbf{k}}\left\{u^{2}(\mathbf{k})n_{A}(\mathbf{k})+v^{2}(\mathbf{k})[1-n_{A}(\mathbf{k})]\right\}\\ +2\epsilon_{c}\sum_{\mathbf{k}}\left\{p^{2}(\mathbf{k})n_{C}(\mathbf{k})+q^{2}(\mathbf{k})[1-n_{C}(\mathbf{k})]\right\}\\ +2gsw\left(1-2n_{B}\right)\sum_{\mathbf{k}}\bigg{[}u(\mathbf{k})v(\mathbf{k})(1-2n_{A}(\mathbf{k}))|V(\mathbf{k})|\bigg{.}\\ \bigg{.}+p(\mathbf{k})q(\mathbf{k})(1-2n_{C}(\mathbf{k}))]|\tilde{V}(\mathbf{k})|\bigg{]},$$

$$n_A(\mathbf{k})=\frac{1}{\exp \frac{{ϵ}_A(\mathbf{k})}{T}+1},$$

$$n_C(\mathbf{k})=\frac{1}{\exp \frac{{ϵ}_C(\mathbf{k})}{T}+1},$$

$${ϵ}_A(\mathbf{k})$$

$${ϵ}_C(\mathbf{k})$$

$$\epsilon_{B}=\frac{g^{2}}{16N}(1-2n_{B})\sum_{\mathbf{k}}\bigg{[}\frac{(1-2n_{A}(\mathbf{k}))}{\epsilon_{A}(\mathbf{k})}|V(\mathbf{k})|^{2}\\ +\frac{(1-2n_{C}(\mathbf{k}))}{\epsilon_{C}(\mathbf{k})}|\tilde{V}(\mathbf{k})|^{2}\bigg{]}.$$

$${ϵ}_A(\mathbf{k})$$

$${ϵ}_C(\mathbf{k})$$

$$(1-2n_B)$$

$$\begin{array}{r}{T}_c=\frac{g^2}{8}[\left(\frac{\exp (|{ϵ}_a|/T_c)-1}{\exp (|{ϵ}_a|/T_c)+1}\right)\frac{1}{|{ϵ}_a|}\\ +\left(\frac{\exp (|{ϵ}_c|/T_c)-1}{\exp (|{ϵ}_c|/T_c)+1}\right)\frac{1}{|{ϵ}_c|}],\end{array}$$

$$\Delta =\sqrt{{ϵ}_a^2+g^2(1-2n_B{)}^2},$$

$${ϵ}_B=\sqrt{{ϵ}_b^2+\frac{g^2}{64}{\tilde{C}}^2},$$

Full text

††preprint: APS/123-QED

I Introduction

Interplay between superconductivity and the onset of electronic spin and charge modulations in cuprate superconductors remains one of the intriguing and unresolved issues in the field of high-temperature superconductivity. Manifestations of electronic modulations are reported in a broad doping range for several families of cuprates — most noticeably around the doping level of 1/8Tranquada et al. (1995); Yamada et al. (1998); Hoffman et al. (2002); McElroy et al. (2003); Vershinin et al. (2004); Hanaguri et al. (2004); Abbamonte et al. (2005); McElroy et al. (2005); Komiya et al. (2005); Wise et al. (2008); da Silva Neto et al. (2014); Comin et al. (2015a); Comin and Damascelli (2016); Peng et al. (2018).

For 1/8-doped lanthanum cuprates, the modulated structure is widely believed to exhibit one-dimensional pattern often referred to as “stripes” Tranquada et al. (1995); Yamada et al. (1998). Yet the principal aspects of the same experimental evidence are also consistent with the possibility of two-dimensional modulations called “checkerboards” Tranquada et al. (1999); Mitsen and Ivanenko (2004); Robertson et al. (2006); Tranquada (2007); Fine (2004, 2007, 2011); Brandenburg and Fine (2013). The experiment-based arguments discriminating between stripes and checkerboards in 1/8-doped lanthanum cuprates are, at present, rather indirect. At the same time, the issue cannot be resolved on purely theoretical grounds, because it requires accuracy of the calculations of the ground-state energy not achievable by first-principles theories. A particularly focused effort to investigate the consequences of the checkerboard scenario was made in Ref.Fine (2004). That analysis was based on a particular kind of checkerboard called “grid”. Later, the grid checkerboard was shown to be inconsistent with the results of spin-polarized neutron scattering experiment of Ref.Christensen et al. (2007). This experiment, however, did not rule out another version of a checkerboard representing a two-dimensional arrangement of spin vortices Fine (2007) shown in Fig.1. Somewhat similar noncollinear spin textures were also considered in Refs. Seibold (1998); Berciu and John (1999); Timm and Bennemann (2000); Koizumi (2008); Wilson (2009); Azzouz et al. (2010). Recently an analogous superstructure called ”spin-vortex crystal” was proposed to exist in iron-based superconductorsAvci et al. (2014); Böhmer et al. (2015); O’Halloran et al. (2017); Meier et al. (2018).

Spin-vortex checkerboard in the context of cuprates was introduced in Ref.Fine (2007) . Its various properties were analyzed in Refs.Fine (2007, 2011); Brandenburg and Fine (2013); Dolgirev and Fine (2017). So far, however, this analysis has not touched the superconducting properties. In Section II of the present article, we introduce a microscopic model aimed at describing superconductivity in the background of spin-vortex checkerboard. The model is a generalization of another one proposed in Ref.Fine (2004) for the grid checkerboard. In Sections III and IV, we find the mean-field solution of the generalized model in the same way as done in Ref.Fine (2004). However, our subsequent analysis in Section V differs from that of Ref.Fine (2004) in one significant respect, namely, we find that both the original model and its generalization made in Section II do not induce non-zero superfluid stiffness. (In this regard, we correct an error made in Ref.Fine (2004).) We further show in Section V.2 that non-zero superfluid stiffness emerges once arbitrarily small hopping terms are included in the modeling.

Experimentally, the onset of static spin modulations in 1/8-doped lanthanum cuprates — spin vortices or stripes — largely suppresses three-dimensional superconductivity but appears to coexist with two-dimensional superconductivity Li et al. (2007); Berg et al. (2007); Tranquada et al. (2008); Li et al. (2018). In Section VI, we show that, similar to the stripe scenario, the suppression of the three dimensional superconductivity in spin-vortex scenario can be explained by the displacement of the modulation patterns in the adjacent CuO2 planes.

As further argued in Section VI, the generic setting in cuprates beyond 1/8-doped lanthanum family, possibly, involves the fluctuating counterpart of the static spin-vortex texture used in the present work. These fluctuations are likely caused by system’s proximity to the threshold of electronic phase separation Fine and Egami (2008). They are expected to couple spin, charge and lattice degrees of freedom — see e.g.Egami et al. (2010). In this respect, our model illustrates the potential of the general two-component scenariosRanninger and Domański (2010); Pawłowski et al. (2010) in the limit of initially localized components for describing the superconductivity in cuprates. In such scenarios, the first component represents unpaired fermions, while the second component represents preformed fermionic pairs. The spin-vortex checkerboard just suggests us the coupling connectivity between the low-energy fermionic states. We assume that in other cuprate families, the inhomogeneous patterns have a two-dimensional character, checkerboard-like or more disordered, where regions of stronger and weaker antiferromagnetic correlations alternate with other, but the size of these regions may be different for different dopings, and the doping dependence of that size may be different for different cuprate families. This would be consistent with rather diverse phenomenology summarized in Ref.Comin and Damascelli (2016). The dimensionality of charge modulations was, in particular, discussed for the yttrium family of curates in Refs.Comin et al. (2015b); Fine (2016); Comin et al. (2016).

Finally, in Section VII, we compare the temperature evolution of the energy gap obtained from the mean-field solution of our model with experiments in bismuth cuprates. Several important but technical calculations behind the reported results are placed in the Appendices.

II Microscopic model: principal terms

The model to be considered has two different kinds of fermionic states physically located in magnetic and non-magnetic parts of the underlying spin texture. The general reasoning for constructing the model is the same as in Ref.Fine (2004). Namely, the entire texture is divided into plaquets having different kinds of spin background, and then, for each plaquet, only one-particle fermionic states closest to the chemical potential are retained. Given that plaquets are rather small, it can be estimatedFine (2004) that the spacing of one-particle energies within each plaquet is of the order of 40 meV, which implies that, for temperatures much smaller than 400K, it is appropriate to retain only the levels closest to the chemical potential.

As shown in Fig.1, spin vortex checkerboard can be represented as a square superlattice of $8d\times 8d$, where $d$ is the period of the underlying square lattice of Cu atoms. We denote the total number of such unit cells in the system as $N$. Each unit cell in Fig.1 is further divided into four spin-polarized plaquets and four spin unpolarized plaquets. We expect that the lowest one-particle states in spin-polarized plaquets are non-spin-degenerate, and hence we include exactly one state per plaquet. We refer to two of the resulting states as “a-states” and to the remaining two as “c-states”. Two different kinds of a-states are distinguished by index $\eta=\pm 1$, and c-states — by index $\zeta=\pm 1$. Two a-states or two c-states with different values of $\eta$ or $\zeta$ respectively are expected to have orthogonal spin wave functions that can be obtained from each other by spin inversion. The lowest-energy states of spin-unpolarized plaquets around the cores of spin vortices are assumed to be spin-degenerate. We, therefore, place two fermionic states on each such plaquet with spins “up” or “down” along any chosen direction. We call them “b-states”. Since the spin texture contains four nonequivalent kinds of spin-vortex cores, we distinguish the corresponding b-states by index $\alpha=1,2,3,4$ and by spin index $\uparrow$ or $\downarrow$ — see Fig. 2.

We now construct a low-energy Hamiltonian similar to that of Ref.Fine (2004). We expect that charge carriers in the background of spin-vortex checkerboard are heavily dressed by antiferromagnetic fluctuations, which should lead to relatively large effective masses. Charge carriers are also supposed to experience strong Coulomb repulsion, which should further prevent them from hopping. Therefore, in the zeroth order, we neglect single-fermion hopping between checkerboard plaquets. (Hopping will be treated in Section V.2 as a small correction.) As far as the interaction terms are concerned, we use the same ”heavy-mass/Coulomb-repulsion” argument to include in the model only those terms that do not change the center-of-mass/center-of-charge positions of the two fermions involved. We also neglect a possible ”on-site” coupling of two $b$-states on the same plaquet — in our preliminary analysis, this term would not change the essential features of the solution, but we still plan to investigate it elsewhere. The above assumptions amount to a relatively crude overall approximation, which should nevertheless allow us to capture the principal behavior of the variational solution and, at the same time, avoid introducing too many adjustable parameters.

We are, finally, left with terms representing on-site energies $\epsilon_{a}$, $\epsilon_{b}$ and $\epsilon_{c}$ (with $\epsilon_{a}=\epsilon_{c}$) and with the following two-kinds of effective interaction terms, namely: two a-states or two c-states adjacent to a given spin vortex core making transitions to the two b-states inside the core or vice versa. The resulting Hamiltonian is:

[TABLE]

where $g$ is the interaction constant, $\epsilon_{a}$, $\epsilon_{b}$ and $\epsilon_{c}$ are on-site energies defined with respect to the chemical potential $\mu$, which we set equal to zero, index $i$ labels unit cells depicted in Fig 2, and indices $\eta$, $\alpha$ and $\zeta$ label the plaquets within the unit cell as illustrated in Fig 2. Following Ref.Fine (2004), whenever the specific value of subscripts $\eta$ or $\zeta$ is fixed, as is the case in the interaction term of Hamiltonian (1), we use subscript “e” for $\eta=1$ or $\zeta=1$ referring to the corresponding plaquets as “even” and subscript “o” for $\eta=-1$ or $\zeta=-1$ referring to the corresponding plaquets as “odd”. Double-subscripts notations such as $a_{j\text{e}_{[i,\alpha]}}$ imply that a-states labeled as $\{j\text{e}\}$ must be adjacent to the b-states labeled as $\{i,\alpha\}$.

If all terms containing c-states are removed from Hamiltonian (1), the result would be exactly equivalent to the Hamiltonian considered in Ref Fine (2004). Since c-states do not directly couple to a-states, and since c-states have the same connectivity with the b-states as a-states (but shifted), the mean-field solutions of the two models are very similar with the only difference being that b-states now experience mean field from both a-states and c-states, which, in turn, makes that mean field two times larger, and, as a result, the value of the superconducting transition temperature becomes modified.

Since the entire mean-field solution has nearly the same structure and logic as that of Ref. Fine (2004), below we only include the formal structure of the derivation and the results, leaving the justification mostly to Ref. Fine (2004).

As explained later in Section V.2, the system described by Hamiltonian (1) will not be superconducting, because it will have zero superfluid stiffness. In order to make it superconducting, arbitrarily small hopping terms will need to be added. However, the model based on Hamiltonian (1) is easier to handle, and its solution already captures the important aspects of the resulting superconducting phase, such as the energy gap.

III Bogoliubov transformations

In the model considered, each of the fermionic states couples to relatively few other states, which makes a mean-field solution rather approximate. We, nevertheless, assume that it gives at least the right qualitative picture of model’s behavior. The first step of this solution is to introduce the Bogoliubov transformation for b-states within the same plaquet:

[TABLE]

where $s$ and $w$ are positive real numbers satisfying a constraint arising from canonical fermionic anticommutation relations

[TABLE]

and $\varphi_{\alpha}$ are the transformation phases, which are to be determined later by minimizing system’s energy.

Substituting Bogoliubov transformation for b-states in (1) and keeping only the thermal averages of terms that do not change the occupations of Bogoliubov $B$-states, we obtain partially averaged Hamiltonian,

[TABLE]

where $\epsilon_{B}$ is the energy of $B$-quasiparticles and

[TABLE]

is their occupation number. We set $k_{{\text{B}}}=1$.

As explained in Ref.Fine (2004), in order to assure proper fermionic anticommutation relations for the Bogoliubov counterparts of a- and c-states, the Bogoliubov transformation for these states should be made in the quasimomentum space. Therefore, we need to rewrite the Hamiltonian in terms of the real-space Fourier transforms for a- and c-operators. To do the Fourier transforms, we first change the notations from $a_{i\eta}$, $c_{i\eta}$ to the notations $a(\mathbf{r})$, $c(\mathbf{r})$, where $\mathbf{r}$ indicates the position of the center of the respective plaquet. We also need to define the following vectors — all in units of underlying crystal lattice period $d$:

[TABLE]

which connects an even a-state with an adjacent even c-state, and

[TABLE]

which connect an even a-state with four adjacent odd a-states. The subscript $\alpha$ in ${\bf R}_{\alpha}$ is chosen such that $\alpha$th *b-*states are located between the pairs a-states connected by vector ${\bf R}_{\alpha}$ originated from an even a-state. Now, we define the position of each unit cell by the position ${\bf r}_{\text{e}}$ of an even a-state within this cell. Therefore, even a-states are located at a set of positions $\{{\bf r}_{\text{e}}\}$, odd a-states at $\{{\bf r}_{\text{e}}+{\bf R}_{1}\}$, even c-states at $\{{\bf r}_{\text{e}}+{\bf L}\}$, and odd c-states at $\{{\bf r}_{\text{e}}+{\bf L}+{\bf R}_{1}\}$. Finally, we rewrite the Hamiltonian (4) as follows,

[TABLE]

We now explicitly write separate Fourier transforms for even and odd a- and c-states as follows:

[TABLE]

Since the superlattice periods for each of the above four kinds of states are the same, the sets of wave vectors $\{\mathbf{k}\}$ are also the same, even though the corresponding states are shifted with respect to each other in real space. Substituting these transformations to (8), we obtain

[TABLE]

where

[TABLE]

and

[TABLE]

Bogoliubov transformations for a- and c-states can now be defined as,

[TABLE]

where $u(\mathbf{k}),$ $v(\mathbf{k})$ and $p(\mathbf{k}),$ $q(\mathbf{k})$ are the real-valued coefficients for a- and c-states respectively, subjected to a constraint arising from the fermionic canonical commutation relations for $A$- and $C$-operators:

[TABLE]

and $\phi_{a}$ and $\phi_{c}$ are complex phases — all to be found by the energy minimization.

We now complete the following steps: (i) substituting the above canonical transformation for a- and c-states (16-19) into Hamiltonian (13), then (ii) obtaining the energy of the system by summing over the thermal averages of the diagonal terms — the result is given in Appendix A by Eq.(76)— and then (iii) minimizing the resulting energy with respect to the choice of phases $\phi_{a}(\mathbf{k})$ and $\phi_{c}(\mathbf{k})$. As explained in Appendix A, these steps lead to conditions:

[TABLE]

were $\phi_{V}(\mathbf{k})$ and $\phi_{\tilde{V}}(\mathbf{k})$ are the complex phases of $V(\mathbf{k})$ and $\tilde{V}(\mathbf{k})$ respectively, which, in turn, depend on phases $\{\varphi_{\alpha}\}$. [Phases $\phi_{a}(\mathbf{k}$) and $\phi_{c}(\mathbf{k})$ do not need to be obtained explicitly, because they will not enter any quantity computed in this paper.] With the above conditions, the expression for the energy of the system becomes:

[TABLE]

where,

[TABLE]

[TABLE]

are the Bogoliubov quasiparticle occupation number and $\epsilon_{A}(\mathbf{k})$, $\epsilon_{C}(\mathbf{k})$ are their energies obtained in the next section.

IV Single-particle excitations, energy gap and the critical temperature

As argued in Ref Fine (2004), the chemical potential of the system is, likely, to coincide with either $\epsilon_{b}$ or $\epsilon_{a}$ (same as $\epsilon_{c}$), which, given our convention $\mu=0$, means that either $\epsilon_{b}=0$ or $\epsilon_{a}=\epsilon_{c}=0$. Below, we treat these two cases separately, referring to them as “Case I” and “Case II” respectively, and also refer to the case of $\epsilon_{a}=\epsilon_{b}=\epsilon_{c}=0$ as “critical”.

The coefficients of the Bogoliubov transformations for both Cases I and II are obtained in Appendices B and C, respectively, by minimizing the total energy (24) at fixed quasiparticle occupation numbersLifshitz and Pitaevskii (1980). We then substitute those coefficients back to Eq (24) and obtain the energy of a Bogoliubov quasiparticle by taking derivative of the total energy (24) with respect to the quasiparticle occupation numbers $n_{A}(\mathbf{k})$, $n_{C}(\mathbf{k})$ or $n_{B}$.

IV.1 Case I: $\epsilon_{b}=0$

In Case I, the above procedure gives the following quasiparticle energies:

[TABLE]

[TABLE]

[TABLE]

The mean-field approach now requires finding a nontrivial solution for $\epsilon_{A}(\mathbf{k})$, $\epsilon_{C}(\mathbf{k})$, and $\epsilon_{B}$ from Eqs.(5, 25, 26,27, 28, 29). In general, it can only be done numerically, but one can also obtain a closed analytical equation for the critical temperature $T_{c}$ using the fact that, near the transition, the ordered state is close to the normal states, which allows one to use the limits

[TABLE]

This gives

[TABLE]

from which the mean-field $T_{c}$ can be obtained numerically.

As explained in Ref.Fine (2004), the energy densities of $A$- and $C$-states described by Eqs. (27) and (28) have Van Hove singularities located in the both cases at the value

[TABLE]

corresponding to $|V(\mathbf{k})|=|\tilde{V}(\mathbf{k})|=2$. We refer to $\Delta$ as the “energy gap”. As $T\rightarrow T_{c}$, $\Delta$ approaches not zero but $|\epsilon_{a}|$, which we associate with the pseudogap.

IV.2 Case II: $\epsilon_{a}=\epsilon_{c}=0$

Following the same procedure as for Case I, we obtain:

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

Detailed calculations can be found in Appendix C. In the grid modelFine (2004), the $c-$states are absent, hence, $C_{c}=0$, and the value of $C_{a}$ is reported to be $0.958$.

The same approach as in Case I now gives the critical temperature

[TABLE]

and the gap parameter

[TABLE]

associated with the Van Hove singularity for $A$- and $B$- states located at $|V(\mathbf{k})|=|\tilde{V}(\mathbf{k})|=2$.

IV.3 Temperature evolution of the energy gap

In Fig. 3, we present temperature dependences of energy gaps for Cases I and II given by Eqs.(31) and (39) respectively. These dependencies were obtained by the numerical solution of the system of equations Eqs.(5, 25, 26,27, 28, 29) for Case I, or Eqs.(5, 25, 26, 32, 33, 34) for Case II.

The families of plots for Cases I and II are connected through the critical case $\epsilon_{a}=\epsilon_{b}=\epsilon_{c}=0$, which is represented by the thick red line. This case corresponds to the ratio $\Delta(0)/T_{c}=2\sqrt{2}\approx 2.82$. Plots above the critical-case line correspond to Case I: at $T=T_{c}$, they all end at nonzero values $\Delta(T_{c})=\epsilon_{a}$. Plots below the critical-case line correspond to Case II: they all have $\Delta(T_{c})=0$, and, moreover approach closely the canonical BCS dependence for $\epsilon_{b}/g\rightarrow\infty$.

Thus, if the assumptions of the present model are valid, the critical-case ratio $\Delta(0)/T_{c}=2\sqrt{2}$ signifies the transition from the conventional behavior $\Delta(T_{c})=0$ for $\Delta(0)/T_{c}<2\sqrt{2}$ to unconventional behavior $\Delta(T_{c})\neq 0$ for $\Delta(0)/T_{c}>2\sqrt{2}$. The value of $\Delta(0)/T_{c}=2\sqrt{2}$ for the critical case makes important quantitative difference from the critical case result $\Delta(0)/T_{c}=4$ for the grid-based model of Refs.Fine (2004, 2005), which involved only a- and b-states. Such a difference was to be expected, because the coupling between b- and c-states in the present model leads to additional energy advantage for the superconducting state and hence higher superconducting transition temperature for the same value of the coupling constant $g$.

V Superconducting properties

V.1 Anomalous correlation functions

Bogoliubov transformations (2,16-19) can be used to obtain the following anomalous correlation functions for $T<T_{c}$:

[TABLE]

where ${\mathbf{r}}_{i\alpha}$ is the position of $\alpha$th b-state in the ith unitcell, and $\delta({\mathbf{r}}_{i\alpha}-{\mathbf{r}}_{j\beta})$ is defined as Kronecker delta on the discrete superlattice. Two different components of the SC order parameter corresponding to a, c, and b-states are the correlation functions (40), (41), and (42) respectively.

The anomalous averages for a and c-components, in real space, can be written as following,

[TABLE]

where $\mathbf{r}_{\text{e}_{[i,\alpha]}}$ and $\mathbf{r}_{\text{o}_{[j,\beta]}}$ are, respectively, the positions of an even ($\eta=1$) and odd ($\eta=-1$) a-states adjacent to the $b$-states labelled by indices $[i,\alpha]$ and $[j,\beta]$, and, likewise, $\mathbf{r}^{\prime}_{\text{e}_{[i,\alpha]}}$ and $\mathbf{r}^{\prime}_{\text{o}_{[j,\beta]}}$ are, respectively, the positions of an even ($\zeta=1$) and odd ($\zeta=-1$) c-states adjacent to the same $b$-states. The anomalous averages given by Eqs.(43) and (44) have non-zero values only when their two arguments correspond to states of different kind (i.e. even and odd). They can be expressed as

[TABLE]

where $\Psi_{a}({\mathbf{k}})$ and $\Psi_{c}({\mathbf{k}})$ are given by Eq.(40) and (41). These anomalous averages also obey the following relations:

[TABLE]

Relations (47) and (48) follow from the fermionic anticommutation rule. Variables $\tilde{\mathbf{r}}_{\text{e}}$, $\tilde{\mathbf{r}}_{\text{o}}$, $\tilde{\mathbf{r}}^{\prime}_{\text{e}}$ and $\tilde{\mathbf{r}}^{\prime}_{\text{o}}$ in Eqs. (49) and (46)represent different positions of the same kind as $\mathbf{r}_{\text{e}}$, $\mathbf{r}_{\text{o}}$, $\mathbf{r}^{\prime}_{\text{e}}$ and $\mathbf{r}^{\prime}_{\text{o}}$ respectively.

The coherence length of the order parameters $\Psi_{a}(\mathbf{r}_{\text{e}}$, $\mathbf{r}_{\text{o}})$ and $\Psi_{c}(\mathbf{r}_{\text{e}}$, $\mathbf{r}_{\text{o}})$ should be inversely proportional to the characteristic $\mathbf{k-}$space scale of $V({\mathbf{k}})$ and $\tilde{V}({\mathbf{k}})$, respectively. The examination of Eqs.(III) and (III) reveals that this characteristic scale is $\pi/l$, where $l=8d$ . Therefore, the coherence length associated with $\Psi_{a}(\mathbf{r}_{\text{e}},\mathbf{r}_{\text{o}})$ and $\Psi_{c}(\mathbf{r}^{\prime}_{\text{e}},\mathbf{r}^{\prime}_{\text{o}})$ can be estimated as the modulation period $l$.

The coherence length associated with $\Psi_{b}$ is equal to zero, which means that only b-states located on the same plaquet form coherent pairs.

We now obtain three quantities that characterize the short-range pair correlations:

[TABLE]

In Case I, the explicit expression for $\Psi_{b(i\alpha)}$ can be obtained by substituting the values of $s$ and $w$ given by Eqs.(85) into Eq.(42) for ${\mathbf{r}}_{i\alpha}={\mathbf{r}}_{j\beta}$, which gives

[TABLE]

One can then obtain both $\Psi_{a(i,\alpha)}$ and $\Psi_{c(i,\alpha)}$ by making use of the fact that

[TABLE]

where $E_{\hbox{\footnotesize int}}^{a}$ and $E_{\hbox{\footnotesize int}}^{c}$ are the interaction parts of the energy (24):

[TABLE]

[TABLE]

After $E_{\hbox{\footnotesize int}}^{a}$ and $E_{\hbox{\footnotesize int}}^{c}$ are evaluated with the help of Eqs.(77-80), one can use Eqs.(54,55,56) to obtain:

[TABLE]

[TABLE]

In Case II, the expressions analogous to (59,60,54) are

[TABLE]

V.2 Superfluid phase stiffness and emergence of superconductivity

The model introduced in Section II is defined in terms of localized fermionic states with a Hamiltonian that does not contain hopping terms or interaction terms changing the center of mass of the particles. As a result, the mean-field solution obtained in the preceding sections has zero superfluid stiffness and hence cannot sustain superconductivityerr . In order to illustrate this, let us consider the following phase transformation of fermionic operators:

[TABLE]

where $\theta(\mathbf{r})$ is the position-dependent quantum phase, subscript $\sigma$ represents the spins of the $b$-states. All other former subscripts “$i$”,“e”,“o”, and “$\alpha$” are unambiguously determined by the position $\mathbf{r}$. Technically, the absence of the superfluid stiffness originates from the fact that the Hamiltonian of our model is invariant under transformation (64-66) with

[TABLE]

where $\mathbf{G}$ is an arbitrary vector. Therefore, there exists a continuous set of mean-field solutions of equal energy, which can be obtained from each other by applying the above set of transformations.

Despite having zero superfluid stiffness, the mean-field solution of the model based on Hamiltonian (1) massively suppresses phase fluctuations of the pair correlations that eventually lead to superconductivity. This happens for two reasons. First, the solution fixes the relative phase of the fermionic pairs occupying $a$- and $c$-states with respect to the $b$-states. Second, each $a$- or $c$-state participates in a fermionic pair with four adjacent states of the same kind but opposite spins, which implies that the phase fluctuations other than the linear ones having form (67) across the entire system do cost energy. The above constraints leave the simultaneous phase fluctuations of $a$-, $b$- and $c$-states with $\theta(\mathbf{r})$ given by Eq.(67) as a rather narrow channel, through which superconductivity is suppressed. As we illustrate below, this channel is removed and the superconductivity is recovered once arbitrarily small hopping terms are added to the model.

In general, the hopping terms can connect two nearest regions ($a\leftrightarrow b$ and $c\leftrightarrow b$), two next-to-nearest regions ($a\leftrightarrow c$ or $b\leftrightarrow b$), two next-to-next-nearest regions, etc. Because our goal in this section is to describe the main qualitative effects of hopping, we would like to consider the terms whose impact is simplest to calculate. Thus we include into consideration the hopping terms between two nearest $a$-states with the same spins and between two nearest $c$-states with the same spins:

[TABLE]

where

[TABLE]

Here $t$ is the hopping parameter, and $\mathbf{r}_{a}$ and $\mathbf{r}_{c}$ run over all $a$- and $c$-plaquets respectively.

For the model with Hamiltonian (68), the variational mean-field solution can be constructed in the same way as in the preceding sections — with only one difference, namely, the energies of unpaired fermions occupying $a$- and $c$-states now form a band with the $\mathbf{k}$-dependence:

[TABLE]

Here we again assume $\epsilon_{a}=\epsilon_{c}$. The Bogoliubov transformations are still defined by Eqs.(2,16,17,18,19). As shown in Appendix D, for small $t$, the $\mathbf{k}$-dependence of bare energies (70) changes the values of the Bogoliubov coefficients $u(\mathbf{k})$, $v(\mathbf{k})$, $p(\mathbf{k})$ and $q(\mathbf{k})$ only very little. The variational wave function of the system can still be expressed as

[TABLE]

and the resulting equation for the critical temperature changes also only slightly. We note, however, that the hopping part of Hamiltonian $H_{t}$ introduces additional constraints on phases $\varphi_{i}$, which are important for computing the superfluid stiffness. These constraints are derived in Appendix E.

In order to calculate the phase stiffness of the above mean-field solution, we introduce an auxiliary wave function $|\tilde{\Psi}\rangle$ obtained from $|\Psi\rangle$ with the help of transformation (64-66). Namely, in the right-hand-side of Eq.(71), we use the inverse Fourier transform for $a$- and $c$-operators defined by Eqs.(9-12), then replace all lattice subscripts with the position argument $\mathbf{r}$, and, finally, substitute the creation operators $a^{\dagger}(\mathbf{r})$, $b^{\dagger}_{\sigma}(\mathbf{r})$, and $c^{\dagger}(\mathbf{r})$ with their respective counterparts $\tilde{a}^{\dagger}(\mathbf{r})$, $\tilde{b}^{\dagger}_{\sigma}(\mathbf{r})$, and $\tilde{c}^{\dagger}(\mathbf{r})$ defined by Eqs.(64-66). The state represented by the wave function $|\tilde{\Psi}\rangle$ carries supercurrent. It is the lowest-energy state of the Hamiltonian $\tilde{H}^{\prime}$ which has the same form as $H^{\prime}$ but in terms of $\tilde{a}(\mathbf{r})$, $\tilde{b}_{\sigma}(\mathbf{r})$ and $\tilde{c}(\mathbf{r})$.

We define the zero-temperature phase stiffness $J$ in the following way:

[TABLE]

where $\theta$ is the phase appearing in transformations (64-66), and the sum in the left-hand side is taken over all $2L\times 2L$ unit cells of the spin vortex checkerboard. (In the continuum limit, this sum is replaced by $\displaystyle\int\frac{J}{2}(\mathbf{\nabla}\theta)^{2}d^{2}{\mathbf{r}}$.)

In Appendix F, we obtain the following expression for the value of $J$ corresponding to the Hamiltonian (68) in the case where $\epsilon_{b}=\mu\equiv 0$:

[TABLE]

We explicitly evaluated this expression in Case I under additional assumption $|t|\ll|g|\ll|\epsilon_{a}|$. In the leading order, the result is:

For $\displaystyle\frac{4|t\epsilon_{a}|}{g^{2}}\leqslant 1$,

[TABLE]

while, for $\displaystyle\frac{4|t\epsilon_{a}|}{g^{2}}\geqslant 1$,

[TABLE]

In this work, we do not attempt to calculate the superconducting properties at finite temperatures. Below we only discuss qualitative considerations about the resulting thermodynamic behavior:

Without hopping, the phase transition that we obtained in the mean-field approximation for Hamiltonian (1) turns into a crossover, whose sharpness is limited by the phase fluctuations of the order parameter. The addition of a small hopping part (69) to Hamiltonian (1) introduces into the system finite superfluid stiffness and the ability to carry current at sufficiently low temperatures. In such a setting, we expect the temperature $T_{ph}$, at which in-plane superfluid stiffness appears, to be significantly smaller than the the mean-field transition temperature $T_{c}$. Strong pair fluctuations should thus be present in the temperature range between $T_{ph}$ and $T_{c}$. However, in a three-dimensional setting of stacked layers and in the presence of Coulomb interaction between charges, the above fluctuations should be somewhat suppressed.

Once the overall phase of the superconducting solution is stabilized, the strong spatial dependence of the anomalous correlation functions described in Section V.1 can lead to an unorthodox interpretation of the experiments sensitive to the sign of the superconducting phase in cuprates — see the discussion in Ref.Fine (2004). It also implies that the standard factorization of the Cooper-pair wave functions in terms of the center-of-mass dependence and the relative-coordinate dependence is not applicable, which also means that the classification of Cooper pairs in terms of spin singlets and spin triplets does not apply. This in turn, may lead to non-trivial spin susceptibility of the superconducting state, which requires a separate investigation extending beyond the scope of the present paper. Here, we would only like to remark that the current experimental knowledge of the spin susceptibility of superconducting cuprates is largely based on the Knight shift nuclear magnetic resonance (NMR) experimentsSlichter (2006), which are subject to a number of assumptions about the chemical shifts and about certain accidental cancellation of hyperfine coefficients. According to a recent reviewHaase et al. (2017), the overall Knight-shift phenomenology in cuprates is still evolving, and this may lead to a re-examination of the standard NMR interpretations adopted in the past.

VI Theoretical scenarios for cuprates

VI.1 1/8-doped lanthanum cuprates

The model presented in this article is primarily relevant to the properties of 1/8-doped lanthanum cuprates where the static spin and charge order is stabilized, such as La1.875Ba0.125CuO4 (1/8-doped LBCO). On cooling, the latter system exhibits static charge order at temperatures below $T_{CO}=54$K. On further cooling at zero magnetic field, 1/8-doped LBCO sample investigated in Ref.Li et al. (2007) exhibits the onset of static spin order and a simultaneous sharp drop of in-plane resistivity by a factor of 10 at temperature $T_{SO}\approx 40$K. At yet lower temperature, the in-plane resistivity decays to unmeasurably small values following the Berezinskii-Kosterlitz-Thouless-like Berezinskii (1971); Kosterlitz and Thouless (1973) scaling with characteristic temperature $T_{BKT}=16$K. Finally, at temperature $T_{c}=4$K, the sample exhibits bulk superconductivity. Another sample reported by the same group in Ref. Li et al. (2018) exhibits similar behavior as far as the initial resistivity drop and the bulk superconducting transitions are concerned, but, in the temperature range between 4K and 40K, the temperature dependence of resistivity appears to be sample-dependent. The latter non-universality, probably, indicates a “fragile” quantum mechanism of resistivity in the interval between the temperature of the initial resistivity drop and the temperature of the onset of bulk superconductivity.

A theoretical interpretation of the above phenomenology in the framework of the stripe scenario was presented in Refs.Berg et al. (2007, 2009). It was based on the idea of strongly suppressed Josephson coupling between the adjacent CuO2 layers caused by the in-plane change of the sign of the SC order parameter between the adjacent stripes and then by the mismatch of that sign between the adjacent layers.

A similar suppression of the Josephson coupling between adjacent layers can exist in the checkerboard scenario. As shown in Ref.Fine (2004) for the grid checkerboard and in the present work for the spin-vortex checkerboard, the sign of the SC order parameters of the $a$- (and $c$-) fermionic components is supposed to change in real space, remaining zero on average - thus falling under the definition of the pair-density wave proposed in Ref.Berg et al. (2009). The sign of SC order for the $b$-component may also change in real space but it may or may not lead the zero real-space average. If layers described by our model are stacked on the top of each other, then, due to the Coulomb repulsion, the $b$-plaquets are supposed to be located over the $a$- or $c$-plaquets - see Fig.4. This mismatch suppresses the Josephson tunneling between the adjacent planes, which can lead to different transition temperatures for the onsets of the two-dimensional and three-dimensional superconductivity. (We note here that the above-mentioned mismatch cannot be perfect, because the adjacent CuO2 planes in lanthanum cuprates are shifted by half a crystal period.)

Apart from the difference between stripes and checkerboards, another crucial difference between our proposal and that of Ref.Berg et al. (2009) is that the model of Ref.Berg et al. (2009) is based on one fermionic component forming Cooper pairs inside non-magnetic regions, while our model involves two fermionic components — one residing in the magnetic regions ($a$-/$c$-states) and the other one in non-magnetic regions ($b$-states), and, moreover, the Hamiltonian of our model is dominated by the center-of-mass-conserving interaction between the two components, which leads to anomalously low in-plane superfluid stiffness. As discussed at the end of the preceding section, our model is then supposed to exhibit a higher-temperature crossover associated with the appearance of anomalous fermionic averages (42)-(44) with fluctuating phases and the lower-temperature transition that leads to establishing the long-range phase coherence. We associate the above higher-temperature crossover with the initial resistivity drop in 1/8-doped LBCO at $T\sim 40$K, while the lower-temperature transition is assigned to the observed BKT-like crossover.

Now we address several issues that naturally arise in the context of the above scenario:

(i) Two-dimensional modeling versus three-dimensional character of real materials. In the present work, we only obtain a mean-field description of a two-dimensional model. On the one hand, it is known that long-range fluctuations and the proliferation of topological defects turn the transitions obtained from such a solution into a crossover. On the other hand, the CuO2 planes described by our model are stacked in a three-dimensional setting and, as a result, coupled to each other via Coulomb interaction, lattice strain, spin exchange and the hopping of charge carriers. This, in turn, suppresses both the long-range fluctuations and the formation of topological defects, thereby, making the mean-field solution adequate at least semi-quantitatively.

(ii) Reason for the resistivity drop to coincide with the mean-field transition temperature in our model. Let us consider Case IA of our model in the limit $g\ll\epsilon_{a}=\epsilon_{c}$ with small hopping between $a$-states and between $c$-states, as also done in Section V.2. (In terms of the concentration of charge carriers, Case IB with $g\ll|\epsilon_{a}|=|\epsilon_{c}|$ is, likely, closer to the situation in 1/8-doped LBCO, but the mathematics of this case is the same as that of Case IA, while the latter is more intuitive to discuss.) In Case IA, $T_{c}\ll\epsilon_{a}$ and, therefore, just above $T_{c}$, the occupations of the $a$- and $c$-states are exponentially small, which means large resistivity. Below $T_{c}$, the anomalous averages of $a$- and $c$- states (43,44) appear, which means that Cooper pairs can now occupy these states. Due to anomalously low superfluid stiffness the phase of these Cooper pairs fluctuates. Hence, instead of carrying a supercurrent, they carry normal current but the resulting state still has resistivity much smaller than that above $T_{c}$.

(iii) Coincidence between the mean-field transition temperature in our model and the onset of static spin modulations. In Ref.Li et al. (2007), the drop of resistivity in zero magnetic field was observed to coincide with the onset static spin modulations. Since we attribute the same resistivity drop to the mean-field transition in our model, the question arises: what can our model have to do with the onset of static spin modulations? In the model, we do not obtain the spin background but rather treat it as imposed externally. We can imagine, however, that, in a more complete theory, spin-background degrees of freedom can be included in the modeling together with the terms considered by us. In such a case, the energy gain due to the onset of superconducting correlations can simultaneously become a factor stabilizing the spin background, thereby making it static. We would then expect that the magnetic field that suppresses the onset of superconducting fluctuations should simultaneously suppress the onset of static spin order. Such a behavior can be contrasted with the weak sensitivity to magnetic fields of spin orders driven by the short-range exchange interactions.

VI.2 Other cuprates

Let us now turn to cuprates beyond doping 1/8 and beyond the lanthanum family. As shown in Ref.Fine and Egami (2008), cuprates are generically close to the threshold of Coulomb-frustrated phase separation into nanoregions of stronger antiferromagnetic correlations and lower density of charge carriers and nearly non-magnetic nanoregions attracting charge carriers. It is, therefore, reasonable to assume that nanoscale charge inhomogeneities are generically present in cuprates, but they are not necessarily completely static or periodic. (Yet the evidence for charge modulations that are both static and periodic was reported for a variety of cupratesHoffman et al. (2002); Vershinin et al. (2004); Hanaguri et al. (2004); Wise et al. (2008); Comin and Damascelli (2016); Peng et al. (2018).) It is further plausible that, in the general case, clusters with stronger antiferromagnetic correlations are not magnetically ordered with respect to each other. In fact, each of them might be in a singlet state as far as the total spin is concerned (see Ref.Lychkovskiy and Fine (2018)). As mentioned earlier, our model in the present form is lacking variables describing the degrees of freedom of the spin background. We expect that, generic situation in cuprates is such that the effective hopping between different nanoregions of inhomogeneous electronic background is larger than that for 1/8-doped lanthanum cuprates. This larger hopping can be sufficient to make $T_{ph}$, $T_{c}$ to approach each other, but still not large enough to dominate over the principal center-of-mass conserving interaction term appearing in Hamiltonian (1). Therefore, apart from suppressing the superconducting phase fluctuations, the model with the dynamic spin background can, quite plausibly, have the variational solution close in its principal features to the one obtained in the present work, in particular, as far as the gap of the fermionic excitations is concerned.

VII Comparison of model predictions with experimentally measured energy gaps

We now demonstrate that the energy gap obtained in Section IV.3 from the variational solution of the model based on Hamiltonian (1) exhibits temperature dependence quite similar to that of the superconducting energy gap experimentally observed in a variety of cuprates. In Figs. 5 and 6, we make comparisons with the available experimental results for break junctions (BJ)Vedeneev et al. (1994); Miyakawa et al. (1998); Akimenko et al. (1999) and the interlayer tunneling (ILT) Suzuki et al. (1999); Suzuki and Watanabe (2000); Krasnov et al. (2000); Krasnov (2002a, b); Yurgens et al. (2003); Yamada et al. (2003); Bae et al. (2008); Suzuki et al. (2012); Ren et al. (2012); Kambara et al. (2013); Krasnov (2016) for the bismuth family of cuprates. Figure  5 includes the data sets reviewed in Ref.Fine (2005), while Fig. 6 covers the experiments done later. The predictions of the grid-based model are also plotted in Figs. 5 and 6.

The model predictions, when limited to Cases I or II only, require two input parameters $\Delta(0)$ and $T_{c}$, which help us to determine $g$ and $|\epsilon_{a}|$ for Case I, or $|\epsilon_{b}|$ for Case II. The choice between Cases I and II is made on the basis of the ratio $\Delta(0)/T_{c}$ being larger or smaller than $2\sqrt{2}$.

All plots in Figs. 5 and 6 demonstrate either very good or satisfactory agreement between the predictions of the present model and the experiment. In comparison with the predictions of the grid-based model, the agreement with experiment has improved overall. It should be remarked here that the experimental data themselves are subject to a number of uncertainties, including, in particular, the overheating effect for the ILT measurements Suzuki et al. (1999); Suzuki and Watanabe (2000); Yamada et al. (2003); Krasnov et al. (2000); Krasnov (2002a, b); Yurgens et al. (2003); Zavaritsky (2004a, b); Yurgens et al. (2004); Krasnov et al. (2005); Kurter et al. (2010); Krasnov (2011); Kurter et al. (2011); Suzuki et al. (2012).

We further remark that, despite the significant experimental difficulty of measuring $\Delta(T)$ close to $T=T_{c}$, the very notion of the existence fo the critical ratio $\Delta(0)/T_{c}$ which separates the dependencies ending with $\Delta(T_{c})=0$ from those ending with $\Delta(T_{c})\neq 0$ appears to be reasonably supported by experiments, and, moreover, the value $2\sqrt{2}$ for such a critical ratio obtained in this work leads to more consistent predictions than the critical value 4 obtained in Ref.Fine (2005) for the grid-based model.

Finally, in Fig. 7, we present the empirical dependence of the model parameters $g$ and $\epsilon_{a}$ on the doping level. We extract this dependence from the ILT experiments reported in Ref. Suzuki et al. (2012), where the values of doping were explicitly indicated. The parameters $g$ and $\epsilon_{a}$ were obtained numerically by solving Eqs.(30) and (31) with $\epsilon_{b}=0$, $n_{B}=0$ and with experimentally determined input parameters $\Delta(0)$ and $T_{c}$. The resulting plot shows that, in the doping range between 0.09 and 0.2, $g$ depends on the doping level rather weakly, while $\epsilon_{a}$ depends strongly — it decreases nearly linearly with increasing doping and, if extrapolated, reaches zero around the doping level 0.25, i.e. close to the level, above which superconductivity disappears. Such a behavior of $\epsilon_{a}$ is consistent with the behavior of the pseudogap expected from other experiments. The eventual disappearance of superconductivity at higher dopings implies, that, if the present model is relevant, then the interaction constant $g$ should steeply decrease to zero for doping levels beyond the plotting range of Fig. 7 (i.e. above 0.2).

VIII Conclusions

In the present paper, we generalized the microscopic model proposed in Ref. Fine (2004) for the grid background to the background formed by the checkerboard of spin vortices. The technical difference is that the former involves two kinds of fermionic states, while the later involves three, even though two of the three are similar. We have shown that the predictions of the grid-based model for the temperature evolution of the energy gap largely remain intact. The most important difference between the spin-vortices-based model and the grid-based model turns out to be the critical ratio $\Delta(0)/T_{c}$ above which the temperature dependence of the energy gap ends at the value $\Delta(T_{c})\neq 0$, which, in turn, is, probably, related to the pseudogap. For spin vortices, this critical value is $2\sqrt{2}$, while, for grid, it is 4. We have further demonstrated that the predictions of the spin-vortices-based model for the temperature evolution of the energy gap exhibit good agreement with experiments, and, moreover, this agreement is somewhat better than that for the grid-based model.

An important difference of the present work from that of Ref.Fine (2004) is the treatment of the superfluid stiffness of the resulting solution. We have shown that nonlocal character of the interaction term of the kind used in Ref.Fine (2004) is not sufficient to induce non-zero superfluid stiffness. The latter appears only once arbitrarily small hopping terms are added to the model. The resulting behavior then implies large superconducting fluctuations of the kind that, possibly, exist in 1/8-doped LBCO.

In the broader context of cuprate superconductivity, the model considered in this work is still rather oversimplified. However, one can use it to develop intuition about more realistic settings that must involve the fluctuations of the spin background, as well as other interactions between fermions.

Acknowledgements.

The authors are grateful to J. Haase, M. Jurkutat, J. Kohlrautz and A. V. Rozhkov for discussions. This project was funded by the Skoltech NGP Program (Skoltech-MIT joint project).

Appendix A Total energy of the system

In this Appendix, we elaborate on the derivation steps (ii) and (iii) mentioned after Eq.(20).

Substituting the canonical transformation for $a$- and $c$-states (16-19) in (13), we obtain

[TABLE]

where all variables are defined in Section (III).

Two interaction terms in the above expression have phase-dependent factors $\cos[\phi_{V}(\mathbf{k})+\phi_{a}(\mathbf{k})]$ and $\cos[\phi_{\tilde{V}}(\mathbf{k})+\phi_{c}(\mathbf{k})]$. Eventually, the variational ground-state energy obtained by finding $u(\mathbf{k})$, $v(\mathbf{k})$, $p(\mathbf{k})$ and $q(\mathbf{k})$ will monotonically decrease with the increasing absolute value of these terms. This implies that the variational energy will be minimized for $|\cos[\phi_{V}(\mathbf{k})+\phi_{a}(\mathbf{k})]|=1$ and $|\cos[\phi_{\tilde{V}}(\mathbf{k})+\phi_{c}(\mathbf{k})]|=1$. Choosing the sign of cosines in these relations is just a matter of sign convention for the Bogoliubov transformation coefficients later converting into the sign of the products $swu(\mathbf{k})v(\mathbf{k})$ and $swp(\mathbf{k})q(\mathbf{k})$. We adopt the convention that the signs of the above products are negative in the state minimizing the total energy.

Appendix B Case I

For $\epsilon_{b}=0$, the Bogoliubov transformation parameters $s$ and $w$ enter the energy (24) only as a term proportional to $sw$. For such a case, given the constraint $s^{2}+w^{2}=1$, the minimization of energy (24) gives, $s=\sqrt{\frac{1}{2}},$ $w=-\sqrt{\frac{1}{2}}$. The relative negative sign of $s$ and $w$ implies later the positive relative sign for the pairs of transformation parameters $\{u(\mathbf{k}),v(\mathbf{k})\}$ and $\{p(\mathbf{k}),q(\mathbf{k})\}$. For $\epsilon_{a}\geq 0$, the minimization of energy (24) with respect to $u(\mathbf{k})$, $v(\mathbf{k})$, $p(\mathbf{k})$, and $q(\mathbf{k})$, finally, gives

[TABLE]

[TABLE]

where,

[TABLE]

[TABLE]

The total energy $E$ of the system in this case can be expressed as

[TABLE]

which is an implicit function of $|V(\mathbf{k})|$ and $|\tilde{V}(\mathbf{k})|.$ Both $|V(\mathbf{k})|$ and $|\tilde{V}(\mathbf{k})|$ are a function of four phases $\varphi_{1}$, $\varphi_{2}$, $\varphi_{3}$ and $\varphi_{4}$ entering Eqs. (III, III). Therefore the energy equation should be further minimized with respect to the values of these phases. Such minimization imposes only one constraintFine (2004)

[TABLE]

Appendix C Case II

In this case, $\epsilon_{a}=\epsilon_{c}=0$ in Eq. (24) and, as a result, the minimization of energy gives $u(\mathbf{k})=v(\mathbf{k})=p(\mathbf{k})=q(\mathbf{k})=1/\sqrt{2}$. The minimization with respect to $s$ subject to condition (3), now gives

[TABLE]

where we introduced the following parameters:

[TABLE]

The parameters $C_{a}$ and $C_{c}$ are defined by Eqs.(36,37). The solution of the biquadratic equation (83) that minimizes the total energy under condition $\epsilon_{b}\geq 0$ is

[TABLE]

We obtain $\epsilon_{B}$ by varying energy (24) with respect to $n_{B}$:

[TABLE]

The substitution of the parameters defined by Eqs.(85,84) now gives

[TABLE]

which, after some manipulations, leads to Eq.(32).

The quasiparticle excitation energies for $A$- and $C$- states are calculated in a similar way — as $\epsilon_{A}(\mathbf{k})\equiv\frac{1}{2}\frac{dE}{dn_{A}(\mathbf{k})}$ and $\epsilon_{C}(\mathbf{k})\equiv\frac{1}{2}\frac{dE}{dn_{C}(\mathbf{k})}$, which gives Eqs.(33) and (34) respectively. The total energy of the system thus becomes

[TABLE]

In this case, phases $\varphi_{\alpha}$ also obey the constraint (82).

Appendix D Effects of non-zero hopping on the variational solution

In this Appendix, we discuss the properties of the superconducting variational solution for Hamiltonian given by Eq. (68). We introduce Bogoliubov transformations in exactly the same way as was done in Section III. Then the energy of the system has the following form:

[TABLE]

where $\mathcal{E}(\mathbf{k})$ is the band energy given by equation (70). The only aspect that is different here in comparison with Eq.(24) is the $\mathbf{k}$-dependence of the band energies for the $a$- and $c$-states.

We are unable to find an analytic expression for the Bogoliubov transformation coefficients in the general case, but we still can do it in Case I ($\epsilon_{b}$=0), where the minimization of energy (87) gives $s=\sqrt{\frac{1}{2}},$ $w=-\sqrt{\frac{1}{2}}$. The relative negative sign of $s$ and $w$ implies later the positive relative sign for the pairs of transformation parameters $\{u(\mathbf{k}),v(\mathbf{k})\}$ and $\{p(\mathbf{k}),q(\mathbf{k})\}$. The minimization of energy (87) with respect to $u(\mathbf{k})$, $v(\mathbf{k})$, $p(\mathbf{k})$, and $q(\mathbf{k})$, leads to the result having the same form as the one given by equations (77)-(78), with the only difference being that the parameters $\mathcal{T}_{a}(\mathbf{k})$, $\mathcal{Q}_{a}(\mathbf{k})$, $\mathcal{T}_{c}(\mathbf{k})$ and $\mathcal{Q}_{c}(\mathbf{k})$ are now defined as:

[TABLE]

[TABLE]

The distinction between Eqs.(88,89) and Eqs.(79, 80) is that the former contain the band energy $\mathcal{E}(\mathbf{k})$ given by Eq.(70), while the latter contain the on-site energies $\epsilon_{a}$ and $\epsilon_{c}$. Equation (30) for the critical temperature is now replaced by the following one:

[TABLE]

Appendix E Additional constraints on phases $\varphi_{i}$ caused by hopping

Energy minimization for Hamiltonian (68) leads to additional constraints on phases $\varphi_{\alpha}$ appearing in Bogoliubov transformations (2). In this Appendix, we derive these constraints for the case $\epsilon_{b}=0$.

For convenience, let us introduce new phases:

[TABLE]

[TABLE]

[TABLE]

and then express $|V(\mathbf{k})|^{2}$ and $|\tilde{V}(\mathbf{k})|^{2}$ from Eqs. (III), (III) as

[TABLE]

and

[TABLE]

In the presence of non-zero hopping $t$, the expression (81) for the total energy of the system is modified as follows:

[TABLE]

where $\mathcal{E}(\mathbf{k})$ is the band energy given by Eq.(70).

Necessary conditions for phases $\alpha$, $\beta$ and $\gamma$ to minimize the total energy (96) are:

[TABLE]

In general, the system of equations (97)-(97) has more than one solution for a given set of parameters $\epsilon_{a}$, $g$ and $t$. The solution representing the global minimum should then be selected by computing the corresponding values of $E$.

We were able to find analytically three types of solutions of the system (97)-(99).

Solutions of type A:

[TABLE]

Solutions of type B:

[TABLE]

Solutions of type C (exist, only if $\displaystyle\frac{4|t\epsilon_{a}|}{g^{2}}\leqslant 1$):

[TABLE]

Here $k$, $m$ and $n$ are arbitrary integer numbers. We cannot prove analytically that the system (97)-(97) has no other solutions, but we checked numerically that the minimum of $E$ corresponds to the solutions of either type A, or B, or C.

Let us now prove that the sets of $\alpha$, $\beta$ and $\gamma$ given by Eqs.(100-103) are, indeed, the solutions of the system (97)-(99).

The substitution of Eq.(96) into Eqs.(97) and (98) gives, respectively,

[TABLE]

and

[TABLE]

We now observe that $\alpha=\pi n$, $\beta=\pi m$ and arbitrary $\gamma$ satisfy both Eq.(104) and Eq.(105), because these equations would then contain only summations of functions that are odd with respect to either $k_{x}$ or $k_{y}$, while the summation regions are symmetric. Thus all sets of $\alpha$, $\beta$ and $\gamma$ entering (100) -(103) satisfy Eqs.(104) and (105), and, hence, Eqs.(97) and (98)

The substitutions of Eq.(96) into Eq.(99) yields:

[TABLE]

The solutions of type A given by Eq.(100) imply that $\sin{\gamma}=0$. Hence they, obviously, satisfy Eq.(106).

Now let us show that the solutions of type B also satisfy Eq.(106). Once we substitute (101) to (106), we get an expression, proportional to

[TABLE]

This expression is zero, because it changes sign under the transformation $(k_{x},k_{y})\longrightarrow(k_{y},k_{x})$.

Finally, let us consider solutions of type C given by (102) and (103) and prove that they satisfy Eq.(106). In order to do it, let us examine the expressions under the square roots in the denominators in Eq.(106):

[TABLE]

and

[TABLE]

When we substitute (102) or (103) into (108)-(109), the terms 4t\epsilon_{a}\big{[}\cos{\left(2k_{x}L\right)}+\cos{\left(2k_{y}L\right)}\big{]} and g^{2}\cos{\gamma}\big{[}\cos{(2k_{x}L+\alpha)}+\cos{(2k_{y}L\pm\beta)}\big{]} cancel each other. For this reason, the sum over $\mathbf{k}$ in Eq.(106) for $\alpha$, $\beta$ and $\gamma$ given by (112) becomes proportional to

[TABLE]

Now we note that, after the simultaneous change of the signs of $\cos{(2k_{x}L)}$ and $\cos{(2k_{y}L)}$, the function under the sum also changes sign. A possible transformation achieving this is :

[TABLE]

It is illustrated in Fig. 8. The fact that transformation (111) maps the Brillouin zone onto itself, while the function summed in Eq.(110) changes sign, means that the sum itself is equal to zero and, hence, Eq.(106) is satisfied. This finishes the proof that Eqs.(100-103) represent the solutions of the system (97-99) and thus correspond to the extrema of $E$.

We have conducted extensive numerical tests of the minima of E in the space of parameters $\epsilon_{a}$, $g$ and $t$. In all these tests, the minimuma corresponded to the solutions of type A, B or C given by Eqs.(100-103). The numerically found “phase diagram” assigning the regions of the parameter space $(\epsilon_{a},t)$ to different types of solutions for $g=1$ is shown in Fig. 9.

In Section V.2 and in Appendix F, we focus on the limit $\displaystyle|t|\ll|g|\ll|\epsilon_{a}|$. In this case, the minima of $E$ are the following:

If $\displaystyle\frac{4|t\epsilon_{a}|}{g^{2}}\leqslant 1$, then

[TABLE]

If $\displaystyle\frac{4|t\epsilon_{a}|}{g^{2}}\geqslant 1$, then

[TABLE]

The transition between the two kinds of minima is illustrated in Fig. 10 through the evolution of $|\cos{\gamma}|$ as a function of $t$ for $\epsilon_{a}=10$ and $g=1$.

Appendix F Calculation of superfluid phase stiffness for Hamiltonian (68)

In this section, we obtain zero-temperature superfluid phase stiffness for Hamiltonian (68) with $\epsilon_{b}=0$.

Let us assume for simplicity that linearly changing in space phase $\theta$ has gradient only along the $x$-direction. We also introduce variable $\theta_{0}$ to denote the change of phase $\theta$ across the $2L\times 2L$ unit cell of the spin-vortex checkerboard. The gradient of $\theta$ can now be expressed as

[TABLE]

and, as a result, the energy (72) of the system associated with the phase gradient can be written as

[TABLE]

Using Eq.(72), we then obtain

[TABLE]

We now observe that $\langle\tilde{\Psi}|H^{\prime}|\tilde{\Psi}\rangle=\langle\Psi|\tilde{H}^{\prime}|\Psi\rangle$, where $\tilde{H}^{\prime}$ is the Hamiltonian obtained by replacing operators $a$, $b$ and $c$ in Hamiltonian (68) with operators $\tilde{a}$, $\tilde{b}$ and $\tilde{c}$ defined by Eqs.(64-66). This allows us to rewrite Eq.(116) as

[TABLE]

where

[TABLE]

For the mean-field ground state in the presence of hopping elements, one can use Eqs.(77) and (78) with ingredients from Eqs.(88) and (89) to obtain

[TABLE]

and

[TABLE]

After substituting Eqs.(119) and (120) into Eq.(117), we, finally, get

[TABLE]

(Here and below, we do not include terms that give zero after the summation.)

Now let us study expression (121) in the limit $|t|\ll|g|\ll|\epsilon_{a}|$ and find its dependence on $t$. Expanding it in the powers of $g/\epsilon_{a}$ and $t/g$, we get, in the leading order,

[TABLE]

For $\alpha$, $\beta$, $\gamma$ corresponding to the energy minima (112) and (113), we obtain from Eqs.(94) and (95) that,

if $\displaystyle\frac{4|t\epsilon_{a}|}{g^{2}}\leqslant 1$:

[TABLE]

and,

if $\displaystyle\frac{4|t\epsilon_{a}|}{g^{2}}\geqslant 1$:

[TABLE]

This leads to the following result:

if $\displaystyle\frac{4|t\epsilon_{a}|}{g^{2}}\leqslant 1$:

[TABLE]

if $\displaystyle\frac{4|t\epsilon_{a}|}{g^{2}}\geqslant 1$:

[TABLE]