Fate of disorder-induced inhomogeneities in strongly correlated d-wave superconductors

We describe the phases of the strongly correlated cuprates by the Hubbard model, which is a minimal lattice model representing repulsively interacting band electrons:

$$\begin{eqnarray}&&{{\mathcal{H}}_{{\rm Hubb}}}=\mathop{\sum }\limits_{ij\sigma }{{t}_{ij}}(c_{i\sigma }^{\dagger }{{c}_{j\sigma }}+h.c.)+U\mathop{\sum }\limits_{j}{{n}_{j\uparrow }}{{n}_{j\downarrow }}.\end{eqnarray} \tag{ 1 }$$

Here, the first term indicates electrons hopping on a two-dimensional square lattice, with $c_{i\sigma }^{\dagger }({{c}_{i\sigma }})$ as the creation (annihilation) operator of the electron on site i and with spin σ. We take ${{t}_{ij}}=-t$, when i and j are nearest neighbors, denoted as $\langle ij\rangle$, and ${{t}_{ij}}=t^{\prime}$, when i and j are next-nearest neighbors, with the notation of $\langle \langle ij\rangle \rangle$. We choose ${{t}_{ij}}=0$ for all other pairs of i and j. U is the onsite interaction energy between electrons occupying the same site and mimics a Coulomb-type repulsion where inter-site contributions are neglected in favor of screening. The local density operator with spin σ at site j is denoted as ${{n}_{j\sigma }}=c_{j\sigma }^{\dagger }{{c}_{j\sigma }}$.

In the strongly correlated limit $U\gg t$, the Schrieffer–Wolff transformation on ${{\mathcal{H}}_{{\rm Hubb}}}$ yields an effective 't–J' model [63]:

$$\begin{eqnarray}\begin{array}{rcl} {{\mathcal{H}}_{t-J}} & = & \mathop{\sum }\limits_{ij\sigma }{{t}_{ij}}(\tilde{c}_{i\sigma }^{\dagger }{{\tilde{c}}_{j\sigma }}+h.c.)+\mathop{\sum }\limits_{ij}{{J}_{ij}}\left( {{\tilde{{\bf S}}}_{i}}.{{\tilde{{\bf S}}}_{j}}-\frac{{{\tilde{n}}_{i}}{{\tilde{n}}_{j}}}{4} \right) \\ {} & {} & -\frac{J}{4}\mathop{\sum }\limits_{\frac{\langle ijm\rangle ,\sigma }{m\ne i}}(\tilde{c}_{i\sigma }^{\dagger }{{\tilde{n}}_{j\bar{\sigma }}}{{\tilde{c}}_{m\sigma }}-\tilde{c}_{i\sigma }^{\dagger }\tilde{c}_{j\bar{\sigma }}^{\dagger }{{\tilde{c}}_{j\sigma }}{{\tilde{c}}_{m\bar{\sigma }}}+h.c.) \\ {} & {} & -\frac{J^{\prime} }{4}\mathop{\sum }\limits_{\frac{\langle \langle ijm\rangle \rangle ,\sigma }{m\ne i}}(\tilde{c}_{i\sigma }^{\dagger }{{\tilde{n}}_{j\bar{\sigma }}}{{\tilde{c}}_{m\sigma }}-\tilde{c}_{i\sigma }^{\dagger }\tilde{c}_{j\bar{\sigma }}^{\dagger }{{\tilde{c}}_{j\sigma }}{{\tilde{c}}_{m\bar{\sigma }}}+h.c.) \\ {} & {} & +\frac{2tt^{\prime} }{U}\mathop{\sum }\limits_{\frac{\sigma ,\langle ij\rangle }{\langle \langle jm\rangle \rangle }}(\tilde{c}_{i\sigma }^{\dagger }{{\tilde{n}}_{j\bar{\sigma }}}{{\tilde{c}}_{m\sigma }}-\tilde{c}_{i\sigma }^{\dagger }\tilde{c}_{j\bar{\sigma }}^{\dagger }{{\tilde{c}}_{j\sigma }}{{\tilde{c}}_{m\bar{\sigma }}}+h.c.), \\ \end{array}\end{eqnarray} \tag{ 2 }$$

where all terms up to $\mathcal{O}({{t}^{2}}/U)$ are kept, and ${{J}_{ij}}=4t_{ij}^{2}/U$. Here, we use the notation ${{J}_{ij}}=J$ for $\langle ij\rangle$, and ${{J}_{ij}}=J^{\prime}$ for $\langle \langle ij\rangle \rangle$, and ${{\tilde{c}}_{i\sigma }}={{c}_{i\sigma }}(1-{{n}_{i\bar{\sigma }}})$ is the annihilation operator in the 'projected space' prohibiting double occupancy at site i. The terms appearing in the second to fourth lines of equation (2) represent the 'three-site terms', which, along with the $J^{\prime}$-term the are often neglected in standard $t-J$ model, even though they are of the same order as the J-term. In fact, these terms are capable of generating additional broken symmetry phases in the Cooper channel, as we discuss in section 2.2. We used the quotes in 't–J' to emphasize that it makes contributions beyond the standard $t-J$ model. The notations $\langle ijm\rangle$ and $\langle \langle ijm\rangle \rangle$ are for the three sites with (i, j) and (j, m) being nearest neighbors and next-nearest neighbors to each other, respectively.

We introduce disorder by redefining $\mathcal{H}$_'t–J' to ${{\mathcal{H}}_{t-J}}+{{\sum }_{i\sigma }}({{V}_{i}}-\mu ){{n}_{i\sigma }}$, where V_i is the (non-magnetic) impurity potential at site i and μ is the chemical potential that fixes the average density of electrons in the system to a desired value. There are different ways of choosing V_iʼs, depending on the specifics of the experimental realizations. It is, however, important to quantify the degree of disorder by a small number of parameters defining the distribution of V_i for an appropriate comparison with the experiments.

For our calculations, we used three models of disorder: (i) box disorder: V_iʼs are drawn from a 'box' distribution, such that ${{V}_{i}}\in [-V/2,V/2]$ uniformly; (ii) concentration disorder (or conc disorder): ${{V}_{i}}={{V}_{0}}$ on a given ${{n}_{{\rm imp}}}$ fraction of the sites, which are randomly chosen, and V_i = 0 for all other sites; and (iii) out-of-plane disorder, as detailed in the supplementary material (available from stacks.iop.org/njp/16/103018/mmedia). Although the box disorder is quantified by a single parameter V, the conc disorder needs two parameters for its characterization: V₀ and ${{n}_{{\rm imp}}}$. We will focus, in particular, on the box disorder, which is better suited to address several conceptual issues, as we will see in section 3.3. All data are averaged over 10–15 independent realizations for a given strength of disorder.

Strong electronic repulsions restrict the Hilbert space of ${{\mathcal{H}}_{t-J}}$ by eliminating all double occupancies. One simple method for implementing such a restriction is the GA, in which the higher-order off-site correlations are neglected [46]. Following reference [54], we consider the GS wave function in the projected space to be $|\psi \rangle =P|{{\psi }_{0}}\rangle$, where $|{{\psi }_{0}}\rangle$ is the GS wave function in the Hilbert space that allows double occupancy and P is the projection operator, $P={{\prod }_{i}}{{P}_{i}}$, with ${{P}_{i}}=\gamma _{i}^{{{n}_{i}}/2}(1-{{n}_{i\uparrow }}{{n}_{i\downarrow }})$. Here, γ_i are the local fugacity factors obtained by demanding conservation of the local electron densities, ${{\rho }_{i\sigma }}=\langle {{n}_{i\sigma }}\rangle ={{\langle {{n}_{i\sigma }}\rangle }_{0}}$. In GA, the effects of projection on impure dSC are absorbed by local Gutzwiller renormalization factors [53]. The expectation value of a general operator $\hat{A}$ in different GS wave functions can be written as ${{\langle \hat{A}\rangle }_{0}}=\langle {{\psi }_{0}}|\hat{A}|{{\psi }_{0}}\rangle /\langle {{\psi }_{0}}|{{\psi }_{0}}\rangle$ and $\langle \hat{A}\rangle =\langle \psi |\hat{A}|\psi \rangle /\langle \psi |\psi \rangle$, and GA amounts to the following simple relations:

