text:
# An Impact Parameter Dipole Saturation Model

## I Introduction

One of the most important observation of the HERA experiments is the rapid rise of the $`\gamma ^{}p`$ cross-section with the $`\gamma ^{}p`$ center of mass energy $`W`$ in the deep inelastic scattering (DIS) region. This observation indicates that over a region of size $`1/Q^2`$ additional gluons are abundantly radiated at small $`xQ^2/W^2`$. The simultaneous observation of a significant number of diffractive events in DIS lead to the dipole saturation model of Golec-Biernat and Wüsthoff (GBW) ref.GBW . The model successfully described the data in the transition region between low and intermediate $`Q^2`$ and indicated the emergence of saturation phenomena in DIS. In spite of its compelling success the model has some obvious shortcomings. First, the treatment of QCD evolution is only rudimentary. Second, the dipole cross-section is integrated over the transverse coordinate although the gluon density is expected to be a strongly varying function of the impact parameter. The evolution was improved by incorporating DGLAP into the dipole cross-section thereby extending the model applicability to the high $`Q^2`$ region ref.BGK .

In this paper we further improve the dipole model by introducing the impact parameter of the proton into the dipole dynamics. For the form of the dipole cross-section we adopt the Glauber-Mueller formula which can be derived under simplifying assumptions ref.GlMu . The transverse profile of the proton can be probed experimentally by measuring the momentum transfer $`t`$ to the proton in diffractive events. Thus diffractive and inclusive $`\gamma ^{}p`$ processes form a tool to determine the gluon density and saturation dynamics at HERA. Whether the experimental data are sufficient to reveal the rich non-linear dynamics of saturation and the Color Glass Condensate CGC  remains an important phenomenological question ref.MSM ; LevinTotal ; LevinJ  which this work addresses.

The Impact Parameter (IP) Saturated dipole model can be extended to DIS interactions in nuclei in a natural way. We study nuclear effects and compare them with data on nuclear shadowing. Finally we estimate the saturation scale in heavy nuclei.

## II Dipole description of DIS

In the dipole picture dip1 ; dip2 , the $`\gamma ^{}p`$ interaction proceeds in three stages: First, the incoming virtual photon fluctuates into a quark-antiquark pair, then the $`q\overline{q}`$ pair elastically scatters on the proton, and finally the $`q\overline{q}`$ pair recombines to form a virtual photon. This is shown schematically in Fig. 1.

The amplitude for the complete process is simply the product of the amplitudes of these three subprocess. We discuss each in turn.

The amplitude for the incoming virtual photon (with polarization $`\lambda =+,,0`$) to fluctuate into a quark-antiquark pair (of flavor $`f`$ with helicities $`h,\overline{h}=\pm 1/2`$) is given by the photon light cone wave-function $`(\psi _{\gamma ,f})_\lambda ^{h\overline{h}}`$ which depends upon $`(Q^2,z,𝐫)`$. Here $`z`$ denotes the longitudinal momentum fraction of the photon carried by the quark and $`𝐫`$ the transverse size of the $`q\overline{q}`$ pair. These wave functions are determined from light cone perturbation theory to leading order in the fermion electromagnetic charge ref.BjKS . Similarly, the amplitude for the $`q\overline{q}`$ to recombine to form a virtual photon is $`(\psi _{\gamma ,f}^{})_\lambda ^{h\overline{h}}`$.

The elastic scattering of the $`q\overline{q}`$ pair with squared momentum transfer $`\mathrm{\Delta }^2=t`$ is described by the elastic scattering amplitude, $`A_{el}^{q\overline{q}}(x,r,\mathrm{\Delta })`$. (The notation in this section follows S. Munier et al ref.MSM .) Here $`A_{el}`$ is defined such that the elastic $`q\overline{q}`$ cross section is

$$\frac{d\sigma _{q\overline{q}}}{dt}=\frac{1}{16\pi }\left|A_{el}(x,r,\mathrm{\Delta })\right|^2.$$
(1)
It is convenient to work in coordinate space and define the S-matrix element at a particular impact parameter $`b`$

$$S(b)=1+\frac{1}{2}d^2\mathrm{\Delta }e^{i𝐛𝚫}A_{el}^{q\overline{q}}(x,r,\mathrm{\Delta }).$$
(2)
This corresponds to ones intuitive notion of impact parameter only when the dipole size is small compared to the size of the proton. The total cross section for the $`q\overline{q}`$ pair is found by taking the imaginary part of $`iA_{el}`$ or in terms of the S-matrix element

$$\sigma _{q\overline{q}}(x,r)=\text{Im}iA_{el}^{q\overline{q}}(x,r,0)=d^2b\mathrm{\hspace{0.17em}2}(1\text{Re}S(b)).$$
(3)
This motivates the definition of the $`q\overline{q}`$ differential cross section

$$\frac{d\sigma _{qq}}{d^2b}=2(1\text{Re}S(b)).$$
(4)
The total $`\gamma ^{}p`$ cross section for transverse and longitudinal virtual photons is then given by the imaginary part of the forward scattering amplitude

$$\sigma _{L,T}^{\gamma ^{}p}(x,Q^2)=d^2r\frac{dz}{4\pi }\underset{f}{}\left(\gamma |\gamma \right)_{L,T}^f\sigma _{q\overline{q}}(x,r),$$
(5)
with

$$\sigma _{qq}(x,r)=d^2b\frac{d\sigma _{qq}}{d^2b}.$$
(6)
$`(\gamma |\gamma )_{L,T}^f`$ denotes the probability for a polarized virtual photon to fluctuate into a $`q\overline{q}`$ pair with flavor $`f`$

