Superfluid stiffness in cuprates: Effect of Mott transition and phase competition

TL;DR

This study investigates how Mott physics and phase competition influence the superfluid stiffness in cuprates, revealing dominant effects of Mott physics on hole-doped and antiferromagnetic competition on electron-doped sides.

Contribution

It provides a detailed computational analysis of superfluid stiffness in cuprates considering Mott transition and phase coexistence using Cellular Dynamical Mean-Field Theory.

Findings

Mott physics dominates superfluid stiffness on hole-doped cuprates.

Antiferromagnetism suppresses superfluid stiffness near half-filling on electron-doped cuprates.

Superfluid stiffness approaches BCS behavior at high doping levels.

Abstract

Superfluid stiffness $\rho_s$ is a defining characteristic of the superconducting state, allowing phase coherence and supercurrent. It is accessible experimentally through the penetration depth. Coexistence of $d$-wave superconductivity with other phases in underdoped cuprates, such as antiferromagnetism (AF) or charge-density waves (CDW), may drastically alter $\rho_s$. To shed light on this physics, the zero-temperature value of $\rho_s=\rho_{zz}$ along the $c$-axis was computed for different values of Hubbard interaction $U$ and different sets of tight-binding parameters describing the high-temperature superconductors YBCO and NCCO. We used Cellular Dynamical Mean-Field Theory for the one-band Hubbard model with exact diagonalization as impurity solver and state-of-the-art bath parametrization. We conclude that Mott physics plays a dominant role in determining the superfluid stiffness on the hole-doped side of the phase diagram. On the electron-doped side, antiferromagnetism wins over superconductivity near half-filling. But upon approaching optimal electron-doping, homogeneous coexistence between superconductivity and antiferromagnetism causes the superfluid stiffness to drop sharply. Hence, on the electron-doped side, it is competition between antiferromagnetism and $d$-wave superconductivity that plays a dominant role in determining the value of $\rho_{zz}$ near half-filling. At large overdoping, $\rho_{zz}$ behaves in a more BCS-like manner in both the electron- and hole-doped cases. We comment on some qualitative implications of these results for the superconducting transition temperature.

Tables

Table 1. Table 1: Tight-binding band parameters

compounds/parameters	$t^{\prime }/t$	$t^{\prime \prime }/t$
YBCO/BSCCO	-0.3	0.2
LSCO/NCCO	-0.17	0.03

Equations

$$\hat{\mathcal{H}}=\sum _{ij,\sigma }{t}_{ij}\left({\hat{c}}_{i,\sigma }^{†}{\hat{c}}_{j,\sigma }+{\hat{c}}_{j,\sigma }^{†}{\hat{c}}_{i,\sigma }\right)+U\sum _i{\hat{n}}_{i,\uparrow }{\hat{n}}_{i,\downarrow },$$

$${\hat{\mathcal{H}}}^{\prime }=$$

$$+\sum _{i\alpha ,\sigma }{\theta }_{i\alpha ,\sigma }\left({\hat{c}}_{i,\sigma }^{†}{\hat{a}}_{\alpha }+\text{H.c.}\right)+\sum _{\alpha ,\sigma }{ϵ}_{\alpha ,\sigma }{\hat{a}}_{\alpha ,\sigma }^{†}{\hat{a}}_{\alpha ,\sigma },$$

$${\mathcal{G}}_{{\mathbf{R}}_i{\mathbf{R}}_j}^{c,\sigma }(i{\omega }_n)={\left[(i{\omega }_n+\mu )\mathbf{I}-{\mathbf{t}}^{\prime }-{\mathbf{\Gamma }}^{\sigma }(i{\omega }_n)-{\mathbf{\Sigma }}_{c,\sigma }(i{\omega }_n)\right]}_{ij}^{-1},$$

$${\mathbf{\Gamma }}_{{\mathbf{R}}_i{\mathbf{R}}_j}^{\sigma }\left(i{\omega }_n\right)=\sum _{\alpha }\frac{{\theta }_{i\alpha ,\sigma }{\theta }_{j\alpha ,\sigma }^{\ast }}{i{\omega }_n-{ϵ}_{\alpha ,\sigma }}.$$

$${\bar{\mathcal{G}}}^{\sigma }(i{\omega }_n)=\frac{N_c}{N}\sum _{\tilde{\mathbf{k}}}\frac{1}{i{\omega }_n+\mu -\mathbf{t}(\tilde{\mathbf{k}})-{\mathbf{\Sigma }}_{c,\sigma }(i{\omega }_n)}$$

$$d=\sum _{\begin{array}{c}{\mathbf{R}}_i{\mathbf{R}}_j\\ i{\omega }_n\le i{\omega }_c\end{array}}\sum _{\sigma }W(i{\omega }_n){\left|{\left({\mathcal{G}}^{c,\sigma }(i{\omega }_n{)}^{-1}-{\bar{\mathcal{G}}}^{\sigma }(i{\omega }_n{)}^{-1}\right)}_{\begin{array}{c}{\mathbf{R}}_i{\mathbf{R}}_j\end{array}}\right|}^2.$$

$$\mathcal{G}(\tilde{\mathbf{k}},i{\omega }_n{)}^{-1}=(i{\omega }_n+\mu )\mathbf{I}-\mathbf{t}(\tilde{\mathbf{k}})-{\mathbf{\Sigma }}_c(i{\omega }_n),$$

$${\mathbf{\Sigma }}_c(i{\omega }_n)=(i{\omega }_n+\mu )\mathbf{I}-{\mathbf{t}}^{\prime }-{\mathcal{G}}^c(i{\omega }_n{)}^{-1}-\mathbf{\Gamma }(i{\omega }_n).$$

$$\begin{array}{rl}& {\hat{\Psi }}_{\tilde{\mathbf{k}}}=\left(\begin{array}{c}{\hat{c}}_{\tilde{\mathbf{k}}\uparrow ,1}\\ {\hat{c}}_{\tilde{\mathbf{k}}\uparrow ,2}\\ ⋮\\ {\hat{c}}_{-\tilde{\mathbf{k}}\downarrow ,N_c-1}^{†}\\ {\hat{c}}_{-\tilde{\mathbf{k}}\downarrow ,N_c}^{†}\end{array}\right)\\ & {\hat{\Psi }}_{\tilde{\mathbf{k}}}^{†}=\left(\begin{array}{ccccc}{\hat{c}}_{\tilde{\mathbf{k}}\uparrow ,1}^{†}& {\hat{c}}_{\tilde{\mathbf{k}}\uparrow ,2}^{†}& \dots & {\hat{c}}_{-\tilde{\mathbf{k}}\downarrow ,N_c-1}& {\hat{c}}_{-\tilde{\mathbf{k}}\downarrow ,N_c}\end{array}\right).\end{array}$$

$$\mathcal{G}(\tilde{\mathbf{k}},\tau )=-⟨{\hat{\mathcal{T}}}_{\tau }\hat{\Psi }(\tau ){\hat{\Psi }}^{†}(0){⟩}_{\hat{\mathcal{H}}}.$$

$$\mathcal{G}(\mathbf{k},i{\omega }_n)=\frac{1}{N_c}\sum _{{\mathbf{R}}_i,{\mathbf{R}}_{\mathbf{j}}}{e}^{-i\mathbf{k}\cdot \left({\mathbf{R}}_i-{\mathbf{R}}_j\right)}{\mathcal{G}}_{{\mathbf{R}}_i{\mathbf{R}}_j}\left(\tilde{\mathbf{k}},i{\omega }_n\right),$$

$${\rho }_{ab}$$

$$\displaystyle=\frac{-1}{V}\int_{0}^{\beta}d\tau\int d(\bm{r}-\bm{r}^{\prime})\frac{\delta^{2}\mathcal{F}[\bm{\mathcal{G}}]}{\delta A_{a}(\bm{r},\tau)\delta A_{b}(\bm{r}^{\prime},0)}\bigg{\rvert}_{\bm{A}=\bm{0}},$$

$$⟨{\hat{J}}_a(\mathbf{q}=0,\omega =0)⟩={\rho }_{ab}{A}_b(\mathbf{q}=0,\omega =0).$$

$${\lambda }_{ab}^{-2}={\rho }_{ab}{\mu }_0.$$

$${\lambda }_{ab}^{-2}=\frac{n_s{e}^2}{m^{\ast }}{\mu }_0,$$

$${\rho }_{ab}$$

$$\displaystyle+\text{tr}\left[\bm{\mathcal{G}}(\bar{k})\lambda^{ab}_{\bar{\bm{k}}}\right]\biggr{)}\bigg{\rvert}_{\bm{A}=\bm{0}},$$

