text:
# HTL quasiparticle models of deconfined QCD at finite chemical potential

## I Introduction

Asymptotic freedom suggests that at sufficiently high temperature and/or quark chemical potential QCD is deconfined, i.e. can be described in terms of the fundamental quark and gluonic degrees of freedom Collins:1975ky . High temperature and small chemical potential are of relevance to the quest for the quark-gluon plasma in heavy-ion experiments; sufficiently high chemical potential (at comparatively low temperature) may be reached in the cores of compact stars. In the latter case, novel color superconducting phases may occur Alford:2002wf , which should however have only comparatively small effects on the equation of state.

The thermodynamic pressure of QCD at high temperature $`T\mathrm{\Lambda }_{\mathrm{QCD}}`$ has been determined in perturbation theory to order $`\alpha _s^{5/2}`$ in QCD Arnold:1995eb ; Zhai:1995ac  and even to $`\alpha ^3\mathrm{log}\alpha `$ in a recent heroic effort Kajantie:2002wa . Strict perturbation theory, however, shows extremely poor convergence for any temperature of practical interest so that further resummations appear to be necessary. In recent years, resummations which are based on the hard-thermal-loop (HTL) effective action Braaten:1992gm  have been proposed, alternatively in form of so-called HTL perturbation theory Andersen:1999fw ; Andersen:2002ey  or based upon the 2-loop $`\mathrm{\Phi }`$-derivable approximation Blaizot:1999ip ; Blaizot:2000fc ; Blaizot:2003tw  (see also Peshier Peshier:2000hx ). The latter approach, which assumes weakly interacting quasiparticles as determined by the HTL propagators and NLO corrections thereof, leads to results which agree remarkably well with lattice data down to about $`3T_c`$ if standard 2-loop running $`\alpha _s`$ is adopted,111This is in fact not the case for 2-loop HTL perturbation theory Andersen:2002ey , but it appears that the source of the difficulty has to do with the necessity of thermal counter-terms and incompletely compensated hard contributions, see Ref. Blaizot:2003iq . and it has been recently shown that these results are in fact consistent with those obtained in high-order perturbation theory if the latter is organized through effective (dimensionally reduced) field theory according to Braaten:1996jr  and effective-field-theory parameters are kept without further expanding in powers of the coupling Blaizot:2003iq , as already advocated in Ref. Kajantie:2002wa .

The behavior of the thermodynamic potentials closer to the transition temperature needs to be studied by fully non-perturbative means. For zero quark chemical potential $`\mu =0`$ lattice QCD calculations are seemingly up to the task of determining thermodynamic quantities for a quark-gluon plasma very accurately by now Karsch:2001cy ; Laermann:2003cv . Although recently there has been important progress also for non-vanishing $`\mu `$ Fodor:2001au ; deForcrand:2002ci ; Fodor:2002km ; Allton:2002zi ; D'Elia:2002gd ; Gavai:2003mf , the case of large $`\mu T`$ (relevant for the equation of state for cold dense matter, which is of importance in astrophysical situations Peshier:1999ww ; Peshier:2002ww ; Fraga:2001xc ; Andersen:2002jz ) is currently beyond the reach of lattice calculations.

As a remedy for this situation, Peshier *et al.* Peshier:1999ww ; Peshier:2002ww  proposed a method which can be used to map the available lattice data for $`\mu =0`$ to finite $`\mu `$ and small temperatures by describing the interacting plasma as a system of massive quasiparticles (QPs); section 2 gives a short review of this technique. While in their simple QP model the (thermal) masses of the quarks and gluons are approximated by the asymptotic limit of the hard-thermal-loop (HTL) self-energies, we consider models based on the HTL-resummed entropy Blaizot:1999ip ; Blaizot:2000fc ; Blaizot:2003tw , which include more of the physically important plasmon effect than in the simple QP model, and (NLO) extensions thereof, which include the full plasmon effect. In perturbation theory it is essentially this effect that spoils convergence, so a priori it was not granted that a model including this effect fully would give physically sensible results. As it turns out, however, the results based on our models (which are presented in section 3 and 4, respectively) turn out to be well-behaved and indeed lie not too far from the simple QP model.

In section 5 we give the 2-parameter fits of the various models to the lattice data AliKhan:2001ek  for 2 flavors and compare the model results for the pressure at small temperatures and large chemical potential in section 6. Furthermore, we compare our results for the isobar emerging at222Here and in the following $`T_c`$ is always to be understood as $`T_c|_{\mu =0}`$. $`T_c`$ to recent lattice results and give our conclusions in section 7.

## II Mapping lattice data to finite chemical potential using QP models

Assuming an SU($`N`$) plasma of gluons and $`N_f`$ light quarks in thermodynamic equilibrium can be described as a weakly interacting gas of massive quasiparticles with residual interaction $`B`$, the pressure of the system is given by Peshier:1999ww

$$P(T,\mu )=\underset{i=g,q}{}p_i(T,\mu _i,m_i^2)B(m_g,m_q),$$
(1)
where the sum runs over gluons (g), quarks and anti-quarks (q) with respective chemical potential $`0,\pm \mu `$; $`p_i`$ are the model dependent QP pressures and $`m_i`$ are the QP masses. The latter are functions of an effective strong coupling $`G^2(T,\mu )`$, which appears only within $`m_i(G^2,T,\mu )`$ in a predefined form.

Using the stationarity of the thermodynamic potential under variation of the self-energies and Maxwell’s relations Peshier:1999ww , one obtains a partial differential equation for $`G^2`$,

$$a_T\frac{G^2}{T}+a_\mu \frac{G^2}{\mu }=b,$$
(2)
where $`a_T`$, $`a_\mu `$ and $`b`$ are coefficients that are given by integrals depending on $`T,\mu `$ and $`G^2`$. Given a valid boundary condition, a solution for $`G^2(T,\mu )`$ is found by solving the above flow equation by the method of characteristics; once $`G^2`$ is thus known in the $`T,\mu `$ plane, the QP pressure is fixed completely. The residual interaction $`B`$ is then given by the integral

$$B=\underset{i}{}\frac{p_i}{m_i^2}\left(\frac{m_i^2}{\mu }d\mu +\frac{m_i^2}{T}dT\right)+B_0,$$
(3)
where $`B_0`$ is an integration constant that has to be fixed by lattice data (usually by requiring $`P(T_c,\mu =0)=P_{lattice}(T_c)`$). Motivated by the fact that at $`\mu =0`$ and $`TT_c`$ the coupling should behave as predicted by the perturbative QCD beta-function, Ref. Peshier:1999ww  used the ansatz

$$G^2(T,0)=\frac{48\pi ^2}{(11N2N_f)\mathrm{ln}\frac{T+T_s}{T_c}\lambda }.$$
(4)
The parameters $`\lambda `$ and $`T_s`$ are determined by fitting the entropy of the model (which is independent of $`B`$) to available lattice data at $`\mu =0`$. Using (4) as boundary condition for (2), the above procedure allows one to map the lattice pressure from $`\mu =0`$ to the whole $`T,\mu `$-plane.

## III Simple and HTL QP model

A simple ansatz for the QP pressure is obtained by taking $`p_i`$ to be the pressure of a free gas of massive particles Peshier:1999ww ; Peshier:2002ww , with the masses defined as the asymptotic (i.e. large momentum) gluonic and fermionic self-energies, respectively Flechsig:1996ju :

$`\widehat{m}_{\mathrm{}}^2`$ $`=`$ $`(2N+N_f){\displaystyle \frac{G^2T^2}{12}}+N_f{\displaystyle \frac{G^2\mu ^2}{4\pi ^2}},`$ (5)
$`\widehat{M}_{\mathrm{}}^2`$ $`=`$ $`{\displaystyle \frac{G^2(N^21)}{8N}}\left(T^2+{\displaystyle \frac{\mu ^2}{\pi ^2}}\right).`$ (6)
Comparing with perturbation theory one finds that while the Stefan-Boltzmann and leading-order interaction terms are correctly reproduced, only $`1/(4\sqrt{2})`$ of the NLO term of the interaction pressure (the plasmon effect $`G^3T^4`$) is present Blaizot:2000fc .

To include more of this physically relevant effect, one can make use of the HTL-resummed entropy Blaizot:1999ip ; Blaizot:2000fc ; Blaizot:2003tw , which takes into account the momentum dependence of the QP excitations as well as Landau-damping effects. The HTL model ansatz for the QP pressure for gluons ($`p_g`$) and quarks $`p_q`$ then reads

$`p_g`$ $`=`$ $`d_g{\displaystyle }{\displaystyle \frac{d^3k}{(2\pi )^3}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{d\omega }{2\pi }}n(\omega )[2\mathrm{I}\mathrm{m}\mathrm{ln}(\omega ^2+\mathrm{k}^2+\widehat{\mathrm{\Pi }}_\mathrm{T})2\mathrm{I}\mathrm{m}\widehat{\mathrm{\Pi }}_\mathrm{T}\mathrm{Re}\widehat{\mathrm{D}}_\mathrm{T}`$ (7)
$`+\mathrm{Im}\mathrm{ln}(\mathrm{k}^2+\widehat{\mathrm{\Pi }}_\mathrm{L})+\mathrm{Im}\widehat{\mathrm{\Pi }}_\mathrm{L}\mathrm{Re}\widehat{\mathrm{D}}_\mathrm{L}]`$
$`p_q`$ $`=`$ $`d_q{\displaystyle }{\displaystyle \frac{d^3k}{(2\pi )^3}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{d\omega }{2\pi }}(f_+(\omega )+f_{}(\omega ))[\mathrm{Im}\mathrm{ln}(\mathrm{k}\omega +\widehat{\mathrm{\Sigma }}_+)`$ (8)
$`\mathrm{Im}\widehat{\mathrm{\Sigma }}_+\mathrm{Re}\widehat{\mathrm{\Delta }}_++\mathrm{Im}\mathrm{ln}(\mathrm{k}+\omega +\widehat{\mathrm{\Sigma }}_{})+\mathrm{Im}\widehat{\mathrm{\Sigma }}_{}\mathrm{Re}\widehat{\mathrm{\Delta }}_{}],`$
where $`d_g=2(N^21)`$, $`d_q=2NN_f`$ for gluons and quarks/anti-quarks, respectively; $`n(\omega )`$ and $`f_\pm (\omega )`$ are the bosonic and fermionic distribution functions. $`\widehat{D}_{T,L}`$, $`\widehat{\mathrm{\Delta }}_\pm `$ are the HTL propagators with $`\widehat{\mathrm{\Pi }}_{T,L}`$ and $`\widehat{\mathrm{\Sigma }}_\pm `$ the corresponding self-energies,

$`\widehat{\mathrm{\Pi }}_L(\omega ,k)`$ $`=`$ $`\widehat{m}_D^2\left[1{\displaystyle \frac{\omega }{2k}}\mathrm{ln}{\displaystyle \frac{\omega +k}{\omega k}}\right]`$ (9)
$`\widehat{\mathrm{\Pi }}_T(\omega ,k)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[\widehat{m}_D^2+{\displaystyle \frac{\omega ^2k^2}{k^2}}\widehat{\mathrm{\Pi }_L}\right]`$ (10)
$`\widehat{\mathrm{\Sigma }}_\pm `$ $`=`$ $`{\displaystyle \frac{\widehat{M}^2}{k}}\left[1{\displaystyle \frac{\omega k}{2k}}\mathrm{ln}{\displaystyle \frac{\omega +k}{\omega k}}\right];`$ (11)
which we have written in terms of the HTL Debye mass squared $`\widehat{m}_D^2=2\widehat{m}_{\mathrm{}}^2`$ and the HTL zero-momentum quark mass squared $`\widehat{M}^2=\widehat{M}_{\mathrm{}}^2/2`$.

The HTL QP model includes $`1/4`$ of the full plasmon effect at $`\mu =0`$, which is a factor of $`\sqrt{2}`$ more than in the simple QP model which only uses the asymptotic HTL masses. The remainder of the plasmon effect is in fact entirely due to NLO corrections to these asymptotic masses at order $`G^3T^2`$ Blaizot:1999ip ; Blaizot:2000fc ; Blaizot:2003tw .