$`(\gamma |\gamma )_L^f`$ $``$ $`{\displaystyle \underset{h\overline{h}}{}}(\psi _{\gamma ,f}^{})_0^{h\overline{h}}(\psi _{\gamma ,f})_0^{h\overline{h}}`$ (7)
$`(\gamma |\gamma )_T^f`$ $``$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\stackrel{h\overline{h}}{\lambda =\pm 1}}{}}(\psi _{\gamma ,f}^{})_\lambda ^{h\overline{h}}(\psi _{\gamma ,f})_\lambda ^{h\overline{h}}.`$ (8)
In the dipole picture elastic diffractive vector meson production appears in a similarly transparent way BrodskyVM . We denote amplitude for a vector meson to fluctuate into a $`q\overline{q}`$ pair by $`(\psi _V)_\lambda ^{h\overline{h}}`$. This vector-meson wave function will be discussed below. The amplitude for elastic diffractive vector meson production with squared momentum transfer $`\mathrm{\Delta }^2=t`$ and transverse (T) or longitudinal polarization (L) is

$$A_{L,T}(\mathrm{\Delta })=d^2r\frac{dz}{4\pi }d^2b(V|\gamma )_{L,T}e^{i𝐛𝚫}\mathrm{\hspace{0.17em}2}(1S(b)).$$
(9)
Assuming that the S-matrix element is predominantly real we may substitute $`2(1S(b))`$ with $`\frac{d^2\sigma _{q\overline{q}}}{d^2b}`$. Then, the elastic diffractive cross section is

$$\frac{d\sigma _{L,T}}{dt}=\frac{1}{16\pi }|d^2𝐫\frac{dz}{4\pi }d^2b(V|\gamma )_{L,T}e^{i𝐛𝚫}\frac{d\sigma _{qq}}{d^2b}|^2.$$
(10)
Note, it is the total dipole cross section at a particular impact parameter which appears in this formula. S. Munier et al ref.MSM  extracted $`S(b)`$ from the diffractive data using the formalism outlined above.

## III Model Description

In order to describe the data it is necessary to make specific assumptions about the dipole cross section and the photon and vector meson wave-functions. First, we give a simple model for the dipole cross section and contrast this model with other models that have appeared previously. Then we define the gluon structure function $`xg(x,\mu ^2)`$ and the proton shape function $`T(b)`$ which enters the dipole cross-section. Then we define and discuss the properties of the vector meson wave-functions in some detail. Finally after the dipole cross-section and the wave-functions are specified, we determine the model parameters by fitting the inclusive total cross section and elastic diffractive vector meson data.

### III.1 The Dipole Cross Section

The total cross section for a small $`q\overline{q}`$ dipole to pass through a dilute gluon cloud is proportional to the dipole area, the strong coupling constant, and the number of gluons in the cloud ref.FRS

$`\sigma _{q\overline{q}}={\displaystyle \frac{\pi ^2}{N_c}}r^2\alpha _s(\mu ^2)xg(x,\mu ^2),`$ (11)
where $`xg(x,\mu ^2)`$ is gluon density at some scale $`\mu ^2`$. Now imagine that the density of gluons in the target is not small. Divide the target into thin slices of thickness $`dz`$ 111$`z`$ denotes the longitudinal direction and not the quark light cone momentum fraction in this paragraph.. The probability that a dipole at impact parameter $`b`$ does not suffer an inelastic interaction passing through one slice of the proton is

$$P(b)=1\frac{\pi ^2}{N_c}r^2\alpha _s(\mu ^2)xg(x,\mu ^2)\rho (b,z)dz.$$
(12)
Here $`\rho (b,z)`$ denotes density of the gluons within a proton and is normalized to one

$`{\displaystyle d^2b𝑑z\rho (b,z)}=1.`$ (13)
Exponentiating this result, the probability that the dipole suffers no inelastic interactions passing through the entire proton is $`|S(b)|^2`$

$`|S(b)|^2`$ $`=`$ $`\mathrm{exp}\left({\displaystyle \frac{\pi ^2}{N_c}}r^2\alpha _s(\mu ^2)xg(x,\mu ^2)T(b)\right).`$ (14)
Here $`T(b)`$ is the thickness function,

$$T(b)=_{\mathrm{}}^{\mathrm{}}𝑑z\rho (b,z).$$
(15)
To obtain the total cross section at a given impact parameter we assume that the S-matrix element is predominantly real. Then the cross section at a given impact parameter $`b`$ is given by $`2(1\text{Re}S(b))`$ or

$$\frac{d\sigma _{q\overline{q}}}{d^2b}=2\left[1\mathrm{exp}\left(\frac{\pi ^2}{2N_C}r^2\alpha _s(\mu ^2)xg(x,\mu ^2)T(b)\right)\right].$$
(16)
This expression is known as the Glauber-Mueller dipole cross section ref.GlMu  and can also be obtained within the McLerran-Venugopalan model ref.McLV . It provides a simple model for the interaction of a $`q\overline{q}`$ probe with a dense target. The function $`T(b)`$ and $`xg(x,\mu ^2)`$ are determined from the fits to the data. Figure 2 shows the dipole cross

section as a function of the impact parameter $`b`$ for various dipole sizes at two distinct $`x`$-values. The properties of the dipole cross-section will be discussed below. Similar dipole cross sections were found in the model of Ref. Steffen . The cross-section is evaluated for the largest dipole sizes only to illustrate saturation in the Glauber-Mueller approach. This dipole size is not contributing to the observed cross-section.

The Glauber-Mueller dipole cross section is a simple model which respects unitarity. As discussed below, in the large $`Q^2`$ region the model matches smoothly onto DGLAP since we will evolve $`xg(x,\mu ^2)`$ with DGLAP. In the small $`x`$ region, the model evolution can be improved. Indeed, evolution based upon the Balitsky-Kochegov equation BALITSKY ; KOV  has been extensively studied LevinTotal ; LevinJ  and provides a reasonable description of the total and total diffractive cross sections. Nevertheless, the Glauber-Mueller dipole cross section provides a simple model which exhibits saturation as a function of impact parameter in a reasonably generic way and is therefore a useful tool in exposing saturation effects in the data.

### III.2 The Proton Shape $`T(b)`$

The impact parameter dependence of the dipole cross-section determines the properties of the $`t`$distributions in elastic diffractive processes. Here the variable $`t`$ denotes the squared four momentum transfer between the incoming and outgoing proton in the diffractive process, $`\gamma ^{}pVp`$. The $`t`$distribution is directly related to the Fourier transform of $`T(b)`$ via Eqs. 16 and 10. Presently, the best data on diffractive $`t`$-distributions are available for vector-meson production. HERA experiments reported ref.HERAVM ; ref.H1JP ; ref.ZEUSJP  that all diffractive vector-meson $`t`$-distributions can be described by a single exponent $`B`$ for $`|t|<1`$ GeV2

$`{\displaystyle \frac{d\sigma _{VM}^{\gamma ^{}p}}{dt}}\mathrm{exp}(B|t|).`$ (17)
Since the $`J/\psi `$ wave-function is better known than the $`\rho `$ or $`\varphi `$ wave-functions, we will concentrate on $`J/\psi `$ production. The value of the coefficient $`B`$ for the $`J/\psi `$ $`t`$-distribution is around 4 GeV-2 and is only weakly dependent $`W^2`$ and $`Q^2`$. The exponential form in Eq. 17 is the Fourier transform of

$$T_G(b)=\frac{1}{2\pi B_G}\mathrm{exp}(b^2/2B_G).$$
(18)
We therefore tried this form for $`T(b)`$. To account for the expected exponential behavior at large impact parameters we also considered the form

$$T_{GY}(b)d^2b^{}\mathrm{exp}((𝐛𝐛^{})^2/2B_{GY})K_0(b^{}/W_{GY}).$$
(19)
Due to the asymptotic behavior of $`K_0(x)\sqrt{\pi /2x}\mathrm{exp}(x)`$, $`T_{GY(}b)`$ falls off exponentially at large impact parameters.

The parameters $`B_G`$, $`B_{GY}`$ and $`W_{GY}`$ are only weakly dependent on the other parameters in the model and were determined iteratively. From the fit to the limited $`t`$ range, $`|t|<0.7`$ GeV2, the following values were found: $`B_G=4.25`$ GeV-2 for $`T_G(b)`$ and $`B_{GY}=3.25`$ GeV-2 and $`W_{GY}=0.66`$ GeV-1 for $`T_{GY}(b)`$. Figure 3 shows a representative subset of data ref.ZEUSJP  together with the results of a fit to the data using both forms of $`T(b)`$. The slope of the data fixes the parameters $`B_G`$ and $`B_{GY}`$ and $`W_{GY}`$. The normalization of the curves in Figure 3 is a determined by the underlying $`J/\psi `$ light cone wave function and is therefore model dependent. We will use three model wavefunctions which are discussed more fully in Sect. III.6 and in the Appendix. The figure is for the CORNELL wavefunction. Unless explicitly stated, the profile function adopted throughout the rest of this work is $`T_{GY}(b)`$ which for practical purposes is equivalent to $`T_G(b)`$.

### III.3 Gluon Structure Function and Dipole Evolution

The dipole cross-section (Eq. 16) requires the gluon density $`xg(x,\mu ^2)`$ for all scales $`\mu ^2`$. As in previous work ref.BGK , the scale $`\mu ^2`$ is related to the dipole size by

$`\mu ^2={\displaystyle \frac{C}{r^2}}+\mu _0^2.`$ (20)
The gluon density is evolved to the scale $`\mu ^2`$ with leading order DGLAP without quarks for the gluon density. The evolution of the gluon density is performed numerically for every dipole size $`r`$ during the integration of Eq. 5. The initial gluon distribution is taken at the scale $`\mu _0^2`$ in the form

$`xg(x,\mu _0^2)=A_gx^{\lambda _g}(1x)^{5.6}.`$ (21)
The parameters of the model $`C`$, $`\mu _0^2`$, $`A_g`$, and $`\lambda _g`$ are determined from a fit to the DIS data. For light quarks the gluon distribution is evaluated at $`x=x_{BJ}=\frac{Q^2}{Q^2+W^2}`$. For the charm quark the gluon structure function is evaluated at $`x=(4m_c^2+Q^2)/(Q^2+W^2)`$, where $`m_c`$ is the charm quark mass Kretzer . This implies that the dipole cross section is flavor dependent when performing the sum over flavors in Eq. 5.

The present model has, in some sense, one less parameter than the model of Ref. ref.BGK . The parameter $`\sigma _0`$ which previously determined the asymptotic size of the dipole cross-section is now replaced by the parameter $`B_G`$ (or $`B_{GY}`$ and $`W_{GY}`$) which is (are) fixed by the shape of the diffractive $`t`$-distribution.

The gluon structure function obtained from the fit is shown in Fig. 4 for various dipole sizes. The dipole size determines the evolution scale $`\mu ^2`$ through Eq. 20. QCD evolution changes the rate of rise of the gluon structure function $`xg(x,\mu ^2)`$. The rate of rise is usually quantified by the exponent $`\lambda _{eff}`$,

$`\lambda _{eff}={\displaystyle \frac{d\mathrm{log}(xg(x,\mu ^2))}{d\mathrm{log}(1/x)}}.`$ (22)
The exponent $`\lambda _{eff}`$ is a measure of the strength of gluon emission process, see e.g. ref.BK . Figure 5 shows a rapid increase of the exponent $`\lambda _{eff}`$ with decreasing $`x`$ for small values of the dipole radius $`r`$.

### III.4 Comparison with the GBW model

In the original GBW model ref.GBW  the dipole cross section $`\sigma _{q\overline{q}}(x,r)`$ was independent of impact parameter and gave a good description of the inclusive total $`\gamma ^{}p`$ cross section. It also described some diffractive processes ref.GBW ; ref.BGK ; ref.CaldSo  without providing any information about the $`t`$-dependence. In the GBW model the dipole cross-section $`\sigma _{qq}`$ is given by

$$\sigma _{q\overline{q}}^{GBW}=\sigma _0\left[1\mathrm{exp}\left(\frac{r^2}{R_0^2}\right)\right],$$
(23)
where $`R_0`$ denotes the saturation radius and $`\sigma _0`$ is a constant. The saturation radius is $`x`$ dependent

$$R_0^2=\frac{1}{\mathrm{GeV}^2}\left(\frac{x}{x_0}\right)^{\lambda _{GBW}}.$$
(24)
The parameters $`\sigma _0=23`$ mb, $`\lambda _{GBW}=0.29`$ and $`x_0=310^4`$ were determined from a fit to the data. The saturation radius $`R_0`$ is analogous to the gluon distribution and determines the growth of the total and diffractive cross sections with decreasing $`x`$. For large size dipoles $`\sigma _{qq}^{GBW}`$ saturates by approaching a constant value $`\sigma _0`$ which is independent of the exponent $`\lambda _{GBW}`$.

Once the impact parameter is included saturation becomes $`b`$ dependent. For small values of $`b`$, the cross-section grows rapidly with $`r`$ until it reaches the saturation plateau seen in Fig. 2, $`d\sigma _{qq}/d^2b=2`$. The extension of the plateau $`b_S`$ can be estimated from the condition $`\frac{\pi ^2}{2N_C}r^2\alpha _s(\mu ^2)xg(x,\mu ^2)T(b_S)=1`$. Using the Gaussian form of $`T(b)`$ and approximating $`xg(x,\mu ^2)(x_0/x)^{\lambda _{eff}}`$ we obtain

$$b_S^2=2B_G\mathrm{log}\left(\frac{\pi }{4N_CB_G}r^2\alpha _s(\mu ^2)(x_0/x)^{\lambda _{eff}}\right).$$
(25)
Thus, the integrated dipole cross-section $`\sigma _{qq}`$ increases logarithmically for large $`r`$ and small $`x`$ due to the growth of the plateau region. The logarithmic growth of the cross-section at small $`x`$ is determined by $`\lambda _{eff}`$ and $`B_G`$. This behavior should be compared to the GBW cross section which becomes a constant $`\sigma _0`$ that is independent of $`\lambda _{GBW}`$. The differences are illustrated in Fig. 6.

### III.5 Photon wave-function

The longitudinal and transverse photon wave-function in the conventions of Ref. ref.BL  are given by

$`(\psi _{\gamma ,f})_0^{h\overline{h}}(r,z)`$ $`=`$ $`e_f\delta _{f\overline{f}}\sqrt{N_c}\left[2z(1z)Q\left(\delta _{h+\overline{h}}+\delta _{h\overline{h}+}\right)\right]{\displaystyle \frac{K_0(ϵr)}{2\pi }}`$
$`(\psi _{\gamma ,f})_{+1}^{h\overline{h}}(r,z)`$ $`=`$ $`e_f\delta _{f\overline{f}}\sqrt{2N_c}\left[ie^{+i\phi _r}\left(z\delta _{h+\overline{h}}(1z)\delta _{h\overline{h}+}\right)_r+m_f\delta _{h+\overline{h}+}\right]{\displaystyle \frac{K_0(ϵr)}{2\pi }}`$ (26)
$`(\psi _{\gamma ,f})_1^{h\overline{h}}(r,z)`$ $`=`$ $`e_f\delta _{f\overline{f}}\sqrt{2N_c}\left[+ie^{i\phi _r}\left((1z)\delta _{h+\overline{h}}z\delta _{h\overline{h}+}\right)_r+m_f\delta _{h\overline{h}}\right]{\displaystyle \frac{K_0(ϵr)}{2\pi }}`$
where $`re^{i\phi _r}=r_x+ir_y`$, $`ϵ^2=z(1z)Q^2+m_f^2`$, $`e_f`$ is the electromagnetic charge of quark flavor $`f`$, and $`_rK_0(ϵr)=ϵK_1(ϵr)`$.

For small dipoles ($`ϵr1`$) the photon wave-function scales as $`K_1(ϵr)1/ϵr`$. For large dipoles ($`ϵr1`$), the wave-function falls off exponentially with the quark mass, $`K_1(ϵr)\sqrt{\pi /2ϵr}\mathrm{exp}(ϵr)`$.

### III.6 Vector Meson Wave-Function

To describe vector meson production, we developed several models of the $`J/\psi `$ light cone wave-function, $`(\psi _V)_\lambda ^{h\overline{h}}`$. These phenomenological wave-functions should satisfy several model independent constraints FS ; DKP . First, the wave-function must satisfy the normalization condition

$`1`$ $`=`$ $`{\displaystyle \frac{d^2kdz}{16\pi ^3}\underset{h\overline{h}}{}\left|(\stackrel{~}{\psi }_V)_\lambda ^{h\overline{h}}(k,z)\right|^2}.`$ (27)
Second, the decay width of the vector meson further constrains the wave-function. Precise definitions and explicit formulas for the decay width and other constraints on the longitudinal and transverse wave functions are given in the Appendix.

As an example, we model the $`J/\psi `$ wave-function as a Gaussian in momentum space in the meson rest frame

$$\stackrel{~}{\psi }_{NR}(p)=N_{L,T}\mathrm{exp}(\frac{p^2}{\omega _{L,T}^2}),$$
(28)
were $`\stackrel{}{p}`$ denotes the charmed quark 3-momentum. In DIS reactions the meson wave-functions are defined on the light-cone and not in the $`J/\psi `$ rest frame. Following Frankfurt and Strikman FS  we boost the $`J/\psi `$ wave-function from the rest frame to the light-cone. The constants $`N_{L,T}`$ and $`\omega _{L,T}`$ were adjusted to satisfy the normalization and decay constraints. This wave-function is denoted Gaus-RF (Rest Frame).

As noted in Ref. FS  this procedure does not lead to the proper asymptotic form of the light cone wave function at small distances, $`\stackrel{~}{\varphi }z(1z)`$. Therefore, we also considered a Gaussian wave-function defined directly on the light cone

$$\stackrel{~}{\varphi }_{L,T}(k,z)=N_Lz(1z)\mathrm{exp}(\frac{k^2}{\omega _L^2}).$$
(29)
Once again, the norm and decay width constrain $`N_L`$ and $`\omega _L`$. This wave-function is denoted Gaus-LC (Light Cone).

$`J/\psi `$ wave-functions were also studied with potential models Eichten ; Buchmueller ; EichQuig  and successfully describe the energy levels of $`c\overline{c}`$ bound states. Potential models were less successful at reproducing the meson decay width which is sensitive to the short distance behavior of the wave-function. This shortcoming may be less important since the short distance behavior of the wave-function is suppressed by the vanishing dipole cross section at small $`r`$. We determined the non-relativistic wave-function by integrating the radial Schrödinger equation with the Cornell potential Eichten . As before we naively boosted this wave-function to the light-cone. Naturally, this wave-function is denoted CORNELL.

In order to visualize the meson wave-function we define the meson photon overlap $`p_{\gamma ^{}V}(r)`$ as

$$p_{\gamma ^{}V}(r)=r\frac{dz}{4\pi }(\gamma |V)_T,$$
(30)
which has the simple interpretation as the amplitude for diffractive vector meson production from the elastic scattering of a dipole of radius $`r`$. The overlap function is shown in Fig. 7 for the transverse component of the three $`J/\psi `$ wave-functions considered in the paper.

### III.7 Dipole size

To interpret the IP Saturation model physically, we determined the typical dipole size contributing to the total cross-section. The distribution of dipole sizes $`\frac{1}{\sigma ^{\gamma ^{}p}}\frac{d\sigma ^{\gamma ^{}p}}{dr}`$ is the integrand of Eq. 5 after the $`z`$ integration is performed. Figure 8 shows this distribution at representative values of $`x`$ and $`Q^2`$. For each illustrated $`x`$ and $`Q^2`$ the arrow indicates the median radius of a dipole which scatters off the proton. Note that at smallest $`Q^2`$ and $`x=10^5`$ the majority of dipoles are smaller than $`r2`$ GeV-1 and thus half of the dipoles are perturbative.

## IV Results for $`\gamma ^{}p`$ Processes

### IV.1 Inclusive total cross-section

The fit to the total $`\gamma ^{}p`$ cross-section data was performed in the range $`x<0.01`$ and $`Q^2>0.25`$ GeV2. The small $`x`$ cut leads to an upper limit on $`Q^2`$, $`Q^2<200`$ GeV2. Figure 9 shows the $`\gamma ^{}p`$ cross section as a

function of $`W^2`$ at the various $`Q^2`$ measured by the HERA experiments ref.sigtot .

Following the rise of the total cross section with $`W^2`$ at fixed $`Q^2`$, we observe a striking change: for small $`Q^2`$ (the top of Fig. 9) the cross section rises slowly with the center of mass energy $`W^2`$, while for large $`Q^2`$ (the bottom of Fig. 9) the cross section rises strongly. If we parameterize the cross section by a power law

$$\sigma _{tot}^{\gamma ^{}p}(W^2)^{\lambda _{tot}(Q^2)}(1/x)^{\lambda _{tot}(Q^2)},$$
(31)
this change in the rate of rise from low to high $`Q^2`$ translates into a $`Q^2`$ dependence of the exponent $`\lambda _{tot}`$. $`\lambda _{tot}(Q^2)`$ is illustrated in Fig. 10 and summarizes one of the principle observations at HERA. At small $`Q^2`$, $`\sigma _{\gamma ^{}p}`$ increases with $`W^2`$ as in all hadronic reactions, $`\lambda _{tot}0.08`$ ref.DonLan . At at large $`Q^2`$ the exponent $`\lambda _{tot}`$ increases substantially.

In the IP Saturation model the increase of the exponent $`\lambda _{tot}`$ with $`Q^2`$ is a consequence of the DGLAP evolution. When the dipole size becomes smaller than $`\frac{2}{\mu _o}`$ DGLAP evolution becomes increasingly effective and the growth of the gluon distribution at small $`x`$ becomes increasingly rapid. Mathematically speaking, $`\lambda _{eff}`$ increases for small size dipoles as seen in Fig. 5. Averaging over dipole sizes with the photon wave function translates $`\lambda _{eff}(r)`$ into $`\lambda _{tot}(Q^2)`$. The large increase in $`\lambda _{tot}`$ with $`Q^2`$ indicates a transition from a predominance of large to a predominance of small dipoles when evaluating the cross-section. The contribution of large dipoles at low $`Q^2`$ is determined by the light quark mass, $`m_{u,d,s}`$.

The parameters $`A_g`$, $`\lambda _g`$ and $`\mu _0^2`$ were found by fitting the total $`\gamma ^{}p`$ cross-section for various light quark masses. The parameter $`C`$ is strongly correlated to the parameters $`A_g`$ and $`\lambda _g`$ and does not affect the fit quality. Therefore, we set $`C=4`$ by fiat. The mass of the charm quark was $`m_c=1.3`$ GeV as required by elastic diffractive $`J/\psi `$ production discussed below. The results of the fits to the ZEUS data are collected in Table 1. Because of the good agreement between the H1 and ZEUS data sets the combined H1-ZEUS fit returns nearly identical parameters. The fits show a clear preference for the light quark mass of $`m_{u,d,s}=50`$ MeV, which is taken as a reference value. Note, the quality of the fit reported here is even better than the fits of Ref. ref.BGK .

As seen in Table 1, when light quark mass is increased, the fit subsequently increases $`\lambda _g`$. DGLAP evolution can only increase the observed rate of rise of the cross section with $`W^2`$, $`\lambda _{tot}>\lambda _g`$. Therefore, as seen in Fig. 9, for larger masses of the light quarks the cross section rises too quickly with $`W^2`$ at low $`Q^2`$. It is certainly known that the strange quark mass is heavier than the $`u`$ and $`d`$ quarks. However, the influence of the strange quark on the overall fit is relatively small since its contribution to $`F_2`$ is suppressed relative to the $`u,d`$ contribution by the ratio of charges squared, $`e_s^2/(e_u^2+e_d^2)=1/5`$. Therefore, we have neglected the difference between the $`s`$ and $`u,d`$ quark masses. Including a heavier strange quark would make the fit only marginally worse. However, the clear dependence of the fit quality on the mass of the light quark is a property of the model. It is due to the interplay of DGLAP evolution and the properties of the photon wave function. It may indicate the shortcomings of the wave function at large distances.

In contrast to the GBW model ref.GBW  and the model of Ref. ref.BGK , saturation is not responsible for the change in $`\lambda _{tot}`$ with $`Q^2`$ once the impact parameter is included. This is because the total dipole cross section continues to grow within the saturation region although the differential dipole cross-section has attained its asymptotic value $`\frac{d\sigma }{d^2b}=2`$. Therefore, for the total cross-section saturation effects are smaller than in the original GBW model. However, for diffractive $`t`$-distributions saturation will appear more clearly.

### IV.2 Charm Production

The dipole model makes a direct prediction for the inclusive charm contribution to $`F_2`$. In the dipole approach the charm quark distribution is calculated from the gluon distribution. Figure 11 shows a comparison between the measured and predicted values of the charm structure function $`F_2^c`$. The results depend weakly on the charm mass. The good agreement with data for both ZEUS and H1 experiments ref.ZEUSCH ; ref.H1CH  confirms the consistency of the model and supports the dipole picture.

### IV.3 Determination of $`F_L`$

As a byproduct of our investigation we determined the longitudinal structure function $`F_L`$. A measurement of $`F_L`$ would be a critical test of our understanding of the gluon distribution. $`F_L`$ is zero in the naive parton model and is non-zero only after gluon radiation is included. In dipole models at small $`x`$, both $`F_L`$ and $`F_2`$ are governed by the gluon distribution and therefore they behave similarly. Figure 12 shows the ratio of $`F_L`$ to $`F_2`$ as a function of $`x`$ for various photon virtualities. The ratio is around 18 $`\%`$ and gently varies with $`x`$. In the standard QCD fits $`F_L/F_2`$ may show a stronger variation with $`x`$ at low $`Q^2`$ RThorn . This is because the shape of the gluon structure function seen in Fig. 13 is very different from the dipole model presented here. The difference in the gluon structure function between the MRST fit MRST  and the dipole model fit is due to the sea quark contribution. In the MRST fit the initial sea quark contribution is independent of the gluon distribution and subsequently influences the gluon evolution. In the dipole fit the sea quark distribution is always derived from the gluon distribution. It is important to measure $`F_L`$ in order to determine the gluon density precisely and to clarify the relative roles of gluons and sea quarks.

### IV.4 $`J/\psi `$ production

The cross sections for elastic diffractive $`J/\psi `$ production were calculated using Eq. 10 with the dipole cross section determined by fitting $`\sigma _{\gamma ^{}p}`$. Figure 14 compares the computed total $`J/\psi `$ cross-section with H1 measurements ref.H1JP  for various photon virtualities at $`W=90`$ GeV. The absolute magnitude of the cross-sections is very sensitive to the charm quark mass. The charm quark mass was adjusted to $`m_c=1.3`$ GeV in order to match the experimental cross section at $`Q^2=0`$. A small shift of 50 MeV in the charm mass changes the cross-section by around 15$`\%`$.

The predicted total $`J/\psi `$ cross-sections for the Gaus-RF and Gaus-LC wave-functions are very similar in the $`Q^2`$ range considered. Therefore, we only show the predictions of the Gaus-RF wave-function in what follows. The CORNELL wave-function, which is more concentrated at the origin leads to a prediction which is in significantly better agreement with data at large $`Q^2`$, see Fig. 14. In contrast to Ref. FS , the scale $`\mu ^2=C/r^2+\mu _0^2`$ is fixed with $`C=4`$. Within the range $`C=110`$ the results were almost independent of the parameter $`C`$ after refitting total $`\gamma ^{}p`$ cross-section.

Figure 15 compares model predictions with the measured $`J/\psi `$ cross-section as a function of energy $`W`$ at $`Q^2=0`$. In this case, both $`J/\psi `$ wave-functions are in good agreement with the data.

Note, there is a $`15\%`$ discrepancy in the normalization of the $`J/\psi `$ cross-sections between ZEUS and H1. This discrepancy was absorbed in the charm quark mass; we take $`m_c=1.25`$ GeV when comparing with ZEUS data and $`m_c=1.3`$ GeV when comparing with H1.

### IV.5 $`t`$ distributions for $`J/\psi `$ production

Figure 16 compares the $`J/\psi `$ diffractive cross-section as a function of $`t`$ for various $`W^2`$ to model predictions. In the measured region, both $`J/\psi `$ wave-functions lead to very similar results although at large $`t`$ they start to differ significantly. In the large $`t`$ region the exponential behavior of the differential cross section is considerably distorted and the IP Saturation model predicts diffractive dips. Figure 17 illustrates the emergence of the dip by displaying the Fourier transform of the differential dipole cross section for various dipole sizes

$$\frac{d\sigma _{qq}}{dt}=\left|d^2be^{i𝐛𝚫}\frac{d\sigma _{qq}}{d^2b}\right|^2.$$
(32)
For small $`r`$, the $`b`$ dependence of the dipole cross section is essentially Gaussian and its Fourier transform is a pure exponential in $`t`$. For large $`r`$, due to saturation, the shape of the dipole cross section becomes more box-like, as seen in Fig. 2. The two dimensional Fourier transform of a box is proportional to the Bessel function $`J_1`$ which leads to dips in the $`t`$ distribution. Although the distortions to the proton shape seen in Fig. 2 are small when displayed in coordinate space, they appear vividly after taking the Fourier transform and examining the $`t`$-distribution. For dipoles with $`r3`$ GeV-1 the $`t`$-distribution is substantially changed.

### IV.6 Saturation versus Non-Saturation

In order to isolate the effects of saturation we compared the results of the IP Saturation model with a model which is not unitarized. This model, denoted IP Non-Sat, replaces the dipole cross-section of Eq. 11 with the first term after expanding of the exponential

$$\frac{d\sigma _{q\overline{q}}}{d^2b}=\frac{\pi ^2}{N_C}r^2\alpha _s(\mu ^2)xg(x,\mu ^2)T(b).$$
(33)
The proton shape $`T(b)`$ was assumed to be the same as in the saturated case. The fit to the total cross section was repeated and a non-saturated gluon structure function was determined. The optimal fit was obtained for a higher quark mass than in the saturated case, $`m_{u,d,s}=150`$ MeV. This is expected since the dipole cross-section $`\sigma _{q\overline{q}}`$ grows faster with $`r`$ than in the saturated case. The fit yields $`A_g=3.5`$, $`\lambda _g=0.13`$, $`\mu _0^2=0.80`$ and $`\chi ^2/N_{DF}=0.8`$. Naturally, the saturated and non-saturated dipole cross sections are very similar since they both fit the inclusive $`\gamma ^{}p`$ cross-sections very well. Consequently, all observable cross-section which only are sensitive to the dipole cross-section integrated over impact parameter are indistinguishable. This was verified explicitly for the total diffractive $`J/\psi `$ production. In particular, this agreement can be seen in Fig. 10 for the $`\lambda _{tot}(Q^2)`$ distribution. However, the gluon distributions are not the same as illustrated in Fig. 18 which compares the saturated and non-saturated models. For small dipole sizes the distributions are similar but they differ significantly as $`r`$ becomes large.

The $`t`$-distribution provides additional insight into the dynamics of saturation in DIS. Figure 16 compares predictions for diffractive $`J/\psi `$ production in the saturated and non-saturated modes . As discussed above saturation distorts the impact parameter distribution and leads unavoidably to diffractive dips in the $`t`$-distribution. Without saturation, as seen in Fig. 16, the exponential decrease continues to large $`t`$. Diffractive dips in the $`t`$ distribution were clearly observed in elastic hadronic reactions Goulianos . In DIS processes at HERA the measured range of the $`t`$-distribution is too small for a dip to appear. In addition the HERA experiments have always measured elastic diffractive events with some small admixture of inelastic diffractive events. Even if inelastic events are a small fraction of the total, these events will dominate the $`t`$-distribution at large $`t`$ since the inelastic $`t`$-distribution falls slowly ref.ZEUSJP ; Goulianos .

In the small $`t`$ range, saturation slightly increases the slope parameter $`B`$ of diffractive vector meson production at large $`W`$. Figure 19 compares

the slope parameter $`B`$ determined by the ZEUS experiment with model predictions. The observed increase of the slope $`B`$ with $`W`$ (called shrinkage) is slightly larger than the model prediction. This could be an indication that the shape $`T(b)`$ is somewhat dependent on $`x`$ due to QCD evolution. Since considerably more data is available than was used in the present analysis we plan further study of the observed shrinkage.

Concluding we note that the $`t`$-distribution provides precise information about the interplay of saturation and evolution phenomena. Therefore, it is crucially important to measure elastic vector meson production with a small systematic uncertainty over a wide range in $`t`$.

## V Saturation Scale $`Q_S^2`$

In the limit of high gluon density called the Color Glass Condensate (CGC)CGC , the interaction probability in DIS becomes large. The probability of no interaction is measured by the square of the S-matrix element, $`S^2(b)`$. Normally, in high energy reactions the interaction probability is small, $`S^21`$. However, in the CGC limit this probability approaches unity, $`S^20`$. Figure 20 shows $`S^2(b=0)`$ as a

function of $`x`$ for various dipole sizes. For small size dipoles the gluon structure function grows strongly with decreasing $`x`$ and the S-matrix diminishes quickly. For large size dipoles, the gluon density density does not grow and the S-matrix is small even at relatively large $`x`$.

To delineate the saturation region in the ($`x`$,$`r^2`$) plane we first define a density profile $`D`$

$$D(b,x,r^2)=\frac{2\pi ^2}{N_C}\alpha _s(\mu ^2(r^2))xg(x,\mu ^2(r^2))T(b).$$
(34)
Then, the saturation radius $`r_S`$ is defined as the dipole size for which proton consists of one interaction length

$$S^2=\mathrm{exp}\left(\frac{D(b,x,r_S^2)r_S^2}{2}\right)=e^1.$$
(35)
This defines an implicit equation for $`r_S^2`$ which can be solved numerically. Now the saturation scale $`Q_S^2`$ is defined as the density profile $`D`$ evaluated at the saturation radius

$$Q_S^2(b,x)=D(b,x,r_S^2)=\frac{2}{r_S^2}.$$
(36)
In the GBW model $`S^2=\mathrm{exp}(r^2/2R_0^2)`$ and we find $`Q_S^2=1/R_0^2`$ in accord with the original definition of the saturation scale. The saturation region is then defined by the condition that $`r>r_S(b,x)`$.

Figure 21 shows the saturation scale $`Q_S^2`$ as a function of $`1/x`$ in the central impact parameter region. We observe that the saturation scale is lower than the GBW value. This is primarily because the charm quark contribution to $`F_2`$ was neglected when the saturation scale was evaluated. Including the charm quark reduces the gluon density and the saturation scale.

We also calculated the S-Matrix element and determined the corresponding effective saturation scale for the non-saturated model. As seen in Fig. 21, the saturation scale is very similar for IP Saturated and IP Non-Saturated models. This is because the corresponding gluon structure functions are very similar for small size dipoles, as seen in Fig. 18. Therefore, our estimates of the saturation scale are only weakly model dependent.

Table 2 records the saturation

radii and the percentage of the total $`\gamma ^{}p`$ cross-section in the saturation region, $`r>r_S`$. At low $`Q^2`$ ($`Q^2=0.4\text{GeV}^2)`$ the cross-section is dominated by dipoles which are larger than the saturation radius. As seen in Fig. 8, a fraction of these dipoles may be considered perturbative since they are smaller than $`r2\text{GeV}^1`$. In the large $`Q^2`$ region saturation is of secondary importance. The region $`Q^24`$ GeV2 is particularly interesting since a significant fraction of the dipoles are perturbative and are approximately the size of the saturation radius. BFKL and saturation dynamics are expected to influence observable cross sections in this kinematic domain.

It is also interesting to define the saturation exponent $`\lambda _S`$

$`\lambda _S={\displaystyle \frac{d\mathrm{log}(Q_S^2)}{d\mathrm{log}(1/x)}}.`$ (37)
Figure 22 shows the saturation exponent $`\lambda _S`$ as a function of $`1/x`$ for various impact parameters. The saturation exponent exhibits a sizable $`x`$ dependence. It varies between $`\lambda _S=0.15`$ at $`x=10^2`$ and $`\lambda _S=0.35`$ at $`x=10^5`$ which should be contrasted with the constant GBW value, $`\lambda _{GBW}=0.29`$. An analysis of BFKL evolution near the saturation boundary found that the exponent $`\lambda _S`$ is nearly constant as a function of $`\mathrm{log}(1/x)`$ in the HERA kinematic window Triant .

## VI Nuclear Dependence

Now that impact parameter distribution of the proton is taken into account, the dipole picture provides a framework for modeling nuclear phenomena in DIS. The distribution of nucleons in the nucleus $`\rho _A(\stackrel{}{r})`$ is parametrized according to the Woods-Saxon distribution Bohr

$$\rho _A(\stackrel{}{r})=\frac{N}{\mathrm{exp}\left(\frac{(rR_A)}{\delta }\right)+1},$$
(38)
where $`A`$ denotes the number of nucleons in the nucleus, $`\delta =0.54\text{fm}`$, $`R_A=(1.12\text{fm})A^{1/3}(0.86\text{fm})A^{1/3}`$ and $`N`$ is adjusted to normalize the distribution to one

$$d^3\stackrel{}{r}\rho _A(\stackrel{}{r})=1.$$
(39)
The transverse distribution is defined as

$`T_A(𝐛)={\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑z\rho _A(\sqrt{𝐛^\mathrm{𝟐}+z^2}).`$ (40)
Figure 23 compares the proton shape $`T_p(b)`$ with

the transverse density $`T_A(b)`$ for several light nuclei.

For a first estimate of nuclear effects we assume that the gluons are distributed uniformly throughout the nucleus. This assumption is realized by replacing the proton shape in Eq. 16 by the integrated Woods-Saxon distribution $`T_A(b)`$ scaled by the number of nucleons,

$$\frac{d\sigma _{q\overline{q}}}{d^2b}=2\left[1\mathrm{exp}\left(r^2\frac{\pi ^2}{2N_c}\alpha _s(\mu ^2)xg(x,\mu ^2)AT_A(b)\right)\right].$$
(41)
This naive assumption, called the smooth nucleus assumption below, cannot be true since Fig. 23 shows that the center of the proton is on average more dense than the carbon nucleus. Nevertheless it is instructive to evaluate $`F_2^A/AF_2^p`$ for this case and to compare it with NMC measurements of shadowing at small-$`x`$ NMC.dat , see Fig. 24.

The figure shows that the smooth nucleus assumption underestimates nuclear shadowing considerably as can be anticipated by comparing the densities in Fig. 23.

The proton profile $`T_p(b)`$ indicates that the gluons are distributed within individual nucleons with a radius of $`0.5\text{fm}`$. The nucleons are then distributed within the nucleus according to a Woods-Saxon distribution. Therefore the gluons form a lumpy distribution within the nucleus. To account for this correlation among the gluons we generalize the model in the following way. The probability of a given configuration of nucleons within the nucleus $`\{(𝐛_\mathrm{𝟏},z_1),(𝐛_\mathrm{𝟐},z_2),\mathrm{},(𝐛_𝐀,z_A)\}`$ is

$$P_{\{(𝐛_\mathrm{𝟏},z_1),(𝐛_\mathrm{𝟐},z_2),\mathrm{},(𝐛_𝐀,z_A)\}}=\mathrm{\Pi }_{i=1}^A\rho (𝐛_𝐢,z_i).$$
(42)
For a given configuration of nucleons the probability that the dipole scatters at impact parameter $`𝐛`$ is

$$\left(\frac{d\sigma _{q\overline{q}}}{d^2b}\right)_{\{(𝐛_\mathrm{𝟏},z_1),(𝐛_\mathrm{𝟐},z_2),\mathrm{},(𝐛_𝐀,z_A)\}}=2\left[1\mathrm{exp}\left(\frac{\pi ^2}{2N_C}r^2\alpha _s(\mu ^2)xg(x,\mu ^2)\underset{i=1}{\overset{A}{}}T_p(𝐛𝐛_𝐢)\right)\right].$$
(43)
Now we may average over all possible configurations of nucleons to determine the cross section

$$\frac{d\sigma _{q\overline{q}}}{d^2b}=d^2b_1𝑑z_1\mathrm{}d^2b_A𝑑z_AP_{\{(𝐛_\mathrm{𝟏},z_1),(𝐛_\mathrm{𝟐},z_2),\mathrm{},(𝐛_𝐀,z_A)\}}\left(\frac{d\sigma _{q\overline{q}}}{d^2b}\right)_{\{(𝐛_\mathrm{𝟏},z_1),(𝐛_\mathrm{𝟐},z_2),\mathrm{}(𝐛_𝐀,z_A)\}}.$$
(44)
Noticing that the integrals factorize we have

$`{\displaystyle \frac{d\sigma _{q\overline{q}}}{d^2b}}`$ $`=`$ $`2[1(1I)^A]`$ (45)
$`I`$ $``$ $`{\displaystyle d^2b^{^{}}T_A(𝐛+𝐛^{^{}})\left[1\mathrm{exp}\left(\frac{\pi ^2}{2N_C}r^2\alpha _s(\mu ^2)xg(x,\mu ^2)T_p(𝐛^{^{}})\right)\right]}.`$ (46)
To simplify this further, we notice that the integral over $`𝐛^{^{}}`$ is different from zero only over the size of the proton which is small compared to the nucleus even for the deuteron. Thus $`T_A(𝐛+𝐛^{^{}})T_A(𝐛)`$. Using the definition of the total cross section of the $`q\overline{q}`$ pair on the proton $`\sigma _{q\overline{q}}^p`$ (Eq. 6 and 16) we have

$`{\displaystyle \frac{d\sigma _{q\overline{q}}}{d^2b}}`$ $`=`$ $`2\left[1(1{\displaystyle \frac{T_A(𝐛)}{2}}\sigma _{q\overline{q}}^p)^A\right].`$ (47)
Only if A is large and $`r`$ is small do we recover the smooth nucleus formula of Eq. 41. Similar considerations were discussed in the context of setting up initial conditions for classical QCD simulations of nucleus-nucleus collisions KNV2 . Figure 24 shows that the lumpy nucleus model slightly overestimates the observed shadowing. However, good overall agreement with data clearly indicates that the IP Saturation model can be used to study of nuclear effects.

To estimate the saturation scale for heavy nuclei we follow the same procedure as for proton. We set the lumpy nucleus S-matrix element squared, $`|S|^2=(1\frac{T_A(𝐛=0)}{2}\sigma _{q\overline{q}}^p)^{2A}`$ equal to $`e^1`$. This is an implicit equation for the saturation radius, $`r_S`$. Once $`r_S`$ is determined numerically the saturation scale is determined from the relation, $`Q_S^2=2/r_S^2`$.

Figure 25 shows the saturation scale for various nuclei calculated in this fashion. It is interesting to observe that the saturation scale is not increasing like $`A^{1/3}`$. Rather, for Au the nuclear enhancement factor is only a factor of $`2`$ at $`x=10^2`$ and $`3`$ at $`x=10^4`$. When comparing the saturation scale to other estimates it is important to remember that this is the saturation scale for quarks. The saturation scale for gluons is $`(Q_s^2)_g=\frac{N_c}{C_F}(Q_s^2)_q2.25(Q_s^2)_q`$ am-qm . Thus for RHIC conditions with $`x10^2`$, $`(Q_s^2)_g`$ is $`1.3\text{GeV}^2`$ which is slightly smaller than was estimated previously by analyzing RHIC multiplicity data am-qm ; KZN ; KNV ; KZL .

## VII Summary

We developed an Impact Parameter dipole saturation model which describes HERA the total $`\gamma ^{}p`$ cross section and differential diffractive $`J/\psi `$ with high precision. The new feature of the model is an impact parameter dependent dipole cross-section which respects unitarity constraints. The impact parameter distribution is determined by measured $`J/\psi `$ $`t`$-distributions to a large extent. For the functional form of the dipole cross-section we adopt the Glauber-Mueller formula ref.GlMu  which can also be obtained within the McLerran-Venugopalan model ref.McLV .

The parameters of the dipole cross-section are determined from a fit to HERA data. The fit depends on the specific form of photon and the $`J/\psi `$ wave function. At short distances, the photon wave-function is well known from QED. For the $`J/\psi `$ wave-function we proposed several ansatzes and verified that the dipole cross-section was almost independent of these alternatives.

QCD evolution effects are taken into account by evolving the gluon structure function with the lowest order DGLAP equations. DGLAP evolution is essential in contrast to the GBW model.

In the dipole approach only the gluon structure function is modeled. All sea-quark contributions to $`F_2`$ are obtained from the photon wave-function and the subsequent dipole-gluon scattering. The validity of this approximation is corroborated by the excellent agreement between the predicted and observed open charm contribution to $`F_2`$.

The impact parameter dependence of the dipole cross section determines the properties of the $`J/\psi `$ total and diffractive cross-sections. The saturation ansatz for the dipole cross section predicts diffractive dips in the $`t`$ distribution which are not necessarily predicted by an ansatz without saturation. For $`J/\psi `$ production the diffractive dips appear in a region of $`t`$ which is not covered by data but which is accessible to future measurements. In addition, the saturation ansatz predicts a slight increase of the diffractive slope $`B`$, which reflects a small increase in the proton size at small $`x`$. The change of the slope $`B`$ with energy, called shrinkage, provides information about the interplay of saturation and evolution phenomena. Since shrinkage was observed for several vector-meson species and since considerably more data exists than was used in the present paper, further study of diffractive vector meson production and their $`t`$-distributions is warranted. We also point out that additional measurements of vector-meson production in a extended $`t`$-range would provide crucial information about the gluon density.

With an accurate determination of the S-matrix element we evaluated the saturation scale $`Q_S^2(x)`$ – the scale at which multiple collisions becomes important. The $`Q_S^2`$ determined in the model is approximately half of the GBW value and is only weakly dependent on the specific form of the dipole cross-section.

In the low $`Q^2`$ region up to $`4`$ GeV2, saturation dynamics contributes significantly to the observed cross-sections. The saturated state is at least partially perturbative since it is frequently initiated by small dipoles. Above $`Q^210`$ GeV2 saturation plays a secondary role. Dipoles are very small and although they are strong gluon emitters the gluon phase-space density remains small.

In order to isolate the effects of saturation, we also studied a dipole model which is not unitarized . This IP Non-Saturated model also obtained a good description of the data. However, even with the non-saturated model the fit returned large gluon densities and unitarity corrections became important in the same region of $`x`$ and $`Q^2`$. Thus, the saturation ansatz is the more consistent approach.

The extension of the model to nuclei is fairly straightforward. We evaluated the nuclear cross-section with two different ansatzes for the distribution of gluons within the nucleus and compared the results to data on nuclear shadowing at low $`x`$. We found that the lumpy nucleus ansatz, in which gluons are distributed around the nucleons in a similar manner as around a free proton, leads to a good description of the data. We then used this ansatz to evaluate the saturation scale in heavy nuclei and found that $`Q_S^2`$ increases more slowly then $`A^{1/3}`$. Nevertheless, for Gold and Uranium $`Q_S^2`$ is enhanced by factor of $`2`$ when compared to the proton. For RHIC conditions the gluon saturation scale was found to be $`Q_s^21.3\text{GeV}^2`$, which is somewhat smaller than was estimated previously from RHIC multiplicity data am-qm .

The success of the model suggests that the saturated state known as the Color Glass Condensate plays an important role at HERA. Saturation is particularly important in the small $`x`$ and low $`Q^2`$ region where the $`\gamma ^{}p`$ total inclusive cross-section exhibits the same universal rate of growth as all known hadronic cross-sections. Revealing the Color Glass Condensate experimentally would represent a fundamental advance towards an understanding of the high energy limit of QCD.

## VIII Acknowledgment

We have profited from helpful discussions with many colleagues: A. Caldwell, J. Jalilian-Marian, D. Kharzeev, D. N. Triantafyllopoulos and R. Venugopalan. Our special thanks go to A. H. Mueller and L. McLerran. One of us (H.K.) thanks the Physics Department of Brookhaven National Laboratory and the Department of Physics of Columbia University for hospitality. Special thanks go to the Nevis Laboratories of Columbia University.

This work was partially supported by the National Science Foundation grant NSF PHY 00-98826 and by the U.S. Department of Energy grant DE-AC02-98CH10886.

## Appendix A The $`J/\psi `$ Wave-Function

The discussion here closely follows DKP  and uses some of the concepts of FS . The light cone wave-function for $`J/\psi `$ can be written as an expansion in the Fock space,

$`|V=\stackrel{~}{\psi }_{c\overline{c}}|c\overline{c}+\stackrel{~}{\psi }_{q\overline{q}g}|q\overline{q}g+\mathrm{}.`$ (48)
Below we neglect higher Fock components in the wave-function and treat the meson as a predominantly a $`c\overline{c}`$ pair. Restoring indices, the $`J/\psi `$ state is given by

$`|V^\lambda ={\displaystyle \frac{d^2kdz}{16\pi ^3}(\stackrel{~}{\psi }_V)_\lambda ^{h\overline{h}}(𝐤,z)\frac{\delta _{ab}}{\sqrt{N_c}}|c^a(𝐤,z,h)\overline{c}^b(𝐤,1z,\overline{h})}`$ (49)
Here $`\lambda =0,\pm 1`$ denotes the polarization of the vector meson, $`h,\overline{h}=\pm `$ denote the helicities of the c-quark and anti-c-quark respectively, and $`a,b`$ are the color indices of the c-quark and anti-c-quark. Here and below we strictly adhere to the light-cone conventions of Ref. ref.BL .

The $`J/\psi `$ wave-function is modeled after the photon wave-function. The longitudinally polarized vector meson ($`\lambda =0`$) is

$`(\stackrel{~}{\psi }_V)_0^{h\overline{h}}(k,z)`$ $`=`$ $`{\displaystyle \frac{\delta _{h+\overline{h}}+\delta _{h\overline{h}+}}{\sqrt{2}}}\stackrel{~}{\varphi }_L(k,z),`$ (50)
while the two transverse polarizations ($`\lambda =\pm 1`$) are given by

$`(\stackrel{~}{\psi }_V)_{+1}^{h\overline{h}}(k,z)`$ $`=`$ $`\left[+{\displaystyle \frac{k}{m_c}}e^{+i\phi _k}(z\delta _{h+\overline{h}}(1z)\delta _{h\overline{h}+})+\delta _{h+\overline{h}+}\right]\stackrel{~}{\varphi }_T(k,z)`$
$`(\stackrel{~}{\psi }_V)_1^{h\overline{h}}(k,z)`$ $`=`$ $`\left[{\displaystyle \frac{k}{m_c}}e^{i\phi _k}((1z)\delta _{h+\overline{h}}z\delta _{h\overline{h}+})+\delta _{h\overline{h}}\right]\stackrel{~}{\varphi }_T(k,z),`$ (51)
where $`ke^{i\phi _k}=k_x+ik_y`$. Once the wave-function in momentum space is specified then Fourier transforms return the wave-function in coordinate space

$$(\psi _V)_\lambda ^{h\overline{h}}(r,z)=\frac{d^2k}{(2\pi )^2}e^{i𝐤𝐫}(\stackrel{~}{\psi }_V)_\lambda ^{h\overline{h}}(k,z).$$
(52)
Note under Fourier transform $`\pm ke^{\pm i\phi _k}\stackrel{~}{\varphi }(k,z)`$ becomes $`ie^{\pm i\phi _r}_r\varphi (r,z)`$. The overlap with the photon wave-function is then

$`(\gamma |V)_L{\displaystyle \underset{h\overline{h}}{}}(\psi _\gamma ^{})_0^{h\overline{h}}(\psi _V)_0^{h\overline{h}}`$ $`=`$ $`e_c\sqrt{2N_c}\left[2z(1z)Q{\displaystyle \frac{K_0(ϵr)}{2\pi }}\varphi _L\right]`$
$`(\gamma |V)_T{\displaystyle \frac{1}{2}}{\displaystyle \underset{\stackrel{h\overline{h}}{\lambda =\pm 1}}{}}(\psi _\gamma ^{})_\lambda ^{h\overline{h}}(\psi _V)_\lambda ^{h\overline{h}}`$ $`=`$ $`e_c\sqrt{2N_c}\left[(z^2+(1z)^2)({\displaystyle \frac{1}{m_c}}{\displaystyle \frac{\varphi _T}{r}}){\displaystyle \frac{ϵK_1(ϵr)}{2\pi }}+\varphi _Tm_c{\displaystyle \frac{K_0(ϵr)}{2\pi }}\right].`$ (53)
Two inputs are used to constrain the wave-function $`\stackrel{~}{\varphi }_{L,T}`$. The first input is the normalization of the wave-function, $`V^\lambda ^{}|V^\lambda =(2\pi )^32p^+\delta (p^+p^+)\delta ^2(𝐩^{}𝐩)\delta ^{\lambda ^{}\lambda }`$. This leads to the requirement that

$`1={\displaystyle \frac{d^2kdz}{16\pi ^3}\underset{h\overline{h}}{}\left|(\stackrel{~}{\psi }_V)_\lambda ^{h\overline{h}}(k,z)\right|^2},`$ (54)
or

$$1=\frac{d^2kdz}{16\pi ^3}\left|\stackrel{~}{\varphi }_L(k,z)\right|^2$$
(55)
for $`\stackrel{~}{\varphi }_L`$ and

$$1=\frac{d^2kdz}{16\pi ^3}\left\{(z^2+(1z)^2)\frac{k^2}{m_c^2}+1\right\}\left|\stackrel{~}{\varphi }_T(k,z)\right|^2$$
(56)
for $`\stackrel{~}{\varphi }_T`$ . Note that the formula (B19b) in DKP  is incorrect by a factor of two for the transverse wave-function. The second input is the leptonic decay width of the vector meson. The decay width of the meson is given by

$`\mathrm{\Gamma }_V`$ $`=`$ $`{\displaystyle \frac{4\pi \alpha _{em}^2}{3}}{\displaystyle \frac{f_V^2}{M_V}},`$ (57)
where $`f_V`$ is given by

$`eM_Vf_Vϵ_\lambda ^\mu `$ $`=`$ $`0\left|J_{em}^\mu \right|V^\lambda .`$ (58)
In this formula $`ϵ_\lambda ^\mu =(ϵ^+,ϵ^{},ϵ_x,ϵ_y)`$ is the polarization vector of the meson and $`e`$ is the charge of the positron. $`ϵ_{\pm 1}^\mu =1/\sqrt{2}(0,0,1,\pm i)`$ for transversely polarized mesons and $`ϵ_0^\mu =(\frac{p^+}{M_V},\frac{M_V}{p^+},0,0)`$ for the longitudinally polarized mesons. Evaluating this matrix element in light cone perturbation with the conventions of ref.BL  gives

$`f_V`$ $`=`$ $`\widehat{e_c}\sqrt{2N_c}\mathrm{\hspace{0.17em}2}{\displaystyle \frac{d^2kdz}{16\pi ^3}\stackrel{~}{\varphi }_L(k,z)}`$
$`f_V`$ $`=`$ $`\widehat{e_c}\sqrt{2N_c}{\displaystyle \frac{m_c}{M_V}}{\displaystyle \frac{d^2kdz}{16\pi ^3}\left\{(z^2+(1z)^2)\frac{k^2}{m_c^2}+1\right\}\frac{\stackrel{~}{\varphi }_T(k,z)}{z(1z)}},`$ (59)
where $`\widehat{e_c}=+\frac{2}{3}`$ (in contrast to $`e_c=\sqrt{4\pi \alpha _{em}}\widehat{e_c}`$ ) is the charge in units of the positron charge.

To go further, we must specify the form of the wave-function. We have studied three model wave-functions, GAUS-RF, CORNELL, and GAUS-LC. Two of the models (GAUS-RF and CORNELL) are motivated by the expected long distance behavior of the wave-function. One of models (GAUS-LC) is motivated by the expected short distance behavior of the wave-function.

For large distances, and for truly heavy vector mesons, the non-relativistic wave-function should give a good description of the light cone wave-function. To make the correspondence between the light cone and non-relativistic wave-functions functions we repeat the discussion in Ref. FS . In the non relativistic limit, terms proportional to $`\left(\frac{k_t}{m_c}\right)^2`$ can be dropped and the the normalization condition Eq. 54 is identical for the longitudinal and transverse wave-functions. Simple kinematic relations give the invariant mass of the the quark antiquark pair $`M_V^2=(p_1+p_2)^2=\frac{m_c^2+k_t^2}{z(1z)}`$. This relation assumes that the quark and antiquark are on shell, which is true for non-relativistic quarks to order $`\left(\frac{k_t}{m_c}\right)^2`$. In the rest frame of the meson (where the quarks have momentum $`\stackrel{}{p}`$ and $`\stackrel{}{p}`$) the invariant mass is $`(p_1+p_2)^2=4(p^2+m_c^2)`$. Comparing these terms gives the identification

$`p^2`$ $``$ $`{\displaystyle \frac{k_t^2+(2z1)^2m_c^2}{4z(1z)}}`$ (60)
$`d^3p`$ $``$ $`{\displaystyle \frac{\sqrt{k_t^2+m_c^2}}{4\left[z(1z)\right]^{3/2}}}d^2k_tdz.`$ (61)
Requiring probability conservation per unit phase space

$$\frac{d^3p}{(2\pi )^3}\left|\stackrel{~}{\psi }_{NR}(p)\right|^2=\frac{d^2k_tdz}{16\pi ^3}\left|\stackrel{~}{\varphi }_{L,T}(k_t,z)\right|^2,$$
(62)
leads to the identification

$$\stackrel{~}{\varphi }_{L,T}(k_t,z)=N_{L,T}\left(\frac{k_t^2+m_c^2}{4z^3(1z)^3}\right)^{\frac{1}{4}}\stackrel{~}{\psi }_{NR}\left(p=\sqrt{\frac{k_t^2+(2z1)^2m_c^2}{4z(1z)}}\right).$$
(63)
Here $`\stackrel{~}{\psi }_{NR}(p)`$ is the non-relativistic wave-function in momentum space

$$\psi _{NR}(\stackrel{}{x})=\frac{d^3p}{(2\pi )^3}e^{i\stackrel{}{p}\stackrel{}{x}}\stackrel{~}{\psi }_{NR}(\stackrel{}{p}).$$
(64)
For the two long distance models (GAUS-RF and CORNELL) $`\stackrel{~}{\psi }_{NR}(p)`$ is specified. Then with the form of $`\stackrel{~}{\varphi }_{L,T}`$ specified, the non-relativistic motivation for the light cone wave-function is abandoned and $`N_{L,T}`$ are constants adjusted to satisfy the normalization and decay conditions (Eqs. 54 and A) including the terms proportional to $`\left(k_t/m_c\right)^2`$. Terms proportional to $`\left(k_t/m_c\right)^2`$ mock up relativistic effects and constitute $`25\%`$ of the total elastic diffractive cross section. For the short distance model (GAUS-LC) $`\stackrel{~}{\varphi }_{L,T}`$ is specified directly. For the numerical results quoted below we take $`M_V=3.096`$GeV and $`m_c=1.3`$GeV.

GAUS-RF. In this model we take

$$\stackrel{~}{\psi }_{NR}(p)=N_{L,T}\mathrm{exp}\left(\frac{p^2}{\omega _{L,T}^2}\right),$$
(65)
where $`\omega _{L,T}`$ and $`N_{L,T}`$ are adjusted to satisfy Eq. A for the decay width and Eq. 54 for the norm. These parameters are $`N_L=1.74`$ GeV-3/2, $`\omega _L=0.66`$ GeV and $`N_T=1.86`$ GeV-3/2 and $`\omega _T=0.61`$ GeV.

CORNELL. In this model we take the form of $`\stackrel{~}{\psi }_{NR}(p)`$ from the wave-function of a non-relativistic potential model of $`J/\psi `$. Using a shooting method NR  we solve the radial Schrödinger equation for the Cornell potential to find $`\psi _{NR}(\stackrel{}{x})`$. The parameters of the potential are taken from Ref. Eichten . Note, the constituent quark mass of the potential model ($`m_c=1.8\text{GeV}`$) differs from the current quark mass ($`m_c=1.3\text{GeV}`$) in Eqs. 50, A, 63 for the the light cone wave-function. In this model the decay condition Eq. A is not satisfied.

GAUS-LC. We now turn to the short distance model GAUS-LC. At short distances the light cone vector meson wave-function is proportional to $`z(1z)`$ ref.BL ; FS . This motivates the following ansatz for the wave-function

$$\stackrel{~}{\varphi }_{L,T}(r,z)=N_{L,T}z(1z)\mathrm{exp}(\frac{k^2}{\omega _{L,T}^2})$$
(66)
$`N_{L,T}`$ and $`\omega _{L,T}`$ are adjusted to reproduce the normalization condition and the decay width. For $`J/\psi `$ we find: $`N_L=1.21`$ GeV-1, $`\omega _L=0.80`$ GeV and $`N_T=1.63`$ GeV-1 and $`\omega _T=0.58`$ GeV.

With the wave-function $`\stackrel{~}{\varphi }_{L,T}(k,z)`$ specified the wave-function $`\varphi _{L,T}(r,z)`$ and its derivative $`_r\varphi (r,z)`$. are found taking Fourier transforms. The final overlap between the photon and vector meson wave-function is given by Eq. A.


format: text

text:
# Strangeness Saturation: Dependence on System-Size, Centrality and Energya

## 1 Introduction

Statistical-thermal models (see  for a review) aim at reproducing hadron multiplicities in relativistic heavy-ion collisions in certain phase space regions by a small number of parameters. As these parameters evolve smoothly with external conditions, thermal models have predictive power. Moreover, the thermal parameters, taken literally, allow for a quantification of the conditions reached in the transient stage of compressed and heated strongly interacting matter. While there is exhaustive literature on the successful use of the thermal models, the foundation of the applicability is far from settled. One should also keep in mind that thermal models address only a subset of the wealth of observables.

The present contribution summarizes our recent work  investigating the dependence of strangeness saturation on system size, centrality and energy, both in full phase space ($`4\pi `$) and at mid-rapidity ($`y0`$).

## 2 Thermal Model

Our model is based on the grand canonical ensemble expression for the abundance of hadron species $`i`$ (Fermions \[+1\] or Bosons \[-1\]) with degeneracy $`g_i`$

$$N_i^{\mathrm{prim}}=Vg_i\frac{d^3p}{(2\pi )^3}𝑑m_i\left[\gamma _s^{\left|S_i\right|}\mathrm{exp}\left\{\frac{E_i\stackrel{}{\mu }_i\stackrel{}{Q}_i}{T}\right\}\pm 1\right]^1\text{BW}(m_i,\mathrm{\Gamma }_i),$$
(1)
where $`T`$ is the temperature, $`\stackrel{}{\mu }`$ a set of chemical potentials corresponding to the conserved charges $`\stackrel{}{Q}`$; $`E_i=\sqrt{\stackrel{}{p}^2+m_i^2}`$ with vacuum mass $`m_i`$.

For small particle numbers one has to turn to a canonical ensemble, e.g. by the projection method . It is important to include resonances (we include meson (baryon) states composed of the three lightest quarks and anti-quarks up to 2.3 (2.6) GeV); their widths are taken into account via the Breit-Wigner parameterization $`BW(m_i,\mathrm{\Gamma }_i)`$, which collapses to a $`\delta `$-function for the respective ground states. Due to feeding and resonance decays, the final hadron multiplicities become $`N_i=N_i^{\mathrm{prim}}+_j\{N_j^{\mathrm{prim}}Br(ji)N_i^{\mathrm{prim}}Br(ij)]\}`$, where $`Br`$ denotes the corresponding branching ratio.

Our focus here is the phenomenological strangeness saturation factor $`\gamma _s`$ in Eq. (1). It is thought  to parameterize possible deviations from chemical equilibrium in the strange sector. Whether one common factor $`\gamma _s`$ is sufficient, with powers $`|S_i|`$ determined by the total strangeness content of hadron $`i`$, needs still to be tested by analyzing data. Also in the non-strange sector such an off-equilibrium parameter should be introduced . Due to the small data samples at our disposal, such an additional parameter cannot be unambiguously fixed; therefore, we are forced here to assume full equilibration in the non-strange sector.

The parameters $`T`$ and $`\stackrel{}{\mu }`$ appear exponentially in Eq. (1), while the sensitivity on variations of $`\gamma _s`$ is weaker. For this reason, we study the mentioned parameter dependencies based on analyses of data sets comprised of the same hadron species. Otherwise variations of $`\gamma _s`$ can be absorbed in small changes of $`T`$ and $`\stackrel{}{\mu }`$. Unfortunately, this shrinks the wealth of available data to quite poor subsets.

Within the Cooper-Frye formalism, dynamical effects factor out of fully-integrated hadron yields , provided that freeze-out occurs on a hyper-surface with constant $`T`$ and $`\stackrel{}{\mu }`$. In this respect, Eq. (1) covers also dynamical situations. Instead of considering hadron ratios, we consider multiplicities $`N_i^{\mathrm{prim}}`$ with a fiducial normalization volume $`V`$ in Eq. (1).

## 3 Data Analysis

### 3.1 CERN-SPS energies

#### 3.1.1 Centrality and System Size Dependence at $`E_{\mathrm{beam}}=158`$ AGeV

$`4\pi `$ multiplicities of $`\pi ^\pm `$, $`K^\pm `$, $`\varphi `$ and $`N_{\mathrm{part}}`$ are at our disposal for (i) central collisions of C + C and Si + Si, and (ii) centrality binned collisions Pb + Pb. More baryon information is desirable, but only $`\overline{p}`$’s are available additionaly for Pb + Pb. The left panel of Fig. 1 summarizes our findings . We observe,
$``$ Pb + Pb: $`\gamma _s`$ increases with centrality, but stays below unity, irrespective of whether $`\overline{p}`$’s are included or not; inclusion of $`\overline{p}`$’s reduces $`\gamma _s`$ somewhat; minimizing $`\chi ^2`$ or the quadratic deviation of data and the model multiplicities, $`q^2`$, yields slightly different results.
$``$ C + C, Si + Si: the extracted values of $`\gamma _s`$ off-set the down-extrapolation for Pb + Pb. This again emphasizes the findings of NA49 : system-size and centrality dependencies differ.
$``$ Most striking is that $`\gamma _s`$ goes as the number of participants which underwent multiple collisions, as extracted from a Glauber model calculation (cf. Fig. 3 in ).

There is some controversy on the use of $`4\pi `$ or mid-rapidity data. Our analysis  of $`\pi ^\pm `$, $`K^\pm `$ and $`p^\pm `$ centrality-binned mid-rapidity data from Pb + Pb at 158 AGeV reveals that $`\gamma _s`$ does indeed differ for $`4\pi `$ and $`y0`$ data (see Fig. 1, middle panel). This is due to the fact that the rapidity distributions of various hadrons differ. Combining the results, even the very restricted data samples point to $`\gamma _s=1\pm 0.25`$ for most central collisions. This suggest that $`\gamma _s=1`$ with slightly renormalized values of $`T`$ and $`\stackrel{}{\mu }`$ is a good choice for the given beam energy and high centrality.

#### 3.1.2 Energy Dependence

Unfortunately, the previous analysis cannot be extended to extract energy dependencies. Rather, the $`4\pi `$ multiplicities of $`\pi ^\pm `$, $`K^\pm `$, $`\mathrm{\Lambda }`$ and $`\overline{\mathrm{\Lambda }}`$ at beam energies of 40, 80 and 158 AGeV are at our disposal as comparable data samples. The results of our analysis  are shown in Fig. 2: $`\gamma _s`$ displays a dip at 80 AGeV, while $`\gamma _s(40\mathrm{A}\mathrm{G}\mathrm{e}\mathrm{V})\gamma _s(158\mathrm{A}\mathrm{G}\mathrm{e}\mathrm{V})`$. (The significance of the drop at 80 AGeV still needs to be clarified.) $`\gamma _s(158\mathrm{A}\mathrm{G}\mathrm{e}\mathrm{V})`$ is compatible with the strangeness saturation factor reported in the previous subsection emerging from a data set suitable for a centrality analysis. The open circle follows from an analysis of mid-rapidity data of $`\pi ^\pm `$, $`K^\pm `$, $`\mathrm{\Lambda }`$ and $`\overline{\mathrm{\Lambda }}`$. The change relative to results reported in the previous subsection highlights the sensitivity on different data sets, if the number of analyzed hadron species is small.

### 3.2 RHIC Energies

Our analysis is restricted to the published data of $`\pi ^\pm `$, $`K^\pm `$ and $`p^\pm `$ at $`y0`$ in centrality-binned collisions of Au + Au at $`\sqrt{s_{NN}}=130`$ GeV . The upcoming data for further hadron species, including excited states, and data at $`\sqrt{s_{NN}}=200`$ GeV will complete our picture on the systematics of hadron multiplicities.

#### 3.2.1 Centrality Dependence

As for CERN-SPS energies, $`\gamma _s`$ increases with centrality, see right panel in Fig. 1. The overall trend follows again the number of multiple participant collisions obtained in a Glauber model calculation . The most central collisions point to strangeness saturation, in particular when feed-down from weak decays is included.

#### 3.2.2 Energy Dependence

Fig. 2 suggests that increasing energy causes increased strangeness saturation. However, this conclusion should be taken with caution, as Fig. 2 displays in one plot analyses of $`4\pi `$ and $`y0`$ data and data sets comprising different hadron species. Assuming that at CERN-SPS the $`4\pi `$ data are sensible, while at RHIC only $`y0`$ data should be analyzed within thermal models, then from Fig. 2 one does indeed derive the mentioned expectation: with increasing energy, $`\gamma _s`$ evolves towards saturation.

One observes in Fig. 2 that the thermal model delivers for $`pp`$ and $`\overline{p}p`$ collisions  values of $`\gamma _s`$ compatible with peripheral collisions.

## 4 Agreement with Data

Up to now we have reported the system-size, centrality and energy dependence of $`\gamma _s`$. The values of $`T`$ coincide approximately with the chiral (deconfinement) transition temperature of $`T165`$ MeV, smoothly increasing from 150 MeV at beam energy of 40 AGeV. $`\mu _B`$ (the baryon potential dominating $`\stackrel{}{\mu }`$) decreases from 350 to 220 MeV at CERN-SPS and drops to about 30 MeV at RHIC. The actual values depend on the data samples analyzed. The agreement of data with multiplicities from the thermal model is extremely good with the following exceptions:
$``$ $`\varphi /K^+`$ is underpredicted for peripheral collisions at beam energy of 158 AGeV, even if the $`\varphi `$ yields are included in the fit;
$``$ $`\mathrm{\Lambda }+\mathrm{\Xi }^0`$ yields follow the data trend but are in absolute normalization below data when using the thermal parameters extracted from $`\pi ^\pm `$, $`K^\pm `$ and $`p^\pm `$ yields at $`\sqrt{s_{NN}}=130`$ GeV; the $`\mathrm{\Xi }^{}`$ yields are at the lower end of the error bars.

To understand which trends in the data drive the change of $`\gamma _s`$ we exhibit in Fig. 3 scaled multiplicities as a function of $`N_{\mathrm{part}}`$. A scaling of $`\overline{p}`$ by $`N_{\mathrm{part}}^{1.15}`$ results in a flat distribution at beam energy of 158 AGeV. The scaled $`\pi ^\pm `$ is slightly decreasing, while $`K^\pm `$ increases. It is the relative decrease of the ”gap” between pion and kaon yields which drives $`\gamma _s`$ towards unity. Similarly, at $`\sqrt{s_{NN}}=130`$ GeV, the $`p^\pm `$ distributions are flat, as also $`K^\pm `$, but the $`\pi ^\pm `$ decrease with increasing $`N_{\mathrm{part}}`$. Again, the decreasing ”gap” between $`\pi ^\pm `$ and $`K^\pm `$ drives the increase of $`\gamma _s`$. Note that the effects are small and hardly visible in a log plot with unscaled yields.

As an aside, we mention that we adjusted the parameters of the color glass model  and obtain a good description of $`\pi ^\pm `$ and $`p^\pm `$ yields, while the model deviates in details from the $`K^\pm `$ data, see right panel of Fig. 3.

## 5 Summary

In summary we report here the dependence of the strangeness saturation factor, $`\gamma _s`$, on system size, centrality and energy. Most striking is the increase of $`\gamma _s`$ with centrality and the number of multiple participant collisions. Details of the energy dependence of $`\gamma _s`$ are puzzling and hampered by restricted data sets. Whenever possible, one should compare the findings with results obtained from the richer data samples for central collisions. Clearly, the present investigations should be extended to include multi-strange baryons and resonances in the analysis of more complete data sets. This will further test the applicability and predictive power of thermal models of statistical hadronization in describing hadron multiplicities.

## Note(s)

* The work is supported by BMBF 06DR921 and 06DR121.
* Visiting Fulbright Professor on leave of absence from the Brookhaven National Laboratory, Upton, NY, USA


text:
# Structural order parameter and itinerant electron magnetism in Cd2Re2O7

## Abstract

We describe the low temperature behaviour of the magnetic susceptibility of Cd2Re2O7 in terms of a Landau theory of structural phase transitions. We calculate the zero temperature Pauli susceptibility using a tight-binding approach to reveal the mechanism of coupling between the structural order parameter and itinerant magnetism.

,

Cd2Re2O7 is the only pyrochlore superconductor . Upon cooling, Cd2Re2O7 undergoes two structural phase transitions (SPT’s) well above the very low temperature superconducting phase transition . The first occurs at 200 K and is of second-order, while the other is first-order at 120 K. Neither SPT involves multiplication of the bcc primitive cell. We will refer to the phases as I, II and III, with corresponding space groups Fd$`\overline{3}`$m$``$I$`\overline{4}`$m$`2`$I$`4_122`$ respectively. The SPT’s are accompanied by anomalous behaviour of physical properties, including the electrical resistivity, magnetic susceptibility, Hall coefficient and thermoelectric power. In this paper we concentrate on the magnetic properties of Cd2Re2O7.

Cd2Re2O7 is a paramagnetic metal in the normal state. Above 400 K, the magnetic susceptibility displays approximate Curie-Weiss behaviour with a large, negative Weiss temperature . On cooling it reaches a broad maximum, and then falls off below 200 K . No anomaly has been found at the second SPT at 120 K. The kind of anomaly observed at 200 K is often associated with an antiferromagnetic phase transition, but Re nuclear quadrupole resonance experiments find no evidence of magnetic ordering in the low temperature phases , thus revealing the itinerant electronic nature of its paramagnetism.

In this article we show theoretically that the magnetic anomaly at 200 K is due to the electrostatic interaction between itinerant electrons and tetragonal distortions of the pyrochlore lattice. Recently, we found a single two-dimensional order parameter (OP) that describes both SPT’s in Cd2Re2O7 which corresponds to an instability of the lattice with respect to a long wavelength phonon mode of $`E_u`$ symmetry .

We introduce the two components of the structural OP $`(\eta _1,\eta _2)`$, which are linear combinations of the displacements of four Re atoms, corresponding to one bcc site, from their ideal positions in phase I, $`(x_m,y_m,z_m),m=1,2,3,4`$ .

$$\eta _1=(XY)/\sqrt{2},\eta _2=(X+Y2Z)/\sqrt{6},$$
(1)
where $`X=(x_1+x_2x_3x_4)/2`$, $`Y=(y_1y_2+y_3y_4)/2`$ and $`Z=(z_1z_2z_3+z_4)/2`$. Phase II is characterized by $`\eta _1=0`$, $`\eta _20`$ and phase III by $`\eta _10`$, $`\eta _2=0`$. The OP spans the $`E_u`$ representation of the cubic point group $`O_h`$. The magnetisation $`𝐌`$ corresponds to the irreducible representation $`F_{1g}`$. Considering the symmetric products $`[E_uE_u]=A_{1g}E_g`$ and $`[F_{1g}F_{1g}]=A_{1g}E_gF_{2g}`$, we find that there are two independent constants describing biquadratic coupling between the OP and $`𝐌`$. In an applied magnetic field $`𝐇`$, the magnetic part of the free energy is

$`F_M`$ $`=`$ $`A𝐌^2\mathrm{𝐌𝐇}+\gamma (\eta _1^2+\eta _2^2)𝐌^2`$
$`+\delta [(\eta _1^2\eta _2^2)(2M_z^2M_x^2M_y^2)`$
$`+2\sqrt{3}\eta _1\eta _2(M_x^2M_y^2)].`$
In all three phases the product $`\eta _1\eta _2`$ is zero, and the inverse magnetic susceptibility is calculated to be

$`\chi _{zz}^1=\chi _0^1+2\gamma (\eta _1^2+\eta _2^2)+4\delta (\eta _1^2\eta _2^2)`$
$`\chi _{xx}^1=\chi _0^1+2\gamma (\eta _1^2+\eta _2^2)2\delta (\eta _1^2\eta _2^2)`$ (3)
where $`\chi _0=(2A)^1`$ is the susceptibility of the high symmetry phase. In order to account for the experimentally observed behaviour at the two SPT’s, we should assume $`\delta \gamma `$ and $`\gamma >0`$ . Note the similarity between the present theory and the antiferromagnetic case, where the anomaly in the susceptibility also arises from biquadratic coupling between $`𝐌`$ and the antiferromagnetic vector.

The phenomenological model (Structural order parameter and itinerant electron magnetism in Cd2Re2O7) is based on group-theoretical arguments only and does not distinguish between localised and itinerant electronic magnetism. We now restrict our attention to the itinerant case. The tight-binding approach describes the dispersion of electronic bands $`ϵ_𝐤`$ resulting from the electrostatic potential of the lattice. Overlap integrals, which constitute matrix elements of the tight-binding Hamiltonian, depend on the radius vector between neighbouring atoms . Due to the relations (1), this dependence can be expressed in terms of the OP, $`ϵ_𝐤(\eta _1,\eta _2)`$.

For simplicity we consider a simple one-band model for small $`𝐤`$ around the center of the Brillouin zone. Group-theoretical considerations yield

$`ϵ_𝐤(\eta _1,\eta _2)`$ $`=`$ $`\mathrm{}^2/2m^{}[𝐤^2+\alpha (\eta _1^2+\eta _2^2)𝐤^2`$ (4)
$`+\beta [(\eta _1^2\eta _2^2)(2k_z^2k_x^2k_y^2)`$
$`+2\sqrt{3}\eta _1\eta _2(k_x^2k_y^2)].`$
Here $`m^{}`$ is the effective mass in phase I, and $`\alpha `$ and $`\beta `$ are constants which may be calculated by diagonalising the tight-binding Hamiltonian. Thus the effective mass components depend on the OP. In all three phases $`\eta _1\eta _2=0`$, therefore $`m_x=m_ym_z`$.

The density of states at the Fermi level (including spin) is

$`D(E_F)`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{}^2}}\left({\displaystyle \frac{3m_x^2m_zN}{\pi ^4V}}\right)^{1/3}`$ (5)
$``$ $`{\displaystyle \frac{m^{}}{\mathrm{}^2}}\left({\displaystyle \frac{3N}{\pi ^4V}}\right)^{1/3}[1\alpha (\eta _1^2+\eta _2^2)],`$
where $`N/V`$ is the number of conduction electrons per unit volume. Then the Pauli susceptibility at zero temperature is

$$\chi _P=\frac{g^2\mu _B^2m^{}}{4\mathrm{}^2}\left(\frac{3N}{\pi ^4V}\right)^{1/3}[1\alpha (\eta _1^2+\eta _2^2)].$$
(6)
We note that up to the second order in the OP, only the isotropic part of the effective mass contributes to the density of states and Pauli susceptibility, but higher order corrections will yield an anisotropic contribution from the OP. Moreover, spin-orbit coupling, according to band structure calculations , has a significant effect on the electronic spectrum, and also gives rise to anisotropy of the $`g`$ factor . Nevertheless, the simple model proposed here provides a mechanism of coupling between itinerant electrons and structural distortions, and the lowest order result (6) supports our assumption about the relative sizes of the phenomenological constants in (Structural order parameter and itinerant electron magnetism in Cd2Re2O7).

This work was supported by NSERC of Canada.


text:
# Study of pure annihilation type decays 𝐵→𝐷_𝑠^∗𝐾

## 1 Introduction

Perturbative QCD (PQCD) method for B decays has been developed for some years . It is successfully applied to exclusive $`B`$ meson decays recently, such as $`B\pi \pi `$ , $`BK\pi `$ , $`B\pi \rho `$  and other channels .

Very recently, pure annihilation type B decays are discussed in the PQCD approach, such as $`B^0D_s^{}K^+`$ and $`B^+D_s^+\overline{K}^0`$  decays and $`B^+D_s^+\varphi `$ decay . It is found that the annihilation type decay $`B^0D_s^{}K^+`$ has a sizable branching ratio of $`10^5`$, which has already been measured by experiments . In this paper, we will continue to compute the branching ratios of similar decays $`B^0D_s^{}K^+`$ and $`B^+D_s^+\overline{K}^0`$.

Because the four valence quarks in the final states $`D_s^{}`$ and $`K(\overline{K})`$ are different from the ones in the $`B`$ meson, the rare decays $`BD_s^{}K`$ are pure annihilation type decays. In the usual factorization approach (FA) , this decay picture is described as $`B`$ meson annihilating into vacuum and $`D_s^{}`$ and $`K`$ meson produced from vacuum afterwards. To calculate this decay in the FA, one needs the $`D_s^{}K`$ form factor at very large timelike momentum transfer $`𝒪(M_B)`$. However the form factor at such a large momentum transfer is not known in FA. The annihilation amplitude is a phenomenological parameter in QCD factorization approach (QCDF) , and the QCDF calculation of these decays is also unreliable.

In this paper, we will calculate these decays in PQCD approach. Similar to the $`BD_SK`$ decays, the $`W`$ boson exchange induce the four quark operator $`\overline{b}d\overline{c}u`$ or $`\overline{b}u\overline{d}c`$, and the $`\overline{s}s`$ quarks included in $`D_s^{}K`$ are produced from a gluon. This gluon attaches to any one of the quarks participating in the four quark operator. In the rest frame of $`B`$ meson, the produced $`s`$ or $`\overline{s}`$ quark included in $`D_s^{}K`$ final states has $`𝒪(M_B/4)`$ momentum, and the gluon producing them has $`q^2𝒪(M_B^2/4)`$. This is a hard gluon, so we can perturbatively treat the process where the four-quark operator exchanges a hard gluon with $`s\overline{s}`$ quark pair. It is a perturbative six quark interaction now.

In PQCD, the decay amplitude is separated into soft ($`\mathrm{\Phi }`$), hard ($`H`$), and harder ($`C`$) dynamics characterized by different scales. It is conceptually written as the convolution,

$$\text{Amplitude}d^4k_1d^4k_2d^4k_3\mathrm{Tr}\left[C(t)\mathrm{\Phi }_B(k_1)\mathrm{\Phi }_{D_s^{}}(k_2)\mathrm{\Phi }_K(k_3)H(k_1,k_2,k_3,t)\right],$$
(1)
where $`k_i`$’s are momenta of light quarks included in each meson. $`C(t)`$ is Wilson coefficient which results from the radiative corrections at short distance. $`\mathrm{\Phi }_M`$ is the wave function which describes hadronization of the quark and anti-quark into the meson $`M`$. $`H`$ describes the four quark operator and the quark pair from the sea connected by a hard gluon whose scale is at the order of $`M_B`$, so the hard part $`H`$ can be perturbatively calculated.

Some analytic formulas for the decay amplitudes of $`BD_s^{}K`$ decays will be given in the next section. In section 3, we give the numerical results and discussion. Finally, we conclude this study in section 4.

## 2 Analytic formulas

We consider the $`B`$ meson at rest for simplicity. In the light-cone coordinate, the $`B`$ meson momentum $`P_1`$, the $`D_s^{}`$ meson momentum $`P_2`$ and $`K`$ meson momentum $`P_3`$ are taken to be:

$$P_1=\frac{M_B}{\sqrt{2}}(1,1,\mathrm{𝟎}_T),P_2=\frac{M_B}{\sqrt{2}}(1,r^2,\mathrm{𝟎}_T),P_3=\frac{M_B}{\sqrt{2}}(0,1r^2,\mathrm{𝟎}_T),$$
(2)
where $`r=M_{D_s^{}}/M_B`$ and we neglect the $`K`$ meson’s mass $`M_K`$. The $`D_s^{}`$ meson’s longitudinal polarization vector $`ϵ`$ is given by $`ϵ=\frac{M_B}{\sqrt{2}M_{D_s^{}}}(1,r^2,\mathrm{𝟎}_T)`$. Denoting the light (anti-)quark momenta in $`B`$, $`D_s^{}`$ and $`K`$ mesons as $`k_1`$, $`k_2`$, and $`k_3`$, respectively, we can choose $`k_1=(x_1P_1^+,0,𝐤_{1T})`$, $`k_2=(x_2P_2^+,0,𝐤_{2T})`$, $`k_3=(0,x_3P_3^{},𝐤_{3T})`$. The decay amplitude in eq.(1) leads to:

$$\begin{array}{c}\text{Amplitude}𝑑x_1𝑑x_2𝑑x_3b_1𝑑b_1b_2𝑑b_2b_3𝑑b_3\hfill \\ \hfill \mathrm{Tr}\left[C(t)\mathrm{\Phi }_B(x_1,b_1)\mathrm{\Phi }_{D_s^{}}(x_2,b_2,ϵ)\mathrm{\Phi }_K(x_3,b_3)H(x_i,b_i,ϵ,t)S_t(x_i)e^{S(t)}\right],\end{array}$$
(3)
where $`b_i`$ is the conjugate space coordinate of $`k_{iT}`$, and $`t`$ is the largest energy scale in $`H`$. The large logarithms $`\mathrm{ln}(m_W/t)`$ coming from QCD radiative corrections to four quark operators are included in the Wilson coefficients $`C(t)`$. The last term, $`e^{S(t)}`$, contains two kinds of logarithms. One of the large logarithms is due to the renormalization of ultra-violet divergence $`\mathrm{ln}tb`$, the other is double logarithm $`\mathrm{ln}^2b`$ from the overlap of collinear and soft gluon corrections. This Sudakov form factor suppresses the soft dynamics effectively . Thus it makes perturbative calculation of the hard part $`H`$ applicable at intermediate scale, i.e.$`M_B`$ scale.

The $`B`$ meson wave function for incoming state and the $`D_s^{}`$ and $`K`$ meson wave function for outgoing state with up to twist-3 terms are written as:

$$\mathrm{\Phi }_{B,\alpha \beta }(x,b)=\frac{i}{\sqrt{2N_c}}\left[(\overline{)}P_1\gamma _5)_{\alpha \beta }+M_B\gamma _{5\alpha \beta }\right]\varphi _B(x,b),$$
(4)
$$\mathrm{\Phi }_{D_s^{},\alpha \beta }(x,b)=\frac{i}{\sqrt{2N_c}}\left[\overline{)}ϵ(\overline{)}P_2+M_{D_s^{}})_{\alpha \beta }\right]\varphi _{D_s^{}}(x,b),$$
(5)
$`\mathrm{\Phi }_{K,\alpha \beta }(x_3,b_3)`$ $`={\displaystyle \frac{i}{\sqrt{2N_c}}}[\gamma _5\overline{)}P_3\varphi _K^A(x_3,b_3)+m_{0K}\gamma _5\varphi _K^P(x_3,b_3)`$
$`+m_{0K}\gamma _5(\mathit{}\mathit{}1)\varphi _K^T(x_3,b_3)]_{\alpha \beta },`$ (6)
where $`v=(0,1,\mathrm{𝟎}_T)P_3`$, $`n=(1,0,\mathrm{𝟎}_T)`$, and $`m_{0K}=M_K^2/(m_u+m_s)`$.

