text:
# Introduction

## Introduction

Present data from B factories on inclusive decays already require precise knowledge of the bottom quark mass parameter with a numerical precision of the order $`50`$ MeV with a reliable estimate of the uncertainty. This will become even more acute in the future when more data becomes available. Using methods based on perturbative QCD, there have been several approaches in the past aiming at uncertainties of less than $`100`$ MeV. The most frequently used method is based on large-$`n`$ ($`n>4`$) moments of the $`b\overline{b}`$ production cross section in $`e^+e^{}`$ annihilation ,

$$P_n=\frac{ds}{s^{n+1}}R_{b\overline{b}}(s),$$
(1)
where $`R_{b\overline{b}}=\sigma (e^+e^{}b\overline{b}+X)/\sigma (e^+e^{}\mu ^+\mu ^{})`$ and the contributions from the virtual Z are neglected. From fitting moments obtained from experimental data to the corresponding theoretical expressions the bottom mass can be determined in threshold schemes (for recent reviews see Refs. ).111 Recently, the bottom mass has been determined in different threshold schemes from moments of inclusive semileptonic and radiative B meson decay spectra with an uncertainty of about $`100`$ MeV, which is dominated by experimental errors. (See also Ref. .)  Large-$`n`$ (“non-relativistic”) moments have the advantage that the badly known experimental continuum $`b\overline{b}`$ cross section above the $`\mathrm{{\rm Y}}(\text{6S})`$ is strongly suppressed in comparison to the rather well known resonance region. However, the order $`\alpha _s^2`$ corrections in the framework of the non-relativistic expansion are generally quite large and small uncertainties below $`100`$ MeV can only be achieved with additional assumptions on higher order corrections.  This indicates that an improvement in the treatment of the non-relativistic bottom quark dynamics might be needed.

A different method uses low-$`n`$ (“relativistic”) moments  where $`n<4`$. In the recent past, relativistic moments have been used less frequently because the badly known continuum region represents a major source of uncertainty that is not reducible without additional assumptions. Theoretically, the usual loop expansion in powers of $`\alpha _s`$ can be employed since for small $`n`$ the bottom dynamics is relativistic. Here, the $`\overline{\text{MS}}`$ mass is an appropriate mass definition to be used and extracted. In contrast to large-$`n`$ moments, the low-$`n`$ moments show a quite good perturbative behavior. A recent analysis by Kühn and Steinhauser  used moments at order $`\alpha _s^2`$. Adopting a theory-driven perspective, the experimental continuum data above the $`\mathrm{{\rm Y}}(\text{6S})`$ was obtained from theoretical results for $`R_{b\overline{b}}`$, basically eliminating uncertainties from the continuum region. In Ref.  the important conclusion was drawn that, using the strategy of Ref. , a substantially more accurate measurement of the $`\mathrm{{\rm Y}}(\text{1S})`$ \- $`\mathrm{{\rm Y}}(\text{6S})`$ region at CLEO  could result in an uncertainty in $`\overline{m}_b(\overline{m}_b)`$ of only $`30`$ MeV.

It is the main purpose of this paper to give a detailed compilation of all sources of uncertainties in the bottom $`\overline{\text{MS}}`$ mass $`\overline{m}_b(\overline{m}_b)`$ obtained from low-$`n`$ moments, including a more conservative treatment of the experimental continuum region. We believe that this compilation can contribute to a more differentiated view on the current uncertainties in $`\overline{m}_b(\overline{m}_b)`$ from low-$`n`$ moments and on the impact of new more precise data in the $`\mathrm{{\rm Y}}`$ resonance region from CLEO. In our analysis we also include the contributions from secondary $`b\overline{b}`$ production from gluon splitting and the effects of the non-zero charm quark mass, which have to our knowledge not been taken into account before. Both effects turn out to be small. Finally, we point out an inconsistency in the way the PDG has treated the original results for the electronic partial widths of $`\mathrm{{\rm Y}}(\text{4S})`$ and $`\mathrm{{\rm Y}}(\text{5S})`$ from CUSB  and CLEO , which leads to a contribution to the experimental moments $`P_n`$ that is smaller than the contributions obtained from the original data in that energy region, both from CUSB and CLEO.

## Theoretical Moments

For the QCD parameters used in this work we adopt the $`\overline{\text{MS}}`$ renormalization scheme and the convention that the bottom quark participates in the running ($`n_f=5`$). The masses of the quarks in the first two generations are set to zero. In terms of $`\overline{m}_b(\mu )`$ the moments in the OPE, including the known perturbative corrections to order $`\alpha _s^2`$ and the contribution from the dimension four gluon condensate, take the form

$`P_n`$ $`=`$ $`{\displaystyle \frac{1}{(\mathrm{\hspace{0.17em}4}\overline{m}_b(\mu ))^n}}\{f_n^0+\left({\displaystyle \frac{\alpha _s(\mu )}{\pi }}\right)(f_n^{10}+f_n^{11}\mathrm{ln}\left({\displaystyle \frac{\overline{m}_b^2(\mu )}{\mu ^2}}\right))`$ (2)
$`+\left({\displaystyle \frac{\alpha _s(\mu )}{\pi }}\right)^2\left(f_n^{20}(r)+f_n^{21}\mathrm{ln}\left({\displaystyle \frac{\overline{m}_b^2(\mu )}{\mu ^2}}\right)+f_n^{22}\mathrm{ln}^2\left({\displaystyle \frac{\overline{m}_b^2(\mu )}{\mu ^2}}\right)\right)`$
$`+{\displaystyle \frac{\frac{\alpha _s}{\pi }𝐆^2}{(\mathrm{\hspace{0.17em}4}\overline{m}_b(\mu ))^2}}[g_n^0+\left({\displaystyle \frac{\alpha _s(\mu )}{\pi }}\right)(g_n^{10}+g_n^{11}\mathrm{ln}\left({\displaystyle \frac{\overline{m}_b^2(\mu )}{\mu ^2}}\right))]\},`$
where $`rm_c/m_b`$. In terms of the more specific choice of $`\overline{m}_b(\overline{m}_b)`$, the moments have the simpler form

$`P_n`$ $`=`$ $`{\displaystyle \frac{1}{(\mathrm{\hspace{0.17em}4}\overline{m}_b(\overline{m}_b))^n}}\{f_n^0+\left({\displaystyle \frac{\alpha _s(\mu )}{\pi }}\right)f_n^{10}`$ (3)
$`+\left({\displaystyle \frac{\alpha _s(\mu )}{\pi }}\right)^2\left(f_n^{20}(r){\displaystyle \frac{1}{4}}\beta _0f_n^{10}\mathrm{ln}\left({\displaystyle \frac{\overline{m}_b^2(\overline{m}_b)}{\mu ^2}}\right)\right)`$
$`+{\displaystyle \frac{\frac{\alpha _s}{\pi }𝐆^2}{(\mathrm{\hspace{0.17em}4}\overline{m}_b(\overline{m}_b))^2}}[g_n^0+\left({\displaystyle \frac{\alpha _s(\mu )}{\pi }}\right)g_n^{10}]\},`$
where $`\beta _0=112/3n_f`$. The Born and order $`\alpha _s`$ terms of the moments are known since a long time  and the order $`\alpha _s^2`$ contributions for primary $`b\overline{b}`$ production for massless light quarks have been determined in Ref. . We have cross-checked these contributions with the explicit expressions for the corresponding contributions to $`R_{b\overline{b}}`$ given in Ref. . The order $`\alpha _s^2`$ contributions to the moments from secondary $`b\overline{b}`$ production, where the $`b\overline{b}`$ pair is produced through gluon radiation off light quarks, has been computed from the corresponding results for the $`R`$-ratio given in Refs. . These contributions only affect the coefficient $`f_n^{20}`$. The coefficients of the gluon condensate have been taken from Ref. . For convenience, the numerical results for the coefficients in Eqs. (2) and (3) for $`n=1,2,3,4`$ and $`m_c=0`$ ($`r=0`$) are collected in Tab. 1.

The numbers for the $`f_n`$’s agree with Ref.  up to a convention dependent factor of $`1/4`$, except for the results for $`f_n^{20}(0)`$, which are slightly larger accounting for the contributions from secondary $`b\overline{b}`$ production. The effects of the non-zero charm quark mass are generated either through virtual gluon self-energy effects or through real primary or secondary associated charm production. The corresponding contributions to $`R_{b\overline{b}}`$ for arbitrary mass constellations have been given in Refs. . Numerical values of the charm quark mass corrections to the coefficient $`f_n^{20}`$ are displayed in Tab. 2 for values of $`r`$ between $`0.1`$ and $`0.5`$.

We note that the numbers given in Tab. 2 also include the non-zero charm mass effects in the bottom quark pole-$`\overline{\text{MS}}`$ mass relation .

## Experimental Moments

For the contributions to the experimental moments from the $`\mathrm{{\rm Y}}(\text{1S})`$, $`\mathrm{{\rm Y}}(\text{2S})`$, $`\mathrm{{\rm Y}}(\text{3S})`$ and $`\mathrm{{\rm Y}}(\text{6S})`$ we use the averages for masses and $`e^+e^{}`$ widths given by the PDG . In Tab. 3 a collection of all contributions to the moments including uncertainties is given. The averages for the $`\mathrm{{\rm Y}}(\text{1S})`$, $`\mathrm{{\rm Y}}(\text{2S})`$ and $`\mathrm{{\rm Y}}(\text{3S})`$ are dominated by data from ARGUS  and the averages for the $`\mathrm{{\rm Y}}(\text{6S})`$ are from results from CUSB  and CLEO .

