text:
# 63Cu NQR Study of the Inhomogeneous Electronic State in La2-xSrxCuO4

## I Introduction

NQR (nuclear quadrupole resonance) and NMR (nuclear magnetic resonance) have proven to give unique information regarding the inhomogeneous electronic state in the CuO2 plane of various high $`T_c`$ cuprates yoshimura ; tou ; hammel ; stat ; fujiyama ; hunt ; singer ; curro ; haase ; julien ; hunt2 ; singerprl ; tobe  including La2-x\[Sr,Ba\]xCuO4, \[La,Nd,Eu\]2-x\[Sr,Ba\]xCuO4 and La2CuO4+δ. In these materials, the high $`T_c`$ superconductivity is achieved by doping holes into the CuO2 plane. In the case of La2-x\[Sr,Ba\]xCuO4 and its Nd or Eu co-doped compound \[La,Nd,Eu\]2-x\[Sr,Ba\]xCuO4, hole doping is achieved by substituting ions of different ionicity, thereby creating an alloy with intrinsic inhomogeneities caused by chemical substitution burgy . Over the past several years, studies using 63Cu NQR and NMR wipeout hunt ; singer ; curro ; julien ; hunt2  have characterised the glassy nature of the slowing down of the stripe inhomogeneity tranquada  in these materials, where the Coulomb potential from the distorted lattice slows down spin and charge density waves. In the case of La2CuO4+δ, hole doping is achieved by super-oxygenation, and it has been shown that the high mobility of the excess oxygen atoms results in electronic phase separation sokol ; kivelson ; castro  between the superconducting and antiferromagnetic phase, as evidenced by 139La NMR hammel ; stat .

No clear picture has emerged which discerns and relates the effects of genuine electronic phase separation, stripe modulation, and the random substitution of donor ions. On the other hand, recent STM (scanning tunnelling microscopy) studies on the surface state of Bi2Sr2CaCu2O8+δ cleaved at low temperature pan  reveal spatial variations of the electronic state on a short length scale $``$nm in the surface plane. Whether such nm modulations are universally observable in the bulk and other high $`T_c`$ cuprates remains to be seen, but the STM results have enhanced the interest and potential impact of the spatial inhomogeniety of the electronic properties in cuprates. One emerging counter example is YBa2Cu3Oy where 89Y NMR alloulpre  measurements indicate that the spatial inhomogeneity in YBa2Cu3Oy is less than in Bi2Sr2CaCu2O8+δ or La2-xSrxCuO4.

Unlike scattering techniques which probe coherent phenomena over length scales larger than tens of nm’s, NMR and NQR are strictly $`local`$ probes which make them an ideal tool for revealing the short length scale physics of the CuO2 plane. Among earlier reports of the inhomogeneous state at short length scales in La2-xSrxCuO4 are: the splitting of the 63Cu NQR and NMR lines due to inequivalent Cu sites known as the ‘A” and “B”-sites resulting from different EFG (electric field gradient) tensors yoshimura ; stat  due to the presence of nearby +2Sr ions, the drastic broadening of the zero field 63Cu line at temperatures below $``$ 4 K tou ; hunt , the 63Cu NMR line broadening from short length scale modulations in orbital shifts haase , and also the 63Cu NQR and NMR wipeout hunt ; curro ; julien  as a result of the glassy slowing down of the stripe inhomogeneity in the temperature region $``$ 100 K.

These evidences clearly support the existence of some sort of short length scale inhomogeneity in the CuO2 planes of La2-xSrxCuO4 which raises questions regarding recent theoretical debates of a “universal electronic phase diagram”, including La2-xSrxCuO4 hwang , which are based on the assumption that hole doping is homogeneous. More recently, we reported evidence of an inhomogeneous electronic state in La2-xSrxCuO4 using 63Cu NQR in 63Cu isotope enriched samples in the range $`0.04x0.16`$ singerprl . As discussed in Ref. singerprl , the two essential ingredients for the spatial variation in local hole concentration $`{}_{}{}^{63}x_{local}^{}(0)`$ are (a) the fact that hole doping is achieved by random substitution of donor ions in the lattice with different ionicity, and (b) the presence of a short electronic length scale $`{}_{}{}^{63}R_{patch}^{}3.04.0`$ nm for the spatial variation. To the best of our knowledge, our results reported in this paper and in Ref. singerprl  are the first of its kind to detect the temperature dependence of the inhomogeneous electronic state in La2-xSrxCuO4, or any other high $`T_c`$ materials with quenched disorder.

Most recently, 17O NMR in high-quality La2-xSrxCuO4 single crystals in the range $`0.035x0.15`$ also reveal a substantial inhomogeneity in the local electronic state tobe . The spatial variation in local hole concentration $`{}_{}{}^{17}x_{local}^{}`$ determined using 17O NMR tobe  is found to reveal a consistent value with $`{}_{}{}^{63}x_{local}^{}`$ reported here and in Ref. singerprl . The spatial variation $`{}_{}{}^{17}x_{local}^{}`$ tobe  is determined through the observed spatial variation in the $`spin`$ susceptibility across the CuO2 plane, whereas $`{}_{}{}^{63}x_{local}^{}`$ is determined through the spatial variation in the EFG which is purely a charge effect. The fact that $`{}_{}{}^{63}x_{local}^{}^{\mathrm{\hspace{0.33em}17}}x_{local}`$ therefore implies that the inhomogeneity in the spin and charge channels are highly correlated.

As presented in Ref. tobe , we also compare the extent of the inhomogeneity in high-quality single crystals tobe  and our poly-crystalline samples (reported here and in Ref. singerprl ) by comparing the extent of the frequency dependence of $`{}_{}{}^{63}1/T_1`$ across the 63Cu NQR spectrum. We find consistent results between single crystal and poly-crystalline samples tobe , which immediately establishes that the inhomogeneity in the electronic state is an $`intrinsic`$ phenomenon in La2-xSrxCuO4. It also rules out certain claims that our poly-crystalline samples are somehow more inhomogeneous than those reported elsewhere. Such claims on the quality of our poly-crystalline samples originated as an excuse to negate the fact that hole doping is inhomogeneous in La2-xSrxCuO4 singerprl . In this paper we report the details of the 63Cu NQR results and analysis originally outlined in Ref. singerprl , and while the 17O NMR results are left to Ref. tobe , we mention the conclusions of the 17O NMR results wherever appropriate.

The dynamics of the electronic state are probed using the spin-lattice relaxation rate $`{}_{}{}^{63}1/T_1`$. $`{}_{}{}^{63}1/T_1`$ is a measure of the local spin fluctuations in the CuO2 plane at the 63Cu NQR frequency $`{}_{}{}^{63}\nu _{Q}^{}`$ (ranging between $`3240`$ MHz). Due to the experimental ease of the measurement, the temperature dependence of $`{}_{}{}^{63}1/T_1`$ is generally only measured at the peak of the spectrum. The most striking feature of our 63Cu NQR data, however, is that $`{}_{}{}^{63}1/T_1`$ shows qualitatively different temperature dependence depending on where along the NQR line it is measured singerprl . This implies, without using any kind of model, that certain regions of the CuO2 plane are more metallic while others are more insulating singerprl .

In order to measure the frequency dependence in $`{}_{}{}^{63}1/T_1`$ and obtain any quantitative information over a wide temperature range, it is $`essential`$ to separate the signal from 63Cu and 65Cu isotopes. Previous work by S. Fujiyama et al. fujiyama  measured $`{}_{}{}^{(63,65)}1/T_1`$ for naturally abundant Cu in which they successfully deduced that a substantial frequency dependence exists across the (63,65)Cu NQR spectrum, however, no statement beyond that could be inferred. Through 63Cu isotope enrichment, together with a systematic study as a function of $`x`$, we can use the extent of the frequency dependence in $`{}_{}{}^{63}1/T_1`$ across the 63Cu NQR line to determine the local hole concentration $`{}_{}{}^{63}x_{local}^{}`$ given by

$${}_{}{}^{63}x_{local}^{}=x\pm ^{63}\mathrm{\Delta }x_{local},$$
(1)
where $`{}_{}{}^{63}\mathrm{\Delta }x_{local}`$ is the characteristic amplitude or extent of the spatial variation $`{}_{}{}^{63}x_{local}^{}`$.

In Fig. 1 we show an example of the frequency dependence in $`{}_{}{}^{63}1/T_1`$ across the 63Cu NQR line at 300 K. We first remark on the overlap between samples with different nominal hole concentration $`x`$, shown in Fig. 1(b). Taking $`x=0.115`$ as an example, we see that the upper (lower) half intensity point of the $`x=0.115`$ A-line (or B-line) roughly coincides with the peaks of the $`x=0.16(0.07)`$, which implies that the characteristic local hole concentration $`{}_{}{}^{63}x_{local}^{}`$ is given by $`{}_{}{}^{63}x_{local}^{}0.16`$ in the more metallic regions, and by $`{}_{}{}^{63}x_{local}^{}0.07`$ in the more insulating regions. Using these limits and Eq. (1), we estimate the upper bound to the amplitude $`{}_{}{}^{63}\mathrm{\Delta }x_{local}0.045`$ for $`x=0.115`$ at 300 K. The reason that only an upper boundary to $`{}_{}{}^{63}x_{local}^{}`$ (and an upper bound to $`{}_{}{}^{63}\mathrm{\Delta }x_{local}`$) can be deduced is that we have assumed that the linewidths in Fig. 1(b) are dominated by $`{}_{}{}^{63}x_{local}^{}(x)`$ alone. In fact, as we calculate in section IV, there is a substantial intrinsic lattice linebroadening (defined as $`\mathrm{\Delta }\nu _{latt}`$) which is independent of the broadening arising from $`{}_{}{}^{63}x_{local}^{}`$. $`\mathrm{\Delta }\nu _{latt}`$ originates from the distribution in local EFG values due to the random positioning of +2Sr ions in the lattice, and we calculate its size by using a point charge lattice summation. In Fig. 2 we illustrate the extent of $`\mathrm{\Delta }\nu _{latt}`$ as the dashed curve, which is comparable to the experimentally observed linebroadening. In order to determine a better estimate of $`{}_{}{}^{63}x_{local}^{}`$ in a model independent way, we must therefore measure $`{}_{}{}^{63}1/T_1`$ at various positions across the line as shown in Fig. 1(a) and Fig. 2.

In order to determine $`{}_{}{}^{63}x_{local}^{}`$ from $`{}_{}{}^{63}1/T_1`$, we first define the values $`{}_{}{}^{63}1/T_{1,A}`$ and $`{}_{}{}^{63}1/T_{1,B}`$ taken at various positions across the NQR lineshape (shown in Fig. 2). $`{}_{}{}^{63}1/T_{1,A}^{(0)}`$ defines $`{}_{}{}^{63}1/T_1`$ measured at the CG (center of gravity) of the A-line, $`{}_{}{}^{63}1/T_{1,A}^{(+)}`$ defines $`{}_{}{}^{63}1/T_1`$ measured at the half intensity of the upper (+) frequency side of the A-line, $`{}_{}{}^{63}1/T_{1,A}^{()}`$ at the lower half intensity side, and $`{}_{}{}^{63}1/T_{1,A}^{(1/10)}`$ at the lower one-tenth intensity of the A-line. We also define the corresponding quantities for the B-line. We shall use these definitions throughout this paper to deduce $`{}_{}{}^{63}x_{local}^{}`$ as a function of temperature and $`x`$. Next, we assume that $`{}_{}{}^{63}1/T_{1,A}^{(0)}`$ and $`{}_{}{}^{63}1/T_{1,B}^{(0)}`$ taken at the CG of each spectrum represents the $`{}_{}{}^{63}1/T_1`$ value for the nominal hole concentration $`x`$, regardless of the underlying linebroadening. We can then compare $`{}_{}{}^{63}1/T_1`$ at the upper and lower half intensity points of the A-line (or B-line) to nominal values $`{}_{}{}^{63}1/T_{1,A}^{(0)}`$ from different $`x`$. As we shall show in detail in section III, $`{}_{}{}^{63}1/T_{1,A}^{()}`$ for $`x=0.115`$ at 300 K is the same as $`{}_{}{}^{63}1/T_{1,A}^{(0)}`$ for $`x=0.08`$, while $`{}_{}{}^{63}1/T_{1,A}^{(+)}`$ for $`x=0.115`$ is the same as $`{}_{}{}^{63}1/T_{1,A}^{(0)}`$ for $`x=0.15`$. Using Eq. (1) we therefore deduce that $`{}_{}{}^{63}\mathrm{\Delta }x_{local}=0.035`$ for $`x=0.115`$ at 300 K. Since we do not have the $`continuous`$ function $`{}_{}{}^{63}1/T_{1,A}^{(0)}`$ as a function of nominal $`x`$ (i.e. we have a limited number of samples), we interpolate $`{}_{}{}^{63}1/T_{1,A}^{(0)}`$ between the discrete set of $`x`$ we do have. The same argument applies to the B-line, and we shall demonstrate that $`{}_{}{}^{63}\mathrm{\Delta }x_{local}`$ deduced using the B-line is consistent with that from the A-line. In order to be consistent with Ref. singerprl , we define the characteristic amplitude $`{}_{}{}^{63}\mathrm{\Delta }x_{local}`$ with respect to the half intensity points, however, following the same proceedure using $`1/T_1^{(1/10)}`$ instead results in an overall 20 - 40 % increase in $`{}_{}{}^{63}\mathrm{\Delta }x_{local}`$ and does not effect our conclusions. Using the frequency dependence of $`{}_{}{}^{63}1/T_1`$, we therefore deduce that within each sample with nominal $`x`$ there exist certain regions in the CuO2 plane with higher and lower local hole concentrations ($`{}_{}{}^{63}x_{local}^{}x`$). Systematic studies as a function of temperature then allow us to follow the temperature dependence of $`{}_{}{}^{63}x_{local}^{}`$.

The next task in this paper is the model analysis of the NQR spectrum. For this task, we develop a static real space picture consisting of a patch-by-patch spatial variation for $`{}_{}{}^{63}x_{local}^{}`$ singerprl . Using a point charge lattice summation of the EFG together with a patch-by-patch spatial variation, we successfully fit the NQR spectrum (including the spectral overlap for different $`x`$ and the B-line) with one free parameter ($`{}_{}{}^{63}R_{patch}^{}`$). As detailed in section IV, the simulation of the EFG consists of two essential ingredients: (a) the random positioning of +2Sr donor ions in the lattice, and (b) the presence of a short length scale $`{}_{}{}^{63}R_{patch}^{}3`$ nm which describes the inhomogeneous electronic state. In our model, the CuO2 plane is sectioned into patches (i.e. circles) of equal radius $`{}_{}{}^{63}R_{patch}^{}(3`$ nm), and each patch $`i`$ has a different local hole concentration $`x_{local}^i(x)`$. The value of $`x_{local}^i`$ is determined by the random number of +2Sr donor ions in the vicinity of patch $`i`$. The holes are placed uniformily within patch $`i`$, and the resonance frequency of the 63Cu nuclei within patch $`i`$ is calculated in a self-consistent way using our pool of NQR data. By sampling a large number of patches ($`10^4`$) and calculating the resonance frequency for each patch, we build up a histogram spectrum of resonance frequencies which constitutes the simulated NQR lineshape. $`{}_{}{}^{63}R_{patch}^{}`$ is the only free parameter in the simulation, and we optimize $`{}_{}{}^{63}R_{patch}^{}`$ to fit the NQR spectrum, as shown in Fig. 2.

Using $`{}_{}{}^{63}R_{patch}^{}=\mathrm{}`$ is equivalent to uniform homogeneous doping where $`{}_{}{}^{63}x_{local}^{}=x`$ throughout the plane (i.e. $`{}_{}{}^{63}\mathrm{\Delta }x_{local}=0`$). This is also equivalent to excluding ingredient (b) listed above and calculating the broadening solely from the random positioning of +2Sr donor ions in the lattice (i.e. $`\mathrm{\Delta }\nu _{latt}`$). Clearly, however, in the case $`{}_{}{}^{63}R_{patch}^{}=\mathrm{}`$ the observed linewidth is underestimated by $`\mathrm{\Delta }\nu _{latt}`$ alone, likewise for the B-line. Using a very short length scale $`{}_{}{}^{63}R_{patch}^{}=1.6`$ nm overestimates the NQR linebroadening, whereas the optimum value $`{}_{}{}^{63}R_{patch}^{}=2.6`$ nm results in the best fit to the data, including the B-line. In Ref. tobe  we also show that $`{}_{}{}^{63}R_{patch}^{}`$ thus deduced is consistent with $`{}_{}{}^{17}R_{patch}^{}`$ determined in a similar way using the inhomogeneous linebroadening of the 17O NMR spectra.

Once $`{}_{}{}^{63}R_{patch}^{}`$ is optimized to fit the NQR spectra, we deduce the spatial variation $`{}_{}{}^{63}x_{patch}^{}`$ given by