Since the wave functions are process independent, one can use the same forms constraint by other decay channels  to make predictions here. Now the hard part $`H`$ in eq.(3) is the only channel dependent part for us to calculate perturbatively.

### 2.1 $`B^0D_s^{}K^+`$ decay

In the decay $`B^0D_s^{}K^+`$, the effective Hamiltonian at the scale lower than $`M_W`$ is given as

$$H_{\mathrm{eff}}=\frac{G_F}{\sqrt{2}}V_{cb}^{}V_{ud}\left[C_1(\mu )O_1(\mu )+C_2(\mu )O_2(\mu )\right],$$
(7)
$$O_1=(\overline{b}d)_{VA}(\overline{u}c)_{VA},O_2=(\overline{b}c)_{VA}(\overline{u}d)_{VA},$$
(8)
where $`C_{1,2}(\mu )`$ are Wilson coefficients at renormalization scale $`\mu `$, and summation in $`\mathrm{SU}(3)_c`$ color’s index $`\alpha `$ and chiral projection, $`_\alpha \overline{q}_\alpha \gamma ^\nu (1\gamma _5)q_\alpha ^{}`$ are abbreviated to $`(\overline{q}q^{})_{VA}`$. The lowest order diagrams of $`B^0D_s^{}K^+`$ are drawn in Fig.1 according to this effective Hamiltonian. There are altogether four quark diagrams. Just as what we said above, the decay $`BD_s^{}K`$ only has annihilation diagrams.