$${\lambda }_{\bar{\mathbf{k}}}^i$$

$${\lambda }_{\bar{\mathbf{k}}}^{ji}$$

$$t_{\perp }^2(\mathbf{k})=t_{\text{bi}}^2{\cos }^2{k}_z{\left(\cos k_x-\cos k_y\right)}^4,$$

$${\rho }_{zz}^{SC}=\frac{e^2}{{\hslash }^2\beta VN}$$

$$\displaystyle\bigg{(}\text{tr}\left[\bm{\mathcal{G}}(k)\bm{\mathcal{G}}(k)\right]-\text{tr}\left[\sigma_{3}\bm{\mathcal{G}}(k)\sigma_{3}\bm{\mathcal{G}}(k)\right]\bigg{)},$$

$${\bar{t}}_{\perp }^2$$

$$=\frac{t_{\text{bi}}^2}{2}(\cos k_x-\cos k_y{)}^4,$$

$$T_{lm}={\sigma }_{\alpha \beta }^l{\tau }_{ab}^m,$$

$${\rho }_{zz}^{AF+SC}$$

$$\displaystyle\bigg{(}\text{tr}\left[\bm{\mathcal{G}}T_{01}\bm{\mathcal{G}}T_{01}\right]-\text{tr}\left[\bm{\mathcal{G}}T_{31}\bm{\mathcal{G}}T_{31}\right]\bigg{)}.$$

$$\mathcal{G}(\tilde{k})=\frac{1}{i{\omega }_n+\mu -{\mathcal{H}}_{\tilde{\mathbf{k}}}^0-{\Sigma }_c(i{\omega }_n)},$$

$${\hat{\Psi }}_{\tilde{\mathbf{k}}}^{†}=\left(\begin{array}{ccccc}{\hat{c}}_{\tilde{\mathbf{k}},\uparrow ,A,1}^{†}& {\hat{c}}_{\tilde{\mathbf{k}},\uparrow ,A,2}^{†}& \dots & {\hat{c}}_{-\tilde{\mathbf{k}},\downarrow ,B,1}& {\hat{c}}_{-\tilde{\mathbf{k}},\downarrow ,B,2}\end{array}\right).$$

$${\mathcal{H}}_{MF}^{\alpha \beta ;ab;rs}(\tilde{\mathbf{k}})=$$

$$\text{medsize}\left(\begin{array}{cccccccc}{\Omega }_{\tilde{\mathbf{k}}}-M& {\zeta }_{\tilde{\mathbf{k}}}& {ϵ}_{\tilde{\mathbf{k}}}& {ϵ}_{\tilde{\mathbf{k}}}& {\Delta }_s& {\Delta }_{p,\tilde{\mathbf{k}}}& {\Delta }_{d,\tilde{\mathbf{k}}}^x& {\Delta }_{d,\tilde{\mathbf{k}}}^y\\ {\zeta }_{\tilde{\mathbf{k}}}& {\Omega }_{\tilde{\mathbf{k}}}-M& {ϵ}_{\tilde{\mathbf{k}}}& {ϵ}_{\tilde{\mathbf{k}}}& {\Delta }_{p,\tilde{\mathbf{k}}}& {\Delta }_s& {\Delta }_{d,\tilde{\mathbf{k}}}^y& {\Delta }_{d,\tilde{\mathbf{k}}}^x\\ {ϵ}_{\tilde{\mathbf{k}}}& {ϵ}_{\tilde{\mathbf{k}}}& {\Omega }_{\tilde{\mathbf{k}}}+M& {\zeta }_{\tilde{\mathbf{k}}}& {\Delta }_{d,\tilde{\mathbf{k}}}^x& {\Delta }_{d,\tilde{\mathbf{k}}}^y& {\Delta }_s& {\Delta }_{p,\tilde{\mathbf{k}}}\\ {ϵ}_{\tilde{\mathbf{k}}}& {ϵ}_{\tilde{\mathbf{k}}}& {\zeta }_{\tilde{\mathbf{k}}}& {\Omega }_{\tilde{\mathbf{k}}}+M& {\Delta }_{d,\tilde{\mathbf{k}}}^y& {\Delta }_{d,\tilde{\mathbf{k}}}^x& {\Delta }_{p,\tilde{\mathbf{k}}}& {\Delta }_s\\ {\Delta }_s^{\ast }& {\Delta }_{p,\tilde{\mathbf{k}}}^{\ast }& {{\Delta }^x}_{d,\tilde{\mathbf{k}}}^{\ast }& {{\Delta }^y}_{d,\tilde{\mathbf{k}}}^{\ast }& -{\Omega }_{-\tilde{\mathbf{k}}}+M& -{\zeta }_{-\tilde{\mathbf{k}}}& -{ϵ}_{-\tilde{\mathbf{k}}}& -{ϵ}_{-\tilde{\mathbf{k}}}\\ {\Delta }_{p,\tilde{\mathbf{k}}}^{\ast }& {\Delta }_s^{\ast }& {{\Delta }^y}_{d,\tilde{\mathbf{k}}}^{\ast }& {{\Delta }^x}_{d,\tilde{\mathbf{k}}}^{\ast }& -{\zeta }_{-\tilde{\mathbf{k}}}& -{\Omega }_{-\tilde{\mathbf{k}}}+M& -{ϵ}_{-\tilde{\mathbf{k}}}& -{ϵ}_{-\tilde{\mathbf{k}}}\\ {{\Delta }^x}_{d,\tilde{\mathbf{k}}}^{\ast }& {{\Delta }^y}_{d,\tilde{\mathbf{k}}}^{\ast }& {\Delta }_s^{\ast }& {\Delta }_{p,\tilde{\mathbf{k}}}^{\ast }& -{ϵ}_{-\tilde{\mathbf{k}}}& -{ϵ}_{-\tilde{\mathbf{k}}}& -{\Omega }_{-\tilde{\mathbf{k}}}-M& -{\zeta }_{-\tilde{\mathbf{k}}}\\ {{\Delta }^y}_{d,\tilde{\mathbf{k}}}^{\ast }& {{\Delta }^x}_{d,\tilde{\mathbf{k}}}^{\ast }& {\Delta }_{p,\tilde{\mathbf{k}}}^{\ast }& {\Delta }_s^{\ast }& -{ϵ}_{-\tilde{\mathbf{k}}}& -{ϵ}_{-\tilde{\mathbf{k}}}& -{\zeta }_{-\tilde{\mathbf{k}}}& -{\Omega }_{-\tilde{\mathbf{k}}}-M\end{array}\right).$$

$${\hat{\mathcal{H}}}_{MF}$$

Full text

I Introduction

Long-range order leads to emergent phenomena, such as generalized rigidities Anderson (2018). The generalized rigidity associated with the superconducting state is superfluid stiffness, which allows phase coherence and supercurrent. Superfluid stiffness is accessible through penetration depth measurements Uemura et al. (1991). In cuprates, one expects that strong interactions will modify the BCS predictions. Mott physics should make $\rho_{s}$ smaller than the BCS value as one approaches half-filling, but what about the effect of a competing order, such as antiferromagnetism (AF) (Mandal et al., 2018; Motoyama et al., 2007; Armitage et al., 2010; Saadaoui et al., 2015) or charge-density waves (Ghiringhelli et al., 2012; Cyr-Choinière et al., 2018)? This is particularly important for the electron-doped cuprates where long-range AF order has been observed far from half-filling, competing with $d$-wave superconductivity ($d$SC). This competition could explain the fall of both $H_{c_{2}}$ and $T_{c}$ in the underdoped cuprates Mandal et al. (2018); Saadaoui et al. (2015); Keimer et al. (2015).

Very few theoretical works have addressed the question of the effect of microscopic homogeneous coexistence between AF and $d$SC on $\rho_{s}$. These works, based on mean-field calculations, have come to the conclusion that microscopic coexistence should decrease $\rho_{s}$ (Tobijaszewska and Micnas, 2005; Sharapov and Carbotte, 2006; Atkinson, 2007; Xiang and Wheatley, 1996; Mallik et al., 2018). Similar conclusions are reached with mean-field equations that use effective interactions generated by the functional renormalization group (Metzner and Yamase, 2019). But all these theoretical works discard the effect of the strong electron-electron interaction and of the Mott transition, while it is known that the cuprates are doped Mott insulators Anderson (1987).