In the 4S–5S region between $`10.5`$ and $`10.95`$ GeV there have been measurements from CUSB  and CLEO . We find it remarkable that both experiments observed an additional resonance-like enhancement between the $`\mathrm{{\rm Y}}(\text{4S})`$ and $`\mathrm{{\rm Y}}(\text{5S})`$ at about $`10.7`$ GeV. Whereas CUSB fitted for $`\mathrm{{\rm Y}}(\text{4S})`$ and $`\mathrm{{\rm Y}}(\text{5S})`$ resonances within in a coupled channel model, the CLEO experiment was fitting an additional resonance, called “$`B^{}`$”, at $`m_B^{}=10.684\pm 0.013`$ GeV with an $`e^+e^{}`$ width of $`\mathrm{\Gamma }_B^{}^{e^+e^{}}=0.20\pm 0.11`$ keV. As a consequence the $`e^+e^{}`$ widths for $`\mathrm{{\rm Y}}(\text{4S})`$ and $`\mathrm{{\rm Y}}(\text{5S})`$ obtained from CLEO were systematically smaller than those from CUSB. In the PDG compilation, however, the existence of the “$`B^{}`$” contribution in the CLEO analysis was ignored when the averages for $`\mathrm{{\rm Y}}(\text{4S})`$ and $`\mathrm{{\rm Y}}(\text{5S})`$ have been determined.222 The PDG number for the electronic width of the $`\mathrm{{\rm Y}}(\text{4S})`$ is a weighted average of the results from CUSB , CLEO  and ARGUS , where the ARGUS result agrees better with the one from CUSB.  As a consequence, the contribution to the moments from the region between $`10.5`$ and $`10.95`$ GeV using the PDG averages is, although compatible within errors, systematically lower than the contributions one obtains using the numbers given in the original CLEO and CUSB publications (see Tab. 3). For our analysis we decided to ignore the PDG averages and to take the average of the original CLEO (including the “$`B^{}`$”) and the CUSB contributions to the moments for this energy region. We adopted the larger CLEO error assuming that the respective uncertainties are correlated. We believe that this conservative treatment is justified as long as the situation is not clarified. We note that the PDG treatment of the 4S–5S region does not affect the results for the hadronic vacuum polarization effects of $`\alpha (M_Z)`$ and $`(g2)_\mu `$ because the corresponding differences are much smaller than the total uncertainties. For the $`\mathrm{{\rm Y}}(\text{6S})`$ we used the PDG averages since it is unlikely that the different treatment of the enhancement observed between $`\mathrm{{\rm Y}}(\text{4S})`$ and $`\mathrm{{\rm Y}}(\text{5S})`$ has affected the fits above the $`\mathrm{{\rm Y}}(\text{5S})`$.

There is no direct experimental data for $`\sigma (e^+e^{}\gamma ^{}b\overline{b}+X)`$ in the region above $`11.1`$ GeV. However, there are measurements of the total hadronic cross section taken by a number of experiments up to energies close to $`M_Z`$ that are compatible with the Standard Model predictions. Furthermore, from measurements of $`R_b`$ at the Z pole by LEP and in the region between $`133`$ and $`207`$ GeV by LEP2333 At the Z pole $`R_b`$ is defined as the ratio of the total b quark partial width of the Z to its total hadronic partial width, $`\mathrm{\Gamma }_{b\overline{b}}/\mathrm{\Gamma }_{\mathrm{had}}`$, and for LEP2 energies it is defined as the ratio of the total $`b\overline{b}`$ cross section to the total hadronic one. , it is known that perturbative QCD agrees with the data for $`b\overline{b}`$ production to about $`1\%`$ at $`M_Z`$ and to about $`10\%`$ in the LEP2 region . It is therefore not unreasonable to estimate the experimental contribution to the moments from above the $`\mathrm{{\rm Y}}(\text{6S})`$ from perturbation theory itself. Although order $`\alpha _s^3`$ corrections to $`R_{b\overline{b}}`$ in the high energy expansion are known, we use the perturbative $`b\overline{b}`$ cross section to order $`\alpha _s^2`$ to estimate the continuum contributions because the theoretical moments are likewise only available to order $`\alpha _s^2`$. A much more subtle question is how to estimate the “experimental” uncertainties in this region. The approach of Ref.  assumes that the experimental data for the $`b\overline{b}`$ cross section including errors lie within the theoretical predictions and uses the small theoretical errors. This approach is quite similar to using finite energy sum rules  where an upper “duality” cutoff $`s_{\mathrm{max}}>11.1`$ GeV is used in the integral in Eq. (1).

In Tab. 3 we have displayed the continuum contributions to the moments to order $`\alpha _s^2`$. The charm quark mass has been set to zero. We have subdivided the continuum contribution into three parts coming from $`11.112.0`$ GeV (region 1), $`12\text{GeV}M_Z`$ (region 2) and $`M_Z\mathrm{}`$ (region 3) in order to visualize the impact of the various energy regions. The theoretical errors shown in Tab. 3 come from varying the strong coupling in the range $`\alpha _s(M_Z)=0.118\pm 0.003`$ and the $`\overline{\text{MS}}`$ bottom mass in the conservative PDG range $`\overline{m}_b(\overline{m}_b)=4.2\pm 0.2`$. The renormalization scale $`\mu `$ has been varied between $`2.5`$ and $`10`$ GeV. For the running of $`\alpha _s`$ and $`\overline{m}_b`$ four-loop renormalization group equations have been used. Our theoretical errors are larger than in Ref.  where also order $`\alpha _s^3`$ contributions have been included . Our central values have been obtained from the average of the respective extremal values. Region 1 has been displayed separately because data for the $`b\overline{b}`$ cross section could potentially be collected there by CLEO  in the near future.

Since the continuum region above $`11.1`$ GeV is unsuppressed and constitutes a sizeable contribution to the experimental moments444 We note that the relative contribution to the low-$`n`$ $`c\overline{c}`$ spectral moments coming from continuum energies above $`4.5`$ GeV is considerably smaller than for the continuum region above $`11.1`$ GeV in $`b\overline{b}`$ spectral sum rules, see e. g. Ref. . , using the small theory errors leads to a considerable model-dependence of the bottom quark mass. For our final error estimate of $`\overline{m}_b(\overline{m}_b)`$ we adopt a $`10\%`$ correlated (relative) error for the continuum regions 2 and 3 and ignore the respective theory errors given in Tab. 3. This choice is, in principle, as arbitrary as using the theoretical errors (or no errors at all), but should reduce the model-dependence to an acceptable level. For region 1 we use the theory error since here the variation of the $`\overline{\text{MS}}`$ mass in the conservative PDG bounds has the largest impact and leads to a variation of more than $`10\%`$.

## Uncertainties in $`\overline{m}_b(\overline{m}_b)`$

For the determination of $`\overline{m}_b(\overline{m}_b)`$ and the uncertainties we have used 4 methods:

1. The bottom mass $`\overline{m}_b(\overline{m}_b)`$ is determined from single moment fits ($`n=1,2,3`$) using Eq. (3).
2. The bottom mass $`\overline{m}_b(\mu )`$ is determined from single moment fits ($`n=1,2,3`$) using Eq. (2) and $`\overline{m}_b(\overline{m}_b)`$ is computed subsequently using renormalization group equations.
3. The bottom mass $`\overline{m}_b(\overline{m}_b)`$ is determined from fits to ratios $`P_n/P_{n+1}`$ ($`n=1,2`$) using Eq. (3).
4. The bottom mass $`\overline{m}_b(\mu )`$ is determined from fits to ratios $`P_n/P_{n+1}`$ ($`n=1,2`$) using Eq. (2) and $`\overline{m}_b(\overline{m}_b)`$ is computed subsequently using renormalization group equations.

For the analysis we employed only moments for $`n=1,2,3`$ to avoid the large higher order contributions $`(\alpha _s\sqrt{n})^k`$ that are characteristic for the large-$`n`$ moments and need to be summed. For method 3 and 4 we did not expand the perturbative series in the theoretical ratios $`P_n/P_{n+1}`$. We checked that expanding the theoretical ratios has only very small effects on the results. We employed four-loop renormalization group equations and used $`\alpha _s(M_Z)=0.118\pm 0.003`$, $`m_c=1.3\pm 0.2`$ GeV, $`\frac{\alpha _s}{\pi }𝐆^2=(0.024\pm 0.024)\text{GeV}^4`$ as theoretical input. The renormalization scale $`\mu `$ was varied between $`2.5`$ and $`10`$ GeV.

In Tab. 4 the results of our analysis for $`\overline{m}_b(\overline{m}_b)`$ are displayed in detail. The table shows the respective central values (in units MeV), which were obtained using the experimental and theoretical central values given before and $`\mu =5`$ GeV. The central values obtained with the four methods are within $`15`$ MeV around $`4.20`$ GeV. All uncertainties (in units of MeV) are presented separately. The experimental errors correspond to the uncertainties given in Tab. 3, where for the 4S-5S region our conservative CUSB-CLEO average has been used. The uncertainty from the 4S-5S region constitutes the largest experimental error from the resonance region. The errors from the continuum regions indicated with a subscript ”th” are obtained from the corresponding theory errors shown Tab. 3. For the continuum regions 2 and 3 also the errors coming from a $`10\%`$ deviation from the theory prediction are displayed (having no subscript). We note that the latter errors scale roughly linearly, i.e. assuming a $`5\%`$ ($`20\%`$) deviation the error decreases (increases) by a factor of two, etc.. This illustrates how strongly the bottom quark mass depends on assumptions for the experimentally unknown $`b\overline{b}`$ continuum cross section above the $`\mathrm{{\rm Y}}(\text{6S})`$. The theoretical errors have been obtained by varying each of the theoretical parameters in the ranges given above while the respective other parameters were fixed to their central values. For method 3 and 4, which are based on fitting ratios $`P_n/P_{n+1}`$, the same theoretical input parameters have been used for the moments in the numerator and those in the denominator. If $`\alpha _s(M_Z)`$ is chosen independently for the moments in the numerator and denominator, the errors for method 3 (4) are 38 (26) MeV for $`n=1`$ and 39 (13) MeV for $`n=2`$. If $`\mu `$ is chosen independently for the moments in the numerator and denominator, the errors for method 3 are 52 MeV for $`n=1`$ and 63 MeV for $`n=2`$. We found that, except for the variations of $`\mu `$, all resulting errors scale linearly with changes of the input parameters. The errors coming from variations of $`\mu `$ have been chosen to be the larger ones of the two deviations obtained in the ranges $`2.5\text{GeV}<\mu <5`$ GeV and $`5\text{GeV}<\mu <10`$ GeV. Note that the overall shift in the central value of $`\overline{m}_b(\overline{m}_b)`$ coming from the gluon condensate contribution is between $`0.1`$ and $`1`$ MeV. The shift caused by the non-zero charm quark mass is between $`1`$ and $`3`$ MeV, which is an order of magnitude smaller than for large-$`n`$ moments, where the charm mass effects are enhanced by a factor $`1/\alpha _s^2n`$ . Using the PDG average for the 4S-5S region instead of the CUSB-CLEO average (see Tab. 3), $`\overline{m}_b(\overline{m}_b)`$ is shifted by $`10`$ to $`15`$ MeV for method 1 and 2 and by $`2`$ to $`6`$ MeV for method 3 and 4. We consider the numbers given in Tab. 4 as the main result of this work.

