text:
# 1 Introduction

## 1 Introduction

A massless up-quark represents a simple and elegant solution to the strong CP problem. Consequently, the question of whether or not $`m_\mathrm{u}`$ is indeed zero has been the subject of much debate over many years (for a review see ref. ). Traditionally the problem is studied in the framework of Chiral Perturbation Theory (ChPT). Although the most recent estimates point to a non-zero value for the ratio $`m_\mathrm{u}/m_\mathrm{d}`$ , the situation is complicated by the presence of a hidden symmetry in the effective chiral Lagrangian . This so-called “Kaplan-Manohar ambiguity” implies that $`m_\mathrm{u}/m_\mathrm{d}`$ can only be constrained after supplementing ChPT with additional theoretical assumptions. Although the validity of the commonly used assumptions is plausible , it is clear that the question whether the up-quark is massless cannot be studied from first principles in ChPT.

More recently attention has focussed on lattice simulations to tackle this problem. A reliable, direct lattice calculation of $`m_\mathrm{u}`$, however, presents considerable difficulties, even on today’s massively parallel computers. It was therefore proposed to use a more indirect approach, based on a combination of lattice QCD and ChPT . The aim of this method is a lattice determination of the so-called “low-energy constants” in the effective chiral Lagrangian, whose precise values are required to constrain $`m_\mathrm{u}/m_\mathrm{d}`$ using chiral symmetry and phenomenological input. A variant of this proposal, which allows for a determination of the low-energy constants with good statistical accuracy was discussed in ref.  and tested in the quenched approximation.

Here we extend the study of  to QCD with two flavours of dynamical quarks. While this addresses the important issue of dynamical quark effects, it still does not correspond to the physical three-flavour case, and thus we are yet unable to give a final answer to the question in the title of this paper. Nevertheless, our study represents an important step in an ultimately realistic treatment of the problem, by studying the dependence of the low-energy constants on the number of flavours. If our results can be taken over to the physical case without large modifications – and there are indications that this is not unreasonable – then the possibility of a massless up-quark is strongly disfavoured.

In Section 2 we briefly review the Kaplan-Manohar ambiguity and its relevance for the value of $`m_\mathrm{u}`$. Section 3 contains details of our lattice simulations, whose results are described in Section 4. In Section 5 we discuss the implications of our findings and present an outlook to future work.

## 2 Low-energy constants and $`m_\mathrm{u}=0`$

In order to make this paper self-contained, we briefly review the implications of the Kaplan-Manohar ambiguity for the ratio $`m_\mathrm{u}/m_\mathrm{d}`$. The strategy to address the problem in lattice simulations will then become clear. A more complete discussion can be found in refs. .

A determination of $`m_\mathrm{u}/m_\mathrm{d}`$ in ChPT which is able to distinguish between a massless and a massive up-quark requires precise knowledge of the first-order mass correction term $`\mathrm{\Delta }_\mathrm{M}`$. At order $`p^4`$ in the chiral Lagrangian it is given by

$$\mathrm{\Delta }_\mathrm{M}=\frac{m_\mathrm{K}^2m_\pi ^2}{(4\pi F_\pi )^2}(2\alpha _8\alpha _5)+\text{chiral logs},$$
(1)
where $`F_\pi =93.3\mathrm{MeV}`$ is the pion decay constant, and $`\alpha _5`$, $`\alpha _8`$ are low-energy constants, whose values have to be determined from phenomenology. Throughout this paper we adopt a convention in which the low-energy constants $`\alpha _i`$ are related to the corresponding constants $`L_i`$ of ref.  through $`\alpha _i=8(4\pi )^2L_i`$. Furthermore, we always quote low-energy constants in the $`\overline{\mathrm{MS}}`$-scheme at scale $`\mu =4\pi F_\pi `$.

The value of $`\alpha _5`$ can be extracted from the ratio of pseudoscalar decay constants, $`F_\mathrm{K}/F_\pi `$, and is given by

$$\alpha _5=0.5\pm 0.6.$$
(2)
By contrast, there is no direct phenomenological information on $`\alpha _8`$ or the combination $`(2\alpha _8\alpha _5)`$. Although $`\alpha _8`$ is contained in the correction to the Gell-Mann–Okubo formula, i.e.

$$\mathrm{\Delta }_{\mathrm{GMO}}=\frac{m_\mathrm{K}^2m_\pi ^2}{(4\pi F_\pi )^2}(\alpha _512\alpha _76\alpha _8)+\text{chiral logs},$$
(3)
its determination requires prior knowledge of $`\alpha _7`$. It is at this point that the Kaplan-Manohar (KM) ambiguity becomes important. It arises from the observation that a simultaneous transformation of the quark masses

$$m_\mathrm{u}m_\mathrm{u}+\lambda m_\mathrm{d}m_\mathrm{s},m_\mathrm{d}m_\mathrm{d}+\lambda m_\mathrm{s}m_\mathrm{u},m_\mathrm{s}m_\mathrm{s}+\lambda m_\mathrm{u}m_\mathrm{d},$$
(4)
and coupling constants according to

$$\alpha _6\alpha _6+\lambda \frac{(4\pi F_0)^2}{4B_0},\alpha _7\alpha _7+\lambda \frac{(4\pi F_0)^2}{4B_0},\alpha _8\alpha _8\lambda \frac{(4\pi F_0)^2}{2B_0},$$
(5)
leaves the effective chiral Lagrangian invariant. Here, $`\lambda `$ is an arbitrary parameter, and $`F_0`$, $`B_0`$ are coupling constants in the lowest-order chiral Lagrangian.222$`F_0`$ coincides with $`F_\pi `$ at lowest order. Thus, chiral symmetry cannot distinguish between different sets of quark masses and coupling constants, which are related through eqs. (4) and (5). Indeed the correction $`\mathrm{\Delta }_{\mathrm{GMO}}`$ is invariant under the above transformations. The value of $`\alpha _7`$ can be fixed by invoking additional theoretical assumptions, such as the validity of large-$`N_c`$ arguments for $`N_c=3`$. In accordance with these assumptions, Leutwyler  constrained the correction term $`\mathrm{\Delta }_\mathrm{M}`$ to be small and positive:

$$0<\mathrm{\Delta }_\mathrm{M}0.13.$$
(6)
This gives $`m_\mathrm{u}/m_\mathrm{d}=0.553(43)`$ and hence a non-zero value for $`m_\mathrm{u}`$. The “standard” values for $`\alpha _7`$ and $`\alpha _8`$ which are compatible with eq. (6) are

$$\alpha _7=0.5\pm 0.25,\alpha _8=0.76\pm 0.4.$$
(7)
In view of the importance of the strong CP problem, one may regard any analysis based on theoretical assumptions beyond chiral symmetry as insufficient. In particular, since the uncertainties in the estimates for $`\alpha _7`$ and $`\alpha _8`$ are quite large, the possibility that $`m_\mathrm{u}=0`$ does not appear to be ruled out completely. A massless up-quark would require

$$\alpha _7=0.25\pm 0.25,\alpha _8=0.9\pm 0.4,$$
(8)
resulting in a large, negative first-order correction $`\mathrm{\Delta }_\mathrm{M}`$. In order to decide which scenario is realised and to pin down the value of $`m_\mathrm{u}`$ one has to replace the theoretical assumptions by a solution of the underlying theory of QCD.

The KM ambiguity implies that the low-energy constants $`\alpha _7`$, $`\alpha _8`$ (and $`\alpha _6`$) can be determined from chiral symmetry and phenomenology only if the physical quark masses are known already. At this point it is important to realise that QCD is not afflicted with the KM ambiguity, and that the formalism of ChPT also holds for unphysical quark masses. Since quark masses are input parameters in lattice simulations of QCD, their relations to hadronic observables need not be known a priori. Hence, the low-energy constants can be determined by studying pseudoscalar meson masses and matrix elements for unphysical quark masses and fitting their quark mass dependence to the expressions found in ChPT. In this way it is possible to determine $`\alpha _5`$ and – more importantly – the combination $`(2\alpha _8\alpha _5)`$ directly in lattice simulations.

## 3 Lattice setup and simulation details

In ref.  it was shown how the low-energy constants can be extracted from lattice data for suitably defined ratios of pseudoscalar masses and matrix elements, $`R_\mathrm{M}`$ and $`R_\mathrm{F}`$. Here we repeat their definitions in order to explain the necessary notation. For more details we refer the reader to the original paper .

The actual determination of the low-energy constants proceeds by studying the mass dependence of $`R_\mathrm{M}`$ and $`R_\mathrm{F}`$ around some reference quark mass $`m_{\mathrm{ref}}`$. As pointed out before, $`m_{\mathrm{ref}}`$ does not have to coincide with a physical quark mass , as long as it is small enough for ChPT to be applicable.

In this paper we will be concerned with simulations of “partially quenched” QCD, where valence and sea quarks can have different masses. Let us therefore consider a pseudoscalar meson with valence quark masses $`m_1`$ and $`m_2`$ at a fixed value of the sea quark mass $`m^{\mathrm{sea}}`$. In  the dimensionless mass parameters

$$y=\frac{B_0(m_1+m_2)}{(4\pi F_0)^2},y_{\mathrm{ref}}=\frac{2B_0m_{\mathrm{ref}}}{(4\pi F_0)^2},x=y/y_{\mathrm{ref}},$$
(9)
were defined, as well as the ratios

$`R_\mathrm{M}(x)`$ $`=`$ $`\left({\displaystyle \frac{F_{\mathrm{PS}}(y)}{G_{\mathrm{PS}}(y)}}\right)/\left({\displaystyle \frac{F_{\mathrm{PS}}(y_{\mathrm{ref}})}{G_{\mathrm{PS}}(y_{\mathrm{ref}})}}\right)`$
$`R_\mathrm{F}(x)`$ $`=`$ $`F_{\mathrm{PS}}(y)/F_{\mathrm{PS}}(y_{\mathrm{ref}}),x=y/y_{\mathrm{ref}}.`$ (10)
Here $`F_{\mathrm{PS}}`$ is the pseudoscalar decay constant, and $`G_{\mathrm{PS}}`$ is the matrix element of the pseudoscalar density between a pseudoscalar state and the vacuum. The parameter $`x`$ denotes the fraction of the reference quark mass at which the ratios $`R_\mathrm{M}`$ and $`R_\mathrm{F}`$ are considered.

As emphasised in , $`R_\mathrm{M}`$ and $`R_\mathrm{F}`$ are well-suited to extract the low-energy constants, since ratios are usually obtained with high statistical accuracy in numerical simulations. Furthermore, any renormalisation factors associated with $`F_{\mathrm{PS}}`$ and $`G_{\mathrm{PS}}`$ drop out, so that $`R_\mathrm{M}`$ and $`R_\mathrm{F}`$ can be readily extrapolated to the continuum for every fixed value of $`x`$. As has been shown in the quenched approximation , discretisation errors in $`R_\mathrm{M}`$ and $`R_\mathrm{F}`$ are very small, so that good control over lattice artefacts is achieved. This issue is important, since the determination of low-energy constants requires an unambiguous separation of the mass dependence from effects of non-zero lattice spacing.

We now give the expressions for $`R_\mathrm{M}`$ and $`R_\mathrm{F}`$ in partially quenched QCD which are relevant for our study. They were obtained using the results of ref.  and are listed in appendix A of . Our reference point $`y_{\mathrm{ref}}`$ was always defined at

$$m_1=m_2=m^{\mathrm{sea}}=m_{\mathrm{ref}},y_{\mathrm{ref}}=\frac{2B_0m_{\mathrm{ref}}}{(4\pi F_0)^2}.$$
(11)
In order to map out the mass dependence of $`R_\mathrm{M}`$ and $`R_\mathrm{F}`$ for a fixed value of $`m^{\mathrm{sea}}`$, we have considered the cases labelled “VV” and “VS1” in ref. . The first uses degenerate valence quarks and is defined by

$$\text{VV:}m_1=m_2=xm_{\mathrm{ref}},m^{\mathrm{sea}}=m_{\mathrm{ref}},$$
(12)
which leads to the expressions