## IV NLO models

The complete (momentum-dependent) NLO corrections to the asymptotic masses of quarks and gluons have not been calculated yet, but as far as the plasmon effect is concerned, these corrections contribute in the averaged form Blaizot:1999ip ; Blaizot:2000fc

$$\overline{\delta }m_{\mathrm{}}^2=\frac{𝑑kkn^{}(k)\mathrm{Re}\delta \mathrm{\Pi }_T(\omega =k)}{𝑑kkn^{}(k)}=\frac{1}{2\pi }G^2NT\widehat{m}_D$$
(12)
and similarly

$$\overline{\delta }M_{\mathrm{}}^2=\frac{𝑑kk(f_+^{}(k)+f_{}^{}(k))\mathrm{Re}\mathrm{\hspace{0.17em}2}k\delta \mathrm{\Sigma }_+(\omega =k)}{𝑑kk(f_+^{}(k)+f_{}^{}(k))}=\frac{(N^21)/(2N)}{2\pi }G^2T\widehat{m}_D.$$
(13)
For the values of the coupling $`G`$ considered here, these corrections are so large that they give tachyonic masses when treated strictly perturbatively. In Ref. Blaizot:2000fc  it has been proposed to incorporate these corrections through a quadratic gap equation which works well as an approximation in the exactly solvable scalar O($`N\mathrm{}`$)-model, where strict perturbation theory would lead to identical difficulties. However, for the fermionic asymptotic masses, in order to have the correct scaling of Casimir factors in the exactly solvable large-$`N_f`$ limit of QCD Moore:2002md , a corresponding gap equation has to remain linear in the fermionic mass squared. Choosing the correction therein to be determined by the solution to the gluonic gap leads to Rebhan:2003fj ; Blaizot:2003tw

$`\overline{m}_{\mathrm{}}^2`$ $`=`$ $`\widehat{m}_{\mathrm{}}^2{\displaystyle \frac{G^2NT}{\sqrt{2}\pi }}\overline{m}_{\mathrm{}}`$ (14)
$`\overline{M}_{\mathrm{}}^2`$ $`=`$ $`\widehat{M}_{\mathrm{}}^2{\displaystyle \frac{G^2(N^21)T}{2\sqrt{2}\pi N}}\overline{m}_{\mathrm{}},`$ (15)
where $`\widehat{m}_{\mathrm{}}^2`$ and $`\widehat{M}_{\mathrm{}}^2`$ are the leading-order gluonic and fermionic asymptotic masses as given in (5) and (6). This in fact avoids tachyonic masses for the fermions as long as $`N_f3`$. 333For $`N_f=3`$ the solutions to the above approximate gap equations happen to coincide with those obtained in the original version of two independent quadratic gap equations of Ref. Blaizot:2000fc . For the case $`N_f=2`$ considered in this paper, the differences are fairly small. For $`N_f>3`$, however, the necessity to avoid tachyonic masses would restrict the range of permissible coupling strength $`G`$.

Finally, since the averaged quantities $`\overline{m}_{\mathrm{}}^2`$ and $`\overline{M}_{\mathrm{}}^2`$ are the effective masses at hard momenta only, a cutoff scale $`\mathrm{\Lambda }=\sqrt{2\pi T\widehat{m}_Dc_\mathrm{\Lambda }}`$ is introduced that separates soft from hard momenta. The QP pressure for this model, which in the following will be referred to as NLA-model, then separates into a soft and a hard component for both gluons and fermions. The soft parts are given by expressions like (7) and (8), but with $`\mathrm{\Lambda }`$ as upper limit for the momentum integration. For the hard parts, the momentum integrations run from $`\mathrm{\Lambda }`$ to $`\mathrm{}`$ and the mass pre-factors in the HTL self-energies (9,10,11) are replaced by their asymptotic counterparts, $`\widehat{m}_D^22\overline{m}_{\mathrm{}}^2`$ and $`\widehat{M}^2\frac{1}{2}\overline{M}_{\mathrm{}}^2`$.

The single free parameter $`c_\mathrm{\Lambda }`$ in $`\mathrm{\Lambda }`$ can be varied around $`1`$ to obtain an idea of the “theoretical error” of the model. In the following we will consider the range $`c_\mathrm{\Lambda }=\frac{1}{4}`$ to $`c_\mathrm{\Lambda }=4`$; note that $`c_\mathrm{\Lambda }=\mathrm{}`$ corresponds to the HTL-model (which has been used as cross-check) since all hard corrections are ignored. On the other hand, $`c_\mathrm{\Lambda }=0`$ would assume that (14) and (15) represent good approximations for the NLO corrections to the spectral properties of soft excitations. However, the few existing results, in particular on NLO corrections to the Debye mass Rebhan:1993az  and the plasma frequency Schulz:1994gf , appear to be rather different so that it seems safer to leave the soft sector unchanged by keeping a finite $`c_\mathrm{\Lambda }`$.

## V Entropy and pressure for $`\mu =0`$

To obtain the input parameters $`T_s`$ and $`\lambda `$ we fitted the entropy expressions from the models under consideration to lattice data AliKhan:2001ek  for $`N_f=2`$, with an estimated continuum extrapolation as used in Ref. Peshier:2002ww . We find

|  | HTL | $`c_\mathrm{\Lambda }=4`$ | $`c_\mathrm{\Lambda }=1`$ | $`c_\mathrm{\Lambda }=1/4`$ |
| --- | --- | --- | --- | --- |
| $`T_s/T_c`$ | -0.89 | -0.89 | -0.84 | -0.61 |
| $`\lambda `$ | 19.4 | 18.64 | 11.43 | 3.43 |

where it can be seen that the NLA model with $`c_\mathrm{\Lambda }=4`$ is very close to the HTL model. The fits to the entropy data (shown in Fig. 2) all lie in a narrow band for all the models considered. The fitted effective coupling $`G^2`$ for the various models is shown in figure 2; for comparison, also the 2-loop perturbative running coupling in $`\overline{\text{MS}}`$ is shown, where the renormalization scale is varied between $`\pi T`$ and $`4\pi T`$ and following Ref. Gupta:2000hr  we have chosen $`T_c=0.49\mathrm{\Lambda }_{\overline{\text{MS}}}`$. As can be seen from the plot, the results for the effective coupling are well within the range of the 2-loop perturbative running coupling (for the case $`c_\mathrm{\Lambda }=1/4`$ and renormalization scale $`\pi T`$ the results even seem to be identical, which is, however, probably only a coincidence). We also show the result for the coupling obtained in the semiclassical approach of Ref. Schneider:2003uz 444For finite $`\mu `$ and constant temperatures, however, the coupling obtained in Schneider:2003uz  rises while our results indicate a decrease of the coupling (which is consistent with the standard QCD running coupling with renormalization scale proportional to $`\sqrt{T^2+(\mu /\pi )^2}`$)..

In general one can see that the effective coupling becomes bigger when $`c_\mathrm{\Lambda }`$ gets smaller; this is because the hard masses (for equal values of the coupling) are smaller than the soft masses which makes the entropy increase when the hard parts become more important. Accordingly, the coupling has to rise in order for the entropy to match the data (therefore, for the extreme case $`c_\mathrm{\Lambda }=0`$ one finds huge values of the effective coupling constant).

Once the effective coupling is known for $`\mu =0`$ one can proceed to evaluate the QP pressure $`p_i`$ and the residual interaction $`B`$ from Eq. (3), shown in Fig. 4; in this figure, the HTL results (not shown) would lie marginally above the upper boundary of the band (cf. Romatschke:2002pb  for a comparison between simple QP and HTL model). A comparison between the full pressure (1) and the lattice data is shown in figure 4, where the integration constant $`B_0`$ was set so that $`P(T_c)=0.536(1)`$:

|  | HTL | $`c_\mathrm{\Lambda }=4`$ | $`c_\mathrm{\Lambda }=1`$ | $`c_\mathrm{\Lambda }=1/4`$ |
| --- | --- | --- | --- | --- |
| $`B_0/T_c^4`$ | 0.82 | 0.73 | 0.6 | 0.47 |

## VI Finite chemical potential

A calculation and subsequent evaluation of the coefficients $`a_\mu `$, $`a_T`$ and $`b`$ in (2) shows that $`b`$, in contrast to $`a_\mu `$ and $`a_T`$, is strongly model-dependent. Solving the flow equation (2) and integrating $`B`$ along the characteristics using Eq. (3) one obtains pressure and effective coupling in the whole $`T,\mu `$ plane. Here we focus on the pressure but the calculation of the other thermodynamic quantities is straightforward once the effective coupling and the pressure have been obtained.

Concerning the shape of the characteristics, it has already been noticed in Romatschke:2002pb  that the crossing of characteristics in the simple QP model Peshier:1999ww  does not occur in the HTL-model; this feature is preserved in the NLA models.

### VI.1 Lines of constant pressure

It is straightforward to extract from our data the line where $`P(T,\mu )=P(T_c,0)`$; for $`N_f=2`$, the results for the simple and HTL QP model as well as for the NLA models for $`c_\mathrm{\Lambda }>1`$ are rather close to each other. For NLA models with $`c_\mathrm{\Lambda }<1`$, however, the slope of the constant pressure line (until $`\mu /T_c=1`$) deviates rapidly from the values found for $`c_\mathrm{\Lambda }>1`$:

|  | HTL | $`c_\mathrm{\Lambda }=4`$ | $`c_\mathrm{\Lambda }=1`$ | $`c_\mathrm{\Lambda }=1/4`$ |
| --- | --- | --- | --- | --- |
| $`T_c\frac{dT}{d\mu ^2}`$ | -0.06818(8) | -0.06810(6) | -0.06329(34) | -0.041(9) |

We have also compared these results with recent lattice data for 2+1 flavors Fodor:2002km  in figure 6. The variations of the lattice data in this figure result from fitting the data with a 2nd order (upper limit) and 4th order polynomial (lower limit) as used in Szabo:2003kg . The comparison shows that our 2-flavor results obtained from simple, HTL and NLA QP models with $`c_\mathrm{\Lambda }>1`$ (indicated in figure 6 by the two lower lying pointed lines) are in very good agreement with the lattice data whereas for NLA models with $`c_\mathrm{\Lambda }<1`$ the slope is much flatter. Also shown are the lattice results for the constant pressure line from Allton:2002zi  for 2 flavors, with a slope of $`T_c\frac{dT}{d\mu ^2}=0.107(22)`$, which is significantly steeper.

### VI.2 Susceptibilities

It been noticed in Fodor:2002km  that the quantity

$$\frac{\mathrm{\Delta }P}{\mathrm{\Delta }P_{SB}}=\frac{P(T,\mu )P(T,0)}{P_{SB}(T,\mu )P_{SB}(T,0)}$$
(16)
is essentially $`\mu `$-independent; we recover a similar scaling behavior for our models. As expected, for small chemical potential this curve is very close to the quark-number susceptibilities $`\chi (T)/\chi _0(T)`$ obtained at $`\mu =0`$. In fact, we found that also for larger $`\mu `$ the pressure can be well approximated by

$$P(T,\mu )P(T,0)=\frac{\chi (T)}{2}(\mu ^2+\frac{\mu ^4}{2\pi ^2T^2}).$$
(17)
To quantify this assertion we consider the subtracted quantity $`\delta P=P(T,\mu )P(T,0)\frac{\chi (T)}{2}(\mu ^2+\frac{\mu ^4}{2\pi ^2T^2})`$, which should be zero, if the pressure at finite chemical potential was fully determined by the susceptibilities at $`\mu =0`$ according to (17); the results (scaled with $`T^4`$) are shown in figure 6 for $`\mu /T_c1`$. As can be seen, dropping the $`\mu ^4`$ term in Eq. (17) (open symbols) deteriorates the result considerably. Trying to determine the higher-order susceptibility, $`\overline{\chi }=\frac{^4P}{\mu ^4}`$ which we parameterized as

$$\overline{\chi }=\frac{\chi (T)}{2T^2}\frac{4!}{x(T)^2},$$