In order to obtain combined errors from the uncertainties of the resonance data we treated one half of each error as correlated (being added linearly) and the other half uncorrelated (being added quadratically) because all data came from $`e^+e^{}`$ machines with common systematic uncertainties and, roughly, systematic and statistical uncertainties were found to have comparable sizes . In Tab. 4 the resulting combined resonance errors are shown in the line below the numbers for the $`\mathrm{{\rm Y}}(\text{6S})`$. The theoretical uncertainties including the errors adopted for the three continuum regions have been combined linearly since they do not have any statistical meaning and, in particular, the division into three continuum regions is completely arbitrary. Note that for regions 2 and 3 the uncertainties from the $`10\%`$ variation of the theory prediction have been adopted and not the smaller theoretical errors. The resulting combined error is displayed in the line below the numbers for $`\delta \mu `$. To obtain our total error (last line in Tab. 4) we added the resonance, the continuum and the theory errors linearly. As expected we find that the uncertainties from the continuum have the largest impact for the moments with $`n=1`$ and $`2`$ and that they are partially canceled when ratios are used for the fitting. However, the theoretical errors are larger for fits with ratios than with single moments, if the theoretical parameters are chosen independently for numerator and denominator. In general, the total error decreases for larger $`n`$. This trend does, however, not continue for higher values $`n>3`$ particularly for methods 2 and 4. Compared to the results of Ref.  our total errors are much larger, particularly for fits involving $`P_1`$ and $`P_2`$. This is mainly because we adopted more conservative errors for the continuum region in the experimental moments, and we combined uncertainties linearly, when they cannot be treated statistically. Since we believe that $`P_3`$ can be computed reliably using Eqs. (2) and (3), we adopt the average of the total errors in the third column of Tab. 4 as our final estimate for the uncertainty and obtain (rounded to units of $`10`$ MeV)

$$\overline{m}_b(\overline{m}_b)=\mathrm{\hspace{0.17em}4.20}\pm \mathrm{\hspace{0.17em}0.09}\text{GeV}.$$
(4)
We do not take into account the small $`50`$ MeV error from $`P_2/P_3`$ for method 3 because the small error from variations of $`\mu `$ only persists if the same $`\mu `$ is chosen for both moments. For an independent choice the error is considerably larger (see text above). If a $`20\%`$ relative uncertainty is assumed for the continuum regions 1, 2 and 3, the final error increases by $`20`$ MeV. Our result in Eq. (4) is compatible with the result from Ref. . Our result is also compatible with results from fitting large-$`n`$ moments at NNLO in the non-relativistic expansion  and with recent results, using different methods, obtained in Refs. .

Based on the numbers given in Tab. 3 and 4 it is straightforward to discuss the impact of improved measurements of the resonance parameters at CLEO. Assuming improved measurements for the electronic widths of $`\mathrm{{\rm Y}}(\text{1S})`$, $`\mathrm{{\rm Y}}(\text{2S})`$ and $`\mathrm{{\rm Y}}(\text{3S})`$  at the level of $`2\%`$, the combined resonance errors shown in Tab. 3 would be reduced by a factor $`2/3`$, which would reduce the total error (obtained from the third column) by about $`10`$ MeV. An improved measurement of the 4S-5S region and the $`\mathrm{{\rm Y}}(\text{6S})`$ with the same precision would result approximately in a reduction of the total error by an additional $`10`$ MeV. A further reduction of the error below about $`70`$ MeV, however, will be difficult to achieve without real experimental data for $`\sigma (e^+e^{}b\overline{b})`$ in the continuum region above the $`\mathrm{{\rm Y}}(\text{6S})`$ with a precision of better than ten percent.

## Conclusion

We have given a detailed analysis of the uncertainties in the $`\overline{\text{MS}}`$ bottom quark mass $`\overline{m}_b(\overline{m}_b)`$ using low-$`n`$ spectral moments of the cross section $`\sigma (e^+e^{}\gamma ^{}b\overline{b})`$. For the experimental moments we employed experimental data for the $`\mathrm{{\rm Y}}(\text{1S})`$ \- $`\mathrm{{\rm Y}}(\text{6S})`$ resonances and the order $`\alpha _s^2`$ QCD predictions for the continuum region above $`11.1`$ GeV. For the 4S–5S region between $`10.5`$ and $`10.95`$ GeV we found that the PDG averages for the electronic partial widths of $`\mathrm{{\rm Y}}(\text{4S})`$ and $`\mathrm{{\rm Y}}(\text{5S})`$ based on data from CUSB, CLEO and ARGUS contain an inconsistency, stemming from the fact CLEO also assumed the existence of an additional resonance at about $`10.7`$ GeV which was ignored in the averaging procedure. For our analysis we therefore used the original CUSB and CLEO results. For the continuum region above the $`\mathrm{{\rm Y}}(\text{6S})`$ we assumed a $`10\%`$ error in our final error estimate. For the theoretical moments we used perturbative results at order $`\alpha _s^2`$ including also the contributions from secondary $`b\overline{b}`$ production and finite charm mass effects. As our final result we get $`\overline{m}_b(\overline{m}_b)=4.20\pm \mathrm{\hspace{0.17em}0.09}`$ GeV and we conclude that more precise data for the electronic partial widths of the $`\mathrm{{\rm Y}}`$ resonances at CLEO could reduce the error by about $`20`$ MeV.

## Acknowledgements

We thank Ed Thorndike for providing details of the analysis in Ref. .


format: text

text:
# A simplified numerical model of coronal energy dissipation based on reduced MHD

## 1 Introduction

One of the main unresolved problems in solar physics concerns the mechanism by which the solar atmosphere is heated from several thousands degrees in the photosphere to millions in the corona. It is now commonly accepted that the ultimate source of energy lies in the convective motions in and below the photosphere, and that a reliable model of coronal heating has to deal with the transfer, the storage and finally the release of this energy into the solar corona. Several conceptual models have been proposed, such as Alfvén waves, electric currents and MHD turbulence – see for a review Zirker (1993) – where in all cases the magnetic field plays a key role in the dynamics.

On the other hand, it is also well established that impulsive events (e.g. solar flares, X-ray bright points) are distributed in the corona over a large range of scales in size, energy and duration (Dennis, 1985; Crosby, Aschwanden & Dennis, 1993; Pearce, Rowe & Yeung, 1993; Krucker & Benz, 1998; Aletti et al., 2000; Aschwanden et al., 2000) and that large events seem to be made up of the superposition of a myriad of smaller unresolved events. It was Parker (1988) who suggested that the active X-ray and UV corona is composed of a swarm of localized impulsive bursts of energy called nanoflares (or picoflares). Much theoretical work has been done to investigate Parker’s conjecture and more generally the statistical nature of the solar coronal heating. They mainly follow two complementary schools which refer on the one hand to the dynamics of complex systems and on the other hand to fluid mechanics.

The statistical properties of flaring activity allow one to view the solar corona as a complex system which can be described with cellular automata (CA) models (Lu & Hamilton, 1991; Lu et al., 1993; Vlahos et al., 1995; Galsgaard, 1996; Georgoulis & Vlahos, 1998; Vassiliadis et al., 1998; Isliker, Anastasiadis & Vlahos, 2000, 2001; Charbonneau et al., 2001; Krasnoselskikh et al., 2001). CA models have become increasingly useful in the study of complex systems because they permit the study of an entire system without ignoring the effects of individual components of the system. There are many natural applications such as substorms in the magnetotail (Takalo et al., 1999), star formation in spiral galaxies (Lejeune & Perdang, 1996) or earthquakes (Carlson & Langer, 1989). A cellular automaton is based upon the idea of the locality of influence : a system is distributed in space, and nearby regions have more influence than those far apart (see for instance MacKinnon & Macpherson (1997) for a study of a nonlocal communication). A grid of cells is used to represent the components of a system, and each cell is given a set of phenomenological rules concerning its surrounding neighbors. The system evolves over several iterations by allowing each cell to interact using the given rules. What makes CA so interesting and useful is that after many iterations they reveal complex structures and arrangements that form across great distances even though each cell only takes into account local information. Self-Organized Criticality (SOC) (Bak et al., 1987; Kadanoff et al., 1989; Hwa & Kardar, 1992; Sornette, 2000) refers to the spontaneous organization of such an externally driven system into a globally stationary state over many scales.

On the other hand, we find fluid models which give the physical description that is missing in CA. Much work has been done on statistical solar flares (Longcope & Sudan, 1994; Walsh, Bell & Hood, 1995; Einaudi et al., 1996; Galsgaard & Nordlund, 1996; Hendrix & Van Hoven, 1996; Dmitruk et al., 1998; Galtier & Pouquet, 1998; Georgoulis et al., 1998; Galtier, 1999; Walsh & Galtier, 2000) but most of them suffer from the fact that statistical simulations of flares studied in the context of forced resistive MHD equations are possible only at the cost of huge computational expenses. Nevertheless it has been possible to show important properties, e.g. that the dissipative events produced exhibit power-law distributions (for total energy, peak of luminosity and duration) in agreement with X-ray observations, but with generally a much smaller “inertial range” than the CA counterpart.

A recent debate about the possible existence of sympathetic flaring, i.e. the correlation in time of two successive events (Pearce, Rowe & Yeung, 1993; Wheatland, Sturrock & McTiernan, 1998; Boffetta et al., 1999; Wheatland, 2000; Lepreti, Carbone & Veltri, 2001; Galtier, 2001), suggests the possibility to dismiss CA as a model of solar flares since standard CA models do not produce correlated events (non-Poissonian statistics such as power-law waiting time distributions). But in fact many CA models exist in the literature like nonconservative models (Christensen & Olami, 1992) which turn out to be able to generate the statistics expected for sympathetic flaring. But the question of the existence of sympathetic flaring in the corona has not yet found an answer since in particular there is still a debate about what we mean by event.

The problem of coronal heating is intimately linked to the existence of nanoflares whose Probability Distribution Function (PDF) in energy is supposed to be a power-law steeper than that for regular flares. Let us assume that the PDF in energy $`E`$ of events is distributed according to a power-law of index $`\zeta `$, i.e. $`\mathrm{Pr}(E)E^\zeta `$. It is then possible to show that there exists a critical slope of index $`\zeta _c=2`$ (Hudson, 1991). Indeed, the total energy released in the corona by events between $`E_{min}`$ and $`E_{max}`$ is $`(E_{max}^{2\zeta }E_{min}^{2\zeta })/(2\zeta )`$ which means that if $`\zeta <2`$ the main contribution comes from high energy events, whereas if $`\zeta >2`$ it comes from smaller events (the swarm of nanoflares). The average power dissipated in a large flare is of the same order of magnitude as the total average power emitted by the corona, $`10^3\mathrm{W}\mathrm{m}^2`$, which proves that regular flares can not account for coronal heating since they are episodic events seen over and above the average coronal background. It is then natural to think that a swarm of very small and still unobservable events may dominate the heating process. One of the main challenges of statistical flare models is to know whether or not it is possible to produce power-law distributions for any relevant quantity such as energy, luminosity or duration, and what the power-law indices are.