$${}_{}{}^{63}x_{patch}^{}=x\pm ^{63}\mathrm{\Delta }x_{patch},$$
(2)
where $`{}_{}{}^{63}\mathrm{\Delta }x_{patch}`$ is the calculated amplitude for $`{}_{}{}^{63}x_{patch}^{}`$. $`{}_{}{}^{63}\mathrm{\Delta }x_{patch}`$ is found to follow the simple relation $`{}_{}{}^{63}\mathrm{\Delta }x_{patch}1/^{63}R_{patch}`$. The results of the EFG simulation ($`{}_{}{}^{63}x_{patch}^{}`$) are found to be consistent with $`{}_{}{}^{63}x_{local}^{}`$ deduced from the frequency dependence in $`{}_{}{}^{63}1/T_1`$. Since $`{}_{}{}^{63}x_{local}^{}`$ is determined in a model independent way, the fact that $`{}_{}{}^{63}x_{patch}^{}^{63}x_{local}`$ justifies the use of our patch-by-patch model. We note, however, that apart from the linebroadening mechanisms (a) and (b) discussed above, we ignore all other possible NQR linebroadening mechanisms in our model. Other possible mechanisms include distributions in local lattice distortions (section V). Our calculation in section IV ignores such possibilities which implies that $`{}_{}{}^{63}R_{patch}^{}`$ is a lower bound, that $`{}_{}{}^{63}\mathrm{\Delta }x_{patch}(1/^{63}R_{patch})`$ is correspondingly an upper bound, and therefore that $`{}_{}{}^{63}x_{patch}^{}`$ is an upper boundary. The fact that $`{}_{}{}^{63}x_{patch}^{}^{\mathrm{\hspace{0.33em}\hspace{0.33em}63}}x_{local}`$, however, indicates that our simulation of the NQR spectra and the patch-by-patch model is appropriate.

In section V we introduce the effects of the local orthorhombic distortions into the simulation. We find that the NQR data separates into two temperature regions, above and below $`T_Q(400)`$ K. For $`T>T_Q`$, our data indicates that the lattice is locally in the tetragonal phase and we use this region to estimate $`{}_{}{}^{63}R_{patch}^{}`$ and $`{}_{}{}^{63}x_{patch}^{}`$ as discussed above. For $`TT_Q`$, however, we find evidence for local orthorhombic distortions which, in the case of $`x0.04`$, start somewhat above the structural phase transition temperature $`T_{st}(`$ 515 K) according to neutron diffraction studies yamada . In the doping region $`x0.04`$, we show the systematic differences between the local and long range orthorhombic order.

The rest of this paper is sectioned as follows: in section II we describe the experimental characteristics including the sample characteristics and NQR background. In section III we present the experimental results which include the frequency dependence of $`{}_{}{}^{63}1/T_1`$ and the determination of $`{}_{}{}^{63}x_{local}^{}`$. We also present the temperature dependence of the NQR spectra. In section IV we describe the analysis of the NQR spectra using a patch-by-patch model for the spatial variation in $`{}_{}{}^{63}x_{local}^{}`$. We present the details of the EFG simulation incorporating randomness effects, and we optimize the length scale $`{}_{}{}^{63}R_{patch}^{}`$ (and $`{}_{}{}^{63}x_{patch}^{}`$) to fit the experimental data. In section V we incorporate the local orthorhombic distortions into the EFG simulation. In section VI we then present our conclusions.

## II EXPERIMENTAL Characteristics

### II.1 Samples

All of the single phase, poly-crystalline samples of La2-xSrxCuO4 used in this study were prepared using conventional solid state reactions takagisusc ; radaelli . We mix pre-dried La2O3 (99.995$`\%`$), SrCO3 (99.995$`\%`$), and 63CuO (99.995$`\%`$) with correct nominal compositions using an agate mortar and pestle until an intimate mixture is obtained. A pre-reaction is carried out for 20 h in a box furnace at $`850^{}`$C followed by repeated grindings and sinterings (also 20 h) at temperatures between $`950^{}`$C and $`1000^{}`$C. During these initial firings the materials are made into low density pellets with a hand press. The samples are then pelletized with the pneumatic press and high temperature annealed in flowing O2 gas at $`1100^{}`$C to $`1150^{}`$C for 24 h to 48 h before a slow and controlled cooling cycle that includes low temperature annealing at $`800^{}`$C (24 h) and $`500^{}`$C (24 h). The long annealing in O2 insures that the oxygen content is uniform and stoichiometric and high annealing temperatures insure fast reaction kinetics. By using a large number of grindings (typically 5-8), we achieve high quality poly-crystalline samples.

The most important test for the quality of our samples was discussed earlier (and presented in Ref. tobe ) where we confirmed that the frequency dependence of $`{}_{}{}^{63}1/T_1`$ across the 63Cu NQR spectrum in our poly-crystalline samples is the $`same`$ as in high-quality single crystals grown in Tokyo tobe . This confirms that the inhomogeneity we report in this paper and in Ref.’s singerprl ; tobe  are intrinsic to La2-xSrxCuO4. Apart from 63Cu NQR itself, we characterise our poly-crystalline samples by measuring the superconducting transition $`T_c`$ and the room temperature lattice parameters. We deduce $`T_c`$ by measuring the Meissner signal using a SQUID magnetometer in the field cooled mode with a constant applied magnetic field of 15 Oe , the results of which are shown in Fig. 3. We find that the onset temperature of the diamagnetic susceptibility, $`T_c`$, and the volume fraction agree with previous measurments in poly-crystalline samples prepared in a similar fashion takagisusc ; radaelli ; takagidelta . We note that the onset of the diamagnetic susceptibility is somewhat tailed in the case of $`x=0.20`$, consistent with takagisusc ; radaelli ; takagidelta . In Ref. takagidelta  it is even argued that bulk superconductivity disappears in the overdoped region, which coincides with the disappearance of LRO (long range order) into the orthorhombic structural phase. As we present in section III, we also find a qualitative change in $`{}_{}{}^{63}1/T_1`$ in the overdoped region. For these reasons we limit our determination of $`{}_{}{}^{63}x_{local}^{}`$ in section III to $`x0.16`$.

We deduce the lattice parameters with an X-ray diffractometer using the Cu $`K\alpha ^1`$ line. The lattice parameters presented in Fig. 4 (the values of which we use in the analysis in section IV) show close agreement with previous results radaelli . In this paper, we define $`a_o`$ and $`b_o`$ to be the in-plane lattice parameters of the orthorhombic cell and $`[a_o,b_o,c]`$ as the vector along the orthorhombic directions. Likewise, we define $`a`$ and $`b(a`$) as the lattice parameters in the tetragonal cell and $`[a,b,c]`$ as the vector along the tetragonal direction. Our typical X-ray linewidths were comparable to the instrumental resolution with a HWHM of $`\delta _{Instr}0.025^o`$ determined using a high quality Si standard. Previous higher resolution X-ray diffraction experiments takagidelta  (where $`\delta _{Instr}0.01^o`$) were able to use the observed X-ray linewidths and overlaps between different samples $`x`$ to deduce a lower bound to the compositional distribution with a HWHM of $`\delta x_{Sr}0.01`$. Only a lower bound to $`\delta x_{Sr}`$ could be inferred because X-ray diffraction takes a spatial average over length scales larger than tens of nm’s takagidelta . They attributed $`\delta x_{Sr}0.01`$ to imperfect mixing inherent to the solid state reaction. Similar conclusions were reached in Ref. radaelli  using the overlap of X-ray diffraction and neutron Bragg peaks between different $`x`$, where $`\delta x_{Sr}0.015`$ was deduced.

We confirmed the same order lower bound $`\delta x_{Sr}0.013`$ in single crystal $`x=0.15`$ singerprl  using EMPA (electron micro-probe analysis) whose length scale is determined by the focus area of the beam $`1`$ $`\mu `$m. This means that even high-quality single crystals whose spatial variation of +2Sr content is as small as $`\delta x_{Sr}0.013`$ ($`^{63}\mathrm{\Delta }x_{local}`$) over $`1`$ $`\mu `$m, exhibit a large variation in the local electronic states. Furthermore, we will demonstrate in section III that $`{}_{}{}^{63}\mathrm{\Delta }x_{local}`$ is temperature dependent and increases with decreasing temperature. This means that the variation in +2Sr content from the solid state reaction alone cannot account for our findings. Instead, a thermodynamic process must be involved.

### II.2 NQR Background

The essential condition for NQR is that the ground state spin $`{}_{}{}^{n}I`$ of the nuclear isotope $`n`$ satisfy $`{}_{}{}^{n}I1`$ abragam ; slichter . In the case of Cu, there exist two stable isotopes $`n=63(65)`$ with natural abundances of 69 (31) % respsectively, both with spin $`{}_{}{}^{(63,65)}I=3/2`$. This is both an advantage and a disadvantage. It is generally a disadvantage in the sense that their quadrupole moments $`{}_{}{}^{n}Q`$ happen to be similar $`{}_{}{}^{63}Q/^{65}Q=1.081`$, and therefore their resonance spectrum tend to coincide. This is the case in La2-xSrxCuO4 for $`x>0.05`$ where the NQR spectrum from each isotope is broad enough that they merge in the temperature range of interest, and therefore 63Cu isotope enrichment is essential to separate multiple lines.

The observed NQR frequency $`{}_{}{}^{63}\nu _{Q}^{}`$ is given by the total EFG (electric field gradient) tensor $`V_{tot}^{(\beta ,\gamma )}`$ (where $`\beta `$ and $`\gamma `$ are any two orthogonal spatial directions) surrounding the resonating 63Cu nucleus, and its quadrupole coupling constant $`{}_{}{}^{63}Q`$. In terms of magnetic eigenstates $`m_I`$, NQR for $`I=3/2`$ corresponds to transitions between $`m_I=|\pm 1/2|`$ and $`m_I^{}=|\pm 3/2|`$ eigenstates abragam ; slichter . The EFG cannot split the $`\pm `$ degeneracy of the magnetic eigenstates since it is of charge origin. The most general expression for the NQR tensor $`{}_{}{}^{63}\nu _{Q}^{(\beta ,\gamma )}`$ is given by abragam ; slichter :

$${}_{}{}^{63}\nu _{Q}^{(\beta ,\gamma )}=\frac{e^{63}Q}{2h}V_{tot}^{(\beta ,\gamma )},$$
(3)
where $`V_{tot}^{(\beta ,\gamma )}`$ is the traceless symmetric EFG tensor. $`V_{tot}^{(\beta ,\gamma )}`$ can be rotated and diagonalised to point towards the principle axes $`\alpha `$ of the EFG. This leads to the following relation

$${}_{}{}^{63}\nu _{Q}^{\alpha }=\frac{e^{63}Q}{2h}V_{tot}^\alpha .$$
(4)
In the tetragonal phase, the principle axes $`\alpha `$ are along the crystal axes $`[a,b,c]`$, and satisfy the general traceless condition $`V^c+V^b+V^a=0`$. The direction of the largest component $`|V^c|`$ (where $`|V^c|>|V^b|>|V^a|`$) defines the main principle axis in the tetragonal phase. The asymmetry parameter $`\eta `$ ($`0\eta 1`$), defined as

$$\eta =\frac{V^aV^b}{V^c},$$
(5)
is typically small $`\eta 0.06`$ for the planar Cu site in the cuprates penn  and for both A and B-sites song . The observed NQR frequency $`{}_{}{}^{63}\nu _{Q}^{}`$ is then given by

$${}_{}{}^{63}\nu _{Q}^{}=\frac{e^{63}Q}{2h}V_{tot}^c\left(1+\eta ^2/3\right)^{1/2}.$$
(6)
In the present case, $`\eta `$ is small $`(\eta 0.06)`$ and to a very good approximation

$${}_{}{}^{63}\nu _{Q}^{}\frac{e^{63}Q}{2h}V_{tot}^c.$$
(7)
The value of $`{}_{}{}^{63}Q`$ used varies between groups, but is generally taken to be either $`{}_{}{}^{63}Q=0.211`$ barns stern  which is based on theoretical calculation, or $`{}_{}{}^{63}Q=0.16`$ barns bleaney  which is based on ESR measurements in cuprate salts. Since we are dealing with experimental results, we choose to use $`{}_{}{}^{63}Q=0.16`$ barns bleaney . As discussed in appendix A, however, our final results of the length scale $`{}_{}{}^{63}R_{patch}^{}`$ is insensitive to the absolute value of $`{}_{}{}^{63}Q`$.

In NQR, the direction along which the spin-lattice relaxation rate $`{}_{}{}^{63}1/T_1`$ is measured is given by the main principle axis of the EFG (i.e. the $`c`$-axis), and cannot be changed externally. The most general expression for $`{}_{}{}^{63}1/T_1`$ in NQR with relaxation by a single magnetic process is given by moriya

$${}_{}{}^{63}\frac{1}{T_1}=\frac{{}_{}{}^{63}\gamma _{}^{2}k_BT}{\mu _B^2}\underset{𝐪}{}2\left|{}_{}{}^{63}A(𝐪)^{}\right|^2\frac{\chi _{}^{\prime \prime }(𝐪,\omega _n)}{\omega _n},$$
(8)
where $`\omega _n/2\pi =^{63}\nu _Q(=3240`$ MHz) is the NQR frequency, $`{}_{}{}^{63}\gamma /2\pi =11.285`$ MHz/Tesla is the nuclear gyro-magnetic ratio, and $`𝐪`$ is the reciprocal lattice vector (in tetragonal notation).

$`{}_{}{}^{63}A(𝐪)^{}`$ is the hyperfine form factor perpendicular to the main principle axis. The two perpendicular directions ($``$) in the present case both lie in the CuO2 plane (i.e. in the $`ab`$ plane). The $`𝐪`$ dependence of the in-plane value $`{}_{}{}^{63}A(𝐪)^{}`$ is given by

$${}_{}{}^{63}A(𝐪)^{}=A_{}+2B\left(\mathrm{cos}(q_xa)+\mathrm{cos}(q_ya)\right),$$
(9)
where $`A_{}38`$ kOe/$`\mu _B`$ and $`B42`$ kOe/$`\mu _B`$ are found to be the same in La2-xSrxCuO4 mmplsco ; imai  and YBa2Cu3Oy takigawa1991  compounds. $`\chi _{}^{\prime \prime }(𝐪,\omega _n)`$ (in emu/\[mol f.u.\]) is the in-plane imaginary part of the dynamic spin susceptibilty of the electrons. According to Eq. (9), there is a large contribution to $`{}_{}{}^{63}1/T_1`$ from the anit-ferromagnetic wavevector $`𝐐_{\mathrm{𝐀𝐅}}=(\pi /a,\pi /a)`$ where $`{}_{}{}^{63}A(𝐐_{\mathrm{𝐀𝐅}})^{}=139`$ kOe$`/\mu _B`$ imai . Since $`\chi _{}^{\prime \prime }(𝐪,\omega _n)`$ is peaked in the vicinity of $`𝐪=𝐐_{\mathrm{𝐀𝐅}}`$, this implies that $`{}_{}{}^{63}1/T_1`$ is a local probe of the low-frequency (i.e. $`\mathrm{}\omega _nk_BT`$) in-plane anti-ferromagnetic spin fluctuations. $`{}_{}{}^{63}\gamma _{}^{2}`$ enters into Eq. (8) because the spin-lattice relaxation process is magnetic in origin, whereas the splitting of the energy levels $`\omega _n`$ is determined by the EFG which is of charge origin.

We measure $`{}_{}{}^{63}1/T_1`$ by applying the following pulse sequence

$$\pi t[\pi /2\tau \pi ]\tau echo,$$
(10)
where an r.f. pulse $`\pi `$ inverts the nuclear magnetization $`M(t)`$ at time $`t`$ prior to the spin-echo sequence in brackets. The time dependence of $`M(t)`$ is determined by taking the integral of the $`echo`$ in the time domain as a function of delay time $`t`$. $`{}_{}{}^{63}1/T_1`$ is then determined by fitting $`M(t)`$ to the standard recovery form narath

$$M(t)=M(\mathrm{})+\left(M(0)M(\mathrm{})\right)\left[\mathrm{exp}\left(\frac{3}{T_1}t\right)\right],$$
(11)
appropriate for NQR with $`{}_{}{}^{n}I=3/2`$.

The results of $`{}_{}{}^{63}1/T_1`$ were independent of the r.f. pulse width $`\pi `$ in Eq. (10). In the case of $`x0.16`$, $`M(t)`$ fit well to Eq. (11), implying that the spin-lattice relaxation rate is dominated by a single relaxation mechanism. In the case of $`x0.16`$ below $`100`$ K, however, detailed measurements of $`M(t)`$ revealed a multi-exponential recovery, implying that $`{}_{}{}^{63}1/T_1`$ is distributed at a fixed frequency. In appendix B we use the multi-exponential recovery of $`M(t)`$ to deduce the intrinsic lattice broadening $`\mathrm{\Delta }\nu _{latt}^{T1}`$ of the NQR line. In the case of $`x=0.16`$, we find that $`\mathrm{\Delta }\nu _{latt}^{T1}=0.62(\pm 0.07)`$ MHz, which is consistent with our independent results of the EFG simulation $`\mathrm{\Delta }\nu _{latt}=0.49`$ MHz presented in section IV. All $`{}_{}{}^{63}1/T_1`$ data presented in this paper (except for appendix B) correspond to results of force-fits (in the case $`x0.16`$) deduced using Eq. (11). The force-fit values from Eq. (11) correspond to the average value of the underlying distribution in spin-lattice relaxation rate.