By calculating the hard part $`H`$ at the first order of $`\alpha _s`$, we get the following analytic formulas. With the meson wave functions, the amplitude for the factorizable annihilation diagrams in Fig.1(a) and (b) results in $`F_a^{(i=2)}`$ as:

$$\begin{array}{c}F_a^{(i)}=16\pi C_FM_B^2_0^1𝑑x_2𝑑x_3_0^{\mathrm{}}b_2𝑑b_2b_3𝑑b_3\varphi _{D_s^{}}(x_2,b_2)\hfill \\ \hfill \times [\{(1r^2)(1(1r^2)x_3)\varphi _K^A(x_3,b_3)+r(1r^2)x_3r_K\varphi _K^P(x_3,b_3)\\ \hfill r(1+r^2)r_K\varphi _K^T(x_3,b_3)\}E_f^i(t_a^1)h_a(x_2,x_3,b_2,b_3)\{(1r^2)x_2\varphi _K^A(x_3,b_3)\\ \hfill +2r(1r^2x_2)r_K\varphi _K^P(x_3,b_3)\}E_f^i(t_a^2)h_a(1x_3,1x_2,b_3,b_2)].\end{array}$$
(9)
In the function, $`C_F=4/3`$ is the group factor of $`\mathrm{SU}(3)_c`$, and $`r_K=m_{0K}/M_B`$. The functions $`E_f^i`$, $`t_a^{1,2}`$, $`h_a`$ will be given in the appendix. The distribution amplitudes $`\varphi _M`$, are given in the next section. The amplitude for the nonfactorizable annihilation diagrams in Fig.1(c) and (d) results in