The aim of this paper is to introduce a hybrid model for a solar magnetic loop which is somewhat intermediate between CA models and full MHD or reduced MHD models. In this model, we will inject and store energy into a coronal loop (our numerical domain) via wave propagation from the photosphere (our numerical boundary). The trigger for an event is determined in a way analogous to conventional CA models, i.e. with a threshold in the current. However, during the subsequent event the current is dissipated and the magnetic field recomputed using Maxwell’s equations. Let us note that this is a minimal consistency requirement for the field evolution which is not always incorporated in CA models. In practice, the model allows current concentrations to form kinematically (advection from the photosphere), but not dynamically (the nonlinear part of the Lorenz force, $`j\times b`$). The model is non-trivial because of the threshold dynamics of the dissipation, which mimics the nonlinear terms, but the model is much simpler to integrate than the full MHD equations, therefore allowing a fast computation of statistics (events sizes and durations), and a comparison both with observations and full numerical simulations.

The paper is organized as follows. In Section 2 we give a detailed description of the CA model and show its basic behavior through some numerical experiments. In Section 3 the results of a parametric study are given and discussed. In Section 4 we summarize the properties of the model, we present a comparison with observations, and we draw some conclusions.

## 2 The model

In the original 3D lattice model developed by Lu & Hamilton (1991) the physical quantity defined on each lattice is the magnetic field. The system is driven from the outside by adding randomly in space a random magnetic field. The process continues until a reconnection instability criterion is satisfied at any point of the 3D lattice, i.e. until the magnetic gradient exceeds a critical value at this point. Then the magnetic field is redistributed (diffused) towards neighboring nodes with the possibility to transfer the instability as well. The redistribution process stops when the system is completely relaxed. Then another random amount of magnetic field is added to the system. An event called avalanche is associated to the rapid diffusion of the magnetic field.

Subsequent models use the magnetic vector potential $`A`$ rather than the magnetic field since the divergence-free condition for the magnetic field is then automatically satisfied. For example in the recent model developed by Isliker, Anastasiadis & Vlahos (2000, 2001) where the 3D lattice represents an ensemble of magnetic loops, the knowledge of $`A`$ allows to reconstruct the magnetic field topology and eventually the structure of the current density. To do so they introduce the notion of derivative. The present simplified model belongs to this class of models but only one typical coronal loop will be considered and simulated. The detailed description of the model is now given.

### 2.1 Basic idea : on-off mechanism and turbulence

A possible reason for the behavior of the corona is that it lies in a turbulent state. A model of coronal loops should therefore allow for the effects of turbulent fluctuations. This is possible with CA models at a very superficial level through an on-off mechanism. The idea behind the threshold dynamics of our model, the on-off mechanism, is the following. The forcing due to the convective granules, although applied on a range of scales, has a typical length scale that is supposed to be far greater than the dissipative scale. The connection between the forcing and the dissipative length scales is made through a turbulent mechanism. During the “off” phase, i.e. the loading phase, the plasma is in a laminar state where the dynamics is essentially governed by the linear terms (and the loading). Because the system is driven slowly, parts of the system or even the entire system can stay in principle in this state for very long periods of time. When sufficient energy is accumulated in the loop some nonlinear instability appears which triggers the rapid generation of small scales. The inertial range of the turbulent energy spectrum extends to small scales and makes the link between the typical forcing length scale and the dissipative length scale. This “on” phase is therefore characterized by a sudden increase of the dissipative terms in the RMHD equations (see section 2.2) leading to a bursty event. Then the system returns immediately to an “off” phase. The nonlinearities of the RMHD equations are therefore included in the model in a very schematic way through a threshold dynamics only.

### 2.2 Description of the model

#### Geometry of the model.

The model describes a coronal loop anchored in the photosphere whose footpoints are randomly moved. The presence of a strong axial magnetic field leads to essentially 2D dynamics, i.e. perpendicular to the mean magnetic field, for which the approximation of the RMHD equations (Strauss, 1976) is a good model. As we can see on Fig. 1, the 3D regular grid is made up of a set of planes distributed along the loop between which the information propagates through Alfvén waves. Therefore each plane will evolve essentially independently from each other. Both boundary planes represent the photospheric footpoints, while the intermediate planes represent the loop itself, as if it were unbent. The curvature of the loop is not taken into account since the width of observed coronal loops (see e.g. observations from the TRACE instrument) appears to be constant along the loop and is much smaller than their length.

#### RMHD equations.

Our aim is to compute the temporal evolution of the velocity field $`v`$ and the magnetic field $`B`$ (or $`bB/\sqrt{\rho _0\mu _0}`$ if we consider only fields with the same physical dimension). We assume the presence of a strong and uniform axial magnetic field along the $`z`$-axis $`b_0=b_0e_z`$ ($`B_0/\sqrt{\rho _0\mu _0}`$) and we consider small perturbations to this field. We separate $`bb_0`$ and $`v`$ into parallel components ($`b_ze_z`$ and $`v_ze_z`$) and orthogonal components ($`b_{}`$ and $`v_{}`$). With the following additional hypotheses, (i) the scales along the $`z`$-axis are larger than the scales in the orthogonal directions (gradients along the $`z`$ axis are negligible), and (ii) the kinetic pressure is negligible compared to the magnetic pressure (i.e. the plasma is cold, $`\beta 1`$), we then obtain the RMHD equations (Strauss, 1976)

$`_tv_{}+(v_{}_{})v_{}`$ $`=`$ $`b_0_zb_{}+\nu \mathrm{\Delta }_{}v_{}`$
$`+(b_{}_{})b_{}(b_{}^2/2),`$
$`_tb_{}+(v_{}_{})b_{}`$ $`=`$ $`b_0_zv_{}+\eta \mathrm{\Delta }_{}b_{}`$
$`+(b_{}_{})v_{},`$
where $`\nu `$ and $`\eta `$ are respectively the kinematic viscosity and the magnetic resistivity. To each grid point are associated two scalar fields $`a^s`$, with $`s=\pm `$, from which the Elsässer fields are derived

$$z^s=v_{}+sb_{}=_{}\times a^se_z.$$
(3)
All other fields (magnetic and velocity fields, current density, vorticity…) are derived from $`a^s`$ by analogy with the standard magnetic and kinematic equations.

#### Initial state and boundaries.

The fields $`a^s`$ are taken to be zero initially. Each cross-sectional plane along the loop is periodic for the $`a^s`$ variables. All numerical computations for each plane are made in the Fourier space. We use inverse Fast Fourier Transforms (FFT) when we temporarily need to know the values of a field at some positions in real space.

#### Alfvén wave propagation.

In the right hand side of equations (2.2-2.2), the first terms correspond to the Alfvén waves propagation along the $`z`$-axis : $`a^+`$ propagates to the bottom of the simulation box (low values of $`z`$), while $`a^{}`$ propagates to the top. These propagations are modeled by specific cellular automaton rules used at each time step $`\delta t`$ of the simulation, corresponding to discretization of the Alfvén waves terms of equations (2.2-2.2) in the form :

$$a^s(z,t+\delta t)=a^s(z+s\delta z,t)$$
(4)
There is no loss of energy (or reflection) during the Alfvén wave propagation since the density is assumed to be constant ; however we assume that there is a total reflection of the waves when they reach the photosphere, i.e. the two boundary planes of the simulation box.

#### Loading.

The action of the photospheric granules on the magnetic footpoints is modeled by a random increment $`\delta \mathrm{\Psi }`$ to the fields on both boundary planes. This increment has random Fourier coefficients but has globally a power-law energy spectrum in $`k^\alpha `$ (the total intensity $`P_{\mathrm{load}}=(k\delta \mathrm{\Psi })^2dk`$ is also a parameter of the simulation). Indeed, observational evidence suggests that the convective layer is in a turbulent state : the photospheric granules exhibit a turbulent power spectrum of velocity consistent with a Kolmogorov energy spectrum in $`k^{5/3}`$ but only for a narrow inertial range of wavenumbers ($`\mathrm{}1/k<3\mathrm{arcsec}`$) (Roudier & Muller, 1986; Chou et al., 1991; Espagnet et al., 1993). We emphasize that there is no loading in the other parts of the loop : energy is solely carried by Alfvén waves. Furthermore, these waves reflect on both photospheric boundary planes.

#### Dissipation criterion.

We assume that dissipation occurs when an instability criterion is satisfied, which is the condition that the magnitude of the current density $`J`$ exceeds a critical value $`J_c`$. As $`J_z=Je_z`$ can be derived from the computed variables $`a^s`$ ($`J_x`$ and $`J_y`$ are negligible), this dissipation criterion is likely to have more physical meaning than the criteria used for example in classical sandpile and in more elaborate models (see Charbonneau et al. (2001) for a review). However, there is still some doubt about the quality of this criterion, as will be discussed later (see section 4.1).

#### Reconfiguration of the field.

When the current density exceeds a given threshold at some real-space grid points in a given plane, the nonlinear terms of equations (2.2-2.2) which are negligible during the loading phase (the off phase) become large and dominate the dynamics of the fields. They are quickly balanced by the dissipative terms when the energy cascade reaches the dissipative scale. This “on” phase (see section 2.1) is modeled by a diffusion-like process for the magnetic and velocity fields which tends to reduce the magnitude of the current density and the vorticity.

The detailed algorithm is an updated version of the one introduced by Einaudi & Velli (1999). At a time $`t`$ for plane $`z_0`$ we compute the current density $`J_z(x,y,z_0)=\mathrm{\Delta }_{}A_z`$, where $`A_z=(a^+a^{})/2`$ is the magnetic vector potential and $`\mathrm{\Delta }_{}`$ denotes the Laplacian operator in the plane $`z_0`$. If at some grid point $`(x,y,z_0)`$ the value of $`|J_z|`$ exceeds the threshold $`J_c`$, $`A_z`$ is updated in the time $`\delta t_c`$ (with $`\delta t_c\delta t`$ ) according to the equation $`A_z(x,y,z_0;t+\delta t_c)=A_z(x,y,z_0;t)\eta \delta t_cJ_z(x,y,z_0;t)`$, which corresponds to current dissipation. The current density $`J_z`$ corresponding to $`A_z(t+\delta t_c)`$ is then computed, and this dissipation process is iterated until $`J_z`$ does not exceed the threshold anywhere in the plane $`z_0`$. However, note that after the first iteration of this process, we take $`CJ_c`$ as a threshold instead of $`J_c`$. The “dissipation efficiency” $`C`$ is a number between $`0`$ and $`1`$ which guarantees that the system is in a relaxed state after the whole dissipation process.

