Critical nematic correlations throughout the superconducting doping range in Bi2−zPbzSr2−yLayCuO6+x

Charge modulations have been widely observed in cuprates, suggesting their centrality for understanding the high-Tc superconductivity in these materials. However, the dimensionality of these modulations remains controversial, including whether their wavevector is unidirectional or bidirectional, and also whether they extend seamlessly from the surface of the material into the bulk. Material disorder presents severe challenges to understanding the charge modulations through bulk scattering techniques. We use a local technique, scanning tunneling microscopy, to image the static charge modulations on Bi2−zPbzSr2−yLayCuO6+x. The ratio of the phase correlation length ξCDW to the orientation correlation length ξorient points to unidirectional charge modulations. By computing new critical exponents at free surfaces including that of the pair connectivity correlation function, we show that these locally 1D charge modulations are actually a bulk effect resulting from classical 3D criticality of the random field Ising model throughout the entire superconducting doping range.

Charge modulations have been widely observed in cuprates, suggesting their centrality for understanding the high-T c superconductivity in these materials. However, the dimensionality of these modulations remains controversial, including whether their wavevector is unidirectional or bidirectional, and also whether they extend seamlessly from the surface of the material into the bulk. Material disorder presents severe challenges to understanding the charge modulations through bulk scattering techniques. We use a local technique, scanning tunneling microscopy, to image the static charge modulations on Bi 2−z Pb z Sr 2−y La y CuO 6+x . The ratio of the phase correlation length ξ CDW to the orientation correlation length ξ orient points to unidirectional charge modulations. By computing new critical exponents at free surfaces including that of the pair connectivity correlation function, we show that these locally 1D charge modulations are actually a bulk effect resulting from classical 3D criticality of the random field Ising model throughout the entire superconducting doping range.
Charge order (CO), long seen on the surface of Bi2212 1-3 , Bi2201 4 , and Na-CCOC 3,5 , has recently been demonstrated in the bulk of many superconducting cuprates by NMR and scattering techniques [6][7][8][9][10][11][12][13][14][15][16][17][18] . Its apparent universality prioritizes its microscopic understanding and the question of its relationship to superconductivity (SC). However, severe material disorder presents both a challenge and an opportunity 19 . The challenge is that material disorder disrupts longrange order and limits macroscopic experimental probes to reporting spatially averaged properties. In particular, while numerous theories rest upon the 1D (stripe) or 2D (checker) nature of the CO, bulk probes may collect signal from multiple domains, obscuring the underlying dimensionality within a single domain.
The opportunity is for local probes to employ disorder as a knob that spatially varies parameters such as doping and strain within a single sample, to test and quantify the relationship of CO to SC. However, this strategy rests on the premise that what is seen on the surface is not merely a surface effect, but is reflective of the bulk of the sample. In Bi2201, while CO has been observed in the bulk (via, e.g., resonant X-ray scattering 9 ) with the same average wavevector as on the surface, it has not yet been demonstrated that they are locally the same phenomenon. That work also left open the question of whether the CO is locally unidirectional ("stripe-like") or bidirectional (a "checkerboard"). Here, we combine a local probe, scanning tunneling microscopy (STM), with a theoretical framework known as cluster analysis 19 , appropriate near a critical point, in order to test whether the surface CO is connected to the bulk CO. We find that the charge modulations in Bi2201 have significant stripe character. By computing new critical exponents at free surfaces including that of the pair connectivity correlation function, we moreover show that these charge modulations pervade the bulk of the sample, and that their spatial correlations are critical throughout the doping range of superconductivity.
We use STM to study the cuprate high-temperature superconductor Bi 2−z Pb z Sr 2−y La y CuO 6+x (Bi2201) at the dopings shown in Fig. 1a, from underdoped to overdoped superconducting samples, as well as an optimally doped sample with superconducting transition temperature T c = 35 K, as a function of hole concentration p. Figure 1b shows a topographic image of slightly underdoped Bi2201 with transition temperature T c = 32 K (UD32K), drift-corrected as described in Ref. 20. La and O doping supplies the holes, while Pb doping suppresses the structural supermodulation, leaving only the atomic corrugations with a periodicity of a 0 = 3.8 Åbetween copper atoms in the Cu-O planes.