$`R_\mathrm{M}^{\mathrm{VV}}(x)`$ $`=`$ $`1{\displaystyle \frac{1}{N_\mathrm{f}}}y_{\mathrm{ref}}\left[(2x1)\mathrm{ln}x+2(x1)\mathrm{ln}y_{\mathrm{ref}}\right]`$ (13)
$`y_{\mathrm{ref}}(x1)\left[(2\alpha _8\alpha _5)+\frac{1}{N_\mathrm{f}}\right]`$
$`R_\mathrm{F}^{\mathrm{VV}}(x)`$ $`=`$ $`1{\displaystyle \frac{N_\mathrm{f}}{4}}y_{\mathrm{ref}}\left[(x+1)\mathrm{ln}\left(\frac{1}{2}(x+1)\right)+(x1)\mathrm{ln}y_{\mathrm{ref}}\right]`$ (14)
$`+y_{\mathrm{ref}}(x1)\frac{1}{2}\alpha _5,`$
where $`N_\mathrm{f}`$ denotes the number of dynamical quark flavours. The case labelled “VS1”, based on non-degenerate valence quarks, is defined by

$$\text{VS1:}m_1=xm_{\mathrm{ref}},m_2=m^{\mathrm{sea}}=m_{\mathrm{ref}}.$$
(15)
According to Table 1 of ref.  the expressions for the ratios are then given by

$`R_\mathrm{M}^{\mathrm{VS1}}(x)`$ $`=`$ $`1{\displaystyle \frac{1}{N_\mathrm{f}}}y_{\mathrm{ref}}\left[x\mathrm{ln}x+(x1)\mathrm{ln}y_{\mathrm{ref}}\right]y_{\mathrm{ref}}(x1)\frac{1}{2}(2\alpha _8\alpha _5)`$ (16)
$`R_\mathrm{F}^{\mathrm{VS1}}(x)`$ $`=`$ $`1{\displaystyle \frac{N_\mathrm{f}}{8}}y_{\mathrm{ref}}\left[(x+1)\mathrm{ln}\left(\frac{1}{2}(x+1)\right)+(x1)\mathrm{ln}y_{\mathrm{ref}}+\frac{2}{N_\mathrm{f}^2}\mathrm{ln}x\right]`$ (17)
$`+y_{\mathrm{ref}}(x1)\frac{1}{4}\left(\alpha _5+\frac{1}{N_\mathrm{f}}\right).`$
Our simulations were performed for $`N_\mathrm{f}=2`$ flavours of dynamical, O($`a`$) improved Wilson fermions. The value of the bare coupling was set to $`\beta =6/g_0^2=5.2`$. The improvement coefficient $`c_{\mathrm{sw}}`$, which multiplies the Sheikholeslami-Wohlert term in the fermionic action, was taken from the interpolating formula of ref. . Here we considered a single value of the sea quark mass, corresponding to a hopping parameter $`\kappa ^{\mathrm{sea}}=0.1355`$. For this choice of parameters we generated 208 dynamical gauge configurations on a lattice of size $`16^332`$. For further details we refer to . Here we only mention that the hadronic radius $`r_0`$ defined through the force between static sources  has been determined as $`r_0/a=5.041(40)`$ . For $`r_0=0.5\mathrm{fm}`$ this implies that the lattice spacing in physical units is $`a=0.099(1)\mathrm{fm}`$.

We have computed quark propagators for valence quarks with hopping parameters $`\kappa ^{\mathrm{val}}=0.1340`$, 0.1345, 0.1350, 0.1355 and 0.1358. In addition to local operators, we have also calculated propagators for fuzzed sinks and/or sources, using the procedure of . In the pseudoscalar channel we employed both the pseudoscalar density and the temporal component of the axial current as interpolating operators. These two types were used together with the different combinations of fuzzed and local propagators to construct a $`4\times 4`$ matrix correlator for pseudoscalar mesons. By performing factorising fits using an ansatz that incorporates the ground state and the first excitation, we were able to extract the pseudoscalar mass $`m_{\mathrm{PS}}`$, as well as the matrix elements of the axial current and the pseudoscalar density, i.e.

$$\zeta _\mathrm{A}=0|A_0|\text{PS},\zeta _\mathrm{P}=0|P|\text{PS}.$$
(18)
In order to be consistent with O($`a`$) improvement the amplitudes $`\zeta _\mathrm{A}`$ and $`\zeta _\mathrm{P}`$ must be related to the matrix elements of the improved currents and densities at non-zero quark mass. Using the definitions of  it is then easy to see that the (unrenormalised) pseudoscalar decay constant is given by

$$F_{\mathrm{PS}}=(1+b_\mathrm{A}am_\mathrm{q})\left\{\frac{\zeta _\mathrm{A}}{m_{\mathrm{PS}}}+c_\mathrm{A}\frac{\zeta _\mathrm{P}}{m_{\mathrm{PS}}}\mathrm{sinh}(am_{\mathrm{PS}})\right\},$$
(19)
whereas $`G_{\mathrm{PS}}`$ is given by

$$G_{\mathrm{PS}}=(1+b_\mathrm{P}am_\mathrm{q})\zeta _\mathrm{P}.$$
(20)
The improvement coefficients $`b_\mathrm{A}`$, $`b_\mathrm{P}`$ and $`c_\mathrm{A}`$ were computed in one-loop perturbation theory  in the bare coupling, since non-perturbative estimates are not available at present. For the improvement coefficient $`c_\mathrm{A}`$ this procedure yields

$$c_\mathrm{A}=0.0087.$$
(21)
For degenerate valence quarks the quark mass $`m_\mathrm{q}`$ is given by

$$am_\mathrm{q}=\frac{1}{2}\left(\frac{1}{\kappa ^{\mathrm{val}}}\frac{1}{\kappa _{\mathrm{crit}}}\right),$$
(22)
where we have inserted $`\kappa _{\mathrm{crit}}=0.13693`$ for the critical hopping parameter, as estimated in . Finally we note that the current quark mass $`m`$ defined through the PCAC relation in the O($`a`$) improved theory is obtained from $`\zeta _\mathrm{A}`$, $`\zeta _\mathrm{P}`$ and $`m_{\mathrm{PS}}`$ via

$$am=\frac{\zeta _\mathrm{A}}{2\zeta _\mathrm{P}}\mathrm{sinh}(am_{\mathrm{PS}})+\frac{1}{2}c_\mathrm{A}\mathrm{sinh}^2(am_{\mathrm{PS}}).$$
(23)
All our statistical errors were obtained using a bootstrap procedure .

## 4 Results

Our results for pseudoscalar masses, matrix elements and current quark masses are listed in Table 1. Compared with  the numbers for the pseudoscalar masses reported here may differ by up to one standard deviation, as a result of using a larger matrix correlator in the fitting procedure.

Following eq. (11) we have defined the mass $`m_{\mathrm{ref}}`$ at the reference point for $`\kappa _1^{\mathrm{val}}=\kappa _2^{\mathrm{val}}=\kappa ^{\mathrm{sea}}=0.1355`$, which corresponds to $`(r_0m_{\mathrm{PS}})^2|_{m=m_{\mathrm{ref}}}=2.092`$. Using the leading-order relations $`m_{\mathrm{PS}}^2=2B_0m`$ and $`F_0=F_\pi =93.3\mathrm{MeV}`$ this implies

$$y_{\mathrm{ref}}=0.2370.$$
(24)
It is instructive to compare these values with those of the previous quenched study , where $`y_{\mathrm{ref}}=0.3398`$, $`(r_0m_{\mathrm{PS}})^2|_{m=m_{\mathrm{ref}}}=3`$, and $`m_{\mathrm{ref}}m_\mathrm{s}`$. Hence, the value of $`m_{\mathrm{ref}}`$ employed in this paper is smaller by 30%, such that $`m_{\mathrm{ref}}0.7m_\mathrm{s}`$. A smaller value of $`m_{\mathrm{ref}}`$ is clearly desirable for our purpose, since the predictions of ChPT are expected to hold more firmly for smaller masses. Nevertheless, our sea quarks are still relatively heavy, as signified by the ratio $`(m_{\mathrm{PS}}/m_\mathrm{V})|_{m=m^{\mathrm{sea}}}0.58`$ , which is to be compared to $`m_\pi /m_\rho =0.169`$.

In the spirit of partially quenched QCD we have considered valence quarks that are lighter than the sea quarks. We were thus able to extend the quark mass range down to about $`m_\mathrm{s}/2`$, which is quite a bit below the smallest mass reached in . This may be an indication that the inclusion of dynamical quarks alleviates the problem of exceptional configurations, which precludes attempts to work at very small quark masses in the quenched approximation. However, our attempts to push to valence quark masses below $`m_\mathrm{s}/2`$ have proved unsuccessful, due to the appearance of exceptional configurations for $`m^{\mathrm{val}}/m^{\mathrm{sea}}\stackrel{<}{}\mathrm{\hspace{0.33em}0.7}`$. More precisely, we observed large statistical fluctuations in hadron correlators computed from quark propagators at $`\kappa ^{\mathrm{val}}=0.1360`$, despite the fact that the inversion algorithm converged.

The results in Table 1 can now be used to compute $`R_\mathrm{M}`$ and $`R_\mathrm{F}`$ through eqs. (20), (19) and (10). In order to map out their quark mass dependence in some detail, we have performed local interpolations of the results for $`F_{\mathrm{PS}}`$ and $`F_{\mathrm{PS}}/G_{\mathrm{PS}}`$ to 20 different values of the dimensionless mass parameter $`x`$ in the range $`0.7x2.6`$, separated by increments of 0.1.

Since we only have one $`\beta `$-value we cannot extrapolate $`R_\mathrm{M}`$, $`R_\mathrm{F}`$ to the continuum limit for fixed $`x`$. Unlike ref. , where such extrapolations could be performed, we must compare our data to ChPT at non-zero lattice spacing. Our estimates for the low-energy constants are therefore subject to an unknown discretisation error. It is reasonable to assume, however, that cutoff effects in $`R_\mathrm{M}`$ and $`R_\mathrm{F}`$ are fairly small, owing to cancellations of lattice artefacts of similar size between numerator and denominator. We will return to this issue below, when we discuss our final estimates for the low-energy constants.

In order to extract the low-energy constants we have restricted the $`x`$-interval to $`0.7x1.1`$, thereby seeking to maximise the overlap with the domain of applicability of ChPT, whilst maintaining a large enough interval to check the stability of our results. Estimates for the low-energy constants were obtained by fitting the data for $`R_\mathrm{M}`$ and $`R_\mathrm{F}`$ to the corresponding expressions for the “VV” and “VS1” cases listed in Section 3. Since $`R_\mathrm{F}`$ is linear in $`\alpha _5`$ one can obtain this low-energy constant also from simple algebraic expressions involving the difference $`R_\mathrm{F}(x_1)R_\mathrm{F}(x_2)`$ for two distinct arguments, $`x_1`$ and $`x_2`$. A similar relation can be used to compute $`(2\alpha _8\alpha _5)`$ from $`R_\mathrm{M}(x_1)R_\mathrm{M}(x_2)`$. We have checked that both methods give consistent results and obtain

VV: $`\alpha _5^{(2)}=1.20\begin{array}{c}+\text{11}\hfill \\ \text{16}\hfill \end{array},(2\alpha _8^{(2)}\alpha _5^{(2)})=0.36\begin{array}{c}+\text{10}\hfill \\ \text{10}\hfill \end{array},`$ (29)
VS1: $`\alpha _5^{(2)}=1.22\begin{array}{c}+\text{11}\hfill \\ \text{16}\hfill \end{array},(2\alpha _8^{(2)}\alpha _5^{(2)})=0.36\begin{array}{c}+\text{10}\hfill \\ \text{12}\hfill \end{array},`$ (34)
where the errors are purely statistical. From here on we also indicate the number of dynamical quark flavours as a superscript, to distinguish these estimates from the corresponding ones in the quenched and three-flavour cases. Our results can now be inserted back into the expressions for $`R_\mathrm{M}`$ and $`R_\mathrm{F}`$. The resulting curves are plotted together with the data in Fig. 1.

It is striking that the data for $`R_\mathrm{F}`$ in the VV case are described remarkably well over the whole mass range, despite the fact that $`\alpha _5^{(2)}`$ has only been determined for $`x1.1`$. The qualitative behaviour of $`R_\mathrm{F}`$ – which features a slight curvature – is thus rather well modelled by eqs. (14) and (17), which include a linear term as well as chiral logarithms. By contrast, there are no logarithmic contributions to $`R_\mathrm{F}`$ in quenched ChPT, and the expected purely linear behaviour has indeed been observed in the data . Of course, higher-order terms in the quark mass could in principle produce a curvature, and therefore these observations do not provide unambiguous evidence for chiral logarithms. Nevertheless, it is remarkable that the clear distinction between the expressions for $`R_\mathrm{F}`$ in partially quenched and quenched ChPT (i.e. the presence, respectively absence of chiral logarithms) is accompanied by corresponding qualitative differences in the numerical data.