As a result of the distribution in the spin-lattice relaxation rate, the value of $`{}_{}{}^{63}1/T_1`$ determined using Eq. (11) showed a slight dependence on the pulse separation time $`\tau `$ in Eq. (10). The most likely explanation is that using a long $`\tau `$ results in suppressed signal intensity from 63Cu nuclei with fast spin-spin relaxation rates ($`{}_{}{}^{63}1/T_2`$). Since 63Cu nuclei with fast $`{}_{}{}^{63}1/T_1`$ also have fast component in $`{}_{}{}^{63}1/T_2`$ red , one would expect a decrease in the average value of $`{}_{}{}^{63}1/T_1`$ with increasing $`\tau `$. In the case of $`x=0.16`$, increasing $`\tau `$ from 12 to 24 $`\mu `$s resulted in a uniform $`10`$ % decrease for $`{}_{}{}^{63}1/T_1`$ across the NQR line. In order to retain the maximum possible signal intensity and measure the most representative value of $`{}_{}{}^{63}1/T_1`$ at each frequency, we used the shortest possible value $`\tau =12`$ $`\mu `$s the experiment allows for all measurements in this paper. Using the same experimental conditions throughout also allowed for a consistent determination of $`{}_{}{}^{63}x_{local}^{}`$.

Finally, we note that there are two relevant time-scales in the NQR experiment. The first is the ”dynamic” time-scale associated with $`{}_{}{}^{63}1/T_1`$. Since $`{}_{}{}^{63}1/T_1`$ measures the spin fluctuations specifically at the NQR frequency $`{}_{}{}^{63}\nu _{Q}^{}`$ (Eq. (8)), the appropriate ”dynamic” time-scale is given by $`1/^{63}\nu _Q0.027`$ $`\mu `$s. The other relevant time-scale is the ”static” one associated with the NQR spectrum. The ”static” time-scale is given by the spin-echo refocusing time $`2\tau =24`$ $`\mu `$s (Eq. (10)). $`2\tau `$ defines the effective exposure time (or equivalently, an effective shutter speed of $`1/(2\tau )`$) for the NQR experiment. All dynamic effects with time-scales shorter than $`2\tau `$ are unobservable through the NQR spectrum.

## III Experimental Results

### III.1 $`{}_{}{}^{63}1/T_1`$ Results

Before presenting the frequency dependence of $`{}_{}{}^{63}1/T_1`$, we first present the temperature dependence of $`{}_{}{}^{63}1/T_{1,A}^{(0)}`$ and $`{}_{}{}^{63}1/T_{1,B}^{(0)}`$ taken at the CG of each spectrum. In Fig. 5 and Fig. 6 we present the temperature dependence of $`{}_{}{}^{63}1/T_{1,A}^{(0)}`$ and $`{}_{}{}^{63}1/T_{1,B}^{(0)}`$, respectively. While the solid curves through the $`{}_{}{}^{63}1/T_{1,A}^{(0)}`$ data in Fig. 5 are guides for the eye, the solid curves in Fig. 6 represent $`\left(ϵ^{63}1/T_{1,A}^{(0)}\right)`$ data taken from Fig. 5. We plot Fig. 6 as such to directly compare the temperature dependence of $`\left(ϵ^{63}1/T_{1,A}^{(0)}\right)`$ with $`{}_{}{}^{63}1/T_{1,B}^{(0)}`$, where $`ϵ`$ is a uniform scaling factor taken to be $`ϵ=[0.87,0.87,0.90,0.84]`$ for $`x=[0.20,0.16,0.115,0.07]`$, respectively. As shown in Fig. 6, we find semi-quantitatively the same temperature dependence between the CG of the A and B-lines, with an overall 10 - 16 % smaller value of $`{}_{}{}^{63}1/T_{1,B}`$ compared to $`{}_{}{}^{63}1/T_{1,A}`$, i.e. $`{}_{}{}^{63}T_{1,A}^{}/^{63}T_{1,B}=0.90.84`$. Previous reports in $`non`$-isotope enriched samples itoh ; itoh2  also found that $`{}_{}{}^{63}T_{1,A}^{}/^{63}T_{1,B}<1`$.

At first glance of Fig. 6, our data for $`x=0.16`$ and $`x=0.20`$ below $`150`$ K may suggest that the ratio $`{}_{}{}^{63}T_{1,A}^{}/^{63}T_{1,B}`$ tends to decrease $`1015`$ % with decreasing temperature, similar to previous trends reported in itoh  for $`x=0.20`$. We point out, however, that caution must be used when deducing any significance from temperature variations in $`{}_{}{}^{63}T_{1,A}^{}/^{63}T_{1,B}`$. The reason for this stems from the fact that spin-lattice relaxation rate is distributed at a fixed frequency, and that finite $`\tau =12`$ $`\mu `$s conditions are used. As discussed in section II, the force-fit value $`{}_{}{}^{63}1/T_1`$ will tend to underestimate the CG of the underlying distribution, the extent of which will depend on $`\tau `$ and the one-to-one correspondence between the spin-spin relaxation rate and spin-lattice relaxation rate for all 63Cu nuclear sites. However, it is known that the spin-spin relaxation rate shows a complicated temperature dependence for both A and B-sites. In particular, the Gaussian-like component of the spin-spin relaxation rate ($`{}_{}{}^{63}1/T_{2G}`$ T2G ) tends to get smaller with decreasing temperature while the Lorenztian- like component ($`{}_{}{}^{63}1/T_{2L}`$) becomes dominant. Below the onset of 63Cu wipeout $`T_{NQR}`$, the spin-spin relaxation rate becomes totally Lorentzian-like hunt ; curro ; singer ; hunt2  at both A and B-sites. Slight differences in the temperature dependence of $`{}_{}{}^{63}1/T_{2G}`$ and $`{}_{}{}^{63}1/T_{2L}`$ between the A and B-sites complicate the one-to-one correspondence with spin-lattice relaxation rate, and result in different underestimations of $`{}_{}{}^{63}1/T_1`$ between the A and B-sites with decreasing temperature. Such factors must be considered when interpreting the temperature dependence in $`{}_{}{}^{63}T_{1,A}^{}/^{63}T_{1,B}`$ for $`150`$ K and $`x0.16`$.

In Fig. 7 we present the temperature dependence of $`{}_{}{}^{63}1/T_{1,A}`$ at various frequencies across the line, labeled according to Fig. 2. In the background of Fig. 7 we also show $`{}_{}{}^{63}1/T_{1,A}^{(0)}`$ taken from Fig. 5 which represent curves of constant $`x`$. The most surprising discovery of the present work is that $`{}_{}{}^{63}1/T_1^{(+)}`$, $`{}_{}{}^{63}1/T_1^{()}`$ and $`{}_{}{}^{63}1/T_1^{(1/10)}`$ all show qualitatively different temperature dependence. For example, $`{}_{}{}^{63}1/T_1^{(+)}`$ for $`x=0.07`$ exhibits semi-quantitatively the same behaviour as $`{}_{}{}^{63}1/T_{1,A}^{(0)}`$ for $`x=0.115`$, while $`{}_{}{}^{63}1/T_1^{()}`$ for $`x=0.07`$ exhibits semi-quantitatively the same behaviour as $`{}_{}{}^{63}1/T_{1,A}^{(0)}`$ for $`x=0.04`$. This is consistent with the fact that the upper (lower) frequency side of the NQR spectrum for $`x=0.07`$ roughly coincides with the peak NQR frequency of $`x=0.115`$ ($`x=0.04`$) respectively, as shown in Fig. 1. We also note that with decreasing temperature, $`{}_{}{}^{63}1/T_1^{(+)}`$, $`{}_{}{}^{63}1/T_1^{()}`$ and $`{}_{}{}^{63}1/T_1^{(1/10)}`$ do not exactly follow $`{}_{}{}^{63}1/T_{1,A}^{(0)}`$ for a given $`x`$, which indicates that $`{}_{}{}^{63}\mathrm{\Delta }x_{local}`$ is gradually growing with decreasing temperature.

In Fig. 8 we illustrate the details of the process used to extract $`{}_{}{}^{63}x_{local}^{}`$ from the data in Fig. 7. We first take the data from Fig. 7 and plot $`{}_{}{}^{63}1/T_{1,A}^{(+)}`$, $`{}_{}{}^{63}1/T_{1,A}^{()}`$ and $`{}_{}{}^{63}1/T_{1,A}^{(0)}`$ for each $`x`$ at a fixed temperature of 300 K. For clarity we connect the $`{}_{}{}^{63}1/T_{1,A}`$ data for fixed $`x`$ by the dashed black lines. We then create a smooth interpolation of $`{}_{}{}^{63}1/T_{1,A}^{(0)}`$ for all $`x`$ shown by the solid grey curve. We show the procedure used to determine $`{}_{}{}^{63}x_{local}^{}`$ in the case of $`x=0.07`$ shown as the solid black vertical and horizontal lines. In the case of $`{}_{}{}^{63}1/T_1^{(+)}`$, the horizontal black line determines what value of $`{}_{}{}^{63}1/T_{1,A}^{(0)}`$ (solid grey curve) corresponds to $`{}_{}{}^{63}1/T_1^{(+)}`$ , then the vertical line yields $`{}_{}{}^{63}x_{local}^{(+)}=0.10`$ $`(\pm 0.01)`$ for the more metallic regions. A similar proceedure for $`{}_{}{}^{63}1/T_1^{()}`$ yields $`{}_{}{}^{63}x_{local}^{()}=0.043`$ $`(\pm 0.008)`$ for the more insulating regions. The upper and lower frequency side therefore yield a consistent deviation $`{}_{}{}^{63}\mathrm{\Delta }x_{local}=0.028(\pm 0.009)`$ for $`x=0.07`$ at 300 K. We then determine $`{}_{}{}^{63}x_{local}^{(+)}`$ and $`{}_{}{}^{63}x_{local}^{()}`$ for all $`x`$ depending on available $`{}_{}{}^{63}1/T_{1,A}^{(0)}`$ data and desired accuracy. We cannot determine $`{}_{}{}^{63}x_{local}^{()}`$ for $`x=0.04`$ since we do not have $`{}_{}{}^{63}1/T_{1,A}^{(0)}`$ data below $`x=0.04`$ at 300 K. The reason for this is that $`{}_{}{}^{63}1/T_1`$ cannot be measured accurately for values $`7`$ (ms)-1. Furthermore, 300 K corresponds to the wipeout temperature $`T_{NQR}`$ for $`x0.035`$ hunt , and we only determine $`{}_{}{}^{63}x_{local}^{()}`$ from our data for $`TT_{NQR}`$ where full Cu signal intensity is observable. As for the limitations of $`{}_{}{}^{63}x_{local}^{(+)}`$, we note that at 300 K, $`{}_{}{}^{63}1/T_{1,A}^{(0)}`$ tends to merge to the same value of $`3`$ (ms)-1 for $`x0.09`$ imai . This implies that one needs to measure $`{}_{}{}^{63}1/T_1^{(+)}`$ beyond the experimental uncertainties in order to get a reliable estimate of $`{}_{}{}^{63}x_{local}^{(+)}`$.

The whole proceedure is repeated at different temperatures, the final results of which are summarized in Fig. 9. At different temperatures we find similar experimental limitations in determining $`{}_{}{}^{63}x_{local}^{(+)}`$ and $`{}_{}{}^{63}x_{local}^{()}`$, however, in cases such as $`x=0.07`$ where both can be determined, we consistently find that $`{}_{}{}^{63}\mathrm{\Delta }x_{local}^{(+)}^{\mathrm{\hspace{0.33em}63}}\mathrm{\Delta }x_{local}^{()}`$. In the overdoped region $`x0.20`$ we find a comparable frequency dependence of $`{}_{}{}^{63}1/T_1`$ across the NQR line to that of $`x=0.16`$. However, as shown in Fig. 8, we also find that $`{}_{}{}^{63}1/T_{1,A}^{(0)}`$ starts to $`increase`$ with increasing $`x`$ for $`x0.20`$ at 300 K, and the same is true at all temperatures. As discussed in section II, we limit our determination of $`{}_{}{}^{63}x_{local}^{}`$ to the region $`x0.16`$ and for clarity we separate the data in Fig. 8 for $`x>0.16`$ by using a dashed grey line.

In Fig. 10 we show the temperature dependence of $`{}_{}{}^{63}1/T_1`$ for $`x=0.16`$ across the superconducting boundary $`T_c`$ = 38 K (Fig. 3) at various positions across the NQR spectrum. Despite the fact that the spatial variation $`{}_{}{}^{63}x_{local}^{}`$ varies as much as $`0.10^{63}x_{local}0.22`$ at $`100`$ K, all values of $`{}_{}{}^{63}1/T_1`$ show a comparable fractional decrease below $`T_c`$ where the superconducting gap opens. In the case of $`{}_{}{}^{63}1/T_1^{()}`$ where $`{}_{}{}^{63}x_{local}^{()}=0.10`$, $`T_c`$ should be around 25 K, however, there is already a large drop in $`{}_{}{}^{63}1/T_1^{()}`$ by 30 K. In the case of $`{}_{}{}^{63}1/T_1^{(1/10)}`$ where $`{}_{}{}^{63}x_{local}^{(1/10)}=0.08`$, $`T_c`$ should be as low as 20 K, however $`{}_{}{}^{63}1/T_1^{(1/10)}`$ shows a drop below 30 K. This shows that the superconducting transition is a genuine bulk phenomenon which effects all patches, regardless of the local hole concentration $`{}_{}{}^{63}x_{local}^{}`$.

### III.2 $`{}_{}{}^{63}\nu _{Q}^{}`$ Results

In Fig. 11(a) and (b) we show the temperature dependence of the resonance frequency at the CG of the NQR spectrum defined as $`{}_{}{}^{63}\nu _{Q}^{k}`$ where $`k=(A,B)`$ for the A and B-lines, respectively. We define the first and second moments of the NQR spectrum as such

$`{}_{}{}^{63}\nu _{Q}^{k}`$ $`=`$ $`{\displaystyle \frac{_j^{63}\nu _{j,Q}^k}{_j}}`$
$`{}_{}{}^{63}\mathrm{\Delta }\nu _Q^k`$ $`=`$ $`\sqrt{\alpha _0}\sqrt{{\displaystyle \frac{_j\left({}_{}{}^{63}\nu _{j,Q}^{k}{}_{}{}^{63}\nu _{Q}^{k}\right)^2}{_j}}},`$ (12)
where $`{}_{}{}^{63}\nu _{j,Q}^{k}`$ corresponds to the $`j`$th data point of the observed NQR spectrum and $`\sqrt{\alpha _0}=\sqrt{2Ln2}=1.177`$. For a Gaussian spectrum (found to be the case for $`x>0.07`$) the use of the prefactor $`\sqrt{\alpha _0}`$ reduces $`{}_{}{}^{63}\mathrm{\Delta }\nu _Q^k`$ to the HWHM (half width at half maximum) exactly.

The $`{}_{}{}^{63}\nu _{Q}^{k}`$ data indicates two temperature regimes. Above $`T_Q`$ (black fit), we find a linear decrease of $`{}_{}{}^{63}\nu _{Q}^{A}`$ with decreasing temperature and a roughly constant value of $`{}_{}{}^{63}\nu _{Q}^{B}`$. In section IV we attribute the temperature dependence above $`T_Q`$ to thermal contraction of the lattice constants within the tetragonal phase, and in appendix A we also use this temperature region to determine the anti-shielding factors needed for the EFG simulation. Below $`T_Q`$ (grey fit), $`{}_{}{}^{63}\nu _{Q}^{k}`$ increases with decreasing temperature for both A and B-lines. The increase in this region is attributed to the increasing local orthorhombic distortions imai  which in the case of $`x0.04`$ set in above the LRO structural phase transition temperature $`T_{st}(515`$ K) according to neutron diffraction yamada . In the case of $`x=0.20`$, we find evidence for local orthorhombic distortion starting as high as $`T_Q`$ 350 K where diffraction results show no sign of LRO. In the section V we fit the data below $`T_Q`$ and deduce the degree of tilting of the CuO6 octahedra away from the $`c`$-axis (defined as $`\theta _{local}`$), and we find that there is a sharp onset temperature for $`\theta _{local}`$ only the case of $`x=0.0`$. For $`x>0.0`$ we see a somewhat rounded transition into the local orthorhombic phase. We note that these local precursive effects to the LRO are consistent with pair distribution function analysis of neutron powder diffraction data bozin  and XAFS (X-ray Absorption Fine Structure) analysis haskel  which find evidence for local lattice distortions in the temperature region above $`T_{st}`$, and also into the overdoped regime $`x0.20`$.