The best way to take Mott physics into account in two dimensions is to use cluster generalizations of dynamical mean-field theory (Maier et al., 2005; Kotliar et al., 2006; Tremblay et al., 2006) for the Hubbard model. The only calculation of superfluid stiffness using these methods was done in the uniform superconducting state (Gull and Millis, 2013), not in a phase where superconductivity coexists microscopically with antiferromagnetism. By microscopic coexistence, which we are interested in, we mean that both order parameters are present simultaneously and homogeneously in the ground state. By contrast, macroscopic coexistence would refer to what happens at a first-order transition where phases coexist in separate macroscopic regions.

In this paper, we address the following questions: (1) What is the effect of the Mott transition on $\rho_{s}(n)$ near half-filling ($n=1$)? (2) Is microscopic coexistence with antiferromagnetism in the underdoped regime even more detrimental to $\rho_{s}$ than the Mott transition? (3) Is there a range of filling where BCS behavior is recovered? (4) What is the effect on $\rho_{zz}$ of the in-plane modulation of interplane hopping? To answer these questions, we compute the $c$-axis superfluid stiffness $\rho_{zz}$ for the one-band two-dimensional Hubbard model with band parameters appropriate to hole- and electron-doped cuprates. Along that direction, vertex corrections can be neglected, as we will discuss. We solve the Hubbard model using cellular dynamical mean-field theory (CDMFT) on a $2\times 2$ plaquette using an exact-diagonalization solver. The sites represent the Cu 3$d_{x^{2}-y^{2}}$ orbitals within the CuO2 planes of the cuprates. We relax symmetries to allow both AF, $d$SC and their coexistence. We call “coexistence” the regime of calculations that allows for both the AF and $d$SC order parameters to coexist microscopically111Thereof, in the CMDFT calculations, the pure $d$SC phase can be emerging from pure or coexistence regimes.. By contrast in the “pure” regime, the only symmetry breaking is $d$SC. The Hubbard model and the method to solve it is presented in section II. In section III, we follow up with the presentation of the formulae for $\rho_{zz}$ in both the pure $d$SC and coexisting AF+$d$SC states. We show the results in section IV. The discussion is in section V. We conclude in section VI. The supplemental material 222See Supplemental Material at [URL will be inserted by publisher] for estimates of the in-plane superfluid stiffness and comments on the implications for $T_{c}$. Additional referenceSchrieffer (2018). contains results for the in-plane superfluid stiffness, neglecting vertex corrections. That allows us to comment briefly in section V on the expected qualitative consequences of our results on the value of the superconducting transition temperature. This work is based on Ref. Simard, 2019 where further details may be found.

II Model and method

The following subsections present in turn the model, the method and the periodization procedure.

II.1 Hubbard model

To simulate interactions affecting electrons in high-Tc cuprates, it was suggested by Anderson (Anderson, 1987) that the Hubbard model

[TABLE]

would encompass key aspects of these strongly correlated materials. Here, $t_{ij}$ are hopping amplitudes, $\sigma\in\{\uparrow,\downarrow\}$ are spin indices, $\hat{c}^{(\dagger)}_{i,\sigma}$ are annihilation (creation) operators in localized Wannier states labeled by $i,j$, while $\hat{n}_{i\sigma}=\hat{c}^{\dagger}_{i,\sigma}\hat{c}_{i,\sigma}$ is the number operator, and $U$ is the local repulsion normalized by the first-neighbor hopping term $t$. The Hubbard model for CuO2 planes of cuprates Jorgensen et al. (1990) is on a square lattice with spacing $a$. We take $c$ for the lattice spacing in the perpendicular $z$ direction. We set $\hbar$, $k_{B}$, electric charge $e$ and lattice spacings $a$, $c$ equal to unity for the figures. Physical units are restored for a few estimates and for some formulae. We used first-, second- and third-neighbor hopping terms to simulate bare electronic dispersion relations. We denote these hoppings respectively by $t$, $t^{\prime}$ and $t^{\prime\prime}$ (see Figs. 5 and 6). The tight-binding band parameters used are displayed in Table 1 (Pavarini et al., 2001; Kyung et al., 2004). YBa2Cu3O7 (YBCO), La2-xSrxCuO4 (LSCO) and Bi2Sr2CaCu2O8+x (BSCCO) are hole-doped compounds while Nd2-xCexCuO4 (NCCO) is electron-doped. Nevertheless, to highlight the physics we consider the whole range of dopings for all sets of parameters.

In this work, we used the Green’s functions obtained in Ref. Foley et al., 2019 using CDMFT with the best available bath parametrization method, as described in the following subection. We also used the definitions given in Ref. Foley et al., 2019 for the superconducting and antiferromagnetic order parameters, respectively denoted $\langle D\rangle$ and $\langle M\rangle$.

II.2 ED-CDMFT

In CDMFT Biroli and Kotliar (2002), a cluster of size $2\times 2$ representing a finite portion of the full lattice is hybridized to a bath of non-interacting electrons to simulate the effect of the environment on the cluster’s electron Green’s function. Hence, the number of orbitals with interactions is $N_{c}=8$ (counting spin degeneracy). The cluster Hamiltonian $\hat{\mathcal{H}}^{\prime}$ including the hybridization to the baths reads Sénéchal (2008); Caffarel and Krauth (1994)

[TABLE]

where $\hat{c}^{(\dagger)}$ annihilates (creates) an electron on the cluster and $\hat{a}^{(\dagger)}$ annihilates (creates) an electron in the bath. The intra-cluster hopping matrix is $t_{ij}$ with $i$ and $j$ labelling the cluster sites and $\sigma\in\{\uparrow,\downarrow\}$. The baths are coupled to the cluster via the hybridization matrix $\theta_{i\alpha,\sigma}$ with $\alpha$ labelling the bath-orbital energy: the $\theta_{i\alpha,\sigma}$ represents the hopping of electrons between the cluster sites and the bath sites while $\epsilon_{\alpha,\sigma}$ is the energy of each orbital. The cluster Green’s function is computed with an ED (impurity) solver based on the Lanczos algorithm Sénéchal (2008); Caffarel and Krauth (1994).

In quantum cluster methods, the position is written as $\bm{r}=\bm{\tilde{r}}+\bm{R}$ where $\bm{\tilde{r}}$ is the base position of the cluster and $\bm{R}$ the position within the cluster. Likewise, a wave vector $\bm{k}$ in the Brillouin zone is decomposed as $\bm{k}=\bm{\tilde{k}}+\bm{K}$ where $\bm{\tilde{k}}$ belongs to the Brillouin zone of the superlattice of clusters (or reduced Brillouin zone) and $\bm{K}$ (which can also be seen as labeling the irreducible representations of the symmetry group of the cluster) belongs to the reciprocal superlattice.

We work on the imaginary axis and the fermionic Matsubara frequencies are $\omega_{n}=2\pi(n+1)/\beta$ where $n\in\mathbb{Z}$ and $\beta$ is the inverse temperature. The fictitious temperature defining the Matsubara grid is $\beta=50/t$. The interacting cluster Green’s function $\bm{\mathcal{G}}^{c,\sigma}_{\bm{R}_{i}\bm{R}_{j}}(i\omega_{n})$ in the cluster-site mixed basis $(\bm{\tilde{k}},\bm{R})$ breaks down as follows

[TABLE]

where $\bm{\Sigma}_{c,\sigma}$ is the cluster self-energy matrix, $\bm{t^{\prime}}$ the intra-cluster hopping matrix, $\mu$ the chemical potential, and $\bm{\Gamma}^{\sigma}$ the hybridization function whose expression can be deduced from Eq. (II.2):

[TABLE]

Each bath site is chosen to be in one of the irreducible representations of the cluster. That determines the symmetries of the $\theta$’s (Foley et al., 2019). In the following, on some occasions, the cluster-site indices and the spin will be left implicit. The $\bm{\Sigma}_{c,\sigma}$ used in our calculations is the one that satisfies the convergence criterion for the hybridization function. More specifically the cluster-projected Green’s function in the cluster-site mixed basis $(\bm{\tilde{k}},\bm{R})$

[TABLE]

and the cluster Green’s function Eq. (3) should be equal within a tolerance that sets the upper bound of the distance function $d$ which we minimize (Foley et al., 2019):

[TABLE]

When a finite number of bath orbitals is used to represent the environment, one can’t expect to obtain $d=0$. Therefore, to capture the important degrees of freedom, one introduces a frequency cutoff $i\omega_{c}$, with $\omega_{c}=2t$, to focus on the low-energy scale. The weight function $W(i\omega_{n})$ is such that $W(i\omega_{n})=1$ if $\omega_{n}<2t$ and $W(i\omega_{n})=0$ otherwise. Further details about the implementation can be found in Ref. Foley et al., 2019. In the equation for the lattice Green’s function Eq. (5), $\bm{t}(\bm{\tilde{k}})=\bm{t^{\prime}}+\delta\bm{t}(\bm{\tilde{k}})$ represents the complete lattice hopping matrix, with $\delta\bm{t}(\bm{\tilde{k}})$ the intercluster hopping amplitude matrix carrying a phase proportional to both $\bm{\tilde{k}}$ and the lattice parameters. $N$ stands for the total number of sites on the full lattice. Once $d$ has been minimized, the full lattice Green’s function $\bm{\mathcal{G}}(\bm{\tilde{k}},i\omega_{n})$, dropping spin indices, reads