$$\begin{array}{c}M_a=\frac{1}{\sqrt{2N_c}}64\pi C_FM_B^2_0^1𝑑x_1𝑑x_2𝑑x_3_0^{\mathrm{}}b_1𝑑b_1b_2𝑑b_2\varphi _B(x_1,b_1)\varphi _{D_s^{}}(x_2,b_2)\hfill \\ \hfill \times [\{(1r^2)((1r^2)(1x_3)r^2x_2)\varphi _K^A(x_3,b_2)\\ \hfill r\left(x_2(1r^2)(1x_3)\right)r_K\varphi _K^P(x_3,b_2)\\ \hfill r(x_2+(1r^2)(1x_3))r_K\varphi _K^T(x_3,b_2)\}E_m(t_m^1)h_a^{(1)}(x_1,x_2,x_3,b_1,b_2)\\ \hfill \{(1r^2)((1r^2)x_2+r^2)\varphi _K^A(x_3,b_2)\\ \hfill +r(x_2(1r^2)(1x_3))r_K\varphi _K^P(x_3,b_2)\\ \hfill +r(x_22r^2(1r^2)(1+x_3))r_K\varphi _K^T(x_3,b_2)\}E_m(t_m^2)h_a^{(2)}(x_1,x_2,x_3,b_1,b_2)],\end{array}$$
(10)
where $`x_1`$ dependence in the numerators of the hard part are neglected by the assumption $`x_1x_2,x_3`$.