The temperature dependence of the observed HWHM $`{}_{}{}^{63}\mathrm{\Delta }\nu _Q^k`$ of the NQR spectra, defined using Eq. (12), is presented in Fig. 12. We also reproduce the temperature regions above and below $`T_Q`$, taken from Fig. 11, as the black and grey curves respectively. In section IV we use the $`{}_{}{}^{63}\mathrm{\Delta }\nu _Q^k`$ data for $`T>T_Q`$ to deduce a lower bound for the patch radius $`{}_{}{}^{63}R_{patch}^{}`$, and the corresponding an upper bound $`{}_{}{}^{63}\mathrm{\Delta }x_{patch}(1/^{63}R_{patch})`$.

## IV analysis

### IV.1 EFG Background

The goal of this section is to determine a static real-space model of the spatial variation in local hole concentration $`{}_{}{}^{63}x_{local}^{}`$. We shall use a point charge lattice summation incorporating randomness effects to calculate the inhomogeneous $`distribution`$ in the EFG which constitutes the observed 63Cu NQR spectrum. To the best of our knowledge, such calculations incorporating randomness effects into the EFG calculation have not been reported. We shall deduce all parameters necessary for this task in a self-consistent way by using our pool of NQR data, and we shall successfully account for the entire 63Cu NQR spectrum, including the B-line, with one adjustable parameter ($`{}_{}{}^{63}R_{patch}^{}`$) which defines the length scale of the spatial variation $`{}_{}{}^{63}x_{local}^{}`$. Once $`{}_{}{}^{63}R_{patch}^{}`$ is optimized to fit the 63Cu NQR spectrum, we deduce an upper bound to the amplitude of the spatial variation across the plane $`{}_{}{}^{63}\mathrm{\Delta }x_{patch}(1/^{63}R_{patch}`$), which we show is consistent with $`{}_{}{}^{63}\mathrm{\Delta }x_{local}`$ determined in a model independent way from section III.

The distribution in the EFG is dominated by three mechanisms (a) the random substitution of donor ions, (b) the variation in local hole concentration $`{}_{}{}^{63}x_{local}^{}`$ over short length scales $`{}_{}{}^{63}R_{patch}^{}3.0`$ nm, and (c) the distribution in local lattice distortions. Without a priori knowing the lattice distortions in (c), we attribute the whole NQR linebroadening to mechanisms (a) and (b). This implies that the deduced $`{}_{}{}^{63}R_{patch}^{}`$ we present is a lower bound and correspondingly $`{}_{}{}^{63}\mathrm{\Delta }x_{patch}(1/^{63}R_{patch})`$ is an upper bound. However we shall demonstrate that $`{}_{}{}^{63}\mathrm{\Delta }x_{patch}^{63}\mathrm{\Delta }x_{local}`$, which implies that using mechanisms (a) and (b) alone is justified in the temperature range $`T>T_Q`$. Below $`T_Q`$ where local orthorhombic distortions set in (Fig. 11), we terminate our analysis of $`{}_{}{}^{63}\mathrm{\Delta }x_{patch}`$. The reason for this is that we also see an increase in the NQR linewidth (Fig. 12) below $`T_Q`$ in the case of $`x0.07`$, which suggests that ignoring mechanism (c) above is no longer justified for $`T<T_Q`$. In section V we make predictions about the magnitude of mechanism (c) by attributing the extra linebraodening for $`T<T_Q`$ to distributions in the lattice distortions. We note, however, that with the exception of section V, we ignore mechanism (c) in our analysis.

There are two approaches for computing the EFG. The first is an $`ab`$ $`initio`$ approach martin ; plib ; husser  which involves quantum chemistry calculations, and the second approach makes use of the experimental data penn ; takigawa1989 ; shimizu ; walstedtprb ; singerprl  to deduce all the necessary parameters in an empirical way. The $`ab`$ $`initio`$ approach has been extensively used to calculate the EFG at the 63Cu site in La2CuO4 martin ; plib ; husser , and gives a consistent value for the observed resonance frequency $`{}_{}{}^{63}\nu _{Q}^{}33`$ MHz imai  to within uncertainties in the quadrupole moment $`{}_{}{}^{63}Q`$. The basic idea behind the ab initio approach is that since the EFG decreases rapidly $`1/r^3`$ away from the origin of the calculation, the most significant contributions should be from local EFG components. One can then justifiably separate a cluster of ions in the immediate vicinity of the Cu nucleus whose contributions are calculated using a full spin-polarized DF (density functional) or HF (Hatree-Fock) calculation, while the rest of the ions in the crystal are treated as point charges. The smallest realistic cluster is typically CuO6/Cu4La10, where the central (CuO$`{}_{6}{}^{})^{10}`$ ionic cluster consists of 23 molecular orbitals for each spin projection made up of linear combinations of 5 $`3d`$ and 18 $`2p`$ atomic orbitals. Such calculations are found to depend on the cluster size used, where the larger clusters are more reliable yet involve increasingly complex calculations. They also depend on whether one uses DF or HF, however both successfully predict that the molecular orbital with the highest energy is the anti-bonding hybridization between the Cu $`3d_{x^2y^2}`$ orbital and the O $`2p_x`$ and $`2p_y`$ orbitals. The DF method, however, predicts husser  more covalent-like bonding where the localized atomic spin density on the Cu is $`\rho _{Cu}=0.67`$, while HF predicts more ionic-like bonding with $`\rho _{Cu}=0.90`$.

Before discussing our new approach for calculating the inhomogeneous distribution of the EFG, let us clarify exactly what inhomogeneous broadening implies in the context of 63Cu NQR. The linebroadening of the 63Cu NQR spectrum for $`x>`$ 0.02 is dominated by the $`inhomogeneous`$ distribution in the EFG. The inhomogeneous broadening may be pictured as such: say there is a Cu nucleus $`j`$ lying in a particular EFG with value $`V_j`$ which resonates at frequency $`\nu _j`$, while a distant nucleus $`k`$ sits in a distinct EFG environment $`V_k`$ and resonates at $`\nu _k`$. As is predominantly the case across the sample, the separation between the two nuclei $`𝐫_{(j,k)}`$ is $`larger`$ than the range $`\xi _{AF}`$3 $`a`$ of the indirect nuclear spin-spin coupling pennslicht  where $`\xi _{AF}`$ is the correlation length of the anti-ferromagnetic fluctuations. In such cases, one may flip the $`j`$ nucleus using an r.f. pulse and observe its echo $`without`$ being effected by what the r.f. pulse does on the $`k`$ nucleus. One then measures the 63Cu NQR spectrum by resonating at $`\nu _j`$ and effectively counting the number of nuclei $`N_j`$ in the sample with EFG values $`V_j`$, then changing the resonance frequency to $`\nu _k`$ and counting the number of nuclei $`N_k`$, thereby building up the $`inhomogeneous`$ NQR spectrum.

As one approaches the undoped limit $`x<0.02`$, the random effects become less and less significant, and the linebroadening becomes predominantly $`homogeneous`$ in nature and dominated by the indirect nuclear spin-spin coupling slichter . As we show in appendix B, in the case of pure homogeneous broadening for $`x=0.0`$, there is no frequency dependence in $`{}_{}{}^{63}1/T_1`$ across the NQR spectrum, even in the orthorhombic phase $`T<T_{st}`$. We note that the two $`distinct`$ length scales $`\xi _{AF}`$ and $`{}_{}{}^{63}R_{patch}^{}`$ which determine the linebroadening of the homogeneous and inhomogeneous NQR spectra, respectively, have qualitatively different temperature dependences. As we shall show, $`{}_{}{}^{63}R_{patch}^{}`$ mildly decreases with decreasing temperature, while $`\xi _{AF}`$ increases with decreasing temperature yamada . Furthermore $`{}_{}{}^{63}R_{patch}^{}\xi _{AF}`$ within the $`x`$ and temperature range of interest.

### IV.2 EFG Simulation

In the experimental approach, the total EFG tensor $`V_{tot}^{(\beta ,\gamma ),k}`$ (see Eq. (7)) along the orthogonal spatial directions $`(\beta ,\gamma `$) for the Cu site $`k=(A,B)`$ is segregated into an onsite $`3d`$ contribution $`V_{3d}^{(\beta ,\gamma )}`$ and an offsite lattice contribution $`V_{latt}^{(\beta ,\gamma ),k}`$ as such

$$V_{tot}^{(\beta ,\gamma ),k}=\zeta _{3d}^kV_{3d}^{(\beta ,\gamma )}+\zeta _{latt}^kV_{latt}^{(\beta ,\gamma ),k},$$
(13)
where we introduce the isotropic anti-shielding factor slichter  $`\zeta _{3d}^k`$ originating from the local $`3d`$ electrons and the isotropic anti-shielding factor $`\zeta _{3d}^k`$ orginating from the lattice. We also leave the possibilty open that they depend on which Cu site $`k=(A,B)`$ is being calculated. In appendix A, we use a self-consistent approach to deduce that $`\zeta _{latt}^A=27.7`$ and $`\zeta _{3d}^A=0.93`$, while $`\zeta _{latt}^B=18.6`$ and $`\zeta _{3d}^B=0.82`$. We also show that the main principle value of the onsite EFG contribution $`V_{3d}^c`$ is negative while the lattice contribution $`V_{latt}^{c,k}`$ is positive, resulting in an overall negative value of $`V_{tot}^{c,k}`$.

The anti-shielding factors must be determined experimentally, while in the $`ab`$ $`initio`$ approach they are, in a sense, already accounted for. The anti-shielding factors are found $`not`$ to vary significantly between different classes of high $`T_c`$ cuprates shimizu . Likewise, the anti-shielding factors deduced in a similar way for the planar oxygen site takigawa1989 ; walstedtprb ; tobe  is also found to be consistent between different classes of cuprates, which supports the use of the experimental approach.

In accord with B. Bleaney and co-workers bleaney , the main principle value of the onsite contribution $`V_{3d}^c`$ arising from the $`3d_{x^2y^2}`$ hole is taken to be

$$V_{3d}^c=\frac{4}{7}e\left(14f_\sigma ^o\right)r_{3d}^3,$$
(14)
where we have allowed for covalency in the form of $`f_\sigma ^o`$ which represents the fractional spin denisty on the four neighbouring planar oxygen $`2p_\sigma `$ orbitals. In order to be consistent with our results of 17O hyperfine coupling analysis we use $`f_\sigma ^o=0.076`$ tobe  which favours more covalent-like bonding. Using ionic bonding with $`f_\sigma ^o=0`$ results in an overall $``$ 14 % increase in $`{}_{}{}^{63}R_{patch}^{}`$ and corresponding $``$ 12 % decrease in our theoretical estimate $`{}_{}{}^{63}\mathrm{\Delta }x_{patch}`$ over covalent-like bonding. We note that $`f_\sigma ^o=0.076`$ deduced from the oxygen hyperfine couplings walstedtprb ; tobe  is consistent with Cu hyperfine coupling analysis using NMR takigawa1991  where $`f_{Cu}=(14f_\sigma ^o)=0.7`$ was found. It is also generally accepted that the $`ab`$ $`initio`$ DF calculation which favors stronger covalency is a more realistic approach than the HF approach.

The lattice contribution $`V_{latt}^{(\beta ,\gamma )}`$ to the total EFG is calculated using a point charge lattice summation with fractional point charges. We use the standard expression for the summation

$$V_{latt}^{(\beta ,\gamma )}=\underset{j}{}q_j\left(3x_{\beta ,j}x_{\gamma ,j}r_j^2\delta _{\beta ,\gamma }\right)/r_j^5,$$
(15)
where the sum over $`j`$ refers to the sum over the surrounding point charges $`q_j`$ in the lattice a distance $`r_j`$ away. The charges $`q_j`$ are assigned as follows

$`q_{Cu}`$ $`=`$ $`+\left(24f_\sigma ^o\right)e`$
$`q_{O_{pl}}`$ $`=`$ $`\left(22f_\sigma ^ox/2\right)e`$
$`q_{O_{ap}}`$ $`=`$ $`2e`$
$`q_{La}`$ $`=`$ $`+3e`$
$`q_{Sr}`$ $`=`$ $`+2e,`$ (16)
where $`x`$ is the nominal hole concentration, Opl stands for the planar oxygen site, Oap stands for the apical oxygen site and $`e(>0)`$ is the electronic charge (in emu). Note that we make use of the full difference in valency between the +3La and +2Sr ions. This implies that depending on the the geometrical distribution of the +2Sr ions in the lattice, $`V_{latt}^{(\beta ,\gamma )}`$ takes on different values. In Fig. 13 we illustrate an example of the random positioning of the +2Sr ions surrounding the central Cu site from which the EFG is calculated.

The standard approach in all earlier works, except for singerprl , has been to set $`q_{La}=q_{Sr}=(3x/2)`$ which effectively bypasses any randomness effects and results in delta function spectra of the EFG. We start the EFG simulation by first taking $`{}_{}{}^{63}R_{patch}^{}=\mathrm{}`$ which uniformly places the donated holes $`x`$ from the +2Sr ions with full mobility onto the planar O sites. The added holes reduce $`|q_{O_{pl}}|`$ in Eq. (16), bearing in mind that there are two planar oxygens per unit cell. Placing the donated holes onto the planar oxygen sites as such is consistent with high energy spectroscopy studies fujimori .

The positions of the ions are allocated using an orthorhombic cell (shown in Fig. 14) which consists of $`4\times `$(La2-xSrxCuO4) formula units, each with the K2NiF4 structure ginsberg . Using an orthorhombic cell will later permit us to incorporate the orthorhombic distortions in section V. The absolute values of the lattice constants $`[a_o,b_o,c]`$ were calibrated at 295 K using our data from Fig. 4. The positions $`𝐫_j(x,T)`$ of the ions as a function of $`x`$ and $`T`$ were taken from a systematic X-ray powder diffraction study by P.G. Radaelli et al. radaelli . The EFG calculation is naturally sensitive to inputs $`𝐫_j(x,T)`$, therefore we use smooth interpolation of the values $`𝐫_j(x,T)`$ for different $`x`$ in order to avoid unnessecary scattering in the output $`V_{latt}^{(\beta ,\gamma )}`$. We also use the same thermal expansion coefficients radaelli ; braden2  $`\alpha _a=+1.4510^5`$ K-1 and $`\alpha _c=+1.4210^5`$ K-1 for all $`x`$, consistent with corresponding thermal coefficients found in La2-xBaxCuO4 suzuki .

The orthorhombic cell is then repeated in all space out to a radius of 50 Å from the origin $`[0,0,0]`$ of the EFG calculation. The procedure for the random placement of the +2Sr ions is carried out as follows: a different random number $`\lambda `$ with a flat probability distribution $`0<\lambda <1`$ is generated at each La site. If $`\lambda >x/2`$ at a particular La site, a La ion with charge $`+3e`$ is assigned to that site, while if $`\lambda <x/2`$, a Sr ion with charge $`+2e`$ is assigned. In equation form this gives the following prescription

$`\lambda >x/2`$ $`()`$ $`+3e`$
$`\lambda <x/2`$ $`()`$ $`+2e`$ (17)
at each La/Sr site. The factor 1/2 in the probability condition arises because there are 2 La/Sr sites per formula unit. Once this is carried out at each La site, we have achieved the particular random configuration $`\kappa `$ of +2Sr ions in the lattice, an illustration of which is given in Fig. 13. In our formulation, we exclude the possibility of any +2Sr - +2Sr clustering beyond probablilty theory, which is supported by the absence of any significant diffuse scattering in neutron diffraction experiments braden .

Once the random configuration $`\kappa `$ is determined, the summation in Eq. (15) is carried out and $`{}_{}{}^{\kappa }V_{latt}^{(\beta ,\gamma )}`$ is diagonalised and the principle value $`{}_{}{}^{\kappa }V_{latt}^{c}`$ is stored into a vector $`P_{latt}`$. The calculation is computed in the tetragonal phase where the main principle axis of $`{}_{}{}^{\kappa }V_{latt}^{(\beta ,\gamma )}`$ is found to lie primarily along the $`c`$-axis (i.e. $`{}_{}{}^{\kappa }V_{latt}^{c}`$ is the main principle value). We note, however, that each configuration $`\kappa `$ causes a slight rotation of the the main principle axis away from the $`c`$-axis. Strictly speaking, for each configuration $`\kappa `$ one should then go one step further and take the projection of the main principle value abragam  along the $`c`$-axis (i.e. along the average main principle axis). In the case of 63Cu NQR, however, the projection only makes $`2`$ % difference in $`{}_{}{}^{\kappa }V_{tot}^{c}`$, and can be neglected.

Next, the whole lattice is re-randomized and a new random configuartion $`\kappa ^{}`$ of +2Sr ions is determined, the lattice summation in Eq. (14) is computed and the new value of $`{}_{}{}^{\kappa ^{}}V_{latt}^{c}`$ is stored into $`P_{latt}`$. This proceedure is re-iterated $`10^4`$ times until the dimension of $`P_{latt}`$ is sufficiently large to create a histogram spectrum of $`P_{latt}`$ with $`50`$ bins across the lineshape, the results of which we show at 600 K in Fig. 15 as the dashed grey lines. Note that we always calculate the EFG for the Cu site at the origin using a different random configuration of +2Sr ions for each run. An equivalent approach is to stick to one random configuration $`\kappa `$ for the whole sample and calculate the EFG at each Cu site in the lattice thereby building up $`P_{latt}`$, however, this latter approach is computationally more intense. Since we are dealing with a purely random system, both methods are equivalent, therefore we use the former method.