we find that – within our numerical errors – the result for $`x(T)`$ is consistent with $`\sqrt{2}\pi =4.44\mathrm{}`$ for $`T>1.5T_c`$ (which was our original assumption in Eq. (17)); therefore, our results are close to but about $`10`$ percent lower than the values found in Ref. Gavai:2003mf .

Clearly, $`\delta P`$ gets smaller with increasing temperature, but also at relatively low temperatures and high chemical potential, $`\delta P`$ is remarkably small (e.g., for the HTL model for $`T=1.5T_c`$ and $`\mu /T_c4`$ we find a $`\delta P/T^4`$ smaller than $`10`$ percent).

In Fig. 8 we compare the results for the susceptibilities of our NLA models with the lattice data for 2 flavors from Ref. Gavai:2001ie  and the “scaling curve” in Fodor:2002km , which provides an approximation to the 2+1 flavor quark-number susceptibility at $`\mu =0`$. We find that all these results agree very well, with the agreement getting better at higher temperatures.

It should be noted however that none of the lattice data we are using has yet been rigorously extrapolated to the continuum limit, so that both the data and the extrapolation by means of QP models are still likely to change somewhat when this will be done eventually.

### VI.3 Very small temperatures

Extending the lines of constant pressure from our models to very small temperatures we obtain a crude estimate of the phase transition line in this region of phase space. Denoting with $`\mu _c`$ the chemical potential where (for vanishing temperature) the pressure equals the lattice pressure at $`\mu =0`$, $`P(0,\mu _c)=P(T_c,0)`$, we find for our models (assuming $`T_c=172`$ MeV)

|  | HTL | $`c_\mathrm{\Lambda }=4`$ | $`c_\mathrm{\Lambda }=1`$ | $`c_\mathrm{\Lambda }=1/4`$ |
| --- | --- | --- | --- | --- |
| $`\mu _c`$ | 533 MeV | 536 MeV | 558 MeV | 584 MeV |
| $`\mu _0`$ | 509 MeV | 511 MeV | 537 MeV | 567 MeV |

Here we have also given the values $`\mu _0`$ where the pressure vanishes, which may be taken as a definite lower bound for the critical chemical potential within the respective models.

In general, these results are in agreement with the estimates for $`\mu _c`$ from Szabo:2003kg ; Peshier:2002ww ; for the QP models considered, the lowest and highest results for $`\mu _c`$ are obtained for the HTL model and the NLA model with $`c_\mathrm{\Lambda }=1/4`$, respectively, while the $`\mu _c`$ of the simple QP model lies between the NLA $`c_\mathrm{\Lambda }=4`$ and $`c_\mathrm{\Lambda }=1`$ model values. However, even our lowest result for $`\mu _c`$ turns out to exceed the value for the critical chemical potential expected in Ref. Fraga:2001id ; Fraga:2001xc .

A comparison of the perturbative pressure Freedman:1977ub  at vanishing temperature with the results for our models at $`T/T_c0.01`$ is shown in figure 8. We plotted the perturbative results converted to $`\overline{\text{MS}}`$ with the standard 2-loop running coupling, $`\mathrm{\Lambda }_{\overline{\text{MS}}}=T_c/0.49`$ Gupta:2000hr  and renormalization scale varied from $`\mu `$ to $`3\mu `$; a comparison of this coupling and our results is shown in figure 9.

### VI.4 Equation of state for cold dense matter

By calculating the number density at $`T/T_c0.01`$ and using our results for the pressure we obtain an equation of state for cold deconfined matter. As is the case for the simple QP model Peshier:2002ww , we find that the energy density $``$ is well fitted by the linear relation

$$(p)=4\stackrel{~}{B}+\alpha p$$

with

|  | HTL | $`c_\mathrm{\Lambda }=4`$ | $`c_\mathrm{\Lambda }=1`$ | $`c_\mathrm{\Lambda }=1/4`$ |
| --- | --- | --- | --- | --- |
| $`4\stackrel{~}{B}/T_c^4`$ | 11.1(8) | 12.3(8) | 14.7(9) | 19.2(1.6) |
| $`\alpha `$ | 3.23(5) | 3.22(4) | 3.22(4) | 3.17(4) |

A value of $`\alpha 3.2`$ for $`N_f=2`$ has also been found in the simple quasiparticle model in Peshier:2002ww  and thus seems to be model independent, in contrast to the bag constant $`\stackrel{~}{B}^{1/4}`$, which varies between $`314`$ and $`360`$ MeV. By using the Tolman-Oppenheimer-Volkov equations and the equation of state one can determine the mass-radius relations of non-rotating quark-stars Andersen:2002jz ; choosing $`\alpha =3.2`$ and taking the bag constant values from above, we find stable star configurations with radii ranging from $`3.6`$ to $`4.9`$ km and masses from $`0.7`$-$`0.8`$ solar masses for NLA, $`c_\mathrm{\Lambda }=1/4`$ and HTL models, respectively. However, it should be kept in mind that the outermost layers of such a quark star are metastable with respect to hadronic matter and therefore the details of the star structure will depend sensitively on the hadronic equation of state Peshier:2002ww .

## VII Conclusions and Outlook

We considered an improvement of simple quasiparticle models Peshier:1999ww  by using the full HTL approximation and certain NLO corrections thereof along the lines of Ref. Blaizot:1999ip ; Blaizot:2000fc ; Blaizot:2003tw  and investigated the respective predictions at finite chemical potential when these models are matched to $`N_f=2`$ lattice data at $`\mu =0`$.

We found that the slope of constant pressure for those models with $`c_\mathrm{\Lambda }>1`$ agrees very well with recent lattice data for 2+1 flavors Fodor:2002km ; Szabo:2003kg , while the lattice data for $`N_f=2`$ Allton:2002zi  indicate significantly steeper slopes.

We also found that for our models the pressure for non-vanishing $`\mu `$ scales with the quark number susceptibilities at $`\mu =0`$, provided that the $`\mu ^4`$ terms (corresponding to a higher order susceptibility) are not dropped but added with the same weight as in the ideal gas limit. We therefore interpreted the “scaling curve” from 2+1 flavor lattice data Fodor:2002km  as susceptibilities at $`\mu =0`$ and found that a comparison with $`N_f=2`$ susceptibilities from another lattice group Gavai:2001ie  as well as those obtained from our models all agree very well.

We finally extended our results to the small temperature, high density regime and obtained an estimate for the critical density at $`T=0`$ which is consistent with earlier results, though in excess of the expectations of Ref. Fraga:2001id ; Fraga:2001xc . Furthermore, we obtained an equation of state for cold, dense matter, which allows for pure quarks stars with masses of $`0.8M_{}`$ and radii of less than $`5`$ km, similar to the results obtained by Fraga:2001xc ; Peshier:2002ww ; Blaschke:1998hy .

###### Acknowledgements.

We want to thank Z. Fodor and K.K. Szabo for kindly providing their lattice data and A. Peshier, C. Schmidt and R.A. Schneider for discussions. This work has been supported by the Austrian Science Foundation FWF, project no. 14632.


format: text

text:
# Evaluation of the two-photon exchange diagrams for the (𝟏𝐬)^𝟐𝟐𝐩_{𝟑/𝟐} electron configuration in Li-like ions

## I Introduction

At present the lowest-lying states in heavy Li-like ions can be investigated very precisely both theoretically and experimentally. One of the most precise experimental results in these systems has been obtained by Beiersdorfer and co-workers beiersdorfer98  for the $`2p_{3/2}`$-$`2s`$ transition energy in Li-like bismuth, which was determined with an accuracy of $`0.04`$ eV. Accurate experimental data are available at present for a number of other elements as well. For the latest high-precision measurements we refer to Refs. bosselmann99 ; feili00 ; brandau00 ; the outline of earlier investigations can be found in Ref. bosselmann99 .

The accuracy reached in experimental investigations provides a promising tool for probing QED corrections in the strong Coulomb field of the nucleus up to second order in the fine structure constant $`\alpha `$. For the $`2p_{1/2}`$-$`2s`$ transition, this project has been carried out in a series of our previous investigation artemyev99 ; yerokhin99 ; yerokhin00 ; yerokhin01pra2p . In Ref. yerokhin01pra2p  we completed the evaluation of all two-electron QED corrections of second order in $`\alpha `$ and obtained most accurate theoretical predictions for the $`2p_{1/2}`$-$`2s`$ splitting within a wide range of nuclear charge numbers $`Z`$. Based on a careful estimate of the uncertainty of the theoretical values, we concluded that already now the comparison of theory and experiment for Li-like uranium provides a test of QED effects of second order in $`\alpha `$ at the level of accuracy of about 17%. For the $`2p_{3/2}`$-$`2s`$ and $`2p_{1/2}`$-$`2s`$ transitions in Li-like bismuth, analogous calculations have been performed recently by Sapirstein and Cheng sapirstein01 . However, in order to match with the experimental accuracy for the $`2p_{3/2}`$-$`2s`$ splitting, rigorous evaluations of second-order QED corrections are required also for other ions than bismuth. The first step in this direction has been performed in our earlier investigation artemyev99  where we have evaluated the vacuum-polarization screening correction for several energy levels of Li-like ions, including the $`(1s)^22p_{3/2}`$ state. The aim of the present work is to calculate the two-photon exchange correction for this state (for extensive calculations of these corrections for the lower states in Li-like ions and for non-mixed low-lying states in He-like ions we refer the reader to Refs. yerokhin00 ; yerokhin01pra2p ; blundell93b ; lindgren95pra ; moh00 ; and01 ; and02 ). After all that, the self-energy screening correction remains the last uncalculated two-electron second-order QED contribution for this state.

This paper is organized as follows. In the next section we present the basic formulas for the two-photon exchange correction for the $`(1s)^22p_{3/2}`$ state. The description of our numerical procedure is given in Sec. III, and the results obtained are discussed in Sec. IV.

Relativistic units ($`\mathrm{}=c=1`$) are used throughout this paper.

## II Basic formulas

The detailed derivation of the two-photon exchange corrections to the $`(1s)^22s`$ and $`(1s)^22p_{1/2}`$ states of Li-like ions can be found in our previous paper yerokhin01pra2p . For the $`(1s)^22p_{3/2}`$ state the derivation is performed along the same lines. Thus, we present mainly the final formulas here. Our derivation is based on the two-time Green function (TTGF) method shabaev90ivf ; shabaev94ttg2 . For the detailed description of the method we refer to the recent review shabaev02rep .

The two-photon exchange corrections to the $`(1s)^22p_{3/2}`$ state of the Li-like ions can be conveniently separated in three parts: the two-photon exchange contribution due to the interaction between two $`1s`$ electrons, the two-photon exchange contribution due to the interaction between the valence electron and one of the $`1s`$ electrons, and the three-electron contribution. The first part coincides with the two-photon exchange correction to the ground-state energy of He-like ions. Its calculation was carried out in blundell93b ; lindgren95pra . This correction does not contribute to the $`2p`$-$`2s`$ splitting in Li-like ions and is not considered here. The remaining two-electron and three-electron corrections are diagrammatically depicted in Fig. 1.

We start from the expression for the second-order correction to the energy shift of the level $`k`$ shabaev02rep ,

$`\mathrm{\Delta }E_k^{(2)}`$ $`=`$ $`{\displaystyle \frac{1}{2\pi i}}{\displaystyle _\mathrm{\Gamma }}𝑑E\mathrm{\Delta }E\mathrm{\Delta }g_{kk}^{(2)}(E)`$ (1)
$``$ $`{\displaystyle \frac{1}{2\pi i}}{\displaystyle _\mathrm{\Gamma }}𝑑E\mathrm{\Delta }E\mathrm{\Delta }g_{kk}^{(1)}(E){\displaystyle \frac{1}{2\pi i}}{\displaystyle _\mathrm{\Gamma }}𝑑E^{}\mathrm{\Delta }g_{kk}^{(1)}(E^{}),`$
where $`\mathrm{\Delta }g_{kk}(E)=g_{kk}(E)g_{kk}^{(0)}(E)`$, $`g_{kk}(E)=u_k|g(E)|u_k`$, $`u_k`$ is the unperturbed wave function, $`\mathrm{\Delta }E=EE_k^{(0)}`$, $`E_k^{(0)}`$ is the unperturbed energy of the state $`k`$, and $`g_{kk}^{(0)}(E)=(\mathrm{\Delta }E)^1`$ is the function $`g_{kk}(E)`$ in the zeroth-order approximation. The function $`g(E)g(E,𝐱_1^{},\mathrm{},𝐱_N^{};𝐱_1,\mathrm{},𝐱_N)`$ is the temporal Fourier transform of the N-electron two-time Green function. Its definition and the corresponding Feynman rules can be found in shabaev02rep . The superscripts in Eq. (1) indicate the order of the contribution in $`\alpha `$.