Comparing with $`B^0D_s^{}K^+`$ , we find that the leading twist contribution, which is proportional to $`\varphi _K^A`$, is almost the same. However the subleading twist contribution, which is proportional to $`rr_K`$, change significantly, especially for the $`\varphi _K^P(x_3,b_2)`$ terms in (9) and (10).

The total decay amplitude for $`B^0D_s^{}K^+`$ decay is given as $`A=f_BF_a^{(2)}+M_a`$. The decay width is then

$$\mathrm{\Gamma }(B^0D_s^{}K^+)=\frac{G_F^2M_B^3}{128\pi }(1r^2)|V_{cb}^{}V_{ud}A|^2.$$
(11)
The charge conjugate decay $`\overline{B}^0D_s^+K^{}`$ is the same value as $`B^0D_s^{}K^+`$ just because $`|V_{cb}^{}V_{ud}|`$ is the same as $`|V_{cb}V_{ud}^{}|`$. Since there is only one kind of CKM phase involved in the decay, there is no CP violation in the standard model for this decay channel.

### 2.2 $`B^+D_s^+\overline{K}^0`$ decay

The effective Hamiltonian related to $`B^+D_s^+\overline{K}^0`$ decay is given as

$$H_{\mathrm{eff}}=\frac{G_F}{\sqrt{2}}V_{ub}^{}V_{cd}\left[C_1(\mu )O_1(\mu )+C_2(\mu )O_2(\mu )\right],$$
(12)
$$O_1=(\overline{b}d)_{VA}(\overline{c}u)_{VA},O_2=(\overline{b}u)_{VA}(\overline{c}d)_{VA}.$$
(13)
The amplitude for the factorizable annihilation diagrams results in $`F_a^{(i=1)}`$. The amplitude for the nonfactorizable annihilation diagrams is given as