[TABLE]

where at each iteration the lattice self-energy is the same as that of the cluster $\bm{\Sigma}_{c}$:

[TABLE]

To account for superconductivity, the lattice Green’s function Eq. (7) is expressed in the following Nambu basis, assuming singlet pairing:

[TABLE]

In imaginary time, the definition is

[TABLE]

The above formulae for ED-CDMFT must be expressed in Nambu space, taking into account that they are no longer diagonal in Nambu indices.

To avoid difficulties associated with the discreteness of the spectrum in ED, a ficticious temperature $\beta$ is introduced to compute $\rho_{zz}$. Since $\rho_{zz}$ converges rapidly with increasing $\beta$, this can be done with minimal effect on the accuracy of the zero-temperature calculation. All the results shown in section IV were computed using $500$ Matsubara frequencies and $\beta=\frac{500}{t}$. In Ref. Simard, 2019, it is shown explicitly that $\rho_{s}$ converges fast with respect to the number of Matsubara frequencies used in the summation ($\rho_{s}\propto\frac{1}{(i\omega_{n})^{4}}$) and with respect to the fictitious temperature.

II.3 Periodization

Once the lattice Green’s function $\bm{\mathcal{G}}(\bm{\tilde{k}},i\omega_{n})$ has been computed, one can periodize the latter to define it over the original Brillouin zone and recover translational invariance. For example, in a AF+$d$SC coexistence phase, periodizing $\bm{\mathcal{G}}(\bm{\tilde{k}},i\omega_{n})$ to extend it over the reduced AF Brillouin zone (AF-BZ) seems natural (see Fig. 1). Doing so, the initially $8\times 8$ cluster Green’s function in the mixed basis shrinks to $4\times 4$. The periodized cluster Green’s function is Sénéchal et al. (2000); Sakai et al. (2012); Sénéchal (2015)

[TABLE]

where $N_{c}$ accounts for the number of cluster sites and $\bm{k}=\bm{\tilde{k}}+\bm{K}$. For periodization in the AF+$d$SC phase, $N_{c}=2$ and $\bm{K}_{i}\in\{(0,0),(\pi,0)\}\ \text{or}\ \{(0,0),(0,\pi)\}$, while in the SC state, $N_{c}=4$ and the reciprocal-superlattice wavevectors are $\bm{K}_{i}\in\{(0,0),(\pi,0),(0,\pi),(\pi,\pi)\}$. In the procedure with coexistence, the two sets of $\bm{K}$ values lead to exactly the same periodized Green’s function, as can be understood with the aid of Fig. 1. Periodizing the Green’s function Eq. (7) using Eq. (11) reduces its dimensionality: for the case where AF and $d$SC coexist, the cluster Green’s function in the reduced AF Brillouin zone shown in Fig. 1 suffices to compute the superfluid stiffness. Eq. (11) is not a unitary transformation, because a unitary transformation would involve off-diagonal reciprocal-superlattice wavevectors, but one would not recover translational invariance.

From now on, we use the four-vector notation $k\equiv(\bm{k},i\omega_{n})$ to lighten the notation, namely $\mathcal{G}(\bm{k},i\omega_{n})\to\mathcal{G}(k)$. Note that cumulant periodization (Stanescu and Kotliar, 2006) gives unphysical results for the superfluid stiffness Simard (2019), especially for YBCO-like tight-binding calculations. Self-energy periodization Biroli and Kotliar (2002) leads to unphysical states in the Mott gap at half-filling Sénéchal (2015), hence we do not consider these here.

III Superfluid Stiffness

In this section we explain the general formalism for computing the superfluid stiffness on the superlattice of clusters and explain why vertex corrections can be neglected (subsection III.1). Then we give the formula for pure $d$-wave superconductivity (subsection III.2), and finally the formula for the regime of coexistence between antiferromagnetism and superconductivity (subsection III.3).

III.1 General formula

The superconducting order parameter is a consequence of spontaneous $U(1)$ symmetry breaking, the global phase being arbitrarily fixed. The phase rigidity, or superfluid stiffness, of the superconducting ground state accounts for the change in free energy when twisting the phase of the order parameter. In the linear response framework, the superfluid stiffness ${\rho}_{ab}$ for the response to a transverse vector potential $\bm{A}$ is related to the current-current correlation function by

[TABLE]

where subscripts $a,b\in\{x,y,z\}$ denote the cartesian axes, $V$ is the volume of a unit cell and $\mathcal{F}$ is the free energy (or energy at $T=0$) of the system. In particular, we evaluate the current $\hat{J}_{z}(\bm{r},\tau)$ along the $c$-axis induced by a magnetic field applied in the transverse direction (in the plane). As we discuss below, this allows a calculation where neglecting vertex corrections can be justified. The vector potential $A_{z}(\bm{r}^{\prime},\tau^{\prime})$ representing the magnetic field is chosen along the $c$-axis as well. The above fomula Eq. (III.1) assumes that we are in the London limit where the kernel of the electromagnetic response can be evaluated in the zero wave vector ($\bm{q=0}$) limit so that linear response theory gives us

[TABLE]

Taking the curl of this equation and using Ampère’s law $\nabla\times\bm{B}=\mu_{0}\bm{J}$ with $\mu_{0}$ the permeability of the vacuum, one finds that the London penetration depth $\lambda$ is related to the superfluid stiffness by

[TABLE]

In the BCS or Ginzburg-Landau formalism, this is written in terms of the superfluid density $n_{s}$

[TABLE]

where $e$ is the electric charge and $m^{\ast}$ the effective mass of the electrons.

On the lattice, coupling of the Bloch electrons to the electromagnetic field is done via the Peierls substitution in the orbital basis ($\tilde{\bm{k}},\bm{R}$) (mixed basis). Since we can work in the $\bm{q=0}$ limit, the vector potential is a constant and the Peierls substitution leads to the replacement $\partial_{A_{i}}\rightarrow-\frac{e}{\hbar}\partial_{\bar{k}_{i}}$, where $\bar{k}_{i}\equiv k_{i}-\frac{e}{\hbar}A_{i}$, as long at the phase difference between atoms in the same unit cell is taken into account in the Fourier transforms (Nourafkan and Tremblay, 2018). Otherwise, the expression for the currents is different (Tomczak and Biermann, 2009). This is discussed further in Appendix A.

When vertex corrections are neglected, the superfluid stiffness is given by

[TABLE]

where the trace $\text{tr}\left[\ldots\right]$ acts in the cluster-site mixed basis and $N$ is the number of unit cells. The first and second terms of Eq. (III.1) are, respectively, the paramagnetic and diamagnetic contributions to $\rho_{ab}$. Since the Nambu formalism involves a particle-hole transformation for the down electrons, we must evaluate the derivative with respect to the vector potential as follows $T_{0(m\times m)}\partial_{A_{i}}=-\frac{e}{\hbar}T_{3(m\times m)}\partial_{\bar{k}_{i}}$, where the tensors are defined by $T_{0(m\times m)}\equiv\sigma^{0}\otimes\mathbb{I}_{m\times m}$ and $T_{3(m\times m)}\equiv\sigma^{3}\otimes\mathbb{I}_{m\times m}$ with $\sigma^{0}$ the $2\times 2$ identity matrix and $\sigma^{3}$ the diagonal Pauli matrix whose components indicate whether we are in the spin up or spin-down part of the Nambu spinor Eq. (9), a minus sign needing to be included in the spin-down part. The identity matrix $\mathbb{I}_{m\times m}$ depicts the remaining components of dimension $m\times m$ of the Nambu space. When vertex corrections are neglected, the partial derivative acts only on the kinetic energy term and not on the self-energy. Hence, we have defined the bare vertices

[TABLE]