For the two-photon exchange correction, Feynman diagrams contributing to $`\mathrm{\Delta }g^{(2)}(E)`$ are presented in Fig. 1. We refer to the corresponding contributions as the ladder (a), the crossed (b), and the three-electron (c) terms. The second term in Eq. (1) is known as the disconnected contribution. It vanishes completely when considered together with the reducible contribution (for details, see yerokhin01pra2p ). In our case, the unperturbed wave function is

$`u_k={\displaystyle \frac{1}{\sqrt{3!}}}{\displaystyle \underset{P}{}}(1)^P\psi _{Pa}(𝐱_1)\psi _{Pb}(𝐱_2)\psi _{Pv}(𝐱_3),`$ (2)
where $`v`$ denotes the valence electron, $`a`$ and $`b`$ are the electrons in $`(1s)^2`$ shell, and $`P`$ is the permutation operator (in the factor $`(1)^P`$, the parity of the permutation is implied by $`P`$). For brevity we will use also the following notations:

$`I(\omega )`$ $`=`$ $`e^2\alpha ^\mu \alpha ^\nu D_{\mu \nu }(\omega ),`$ (3)
$`I_{abcd}(\omega )`$ $`=`$ $`ab|I(\omega )|cd,`$ (4)
$`I_{ab;cd}`$ $`=`$ $`I_{abcd}(\mathrm{\Delta }_{bd})I_{bacd}(\mathrm{\Delta }_{ad}),`$ (5)
$`I^{}(\omega )`$ $`=`$ $`{\displaystyle \frac{dI(\omega )}{d\omega }},`$ (6)
where $`\mathrm{\Delta }_{ab}=\epsilon _a\epsilon _b`$, $`\alpha ^\mu =(1,𝜶)`$ are the Dirac matrices, and $`D_{\mu \nu }(\omega )`$ is the photon propagator.

We separate the contributions of the diagrams under consideration into two parts: the reducible, with the energy of the intermediate state coinciding with the energy of the initial (final) state, and the irreducible, for the remainder, respectively. Omitting the derivation similar to that in Ref. yerokhin01pra2p , we present here only the final expressions for the energy shift. The reducible (”red”) and irreducible (”ir”) three-electron contributions read

$`\mathrm{\Delta }E_{\mathrm{ir}}^{3\mathrm{e}\mathrm{l}}`$ $`=`$ $`{\displaystyle \underset{PQ}{}}(1)^{P+Q}`$ (7)
$`\times `$ $`{\displaystyle \underset{n}{\overset{}{}}}{\displaystyle \frac{I_{P2P3nQ3}(\mathrm{\Delta }_{P3Q3})I_{P1nQ1Q2}(\mathrm{\Delta }_{Q1P1})}{\epsilon _{Q1}+\epsilon _{Q2}\epsilon _{P1}\epsilon _n}},`$
$`\mathrm{\Delta }E_{\mathrm{red}}^{3\mathrm{e}\mathrm{l}}`$ $`=`$ $`{\displaystyle \underset{\mu _a}{}}[I_{vaav}^{}(\mathrm{\Delta }_{va})(I_{ab;ab}I_{bv;bv})`$ (8)
$`+`$ $`{\displaystyle \frac{1}{2}}I_{av\stackrel{~}{v}b}^{}(\mathrm{\Delta }_{va})I_{b\stackrel{~}{v};av}+{\displaystyle \frac{1}{2}}I_{b\stackrel{~}{v}va}^{}(\mathrm{\Delta }_{va})I_{va;\stackrel{~}{v}b}],`$
where $`P`$ and $`Q`$ are the permutation operators, and the prime in the sum in Eq. (7) indicates that terms with the vanishing denominator should be omitted in the summation. In Eq. (8) $`a`$ and $`b`$ denote $`1s`$ electrons with opposite angular-momentum projections $`\mu _a=\mu _b`$, $`v`$ stands for the valence state with the angular-momentum projection $`\mu _v`$, and $`\stackrel{~}{v}`$ is the valence state with $`\mu _{\stackrel{~}{v}}=2\mu _a+\mu _v`$ (the corresponding contribution is assumed to be zero when $`\mu _{\stackrel{~}{v}}`$ is out of the range $`j_v,\mathrm{},j_v`$).

The irreducible two-electron contribution is

$`\mathrm{\Delta }E_{\mathrm{"}\mathrm{ir}\mathrm{"}}^{2\mathrm{e}\mathrm{l}}`$ $`=`$ $`\mathrm{\Delta }E_{\mathrm{dir}}^{\mathrm{lad}}+\mathrm{\Delta }E_{\mathrm{exch}}^{\mathrm{lad}}+\mathrm{\Delta }E_{\mathrm{dir}}^{\mathrm{cr}}+\mathrm{\Delta }E_{\mathrm{exch}}^{\mathrm{cr}},`$ (9)
$`\mathrm{\Delta }E_{\mathrm{dir}}^{\mathrm{lad}}`$ $`=`$ $`{\displaystyle \underset{n_1n_2}{}}{}_{}{}^{}{\displaystyle \frac{i}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\omega `$ (10)
$`\times `$ $`{\displaystyle \frac{F_{\mathrm{dir}}^{\mathrm{lad}}(\omega ,n_1n_2)}{(\epsilon _c\omega \epsilon _{n_1}u)(\epsilon _v+\omega \epsilon _{n_2}u)}},`$
$`\mathrm{\Delta }E_{\mathrm{exch}}^{\mathrm{lad}}`$ $`=`$ $`{\displaystyle \underset{n_1n_2}{}}{}_{}{}^{}{\displaystyle \frac{i}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\omega `$ (11)
$`\times `$ $`{\displaystyle \frac{F_{\mathrm{exch}}^{\mathrm{lad}}(\omega ,n_1n_2)}{(\epsilon _v\omega \epsilon _{n_1}u)(\epsilon _c+\omega \epsilon _{n_2}u)}},`$
$`\mathrm{\Delta }E_{\mathrm{dir}}^{\mathrm{cr}}`$ $`=`$ $`{\displaystyle \underset{n_1n_2}{}}{}_{}{}^{}{\displaystyle \frac{i}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\omega `$ (12)
$`\times `$ $`{\displaystyle \frac{F_{\mathrm{dir}}^{\mathrm{cr}}(\omega ,n_1n_2)}{(\epsilon _c\omega \epsilon _{n_1}u)(\epsilon _v\omega \epsilon _{n_2}u)}},`$
$`\mathrm{\Delta }E_{\mathrm{exch}}^{\mathrm{cr}}`$ $`=`$ $`{\displaystyle \underset{n_1n_2}{}}{}_{}{}^{}{\displaystyle \frac{i}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\omega `$ (13)
$`\times `$ $`{\displaystyle \frac{F_{\mathrm{exch}}^{\mathrm{cr}}(\omega ,n_1n_2)}{(\epsilon _v\omega \epsilon _{n_1}u)(\epsilon _v\omega \epsilon _{n_2}u)}}.`$
Here we introduced the labels ”lad” and ”cr” for the ladder and the crossed diagram, and ”dir” and ”exch” for the direct and the exchange parts. The other notations are:

$`F_{\mathrm{dir}}^{\mathrm{lad}}(\omega ,n_1n_2)`$ $`=`$ $`{\displaystyle \underset{\mu _c\mu _{n_1}\mu _{n_2}}{}}I_{cvn_1n_2}(\omega )I_{n_1n_2cv}(\omega ),`$ (14)
$`F_{\mathrm{exch}}^{\mathrm{lad}}(\omega ,n_1n_2)`$ $`=`$ $`{\displaystyle \underset{\mu _c\mu _{n_1}\mu _{n_2}}{}}I_{vcn_1n_2}(\omega )I_{n_1n_2cv}(\omega \mathrm{\Delta }_{vc}),`$ (15)
$`F_{\mathrm{dir}}^{\mathrm{cr}}(\omega ,n_1n_2)`$ $`=`$ $`{\displaystyle \underset{\mu _c\mu _{n_1}\mu _{n_2}}{}}I_{cn_2n_1v}(\omega )I_{n_1vcn_2}(\omega ),`$ (16)
$`F_{\mathrm{exch}}^{\mathrm{cr}}(\omega ,n_1n_2)`$ $`=`$ $`{\displaystyle \underset{\mu _c\mu _{n_1}\mu _{n_2}}{}}I_{vn_2n_1v}(\omega )I_{n_1ccn_2}(\omega \mathrm{\Delta }_{vc}),`$ (17)
$`u=(1i0)`$, and the prime on the sum indicates that some terms are excluded from the summation. First of all, we omit the reducible contribution, i.e. the terms for which the intermediate two-electron energy $`\epsilon _{n_1}+\epsilon _{n_2}`$ equals the energy of the initial two-electron state $`\epsilon _v+\epsilon _c`$. Those are: $`(\epsilon _{n_1}\epsilon _{n_2})=(\epsilon _c\epsilon _v)`$ and $`(\epsilon _v\epsilon _c)`$. In addition, we exclude also the infrared-divergent terms (see ShabaevFokeeva94 ; yerokhin01pra2p  for details), namely those with $`(\epsilon _{n_1}\epsilon _{n_2})=(\epsilon _c\epsilon _v)`$ in the direct crossed part and with $`(\epsilon _{n_1}\epsilon _{n_2})=(\epsilon _c\epsilon _c)`$ and $`(\epsilon _v\epsilon _v)`$ in the exchange crossed part. These terms should be considered together with the reducible contribution. Their sum can be shown to be infrared finite. We employ the notations $`\mathrm{\Delta }E_{\mathrm{"}\mathrm{ir}\mathrm{"}}^{2\mathrm{e}\mathrm{l}}`$ and $`\mathrm{\Delta }E_{\mathrm{"}\mathrm{red}\mathrm{"}}^{2\mathrm{e}\mathrm{l}}`$ in order to emphasize that the corresponding terms are not ”pure” irreducible and reducible contributions.

We mention that the case under consideration differs from the cases of the $`2s`$ or $`2p_{1/2}`$ valence electrons considered previously in yerokhin01pra2p  by the fact that for the $`2p_{3/2}`$ Dirac state there is no adjoining state separated only by the finite-nuclear-size effect. Consequently, there is no need to exclude any further terms from the crossed contribution, as we had to proceed in Ref. yerokhin01pra2p  in the case of the $`2s`$ and the $`2p_{1/2}`$ valence electron.

Finally, we note the ”reducible” contribution

$`\mathrm{\Delta }E_{\mathrm{"}\mathrm{red}\mathrm{"}}^{2\mathrm{e}\mathrm{l}}`$ $`=`$ $`{\displaystyle \frac{i}{4\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}d\omega {\displaystyle \frac{1}{(\omega +i0)^2}}[2F_{\mathrm{exch}}^{\mathrm{cr}}(\omega +\mathrm{\Delta }_{vc},cc)`$ (18)
$`+`$ $`2F_{\mathrm{exch}}^{\mathrm{cr}}(\omega ,vv)`$
$``$ $`F_{\mathrm{exch}}^{\mathrm{lad}}(\omega +\mathrm{\Delta }_{vc},cv)F_{\mathrm{exch}}^{\mathrm{lad}}(\omega +\mathrm{\Delta }_{vc},cv)`$
$``$ $`F_{\mathrm{dir}}^{\mathrm{lad}}(\omega \mathrm{\Delta }_{vc},vc)F_{\mathrm{dir}}^{\mathrm{lad}}(\omega \mathrm{\Delta }_{vc},vc)`$
$``$ $`F_{\mathrm{exch}}^{\mathrm{lad}}(\omega ,vc)F_{\mathrm{exch}}^{\mathrm{lad}}(\omega ,vc)].`$
## III Numerical evaluation

