text:
# Effective potentials for atom-atom interaction at low temperatures

\]http://bgaowww.physics.utoledo.edu

(August 8, 2002)

## Abstract

We discuss the concept and design of effective atom-atom potentials that accurately describe any physical processes involving only states around the threshold. The existence of such potentials gives hope to a quantitative, and systematic, understanding of quantum few-atom and quantum many-atom systems at relatively low temperatures.

The concept of model potential has played an important role in many branches of physics. The well known examples include the Morse potential Morse (1929) and the Lennard-Jones potential for molecular systems, and the hard-sphere potential and the delta-function pseudopotential Huang and Yang (1957) widely used in many-body theories, including theories for Bose-Einstein condensates (BEC) Leggett (2001).

There are many reasons why one uses a model potential instead of the “real” potential. Here, we only emphasize that a problem can simply become unmanageable if the “real” potential is used. For quantum systems, this statement quickly becomes true for three or more atoms. This explains, in part, that despite its many limitations, it has proven difficult to go substantially beyond the Gross-Pitaveskii theory for BEC that is based on the delta-function pseudopotential Leggett (2001).

Our goal here is to discuss the concept and design of model potentials that better reflect the reality of atom-atom interaction than either the hard-sphere potential or the delta-function pseudopotential, yet simple enough to allow for applications in quantum few-atom and quantum many-atom systems. Here, being “simple”, to a large extent, means being shallow, as it is the depth of a potential, which can be measured by the number of bound states it supports, that determines the complexity of the resulting quantum few-atom and quantum many-atom problems.

One of our key conclusions is the following. For N-atom states around the N-atom threshold (such as the BEC state Leggett (2001)), or physical processes that involve only states around the threshold (such as the three-body recombination process Nielsen et al. (2001); Suno et al. (2002)), the interaction potential between a pair of atoms, no matter how deep it might be, can be replaced by an effective potential supporting only one or a few bound states. Furthermore, because different partial waves are described by the same effective potential, it can be used in precisely the same manner as any “real” potential. And in doing so, one reduces the complexity of the resulting quantum few-atom and quantum many-atom problems to a level comparable to that for He, a level that we are quickly learning to deal with Blume and Greene (2000); Suno et al. (2002).

If one thinks of the delta-function pseudopotential Huang and Yang (1957) as describing the atomic interaction at the longest length scale in the zero-energy limit, $`2\pi /k`$, the natural next step is the description of atomic interaction at the next, shorter, length scale. This scale is $`\beta _n=(2\mu C_n/\mathrm{}^2)^{1/(n2)}`$, which characterizes the long-range atomic interaction of the form of $`C_n/r^n`$ ($`n>2`$). The angular-momentum-insensitive quantum-defect theory (AQDT) Gao (2001, 2000) provides a systematic description of atom-atom interaction at this scale, and is the basis of our discussion.

Reference Gao (2001) focused on two-atom systems with $`V(r)C_6/r^6`$, but the same concepts and formulation are readily generalized to any $`n>2`$. In this general formulation, a two-atom system with an asymptotic potential of the form of $`C_n/r^n`$ ($`n>2`$) is described by a dimensionless K matrix $`K^c(ϵ,l)`$ and a set of universal functions that are determined from the solutions of

$$\left[\frac{d^2}{dr^2}\frac{l(l+1)}{r^2}+\frac{\beta _n^{n2}}{r^n}+\overline{ϵ}\right]u_{ϵl}(r)=0,$$
(1)
where $`\overline{ϵ}2\mu ϵ/\mathrm{}^2`$. Specifically, $`K^c(ϵ,l)`$ is defined by writing the wave function at large distances as a linear superposition of a pair of reference solutions of Eq. (1):

$$u_{ϵl}(r)=A_{ϵl}[f_{ϵl}^c(r)K^c(ϵ,l)g_{ϵl}^c(r)],$$
(2)
where $`f^c`$ and $`g^c`$ are purposely chosen to have the behavior

$`f_{ϵl}^c`$ $`\stackrel{r\beta _n}{}`$ $`(2/\pi )^{1/2}(r/\beta _n)^{n/4}\mathrm{cos}\left(y\pi /4\right),`$ (3)
$`g_{ϵl}^c`$ $`\stackrel{r\beta _n}{}`$ $`(2/\pi )^{1/2}(r/\beta _n)^{n/4}\mathrm{sin}\left(y\pi /4\right),`$ (4)
for all energies Gao (2001). Here $`y=[2/(n2)](\beta _n/r)^{(n2)/2}`$.

AQDT asserts that $`K^c(ϵ,l)`$ is approximately a constant that is independent of both $`ϵ`$ and $`l`$, provided that $`\beta _n`$ is greater than other, energy-independent, length scales present in the system Gao (2001). For our purposes here, the most important conclusion of AQDT is the following. To the extent of $`K^c(ϵ,l)`$ being energy and angular-momentum independent, potentials with the same type of long-range behavior (namely the same $`n`$) and the same $`K^c`$ have, on a scaled energy basis, the same bound spectra and scattering properties around the threshold Gao (2001). Here the bound spectra and scattering properties include all angular momentum states for which $`K^c`$ remains approximately $`l`$-independent.

Ignoring scaling relations implied by this statement, which will discussed elsewhere Gao , AQDT gives the following simple prescription for designing an effective potential that has the same physical properties around the threshold as the system of interest. First, choose a model potential. The only restriction is that it should have the right asymptotic behavior. Second, adjust the short range parameters of the model potential so that

$$K_{\mathrm{mod}}^c(0,l)=K^c(0,l),$$
(5)
for one particular $`l`$. These two conditions do not uniquely determine an effective potential. Another auxiliary condition, which gives a convenient characterization of the depth of a potential, is the number of bound levels supported by a model potential for a particular $`l`$, $`N_l`$. For classes of model potentials discussed below, these conditions uniquely determine a model potential.

We stress here two classes of model potentials for which $`K^c(0,l)`$, and the number of bound levels for each $`l`$, $`N_l`$, can be found analytically. One class is of the type of a hard-sphere with an attractive tail:

$$V_{HST}(r)=\{\begin{array}{ccc}\mathrm{}\hfill & ,& rr_0\hfill \\ C_n/r^n\hfill & ,& r>r_0\hfill \end{array},$$
(6)
which will be denoted by HST. The other class is of the type of Lennard-Jones LJ$`(n,2n2)`$:

$$V_{LJn}(r)=C_n/r^n+C_{2n2}/r^{2n2},$$
(7)
which will be denoted by LJn. In particular, this potential corresponds to a LJ(6,10) potential for $`n=6`$.

For HST potentials, it is not difficult to show that the $`K^c`$ parameter at zero energy is given by Gao

$$K_{\text{HST}}^c(0,l)=\frac{J_{\nu _0}(y_0)\mathrm{cos}(\pi \nu _0/2)Y_{\nu _0}(y_0)\mathrm{sin}(\pi \nu _0/2)}{J_{\nu _0}(y_0)\mathrm{sin}(\pi \nu _0/2)+Y_{\nu _0}(y_0)\mathrm{cos}(\pi \nu _0/2)},$$
(8)
where $`\nu _0=(2l+1)/(n2)`$, $`J`$ and $`Y`$ are the Bessel functions Abramowitz and Stegun (1964), and $`y_0=[2/(n2)](\beta _n/r_0)^{(n2)/2}`$. The number of bound levels for angular momentum $`l`$ is given by

$$N_{\text{HST}}(l)=\{\begin{array}{ccc}m\hfill & ,& j_{\nu _0,m}y_0<j_{\nu _0,m+1}\hfill \\ 0\hfill & ,& y_0<j_{\nu _0,1}\hfill \end{array},$$
(9)
where $`j_{\nu _0,m}`$ ($`m1`$), is the m-th zero of the Bessel function $`J_{\nu _0}(x)`$ Abramowitz and Stegun (1964).

For LJn potentials, the following results can be derived with the help of a local scaling transformation Gao .

$$K_{\text{LJn}}^c(0,l)=\mathrm{tan}(\pi \nu _0/2)[1+h_l(z_0)][1h_l(z_0)]^1,$$
(10)
where $`z_0=(\beta _n/\beta _{2n2})^{n2}/[2(n2)]`$,

$$h_l(z_0)=z_0^{\nu _0}\frac{\mathrm{sin}\pi (z_0+1/2\nu _0/2)\mathrm{\Gamma }(z_0+1/2\nu _0/2)}{\mathrm{sin}\pi (z_0+1/2+\nu _0/2)\mathrm{\Gamma }(z_0+1/2+\nu _0/2)},$$
(11)
and $`\nu _0=(2l+1)/(n2)`$. The number of bound levels for any $`l`$ is given by

$$N_{\text{LJn}}(l)=\{\begin{array}{ccc}\left[z_0+\frac{1}{2}\frac{\nu _0}{2}\right]\hfill & ,& z_0(\nu _0+1)/2\hfill \\ 0\hfill & ,& z_0<(\nu _0+1)/2\hfill \end{array},$$
(12)
where $`[x]`$ means the greatest integer less than or equal to $`x`$. We note that scattering lengths, which are much more restrictive concepts Levy and Keller (1963); Gao (2000), can be derived from $`K^c(0,l)`$. For example, the $`s`$ wave scattering length, which is defined only for $`n>3`$, can be obtained from $`K^c(0,l=0)`$:

$$a_{l=0}/\beta _n=\left[b^{2b}\frac{\mathrm{\Gamma }(1b)}{\mathrm{\Gamma }(1+b)}\right]\frac{K^c(0,0)+\mathrm{tan}(\pi b/2)}{K^c(0,0)\mathrm{tan}(\pi b/2)},$$
(13)
where $`b=1/(n2)`$.

From these results, the HST and LJn types of potentials can be readily designed according to Eq. (5) to have the desired $`K^c(0,l)`$, and the desired $`N_l`$. Table 1 gives a selected set of designs for the triplet state of a 23Na dimer. Here the potentials are designed to have $`K^c(0,l=0)=13.57`$, which is found numerically using the latest potential for a sodium dimer Laue et al. (2002). From Eq. (13), this corresponds to an $`s`$ wave scattering length of 64.57 a.u. The number of bound $`s`$ wave levels supported by this “real” potential is found numerically to be 16.

Figure 1 shows the comparison of $`s`$ and $`d`$ wave partial cross sections of effective LJ(6,10) and effective HST potentials, both designed to support 16 $`s`$ wave bound levels and have a $`K^c(0,l=0)=13.57`$, with the Na-Na partial cross sections computed from the “real” potential Laue et al. (2002). The results are hardly distinguishable over a wide range of energies. In comparison, the hard-sphere potential (HS) fails quickly away from the threshold for the $`s`$ wave, and gives completely wrong results for the $`d`$ wave.

This result confirms the concept of effective potential based on AQDT. For our purposes here, however, what is more important is how shallow the effective potentials can be while still maintain a good description of low-energy characteristics of a real system. Figure 2 shows the comparison of $`s`$ and $`d`$ wave partial cross sections of an effective LJ(6,10) potential, in this case designed to support only a single $`s`$ wave bound state, with the “real” Na-Na results Na . The agreements remain excellent.

Figure 3 shows similar results for a HST potential. In this case, the HST potential that supports only a single $`s`$ wave bound state does not do as well near the $`d`$ wave shape resonance. But it quickly improves, monotonically, as the number of bound levels it is designed to support increases. By $`N_{l=0}=4`$, a good agreement is achieved.

The robustness of these designs is not limited to the description of scattering properties, it also applies to the energies of bound states that are close to the dissociation threshold, and to the wave functions. For example, for effective potentials supporting four $`s`$ wave bound levels, the HST gives a binding energy of 0.2027 GHz for the least-bound $`s`$ state, the LJ(6,10) gives 0.2003 GHz. Both in good agreement with the result for the “real” potential, which gives 0.2044 GHz Na . Figure 4 shows that for effective potentials supporting four $`s`$ wave bound levels, the wave functions are well represented down to $`r=20`$ a.u., covering basically all regions of space in which there is an appreciable amplitude.

We stress that while only the results for sodium are presented here. They are used to illustrate much more general concepts. As the number of bound levels supported by an effective potential increases, all physical properties of states around the threshold converge to the same results (see Figs. 1-4). This is the shape-independence at length scale $`\beta _n`$. The converged results, properly scaled, represent a set of universal properties shared by all quantum systems with the same type of long-range potential, and characterized by the same, $`l`$-independent constant $`K^c=lim_{\beta _x/\beta _n0}K^c(0,l)`$ Gao (2001), where $`\beta _x`$ represent the next shorter length scale present in the system. The examples presented here, Figs. 2-4, show how quickly this set of universal properties are approached as one increases the number of bound levels supported by an effective potential. This quick convergence is due to the fact that deviations from the universal behavior depend on a high power of $`\beta _x/\beta _n`$ Gao . Other quantum systems differ from Na primarily in $`K^c`$ Gao (2001), which does not effect this rapid convergence.

Note that we did not make any distinction between two-atom and N-atom quantum systems in the statements above, because the same apply to a N-atom system. A short argument is simply that diffuse states, in which atoms are mostly at large distances relative to each other, only couple coherently to other diffuse states. A longer argument can proceed as follows. The correct $`K^c(0,0)`$, and therefore $`a_{l=0}`$, ensures the correct results at the mean-field level Leggett (2001). The correct two-atom wave function ensures the correct two-atom correlation. It also ensures the correct three-atom correlation, as follows. Think of a three-atom as a two-atom perturbed by another. Frank-Condon considerations tell one that only two-atom states around the threshold are significantly coupled. This means by having the correct two-atom wave functions around the threshold, one has also the correct three-atom correlation around the threshold. …

To illustrate the savings of computer resources as a result of using an effective potential, consider the problem of N interacting atoms in a symmetric trap of frequency $`\omega `$. If a real potential is used, the fact we need to represent the length scale of $`(2\mu D_e/\mathrm{}^2)^{1/2}`$ means we need roughly $`D_e/\mathrm{}\omega `$ number of harmonic oscillator states for each atom, for a total of $`(D_e/\mathrm{}\omega )^N`$ number of states (ignoring statistics). If an effective potential is used, the corresponding number is $`(D_{e,\text{eff}}/\mathrm{}\omega )^N`$. Thus the saving in the size of basis set is characterized by the factor $`(D_{e,\text{eff}}/D_e)^N`$. For the triplet state of Na, $`D_{e,\text{eff}}/D_e`$ is of the order of $`10^3`$ if an effective potential with two $`s`$ wave bound levels is used (see Table 1). This corresponds to a saving in the size of basis set of $`10^9`$ fold just for a three-atom problem. Even greater savings are achieved for deeper potentials or for more atoms. From another angle, for effective potentials with $`N_{l=0}1`$, all length scales shorter than $`\beta _n`$ have effectively been eliminated. It is this elimination of short length scales that makes a complex problem more manageable.

On the other hand, if a good description over an even wider range of energies around the threshold is desired, the same methodology can be carried to scales shorter than $`\beta _n`$ (e.g., $`\beta _8`$ for atoms with a $`C_6/r^6`$ long-range interaction and a $`C_8/r^8`$ correction) Gao . However, because the ratio $`\beta _8/\beta _6`$ is different for different systems, the results become dependent upon one more system-specific parameter. At this stage, going to shorter length scale seems useful only in specialized two-atom applications Gao . We also point out that it is around the threshold that the quantum effects are most important Gao (1999).

In conclusion, we have established the concept and design of effective potentials describing atomic interaction at the length scale of $`\beta _n`$. It is the scale that one has to deal with in studying quantum few-atom Nielsen et al. (2001) and quantum many-atom systems Leggett (2001) at finite temperature, of high density, or under strong confinement. We expect the method presented to play a role in our understanding of some of the more complex systems and processes at low temperatures, such as the three-body recombination process Nielsen et al. (2001); Suno et al. (2002), excited clusters states Blume and Greene (2000), and quantum liquids Leggett (2001). In all cases, one can look for, and verify universal properties at the scale of $`\beta _n`$, by comparing results from different designs, such as HST and LJn, and by checking convergences as one relaxes an effective potential towards more bound levels.

###### Acknowledgements.

I thank Eite Tiesinga, Mike Cavagnero, and Brett Esry for helpful discussions. Special thanks goes to Eite for providing the potential for sodium dimer. This work was supported in part by the National Natural Science Foundation of China (No. 19834060) and the Key Project of Knowledge Innovation Program of Chinese Academy of Science (No. KJCX2-W7), and in part by the US National Science Foundation.


format: text

text:
# Quasiparticle bandstructure effects on the d hole lifetimes of copper within the 𝐺𝑊 approximation

## Abstract

We investigate the lifetime of d holes in copper within a first–principle $`GW`$ approximation. At the $`G_0W_0`$ level the lifetime of the topmost d bands are in agreement with the experimental results and are four times smaller than those obtained in the “on–shell” calculations commonly used in literature. The theoretical lifetimes and bandstructure, however, worsen when further iterative steps of self-consistency are included in the calculation, pointing to a delicate interplay between self–consistency and the inclusion of vertex corrections. We show that the $`G_0W_0`$ success in the lifetimes calculation is due to the opening of new “intraband” decay channels that disappear at self–consistency.

Very recent experimental and theoretical results on the quasiparticle lifetimes in noble campillo ferdi  and simple silkin  metals show that our present understanding of the electron dynamics in real solids is far from being complete angelrev . For electrons above the Fermi level (hot electrons) time–resolved two–photon photoemission experiments show a non quadratic behavior of the lifetime, in contrast with the Fermi liquid theory prediction graphite . For occupied states, the number of possible scattering events increases rapidly as their energies decrease below the Fermi level, and the agreement between photoemission experiments and theoretical results worsens angel . An open question is whether this discrepancy comes from effects beyond the approximation used in the calculations (i.e., beyond $`GW`$, where vertex corrections are neglected), or it comes from the way used to solve the quasiparticle equation for a given self–energy. In the present work we address quantitatively the latter point, the inclusion of vertex corrections being beyond the scope of this paper.

Density Functional Theory DFT  (DFT) has become the state–of–the–art approach to study ground state properties of a large number of systems, going from molecules, to surfaces, to complex solids dft-app . The success of DFT is based on the idea that the spatial density of a system of interacting particles can be exactly described by a non–interacting gas of Kohn-Sham (KS) independent particles, moving under the action of an effective potential which includes the exchange-correlation potential $`V_{xc}`$. Compared with experiment, the usual local density (LDA) LDA  approximation to the DFT yields semiconductor bandstructure which systematically underestimate the band gap, while for noble metals the discrepancies are both k–point and band dependent noi . Another important drawback of the use of DFT eigenvalues as excitations (bandstructure) energies is that they are by construction real; no lifetime effects are included. An alternative approach is Time Dependent DFT where all neutral excitations are, in principle, exactly described rmp . However, quasiparticle lifetimes have not been considered so far within this approach.

Many–body perturbation theory allows one to obtain band energies and lifetimes in a rigorous way, i.e. as the poles of the one–particle Green’s function $`G`$ rmp . Those are determined by the solution of a Dyson-like equation of the form hedin :

$$\begin{array}{c}\left[\frac{\mathrm{}^2}{2m}_𝐫^2+V_{external}\left(𝐫\right)+V_{Hartree}\left(𝐫\right)\right]\psi _{n𝐤}(𝐫,\omega )+\hfill \\ \hfill 𝑑𝐫^{}\mathrm{\Sigma }(𝐫𝐫^{},\omega )\psi _{n𝐤}(𝐫^{},\omega )=E_{n𝐤}\left(\omega \right)\psi _{n𝐤}(𝐫,\omega ),\end{array}$$
(1)
containing the non-local, generally complex and non-hermitian, frequency dependent self–energy operator $`\mathrm{\Sigma }`$. The poles of G are the quasiparticle (QP) energies $`ϵ_{n𝐤}^{QP}`$, that from Eq. (1) correspond to the generally complex solutions of the equation: $`ϵ_{n𝐤}^{QP}=E_{n𝐤}\left(ϵ_{n𝐤}^{QP}\right)`$. The real part of $`ϵ_{n𝐤}^{QP}`$ gives the quasiparticle bandstructure; whereas the imaginary part yields the inverse quasiparticle lifetime. In the present work $`\mathrm{\Sigma }`$ is evaluated according to the so-called GW approximation, derived by Hedin in 1965 hedin , which is based on an expansion in terms of the dynamically screened Coulomb interaction $`W\left(\omega \right)`$. In the first iteration $`G_0`$ and $`W_0`$ from DFT calculations are used to compute $`\mathrm{\Sigma }`$ as $`\mathrm{\Sigma }_0=iG_0W_0`$. Unlike semiconductors, the case of noble metals has been studied only recently campillo ; noi ; campillo2 . The lifetime of hot–electrons in copper has been calculated with the $`GW`$ self–energy evaluated at the DFT zero–order energies campillo  (namely “on mass-shell” approximation). This approach is based on the assumption of vanishing QP corrections of the DFT states while, very recently, large QP effects on the occupied bands of copper have been found noi . The “on mass-shell” approach applied to the hole lifetimes campillo2  indicates that d–holes in copper exhibit a longer lifetime than excited s/p electrons. The quantitative comparison with experiment angel , however, has shown a large overestimation of the experimental lifetimes measured by means of photoemission spectroscopy. In this paper we go beyond the “on mass-shell” approximation, and calculate the lifetimes of d bands in copper fully solving the QP Eq. (1) in the complex plane without relying on any analytic continuation. The convergence of the results is carefully checked. Our results are significantly different from those obtained within a DFT self–energy based approach and are in good agreement with experiments. The careful analysis of the physics underlying the GW decay of quasiholes will shed light into the quantum–mechanical mechanism behind the electron dynamics in noble metals.

In our approach, we start by solving Eq. (1) on the real axis, as it is commonly done in QP bandstructure calculations gwrev . In this way we obtain $`ϵ_{n𝐤}^{QP,0}`$, a first guess for the real QP energies:

$`ϵ_{n𝐤}^{QP,0}=ϵ_{n𝐤}^{DFT}+Re\left[\mathrm{\Sigma }_{n𝐤}\left(ϵ_{n𝐤}^{QP,0}\right)\right]V_{xc}^{n𝐤},`$ (2)
where $`\mathrm{\Sigma }_{n𝐤}\left(\omega \right)n𝐤|\mathrm{\Sigma }(𝐫𝐫^{},\omega )|n𝐤`$ and $`V_{xc}^{n𝐤}n𝐤|V_{xc}\left(𝐫\right)|n𝐤`$. Even if Eq. (1) in principle requires a diagonalization with respect to the band index $`n`$, it can be reduced to the form of Eq. (2) because the computed off–diagonal matrix elements of $`\mathrm{\Sigma }`$ are at least two orders of magnitude smaller than the diagonal ones ($`\left|\mathrm{\Sigma }_{n,n^{}}\right|<<\left|\mathrm{\Sigma }_{n,n}\right|`$ $`n,n^{}`$). The difference between the requested exact quasiparticle energy $`ϵ_{n𝐤}^{QP}`$ and the first guess defined in Eq. (2) is due to the fact that the self–energy has an imaginary part, namely:

$`ϵ_{n𝐤}^{QP}ϵ_{n𝐤}^{QP,0}=iZ_{n𝐤}Im\left[\mathrm{\Sigma }_{n𝐤}\left(ϵ_{n𝐤}^{QP,0}\right)\right],`$ (3)
and

$`Z_{n𝐤}\left[1{\displaystyle \frac{d\mathrm{\Sigma }_{n𝐤}\left(\omega \right)}{d\omega }}|_{\omega =ϵ_{n𝐤}^{QP,0}}\right]^1.`$ (4)
Since the quasiparticle concept holds when $`Im\left[\mathrm{\Sigma }_{n𝐤}\left(ϵ_{n𝐤}^{QP,0}\right)\right]`$ is small $`ϵ_{n𝐤}^{QP,0}`$ is a natural starting point to get the QP excitation. Note that being $`Z_{n𝐤}`$ complex Eq. (3) will also slightly modify the real part of $`ϵ_{n𝐤}^{QP,0}`$. Thus the corresponding linewidth ($`\mathrm{\Gamma }_{n𝐤}`$) and lifetime ($`\tau _{n𝐤}`$) are given by

$`\tau _{n𝐤}^12\mathrm{\Gamma }_{n𝐤}2Re\left[Z_{n𝐤}\right]Im\left[\mathrm{\Sigma }_{n𝐤}\left(ϵ_{n𝐤}^{QP,0}\right)\right].`$ (5)
A remarkable property of Eq. (5) is that, being $`ϵ_{n𝐤}^{QP,0}`$ solution of Eq. (3), only the $`Re\left[Z_{n𝐤}\right]`$ is needed to define the QP lifetime in contrast with the common expansion of $`\mathrm{\Sigma }_{n𝐤}\left(\omega \right)`$ around the DFT energy $`ϵ_{n𝐤}^{DFT}`$. $`Z_{n𝐤}`$ is the usual renormalization factor, referred to the initial QP guess instead of to the DFT eigenvalue. A first approximate solution of Eqs. (12) can be obtained by fully neglecting the QP correction $`Re\left[\mathrm{\Sigma }_{n𝐤}\left(ϵ_{n𝐤}^{QP,0}\right)\right]`$, i.e., by assuming $`ϵ_{n𝐤}^{QP,0}ϵ_{n𝐤}^{DFT}`$ and hence $`Z_{n𝐤}=1`$. The corresponding lifetime is $`\tau _{n𝐤}\left\{2Im\left[\mathrm{\Sigma }_{n𝐤}\left(ϵ_{n𝐤}^{DFT}\right)\right]\right\}^1`$ and is usually referred in the literature as the “on mass-shell” $`G_0W_0`$ lifetime quinn ; angelrev , because the input energies of the QP equation are supposed to remain constant, and are subsequently used to calculate the lifetimes. Another approach uses Eq. (3) with real $`Z_{n𝐤}`$, and $`ϵ_{n𝐤}^{QP,0}`$ corresponding to a LMTO–$`GW`$ bandstructure ferdi . This method has confirmed the overestimation of the top d–bands lifetimes found in the “on mass-shell” approximation campillo2 . However, we stress that the approximated real QP energies $`ϵ_{n𝐤}^{QP,0}`$ used in Ref. ferdi  are not solution of Eq. (2). As a consequence, $`\mathrm{\Gamma }_{n𝐤}`$, obtained from Eq. (3), contains and additional term proportional to $`Im\left[Z_{n𝐤}\right]`$ that usually is neglected ferdi  or, it is very small. However for the case of noble metals it can be as large as $`20`$ meV for the top d–bands (see Fig. 1), accounting for 10–40% of the total electronic linewidth. This term reduces the lifetimes, in agreement with the experiment and with the results presented below.

Al the present calculations of Green’s function and screened interaction have been performed using a plane wave basis. We have used norm–conserving pseudopotentials MT ; CUPRB  in the DFT–LDA calculation, with a 60 Ry energy cutoff, corresponding to $`\mathrm{\hspace{0.17em}800}`$ plane waves notaconv . A particular care has been devoted to the check of the effect of the Lorentzian broadening $`\eta `$ included in the screening function $`W\left(\omega \right)`$ for numerical reasons (see Ref. gwrev  for details). Calculations with different broadening values have shown an increase up to $`50\%`$ of the lifetimes when $`\eta `$ is reduced from $`0.1`$ to $`0.005`$ eV. This dramatic numerical effect can be avoided by using $`\eta =0.005`$ eV, yielding results which coincides with those extrapolated at $`\eta =0`$ eV. This is an important point since most calculations presented so far have been done using dampings of 0.1 eV or larger.

In Fig. 1 we present the “on mass-shell” $`G_0W_0`$ result compared with the full solution of Eqs. (23) (simply referred as to $`G_0W_0`$), and with the experimental photoemission results angel . The “on mass-shell” $`G_0W_0`$ yields lifetimes which are four times too large at the d–bands top, and three times too small at the d–bands bottom (not shown in Fig. 1). The $`G_0W_0`$ results, instead, are in rather good agreement with experiments. Actually, they are systematically above the experimental values, as expected for the contribution of higher–order electron–electron and residual phononic and impurity contributions. Moreover, as shown in Tab. 1, the energy positions of the quasiparticle peaks are well reproduced in $`G_0W_0`$, while in the “on mass-shell” calculation the eigenvalues (coincident with DFT ones) span a larger range of energies, reflecting the known DFT overestimation of the linewidth of d–bands noi . The origin of the big difference between the two lifetime calculations stems from the large QP self–energy corrections to the d–bands of copper, completely neglected in “on mass-shell” calculations.

The proper inclusion of these non trivial self–energy corrections change the electronic decay channels of the d levels. In particular, we expect a strong effect on the topmost bands. Let us develop this idea further, by looking at the $`G_0W_0`$ expression for the self–energy:

$`\mathrm{\Sigma }_0(𝐫𝐫^{},\omega )={\displaystyle \frac{id\omega ^{}}{2\pi }G_0(𝐫𝐫^{},\omega ^{})W_0(𝐫𝐫^{},\omega \omega ^{})}.`$ (6)
$`W_0(𝐫,𝐫^{};\omega ^{})`$ is the dynamically screened potential (convolution of the inverse dielectric function $`ϵ^1(𝐫,𝐫^{};\omega ^{})`$ with the bare Coulomb potential). We can write explicitly Eq. (3) in terms of a summation over the DFT states embodied in $`G_0`$ as:

$$\begin{array}{c}\mathrm{\Gamma }_{n𝐤}\underset{n^{}}{}\underset{𝐪}{}Im\left[W_{n𝐤n^{}\left(𝐤𝐪\right)}\left(ϵ_{n^{}\left(𝐤𝐪\right)}^{DFT}ϵ_{n𝐤}^{QP,0}\right)\right]\hfill \\ \hfill \theta \left(ϵ_{n^{}\left(𝐤𝐪\right)}^{DFT}ϵ_{n𝐤}^{QP,0}\right)f_{n^{}\left(𝐤𝐪\right)},\end{array}$$
(7)
where

$$\begin{array}{c}W_{n𝐤n^{}\left(𝐤𝐪\right)}\left(\omega \right)𝑑𝐫𝑑𝐫^{}\psi _{n𝐤}^{}\left(𝐫\right)\psi _{n^{}\left(𝐤𝐪\right)}\left(𝐫\right)\hfill \\ \hfill W(𝐫𝐫^{},\omega )\psi _{n^{}\left(𝐤𝐪\right)}^{}\left(𝐫^{}\right)\psi _{n𝐤}\left(𝐫^{}\right),\end{array}$$
(8)
being $`f_{n𝐤}`$ the Fermi occupations and $`\theta \left(\omega \right)`$ the step function. In Eq. (7) a quasihole in the state $`|n𝐤`$ looses energy via transitions to all the possible occupied states with (higher) energy $`ϵ_{n^{}\left(𝐤𝐪\right)}^{DFT}`$; the energy difference is dissipated by the the screening cloud (described by $`W\left(\omega \right)`$) that surrounds the DFT hole $`|n𝐤`$. Thus the $`G_0W_0`$ transitions contributing to the hole linewidth look like single particle transitions from a quasihole to a DFT hole. In higher order of Hedin’s equations, i.e. including vertex corrections, this interpretation of transitions looses meaning. This is due to self–energy effects on the Green’s function of Eq. (6) and interactions of the screening cloud with the state $`|n^{}\left(𝐤𝐪\right)`$ (vertex corrections hedin ) that in Eq. (7) are neglected.

For the top of the occupied d bands the $`G_0W_0`$ calculation yields negative QP corrections noi ; this means that Eq. (7) contains also contributions coming from the decay of the QP d band to exactly the same DFT band. We will refer to these transitions as “intraband” decay channels. These contributions are important because the d–bands of copper are flat, hence the corresponding density of states is large. Moreover the screened interaction between d–bands is strong, as the d–states are spatially localized and screening is less effective at small distances. In the “on mass-shell” calculation the states $`|n^{}\left(𝐤𝐪\right)`$ appearing in Eq. (7) and the quasiparticle states correspond to the same DFT eigenvalues; this means that the “intraband” decay channels occur at zero energy, where the low–energy Drude tail of the dielectric function dominates (and hence the screened interaction is vanishing). This leads to the usual interpretation of the long calculated lifetimes at the d bands top as due to the fact that these d states can only decay to s/p states. These matrix elements are smaller than those involving states with the same l–character.

Even if the results of the $`G_0W_0`$ calculation are in rather good agreement with the experiments, a natural question about the physical meaning of the “intraband” decay channels arises. Being not at self-consistency, the system is described within $`G_0W_0`$ on the basis of quasiparticle states and DFT states (those involved in the hole decay) with different energies. To remove this inconsistency Eq. (1) should be solved iteratively, until converged QP energies are obtained. However, as we show below, iterations beyond the usual $`G_0W_0`$ level worsen both the imaginary and real parts of the calculated QP energies.

So far fully self-consistent $`GW`$ calculations have been performed only for the homogeneous electron gas holm  and for simple semiconductors and metals eguiluz , yielding worse spectral properties than those obtained in the non self-consistent $`G_0W_0`$. The construction of a self-consistent $`GW`$ self–energy is a formidable task even for the simple systems mentioned above. In copper, already the update of the screening function is rather demanding, due to the presence of localized d orbitals that imply a large cut–off in the plane wave expansion. To test the effect of self-consistency on the QP energies and lifetimes we use a simplified $`G_iW_0`$ method, where the self–energy operator is defined as

$`\mathrm{\Sigma }_i(𝐫𝐫^{},\omega ){\displaystyle \frac{id\omega ^{}}{2\pi }G_i(𝐫𝐫^{},\omega ^{})W_{\left(0\right)}(𝐫𝐫^{},\omega \omega ^{})},`$ (9)
being $`i`$ is the iteration number. $`G_i`$ involves the QP energies obtained from $`\mathrm{\Sigma }_{\left(i1\right)}`$, without considering renormalization factors, lifetimes, and energy structures beyond QP peaks. As the QP bandstructure resulting from the first iteration is already in excellent agreement with experiment, our next step is to perform an “on mass-shell” $`G_1W_0`$ calculation. The resulting lifetimes are compared with experiment in Fig.1. One sees that forcing QP energies to appear also in the states involved in Eq. (7), i.e. the states forward which the quasihole decays, yields lifetime results similar to those of the “on mass-shell” $`G_0W_0`$ method. The same overestimation of the lifetimes of the top of the d bands is observed, confirming that a good agreement with experiment depends on the inclusion of the “intraband” decay channels described by Eq. (7).

A further question could now be addressed and is how self-consistency affects the quasiparticle bandstructure obtained within $`G_0W_0`$. As shown in Tab. 1, at higher iteration orders of the quasiparticle $`GW`$ equation the resulting energies worsen. The d–bands width decreases, reducing the agreement with experiment.

These results show that $`G_0W_0`$ describes correctly the d–holes lifetimes as far as “intraband” decay channels between d–like states are included. Those transitions, however, imply an inconsistency between the quasiparticle initial states and the DFT final states of the hole decay, as shown in Eq. (7). A self–consistent solution of Dyson equation removes this inconsistency, worsening, however, the agreement with the experimental results. This is a clear indication of the need of including vertex corrections in the self–consistent procedure. As shown for simple metals shirley , vertex corrections would partially cancel the dressing of the $`G_0`$ Green’s function of Eq. (6), restoring the “intraband” decay channels and, consequently, the $`G_0W_0`$ results.

In conclusion, we have shown that the lifetimes of d–holes in copper, calculated within the $`G_0W_0`$ method, are in good agreement with the experimental results, and can be obtained within the very same scheme which yields a good quasiparticle bandstructure. In contrast, further iterations of the QP equation beyond the $`G_0W_0`$ level yield worse results, for both the real and imaginary parts of self–energy. This can be explained by the need of including also vertex corrections together with self–consistency. This result is quite general and should apply to all metals in which one has two or more sets of electronic states with different degrees of spatial localization.

This work has been supported by the INFM PRA project “1MESS”, MURST-COFIN 99, Basque Country University, Iberdrola S.A. and DGES. and by the EU through the NANOPHASE Research Training Network (Contract No. HPRN-CT-2000-00167). We thank F. Aryasetiawan and P.M. Echenique for helpfull discussions.


text:
# The phase transition in QCD with broken 𝑆𝑈(2) flavour symmetry

## Abstract

We report the first investigation of the QCD transition temperature, $`T_c`$, for two flavours of staggered quarks with unequal masses at lattice spacings of $`1/4T`$. On changing the $`u`$/$`d`$ quark mass ratio in such a way that $`m_{\pi ^0}^2/m_{\pi ^\pm }^2`$ changes from 1 to 0.78, thus bracketing the physical value of this ratio, we find that $`T_c`$ remains unchanged in units of both $`m_\rho `$ and $`\mathrm{\Lambda }_{\overline{MS}}`$.

preprint: TIFR/TH/02-27, hep-lat/0208019

Lattice simulations of QCD with dynamical quarks aim to generate weights for the Euclidean path integral

$$Z=𝒟U\mathrm{exp}[S]\mathrm{Det}M_u\mathrm{Det}M_d,$$
(1)
discretised on a space-time lattice. Here $`U`$ are gauge fields which enter the action $`S`$ and the determinants of the Dirac operator $`M`$ for both $`u`$ and $`d`$ flavours. This leads to the well-known doubling problem in the chiral limit. Several solutions are known and have been used fairly extensively. One is to work with Wilson quarks, for which chiral symmetry is broken through an irrelevant operator at finite lattice spacing and recovered in the continuum. Among the new solutions is the domain wall definition of quarks dwf  in which the theory is extended to five dimensions, and by suitable tweaking of this extended action chiral symmetry is obtained on a four dimensional slice through the lattice when the length of the fifth dimension is sent to infinity. Another recently discovered solution is the overlap definition overlap , in which chiral symmetry is intact but its generator and the Dirac operator are non-local at finite lattice spacing, but become local in the continuum.

In this study we used the solution popular in finite temperature studies— staggered quarks. These have an exact continuous chiral symmetry at all lattice spacings, $`a`$, and in the continuum limit give 4 degenerate flavours in 4-dimensions. Two degenerate flavours are obtained by the prescription

$$\mathrm{Det}M_u\mathrm{Det}M_d=\left(\mathrm{Det}M_{stag}\right)^{1/2},$$
(2)
where $`M_{stag}`$ is the determinant for a single staggered quark field.

In this paper we introduce two staggered quark fields and define two flavours by the prescription

$$\mathrm{Det}M_u\mathrm{Det}M_d=\left(\mathrm{Det}M_{stag(u)}\mathrm{Det}M_{stag(d)}\right)^{1/4},$$
(3)
which we call 1+1 flavours. When the two quark masses are degenerate, this definition gives the same weight in the partition function of eq. (1) as the definition in eq. (2). The new prescription allows us to handle problems involving breaking of vector $`SU(2)`$ flavour symmetry, such as unequal masses for the $`u`$ and $`d`$ quarks (which we explore in this paper) and putting isovector chemical potential on the lattice isovector .

Our strategy here is to proceed in two steps— first to compare 2 flavours with 1+1 flavours of degenerate light quarks, and next to lift the degeneracy by giving unequal masses to the quarks (we shall use the notation $`\overline{m}=(m_u+m_d)/2`$). In each case we simulate the finite temperature theory on $`4\times 8^3`$ lattices for several $`\overline{m}`$. We study the real part of the Wilson line, $`\mathrm{Re}L`$, the quark condensate, $`\overline{\psi }\psi =(\overline{u}u+\overline{d}d)/2`$ and the related susceptibilities

$`\chi _L`$ $`=`$ $`V\left((ReL)^2ReL^2\right),`$
$`\chi _F`$ $`=`$ $`V\left((\mathrm{Tr}M^1)^2\mathrm{Tr}M^1^2+(\mathrm{Tr}M^2)\right)=\chi _m+\chi _\pi ,`$ (4)
where $`\chi _m`$ is the chiral susceptibility defined in bi  and $`\chi _\pi `$ is the pion susceptibility of sg . $`\chi _F`$ has the simple stochastic estimator pushan

$$\chi _F/V=\overline{(r^{}M^1r)^2}\overline{r^{}M^1r}^2,$$
(5)
where each $`r`$ is a complex random vector with each component drawn from a Gaussian distribution of unit width and the bar denotes an average over such an ensemble. We also study the integrated autocorrelation time, $`\tau `$ auto , of $`\mathrm{Re}L`$ and $`\overline{\psi }\psi `$, measured self-consistently over runs of length of at least $`30\tau `$. $`\beta _c`$ can be defined by peaks in $`\chi _L`$, $`\chi _F`$ and $`\tau `$. Within the accuracy of our determinations, they agree. We also make zero temperature measurements of the plaquette, $`P`$, and the masses of the pion, $`m_\pi `$, and the rho meson, $`m_\rho `$, at $`\beta _c`$ on $`8^4`$ lattices (we confirm that the lattice is large enough for these measurements by checking that the zero temperature measurements are independent of boundary conditions). The measurement of $`P`$ at $`T=0`$ also allows us to extract $`\alpha _S`$ blm  and hence to convert $`a`$ into a physical length scale in order to extract $`T_c/\mathrm{\Lambda }_{\overline{MS}}`$.

Our measurements have been taken at bare quark masses, $`a\overline{m}=0.15`$, $`0.1`$, $`0.05`$ and $`0.025`$ for degenerate flavours. For non-degenerate flavours with $`m_d=2m_u`$ we have only worked with the first three values of $`\overline{m}`$. We have also investigated $`m_d=10m_u`$ with $`a\overline{m}=0.15`$. The simulations have been performed with the Hybrid-R algorithm ralg . This algorithm is based on a molecular dynamics (MD) evolution which treats degenerate 1+1 and 2 flavours differently. Hence, the agreement in the results of these two cases, which we show later, is a good test that the numerical treatment of the discretised MD evolution is free of errors at the level of accuracy we achieve. The trajectory lengths have been taken to be one unit of molecular dynamics time, and the MD equations have been integrated over this range in 100 steps. Some typical run time histories of $`\mathrm{Re}L`$ are shown in Figure 1. Run time histories of other quantities are similar. The slow but large fluctuations visible in some of these runs are typical of critical slowing down, and therefore yield estimates of $`\beta _c`$.

Our measurements of $`\mathrm{Re}L`$ and $`\overline{\psi }\psi `$ are collected in Figure 2. Clearly there are no statistically significant differences between the measurements with degenerate quark masses and for $`m_d/m_u=2`$. From these figures it is also clear that $`\beta _c`$ is not sensitive to variation of $`m_d/m_u`$ in the range 1–2. From our measurements of masses, detailed later, it turns out this range of bare quark mass ratios includes the range of $`m_{\pi ^0}^2/m_{\pi ^+}^2`$ between 1 and roughly 0.8. Since the physical pion masses pdg  yield $`m_{\pi ^0}^2/m_{\pi ^+}^2=0.935`$, which is inside this range, we may conclude that the shift in $`\beta _c`$ due to realistic vector $`SU(2)`$ flavour symmetry breaking is negligible at the level of accuracy we can reach. A significant downward shift in $`\beta _c`$ is seen when the ratio $`m_d/m_u`$ is made as large as 10,

The measurements of autocorrelations yield strong peaks, sometimes as high as 300 MD time units, at couplings which we identify as $`\beta _c`$. These runs were extremely time consuming, since some of them required more than $`10^4`$ trajectories for reliable estimates of the susceptibilities shown in Figures 3 and 4. These also peaked at the same values of $`\beta `$, and are our primary means for the identification of $`\beta _c`$. Very good agreement between measurements of $`\chi _L`$ for $`m_d/m_u=2`$ and 1 indicate that there is no shift in $`\beta _c`$. The only exception is the value of $`\chi _L`$ near $`\beta _c`$ for $`a\overline{m}=0.15`$, where the difference can be attributed to the factor two difference in statistics between the runs with $`m_d/m_u=2`$ and 1— the former have both larger statistics and $`\chi _L`$. In Figure 4 we show $`\chi _F`$, which also peak at the same $`\beta _c`$, but unlike $`\chi _L`$ have stronger quark mass dependence.

Our estimates of $`\beta _c`$ for 1+1 flavours are shown in Table 2. Measurements of $`\beta _c`$ with the 2 flavour definition of staggered quarks have been made earlier in gott  on $`4\times 8^3`$ lattices with $`a\overline{m}=0.1`$, $`0.05`$ and $`0.025`$. Estimates of $`\beta _c`$ have also been made on lattices with $`N_t=4`$ for $`a\overline{m}=0.025`$ and $`0.0125`$ in fuku . These earlier results are fully in agreement with our measurements. Recall that studies on fixed volumes cannot decide the order of the transition, or indeed whether the peaking of the susceptibilities are due to a phase transition or a cross-over. That requires a study with varying volumes— a work that is outside the scope of this investigation and left for later.

We turn now to measurements at zero temperature. For this part of the work we generated configurations on $`8^4`$ lattices at the coupling $`\beta _c`$ for the value of $`a\overline{m}`$ under study. The autocorrelations of $`\overline{\psi }\psi `$ were found to be less than 5 MD time units in all cases, and the thermalisation time was less than 50 trajectories. We discarded the first 50 trajectories and stored 50 configurations separated by 5 trajectories in each of these simulations.

Meson correlators were computed on the stored configurations by inverting the staggered Dirac operator with appropriate masses. For $`m_um_d`$, the flavour combinations which give the meson propagators are—

$`C_{\pi ^+}(r)`$ $`=`$ $`M_u^1(0,r)M_d^1(r,0)`$
$`C_{\pi ^0}(r)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[M_u^1(0,r)M_u^1(r,0)+M_d^1(0,r)M_d^1(r,0)\right],`$ (6)
where $`r`$ labels a point of the lattice with respect to the chosen origin. The correlator for $`\pi ^{}`$ is obtained by flipping the roles of the $`u`$ and $`d`$ quarks in the expectation value on the right. The analogous combinations for the $`\rho `$ have the usual staggered fermion phase factor. The correlators, and hence the masses, for opposite charge states are identical. For degenerate quark masses, the propagators in all the charge states are identical. As usual, these correlators are summed over spatial slices to give the zero momentum propagators.

Meson masses were obtained by fitting these zero momentum correlators to the form

$$C(t)=A\mathrm{cosh}\left[m\left(\frac{L}{2}t\right)\right]+A^{}\mathrm{cosh}\left[m^{}\left(\frac{L}{2}t\right)\right],$$
(7)
or the corresponding single mass formula obtained by dropping the second term. We explored the stability of the fits by changing the range of $`t`$ over which the single mass form was fitted, and by comparing the result with the two mass fits.

The detailed comparisons are given in Table 1. For the flavour symmetric pion and the $`\pi ^+`$, there is a clear stable mass for fits over the range $`3at13`$. For the $`\pi ^0`$ the mass extracted from the two mass fit lies below that obtained from single mass fits. In all cases we took the lowest estimate to be the value of the mass. For the various $`\rho `$, the fit errors are larger and we accept the single mass fit over the full range of $`t`$ as our best estimate of the mass. The 1+1 flavour masses obtained with $`m_u=m_d`$ are fully consistent with previous measurements for 2 flavours gott .

We draw attention to the fact that the mass splitting between the charged and uncharged $`\rho `$ is not visible within measurement errors, and both these masses are equal to that obtained with degenerate quarks. However, the splitting between the neutral and charged pions is clearly visible. Interestingly, the latter are statistically indistinguishable from the pion mass measured with equal quark masses. Breaking $`SU(2)`$ flavour symmetry only results in lowering the $`\pi ^0`$ mass.

From our extraction of $`am_\rho `$, we extracted the ratio $`T_c/m_\rho =1/N_tam_\rho `$. These estimates are collected in Table 2. The ratio $`T_c/\mathrm{\Lambda }_{\overline{MS}}`$ increases marginally as $`a\overline{m}`$ decreases, due to the small decrease in the measured values of $`am_\rho `$, We have extrapolated these estimates to the values of $`a\overline{m}`$ for which the physical value of $`m_\pi /m_\rho `$ is obtained, through the form—

$$m_\rho =a+b\overline{m}+\mathrm{}.$$
(8)
For the value of $`a\overline{m}`$ where the physical value of $`m_\pi /m_\rho `$ is obtained, we find $`T_c/m_\rho =0.227\pm 0.008`$. This gives $`T_c=175\pm 6`$ MeV. Mutually consistent values of $`T_c`$ are obtained using the 2-flavour data of gott , and our 1+1 flavour data with $`m_d/m_u=1`$ as well as 2.

Alternatively, the lattice spacing can be traded for $`\mathrm{\Lambda }_{\overline{MS}}`$ by using our measurements of the plaquette, $`P`$, to extract the running coupling, $`\alpha _S`$, and using this in the 2-loop $`\beta `$ function. This procedure was introduced in blm  and later used to determine $`T_c/\mathrm{\Lambda }_{\overline{MS}}`$ from the 2 flavour data for realistic quark masses in tclam . With such an analysis, the values of $`T_c/\mathrm{\Lambda }_{\overline{MS}}`$ extracted from our measurements in the 1+1 flavour case agree completely with the previous values at the same $`a\overline{m}`$. Extrapolation of our 1+1 flavour results to quark masses which give the correct value of $`m_\rho /\mathrm{\Lambda }_{\overline{MS}}`$, yield $`T_c/\mathrm{\Lambda }_{\overline{MS}}=0.49\pm 0.02`$, as in tclam . Then, using the 2-loop value of $`\mathrm{\Lambda }_{\overline{MS}}=343_{28}^{+31}`$, as appropriate below the charm quark threshold pdg , we obtain $`T_c=167\pm 9_{14}^{+15}`$ MeV (where the first error comes from $`T_c/\mathrm{\Lambda }_{\overline{MS}}`$ and the second from $`\mathrm{\Lambda }_{\overline{MS}}`$). The two estimates of $`T_c`$ are totally compatible with each other, and also from other estimates using improved actions improved .

The above method depends on the expansion of $`P`$ at $`T=0`$ in a series in $`\alpha _S`$. Since this is scheme dependent schemes , the extraction of $`T_c/\mathrm{\Lambda }_{\overline{MS}}`$ can have large uncertainties if the cutoff is large. To estimate these, we have extrapolated the temperature scale to smaller lattice spacings and compared to direct measurements of the scale. At the smallest $`\overline{m}=0.025/a`$ we have made a few $`T=0`$ runs at larger $`\beta `$ to estimate the lattice spacing through a measurement of $`P`$. Since the lattice spacing, $`a`$, is an outcome of this computation, a series of runs is needed to tune $`\overline{m}`$. The scale determined in this way can be converted into a temperature scale for simulations with $`N_t=4`$ at the corresponding $`\beta `$. This is shown in Figure 5. The temperature scale for $`N_t=4`$ can also be calibrated by direct measurements of $`\beta _c`$ at larger $`N_t`$. For 2 flavours of quarks and $`\overline{m}=0.025/a`$, such measurements have been performed milc . These results are also plotted in Figure 5. The good agreement between the two methods of setting the scale implies that the lattice spacings are fine enough for 2-loop scaling to work. As a result, we expect that $`T_c/\mathrm{\Lambda }_{\overline{MS}}`$ obtained here are relevant to the continuum limit.

The major remaining uncertainties are in the extrapolation to zero quark mass, and in possible power corrections in $`a`$ to various quantities we have measured. The question of the order of the phase transition needs a finite size scaling study and has not been addressed here. A detailed study of these issues lies outside the scope of this paper, and is left to the future.


text:
# Bulk magnetic properties and phase diagram of Li-doped La2CuO4: Common magnetic response of hole-doped CuO2 planes

## I INTRODUCTION

Depending on the nature and concentration of dopants, La2CuO4 (LCO) displays a wide variety of phenomena such as antiferromagnetism, spin glass (SG) behavior, an anomalous metallic response, and high-temperature superconductivity (HTS). Undoped LCO contains weakly coupled CuO2 planes and exhibits antiferromagnetic (AF) order of the Cu2+ spin-1/2 moments below $`T_N325`$ K. Replacement of Cu2+ by nonmagnetic Zn2+ or Mg2+ models random spin dilution, leading the system into a disordered state at doping concentrations above $`\mathrm{\hspace{0.17em}40}\%`$. On the other hand, substitution of divalent alkaline earth cations for La3+ or the introduction of excess interstitial oxygen introduces hole charge carriers into the CuO2 planes which frustrate the spin system. A concentration of $`x=0.02`$ in La2-xSrxCuO4 (Sr-LCO) is enough to destroy the AF order; this value is one order of magnitude smaller than for the spin-dilution case. Spin freezing has been found at low temperatures in the AF doping regime ($`x<0.02`$), with recent direct evidence for electronic phase separation. Further doping ($`x>0.02`$) leads to the emergence of a SG phase, followed by superconductivity for $`x0.060.25`$. SG order is found to coexist with superconductivity, with no apparent anomaly in the SG temperature at the doping level $`x0.06`$ at which superconductivity first occurs. The phase diagram of the substitutionally-doped bilayer material Y1-xCaxBa2Cu3O6 (Ca-YBCO) closely resembles that of Sr-LCO, and spin freezing has also been found in LCO (Ref. 11) and YBCO (Ref. 12) doped with excess oxygen.

The extrapolated disappearance of SG order in Sr-LCO and Ca-YBCO appears to coincide with the doping level at which the normal state pseudogap extrapolates to zero, and it has been suggested that this might be consistent with predictions involving quantum criticality. In the $`d`$-density wave picture of HTS, the reason for the lack of a genuine phase transition at the pseudogap temperature is that the disorder present in all existing cuprates corrupts the $`d`$-density wave order and transforms the transition into the low-temperature SG transition. While the freezing of $`d`$-density wave fluctuations is one proposal for the origin of the SG in doped cuprates, there also exist other interpretations. For example, it has been argued that the glassiness found in doped Mott insulators may be self-generated, due to the competition between interactions on different length scales, and that quenched disorder may merely further stabilize SG order. The most discussed scenario has been the so-called cluster SG, with holes on the cluster boundaries and in the clusters giving rise to SG physics and to the experimentally observed incommensurate spin correlations.

Given the enormous interest in the connection between magnetic correlations and HTS, it would be valuable to investigate the detailed magnetic properties of related, nonmetallic materials such as La2Cu1-xLixO4 (Li-LCO). Since Li+ not only provides one hole carrier, but also removes a Cu2+ spin, this system experiences the dual effects of spin dilution and frustration. Unlike Sr-LCO, Li-LCO does not superconduct and shows no evidence of incommensurate AF correlations. However, early work on polycrystalline samples reported a rapid suppression of AF order with doping, similar to Sr doping rather than Zn doping, and evidence for spin freezing in the Néel state. The remarkable similarity of magnetic properties of lightly doped Li-LCO with those of Sr-LCO has been interpreted as due to new collective behavior of the holes. Therefore, in connection with HTS, the detailed magnetic properties and phase diagram of Li-LCO are of considerable interest. Experimentally, spin freezing is observed at a different temperature, depending on the time scale of the probe. Unlike NQR, $`\mu `$SR, and neutron scattering, magnetometry using a superconducing quantum interference device (SQUID) is essentially a static probe, allowing a more accurate extraction of the SG transition temperature $`T_{\mathrm{sg}}`$. Our results for lightly-doped samples demonstrate the feasibility of using bulk magnetometry to extract spin freezing temperatures in the Néel regime. The present magnetometry study establishes the existence of a nearly quantitative agreement of the bulk magnetic phase diagram of Li-LCO with that of Sr-LCO. Since the strength and type of the disorder as well as the charge transport differ significantly in these two materials, this finding places constraints on the origin of the SG degrees of freedom in hole-doped cuprates.

## II EXPERIMENTAL

Using the traveling-solvent floating-zone method, we have succeeded in growing large single crystals of La2Cu1-xLixO4 (0 $``$ $`x`$ $``$ 0.07). The crystal axes were precisely determined by the x-ray Laue backscattering technique. In order to eliminate possible hole doping by excess oxygen, the crystals were carefully heat treated under reducing conditions. The Li concentrations were estimated to within $`\pm 0.003`$ from x-ray diffraction measurements of the lattice constants. For a few samples, we confirmed this estimate using neutron diffraction to determine structural and Néel transition temperatures. All magnetometry data reported here were taken with a commercial SQUID magnetometer with the magnetic field along the tetragonal $`a`$-axis.

## III RESULTS AND DISCUSSIONS

Figure 1 summarizes the temperature and doping dependence of the magnetic susceptibility of La2Cu1-xLixO4, taken with a field of $`H`$ = 0.1 kOe along the tetragonal $`a`$-axis. The data exhibit several features which systematically shift with doping up to $`x`$ = 0.03, and change their nature above this doping level. Earlier powder studies had demonstrated a rapid suppression of the AF order with Li doping, and the Li concentration at which Néel order is lost corresponds to $`x`$ 0.03. The doping dependence of the Néel temperature determined by local and bulk magnetic probes agree quite well, and until now this has been the only magnetic phase boundary known in the Li system. However, the present results indicate that the phase diagram in Li-LCO is more complicated than what has been determined from powder samples.

$`\mu `$SR results indicate magnetic order at low temperatures even up to $`x`$ 0.10, but bulk magnetization measurements in powder samples only show spin paramagnetism for $`x>0.03`$. The temperature dependent susceptibility for $`x>0.03`$ can therefore be expressed by the extended Curie-Weiss formula:

$$\chi (T)=\chi _0+\frac{C}{(T\Theta )},$$
(1)
where $`\chi _0`$, $`C`$, and $`\Theta `$ are the $`T`$-independent susceptibility, Curie constant, and Curie-Weiss temperature, respectively. Figure 2 shows $`\chi (T)`$ for $`x=0.05`$. As in the earlier powder report, paramagnetism is observed in the high temperature regime. A fit to Eq. (1) (solid line in Fig. 2) resulted in $`\chi _0`$ = 2.35 $`\times `$ 10-7 emu/g, $`C`$ = 7.18 $`\times `$ 10-6 emu$``$K/g, and $`\Theta <`$ 0.5 K. We note that these parameters are roughly comparable with those Sr-LCO: $`\chi _0`$ = 0.4 $`\times `$ 10-7 emu/g for $`x`$ = 0.04, $`C`$ = 2$``$5 $`\times `$ 10-6 emu K/g, and $`\Theta `$ = 0 K for $`x`$ = 0.03$``$0.05.

At low temperatures, $`\chi (T)`$ deviates significantly from a Curie-Weiss law and shows signatures of a SG. As described in detail below by the appropriate limits of the scaling function, the SG order parameter in the zero-field limit increases from zero upon cooling below $`T_{\mathrm{sg}}`$, which lead to a decrease of $`\chi (T)`$ for $`T<T_{\mathrm{sg}}`$. Furthermore, the resulting peak becomes broader as the applied magnetic field is increased due to the enhancement of the order parameter both above and below $`T_{\mathrm{sg}}`$. This behavior is demonstrated in the inset of Fig. 2. From the peak in the lowest field ($`H=0.1`$ kOe) we determined $`T_{\mathrm{sg}}=6.2(1)`$ K, which should be compared to $`T_{\mathrm{sg}}=5.0(5)`$ K for Sr-LCO with $`x=0.05`$. Hysteresis below $`T_{\mathrm{sg}}`$, observed between zero-field-cooled and field-cooled curves, is indicated with the help of arrows in the inset of Fig. 2. Such a behavior is also characteristic of a SG. Observing these features in either small crystals or polycrystalline samples would be difficult because they become obscured as the field is increased, or when the field direction is canted away from the CuO2 plane.

The large oriented single crystals used in the present study have enabled us to further characterize the SG state by means of a critical scaling analysis of the SG order parameter $`q`$, which is experimentally associated with the deviation of the observed equilibrium (field-cooled) susceptibility $`\chi (T,H)`$ from Curie behavior. Normalization to satisfy 0 $`q`$ 1 yields

$$q(T,H)=\left[\left(\chi _0+\frac{C}{T}\right)\chi (T,H)\right]/\left(\frac{C}{T}\right).$$
(2)
Theoretically, the SG transition should obey a scaling relation as observed in other critical phenomena. Using the reduced temperature $`t`$ = ($`T`$ $``$ $`T_{\mathrm{sg}}`$)/$`T_{\mathrm{sg}}`$ and a scaling function $`F_\pm (z)`$, the SG order parameter is defined as

$$q(T,H)=|t|^\beta F_\pm (H^2/|t|^\varphi ),$$
(3)
where $`\beta `$ and $`\gamma =\varphi \beta `$ are the critical exponents characterizing the SG state. The behavior of the scaling function $`F_\pm (z)`$ is well known in the following three limits: (i) $`F_+(z0)`$ = 0 ($`t>0`$), (ii) $`F_{}(z0)`$ = const ($`t<0`$), and (iii) $`F_\pm (z\mathrm{})`$ = $`z^{\beta /\varphi }`$. Because a relation $`qt^\beta `$ is immediately found from limit (ii), the exponent $`\beta `$ can be estimated from the slope in a log-log plot of $`q`$ (below $`T_{\mathrm{sg}}`$) versus $`t`$ in the limit of zero field. Such an analysis is performed in Fig. 3(a), and an estimate of $`\beta `$ = 0.7(1) is obtained. Similarly, using the limit (iii), the other exponent, $`\varphi `$, can be determined from the dependence of $`q`$ on $`H^2`$ at $`T_{\mathrm{sg}}`$ from the relation $`q(T_{\mathrm{sg}},H)`$ = $`(H^2)^{\beta /\varphi }`$. The slope in Fig. 3(b) then gives $`\beta /\varphi `$ = 0.163(5) or $`\varphi `$ = 4.3(8).

Now that we have a good estimate of the critical exponents, we can directly test the scaling relation described by Eq. (3). As shown in Fig. 3(c), we obtain excellent scaling. Through the process of getting the best scaling result, shown in Fig. 3(c), the critical exponents have been optimized further: $`\beta =0.78(5)`$ and $`\gamma =4.1(5)`$.

The same SG features as for $`x`$ = 0.05 were observed for $`x`$ = 0.06 and 0.07. The critical exponents obtained in this manner do not depend on $`x`$, although $`T_{\mathrm{sg}}`$ decreases with $`x`$. We note that, in addition to the comparable values of $`T_{\mathrm{sg}}`$, the critical exponents are found to be almost identical to those of Sr-LCO, suggesting the existence of the same SG state in both systems. While the observed exponents are consistent with those of canonical SG materials, we note that recent results for untwinned Sr-LCO crystals reveal unconventional anisotropic behavior.

For $`x<`$ 0.03, the magnetic susceptibility differs from that for $`x>`$ 0.03. Data for $`x`$ = 0.02 are shown in Fig. 4(a). Two magnetic anomalies are present. One is a high-temperature cusp associated with the AF transition which has already been observed in powder samples. As shown in Fig. 4(b), neutron data complement the SQUID observations, providing unambiguous evidence of a well-defined Néel transition around this high-temperature bump. We have confirmed the development of the antiferromagnetic peak at the (1,0,0) position (orthorhombic notation) just below the susceptibility anomaly, and determined the Néel temperature to be $`T_\mathrm{N}`$ = 135 K. The second anomaly, recognized well below $`T_\mathrm{N}`$, appears somewhat analogous to the SG behavior found for $`x>`$ 0.03 in that we find hysteresis below the anomaly and the onset temperature of magnetic irreversibility is comparable to $`T_{\mathrm{sg}}`$ for $`x>`$ 0.03. Neutron characterization furthermore revealed a sharp structural transition from the tetragonal to orthorhombic phase at 490 K, consistent with a doping level of $`x`$ = 0.02. This indicates that our samples have a high degree of chemical homogeneity, and that the observed magnetic anomaly is intrinsic to the $`x`$ = 0.02 phase. Local-probe 139La-NQR and $`\mu `$SR experiments have reported a corresponding spin-freezing within the AF state in Sr-LCO (Refs. 3 and 9) and Li-LCO. This has been ascribed to the continuous freezing of the spins of the doped holes on the antiferromagnetic background, with independent ordering of Cu2+ spins and doped holes, or due to a collective hole behavior. However, to the best of our knowledge, this spin freezing has not been reported thus far from bulk magnetometry for the cuprates. Because the SG behavior is very sensitive to the dopant concentration ($`x`$-dependence will be discussed later) and to the field direction, the success of the present observations depended greatly on high-quality single crystals as well as extended measurements in the low-temperature region.

Although the second anomaly appears to be a SG transition, further characterization in terms of critical scaling analysis is not possible due to the ambiguity in the definition of the SG order parameter. Nevertheless, it would be worthwhile to determine a characteristic freezing temperature, $`T_f`$. We define $`T_f`$ as the lower-temperature “shoulder” of the zero-field-cooled data at low $`H`$, which corresponds to the onset of magnetic hysteresis for $`x`$ = 0.02, as seen in the inset of Fig. 4(a). We note that, although not shown here, a change in the imaginary part of the ac susceptibility was also observed at this temperature, which can be associated with the drastic change of the spin response due to the onset of the spin freezing.

Our results yield the phase diagram summarized in Fig. 5. NQR ($`x<0.02`$) (Ref. 3) and magnetometry ($`x>0.02`$) (Ref. 7) data from Sr-LCO single crystals are included for comparison. Despite the different nature of the dopants, the phase boundaries for these two compounds are almost identical. Figure 5 demonstrates that $`T_\mathrm{f}x`$ as in Sr-LCO, opposite to the trend with doping for either $`T_N`$ or $`T_{\mathrm{sg}}`$. We find $`T_\mathrm{f}`$ = (339 (14) K)$`\times x`$, as compared to $`T_\mathrm{f}`$ (815 K)$`\times x`$ for Sr-LCO from 139La NQR.

The dopant concentrations at which Néel order is lost in the two materials differ by a factor of $``$ 2/3. Qualitatively, this shift of the phase boundaries can be understood to result from reduced magnetic frustration in the Li-doped system. As shown in the inset of Fig. 5, when normalized by that factor, $`T_{\mathrm{sg}}`$ for Li-LCO ($`x>0.03`$) $`quantitatively`$ agrees with the SG temperature for Sr-LCO. Since NQR and magnetometry probe very different time scales, systematic Sr-LCO magnetometry data for $`x<0.02`$ would be very desirable for a proper quantitative comparison at low doping.

Recent neutron scattering results reveal direct evidence for electronic phase separation in Sr-LCO ($`x<0.02`$) into regions with hole concentrations $`0`$ and $`0.02`$. The latter phase exhibits diagonal stripes. Spin freezing occurs below the doping-independent phase separation temperature $`T_{ps}30`$ K determined from neutron scattering, and consistent with previous NQR results. For Li-LCO, on the other hand, commensurate AF correlations have been both predicted and observed, consistent with the expected stronger pinning potential of the in-plane dopant Li+. Consequently, the cluster model, originally proposed for Sr-LCO, but which predicts commensurate AF correlations, might more accurately describe the physics of the Li-doped variant. This model predicts $`k_\mathrm{B}T_\mathrm{f}J_{\mathrm{eff}}x`$, where $`J_{\mathrm{eff}}`$ is the effective in-plane exchange coupling constant. While this linear doping dependence is indeed consistent with our observations, the same behavior is found in Sr-LCO. On the other hand, it has been speculated that $`T_\mathrm{f}`$ in Sr-LCO might depend linearly on the volume fraction of the SG phase. The close analogy between the NQR (Refs. 3 and 22) and bulk magnetic properties of the two materials suggests that Li-LCO might phase separate as well.

At higher hole concentrations, the behavior of the carriers differs significantly between Li-LCO and Sr-LCO. The latter material shows an insulator to metal transition and HTS, while the former remains insulating. What is particularly interesting from our present observations is that, in spite of the great difference in charge dynamics for $`x>`$ 0.03, the spin degrees of freedom measured via magnetometry are remarkably similar in both compounds. Both materials show a spin glass transition, with comparable transition temperatures at the same doping value and with the same critical exponents. The values of $`T_{\mathrm{sg}}(x)`$ for the more disordered material Li-LCO lie above those of Sr-LCO. The relative shift of the SG phase boundaries likely results from an effective decrease of magnetic frustration due to the presence of nonmagnetic Li+ in Li-LCO. In this context, it is worth noting that the magnetic phase boundaries in the double-layer cuprates lie above those of both Li-LCO and Sr-LCO. While differences in the strength and type of disorder might play a subtle role in setting the temperature scales for Néel and SG order, we believe that the predominant effect is the difference in the effective three-dimensional AF coupling between single- and double-layer materials. Specifically, while the interplanar AF coupling is nearly frustrated in doped LCO, this is not the case for the double-layer materials. On the other hand, the relative insensitivity of the low-temperature magnetic phase boundaries of the structurally identical materials Li-LCO and Sr-LCO to the type and strength of the quenched disorder is consistent with the notion that the glassiness in these doped Mott insulators is primarily self-generated.

## IV CONCLUSION

The bulk magnetic properties of Li-LCO single crystals are found to be richer than previously reported from powder samples. Our results demonstrate the feasibility of extracting low-temperature freezing temperatures in the Néel state from magnetometry, and call for similar measurements in other cuprates for a proper quantitative comparison. We find that $`T_\mathrm{f}x`$, which also has been reported for Sr-LCO. Since the spin freezing in the latter material occurs in a phase separated state, the close similarity between the doping dependence of the freezing temperature suggests that lightly doped Li-LCO might phase separate as well, despite the stronger pinning potential of the in-plane dopant Li+. Outside of the Néel phase, the spin glass phase boundaries differ by a mere scale factor for the effective doping level. The experimental results for Li-LCO obtained here extend the universality of the bulk magnetic phase diagram to hole-doped but non-superconducting cuprates.

## ACKNOWLEDGMENTS

We thank D. Fisher, H. Eisaki, K. Kishio, R.S. Markiewicz, Ch. Niedermayer, C. Panagopoulos, and H. Takagi for fruitful discussions. P.K.M., O.P.V., and M.G. would like to thank the staff of the NIST Center for Neutron Research for their warm hospitality. The work at Stanford was supported by the U.S. DOE under Contract Nos. DE-FG03-99ER45773 and DE-AC03-76SF00515, by NSF CAREER Award No. DMR-9985067, and by the A.P. Sloan Foundation. T.S. was partially supported by JSPS.


text:
# Saturation and Pion Production in Proton-Nucleus Collisions

\[

## Abstract

We study the effects of gluon saturation on pion production in high energy proton-nucleus collisions using the color glass condensate model. At high $`p_{}`$, we show that the $`p_{}`$-distribution of gluons behaves as $`1/p_{}^7`$ in accordance with both conventional perturbative QCD calculations and experiment. Fragmentation of gluons into pions leads to a rapidity dependent depletion of pions relative to the conventional perturbative QCD predictions. We argue that these clear and systematic differences provide a signal for the onset of gluon saturation which is accessible in upcoming experiments.

NORDITA-2002-35 HE, hep-ph/0208007

\]

Since the pioneering work of Ref. , the phenomenon of gluon saturation, occurring at very high gluon density, has been under intense study. One hopes that saturation would naturally tame the rise of the gluon distribution function at small values of Bjorken $`x`$ and circumvent the subsequent violation of unitarity.

A promising attempt to address the problem of saturation involves replacing the partonic fields by coherent classical fields. The justification is that for large occupation numbers the dynamics are essentially classical. This program was initiated in Ref.  and developed further in e.g. Refs. . Within this approach one is able to show that an intrinsic scale, the saturation scale $`Q_{S,A}(x)`$, is generated. This scale depends on both Bjorken $`x`$ and atomic number, $`A`$, and roughly delineates the dense and dilute regimes. Throughout this work we shall refer to this approach as the color glass condensate (CGC) framework.

The parametric dependence of the saturation scale has been extensively studied. Phenomenological models based on the idea of parton saturation have been successful in describing both the inclusive and diffractive HERA data for $`F_2`$ . According to this model, the dependence of the saturation scale on Bjorken $`x`$ is

$$Q_{S,A}^2(x)=(1\mathrm{G}\mathrm{e}\mathrm{V})^2\left(\frac{x_0}{x}\right)^\lambda A^{1/3},$$
(1)
where $`\lambda =0.25\mathrm{}0.3`$ and $`x_0=310^4`$. For the numerical calculations, we will use $`\lambda =0.3`$. Since $`Q_{S,A}^2`$ can be interpreted as the transverse density of partons, it is enhanced by considering larger nuclei at any given $`x`$ by a factor of $`A^{1/3}`$. As a result, one may expect to observe saturation effects better in nuclear collisions at currently available collider energies. Indeed, the present data from BNL–RHIC  seem to be compatible with saturation based phenomenology . These results, however, are far from being conclusive on account of numerous model dependencies related to our poor knowledge of the details of the spacetime evolution of the produced strongly interacting matter.

The onset of saturation has important qualitative and quantitative consequences for the $`p_{}`$-distribution of the produced gluons. This was already pointed out in Ref.  in the context of proton–proton collisions where an asymmetry in the saturation scales of the two protons appears away from the mid-rapidity. As one hadron moves deeper into the saturation region, towards smaller $`x`$, its saturation scale increases, while the saturation scale of the second hadron decreases as is evident from Eq. (1). In proton–nucleus collisions, a similar situation is realized already at mid-rapidity due to the atomic number asymmetry. In particular, one expects three different regions in the transverse momentum distribution of the produced gluons on account of the asymmetries in the saturation momenta, $`Q_{S,A}>Q_{S,p}`$. There is the region $`p_{}<Q_{S,p}<Q_{S,A}`$ in which both the proton and the nucleus are saturated. The second possibility is that the nucleus is saturated but the proton is in the linear regime ($`Q_{S,p}<p_{}<Q_{S,A}`$). The final region is the usual perturbative domain in which both gluon distributions are entirely governed by the linear evolution ($`Q_{S,p}<Q_{S,A}<p_{}`$). Consequently, proton–nucleus collisions offer the best possibility of experimentally probing the systematics of gluon saturation. The problem of evaluating gluon production in proton–nucleus collisions has been found to be analytically tractable above the smaller saturation scale and in Ref.  this was qualitatively studied in the CGC framework.

In this Letter, we undertake the quantitative computation of the produced hadronic spectra based on the framework of Ref.  and compare it to the conventional perturbative QCD (pQCD) calculation which does not include any saturation effects. The main emphasis is on a numerical calculation of the fragmentation of the gluon distributions into pions, and on the study of how saturation effects manifest themselves in the $`p_{}`$-distributions of the measured hadrons. We shall see that although the spectra at the gluonic level are qualitatively very different, this difference is made less striking by the fragmentation of the gluons into pions.

In pQCD, the inclusive cross section for gluon production is given as a convolution of the hard partonic scattering cross section with the gluon distribution functions. The hard cross section behaves as $`1/p_{}^4`$, but the final result receives additional contributions from the two gluon distributions and from the running of the strong coupling constant. As a result, the slope of the transverse momentum distribution actually falls much more rapidly with a scaling closer to $`1/p_{}^7`$. Any realistic model of particle production must reproduce this scaling at large transverse momenta.

We consider a CGC model obtained by combining two separate limits of the results of Ref. . The gluon production cross section in the perturbative domain is given by

$$\frac{d\sigma _g^{\mathrm{pert}}}{d^2p_{}dy}=\frac{8N_c(N_c^21)}{\pi }\frac{\alpha _s^3}{p_{}^4}\chi _p(x_p,p_{}^2)\chi _A(x_A,p_{}^2),$$
(2)
where

$`\chi _i(x_i,p_{}^2)`$ $`=`$ $`{\displaystyle \frac{1}{\pi R_i}}{\displaystyle \frac{N_c}{N_c^21}}`$ (3)
$`\times `$ $`\left({\displaystyle _{x_i}^1}𝑑xg_i(x,p_{}^2)+{\displaystyle \frac{C_F}{N_c}}{\displaystyle _{x_i}^1}𝑑xq_i(x,p_{}^2)\right),`$ (4)
with $`i=A,p`$. $`R_i`$ is the transverse radius, $`x_{A,p}=p_{}e^y/\sqrt{s}`$, $`C_F=(N_c^21)/(2N_c)`$ and the nucleus is chosen to have negative beam rapidity. Here $`q_i`$ and $`g_i`$ denote the valence quark and gluon distributions in the nucleus or in the proton. One needs to evaluate the color charge functions $`\chi _i`$ only in the weak field regime, and therefore the use of DGLAP evolved parameterizations is appropriate.

In the saturation domain of the nucleus, the cross section is given by

$$\frac{d\sigma _g^{\mathrm{sat}}}{d^2p_{}dy}=\frac{C(N_c^21)}{2\pi ^2}\frac{\alpha _s}{p_{}^2}\chi _p(x_p,p_{}^2)\pi R_A^2,$$
(5)
The constant $`C`$ is chosen so that the sub-cross section is continuous at $`Q_{S,A}`$. We have neglected the logarithms appearing in the results of Ref.  since one expects the overall behavior to be driven by the strong powers of $`p_{}`$ and the constant under the logarithm is in any case difficult to fix. The overall normalization is unimportant as we will discuss later.

From these equations, it is clear that retaining the DGLAP evolved parton distributions for the weak fields and letting the strong coupling run will affect the qualitative results of Ref.  which display the expected $`1/p_{}^4`$ and $`1/p_{}^2`$ behavior in the perturbative domain and in the saturation domain of the nucleus, respectively. The inclusion of these effects is essential for making the proper contact with conventional result for high $`p_{}`$ mentioned above.

In order to compute the produced hadronic spectra, we convolute the above inclusive gluon production cross section with the appropriate fragmentation functions. For the sake of simplicity, we do not attempt to account for the Cronin effect , nuclear shadowing  or nuclear modifications of the fragmentation functions . So, $`xg_A(x,Q^2)=Axg_p(x,Q^2)`$ in the above equations and similarly for quarks. The generalization to nucleus–nucleus collisions with $`A_1A_2`$ is trivial through replacements $`pA_1`$, $`AA_2`$. We note that the overall normalization of the results is difficult to fix. In Ref.  it was found that in order to reproduce the data from hadronic collisions above a cutoff momentum of $`p_01\mathrm{}3`$ GeV using pQCD one needs a $`\sqrt{s}`$-dependent $`K`$-factor to account for the higher order corrections. The relative normalization between the pQCD and CGC results, however, is almost fixed since both results should match at high $`p_{}`$. Therefore, we do not attempt to fine-tune the overall normalization, but rather consider a leading order pQCD calculation, including only gluons, without any $`K`$-factors and point out the differences relative to the CGC calculation which will be normalized at high $`p_{}`$ to the pQCD result. Our central results, the slopes of the $`p_{}`$\- and $`y`$-distributions, are not sensitive to the normalization ambiguities. One should note, however, that the $`\sqrt{s}`$ growth of the gluon multiplicity in the CGC model is driven by the small $`x`$-growth of the gluon distributions, and is therefore of the order of $`dN/dy\sqrt{s}^\lambda `$. We will focus exclusively on pions since this distribution approximates well the distribution of all hadrons.

To compute the fragmentation of gluons into pions, we assume that the pion is produced collinearly with its parent gluon, $`\eta _g=\eta _\pi `$, and carries a fraction $`z`$ of the parent gluon’s energy. The cross section for the production of pions is

$$\frac{d\sigma ^{pA\pi X}}{d^2q_{}dy}=J(m_{},y)\frac{dz}{z^2}D_g^\pi (z,q_{}^2)\frac{d\sigma ^{pAgX}}{d^2p_{}dy_g},$$
(6)
where

$`p_{}={\displaystyle \frac{q_{}}{z}}J(m_{},y),y_g=\mathrm{sinh}^1\left({\displaystyle \frac{m_{}}{q_{}}}\mathrm{sinh}y\right),`$ (7)
$`J(m_{},y)=\left(1{\displaystyle \frac{m_\pi ^2}{m_{}^2\mathrm{cosh}^2y}}\right)^{1/2}`$ (8)
and $`m_{}^2=q_{}^2+m_\pi ^2`$. The integration over $`z`$ is limited by the maximal energy the gluon can carry and by the lower bound on $`p_{}`$. In pQCD, this latter scale is a fixed cutoff, while in the CGC calculation it is determined by the saturation scale of the proton,

$$\frac{am_{}}{\sqrt{s}}\mathrm{cosh}yz\mathrm{min}(1,\frac{q_{}}{p_{,min}}J(m_{},y)),$$
(9)
where $`a=1,2`$ for the CGC and pQCD calculations respectively. This difference arises from the underlying parton kinematics. In LO pQCD, one produces two minijets, back-to-back in the transverse plane, carrying at most $`E=\sqrt{s}/2`$. In the CGC calculation, the scattering is a BFKL-type $`21`$ fusion, and the produced gluon can have energy up to $`E=\sqrt{s}`$. We use the CTEQ5 parameterization of parton distributions  and KKP parameterization of the fragmentation functions .

In Fig. 1, we plot the $`p_{}`$-distributions of the produced gluons at different rapidities. The topmost three curves are for $`y=0`$ and the lower two are for $`y=3`$. The dashed curve is the CGC result including only gluons in the sources $`\chi _i`$. Note that at $`y=0`$ the slopes of the CGC and pQCD calculations almost match at high $`p_{}`$. Below the saturation scale of the nucleus, the two results strongly deviate. In a full computation of the gluon distribution, this sharp bend at $`Q_{S,A}(y)`$ will become smooth. The essential point is that there is a marked depletion of gluons for $`p_{}<Q_{S,A}(y)`$.

If this spectrum were measured, one should be able to determine whether the parton dynamics in the collision are better described by pQCD or by the CGC model. To obtain a measurable spectum one needs to fragment each of the distributions in Fig. 1 to hadrons.

In Fig. 2, we show the $`p_{}`$-distribution of the pions produced from the two models. The fragmentation changes the qualitative behavior of the spectrum of the CGC model dramatically with the sharp bend disappearing completely on account of the redistribution of momentum from the gluons to the pions. We note that there are significantly fewer produced pions at lower transverse momentum in the CGC calculation than in the conventional quantum pQCD one and that this deficit is strongly rapidity dependent. We have checked using a smooth fit interpolating between the two limiting forms of the gluon distribution that the qualitative, and to some extent even quantitative, results for the pion spectrum in the CGC model are not sensitive to the sharpness of the knee in the gluon distribution.

To obtain an useful simple parameterization of the pion spectrum, we consider the form

$$\frac{d\sigma }{dydp_{}^2}=\frac{C(y)}{(p_{}^2+\mu (y)^2)^{b(y)}},$$
(10)
where the emergence of a rapidity dependent scale, $`\mu (y)`$, is a natural consequence of the rapidity dependence of the saturation scale. The rapidity dependence of the effective slope $`b(y)`$ is expected since the slope of the gluon distribution already depends on rapidity. We find that $`\mu (y)^2=0.50e^{\lambda y}`$, where $`\lambda =0.3`$. The values of the fit parameters $`C(y)`$ and $`b(y)`$ at various rapidities are collected in Table I. The inclusion of the quarks plays a role for the overall normalization and the slopes are left almost unchanged as can be seen already from Fig. 1

Finally, we consider the rapidity distributions at fixed $`p_{}`$. One cannot approach the beam rapidity of the nucleus in this framework, since at such large negative rapidities the saturation scales of the nucleus and the proton become comparable and the underlying model breaks down. Near midrapidity, the CGC model gluon spectrum in perturbative domain is characterized by constant behavior in rapidity, as is the case for pQCD spectrum. In the saturation domain, the CGC gluon spectrum has a nonzero slope also at midrapidity. The CGC model gluon spectra are shown in Fig. 3 with dotted curves for $`p_{}=3`$ and 1.4 GeV, and the corresponding pion spectrum by solid lines for $`p_{}=0.4,0.60.81.0`$ and 1.7 GeV. The pion spectrum from pQCD, shown by dashed lines for 0.4 and 1.0 GeV has the flat behavior at midrapidity over the whole $`p_{}`$ range. We observe that the small $`p_{}`$ pions might reveal the characteristic form carried by the gluon spectrum in the saturation domain. This result is easy to understand, since the pion spectrum at scale $`q_{}`$ effectively probes the gluon spectrum at some higher scale $`q_{}/\overline{z}`$. For the perturbative distribution one finds that $`\overline{z}0.5\mathrm{}0.6`$. This is clearly seen by comparing pions at $`p_{}=1.7`$ GeV and gluons at $`p_{}=3`$ GeV. As the slope of the gluon distribution decreases, $`\overline{z}`$ moves towards smaller values as well. This causes the gluons in saturation domain to fragment mainly to much smaller values of $`p_{}`$ than the perturbative gluons would. The situation improves at LHC energies where the saturation domain extends already at midrapidity to $`2`$ GeV. At large rapidities it is difficult to distinguish between saturation and perturbative spectra, and presumably different dynamics which are not contained in either of these models become dominant towards the fragmentation regions.

We have investigated the quantitative behavior of the $`p_{}`$\- and $`y`$-distributions of gluons and pions in collisions of nuclei with a large atomic number asymmetry. We have shown that within the color glass condensate framework one indeed reproduces the behavior of the conventional pQCD at high $`p_{}`$ which significantly deviates from commonly cited expectation of $`1/p_{}^4`$ in the literature on saturation. The fragmentation of gluons into pions was shown to significantly change the qualitative shape of the spectrum. However, a clear difference from the spectrum obtained from the conventional pQCD calculation in the form of a rapidity dependent pion deficit can still be extracted by measuring the rapidity dependence of the $`p_{}`$-distributions. The rapidity dependence of the slopes were found to be different at small $`p_{}`$, and a change in the form of the rapidity distribution of pions as a function of $`p_{}`$ could be used to determine the onset of saturation.

Acknowledgements: We would like to thank A. Dumitru, K.J. Eskola, K. Kajantie, J. Jalilian–Marian, L. McLerran and W. Vogelsang for useful discussions. J.T.L. thanks the Nuclear Theory Group at Brookhaven National Laboratory, where part of this work was done, for hospitality.