We notice three immediate features about the results in Fig. 15. The first is the overall decrease in CG $`V_{latt}^c`$ with increasing $`x`$. This is a consequence of the change in lattice paramters and the decrease in $`|q_{O_{pl}}|`$ with increasing $`x`$. The second feature is the change in HWHM which is found to increase with increasing $`x`$ as a result of the increased amount of quenched disorder in the lattice. The third feature in our lattice summation is the presence of the secondary peak known as the B-line.

### IV.3 B-line

The 63Cu NQR spectrum is known to show a secondary peak known as the B-line yoshimura . Generally, the origin of the B-line has been attributed to structural effects from the dopant ions yoshimura ; song , but the limited range of $`x`$ investigated in earlier studies did not allow unambiguous identification of the origin of the B-line. In fact, there has been a persistent claim that the origin of the B-line is more exotic in nature hammel ; martin ; stat  and originates from an intrinsic response of the material to the presence of doped holes. Such a conjecture that the B-line is an electronic rather than a structural effect was based on the observed similarity between the A-B line splitting $`3`$ MHz for both La2-xSrxCuO4 and La2CuO4+δ systems hammel ; martin ; stat . It was thought that the same A-B line splitting between such distinct systems was clear evidence that the second site was a result of the presence of the doped holes themselves, independent of the means of doping. If this is the case, however, then it should equally well be observed in other 214 compounds such as La2-xBaxCuO4. However, a much larger A-B line splitting of $`6`$ MHz yoshimura ; hunt2  is observed in La2-xBaxCuO4 where the +2Ba ion is known to cause much larger ionic size effect disorder bozin  than the +2Sr ion, and which also stablizes an additional structural phase transition at low temperatures suzuki . We also recall that the 63Cu NQR spectrum in La2-x-yEuySrxCuO4 for 63Cu isotope enriched samples hunt2  shows a third structural peak (the C-line) whose fractional intensity is observed equal $`y`$. The C-line corresponds to a Cu nucleus directly above or below a Eu3+ ion.

In Fig. 16 we present the experimental value of the relative intensity ($`f_B^{exp}`$) of the B-line ($`I_B`$) to the total ($`I_A+I_B`$) defined as such

$$f_B^{exp}=\frac{I_B}{I_A+I_B},$$
(18)
where both $`I_A`$ and $`I_B`$ have been corrected for differences in spin-spin relaxation rates. Typically, $`{}_{}{}^{63}T_{2,A}^{}/^{63}T_{2,B}0.85`$ for both components of the spin-spin relaxation rate T2G , similar to the ratio of $`{}_{}{}^{63}1/T_1`$. Our data in Fig. 16 shows that $`f_B^{exp}x`$ over a large range of $`0.04x0.20`$, which is consistent with the EFG simulation where $`f_B=x`$. According to the EFG simulation, the Cu B-sites (illustrated in Fig. 14) correspond to Cu nuclei located at $`𝐫=[n_aa,n_ba,n_cc]`$ (where $`[n_a,n_b,n_c]`$ are $`\pm `$ integers ranging from zero to $`\mathrm{}`$), which have distinct EFG values due to the presence of +2Sr ions located at positions $`𝐫=[n_aa,n_ba,(n_c\pm 0.361)c]`$. In this interpretation, the relative intensity of such Cu sites ($`f_B`$) is then given by the concentration of +2Sr ions $`x`$, i.e. $`f_B=x`$. This site assignment for the B-site Cu nuclei can naturally account for the observed relative intensity $`f_B^{exp}x`$ (Fig. 16).

Other evidence that the B-line is primarily structural in origin is the similarity in the frequency dependence of $`{}_{}{}^{63}1/T_{1,B}`$ to $`{}_{}{}^{63}1/T_{1,A}`$ (Fig. 1). The similarity between in the frequency dependence of $`{}_{}{}^{63}1/T_{1,k}`$ indicates that the length scale of the spatial variation in $`{}_{}{}^{63}x_{local}^{}`$ is $`larger`$ than the average +2Sr - +2Sr distance $`l_{Sr}=a/\sqrt{x}`$, or equivalently the average B-site to B-site distance. Later we show that the deduced length scale $`{}_{}{}^{63}R_{patch}^{}3`$ nm is larger than $`l_{Sr}`$ in the region of interest $`x>0.02`$. Furthermore, as shown in Fig. 9, $`{}_{}{}^{63}x_{local}^{}`$ across the B-line is consistent with $`{}_{}{}^{63}x_{local}^{}`$ across the A-line.

We also note that in Fig. 6, $`{}_{}{}^{63}T_{1,A}^{}/^{63}T_{1,B}`$ = 0.9-0.84, i.e. that $`{}_{}{}^{63}1/T_{1,B}`$ is uniformly suppressed compared with $`{}_{}{}^{63}1/T_{1,A}`$. Our systmetic study up to $`x=0.20`$ showed that in the overdoped region $`x0.20`$, $`{}_{}{}^{63}1/T_1`$ increases with increasing $`x`$. This implies that if there is hole localisation in the vicinity of a Cu B-site as claimed in Ref. hammel ; martin ; stat , then $`{}_{}{}^{63}1/T_{1,B}`$ at the B-site should be $`larger`$ than $`{}_{}{}^{63}1/T_{1,A}`$ in the case of $`x=0.16`$. Fig. 6 clearly shows that this is not the case, and further supports our interpretation that the B-line is primarily structural in origin. The most likely explanation for $`{}_{}{}^{63}T_{1,A}^{}/^{63}T_{1,B}=0.90.84`$ (and $`{}_{}{}^{63}T_{2,A}^{}/^{63}T_{2,B}`$) is that the hyperfine coupling constant $`|^{63}A(𝐐_{\mathrm{𝐀𝐅}})^{}|`$ is 5-10 % lower at the B-site due to the increased lattice distortions in the vicinity of a +2Sr ion. It is known from XAFS haskelrapid  that the +2Sr ion causes large structural distortions in its vicinity, which could well account for such changes in the hyperfine couplings.

### IV.4 EFG $`^{63}\nu _Q`$

We now proceed with our simulation taking the B-site assignment shown in Fig. 14. We can separate the A and B-sites within the calculation by separating any randomly generated $`{}_{}{}^{\kappa }V_{latt}^{c}`$ value with a +2Sr ion in a position illustrated in Fig. 14, and placing it into a new vector $`P_{latt}^B`$ , while placing all other generated $`{}_{}{}^{\kappa }V_{latt}^{c}`$ values into another vector $`P_{latt}^A`$. This allows us to separate the A and B-lines as shown by the solid black lines in Fig. 15, and we attach the superscript $`k=(A,B)`$ to $`{}_{}{}^{\kappa }V_{latt}^{c,k}`$. The first and second moments of the lattice contribution to the EFG are then computed as such

$`V_{latt}^{c,k}`$ $`=`$ $`{\displaystyle \frac{_jP_{j,latt}^k}{_j}}`$
$`\mathrm{\Delta }V_{latt}^c`$ $`=`$ $`\sqrt{\alpha _0}\sqrt{{\displaystyle \frac{_j\left(P_{j,latt}^kV_{latt}^{c,k}\right)^2}{_j}}}.`$ (19)
In order to determine the total EFG defined in Eq. (13) , we now determine the anti-shielding factors for both the A and B-lines from the CG data as such

$${}_{}{}^{63}\nu _{Q}^{k}=\frac{e^{63}Q}{2h}\left(\zeta _{3d}^kV_{3d}^c+\zeta _{latt}^kV_{latt}^{c,k}\right),$$
(20)
where we have used Eq. (13) and Eq. (19), along with the experimental data $`{}_{}{}^{63}\nu _{Q}^{k}`$ (see appendix A for more details). Eq. (20) and the anti-shielding factors convert the CG of the EFG simulation to the CG of the experimentally observed data.

We are now in a position to convert each EFG value ($`{}_{}{}^{\kappa }V_{latt}^{c,k}`$) stored in $`P_{latt}^k`$ into a frequency as such

$${}_{}{}^{\kappa }\nu _{}^{k}=\frac{e^{63}Q}{2h}\left(\zeta _{3d}^kV_{3d}^c+\zeta _{latt}^k{}_{}{}^{\kappa }V_{latt}^{c,k}\right).$$
(21)
Using Eq. (20) and Eq. (21), we can finally predict the intrinsic lattice linebroadening ($`\mathrm{\Delta }\nu _{latt}^k`$) as such

$$\mathrm{\Delta }\nu _{latt}^k=\frac{e|^{63}Q|}{2h}\zeta _{latt}^k\mathrm{\Delta }V_{latt}^c.$$
(22)
First, note that in Eq. (22) and Eq. (19) we have $`not`$ attach a superscript $`k=(A,B)`$ to $`\mathrm{\Delta }V_{latt}^c`$. This is because the EFG simulation predicts the $`same`$ width $`\mathrm{\Delta }V_{latt}^c`$ for both the A and B-lines. Second, note that we have assumed that the total width from the calculation $`\mathrm{\Delta }\nu _{latt}^k`$ is proportional to the intrinsic lattice width $`\mathrm{\Delta }V_{latt}^c`$ with NO corresponding width from the onsite contribution, i.e. $`\mathrm{\Delta }V_{3d}^c=0`$. We justify this by the fact that the ratio of widths between the A and B-lines is experimentally found to be $`{}_{}{}^{63}\mathrm{\Delta }\nu _Q^B/^{63}\mathrm{\Delta }\nu _Q^A0.63`$ in the temperature region $`T>T_Q`$. According to appendix A, this ratio is closer to the ratio of the anti-shielding factors for the lattice contribution $`\zeta _{latt}^B/\zeta _{latt}^A=0.67`$ than for the onsite contribution $`\zeta _{3d}^B/\zeta _{3d}^A=0.89`$. Note that the anti-shielding factors are deduced without any information about the widths, therefore taking $`\mathrm{\Delta }V_{3d}^c=0`$ is justified.

Up to this point in the EFG simulation, we have taken care of the linebroadening due to the random placement of +2Sr ions in the lattice (mechanism (a)), however, we are still assuming that the holes are uniformily distributed across the plane, i.e. $`{}_{}{}^{63}R_{patch}^{}=\mathrm{}`$. Clearly this is not the whole picture since our calculation of $`\mathrm{\Delta }\nu _{latt}^k`$ underestimates the observed linewidth $`{}_{}{}^{63}\mathrm{\Delta }\nu _Q^k`$ by a factor $`2`$, as can be seen in Fig. 2.

### IV.5 $`{}_{}{}^{63}R_{patch}^{}`$ Simulation

The introduction of the spatial variation in local hole concentration $`{}_{}{}^{63}x_{local}^{}`$ into the EFG calculation is an extension of what we have described in arriving at the intrinsic width $`\mathrm{\Delta }\nu _{latt}^k`$. We start by randomly positioning the +2Sr donor ions into the +3La sites as described earlier. Next, we define a patch (i.e. a circle) of radius $`{}_{}{}^{63}R_{patch}^{}`$ around the Cu site at the origin, as shown in Fig. 13. We then determine the total number of Cu sites $`N^{R_{patch}}`$ within $`{}_{}{}^{63}R_{patch}^{}`$ by multipying the area ($`\pi ^{63}R_{patch}^2`$) the patch covers by the areal Cu density ($`1/a^2`$) as such

$$N^{R_{patch}}=\frac{\pi ^{63}R_{patch}^2}{a^2},$$
(23)
where we round $`N^{R_{patch}}`$ to the nearest integer. These $`N^{R_{patch}}`$ Cu sites within the radius $`{}_{}{}^{63}R_{patch}^{}`$ will uniformly share the donated holes from neighbouring +2Sr ions. Next, we determine how many donated holes $`N_{Sr}^\kappa `$ are to be shared within $`{}_{}{}^{63}R_{patch}^{}`$, where $`N_{Sr}^\kappa `$ itself is the random number to be calculated. As shown in Fig. 14, there are two distinct Sr/La sites (labeled $`d_1`$ and $`d_2`$), from which we can donate holes. This comes about because there are two distinct Sr/La sites within the unit cell. For donation of holes from $`d_1`$, we accept a donor hole into $`{}_{}{}^{63}R_{patch}^{}`$ provided there is a +2Sr ion which satisfies the position (in tetragonal notation)

$`𝐫^{d_1}`$ $`=`$ $`[n_aa,n_ba,\pm 0.361c]`$
$`\sqrt{n_a^2+n_b^2}`$ $`<`$ $`{}_{}{}^{63}R_{patch}^{}/a`$ (24)
relative to the origin. For donation of holes from $`d_2`$, we accept a donor hole into $`{}_{}{}^{63}R_{patch}^{}`$ provided there is a +2Sr ion which satisfies the position

$`𝐫^{d_2}`$ $`=`$ $`[(n_a+1/2)a,(n_b+1/2)a,\pm 0.139c]`$
$`\sqrt{n_a^2+n_b^2}`$ $`<`$ $`{}_{}{}^{63}R_{patch}^{}/a,`$ (25)
and in both cases, $`n_a`$ and $`n_b`$ are $`\pm `$ integers ranging from zero to $`\mathrm{}`$. We do not a priori know which site $`d_1`$ or $`d_2`$ are the donor sites, or if the donors come from a combination of $`d_1`$ and $`d_2`$. We therefore treat both or combinations of both as equal possibilities. Technically, the n.n. (nearest neighbour) $`d_2`$ is the closest hole donor to the central Cu nucleus when occupied by a +2Sr ion. The n.n. $`d_1`$ site on the other hand has more bonding with the central Cu via the apical oxygens. The n.n $`d_1`$ site has the largest effect on $`{}_{}{}^{63}\nu _{Q}^{}`$ when occupied by a +2Sr ion and in particular gives rise to the B-site.

Once we have counted all the $`N_{Sr}^\kappa `$ satisfying either Eq. (24) or Eq. (25) for a particular random configuration of +2Sr ions $`\kappa `$, we deduce the local hole concentration $`{}_{}{}^{63}x_{local}^{\kappa }`$ as such

$${}_{}{}^{63}x_{local}^{\kappa }=N_{Sr}^\kappa /N^{R_{patch}}.$$
(26)
We now carryout the EFG summation in Eq. (14) for the particular configuration of +2Sr ions $`\kappa `$, where the planar oxygen charge $`q_{O_{pl}}`$ in Eq. (16) is replaced by

$$q_{O_{pl}}^\kappa =\left(22f_\sigma ^o^{63}x_{local}^\kappa /2\right)e$$
(27)
for all planar oxygens within the patch radius $`{}_{}{}^{63}R_{patch}^{}`$. All other planar oxygen sites outside the radius $`{}_{}{}^{63}R_{patch}^{}`$, including those on other CuO2 planes, are assigned the uniform hole concentration $`x`$. Using this scheme effectively ignores 3-d correlations of $`{}_{}{}^{63}x_{local}^{\kappa }`$ from neighbouring CuO2 planes, however, such effects are small and of higher order.

As before, we then convert $`{}_{}{}^{\kappa }V_{latt}^{c,k}`$ into a frequency using Eq. (21) and we store $`{}_{}{}^{\kappa }\nu _{}^{k}`$ into a vector $`P_R^k`$. The position of the +2Sr ions are then re-randomized to obtain a new random configuration $`\kappa ^{}`$, then a new value of $`{}_{}{}^{\kappa ^{}}\nu _{}^{k}`$ is deduced and stored into $`P_R^k`$. This procedure is repeated $`10^4`$ times until a sufficient number of data points exist in $`P_R^k`$ so as to create a histogram spectrum of the lineshape over $``$50 bins. The CG of the spectra $`P_R^k`$ coincide with the experimental data $`{}_{}{}^{63}\nu _{Q}^{k}`$ (Eq. (20)), while a new HWHM ($`\mathrm{\Delta }\nu _{R_{patch}}^k`$) is determined as such

$$\mathrm{\Delta }\nu _{R_{patch}}^k=\sqrt{\alpha _0}\sqrt{\frac{_j\left(P_{j,R}^k{}_{}{}^{63}\nu _{Q}^{k}\right)^2}{_j}}.$$
(28)
We then repeat this whole procedure and calculate $`\mathrm{\Delta }\nu _{R_{patch}}^k`$ for various values of $`{}_{}{}^{63}R_{patch}^{}`$, the results of which are shown in Fig. 17 for the A-line. The results for the B-line are the same but with a uniform decrease $`\mathrm{\Delta }\nu _{R_{patch}}^B/\mathrm{\Delta }\nu _{R_{patch}}^A=\zeta _{latt}^B/\zeta _{latt}^A(=0.67`$) for all $`{}_{}{}^{63}R_{patch}^{}`$. Also shown in Fig. 17 is how to deduce the best fit value of $`{}_{}{}^{63}R_{patch}^{}`$ to the lineshape at 600 K. We match experimental values of $`\mathrm{\Delta }\nu _Q^k`$ (shown by the grey horizontal line) to the calculation shown by the solid black curves, and extrapolate $`{}_{}{}^{63}R_{patch}^{}`$ shown by the grey vertical arrows.