Unlike $`R_\mathrm{F}`$ the ratio $`R_\mathrm{M}`$ is not described well for larger masses, which may signal a breakdown of the chiral expansion for this quantity for masses not much larger than $`m_\mathrm{s}`$. It is therefore conceivable that higher orders in ChPT affect the extraction of $`(2\alpha _8^{(2)}\alpha _5^{(2)})`$. However, without access to quark masses that are substantially lower than our simulated ones, it is not easy to quantify reliably any uncertainty due to neglecting higher orders.

One way to examine the influence of higher orders is to extract the low-energy constants from a mass interval of fixed length, $`(x_{\mathrm{max}}x_{\mathrm{min}})=0.4`$, which is then shifted inside an extended range of $`0.7x1.5`$. The spread of results so obtained then serves as an estimate of the systematic errors incurred by neglecting higher orders. We note that $`x=1.5`$ corresponds to a quark mass slightly larger than $`m_\mathrm{s}`$. For $`\alpha _5^{(2)}`$ such a procedure yields only a small variation of $`\pm 0.05`$. This is not surprising, since $`R_\mathrm{F}`$ is modelled very well over the entire mass range. By contrast, the spread of results obtained for $`(2\alpha _8^{(2)}\alpha _5^{(2)})`$ is as large as $`\pm 0.15`$. Whether or not these numbers represent realistic estimates of the actual uncertainty cannot be decided at this stage. In order to be more conservative we have decided to quote a systematic error of $`\pm 0.2`$ for all low-energy constants. We note that this level of uncertainty due to neglecting higher orders was also quoted in the quenched case , where quark masses were slightly larger.

Since we do not have enough data to extrapolate $`R_\mathrm{F}`$ and $`R_\mathrm{M}`$ to the continuum limit we also have to estimate a systematic error due to cutoff effects. As explained above, however, we expect such effects to be small. In order to get an idea of the typical size of discretisation errors we have looked again at quenched data obtained at $`\beta =5.93`$  and 6.0 , for which the lattice spacing in physical units is roughly the same as in our dynamical simulations ($`a0.1\mathrm{fm}`$). For both $`\beta =5.93`$ and 6.0 the results for $`R_\mathrm{M}`$ and $`R_\mathrm{F}`$ are mostly consistent within errors with the corresponding values in the continuum limit (see, for instance, Fig. 1 in ). Furthermore, low-energy constants extracted for $`a0.1\mathrm{fm}`$ differ from the results in the continuum limit by less than one standard deviation. Although these findings cannot be taken over literally to the dynamical case without direct verification, they nevertheless indicate that lattice artefacts are small enough such that does not have to expect large distortions in our estimates for the low-energy constants. In order to take account of these observations we have decided to quote an additional systematic error due to lattice artefacts, which is as large as the statistical error.

Since non-perturbative estimates for the improvement coefficients $`b_\mathrm{A}`$, $`b_\mathrm{P}`$ and $`c_\mathrm{A}`$ are not available for $`N_\mathrm{f}=2`$, one may be worried that there are large uncancelled lattice artefacts of order $`a`$ in our data. We have addressed this issue by studying the influence of different choices for improvement coefficients on our results. To this end we have repeated the complete analysis using non-perturbative values for $`c_\mathrm{A}`$ and the combination $`b_\mathrm{A}b_\mathrm{P}`$ obtained in the quenched approximation  at a similar value of the lattice spacing, $`a0.1\mathrm{fm}`$. We found that the resulting variation in the estimates for $`\alpha _5^{(2)}`$ and $`(2\alpha _8^{(2)}\alpha _5^{(2)})`$ is typically a factor 10 smaller than the statistical error. Thus we conclude that the influence of improvement coefficients on our results is very weak indeed.

As a final comment we point out that we have not taken finite volume effects into account in our error estimates, because different lattice sizes were not considered in our study (unlike in earlier simulations ). However, since $`Lm_{\mathrm{PS}}=4.59`$ at the reference point and $`Lm_{\mathrm{PS}}=3.68`$ at the lightest valence quark mass, one may not be totally convinced that such effects may be entirely neglected. We stress, though, the the definition of $`R_\mathrm{M}`$ and $`R_\mathrm{F}`$ implies that only the relative finite-size effects between hadronic quantities is relevant. Thus, as long as the mass parameter $`x`$ does not differ too much from unity, one can reasonably expect that finite-volume effects largely cancel in the ratios used to determine the low-energy constants. Finite-size effects in pseudoscalar masses and decay constants have also been studied in ChPT . These calculations indicate that the typical relative finite-volume effect in $`m_{\mathrm{PS}}`$ and $`F_{\mathrm{PS}}`$ between our reference point and the smallest quark mass is less than 1%.

## 5 Discussion and outlook

After combining the different systematic errors in quadrature we obtain as our final results in two-flavour QCD:

$`\alpha _5^{(2)}`$ $`=`$ $`1.22_{0.16}^{+0.11}\text{(stat)}_{0.26}^{+0.23}\text{(syst)}`$ (35)
$`(2\alpha _8^{(2)}\alpha _5^{(2)})`$ $`=`$ $`0.36\pm 0.10\text{(stat)}\pm 0.22\text{(syst)}`$ (36)
$`\alpha _8^{(2)}`$ $`=`$ $`0.79_{0.07}^{+0.05}\text{(stat)}\pm 0.21\text{(syst)},`$ (37)
where eqs. (35) and (36) have been combined to produce the result for $`\alpha _8^{(2)}`$.

We can now investigate the dependence of the low-energy constants on the number of dynamical quark flavours. In the quenched approximation  (i.e. for $`N_\mathrm{f}=0`$) it was found that333When extracting the result for $`\alpha _8^{(0)}`$ it was assumed that the coefficients multiplying quenched chiral logarithms were set to $`\delta =0.12,\alpha _\mathrm{\Phi }=0.0`$.

$`\alpha _5^{(0)}`$ $`=`$ $`0.99\pm 0.06\text{(stat)}\pm 0.2\text{(syst)}`$ (38)
$`\alpha _8^{(0)}`$ $`=`$ $`0.67\pm 0.04\text{(stat)}\pm 0.2\text{(syst)}.`$ (39)
A comparison with eq. (35) and (37) then shows that $`\alpha _5^{(2)}`$ and $`\alpha _8^{(2)}`$ are larger than their quenched counterparts by 23% and 18% respectively. We can thus conclude that the $`N_\mathrm{f}`$-dependence of the low-energy constants is fairly weak: variations between the quenched and two-flavour theories are about as large as the error due to neglecting higher orders.

Although there is a priori no reason why the weak $`N_\mathrm{f}`$-dependence should extend to the physical three-flavour case, it is still instructive to compare eqs. (35) and (37) with phenomenological values of the low-energy constants. It then becomes obvious from eqs. (2) and (7) that our results for $`\alpha _5^{(2)}`$ and $`\alpha _8^{(2)}`$ are compatible with the standard estimates found in the literature. By contrast, our numerical data for $`R_\mathrm{F}`$ and $`R_\mathrm{M}`$ suggest that a large negative value for $`\alpha _8`$, which is required for the scenario of $`m_\mathrm{u}=0`$ (see eq. (8)), is practically ruled out. Thus, provided that the quark mass behaviour in the physical three-flavour case is not fundamentally different, the possibility of a massless up-quark is strongly disfavoured.

By how much does one expect the mass dependence of $`R_\mathrm{M}`$ and $`R_\mathrm{F}`$ to differ between $`N_\mathrm{f}=2`$ and 3 ? Ultimately this must be answered by a direct simulation of the three-flavour case. For the time being we have to be content with the following gedanken simulation. Suppose that we had analysed our $`N_\mathrm{f}=2`$ data under the erroneous assumption that they had been obtained in the physical three-flavour case. We would then have set $`N_\mathrm{f}=3`$ in eqs. (14) and (17) to extract $`\alpha _5^{(3)}`$, giving $`\alpha _5^{(3)}=0.98\begin{array}{c}+\text{11}\hfill \\ \text{16}\hfill \end{array}\begin{array}{c}+\text{23}\hfill \\ \text{26}\hfill \end{array}`$. This value can be inserted into the expression for $`F_\mathrm{K}/F_\pi `$ in the physical theory , to yield

$$\frac{F_\mathrm{K}}{F_\pi }=1.247_{0.013}^{+0.009}\text{(stat)}_{0.021}^{+0.019}\text{(syst)},$$
(40)
which is in fair agreement with the experimental result $`F_\mathrm{K}/F_\pi =1.22\pm 0.01`$. This shows that the experimental value can only be reproduced if the quark mass dependence of $`R_\mathrm{F}`$ in the physical case is not much different from that encountered in our simulations for $`N_\mathrm{f}=2`$. In other words, it is reasonable to assume that the mass dependence of $`R_\mathrm{F}`$ is only weakly distorted by neglecting the dynamical quark effects due to a third flavour.

Similarly we can apply the expressions for $`R_\mathrm{M}`$ for $`N_\mathrm{f}=3`$ to our data, which gives $`(2\alpha _8^{(3)}\alpha _5^{(3)})=0.20\begin{array}{c}+\text{10}\hfill \\ \text{12}\hfill \end{array}\begin{array}{c}+\text{22}\hfill \\ \text{23}\hfill \end{array}`$. The corresponding curve is shown in Fig. 2. The first-order mass correction $`\mathrm{\Delta }_\mathrm{M}`$ is then obtained as

$$\mathrm{\Delta }_\mathrm{M}=0.04_{0.06}^{+0.05}\text{(stat)}\pm 0.11\text{(syst)},$$
(41)
which is consistent with Leutwyler’s estimate (see eq. (6)). We emphasise that this does not represent a reliable result for $`\mathrm{\Delta }_\mathrm{M}`$ derived from first principles. Nevertheless, the above discussion shows that there are examples which support the idea that the gross features of the mass dependence do not differ substantially in the two- and three-flavour cases. On the basis of this assumption one may conclude that the correction factor $`\mathrm{\Delta }_\mathrm{M}`$ is indeed small, ruling out the scenario of a massless up-quark. As a further illustration we have included in Fig. 2 the curve which one would expect if $`m_\mathrm{u}=0`$, by taking the central values for $`\alpha _5`$ and $`\alpha _8`$ from eqs. (2) and (8).

The first priority for future work is undoubtedly the application of the method to simulations employing $`N_\mathrm{f}=3`$ flavours of dynamical quarks, and the extension of the quark mass range towards the chiral regime. While efficient simulations of QCD with odd $`N_\mathrm{f}`$ and light dynamical quarks represent an algorithmic challenge, some efforts in this direction have already been made . It would also be interesting to extend applications to the case of flavour singlets, which allow a determination of $`\alpha _7`$ , i.e. another low-energy constant afflicted with the KM ambiguity. Methods to improve the notoriously bad signal/noise ratio in flavour-singlet correlators have been developed , so that there are good prospects for a successful implementation.

## Acknowledgements

We acknowledge the support of the Particle Physics & Astronomy Research Council under grants GR/L22744 and PPA/G/O/1998/00777. We are grateful to the staff of the Edinburgh Parallel Computing Centre for maintaining service on the Cray T3E.


format: text

text:
# The onset of magnetic order in fcc-Fe films on Cu(100)

## Abstract

On the basis of a first-principles electronic structure theory of finite temperature metallic magnetism in layered materials, we investigate the onset of magnetic order in thin (2-8 layers) fcc-Fe films on Cu(100) substrates. The nature of this ordering is altered when the systems are capped with copper. Indeed we find an oscillatory dependence of the Curie temperatures as a function of Cu-cap thickness, in excellent agreement with experimental data. The thermally induced spin-fluctuations are treated within a mean-field disordered local moment (DLM) picture and give rise to layer-dependent ‘local exchange splittings’ in the electronic structure even in the paramagnetic phase. These features determine the magnetic intra- and interlayer interactions which are strongly influenced by the presence and extent of the Cu cap.

The Disordered Local Moment (DLM) picture , when implemented by the KKR-CPA method , is a first-principles, that is to say material specific and yet parameter-free, mean field theory of metallic magnetism in the limit where good local moments form. It works well for bulk $`Fe`$ and $`Co`$ and their alloys but appears to be inadequate for $`Ni`$ . Under these circumstances it is of fundamental interest to test its applicability by using it in different, but nevertheless similar, physical situations such as in a thin films whose thickness, $`d`$, can act as a tuneable, new thermodynamic variable. With this in mind we deployed the theory for $`Fe`$ on and embedded in $`Cu`$. In this letter we report some of our most significant findings.

