Disorder-robust phase crystal in high-temperature superconductors stabilized by strong correlations

The simultaneous interplay of strong electron-electron correlations, topological zero-energy states, and disorder is yet an unexplored territory but of immense interest due to their inevitable presence in many materials. Copper oxide high-temperature superconductors (cuprates) with pair breaking edges host a flat band of topological zero-energy states, making them an ideal playground where strong correlations, topology, and disorder are strongly intertwined. Here we show that this interplay in cuprates generates a new phase of matter: a fully gapped ``phase crystal"state that breaks both translational and time-reversal invariance, characterized by a modulation of the $d$-wave superconducting phase co-existing with a modulating extended $s$-wave superconducting order. In contrast to conventional wisdom, we find that this phase crystal state is remarkably robust to omnipresent disorder, but only in the presence of strong correlations, thus giving a clear route to its experimental realization.

Different phases of matter appear in condensed matter physics due to the presence of strong correlations, 1,2 disorder, 3,4 or topology. 5,6 While they have all been intensively studied individually and also to some extent pairwise, [7][8][9] the combination of all three effects is one of the currently most challenging problem in physics. Cuprates with pair breaking edges provide an interesting platform where this three-way interplay can become manifest. First, the d-wave pairing symmetry in cuprates leads to Andreev bound states at the pair breaking [110] edges forming flat bands of zeroenergy states, 10,11 protected by the bulk topology. 12,13 Secondly, electron correlations are exceptionally strong in the cuprates, with the undoped parent compound even being a Mott insulator. 1 Finally, disorder, both intrinsic or extrinsic, is prominent in all cuprates. 14,15 It has already been shown that the flat band of zero-energy states is thermodynamically unstable due to the extensive ground-state degeneracy and is hence highly susceptible to electronic correlations. 16 However, the question of how these zero-energy states respond to a simultaneous presence of strong correlations and disorder is still an unsolved problem.
The thermodynamic instability of the zero-energy states opens for a possible appearance of competing phases, breaking the topological protection. [16][17][18][19][20][21] The likely most promising intrinsic scenario involves an exotic state called "phase crystal", which spontaneously breaks both time-reversal and translational symmetries along the [110] edge. [22][23][24][25] In this state, the phase of the d-wave superconducting order parameter is modulated along the edge, creating superflow patterns which lower the free energy by Doppler shifting some, but not all, states away from zero energy. This state was originally found in quasiclassical calculations, 22 but subsequently also in a weak-coupling tight-binding model. 25 However, these and most other studies of competing phases, do not take into account two of the most essential ingredients of cuprates: strong correlations and disorder. It is both unknown if the phase crystal even survives in the presence of strong correlations, and, most importantly, disorder is known to essentially eliminate the ground state degeneracy, 26,27 thus making any competing phase seemingly unlikely when real-world disorder is present.
Here we perform fully self-consistent calculations of the superconducting state in cuprate superconductors and find that strong correlations stabilize a phase crystal state along the [110] edge that develops a full energy gap due to a co-existing s-wave superconducting order. We first show how strong correlations increase the number of zero-energy states for any uniform superconducting phase, contradicting simple topological arguments. Then, when allowing for a non-uniform solution, we find a cascade of phase transitions occurring at different temperatures: d-wave superconductivity occurs below a transition temperature T c , the phase crystal appears further below at temperature T * ∼ 0.2T c , and finally an additional extended s-wave order, with the same spatial modulations as the phase crystal, is generated below T s , generating a full energy gap. Taken together, these phase transitions explain a set of so far seemingly contradictory experimental results. [28][29][30][31][32] Furthermore, we find that the phase crystal state is unexpectedly very robust to disorder, but notably only in the presence of strong correlations. Our results show that the combined effects of strong correlations and topology lead to the emergence of novel phases of matter that survives strong disorder in a highly non-intuitive manner. Model and approach. The d-wave superconducting state in the strongly correlated cuprates has long been assumed to be described by the strong coupling repulsive Hubbard-U model 33 and the equivalent t-J model: 34 Here c † iσ (c iσ ) is the creation (annihilation) operator of an electron with spin σ at lattice site i in a two-dimensional square lattice, S i and n i are the spin and electron density operators, respectively, µ is the chemical potential fixing the average electron density, ij denotes nearest neighbor bonds, J = 4t 2 /U is the super-exchange interaction with U being the onsite Hubbard repulsion strength, and t ij is the hopping amplitude for an electron from site i to j, where we consider arXiv:2103.12756v2 [cond-mat.supr-con] 10 Jan 2023 t ij = −t for nearest neighbor bonds and t ij = t for nextnearest neighbor bonds. For simplicity we express energies in units of t, lengths in units of the lattice spacing a, and set = k B = 1. Furthermore, P is the projection operator which prohibits double occupancies on each lattice site due to the strong onsite repulsive U . The effects of the projection operator can be implemented by the Gutzwiller approximation, where the Hamiltonian parameters are renormalized by statistical factors ensuring no double occupancy. 35 The Gutzwiller approximation has been verified to agree well with variational Monte Carlo calculations, which treat the effects of the projection exactly, 36 in both homogeneous systems 37,38 and for bulk impurities. 39 In the presence of inhomogeneities, such as edges or disorder, the Gutzwiller approximation should be implemented with renormalization factors determined by the local electron density. 40 In this way, the Gutzwiller approximation becomes a powerful theoretical tool encapsulating the simultaneous effects of strong correlations and inhomo-geneities.
We perform a Gutzwiller inhomogeneous mean-field theory (GIMT) treatment of the Hamiltonian in Eq. (1). 40,41 The super-exchange J term gives rise to superconductivity with spin-singlet Cooper pairs living on nearest neighbor bonds, which can be represented by the local dand extended s-wave pairing amplitudes: 40 Here ∆ δ i is the pairing amplitude on the nearest neighbor bond in direction δ. We are further interested in the phase θ of the d-wave pairing amplitude characterized through ∆ d (i) = |∆ d (i)| exp(iθ). To quantify the effects of strong correlations, we compare our GIMT results with the results of standard weak-coupling Bogoliubov-de Gennes (BdG) calculations, 42 where the important projection P in Eq. (1) is ignored and the only effect of the correlations is to create superconducting pairing. Finally, we introduce generic non-magnetic disorder by studying where V i is a site-dependent non-magnetic impurity potential drawn from a random distribution, such that V i ∈ [−V /2, V /2] uniformly, also known as Anderson disorder. Details of the GIMT and BdG methods, system geometry, and choice of the parameters are given in the Methods section. Clean superconductor edge. We begin by looking at a clean superconductor with a pair breaking [110] edge. In Fig. 1(a,b) we show the phase θ of the d-wave pairing amplitude of the ground state at a low temperature, obtained with strong correlations within GIMT and as a comparison without strong correlations within BdG. In both cases θ acquires distinct and modulating non-zero values near the pair breaking [110] edge. These non-zero values of the phase of the d-wave pairing amplitude near the edge imply that this phase cannot be gauged away and, consequently, time-reversal symmetry is broken. 43 Additionally, the modulating nature of this phase shows that the translational symmetry along the [110] edge is also broken. These modulations define a phase crystal state, 24 previously found in the absence of strong correlations in both quasi-classical theory 22,44 and BdG. 25 Clearly, the phase crystal also thrives in the presence of the strong correlations found in the cuprates, and remarkably here also features multiple new distinct properties, as we report below. Before proceeding we also note that the modulating phase of the superconducting order additionally results in circulating orbital currents which we demonstrate in the Supplementary Material (SM).
To start characterizing these new properties, we investigate in Fig. 1(c,d) the behavior of the relative change of the magnitude of the d-wave pairing amplitude, δ |∆ d |, and its phase, θ, along x ⊥ , i.e. perpendicular to the pair breaking edge. As seen, the length scale λ ⊥ , at which the phase θ decays to its zero bulk value, is the same as the healing length of |∆ d |, with and without strong correlations. But notably, λ ⊥ is much shorter when we appropriately include the strong correlations. The reason for this dramatically shorter λ ⊥ is two-fold: First, strong correlations renormalize the hopping amplitudes t ij and the super-exchange interaction J through Gutzwiller factors in the bulk. Hence the bulk superconducting coherence length is reduced, which automatically leads to a decreasing healing length of |∆ d |. We however note that these factors are not significantly different between the bulk and the edge. Secondly, strong electronic repulsion suppresses the charge fluctuations formed at the edge due an increased proximity to the Mott insulating normal state, as previously also established for local impurities. 45 This suppression of the charge fluctuation is mimicked by the local Hartree shift (see Methods), which is notably different on the edge compared to the bulk. As a result, the charge density heals to its bulk value over a short distance and also the healing length of |∆ d | is closely tied to the healing length of the charge density in a strongly correlated state. 41,45 In fact, even near the edge we find that |∆ d | is closely following the relative increase δρ of the charge density, as seen in Fig. 1(c). To summarize, both a short bulk superconducting coherence length and a short charge density healing length result into a short healing length of |∆ d | into the bulk, and consequently a short λ ⊥ . In contrast, λ ⊥ in BdG 25 (and also quasiclassical results 24 ) is only related to the bulk superconducting coherence length, with no relation to the charge density fluctuations, as also clearly seen in Fig. 1(d). The strikingly short λ ⊥ with strong correlations has important consequences for the phase crystal as its modulations along the edge is set by an intricate energy balance: A spatially varying θ in the bulk costs kinetic energy, whereas such variations in θ along the edge instead lower the free energy through a Doppler shift of the zero-energy states. 23,24 Hence, the very short λ ⊥ in GIMT automatically gives a very short modulation wavelength of the phase crystal along the edge, as is evident in Fig. 1(a).
Another striking feature of strong correlations is the appearance of an extended s-wave pairing amplitude |∆ s | at the edge with a modulation length scale set by that of the phase crystal, as seen in Fig. 1(e), while the s-wave amplitude is always negligible in the bulk. Although not shown in Fig. 1(e), the phase of the s-wave pairing amplitude is nearly constant in all regions where its magnitude is finite, with a value of π/2 relative to the phase of the d-wave amplitude. We attribute the existence of this s-wave component to the very short modulation wavelength of the phase crystal as it gives a large number of phase nodes. These nodes lead to some zero-energy states still persisting in the phase crystal, although their degeneracy is reduced. At low enough temperatures, these remaining zeroenergy states make the system thermodynamically unstable, with the consequence that the extended s-wave state is developed, gapping out the last states and generating a full energy gap. Ignoring strong correlations, the phase crystal modulation wavelength is instead longer, no extended s-wave pairing amplitude emerges, and consequently the energy spectrum is never fully gapped. 22,25 Here we note that the time-reversal symmetry breaking due to the formation of the phase crystal state can technically also generate subdominant odd parity pairing amplitudes, as in superconductor-magnet hybrid systems. 46,47 However, with the interaction strength J limited to the spin-singlet channel in the 't-J' model and the cuprate superconductors, such pairing amplitudes never appear in the Hamiltonian and thus cannot affect the energetics. 48,49 Similarly, no longer range spin-singlet even parity pairing amplitudes, beyond the considered nearest neighbor dand s-wave pairing amplitudes, are important as they do not either appear in the Hamiltonian in Eq. (1).
Having found a phase crystal state with accompanied swave pairing at low temperature in the cuprate superconductors, we next investigate how this exotic state develops with temperature and, in particular, how it compares with a uniform state (say θ = 0) throughout the superconductor, as that is also a viable superconducting solution to Eq. (1). In order to investigate the thermodynamic properties, we calculate the free energy Ω = E g − T S, where E g is the internal energy corresponding of the Hamiltonian Eq. (1) and S is the entropy (see Methods). In Fig. 2 we plot Ω for both the phase crystal state and the uniform phase state. The phase crystal state has a lower free energy than the uniform phase state for Thermodynamics for clean sample. Free energy Ω calculated in GIMT for phase crystal and uniform phase states as a function of scaled temperature T /Tc, where Tc is the d-wave superconducting transition temperature. Green arrow indicates the phase crystal transition temperature T * . Upper inset shows the phase transition temperatures T * and Ts obtained by plotting the phase sin(θ) rms and the magnitude of the extended s-wave pairing amplitude |∆s| rms, respectively, averaged along the [110] edge (a very small temperature-independent value of ∆s of the uniform phase state has been subtracted to extract Ts). Lower inset shows difference in free energy between the phase crystal state, Ωpc, and the metastable uniform phase state, Ωms, within GIMT and BdG. The energy jump in GIMT near T = T * is mainly from the internal energy Eg.
all temperatures below T * . Thus, T * defines the transition temperature of the phase crystal state, with its spontaneous breaking of both time-reversal and translational symmetry. At an even lower T = T s , the phase crystal free energy shows a marked downturn. This temperature corresponds to the appearance of the modulating s-wave pairing inside the phase crystal state, see upper inset in Fig. 2. Finally, to understand the effects of strong correlations on the stability of the phase crystal state, we show the free energy difference between the phase crystal state and the (metastable) uniform phase state in the lower inset of Fig. 2 for both GIMT and BdG calculations. The energy gain due to the formation of the phase crystals is much larger in GIMT compared to BdG. We further expect the enhancement of the energy gain in GIMT to be even more pronounced for a lattice with larger pair breaking [110] edges than the one considered here, since phase crystal is only occurring at the pair breaking edge. Also, as indicated by arrows, T * (GIMT) > T * (BdG). Both the large free energy gain and the higher T * demonstrate that the strong correlations give a notably enhanced thermodynamic stability of the phase crystal state.
The remarkably increased thermodynamic stability of the phase crystal state with strong correlations can be understood by looking at the eigenvalues of the Hamiltonian Eq. (1). In Fig. 3 we plot the distribution of the low-energy eigenvalues Eigenvalues En within GIMT and BdG in the uniform phase (a) and phase crystal (b) states at a low temperature, T = 0.08Tc < Ts < T * . Strong correlations generate a higher number of zeroenergy states in the uniform phase in (a). The phase crystal state shifts most zero-energy states to finite energies within both BdG and GIMT, but GIMT displays larger shifts and a full energy gap (min{En} = 0.006) due to the formation of co-existing s-wave state. Here, the eigenindex n is shifted by the total number of lattice sites N . E n , where n is the eigenindex, at a low temperature, for both the phase crystal and metastable uniform phase states and with and without strong correlations.
In the uniform-phase state, zero-energy Andreev bound states forms at [110] edges. 10,11 The origin of these zeroenergy states lies in the change of the sign of the d-wave superconducting order parameter experienced by quasiparticles Andreev scattering off the [110] edge. 10 Such a sign change does not occur at edges parallel to the crystallographic axes, such as the [100] edge, and hence the zero-energy states are only existing at the [110] edge in our sample. Moreover, since scattering trajectories are dependent on the quasiparticle band structure and its superconducting nodes, there exists an inherent band structure effect. Alternatively, and completely equivalently, the zero-energy states at [110] edges in a d-wave superconductor can be seen as having a topological origin set by the nodal structure of the quasiparticle band structure. 12,13 More specifically, the number of zero-energy edge states is equal to the number of states in the region of the edge Brillouin zone bounded by the projections of the superconducting nodes in the bulk band structure. 16 It is here important to emphasize that the topological nature of the edge states thus depends both on the orientation of the edge and the nodal structure band structure, which distinguishes the behavior topological edge states of a d-wave superconductor from the edge states of a topological insulator. As expected, we find zeroenergy states in the uniform phase of both GIMT and BdG, as shown in Fig. 3(a). However, the number of zero-energy states is significantly increased with strong correlations in GIMT.
We can relate this increase of the number of zero-energy states to a notable increase in the charge density at the edge always found when including strong correlations within GIMT, but not in BdG. This increase of the charge density at the edge can in turn be understood using a simple argument: Strong electronic repulsion tries to keep the electrons as far apart from each other as possible, or in other words, strong correlations smear out any charge accumulation. However, at the edge, this smearing effect is always reduced due to a reduced number of neighbors, and, as a result, the electrons accumulate at the edge compared to the bulk. In a sense, the enhancement of charge density at the edge can thus be thought of as a 'surface tension' caused by strong correlations. Notably, this result asks for a reconsideration of the straightforward use of the bulk-boundary correspondence in nodal superconductors, which is used to build a topological relation between the gap nodes and the number of zero-energy edge states. 13,50 A straightforward use of the bulk-boundary correspondence for the [110] edge of a d-wave superconductor would mean that the number of zero-energy edge states is set by the projections of the gap nodes in the bulk band structure onto the edge Brillouin zone. 16 However, the change in the edge charge density when including strong correlations results in a different effective band structure near the edge compared to that of the bulk. As a consequence, a reasonable understanding of the number zero-energy edge states can only be established if bands corresponding to the edge charge density are considered, instead of those corresponding to the bulk charge density, see SM for additional details Thus, the bulk-boundary correspondence should be applied using the modified edge band structure as the new "bulk" part in the correspondence. This is a clear example where the topological relation between the superconducting gap nodes and the zero-energy edge states in nodal superconductors needs to be reformulated in the presence of strong correlations.
In the phase crystal state, eigenvalues are also shifted to finite energies, see Fig. 3(b), making it the ground state below T * in both GIMT and BdG. Still, there exist three crucial ef-fects of strong correlations. First, the number of zero-energy states is significantly larger in the GIMT uniform phase than in the BdG uniform phase and consequently the energy gain in the phase crystal state, due to the shift of these states to finite energies, is much larger in GIMT. Secondly, due to the presence of the extended s-wave order, the GIMT phase crystal state features a full energy gap below T s , while there are always zero-energy states present in BdG calculations. Finally, even outside the full energy gap, the eigenvalues of the negative energy states are found at lower energies in GIMT compared to in BdG, and also display a moderate temperature dependence, see SM for additional data. These effects all contribute to making the phase crystal state more heavily preferred in the presence of strong correlations. We thus find that strong electronic correlations are crucial to generate both the correct stability and the full energy gap of the phase crystal state. Dirty superconductor edge. After having established a fully gapped phase crystal state at pair breaking edges of clean d-wave superconductors when accounting for strong correlations, we turn to another important aspect of cuprate superconductors, that of disorder. The bulk both possesses intrinsic disorder 14 and disorder is generated with chemical substitution and doping. 15 For systems with edges, disorder is even more important as it inevitably appears while preparing any edge. Since we are here primarily interested in edge properties, we only consider non-magnetic disorder for which the bulk d-wave superconducting state has been shown to be robust when including strong correlation effects. 40,41,51 In Fig. 4(a,b) we show the evolution of the phase of the d-wave pairing amplitude at a temperature within the phase crystal state for different disorder strengths V . The phase modulations near the pair breaking edge persist in the whole range from weak to strong disorder strength V = 1.5. This is also seen in the line plots of the phase at the pair breaking edge in Fig. 4(c). Moreover, the disorder averaged values of the phase crystal transition temperature T * and the accompanying extended s-wave transition temperature T s also remain constant with disorder, see Fig. 4(d). Thus, the phase crystal state is remarkably robust to disorder when strong correlations are appropriately included. The almost complete disorder insensitivity of the phase crystal is very surprising. It is true that strong correlations have previously been seen to be the reason for the robustness of d-wave superconductors to disorder, 40,41,51 but this robustness has always involved the magnitude of the pairing amplitude. In contrast, the phase crystal state exists only near pair breaking edges where the magnitude is suppressed and it is here instead the superconducting phase that is completely robust to disorder.
To further imprint the importance of strong correlations for the disorder robustness, we plot in Fig. 4(e,f) the phase of the d-wave pairing amplitude obtained within BdG for V = 0.5. Even for this relatively weak V = 0.5, where bulk dwave superconductivity clearly still thrives, 41 the edge phase modulations are almost completely disrupted. Thus, in the absence of strong correlations, the phase crystal state is clearly extremely sensitive to disorder, much more so than the bulk d-wave state. This further emphasizes the remarkable role of strong correlations in the stability of the phase crystal state. We have verified that these results hold also for other models of disorder, see SM.
The disorder robustness of the phase crystal state is even more surprising if we consider the energetic origin of the phase crystal state: a large degeneracy of zero-energy states in the uniform phase state. It has recently been shown using both T-matrix calculations 26 and topological arguments 27 that disorder generally weakens zero-energy state degeneracies. This we also see in the uniform phase state in Fig. 5(a), where we plot the low-energy eigenvalues E n . In the uniform state, the spectrum of eigenvalues is changed by disorder and the number of zero-energy states is indeed rather dramatically reduced. Based on this, we would naively expect that the phase crystal state would quickly disappear with increasing disorder. Instead, we see that disorder has very little effect on the phase crystal states and its low-energy states, when we appropriately include strong correlation effects in Fig. 5(a). As a consequence, the energy gain due to the formation of the phase crystal state is still higher than the disorder-induced shift of the energy spectrum in the uniform phase state. This results in the phase crystal state being robust to disorder in the presence of strong correlations, even though the degeneracy of the zero-energy states in the uniform state is nearly eliminated. In contrast to the disorder robustness when strong correlations are appropriately included, we find in Fig. 5(b) that the BdG eigenspectra of the uniform and phase crystal states start to overlap even for small disorder. As a result, the formation of a phase crystal state does not reduce the total free energy and thus the phase crystal state is destroyed even at weak disorder when strong correlations are ignored.
The striking disorder robustness of the phase crystal state in the presence of strong correlations has two underlying reasons. First, the enhanced number of zero-energy states in the uniform phase state of GIMT results in a larger energy gain of the phase crystal state for the clean sample, compared to the non-correlated case. This energy gain is clear from the lower inset of Fig. 2. All by itself, this large energy gain already suggests that the phase crystal within GIMT can withstand disorder up to a larger disorder strength compared to BdG. Second, strong electronic repulsion smears out the charge inhomogeneities caused by the presence of disorder. This results in a strongly weakened effective disorder. 41 These two features cumulatively make the phase crystal state robust to disorder in the presence of strong correlations.
Our results demonstrate that the combined effects of strong correlations, topology, and disorder can be very non-intuitive. The phase crystal state, which initially is formed due to a large degeneracy of topologically protected zero-energy states, is extremely robust to disorder but only in the presence of strong correlations. This result holds even though the original degeneracy of the zero-energy states is drastically reduced by disorder.
Comparison with experiments. Above we established the existence of a disorder robust phase crystal state with a full gap at the lowest temperatures at edges of strongly correlated high-temperature cuprate superconductors. These findings have direct experimental application in superconducting cuprate devices since the presence of zero-energy edge states without the phase crystal state should, if indeed present, significantly affect transport properties. In fact, different experiments on cuprate devices have for a long time given contradictory results, with no previous viable theoretical explanation. At intermediate temperatures, tunneling experiments across cuprate-normal metal junctions unanimously see a large zerobias conductance peak, 10,52 consistent with the existence of many zero-energy edge states. However, at lower temperatures, some experiments reveal a full gap in the tunneling spectra, [28][29][30] whereas others only see a so-called temperatureindependent broadening, 31,32 i.e. broadening beyond standard thermal broadening through the Fermi-Dirac distribution, of the zero-bias conductance peak. This dichotomy at low temperatures has previously been notoriously hard to explain theoretically. A set of theories, based on either spontaneous formation of subdominant s-wave superconducting order, 19,20 magnetic order, 16,17 or supercurrents 18 all give only a split in the zero-energy density of states at low temperatures. On the other hand, a phase crystal state in the absence of strong correlations only explains the temperature-independent broadening of the zero-energy states 22 but not the full gap. By including strong correlation effects we remedy these theory shortcomings and can now explain all experimental data.
To connect to the experiments, in Fig. 6 we show the spa- tially averaged low-energy density of states N (E) obtained within GIMT at different temperatures in clean (a) and disordered (b) samples. To numerically evaluate N (E), we choose a small and temperature independent energy smearing in N (E) (see Methods for details). We do not expect broadening due to any other scattering processes to be temperature dependent for the low temperatures considered in this work, as the energy scales associated with these temperatures are much smaller than the superconducting gap. We find that for both clean and disordered samples, at high temperatures, T > T * , the zero-energy states of the uniform phase state persist, giving a pronounced zero-bias peak in the N (E). Below a lower temperature, T * , the phase crystal state develops, which induces a broadening of the lowest energy peak. This is most clearly seen in the temperature evolution of the full width at half maximum (FWHM) of the zero-energy peak, which we plot in the insets of Fig. 6. Most importantly, the zero-energy peak broadening increases with decreasing temperature and exists independently of the thermal broadening (which is not included), therefore being a temperature-independent effect. It is caused by an increase in the Doppler shift of the low energy states in the phase crystal state with lowering temper-atures, see also SM for the explicit temperature dependence of the eigenvalues. Finally, at an even lower temperature, T s , an additional extended s-wave component appears causing a full energy gap to develop. Thus, it is very feasible to experimentally measure, across different devices with their different T * and T s due to microscopic details, both a full energy gap and a temperature-independent broadening of the zerobias conductance peak at low temperatures, while at intermediate temperatures consistently see the large zero-bias peak of the uniform state. This illustrates how our findings likely resolve the long-standing experimental dichotomy and provide future guidelines for the design of superconducting devices. Discussion. The strong disorder robustness of the phase crystal state in strongly correlated superconductors makes it a very promising candidate for finally explaining the physics of boundaries and edges in the cuprate superconductors. The disorder robustness of the phase crystal state is also remarkable from the point of view that it concerns the phase of the superconducting pairing, not the magnitude. Conventional s-wave superconductors have for a long time been known to be highly robust against disorder thanks to the celebrated Anderson's theorem. 53 In the presence of strong correlations, an analogy of the Anderson's theorem has recently been constructed for d-wave superconductors. 54 However, these results both only concern the magnitude of the superconducting pairing. The phase crystal state on the other hand is characterized by its strong modulation of the phase of the superconducting pairing and it even only exists in regions where the magnitude is suppressed. Thus, there exists no established reason to expect a phase crystal state to be anything but unstable towards disorder. This disorder sensitivity is even verified in the absence of strong correlations, where the phase crystal is very easily destroyed by disorder. It is only when including strong correlation effects that the phase crystal state is stabilized to the degree that it can be present in a real material. Thus our results establish the importance of strong correlations for disorder robustness of superconducting phase modulations. Since disorder robustness is often used to predict the underlying pairing symmetry of various superconductors, our findings can additionally play a pivotal role in determining the pairing symmetry of newly discovered superconductors such as nickelates or twisted bilayer graphene, where electronic correlations are considered to be strong. 55,56 Methods In this section, we first provide details of the two methods, GIMT and BdG, implemented to solve the Hamiltonian in Eq.
(2) to extract the self-consistent superconducting pairing amplitude, eigenvalues, and thermodynamic quantities. Then we provide values and motivations for the lattice geometry and different parameters used to obtain the results.
(2) requires handling the projection operators that prohibit formation of double occupancy on any lattice site. Within the Gutzwiller inhomogeneous mean-field theory (GIMT) the effects of this projection is treated locally, i.e. at each site, using the Gutzwiller approximation. To implement the Gutzwiller approximation, we first consider the Gutzwiller wave function 34 (2) are assumed to be related by 35,58,59 where g t ij and g J ij are known as the Gutzwiller renormalization factors. In the absence of magnetic order and ignoring intersite correlations, g t ij and g J ij are given in terms of the local hole doping xi: 38,40 g t ij = 4xixj (1 + xi)(1 + xj) and g J ij = 4 (1 + xi)(1 + xj) . (4) Using the relations obtained by the Gutzwiller approximation in Eq.
(3), the internal energy of the Hamiltonian in Eq.
(2), Eg = ψ0|H|ψ0 , can be expressed as where ρi = σ niσ = σ niσ 0 is the local electron density, ∆ij = c j↓ c i↑ 0 − c j↑ c i↓ 0 is the spin-singlet Cooper pairing amplitude on each nearest neighbor bond, and τij = c † i↓ c j↓ 0 = c † i↑ c j↑ 0 is the particle-hole bond expectation value. The local doping is given by xi = 1 − ρi. In writing Eq. (5), we assume full spin rotational symmetry (due to no magnetic order) and no spin-triplet superconducting order, both appropriate for cuprate superconductors. Finally, the extremum of Eg with respect to the different mean-field variables gives the GIMT Hamiltonian 39,41,60-62 which, when combining Eqs. (5) and (6), can be written as where now nearest neighbor sites of i are denoted as j = i + δ, δ = ±x, ±y, and the next-nearest neighbor sites as j = i +δ, δ = ±(x±y). Consequently, ∆ij, which resides on nearest neighbor bonds, is denoted as ∆ δ i , and similarly for τ . Here, W FS iδ and µ HS i are the Fock and Hartree shifts, respectively, given by The derivatives of the Gutzwiller factors in Eq. (8) are calculated analytically using the expressions in Eq. (4). Notably, these derivatives suppress any change in the charge density between the site i and its neighbors j and are thus responsible for reducing charge accumulation within GIMT. The mean-field Hamiltonian HGIMT in Eq. (7) is diagonalized using the Bogoliubov-de Gennes transformations, 42 where γ † nσ and γnσ are the creation and annihilation operators of the Bogoliubov quasiparticles and ui,n and v * i,n the eigenfunctions with eigenvalues En,. The resulting eigensystem is then solved self-consistently for all independent local variables: and ρi. Finally, in order to analyze the pairing amplitude we project ∆ δ i on its local dand swave components ∆ d,s (i) as outlined in the main text. In this work we always report on the superconducting pairing amplitudes ∆ δ i , as is commonly done in GIMT calculations, 40,41,61,62 while technically the superconducting order parameter is given by g t i,i+δ ∆ δ i . 35,63 With both ∆ δ i and τ δ i possibly being complex-valued, we obtain self-consistency on both the real and the imaginary parts. A crucial aspect in this process is the initial guess of the self-consistent variables. With real-valued inputs for ∆ δ i and τ δ i , the solution always converges to a real ∆ δ i . Thus, a purely real ∆ δ i is always a solution of the Hamiltonian, albeit it might not be the ground state. In order to converge to the true ground state, that might also include complex values of ∆ δ i , complex-valued initial guesses of either ∆ δ i or τ δ i are essential. We have tried many different inputs as starting guesses. All results for the phase crystal state reported in this work are obtained using completely random complex initial inputs for both the amplitude and phase of ∆ δ i , while the uniform phase state is obtained using completely random real-valued amplitude inputs with zero phase.
We also extract several different physical properties. In order to compare the thermodynamic stability of different solutions and determine the ground state, we calculate the free energy Ω = Eg − T S, where Eg is given by Eq. (5) and S is the entropy given by where f is the Fermi-Dirac distribution function.
We also calculate the spatially averaged density of states N (E) = 1/N i,n g t i,i (|ui,n| 2 δ(E−En)+|vi,n| 2 δ(E+En)). Here N is the total number of lattice sites. To numerically evaluate N (E), we use a Lorentzian with fixed width 0.0015 to calculate the delta-function and implement a small and temperature independent energy smearing. The fixed width explicitly excludes the standard thermal broadening effect caused by the Fermi-Dirac distribution of electrons, allowing us to isolate effects beyond the standard thermal broadening. BdG. To identify the effects of strong electronic correlations, we compare the results of GIMT with a standard Bogoliubov-de Gennes (BdG) calculation, 42 where the sole effect of the correlations becomes the creation of the superconducting pairing amplitude. Within the BdG method, the (disordered) 't-J' model Eq.
(2) is thus solved by ignoring the effects of projection, i.e., the Hilbert space does not have any restriction on the formation of double occupancy, and mean-field decoupling of the super-exchange interaction term is only performed in the Cooper pairing channel. As a consequence, the BdG method does not capture the no-double occupancy constraint imposed by the strong correlations, nor does it include the Hartree (µ HS i ) and Fock (W FS iδ ) shifts also given by the strong correlations. However, the lattice scale effects (short wavelength fluctuations) that are relevant for high-temperature superconductors with short superconducting coherence lengths are still included in a BdG solution. The resulting BdG mean-field Hamiltonian is given by Practically this BdG Hamiltonian is obtained by setting the Gutzwiller factors g t ij and g J ij to unity and ignoring the Hartree (µ HS i ) and Fock (W FS iδ ) shifts in the GIMT Hamiltonian in Eq. (7). Then, we follow the same iterative self-consistency treatment as in GIMT, but here self-consistency is only needed for ∆ δ i . Lattice geometry. We primarily consider the lattice geometry shown in Fig. 7, using N = 1926 lattice points. This particular lattice geometry gives us an optimum length for the [110] edge (dashed line in Fig. 7) with minimum interference from other edges. This is important since the [110] edge is pair breaking for d-wave pairing on nearest neighbor bonds (technically d x 2 −y 2 symmetry, but here abbreviated as d-wave) as the Andreev reflection at this edge mixes opposite signs of the superconducting pairing amplitudes. We have also performed calculations for other system sizes with similar qualitative outcomes, provided the [110] edge is long enough to host multiple modulations of the phase crystal state. Parameters. Within GIMT we set t = 1 and t = 0.5 to mimic a normal-state band with Fermi velocities at the anti-nodes being very similar to the nodal ones, which is the case for cuprate materials such as Bi2Sr2CaCu2O8+x (BSCCO), 64 while only using a minimal number of hopping parameters. We fix the average density to be ρav = 0.8 (i.e. average hole doping 0.2), which gives a doping slightly below the optimal doping for superconductivity for the chosen parameters. Our choice of ρav allows us to ignore any competing magnetic order (which appears near the extremely underdoped regime) or competing extended s-wave superconducting order (which appears in the extremely overdoped regime) in the bulk. Further, we set the super-exchange interaction to J = 0.33, a value close to the one obtained in scattering experiments on cuprates. 65,66 The chosen value of J corresponds to a large value of U = 12, mimicking the strong Coulomb interaction in cuprates. In the bulk this gives rise to a bulk superconducting Tc = 0.148 and |∆ d | = 0.236. We have also checked that the qualitative findings of this work do not change for other lower values of J. In terms of the Gutzwiller renormalization factors, we find that they are all significant but that both the hopping amplitudes tij and the super-exchange interaction J are essentially the same in the bulk as on the edge. However, through the derivatives in Eq. (8), the local Hartree shift of the chemical potential differs substantially, being around 1.5 in the bulk but 1.15 at the edge.
For the BdG calculations we choose the same parameters as used in Ref. [25] to study the phase crystal state. In particular, within BdG we set t = 1, t = 0.25, J = 1.4, and µ = 0.0, giving an average density to be ρav = 1.18 and the bulk superconducting Tc = 0.122 and |∆ d | = 0.108. These parameters both reproduce earlier data in the literature for BdG calculations and also generate normal-state band parameters that are very close to those of the GIMT calculations. We note however that both the density and transition temperatures are different compared to in GIMT, but we can still compare the GIMT and BdG solutions by reporting properties rescaled by Tc. In the SM we provide an alternative BdG (alt-BdG) calculation where we choose the BdG parameters explicitly such that both the bulk average density ρav and pairing amplitude ∆ d (i) are matching those of the GIMT calculations. This allows us to compare the number of zero-energy states in the uniform phase state, which is the purpose of the SM work. However, the phase crystal state become extremely weak in this alt-BdG setup, with also a large s-wave component, making the edge region being almost instead in a d + is-wave state. This result in fact illustrates a common feature of our BdG calculations, that they are surprisingly hard to converge into a phase crystal state. Including strong correlations, this problem completely disappears and the phase crystal is stable in a large parameter regime.
Supplementary Material for "Disorder-robust phase crystal in high-temperature superconductors stabilized by strong correlations" In this Supplementary Material (SM) we provide additional details to further support some of the results in the main text. In particular, we first show the circulating currents near the pair breaking edge obtained within GIMT, then we show the temperature dependence of the low-energy eigenvalues of the ground state for a clean superconductor, then we provide a detailed discussion on the relation between the number of zero-energy states and the local charge density in the uniform phase state of a clean superconductor, and finally we illustrate the disorder robustness of phase crystals in GIMT with a different model for the disorder. In the main text, we state that time-reversal symmetry is broken for T < T * due to the formation of the phase crystal with its modulating phase. In this section, we show that the broken time-reversal symmetry also leads to the appearance of circulating currents. Using the charge continuity equation together with the Heisenberg equation for the particle density (see e.g. Refs. 67 and 68), the total particle (electron) current on a bond between site i and i + δ for the Hamiltonian in Eqs. (1) and (7) is given by With J i,i+δ only present on bonds, we also define for visualization purposes a vector field where δ = x, y, y + x, y − x andδ is the unit vector along the direction of δ. In Fig. S1 we illustrate the spontaneous currents in a clean sample within GIMT at two different temperatures. Here we use a color density plot for the magnitude of the current vector field J i along with a vector plot for the direction of the current vector field. As seen in Fig. S1(a) for T s < T < T * , i.e. in the phase crystal state but without extended s-wave pairing, we find alternating circulating currents near the pair breaking edge. This result is very similar to the findings of quasi-classical theory 22,24 and BdG 25 for the phase crystal state. However, for T < T s , where extended s-wave pairing also appears near the nodes of the phase crystal, this pictures changes. Within this temperature regime, local d + is-wave regions are formed and they add additional currents 19 on top of the circulating currents from the phase crystal. The result is the current pattern in Fig. S1(b), where the current is clearly finite but does not show a simple circulating pattern. The temperature evolution of En of the ground state (the state with lowest free energy: the phase crystal state for T < T * and the uniform phase state for T > T * ) are shown for GIMT (a) and BdG (b). The full gap in GIMT disappears for T > Ts with the slope of the curve increasing very slowly as T * is approached. In contrast, the slopes of the curves change monotonically with temperature in BdG. Note that En for T /Tc = 0.2 and T /Tc = 0.27 are not shown in (b), as the distribution does not change with temperature for T * < T < Tc.
In Fig. 3 of the main text we show the distribution of the low-energy eigenvalues of a clean sample but only for one specific low temperature. To offer complementary data, we plot in Fig. S2 the temperature evolution of the eigenvalues E n of the ground state, i.e. the phase crystal state for T < T * and the uniform phase state for T > T * . Within the GIMT, where strong correlation effects are included, we see in Fig. S2(a) that there are two temperature regimes for T < T * . First, for T ≤ T s , where T s marks the emergence of an extended s-wave component at the edge, the full gap in E n around zero energy reduces with increased temperature. Furthermore, with increasing temperature, more states move toward zero energy, as seen from the increased curve slopes. The rate of slope increase is strong in this temperature regime. Secondly, for T s < T ≤ T * the rate at which the states move towards zero energy is much slower with increasing temperature and eventually the distribution E n reaches that of the uniform phase state found at T = T * . We note here that these temperature dependencies found below T * are affecting the eigenvalues of the system, and are as such a temperature effect that goes beyond the standard thermal broadening of physical properties caused by the Fermi-Dirac distribution of electrons. As such, they will cause temperature effects, such as conductance peak broadening, that cannot be related to standard Fermi-Dirac temperature broadening. In the uniform phase state the lowenergy E n are instead insensitive to temperature. In Fig. S2(b) we present the corresponding plot within BdG where strong correlations effects are ignored. Here, the rate at which the E n distribution approaches the uniform phase state is essentially monotonic in temperature and there is also no emerging s-wave component leaving the spectrum ungapped at all temperatures. This contrasting temperature evolution of the eigenvalues is also responsible for a larger T * in GIMT.

