Spatial complexity due to bulk electronic nematicity in a superconducting underdoped cuprate

Abstract

Surface probes such as scanning tunnelling microscopy have detected complex electronic patterns at the nanoscale in many high-temperature superconductors. In cuprates, the pattern formation is associated with the pseudogap phase, a precursor to the high-temperature superconducting state. Rotational symmetry breaking of the host crystal in the form of electronic nematicity has recently been proposed as a unifying theme of the pseudogap phase. However, the fundamental physics governing the nanoscale pattern formation has not yet been identifed. Here we introduce a new set of methods for analysing strongly correlated electronic systems, including the effects of both disorder and broken symmetry. We use universal cluster properties extracted from scanning tunnelling microscopy studies of cuprate superconductors to identify the fundamental physics controlling the complex pattern formation. Because of a delicate balance between disorder, interactions, and material anisotropy, we find that the electron nematic is fractal in nature, and that it extends throughout the bulk of the material.

ARTICLE

Received 14 Oct 2011 | Accepted 18 May 2012 | Published 26 Jun 2012 DOI: 10.1038/ncomms1920

Spatial complexity due to bulk electronic nematicity in a superconducting underdoped cuprate

B. Phillabaum1, E.W. Carlson1 & K.A. Dahmen2

Surface probes such as scanning tunnelling microscopy have detected complex electronic patterns at the nanoscale in many high-temperature superconductors. In cuprates, the pattern formation is associated with the pseudogap phase, a precursor to the high-temperature superconducting state. Rotational symmetry breaking of the host crystal in the form of electronic nematicity has recently been proposed as a unifying theme of the pseudogap phase. However, the fundamental physics governing the nanoscale pattern formation has not yet been identifed. Here we introduce a new set of methods for analysing strongly correlated electronic systems, including the effects of both disorder and broken symmetry. We use universal cluster properties extracted from scanning tunnelling microscopy studies of cuprate superconductors to identify the fundamental physics controlling the complex pattern formation. Because of a delicate balance between disorder, interactions, and material anisotropy, we nd that the electron nematic is fractal in nature, and that it extends throughout the bulk of the material.

1 Department of Physics, Purdue University, West Lafayette, Indiana 47907, USA.2 Department of Physics, University of Illinois , Urbana-Champaign , Ilinois

61801 , USA . Correspondence and requests for materials should be addressed to E.W.C (email: [email protected] ) .

1

NATURE COMMUNICATIONS | 3:915 | DOI: 10.1038/ncomms1920 | www.nature.com/naturecommunications

2012 Macmillan Publishers Limited. All rights reserved.

ARTICLE

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms1920

Whereas lanthanum-based cuprate superconductors show striking evidence of electronic liquid crystal phases in the pseudogap regime 1 , the issue is less clear in the

higher transition temperature compounds Yba 2Cu 3O6+x (YBCO) and Bi 2Sr 2CaCu 2O8+ (BSCCO). Recent experiments on YBCO report an electron nematic in the pseudogap regime, detected via

Nernst eect 2, transport3, and neutron scattering 4, and evidence of time-reversal symmetry breaking detected via neutron scattering 5 and the Kerr eect 6. Th e detection via scanning tunnelling microscopy (STM) of a glass of unidirectional domains in BSCCO 79 and in Ca 2x Na x CuO2Cl2 and Bi 2Sr 2Dy x Ca 1x Cu 2O8+ (Dy-Bi2212)10 has now been followed by the dramatic demonstration of a large (>40 nm) electron nematic domain at the surface of BSCCO 11.

However, it is not known whether such structures extend into the bulk of the system, or whether they are merely surface eects in the Bi-based compounds. Because such electronic liquid crystals have been proposed as a unifying theme of the pseudogap phase 12,11 (a precursor to the superconducting phase), it is important to understand to what extent they inuence bulk properties such as superconductivity. In addition, the issue of how to classify broken symmetry states, in the presence of disorder in strongly correlated electronic matter, is important but poorly understood.

Here we introduce new methods for analysing strongly correlated electronic systems including the eects of both disorder and broken symmetry. We use universal cluster techniques to show that the quantitative geometrical properties displayed at the surface of Dy-Bi2212 (ref. 10) indicate that the electron nematic, while not necessarily long-range ordered, possesses large clusters throughout the bulk of the material, with important implications for the coexisting superconductivity.

ResultsDening the model. Th e physics of the orientational degree of freedom of the electron nematic in the cuprates in the presence of quenched disorder maps to a disordered Ising model 13.When an electron nematic forms, there is a preferred orientation to the electronic degrees of freedom, leading to rotational symmetry breaking of the host crystal. We consider Cu O planes that are C4 symmetric, with an incipient electron nematic that breaks the local rotational symmetry of the host crystal from C4 to C2 , leading to two possible nematic orientations. We coarse grain the system, and dene a local nematic order parameter by

= 1,corresponding to the two allowed orientations. Th e tendency for neighbouring nematic regions to align is modelled as a ferromagnetic nearest-neighbour interaction.

Material disorder (such as dopants) competes with the ferromagnetic coupling between local nematic directors. We consider a general model encompassing both disorder in the coupling strengths, as well as random eld disorder that couples linearly to the nematic director.

H J J J J h h

ij

ij

i j

the system. In addition, the local amplitude of the nematic order parameter can vary spatially 11, an eect that may be subsumed into randomness in the bond strengths in an order parameter description. Local energy density disorder has other microscopic realizations, including site dilution, but all types of non-frustrating local energy density disorder belong to the same universality class as random bond disorder 14,15. Frustrating disorder occurs when there are variations in the sign of J, in which case xed points associated with spin glass behaviour can arise. While not forbidden to occur, we have not included this extension of equation 1 because it is physically unlikely in the present system. In equation 1 , dJij|| represents in-plane bond disorder in the coupling strength, and dJij represents bond disorder in the interplane coupling strength.

Th e other class of disorder, random eld disorder, arises because any local pattern of disorder (such as dopant atoms) breaks rotational symmetry, thus favouring one or the other orientation of the nematic director in that region. Th is type of disorder couples linearly to the nematic director. Th e random eld h i is chosen from a gaussian probability distribution centred about zero, with width , which we call random eld strength.