For complexity in magnetic systems ultrathin $`Fe`$ films on $`Cu(100)`$ are unsurpassed . How the films grow, their structure and morphology are profoundly intertwined with their magnetic properties and a satisfactory description of these systems has become a benchmark for theories of thin film magnetism. After much extensive experimental work on thermally deposited films it is now clear that below a critical thickness of 10-12 monolayers (ML) the films take on the fcc structure of the substrate while thicker films revert to the bcc structure of bulk iron. In the ultra thin regime when there are fewer than 4-5 ML a ferromagnetic (FM) phase is observed whilst the thicker films of 6-11 ML seem to be antiferromagnetic (AF) with a net moment across the film. Adding covering layers of $`Cu`$ to the films has a marked effect upon their magnetic properties. A single $`Cu`$ monolayer suppresses the magnetic ordering temperatures $`T_c`$ whilst the $`T_c`$’s of 2ML Cu-capped films are partly restored. In fact $`T_c`$ oscillates weakly as further Cu layers are added .

The experimental data on the uncapped films has been interpreted  in terms of models in which the top two surface iron layers are coupled together ferromagnetically in an otherwise AF $`Fe`$ film. This picture has been prompted by a series of ab-initio $`T=0`$K electronic structure total energy calculations . Pajda et al. have calculated effective exchange interactions in a $`Cu`$-capped $`Fe`$ monolayer on $`Cu(100)`$ and obtained an oscillatory variation of $`T_c`$ with cap thickness. Here we describe the full DLM picture of both electronic and magnetic structure pertinent to the paramagnetic states of $`Fe`$-films on $`Cu(100)`$ of up to 8ML thickness at finite temperatures and investigate both the growth of magnetic correlations and onset of magnetic order. Moreover we are able to extract intra- and interlayer ‘exchange’ interactions.

The DLM picture is based on assuming a separation between fast and slow degrees of freedom in the interacting many-electron system. For times, $`\tau `$, long in comparison with an electronic hopping times,$`\mathrm{}/W`$ ($`10^{15}`$ seconds), where $`W`$ represents a relevant band-width, but short when compared with typical spin fluctuation times, the spin orientations of the electrons leaving an atomic site are sufficiently correlated with those arriving that the magnetisation integrated over a unit cell and averaged over $`\tau `$ is non-zero. These are the ‘local moments’ which can change their orientations, described by a set of unit vectors $`\{\widehat{e}_i\}`$ labelled by the site index $`i`$, on a time scale longer than $`\tau `$ while their magnitudes fluctuate rapidly on the time scale $`\tau `$. This timescale demarcation has been invoked in much recent work  on bulk systems but has not been examined and tested in thin film systems. This is the issue we address in this letter.

For an ab-initio implementation of this picture, standard Spin Density Functional theory for studying electrons in spin-polarised metals is adapted to describe the states of the system for each orientational configuration $`\{\widehat{e}_i\}`$. At the heart of the theory is the generalised electronic Grand Potential $`\mathrm{\Omega }\{\widehat{e}_i\}`$ of the system so constrained. A long time average can be replaced by an ensemble average with respect to the Gibbsian measure $`P\{\widehat{e}_i\}=Z^1\mathrm{exp}\beta \mathrm{\Omega }\{\widehat{e}_i\}`$, where the partition function $`Z=_i𝑑\widehat{e}_i\mathrm{exp}\beta \mathrm{\Omega }\{\widehat{e}_i\}`$. The thermodynamic free energy, which accounts for the entropy associated with the orientational fluctuations as well as creation of electron-hole pairs, is given by $`F=k_BT\mathrm{log}Z`$. Evidently $`\mathrm{\Omega }\{\widehat{e}_i\}`$ plays the role of a classical ‘spin’ (local moment) Hamiltonian, albeit a highly complicated one. By choosing a suitable reference single-site ‘spin’ Hamiltonian $`\mathrm{\Omega }_0\{\widehat{e}_i\}=_i\omega _i(\widehat{e}_i)`$ and using the Feynman-Peierls’ inequality  a mean field theory is constructed . This is a ‘first principles’ formulation of the Disordered Local Moment (DLM) picture which can be implemented by an adaptation of SCF-KKR-CPA method ideally suited for calculating the partial averages $`\omega _i(\widehat{e}_i)`$. The approach can be further improved via the construct of a generalised Onsager cavity field .

When the above procedure was used to study bulk magnetic metals, it was found that the local electronic structure can possess a ‘local exchange’ splitting even in the paramagnetic state . This means that an electron spin-polarised parallel to a local moment will have a different density of states to that polarised anti-parallel. When all orientations of the moments are averaged over in the paramagnetic state the electronic structure is inevitably unpolarised but consequences from the presence of local moments can still be identified. (Moreover the magnitude of this splitting is expected to vary sharply as a function of wave-vector and energy. At some points, if the ‘bands’ are flat in a region of wave-vector and energy space, the local exchange splitting will be roughly of the same size as the rigid exchange splitting of the majority and minority-spin bands of the ferromagnetic state whereas at other points, where the ‘bands’ have greater dispersion, the splitting will vanish entirely.) This local exchange-splitting is the cause of the establishment of a local moment. Photoemission  and inverse photoemission (IPES) experiments  have found these qualitative features in bcc $`Fe`$. In our generalisation of the DLM theory to layered systems we find that the layer dependence of similar features in the electronic structure of paramagnetic $`Fe`$ grown on a $`Cu(100)`$ substrate drives the onset of magnetic order and the form of the magnetic interactions between $`Fe`$ layers.

To probe the content of the theory consider the response of the paramagnetic DLM state to the application of a small, external spin-only magnetic field, $`\{𝐡_𝐢\}`$, varying from site to site, layer to layer in both orientation and magnitude. The induced magnetisation is predominantly caused by the local moments changing their orientations to align with the field causing a change in the single-site probabilities, $`\delta P_i(𝐞_i)`$. In short the site by site paramagnetic spin susceptibility  is the solution of the matrix equation,

$$\chi ^{ij}=\frac{\beta }{3}\mu _i^2\delta _{ij}+\frac{\beta }{3}\underset{k}{}S_{ik}^{(2)}\chi ^{kj}$$
(1)
where $`S_{ik}^{(2)}`$ is the corresponding direct correlation function and $`\mu _i`$ is the magnitude of the local magnetic moment on the $`i`$-th site. Expressions for $`S_{ij}^{(2)}`$ involving the electronic structure of the paramagnetic state and techniques for calculating them, using the KKR-CPA, for bulk systems have been given elsewhere . For layered systems with two-dimensional translational symmetry the magnitudes of the ‘local moments’, $`\{\mu _i\}`$, assume one value $`\mu _P`$ per site in a given layer labelled $`P`$ but vary from layer to layer. Taking a 2D lattice Fourier transform over sites within each layer equation (1) can be rewritten as,

$`\chi ^{PQ}(𝐪_{})`$ $`=`$ $`{\displaystyle \underset{ij}{}}\chi ^{PiQj}\mathrm{exp}[i𝐪_{}(𝐑_i𝐑_j)]`$ (2)
$`=`$ $`{\displaystyle \frac{\beta }{3}}\mu _P^2\delta _{PQ}+{\displaystyle \frac{\beta }{3}}{\displaystyle \underset{S}{}}S_{PS}^{(2)}(𝐪_{})\chi ^{SQ}(𝐪_{})`$ (3)
where, $`𝐪_{}`$ is a wave-vector in the 2D layer Brilliouin Zone. Once the $`S_{PS}^{(2)}(𝐪_{})`$’s have been obtained from the SCF-KKR-CPA calculations and loaded into a $`n\times n`$ matrix, $`S^{(2)}(𝐪_{})`$ (where $`n`$ is the number of layers in the film), this set of equations (3) is solved by a simple matrix inversion i.e.$`\chi ^{PQ}(𝐪_{})=\left[3k_BTIS^{(2)}(𝐪_{})\right]_{PQ}^1\mu _Q^2`$ ($`I`$ is a unit matrix). For films in which the intralayer exchange coupling is ferromagnetic, as for the $`Fe`$ films on $`Cu`$, the Curie temperature $`T_c`$ is specified by the condition $`3k_BT_CIS^{(2)}(𝐪_{}=0)=0`$ or in terms of the largest eigenvalue of $`S^{(2)}(𝐪_{}=0)`$. Full technical details on how $`S^{(2)}(𝐪_{})`$ is calculated for layered systems will be provided elsewhere.

We use the spin-polarised screened KKR method  for layered systems adapted for the DLM picture. The lattice constants of the layers are taken to be the same as that of the substrate (6.83 a.u.) so that the effects of lattice strain on the electronic structure are neglected. For each $`n`$, the electronic structure of the DLM state is calculated self-consistently using 78 $`k_{}`$ points in the irreducible part of the surface Brillouin zone. A buffer of three layers of the substrate as well as a buffer of at least three layers of vacuum is calculated self-consistently along with the potentials on each layer. The self-consistent layer-resolved potentials are then used to calculate the magnetic moments and the effective ’Weiss’ field on each layer as described above. By introducing a small change in the average magnetisation on one particular layer and calculating the effective ’Weiss’ field on each layer we determine $`S_{PQ}^{(2)}(𝐪_{}=0)`$, and hence the static paramagnetic spin susceptibility, $`\chi ^{PQ}(𝐪_{}=0)`$ together with the Curie temperatures $`T_c`$’s for various thicknesses of both the $`Fe`$ film and the $`Cu`$ cap.

Table I shows our results for the layer-dependent ‘local’ magnetic moments in uncapped Fen/Cu(100) films together with the dependence of the Curie temperatures on Fe-film thickness. The moments on the topmost layer are always the largest at around 2.5$`\mu _B`$ followed by the moments on the Fe layer closest to the substrate next at roughly 2.2$`\mu _B`$. In the interior of the film the moments are reduced to 1.7$`\mu _B`$ close to what was found for fcc-Fe for this lattice spacing . The magnetic ordering temperatures reduce monotonically to a constant value of 485 K for 6 MLs onwards. Table II shows the effective ‘exchange’ interactions, that is to say the direct correlation function,$`S_{PQ}^{(2)}`$ which lead to these transition temperatures. Examples of a 7ML and a 3ML Fe layer systems are shown to bring out the dominant features. The intralayer couplings are ferromagnetic throughout. The layer at the top of film and the layer nearest the Cu substrate have the largest values owing to their reduced effective coordination. It is the interlayer couplings that show an important and interesting trend. The top two layers are strongly ferromagnetically (FM) coupled (positive values) whereas the coupling to nearest neighbor layers within the films is anti-ferromagnetic (AF). There is also significant coupling to subsequent layers which alternates between FM and AF.

Adding a copper cap to these films alters these results profoundly. Figure 1 highlights this for 3 and 7 monolayer Fe/Cu(100) films and shows the magnetic ordering temperature $`T_c`$ to have an oscillatory dependence on cap thickness. A single Cu monolayer suppresses $`T_c`$ by some 50K with the deficit being restored by a second layer. A third layer reduces $`T_c`$ once again and further layers cause minor oscillations about this lower value. This behavior is in excellent agreement with experimental data on the same systems by Vollmer et al. (see their figure 6). In figure 2 the dominant changes to intra- and interlayer magnetic interactions are shown. These occur in the top two Fe layers nearest the cap. A single Cu ML cap switches the coupling between these layers from FM to AF whilst further Cu layers strengthen the coupling within the topmost Fe layer.

Calculations for bulk fcc-Fe have shown how the tendency to magnetic order changes as the lattice spacing is varied . On an expanded lattice, a ferromagnetically ordered state is stable at low temperatures and the coupling between ‘local moments’ in the paramagnetic DLM state mirrors this aspect. As the lattice is contracted the magnetic ordering tendencies become anti-ferromagnetic. Capping fcc Fe films with Cu causes a similar effect. Figure 3 shows that the electronic density of states (DOS) of an uncapped film is ‘exchange split’ over a large energy range. These states are nearly fully occupied and FM coupling between neighboring sites results. The reduced coordination of sites within this layer produce an effect akin to that of increasing the lattice spacing of bulk fcc Fe. There is hybridisation of these Fe states with quantum well states in the Cu cap formed between the vacuum and the Fe film. This causes the exchange splitting to collapse in the region of this hybridisation just as in fcc Fe for smaller lattice spacings. The remaining exchange-split states are roughly half-filled which leads an AF coupling along the (100) direction as in bulk fcc Fe. Figure 3 shows these effects in the DOS in the top layer of a Fe7/Cu(100) system both uncapped and capped by a single Cu ML. The narrow band of states in the Cu layer are also shown.