The neglect of vertex corrections for $c-$axis superfluid stiffness is justified as follows. In the Green’s function, the small $c$-axis hopping amplitude $t_{\perp}$ comes in the DMFT self-consistency equation only through the lattice Green’s function, where it can be neglected compared with hopping in the plane. On the other hand, the current vertices coming from $t_{\perp}$ mean, in space-time, that the Green’s functions entering the particle-hole bubble for the superfluid stiffness are initially in different planes; these two Green’s functions do not allow hopping back to the same planes when $t_{\perp}$ is neglected. Since interactions are purely local in the Hubbard model, they can’t act on Green’s functions for electrons propagating in planes that are different. So the current vertex corrections for stiffness along the $c$-axis can be dropped out. The above argument also shows that the vertex corrections are of order $t_{\perp}^{2}$ compared with the leading terms. They can thus be neglected. Note that for the longitudinal response, which obeys the $f$-sum rule unlike the transverse response, the vector potential must be frequency dependent and, in addition, vertex corrections cannot be neglected (Schrieffer, 2018).

The $c$-axis hopping amplitude branches out into many different forms depending on the class of cuprates studied. This is because the matrix element for hopping between planes depends a lot on which orbitals overlap. We defer for the details to Refs. Chakravarty et al., 1993; Xiang and Wheatley, 1996; Panagopoulos et al., 1997; Andersen et al., 1995; Markiewicz et al., 2005; Gull and Millis, 2013. We chose a generic form describing $t_{\perp}$:

[TABLE]

where $t_{\text{bi}}\sim\frac{t}{25}=10$meV (Markiewicz et al., 2005; Panagopoulos et al., 1997). For the figures, we take $t_{\text{bi}}=1$, except when we show values for the penetration depth in physical units. In momentum space, from ARPES experiments at temperatures between pseudogap crossover $T^{\ast}$ and $T_{c}$, the structure of the pseudogap appears to mimic the essential features of the $d$-wave superconducting gap Loeser et al. (1996); Hoffman et al. (2002): the pseudogap is apparent only in the antinodal regions of the Brillouin zone where the $d$-wave gap is largest. Hence, the momentum dependence of $t_{\perp}$, of the form $(\cos{k_{x}}-\cos{k_{y}})^{2}$, suggests that the opening of the pseudogap in the CuO2 plane will lead to a large effect on the superfluid stiffness. The current vertices $\lambda_{\bm{k}}^{i}$ in Eq. (III.1) are obtained from the partial derivative along $z$ of $t_{\perp}(\bm{k})$.

To compute the London penetration depth $\lambda_{c}\equiv\lambda_{zz}$ along the $c$-axis in physical units, we set nearest-neighbor in-plane hopping to $t\sim 250$meV, lattice constants to $a=b\simeq 3.8\mbox{\AA}$ and $c\simeq 11.7\mbox{\AA}$ for the YBCO-like results, $a=b\simeq 3.8\mbox{\AA}$ and $c\simeq 13.2\mbox{\AA}$ for the NCCO-like results with $t_{\text{bi}}\sim 10$meV.

III.2 $d$SC regime

The superfluid stiffness without current vertex corrections comprising only $d$-wave superconductivity ($d$SC) reads (Gull and Millis, 2013; Schrieffer, 2018)

[TABLE]

where $\sigma_{3}$ is the diagonal Pauli matrix. The trace $\text{tr}\left[\ldots\right]$ operates on Nambu space $\hat{\Psi}_{\bm{k}}=\left(\hat{c}_{\bm{k},\uparrow}\ \hat{c}^{\dagger}_{-\bm{k},\downarrow}\right)^{\intercal}$. The current vertices give a contribution

[TABLE]

where the integral over $k_{z}$ can be performed because $t{{}_{\text{bi}}}$ is neglected in the Green’s functions. To compute $\bm{\rho}^{SC}_{zz}$ with the above formula, we first periodize the cluster Green’s function $\bm{\mathcal{G}}(k)$ using the full set of superlattice reciprocal wavevectors $\bm{K}_{i}$ (see Fig. 1). The periodized Green’s function is of size $2\times 2$.

III.3 Coexistence regime $d$SC + AF

We derived a formula to compute the superfluid stiffness in the regime where $d$-wave superconductivity and antiferromagnetism coexist homogeneously. This formula requires that one periodizes the cluster Green’s function to map onto the reduced AF Brillouin zone (AF-BZ).

First we define

[TABLE]

where $\sigma$ and $\tau$ are Pauli matrices, $\sigma^{l}$ acting in Nambu space and the $\tau^{m}$ acting in the AF sublattice space spanned by sublattices $A$ and $B$. We define $\sigma_{0}$ and $\tau_{0}$ as the identity matrix $\mathbb{I}_{2\times 2}$. The superfluid stiffness when AF and $d$SC coexist then reads ($\bm{\mathcal{G}}(k)\to\bm{\mathcal{G}}$):

[TABLE]

A detailed derivation of Eq. (III.3) is given in Appendix A. It can be extended to any phase coexistence scenario. In the above equation, the Green’s functions extracted from the CDMFT procedure are periodized using Eq. (11) with either $\bm{K}_{y}=\{(0,0),(0,\pi)\}$ or $\bm{K}_{x}=\{(0,0),(\pi,0)\}$ as the set of superlattice wavevectors (cf. Fig.1). The Green’s functions are then of dimension $4\times 4$ instead of $8\times 8$. Note that the dimension of the Green’s function in Eq. (9) is $8\times 8$ for unperiodized Green’s functions.

It is important to stress that the Green’s functions must be periodized in the AF-BZ prior to using Eq. (III.3) whether they have converged to a pure $d$SC state, a pure AF state, or a microscopic AF+$d$SC state. When the solution converges to a pure $d$SC state instead of microscopic AF+$d$SC, the superfluid stiffness obtained with either periodizations, namely Eq. (III.2) or (III.3), are indistiguishable on the plots.

IV Results

We study the superfluid stiffness $\rho_{zz}$ for a variety of parameters within the one-band Hubbard model Eq. (1), both with and without homogeneous microscopic AF+$d$SC coexistence. We find such coexistence in the CDMFT solutions of the cluster Green’s function only on the electron-doped side $(n>1)$ . The hole-doped side corresponds to band filling $n<1$. Whether antiferromagnetism is present or not, superconductivity is supressed at half-filling when the Hubbard interaction $U$ becomes larger than the value $U_{c}~{}\sim 6$ that leads to a Mott insulator (see Figs. 2 a), 2 b) and 3 a)). Overdoping means small $n$ for $n<1$ and large $n$ for $n>1$. In both cases, underdoping is near $n=1$.

We consider in turn band parameters that are close to those of YBCO and those of NCCO. The last subsection will show the effect of the $\bm{k}_{\parallel}$-dependence of $t_{\perp}$ in Eq. 19, giving us some insight on the parts of the Fermi surface that are most relevant for superconductivity.

IV.1 YBCO-like band parameters

Figures 2 a) and 2 b) illustrate both $\rho^{SC}_{zz}$ and $\rho^{AF+SC}_{zz}$ with respect to band filling $n$ per Cu 3$d_{x^{2}-y^{2}}$ orbital for the YBCO tight-binding parameters at $U=12t$ and $U=8t$, respectively.

The superfluid stiffness for both values of $U$ and for both hole- and electron-doping falls abruptly to zero in the overdoped regimes, where there is no coexistence. This suggests that in this limit, the system eventually reaches BCS-like behavior where at $T=0$ that sudden drop is expected. Finite resolution in the distance function, which contains an artifical temperature, probably explains why that drop is not perfectly discontinuous.

Contrast this BCS-like behavior with the behavior near half-filling for $U=12$ in Fig. 2 a) where the fall is much more gradual, as has been observed experimentally both along the $c$-axis and in the plane Panagopoulos et al. (1997); Uemura et al. (1989); Tallon et al. (2003). This is clearly the effect of the Mott transition since it does not appear when $U$ is not large enough ($U<U_{c}$), as can be seen in Fig. 3 b). The gradual fall of the superfluid stiffness has been interpreted as indicating that the superfluid density, as defined by the penetration depth Eq. (15), vanishes at half-filling and increases roughly proportionally to the doping, as if the number of carriers had to be measured with respect to half-filling.

Let us move to the effect of the competition with antiferromagnetism. Although there is no coexistence on the hole-doped side, antiferromagnetism is detrimental to superconductivity for $U=8$, as can be seen in Fig. 2 b) where the superconducting order parameter vanishes before half-filling is reached.

Comparing figures 2 a) and 2 b) in the region where there is coexistence, namely on the electron-doped side, we see that as $U$ increases, the domain of dopings where $d$SC and AF coexist grows. With increasing $U$, coexistence also ends at larger dopings when it reaches the pure superconducting phase near optimal doping. The jump in superfluid stiffness at this point is quite remarkable, but it may just reflect that the transition between the pure and coexisting phases is a first-order transition. The numerical values of $c$-axis superfluid stiffness given on the plots in physical units are of the correct order of magnitude compared with experimental measurements in cuprates Panagopoulos et al. (1997). Another remarkable property of superfluid stiffness in the coexistence region is that it is small and it decreases extremely rapidly as half-filling is approached. Coexistence ends relatively far away from half-filling.