III. RELATION BETWEEN NUMBER OF ZERO-ENERGY STATES AND LOCAL CHARGE DENSITY
In this section, we discuss in more detail the reasons behind the increased number of zero-energy states at the edge in the presence of strong correlations, as shown in Fig. 3(a) of the main text, and why an application of the bulk-boundary correspondence in nodal superconductors using bulk properties is unattainable when strong correlations effects are included.
Strong electron correlations modify the properties of a system in primarily two ways. First, strong correlations renormalize the hopping amplitudes and the super-exchange interaction through the Gutzwiller factors, even in the absence of any inhomogeneity. Secondly, these renormalization factors, g t ij and g J ij , depend on the local charge density and consequently alter the effects of any inhomogeneity. In order to probe the importance of strong correlations for the zero-energy edge states, we need to focus on the inhomogeneous effects caused by the edge. We can do this by comparing our strongly correlated results with an alternative BdG (alt-BdG) calculation where we choose parameters such that the bulk solution in the absence of any disorder is exactly the same as in GIMT. Thus, in the alt-BdG calculation we choose t, t , J and µ such that the bulk homogeneous solution, both charge density and superconducting pairing amplitude, of Eq. (7) and Eq. (9) match exactly. In Fig. S3(a,b) we show the resulting spatial color maps of the local charge density ρ i in the uniform phase state in a clean sample obtained in GIMT and alt-BdG, respectively. By construction, the bulk charge densities are exactly the same in both plots. However, the edge density shows alt-BdG contrasting features: within GIMT ρ edge > ρ bulk , while within alt-BdG ρ edge < ρ bulk . The increase in the charge density at the edge in GIMT is mainly due to the local Hartree shift µ HS i given in Eq. (8). The most significant contribution to µ HS i comes from the derivatives of g t i,i+δ on the nearest neighbor bonds, i.e. the first term in Eq. (8). Since ∂g t i,i+δ /∂ρ i < 0 and the sites on the edge have only two nearest neighbors, µ HS edge < µ HS bulk . As a consequence, within GIMT edges very generally host a larger charge density than the bulk, not just for the system parameters used to produce Fig. S3. In contrast, the charge density at system edges in (alt-)BdG, where strong correlation effects are ignored, is either lower (as in Fig. S3(b)) or very similar to the bulk density (as in Fig. 1(d)), depending on the bulk band structure. To summarize, the charge density behavior at the edges in Fig. S3(a,b) is very generic for GIMT and BdG solutions, respectively. Now that we have fixed the bulk to be the same in GIMT and alt-BdG, we can directly compare the number of zero-energy edge states between the two methods and check how they compare and also how that compares with the prediction of the bulkboundary correspondence. In Fig. S3(c) we plot the distribution of the low-energy eigenvalues for the uniform phase state in a clean superconductor. We find that the number of zero-energy states is much larger in GIMT than in alt-BdG, quite similar to the results in the main text where we instead compared GIMT to BdG. Most notably, we find a very different number of zero-energy states in GIMT and alt-BdG, despite exactly the same bulk properties. Thus we have explicitly shown that applying the bulk-boundary correspondence using the bulk properties cannot possibly generate the correct number of edge states in both cases To explain this discrepancy between number of edge states and the bulk properties, we plot in Fig. S3(d) the Fermi surface for three different cases: (1) a band giving the bulk charge density (yellow), (2) a band giving a bulk density equal to the edge density of GIMT (green), and (3) a band giving the bulk density equal to the edge density of alt-BdG (blue). Due to translational invariance along x || , the edge states of the [110] edge form one-dimensional flat bands at zero energy that terminate at the bulk superconducting nodal points, 16,69 indicated for each case by the red dots in Fig. S3(d). The number of zero-energy edge states is thus proportional to the distance between these bulk nodal points, indicated by the black double-arrows and Λ i labels for the three bands i = 1, 2, 3. As a consequence, the number of zero-energy edge states is directly related to the normal-state band structure and its corresponding charge density. However, the problem in a real system is to know which charge density to use when it is changing between the edge and the bulk? With the edge states physically most connected to the edge charge density, it is natural to assume that the number of zero-energy edge states is more related to the edge charge density than that of the bulk. Based on this assumption connecting the edge charge density with number of zero-energy states, the number of zero-energy states becomes proportional to Λ 3 in alt-BdG, but Λ 2 in GIMT, as shown in Fig. S3(c,d). Since Λ 2 > Λ 3 and also Λ 2 > Λ 1 , we always get more zero-energy states in the GIMT solution than in any BdG solution, even if the bulk conditions are exactly the same, as we also see in our numerical results. We thus conclude that a simple-minded application of the bulk-boundary correspondence in nodal superconductors to extract the number of zero-energy states from the bulk band structure and its charge density is not correct, especially not in the presence of strong correlations where the charge density at the edge always increases, irrespective of the normal state in the bulk. Additionally, our results show that including strong correlations always leads to more zero-energy edge states due to an enhanced edge charge density, which enhances the stability of the phase crystal state in the presence of strong correlations.

IV. DISORDER ROBUSTNESS OF PHASE CRYSTAL STATE FOR CONCENTRATION DISORDER
In the main text, we show the results for 'Anderson type' disorder, where a random potential V i is added to each lattice site. To verify that our results are not sensitive to this particular choice of disorder, we also consider another model of disorder, referred to as concentration disorder, where V i = V 0 = 1 only on a given n imp fraction of randomly chosen lattice sites, while V i = 0 for all other sites. Thus concentration disorder contains fewer but stronger scatterers, while Anderson disorder instead generates a continuum of relatively weaker scatterers. The results when using concentration disorder are shown in Fig. S4, where we plot the color density maps of sin(θ) in both GIMT and BdG for different concentrations n imp . The remarkable robustness of the phase crystal state in GIMT persists even up to the very high disorder concentration n imp = 20%, which only induces minor irregularities in the phase crystal modulations. In contrast, the phase crystal state in BdG is extremely sensitive to concentration disorder, with the modulations getting disrupted even for a low concentration n imp = 5%. We thus conclude that the disorder results obtained in the main text are not sensitive to the type of the disorder and the phase crystal state is robust to all types of (non-magnetic) disorder when strong correlations are appropriately included.