In conclusion we note that the dramatic variations of interaction between local moments at different distances from the various surfaces and interfaces are quintessentially itinerant effects. As the much studied variation of the local moments near such planar defects show, they are due to the motion of the electrons not limited to an atom but confined by the geometry of the sample. In this respect they are manifestations of the same physics as is at work in the oscillatory magnetic coupling in metallic multilayers . Consequently, the detailed quantitative account of the complex experimental data by our calculations can be taken as significant new evidence that the first-principles (KKR-CPA) implementation of the DLM picture correctly captures the essential physics of magnetic order, due to mobile electrons, in the limit where local moments form.

Acknowledgments: This work is a result of a collaboration under the TMR-Network (contract no. ERBFMRXCT 96-0089) on “Ab initio calculations of magnetic properties of surfaces, interfaces, and multilayers” and also the RTN ’Computational Magnetoelectronics’ network (contract number HPRN-CT-2000-00143) and is supported by the Engineering and Physical Sciences Research Council (UK) and the Hungarian National Science Foundation (contracts OTKA T030240 and T029813).


text:
# 𝐽/Ψ suppression in an expanding hadron gas Talk based on the work done in collaboration with L.Turko [1]. Presented at the XLI Cracow School of Theoretical Physics ”Fundamental Interactions”, Zakopane, Poland, June 2-11 2001.

## I Introduction

On the 10th of February 2000 CERN announced the discovery of the novel state of matter — the so-called quark-gluon plasma (QGP). The existence of the QGP has been predicted upon lattice QCD calculations (for a review see  and references therein) and the critical temperature $`T_c`$ for the ordinary hadron matter-QGP phase transition has been estimated in the range of 150 - 270 MeV (this corresponds to the broad range of the critical energy density $`ϵ_c0.265.5`$ GeV/fm-3). The upper limit belongs to a pure $`SU(3)`$ theory, whereas adding quarks causes lowering of $`T_c`$ even to about 150 MeV ($`ϵ_c0.26`$ GeV/fm-3 respectively). Since NA50 Collaboration estimates for the energy density obtained in the central rapidity region (CRR) give the value of 3.5 GeV/fm-3 for the most central Pb-Pb data point , it is argued that the region of the existence of the QGP has been reached in Pb-Pb collisions at CERN SPS.

Soon after the announcement, the paper entitled ”Evidence for deconfinement of quarks and gluons from the $`J/\mathrm{\Psi }`$ suppression pattern measured in Pb-Pb collisions at the CERN-SPS” appeared . As the title suggests, the main argument for the QGP creation during Pb-Pb collisions at the CERN SPS is the observation of the suppression of $`J/\mathrm{\Psi }`$ relative yield. $`J/\mathrm{\Psi }`$ suppression as a signal for the QGP formation was originally proposed by Matsui and Satz . The clue point of the NA50 Collaboration paper is the figure (denoted as Fig.6 there), where experimental data for Pb-Pb collision values of $`\frac{B_{\mu \mu }\sigma _{J/\psi }}{\sigma _{DY}}`$ (the ratio of the $`J/\mathrm{\Psi }`$ to the Drell-Yan production cross-section times the branching ratio of the $`J/\mathrm{\Psi }`$ into a muon pair) are presented together with some conventional predictions. Here, ”conventional” means that $`J/\mathrm{\Psi }`$ suppression is due to $`J/\mathrm{\Psi }`$ absorption in ordinary hadron matter. Since all those conventional curves saturate at high transverse energy $`E_T`$, but the experimental data fall from $`E_T`$ 90 GeV much lower and this behaviour could be reproduced on the base of $`J/\mathrm{\Psi }`$ disintegration in QGP , it is argued that QGP had to appear during most central Pb-Pb collisions. In general, the immediate reservation about such reasoning is that besides those already known (see e.g.  and references \[12-15\] in ), there could be a huge number of different conventional models, so until this subject is cleared up completely, no one is legitimate to ruled out the absorption picture of $`J/\mathrm{\Psi }`$ suppression.

In the following talk, the more systematic and general description of $`J/\mathrm{\Psi }`$ absorption in the framework of statistical analysis will be presented. The main features of the model are :

* 1. a multi-component non-interacting hadron gas appears in the CRR instead of the QGP;
* 2. the gas expands longitudinally and transversely;
* 3. $`J/\mathrm{\Psi }`$ suppression is the result of inelastic scattering on constituents of the gas and on nucleons of colliding ions.

## II The timetable of events in the CRR

For a given A-B collision $`t=0`$ is fixed at the moment of the maximal overlap of the nuclei (for more details see e.g. ). As nuclei pass each other charmonium states are produced as the result of gluon fusion. After half of the time the nuclei need to cross each other ($`t`$ 0.5 fm), matter appears in the CRR. It is assumed that the matter thermalizes almost immediately and the moment of thermalization, $`t_0`$, is estimated at about 1 fm . Then the matter begins its expansion and cooling and after reaching the freeze-out temperature, $`T_{f.o.}`$, it ceases as a thermodynamical system. The moment when the temperature has decreased to $`T_{f.o.}`$ is denoted as $`t_{f.o.}`$. Since the matter under consideration is the hadron gas, any phase transition does not take place during cooling here.

For the description of the evolution of the matter, relativistic hydrodynamic is explored. The longitudinal component of the solution of hydrodynamic equations (the exact analytic solution for an (1+1)-dimensional case) reads (for details see e.g. )

$$s(\tau )=\frac{s_0\tau _0}{\tau },n_B(\tau )=\frac{n_B^0\tau _0}{\tau },v_z=\frac{z}{t}$$
(1)
where $`\tau =\sqrt{t^2z^2}`$ is a local proper time, $`v_z`$ is the $`z`$-component of fluid velocity ($`z`$ is a collision axis) and $`s_0`$ and $`n_B^0`$ are initial densities of the entropy and the baryon number respectively. For $`n_B=0`$ and the uniform initial temperature distribution with the sharp edge at the border established by nuclei radii, the full solution of (3+1)-dimensional hydrodynamic equations is known . The evolution derived is the decomposition of the longitudinal expansion inside a slice bordered by the front of the rarefaction wave and the transverse expansion which is superimposed outside of the wave. Because small but nevertheless non-zero baryon number densities are considered here, the above-mentioned description of the evolution has to be treated as an assumption in the presented model. The rarefaction wave moves radially inward with a sound velocity $`c_s`$. The sound velocity squared is given by $`c_s^2=\frac{P}{ϵ}`$ and can be evaluated numerically .

## III  The multi-component hadron gas

For an ideal hadron gas in thermal and chemical equilibrium, which consists of $`l`$ species of particles, energy density $`ϵ`$, baryon number density $`n_B`$, strangeness density $`n_S`$ and entropy density $`s`$ read ($`\mathrm{}=c=1`$ always)

$$ϵ=\frac{1}{2\pi ^2}\underset{i=1}{\overset{l}{}}(2s_i+1)_0^{\mathrm{}}𝑑p\frac{p^2E_i}{\mathrm{exp}\left\{\frac{E_i\mu _i}{T}\right\}+g_i},$$
(3)
$$n_B=\frac{1}{2\pi ^2}\underset{i=1}{\overset{l}{}}(2s_i+1)_0^{\mathrm{}}𝑑p\frac{p^2B_i}{\mathrm{exp}\left\{\frac{E_i\mu _i}{T}\right\}+g_i},$$
(4)
$$n_S=\frac{1}{2\pi ^2}\underset{i=1}{\overset{l}{}}(2s_i+1)_0^{\mathrm{}}𝑑p\frac{p^2S_i}{\mathrm{exp}\left\{\frac{E_i\mu _i}{T}\right\}+g_i},$$
(5)
$$s=\frac{1}{6\pi ^2T^2}\underset{i=1}{\overset{l}{}}(2s_i+1)_0^{\mathrm{}}𝑑p\frac{p^4}{E_i}\frac{(E_i\mu _i)\mathrm{exp}\left\{\frac{E_i\mu _i}{T}\right\}}{\left(\mathrm{exp}\left\{\frac{E_i\mu _i}{T}\right\}+g_i\right)^2},$$
(6)
where $`E_i=(m_i^2+p^2)^{1/2}`$ and $`m_i`$, $`B_i`$, $`S_i`$, $`\mu _i`$, $`s_i`$ and $`g_i`$ are the mass, baryon number, strangeness, chemical potential, spin and a statistical factor of specie $`i`$ respectively (an antiparticle is treated as a different specie). And $`\mu _i=B_i\mu _B+S_i\mu _S`$, where $`\mu _B`$ and $`\mu _S`$ are overall baryon number and strangeness chemical potentials respectively.

To obtain the time dependence of temperature and baryon number and strangeness chemical potentials one has to solve numerically equations (4 \- 6) with $`s`$, $`n_B`$ and $`n_S`$ given as time dependent quantities. For $`s(\tau )`$ and $`n_B(\tau )`$ expressions (1) are taken and $`n_S=0`$ since the overall strangeness equals zero during all the evolution (for more details see ).

## IV $`J/\mathrm{\Psi }`$ absorption in the expanding hadronic gas

As it has been already mentioned in Sect. II charmonium states are produced in the beginning of the collision, when nuclei overlap. For simplicity, it is assumed that the production of $`c\overline{c}`$ states takes place at $`t=0`$. To describe $`J/\mathrm{\Psi }`$ absorption quantitatively, the idea of  is generalized to the case of the multi-component massive gas, here. Since in the CRR longitudinal momenta of particles are much lower than transverse ones (in the c.m.s. frame of nuclei), $`J/\mathrm{\Psi }`$ longitudinal momentum is put at zero. Additionally, only the plane $`z=0`$ is under consideration, here. For the simplicity of the model, it is assumed that all charmonium states are completely formed and can be absorbed by the constituents of a surrounding medium from the moment of creation. The absorption is the result of a $`c\overline{c}`$ state inelastic scattering on the constituents of the hadron gas through interactions of the type

$$c\overline{c}+hD+\overline{D}+X,$$
(7)
where $`h`$ denotes a hadron, $`D`$ is a charm meson and $`X`$ means a particle which is necessary to conserve the charge, baryon number or strangeness.

Since the most crucial for the problem of QGP existence are Pb-Pb collisions (see remarks in Sect. I), the further considerations are done for this case. So, for the Pb-Pb collision at impact parameter $`b`$, the situation in the plane $`z=0`$ is presented in Fig. 1, where $`S_{eff}`$ means the area of the overlap of the colliding nuclei.

Additionally, it is assumed that the hadron gas, which appears in the space between the nuclei after they have crossed each other, also has the shape of $`S_{eff}`$ at $`t_0`$ in the plane $`z=0`$. Then, the transverse expansion starts in the form of the rarefaction wave moving inward $`S_{eff}`$ at $`t_0`$.

From the considerations based on the relativistic kinetic equation (for details see ), the survival fraction of $`J/\mathrm{\Psi }`$ in the hadron gas as a function of the initial energy density $`ϵ_0`$ in the CRR is obtained:

$$𝒩_{h.g.}(ϵ_0)=𝑑p_Tg(p_T,ϵ_0)\mathrm{exp}\left\{_{t_0}^{t_{final}}𝑑t\underset{i=1}{\overset{l}{}}\frac{d^3\stackrel{}{q}}{(2\pi )^3}f_i(\stackrel{}{q},t)\sigma _iv_{rel,i}\frac{p_\nu q_i^\nu }{EE_i^{}}\right\},$$
(8)
where the sum is over all taken species of scatters (hadrons), $`p^\nu =(E,\stackrel{}{p}_T)`$ and $`q_i^\nu =(E_i^{},\stackrel{}{q})`$ are four momenta of $`J/\mathrm{\Psi }`$ and hadron specie $`i`$ respectively, $`\sigma _i`$ states for the absorption cross-section of $`J/\mathrm{\Psi }h_i`$ scattering, $`v_{rel,i}`$ is the relative velocity of $`h_i`$ hadron with respect to $`J/\mathrm{\Psi }`$ and $`M`$ and $`m_i`$ denote $`J/\mathrm{\Psi }`$ and $`h_i`$ masses respectively ($`M=3097`$ MeV). The function $`g(p_T,ϵ_0)`$ is the $`J/\mathrm{\Psi }`$ initial momentum distribution. It has a gaussian form and reflects gluon multiple elastic scatterings on nucleons before their fusion into a $`J/\mathrm{\Psi }`$ in the first stage of the collision . The upper limit of the integration over time in (8), namely $`t_{final}`$ is the minimal value of $`t_{esc}`$ and $`t_{f.o.}`$. The quantity $`t_{esc}`$ is the average time of the escape of $`J/\mathrm{\Psi }`$’s from the hadron medium for given values of $`b`$ and $`J/\mathrm{\Psi }`$ velocity $`\stackrel{}{v}=\stackrel{}{p}_T/E`$. Note that the average is taken with the weight