In the coexistence regime, the converged cluster Green’s functions were periodized in the AF-BZ and used in Eq. (III.3). On the other hand, the superfluid stiffness in the pure $d$SC phase was computed by periodizing the cluster Green’s function in the BZ and used in Eq. (III.2). The agreement between both formulae when there is only a pure $d$SC phase (not shown on the figures) is non trivial and suggests that the superfluid stiffness formulae and the methods are consistent. This correspondence is also observed for calculations with the NCCO-like parameters.

Electron-doped material generally do not have band parameters close to those of YBCO. Electron-doped NCCO-like band parameters are explored in the next subsection.

IV.2 NCCO-like band parameters

Comparing calculations with experiments suggests that electron-doped cuprates, such as NCCO, are described by a Hubbard model with a value of $U$ in the vicinity of the Mott transition Sénéchal and Tremblay (2004); Kyung et al. (2004); Tremblay et al. (2006); Weber et al. (2010). The results for $\rho_{zz}$ appear in Figs. 3 a) and 3 b). Contrary to above, the discontinuity in $\rho_{zz}$ when antiferromagnetism appears near optimal doping has disappeared. The values of $U$ are quite close for the two plots, $U=6.55t$ in Fig. 3 a) and $U=5t$ in Fig. 3 b), leading to values of $\rho_{zz}$ that are quite close on the electron-doped side near optimal doping. While $\rho_{zz}$ looks continuous as a function of $n$ in Fig. 3 b) when antiferromagnetism appears upon decreasing doping, in the doped Mott insulator regime (see Fig. 3 a)), there is a rapid change in slope as a function of $n$ when antiferromagnetism appears.

Even though the values of $U$ in Figs. 3 a) and 3 b) are quite close, the difference between the two is quite striking. The case $U=5t$ in Fig. 3 b) is below $U_{c}$ for the Mott transition. This allows superconductivity to survive at half-filling when we do not allow antiferromagnetism to set in. The fall of $\rho_{zz}$ in the two overdoped regimes is abrupt, in BCS-like fashion. BCS would predict that $\rho_{zz}$ is proportional to band filling. Since the non-interacting Fermi surface is hole-like even on the electron-doped side, this is consistent with the increase in superfluid stiffness as $n$ decreases or hole density $\left|1-n\right|$ increases. The non-interacting van-Hove singularity where the Fermi surface becomes electron-like is at $n=0.8$, but this is shifted by interactions.

IV.3 Effect of the ${k}_{\parallel}$-dependence of $t_{\perp}$

Since $t_{\perp}$ is maximum at the $(\pi,0)$, $(0,\pi)$ points, as can be seen from $t_{\perp}^{2}(\bm{k})=t_{\text{bi}}^{2}\cos^{2}{k_{z}}\left(\cos{k_{x}}-\cos{k_{y}}\right)^{4}$ (Eq. (19)), an interesting question arises. Since both the pseudogap and the antiferromagnetic gap vary along the Fermi surface in the plane, their effect on $c$-axis superfluid stiffness $\rho_{zz}$ should be influenced by modulations of the $c$-axis hopping integral in the plane. What is the net effect of this modulation? The answer is in Figs. 4 a), 4 b).

We have computed $\rho_{zz}$ with and without the $\bm{k}_{\parallel}$-dependence of $t_{\perp}$. In the plots, by “no $t_{\perp}$”, we mean “in the absence of the $\bm{k}_{\parallel}$-dependence of the bilayer hopping term $t_{\perp}$”. In other words, we have replaced the in-plane modulation of perpendicular hopping $(\cos{k_{x}}-\cos{k_{y}})^{4}$ by $9/8$ since this is its average over the AF Brillouin zone. Figure 4 a) shows the effect of the $\bm{k}_{\parallel}$-dependence on $\rho_{zz}$ for YBCO band parameters, $U=12t$ and $n>1$. Figure 4 b) shows the same for NCCO band parameters, $U=6.55t$ and $n>1$. The results are qualitatively similar for the two sets of parameters.

The $\bm{k}_{\parallel}$-dependent term of $t_{\perp}$ (Eq. (19)) takes its maximum values in the portions of the Brillouin zone where $\bm{k}=(0,\pi)$ or $(\pi,0)$. These portions of the BZ are the antinodal regions where the $d$SC gap is the largest. Hence this is the region of the Brillouin zone that contributes most to the superfluid stiffness. In the coexistence region, which is electron-doped, the AF Fermi surface still has weight where the superconducting gap is largest. Hence, increasing the importance of these regions makes the superfluid stiffness larger. Also, the $\mathbf{k}_{\parallel}$-dependence of $t_{\perp}$ increases the contribution to $\rho_{zz}$ of the states in the vicinity of the van-Hove singularity on the hole-doped side (not shown). The latter can be checked by means of simple mean-field calculations.

By contrast, when superconductivity gaps the pseudogap normal state without coexisting antiferromagnetism, the situation is different. The pseudogap in the normal state is near $(\pm\pi/2,\pm\pi/2)$ in the electron-doped case. The superconductivity in that region is effective in lowering the energy because it replaces the pseudogap by quasiparticles. Hence, a more uniform weighting of the contributions across the Brillouin zone is more favourable in this case (not shown). This is also why the superfluid stiffness becomes larger without the $\mathbf{k}_{\parallel}$-dependence for $n>1.2$ in Fig. 4 a).

V Discussion

The $c$-axis striffness $\rho_{zz}$ has been calculated in Ref. Gull and Millis, 2013 using 8-site DCA for $U=6t>U_{c}$ and $\beta=60/t$ for particle-hole symmetric in-plane nearest-neighbor hopping. Their conclusions are qualitatively similar to the ones shown in Fig. 2 b) and 3 a): Mott physics suffices to lead to a vanishing superfluid stiffness as half-filling is approached and the fall of $\rho_{zz}$ when superconductivity disappears is more BCS-like in the overdoped regime. The authors noted that finite-temperature effects were likely to influence the results in the latter case, as also suggested in Ref. Metzner and Yamase, 2019.

As noted in the introduction, due to the finiteness of the cluster, the order parameter $\langle D\rangle$ indicates the formation of local Cooper pairs and the order parameter $\langle M\rangle$ indicates the formation of $(\pi,\pi)$ particle-hole bound states, or equivalently, local AF spin correlations.

V.1 Relation between $\rho_{s}$ and $T_{c}$

Before we discuss our results further, we comment on the possible relation between $\rho_{s}$ and $T_{c}$. Since cuprate superconductors are layered highly anisotropic superconductors, one expects Kosterlitz-Thouless physics (Nelson and Kosterlitz, 1977; Kosterlitz and Thouless, 1973) to play an important role in determining the actual superconducting transition temperature $T_{c}$ (Emery and Kivelson, 1995; Lee et al., 2008). This is brought out by experiments (Boschini et al., 2018; Yong et al., 2012; Hetel et al., 2007) even though there are contradictory views (Grbić et al., 2011; Chang et al., 2012; Tafti et al., 2014). The relevance of Kosterlitz-Thouless physics comes out clearly in the experimentally-determined proportionality between $T_{c}$ and $\rho_{s}$ in underdoped cuprates, both hole-doping(Uemura et al., 1989; Tallon et al., 2003) and electron-doping (Li et al., 2018).

Finite-temperature DCA studies with eight sites (Gull and Millis, 2013) and twelve sites (Maier and Scalapino, 2018) show that in the finite-temperature underdoped regime, increased phase fluctuations improve the agreement between the calculated and the observed shape of the superconducting transition-temperature dome. These theoretical results are consistent with the importance of long wave-length phase fluctuations, as in the Kosterlitz-Thouless theory.

When Kosterlitz-Thouless physics applies, $T_{c}$ can be bounded (Emery and Kivelson, 1995; Esterlis et al., 2018) using the zero-temperature value of $\rho_{s}$. However, it is the in-plane $\rho_{s}$ that counts. In the supplemental material (Note, 2), we compute that quantity, assuming that vertex corrections can be neglected. We find that the zero-temperature in plane $\rho_{s}$ scales with doping in a manner very similar to the $c$-axis stiffness. Further finite-temperature calculations and their relation to $T_{c}$ will be the subject of future work.

V.2 Hole-doped cuprates