#### Energy release

During this relaxation process, magnetic and kinetic energies are released. The energy release in each plane can easily be computed from the variations of $`b_{}^2`$ and $`v_{}^2`$ in the plane. It is the primary variable for our statistics. Note that topological modifications of the magnetic field may be expected : the connectivity of the magnetic field lines is modified because of the field diffusion. One of the possible interpretations of this phenomenon is magnetic reconnection (see however the discussion in section 4.1).

### 2.3 Time and space scales.

Let $`\delta x`$ and $`\delta z`$ be the distance between grid points in the $`x`$ (or $`y`$) and $`z`$ directions respectively, and let $`\delta t`$ be the time step. If we assume that the loop length $`L`$ is $`1`$ to $`100\mathrm{Mm}`$, then $`\delta z`$ is $`30\mathrm{km}`$ to $`3\mathrm{Mm}`$ for a typical resolution of $`N_L=30`$ points (i.e. 30 planes along $`z`$). The analog assumptions for a loop width $`\mathrm{}`$ ($`=L/10`$) of $`0.1`$ to $`10\mathrm{Mm}`$ give $`\delta x`$ between $`1.5`$ and $`150\mathrm{km}`$ for a typical resolution of $`N_{\mathrm{}}=64`$.

We can also determine time scales for the model : as the Alfvén speed is one in the model units, i.e. in units of $`\delta z/\delta t`$, we have $`\delta t=\delta z/v_A`$ : the time step is the time needed by the Alfvén wave to propagate from one plane to its neighbors. For $`B_0=10^3`$ to $`10^2\mathrm{T}`$ (i.e. $`v_A1`$ to $`10\mathrm{Mm}\mathrm{s}^1`$ with density $`\rho _010^{12}\mathrm{kg}\mathrm{m}^3`$), this yields $`\delta t`$ between $`3\times 10^3\mathrm{s}`$ and $`3\mathrm{s}`$. Another time scale in the system is the coherence time of photospheric loading $`\delta t_l`$. It is modeled by a periodic re-initialization of the coefficients of the loading increment $`d\mathrm{\Psi }`$, which occurs every $`200`$ time steps $`\delta t`$, or $`0.6`$ to $`600\mathrm{s}`$. This is small compared to observational evaluations of the photospheric coherence time, but the relevant point is the good separation between time scales ; besides, larger values of the photospheric loading coherence time do not alter the results of the model. When no other indication is given, the time step $`\delta t`$ is the unit of time ; for example, on Fig. 2, the $`x`$-axis range maximum is $`\mathrm{10\hspace{0.17em}000}\delta t`$, i.e. between approximately $`30\mathrm{s}`$ and $`8\mathrm{h}`$. At last, the shortest time scale in the model is the cascade time scale $`\delta t_c`$, which is the time step for dissipations within a cascade, and which is analogous to the non-linear time scale of MHD models. It is assumed to be completely separated from the other time scales, i.e. $`\delta t_c\delta t\delta t_l`$.

A direct consequence of the separation between the cascade time scale $`\delta t_c`$ and the time $`\delta t`$ of wave propagation between planes is that the fields of neighboring planes are expected to be uncorrelated.

## 3 Numerical results of the model

### 3.1 Model behavior

The simulations presented in this paper have been performed on a local quadri-RS/6000 IBM workstation at IAS. A typical run of $`\mathrm{200\hspace{0.17em}000}`$ time steps with a resolution of $`N_L\times N_{\mathrm{}}^2=30\times 64^2`$ takes between 2 days and 2 weeks for one CPU, depending on the parameters.

#### Initial growths of energy and dissipation.

As the initial fields in the simulation box are zero, the initial kinetic and magnetic energies are zero. The loading phase inputs energy into the system at each time step $`\delta t`$ which gives a growth of the total energy of the system as shown on Fig. 2. Then the current density threshold can be reached at some points and dissipation occurs, which slows down the initial energy growth. At the same time, the average rate of dissipation increases until a stationary state is reached.

#### Histograms and fitting methodology.

The heights of the bars of the histograms we plot are normalized by their width and they are are divided by the number of events, i.e. our histograms are empirical PDFs. A least-squares linear fit is then performed in bi-logarithmic axes, on a range determined by visual inspection of the histogram (see Figure 3b). This gives error bars on the slope of the linear fit, which is the slope of the expected histogram power-law. However, we should keep in mind that the choice of the fitting range often introduces much larger error bars (typically $`\pm 0.1`$ to $`0.2`$) than the error bars of the least-squares linear fit of the slope (typically $`\pm 0.01`$ to $`0.05`$). The error bars we give from now are conservative estimates taking into account the fitting range uncertainty.

#### Choice of the variable used for the statistics and general shape of the PDF.

Former studies (Aletti et al. (2000); Aletti (2001)) plotted the histograms or PDFs of the magnetic energy dissipation $`\mathrm{\Delta }E_{\text{tot}}`$ calculated in the whole simulation box (Fig. 3). The global shape of the PDF was a Gaussian. A power-law $`\mathrm{Pr}(\mathrm{\Delta }E_{\text{tot}})\mathrm{\Delta }E_{\text{tot}}^6`$ seemed to appear as a deviation from the Gaussianity in the tail of the distribution, but it only spanned half a decade, which makes it perhaps not so relevant. On the contrary, we choose to plot the PDF of the magnetic energy $`\mathrm{\Delta }E_i`$ dissipated in one given plane $`i[1,N_L]`$ (Fig. 3b). As the computations are done in the Fourier space, this is our primary variable. The power-law that can be fitted to the PDF of this variable has a much wider range (more than $`2`$ decades) and is much less steep (the index is between $`1`$ and $`2`$) than in the former case.

This can be explained by a Central Limit Theorem after remarking that $`\mathrm{\Delta }E_{\text{tot}}=_i\mathrm{\Delta }E_i`$ and that the $`E_i`$’s are quasi-independent (as expected, the correlation between fields in neighboring planes is very low), thus the PDF of $`\mathrm{\Delta }E_{\text{tot}}`$ is the convolution of the PDFs of all $`\mathrm{\Delta }E_i`$ for $`i[1,N_L]`$. The difference between the PDFs in both cases stresses the importance of the choice of the variable used for the statistics. It also emphasizes that in the case of the statistics of observational data, we have to be careful about the definition of an “event”.

Both distributions of $`\mathrm{\Delta }E_i`$ and of $`\mathrm{\Delta }E`$ show a maximum. In the case of $`\mathrm{\Delta }E_i`$, it is a consequence of the finite range of the power-law distribution ; the position of the maximum depends the average event size and on the slope. In the case of $`\mathrm{\Delta }E`$, knowing the distributions $`\mathrm{\Delta }E_i`$ in all planes $`i`$, it can be seen as a simple consequence of the Central Limit Theorem.

#### Effect of initial growth on statistics.

During the initial energy growth, the PDF of the energy of events is different than during the stationary state. In particular, it is shifted to the left, i.e. the events are smaller. As a result, the left part of the PDF of events energy gets higher than what it would be if stationary state events only were taken into account, as seen on Fig. 4. As we are interested in stationary state events, we choose to exclude events occuring during the initial energy growth from the statistics.

#### Typical fields.

As the model is built on phenomenological evolution rules, it is not expected that the fields produced by the simulation coincide with any real picture. However, as we have tried to be as close as possible to the original MHD equations it is interesting to see how far the fields are realistic and what the limits of this phenomenological model could be. Typical magnetic and current density fields for $`\alpha =2`$ and $`\alpha =4`$ are shown on Fig. 5. On both samples but especially for high values of $`\alpha `$, we can notice that high current densities occur in magnetic “islands” and in regions where magnetic field densities are high. We do not observe many structures such as current sheets or possible reconnection sites, as will be discussed in section 4.1, although they are more present for small values of $`\alpha `$. Large-scale photospheric forcing (large $`\alpha `$) leads to large-scale structures.

### 3.2 Parametric study of event energy PDFs

A parametric study is performed in order to explore the influence of the simulation parameters on the magnetic energy dissipations PDFs. A reference set of parameters, called $`\text{a}`$, is chosen (see Table 1), and it gives the PDF shown on Fig. 3 . The PDFs obtained for other sets of parameters will be compared to the PDF obtained for $`\text{a}`$. Most of the sets of parameters correspond to the modification with respect to $`\text{a}`$ of one parameter (dissipation efficiency $`C`$, magnetic resistivity $`\eta `$, loading spectrum index), which is in italic in Table 1. All simulations were performed on $`\mathrm{200\hspace{0.17em}000}`$ timesteps, which seems sufficient to achieve a long stationary state after the initial energy growth phase, and to achieve good statistics for the PDFs. Histograms were done with data from the $`\mathrm{100\hspace{0.17em}000}`$ last timesteps.

Other parameters, which are not changed during the study, include the grid size (see above) and the current density threshold $`J_c=300`$. A higher grid resolution would have been interesting so as to get a broader wavelength range, but it would need a rescaling of other parameters and longer computation times. The current density threshold fixes the scale of current density, so its value has no intrinsic meaning.

#### Dissipation efficiency.

Dissipation efficiency tells how much the system gets relaxed after a series of iterative dissipations : the current density threshold $`J_c`$ is replaced by a new threshold $`CJ_c`$ after the first dissipation. With respect to the value $`C=0.5`$ used in parameters set $`\text{a}`$, the dissipation efficiency can be set to almost any value of its range $`[0,1]`$ of valid values with almost no visible change in the PDFs.

#### Magnetic resistivity.

Magnetic resistivity $`\eta `$ could vary in the range $`[310^4,310^3]`$. A numerical stability analysis, given the time step fixed to one and the wavenumber range, shows that higher values of $`\eta `$ would result in numerical instability. On the other hand, lower values of $`\eta `$ would lead to longer computational time. However, $`\eta `$ has mainly an influence only on the dissipation process length; a change in the value of $`\eta `$ has little influence on the histograms of dissipated energies.

#### Loading spectrum.