$$p_{J/\mathrm{\Psi }}(\stackrel{}{r})=\frac{T_A(\stackrel{}{r})T_B(\stackrel{}{r}\stackrel{}{b})}{T_{AB}(b)},$$
(9)
where $`T_A(\stackrel{}{s})=𝑑z\rho _A(\stackrel{}{s},z)`$ and $`\rho _A(\stackrel{}{s},z)`$ is the nuclear matter density distribution. For the last quantity, the Woods-Saxon nuclear matter density distribution with parameters from  is taken. In the integration over hadron momentum in (8) the threshold for the reaction (7) is included, i.e. $`\sigma _i`$ equals zero for $`(p^\nu +q_i^\nu )^2<(2m_D+m_X)^2`$ and is constant elsewhere ($`m_D`$ is a charm meson mass, $`m_D=1867`$ MeV). Also the usual Bose-Einstein or Fermi-Dirac distribution for hadron specie $`i`$ is used in (8)

$$f_i(\stackrel{}{q},t)=f_i(q,t)=\frac{2s_i+1}{\mathrm{exp}\left\{\frac{E_i^{}\mu _i(t)}{T(t)}\right\}+g_i}.$$
(10)
As far as $`\sigma _i`$ is concerned, there are no data for every particular $`J/\mathrm{\Psi }h_i`$ scattering. Therefore, only two types of the cross-section, the first, $`\sigma _b`$, for $`J/\mathrm{\Psi }baryon`$ scattering and the second, $`\sigma _m`$, for $`J/\mathrm{\Psi }meson`$ scattering are assumed, here. In the model presented $`\sigma _b`$ is put at $`\sigma _{J/\psi N}`$ — $`J/\mathrm{\Psi }Nucleon`$ absorption cross-section. And values of $`\sigma _{J/\psi N}`$ in the range of 3-5 mb are taken from p-A data . The $`J/\mathrm{\Psi }meson`$ cross-section $`\sigma _m`$ is assumed to be 2/3 of $`\sigma _b`$, which is due to the quark counting.

As it has been already suggested  also $`J/\mathrm{\Psi }`$ scattering in nuclear matter should be included in any $`J/\mathrm{\Psi }`$ absorption model. This could be done with the introduction of $`J/\mathrm{\Psi }`$ survival factor in nuclear matter

$$𝒩_{n.m.}(ϵ_0)\mathrm{exp}\left\{\sigma _{J/\psi N}\rho _0L\right\},$$
(11)
where $`\rho _0`$ is the nuclear matter density and $`L`$ the mean path length of the $`J/\mathrm{\Psi }`$ through the colliding nuclei. The length $`L`$ is expressed by the following formula :

$$\rho _0L(b)=\frac{1}{2T_{AB}}d^2\stackrel{}{s}T_A(\stackrel{}{s})T_B(\stackrel{}{s}\stackrel{}{b})\left[\frac{A1}{A}T_A(\stackrel{}{s})+\frac{B1}{B}T_B(\stackrel{}{s}\stackrel{}{b})\right],$$
(12)
where $`T_{AB}(b)=d^2\stackrel{}{s}T_A(\stackrel{}{s})T_B(\stackrel{}{s}\stackrel{}{b})`$.

Since $`J/\mathrm{\Psi }`$ absorptions: in nuclear matter and in the hadron gas, are separate in time, $`J/\mathrm{\Psi }`$ survival factor for a heavy-ion collision with the initial energy density $`ϵ_0`$, could be defined as

$$𝒩(ϵ_0)=𝒩_{n.m.}(ϵ_0)𝒩_{h.g.}(ϵ_0).$$
(13)
Note that since right sides of (8) and (11) include parts which depend on impact parameter $`b`$ and the left sides are functions of $`ϵ_0`$ only, the expression converting the first quantity to the second (or reverse) should be given. This is done with the use of the dependence of $`ϵ_0`$ on the transverse energy $`E_T`$ extracted from NA50 data  (for details see ).

To make the model as much realistic as possible, one should keep in mind that only about $`60\%`$ of $`J/\mathrm{\Psi }`$’s measured are directly produced during collision. The rest is the result of $`\chi `$ ($`30\%`$) and $`\psi ^{}`$ ($`10\%`$) decay . Therefore the realistic $`J/\mathrm{\Psi }`$ survival factor could be expressed as

$$𝒩(ϵ_0)=0.6𝒩_{J/\psi }(ϵ_0)+0.3𝒩_\chi (ϵ_0)+0.1𝒩_\psi ^{}(ϵ_0),$$
(14)
where $`𝒩_{J/\psi }(ϵ_0)`$, $`𝒩_\chi (ϵ_0)`$ and $`𝒩_\psi ^{}(ϵ_0)`$ are given also by formulae (8), (11) and (13) but with $`g(p_T,ϵ_0)`$, $`\sigma _{J/\psi N}`$ and $`M`$ changed appropriately (for details see ).

To complete, also values of cross-sections for $`\chi baryon`$ and $`\psi ^{}baryon`$ scatterings are needed. For simplicity, it is assumed that both these cross-sections are equal to $`J/\mathrm{\Psi }baryon`$ one. Since $`J/\mathrm{\Psi }`$ is smaller than $`\chi `$ or $`\psi ^{}`$, this means that $`J/\mathrm{\Psi }`$ suppression evaluated according to (14) is underestimated here.

## V  Results

To calculate formula (8) initial values $`s_0`$ and $`n_B^0`$ are needed. For the estimation of the last quantity experimental results for S-S  and Au-Au  collisions are explored. The S-S data give $`n_B^00.25`$ fm-3 and Au-Au data $`n_B^00.65`$ fm-3. These values are for central collisions, for more peripheral ones the initial baryon number density should be lower. To simplify numerical calculations, $`n_B^0`$ is kept constant over the whole range of $`b`$. But, to check the possible dependence on $`n_B^0`$, also the much lower value $`n_B^0=0.05`$ fm-3 is examined. To find $`s_0`$, (3 \- 5) with respect to $`T`$, $`\mu _S`$ and $`\mu _B`$ are solved. Then, having put $`T_0`$, $`\mu _S^0`$ and $`\mu _B^0`$ into (6), $`s_0`$ is obtained. The last parameter of the model is the freeze-out temperature. Two values $`T_{f.o.}=100,140`$ MeV are taken here and they agree well with estimates based on hadron yields .

Firstly, formula (14) with the use of (8) only instead of (13) is calculated for the case of (1+1)-dimensional expansion and results are presented in Figs. 3-3.

For comparison, also the experimental data are shown in Figs. 3-3. The experimental survival factor is defined as

$$𝒩_{exp}=\frac{\frac{B_{\mu \mu }\sigma _{J/\psi }^{AB}}{\sigma _{DY}^{AB}}}{\frac{B_{\mu \mu }\sigma _{J/\psi }^{pp}}{\sigma _{DY}^{pp}}},$$
(15)
where $`\frac{B_{\mu \mu }\sigma _{J/\psi }^{AB(pp)}}{\sigma _{DY}^{AB(pp)}}`$ is the ratio of the $`J/\mathrm{\Psi }`$ to the Drell-Yan production cross-section in A-B(p-p) interactions times the branching ratio of the $`J/\mathrm{\Psi }`$ into a muon pair. The values of the ratio for p-p, S-U and Pb-Pb are taken from . The results of numerical evaluations of (14) agree with the data well qualitatively for greater baryonic cross-section $`\sigma _b`$ and (or) for the lower freeze-out temperature $`T_{f.o.}`$, except the last point measured for the highest $`E_T`$ bin.

Now the full (3+1)-dimensional hydrodynamic evolution (in the form described in Sect. II) and $`J/\mathrm{\Psi }`$ absorption in nuclear matter will be taken into account and results are depicted in Figs. 5-5. To draw also S-U data together with Pb-Pb ones, instead of multiplying $`𝒩_{h.g.}`$ by $`𝒩_{n.m.}`$, $`𝒩_{exp}`$ is divided by $`𝒩_{n.m.}`$, i.e. ”the experimental $`J/\mathrm{\Psi }`$ hadron gas survival factor” is defined as

$$\stackrel{~}{𝒩}_{exp}=\mathrm{exp}\left\{\sigma _{J/\psi N}\rho _0L\right\}𝒩_{exp},$$
(16)
and values of this factor are drawn in Figs. 5-5 as the experimental data. The quantity $`r_0`$ is a parameter from the expression for a nucleus radius $`R_A=r_0A^{\frac{1}{3}}`$. The value of $`R_A`$ is used to fix the initial position of the rarefaction wave in the evaluation of $`t_{esc}`$ in (8). It has turned out that the case of $`n_B^0=0.05`$ fm-3 does not differ substantially from that of $`n_B^0=0.25`$ fm-3, so curves for $`n_B^0=0.05`$ fm-3 are not depicted in Figs. 5-5.

Having compared Figs. 5-5 with Figs. 3-3, one can see that adding the transverse expansion changes the final (theoretical) pattern of $`J/\mathrm{\Psi }`$ suppression qualitatively. Now the curves for the case including the transverse expansion are not convex, in opposite to the case with the longitudinal expansion only, where the curves are. As far as the pattern of suppression is concerned, theoretical curves do not fall steep enough at high $`ϵ_0`$ to cover the data area. But for some choice of parameters, namely for $`\sigma _b`$ somewhere between 4 and 5 mb and for $`r_0=1.12`$ fm, a quite satisfactory curve could have been obtained. Precisely, again only the highest $`E_T`$ bin point falls down of the range of theoretical estimates completely. But the error bar of this point is very wide. And also the contradiction in positions of the last three points of the 1998 data can be seen. This means that the high $`E_T`$ region should be measured once more with the better accuracy to state definitely whether the abrupt fall of the experimental survival factor takes place or not there.

To support the conclusion, main results from Figs. 5-5 are presented in Fig. 6 again. The original data  for $`\frac{B_{\mu \mu }\sigma _{J/\psi }^{PbPb}}{\sigma _{DY}^{PbPb}}`$ and $`J/\mathrm{\Psi }`$ survival factors given by (14) multiplied by $`\frac{B_{\mu \mu }\sigma _{J/\psi }^{pp}}{\sigma _{DY}^{pp}}`$ and now as functions of $`E_T`$ are depicted there (for details see ).

As it has been already mentioned, the main disagreement with the data appears in the last experimental point of the 1998 analysis.

The $`charmoniumbaryon`$ inelastic cross-section is the most crucial parameter in the model presented. So to be sure what is the exact result of $`J/\mathrm{\Psi }`$ absorption, one should know how this cross-section behaves in the hot hadron environment. The newest estimations of $`\pi +J/\mathrm{\Psi }`$, $`\rho +J/\mathrm{\Psi }`$ and $`J/\mathrm{\Psi }+N`$ cross-sections at high invariant collision energies  give values of $`\sigma _b`$ and $`\sigma _m`$ of the same order as assumed here. Also it should be stressed that the $`charmoniumhadron`$ inelastic cross-sections are considered as constant quantities in the model presented. But, as the results of just mentioned papers suggest, they are not constant at all. The cross-sections are growing functions of the invariant collision energy $`\sqrt{s}`$. Therefore, one could think naively that the increase of $`ϵ_0`$ (or in other words $`E_T`$) causes the increase of the invariant collision energy $`\sqrt{s}`$ on the average and further the increase of the $`charmoniumhadron`$ inelastic cross-sections. In this way, the line describing $`J/\mathrm{\Psi }`$ survival factor could be placed close to the solid curve of Fig. 6 for low $`ϵ_0`$ ($`E_T`$), and then, as the $`charmoniumhadron`$ inelastic cross-sections increase, this line would go closer to the dashed curve of Fig. 6 for high $`ϵ_0`$ ($`E_T`$). So, the experimental pattern of $`J/\mathrm{\Psi }`$ suppression could be recovered.

As the last remark, the author would like to call reader’s attention to Fig.5 of . The solid and dotted curves in that figure are exactly the same as curves just presented in Fig. 6. But the distortion in the QGP is the main reason for $`J/\mathrm{\Psi }`$ suppression in ! It is also stated there, that results shown ”provide evidence for the production of the quark-gluon plasma in central high-energy Pb-Pb collisions”.

