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


INTRODUCTION
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 pair-wise [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 zero-energy 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 dwave 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 realworld 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 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 y 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 t ij = −t for nearest-neighbor bonds and t ij ¼ t 0 for next-nearest-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 inhomogeneities.
We perform a Gutzwiller inhomogeneous mean-field theory (GIMT) treatment of the Hamiltonian in Eq. (1) 40,41 . The superexchange J term gives rise to superconductivity with spin-singlet Cooper pairs living on nearest-neighbor bonds, which can be represented by the local d-and 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Þ j jexpð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 nonmagnetic 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. 1a, 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 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 Note 1.
To start characterizing these distinct properties, we investigate in Fig. 1c, d the behavior of the relative change of the magnitude of the d-wave pairing amplitude, δ Δ d j j, 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 j j, 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 j j. 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 the "Methods" section), 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 j j 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 j j is closely following the relative increase δρ of the charge density, as seen in Fig. 1c. To summarize, both a short bulk superconducting coherence length and a short charge density healing length result into a short healing length of Δ d j j 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. 1d. 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 zeroenergy 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. 1a.
Another striking feature of strong correlations is the appearance of an extended s-wave pairing amplitude Δ s j j at the edge with a modulation length scale set by that of the phase crystal, as seen in Fig. 1e, while the s-wave amplitude is always negligible in the bulk. Although not shown in Fig. 1e, 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 zero-energy 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 timereversal 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 D. Chakraborty et al. 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 d-and 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 s-wave 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 −TS, where E g is the internal energy corresponding of the Hamiltonian Eq. (1) and S is the entropy (see the "Methods" section). 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 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 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 form at [110] edges 10,11 . The origin of these zero-energy states lies i rms , and relative change in the charge density defined as δρ = 〈ρ〉 rms − ρ(bulk) plotted as a function of x ⊥ for GIMT and BdG, respectively. Here 〈. . . 〉 rms denotes the root mean square average over all the sites in the x ∥ direction for a fixed x ⊥ . e Line plots of the phase modulations together with the magnitude of the extended s-wave pairing amplitude Δ s j j at the pair breaking edge. See the "Methods" section for remaining parameters.
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 band structure, which distinguishes the behavior of topological edge states of a d-wave superconductor from the edge states of a topological insulator. As expected, we find zero-energy states in the uniform phase of both GIMT and BdG, as shown in Fig. 3a. 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 of zeroenergy states can only be established if the bands corresponding to the edge charge density are considered, instead of those corresponding to the bulk charge density, see Supplementary Note 3 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. 3b), making it the ground state below T * in both GIMT and BdG. Still, there exist three crucial effects 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 E n within GIMT and BdG in the uniform phase (a) and phase crystal (b) states at a low temperature, T = 0.08T c < T s < T * . Strong correlations generate a higher number of zero-energy states in the uniform phase in (a). The phase crystal state shifts most zeroenergy states to finite energies within both BdG and GIMT, but GIMT displays larger shifts and a full energy gap (minfE n g ¼ 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.
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 Supplementary Note 2 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. 4a, 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. 4c. 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. 4d).
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. 4e, 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 d-wave 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 Supplementary Note 4).
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. 5a, 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. 5a. 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. 5b 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 noncorrelated 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 zero-bias 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 temperature-independent 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 spatially 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 the "Methods" section 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 zeroenergy 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 temperatures, see also Supplementary Note 2 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 zero-bias 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 recently 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 selfconsistent 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.

GIMT
Solving the (disordered) 't−J' Hamiltonian in Eq. (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,36,57 ψ j i ¼ P ψ 0 j i, where ψ 0 j i is the ground state wave function in the Hilbert space that allows double-occupancy and P is the projection operator: (2) are assumed to be related by 35,58,59 hc y iσ c jσ i % g t ij hc y iσ c jσ i 0 ; 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 x i

38,40
: Using the relations obtained by the Gutzwiller approximation in Eq. (3), the internal energy of the Hamiltonian in Eq. (2), E g ¼ ψ 0 h jH ψ 0 j i, can be expressed as where ρ i ¼ P σ hn iσ i ¼ P σ hn iσ i 0 is the local electron density, Δ ij ¼ hc j# c i" i 0 À hc j" c i# i 0 is the spin-singlet Cooper pairing amplitude on each nearest-neighbor bond, and τ ij ¼ hc y i# c j# i 0 ¼ hc y i" c j" i 0 is the particle-hole bond expectation value. The local doping is given by x i = 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 E g 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, D. Chakraborty et al.
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 H GIMT in Eq. (7) is diagonalized using the Bogoliubov-de Gennes transformations 42 , c iσ ¼ P n ðγ nσ u i;n À σγ y nσ v Ã i;n Þ, where γ y nσ and γ nσ are the creation and annihilation operators of the Bogoliubov quasiparticles and u i,n and v Ã i;n the eigenfunctions with eigenvalues E n ,. The resulting eigen-system is then solved selfconsistently for all independent local variables: Finally, in order to analyze the pairing amplitude we project Δ δ i on its local d-and s-wave 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 selfconsistency on both the real and the imaginary parts. A crucial aspect in this process is the initial guess of the self-consistent variables. With realvalued 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 Ω = E g −TS, where E g is given by Eq. (5) and S is the entropy given by S ¼ À2 P n f ðE n Þ ln f ðE n Þ þ ð1 À f ðE n ÞÞ lnð1 À f ðE n ÞÞ ½ , where f is the Fermi-Dirac distribution function. We also calculate the spatially averaged density of states NðEÞ ¼ 1=N P i;n g t i;i ðju i;n j 2 δðE À E n Þ þ jv i;n j 2 δðE þ E n ÞÞ. 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 hightemperature superconductors with short superconducting coherence lengths are still included in a BdG solution. The resulting BdG meanfield 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

Parameters
Within GIMT we set t = 1 and t 0 ¼ 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 Bi 2 Sr 2 CaCu 2 O 8 + 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 T c = 0.148 and Δ d j j ¼ 0:236. We have also checked that the In terms of the Gutzwiller renormalization factors, we find that they are all significant but that both the hopping amplitudes t ij and the superexchange 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 ¼ 0:25, J = 1.4, and μ = 0.0, giving an average density to be ρ av = 1.18 and the bulk superconducting T c = 0.122 and Δ d j j ¼ 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 T c . In the Supplementary Note 3 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 Supplementary Note 3. 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.

DATA AVAILABILITY
The data are available from the corresponding author upon reasonable request.

CODE AVAILABILITY
The codes are available from the corresponding author upon reasonable request.