Note that according to the arrows, $`{}_{}{}^{63}R_{patch}^{d_2}`$ deduced using dopants from +2Sr ions at $`d_2`$ is 15 % smaller than $`{}_{}{}^{63}R_{patch}^{d_1}`$ deduced using dopants from +2Sr ions at $`d_1`$. Put in another way, when we do the calculation for fixed $`{}_{}{}^{63}R_{patch}^{}`$, we deduce a $``$ 15 % smaller value of $`\mathrm{\Delta }\nu _{R_{patch}^{d_2}}^k`$ compared to $`\mathrm{\Delta }\nu _{R_{patch}^{d_1}}^k`$. The reason for $`\mathrm{\Delta }\nu _{R_{patch}^{d_2}}^k<\mathrm{\Delta }\nu _{R_{patch}^{d_1}}^k`$ can be explained by a correlation between the local hole concentration $`{}_{}{}^{63}x_{local}^{\kappa }`$ on the planar oxygens and the +2Sr donor ions in n.n. positions to the origin. The result of adding one hole from a n.n. +2Sr ion amounts to two separate shifts to the resonance frequency. The first shift is from the added hole itself which decreases of magnitude of the surrounding planar oxygen charge $`|q_{O_{pl}}|`$ and results in a large $`increase`$ in the resonance frequency, regardless of the hole’s origin. We shall define this shift as $`\delta \nu ^{pl}(>0)`$. The second effect of adding one hole from a n.n. +2Sr donor ion is from the n.n. +2Sr ion itself which results in a shift $`\delta \nu ^{Sr}`$, however, the sign of shift $`\delta \nu ^{Sr}`$ depends on whether it is at $`d_2`$ or $`d_1`$.

Let us imagine scenario (1) where the added hole specifically comes from a +2Sr ion in the n.n.n. (next nearest neighbour) $`d_1`$ site. We do not use the n.n. $`d_1`$ site in this case since this is the B-line which, according to the simulation, shows the $`same`$ difference $`\mathrm{\Delta }\nu _{R_{patch}^{d_2}}^B/\mathrm{\Delta }\nu _{R_{patch}^{d_1}}^B(0.85`$) as the A-line. It can be shown that a +2Sr ion in the n.n.n. $`d_1`$ site causes a small positive shift $`\delta \nu ^{d_1}>0`$, therefore the total shift $`\delta \nu ^{(1)}`$ in scenario (1) is given by $`\left(\delta \nu ^{pl}+\delta \nu ^{d_1}\right)=\delta \nu ^{(1)}(>0)`$ where both components are positive and enhance the total effect. Next let us imagine scenario (2) where the +2Sr ion from which the added hole came from is specifically located at a n.n. $`d_2`$ position to the origin. It can be shown that placing a +2Sr ion at the n.n. $`d_2`$ position causes a small $`negative`$ shift $`\delta \nu ^{d_2}<0`$ to the resonance frequency. The total shift $`\delta \nu ^{(2)}`$ in scenario (2) is therefore given by $`\left(\delta \nu ^{pl}+\delta \nu ^{d_2}\right)=\delta \nu ^{(2)}(>0)`$ where the individual components counteract each other. Therefore, in scenario (1) the overall positive shift to the resonance frequency is enhanced while in scenario (2) it is reduced. After many random runs which include the independent random effects from distant +2Sr ions, scenario (1) will lead to a larger overall spread, therefore $`\mathrm{\Delta }\nu _{R_{patch}^{d_2}}^k<\mathrm{\Delta }\nu _{R_{patch}^{d_1}}^k`$.

### IV.6 $`{}_{}{}^{63}R_{patch}^{}`$ Results

In Fig. 18 we show the results for the best fit to the lineshape data at 600 K. Note that the observed asymmetry of the A and B-lines for $`x0.04`$ are well reproduced within the model without any additional parameters. The fits shown in Fig. 18 are deduced using $`d_1`$ dopants with $`{}_{}{}^{63}R_{patch}^{d_1}`$, however the same fits can be deduced using $`d_2`$ dopants with $`{}_{}{}^{63}R_{patch}^{d_2}`$. The best fit values of $`{}_{}{}^{63}R_{patch}^{d_1}`$ and $`{}_{}{}^{63}R_{patch}^{d_2}`$ are shown in Fig. 19(a). In Fig. 18 we also show the data and fit for $`x=0.02`$. As can be seen, however, the asymmetry in the experimentally observed line-profile for $`x=0.02`$ is somewhat larger that our patch model can predict, hence we show the fit as a dashed line. It is interesting to note that in the range $`x0.02`$, $`{}_{}{}^{63}R_{patch}^{d_1}`$ and $`{}_{}{}^{63}R_{patch}^{d_2}`$ tend towards the average donor separation distance $`l_{Sr}=a/\sqrt{x}`$, as shown in Fig. 19. If $`{}_{}{}^{63}R_{patch}^{}l_{Sr}`$, our model implies that each patch only covers one +2Sr ion on average. In such a senario we no longer expect our model to account for the data. In fact, the critical region $`x=0.02`$ which separates the spin-glass phase and the anti-ferromagnetic phase for $`x0.02`$ cho  is known to show complex behaviour rigamonti .

We also note in Fig. 18 that the intensity of the B-line is slightly underestimated by the calculation compared with the experimental data. This is simply due to the fact that the spin-spin relaxation rate is experimentally $``$ 15 % smaller at the B-site together with the fact that we are measuring the data points using $`\tau =12`$ $`\mu `$s. After correcting for differences in the spin-spin relaxation rate, the observed relative intensity of the B-line decreases slightly and is found to go as $`f_B^{exp}x`$ (Fig. 16), which is then consistent with the simulation. In order to fit the NQR spectrum in Fig. 18, we have attributed all of the extra line broadening to $`{}_{}{}^{63}R_{patch}^{}`$. In Fig. 19(b) we compare the experimentally observed width $`{}_{}{}^{63}\mathrm{\Delta }\nu _Q^A`$ with the calculated intrinsic width $`\mathrm{\Delta }\nu _{latt}^A`$.

In Fig. 20 we show the same plot as Fig. 19 except at 300 K, where the $`{}_{}{}^{63}R_{patch}^{}`$ are deduced in a similar way to 600 K shown in Fig. 17. At 300 K, we see an effect not present at 600 K, namely that the observed linewidth $`{}_{}{}^{63}\mathrm{\Delta }\nu _Q^A`$ for $`x0.115`$ in Fig. 20(b) no longer decreases with decreasing $`x`$ as it does at 600 K. Also shown in Fig. 20 is the onset of the LRO orthorhombic phase yamada  shown as the dashed lines through the data points. The result of the observed enhancement of the linewidths $`{}_{}{}^{63}\mathrm{\Delta }\nu _Q^A`$ for $`x0.115`$ in Fig. 20(b) results in a large drop in the calculated $`{}_{}{}^{63}R_{patch}^{}`$ shown in Fig. 20(a). Comparing the $`x`$ dependence in $`{}_{}{}^{63}R_{patch}^{}`$ within the tetragonal phase at 600 K and at 300 K for $`x0.115`$ seems to suggest that $`{}_{}{}^{63}R_{patch}^{}`$ is more or less $`x`$ independent $`provided`$ that the local lattice distortions do not dramatically enhance the observed NQR linebroadening.

The effects of the orthorhombic distortions can also be seen in Fig. 21 where we plot the temperature dependence of the average $`{}_{}{}^{63}R_{patch}^{}`$ defined as such

$${}_{}{}^{63}R_{patch}^{}=\left({}_{}{}^{63}R_{patch}^{d_1}+^{63}R_{patch}^{d_2}\right)/2.$$
(29)
We chose to plot the average value $`{}_{}{}^{63}R_{patch}^{}`$ alone for clarity. This corresponds to a simulation where half of the donors come from $`d_1`$ sites and half come from $`d_2`$ sites. In Fig. 21 the black and grey lines correspond to regions above and below $`T_Q`$ taken from Fig. 11, while the dashed black line shows the onset of LRO $`T_{st}`$ yamada . In the case of $`x=0.04`$ we clearly see a dip in $`{}_{}{}^{63}R_{patch}^{}`$ below $`T_QT_{st}`$ which is most likely due to the fact that we have neglected distributions in the local orthorhombic distortions. In section V we make the reasonable assumption that $`{}_{}{}^{63}R_{patch}^{}`$ for $`x=0.04`$ continues to decrease lineaarly (defined as $`{}_{}{}^{63}R_{patch}^{}`$) with decreasing temperature below $`T_Q`$, similar to $`x0.115`$ samples. Assuming that $`{}_{}{}^{63}R_{patch}^{}`$ is the underlying length scale will allow us to estimate the distribution in local orthorhombic distortions.

We also remark that similar patch-by-patch models of disordered doping in the CuO2 plane have been independently used to account for the smeared density of states found in photoemission peaks as a result of a distribution in the superconducting gap bala . In such models the only free parameter used to account for the smoothness of the density of states is the length scale of the model which is determined by the in plane superconducting coherence length $`\xi _{ab}2.03.0`$ nm. It is interesting to note that in the case of $`x=0.16`$, the lower bound $`{}_{}{}^{63}R_{patch}^{}2.5\pm 0.5`$ nm in the region of $`T_c`$ is comparable to $`\xi _{ab}3.3`$ nm in thin films suzukicoh .

### IV.7 $`{}_{}{}^{63}\mathrm{\Delta }x_{patch}`$ Results

In the temperature region above $`T_Q`$, we can also deduce an upper boundary $`{}_{}{}^{63}x_{patch}^{}`$ for the spatial variation $`{}_{}{}^{63}x_{local}^{}`$ without any further computation. We proceed as described above, but for each random lattice configuration $`\kappa `$, we also store each value $`{}_{}{}^{63}x_{local}^{\kappa }`$ (used to compute $`{}_{}{}^{\kappa }\nu _{Q}^{k}`$) into a new vector $`X_{local}^k`$. We then build up a vector $`X_{local}^k`$ of randomly distributed local hole concentrations for each lattice configuration $`\kappa `$. The mean of the set $`X_{local}^k`$ is just the nominal value $`X_{local}^k=x`$, while the effective HWHM $`{}_{}{}^{63}\mathrm{\Delta }x_{patch}`$ (defined in an analogous manner to Eq. (19)) is found to be independent of which site $`k=(A,B)`$ is computed. By deducing $`{}_{}{}^{63}\mathrm{\Delta }x_{patch}`$ over a large range of the free parameters $`\left[x,^{63}R_{patch}\right]`$, we find our calculated values of $`{}_{}{}^{63}\mathrm{\Delta }x_{patch}`$ fit well to the Binomial theorem

$`{}_{}{}^{63}\mathrm{\Delta }x_{patch}`$ $`=`$ $`\sqrt{\alpha _0}\sqrt{{\displaystyle \frac{x(1x)}{N^{R_{patch}}}}}`$ (30)
$`=`$ $`{\displaystyle \frac{a}{{}_{}{}^{63}R_{patch}^{}}}\sqrt{x(1x)}\sqrt{{\displaystyle \frac{\alpha _0}{\pi }}},`$
where we have used the relation in Eq. (23) for the second line. In Fig. 9 we show the temperature dependence of $`{}_{}{}^{63}x_{patch}^{}`$ using Eq. (2) and Eq. (30) for $`{}_{}{}^{63}R_{patch}^{}`$ defined in Eq. (29). Using $`{}_{}{}^{63}R_{patch}^{d_1}`$($`{}_{}{}^{63}R_{patch}^{d_2}`$) instead of Eq. (29) results in a small $`(+)7`$ % variation in $`{}_{}{}^{63}\mathrm{\Delta }x_{patch}`$, respectively. The fact that $`{}_{}{}^{63}x_{patch}^{}^{63}x_{local}`$ suggests that (a) our patch-by-patch model is appropriate in describing the spatial variation in $`{}_{}{}^{63}x_{local}^{}`$, (b) our estimate of the intrinsic lattice broadening $`{}_{}{}^{63}\mathrm{\Delta }\nu _{latt}`$ is correct in the temperature region above $`T_Q`$, and (c) the static value $`{}_{}{}^{63}x_{patch}^{}`$ deduced from the lineshape analysis is consistent with the low frequency dynamic value $`{}_{}{}^{63}x_{local}^{}`$ deduced using $`{}_{}{}^{63}1/T_1`$.

## V Orthorhombic Lattice Distortions

### V.1 Local Order Parameter

In the orthorhombic phase, the CuO6 octahedra are tilted by an angle $`\theta _c5^o`$ radaelli  away from the $`c`$-axis towards the $`b_o`$-axis, and the primitive lattice constants in the plane are split ($`b_o>a_o`$). We can calculate the effect of these distortions to the EFG by placing the ions in an orthorhombic cell (see Fig. 14) with appropriate orthorhombic modulations ginsberg ; radaelli  and compute the changes to the $`{}_{}{}^{63}\nu _{Q}^{}`$. However, we already know from Fig. 11 that $`{}_{}{}^{63}\nu _{Q}^{}`$ starts to increase above the structural phase transition temperature $`T_{st}`$ deduced from neutron and X-ray diffraction results radaelli . These bulk probes detect LRO, whereas NQR is a local probe, therefore it is more appropriate to deduce the local orthorhombic distortions from our $`{}_{}{}^{63}\nu _{Q}^{}`$ data, and then compare them with LRO results.

The orthorhombic distortions cannot a priori be deduced from the NQR data alone. We must first take certain know relations from LRO results and apply them to NQR. In particular, we shall use the neutron diffraction results reported by P.G. Radaelli et al radaelli  which demonstrate that certain measured orthorhombic parameters such as the orthorhombic splitting $`(b_oa_o)`$, the octahedron tilting angle $`\theta _c(5^o`$) which measures the departure angle of the Oap-Cu-Oap bond from the $`c`$-axis, and the “scissors” angle $`\theta _{scis}(90^o`$) which measures the Opl-Cu-Opl angle, are all related to powers of the underlying orthorhombic order parameter we define as $`\theta _{LRO}(x,T)`$.

At low temperatures ($`70`$ K) where lattice fluctuations are small, mean-field theory gives a good account of the observed lattice distortions radaelli . In light of this, we fit the data in Ref. radaelli  to a mean-field like form and deduced

$`\theta _c`$ $``$ $`\theta _{LRO}(x,T)`$
$`\theta _{scis}`$ $``$ $`90410^2\left[\theta _{LRO}(x,T)\right]^2`$
$`(b_oa_o)`$ $``$ $`610^4{\displaystyle \frac{(b_o+a_o)}{2}}\left[\theta _{LRO}(x,T)\right]^2`$ (31)
where all angles are in degrees. These set of equations provide a constraint for the orthorhombic distortions which determine the positions $`𝐫_j(x,T)`$ of the $`j`$th ion in the lattice. We now assume that the LRO constraints in Eq. (31) also apply to the $`local`$ constraints which determine the local position $`𝐫_j(x,T)`$ of the $`j`$th ion in the lattice (i.e. $`\theta _{LRO}(x,T)\theta _{local}(x,T)`$ in Eq. (31)).

The EFG lattice summation in Eq. (15) is now carried out in a similar way as described in section IV but with the orthorhombic distortions incorporated into the positions of the ions. The positions $`𝐫_j(x,T)`$ of the ions are deduced using the absolute values of the lattice constants $`a_o`$ and $`c`$ at 295 K from Fig. 4 along with their corresponding thermal coefficients $`\alpha _a`$ and $`\alpha _c`$ as before, however, we allow $`b_o`$ to vary according to the constraints in Eq. (31) and the unknown local structural parameter $`\theta _{local}(x,T)`$. We position the La/Sr ions and the apical oxygen ions according to the tilting angle $`\theta _c=\theta _{local}(x,T)`$ towards the $`b_o`$ direction, while the planar oxygens are positioned according to $`\theta _{sc}`$ in Eq. (31) and $`\theta _{local}(x,T)`$. We now have all the ionic positions $`𝐫_j(x,T)`$ as a function of $`\theta _{local}(x,T)`$. Just as in section IV, we then compute the EFG lattice summation and resonance frequency $`{}_{}{}^{\kappa }\nu _{}^{k}`$ (Eq. (21)) for the random configuration $`\kappa `$ and patch radius $`{}_{}{}^{63}R_{patch}^{}`$, however, the EFG lattice summation is taken using the tilted ionic positions specified by $`\theta _{local}(x,T)`$. We store $`{}_{}{}^{\kappa }\nu _{}^{k}`$ into a vector $`P_\theta ^k`$ and re-randomize the lattice until $`P_\theta ^k`$ is sufficiently large to build a histogram spectrum. We define the new CG $`P_\theta ^k`$ incorporating finite $`\theta _{local}`$ as $`\nu _\theta ^k`$.