The eld

h=hint+hext represents an orienting eld that breaks rotational symmetry. Th e external contribution hext may be achieved by the application of, for example, magnetic elds, uniaxial pressure, high currents, or other symmetry-breaking external perturbations 13,16 . In the material considered here, data was taken in the absence of applied elds, hext=0. Another source of nite h may be internal crystal eects hint, such as the chains in YBCO may present.

Such issues do not arise in Dy-Bi2212, and so we set hint=0.

Continuous phase transitions and power-law scaling. The equilibrium behaviour of the model in equation 1 is shown in Fig. 1 . Solid regions denote the ordered nematic phase, from two dimensions (2D, yellow region) to three dimensions (3D, orange region). Solid black lines represent continuous phase transitions in which the nematic order parameter rises continuously from zero on entering the ordered phase. Th e blue arrow on each phase transition line points to the solid green circle representing the xed point that determines the universal properties (such as the power-law scaling behaviour discussed below) of that entire phase transition line.

When the system is at a continuous phase transition (solid black lines in Fig. 1 ), its correlation length is innite, and certain physical properties display power law behaviour, with characteristic critical exponents that are determined by the xed point controlling that phase transition. As the system moves away from the critical point, the correlation length decreases, and the power law behaviour (also known as scaling) is only observed for length scales below the correlation length. Th erefore, one measure of proximity to a continuous phase transition is the number of decades of power-law scaling observed in the physical quantities of the system. We show below that certain STM properties of Dy-Bi2212 display power-law scaling, consistent with the system being near a critical point associated with a continuous phase transition, that is, that the physical system is near one of the black phase transition lines in Fig. 1.

Mapping unidirectional domains to Ising variables. Figure 2a shows scanning tunnelling microscopy data on Bi 2Sr 2Dy 0.2Ca 0.8Cu 2O8+

(Dy-Bi2212) from Supplementary Fig. S3 of ref. 10, reported as an R-map, where at each position [arrowrightnosp]

r,

d s s d s s ssi.

Here J|| sets the overall strength of the in-plane ferromagnetic coupling between nearest-neighbour Ising nematic variables, and J represents the overall coupling strength between Ising variables in neighbouring planes. Th e ratio J /J|| is set by the anisotropy of the material. Th e order parameter n=(1/N) i

i describes the degree of orientational order in the system.

Material disorder can disrupt the coupling between local nematic directors. Th ere are two broad classes of disorder that present themselves at the order parameter level: local energy density disorder (which includes random T c disorder), and random eld disorder 14.

Local energy density disorder may arise in the form of, for example, random bond disorder, in which the strength of the ferromagnetic coupling varies from coarse-grained site to coarse-grained site in

= (1 ) (1 ) ( )

+ + +

|| ||

||

ij i j

i

i

(1)(1)

ij

[arrowrightnosp][arrowrightnosp] [arrowrightnosp]

+ is the ratio of the tunnelling current at positive voltage to that at negative voltage 10,17,18. On the basis of the R-map , the presence of local, unidirectional domains of width 4 ao were noted 10,correspondingto the distance between legs of the 4 ao -wide ladders, where ao is the

Cu Cu distance within the Cu O planes. Th ere is also a coexisting local density wave near ao , corresponding to the distance between rungs on the ladders.

R r V I r V I r V

( , ) = ( , )/ ( , )

NATURE COMMUNICATIONS | 3:915 | DOI: 10.1038/ncomms1920 | www.nature.com/naturecommunications

2012 Macmillan Publishers Limited. All rights reserved.

2

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms1920

ARTICLE

a b

C-3D C-3D

RF-3D

RB-3D

P-3D

C-2D C-2D

3D limit 3D limit

Dy-Bi2212 2D limit 2D limit

Temperature

J /J ||

More layered

More isotropic

Temperature

J /J ||

More layered

More isotropic

P-2D

RF-2D

Random field strength Random bond strength

Figure 1 | Equilibrium phase diagrams and xed points. ( a) Random eld Ising models. ( b) Random bond Ising models. Solid regions denote the ordered nematic phase, from 2D (yellow region) to 3D (orange region). Blue arrows represent how effective model parameters change with increasing coarse graining. Phase transitions are denoted by solid black lines. The corresponding critical exponents are determined by the xed points to which the blue arrows point, denoted by solid green circles. The blue square in ( a) denotes the approximate location of the incipient electron nematic in Dy-Bi2212. While it is not long-range ordered, large planar nematic clusters are present throughout the bulk.

Figure 2b shows our mapping of this dataset to Ising variables. To eect the mapping, we have done the following. In any given region, we take a local spatial Fourier transform (FT) of the R-map, and then focus on k-space intensity near 2 /ao. Th ere are two avours to the unidirectional domains, either oriented along the a axis, or along the b axis. We assign a local Ising variable based on the relative weight of the 2

/ao peak in the a direction to that in the b direction. The quantitative geometric properties of the clusters thus derived below are insensitive to changes in details such as the size of the FT window, the Ising lattice spacing, and the threshhold for which a cluster is counted as

= +1 or

= 1. To eect the mapping, we use the full eld of view of Supplementary Fig. S3 of ref. 10. We use a rolling square FT window, of size 1.42 nm 1.42 nm, which corresponds to 16 16 pixels in the original gure. We sum the integrated intensity about the two peaks at Qx =(2 /ao,0), and subtract the integrated intensity about the two peaks at Qy =(0,2

/ao ) using a square integrationwindow of size1.3251.325centredabout the Qx and Qy . If the result of the subtraction is positive (that is, Q x

is dominant), we assign an Ising variable

= +1,coloured orange in Fig. 2b , and if it is negative (that is, Q y is dominant), we assign an Ising variable

= 1,colouredblueinFig.2b.The distance between centres of FT windows is 4 pixels=3.56, which denes the Ising lattice spacing in Fig. 2b . Although the Ising lattice spacing of panel (b) may appear small, the Ising variable at each site is derived by incorporating information from the surrounding 16 Ising lattice sites.