Stripes vs. checkers
To identify the nature of the charge modulations, we focus on the Rmap, where R(r, V) = I(r, V)/I(r, −V), and I(r, ±V) represents the STM tunneling current at ±V as a function of position r on the surface of the sample 3 . The R-map has the advantage that it cancels out certain unmeasurable quantities, such as the tunneling matrix element and tunnel barrier height. Figure 1c shows the R-map with V = 100 mV in the same field of view (FOV) as Fig. 1b. A local modulation with period near 4a 0 is readily apparent, as confirmed by the two-dimensional Fourier transform (FT) of the R-map in Fig. 1d, showing peaks at Q ** x ∼ ð ±3=4,0Þ2π=a 0 and Q ** y ∼ ð0, ±3=4Þ2π=a 0 . The Q ** peaks carry information about the same charge modulation as the peaks at Q * x ∼ ð±1=4,0Þ2π=a 0 and Q * y ∼ ð0, ±1=4Þ2π=a 0 21,22 , and because they are well-separated from the central broad FT peak, there is less measurement error associated with tracking Q ** . We therefore focus on the Q ** peaks.
There has been experimental evidence in several families of cuprate superconductors for both stripe order (unidirectional CDW) 2,23-25 and checkerboard order (bidirectional CDW) 1,[26][27][28] . Recent experiments on YBCO show that the issue of dimensionality of the CO in cuprates can be quite complex: while the zero field CDW is 2D correlated, the field-induced CDW is both 3D correlated 29 and unidirectional 30 . Recent experiments on LSCO point to a multi-faceted relationship between CO and SC: while short-range CO exists above the superconducting dome, the CO is unidirectional in the superconducting regime 31 . It is difficult to discern from direct observation which tendency (stripes or checkerboards) would dominate in a hypothetical zero disorder limit 32,33 , because the quenched disorder that is always present in real materials can favor the appearance of stripe correlations 33 . One metric for distinguishing whether the underlying electronic tendency favors stripes or checkerboards is to compare the correlation length of the periodic density modulations ξ CDW with the correlation length of the orientation of the modulations ξ orient . Two different theoretical approaches 32,33 predict that ξ orient > 1 2 ξ CDW when the underlying tendency is toward stripes rather than checkerboard modulations.
In order to infer whether the charge modulations would tend toward stripes or checkerboards in Bi2201 in a hypothetical zero disorder limit, we construct the local Fourier components of the R-map at  . b Constant-current STM topography of UD32K sample acquired at I = 400 pA and V s = − 200 mV, over a 30 nm × 30 nm field of view. The arrows corresponds to the two orthogonal Cu-O bond directions throughout this paper. c A typical tunneling asymmetry R-map taken at 100 mV (i.e. R(r, 100 mV) = I(r, 100 mV)/I(r, − 100 mV)) in the same field of view shown in (b). A disordered charge modulation with a period of~4a 0 is evident in real space. d Fourier transform (FT) of R-map over the entire FOV in (c), with Bragg vectors (±1, 0)2π/a 0 and (0, ±1)2π/a 0 marked by black circles. The wavevectors Q ** x ∼ ð3=4,0Þ2π=a 0 and Q ** y ∼ ð0,3=4Þ2π=a 0 from the charge modulation are identified by dashed red and blue circles. e ξ orient /ξ CDW extracted from two different definitions in Refs. 32 (circles) and 33 (diamonds). Ratios of ξ orient /ξ CDW appear >0.5 for all samples, consistent with a striped nature of the charge order. The purple and yellow regions indicate the stripe and checkerboard phases, respectively. x dominates) and blue (σ = À 1, Q ** y dominates) to indicate the local Q ** unidirectional orientations. The unidirectional domains are derived from peaks around Q ** x and Q ** y with details in the text. The R-map has been Fourier-filtered to include only the power spectral density surrounding the four Q ** peaks (dashed circles in Fig. 1d). b, c Fourier transform of red-and blue-colored regions of non-filtered R-map in (a), respectively. Q ** x dominates in the red regions (b), whereas Q ** y in the blue regions (c). Aðq,rÞ = 1 Throughout the paper, we use L = 0.6a 0 for all R-maps. The correlation lengths ξ CDW and ξ orient are then formed from the scalar fields A x ðrÞ = AðQ ** x ,rÞ and A y ðrÞ = AðQ ** y ,rÞ using two different methods as described in Refs. 32,33. Figure 1e summarizes the ratio of ξ orient /ξ CDW obtained from each dataset. In every sample, both methods yield ξ orient /ξ CDW > 0.5, revealing that Bi2201 tends more towards stripes than checkerboards, and that stripe order likely would also be present in the zero disorder limit. Regardless of whether the tendency to stripe modulations survives the hypothetical zero disorder limit, in the Bi2201 samples under consideration, our analyses show that there are local stripe domains present.