$$\begin{array}{c}M_a^{}=\frac{1}{\sqrt{2N_c}}64\pi C_FM_B^2_0^1𝑑x_1𝑑x_2𝑑x_3_0^{\mathrm{}}b_1𝑑b_1b_2𝑑b_2\varphi _B(x_1,b_1)\varphi _{D_s^{}}(x_2,b_2)\hfill \\ \hfill \times [\{(12r^2)x_2\varphi _K^A(x_3,b_2)+r(x_2(1r^2)(1x_3))r_K\varphi _K^P(x_3,b_2)\\ \hfill +r(x_2+(1r^2)(1x_3))r_K\varphi _K^T(x_3,b_2)\}E_m^{}(t_m^1)h_a^{(1)}(x_1,x_2,x_3,b_1,b_2)\\ \hfill \{(1r^2)(1x_2r^2(1r^2)x_3)\varphi _K^A(x_3,b_2)\\ \hfill r\left(x_2(1r^2)(1x_3)\right)r_K\varphi _K^P(x_3,b_2)\\ \hfill r(x_22r^2(1r^2)(1x_3))r_K\varphi _K^T(x_3,b_2)\}E_m^{}(t_m^2)h_a^{(2)}(x_1,x_2,x_3,b_1,b_2)].\end{array}$$
(14)
Comparing with the decay to two pseudo-scaler final states $`B^+D_s^+\overline{K}^0`$ , the change only occur in $`r^2`$ terms. Therefore, the non-factorizable contribution in $`B^+`$ decay does not change much.