It should be noted that there are many possible ways to relate the local nematic order to specic observable quantities in the material, and the choice of which to use is more of a practical question than one of the fundamental physics of nematics. For example, dierent techniques have been used in ref. 11 to extract a nematic order parameter associated with symmetry breaking within the crystal unit cell. Our method focusses on nematicity associated with the ao periodicity itself and has the advantage that it does not depend sensitively on the phase of the complex Fourier transform. For example, it is not necessary with our method to have detailed information about the exact location of each atom, but the tradeo is that we cannot say anything about the degree of nematicity within each unit cell.

Cluster techniques and extraction of exponents. The extracted cluster maps are reminiscent of cluster patterns observed in numerical studies of equation 1 , consistent with the idea of

mapping disordered electron nematics 13 to disordered Ising models. We apply quantitative cluster analysis methods 1923 from the statistical mechanics of disordered systems to identify the fundamental physics controlling the pattern formation. We trackthegeometricclustersdened as connected sets of nearest-neighbour domains in which the Ising variables are oriented in the same direction. Th e statistics of the shapes and sizes of these geometric clusters can be quantied, and used to identify the cause of the complex pattern formation. In particular, we have discovered that quantitative measures of the cluster shapes as well as the distribution of cluster sizes reveal power-law behaviour over multiple decades.

Several quantitative properties can be extracted from the spatial conguration of the clusters in Fig. 2 , each described in more detail below:

characterizes the distribution of cluster sizes 20 and the spin spin correlation function yields the combination d2+

||

(refs 24,25), where d is the dimension of the phenomenon being studied and

|| istheanomalousdimensionat thesurfaceof the material. In addition, we dene below an eective ratio of the hull fractal dimension 25,26 to the volume fractal dimension 25,26 of the clusters, d h */d v * . Near continuous phase transitions, these physical properties display scaling behaviour with exponents that assume universal values set by the corresponding xed point, allowing us to distinguish which of the xed points (solid green circles in Fig. 1 ) is responsible for the observed pattern formation.

Th e exponent

characterizes the distribution of cluster sizes.

Figure 2b shows the geometric clusters identied from Dy-Bi2212. Figure 3a shows the cluster size distribution, a histogram of cluster sizes. Th e results are noisier at larger cluster sizes, owing to the nite eld of view. Figure 3b shows the same distribution with logarithmic binning, a standard technique for analysing power law behaviour 27. Near a critical point, the cluster size distribution D(A) is of the form D(A)A

, where A is the size of a cluster. Using a straightforward t to this form, we nd

=1.710.07. Th ere is one large spanning cluster, represented by the last point in Fig. 3b , which is excluded from the t. Not including the spanning cluster, 2.5 decades of scaling are evident. In this case, the spanning cluster also lies near the regime of scaling. If the spanning cluster is also included, the data display 4 decades of scaling. However, a larger eld of view is necessary to determine whether the spanning cluster is truly in the regime of scaling.

Cluster volumes and surfaces (hulls) become fractal near certain critical points, in the sense that the volume V and hull H scale

NATURE COMMUNICATIONS | 3:915 | DOI: 10.1038/ncomms1920 | www.nature.com/naturecommunications

2012 Macmillan Publishers Limited. All rights reserved.

3

ARTICLE

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms1920

a b

1.41.21.00.8

c

d

Figure 2 | Mapping of Ising nematic domains from tunnelling asymmetry maps in Dy-Bi2212. ( a) A subset of Supplementary Fig. S3 of ref. 10 of STM on Dy-Bi2212, showing the R-map taken at 0.15 V a function of position [arrowrightnosp]

r , where

[arrowrightnosp][arrowrightnosp] [arrowrightnosp]

+ is the ratio of the tunnelling current at positive voltage to that at negative voltage 10,17,18. In the colour scale, brighter corresponds to larger R. (Reproduced with permission from AAAS.) ( b) Ising domains derived from peaks at Q x and Q y . Notice the large, percolating cluster (orange), with isolated ipped domains inside (blue). ( c) The image from( a), masked by the Ising map ( b) so as to show only the blue domains. ( d) The same image, masked by the Ising map ( b) so as to show only the orange

domains. The sample has superconducting transition temperature T c ~ 45 K, and data were taken at T=4.2 (ref. 10).

R r V I r V I r V

( , )= ( , )/ ( , )