Therefore, the final conclusion is the following: the NA50 PB-PB data provides evidence for the production of the thermalized and confined hadron gas in the central rapidity region of a PB-PB collision. This implies that $`J/\mathrm{\Psi }`$ suppression is not a good signal for the quark-gluon plasma appearance at a central heavy-ion collision.


text:
# Ge-substitutional defects and the √3×√3 ↔ 3×3 transition in 𝛼–Sn/Ge(111)

\[

## Abstract

The structure and energetics of Ge substitutional defects on the $`\alpha `$-Sn/Ge(111) surface are analyzed using Density Functional Theory (DFT) molecular dynamics (MD) simulations. An isolated Ge defect induces a very local distortion of the 3$`\times `$3 reconstruction, confined to a significant downwards displacement (-0.31 Å) at the defect site and a modest upward displacement (0.05 Å) of the three Sn nearest neighbours with partially occupied dangling bonds. Dynamical fluctuations between the two degenerate ground states yield the six-fold symmetry observed around a defect in the experiments at room temperature. Defect-defect interactions are controlled by the energetics of the deformation of the 3$`\times `$3 structure: They are negligible for defects on the honeycomb lattice and quite large for a third defect on the hexagonal lattice, explaining the low temperature defect ordering.

\]

The driving force behind the T-induced reversible $`\sqrt{3}\times \sqrt{3}3\times 3`$ phase transition observed for 1/3-ML coverage of Sn or Pb on Ge(111) has been under intense debate since its discovery in Pb/Ge(111). After the early suggestions of surface Fermi wavevector nesting or electron correlations, leading to the formation of a Charge Density Wave  at low T, recent theoretical work indicates that the softening of a 3$`\times `$3 surface phonon, associated with the electronic energy gain due to the Sn dangling bond (DB) rehybridization, plays a major role on this transition. The dynamical fluctuations model associated with this soft mode provides a consistent account of the main experimental facts, in particular the similarities and differences between the isovalent Sn/Ge and Sn/Si(111) systems.

DFT calculations for the $`\alpha `$–Sn/Ge(111) surface predict that the ground state is a 3$`\times `$3 reconstruction in which 1/3 of the Sn adatoms are placed $``$ 0.3 Å higher than the other Sn adatoms . In other words, out of the three 3$`\times `$3 sublattices that form the $`\sqrt{3}\times \sqrt{3}`$ lattice , “up” Sn adatoms ($`Sn_u`$) occupy one of these sublattices (say $`a`$), forming a 3$`\times `$3 hexagonal lattice, while sublattices $`b,c`$ are occupied by ‘down’ Sn adatoms ($`Sn_d`$) (forming a honeycomb pattern). These calculations also show that this reconstruction is correlated with a change in the surface electronic structure: the $`Sn_u`$–DB is now fully occupied while the other two $`Sn_d`$–DBs share the remaining electron per 3$`\times `$3 unit cell. STM experiments are naturally explained in terms of these different DB occupancies. At LT, below $`T_c`$, filled-state STM images show the 3$`\times `$3 phase : bright spots occupy the 3$`\times `$3 sublattice associated with the $`Sn_u`$–DBs, while the other two sublattices, associated with the $`Sn_d`$–DBs, are occupied by darker spots. At RT,however, dynamical fluctuations of the Sn adatom heights, fast with respect to the STM scan speed, induce a $`\sqrt{3}\times \sqrt{3}`$ periodicity in the STM images.

While these models assume an ideal 1/3-ML Sn surface coverage, point defects, with an average concentration of 3$`\pm `$1 % , are always present. 90 % of these defects are Ge substitutionals (called Ge defects in the rest of the paper), while vacancies account for the other 10 %. STM, a real-space experimental technique, has been the method of choice to get non-averaged information on the structure and distribution of these defects on Sn/Ge(111) across a wide temperature range, covering the transition temperature $`T_c`$. Three are the main experimental results: (i) At RT Ge defects are imaged as dark spots in the $`\sqrt{3}\times \sqrt{3}`$ lattice surrounded by six features (brighter than Sn adatoms away from defects) that correspond to the six nearest-neighbours (n.n.) Sn atoms. (ii) At LT most of the Ge defects occupy positions corresponding to the two honeycomb $`Sn_d`$–sublattices $`b,c`$, with equal probability, and are therefore surrounded by three bright and three darker Sn n.n. (iii) At RT Ge defects are randomly distributed over the three sublattices $`a,b,c`$.

These results have been considered as supporting evidence for an interpretation of the $`\sqrt{3}\times \sqrt{3}`$ 3$`\times `$3 phase transition in terms of the interaction of the Ge defects. One of the main ingredients in this argument is the six-fold symmetry assumed for the perturbation induced by an isolated defect. This six-fold symmetry (honeycomb, with one atom dark and two atoms bright in the three possible sublattices), suggested by the RT STM images, is, however, in contradiction with the hexagonal pattern (two atoms dark and one atom bright, with three-fold symmetry) observed in the filled-state images below $`T_c`$. In order to explain the LT experimental results, a defect-defect interaction, with a temperature dependent decay length $`l(T)`$ (e.g. $`l(300K)`$ 11 Å, $`l(100K)`$ 100 Å) has to be introduced. This long–range interaction would be responsible for the ordering of defects at LT on two of the three possible sublattices. Finally, the superposition of the perturbations induced by defects in these two sublattices would give rise to the 3$`\times `$3 phase.

In this letter, we show that ab initio DFT-LDA MD simulations of Ge substitutional defects on the Sn/Ge(111) surface support a completely different point of view: Both the structure and the ordering of the Ge defects at LT are induced by the ground state 3$`\times `$3 reconstruction. In particular, we have studied the atomic geometry, the electronic structure and the energetics of one and two Ge defects in a large 9$`\times `$9 unit cell. In the case of two defects, we have considered three cases, corresponding to the defects on first, second and third n.n. positions. Our analysis shows that Ge defects at LT always create a three-fold distortion around them, the simulated STM filled-states image showing three bright and three darker spots on the n.n. Sn atoms. The distortion is very local, affecting only to the Sn atoms which are nearest neighbours of the defects. We have also performed first-principles MD simulations of a 3$`\times `$3 unit cell with one Ge defect. Our MD simulations indicate that the six-fold pattern observed around Ge defects in the RT STM images should be understood in terms of the RT dynamical fluctuations of the six n.n. Sn adatoms. These dynamical effects contribute, together with structural and electronic effects, to the enhanced brightness of the Sn atoms around the defect.

In our approach we calculate different defects using a Sn/Ge(111)-9$`\times `$9 unit cell with a slab of 4 Ge layers and a layer of H atoms saturating the dangling bonds of the deeper Ge layer (432 atoms in the unit cell). Both the H layer and the last Ge layer are kept fixed during the calculations to mimic the bulk continuation. In our calculations we have used an efficient MD technique, the Fireball96 code , that allows us to analyze the large 9$`\times `$9 unit cell we have considered. We have checked that this code yields accurate surface geometries by comparing it with a plane wave code (CASTEP) for the case of the Sn/Ge(111)-3$`\times `$3 reconstruction and for the case of one Ge defect on a 3$`\times `$3 surface unit cell (see Table I). These results confirm the quality of the calculations provided by the Fireball96 code and justify the wide use we are going to make of it in the rest of the paper. STM images are calculated from the ab-initio electronic structure of the surface using the Keldysh-Green function formalism and considering a W tip (for details on the method and the application to the Sn/Ge(111)-3$`\times `$3 system see ref. ).

Figure 1a shows the relaxed geometry obtained in our calculations for the $`\alpha `$–Sn/Ge(111) surface with one Ge defect per 9$`\times `$9 unit cell (3.7 % concentration). The surface presents the 3$`\times `$3 reconstruction with the Ge adatom on one of the two 3$`\times `$3–sublattices (say $`b`$) that form the honeycomb pattern associated with $`Sn_d`$ adatoms. Our simulations yield for the Ge substitutional a downwards displacement of 0.30 Å with respect to the original $`Sn_d`$ position. The perturbation induced by the Ge defect is well-localized, affecting only to the three n.n. $`Sn_d`$ adatoms ($`c`$ sublattice). These adatoms are displaced upwards $``$ 0.05 Å due to the transfer of charge from the Ge defect: the Ge–DB is empty and each n.n. $`Sn_d`$–DB gains $``$ 1/6 electrons. The occupancy of the n.n. $`Sn_u`$–DBs ($`a`$ sublattice) does not change (they are already fully occupied), and therefore their position is practically the same as in the 3$`\times `$3 reconstruction free of defects. Figure 1b shows the calculated filled-state STM image for this geometry. Probably the most remarkable feature is the three-fold symmetry of the image around the defect (only three of the surrounding Sn atoms are imaged as protrusions), in contrast with the six-fold symmetry assumed for the perturbation induced by the defect from the experimental RT STM images.

Notice that the solution shown in fig. 1 is two-fold degenerate: an equivalent solution is obtained applying a symmetry operation that exchanges the two sublattices $`(a,c)`$ occupied by the up and down Sn. This degeneracy will be important when discussing the effect of temperature on the surface geometry and its corresponding STM image.

Figures 2a,b,c display the relaxed adatom positions we have obtained for the $`\alpha `$–Sn/Ge(111)–9$`\times `$9 surface with two defects on first, second and third n.n. positions. In all the cases the 3$`\times `$3 reconstruction adapts its particular geometry (out of a three-fold degenerate case) to the defect positions, in such a way that the two defects are placed on the two $`Sn_d`$–sublattices forming a honeycomb pattern. When the Ge defects are second n.n. (fig. 2b) they are on the same sublattice (say $`b`$), and the ground state is two-fold degenerate: an exchange of the $`a`$ and $`c`$ sublattices gives a completely equivalent solution. In fig. 2a and 2c the Ge defects occupy different sublattices and the solution is non-degenerate. As in the case of an isolated defect, the perturbation induced by the defects is localized essentially on the Sn-adatoms n.n. to the defects. The total energy for the three cases shown in fig. 2 is, within the accuracy of these calculations ($`20`$ meV/9$`\times `$9 unit cell), the same, confirming the localized nature of the distortion induced by the Ge defects. It is remarkable that the sum of the upward displacements of all the neighbouring Sn-adatoms is, in the three cases, very similar to the sum of the displacements induced by two independent defects ($``$ 0.30 Å, see fig. 1a). This result suggests a roughly linear relationship between the upwards displacement of a Sn adatom and the transfer of charge to its DB.

Coming back to the experimental STM data for the Ge defects on the Sn/Ge(111) surface, we first notice that our solution for a single defect (fig, 1b) has a three-fold and not a six-fold symmetry, with the Sn atoms around the defect presenting two different types of displacements perpendicular to the surface.This apparent contradiction with the experimental data at RT can be resolved, however, considering the dynamics of the atoms around the defect. To this end, we have considered a 3$`\times `$3 unit cell with a single substitutional and performed MD simulations at two different temperatures, T= 180 K and T= 270 K. The initial conditions correspond to the positions of the relaxed structure of the defect (see Table I) and random velocities following a Maxwell-Boltzmann distribution. Then we follow the time evolution of the system integrating the equations of motion with the forces obtained from the DFT-LDA calculations. Figure 3 shows the time evolution of the $`z`$-coordinate of the three upper atoms (two Sn atoms and the Ge susbstitutional). These results show that Ge remains around its equilibrium site, while the two Sn atoms tend to exchange heights, very much in similarity with what we found for the $`\sqrt{3}\times \sqrt{3}`$ structure . Therefore, the RT STM images of the six Sn–adatoms n.n. to a Ge defect must represent a 50 % average of $`Sn_u`$ and $`Sn_d`$ (displaced upwards $``$ 0.05 Å w.r.t. normal $`Sn_d`$–adatoms, and containing 1/6 electrons more), while the RT STM images of Sn–adatoms away from defects represent a 1/3 $`Sn_u`$ \+ 2/3 $`Sn_d`$ average. This dynamical (plus charge transfer) effect explains in a simple way the different pattern observed at RT and LT on the STM images around a Ge defect (six-fold vs. three-fold symmetry) and also explains why Sn–adatoms n.n. to a Ge defect appear on the RT STM images brighter than Sn–adatoms in regions free of defects.