The reference parameters set $`\text{a}`$ has a loading spectrum index $`\alpha =5/3`$, corresponding to the spectrum in the inertial range of Kolmogorov turbulence. In parameters sets $`\text{a}`$ and $`\text{h}`$ to $`\text{o}`$, $`\alpha `$ varies from $`3/2`$ to $`4`$ by a maximal interval of $`1/3`$. Another series of simulations was performed with lower loading power values to explore the influence of the ratio $`P_{\mathrm{load}}/J_c`$. For both series, for high values of $`\alpha `$, the power-law is well defined (3 to 4 orders of magnitude wide). Its slope index is approximately $`1.6`$, and this value depends neither on the loading spectrum index $`\alpha `$ (as seen in Table 2 and Figure 6, and as will be discussed in section 4.1) nor on the loading intensity $`P_{\mathrm{load}}`$. For low values of $`\alpha `$, however, power-laws were difficult to obtain, and their slopes were sensitive to both loading spectrum index and loading intensity.

### 3.3 Statistics of durations of events

Durations of events extend over two decades. They are indeed a discrete variable, which is a multiple of the cascade time step $`\delta t_c`$, and their maximum value is a few hundreds times $`\delta t_c`$. Histograms can be obtained, although their width is too narrow to perform relevant power-law fitting (see figure 7).

The duration of an event is correlated with its energy, like $`dE_idt_i^{1.76}`$, as seen on the scatter-plot on figure 8. Another way to visualize this correlation is to select events according to their duration, and to plot the histogram of energies of events from this population, as shown on figure 9. One possible observational consequence could be that missing long-duration events, due for example to short observation times, can lead to energy histograms with narrower ranges and steeper slopes.

## 4 Discussion and Conclusions

### 4.1 Properties of the model

The CA model presented here differs from previously presented work in two features. On the one hand, the energy pumping due to photospheric motion is known quite accurately via solution of the Alfvén wave propagation equation with reflecting boundary conditions. Each cell therefore receives and sends energy to neighboring cells along the loop axis via a wave equation. Energy redistribution to cells on the same loop plane however occurs using an instability and redistribution criterion. This criterion is a threshold in current, and satisfies the basic requirements for magnetic fields of divergence-free conditions and realistic current redistribution.

#### Power-law slope of energy distributions as a function of parameters.

When parameters sets lead to wide and robust power-law statistics of events energies, it seems that the slope of these power-laws is quite the same and takes a value close to $`\zeta =1.6`$. If this behavior is confirmed by further parametric studies, it could be interpreted as the “universal” behavior of SOC systems, described for example in Bak et al. (1987).

It is however interesting to note that Georgoulis & Vlahos (1998) observe in their model a variability of the slope of event energy (and also peak luminosity and duration) as a function of the loading : the slope is steeper when the amplitude distribution of the loading increments has a steeper power-law. However, this discrepancy is due to the differences in the way the system is driven : their loading increments are discrete in space whereas our loading has a varying spatial spectrum and no power-law amplitude distribution variation.

#### Reconnection, and the quality of the dissipation criterion.

Snapshots of the superposition of the magnetic lines with the current densities (see figure 5) reveal that the sites of dissipation often do not correspond to reconnection of magnetic field lines ; the reason is that one may have intense currents (and our dissipation criterion is based precisely on current intensity) without having a topology of the field where reconnection occurs (instabilities such as resistive kinks and/or tearing modes are triggered by a combination of currents and current gradients). Though the typical field structures observed resemble fields from turbulence simulations, one sees that our simplified model distributes dissipation in a different, more homogeneous way. Of course, a cellular automaton model, which models the non-linear terms of the equations through simple threshold dynamics, is not supposed to generate such structures, but we can try to understand what can be done to improve the model. The physical quantities available in our model make it for example possible to use a more elaborate dissipation criterion which would model more accurately the reconnection instability threshold, such as for example introducing a combined criterion on current and current gradients as a trigger for relaxation.

### 4.2 Comparison with observations

Impulsive coronal events are statistically distributed over an energy range of some eight orders of magnitude. Since Parker’s idea of the existence of nanoflares, it has been thought that at some point in the quest of small scale coronal structures we will break the self-similarity of the solar corona by the observation of a steepening of the flare distribution with finally a power-law index $`\zeta `$ greater than the critical value of $`2`$. Recent data analyses (Krucker & Benz, 1998; Aschwanden et al., 2000; Parnel & Jupp, 2000) seem to show this behavior with observed values of $`\zeta `$ going up to $`2.6`$. Empirical formulas have been used here to determine flares energies from observed luminosities but new analyses of the data (Aschwanden & Charbonneau, 2002; Benz & Krucker, 2002) reveal in fact the existence of a bias due to the finite range of temperature on which the observations are made. The correction of this temperature bias leads eventually to a value of $`\zeta `$ close to $`1.6`$ valid for the whole range of energies, from unresolved, X-ray observations with the Solar Maximum Mission (SMM) to Extreme Ultra-Violet observations with the Extreme-Ultraviolet Imaging Telescope on the Solar and Heliospheric Observatory (SoHO/EIT) and the Transition Region And Coronal Explorer (TRACE) – see Aschwanden et al. (2000). Such an observational bias, as well as the one described in the end of paragraph 3.3, show the importance of defining well what an event is and what its characteristics are, when observing the corona but also when using statistical flare models.

The comparison between model distributions and observed distributions is actually not an easy task, even if it is tempting to compare the power-law slope $`1.6`$ obtained by our model to the observed global slope of $`1.6`$. The first pitfall for comparison between statistical results of observations and models may be linked to the spatial resolution of observations ($`100\mathrm{km}`$ at best) compared to the dissipative scales of the system ($`100\mathrm{m}`$). It is clearly shown in this paper that the plot of the PDF of the magnetic energy dissipated in a given plane of the model has a shape quite different from the PDF of the same observable but for the whole simulation box. The bias here consists of a steepening of the power-law slope with an index greater than $`2`$. This result suggests that the limited instrumental resolution may be a source of error as well. The SOHO/EIT observations analyzed by Aletti et al. (2000) seem to illustrate quite well this interpretation : the pixel intensity distribution power-laws are steeper (with indices going up to $`5.6`$) at lower resolution. Most of the coronal structures in the quiet Sun are indeed smaller than the spatial resolution of EIT. Besides, the domain where the power-law is fitted on the PDF is reduced which leads to larger error on the value of the index. Reliable statistical results could be accessible with higher instrumental resolution but also by using mathematical tools like for example Pearson’s method (Podladchikova, 2002).

However, we think that resolution has not always such a dramatic effect on the slope of the PDFs, and that it is still interesting to model statistics of individual events (i.e. in one plane, in the case of our model). The convergence to a Gaussian when summing the energies of events before doing statistics, predicted by an argument lying on the Central Limit Theorem, may indeed be much slower in the case of observed, real micro-events than in the ideal case of independent events with low-moments distributions: the distributions of real events energies are much wider than the modeled distributions and their moments may be greater. As a result, depending on the observational conditions, there could still exist a quite wide power-law after summation (due to lack of resolution) of a large but not too large number of events, and the slope of this distribution may be still close to the slope of the original distribution of individual events. In this context the slope $`\zeta =1.6`$ we obtain is rather encouraging.

Another prediction of our model is that durations and energies of events scale like $`dE_idt_i^{1.76}`$, or, equivalently, $`dt_idE_i^{1/1.76}=dE_i^{0.57}`$. This exponent $`0.57`$ is in quite good agreement with Berghmans et al. (1998), who report that observed events durations scale like their radiative loss at the power $`0.5`$.

We should however emphasize that statistical flare models usually give energy dissipations, while observations give luminosities at some given wavelengths and the infered energies depend on models. It is therefore crucial in the future to develop models including the production of *observable* quantities, in order to provide stronger links between models and observations but also to quantify more precisely the weight of observational biases. The first agreements obtained during the last decade between statistical predictions made by theoretical models and observations are however very promising.

### 4.3 Summary and conclusion

In this paper we have presented a three-dimensional simplified model inspired by the RMHD equations whose first version was introduced by Einaudi & Velli (1999). This model mimics a coronal magnetic loop anchored in the photosphere whose footpoints are driven randomly by convective motions. The slow driving of the magnetic footpoints leads to storage of energy along the coronal loop and eventually to dissipation through impulsive events. The characteristics of the model are the following : (i) the model describes Alfvén wave propagation along a loop exactly ; the internal structure of the loop is described by a set of planes distributed along the axis and ending in the photosphere from which the information propagates ; (ii) the external forcing applied to the two boundary planes is expressed as a turbulent spectrum in Fourier space ; (iii) when the criterion of instability is satisfied, the dissipation of current density and vorticity takes place non-locally in Fourier space but still locally in physical space ; (iv) Fast Fourier Transforms (FFT) are implemented in the numerical code to use the dual physical/Fourier space.

A numerical study has allowed to quantify the role of the parameters, especially forcing, in the behavior of this model and in statistical properties of coronal events. The slope of event energy histograms was found to be almost constant, in accordance with a SOC-like “universal” behavior, and consistent with the values given by observations. Event durations statistics were performed, and correlations with event energies are also compatible with observations. Different possible observational biases were pointed out, all of them resulting in a narrower power-law range on histograms and in a steeper slope than in the statistics of all elementary events : a bias due to the limited spatial resolution, which gives a possible interpretation of recent observations made with the instrument EIT on board SoHO, and a bias due to limited observation durations.

Improvements are always possible to provide a better description of the physics. One can imagine some ad-hoc rules to obtain for example a correct picture of the reconnection process. But the probably most interesting (and difficult) study is about the incorporation of the non-linear dynamics in a more realistic way than the simple on-off mechanism of CA models, but without going directly to the MHD equations. From a pure observational point of view it seems crucial to have as precise as possible an estimate of the possible biases to determine the effective value of the power-law index $`\zeta `$ of the energy distribution. Indeed the confirmation of the sub-critical value of $`\zeta `$ could be a serious challenge to Parker’s hypothesis of coronal heating by a swarm of nanoflares.

###### Acknowledgements.

The authors acknowledge partial financial support from the PNST (Programme National Soleil–Terre) program of INSU (CNRS). M. Velli and G. Einaudi acknowledge support from MIUR contract MM02242342\_002. E. Buchlin thanks the Scuola Normale Superiore of Pisa, Italy, for support and accommodation. They thank Dr. M. Georgoulis for providing the bottom panel of figure 5. They appreciate the useful criticism of the anonymous referee.


text:
# The ground state of heavily-overdoped non-superconducting La2-xSrxCuO4