Ising domains
Having identified the stripe nature of the local charge modulations in Bi2201, we map out where in the sample there are locally x-oriented domains, and where there are locally y-oriented domains. In Fig. 2, we show this mapping for the UD32K sample, constructed as follows: At each position r, the local FT A(q, r) is calculated according to Eqn. (1), which employs a Gaussian window of width L, with L optimized as described in the Supplementary information. We then integrate the FT intensity in a 2D gaussian window centered on q = ±Q ** x , and divide it by the integrated FT intensity around q = ±Q ** y . If this ratio is greater than a threshold f~1 (i.e. Q ** x is dominant), the region is colored red in Fig. 2a, otherwise the region is colored blue. The pattern thus derived in Fig. 2a is largely insensitive to changes in detail such as the exact center of the integration window, the size of the integration window, We analyze the pattern formation under the assumption that it is driven by a critical point under the superconducting dome. At the critical point of a second order phase transition, a system exhibits correlated fluctuations on all length scales, resulting in power law behavior for measurable quantities, with a different "critical exponent" controlling the power law of each quantity.
If the complex pattern formation shown in Fig. 2 is due to proximity to a critical point, then the critical exponents would be encoded in the geometric pattern, and the quantitative characteristics of the clusters would act like a fingerprint to identify the critical point controlling the pattern formation. This reveals information such as the relative importance of disorder and interactions. Because critical exponents are particularly sensitive to dimension, this analysis can also reveal whether the clusters form only on the surface of the material (like frost on a window), or whether they extend seamlessly from the surface into the bulk (like a tree whose roots reach deep underground). Unless the structures seen on the surface pervade the bulk of the material, they cannot be responsible for the bulk superconductivity.