$$\begin{eqnarray}&&\langle c_{i\sigma }^{\dagger }{{c}_{j\sigma }}\rangle \approx g_{ij}^{t}{{\langle c_{i\sigma }^{\dagger }{{c}_{j\sigma }}\rangle }_{0}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\langle c_{i\sigma }^{\dagger }{{n}_{j\bar{\sigma }}}{{c}_{m\sigma }}\rangle \approx g_{ijm}^{31}{{\langle c_{i\sigma }^{\dagger }{{n}_{j\bar{\sigma }}}{{c}_{m\sigma }}\rangle }_{0}},\;\langle c_{i\sigma }^{\dagger }c_{j\bar{\sigma }}^{\dagger }{{c}_{j\sigma }}{{c}_{m\bar{\sigma }}}\rangle \approx g_{ijm}^{32}{{\langle c_{i\sigma }^{\dagger }c_{j\bar{\sigma }}^{\dagger }{{c}_{j\sigma }}{{c}_{m\bar{\sigma }}}\rangle }_{0}},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\langle {{{\bf S}}_{i}}.{{{\bf S}}_{j}}\rangle \approx g_{ij}^{xy}{{\langle {{{\bf S}}_{i}}.{{{\bf S}}_{j}}\rangle }_{0}},\;{\rm however},\;\langle {{n}_{i}}{{n}_{j}}\rangle \approx {{\langle {{n}_{i}}{{n}_{j}}\rangle }_{0}}.\end{eqnarray} \tag{ 5 }$$

Here g^t_ij, g³¹_ijm, g³²_ijm, and g^xy_ij are the Gutzwiller renormalization factors (GRF) which depend on an appropriate combinations of local doping ${{x}_{i}},{{x}_{j}},{{x}_{m}}$ (here ${{x}_{i}}=1-{{\rho }_{i}}$). The GRFs for our system (without any magnetic order) are given by

$$\begin{eqnarray}&&g_{ij}^{t}=\sqrt{\frac{4{{x}_{i}}{{x}_{j}}}{(1+{{x}_{i}})(1+{{x}_{j}})}},\;\;g_{ij}^{xy}=\frac{4}{(1+{{x}_{i}})(1+{{x}_{j}})},\;\;g_{ijm}^{31}=g_{im}^{t},\;\;g_{ijm}^{32}=\frac{2\;g_{im}^{t}}{(1+{{x}_{j}})}.\end{eqnarray} \tag{ 6 }$$

We emphasize that since the GRFs are different for ${{{\bf S}}_{i}}.{{{\bf S}}_{j}}$ (renormalized by g^xy_ij) and ${{n}_{i}}{{n}_{j}}$ (not renormalized) terms in ${{\mathcal{H}}_{t-J}}$, a flat renormalization of the exchange coupling, J, would not suffice.

Upon complete removal of the double occupancy by GA, the total energy of the system can be calculated in terms of the relevant local-order parameters: ${{\rho }_{i}}={{\sum }_{\sigma }}\langle {{n}_{i\sigma }}\rangle \equiv {{\sum }_{\sigma }}{{\langle {{n}_{i\sigma }}\rangle }_{0}}$, ${{\Delta }_{ij}}\equiv {{\langle {{c}_{j\downarrow }}{{c}_{i\uparrow }}\rangle }_{0}}+{{\langle {{c}_{i\downarrow }}{{c}_{j\uparrow }}\rangle }_{0}}$, and ${{\tau }_{ij}}\equiv {{\langle c_{i\downarrow }^{\dagger }{{c}_{j\downarrow }}\rangle }_{0}}\equiv {{\langle c_{i\uparrow }^{\dagger }{{c}_{j\uparrow }}\rangle }_{0}}$. We minimize $\langle {{\psi }_{0}}|{{\mathcal{H}}_{t-J}}|{{\psi }_{0}}\rangle$, with respect to $|{{\psi }_{0}}\rangle$ with additional constraint $\langle {{\psi }_{0}}|{{\psi }_{0}}\rangle =1$ and ${{\sum }_{i}}{{\rho }_{i}}=\rho$ (ρ being the chosen total density), leading to the minimization of the functional W, where

$$\begin{eqnarray}&&W=\langle {{\psi }_{0}}|{{\mathcal{H}}_{t-J}}|{{\psi }_{0}}\rangle -\beta (\langle {{\psi }_{0}}|{{\psi }_{0}}\rangle -1)-\mu (\mathop{\sum }\limits_{i}{{\rho }_{i}}-\rho ),\end{eqnarray} \tag{ 7 }$$

employing the variational relation $(\delta W)/(\delta \langle {{\psi }_{0}}|)=0$. Here β and μ are the Lagrangian multipliers for the constraints mentioned previously. Upon minimizing W, we obtain the following BdG mean-field Hamiltonian [47, 57, 58, 60]

$$\begin{eqnarray}&&{{H}_{{\rm MF}}}=\mathop{\sum }\limits_{ij\sigma }\frac{\partial W}{\partial {{\tau }_{ij}}}c_{i\sigma }^{\dagger }{{c}_{j\sigma }}+\mathop{\sum }\limits_{ij}\frac{\partial W}{\partial {{\Delta }_{ij}}}(c_{i\uparrow }^{\dagger }c_{j\downarrow }^{\dagger }+c_{j\uparrow }^{\dagger }c_{i\downarrow }^{\dagger })+\mathop{\sum }\limits_{i\sigma }\frac{\partial W}{\partial {{\rho }_{i\sigma }}}{{n}_{i\sigma }},\end{eqnarray} \tag{ 8 }$$

satisfying the eigenvalue equation: ${{H}_{{\rm MF}}}|{{\psi }_{0}}\rangle =\beta |{{\psi }_{0}}\rangle$. Upon expanding the derivatives in equation (8), we obtain the explicit form of the mean-field Hamiltonian:

$$\begin{eqnarray}\begin{array}{rcl} {{\mathcal{H}}_{{\rm MF}}} & = & {{\mathop{\sum }\limits_{i,\delta ,\sigma }}^{\prime }}\{{{t}_{i\delta }}g_{i,i+\delta }^{t}-W_{i\delta }^{{\rm FS}}\}\;c_{i\sigma }^{\dagger }{{c}_{i+\delta \sigma }}+{{\mathop{\sum }\limits_{i,\delta }}^{\prime }}\{G_{1}^{i,\delta }\Delta _{i}^{\delta }(c_{i\uparrow }^{\dagger }c_{i+\delta \downarrow }^{\dagger }+c_{i+\delta \uparrow }^{\dagger }c_{i\downarrow }^{\dagger })+h.c.\} \\ {} & {} & +{{\mathop{\sum }\limits_{i,\delta ,\delta ^{\prime} }}^{\prime }}\{G_{2}^{i,\delta ,\delta ^{\prime} }\Delta _{i}^{\delta }(c_{i\uparrow }^{\dagger }c_{i+\delta ^{\prime} \downarrow }^{\dagger }+c_{i+\delta ^{\prime} \uparrow }^{\dagger }c_{i\downarrow }^{\dagger })+h.c.\} \\ {} & {} & +\mathop{\sum }\limits_{i,\delta ,\tilde{\delta }}\{G_{3}^{i,\delta ,\tilde{\delta }}\Delta _{i}^{\tilde{\delta }}(c_{i\uparrow }^{\dagger }c_{i+\delta \downarrow }^{\dagger }+c_{i+\delta \uparrow }^{\dagger }c_{i\downarrow }^{\dagger })+h.c.\} \\ {} & {} & +\mathop{\sum }\limits_{i,\delta ,\tilde{\delta }}\{G_{3}^{i,\delta ,\tilde{\delta }}\Delta _{i}^{\delta }(c_{i\uparrow }^{\dagger }c_{i+\tilde{\delta }\downarrow }^{\dagger }+c_{i+\tilde{\delta }\uparrow }^{\dagger }c_{i\downarrow }^{\dagger })+h.c.\} \\ {} & {} & +{{\mathop{\sum }\limits_{\frac{i,\delta }{\delta ^{\prime} ,\sigma }}}^{\prime }}G_{4}^{i,\delta ,\delta ^{\prime} }c_{i+\delta \sigma }^{\dagger }{{c}_{i+\delta ^{\prime} \sigma }}+\mathop{\sum }\limits_{i,\sigma }({{V}_{i}}-\mu +\mu _{i}^{{\rm HS}}){{n}_{i\sigma }} \\ {} & {} & +\mathop{\sum }\limits_{i,\delta ,\tilde{\delta },\sigma }G_{5}^{i,\delta ,\tilde{\delta }}\{c_{i+\delta \sigma }^{\dagger }{{c}_{i+\tilde{\delta }\sigma }}+h.c.\}. \\ \end{array}\end{eqnarray} \tag{ 9 }$$

We denote the Cooper-channel order parameters Δ_ij on the bond connecting sites i and j as $\Delta _{i}^{\delta }$, where $j=i+\delta$, with $\delta =\pm x\;{\rm or}\pm y$, $\tilde{\delta }=\pm (x\pm y)$, and the primed summation means.

$$\begin{eqnarray}&&{{\mathop{\sum }\limits_{i,\delta ,\sigma }}^{\prime }}\equiv \mathop{\sum }\limits_{i,\delta ,\sigma }+\mathop{\sum }\limits_{i,\tilde{\delta },\sigma }\;{\rm and}\;{{\mathop{\sum }\limits_{i,\delta ,\delta ^{\prime} }}^{\prime }}\equiv \mathop{\sum }\limits_{i,\delta ,\delta ^{\prime} }+\mathop{\sum }\limits_{i,\tilde{\delta },\tilde{\delta }^{\prime} }\end{eqnarray} \tag{ 10 }$$

The coefficients G₁ to G₅ in equation (9) are built with the GRFs defined in equation (6) and their exact expressions are included in appendix A. $\mu _{i}^{{\rm HS}}$ and $W_{i\delta }^{{\rm FS}}$ are the Hartree- and Fock-shift terms arising from the mean-field decomposition of interactions and are also given in appendix A. We found that the major contributions of $\mu _{i}^{{\rm HS}}$ come from the derivatives of the GRF for the parameters describing cuprates.

We diagonalize ${{\mathcal{H}}_{{\rm MF}}}$ by using the BdG transformations [64], ${{c}_{i\sigma }}={{\sum }_{n}}({{\gamma }_{n,\sigma }}{{u}_{i,n}}-\sigma \gamma _{n,-\sigma }^{\dagger }v_{i,n}^{*})$, where $\gamma _{n,\sigma }^{\dagger }$ and ${{\gamma }_{n,\sigma }}$ are the creation and annihilation operators of the Bogoliubov quasiparticles. The resulting eigen-system is then solved self-consistently for all the local-order parameters. Inclusion of $J^{\prime}$ and the three-site terms in the Hamiltonian allows additional Cooper-channel orders in ${{\mathcal{H}}_{{\rm MF}}}$, namely, s_xy and d_xy, in addition to the standard d-wave (${{d}_{{{x}^{2}}-{{y}^{2}}}}$) and the extended s-wave (s_xs), which take the following forms in the clean systems:

$$\begin{eqnarray}&&{{\Delta }_{{{d}_{{{x}^{2}}-{{y}^{2}}}}}}(\vec{k})=\Delta _{{{d}_{{{x}^{2}}-{{y}^{2}}}}}^{(0)}({\rm cos} {{k}_{x}}-{\rm cos} {{k}_{y}})/2,\;{{\Delta }_{{{s}_{xs}}}}(\vec{k})=\Delta _{{{s}_{xs}}}^{(0)}({\rm cos} {{k}_{x}}+{\rm cos} {{k}_{y}})/2\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{{\Delta }_{{{s}_{xy}}}}(\vec{k})=\Delta _{{{s}_{xy}}}^{(0)}({\rm cos} {{k}_{x}}{\rm cos} {{k}_{y}}),\;{{\Delta }_{{{d}_{xy}}}}(\vec{k})=\Delta _{{{d}_{xy}}}^{(0)}({\rm sin} {{k}_{x}}{\rm sin} {{k}_{y}}),\end{eqnarray} \tag{ 12 }$$

where $\Delta _{{{d}_{{{x}^{2}}-{{y}^{2}}}}}^{(0)}$, $\Delta _{{{s}_{xs}}}^{(0)}$, $\Delta _{{{s}_{xy}}}^{(0)}$, and $\Delta _{{{d}_{xy}}}^{(0)}$ are the uniform values of the pairing amplitudes with respective symmetries. For notational simplification, we denote ${{\Delta }_{{{d}_{{{x}^{2}}-{{y}^{2}}}}}}$ by Δ_d for the rest of the paper. Starting with initial guesses for all the local order parameters, we numerically establish the eigenvalues and eigenfunctions of the BdG Hamiltonian. With these, we recalculate the order parameters and compare them with the initial guesses. If the two do not match on all the sites, the whole process is iterated with a new choice of the order parameters until self-consistency is achieved.

We investigated the physics of ${{\mathcal{H}}_{{\rm MF}}}$ within the framework of GIMT for a wide range of parameters. In the next section, we will describe our results for T = 0 for U = 12 and $t^{\prime} =1/4$ in the units of t for all our GIMT calculations, which are believed to describe a typical cuprate [65]. All other energies are also expressed in units of t. Typical size of the unit cell for our calculations has been 30 × 30, and we focus on the results for the average density of electrons, $\rho =0.8$, which coincides with near, optimal doping in the cuprates. In general, HTSC materials show a low temperature symmetry breaking primarily in the Cooper channel at the optimal doping and hence we do not consider other orders [66–70].

For the box disorder, we tune the disorder strength up to V = 3, which is sufficient to destroy dSC in the IMT calculations, as we will see in section 3.2. The chosen strength of the conc disorder is ${{V}_{0}}=1$ for most of our calculations, setting them to the 'Born scattering limit'. Although we report some results even for ${{V}_{0}}\sim 5-10$ to address the effect of strong scatterers, we avoid the unitary limit [57], e.g., in the case of substitutionary zinc impurities [71]. These impurities arguably produce a local antiferromagnetic surrounding [57] and can lead to effects beyond the scope of our study.

To carry out an appropriate comparison, a value of U = 3.69 is chosen for the IMT calculation which yield the same homogeneous (V = 0) d-wave gap from the GIMT scheme with the cuprate parameters. An alternative IMT scheme could also be set up and is discussed in the supplementary material. The ${{\mathcal{H}}_{{\rm MF}}}$ for the plain IMT is recovered by setting the GRF in equation (6) to unity; however, it presents a technical problem. The three-site terms deplete the ${{d}_{{{x}^{2}}-{{y}^{2}}}}$ order heavily, leading to a cancellation of much of the d-wave pairing amplitudes. The GIMT calculations, however, remain free from such a deficit. Further, the s_xy and d_xy orders from the $J^{\prime}$ term were found to be insignificant in the GIMT calculations, as is illustrated in figure 1(c). Thus, we present the bulk of our results by using just the $t-t^{\prime} -J$ model for the GIMT calculations with the box disorder and for the IMT calculations with both the box and the conc disorder, dropping the three-site terms and the $J^{\prime}$ term in favor of simplicity. Our GIMT results with the conc disorder allow all the Cooper-channel orders since they use the full 't – J' model.

Zoom In Zoom Out Reset image size

We plot along the y-axis of figures 1(a) and (b) the number of times a value of ρ occurs anywhere in the system against ρ itself along the x-axis from the GIMT and IMT findings. The resulting normalized distribution $P(\rho )$ widens with V; however, the disorder dependence of $P(\rho )$ is markedly different in the two calculations. For example, at V = 2.5, ${{P}_{{\rm GIMT}}}(\rho )$ is finite only for $\rho \in [0.6,0.95]$, whereas ${{P}_{{\rm IMT}}}(\rho )$ is non-zero over a much wider range of $\rho =0.1$ to $\rho =1.45$. The Gutzwiller projection ensures ${{\rho }_{i}}\leqslant 1$. It is interesting to note that ${{P}_{{\rm GIMT}}}(\rho )$ develops an asymmetry with increasing V with a larger weight towards larger values of density. Such asymmetry is absent in ${{P}_{{\rm IMT}}}(\rho )$. The explanation for this feature is discussed in section 3.3. Results from the conc disorder (which includes three-site terms and $J^{\prime}$ terms in GIMT) are broadly similar to those for the box disorder with respect to the inhomogeneities and are presented as the insets in figures 1(a) and (b). As expected, $P(\rho )$ develops finite weight for depleted values of ρ on the disorder sites (at $\rho \sim 0.7$ for GIMT and $\rho \sim 0.4$ for the IMT calculations) in addition to the large contribution around its average, $\rho \approx 0.8$. It is interesting that the electron population on the impurities is much larger in GIMT than in plain IMT. This ensures that the electronic correlations prohibit strong inhomogeneities. Further, it reduces the effective disorder strength, if we were to define this strength locally in terms of depleting the population of electrons in our context of repulsive scatterers. The reduction of the effective disorder has serious consequences for the physical properties of the superconductor and is addressed in section 3.4.

The local density plays a crucial role in determining the nature of the local d-wave order parameter. The double occupant sites (and empty sites) in the IMT calculation deplete the ${{d}_{{{x}^{2}}-{{y}^{2}}}}$-wave pairing amplitude locally, defined by ${{\Delta }_{d}}(i)=\frac{1}{4}(\Delta _{i}^{+x}-\Delta _{i}^{+y}+\Delta _{i}^{-x}-\Delta _{i}^{-y})$. This is because these sites have extreme values of local density (namely, ${{\rho }_{i}}=0$ or 2) and do not participate in the particle-hole mixing, which is essential for a non-zero value to ${{\Delta }_{d}}(i)$. As a result, the d-wave pairing amplitude from the plain IMT reduces with increasing V (box disorder) or with ${{n}_{{\rm imp}}}$ (conc disorder). The corresponding distribution, ${{P}_{{\rm IMT}}}({{\Delta }_{d}})$, as shown in figure 1(d), broadens with increasing V with a large amount of weight progressively shifting towards a smaller value. Such an evolution finally leads to a rather skewed ${{P}_{{\rm IMT}}}({{\Delta }_{d}})$ at large V, where a large number of sites have ${{\Delta }_{d}}(i)\approx 0$, as shown in figure 1(d). In fact, this same physical picture [59, 72] yielded a very similar disorder dependence of ${{P}_{{\rm IMT}}}({{\Delta }_{s}})$ for a disordered sSC. In contrast, the d-wave pairing amplitude in the GIMT calculation, which is free of doubly occupied sites (and hence, free of near-empty sites, too, thereby satisfying the constraint of maintaining $\rho =0.8$), stays around its homogeneous value (see figure 1(c)).

Of course, there are inhomogeneities in the GIMT results, which give the increasing width in ${{P}_{{\rm GIMT}}}({{\Delta }_{d}})$ with V. However, the fluctuations in the pairing amplitude are weak (at least for moderate V) and symmetric about $\Delta _{d}^{(0)}$, yielding a ${{P}_{{\rm GIMT}}}({{\Delta }_{d}})$ nearly Gaussian in shape for all V. Although $P(\rho )$ indicates that the degree of inhomogeneity is different in the GIMT and IMT schemes, the results for $P({{\Delta }_{d}})$ also establish that the nature of inhomogeneity is substantially different in the two calculations. Such a disorder evolution of ${{P}_{{\rm GIMT}}}({{\Delta }_{d}})$ is strikingly different from the predictions of the simple Abrikosov–Gorkov theory or its extension in the form of a self-consistent T-matrix approximation [15, 73], in which Δ_d (assumed homogeneous) depletes monotonically with V.

It is interesting to note that the GIMT calculation yields ${{\Delta }_{d}}(i)\geqslant \Delta _{d}^{(0)}$ at certain locations, particularly at large V, whereas the IMT value of ${{\Delta }_{d}}(i)$ is more or less limited by $\Delta _{d}^{(0)}$. The similarity of our ${{P}_{{\rm GIMT}}}({{\Delta }_{d}})$ and that found from the nanoscale inhomogeneity [19, 20] in ${\rm B}{{{\rm i}}_{2}}{\rm S}{{{\rm R}}_{2}}{\rm CaC}{{{\rm u}}_{2}}{{{\rm O}}_{8+x}}$ (BSCCO) is heartening, and the nature of ${{P}_{{\rm IMT}}}({{\Delta }_{d}})$ is not found in cuprates. We draw readers attention to the fact that in plotting $P({{\Delta }_{d}})$ in figure 1, we introduced an additional prefactor $J(3\;g_{ij}^{xy}+1)/8$ such that the GIMT and IMT calculations produce the same Δ_d at V = 0 ($g_{ij}^{xy}=1$ for IMT). We also present the GIMT distribution of one sub-dominant Cooper-channel order for the conc disorder, ${{\Delta }_{{{d}_{xy}}}}$, as shown in the right inset of figure 1(c). The tiny width of $P({{\Delta }_{{{d}_{xy}}}})$ and also the smallness of ${\rm max} \{{{\Delta }_{{{d}_{xy}}}}(i)\}$ at all V illustrate the insignificance of ${{\Delta }_{{{d}_{xy}}}}$ compared with Δ_d. Their redundancy is further verified by repeating the GIMT calculations where ${{\Delta }_{{{d}_{xy}}}}$, ${{\Delta }_{{{s}_{xy}}}}$, and ${{\Delta }_{{{s}_{xs}}}}$ are suppressed by hand, without showing any change in the results. The disorder dependences of $P({{\Delta }_{{{s}_{xy}}}})$ and $P({{\Delta }_{{{s}_{xs}}}})$ are essentially the same as that of $P({{\Delta }_{{{d}_{xy}}}})$. Upon establishing their insignificance, we drop all the sub dominant orders for the box disorder calculation in favor of simplicity and clarity.

The qualitative as well as quantitative differences in the GIMT and plain IMT results for $P(\rho )$ and $P({{\Delta }_{d}})$ raise the question: how do the presence or absence of inhomogeneities in the presence of strong correlations affect the disorder dependence of various physical observables? We address this issue in figures 2 and 3 where we show the evolution of ODLRO, superfluid density, and average DOS with V. The results bring out a striking contrast in the V-dependence of these quantities with or without the strong correlations.

Zoom In Zoom Out Reset image size

3.2.1. ODLRO

First, we discuss the superconducting ODLRO defined as:

$$\begin{eqnarray}&&\Delta _{{\rm OP}}^{2}=\mathop{{\rm lim} }\limits_{|i-j|\to \infty }{{F}_{\delta ,\delta ^{\prime} }}(i-j)\end{eqnarray} \tag{ 13 }$$

where ${{F}_{\delta ,\delta ^{\prime} }}(i-j)=\langle B_{i\delta }^{\dagger }{{B}_{j\delta ^{\prime} }}\rangle$ and $B_{i\delta }^{\dagger }=(c_{i\uparrow }^{\dagger }c_{i+\delta \downarrow }^{\dagger }+c_{i+\delta \uparrow }^{\dagger }c_{i\downarrow }^{\dagger })$ is the singlet Cooper-pair creation operator on the links connecting the neighbouring sites at i and $i+\delta$. ${{F}_{\delta ,\delta ^{\prime} }}(i-j)$ is independent of δ and $\delta ^{\prime}$ in the limit of large $|i-j|$. Since ${{F}_{\delta ,\delta ^{\prime} }}(i-j)$ can be interpreted as simultaneous hopping of a singlet Cooper-pair on a link, the GRF corresponding to this process is $g_{i,j}^{t}g_{i+\delta ,j+\delta ^{\prime} }^{t}$. The evolution of ${{\Delta }_{{\rm OP}}}$ (normalized by its value $\Delta _{{\rm OP}}^{(0)}$ at V = 0) with the box disorder is presented in figure 2(a). With increasing V, ${{\Delta }_{{\rm OP}}}$ from the IMT results decreases monotonically, initially with a gradual fall that is followed by a rapid decrease beyond $V\geqslant 1$, causing it to almost vanish by $V\sim 2.5$. Such behavior is consistent with the V-dependence of ${{P}_{{\rm IMT}}}({{\Delta }_{d}})$ in figure 1(d). The strong correlations make the superconductor robust against disorder. There is hardly any significant fall of the GIMT value of ${{\Delta }_{{\rm OP}}}$ at V = 2.5, a strength by which the transition to a non-superconducting state occurs within the IMT scheme. The Gaussian-type shape of ${{P}_{{\rm GIMT}}}({{\Delta }_{d}})$ with minimally depleted centroid for all V is crucial to ascertain that ${{\Delta }_{{\rm OP}}}$ changes very little in the whole range of V within the GIMT calculation. The comparison of the GIMT and IMT results for ${{\Delta }_{{\rm OP}}}$ from the conc disorder, as presented in the inset of figure 2(a), shows a trend similar to the box disorder results. A numerically inexpensive scheme for approximating ${{\Delta }_{{\rm OP}}}$ is further discussed in appendix B.

The length scale of decay of ${{F}_{\delta ,\delta ^{\prime} }}(i-j)$ with $|i-j|$ defines the coherence length (ξ) of a superconductor. The evolution of ${{F}_{\delta ,\delta ^{\prime} }}(i-j)$ for several V, as presented in the lower inset of figure 2(a) from the GIMT method, implies that there is hardly any V-dependence of ξ. Given that ξ is already small (only a few lattice spacings) in the clean cuprates, it is not expected to show any significant reduction with disorder. A similar insensitivity of ξ on V is also found in the IMT results.

3.2.2. Superfluid stiffness

The defining characteristic of a superconductor lies in the Meissner [74] effect, which is quantified by the stiffness of the BCS–GS wave function to an externally applied phase twist. This stiffness translates into its finite superfluid stiffness, which is proportional to the superfluid density. Within the Kubo formalism, the superfluid stiffness, denoted as D_s, is measured by

$$\begin{eqnarray}&&\frac{{{D}_{s}}}{\pi }=\langle -{{k}_{x}}\rangle -{{\Lambda }_{xx}}({{q}_{x}}=0,{{q}_{y}}\to 0,\omega =0),\end{eqnarray} \tag{ 14 }$$

where k_x is the kinetic energy along the x-direction and Λ_xx is the long wavelength limit of transverse (static) current–current correlation function [75]. It is calculated by Fourier transforming the impurity-averaged Matsubara Greenʼs function;

$$\begin{eqnarray}&&{{\Lambda }_{xx}}({\bf q},i{{\omega }_{n}})=\frac{1}{N}\int _{0}^{1/T}d\tau {{e}^{i{{\omega }_{n}}\tau }}\langle j_{x}^{p}({\bf q},\tau )j_{x}^{p}(-{\bf q},0)\rangle ,\end{eqnarray} \tag{ 15 }$$

where $j_{x}^{p}({\bf q})$ is the paramagnetic current and ${{\omega }_{n}}=2\pi nT$ (n is a positive integer). The GRF for $\langle {{k}_{x}}\rangle$ is $g_{i,i+x}^{t}$ and it is $g_{i,i+x}^{t}g_{j,j+x}^{t}$ for ${{\Lambda }_{xx}}({{{\bf r}}_{i}},{{{\bf r}}_{j}},\tau )$. To obtain a good ${\bf q}$-resolution of Λ_xx from a finite simulation box, we used an effectively larger system employing a 'repeated zone scheme' (RZS) (see appendix C) for a single realization of disorder. The final ${{q}_{y}}\to 0$ extrapolation of the disorder averaged Λ_xx obtained from a single unit cell was guided by $\Lambda _{xx}^{{\rm RZS}}({{q}_{y}})$. The disorder dependence of ${{\Delta }_{{\rm OP}}}$ and D_s are quite similar, whereas the IMT causes them to almost vanish by $V\sim 2.5$; there is hardly any appreciable reduction of them in the GIMT results.

The inset of figure 2(b) shows the individual disorder dependence of $\langle -{{k}_{x}}\rangle$ and ${{\Lambda }_{xx}}({{q}_{y}}\to 0)$. The $\langle -{{k}_{x}}\rangle$ from the GIMT method are smaller than the IMT values due to a complete removal of double occupancy. The kinetic energy, which is the diamagnetic contribution to D_s, decreases slowly with disorder and the rate of decrease is similar for both the IMT and GIMT calculations because it is proportional to the total density. For V = 0, Λ_xx is zero in the absence of any quantum fluctuations, making D_s equal to $\langle -{{k}_{x}}\rangle$. With increasing V, Λ_xx grows in the IMT scheme; it almost equals $\langle -{{k}_{x}}\rangle$ by $V\sim 2.5$ and continues to reduce their differences for larger V. However, Λ_xx remains insensitive to V within the GIMT scheme, making D_s robust to impurities.

Such a strong insensitivity of D_s to V is surprising although consistent with the broad experimental trends [22]. There are recent data on the temperature and doping dependences of D_s [76] for the cuprates. We hope similar experiments will be extended to careful disorder dependence for an appropriate comparison with our results.

3.2.3. DOS

We now turn towards the discussion of the disorder-dependence of the site-averaged DOS, given by:

$$\begin{eqnarray}&&N(\omega )=\frac{1}{N}\mathop{\sum }\limits_{i,n}g_{ii}^{t}\left\{ |{{u}_{i,n}}{{|}^{2}}\delta (\omega -{{E}_{n}})+|{{v}_{i,n}}{{|}^{2}}\delta (\omega +{{E}_{n}}) \right\}\end{eqnarray} \tag{ 16 }$$

where E_ns are the BdG eigenvalues. We focus on $N(\omega )$ for ω up to the coherence peak energy. Although $N(\omega )$ consists only of δ-function peaks at T = 0, it is necessary to obtain a smooth trace such that the relevant features are truly identified and free of any artifact of smoothing. To achieve this, we generate a denser spectrum by using RZS (see appendix C). We broaden each of the δ-function in $N(\omega )$ by an amount comparable to the average spacing of the E_nʼs, and the final $N(\omega )$ is presented in figures 3(a) and (b). After establishing that the additional Cooper-channel orders, e.g., ${{\Delta }_{{{s}_{xs}}}}$, ${{\Delta }_{{{s}_{xy}}}}$, ${{\Delta }_{{{d}_{xy}}}}$ make insignificant contributions to $N(\omega )$, just as in all the other observables, we restrict the DOS calculations implemented with RZS to only the $t-t^{\prime} -J$ model for all classes of disorder.

Zoom In Zoom Out Reset image size

Although the average $N(\omega )$ calculated within the plain IMT shows significant gap filling in figure 3(b) for large V, results from the GIMT scheme in figure 3(a) show that the low-lying DOS remains unaffected with impurities even up to V = 3. Strong correlations thus appear to prohibit the formation of the low-energy states in the GIMT results. Although the coherence peak height is reduced with V in both methods, the rate of this reduction is significantly weaker in the GIMT. Our results for the conc disorder are similar to [29]. The strong correlations thus seem to protect the V-shape of the $N(|\omega |\leqslant \Delta _{d}^{(0)})$, a feature established in the STS-experiments [21, 27].

To develop a deeper understanding of the robustness of low-energy DOS in GIMT, we calculated the LDOS for two regions with large and small Δ_d. We average the LDOS for these regions with a weak ${{\Delta }_{d}}(i)\in [0.18,0.21]$ and with a strong ${{\Delta }_{d}}(i)\in [0.25,0.28]$ in figure 3(c) and find that LDOS in the two regions behave quite differently even for $N(|\omega |\leqslant \Delta _{d}^{(0)})$. We also show in section 3.3 that the regions of large (small) ${{\Delta }_{d}}(i)$ corresponds also to the regions of large (small) ρ_i within the GIMT calculations so that the inhomogeneity in the local density is also reflected in the same LDOS. Although the anti-correlation between the coherence peak heights and the local gap values is clearly identified as in the experiments noted in [21], we did not find the expected homogeneity in LDOS for $|\omega |\lt \Delta _{d}^{(0)}$. Such homogeneity is recovered in the site-averaged $N(\omega )$, which can be easily seen from the $P({{\Delta }_{d}})$ in figure 1(c). The peak remains symmetric for all V centered around $\Delta _{d}^{(0)}$. For the conc disorder, the homogeneity of the low-lying LDOS is maintained fairly well [29] up to $\omega \sim \Delta _{d}^{(0)}/2$. This is shown in the inset presenting the LDOS for sites with ${{\Delta }_{d}}(i)\in [0.14,0.17]$ (blue curve) and for sites with ${{\Delta }_{d}}(i)\in [0.19,0.22]$ (red curve).

3.2.4. Fourier transformed local density of states (FT-LDOS)

Another occurrence that has recently received attention is the energy-resolved FT-LDOS, defined by $N({\bf q},\omega )={{\sum }_{i}}{{e}^{-i{\bf q}.{{{\bf r}}_{i}}}}{{N}_{i}}(\omega )$, where ${{N}_{i}}(\omega )$ is the LDOS at site i, obtained from equation (16) without performing the summation over i. FT-LDOS is extracted from the STS data indicating the existence of a periodic and energy-resolved modulation of the local density. An interpretation of the STS data comes from the simple 'octet'-model [77, 78], based on the scattering of the d-wave quasiparticles from a single impurity. The resulting theory yields an octet of wave-vectors ${\bf q}$ connecting the constant quasiparticle-energy contours for which the scattering probability is maximal. The Fermi surface obtained from the dispersion of these octet members agrees well with its direct measurement from the angle-resolved photo-emission spectroscopy [79]. The extension of the octet analysis to many impurities is a subtle issue. Without invoking strong correlations, results from many impurity calculations are sensitively dependent on the details of the disorder [80] for a reasonable fit to the STS data, and the octet peaks are often masked by a large noise [81]. Further, non-dispersive FT-LDOS peaks have been attributed to competing orders [38] for under-doped samples. However, to the best of our knowledge, the role of strong correlations on the FT-LDOS has never been addressed in a calculation.

We present the FT-LDOS from the GIMT calculations in figure 4. Each density plot shows the power spectrum $P({\bf q},\omega )=\overline{|N({\bf q},\omega ){{|}^{2}}/N}$ for a given ω on the first Brillouin zone (BZ): ${{q}_{x}},{{q}_{y}}\in [-\pi ,\pi ]$, disorder averaged over 10 independent realizations of the conc disorder with ${{n}_{{\rm imp}}}=0.09$. We apply RZS for smooth data from a larger system (see appendix C). Octet peaks can easily be recognized, such as ${{{\bf q}}_{1}}\sim (0,\pm \pi /{{\mathcal{L}}_{1}})$ or $(\pm \pi /{{\mathcal{L}}_{1}},0)$. The dispersion of these peaks results in a continuous change in ${{\mathcal{L}}_{1}}\sim 1.36$ for $\omega =0.05$ to ${{\mathcal{L}}_{1}}\sim 7.48$ for $\omega =0.2$. Another (dispersing) peak, likely the ${{{\bf q}}_{7}}$ in the octet terminology [77], can also be discerned along the diagonal direction of the BZ (better resolved for negative ω). Our results replicate many of the experimental trends, e.g., the shape and location of the FT-LDOS peaks and their dispersion (except perhaps for ${{{\bf q}}_{5}}$), among other things. However, the asymmetry of the FT-LDOS profiles for $\pm \omega$ is found to be much stronger than indicated by the experimental data. The qualitative features of figure 4 remain unaltered with a change of ${{n}_{{\rm imp}}}$ to 0.05 from 0.09. A similar calculation of FT-LDOS within the IMT scheme (shown in the supplementary material) was found to capture the key points of figure 4; however, appropriate choice of the model parameters for the IMT calculations is crucial for an appropriate comparison.

Zoom In Zoom Out Reset image size

Before closing the discussion on the disorder dependence of the observables, we note that the degree of inhomogeneity in Δ_d reduces upon including strong electronic repulsions (see figure 1(c)), at least for the moderate disorder ($0.5\leqslant V\leqslant 2.0$). However, the behaviour of the observables remains very different compared with that of the self-consistent T-matrix predictions [15] based on the assumption of homogeneous but renormalized Δ_d for any V. The inadequacy of self-consistent T-matrix approximations for describing disordered cuprates is discussed in [33], which is a weak coupling theory built on a simple metallic normal state. The strong interactions amplify such differences further due to modifications of the normal states, as we see in section 3.5.

Having seen the dramatic robustness of the observables to impurities in the GIMT scheme, it appears natural to look for the relationship between such insensitivity and the inhomogeneity compared with the plain IMT. We develop an important insight by studying the spatial distribution of the different order parameters on the lattice, as explained in the following discussion.

We present the GIMT and IMT results for ${{\Delta }_{d}}(i)$ on the left and right columns of figure 5(a), respectively. Similarly, the local density, ρ_i, from the GIMT and IMT calculations is shown in the left and right columns of figure 5(b). All results are shown for a particular realization of the box disorder. There are important points to note here: first, the spatial inhomogeneity in density from the GIMT calculation does not evolve with V (not even in its magnitude), as seen from figures 5(b1), (b3) and (b5). Although this observation reinforces the finding of figure 1(a), the additional information here is on the spatial structures, as explained in the following discussion. The profiles of ρ_i in plain IMT, shown in figures 5(b2), (b4), and (b6), are similar to those from the GIMT in figures 5(b1), (b3), and (b5) respectively, although its magnitude differs significantly (See also figure 1(b)). Second, the evolution of the spatial profile in ${{\Delta }_{d}}(i)$ from GIMT with V is also very weak as shown in figure 5(a1), (a3), and (a5), compared with the same spatial profile calculated in the plain IMT scheme as given in figures 5(a2), (a4) and (a6). Third, an impressive spatial correlation is found between the GIMT profiles of ρ_i and ${{\Delta }_{d}}(i)$ for all V, which implies that both are large (or small) approximately on the similar region of space (further illustrated in figure 7(b)), whereas such spatial correlations are completely absent in the IMT results. The IMT data in figure 5(a6) illustrates the formation of 'islands' with large ${{\Delta }_{d}}(i)$, separated by regions where ${{\Delta }_{d}}(i)\approx 0$. We verify that the sites with large Δ_d are indeed the sites where $|{{V}_{i}}-\mu _{i}^{{\rm HS}}|\approx 0$ for large V and correspond to the spatial regions, allowing maximum particle-hole mixing in the plain IMT.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Clearly, the origin of the inhomogeneities leading to the island formation in ${{\Delta }_{d}}(i)$ within the IMT scheme is essentially the same as in the case of an sSC [59, 72]. The lack of the evolution of the GIMT spatial structures in ${{\Delta }_{d}}(i)$, on the other hand, illustrates that there are no well defined 'islands'. Primarily, those sites lying in the deepest valleys of the disorder potential become nearly doubly occupied in the plain IMT. As a result, no significant spatial correlation between ${{\Delta }_{d}}(i)$ and ρ_i is expected at large V, which is clear from our data. Strong electronic correlations in the GIMT, however, ensure that the local density never goes beyond half-filling anywhere in the system which leads to the observation that the sites with large ρ_i (close to half filling) also satisfy $|{{V}_{i}}-\mu _{i}^{{\rm HS}}|\approx 0$ approximately (at large V) for $\rho =0.8$. It is thus not surprising that we see a strong spatial correlation between local density and pairing amplitude in the GIMT results in the left columns of figure 5(a) and 5(b) (see also figure 7(b1) and (b2)), quite in contrast with the IMT findings. This is further illustrated in the supplementary material for out-of-plane disorders where the contrast is even more prominent.

This observation raises a conceptual question: Inhomogeneities in ${{\Delta }_{d}}(i)$ within the IMT calculations arise with a length-scale of ξ, whereas the fluctuations in ρ follow (uncorrelated) disorder with the natural length-scale $k_{F}^{-1}$, the Fermi wavelength. This results in a clear separation of the two length-scales. But such a distinction is difficult because the spatial profile of one follows the other in the GIMT results. This situation is actually not so crucial for the cuprates, in which both these scales are of the order of a few lattice spacings. In fact, a recent extraction [82] of the length scales associated with disorder seems to validate the GIMT picture.

The mean-field Hamiltonian of equation (9) can be cast in the following effective form:

$$\begin{eqnarray}&&{{\mathcal{H}}_{{\rm MF}}}=\mathop{\sum }\limits_{i,\delta ,\sigma }{{t}_{{\rm eff}}}(i,\delta )\;c_{i\sigma }^{\dagger }{{c}_{i+\delta \sigma }}+\mathop{\sum }\limits_{i,\sigma }{{V}_{{\rm eff}}}(i){{n}_{i\sigma }}+\left\{ {\rm pairing}\;{\rm terms}\;{\rm involving}\;\;\Delta _{i}^{\delta } \right\}\end{eqnarray} \tag{ 17 }$$

where the explicit expressions for ${{t}_{{\rm eff}}}$ and ${{V}_{{\rm eff}}}$ can be extracted to compare equations (9) and (17), both of which become inhomogeneous in the presence of the disorder. However, inhomogeneities in ${{t}_{{\rm eff}}}$ from the IMT calculation are very weak and that in ${{V}_{{\rm eff}}}$ is largely similar to the bare disorder. Obviously, the physics of the strong correlations within GIMT is also contained in equation (17) through the local GRFs and their derivatives with respect to ρ_i. The physical origin for the differences between the IMT and GIMT results can be explained by using the following argument. A strong repulsive interaction smears out any local accumulation of the charge carriers. Such a smearing plays an important role in the GIMT but is missing from plain IMT scheme. As a result, the spatial distribution of electrons ought to be different between the two calculations. Such reorganization of the electronic density must modify all other local-order parameters in a self-consistent manner and hence should also alter all the physical observables. The question remains: how is this correlation, driven physics incorporated in the ${{\mathcal{H}}_{{\rm MF}}}$? A significant insight is obtained by studying the spatial fluctuations of ${{t}_{{\rm eff}}}$ and ${{V}_{{\rm eff}}}$ individually, as presented in figure 6. We plot ${{t}_{{\rm eff}}}(i)={{\sum }_{\delta }}{{t}_{{\rm eff}}}(i,\delta )$ and V_eff(i), both from the IMT and GMT calculations, on the lattice. Both the t_eff(i) and the V_eff(i) show qualitatively similar spatial variations in the plain IMT (figures 6(a) and (b)). In contrast, a strong spatial anti-correlation develops in the fluctuating parts of ${{t}_{{\rm eff}}}$ and ${{V}_{{\rm eff}}}$ within the GIMT scheme as demonstrated in figures 7(b3) and (b4) using scatter plots. Comparing their spatial profiles in figures 6(c) and (d) shows that t_eff(i) is small precisely in those regions where V_eff(i) is large and vice versa. Such a competing behavior [83] is a result of the strong repulsive interactions; we elaborate in the following discussion.

Zoom In Zoom Out Reset image size

According to equation (17), t_eff(i) is the probability that an electron will hop onto the site i from the neighbours. A site with ${{V}_{{\rm eff}}}(i)\gt {{V}_{{\rm eff}}}(i+\delta )$ will support ${{\rho }_{i}}\lt {{\rho }_{i+\delta }}$ leading to a larger $g_{i,i+\delta }^{t}$ (compared with when ${{V}_{{\rm eff}}}$ were the same for i and $i+\delta$). This, in turn, enhances t_eff(i), increasing the final ρ_i compared with what V_eff(i) alone would have predicted in the first place. So the high 'hills' of the disorder profile do not remain as sparsely populated by electrons in the GIMT as expected in the IMT calculation. Exactly the opposite mechanism smears out the charge accumulation in the GIMT from sites with deep 'valleys' of disorder potential.

In contrast, t_eff(i) is largely independent of ρ_i within the plain IMT framework, although it develops a very weak density dependence through the Fock shifts. In fact, t_eff(i) becomes weaker in regions with high hills and deep valleys of the disorder profile within the plain IMT. Such sites tend to get decoupled from hopping in the absence of the 'feedback-loop' in terms of $g_{i,i+\delta }^{t}$, quite in contrast to the GIMT scenario discussed previously. This picture is substantiated through our results in figures 6(e) and (f), where we demonstrate the competing effects of t_eff(i) and V_eff(i) on ρ_i, if the two were to work individually. It is thus no surprise that the GIMT local density (left column, figure 5(b), incorporating both ${{t}_{{\rm eff}}}$ and ${{V}_{{\rm eff}}}$) becomes less inhomogeneous than the IMT finding.

Finally, we note that one of the important effects of the strong correlation lies in the significant reduction [58] of the magnitude of V_eff(i) from the bare V_i. We demonstrate this by presenting the distribution of $P({{V}_{{\rm eff}}})$ in figures 8(a) and (b) both from the GIMT and plain IMT calculations. The homogeneous component of ${{V}_{{\rm eff}}}$ (the contribution at V = 0) from the Hartree shift is already subtracted in this analysis. We see that whereas in the IMT calculations the range of V_eff(i) actually increases a bit from the bare value of $\pm V/2$, strong electronic correlations indeed cause its drastic reduction as seen from the resulting ${{P}_{{\rm GIMT}}}({{V}_{{\rm eff}}})$. Interestingly, ${{P}_{{\rm GIMT}}}({{V}_{{\rm eff}}})$ becomes more asymmetric with increasing V with larger weight accumulating for negative V_eff(i). In contrast, the ${{P}_{{\rm IMT}}}({{V}_{{\rm eff}}})$ maintains its symmetry around zero just as in the bare $P({{V}_{i}})$ by construction. Nevertheless, we found that the asymmetric ${{P}_{{\rm GIMT}}}({{V}_{{\rm eff}}})$ still produces $\langle {{V}_{{\rm eff}}}\rangle =0$ for all V. In fact, it is this asymmetry that translates into an oppositely asymmetric ${{P}_{{\rm GIMT}}}({{\rho }_{i}})$ in figure 1(a). We explain this intriguing feature in the following way: For a non-interacting problem at large V, the local density ranges from zero to double occupancy on the highest hill and the deepest valley, respectively. The GA, however, restricts ${{\rho }_{i}}\leqslant 1$, breaking the symmetry of the non-interacting problem. A self-consistent interplay of the strong correlations and the disorder now ensures that the effective disorder has more 'attractive' points such that the maximal ρ_i is limited by half-filling everywhere, still maintaining the global density at the desired value ($\rho =0.8$, in our case). An evolution of spatial correlation in V_eff(i) with disorder, producing larger patches of attractive points, is further illustrated by using colour-density plots in figure 7(a), indicating a growing spatial correlation of such attractive points. Let us emphasize that we verified the crucial physical insights, emerging from the study of the spatial distribution of the local observables, as expressed in all the permutations of disorder we considered.

Zoom In Zoom Out Reset image size

The intriguing evolution of the spatially correlated effective disorder in the GIMT calculations immediately tells us that they must also modify the properties of the underlying normal state (NS). Here, the NS is defined as the solution of the same ${{\mathcal{H}}_{{\rm MF}}}$ of equation (17) but suppressing all the Cooper-channel order parameters $\Delta _{i}^{\delta }$, for all i and δ. We caution readers that our usage of the term 'normal state' might lead to confusion because it is frequently used to identify the $T\gt {{T}_{c}}$ phase of a HTSC, where the d-wave order is already destroyed by thermal fluctuations. Our NS, on the other hand, is a true GS at T = 0, had there been no pairing. Because the locally renormalized GRF leads to spatially correlated disorder, the ${\rm N}{{{\rm S}}_{{\rm GIMT}}}$ must differ from the ${\rm N}{{{\rm S}}_{{\rm IMT}}}$, where the later is essentially the solution of the Anderson model [84] of non-interacting electrons in a disordered background. The ${\rm N}{{{\rm S}}_{{\rm GIMT}}}$ plays a crucial role in explaining our GIMT results, which suggests a physical picture for the dirty cuprates, as we explain in the following discussion.

The robustness of the observables to the impurities, as demonstrated in section 3.2, is actually reminiscent of Andersonʼs theorem [30]. Although the theorem does not apply to a weak-coupling d-wave superconductor [14], our results open up such a possibility for strongly correlated cuprates. To support this scenario, we first argue that our GIMT results are consistent with Andersonʼs mechanism provided we extend the scope of pairing between all the 'correlated' single-particle normal states (${\rm N}{{{\rm S}}_{{\rm GIMT}}}$). To establish this, we solve for $\mathcal{H}_{{\rm MF}}^{{\rm Normal}}$ (equation (17) with all $\Delta _{i}^{\delta }$ terms suppressed) to obtain the ${\rm N}{{{\rm S}}_{{\rm GIMT}}}$ parameters t_eff(i) and V_eff(i) in a self-consistent manner in the first step. These parameters carry information on the local orders in Hartree- and Fock-channels. In the next step, we set up a BdG-type matrix ${{\tilde{\mathcal{H}}}_{{\rm BdG}}}$ in terms of these ${\rm N}{{{\rm S}}_{{\rm GIMT}}}$ parameters in the diagonal blocks, and with $\Delta _{i}^{\delta }$ in the off-diagonal blocks, where $\Delta _{i}^{\delta }$ are obtained from the corresponding GIMT solution of the full ${{\mathcal{H}}_{{\rm MF}}}$, including the dSC order. This ${{\tilde{\mathcal{H}}}_{{\rm BdG}}}$ is then diagonalized only once in the final step, and we recalculate the new (non-self-consistent) local pairing amplitudes $\tilde{\Delta }_{\delta }^{i}$ from this sole diagonalization. We find impressive agreement between $\Delta _{\delta }^{i}$ and $\tilde{\Delta }_{\delta }^{i}$ for all i and δ. Such a non self-consistent calculation is equivalent to an Anderson-type idea (though the pairing in this case is not limited to the time-reversed states alone). This interpretation works because the effective attraction does not seem to renormalize the spatially fluctuating NS orders in the diagonal blocks of ${{\tilde{\mathcal{H}}}_{{\rm BdG}}}$. If instead we had any significant mismatch between $\Delta _{\delta }^{i}$ and $\tilde{\Delta }_{\delta }^{i}$, self-consistency would have required a feedback loop, which could, in general, update NS order parameters as well, implying that the pairing modifies single-particle states. Such a conclusion goes against Andersonʼs philosophy, according to which pairing from effective attraction is independent of disorder. This is because the impurities are already consumed in producing the corresponding 'exact eigenstates'. Note that in the case of the sSC, the normal states are essentially the exact eigenstates of the effective single-particle Hamiltonian because there are no strong correlations.

It is not only the pairing between the normal states which is responsible for the insensitivity of superconductivity to the impurities. The extent of localization of these states is also crucial for the inferences made by Andersonʼs theorem. Although such a pairing theory for a sSC can be formulated in principle for any strength of disorder by tailoring the disorder-induced inhomogeneities it contains, the robustness of s-wave superconductivity breaks down at large disorders ($V\geqslant t$) due to strong localization [59]. It is interesting to note that the insensitivity of observables for dirty cuprates actually extends even beyond $V\gt 2.5$ (see figures 2 and 3), a strength by which the traditional Andersonʼs theorem (for sSC) weakens drastically. To settle this apparent conflict, we study the differences between the NS from the GIMT and IMT schemes, focusing on their localization properties. We quantify the extent of localization by calculating the participation ratio (PR) for the effective one-particle NS near μ, defined as

$$\begin{eqnarray}&&{{(PR)}_{n}}=\mathop{\sum }\limits_{i}{{\left| {{\psi }_{n}}(i) \right|}^{4}}\end{eqnarray} \tag{ 18 }$$

where ${{\psi }_{n}}(i)$ is the normalized eigenfunction of the NS for the eigenvalue ${{\varepsilon }_{n}}\approx \mu$. For better statistics, we actually average PR over all those nʼs such that $|{{\varepsilon }_{n}}-\mu |\leqslant 0.03$. The PR measures the fraction of sites of the system contributing to the nth wave function and is a good estimate of the localization area in two dimensions. The behavior of PR with the disorder strength V in figure 8(c) shows that the one-particle ${\rm N}{{{\rm S}}_{{\rm GIMT}}}$ (with ${{\varepsilon }_{n}}\approx \mu$) remains fairly extended up to V = 1.75 on a scale of the system size. In contrast, the ${\rm N}{{{\rm S}}_{{\rm IMT}}}$, pairing of which leads to Andersonʼs theorem for sSC, localizes relatively fast. The normal state DOS, ${{N}_{{\rm NS}}}(\omega \approx \mu )$, that incorporates the GA was found to be insensitive to V (although it shows strong sample-to-sample fluctuations) [86], consistent with the extended nature of the ${\rm N}{{{\rm S}}_{{\rm GIMT}}}$. The robustness of observables for dirty superconductors follows as a consequence of Andersonʼs traditional pairing philosophy, as long as the exact eigenstates (or NS in our terminology) remain fairly extended in the scale of the system size. A similar picture emerges for the dirty cuprates from the preceding discussions: (a) the GIMT findings can be thought of as a result of pairing between correlated normal states. (b) Such normal states remain delocalized up to a large disorder.

An important consequence of the conventional Andersonʼs theorem [30] is the irrelevance of any detailed spatial information about the system. Such irrelevance for our case is seen by calculating $\tilde{N}(\omega )$ as follows [85]:

$$\begin{eqnarray}&&\tilde{N}(\omega )=\int d{{\Delta }_{d}}P({{\Delta }_{d}}){{N}_{{\rm BCS}}}(\omega ,{{\Delta }_{d}})\end{eqnarray} \tag{ 19 }$$

where ${{N}_{{\rm BCS}}}$ is the density of states of a homogeneous d-wave BCS superconductor and is given by [13]

$$\begin{eqnarray}\begin{array}{rcl} {{N}_{{\rm BCS}}}(\omega ,{{\Delta }_{d}}) & = & \frac{2}{\pi }\frac{|\omega |}{{{\Delta }_{d}}}\mathcal{K}\left( \frac{|\omega |}{{{\Delta }_{d}}} \right)\;{\rm for}\;\frac{|\omega |}{{{\Delta }_{d}}}\leqslant 1, \\ {} & = & \frac{2}{\pi }\mathcal{K}\left( \frac{{{\Delta }_{d}}}{|\omega |} \right)\;{\rm for}\;\frac{|\omega |}{{{\Delta }_{d}}}\gt 1. \\ \end{array}\end{eqnarray} \tag{ 20 }$$

Here, $\mathcal{K}(x)$ is the complete elliptic integral and Δ_d is the d-wave gap at T = 0. Such a definition keeps track only of the average fluctuations in the pairing amplitude; no information about the spatial structures of inhomogeneities is accounted for. The resulting $\tilde{N}(\omega )$, using $P({{\Delta }_{d}})$ from our GIMT calculation, is presented in figure 8(d). We show $\tilde{N}(\omega )$ only for the positive ω because it is manifestly symmetric about $\omega =0$. The similarity of $\tilde{N}(\omega )$ with ${{N}_{{\rm GIMT}}}(\omega )$ in figure 3(a) illustrates the redundancy of any spatial correlation of inhomogeneities in ${{N}_{{\rm GIMT}}}(\omega )$.

Before closing the current discussion, we must mention that the preceding arguments supporting the robustness of cuprate superconductors as being similar to the conventional Andersonʼs theorem are all indirect arguments, although they add up to a significant claim. It will be interesting to substantiate these claims also by direct demonstration. Such a formalism is beyond the scope of this study and will be addressed elsewhere.

It is important to address the fate of the cuprates as the disorder strength is increased beyond what is reported here. It is easy to see that the GIMT scheme, which removes all double-occupancy by enforcing ${{\rho }_{i}}\leqslant 1$ (for all i), would not work on the sites with ${{V}_{i}}\geqslant U$. However, this is not a deficit of the GA. It would require the modification of allowing occasional ${{\rho }_{i}}\gt 1$ based on the relative strength of U and V_i. Such calculations, however, are beyond the scope of our study and will be addressed elsewhere. Nevertheless, we found that the GIMT self-consistency falters technically even for $V\gt 3.0$ (box disorder). In fact, making it work for large V requires significant improvements in the iterative self-consistency in terms of applying accelerated convergence methods. Our usage of Broyden [87] and modified Broyden [88] methods was crucial in pushing the self-consistency to larger disorder strengths.

It turns out that conc disorder helps to push the self-consistent calculations towards a higher V₀, because we can tune ${{n}_{{\rm imp}}}$ independently to keep iterative self-consistency under control. The average DOS calculated with ${{V}_{0}}=5$ and ${{n}_{{\rm imp}}}=0.03$, 0.05 are presented in figure 9(a) which shows the first evidence of a gap-filling within the GIMT scheme [89]. Figures 9(b) and (c) present the LDOS from the GIMT calculation around a single impurity, where the strength of the sole impurity is ${{V}_{0}}=7$ and ${{V}_{0}}=10$, an energy scale larger than t, $t^{\prime}$, J, and $J^{\prime}$. Results for both ${{V}_{0}}=7$ and ${{V}_{0}}=10$ show that growing V₀ produces larger $N({\bf r},\omega \simeq 0)$, where ${\bf r}$ represents the nearest neighbours of the impurity site. A simple picture emerges from repeating the calculations for several combinations of V₀ and ${{n}_{{\rm imp}}}$: the gap filling begins when V₀ grows beyond the bandwidth of the homogeneous, Gutzwiller-renormalized normal state: i.e., ${{V}_{0}}\gt 8\;{{g}^{t}}t$, where ${{g}^{t}}=(1-\rho )/(1-\rho /2)$. This observation motivates the formalism of an alternate IMT scheme elaborated further in the supplementary material.

Zoom In Zoom Out Reset image size