The results of the simulation suggest that the increase in the CG due to $`\theta _{local}`$ $`alone`$ goes as

$$\nu _\theta ^k\nu _{\theta =0}^kc^k\theta _{local}^2,$$
(32)
where $`c^A`$ \+ 0.14 MHz/\[deg.2\] and $`c^B(\frac{\zeta _{latt}^B}{\zeta _{latt}^A}c^A)`$ \+ 0.08 MHz/\[deg.2\] are $`x`$ and $`T`$ independent within $`\pm 15`$ %. Note that NQR does not distinguish between $`\pm \theta _{local}`$, but only the average square magnitude $`\theta _{local}^2`$ across the whole sample. Another conclusion we deduce from the calculation is that after diagonalizing the EFG tensor, the main principle axis is found to be tilted an angle $`|\theta _{NQR}|`$ away from the $`c`$-axis (where $`\theta _{NQR}0.9\theta _{local}`$), while the EFG asymmetry parameter $`\eta `$ remains small $`\eta <0.06`$.

We are now in a position to deduce the average $`\theta _{local}^2`$ from our experimental data $`{}_{}{}^{63}\nu _{Q}^{k}`$ (Fig. 11). In the spirit of the LRO results, we deduced the temperature dependence of $`\theta _{local}^2`$ (presented in Fig. 22) from the $`{}_{}{}^{63}\nu _{Q}^{k}`$ data assuming the paramteric form

$$\theta _{local}^2=\theta _{local}^2^0\left(\frac{T_QT}{T_Q}\right)^{2\beta _{local}},$$
(33)
where $`\theta _{local}^2^0`$ is the square averaged tilting angle at $`T=0`$ K, $`\beta _{local}`$ is the reduced temperature coefficient, and $`T_Q`$ is local onset temperature for the orhorhombic distortions. All three free parameters in Eq. (33) depend on $`x`$ and are optimized to fit $`{}_{}{}^{63}\nu _{Q}^{k}`$ (in Fig. 11). The best fit to the $`{}_{}{}^{63}\nu _{Q}^{k}`$ are shown as the black and grey curves in Fig. 11, and the corresponding temperature dependence of $`\theta _{local}^2`$ is shown in Fig. 22. In Fig. 23 we summarize the $`x`$ dependence of best fit values $`\theta _{local}^2^0`$, $`\beta _{local}`$ and $`T_Q`$. We compare these parameters with equivalent mean field parameters $`\theta _{LRO}^0^2`$, $`\beta `$ and $`T_{st}`$ according to LRO radaelli ; yamada , shown as the solid grey lines in Fig. 23.

In the case of $`x=0.0`$, all three parameters deduced by NQR are consistent with equivalent LRO parameters deduced by neutron and X-ray diffraction. We find that $`\beta _{local}0.34`$ according to NQR for $`x=0.0`$, which is consistent with the LRO critical exponent $`\beta 0.30`$ according to neutron and X-ray scattering for all $`x`$ beta ; ginsberg . A similar critical exponent is also found in La2-xBaxCuO4 where $`\beta =0.33`$ suzuki , which are all consistent with the theoretical predictions from a 3d $`XY`$ model.

For $`x0.04`$ we find consistent values between $`\theta _{local}^2^0`$ and $`\theta _{LRO}^0^2`$. In the case of $`\beta _{local}`$ and $`T_Q`$ for $`x>0.04`$, however, we find qualitative differences with the LRO parameters $`\beta `$ and $`T_{st}`$, respectively. The fact that $`\beta _{local}1`$ for $`x>0.04`$ indicates that locally, the transition into the orthorhombic phase is smeared out in temperature. The fact that $`T_Q400`$ K indicates that the local onset temperature remains large and roughly $`x`$ independent.

In order to account for the differences between local and LRO structural results one must first recall that our NQR results measure the local square deviation $`\theta _{local}^2`$, while the diffraction results measure the coherent average $`\theta _{LRO}`$ across the sample. One must also leave open the possibility that the orthorhombic structural transition for $`x>0.04`$ is partly order-disorder and partly displacive in nature. In an order-disorder dominated structural transition (as proposed by other local probes bozin ; haskel  in the region $`x>0.15`$), $`\theta _{local}`$ is argued to be finite and temperature independent in the LRO tetragonal phase ($`T>T_{st}`$), while $`\theta _{LRO}=0`$ when averaged over large ($`10`$ nm) length scales. In such a scenario, one can picture certain short length scale regions ($`10`$ nm) which align along the positive $`b_o`$ direction where $`\theta _{local}>0`$, while an equal number of other regions align along the negative $`b_o`$ direction with $`\theta _{local}<0`$, therefore $`\theta _{LRO}=0`$. The temperature dependence in an order-disorder type transition does not involve a $`T`$ dependence in the magnitude of $`|\theta _{local}|`$ itself. Instead, the $`T`$ dependence lies in the coherent orientation of the local regions towards a particular direction, for arguments sake say along positive $`b_o`$ direction. The coherent average of the first moment $`\theta _{LRO}`$ will then be finite and positive ($`\theta _{LRO}>0`$), and will appear to increase with decreasing temperature as more and more local regions in the sample choose to align along the positive $`b_o`$ direction.

In a displacive dominated structural transition (as argued for $`x0.04`$), the magnitude of $`|\theta _{local}|`$ itself increases with decreasing temperature below $`T_{st}`$ and $`all`$ the local regions in the sample align along the same $`b_o`$ direction. Our new NQR results combined with LRO results support a displacive structural transition for $`x0.04`$, while a mixture of displacive and order-disorder type transition occurs in the region $`x>0.04`$. For $`x>0.04`$ below $`T<T_Q`$, NQR indicates that locally the octahedra start to tilt in a displacive manner (i.e. the tilting angle of each octahedra increases with decreasing temperature below $`T<T_Q`$ and $`\theta _{local}^2>0`$). In the region $`T_{st}<T<T_Q`$, $`\theta _{local}`$ continues to increase yet there continue to be an equal number of regions with positive and negative tilting with respect to the $`b_o`$ direction, and therefore $`\theta _{LRO}=0`$. In the region $`T<T_{st}`$, $`\theta _{local}`$ continues to increase with decreasing temperature, however, the local regions start to coherently align along one direction (say along the positive $`b_o`$ direction), and $`\theta _{LRO}>0`$. This interpretation is entirely consistent with other local probes bozin ; haskel , but in addition we have deduced the onset temperature $`T_Q`$ for the local orthorhombic distortions. In the case of $`x=0.20`$ we find that $`T_Q350`$ K even though the LRO structural transition disappears yamada .

### V.2 Distribution in Local Parameter

In all of our analysis of the 63Cu NQR linebroadening thus far, we have ignored potential linebroadening from local lattice distortions. In the previous part of this section, we deduced the temperature dependence of the square $`average`$ distortions $`\theta _{local}^2`$ across the sample using the CG data $`{}_{}{}^{63}\nu _{Q}^{A}`$. We shall now estimate the corresponding $`distribution`$ of $`\theta _{local}`$ using the HWHM of the NQR spectrum, $`{}_{}{}^{63}\mathrm{\Delta }\nu _Q^A`$. As discussed earlier, the three linebroadening mechanisms for the 63Cu NQR lineshape are (a) the random placement of +2Sr ions in the lattice (b) the patch-by-patch variation in hole concentration $`{}_{}{}^{63}x_{local}^{}`$, and (c) the local lattice distortions. In section IV we used mechanisms (a) and (b) to account for the experimental data, and deduced the lower bound $`{}_{}{}^{63}R_{patch}^{}`$ for the spatial variation $`{}_{}{}^{63}x_{local}^{}`$. We shall now use all three mechanisms (a),(b) and (c) to account for the observed HWHM data $`{}_{}{}^{63}\mathrm{\Delta }\nu _Q^A`$ in the case of $`x=0.04`$. Note that we cannot a priori separate (a),(b) and (c), therefore in order to obtain the contribution from (c), we must first make reasonable assumptions for (a) and (b).

In Fig. 24 we plot the temperature dependence of the observed data $`{}_{}{}^{63}\mathrm{\Delta }\nu _Q^A`$ for $`x=0.04`$ above and below $`T_QT_{st}450`$ K taken from Fig. 12. For $`T>T_Q`$, we attribute all of the observed linebroadening to mechanism (a) and (b) given by $`\mathrm{\Delta }\nu _{R_{patch}}^A`$. For $`T<T_Q`$, however, we make the ansatz that the $`\mathrm{\Delta }\nu _{R_{patch}^{}}^A`$ continues to increase only linearly, as shown in Fig. 24. $`\mathrm{\Delta }\nu _{R_{patch}^{}}^A`$ in Fig. 24 results from using a different estimate $`{}_{}{}^{63}R_{patch}^{}`$ for the patch radius, shown as the grey dashed line in Fig. 21 for $`x=0.04`$. The motivation for using $`{}_{}{}^{63}R_{patch}^{}`$ for $`x=0.04`$ is that the temperature dependence follows that of $`x0.115`$ where the overall amplitude of the orthorhombic distortions ($`\theta _{local}^2^0`$ in Eq. (33)) are smaller.

In order to simulated the NQR spectrum, we now proceed through the simulation as in section IV, however, we calculate the local hole concentration $`{}_{}{}^{63}x_{local}^{\kappa }`$ for a random configuration $`\kappa `$ using the patch radius $`{}_{}{}^{63}R_{patch}^{}`$ (Fig. 21), instead of $`{}_{}{}^{63}R_{patch}^{}`$. Along with each simulation run $`\kappa `$, we also independently choose a random value of $`|\theta _{local}^\kappa |`$ taken from a trial PDF (probability distribution function), examples of which are shown in Fig. 25. The positions of all the ions in the random lattice configuration $`\kappa `$ are calculated given the value $`|\theta _{local}^\kappa |`$ and Eq. (31) (with $`\theta _{LRO}\theta _{local}^\kappa `$). We then take the lattice summation in Eq. (15) using the the random positions of the ions with $`|\theta _{local}^\kappa |`$, together with the random local hole concentration $`{}_{}{}^{63}x_{local}^{\kappa }`$. The main principle value of the diagonalized resonance frequency $`{}_{}{}^{\kappa }\nu _{}^{k}`$ is then computed and stored into a vector $`P_\theta `$, then a new random lattice $`\kappa ^{}`$ is re-generated along with a new value of $`|\theta _{local}^\kappa ^{}|`$ according to the trial PDF. The corresponding resonance frequency $`{}_{}{}^{\kappa ^{}}\nu _{}^{k}`$ is calculated and stored into $`P_\theta `$. As in section IV, taking account of the random fluctuations of the main principle axis around the average NQR axis $`\theta _{NQR}`$ ($`0.9\theta _{local}`$) makes only minor $`2`$ % corrections to $`{}_{}{}^{\kappa }\nu _{}^{k}`$. As before, we repeat this proceedure $`10^4`$ times until $`P_\theta `$ is sufficiently large to produce a histogram NQR spectrum with $``$ 50 bins. Since $`{}_{}{}^{63}R_{patch}^{}`$ (see Fig. 21) is already determined, the remaining free parameter in the simulation is now the PDF for $`|\theta _{local}|`$.

The simulation is repeated using various trial PDF’s of $`|\theta _{local}|`$ until the the $`same`$ simulation of the NQR spectrum using $`{}_{}{}^{63}R_{patch}^{}`$ from section IV is obtained. A good starting point for the PDF turns out to be a Gaussian distribution in the squared variable $`\theta _{local}^2`$ centered at the values $`\theta _{local}^2`$ taken from Fig. 22 for $`x=0.04`$. The width of the Gaussian distribution can then be optimized to best reproduce the observed NQR spectrum. After various trials, the final results of the optimized PDF’s for $`\left|\theta _{local}\right|(\sqrt{\theta _{local}^2})`$ are shown in Fig. 25 at various temperatures $`TT_Q`$.

## VI Conclusions

We have presented a detailed systematic study of the spatial variation in electronic states using 63Cu NQR in 63Cu isotope enriched poly-crystalline samples of La2-xSrxCuO4 for $`0.04x0.16`$. By analysing the extent of the frequency dependence of $`{}_{}{}^{63}1/T_1`$ across the inhomogeneous linebroadening of the 63Cu NQR spectrum, we determined the spatial variation in hole concentration $`{}_{}{}^{63}x_{local}^{}`$ given by $`{}_{}{}^{63}x_{local}^{}=x\pm ^{63}\mathrm{\Delta }x_{local}`$, where $`{}_{}{}^{63}\mathrm{\Delta }x_{local}`$ was defined as the characteristic amplitude or extent of the spatial variaiton away from nominal hole concentration $`x`$. We showed that $`{}_{}{}^{63}x_{local}^{}(x)`$ is a thermodynamic effect whose extent increases with decreasing temperature below 500-600 K and reaches values as large as $`{}_{}{}^{63}\mathrm{\Delta }x_{local}/x0.5`$ in the temperature region $`150`$ K. Furthermore, we showed that the spatial variation $`{}_{}{}^{63}x_{local}^{}`$ is an $`intrinsic`$ effect in La2-xSrxCuO4 by confirming the $`same`$ extent of the inhomogeneity in our poly-crystalline samples to that in high-quality single crystals tobe . To the best of our knowledge these reports, together with those in Ref.’s singerprl ; tobe , are the first of its kind to show the temperature dependence of the intrinsic inhomogeneity in La2-xSrxCuO4, or any other system with quenched disorder. In Ref. tobe  the results of the inhomogeneous electronic state $`{}_{}{}^{17}x_{local}^{}`$ determined using 17O NMR in high-quality La2-xSrxCuO4 crystals are shown to be consistent with the 63Cu NQR result $`{}_{}{}^{63}x_{local}^{}`$. We believe that this is a powerful demonstration of the validity of our analysis in this work, as well as in singerprl ; tobe . Moreover, consistency between 17O NMR and 63Cu NQR results indicate that the spatial variation in the spin and charge channels are correlated across most of the Brillouin zone.

We have argued that there are two essential ingredients for the spatial variation $`{}_{}{}^{63}x_{local}^{}(x)`$. They are as follows: (a) that the +2Sr donor ions are positioned with random probability throughout the lattice and (b) that there is a short length scale for the spatial variation $`{}_{}{}^{63}x_{local}^{}`$ of the order $`{}_{}{}^{63}R_{patch}^{}3.04.0`$ nm. Using (a) and (b) above we succesfully fit the inhomogeneous 63Cu NQR spectrum (including the B-line) using a patch-by-patch model for the spatial variation $`{}_{}{}^{63}x_{local}^{}`$ with the patch radius $`{}_{}{}^{63}R_{patch}^{}`$ as the only adjustable parameter. We assumed that no extra 63Cu NQR linebroadening from the lattice distortions existed, therefore the calculation for $`{}_{}{}^{63}R_{patch}^{}`$ resulted in a lower bound. As shown in Ref. tobe , however, our lower bound estimate $`{}_{}{}^{63}R_{patch}^{}`$ is shown to be consistent with $`{}_{}{}^{17}R_{patch}^{}`$ deduced from the inhomogeneous 17O NMR spectrum, which justifies our assumptions used to determine $`{}_{}{}^{63}R_{patch}^{}`$ in the temperature region $`T>T_Q`$. Within the EFG simulation, we also deduced an upper boundary to the spatial variation in electronic states $`{}_{}{}^{63}x_{patch}^{}`$ given by $`{}_{}{}^{63}x_{patch}^{}=x\pm ^{63}\mathrm{\Delta }x_{patch}`$, where the extent of the spatial variation $`{}_{}{}^{63}\mathrm{\Delta }x_{patch}`$ was shown to depend on $`{}_{}{}^{63}R_{patch}^{}`$ as such $`{}_{}{}^{63}\mathrm{\Delta }x_{patch}1/^{63}R_{patch}`$ (i.e. $`{}_{}{}^{63}R_{patch}^{}=\mathrm{}`$ implied the homogenous doping limit $`{}_{}{}^{63}x_{patch}^{}=x`$). Within the temperature region $`T>T_Q`$ we found the consistent result $`{}_{}{}^{63}x_{patch}^{}^{63}x_{local}`$, which justified our EFG simulation and short length scale model of the inhomogeneous electronic state.

$`{}_{}{}^{63}R_{patch}^{}`$ was deduced assuming no clustering of the +2Sr ions in the lattice beyond probability theory, however, we note that even in the case of clustering over length scales 3-4 nm, the presence of a short length scale is still needed to account for the observed spatial variation in electronic states. We also argued that the short length scale implied that any experimental probe which averages over spatial extents larger than $`{}_{}{}^{63}R_{patch}^{}3.04.0`$ nm results in an $`underestimate`$ of the extent of the inhomogeneous electronic state. NQR, however, is a local probe, which makes it $`ideal`$ for detecting short $``$nm variations in the electronic state. Our findings have a strong impact on the general, but now incorrect, view that doped holes are homogeneously distributed in the CuO2 planes of La2-xSrxCuO4.