Critical exponents
Near a critical point, the number of clusters D of a particular size s is power-law distributed, D(s) ∝ s −τ , where s is the number of sites in the cluster and τ is the Fisher critical exponent 34 . Figure 3a shows the first 〈s〉, second 〈s 2 〉, and third 〈s 3 〉 moments of the cluster size distribution as a function of the window size W, where s is the observed area of each cluster. Consistent with a system near criticality, the behavior of the moments vs. window size W displays robust power law behavior. The cluster moments are related to critical exponents by hs n i / W ðn + 1ÀτÞd * v , where d * v is the effective volume fractal dimension. Since the first moment hsi depends only weakly on W (leading to larger error in the estimate of the power law), we combine the information from 〈s 2 〉 and 〈s 3 〉 to derive τ. In the UD32K sample (Fig. 3a), we find τ = 1.97 ± 0.08.
The boundaries of clusters become fractal in the vicinity of a critical point, scaling as H / R d h where H is the size of each cluster's hull (outer perimeter), R is the radius of gyration of each cluster, and d h is the fractal dimension of the hull. The interiors of the clusters also become fractal, scaling as V / R d v , where d v is the volume fractal dimension of the clusters. Since STM probes only the sample surface, the observable quantities are the area A / R d * v and perimeter P / R d * h of each cluster, where d * v and d * h represent the effective volume and hull fractal dimension, respectively. In Fig. 3b and c, the cluster properties A and P are plotted vs. R, revealing a robust power law spanning 2.5 decades for both d We turn now to the orientation-orientation correlation function G orient (r), which is the analog of the spin-spin correlation function familiar from Ising models, G orient (r) = G spin (r) = 〈S(r)S(0)〉 − 〈S(0)〉 2 , where r = |R i − R j | is the distance between (x, y) positions, here measured only on the surface. Near criticality, this function displays power law behavior as G orient ðrÞ / 1=r dÀ2 + η || , where η || is the anomalous dimension as measured at the surface, and d is the physical dimension of the phenomenon being studied, whether d = 2 for a surface phenomenon, or d = 3 for physics arising from the bulk interior of the material. Figure 3d shows G orient (r) for UD32K (triangles). For the UD32K sample as well as the other samples studied, G orient (r) does not have the standard power-law behavior expected near a critical point, but instead it decays more quickly with r.
Whereas the orientation-orientation correlation function is not power law in the data, the pair connectivity function, which is the probability that two aligned regions a distance r apart are in the same connected cluster 35 , does display robust power law behavior in the data, with G conn ðrÞ / r ÀðdÀ2 + η conn Þ , with d − 2 + η conn = 0.29 +/− 0.036, as shown in Fig. 3d (circles).
While the pair connectivity function has been widely discussed for uncorrelated percolation fixed points 35 , where it is a power law, it has not previously been characterized at other fixed points. Our simulations of both the clean and random field Ising models show that the pair connectivity function is also a power law at the 2D clean Ising (C-2D) and the 2D random field Ising (RF-2D) fixed points (see Supplementary information). We find that it also displays power law behavior on interior 2D slices 36 and at a free surface for the 3D clean Ising (C-3D) and 3D random field Ising (RF-3D) fixed points. In addition, our simulations of the clean and random field models close to but not at criticality show that there is a regime in which a short correlation length ξ spin is evident in the spin-spin correlation function, in conjunction with robust power law behavior with a long correlation length ξ cluster in the pair connectivity function, consistent with this dataset (see Supplementary information). Figure 4 shows the experimentally determined critical exponents for five dopings spanning from underdoped to overdoped, using cluster maps based on both Q ** and Q * . We find similar results at all dopings, which also show robust power laws, with the same exponents within error bars as the UD32K sample, and a cluster correlation length ξ cluster which exceeds the FOV (see Supplementary information). Figure 5 shows a comparison between the theoretical critical exponents of Eqn. (2) and the experimentally determined values averaged over both Q** and Q* maps and over all dopings. The data-derived value of d * h is inconsistent with the 2D percolation (P-2D) fixed point, indicating that interactions between stripe orientations must be present. In addition, the data-derived value for d − 2 + η conn is inconsistent with that of the C-2D fixed point, and the data-derived value of d * h is inconsistent with the RF-2D fixed point. The remaining candidate fixed points controlling the power law order of stripe orientations are C-3D and RF-3D (denoted C-3Ds and RF-3Ds, respectively, in the figure because we report theoretical values of the exponents at a free surface of the 3D models). Therefore, we find that the data-derived exponents are consistent with those of a layered clean or random field Ising model with J ⊥ > 0 near criticality.
This shows that the fractal patterns observed here via STM are not confined to the surface, like frost growing on a window pane. Rather, these fractal stripe clusters fill the bulk of the material, more like a tree viewed through a 2D window. In the same way that transverse stripe fluctuations help electron pairs condense into a superconducting state rather than into the competing (insulating) pair crystal phase [37][38][39] , the stripe orientation fluctuations observed here could also have a profound effect on superconductivity, since orientation fluctuations of stripes also frustrate the pair crystal.
The doping independence evident in Fig. 4 is surprising, since one would expect there to be a phase transition from ordered to disordered stripe orientations as doping (a source of quenched disorder) is increased, with the critical, power law behavior observed here limited to the vicinity of the phase transition. While a broad region of critical behavior like that observed here is not natural near the C-3D fixed point, a broad region of critical behavior is characteristic of the RF-3D fixed point: for example, the cluster size distribution D(S) displays 2 decades of scaling, 50% away from the RF-3D critical point 40 . Therefore, the overall phenomenology, taking into account the cluster exponents, the static nature of the correlations, and the broad critical region, is consistent with the critical nematic correlations being controlled by the RF-3D fixed point.