Assuming that $\rho_{zz}$ scales with doping in the same way as the in-plane superfluid stiffness (Note, 2), our results on the hole-doped side of the phase diagram in Figs. 2 a) and 2 b) are consistent with the experimental drop of $\rho_{zz}(T=0)$ upon approaching half-filling in cuprates (Uemura et al., 1989, 1991). The quadratic component of the doping dependence that we found even seems consistent with the experimentally-inferred doping dependence of $\rho_{zz}(T=0)$ Hosseini et al. (2004). Earlier calculations of the in-plane superfluid stiffness that do not take into account antiferromagnetism explicitly have also found that superfluid stiffness vanishes near half filling when the Mott transition is taken into account. These calculations were done, for example, with slave-particles (Lee and Wen, 1997; Mallik et al., 2018) for the $t-J$ model or with variational Monte Carlo for the Hubbard model (Paramekanti et al., 2004).

The BCS-like drop in the superfluid stiffness that we find on the highly overdoped side is, however, not consistent with the linear doping dependence found experimentally in Refs. Božović et al., 2016; Mahmood et al., 2019 if we assume that the in-plane superfluid stiffness measured in these experiments behaves in the same way as $\rho_{zz}$ calculated here. It has however been argued theoretically that the behavior of the superfluid stiffness on the overdoped side is consistent with BCS dirty $d$-wave behavior (Lee-Hone et al., 2017, 2018).

At intermediate values of $U$, for example $U=8t$ in Fig. 2 b), antiferromagnetism plays an important role in making the superfluid stiffness vanish before half-filling. For larger clusters, it was found that superconductivity begins at a finite doping away from half-filling (Gull et al., 2013), even in the absence of antiferromagnetism. Nevertheless, comparing Fig. 3 b) for $U$ below the critical $U_{c}$ for the Mott transition with Figs. 2 a) and 2 b) for $U$ larger than $U_{c}$, it is clear that over most of the doping range the much smaller value of $\rho_{zz}$ and its doping dependence at large $U$ is controlled by Mott physics, not by competition with antiferromagnetism since antiferromagnetism appears only close to half-filling.

Note however that our cluster can’t accomodate long-period or incommensurate spin-density waves. These are seen both in experiments (Yamada et al., 1998; Wakimoto et al., 2000, 1999; Fujita et al., 2002; Haug et al., 2010) and in infinite-lattice calculations using methods that are valid for weak- (Yamase et al., 2016; Eberlein et al., 2016; Schulz, 1990) to intermediate-strength interaction (Y.M. Vilk and A.-M.S. Tremblay, 1997). A preprint that appeared as this paper was prepared (Metzner and Yamase, 2019) obtains results similar to ours in the hole-doped regime using mean-field parameters obtained from functional renormalization group. Even though the superfluid stiffness is similar to ours, its fall towards half-filling is caused by coexistence with commensurate antiferromagnetism. Results in the incommensurate regime were not presented. For $U>U_{c}$ our results suggest, but do not prove, that it is possible for superfluid stiffness to control $T_{c}$ in the underdoped regime even when there is no coexisting antiferromagnetism, contrary to the results for weak interaction strength Metzner and Yamase (2019).

V.3 Electron-doped cuprates

It is in electron-doped cuprates that competition with antiferromagnetism is strongest and it is there also that coexistence occurs in our calculations. Even though electronic-structure calculations Weber et al. (2010) and comparisons of theory Sénéchal and Tremblay (2004); Kyung et al. (2004) with photoemission Armitage et al. (2001, 2002) and neutron experiment Motoyama et al. (2007) show that the value of $U$ should be in close vicinity to the Mott transition, this is not crucial for qualitative features of $\rho_{zz}$ as a function of doping. They are quite similar in the case of electron doping for $U=12$ in Fig. 2 a), $U=8$ in Fig. 2 b) and $U=6.55$ in Fig. 3 a), which are all in the doped Mott insulator regime. In all cases: a) there is a small coexistence region where superfluid stiffness decreases rapidly compared with the value it would have in a pure superconducting state, b) antiferromagnetism overcomes completely superconductivity at a doping that is distinctly away from half-filling, as found in experiments Luke et al. (1990); Yu et al. (2007); Motoyama et al. (2007); Saadaoui et al. (2015) (see also note 333When antiferromagnetism is destroyed in the annealing process Horio et al. (2016), superconductivity ends very close to half-filling, which confirms that competition with antiferromagnetism is crucial in electron-doped superconductors), c) as one decreases doping, antiferromagnetism starts to coexist with superconductivity close to the doping where $\rho_{zz}$ reaches its maximum and d) the superfluid stiffness has a jump, or a rapid change in slope at smaller $U$, when one enters the coexistence phase from the pure superconducting phase at large doping. A similar jump was found in Ref. Metzner and Yamase, 2019.

One of the difficulties encountered by the one-band Hubbard or $t-J$ models has been that at zero temperature, when competition with long-range antiferromagnetic order is not allowed (pure regime), the size of the superconducting order parameter is larger on the electron-doped side of the phase diagram (Kancharla et al., 2008; White and Scalapino, 1999). This is quite clear in three bottom subfigures of Figs. 2 a) and 2 b). This was interpreted as a prediction that the transition temperature for electron-doped systems is larger than for hole-doped curates, contrary to observation. But one should not confuse the size of the superconducting order parameter with the value of the transition temperature. The correct value of the order parameters shows that competition with antiferromagnetic long-range order leads to the disappearance of superconductivity near half-filling, which decreases considerably the maximum value that the superconducting order-parameter can reach on the electron-doped side. Assuming that the in-plane superfluid stiffness scales with doping in a way similar to the $c$-axis results, one notices that the superfluid stiffness at optimal doping, that can be dominant in determining the value of $T_{c}$, is in all cases smaller on the electron-doped than on the hole-doped side of the phase diagram. In addition, in the actual materials, the value of $U$ should be somewhat smaller for electron-doped cuprates, as mentioned above.

VI Conclusion

We computed the $c$-axis superfluid stiffness at zero temperature for the one-band two-dimensional square-lattice Hubbard model. We solved the model on a $2\times 2$ plaquette using ED-CDMFT for model parameters appropriate for the cuprates. In finite-temperature $2\times 2$ plaquette calculations, the value of the superconducting transition temperature (Fratino et al., 2016a) indicates the formation of local pairs, not necessarily the actual transition temperature, which, as our calculations suggest, is controlled by superfluid stiffness in the underdoped regime.

On the hole-doped side, for YBCO band parameters and $U$ larger than the critical value for the Mott transition, it is mostly Mott physics that controls the value of the superfluid stiffness $\rho_{zz}$ near half-filling, although competition with antiferromagnetism does play a role just before half-filling. Superfluid stiffness along the $c$-axis increases with hole doping with linear plus quadratic dependence on doping, in qualitative agreement with experiment (Hosseini et al., 2004).

On the electron-doped side, our results suggest that it is the competition between AF and $d$SC that is most important even near optimal doping. This is suggested both by the value of the superconducting order parameter and by the superfluid stiffness $\rho_{zz}$ that jumps down Metzner and Yamase (2019) and then drops precipitously as soon as antiferromagnetism starts to coexist with superconductivity, a prediction for experiment. Just above that doping, $\rho_{zz}$ takes its largest value. This drop in $\rho_{zz}$ strongly depends on the electronic structure and on the value of the interaction $U$. The drop in $\rho_{zz}$ is more prominent for $U=8t$ and YBCO-like parameters. The reduction of $\rho_{zz}$ in the underdoped regime would increase the phase fluctuations of the superconducting order parameter. Hence, phase competition could be, according to the Uemura scaling relation (Uemura et al., 1989), an important factor in the depletion of $T_{c}$ in the underdoped regime for electron-doped cuprates as well. The disappearance of superconductivity closer to half-filling, however, comes from the fact that antiferromagnetism wins the competition with superconductivity in electron-doped cuprates.

For both hole- and electron-doping at large $U$, the superfluid stiffness jumps extremely quickly to zero when the system becomes normal in the overdoped regime, in qualitative agreement with the expected BCS behavior.

The effect of the in-plane modulation of the hopping amplitude along the $c$-axis is important: in the electron-doped case, at large $U$ on the electron-doped side it increases $\rho_{zz}$ in the regime where only superconductivity exists while it decreases it when there is coexistence with antiferromagnetism. This is understood in terms of where the $d$-wave superconducting gap is important compared with the underlying state.

We expect that competition with other types of order could have an effect on $\rho_{zz}$ similar to competition with antiferromagnetism. In future work, we plan to perform finite-temperature calculations to understand some of the unusual features of the superfluid-stiffness (Hosseini et al., 2004) and its more precise role in determining the transition temperature. Even though it has a phase diagram very similar to the one-band model Fratino et al. (2016b), we also plan to study the charge-transfer three-band model.