(Submitted to Physical Review Letters, 9 December 2002)

## Abstract

We report detailed thermodynamic and transport measurements for non-superconducting La1.7Sr0.3CuO4. Collectively, these data reveal that a highly-correlated Fermi-liquid ground state exists in La2-xSrxCuO4 beyond the superconducting dome, and confirm that charge transport in the cuprates is dominated at finite temperatures by intense electron-electron scattering.

The high-$`T_c`$ cuprates (HTC) have emerged as one of the most formidable challenges facing the theory of correlated electrons in solids orenstein00 . In particular, transport properties of the normal state present a number of uneasy paradoxes for the Fermi-liquid (FL) picture. Whilst it has been argued that conventional electron-phonon ($`e`$-$`ph`$) scattering might still account for the $`T`$-linear resistivity in optimally doped cuprates allen01 , it is generally assumed that proximity to a Mott insulator, and thus the unusual strength of on-site electronic repulsion, is the fundamental reason behind the more unusual aspects of HTC physics. This general assumption is backed up by a host of experimental observations which indicate that the anomalous behavior becomes even more pronounced with decreasing the density of carriers, that is, by moving towards the underdoped side. At sufficiently high carrier concentrations however, it has often been assumed that HTC eventually evolve into a conventional FL as the electron correlations become weaker and the system becomes more three-dimensional (3D).

Ironically, the persistence of robust superconductivity on the overdoped (OD) side of the phase diagram has been a major obstacle in the exploration of the metallic non-superconducting ground state in HTC. Indeed, supporting evidence for a FL ground state has only surfaced very recently with the experimental verification of the Wiedemann-Franz (WF) law in OD Tl2Ba2CuO6+δ (Tl2201) ($`T_c`$ $``$ 15K) proust02 . By suppressing superconductivity in a large magnetic field, Proust et al. observed the precise WF ratio $`\kappa _{ab}`$/$`\sigma _{ab}T`$ = $`L_0`$, where $`\kappa _{ab}`$ and $`\sigma _{ab}`$ are the in-plane thermal and electrical conductivities and the Lorenz number $`L_0`$ = 2.44 $`\times `$ 10-8 W$`\mathrm{\Omega }`$/K-2. Surprisingly however, and at odds with a conventional FL picture, the WF relation was found to co-exist with a large linear resistivity term extending down to 0K. This dichotomy raises the question whether the field-induced ’normal state’ in OD HTC, i.e. beyond $`H_{c2}`$, is the same as the ground state that would exist in the absence of a magnetic field, as it does in more conventional superconductors. Moreover, a clear understanding of the experimental situation in OD cuprates has often been compounded by their tendency to undergo phase separation.

In this Letter, we present in- ($`\rho _{ab}`$) and out-of-plane ($`\rho _c`$) resistivity, in- ($`\kappa _{ab}`$) and out-of-plane ($`\kappa _c`$) thermal conductivity, specific heat ($`C`$) and magnetic susceptibility ($`\chi `$) measurements on single-phase non-superconducting La2-xSrxCuO4 (LSCO) crystals ($`x`$ = 0.30) in which these various concerns are removed. Collectively, the data provide a consistent picture of La1.7Sr0.3CuO4 as a highly-correlated FL. The WF law is verified to within our experimental resolution, though in contrast to OD Tl2201, both $`\rho _{ab}`$ and $`\rho _c`$ exhibit strictly $`T`$2 behavior below 50K, with no additional linear term. Significantly, the Kadowaki-Woods ratio, linking the coefficient $`A`$ of the (in-plane) $`T`$2 resistivity and the square of the linear specific heat coefficient $`\gamma _0`$, is found to be anomalously enhanced, even compared with other strongly correlated metals. This latter observation reveals that intense $`e`$-$`e`$ interactions persist beyond the superconducting dome and sheds important new light on the evolution of $`\rho _{ab}`$($`T`$) across the HTC phase diagram.

Seven bar-shaped samples (typical dimensions 3 mm x 0.6 mm x 0.6 mm) were prepared for either $`ab`$-plane (A1-4) or $`c`$-axis (C1-3) measurements from a large La1.7Sr0.3CuO4 single crystal grown in an infra-red image furnace. The individual ingots were post-annealed together with the remaining boule under extremely high partial pressures of oxygen (400 atm) for 2 weeks at 900C to minimize oxygen deficiencies and to ensure good homogeneity within each crystal. Subsequent X-ray analysis revealed good crystallinity and no trace of superconductivity could be detected (resistively) down to 95mK (see top inset to Fig. 1), confirming that these crystals were indeed single-phase. $`\chi `$($`T`$) of the boule was measured with a commercial magnetometer, whilst $`C`$($`T`$) was measured between 0.6K and 10K in a relaxation calorimeter. In both cases, addenda contributions were measured independently and subtracted from the raw data. $`\kappa `$($`T`$) of A1, A2 and C1 were measured in a dilution refrigerator using three RuO2 chips employed as one heater and two thermometers. $`\rho `$($`T`$) measurements were made on all seven samples using a conventional ac four-probe method. For all transport measurements, current and voltage contacts were painted onto the crystals so as to short out any spurious voltage drops from orthogonal components of the conductivity tensor. Uncertainty in the geometrical factor was estimated to be $``$ 20$`\%`$. Finally, $`c`$-axis magnetoresistance $`\mathrm{\Delta }\rho _c`$/$`\rho _c`$ data were taken on C1 at the NHMFL in Florida.

Fig. 1 shows $`\rho `$($`T`$) of A1, A2 and C1 below 300K. (All crystals reported here showed similar behavior). Note the similar metallic $`T`$-dependencies observed along both orthogonal directions, the different residual resistivity ($`\rho _0`$) values for A1 and A2, and the strong upward curvature across the entire temperature range. The resistivity ratio $`\rho _c`$(C1)/$`\rho _{ab}`$(A1) rises from $``$ 65 at 300K to $``$ 80 as $`T`$ 0, presumably due to slight differences in their $`\rho _0`$ values. The limiting low-$`T`$ resistivities of C1 and A3 are reproduced in the top inset to Fig. 1. No current dependence was observed in $`\rho `$($`T`$) down to 0.1A/cm2, well below the typical critical current densities found in HTC. The lower inset shows $`\mathrm{\Delta }\rho _c`$/$`\rho _c`$ for C1 ($`Bab`$) plotted versus $`B`$2 at $`T`$ = 0.5K. Inserting the fit to the low-field data (see caption) into the Boltzmann transport equation for a quasi-2D FL (see, e.g. hussey96 ), we obtain an estimate of the in-plane electronic mean-free-path $`\mathrm{}_{ab}`$ 145 $`\AA `$. Finally, using $`k_F`$ = 0.55$`\AA ^1`$ for La1.7Sr0.3CuO4 ino02  and the ”isotropic-$`\mathrm{}`$” approximation mackenzie96a , we obtain $`\rho _{ab0}`$ = 21$`\mu \mathrm{\Omega }`$cm, in good agreement with $`\rho _0`$(A1).

In the top panel of Fig. 2, $`\rho `$($`T`$) of A3, A4 and C2 are plotted versus $`T`$1.6. In a previous study takagi92 , such non-integer power-law resistivities were shown to extend from the lowest temperatures up to 1000K. Whilst $`\rho `$($`T`$) of our crystals follows very closely a $`T`$1.6 dependence at elevated temperatures, the $`T`$-dependence clearly becomes stronger than $`T`$1.6 as $`T`$ drops below around 100K. Indeed, as shown in the bottom panel in Fig. 2, ($`\rho `$ \- $`\rho _0`$)($`T`$) = $`A`$$`T`$2 up to at least $`T`$ = 55K for both current directions (($`\rho `$ \- $`\rho _0`$)($`T`$) for A4 is shown in an inset for clarity). Significantly, $`A`$ has the same magnitude ($``$ 2.5 $`\pm `$ 0.1 $`n\mathrm{\Omega }`$cm/K2) for all in-plane crystals, even though $`\rho _0`$(A2,A3) $``$ 2 $`\rho _0`$(A1,A4). To our knowledge, this is the first time a pure $`T`$2 resistivity has been resolved in a hole-doped cuprate (recall that in OD Tl2201, ($`\rho \rho _0`$)($`T`$) = $`\alpha `$$`T`$ \+ $`AT`$2 proust02 ) and implies that in La1.7Sr0.3CuO4, charge transport can be understood in terms of an anisotropic 3D FL 2DFL .

In Fig. 3, the low-$`T`$ thermal conductivities of A1, A2 and C1 are plotted as $`\kappa `$/$`T`$ versus $`T`$2. Note that A1 and A2 are better conductors of heat, as expected since $`\sigma _c`$ $``$ $`\sigma _{ab}`$, and that $`\kappa `$(A1) $`>`$ $`\kappa `$(A2), reflecting their respective $`\rho _0`$ values. The sizeable phonon term ($`\kappa _{ph}`$) in HTC has often been a major obstacle to the determination of $`\kappa _{el}`$, the electronic contribution to $`\kappa `$($`T`$) and the chief quantity of interest. Typically, $`\kappa _{ph}`$ can only be unambiguously determined in the ballistic limit, below say 0.2K, where $`\kappa _{ph}`$ is simply proportional to $`\beta _3`$$`T`$3, the low-$`T`$ lattice heat capacity, as found here for C1 and indicated in the figure by a solid line. If $`\kappa _{ph}`$ were isotropic, a more insightful approach might be to measure both $`\kappa _{ab}`$ and $`\kappa _c`$ on identical crystals and make use of the large $`\sigma _{ab}`$/$`\sigma _c`$ anisotropy in HTC to delineate $`\kappa _{el}`$ from $`\kappa _{ab}`$ across the entire temperature range.