We have also compared our local orthorhombic structural distortions results to LRO structural results taken from neutron and X-ray diffraction radaelli ; yamada  techniques which probe coherence phenomena over length scales larger than ten’s of nm takagidelta . In the case of $`x=0.0`$ we found good agreement between NQR and LRO results, however, for $`x0.04`$ we presented evidence that the local orthorhombic distortions started at elevated temperatures $`T_Q(400`$ K) $`above`$ the onset of LRO $`T_{st}`$. We also found that the local distortions were somewhat smeared below $`T_Q`$ compared with $`x=0.0`$. Our local structural results were found to be consistent with a mixture of order-disorder and displacive transitions in the region $`x0.04`$, consistent with other local probes bozin ; haskel . Our NQR studies at elevated temperatures have now quantified the onset of the local orthorhombic distortions at $`T_Q(400`$ K) in the region $`0.04x0.20`$.

It is known that in the $`x<1/8`$ region of the phase diagram, the temperature dependence of the wipeout fraction $`{}_{}{}^{63}F(T)`$ (defined as the fracton of unobservable Cu nuclei that are wiped-out) is somewhat tailed hunt ; hunt2 ; singerprl  below $`T_{NQR}`$ in both La2-xSrxCuO4 and its Eu and Nd co-doped materials. In the case of La2-xSrxCuO4 for $`x=0.04`$, the onset temperature for Cu wipeout is rather high at $`T_{NQR}300`$ K, yet half the Cu nuclei are still observable (i.e. $`{}_{}{}^{63}F(T)1/2`$) at $`T_{NQR}/2`$, implying that $`T_{NQR}`$ does not signify the onset of a $`collective`$ phenomenon for $`x<1/8`$. This may indicate the formation of a small number of local moments in the CuO2 plane which cause the wipeout of the NQR signal for neighbouring Cu sites, as in a classic spin-glass.

Rather than invoking a spin-glass type model, however, we can now naturally account for the tailed wipeout for $`x<1/8`$ using the temperature dependence of the growing spatial variation in electronic states $`{}_{}{}^{63}x_{local}^{}`$. In the temperature region $`T_{NQR}300`$ K for $`x=0.04`$, the characteristic local hole concentration $`{}_{}{}^{63}x_{local}^{}`$ for insulating patches is inferred to be as low as $`{}_{}{}^{63}x_{local}^{}0.015`$ (Fig. 9). Recall that $`{}_{}{}^{63}1/T_1`$ for $`x=0.02`$ at $`300`$ K is already too large for the Cu signal to be observable, therefore the more insulating patches for $`x=0.04`$ at 300 K begin to wipe-out. As one cools the $`x=0.04`$ sample below 300 K, there are a growing number of patches with lower and lower hole concentration which are also wiped-out, leading to an gradual increase in the wipeout fraction $`{}_{}{}^{63}F(T)`$ with decreasing temperature. The fact that $`{}_{}{}^{63}x_{local}^{}`$ only gradually fans out with decreasing temperature naturally explains the tailed nature of the wipeout fraction $`{}_{}{}^{63}F(T)`$ below $`T_{NQR}`$ in the underdoped region $`x<1/8`$.

In Eu and Nd co-doped La2-xSrxCuO4 for $`x<1/8`$, the same tailed feature in $`{}_{}{}^{63}F(T)`$ is found together with an inflection point $`T_{NQR}^{inflec}(<T_{NQR})`$ below which $`{}_{}{}^{63}F(T)`$ increases rapidly with decreasing temperature. In the case of Eu and Nd co-doping with $`x<1/8`$, $`T_{NQR}^{inflec}`$ coincides with the onset temperature $`T_{LTT}`$ of an additional structural phase transition with tetragonal symmetry known as the LTT (low temperature tetragonal) phase $`and`$ the onset of charge ordering $`T_{charge}`$. In the case of La1.6-xNd0.4SrxCuO4, $`T_{charge}T_{LTT}70`$ K ichi . This coincidence suggests that the glassy slowing of the stripe modulation below $`T_{NQR}^{inflec}(T_{charge})`$ is accelerated by the tetragonal symmetry compatible with the charge stripe symmetry tranquada . 63Cu isotope enriched samples of Eu and Nd co-doped La2-xSrxCuO4 are found to have comparable frequency dependence of $`{}_{}{}^{63}1/T_1`$ across its inhomogeneous 63Cu NQR spectrum. In the Eu and Nd co-doped cases, therefore, we can naturally account for the tailed wipeout in the temperature region $`T_{NQR}^{inflec}<T<T_{NQR}`$ by the growing inhomogeneity in electronic states with decreasing temperature, while for $`T<T_{NQR}^{inflec}(T_{charge}`$), the wipeout is dramatically accelerated due to the pinning of charge stripes by the lattice symmetry. Although there is no LRO into a LTT phase in La2-xSrxCuO4 which can stablize the charge stripes, there is evidence for an incipient structural transition into the LTT phase kimura , which may suggest a similar phenomenon is occuring in La2-xSrxCuO4.

## Appendix A Anti-shielding factors

In this appendix we deduce the anti-shielding factors slichter  $`\zeta _{latt}^k`$ and $`\zeta _{3d}^k`$, where $`k=(A,B)`$. We first assume that both anti-shielding factors are isotropic and independent of $`x`$ and $`T`$, but are allowed to vary between the A and B-lines. We produce a cross-plot (shown in Fig. 26(a)) of the resonance frequency at the CG of the A-line $`{}_{}{}^{63}\nu _{Q}^{A}(x,T)`$ (Fig. 11(a)) against the CG of the lattice contribution to the EFG $`V_{latt}^{c,A}(x,T)`$ (Fig. 15). Likewise, we produce the same cross-plot for the B-line shown Fig. 26(b). Note that we only use data in temperature region above $`T_Q`$ (Fig. 11) where local orthorhombic distortions are absent.

Then we deduce the linear fit using the standard empirical form penn ; shimizu :

$${}_{}{}^{63}\nu _{Q}^{k}=\mathrm{\Lambda }_{3d}^k\mathrm{\Lambda }_{latt}^kV_{latt}^{c,k}$$
(34)
where $`k=(A,B)`$ and all quantities are overall positive. The best fit to the data shown by the solid black lines in Fig. 26 indicates that $`\mathrm{\Lambda }_{3d}^A`$ = 75.3 MHz and $`\mathrm{\Lambda }_{latt}^A=16.0`$ MHz/\[emu$`\times 10^{14}`$\] for the A-line, and that $`\mathrm{\Lambda }_{3d}^B`$ = 64.6 MHz and $`\mathrm{\Lambda }_{latt}^B=10.7`$ MHz/\[emu$`\times 10^{14}`$\] for the B-line.

Using Eq. (20) we express the $`\mathrm{\Lambda }`$’s in terms of the unknown anti-shielding factors as such

$`\mathrm{\Lambda }_{3d}^k`$ $`=`$ $`{\displaystyle \frac{4e^2|^{63}Q|}{14h}}\left(14f_\sigma ^o\right)r_{3d}^3\zeta _{3d}^k`$
$`\mathrm{\Lambda }_{latt}^k`$ $`=`$ $`{\displaystyle \frac{e|^{63}Q|}{2h}}\zeta _{latt}^k`$ (35)
Using $`f_\sigma ^o=0.076`$ tobe  , $`{}_{}{}^{63}Q=0.16`$ barns bleaney , and the bare value $`r_{3d}^3=7.5`$ a.u. bleaney  for free cupric ions, we find that for the A-line $`\zeta _{latt}^A=27.7`$ and $`\zeta _{3d}^A=0.93`$, while for the B-line $`\zeta _{latt}^B=18.6`$ and $`\zeta _{3d}^B=0.82`$. These set of results are consistent with previous estimates penn ; shimizu .

Converting the above analysis using $`{}_{}{}^{63}Q=0.211`$ barns stern  instead of $`{}_{}{}^{63}Q=0.16`$ barns bleaney  results in a uniform 24.2 % decrease in our estimates of the $`\zeta `$’s. However, note that what really matters in converting from calculation to experiment are the $`\mathrm{\Lambda }`$’s in Eq. (34). Likewise, when accounting for the linebroadening in section IV to deduce $`{}_{}{}^{63}\mathrm{\Delta }x_{patch}`$, it is also the $`\mathrm{\Lambda }`$’s which ultimately matter, not the $`\zeta `$’s, therefore our end result for $`{}_{}{}^{63}\mathrm{\Delta }x_{patch}`$ and $`{}_{}{}^{63}R_{patch}^{}`$ are independent of what value of $`{}_{}{}^{63}Q`$ we use.

## Appendix B Magnetization Recovery

In this appendix we quantify the observed multi-exponential recovery of the magnetization $`M(t)`$ in the case of $`x=0.16`$. The fact that the recovery is multi-exponential implies that $`{}_{}{}^{63}1/T_1`$ is distributed at a fixed frequency on the NQR line. We can naturally account for the distribution in $`{}_{}{}^{63}1/T_1`$ at fixed frequency by recalling that the 63Cu NQR line has an underlying intrinsic lattice broadening ($`\mathrm{\Delta }\nu _{latt}^k`$) which will tend to smear out the frequency dependence of $`{}_{}{}^{63}1/T_1`$ across the NQR line. Smearing out the frequency dependence of $`{}_{}{}^{63}1/T_1`$ results in a distribution of $`{}_{}{}^{63}1/T_1`$ values at a fixed frequency. Using the the multi-exponential recovery data, we estimate the extent of the underlying lattice broadening $`\mathrm{\Delta }\nu _{latt}^{T1}=0.62(\pm 0.07)`$ MHz in the case of $`x=0.16`$ at the A-line. We find that $`\mathrm{\Delta }\nu _{latt}^{T1}`$ is consistent with the lattice broadening $`\mathrm{\Delta }\nu _{latt}^A=0.49`$ MHz found independently using our point charge calculation (section IV), which also adds weight to our EFG simulation of the random lattice.

The experimentally observed nuclear recovery $`M(t)`$ can be fit to Eq. (11), where $`M(0)`$, $`M(\mathrm{})`$ and $`{}_{}{}^{63}1/T_1`$ are free parameters of the fit. The bare recovery is then defined as $`I(t)`$ where

$$I(t)=\frac{M(t)M(\mathrm{})}{M(0)M(\mathrm{})}.$$
(36)
In the case of single exponential recovery, $`I(t)`$ is then found to be a straight line on a semi-log plot, an example of which is shown in Fig. 27 where we present $`I(t)`$ in the case $`x=0.0`$ at 475 K for various positions across its NQR spectrum. As discussed earlier, the $`x=0.0`$ spectrum is homogeneously broadened, and therefore shows no frequency dependence of $`{}_{}{}^{63}1/T_1`$ across its line (i.e. $`I(t)`$ has the same slope at different frequencies), nor does it show any sign of having a distribution at a fixed frequency. We also note that the data in Fig. 27 is taken at 475 K which is the orthorhombic phase, therefore we can also rule out the possibility that the orthorhombic distortions are the $`direct`$ cause of observed frequency dependence in $`{}_{}{}^{63}1/T_1`$.

In Fig. 28 we show the bare recovery $`I(t)`$ at similar positions across the A-line for $`x=0.16`$ at 100 K. First note that the slopes in the recovery get steeper towards lower frequency. This corresponds to the frequency dependence in $`{}_{}{}^{63}1/T_1`$ where $`{}_{}{}^{63}1/T_{1,A}^{()}>^{63}1/T_{1,A}^{(0)}>^{63}1/T_{1,A}^{(+)}`$. Second, note that $`I(t)`$ shows finite curvature on a semi-log plot. This implies that there is a distribution in recoveries $`I(t)`$ which can be fit to the general form

$$I(t)=\frac{1}{_ja_j}\underset{j}{}a_j\mathrm{exp}\left(\frac{3}{T_{1,j}}t\right).$$
(37)
Fitting the recoveries to a single exponential form (i.e. taking $`a_j`$ =1 and $`1/T_{1,j}=1/T_1`$ in Eq. (37)) results in the force fit values $`{}_{}{}^{63}1/T_{1,A}^{()}`$, $`{}_{}{}^{63}1/T_{1,A}^{(0)}`$ and $`{}_{}{}^{63}1/T_{1,A}^{(+)}`$. These force fit values represent the average of the underlying distribution in $`{}_{}{}^{63}1/T_{1,A}`$.

The dashed lines in Fig. 28 show the best fit result from a sum of exponentials $`1/T_{1,j}`$ whose distribution coefficients $`a_j`$ are presented as the dashed lines in Fig. 29. The distribution coefficients $`a_j`$ in Fig. 29 are determined as such

$`a_j^{()}`$ $`=`$ $`\mathrm{exp}\left({\displaystyle \frac{(P_{T1}1/T_{1,j})^2}{(\mathrm{\Delta }_P^{()})^2}}\right)`$ (38)
$`a_j^{(+)}`$ $`=`$ $`\mathrm{exp}\left({\displaystyle \frac{(P_{T1}1/T_{1,j})^2}{(\mathrm{\Delta }_P^{(+)})^2}}\right),`$ (39)
where Eq. (38) is used for $`1/T_{1,j}<P_{T1}`$ and Eq. (39) is used for $`1/T_{1,j}>P_{T1}`$. The three free parameters used to fit $`I(t)`$ are therefore the peak of the distribution $`P_{T1}`$ and the two widths $`\mathrm{\Delta }_P^{()}`$ and $`\mathrm{\Delta }_P^{(+)}`$. When optimizing the three parameters in Eq. (39), we also use the constraint that the CG of the distribution (shown as the solid vertical lines in Fig. 29) coincide with the force-fit value $`{}_{}{}^{63}1/T_{1,A}`$ (presented in Fig. 30) deduced with a single exponential fit.

We then deduce the second moment of the distributions in Fig. 29 in an analogous manner to Eq. (12) and Eq. (19), and we find that the HWHM of the distributions are $`\mathrm{\Delta }1/T_{1,A}^{()}=0.33`$ (ms)-1, $`\mathrm{\Delta }1/T_{1,A}^{(0)}=0.27`$ (ms)-1 and $`\mathrm{\Delta }1/T_{1,A}^{(+})=0.20`$ (ms)-1. Note that the HWHM of the distributions are $``$ 15 % of average value $`{}_{}{}^{63}1/T_{1,A}`$, which implies that accurate recovery data $`I(t)`$ is needed down to $`I(t)0.001`$.

We can now get an estimate of the underlying intrinsic lattice broadening $`\mathrm{\Delta }\nu _{latt}^{T1}`$ by converting the HWHM $`\mathrm{\Delta }1/T_{1,A}`$ deduced above into equivalent frequency widths. In order to do so we take a local derivative $`\beta _\nu `$ of the underlying frequency dependence of $`{}_{}{}^{63}1/T_{1,A}`$ at a frequency $`\nu `$ defined as such

$$\beta ^\nu =\left|\frac{\delta 1/T_{1,A}}{\delta \nu }\right|_\nu ,$$
(40)
where $`1/T_{1,A}`$ is the interpolated value of $`{}_{}{}^{63}1/T_{1,A}`$ shown in Fig. 30 as the solid line. Using $`{}_{}{}^{63}1/T_{1,A}`$ in units of (ms)-1 and $`\nu `$ in units of MHz, we find that $`\beta ^{()}=0.62`$, $`\beta ^{(0)}=0.39`$, and $`\beta ^{(+)}=0.29`$, from which we deduce the HWHM (in MHz)

$`\mathrm{\Delta }\nu _{latt}^{()}`$ $`=`$ $`\mathrm{\Delta }1/T_{1,A}^{()}/\beta ^{()}=0.53`$
$`\mathrm{\Delta }\nu _{latt}^{(0)}`$ $`=`$ $`\mathrm{\Delta }1/T_{1,A}^{(0)}/\beta ^{(0)}=0.69`$
$`\mathrm{\Delta }\nu _{latt}^{(+)}`$ $`=`$ $`\mathrm{\Delta }1/T_{1,A}^{(+)}/\beta ^{(+)}=0.68`$ (41)
for the lower half intensity, center and upper half intensity of the spectrum, respectively. Finally, we deduce the average value $`\mathrm{\Delta }\nu _{latt}^{T1}`$ of the three estimates in Eq. (41) to be $`\mathrm{\Delta }\nu _{latt}^{T1}=0.62(\pm 0.07)`$ MHz. This is consistent with $`\mathrm{\Delta }\nu _{latt}^A=0.49`$ MHz deduced independently using the point charge calculation, as shown in Fig. 30.

We note that our whole analysis in this appendix is derived using fixed $`\tau =12`$ $`\mu `$s conditions. As discussed in section II and III, however, there exists a small $`\tau `$ dependence on the underlying distribution in the spin-lattice relaxation rate. Such effects could account for the 25 % difference between the lattice width $`\mathrm{\Delta }\nu _{latt}^{T1}`$ derived from Eq. (41) and the point charge calculation $`\mathrm{\Delta }\nu _{latt}^k`$ in section IV.

###### Acknowledgements.

The work at M.I.T. was supported by NSF DMR 98-08941 and 99-71264.