We finally address the following question: why do Ge defects at LT occupy preferentially two 3$`\times `$3–sublattices, with equal probability, and not the three of them or only one?. Our results show that, on the $`\alpha `$–Sn/Ge(111) surface, Ge substitutional defects transfer their DB charge to neighbouring Sn-adatoms DBs and do not participate in the up/down dynamical fluctuations of the Sn–adatoms at RT. Our results also show that at LT, when these fluctuations are supressed, the stable configuration of the surface is a 3$`\times `$3 reconstruction in which Ge defects occupy $`Sn_d`$ positions; the defect-defect interaction is negligible (as long as they stay on $`Sn_d`$ sublattices) and Ge defects occupy with equal probability any of the two $`Sn_d`$ sublattices.

A third Ge defect initially located in the other sublattice (hexagonal) should be forced, when the temperature is lowered, either to present an upward displacement (contrary to what it does when it is isolated) or to deform dramatically its environment in order to move downwards. Any of these solutions would imply a significant energy cost (160 meV for the defect in an upper position according to our calculations), that the system prefers not to pay, locking in all the substitutionals in a honeycomb lattice. We conclude that the ordering of Ge defects at LT is induced by the ground state 3$`\times `$3 reconstruction and not the other way around: Ge substitutional defects are not the driving force of the $`\sqrt{3}\times \sqrt{3}`$ $``$ 3$`\times `$3 transition in the $`\alpha `$–Sn/Ge(111) surface.

In conclusion, our DFT-LDA total-energy and MD calculations show that the ground state geometry around a Ge substitutional defect corresponds to a very local distortion of the 3$`\times `$3 reconstruction of the defect-free surface, confined to the defect site and the Sn nearest neighbours. The lattice deformation can be characterized by a downwards displacement of the Ge defect and upward displacement of the neighbouring Sn atoms with partially occupied DBs, with respect to the ideal 3$`\times `$3 structure. This structural effect is correlated with the charge transfer from the defect DB to the DBs of these Sn atoms. This ground state configuration provides, at LT, an STM image with a three-fold symmetry, while at high T, the dynamical fluctuations between the two degenerate states (with different Sn atoms in the upper positions) would yield the six-fold symmetry observed in the experiments. These dynamical effects, together with the upward displacement and the increased DB occupancy, are responsible of the enhanced brightness of the Sn atoms around the defect. Finally, the ordering of Ge defects at LT is naturally explained as an effect of the 3$`\times `$3 reconstruction.

This work has been partly funded by the spanish CICYT under contract No.PB-97-0028. Part of the calculations have been performed in the CCC-UAM.


text:
# Evidence for a Nodeless Gap from the Superfluid Density of Optimally Doped Pr1.855Ce0.145CuO4-y Films

## Abstract

We present measurements of the ab-plane magnetic penetration depth, $`\lambda (T)`$, in five optimally doped Pr1.855Ce0.145CuO4-y films for $`1.6`$ K $`TT_c24`$ K. Low resistivities, high superfluid densities $`n_s(T)\lambda ^2(T)`$, high $`T_c`$’s, and small transition widths are reproducible and indicative of excellent film quality. For all five films, $`\lambda ^2(T)/\lambda ^2(0)`$ at low $`T`$ is well fitted by an exponential temperature dependence with a gap, $`\mathrm{\Delta }_{min}`$, of 0.85 $`k_BT_c`$. This behavior is consistent with a nodeless gap and is incompatible with d-wave superconductivity.

It is widely accepted that pairing symmetry in the hole-doped cuprates is predominantly $`d_{x^2y^2}`$, at least near optimal doping where most phase-sensitive measurements have been made . Essentially all other experimental results agree with this view. The developing situation in the electron-doped cuprates is not so clear. While recent phase-sensitive measurements on optimally doped Nd1.85Ce0.15CuO4-y (NCCO) and Pr1.85Ce0.15CuO4-y (PCCO) films are consistent with a $`d_{x^2y^2}`$ energy gap , there is substantial evidence for a nodeless gap. The penetration depth, $`\lambda (T)`$, measured by Alff et al.  via field modulation of Josephson junctions fabricated in the ab-plane of optimally doped PCCO films is exponentially flat at low $`T`$, suggesting a gapped superconducting state. The zero-bias peak seen in tunnelling data from hole-doped cuprates – associated with Andreev bound states of a d-wave order parameter – is absent from tunnel junctions in NCCO , suggesting s-wave superconductivity. In light of persuasive evidence for $`d_{x^2y^2}`$ superconductivity from phase-sensitive  and angle resolved photoemission spectroscopy  measurements on electron-doped samples, one might suspect some problem with measurements on in-plane tunnel junctions, or with the films themselves. This suspicion may be strengthened by reports of quadratic behavior of $`\lambda (T)\lambda (0)`$ at low $`T`$ in optimally doped PCCO crystals . Thus the importance of penetration depth measurements that confirm the findings of Alff et al. and demonstrate high sample quality and sample-to-sample reproducibility.

We present measurements of the ab-plane superfluid density, $`n_s(T)\lambda ^2(T)`$, in five optimally doped Pr1.855Ce0.145CuO4-y films. $`\lambda ^2(T)`$ at low $`T`$ in hole-doped cuprates is consistently linear  or quadratic  in temperature. Theory has not found a scenario in which behavior flatter than $`T^2`$ is predicted for $`d_{x^2y^2}`$ superconductors . On the other hand, $`\lambda ^2(T)`$ in gapped – e.g., s-wave – superconductors is exponentially flat at low temperatures. As shown below, $`\lambda ^2(T)`$ in our optimally doped PCCO films is exponential in $`\mathrm{\Delta }_{min}/T`$, with the *same* value of the minimum gap, $`\mathrm{\Delta }_{min}=0.85\pm 0.05`$ $`k_BT_c`$.

Our films were prepared by molecular-beam epitaxy on 12.7 mm $`\times `$ 12.7 mm $`\times `$ 0.35 mm SrTiO3 substrates as detailed elsewhere . Ce content, $`x`$, is measured using inductively coupled plasma spectroscopy and is known to better than $`\pm 0.005`$. Films are highly oriented with their c-axes perpendicular to the substrate. Table I summarizes film properties. Thicknesses vary from 750 Å to 1250 Å, and one film (P6) was grown on a 250 Å buffer layer of insulating Pr2CuO4. ab-plane resistivities, $`\rho (T)`$ (Fig. 1), are reproducible to $`\pm 15\%`$. $`T_c=23.7\pm 0.5`$ K and transition width, $`\mathrm{\Delta }T_c=0.9\pm 0.2`$ K meet or exceed values reported for nominally identical crystals . As noted previously , residual resistivities of PCCO films are smaller than those of PCCO crystals . Curvature in $`\lambda (T)\lambda (0)`$ varies by more than a factor of two from crystal to crystal . These comparisons suggest that film quality is superior to crystal quality at this time. $`\rho (T)`$ in our films at 40 K is a factor of two lower than that of high-quality, comparably doped crystals of La2-xSrxCuO4 (LSCO) , the hole-doped cousin of PCCO. $`\lambda (0)=1800\pm 300`$ Å in our films matches the lowest value, 1930 Å, reported in LSCO near optimal doping , and it is only slightly larger than the a-axis penetration depth, 1600 Å, in YBCO . We therefore conclude that our films have very little disorder.

We measure $`\lambda ^2(T)`$ with a low frequency two-coil mutual inductance technique described in detail elsewhere . A film is centered between two small coils, and a current at about 50 kHz in one coil induces eddy currents in the film. Currents are approximately uniform through the film thickness. Data have been measured to be independent of frequency for $`10`$ kHz $`f100`$ kHz. Magnetic fields from the primary coil and the film are measured as a voltage across the secondary coil. We have checked that the typical excitation field (100 $`\mu `$Tesla $``$ to film) is too small to create vortices in the film. Because the coils are much smaller than the film, the applied field is concentrated near the film’s center and demagnetizing effects at the film perimeter are irrelevant. All data presented here are in the linear response regime.

Because film thicknesses, $`d`$, are less than $`\lambda `$, it is the sheet conductivity $`\sigma (\omega ,T)d=\sigma _1(\omega ,T)di\sigma _2(\omega ,T)d`$ that is measured, with an estimated accuracy of 5%. $`\lambda ^2(T)`$ is obtained from $`\sigma _2`$: $`\lambda ^2(T)\mu _0\omega \sigma _2(T)`$, where $`\mu _0`$ is the magnetic permeability of vacuum, and accuracy is limited by 5% uncertainty in $`d`$. The temperature dependence of $`\lambda ^2(T)/\lambda ^2(0)`$ is unaffected by uncertainty in $`d`$. The origin of 0.2% “wiggles” in $`\lambda ^2(T)`$ at low temperatures is uncertain, but we know that at least some of the effect arises from slow drift in amplifier gain.

$`T_c`$ and $`\mathrm{\Delta }T_c`$ (Table I) are defined to be the position and full-width of the fluctuation peak in $`\sigma _1(T)`$ measured at 50 kHz (Fig. 2). $`T_c`$ from resistivity (Fig. 1) and penetration depth measurements (Figs. 2 and 3) are identical. Structure in $`\sigma _1(T)`$ is due to slight variation in $`T_c`$ through the film thickness. Our measurement technique reveals transitions in all layers of the sample, so $`\mathrm{\Delta }T_c1`$ K indicates excellent film homogeneity.

Figure 3 displays $`\lambda ^2(T)`$ in all five films (thick solid lines). The spread in $`\lambda ^2(0)`$ is larger than expected from uncertainty in $`d`$ . Slight upward curvature in $`\lambda ^2`$ near $`T_c`$ is reproducible. Most importantly, the flatness of the low-temperature data is highly reproducible.

To determine the behavior of $`\lambda ^2(T)`$ at low temperatures, we fit the first $`5`$% drop in $`\lambda ^2(T)/\lambda ^2(0)`$ to

$$\lambda ^2(T)/\lambda ^2(0)1C_{\mathrm{}}e^{D/t},$$
(1)
with $`tT/T_c`$. $`D`$ is roughly the minimum gap value, $`\mathrm{\Delta }_{min}`$, normalized to $`T_c`$ and is fixed at 0.85 for all films to emphasize sample-to-sample reproducibility. $`C_{\mathrm{}}`$ and $`\lambda ^2(0)`$ (Table I) are adjusted for each film to yield the fits (dotted curves) in Fig. 4. The exponential curves in Fig. 4 *are not best fits*, but their $`\chi ^2`$ values (Table I) lie within the estimated experimental noise of $`\pm 0.2`$% (error bars in Fig. 4). The worst fits are improved by a factor of two by adjusting $`D`$ by 10%. We note that a more thorough analysis finds that the peak in the density of states is near 2.2 $`k_BT_c`$ in optimally doped PCCO .

For comparison, quadratic and cubic fits of the form $`\lambda ^2(T)\lambda ^2(0)\left[1\left(T/T_0\right)^x\right]`$ with $`x=2`$ or 3, were also performed. They are shown for film P2 in Fig. 4. For all films, the best quadratic fits are visibly poorer than cubic or exponential fits. $`\chi ^2`$ values for quadratic fits are typically twice as large as for exponential fits. We note that cubic fits possess the smallest $`\chi ^2`$ values, leaving open the possibility that the superconducting density of states increases as energy cubed at low energies rather than being fully gapped . No d-wave theory predicts such behavior.

Figure 3 displays $`\lambda ^2(T)`$ from 0 K to $`T_c`$ along with low-$`T`$ exponential fits from Fig. 4 and low-$`T`$ quadratic and cubic fits extrapolated to higher temperatures. Exponential fits are cut off at 8 K, since the temperature dependence of the gap is neglected in the low-$`T`$ expression, Eq. (1). We have already demonstrated that low-$`T`$ quadratic fits are unacceptably poor. Even if one believes that the low-temperature data are quadratic in $`T`$, there is no sign of a crossover to linear behavior , and the extrapolated low-$`T`$ quadratics go to zero well above $`T_c`$, so the curvature is too weak to be compatible with a dirty d-wave interpretation . We note that previous observations of quadratic behavior in $`\lambda (T)`$  do not find a crossover from quadratic to linear, weakening their support for d-wave superconductivity.

In conclusion, we have presented high-precision measurements of $`\lambda ^2(T)`$ in five optimally doped PCCO films. Resistivities, $`T_c`$’s, transition widths, and superfluid densities are reproducible and indicative of superior film quality. In every film, $`\lambda ^2(T)`$ at low temperatures is flatter than $`T^2`$ and is well fitted by an exponential temperature dependence with the same minimum gap value of 0.85 $`k_BT_c`$, indicating a nodeless gap in optimally doped PCCO.