Discussion
While our findings suggest a prominent role for criticality in the phase diagram of cuprate superconductors, the spatial structures reported here are inconsistent with quantum criticality because these correlations are static on the timescales of several seconds, whereas quantum critical correlations fluctuate in time. In addition, for a quantum critical point tuned by doping, quantum critical scaling is confined to a narrow region close to the critical doping, in a "wedge" emanating from the critical doping and extending up in temperature. By contrast, we find critical, power law correlations at low temperature throughout the entire doping range measured. Whereas the lack of detectable doping dependence in finite FOV's is inconsistent with quantum criticality, it is natural in the classical, three-dimensional random field Ising model. What the RF-3D fixed point shares in common with quantum criticality is that it is also a zero temperature critical point. However, it is tuned by disorder rather than by quantum fluctuations.
In addition, RF-3D is notoriously difficult to equilibrate in the vicinity of the critical point, since the relaxation time scales exponentially with the spin-spin correlation length: τ rel ∼ exp½ξ θ spin , where the violation of hyperscaling exponent θ = 1.5. If the spin-spin correlation length reaches even 10 unit cells, the relaxation time will be exp½10 1:5 ≈ 5 × 10 13 times any bare microscopic timescale. Compare this with critical scaling near the C-3D fixed point, where for a correlation length of 10 unit cells, the relaxation time is on the order of τ rel / ξ z spin ≈ 10 times the bare timescale (where z is the dynamical critical exponent, which is of order 1).
As temperature is lowered on any given sample, it falls out of equilibrium if the relaxation time exceeds a timescale t 0 which is set by the cooling protocol, τ rel / exp½ξ θ spin ≳ t 0 . Thus the orientational correlation length depends on the cooling protocol, rather than on doping, when approaching an ordered ground state (i.e. for doping p < p c ), where the correlation length scales as ξ spin ∝ 1/|T − T c | ν . Beyond that doping, the correlation length scales as ξ spin ∝ 1/|p − p c | ν at low temperature.
Our model unifies the observed power law scaling of stripe orientations in this material with the four decades of scaling seen in stripe orientations in NCCOC 19 . The fact that we see scaling all the way out to the FOV among four different critical exponents that are consistent with each other is strong evidence that the effect is due to criticality. If the scaling were to arise from some other mechanism, it is highly unlikely that all four exponents would line up with any theoretical model. Further tests of our model include the following: (1) At larger FOV, the cluster correlation length diverges as the critical doping is approached from large doping, as ξ cluster / |p À p c | Àν cluster . Such a finding would also serve to identify the critical doping concentration for the vestigial nematic 41 . (2) The pair connectivity function will show a scaling collapse 40 as a function of doping in the vicinity of the critical doping. (3) As a function of uniaxial in-plane strain, stripe orientations will experience switching events (avalanches) displaying Barkhausen noise 40 in the vicinity of the critical point, typical of 3D RFIM criticality. (4) Finally, more fully aligned stripe orientations can be trained into the sample by cooling in the presence of in-plane strain 19 .
Our discovery that the charge modulations observed at the surface of Bi 2−z Pb z Sr 2−y La y CuO 6+x are locally one dimensional and also extend throughout the bulk of the material has important implications for the mechanism of superconductivity in cuprate superconductors. The fractal stripe clusters may have a profound effect on superconductivity, by frustrating competing orders like the pair crystal. The cluster analysis framework demonstrated here extends the capability of all surface probes used to study quantum materials to distinguish surface from bulk behavior. Furthermore, our finding that fractal stripe patterns both permeate the bulk of a cuprate superconductor and that they share universal features throughout the superconducting dome, raises important questions. Because doping naturally introduces disorder, a disorder-driven, zero temperature critical point for electronic nematicity is a very real possibility in other cuprates as well. Indeed, the fact that the CO in La 2−x Sr x CuO 4 is relatively unaffected by the onset of superconductivity 31 may indicate that the superconducting dome in LSCO is far from a quantum critical point, as argued in Ref. 42, making it a good candidate to look for disorder-driven criticality as discussed here. More work is needed to further elucidate the connection between fractal electronic textures and superconductivity. For example, the connectivity correlation length of the stripes exceeds the field of view of our experiments throughout the doping range. An important open question for future studies is to establish the relationship between this correlation length and the optimal superconducting transition temperature.