The three-electron contribution to the energy of $`(1s)^22s`$, $`(1s)^22p_{1/2}`$, and $`(1s)^22p_{3/2}`$ levels of Li-like ions has been calculated in our recent investigation sysak02 . This evaluation appears as relatively simple since the corresponding expressions (7), (8) contain at most one summation over the Dirac spectrum and no integrations over the virtual-photon energy. Thus we focus here on the calculation of the two-electron contribution.

The summation over magnetic substates in Eqs. (10)-(13), (18) was performed by means of standard techniques. The resulting expressions can be found in yerokhin01pra2p . As an independent check we employed also the direct numerical summation of Clebsch-Gordan coefficients.

To calculate infinite summations over the spectrum of the Dirac equation in Eqs. (10)-(13), we employed the method of the B-spline basis set for the Dirac equation johnson88 . Typical basis sets contained 50 positive and 50 negative-energy eigenstates for each value of the angular-momentum quantum number $`\kappa `$. The finite size of the nucleus has been taken into account employing the homogeneously-charged sphere model for the nuclear-charge distribution. The values of the rms radii used in this work are the same as in yerokhin01pra2p . Infinite summations over $`\kappa `$ were truncated typically at $`|\kappa |=10`$. Partial sums of the expansion over $`|\kappa |`$ were fitted to the form

$`S_{|\kappa |}=a_0+{\displaystyle \underset{n=2}{\overset{N}{}}}{\displaystyle \frac{a_n}{|\kappa |^n}}`$ (19)
using the least squares method. The coefficient $`a_0`$ yields the extrapolated value for the sum of the expansion. We found that different fits with $`N=`$ 4-6 yield the same result with an accuracy of at least 5 digits.

The integration over the energy of the virtual photon $`\omega `$ in Eqs. (10)-(13) represents the most difficult part of the calculation. To avoid strong oscillations for large values of $`\omega `$, we performed the Wick rotation of the integration contour. Deforming the contour, one should take care about the poles and the branch cuts of the integrand. The analytic structure of the integrand for Eqs. (11)-(13) is shown in Figs. 2-5. These graphs are very similar to those for the $`2s`$\- and $`2p_{1/2}`$-valence electrons in Ref. yerokhin01pra2p . The only difference is that now three Dirac energy levels occur which are more deeply bound than the valence state: $`1s`$, $`2s`$, and $`2p_{1/2}`$. The terms in Eqs. (11) and (13) containing these states and the valence state as intermediate were treated in a different way than the remainder, as is discussed below.

For the evaluation of the direct parts of the ladder and crossed contributions, we perform the Wick rotation of the integration contours separating the corresponding pole contributions, as shown in Figs. 2 and 3. In the direct part of the reducible contribution, we also perform a Wick rotation and then integrate by parts. This yields the following expression which can be evaluated directly,