Acknowledgements.

We are grateful to Steve Kivelson, Giovanni Sordi and Simon Verret for discussions and to an anonymous referee for detailed constructive criticism. This work has been supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) under grant RGPIN-2014-04584, the Canada First Research Excellence Fund and by the Research Chair in the Theory of Quantum Materials. Simulations were performed on computers provided by the Canadian Foundation for Innovation, the Ministère de l’Éducation des Loisirs et du Sport (Québec), Calcul Québec, and Compute Canada.

Appendix A Derivation of $\rho_{ab}^{AF+SC}$

In this Appendix, we give further details on the calculation of the general superfluid stiffness $\rho_{zz}$ (III.1) for the CDMFT calculation and for the state with AF+dSC coexistence Eq. (III.3). We also explain the expression for the vertices (17) and how they are calculated when the perpendicular hopping amplitude depends on in-plane wave vectors.

A.1 General superfluid stiffness for CDMFT cluster

To understand how the tensors $T_{lmn}$ occur in our calculations, it is simpler to start from a mean-field-like Hamiltonian where hybridization functions are replaced by order parameters.

We start from the following lattice Green’s function

[TABLE]

where the 4-vector is defined by $\tilde{k}\equiv(\tilde{\bm{k}},i\omega_{n})$, $\mathcal{H}^{0}_{\tilde{\bm{k}}}$ is the mean-field-like Hamiltonian, i.e quadratic in field operators, and $\Sigma_{c}$ is the cluster self-energy.

In the CDMFT calculations, we consider a cluster consisting of 4 sites, therefore $N_{c}=4$ in the expression for the spinor that we use as a basis (Eq. (9)):

[TABLE]

Using the definitions in Fig. 5, the mean-field Nambu Hamiltonian then would be

[TABLE]

The superscripts in $\mathcal{H}^{\alpha\beta;ab;rs}_{{MF}}(\tilde{\mathbf{k}})$ take their meaning when the Hamiltonian is written as follows, taking advantage of the tensor-product form of the states on which the creation-annihilation operators (25) act:

[TABLE]

where $\sigma^{l}_{\alpha\beta}$, $\tau_{ab}^{m}$ and $\tilde{\sigma}_{rs}^{n}$ are Pauli and identity matrices and $A$ and $B$ are order parameter tensors when mean-field is used. This structure of the Hamiltonian allowed us to introduce for short-hand in subsection III.1 the tensor

[TABLE]

Equation (A.1) represents the Hamiltonian before periodization to the AF-BZ. This is why there is a Pauli matrix $\tilde{\sigma}$. It is always diagonal in our case.

We stress that we do not do mean-field theory. The effects of long-range order are all contained in the self-energy and hybridization function, not in the cluster Hamiltonian.

Following the linear response procedure in Eq. (III.1) using the Green’s function  (24), the formula obtained for the superfluid stiffness is

[TABLE]

The derivation will become clearer below when we consider the AF+dSC mean-field state.

The current vertices Eq. (17) are:

[TABLE]

They can be obtained from the gradient of the kinetic-energy part of the Hamiltonian because the phase of the Fourier transform within a unit cell was taken into account when writing the Hamiltonian Eq. (A.1) in the orbital basis Nourafkan and Tremblay (2018). The kinetic-energy part of the Hamiltonian, $\mathcal{H}_{\bar{\bm{k}},\sigma}^{0}$, is in the two 4$\times$4 diagonal blocks of Eq. (A.1). Because of the particle-hole transformation of down spins in the Nambu representation, we had to introduce a sign change through $\delta_{A_{i}}T_{000}=-\frac{e}{\hbar}\delta_{\bar{k}_{i}}T_{300}$ in Eq. (A.1) for the superfluid stiffness. In the CDMFT calculations, all off-diagonal terms are contained in the self-energy.

A.2 Superfluid stiffness in the AF+dSC coexistence state

Once again, the Pauli matrices $\sigma^{l}$ span the spin (Nambu) space, $\tau^{m}$ the AF sublattice space and $\tilde{\sigma}$ the sublattice spin projection degrees of freedom. In the infinite lattice with coexisting AF+dSC, the subspace spanned by $\tilde{\sigma}$ is not relevant. When we periodize the cluster Green’s function onto the AF-BZ the superfluous $\tilde{\sigma}$ subspace disappears.

Let us go in more details through the derivation of $\rho_{zz}$ for a mean-field with AF+dSC microscopic coexistence. It will be clear how to handle the case of the periodized CDMFT Green’s function Eq. (7) mapped onto the AF-BZ (Fig. 1).

The mean-field AF+dSC Hamiltonian Eq. (A.1) would read, with $\{i,j\}=\{A,B\}$ and $\{\alpha,\beta\}=\{\uparrow,\downarrow\}$:

[TABLE]

where $t_{ab}$ is the hopping matrix between different AF sublattices, $\bm{Q}=(\pi,\pi)$ is the AF nesting wavevector and $\Delta_{ab}=\Delta$ if $\bm{r}_{a}-\bm{r}_{b}=\pm\bm{e}_{x}$, and $\Delta_{ab}=-\Delta$ if $\bm{r}_{a}-\bm{r}_{b}=\pm\bm{e}_{y}$, corresponding to $d_{x^{2}-y^{2}}$ pairing. The momentum vector $\bm{k}$ is defined in the rBZ.

In orbital Nambu basis, the matrix form of the mean-field hamiltonian Eq. (A.2) is such that it can be written in terms of the SU$(2)\otimes$SU$(2)$ matrices $T_{lm}$ defined in Eq. (22):

[TABLE]

For the current $\langle j_{a}\rangle=-\frac{1}{V}\frac{\delta\mathcal{F}}{\delta A_{a}}$, one first needs

[TABLE]

The bar over $\bm{k}$ reminds us that the vector potential is contained in the wave-vector with a sign that differs between up and down spins. The current then can be written as

[TABLE]

where one can use either the mean-field or the periodized CDMFT Green’s function and where $\text{tr}[\cdots]$ operates in the 4$\times$4 Nambu space. We have supposed that the system is invariant under inversion ($\bm{k}=-\bm{k}$).

The periodized CDMFT Green’s function takes the same form as in Eq. 24 except that $\bm{\mathcal{H}}^{0}_{\bm{k}}$ depends on $\bm{k}$ instead of $\bm{\tilde{k}}$ and has a smaller size since it contains only the $T_{30}$ and $T_{31}$ parts of the mean-field Hamiltonian Eq. (33). All off-diagonal pieces are in the self-consistent off-diagonal self-energies.

Inserting either the mean-field or CDMFT periodized Green’s function, neglecting vertex corrections (i.e. the self-energy dependence of the vector potential $\bm{A}$), and using

[TABLE]

one can compute \bm{\rho}_{ab}=-\frac{\delta\langle j_{a}\rangle}{\delta A_{b}}\big{\rvert}_{\bm{A}=\bm{0}}:

[TABLE]

The second term, so-called paramagnetic term, was obtained from the derivative of the Green’s function Eq. (36). Once the partial derivatives have acted, we set $\bm{A}\to\bm{0}$.

It is convenient to use the periodicity of the Brillouin zone to do a partial integration of the diamagnetic components

[TABLE]

That allows a grouping of the diamagnetic and paramagnetic terms of Eq. (A.2), Indeed, benefiting again from Eq. (36) for the derivative of a Green’s function (with the replacement $A_{b}\to k_{b}$), the final expression for the superfluid stiffness in the AF-$d$SC coexistence regime takes the form ($\bm{\mathcal{G}}(k)\to\bm{\mathcal{G}}$)

[TABLE]

Equation (A.2) is general if vertex corrections are neglected and $a,b\in\{x,y,z\}$. The Green’s functions obtained from periodizing the CDMFT solutions can be introduced where $\bm{\mathcal{G}}$ stands in Eq. (A.2). If one does not allow for symmetry breaking associated with antiferromagnetism, one retrieves the superfluid stiffness formula for pure superconducting systems Eq. (III.2).

A.3 Vertices for $\rho_{zz}$ when there is a dependence on in-plane wave vectors

Since we compute $\rho_{zz}$, we do not need all the terms of Eq. (A.2). Fourier transforming the perpendicular hopping Eq. (19) back to lattice coordinates, one can see that there are three different interlayer hopping terms involved in Eq. (19) and they all make the electrons hop from one AF sublattice to the other, as can be seen from Fig. 6. Hence, only the last term of Eq. (A.2) remains after setting $a=b=z$.