STM measurements
Two different home-built STMs were used to acquire the data in this paper, both in cryogenic ultra-high vacuum. The samples were cleaved in situ at~25 K and inserted immediately into the STM sample stage for imaging at 6 K. A mechanically cut polycrystalline PtIr tip was firstly calibrated in Au single crystals to eliminate large tip anisotropy. To obtain a tunneling current, we applied a bias to the sample while the tip was held at virtual ground. All tunneling spectra, which are proportional to the local density of states at given sample voltage, were measured using a standard lock-in technique.

Theoretical models
Because there are only two orientations of the unidirectional domains, we can map the orientations to an Ising variable 19,43,44 , σ = ±1, where the + sign corresponds to red regions in Fig. 2, and the − sign corresponds to the blue regions. We model the tendency of neighboring unidirectional regions to align by a ferromagnetic interaction within each plane J || as well as an interlayer coupling J ⊥ Any net orienting field, whether applied or intrinsic to the crystal, contributes to the bulk orienting field h 19 . In any given region, the local pattern of quenched disorder breaks the rotational symmetry of the host crystal, corresponding to random field disorder h i in the Ising model. Quenched disorder can also introduce randomness in the couplings J, also known as random bond disorder. In the presence of both random bond and random field disorder, the critical behavior is controlled by the random field fixed point. In the model, h i is chosen from a gaussian distribution of width Δ, which quantifies the disorder strength.

Simulation methods
When comparing to a 2D model, the effective fractal dimensions observed at the surface can be compared directly with those of the model, d When comparing to a 3D model, we have calculated the cluster critical exponents of the model at a free surface, denoted C-3Ds and RF-3Ds in Fig. 5.
For the clean Ising model in three dimensions, because the fixed point (C-3D) controlling the continuous phase transition is at finite temperature, we use Monte Carlo simulations to generate stripe orientation configurations. To calculate the critical exponents of the 3D clean Ising model at a free surface, a 840 × 840 × 840 3D clean Ising model with periodic boundary conditions in the x and y direction and open boundary conditions in the z direction was simulated at T c = 4.51J with 20,000 steps of the parallel Metropolis algorithm. To compare with the finite field of view of the experiments, we average over nine windows of size 256 × 256, taken from a free surface of the final spin configuration. The averages of these critical exponents are shown in Fig. 5. The standard deviations are smaller than the symbol size.
For the random field Ising model in three dimensions, the fixed point (RF-3D) controlling the continuous phase transition is at zero temperature, and we use a mapping to the max-flow min-cut algorithm [45][46][47] to calculate exact ground state spin orientation configurations. The critical point of the 3D random field Ising model occurs at zero temperature. To calculate the 3D random field Ising model surface exponents, ground states were computed for 10 different disorder configurations of a 512 × 512 × 512 3D RFIM with open boundary conditions in the z direction and periodic boundary conditions in the x and y directions with R = 3, using a mapping between RFIM and the max-flow/min-cut problem [45][46][47] . To compare with the finite field of view of the experiments, the top surface of these ground states was windowed to system size 256 × 256 and critical exponents were extracted from the corresponding windows. The averages of these critical exponents are shown in Fig. 5. The standard deviations are smaller than the symbol size.

Cluster methods
While the exponent τ derived from STM data is close to the narrow range allowed by the theoretical models, it is slightly below this range. Estimates of this exponent derived from a finite FOV are known to be skewed toward low values due to a bump in the scaling function, especially in the presence of random field effects 40 . To mitigate this effect, we perform finite-size scaling by analyzing the data as a function of window size W.
To mitigate possible window effects associated with a finite FOV in deriving the fractal dimensions d h and d v , only the internal clusters that touch no edge of the Ising map have been included in the analyses of experimental data as well as simulation results. To extract the effective fractal dimensions, we adopt a standard logarithmic binning technique for analyzing power-law behavior 48 . (See also Supplementary Information and Refs. 51-67 therein.)