A quick inspection of Fig. 3 however reveals that $`\kappa _{ab}`$/$`T`$ and $`\kappa _c`$/$`T`$ are not simply offset from one another, implying either a $`T`$-dependent $`\kappa _{el}`$/$`T`$ at low-$`T`$ or anisotropy in $`\kappa _{ph}`$($`T`$) for in- and out-of-plane heat flow, due, for example, to anisotropy in the phonon velocities or $`e`$-$`ph`$ coupling. The fact though that $`\kappa `$/$`T`$(A1) and $`\kappa `$/$`T`$(A2) are shifted relative to one another (except at the very lowest temperatures) suggests the latter, and simply by scaling $`\kappa _c`$($`T`$) by a constant factor ($`f`$ = 1.6) and subtracting it from $`\kappa _{ab}`$($`T`$), we find, as shown in the inset to Fig. 3, that the resultant $`\kappa _{el}`$/$`T`$ is indeed constant for both crystals and equal to 1.07 and 0.50 mW/cm.K2 for A1 and A2 respectively (the dotted lines). This gives $`\rho _0`$ = $`L_0`$$`T`$/$`\kappa _{el}`$ = 23 (A1) and 49 (A2) $`\mu \mathrm{\Omega }`$cm, compared with the measured values of 19 and 42 $`\mu \mathrm{\Omega }`$cm. Given our somewhat simplistic way of extracting $`\kappa _{ph}`$($`T`$) from $`\kappa _{ab}`$($`T`$), it is encouraging to find that the WF law is satisfied for both crystals to within our experimental (geometrical) error. It should be stressed here that $`\sigma _{ab0}`$(LSCO) is almost an order of magnitude smaller than in OD Tl2201 proust02 , thus making determination of $`\kappa _{el}`$/$`T`$ more difficult.

The surprising feature of Fig. 3 is of course the strong downturn in $`\kappa `$(A1) as $`T`$ 0. (A slight downturn is also observed for A2 below 0.1K, though nowhere near as dramatic.) While slight deviations from the WF law can be attributed to additional anisotropy in $`\kappa _{ph}`$ not yet accounted for, there is no assumption on $`\kappa _{ph}`$($`T`$) which can erase the sudden vanishing of $`\kappa _{el}`$ below 0.2K for A1. Note that this downturn is robust both to changes in the thermal contacts and the crystal dimensions. Intriguingly, a vanishing of $`\kappa _{el}`$ has also been observed in the electron-doped cuprate Pr2-xCexCuO4 (PCCO) hill01 . By suppressing superconductivity with a high magnetic field, Hill et al. found that $`\kappa _{el}`$/$`T`$, whilst violating the WF law for $`T`$ $`>`$ 0.3K (in the sense that it clearly exceeds $`L_0`$/$`\rho `$($`T`$)), rapidly vanishes as $`T`$ falls below 0.3K. These two anomalous observations appear to have a common origin. The sample-dependent aspect of the downturn, however, suggests that it may be extrinsic in origin and unrelated to any exotic electronic behavior, particularly in view of the other, more conventional results obtained here for La1.7Sr0.3CuO4 above 0.2K. We reserve a detailed discussion on the origin of this anomalous downturn to a more complete and systematic investigation paglione .

The absence of superconductivity in La1.7Sr0.3CuO4, coupled with the emergence of a $`T`$2 resistivity at low $`T`$, gives us a unique opportunity to compare experimental manifestations of e-e correlations in La1.7Sr0.3CuO4 to other strongly correlated systems. Electronic correlations in a FL are known to lead to an enhancement in the quasi-particle effective mass. This effect can be detected through the resistivity ($`AT`$2), specific heat ($`\gamma _0`$$`T`$) and magnetic susceptibility ($`\chi _p`$). Empirical relationships that correlate these physical parameters have been found in a wide range of strongly-correlated metals; namely, the Kadowaki-Woods (KW) ratio ($`A/\gamma _0^2`$ = $`a_0`$ = 10$`{}_{}{}^{5}\mu \mathrm{\Omega }`$cm/(mJ/Kmol)2 kadowaki86 ; miyake89 ) and the Wilson ratio ($`R_w`$ = $`(\frac{\pi k_B^2}{3\mu _B^2}`$) $`\chi _p`$/$`\gamma _0`$, with $`R_w`$ 1 for a free electron gas and $``$ 2 for strongly-correlated fermions maeno97 ).

The lower inset to Fig. 4 shows $`C`$($`T`$) of the large crystal boule plotted as $`C`$/$`T`$ versus $`T`$2. $`C`$($`T`$) can be fitted over the whole temperature range (0.6K $`<`$ $`T`$ $`<`$ 10K) to the expression $`C`$ = $`\alpha T`$-2 \+ $`\gamma _0T`$ \+ $`\beta _3T`$3 \+ $`\beta _5T`$5 with $`\alpha `$ = 76 $`\mu `$J.K/mol., $`\gamma _0`$ = 6.9 mJ/mol.K2, $`\beta _3`$ = 0.1 mJ/mol.K4 and $`\beta _5`$ = 0.7 $`\mu `$J/mol.K6. The $`\alpha T`$-2 term represents the high temperature tail of a Schottky anomaly that invariably develops in HTC samples at low $`T`$ due to either a small concentration of isolated paramagnetic impurities, or to nuclear hyperfine or quadrupole splitting. The phononic $`\beta _3T^3`$ term is comparatively small and gives rise to an elevated $`\mathrm{\Theta }_D`$ = 550K, though $`\beta _3`$ may be made artificially low by a $`T`$-dependent electronic term cooper96 .

The magnitude of $`\gamma _0`$, though, is very robust in our measurements, is comparable to previous reports at this doping level cooper96  and gives a corresponding KW ratio, $`A/\gamma _0^2`$ 5$`a_0`$. This puts non-superconducting LSCO well off the so-called ’universal’ line for strongly correlated metals indicated in Fig. 4. A similarly enhanced KW ratio, i.e. $`A/\gamma _0^2`$ 5$`a_0`$, can also be inferred for OD Tl2201 ($`A`$ = 5.4 n$`\mathrm{\Omega }`$cmK-2 proust02 , $`\gamma _0`$ $``$ 10 mJ/mol.K2 wade95 ), and indirectly for OD PCCO, ($`A`$ $``$ 4 n$`\mathrm{\Omega }`$cmK-2 fournier98 , $`\gamma _0`$ $``$ 6.7 mJ/mol.K2 balci02 ). This observation reveals a new aspect of the physics of the cuprates and is the central result of this Letter. Whilst various explanations for the magnitude of the KW ratio in heavy fermion systems have been proposed, including a strong frequency dependence of the quasi-particle self-energy miyake89  and proximity to an antiferromagnetic instability takimoto96 , none have so far predicted deviations from $`a_0`$ as large as those found in the HTC. It is noted, however, that the effect of a strong momentum dependence in the scattering amplitude, thought to be a key feature of HTC hussey96 , has yet to be taken into account okabe98  and efforts to include such an effect could prove highly illuminating.

As one moves across the HTC phase diagram towards the Mott insulator at half-filling, one expects $`e`$-$`e`$ scattering to become even more intense. Thus, the gradual evolution from quadratic to linear resistivity in HTC must somehow reflect the growing strength of $`e`$-$`e`$ interactions as one approaches the Mott insulator from the metallic side. As an indication of how rapidly the $`e`$-$`e`$ scattering intensity might grow, we note here that $`\rho _{ab}`$ of La1.85Sr0.15CuO4 is $``$ 2 - 3 times larger than that of La1.7Sr0.3CuO4 at 300K, even though the Luttinger sum rule dictates that their carrier densities differ by only a few percent. Moreover, given the relatively minor changes observed in the phonon spectra as a function of doping, one does not expect the strength of the $`e`$-$`ph`$ interaction to change forcibly across the phase diagram. Thus, since $`e`$-$`ph`$ scattering appears to give a negligible contribution to $`\rho _{ab}`$($`T`$) in La1.7Sr0.3CuO4, it seems highly unlikely now that $`e`$-$`ph`$ scattering can in any way account for the linear resistivity that appears (in a narrow composition range) around optimal doping. In the light of these results, we conclude that chronic (Umklapp) $`e`$-$`e`$ scattering processes must dominate the normal state transport behavior of the cuprates across the *entire* accessible doping range. Such a scenario is consistent with the large increase in the quasi-particle lifetime below $`T_c`$ observed in thermal and electrical conductivity measurements yu92 ; bonn93 .

Finally, let us comment briefly on the Wilson ratio $`R_W`$. As shown in the upper inset to Fig. 4, $`\chi `$($`T`$) displays significant anisotropy with respect to the field orientation (believed to arise from anisotropic $`g`$-values terasaki92 ) and a strong enhancement at low $`T`$. Assuming a constant Pauli susceptibility $`\chi _p`$ vv  and an isotropic Curie-Weiss term (with $`\mathrm{\Theta }_C`$ $``$ 7.5K), we obtain an average $`\chi _p`$ $``$ 1.75$`\times `$ 10-4 emu/mol and hence, $`R_W`$ 2.0, consistent with values for other strongly correlated metals maeno97 . If however, all the enhancement in $`\chi `$($`T`$) is assumed to be intrinsic, i.e. due to an enhanced $`\chi _p`$($`T`$), we obtain $`R_W`$ $``$ 3.5. We note here that whilst a large enhancement in $`\chi _p`$($`T`$) at low-$`T`$ would imply a similar enhancement in $`\gamma `$($`T`$), such behavior would be difficult to discern from our $`C`$($`T`$) data due to the dominant phonon contribution.

In summary, detailed low-$`T`$ transport and thermodynamic measurements in La1.7Sr0.3CuO4 have revealed for the first time that LSCO develops a highly correlated FL ground state beyond the superconducting dome. The observed enhancement in the KW ratio suggests that $`e`$ scattering in non-superconducting LSCO is much more intense than previously imagined and raises serious question marks over previous claims that the ubiquitous linear resistivity observed in optimally-doped cuprates is a signature of (strong) $`e`$-$`ph`$ coupling.

We acknowledge stimulating discussions with Y. Ando, J. W. Loram, L. Taillefer and J. A. Wilson and experimental assistance from L. Balicas. This work was supported by the Alliance: Franco-British Partnership Program.


text:
# Estimate of {𝐵("→Kπν¯ν ")|}_{𝑆𝑀} from Standard Model fits to 𝜆_𝑡

## I Acknowledgments

We would like to thank P. Cooper, G. Isidori, D.E. Jaffe, P. Mackenzie, U. Nierste, L. Okun and A. Soni for useful discussions and comments and W. Marciano for the original stimulation to perform this type of calculation. This work was in part supported by the U.S. Department of Energy through contract #DE-AC02-98CH10886, in part by the Fermilab Particle Physics Division and in part by IHEP. L.G.L is grateful to the Fermilab Administration and Particle Physics Division for their hospitality and support for his stay at Fermilab during the preparation of this work.

## II Note

During the final preparation of this work for publication we found that Reference considered fitting for the apex of the UT from the CP-violating data only ($`|\epsilon _K|`$ and $`a_{\psi K}`$ ), as we do. However, Reference used ($`\overline{\rho },\overline{\eta }`$), which is dependent on $`|V_{cb}|`$ and is not as suitable for analysis of $`K\pi \nu \overline{\nu }`$ .