$`\mathrm{\Delta }E_{\mathrm{"}\mathrm{red}\mathrm{"},\mathrm{dir}}^{2\mathrm{e}\mathrm{l}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[F_{\mathrm{dir}}^{\mathrm{lad}}(\mathrm{\Delta }_{vc},vc)\right]^{}`$ (20)
$``$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle _0^{\mathrm{}}}𝑑\omega {\displaystyle \frac{\omega }{\mathrm{\Delta }_{vc}^2+\omega ^2}}{\displaystyle \frac{d}{d\omega }}F_{\mathrm{dir}}^{\mathrm{lad}}(i\omega ,vc),`$
where $`F^{}(\mathrm{\Delta })=(dF/d\omega )_{\omega =\mathrm{\Delta }}`$.

Let us now turn to the exchange contribution. As one can see from Figs. 4 or 5, in this case the integration contour is squeezed between two branch cuts of the photon propagators on the interval $`[0,\mathrm{\Delta }_{vc}]`$. Therefore, the standard Wick rotation of the contour is not possible. It is convenient to divide the contributions of Eqs. (11) and (13) into two parts. The first one accounts for the poles of the integrand on the interval $`[0,\mathrm{\Delta }_{vc}]`$ and is referred to as the irregular part. The remainder is denoted as the regular part. This contribution does not possess any poles close to the squeezed part of the contour, which simplifies its numerical evaluation. However, it turns out as is the most time-consuming part of the calculation. One of the integration contours $`C_{\mathrm{reg}}`$ used for the evaluation of the regular part is depicted in Fig. 5.

The evaluation of the irregular part is less time consuming, but its structure is more difficult. In this case we need to take care of single and double poles of the integrand that are located close to the integration contour. The potential occurrences of one or two single poles and one double pole within the interval $`[0,\mathrm{\Delta }_{vc}]`$ were treated by means of the following identities:

$`{\displaystyle _{\omega _1}^{\omega _2}}𝑑\omega {\displaystyle \frac{f(\omega )}{x_0\omega \pm i0}}`$ $`=`$ $`P{\displaystyle _{\omega _1}^{\omega _2}}𝑑\omega {\displaystyle \frac{f(\omega )}{x_0\omega }}i\pi f(x_0),`$ (21)
$`{\displaystyle _{\omega _1}^{\omega _2}}𝑑\omega {\displaystyle \frac{f(\omega )}{(x_0\omega \pm i0)^2}}`$ $`=`$ $`\pm i\pi f^{}(x_0)+{\displaystyle \frac{f(\omega _2)}{x_0\omega _2}}`$ (22)
$`{\displaystyle \frac{f(\omega _1)}{x_0\omega _1}}P{\displaystyle _{\omega _1}^{\omega _2}}𝑑\omega {\displaystyle \frac{f^{}(\omega )}{x_0\omega }},`$
$`{\displaystyle _{\omega _1}^{\omega _2}}𝑑\omega {\displaystyle \frac{f(\omega )}{(x_0\omega \pm i0)(x_1\omega \pm i0)}}`$ $`=`$ $`{\displaystyle \frac{1}{x_1x_0}}[P{\displaystyle _{\omega _1}^{\omega _2}}d\omega {\displaystyle \frac{f(\omega )}{x_0\omega }}`$ (23)
$`P{\displaystyle _{\omega _1}^{\omega _2}}d\omega {\displaystyle \frac{f(\omega )}{x_1\omega }}i\pi f(x_0)\pm i\pi f(x_1)],`$
where $`P`$ indicates the principal value of the integral. In Eq. (23) the choice of the sign before $`i\pi f(x_0)`$ and $`i\pi f(x_1)`$ is determined by the sign of the infinitesimal addition $`\pm i0`$ in the first and the second denominator, respectively. For the numerical evaluation of the irregular contribution we employed the integration contour $`C_{\mathrm{irr}}`$ shown in Fig. 4. It consists of 3 parts: $`[i\mathrm{}ϵ,ϵ]`$, $`[ϵ,\mathrm{\Delta }_{vc}+ϵ]`$, and $`[\mathrm{\Delta }_{vc}+ϵ,\mathrm{\Delta }_{vc}+ϵ+i\mathrm{}]`$. A small positive constant $`ϵ`$ was introduced in order to facilitate the numerical evaluation of the principal-value integrals.

After integration by parts, the exchange contribution of the reducible part can be written as

$`\mathrm{\Delta }E_{\mathrm{"}\mathrm{red}\mathrm{"},\mathrm{exch}}^{2\mathrm{e}\mathrm{l}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[F_{\mathrm{exch}}^{\mathrm{cr}}(\mathrm{\Delta }_vc,cc)+F_{\mathrm{exch}}^{\mathrm{cr}}(0,vv)\right]^{}`$ (24)
$`+{\displaystyle \frac{1}{2\pi i}}P{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{d\omega }{\omega }}{\displaystyle \frac{d}{d\omega }}[F_{\mathrm{exch}}^{\mathrm{cr}}(\mathrm{\Delta }_vc+\omega ,cc)`$
$`+F_{\mathrm{exch}}^{\mathrm{cr}}(\omega ,vv)2F_{\mathrm{exch}}^{\mathrm{lad}}(\omega ,vc)].`$
It is again worth mentioning that the integral in Eq. (24) exists only if the sum of all 3 terms in the brackets is considered. For the each single term, the integral is infrared divergent.

## IV Numerical results and discussion

The results of our calculations are presented in Table 1, where the direct, the exchange, and the three-electron contribution to the two-photon exchange correction of the valence $`2p_{3/2}`$ electron with the $`(1s)^2`$ shell are listed separately. The evaluation was performed within the Feynman gauge. We estimate the numerical uncertainty of our results to be less than $`5\times 10^5`$ a.u. For bismuth, our results can be compared with the calculation by Sapirstein and Cheng sapirstein01 . They report $`6.529`$ and $`6.670`$ eV for the two-electron and the three-electron contribution, respectively. This agrees well with our corresponding results of $`6.5330`$ and $`6.6698`$ eV, respectively.

It is interesting to compare the results of the rigorous QED treatment with approximations evaluations based on relativistic many-body perturbation theory (MBPT). The difference between the QED and MBPT results can be conventionally regarded as a ”nontrivial” QED contribution. In order to deduce the two-photon exchange correction within the framework of MBPT, we should introduce the following changes in our basic formulas: all summations over intermediate states should be restricted to positive-energy states only, the calculation should be performed within Coulomb gauge, and the virtual-photon energy in the photon propagator should be set equal to zero. Within this approximation, all reducible parts vanish, and the integration over the energy of the virtual photon can be carried out employing Cauchy’s theorem. This yields zero for the crossed contribution, and finally we are left with the following expression for the total two-photon exchange correction within the MBPT approximation:

$`\mathrm{\Delta }E_{\mathrm{MBPT}}^{3\mathrm{e}\mathrm{l}}`$ $`=`$ $`{\displaystyle \underset{PQ}{}}(1)^{P+Q}{\displaystyle \underset{\epsilon _n>0}{\overset{}{}}}{\displaystyle \frac{I_{P2P3nQ3}(0)I_{P1nQ1Q2}(0)}{\epsilon _{Q1}+\epsilon _{Q2}\epsilon _{P1}\epsilon _n}},`$ (25)
$`\mathrm{\Delta }E_{\mathrm{MBPT}}^{2\mathrm{e}\mathrm{l}}`$ $`=`$ $`{\displaystyle \underset{\mu _c}{}}{\displaystyle \underset{\epsilon _{n_1}\epsilon _{n_2}>0}{\overset{}{}}}{\displaystyle \frac{\left[I_{cvn_1n_2}(0)I_{vcn_1n_2}(0)\right]I_{n_1n_2cv}(0)}{\epsilon _c+\epsilon _v\epsilon _{n_1}\epsilon _{n_2}}},`$ (26)
where the photon propagators should be taken in the Coulomb gauge and the prime indicates that terms with vanishing denominator should be omitted. We mention that Eqs. (25) and (26) include the contribution due to the exchange by two Breit photons (the $`B\times B`$ term). Strictly speaking, this term is of higher order than the level of validity of the Breit approximation, and, therefore, it appears to be inconsistent to include it within the MBPT scheme.

In Table 2 and in Fig. 6 we compare the results of the rigorous QED treatment of the two-photon exchange correction to the $`2p_{3/2}`$-$`2s`$ splitting with the complete MBPT result \[Eqs. (25) and (26)\], and the MBPT result dropping the $`B\times B`$ term. A similar analysis for the $`2p_{1/2}`$-$`2s`$ splitting has been presented in our previous investigation yerokhin01pra2p . Our results show that in the case under consideration the nontrivial QED contribution is essentially larger than that for the $`2p_{1/2}`$-$`2s`$ transition. E.g., for uranium it yields $`0.011`$ a.u., while the corresponding contribution to the $`2p_{1/2}`$-$`2s`$ splitting is two times smaller, of about $`0.006`$ a.u. Moreover, we see that in our case the total correction changes its sign in the region between $`Z=92`$ and 100. As a result, the MBPT result becomes incorrect by more than 50% at very high values of $`Z`$. A further conclusion that can be drawn from our comparison is that for the $`2p_{3/2}`$-$`2s`$ splitting the $`B\times B`$ term is of the same sign and magnitude as the nontrivial QED contribution. Thus, its inclusion improves the agreement between the MBPT and the QED result. This situation is contrary to the one for the $`2p_{1/2}`$-$`2s`$ splitting, where the $`B\times B`$ term turns out to be of the same order of magnitude, but of different sign than the nontrivial QED contribution.

To summarize this investigation we presented a rigorous QED evaluation of the two-photon exchange correction for the $`(1s)^22p_{3/2}`$ state of Li-like ions. Combining these results with the data for the $`(1s)^22s`$ state from our previous study yerokhin01pra2p , we obtained the two-photon exchange correction for the $`2p_{3/2}`$-$`2s`$ splitting. This is an important step towards the final goal consisting in the evaluation of all two-electron second-order QED corrections to the $`2p_{3/2}`$-$`2s`$ transition energy for the Li isoelectronic sequence.

## V Acknowledgments

This work was supported by the Russian Foundation for Basic Research (Grant No. 01-02-17248), by the Russian Ministry of Education (Grant No. E02-3.1-49), and by the program ”Russian Universities” (Grant No. UR.01.01.072). The work of A.N.A. and V.M.S. was supported by the joint grant of the Russian Ministry of Education and the Administration of Saint Petersburg (Grant No. PD02-1.2-79). V.A.Y. acknowledges the support of the foundation ”Dynasty” and the International Center for Fundamental Physics. The work of V.M.S. was supported by the Alexander von Humboldt Stiftung. We also acknowledge support from BMBF, DFG, DAAD, and GSI.


text:
# 1 Introduction

## 1 Introduction

The study of nuclear structure functions has become a fashionable subject. Apart from its intrinsic interest, such analysis has a great impact on the interpretation of results from heavy ion experiment. At small values of the Bjorken variable $`x`$ ($`<0.01`$, shadowing region), the structure function $`F_2`$ per nucleon turns out to be smaller in nuclei than in a free nucleon . Several explanations to this shadowing have been proposed.

On the one hand, some models use the fact that, in the rest frame of the nucleus the incoming photon splits into a $`q\overline{q}`$ pair long before reaching the nucleus, and this $`q\overline{q}`$ pair interacts with it with typical hadronic cross sections, which results in absorption . Thus nuclear shadowing is a consequence of multiple scattering which in turn is related to diffraction . This relationship will be developed in this paper. Equivalently, in a frame in which the nucleus is moving fast, gluon recombination due to the overlap of the gluon clouds from different nucleons reduces the gluon density in the nucleus . These studies have received great theoretical impulse with the development of semiclassical ideas in QCD and the appearance of non-linear equations for evolution in $`x`$ in this framework (see  and references therein; also  for a simple geometrical approach in this framework).

On the other hand, other approaches  do not address the origin of shadowing but only its evolution with $`\mathrm{ln}Q^2`$: parton densities inside the nucleus are parametrized at some scale $`Q_0^2`$ and then evolved using the DGLAP  evolution equations.

The results from different models usually depend on phenomenological assumptions and their predictions (notably for small values of $`x`$ which are of uttermost importance to compute particle production at RHIC and LHC) turn out to be very different. For example, concerning the $`Q^2`$-dependence of shadowing, it can be either constant , or die out logarithmically  or behave as a higher-twist .

In this paper we will use the relation of diffraction to nuclear shadowing which arises from Gribov Theory , Reggeon Calculus  and the AGK rules . In this way we obtain a parameter-free description of nuclear structure functions in the shadowing region valid for $`x<0.01`$ and $`Q^2<10`$ GeV2, using a model for $`F_2`$ and $`F_{2𝒟}`$ . The same strategy has been used in , but our extrapolation to smaller $`x`$ or higher $`W^2`$ is more reliable than that of  due to the model employed for the nucleon; besides, our description is valid for small $`Q^2`$ while that of  applies to $`Q^24`$ GeV2. In Section 2 the model will be described. In Section 3 numerical results will be presented together with comparisons with experimental data and with other models. In Section 4 the model will be applied to calculate the multiplicity reduction factors  relevant to compute particle production in heavy ion collisions at RHIC and LHC. Finally, the last Section will contain our conclusions.

## 2 Description of the model

We assume that nuclei are made of nucleons in the spirit of the Glauber model. In order to relate diffraction on nucleons with nuclear shadowing, we will follow the procedure explained in . The $`\gamma ^{}`$-nucleus cross section can be expanded in a multiple scattering series containing the contribution from 1, 2,$`\mathrm{}`$ scatterings:

$$\sigma _A=\sigma _A^{(1)}+\sigma _A^{(2)}+\mathrm{}.$$
(1)
$`\sigma _A^{(1)}`$ is simply equal to $`A\sigma _{\mathrm{nucleon}}`$. Let us consider now the first correction to the non-additivity of cross sections which comes from the second-order rescattering $`\sigma _A^{(2)}`$. In Fig. 1 diffractive DIS is shown in both the infinite momentum frame and in the rest frame of the nucleon. In Fig. 2 it becomes clear that the square of such contribution is equivalent to a double exchange with a cut between the exchanged amplitudes, a so-called diffractive cut. To compute the first contribution to nuclear shadowing $`\sigma _A^{(2)}`$, which comes from these two exchanges, we need its total contribution to the $`\gamma ^{}`$-nucleon cross section, which arises from cutting the two-exchange amplitude in all possible ways (between the amplitudes and the amplitudes themselves in all possible manners). For purely imaginary amplitudes, it can be shown  that this total contribution is identical to minus the contribution from the diffractive cut. Thus diffractive DIS is directly related to the first contribution to nuclear shadowing. The final expression reads

$$\sigma _A^{(2)}=4\pi A(A1)d^2bT_A^2(b)_{M_{min}^2}^{M_{max}^2}𝑑M^2\frac{d\sigma _{\gamma ^{}\mathrm{p}}^𝒟}{dM^2dt}|_{t=0}F_A^2(t_{min}),$$
(2)
with $`T_A(b)=_{\mathrm{}}^+\mathrm{}𝑑z\rho _A(\stackrel{}{b},z)`$ the nuclear profile function normalized to 1, $`d^2bT_A(\stackrel{}{b})=1`$, and $`M^2`$ the mass of the diffractively produced system. The usual variables for diffractive DIS: $`Q^2`$, $`x`$, $`M^2`$ and $`t`$, or $`x_P=x/\beta `$, $`\beta =\frac{Q^2}{Q^2+M^2}`$, are shown in Fig. 1.

Coherence effects, i.e. the coherence length of the $`q\overline{q}`$ fluctuation of the incoming virtual photon, is taken into account through

$$F_A(t_{min})=d^2bJ_0(b\sqrt{t_{min}})T_A(b),$$
(3)
with $`t_{min}=m_N^2x_P^2`$ and $`m_N`$ the nucleon mass. This function is equal to 1 at $`x0`$ and decreases with increasing $`x`$ due to the loss of coherence for $`x>x_{crit}(m_NR_A)^1`$.

Let us briefly examine (2). Here the real part of the Pomeron amplitude, which is small for the value of the intercept which will be used , $`\mathrm{\Delta }=\alpha _P(t=0)1=0.2`$, has not been taken into account. Also it has been deduced under the approximation $`R_A^2R_N^2`$, so the $`t`$-dependence of the $`\gamma ^{}`$-nucleon cross section has been neglected.

For $`A>20`$ a nuclear density in the form of a 3-parameter Fermi distribution with parameters taken from  will be employed to compute both $`T_A(b)`$ and (3). For $`2<A20`$ a Gaussian profile function is used :

$$T_A(b)=\frac{3}{2\pi R_A^2}\mathrm{exp}\left(\frac{3b^2}{2R_A^2}\right),R_A=0.82A^{1/3}+0.58\mathrm{fm},$$
(4)
but, in order to take into account the $`t`$-dependence for these light nuclei, we make in the computation of form factors (3) the substitution

$$R_A^2R_A^2+R_N^2,R_N=0.8\mathrm{fm}.$$
(5)
Finally, for deuteron the double rescattering contribution has the form

$$\sigma _A^{(2)}=2_{\mathrm{}}^{t_{min}}𝑑t_{M_{min}^2}^{M_{max}^2}𝑑M^2\frac{d\sigma _{\gamma ^{}\mathrm{nucleon}}^𝒟}{dM^2dt}|_{t=0}F_D(t),$$
(6)
where $`F_D(t)=e^{at}`$, $`a=40`$ GeV-2.

The lower integration limit in (2) and (6) is $`M_{min}^2=4m_\pi ^2=0.08`$ GeV2, while the upper one is taken from the condition:

$$x_P=x\left(\frac{M^2+Q^2}{Q^2}\right)x_{Pmax}M_{max}^2=Q^2\left(\frac{x_{Pmax}}{x}1\right),$$
(7)
with $`x_{Pmax}=0.1`$; this value was used in  motivated by the fact that the model is only valid for $`M^2W^2`$ or $`x_P1`$, i.e. a large rapidity gap is required. In our case, variations of $`x_{Pmax}`$ by a factor 2 do not affect the description of nuclear shadowing at $`x<0.01`$ but the choice $`x_{Pmax}=0.1`$ is convenient as it guarantees the disappearance of nuclear shadowing at $`x0.1`$ (see below) as in the experimental data.

The relation between $`\frac{d\sigma _{\gamma ^{}\mathrm{p}}^𝒟}{dM^2dt}|_{t=0}`$ and $`x_PF_{2𝒟}^{(3)}(Q^2,x_P,\beta )`$ is provided by the model :

$`x_PF_{2𝒟}^{(3)}(Q^2,x_P,\beta )`$ $`=x_P{\displaystyle \frac{Q^2}{4\pi ^2\alpha _{em}}}{\displaystyle _{\mathrm{}}^0}𝑑t{\displaystyle \frac{d\sigma _{\gamma ^{}\mathrm{p}}^𝒟(Q^2,x_P,\beta ,t)}{dx_Pdt}}`$ (8)
$``$ $`{\displaystyle \frac{d\sigma _{\gamma ^{}\mathrm{p}}^𝒟(Q^2,x_P,\beta )}{dM^2dt}}|_{t=0}={\displaystyle \frac{4\pi ^2\alpha _{em}B}{Q^2(Q^2+M^2)}}x_PF_{2𝒟}^{(3)}(Q^2,x_P,\beta ),`$
where the usual factorization has been assumed:

$$\frac{d\sigma _{\gamma ^{}\mathrm{p}}^𝒟(x,Q^2,M^2,t)}{dM^2dt}=\frac{d\sigma _{\gamma ^{}\mathrm{p}}^𝒟(x,Q^2,M^2)}{dM^2dt}|_{t=0}e^{Bt},$$
(9)
with $`B=6`$ GeV-2 (as in , see the discussion there; this value is slightly smaller than the experimental values $`7.2\pm 1.1(\mathrm{stat}.)_{0.9}^{+0.7}(\mathrm{syst}.)`$ GeV-2  at $`Q^2=8`$ GeV2 and $`6.8\pm 0.9(\mathrm{stat}.)_{1.1}^{+1.2}(\mathrm{syst}.)`$ GeV-2  for photoproduction). Note that $`\frac{d\sigma _{\gamma ^{}\mathrm{p}}^𝒟(x,Q^2,M^2)}{dM^2dt}|_{t=0}`$ can be obtained directly from $`\sigma _{tot}`$. However, the model for diffraction we are using  has mainly been tested after integration in $`t`$ (most available data are integrated in $`t`$). For this reason, we use the integrated expression together with the experimental value of $`B`$. While this is legitimate at present values of $`x`$, it can lead to an underestimation of shadowing at very small $`x`$, due to the increase of $`B`$ with energy111Nevertheless, the effect is not too large: we have checked that an increase of $`B`$ from 6 to 7.2 GeV-2 produces an increase of shadowing for Pb of at most 10 % at $`x=10^7`$. As estimates indicate an increase $`<50`$ % in $`B`$ for the smallest $`x`$ we have studied, $`x=10^7`$, the increase of shadowing due to this effect would be at most $`25`$ % for these values of $`x`$..

The model in  is based on the dipole picture of the photon and contains two components. The small-distance ($`S`$) component corresponds to transverse distances $`r`$ between the $`q`$ and the $`\overline{q}`$ of the dipole such that $`r<r_0`$, and the large-distance ($`L`$) component to $`r>r_0`$, with $`r_0=0.2`$ fm. In each component a quasi-eikonal iteration is introduced in order to enforce unitarity. Reggeon and Pomeron exchanges are allowed. For diffraction, a third component is used, namely a contribution from the triple interaction of Reggeons and Pomerons. This model has been designed to describe the small $`x<10^2`$, small or moderate $`Q^2<10`$ GeV2 region, and it contains the basic ingredients which allow its safe extrapolation222In order to use the model for larger $`x`$, $`0.01<x<0.1`$, we have made some modifications in : there, in Eq. (26) $`\beta _{min}`$ in the normalization denominators has been set to 0, and in Eq. (25) the Reggeon-Reggeon contribution has been ignored. These two changes slightly modify the description of diffraction but we have checked that the agreement with experimental data is as good as in the original version of the model. to very small $`x`$ or high $`W^2`$.

Eq. (2) corresponds to the case with only two scatterings. Its extension to include higher order rescatterings is model-dependent. We will use two models: a Schwimmer unitarization  which is obtained from a summation of fan diagrams with triple Pomeron interactions,

$$\sigma _{\gamma ^{}A}^{Sch}=\sigma _{\gamma ^{}\mathrm{nucleon}}d^2b\frac{AT_A(b)}{1+(A1)f(x,Q^2)T_A(b)},$$
(10)
and an eikonal unitarization,

$$\sigma _{\gamma ^{}A}^{eik}=\sigma _{\gamma ^{}\mathrm{nucleon}}d^2b\frac{A}{2(A1)f(x,Q^2)}\left\{1\mathrm{exp}\left[2(A1)T_A(b)f(x,Q^2)\right]\right\},$$
(11)
where we use the relation $`\sigma _{\gamma ^{}\mathrm{nucleon}}=\frac{4\pi ^2\alpha _{em}}{Q^2}F_2(x,Q^2)`$ valid at small $`x`$. Here, $`F_2(x,Q^2)`$ is the nucleon structure function, taken from . Both expressions (10) and (11), expanded to the first non-trivial order, reproduce the second rescattering result (2). Eikonal unitarization will produce larger shadowing than Schwimmer, as can be expected by comparing the second non-trivial order in the expansion of both expressions. Finally,

$$f(x,Q^2)=\frac{4\pi }{\sigma _{\gamma ^{}\mathrm{nucleon}}}_{M_{min}^2}^{M_{max}^2}𝑑M^2\frac{d\sigma _{\gamma ^{}\mathrm{p}}^𝒟}{dM^2dt}|_{t=0}F_A^2(t_{min})$$
(12)
as required by consistency with (2).

The shadowing in nuclei is usually studied through the ratios of cross sections per nucleon for different nuclei, defined as

$$R(A/B)=\frac{B}{A}\frac{\sigma _{\gamma ^{}A}}{\sigma _{\gamma ^{}B}}.$$
(13)
In the simplest case of the ratio over nucleon (equivalent to proton at small $`x`$ where the valence contribution can be neglected), we get:

$$R^{Sch}(A/\mathrm{nucleon})=d^2b\frac{T_A(b)}{1+(A1)f(x,Q^2)T_A(b)},$$
(14)
$$R^{eik}(A/\mathrm{nucleon})=d^2b\frac{1}{2(A1)f(x,Q^2)}\left\{1\mathrm{exp}\left[2(A1)T_A(b)f(x,Q^2)\right]\right\}.$$
(15)
To calculate shadowing in photoproduction, $`x`$ is no longer a relevant kinematical variable. Instead we use the $`\gamma ^{}`$-nucleon center of mass energy $`W^2`$.

In our framework shadowing can also be studied as a function of the impact parameter $`b`$:

$$R(A/\mathrm{nucleon})^{Sch}(b)=\frac{1}{1+(A1)f(x,Q^2)T_A(b)},$$
(16)
$$R(A/\mathrm{nucleon})^{eik}(b)=\frac{1}{2(A1)T_A(b)f(x,Q^2)}\left\{1\mathrm{exp}\left[2(A1)T_A(b)f(x,Q^2)\right]\right\}.$$
(17)
Finally, the region of applicability of our model is the same as that of the model for diffraction on nucleon , i.e. small $`x<0.01`$ and small or moderate $`Q^2<10`$ GeV2, including photoproduction.

## 3 Numerical results

In our model and in  we work in the small $`x`$ region and thus no distinction is made between protons and neutrons. Although usually joined with straight lines, our results are computed at the same $`x`$ and $`Q^2`$ as the experimental data. For the latter, inner error bars show statistical errors, and outer error bars correspond to statistical and systematical errors added in quadrature.

In Figs. 3-6 a comparison with experimental data at small $`x`$ from E665  and NMC  is presented. As expected, eikonal unitarization produces larger shadowing than Schwimmer. The agreement with experimental data is quite reasonable taking into account that no parameters have been fitted to reproduce the data. Two comments are in order: First, for C/D and Ca/D in Fig. 3 which shows the comparison with E665 data, shadowing looks overestimated for $`x0.01`$, while in Fig. 5 which shows the comparison with NMC data, it looks underestimated. This corresponds to the known difference between the results of both experiments for ratios over D, while the compatibility is restored  when ratios are computed over C. Second, from Fig. 6 it becomes clear that the evolution with $`Q^2`$ in the model is too slow at $`x0.01`$, a problem related with the lack of DGLAP evolution in the model  (see  for an application of DGLAP evolution to initial conditions).

In Fig. 7 a comparison of the results of our model with those of others is shown, for $`Q^2=3`$ GeV2 (except the results of  which are at $`Q^2=4`$ GeV2). It can be seen that the results of different models agree within 15 % at $`x0.01`$ where experimental data exist, while they differ up to a factor 0.6 at $`x=10^5`$. At this $`x`$, our results are the lowest ones but roughly agree with those of  and with one set of , while the results from  are the highest ones, and those of  and the second set of  lie in between. Let us briefly comment on these models: In  an initial condition is parameterized at some $`Q_0^2`$ and then evolved using DGLAP; the initial condition is fitted from the comparison of the evolved results with experimental data (see  for a comparison between these two models).  is a model which uses a saturating ansatz for the total $`\gamma ^{}`$-nucleon cross section in the proton, which is introduced in a Glauber expression for its extension to the nuclear case. In  some parameterization of hard diffraction at $`Q_0^2`$, which as in the present work gives nuclear shadowing through Gribov’s Reggeon Calculus, is employed; this nuclear shadowing computed at $`Q_0^2`$ is used as initial condition for DGLAP evolution. In  a Glauber ansatz provides with the initial condition for DGLAP evolution. Finally, in  a non-linear equation for small $`x`$ evolution is numerically solved  and used in the nuclear case. In view of the differences at small $`x`$ among different models, a measurement of $`F_2`$ in nuclei with $`10`$ % precision would be a sensitive test to discriminate among them. Lepton-ion colliders  could provide with such data.

In Fig. 8 our predictions for the ratios D, He, Li, C, Ca, Sn and Pb over nucleon for $`Q^2=0.5`$, 2 and 5 GeV2 are given for $`x>10^8`$. Let us notice that our model is designed for the small $`x`$ region and that no antishadowing or any other effects relevant for $`x>0.1`$ have been introduced. The disappearance of shadowing at $`x0.1`$ is a consequence of both the coherence effects in (3) and the vanishing integration domain in (2), see (7). In Fig. 9 results in photoproduction for the same ratios as in Fig. 8 are given for $`W^2<10^5`$ GeV2, together with predictions for the evolution of the ratios C and Pb over nucleon with impact parameter $`b`$. Values as low as 0.3 are reached for central Pb/nucleon. This evolution with centrality is very important to compute the corresponding evolution of particle production in nuclear collisions, and could also be measured in lepton-ion colliders .

As a last comment in this Section, let us discuss about the twist structure of the model (i.e. its structure in powers of $`1/Q^2`$). In the model of  the unitarity corrections to the $`L`$-component are all of order $`1/Q^2`$. On the contrary, in the $`S`$-component the unitarity corrections are higher-twist (they can be expanded as a sum of terms, each one containing an additional factor $`1/Q^2`$ as compared to the previous one). The fact that diffraction is related to the unitarity corrections allows to study the $`1/Q^2`$ behavior of shadowing in this model. In order to keep only the leading-twist contribution (terms $`1/Q^2`$ in the cross section) we ignore the higher-twist contribution of the $`S`$-component to the diffractive cross section333Concretely, we ignore the $`S`$-component in (17) of , and in (20) of  we set the exponential containing $`\chi _S`$ to 1.. The results are given in Fig. 10. One can see that neglecting these terms introduces only a small difference. The fact that nuclear shadowing corrections are predominantly leading-twist is not unexpected, as the diffractive cross section is also leading-twist for the relevant kinematical region (indeed, in the model of  the $`S`$-component diffraction is almost negligible for small $`Q^2`$ and/or large $`M^2`$). This is also seen in the fact that the ratio of diffractive to inclusive cross sections does not show any strong $`Q^2`$-dependence for large $`M^2`$ . Here a comment is in order: in  a parameter $`s_0`$ is introduced in $`x`$ and $`\beta `$ to control the limit $`Q^20`$, so that all the equations are written for $`\overline{x}=x+s_0/(W^2+Q^2)`$, $`\overline{\beta }=\beta +s_0/(M^2+Q^2)`$. These terms could mimic higher-twist corrections. In Fig. 10 we check that the effect of varying this parameter444And setting $`c=0`$ in Eq. (27) of . from the original $`s_0=0.79`$ GeV2 to $`s_0=0.2`$ GeV2 is also very small. So, we can conclude that the contribution from higher-twist terms to the shadowing of $`F_2`$ is small. In contrast, in  a large higher-twist correction for the shadowing is claimed. The approach in this reference is very similar to ours: the authors also compute shadowing from the diffractive cross section, but using the H1 parametrization . The shadowing obtained in this way for $`Q_0^2=4`$ GeV2 is then employed as initial condition for DGLAP evolution, taking the shadowing for valence quarks from . When evolved downwards to $`Q^2=3`$ GeV2 a disagreement is found with experimental data on the ratio Ca over D. This disagreement is attributed to higher-twist contributions.

## 4 Multiplicity reduction in nucleus-nucleus collisions

In this framework it is also possible to study the reduction of multiplicities in nucleus-nucleus collisions . We will denote the rapidity of the produced system in the center of mass frame by $`y^{}`$. Shadowing as a function of the rapidity of the produced particle can be computed taking into account the general relation with the diffractive variables:

$$y=\mathrm{ln}\left(\frac{1}{x_P}\right)=\mathrm{ln}\left(\frac{s}{M^2}\right).$$
(18)
Then the factor for reduction of multiplicities at fixed impact parameter $`b`$ is

$$R_{AB}(b)=\frac{d^2sR_A(\stackrel{}{s})R_B(\stackrel{}{b}\stackrel{}{s})}{T_{AB}(b)}.$$
(19)
$`R_{A(B)}(b)`$ is given by the r.h.s. of (16) multiplied by $`T_{A(B)}(b)`$ and with $`f(x,Q^2)`$ substituted by $`F(s,y^{})`$ (see below), and

$$T_{AB}(b)=d^2sT_A(\stackrel{}{s})T_B(\stackrel{}{b}\stackrel{}{s}).$$
(20)
(19) takes into account the summation of Schwimmer’s fan-like diagrams for the projectile and target, which are joined by a single Pomeron whose cut gives rise to the produced particle (Fig. 11). It follows from AGK cancellation  that this is the only contribution of this type (more complicated diagrams with lines joining upper and lower parts of the diagram, cancel). This provides the justification for the factorized expression (19), which is true even if more general rescattering diagrams are taken into account.

The reduction factor as a function of the rapidity of the produced particles $`F(s,y^{})`$, can be calculated in several ways. The first one is using (12), but with the integration limits inspired by the parton model for hard processes: for projectile $`A`$ (target $`B`$),

$$x_{A(B)}=\frac{m_T}{\sqrt{s}}e^{\pm y^{}},$$
(21)
with $`y^{}>0`$ for the projectile hemisphere and $`y^{}<0`$ for the target one, and $`m_T=\sqrt{m^2+p_T^2}`$ the transverse mass of the produced particle. Substituting in the general relation for $`M_{max}^2`$, (7), we get

$$M_{max}^{2(A(B))}=Q^2\left(\frac{x_{Pmax}}{x_{A(B)}}1\right)=Q^2\left(\frac{x_{Pmax}\sqrt{s}}{m_T}e^y^{}1\right),$$
(22)
while $`M_{min}^2`$ remains fixed and equal to 0.08 GeV2, and $`Q^2=m_T^2`$.

On the other hand, we can also compute the reduction factor from the formulas

$$F(s,y^{})=4\pi _{y_{min}}^{y_{max}}𝑑y\frac{1}{\sigma _P(s)}\frac{d\sigma ^{PPP}}{dydt}|_{t=0}F_A^2(t_{min}),$$
(23)
where $`\sigma _P(s)`$ is the single Pomeron exchange cross section and $`\frac{d\sigma ^{PPP}}{dydt}`$ the triple Pomeron cross section. Using the standard triple Pomeron formula for the latter, we get

$$\frac{1}{\sigma _P(s)}\frac{d\sigma ^{PPP}}{dydt}|_{t=0}=C\mathrm{\Delta }\mathrm{exp}(\mathrm{\Delta }y),$$
(24)
with $`C=\frac{g_{pp}^P(0)r_{PPP}(0)}{4\mathrm{\Delta }}`$, $`g_{pp}^P(0)`$ the Pomeron-proton coupling and $`r_{PPP}(0)`$ the triple Pomeron coupling, both evaluated at $`t=0`$. In this case, the same integration limits used above correspond to:

$$y_{min}^{(A(B))}=\mathrm{ln}\left(\frac{s}{M_{max}^{2(A(B))}}\right),$$
(25)
with $`M_{max}^{2(A(B))}`$ given by (22), and

$$y_{max}^{(A(B))}=\frac{1}{2}\mathrm{ln}\left(\frac{s}{m_T^2}\right)y^{}.$$
(26)
In the calculations we have used \[$`C=0.31`$ fm2,$`\mathrm{\Delta }=0.13`$\] taken from  (used in ). A value $`m_T=0.4`$ GeV is employed by default (in  the nucleon mass $`m_N`$ was used). The sensibility of our results to variations in $`m_T`$ will be examined.

In Fig. 12 our results at $`y^{}=0`$ are presented for AuAu at RHIC energies and for PbPb collisions at the LHC, versus impact parameter. Reductions of multiplicities at $`b=0`$ by factors $`1/2`$ for RHIC and $`1/3`$ for LHC are found, with a clear increase of the suppression with increasing energy. In Fig. 13 results are presented for AuAu at RHIC and PbPb at LHC for different $`y^{}`$. Finally, in Fig. 14 the variation with $`m_T`$ of the results at $`y^{}=0`$ for AuAu at RHIC is studied. A reduction of the suppression with increasing $`m_T`$ is seen, as expected. Let us make two comments: First, our results for the reduction factors are very similar to the ones estimated in . It has been shown in  that, when these reduction factors are used to correct the results of the Dual Parton Model, one obtains a good description of the RHIC data on multiplicities and their evolution with centrality. Thus, our results provide a detailed calculation of these reduction factors which confirms the estimations in . Second, our results are important in studying particle production in heavy ion collisions. In particular, the dependence of the reduction factors on $`m_T`$ gives the variation of shadowing corrections with the $`p_T`$ of the produced particle555For reduction factors based on other mechanisms, see ..

## 5 Conclusions

In this paper, we have used the relation which arises from Reggeon Calculus and the AGK rules, between the diffractive cross section measured in DIS on nucleons and the first contribution (i.e. double scattering) to nuclear shadowing. The next contributions have been estimated using two different methods for unitarization. In this way we have obtained a description of nuclear shadowing, based on the model of  for diffraction, which agrees with the existing experimental data without any fitted parameter. The model is designed for the region of $`x<0.01`$ and $`Q^2<10`$ GeV2, i.e. small $`x`$ and small or moderate $`Q^2`$.

The same method has been applied in . In , a model for diffraction  has been used that takes into account unitarization effects in an effective manner, so the extrapolation to smaller $`x`$ or larger $`W^2`$ is not so reliable as in the full unitarization program followed in ; furthermore, the description of diffraction in the model we use is substantially better due to the inclusion in the fits of new, more precise experimental data. In , a model for diffraction is used in order to obtain an initial condition for DGLAP evolution at $`Q_0^2=4`$ GeV2, so their leading-twist description for nuclear shadowing is not valid at small $`Q^2`$. On the contrary, we develop a model valid for the full low $`Q^2`$ region which does not correspond to any definite twist but contains contributions from all twist orders. Nevertheless, it turns out that, as discussed at the end of Section 3, the leading-twist contribution is the dominant one in our model, which is in reasonable agreement with the existing experimental data. Precise data on the $`Q^2`$-dependence of nuclear structure functions should disentangle between these two possibilities. The existing data from NMC  can be well reproduced within the leading-twist DGLAP evolution  with an appropriate set of initial conditions. An extension of our results using DGLAP evolution for large values of $`Q^2`$ is thus a natural continuation of our work .

In this framework we have also obtained the factor for multiplicity reduction in nucleus-nucleus collisions. This factor reaches values $`1/2`$ and $`1/3`$ for central AuAu and PbPb collisions at RHIC and LHC respectively. It is therefore a very important ingredient for the computation of particle production at these energies which should be taken into account together with other possible effects.

Comparison among models shows differences of a factor 0.6 for the ratio of structure functions Pb/nucleon at $`x=10^5`$ and $`Q^2=3`$ GeV2. These differences have a large impact in the computation of particle production in nuclear collisions at energies of RHIC and LHC. They should be testable in future lepton-ion colliders .

To conclude, the method which we have followed offers a natural link between the measurements of nucleon diffractive structure functions and nuclear shadowing, and between the latter and the suppression of particle production in nuclear collisions. In this way the study of Low $`x`$ Physics at HERA is linked to that of nuclear structure functions at future lepton-ion colliders and with Heavy Ion Physics at RHIC and LHC .

Acknowledgments: The authors express their gratitude to B. Badelek, E. G. Ferreiro, B. Z. Kopeliovich, J. Raufeisen and M. Strikman for useful comments and discussions. They also thank M. Lublinsky for discussions and for providing the results , and C. Pajares for a reading of the manuscript. A. C., A. B. K. and C. A. S. acknowledge financial support by grant INTAS 00-00366, and A. B. K. by grants RFBR 00-15-96786, RFBR 01-02-17383 and DFG 436 RUS 113/721/0-1. N. A. thanks the Institute for Nuclear Theory at the University of Washington for its hospitality and the US Department of Energy for partial support during the completion of this work. J. L.-A. thanks CERN Theory Division for kind hospitality, and Universidad de Córdoba and Ministerio de Cultura y Deporte of Spain (grant AP2001-3333) for financial support. C. A. S. is supported by a Marie Curie Fellowship of the European Community programme TMR (Training and Mobility of Researchers), under the contract number HPMF-CT-2000-01025.

## Figure captions:

Fig. 1. Diagram showing diffractive DIS with the corresponding kinematical variables in the infinite momentum frame (left) and its equivalence in the rest frame of the nucleon (right).

Fig. 2. Diagram showing the equivalence between diffractive DIS and two exchanged amplitudes with a cut between the amplitudes.

Fig. 3. Results of the model using Schwimmer (solid lines) and eikonal (dashed lines) unitarization compared with experimental data versus $`x`$, for the ratios C/D, Ca/D, Pb/D  and Xe/D  (filled circles correspond to the analysis with hadron requirement and open circles to that with electromagnetic cuts, see the experimental paper for more details).

Fig. 4. Results of the model using Schwimmer (open circles) and eikonal (open triangles) unitarization compared with experimental data versus $`A`$, for the ratios Be/C, Al/C, Ca/C, Fe/C, Sn/C and Pb/C  at two fixed values of $`x`$.

Fig. 5. Id. to Fig. 3 but for the ratios He/D, C/D and Ca/D .

Fig. 6. Results of the model using Schwimmer (solid lines) and eikonal (dashed lines) unitarization compared with experimental data versus $`Q^2`$, for the ratio Sn/C  at two fixed values of $`x`$.

Fig. 7. Comparison of the results of our model using Schwimmer (solid lines) and eikonal (dashed lines) unitarization for the ratio Pb/nucleon with other models, versus $`x`$ at fixed $`Q^2=3`$ GeV2. HKM are the results from , Sarcevic from , Bartels from , Frankfurt from  (at $`Q^2=4`$ GeV2), Armesto from  and EKS98 from .

Fig. 8. Results of the model using Schwimmer (solid lines) and eikonal (dashed lines) unitarization for the ratios D/nucleon, He/nucleon, Li/nucleon, C/nucleon, Ca/nucleon, Sn/nucleon and Pb/nucleon versus $`x`$ at $`Q^2=0.5`$, 2 and 5 GeV2.

Fig. 9. Results of the model for $`Q^2=0`$ using Schwimmer (solid lines) and eikonal (dashed lines) unitarization for the ratios D/nucleon, He/nucleon, Li/nucleon, C/nucleon, Ca/nucleon, Sn/nucleon and Pb/nucleon (upper plot), and for different impact parameters $`b`$ for the ratios C/nucleon (plot in the middle) and Pb/nucleon (lower plot), versus $`W^2`$.

Fig. 10. Results for $`Q^2=0.5`$ (upper plot) and 5 (lower plot) GeV2 using Schwimmer unitarization for the ratio Pb/nucleon versus $`x`$, of the model without modifications (solid lines), without the higher-twist contribution in the short-distance component (dashed lines), and without the higher-twist contribution in the short-distance component plus some modification in parameters (dotted lines) to check the sensibility of the results, see text.

Fig. 11. Diagram showing the contribution to particle production in the central region in AB collisions.

Fig. 12. Results of the model for the multiplicity reduction factor versus impact parameter $`b`$ at $`y^{}=0`$, for AuAu collisions at $`\sqrt{s}=19`$, 130 and 200 GeV per nucleon, and for PbPb collisions at $`\sqrt{s}=5500`$ GeV per nucleon, in the parton model-like realization (solid lines) and for \[$`C=0.31`$ fm2,$`\mathrm{\Delta }=0.13`$\] (dashed lines).

Fig. 13. Id. to Fig. 12 but for AuAu collisions at $`\sqrt{s}=200`$ GeV per nucleon and for PbPb collisions at $`\sqrt{s}=5500`$ GeV per nucleon, for $`y^{}=1`$, 2 and 3.

Fig. 14. Results of the model for the multiplicity reduction factor versus impact parameter $`b`$ at $`y^{}=0`$, for AuAu collisions at $`\sqrt{s}=200`$ GeV per nucleon, in the parton model-like realization (upper plot) and for \[$`C=0.31`$ fm2,$`\mathrm{\Delta }=0.13`$\] (lower plot). In each plot, lines from bottom to top correspond to $`m_T^2=0.16`$, 1, 2, 3, 4 and 5 GeV2.

Figures:

Fig. 1

Fig. 2

Fig. 3

Fig. 4

Fig. 5

Fig. 6

Fig. 7

Fig. 8

Fig. 9

Fig. 10

Fig. 11

Fig. 12

Fig. 13

Fig. 14