like a fractional power of the length scale of the cluster. Using, for example, the radius of gyration R (ref. 25 as a measure of the length scale of a cluster, this implies V Rdv

and H Rdh

. Combining the two relations to eliminate R , one obtains H Vdh dv

/ . Because STM is a surface probe, the available information is the observed area A Rdv

* and the observed perimeter P Rdh

* of each cluster, and

/ and

. Larger elds of view will enable a more reliable determination of this quantity.

Fixed points of the model. Having extracted quantitative measures of the critical exponents that are available from geometric clusters, we are now in a position to decide which of the candidate xed points of Fig. 1 (if any) is controlling the observed pattern formation and is therefore responsible for the observed power laws. We begin by reviewing the xed points, then compare observations with theory. Th e equilibrium phase diagrams and corresponding xed points of the random bond and random eld Ising models are shown in Fig. 1 . In the gure, solid black lines denote phase transitions between an ordered electron nematic and a disordered phase. Th e exponents of each phase transition line are set by the xed point (solid green circle) to which the blue arrow points. In the clean limit of zero random eld strength and no bond disorder, the model has a nite-temperature continuous phase transition from a disordered phase to an ordered electron nematic, that is, with long-range orientational order n0. Th e universal properties of that transition are controlled by the xed point labelled C-2D in strictly 2D systems, and controlled by the C-3D xed point in layered or 3D systems.

In 2D, long-range orientational order is forbidden at any nite random eld strength ( Fig. 1a ), and universal properties near this

therefore P Adh dv

*/ * . For a 2D model, the observed fractal dimensions are the true fractal dimensions of the model, d d

v v

* = and

d d

h h

* = , and therefore P A A

dh dv dh dv

*/ * / . When comparing with a 3D model, we make the reasonable assumption that, at the surface, one is observing a 2D cross-section of a 3D cluster. In this case, the observed fractal dimensions are related to the true fractal dimensions of the model by geometric factors, d d

v v

* = 2 /3 and

d d

h h

* = /2,

yielding P A A

dh dv dh dv

*/ * 3 /4 (see Methods). Using a straightforward t over the 2.5 decades of power-law scaling of P Adh dv

* / * gives an eective fractal dimension ratio of d d

h v

* *

/ = 0.78 0.01

, excluding the spanning cluster. It is evident from Fig. 3c that the spanning cluster is also near the scaling regime for this measure. This extraction of the exponent d d

h v

* *

/ can be considered quite reliable because multiple decades of scaling are observed.

Within the Ising description, the spin spin correlation function of the Ising pseudospin variables is G(r)=S(r)S(0) S(r)S(0)|r|( d 2+ ) . When measured for two points on the surface of a material, G(r)|r|( d 2+ ||). In the eld of view available, we have averaged over all sites to obtain the spin spin correlation function plotted in Fig. 3d . For this measure, scaling is expected at long distances, whereas the correlation function at short distances

can be dominated by nonuniversal eects such as the spanning cluster in this case. We therefore exclude the rst three points in Fig. 3d , whereupon a straightforward t yields d2+

||=0.80.3. However, with less than one decade of scaling, this measure is less reliable than our values for d d

h v

* *

4

NATURE COMMUNICATIONS | 3:915 | DOI: 10.1038/ncomms1920 | www.nature.com/naturecommunications

2012 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms1920

ARTICLE

a

b

2.5

2

2.0

1.5

1

Log 10D(A)

Log 10G(r)Log 10D(A)

0

1

1.0

0.5

0.0

4.0

2

3

4

0

1

2

3 4

0 1 2 3 4

Log10A

Log10A

c d

0.0

0.5

Log 10perimeter

3.53.02.53.01.51.00.5

1.0

1.5

2.0

2.5

3.0

0 1 2 3 4

0.5 1.0

0.0 1.5

Log10area

Log10r

Figure 3 | Cluster size distribution and critical exponents. ( a) Raw cluster size distribution. ( b) Cluster size distribution after logarithmic binning 27, used to calculate the critical exponent

. (c) The effective fractal dimension ratio, d h */d v * , relates the perimeter of clusters to their area. Logarithmic binning has been used. ( d) Logarithmically binned data for calculating the critical exponent d 2+ || .

limit are set by the RF-2D xed point. In 3D, electron nematic order is allowed below the critical disorder strength of the random eld, cD T J

3 ( = 0) = 2.27 [parallelto] (ref. 19 ). In strongly layered systems, the critical disorder strength is nite, but signicantly reduced from the corresponding value in 3D systems 28 . Transitions in 3D and layered systems are controlled by the RF-3D xed point. Because the critical region is large in the random eld Ising model, it is possible to observe the scaling in a broad range of parameters near the critical point 29. For example, two decades of scaling appear in the avalanche size distribution of the nonequilibrium 3D RFIM even at a disorder strength of =4 J, which is 85 % larger than the nonequilibrium critical disorder strength cD J

3 = 2.16 (ref. 29).

In contrast with random eld models, weak random bond disorder ( Fig. 1b ) does not forbid nematic order in 2D. Furthermore, random bond disorder is irrelevant in the renormalization group sense in 2D, and the phase transition is therefore governed by clean Ising model exponents set by the C-2D xed point. In layered and 3D systems, the phase transition of the random bond model is controlled by the disordered xed point, which we have labelled RB-3D, also known in the literature as R 15. In the presence of both random bond and random eld disorder, the universality class is that of the random eld model.

Which xed point is the system nearest? We trackthegeometric clusters dened as connected sets of nearest-neighbour domains in which the Ising variables are oriented in the same direction. Although power law behaviour is generically associated with a continuous phase transition, the geometric clusters do not always display power law behaviour at the continuous Ising nematic to disordered phase transition 30. Th e geometric clusters do exhibit power law behaviour at the percolation xed points 25 (both P-2D and P-3D), at the random eld xed points 19,31 (both RF-2D and RF-3D), and at the 2D clean Ising xed point (C-2D) 32 . At the other xed points, including the 3D clean Ising xed point (C-3D) 32 and the 3D random bond xed point (RB-3D) 15 , the geometric clusters do not display power law behaviour 30 . Rather, it is the Fortuin-Kasteleyn (FK) clusters that exhibit power law behaviour at the C-3D and

RB-3D xed points. FK clusters are related to the geometric clusters we study here by assigning a temperature-dependent probability of breaking each bond in the geometric cluster. The new (smaller) connected clusters thus generated are the FK clusters. We focus on the geometric clusters, because they can be directly extracted from the data without need of further ansatzes, and it is the geometric clusters themselves that exhibit power law behaviour in the data.

At the 3D clean Ising xed point (C-3D) 32 and the 3D random bond xed point (RB-3D) 15 , the geometric clusters do not display power law behaviour 30 . Rather, they display power law behaviour at the temperature at which minority spin clusters rst percolate, T p , which is less than the temperature T c at which the nematic-to-disordered phase transition takes place, T p <T c (ref. 32). In these cases, the power law behaviour of the geometric clusters is controlled not by the phase transition at T c , but by the 3D percolation xed point, P-3D. Because the geometric clusters that we track do not exhibit power law behaviour at the C-3D and RB-3D xed points, these two xed points cannot be the cause of the observed power law cluster size distribution, and so they are excluded from the comparison in Fig. 4 .

Th e 3D percolation xed point can also be ruled out as the cause of the power law behaviour, whether as merely an uncorrelated percolation phenomenon, or at T p <T c in the 3D clean and random bond Ising cases. In the uncorrelated case, the percolation thresh-hold is p c =0.311in3D.Near thisratioof verticalto horizontal

nematic patches, there would be a detectable net nematicity in the system, which has not been observed in bulk measurements. Therefore, we can rule out P-3D in the uncorrelated case. In the correlated cases of the 3D clean and random bond Ising models, although clusters do exhibit power law behaviour (controlled by the P-3D xed point) at T p <T c , this happens deep in the ordered phase. We can thus rule out these cases as well, because bulk experiments to date have not shown the material to possess a 3D-ordered nematic phase. Th erefore, we also exclude P-3D from the comparison in Fig. 4 .

Th e percolation threshold on a square lattice is p c =0.59.If the power law behaviour were due to an uncorrelated purely geometric percolation phenomenon occurring independently within each

5

NATURE COMMUNICATIONS | 3:915 | DOI: 10.1038/ncomms1920 | www.nature.com/naturecommunications

2012 Macmillan Publishers Limited. All rights reserved.

ARTICLE

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms1920

a b

Clean lsing model

Percolation model

Random field lsing model

Percolation model

Random field lsing model

indicating that large planar nematic clusters permeate the bulk of the system. Ultimately at long length scales, such systems ow to the universality class of the higher dimension. Our model therefore predicts that the measured value of d d

h v

* *

2.5 1.0

Clean lsing model

2.01.51.00.5 0

0.90.80.70.6

0.5 C-2D P-2D RF-3D

RF-2D

/ in this material will dri towards the 3D random eld (RF-3D) exponent as larger elds of view are measured.

Th e value of d2+

|| also shows broad variation among the theoretical xed points considered in Fig. 4 , and can be a good measure for distinguishing among the candidate xed points. As with d d

h v

* *

/ ,

C-2D P-2D RF-3D

RF-2D

Figure 4 | Critical exponents. ( a) Critical exponents

the value of d2+

|| extracted from experiment is consistent with a layered random eld Ising model. It is also consistent with a 2D random eld Ising model within the error bars, whereas it is inconsistent with the other candidate xed points. Our model predicts that the measured value of d2+

|| will dri closer to the RF-3D value with larger elds of view.

As can be seen from Fig. 4 , there is a narrow range of values of displayed by the xed points of equation 1 , and therefore this exponent is not a very good way to distinguish among the xed points. One reason for the narrow range of

is that there is a constraint on

and d2+

||.

/ . The horizontal lines in red, green, and purple are our results for the spatial characteristics of clusters in Dy-Bi2212as reported in Supplementary Fig. S3 of ref. 10. Solid circles represent theoretical values reported in the literature for the xed points of equation 1 , summarized in Supplementary Tables S1-S6 . Red circles represent

,

green circles represent d2+

||, and purple circles represent d d

h v

* *

( b) Critical exponent d d

h v

* *

/ . When

/ = / , shown by the solid circles. When comparing with 3D models, we have assumed that, at the surface, one observes a 2D cross-section of a cluster embedded in 3D, which implies d d d d

h v h v

* *

/ =3 /(4 ) . This value is represented by the open purple circle. Thick lines connecting symbols represent putative crossover of exponents from 2D to 3D behaviour in a layered material.

comparing with 2D models, d d d d

* *

h v h v

this exponent: 2<

plane, then one would expect there to be, on average, an equal numberof verticalandhorizontalnematicpatches,corresponding to a percolation density of p=0.5. While this is not equal to the percolation threshold of p c =0.59 at P-2D, it is close enough that we cannot a priori rule out P-2D as the source of the power law behaviour, and we include it in Fig. 4 . Th ere remain only four xed points to consider: C-2D, P-2D, RF-2D and RF-3D. Exponents characterizing these xed points are charted in Fig. 4 , and compared with the observations in Dy-Bi2212.

Th e material of interest, Dy-Bi2212, is a strongly layered material, expected to have weak coupling between planes. In a layered system, the measured critical exponents can be observed to dri from the 2D limit to the 3D limit as observations are extended to larger length scales. We do not consider the possibility of a dri in exponents from the P-2D to the P-3D xed point, as an intermediate percolation dimension has no physical meaning in a crystal. However, it is possible to observe a dri from the RF-2D xed point to the RF-3D xed point. This dri is controlled by the ratio of the interlayer coupling to the in-plane coupling J /J|| in equation 1.

Th erefore, in Fig. 4 , we have included the possibility of this dri in exponents from the RF-2D to the RF-3D xed point, as indicated by the thick lines connecting xed point circles in the gure.

Discussion

Our results using published STM data 10 on Dy-Bi2212, a material that displays evidence of local electron nematic behaviour at the surface 10 , reveal a strong power law behaviour in these measures, spanning at least 2.5 decades for both

and d d

h v

* *

<3 (ref. 33). Our extracted value of

=1.710.07 is 15 % lower than the models shown, and, furthermore, it violates the above constraint. However, the scaling function of the cluster size distribution has a bump in the random eld case 29 , as well as the clean 2D case, which at small elds of view skews the estimate of

to lower values than the true value, as observed here as well as in numerical simulations 34 . Although 2.5 decades of scaling is remarkable, a larger eld of view is necessary to get a more accurate measure of this exponent because of the form of the scaling function itself. However, because of the narrow range of

displayed by the range of xed points available to equation 1 , it will be diffi cult to distinguish among the various xed points based on the value of

,

even with larger elds of view.

We thus conclude that the electron nematic in this material is fractal 35 , because it is near criticality owing to a competition between disorder and interactions. However, because long-range-ordered electronic nematicity has not been detected in Dy-Bi2212 via a bulk probe, the material is likely in the disordered phase. It is also possible that the material is out of equilibrium 16 owing to glassy behaviour induced by disorder, or even that the broken symmetry is suffi ciently weak that it has thus far escaped detection. Regardless, the evidence does not point to the bulk of the material being in the ordered nematic phase, and it is consistent with the bulk of the material being in the disordered phase. The quantitative properties studied here indicate that the material must then be just on the disordered side, in the regime of intermediate random eld disorder 28, where the disorder strength is small within a plane, but strong between planes, J||[greatermuch][greatermuch]J , as indicated by the blue square in Fig. 1a . Although long-range electron nematic order is not present in this regime, the quantitative characteristics of the fractal geometry and other cluster properties studied here indicate that the nematic behaviour extends throughout the bulk of the material, and clusters within each plane have a long correlation length 28, with potentially important implications for the superconductivity arising in these materials.

Although long-range-ordered nematicity has not yet been observed in bismuth-based copper oxide superconductors, we conclude here that large (>400 ) nematic clusters persist into the bulk of the crystal. Th ese clusters are suffi ciently larger than the superconducting coherence length (

coh~10) 36 to support a theory of high-temperature superconductivity based on electronic liquid crystals such as nematics 12 . In such a theory, the long, straight stripes (or ladders) that make up the nematogens are favourable for superconducting pairing 37. However, in a quasi-1D theory, superconducting pairing also leads to a charge density wave (CDW) instability, which is a tendency for the pairs to form a density wave along the long direction of the stripes or ladders. In a system of long-range-ordered

/ , as

shown in Fig. 3 . On the basis of these power laws, the extracted measures are as follows:

= 1.71 0.07,

d d

h v

* *

/ = 0.78 0.01

, and

d 2 +

|| = 0.8 0.3. Fig. 4 charts a direct comparison between our values extracted from data, and the theoretical exponents of disordered Ising models that are consistent with power-law scaling in the cluster properties.

The eective fractal dimension ratio, d d

h v

* *

/ , shows broad variation among the xed points shown in Fig. 4 and, therefore, it is a good way to distinguish among the xed points. Th e value of d d

h v

* *

/ extracted from experiment falls within the theoretical RFIM values, whereas it is inconsistent with the other candidate models. The value of d d

h v

* *

/ lies between the 2D and 3D values of the RFIM, consistent with the behaviour of a layered system (as expected for Dy-Bi2212),

6

NATURE COMMUNICATIONS | 3:915 | DOI: 10.1038/ncomms1920 | www.nature.com/naturecommunications

2012 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms1920

ARTICLE

stripes, although the superconducting pairing is large, ultimately Coulomb repulsion between pairs causes them to crystallize into a high-temperature insulator (the CDW state). On the other hand, if the stripes uctuate enough (either in time or in space), the CDW becomes so disrupted by defects that superconductivity becomes the true ground state 12 . Within a stripes-based theory of high-temperature superconductivity, it is in this sense that long-range-ordered stripes compete with superconductivity. Conversely, shorter, disordered stripes such as those discussed here encourage pairing to a lesser degree, but are better for phase coherence of the superconductivity 8,12,36.

In summary, we have introduced new methods for analysing strongly correlated electronic systems including the eects of both disorder and broken symmetry, using universal cluster techniques from the eld of disordered statistical mechanics. We have used this new method of analysis to show that there is robust power-law scaling of the nematic clusters at the surface of Dy-Bi2212, extending over the entire eld of view. Th is remarkable property is generically associated with critical behaviour. Of the possible models that may explain such power law behaviour in the long-range behaviour of a discrete nematic, the values are consistent with the material being near a phase transition in a layered random eld Ising model, in the regime of intermediate disorder strength of weak disorder within each plane, but strong disorder between planes. In this regime, planar nematic clusters large enough to support a stripes-based mechanism of high-temperature superconductivity are present throughout the bulk of the material, whether or not true long-range nematic order is present.

Methods

Theoretical values of critical exponents. Th e cluster critical exponents we require to compare with known theoretical results of Ising models are not necessarily directly reported in the literature . However, scaling relations allow us to derive the values of the needed exponents from other exponents available in the literature (See Supplementary Tables S1S6). Th is method was used to derive

implies that the Harris criterion reduces to <0. In the 3D clean Ising model,

~ 0.1 (ref. 14), and randomness in the local energy density is relevant. For the 3D case with weak bond disorder, there is a disordered xed point with new exponents, labelled RB-3D in Fig. 1b . In Fig. 4 and

Supplementary Tables S1-S6 , the 3D random bond Ising model exponents quoted are those associated with the xed point labelled RB-3D in Fig. 1 , also knownas R in the literature. In layered systems (which interpolate between the 2D and 3D limits via the coupling ratio J /J|| from equation 1 ), increasing dimension is relevant, so although the ordered nematic phase occupies less and less of the phase diagram as the 2D limit is approached, the critical exponents of the transitionare those of the 3D model for any nite coupling between planes J /J|| >0. In the 2D clean Ising model,

,

d2+

, and

=0 and such randomness is marginal. In the presence of weak bond disorder, it has been shown that the 2D system ows towards the clean model 46,47. For this reason, the 2D random bond Ising exponents quoted in Supplementary Tables S1 S6 are those of the 2D clean Ising model, controlled by the xed point labelled C-2D in Fig. 1.

The Harris criterion does not apply to random field disorder. However, the equilibrium phase diagram and fixed points in two and three dimensions are known 24,19 and summarized in Fig. 1 . In three dimensions, there is a finite critical disorder strength, which at zero temperature goes to cD J

|| . In all cases, d h is quoted directly from the literature. Because geometric clusters do not exhibit power law behaviour at the C-3D and RB-3D xed points, clusters are not fractal at these xed points, and have no well-dened d h .

To derive

, we use the scaling relation 14

d2+

3 = 2.27 ||19.

In strongly layered systems, the critical disorder strength is finite, but significantly reduced from the 3D value 28. Phase transitions in these cases are controlled by the fixed point labelled RF-3D in Fig. 1. In two dimensions,the critical disorder strength is zero, and for any finite disorder strength, the nematic is forbidden at all temperatures. However, there is an unstable fixed point at T = 0 and zero disorder, RF-2D24. Supplementary Table S1 and Fig. 4 report values from the literature for this fixed point, as it may affect the scaling in some regimes.

References

1. Tranquada, J. et al. Quantum magnetic excitations from stripes in copper oxide superconductors.Nature 429, 534538 (2004).

2. Daou, R. et al. Broken rotational symmetry in the pseudogap phase of a high-T c superconductor. Nature 463, 519522 (2010).

3. Ando, Y., Segawa, K., Komiya, S. & Lavrov, A. Electrical resistivity anisotropy from self-organized one dimensionality in high-temperature superconductors . Phys. Rev. Lett. 88, 137005 (2002).

4. Hinkov, V. et al. Electronic liquid crystal state in the high-temperature superconductor YBa2Cu3O6. 45 . Science 319, 597600 (2008).5. Fauqu, B. et al. Magnetic order in the pseudogap phase of High-T c superconductors.Phys. Rev. Let 96, 197001 (2006).

6. Xia, J. et al. Polar kerr-eect measurements of the high-temperature Yba 2Cu 3O6+ x superconductor: Evidence for broken symmetry near the pseudogap temperature.Phys. Rev. Lett. 100, 127002 (2008).

7. Howald, C., Eisaki, H., Kaneko, N. & Kapitulnik, A. Coexistence of periodic modulation of quasiparticle states and superconductivity in Bi 2Sr 2CaCu 2O8+ .

Proc. Natl Acad. Sci. USA 100, 97059709 (2003).8. Kivelson, S. et al. How to detect uctuating order in the high-temperature superconductors.Rev. Mod. Phys. 7, 12011241 (2003).

9. Robertson, J. A., Kivelson, S. A., Fradkin, E., Fang, A. C. & Kapitulnik, A.

Distinguishing patterns of charge order: Stripes or checkerboards . Phys. Rev. B 74, 134507 (2006).10. Kohsaka, Y. et al. An intrinsic bond-centered electronic glass with unidirectional domains in underdoped cuprates . Science 315, 13801385 (2007).

=(2

)/(

+

)+1=(2

)/

(2

)+1. Th e anomalous dimension

is derived through d2+

=2

/v

(ref. 38). Th e surface critical exponent

|| may be derived through d2+

||=2

1/v

(refs 39,40). For 2D models,

|| =

, whereas for 3D models,

||

.

As discussed in the Results section, cluster volumes and surfaces (hulls) become fractal near certain critical points, in the sense that the volume V and hull H scale like a fractional power of the length scale R of the cluster, V Rdv

and

. (For

R, one may use, for example, the radius of gyration 25). In comparing with a surface probe such as STM, the observed eective fractal dimensions dv

* and

H Rdh

dh

* may be extracted from the relation of the observed area A of each cluster to its perimeter P as: P Adh dv

* / * . For 2D models, a direct comparison can be made, and d d d d

h v h v

* *

/ = / .

Because a surface probe only provides 2D information, when comparing STM or any other surface probe to 3D models, we assume that the observation is about a 2D cross-section of putative 3D clusters. Th e hull fractal dimension of an object embedded in d-dimensional Euclidian space, when generated by an isotropic model, scales as d h (d 1). For the case of a fractal cluster embedded in 3D

Euclidian space, the fractal dimension of the cluster hull (that is, its surface) must be a number between 2 and 3, 2 d h 3. For a fractal cluster that is generated in an isotropic manner (for example, from an isotropic 3D model with J J|| in equation 1 ), taking a 2D cross-section of the cluster results in a new surface (the surface of the 2D cross-section) with eective hull fractal dimension d d

h h

* = /2, with the constraint that d h / 2 must be a number between 1 and 3 / 2.

Likewise, the volume fractal dimension of an object embedded in d-dimensional Euclidian space, when generated by an isotropic model, scales as d v d.

For a fractal cluster that is generated in an isotropic manner (for example, from an isotropic 3D model with J J|| in equation 1 ), taking a 2D cross-section of the cluster volume results in a cross-sectional area, with eective volume fractal dimension d d

v v

* = 2 /3 . One can verify that in the limit of compact clusters, d v 3, and the correct result dv* 2

is obtained. Using the eective fractal dimensions as dened above, when comparing with a 3D model, we have that

P A A

dh dv dh dv

* / * 3 /4 .

Theory of surface criticality. Ultimately, critical behaviour observed at a surface is controlled by the theory of surface critical phenomena. A transition taking place

only on the surface corresponds directly to the 2D xed points we have considered. On the other hand, surface criticality can also arise as the bulk orders, giving rise to surface critical exponents at the bulk transition 39,40. In the case of an ordinary surface transition, the bulk orders without a pre-existing surface transition, and d2+

||=2.54,2.59 and 0.336 for the clean41,42, random bond43, and random eld cases 44, respectively. Here

|| denotes a correlation function measured solely at the surface, applicable to our case. For this reason, d2+

|| may be directly compared

with our value of 0.8 0.3.

Th is is consistent with an ordinary surface transition, in which the exponent is dri ing from from RF-2D to RF-3D. However, further theoretical developments are needed to predict the surface exponents corresponding to the cluster size distribution exponent

and the fractal dimension of cluster surfaces d h for the models discussed here, and also to predict the correlation function surface critical exponent

|| in the case of extraordinary transitions 45, where the bulk orders at a lower temperature than the surface. However, the presence of random eld disorder in the system precludes an extraordinary transition in this material, because a 2D Ising system cannot order in the presence of random eld disorder.

Harris criterion considerations. Th e Harris criterion states that local energy density disorder (such as random T c disorder and random bond disorder) is irrelevant if dv>2, where the exponents refer to the clean model. If disorder is irrelevant, then disorder will reduce the overall transition temperature of the phase transition,but the critical exponents associated with the transition remain those of the clean model (that is, with no disorder). In the presence of hyperscaling (obeyed by the clean and random bond cases), dv=2

7

NATURE COMMUNICATIONS | 3:915 | DOI: 10.1038/ncomms1920 | www.nature.com/naturecommunications

2012 Macmillan Publishers Limited. All rights reserved.

ARTICLE

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms1920

11. Lawler, M. J. et al. Intra-unit-cell electronic nematicity of the high-T c copper-oxide pseudogap states . Nature 466, 347351 (2010).

12.Kivelson, S. A., Fradkin, E. & Emery, V. J. Electronic liquid-crystal phases of a doped mott insulator . Nature 393, 550553 (1998).

13.Carlson, E., Dahmen, K., Fradkin, E. & Kivelson, S. Hysteresis and noise from electronic nematicity in high-temperature superconductors . Phys. Rev. Lett. 96, 097003 (2006).

14. Cardy, J. Scaling and Renormalization in Statistical Physics (Cambridge University Press, 1996) Chapter 8.

15. Berche, P. E., Chatelain, C., Berche, B. & Janke, W. Bond dilution in the 3d ising model: a monte carlo study . Eur. Phys. J. B 38, 463474 (2004).

16.Carlson, E. W. & Dahmen, K. A. Using disorder to detect locally ordered electron nematics via hysteresis . Nat. Commun.. 2, 379 (2011).

17. Anderson, P. & Ong, N. Th eory of asymmetric tunneling in the cuprate superconductors. J. Phys. Chem. Solids 67, 15 (2006).

18. Randeria, M., Sensarma, R., Trivedi, N. & Zhang, F. C. Particle-hole asymmetry in doped mott insulators: implications for tunneling and photoemission spectroscopies.Phys. Rev. Lett. 95, 137001 (2005).

19. Middleton, A. & Fisher, D. Three-dimensional random-eld ising magnet: interfaces, scaling, and the nature of states . Phys. Rev. B 65, 134411 (2002).

20.Perkovic, O., Dahmen, K. & Sethna, J. Disorder-induced critical phenomenain hysteresis: numerical scaling in three and higher dimensions . Phys. Rev. B 59, 61066119 (1999).

21.Liu, Y. & Dahmen, K. Unexpected universality in static and dynamic avalanches.Phys. Rev. E 79, 061124 (2009).

22. Liu, Y. & Dahmen, K. Random-eld ising model in and out of equilibrium . Europhys. Lett. 86, 56003 (2009).

23.Prez-Reche, F. & Vives, E. Finite-size scaling analysis of the avalanches in the three-dimensional gaussian random-eld ising model with metastable dynamics.Phys. Rev. B 67, 134421 (2003).

24. Bray, A. & Moore, M. Scaling theory of the random-eld ising model . J. Phys. C Solid State Phys, 18, L927L933 (1985).

25. Stauer, D. & Aharony, A. Introduction to Percolation Theory (Taylor & Francis, 1991).

26. Mandelbrot, B. B. Th e Fractal Geometry of Nature (W.H. Freeman, 1982).27.Newman, M. E. J. Power laws, pareto distributions, and zipfs law. Contemporary Physics 46, 323351 (2005).

28.Zachar, O. & Zaliznyak, I. Dimensional crossover and charge order in half-doped manganites and cobaltites . Phys. Rev. Lett. 91, 036401 (2003).

29. Perkovic, O., Dahmen, K. & Sethna, J. Avalanches, barkhausen noise, and plain old criticality.Phys. Rev. Lett. 75, 45284531 (1995).

30. Dotsenko, V. et al. Self-avoiding surfaces in the 3d ising model.Nucl. Phys. B

448, 577620 (1995).

31.Seppl, E., Petj, V. & Alava, M. Disorder, order, and domain wall roughening in the two-dimensional random eld ising model.Phys. Rev. E 58, 52175220 (1998).

32. Coniglio, A., Nappi, C., Peruggi, F. & Russo, L. Percolation points and critical point in the ising model . J. Phys. A Math. Gen. 10, 205218 (1977).

33. Stauer, D. Scaling theory of percolation clusters.Phys. Rep. 54, 174 (1979).

34.Alava, M. & Rieger, H. Chaos in the random eld ising model.Phys. Rev. E 58,

42844287 (1998).

35. Fratini, M. et al. Scale-free structural organization of oxygen interstitials in La2CuO4+y.Nature 466, 841844 (2010).

36. Carlson, E. W., Emery, V. J., Kivelson, S. A. & Orgad, D. In: Th e Physics of Superconductors Vol. II: (eds Ketterson, J. & Benneman, K.) 275452 (Springer-Verlag, 2004).

37.Emery, V. J., Kivelson, S. A. & Zachar, O. Spin-gap proximity eect mechanism of high-temperature superconductivity.Phys. Rev. B 56, 61206147 (1997).38. Chaikin, P. & Lubensky, T. Principles of Condensed Matter Physics (Cambridge University Press, 1995).

39. Binder, K. In: Phase Transition and Critical Phenomena vol. 8: (eds Domb, C. & Lebowitz, J. L.) 1144 (Academic Press, 1983).

40. Diehl, H. W. In: Phase Transition and Critical Phenomena vol. 10: (eds Domb,C. & Lebowitz, J. L.) 75267 (Academic Press, 1986).41.Diehl, H. & Shpot, M. Surface critical behavior in xed dimensions d<4: nonanalyticity of critical surface enhancement and massive eld theory approach.Phys. Rev. Lett. 73, 34313434 (1994).

42.Pleimling, M. Critical phenomena at perfect and non-perfect surfaces.J. Phys. A Math. Gen. 37, R79R115 (2004).

43.Usatenko, Z., Shpot, M. & Hu, C. K. Surface critical behavior of random systems: Ordinary transition . Phys. Rev. E 63, 056102 (2001).

44.Laurson, L. & Alava, M. Surface criticality in random eld magnets.Phys. Rev. B. 72, 214416 (2005).45.Brown, S., Fradkin, E. & Kivelson, S. Surface pinning of uctuating charge order: an extraordinary surface phase transition . Phys. Rev. B 71, 224512 (2005).

46.Shalaev, B. Critical behavior of the two-dimensional ising model with random bonds.Phys. Rep. 237, 129188 (1994).

47.Picco, M., Honecker, A. & Pujol, P. Strong disorder xed points in the two-dimensional random-bond ising model . J. Stat. Mech. Theor. Exp. 2006, P09006 (2006).

Acknowledgements

We thank B. Brinkman, E. Fradkin, E. Goins, J. Homan, S. Kivelson, Y. Loh, E. Main, H. Nakanishi, and E. Sepp l for conversations. B.P. and E.W.C. acknowledge support from Research Corporation for Science Advancement and NSF Grant No. DMR 11-06187. K.A.D. acknowledges support from NSF Grant No. DMR 10-05209 and from the Materials Computation Center [NSF DMR 03-25939 (MCC)]. E.W.C. thanks cole Sup rieure de Physique et de Chimie Industrielles (ESPCI) for hospitality.

Author contributions

E.W.C. and K.A.D. conceived the research. B.P. carried out the calculations and analyses. E.W.C. prepared the manuscript. Th e manuscript reects the contribution of all authors.

Additional information

Supplementary Information accompanies this paper at http://www.nature.com/ naturecommunications

Competing nancial interests: Th e authors declare no competing nancial interests.

Reprints and permission information is available online at http://npg.nature.com/ reprintsandpermissions/

How to cite this article: Phillabaum, B. et al. Spatial complexity due to bulk electronic nematicity in a superconducting underdoped cuprate. Nat. Commun. 3:915 doi: 10.1038/ncomms1920 (2012).

8

NATURE COMMUNICATIONS | 3:915 | DOI: 10.1038/ncomms1920 | www.nature.com/naturecommunications

2012 Macmillan Publishers Limited. All rights reserved.