The total decay amplitude $`A^{}`$ and decay width $`\mathrm{\Gamma }`$ for $`B^+D_s^+\overline{K}^0`$ decay are given as

$$A^{}=f_BF_a^{(1)}+M_a^{},$$
(15)
$$\mathrm{\Gamma }(B^+D_s^+\overline{K}^0)=\frac{G_F^2M_B^3}{128\pi }(1r^2)|V_{ub}^{}V_{cd}A^{}|^2.$$
(16)
The decay width for CP conjugated mode, $`B^{}D_s^{}K^0`$, is the same value as $`B^+D_s^+\overline{K}^0`$. Similar to the $`B^0`$ decay, there is also no CP violation in this decay within standard model.

## 3 Numerical evaluation

We use the wave function of $`B`$ meson written as

$$\varphi _B(x,b)=N_Bx^2(1x)^2\mathrm{exp}\left[\frac{M_B^2x^2}{2\omega _b^2}\frac{1}{2}(\omega _bb)^2\right],$$
(17)
where $`N_B`$ is a normalization factor. Because the mass difference between $`D_s`$ and $`D_s^{}`$ is not large, for simplicity, their wave functions are chosen to be the same :

$$\varphi _{D_s^{}}(x,b)=\frac{3}{\sqrt{2N_c}}f_{D_s^{}}x(1x)\{1+a_{D_s^{}}(12x)\}\mathrm{exp}\left[\frac{1}{2}(\omega _{D_s^{}}b)^2\right].$$
(18)
Since $`c`$ quark is much heavier than $`s`$ quark, this function is peaked at $`c`$ quark side, i.e. small $`x`$ region.

The $`K`$ meson wave functions are given as

$`\varphi _K^A(x)`$ $`=`$ $`{\displaystyle \frac{f_K}{2\sqrt{2N_c}}}6x(1x)\left\{1a_1^K3\xi +a_2^K{\displaystyle \frac{3}{2}}(1+5\xi ^2)\right\},`$ (19)
$`\varphi _K^P(x)`$ $`=`$ $`{\displaystyle \frac{f_K}{2\sqrt{2N_c}}}\left\{1+a_{p1}^K{\displaystyle \frac{1}{2}}(1+3\xi ^2)+a_{p2}^K{\displaystyle \frac{1}{8}}(330\xi ^2+35\xi ^4)\right\},`$ (20)
$`\varphi _K^T(x)`$ $`=`$ $`{\displaystyle \frac{f_K}{2\sqrt{2N_c}}}(12x)\left\{1+a_T^K3(3+5\xi ^2)\right\},`$ (21)
where $`\xi =2x1`$. The parameters of these distribution amplitudes calculated from QCD sum rule  are given as

$$a_1^K=0.17,a_2^K=0.2,a_{p1}^K=0.212,a_{p2}^K=0.148,a_T^K=0.0527,$$
(22)
for $`m_{0K}=1.6\text{ GeV}`$. In addition, we use the following input parameters :

$$f_B=190\text{ MeV},f_K=160\text{ MeV},f_{D_s^{}}=241\text{ MeV},$$
(23)
$$m_{0K}=1.6\text{ GeV},\omega _b=0.4\text{ GeV},a_{D_s^{}}=0.3,\omega _{D_s^{}}=0.3\text{ GeV}.$$
(24)
For branching ratio estimation, we use the CKM matrix elements and the lifetimes of B mesons as following ,

$$|V_{ud}|=0.9734,|V_{ub}|=3.6\times 10^3,$$
(25)
$$|V_{cb}|=41.2\times 10^3,|V_{cd}|=0.224,$$
(26)
$$\tau _{B^\pm }=1.67\times 10^{12}\text{ s,}\tau _{B^0}=1.54\times 10^{12}\text{ s}.$$
(27)
The predicted branching ratios are

$$\mathrm{Br}(B^0D_s^{}K^+)=2.7\times 10^5,$$
(28)
$$\mathrm{Br}(B^+D_s^+\overline{K}^0)=4.0\times 10^8.$$
(29)
The branching ratio of $`B^+D_s^+\overline{K}^0`$ is much smaller than that of $`B^0D_s^{}K^+`$, due to the suppression from CKM matrix element $`V_{ub}^{}V_{cd}<<|V_{cb}^{}V_{ud}|`$. C-H Chen has also given the branching ratio of $`B^0D_s^{}K^+`$ , our result agrees with his.

In these processes of $`BD_s^{}K`$, only the longitudinal polarization of $`D_s^{}`$ has contribution. If we ignore the difference between $`f_{D_s}`$ and $`f_{D_s^{}}`$ and the difference between $`M_{D_s}`$ and $`M_{D_s^{}}`$, the branch ratios are thought to be the same as the corresponding $`BD_sK`$ decays. But in fact the branching ratio of $`B^0D_s^{}K^+`$ is a little smaller than that of $`B^0D_s^{}K^+`$, because the contribution of twist-3 wave function in $`B^0D_s^{}K^+`$ decay become negative. Although the larger $`f_{D_s^{}}`$ makes the branching ratio larger, a larger $`M_{D_s^{}}`$ leads to a smaller branching ratio. The effect of $`M_{D_s^{}}`$ is more dominant than that of $`f_{D_s^{}}`$. In decays of charged $`B^+`$, $`f_{D_s^{}}`$ makes branching ratios larger, and $`M_{D_s^{}}`$ has little effect, such that the $`B^+D_s^+\overline{K}^0`$ decay in this case has a little larger branching ratio than $`B^+D_s^+\overline{K}^0`$.

For the experimental side, there are only upper limits given at $`90`$% confidence level ,

$`\mathrm{Br}(B^0D_s^{}K^+)`$ $`<1.7\times 10^4,`$ (30)
$`\mathrm{Br}(B^+D_s^+\overline{K}^0)`$ $`<1.1\times 10^3.`$ (31)
Obviously, our results are consistent with the data.

In addition to the perturbative annihilation contributions, there is also a hadronic picture for the $`B^0D_s^{}K^+`$. The $`B`$ meson decays into $`D^{}`$ and $`\rho ^+`$, the secondary particles then exchanging a $`K^0`$, scatter into $`D_s^{}`$, $`K^+`$ through final state interaction afterwards. For charged decay, $`B^+`$ meson decays into $`D^0\rho ^+`$ then scatter into $`D_s^+`$ and $`\overline{K}^0`$ by exchanging a $`K^+`$. But this picture can not be calculated accurately. In , the results from PQCD approach for $`B^0D_s^{}K^+`$ decay were consistent with the experiment, which shows that the soft final state interaction may not be important.

The calculated branching ratios in PQCD are sensitive to various parameters, such as parameters in eqs.(22-24). It is necessary to give the sensitivity of the branching ratios when we choose the parameters to some extent. Table 1 shows the sensitivity of the branch ratios to $`30`$% change of $`a_{D_s^{}}`$, $`\omega _{D_s^{}}`$ and $`\omega _b`$. Here we don’t present the sensitivity to $`a_{p1,p2,T}^K`$ because the branching ratios are insensitive to them. It is found that the uncertainty of the predictions on PQCD is mainly due to $`\omega _{D_s^{}}`$, which characterizes the shape of $`D_s^{}`$ wave function.

From the above discussions, we can derive the uncertainties of the branching ratios within the suitable ranges on $`\omega _{D_s^{}}`$, $`a_{D_s^{}}`$ and $`\omega _b`$. The branching ratios normalized by the decay constants and the CKM matrix elements result in

$$\mathrm{Br}(B^0D_s^{}K^+)=(2.7\pm 0.6)\times 10^5\left(\frac{f_Bf_{D_s^{}}}{190\text{ MeV}240\text{ MeV}}\right)^2\left(\frac{|V_{cb}^{}V_{ud}|}{0.04120.9734}\right)^2,$$
(32)
$$\mathrm{Br}(B^+D_s^+\overline{K}^0)=(4.0\pm 0.8)\times 10^8\left(\frac{f_Bf_{D_s^{}}}{190\text{ MeV}240\text{ MeV}}\right)^2\left(\frac{|V_{ub}^{}V_{cd}|}{0.00360.224}\right)^2.$$
(33)
## 4 Conclusion

In this paper, we calculate the $`B^0D_s^{}K^+`$ and $`B^+D_s^+\overline{K}^0`$ decays in PQCD approach. These two decays occur purely via annihilation type diagrams because the four quarks in final states are not the same as the ones in $`B`$ meson. We argue that the soft final state interaction may be small in B decays. The PQCD calculation of annihilation decays is reliable. The branching ratio for $`B^0D_s^{}K^+`$ decay is sizable at order of $`10^5`$, which can be measured in the current $`B`$ factories Belle and BABAR. The small branching ratio of $`B^+D_s^+\overline{K}^0`$ ($`10^8`$) predicted in the SM, makes it sensitive to new physics contributions, which may be studied in the future LHC-B experiment.

## Acknowledgments

This work is partly supported by National Science Foundation of China under Grant (No. 90103013 and 10135060).

## Appendix A Some formulas used in the text

The function $`E_f^i`$, $`E_m`$, and $`E_m^{}`$ including Wilson coefficients are defined as

$$E_f^i(t)=a_i(t)\alpha _s(t)e^{S_D(t)S_K(t)},$$
(34)
$$E_m(t)=C_2(t)\alpha _s(t)e^{S_B(t)S_D(t)S_K(t)},$$
(35)
$$E_m^{}(t)=C_1(t)\alpha _s(t)e^{S_B(t)S_D(t)S_K(t)},$$
(36)
where

$$a_1(t)=C_1(t)/3+C_2(t),a_2(t)=C_1(t)+C_2(t)/3,$$
(37)
and $`S_B`$, $`S_D`$, and $`S_K`$ result from summing both double logarithms caused by soft gluon corrections and single ones due to the renormalization of ultra-violet divergence. Those factors are given in ref. . In the numerical analysis we use leading logarithms expressions for Wilson coefficients $`C_{1,2}`$ presented in ref..

The functions $`h_a`$, $`h_a^{(1)}`$, and $`h_a^{(2)}`$ in the decay amplitudes consist of two parts: one is the jet function $`S_t(x_i)`$ derived by the threshold resummation , the other is the propagator of virtual quark and gluon. They are defined by

$`h_a(x_2,x_3,b_2,b_3)=S_t(1x_3)\left({\displaystyle \frac{\pi i}{2}}\right)^2H_0^{(1)}(M_B\sqrt{(1r^2)x_2(1x_3)}b_2)`$
$`\times \{H_0^{(1)}(M_B\sqrt{(1r^2)(1x_3)}b_2)J_0(M_B\sqrt{(1r^2)(1x_3)}b_3)\theta (b_2b_3)+(b_2b_3)\},`$ (38)
$`h_a^{(j)}(x_1,x_2,x_3,b_1,b_2)=`$
$`\{{\displaystyle \frac{\pi i}{2}}\mathrm{H}_0^{(1)}(M_B\sqrt{(1r^2)x_2(1x_3)}b_1)\mathrm{J}_0(M_B\sqrt{(1r^2)x_2(1x_3)}b_2)\theta (b_1b_2)`$
$`+(b_1b_2)\}\times (\begin{array}{cc}\mathrm{K}_0(M_BF_{(j)}b_1),& \text{for}F_{(j)}^2>0\\ \frac{\pi i}{2}\mathrm{H}_0^{(1)}(M_B\sqrt{|F_{(j)}^2|}b_1),& \text{for}F_{(j)}^2<0\end{array}),`$ (39)
where $`\mathrm{H}_0^{(1)}(z)=\mathrm{J}_0(z)+i\mathrm{Y}_0(z)`$, and $`F_{(j)}`$s are defined by

$$F_{(1)}^2=(1r^2)(x_1x_2)(1x_3),F_{(2)}^2=x_1+x_2+(1r^2)(1x_1x_2)(1x_3).$$
(40)
We adopt the parametrization for $`S_t(x)`$ of the factorizable contributions,

$$S_t(x)=\frac{2^{1+2c}\mathrm{\Gamma }(3/2+c)}{\sqrt{\pi }\mathrm{\Gamma }(1+c)}[x(1x)]^c,c=0.3,$$
(41)
which is proposed in ref. . The hard scale $`t`$’s in the amplitudes are taken as the largest energy scale in the hard part $`H`$ to kill the large logarithmic radiative corrections:

$$t_a^1=\mathrm{max}(M_B\sqrt{(1r^2)(1x_3)},1/b_2,1/b_3),$$
(42)
$$t_a^2=\mathrm{max}(M_B\sqrt{(1r^2)x_2},1/b_2,1/b_3),$$
(43)
$$t_m^j=\mathrm{max}(M_B\sqrt{|F_{(j)}^2|},M_B\sqrt{(1r^2)x_2(1x_3)},1/b_1,1/b_2).$$
(44)