text:
# Antideuteron production in proton-proton and proton-nucleus collisions

## 1 Introduction

The interest in the study of production of light antinuclei, in particular, antideuterons in proton-proton and proton-nucleus collisions has recently intensified. There are, at least, two major reasons for this. Firstly, studies of antideuteron production in cosmic space can be a very powerful tool to search for antimatter in the Universe Salati . It is planed to measure the antideuteron flux in the future AMS AMS  and PAMELA PAMELA  experiments.

Secondly, the possibility to make experiments with antideuteron beams was discussed recently (see Iazzi99  and references therein).

The theoretical estimation of antideuteron production preformed in the article Salati  are based on the well-known coalescence model CoalModel  which supposes that two nucleons fuse into a deuteron if the momentum of their relative motion is smaller than a certain quantity $`p_0`$, the coalescence radius in momentum space. This momentum $`p_0`$ is considered as a free parameter to be fixed from the experimental data.

More than ten years ago, a quite simple diagrammatic approach to the coalescence model provided a microscopical basis for the coalescence model and expressed the parameter $`p_0`$ in terms of the slope parameter of the inclusive nucleon spectrum and the wave function of the produced nucleus Kolyb . Within this approach, it appears to be possible to explain the empirical fact of approximate equality of the values of the coalescence radii for the description of the yields of various light fragments under similar kinematical conditions.

The aims of this article are to generalize this diagrammatic approach to antinuclei production (by introduction of threshold effects and by taking account of anisotropy of angular distributions) and to apply this model to the antideuteron production in proton-proton and proton-nucleus collisions. Particular attention will be paid to proton-proton collisions due to their interest for astrophysics. We will show that, in the cases where the inclusive antiproton production cross-section and deuteron wave function are well-known, this approach can describe quite well the inclusive antideuteron production cross-section without any additional parameter.

Note, that there are few articles in which the antideuteron production is discussed within different approaches Mro ; Dover ; LeuHeinz ; SchHeinz . Some of them require additional parameters to describe the experimental data and none of them describes the whole ensemble of experimental data.

The article is organized as follows. In section 2 the main ideas of the diagrammatic approach to the coalescence model are described and this approach is generalized to the case of antideuteron production. In Section 3, the description of the experimental data are presented. Finally, we provide a brief summary of the results.

## 2 Diagrammatic approach to the coalescence model

Let us remind the reader the main ideas of the diagrammatic approach to the coalescence model Kolyb . As a basis for the coalescence model, the simplest Feynman diagram of Fig. 1, corresponding to fusion of two nucleons is considered. Here the symbol $`f`$ denotes the state of all other particles except the nucleons 1 and 2 which form the deuteron. The physical picture behind this diagram is quite simple: the nucleons produced in a collision (block A) are slightly virtual and can fuse without further interaction with the nuclear field. This simplest diagram is not the only possible contribution. However, as was shown in Braun , there are mutual cancellations of a number of other diagrams and, as a result, the diagram of Fig. 1 is the dominant one.

Let us remind briefly how to calculate this diagram by using the nonrelativistic technics developed in Shapiro . The amplitude $`M`$ of this processus can be written as

$`M`$ $`=`$ $`{\displaystyle \frac{d^3p_1dE_1}{(2\pi )^4}\frac{d^3p_2dE_2}{(2\pi )^4}M_A}`$
$`{\displaystyle \frac{2im_\text{p}}{\text{p}_1^22m_\text{p}E_1i0}}{\displaystyle \frac{2im_\text{p}}{\text{p}_2^22m_\text{p}E_2i0}}M_\text{d}`$
$`(2\pi )^4\delta ^3(\text{p}_1+\text{p}_2\text{P})\delta (E_1+E_2E\epsilon )`$
Here $`M_A`$ is the diagram corresponding to the block A (production of nucleons 1 and 2 and other particles in the final state $`f`$), $`M_\text{d}`$ the vertex of coalescence to the deuteron, $`m_\text{p}`$ the nucleon mass. Two fractions are propagators of nucleons 1 and 2, the integrals are done over energies and momenta of these virtual particles. The last delta functions reflect energy-momentum conservation in the deuteron vertex ($`\text{P}=\text{p}_1+\text{p}_2`$ is deuteron momentum, $`E=\text{P}^2/4m_\text{p}`$ its kinetic energy, $`\epsilon 2.2`$ MeV its binding energy).

After trivial integration over $`\text{p}_2`$ and $`E_2`$ and introduction of relative momentum $`\text{q}=\frac{1}{2}(\text{p}_1\text{p}_2)`$, one obtains the following expression:

$`M`$ $`=`$ $`{\displaystyle \frac{d^3qdE_1}{(2\pi )^4}M_A(\text{P},\text{q})M_\text{d}\frac{2im_\text{p}}{(\text{q}+\text{P}/2)^22m_\text{p}E_1i0}}`$
$`{\displaystyle \frac{2im_\text{p}}{(\text{q}\text{P}/2)^22m_\text{p}(\text{P}^2/4m_\text{p}E_1+\epsilon )i0}}`$
The final integration over $`dE_1`$ (one supposes as usually that the only singularities in the complex $`E_1`$ plane are those of propagators) gives

$`M=i{\displaystyle \frac{d^3q}{(2\pi )^3}M_A(\text{P},\text{q})\phi _d(\text{q})}`$ (1)
where the deuteron wave function

$`\phi _d(\text{q})={\displaystyle \frac{m_\text{p}M_\text{d}}{\text{q}^2+m_\text{p}|\epsilon |}}`$
is normalized by the condition

$`{\displaystyle \frac{d^3q}{(2\pi )^3}|\phi _d(\text{q})|^2}=1.`$
To perform further calculations one needs to make an assumption about the dependance of the amplitude $`M_A`$ corresponding to the block A on its variables (the particle momenta). It can be shown Kolyb , that the simplest hypothesis that this amplitude is constant gives rise to a wrong result: the production cross-section appears to be zero. One can see it from (1). When $`M_A`$ does not depend on momenta the transition amplitude $`M`$ becomes proportional to

$`M=iM_A{\displaystyle \frac{d^3q}{(2\pi )^3}\phi _d(\text{q})}\phi _d(\text{r}=0)=0,`$
i.e. the deuteron wave function at the origin which, for realistic potentials, is equal to zero.

Therefore, the momentum dependance of $`M_A`$ has to be introduced, for instance, in a “minimal” way: the inclusive nucleon spectra usually have a decreasing form and can be parameterized by a Gaussian function in rather wide parameter regions:

$`E_\text{p}{\displaystyle \frac{d^3\sigma _\text{p}}{dp_\text{p}^3}}\mathrm{exp}\left(\text{p}_\text{p}^2/Q^2\right),`$ (2)
where $`Q`$ is related to the slope parameter. Accordingly, the amplitude $`M_A`$ can be written in the following way:

$`M_A=C\mathrm{exp}\left({\displaystyle \frac{\text{p}_1^2+\text{p}_2^2}{2Q^2}}\right)=C\mathrm{exp}\left({\displaystyle \frac{\text{P}^2}{4Q^2}}{\displaystyle \frac{\text{q}^2}{Q^2}}\right),`$ (3)
where the center-of-mass motion of the two nucleons is separated from their relative motion. The amplitude $`M_A`$ determines the cross-section for simultaneous production of two nucleons which can be expressed in a standard way which supposes statistical independence in production of the two nucleons as a product of inclusive cross-sections:

$`{\displaystyle \frac{d^6\sigma _{\text{pn}}}{dp_\text{p}^3dp_\text{n}^3}}={\displaystyle \frac{1}{\sigma _{inel}}}{\displaystyle \frac{d^3\sigma _\text{p}}{dp_\text{p}^3}}{\displaystyle \frac{d^3\sigma _\text{n}}{dp_\text{n}^3}},`$ (4)
where $`\sigma _{inel}`$ is the cross-section of inelastic interaction of initial particles.

After the substitution of (3) into the expression for the diagram of Fig. 1 and taking into account (4), the cross-section for the formation of deuterons takes the form

$`E_\text{d}{\displaystyle \frac{d^3\sigma _\text{d}}{dp_\text{d}^3}}=12\pi ^3|S|^2{\displaystyle \frac{1}{m_\text{p}\sigma _{inel}}}E_\text{p}{\displaystyle \frac{d^3\sigma _\text{p}}{dp_\text{p}^3}}E_\text{n}{\displaystyle \frac{d^3\sigma _\text{n}}{dp_\text{n}^3}},`$ (5)
where $`\text{p}_\text{d}2\text{p}_\text{p}`$, and

$`S={\displaystyle \frac{d^3q}{(2\pi )^3}\phi _d(\text{q})\mathrm{exp}\left(\frac{\text{q}^2}{Q^2}\right)}.`$ (6)
The structure of (4) is exactly the same as that of the coalescence model and one can obtain easily the following expression of the coalescence radius in momentum space111In this definition of $`p_0`$, the spin of outgoing particles is not taken into account (see Kolyb  for discussion). Note also that the definition of the parameter $`p_0`$ given in Salati  differs from the standard one by a factor of 2.

$`p_0^3=36\pi ^2\left|{\displaystyle \frac{d^3q}{(2\pi )^3}\phi _d(\text{q})\mathrm{exp}\left(\frac{\text{q}^2}{Q^2}\right)}\right|^2.`$ (7)
Thus, the approach based on the diagram of Fig. 1 reproduces the coalescence model with $`p_0`$ which is no more a free parameter but it is determined by the inclusive proton spectrum and by the deuteron wave function. As an example, in Fig. 2 the values of $`p_0`$ as a function of $`Q`$ (equation (8)) are presented for different nucleon-nucleon potentials (Paris Paris , Bonn Bonn , and different versions of Nijmegen potential Nijmegen ).

One can see that all potentials give the same result up to $`Q300`$ MeV/$`c`$ where the deuteron wave function in momentum space is quite well known. For higher values of $`Q`$ the difference between the predictions of different nucleon-nucleon potentials can be very important (taking into account the fact that the cross-section is proportional to the third power of the coalescence radius in momentum space $`p_0`$). In this work, the Paris potential was chosen for further calculations.

The isotropic angular dependence supposed in (2) is quite frequently used in nonrelativistic collisions. In the relativistic case, the dependencies on transversal and longitudinal momentum can be very different. However, the formulae obtained within this approach can be easily generalized to any angular dependence.

If the inclusive nucleon cross-section is parameterized by an amplitude $`M_{\text{p}_1}`$

$`E_1{\displaystyle \frac{d^3\sigma _1}{dp_1^3}}=|M_{\text{p}_1}|^2,`$ (8)
the cross-section for the deuteron formation can be written as (see (1))

$`E_\text{d}{\displaystyle \frac{d^3\sigma _\text{d}}{dp_\text{d}^3}}={\displaystyle \frac{12\pi ^3}{\sigma _{inel}m_\text{p}}}\left|{\displaystyle \frac{d^3q}{(2\pi )^3}M_{\text{p}_1}M_{\text{p}_2}\phi _d(\text{q})}\right|^2.`$ (9)
It is clear that this model can be practically directly used to describe the production of antideuterons. The only problem is a presence of the threshold in the antiparticle production cross-sections. The coalescence model and the approach used here (10) are not valid in the near threshold region. Therefore one needs to propose a phenomenological procedure to describe experimental data near the threshold. For proton-proton collisions, the authors of Salati  have proposed a quite simple prescription: the center of mass energy available for the production of the second antinucleon has to be reduced by twice the energy carried away by the first antinucleon $`E_{\text{}}`$. In other words, the two antinucleons are supposed to be produced at different energies: $`\sqrt{s}`$ and $`\sqrt{s}2E_{\text{}}`$.

In this article, the antideuteron production threshold is taken into account in a slightly different way. In proton-proton collisions, the main reaction giving antideuterons is $`\text{pp}\text{d̄pppn}`$. Near the threshold of this reaction, the energy dependence of the antideuteron production cross-section is mostly determined by the phase space of four nucleons $`\mathrm{\Phi }(\sqrt{s}E_{\text{}};m,m,m,m)`$,

$`E_{\text{}}{\displaystyle \frac{d^3\sigma _{\text{}}}{dp_{\text{}}^3}}\mathrm{\Phi }(\sqrt{s}E_{\text{}};m,m,m,m).`$ (10)
The phase space $`\mathrm{\Phi }`$ for $`n`$ particles with masses, momenta and energies, respectively, $`m_i`$, $`\text{p}_i`$, $`E_i`$ is defined in usual way (in cms)

$`\mathrm{\Phi }(\sqrt{s};m_1,m_2,\mathrm{}m_n)=`$
$`{\displaystyle \underset{i=1}{\overset{n}{}}}{\displaystyle \frac{1}{(2\pi )^3}}{\displaystyle \frac{d^3p_i}{2E_i}}\delta ^3\left({\displaystyle \underset{i=1}{\overset{n}{}}}\text{p}_i\right)\delta \left({\displaystyle \underset{i=1}{\overset{n}{}}}\text{E}_i\sqrt{s}\right)`$
and is calculated by using standard CERN library program CERNlib . $`\sqrt{s}`$ is the total energy available for these $`n`$ particles in the center of mass system.

Therefore one can introduce a phenomenological correction factor $`R`$ to the formula (10) defined as

$`R(\sqrt{s}E_{\text{}})={\displaystyle \frac{\mathrm{\Phi }(\sqrt{s}E_{\text{}};m,m,m,m)}{\mathrm{\Phi }(\sqrt{s}E_{\text{}};0,0,0,0)}},`$ (11)
where the denominator contains the ultrarelativistic phase space to ensure $`R`$ to be dimensionless and to have correct behavior at high energies ($`R1`$). The behavior of $`R(x)`$ is presented in Fig. 3.

There are, at least, two advantages with respect to the prescription chosen in Salati . Firstly, one makes no assumption about mechanism of the production. Secondly, this correction factor has correct kinematic behavior both near the threshold and at high energies. We’ll discuss the role of this factor later.

## 3 Description of experimental data

Before presenting of the results let us make some preliminary remarks about existing experimental data.

* The experimental data on antideuteron production are not abundant and are much less informative than that on deuteron production. There are only a few different experimental observations of antideuteron production in proton-proton ISRdbar73 ; ISRdbar78 ; ISRdbar75 ; ISRout ; IHEPpp , proton-nucleus IHEPpp ; IHEPpAl69 ; IHEPpAl71 ; FL , and nucleus-nucleus collisions Aoki ; AGSAuPb ; NA52 . Some experimental results cannot be analyzed within our approach because the experimental information is not complete. For instance, in Dorfan  there is no data on antiproton production for corresponding energies; in Bozzoli  only relative spectra (antiprotons to $`\pi ^{}`$ and antideuterons to $`\pi ^{}`$) were measured.
* To obtain a reasonable description within the present approach one needs to know the cross-section of antiproton production for the antiproton momentum equal to a half of the antideuteron one. Unfortunately, this condition is rarely satisfied: in most experiments, the differential cross-sections of the antiproton and antideuteron production being measured for approximately the same momentum. Therefore, to apply the method one has to extrapolate the antiproton data to another kinematical region. This procedure, of course, introduces an additional error.
* In principle, the inclusive cross-sections discussed here are functions of two kinematical variables (for instance, transversal and longitudinal momentum) and one has to present the results in three-dimensional form. However, in each experiment one has only a few experimental points and the results are presented as a function of one variable (either total or transversal momentum).

The total inelastic cross-section was taken from the PDG data PDG  (for proton-proton collisions) or described by the well-known parameterization Sigmainel  (for proton-nucleus and nucleus-nucleus collisions).

### 3.1 Proton-proton collisions

Let us begin the analysis with the most informative experiment performed on the ISR at CERN. In a few experiments, the spectra of antiprotons ISRpbar  and antideuterons ISRdbar73 ; ISRdbar78  were measured in pp-collisions at $`\sqrt{s}=53`$ GeV. The detector was situated at 90 degrees (in this geometry, the total momentum of outgoing particles coincides with the transversal one). Inclusive antiproton cross-sections over quite large regions of momentum and $`\sqrt{s}`$ was obtained in ISRpbar . An example of an experimental distribution of antiprotons for $`\sqrt{s}=53`$ GeV is presented in Fig. 4 in comparison with two different parameterizations of the data.

The first parameterization (solid line) proposed by the authors of the experiment ISRpbar :

$`E_{\text{}}{\displaystyle \frac{d^3\sigma _{\text{}}}{dp_{\text{}}^3}}=As^\alpha \mathrm{exp}[Bp_t],`$ (12)
with parameters $`A=0.195`$ mbarn/(GeV/$`c`$)2, $`\alpha =0.310`$, $`B=2.49`$ (GeV/$`c`$)-1 describes perfectly these data. The second one (dashed line)

$`E_{\text{}}{\displaystyle \frac{d^3\sigma _{\text{}}}{dp_{\text{}}^3}}=f\mathrm{exp}[(Ap_t+Bp_t^2)],`$ (13)
is the frequently used parameterization proposed by Tan and Ng TanNg  which works quite well in a wide region of $`p_t`$ and $`\sqrt{s}`$. Here $`f=f(E^{},\sqrt{s})`$, $`A=A(E^{},\sqrt{s})`$ and $`B=B(E^{},\sqrt{s})`$ are known functions of $`\sqrt{s}`$ and of the antiproton energy in center-of-mass system $`E^{}`$. This formula gives here reasonable values of the cross-section but the trend is not well reproduced. We present here both quite close parameterizations to demonstrate the difference in description of the data on antideuteron production. The corresponding cross-section measured in this experiment ISRdbar73 ; ISRdbar78  are given in Fig. 5 and compared with different calculations.

Comparison between the two parameterizations shows that the description of the antideuteron production is quite sensitive to the antiproton production cross-sections: 20-30$`\%`$ difference in description of the p̄ data can result in a factor of 2 for the d̄. This difference can be even more pronounced if one has to extrapolate a chosen parameterization. As we mentioned previously $`p_{\text{}}2p_{\text{}}`$. This condition is satisfied for these ISR data. In some other experiments presented hereafter, it is not the case.

In this figure, one can see also importance of the anisotropy of angular distributions. An exponential parameterization (13) of the antiproton production cross-section can be seen in two ways: as a function of the total antiproton momentum or of the transversal one (as stated previously, in the particular geometry of this experiment, they are equal to each other). However in the integral (10), all directions (and not only transverse one) are presented and the difference between total and transversal momenta can be quite important. Thus for the data under consideration the solid line represents the results with parameterization (13) whereas the dotted line corresponds to the same parameterization but with total momentum instead of transversal one (antiproton production cross-section supposed to be isotropic).

By using the picture of coalescence model, it is quite easy to see when anisotropy can be important in description of experimental data. If the total antideuteron momentum $`p_{cm}`$ is very high with respect to the coalescence radius in momentum space $`p_0`$, antiproton and antineutron are produced in approximately the same direction (the direction of antideuteron momentum) and an anisotropy of antinucleon angular distributions plays no role. If $`p_0p_{cm}`$, the two antinucleons can propagate in quite different directions before the coalescence and it is necessary to take anisotropy into account correctly.

For these ISR data, the total antideuteron momentum $`p_{cm}`$ is of the order of the coalescence momentum $`p_0`$ and the anisotropy effect is seen clearly. In all other experimental data discussed hereafter, $`p_{cm}`$ is very high with respect to $`p_0`$ and the anisotropy effect is not so important.

On can thus understand easily that if one uses anisotropic cross-sections the effect of the D-wave in the deuteron wave function can be quite important (for isotropic Gaussian parametrization (1), the D-wave contribution is explicitly equal to 0). For these ISR data, the introduction of D-wave contribution into the deuteron wave function divides the value of the cross-section by a factor of 2.

In all calculations presented in this article, the Paris wave function is used. For most experimental data analyzed in this article, the choice of the deuteron wave function is not crucial: if characteristic slope parameter of the inclusive antinucleon spectrum is less than approximately 0.5 GeV/$`c`$, all potential models give close values of $`p_0`$ (see Fig. 2). However, for some sets of the data it is not the case (we will mention it where necessary).

Another set of experimental data for the same $`\sqrt{s}`$ but with very large longitudinal component of antideuteron momentum $`p_l57`$ GeV/$`c`$ was measured in ISRdbar75  and is presented in Fig. 6. The antiproton spectrum is taken from an experiment performed by this group pbarforISR75 . Unfortunately, the parameterization of Tang and Ng does not work well here (its prediction exceeds systematically the data by a factor of 2). Therefore, we fitted the data by a Gaussian (proposed also by the authors of pbarforISR75 ) and by an exponential function of $`p_t`$. Both parameterizations give quite good antideuteron production cross section.

However, it appears to be impossible to describe reasonably another set of ISR data measured at very big transfer momentum ISRout . The parameterization of Tang and Ng cannot reproduce the antiproton production cross-section: the order of magnitude is correct but not a general trend which has quite unusual form (the cross-section increases with $`p_t`$ increasing! Thus the absence of the p̄ cross-section parameterization does not allow us to obtain a reasonable description of the d̄ production data (the use of Tan and Ng parameterization as input gives a result 4–10 times higher than the experimental data).

The last measurement of d̄ production in pp-collisions was performed ten years later at the IHEP machine IHEPpp  at lower cms energies ($`\sqrt{s}=11,5`$ GeV) and with different geometry (with fixed target). The data were taken at very high $`p_t`$ and $`p_l`$ and are not very rich (two points both for p̄ and d̄). Tang and Ng parameterization fails to describe the data (by a factor of ten for the highest momentum) and the characteristic slope parameter of the inclusive p̄ spectrum is quite high (where different potential models give quite different (by a factor of 2) predictions). However, we decided to present this instructive example because one can estimate here the role of the threshold effect by the procedure proposed in (12). One can fit the antiproton data (two points) both by a Gaussian and an exponential function of $`p_t`$. The corresponding predictions for the d̄ production cross-section are given in Fig. 7 by dotted and solid line respectively (without threshold effect). Once the threshold is taken into account by introduction of the factor, $`R`$ (12), the agreement with the experimental data is improved significantly (dashed line represents the calculations with an exponential parameterization and the threshold effect).

### 3.2 Proton-nucleus collisions

Unfortunately, there are no more exploitable data on the antideuteron production in proton-proton collisions. To test further the model, we analyzed available data on the d̄ production in proton-nucleus collisions. Here, one has no more general parameterization like TanNg  and, for each set of data, a different parameterization is used 222Note that in high energy proton-nucleus and nucleus-nucleus collisions, the physical center-of-mass system is that of nucleon (from the target) - nucleon (from the beam) one..

In the same IHEP experiment IHEPpp , the production of antideuterons was measured on Be and Pb targets also.

For these data, one can make the same remarks as for the data obtained in the proton-proton collisions (only a few experimental points, very high momenta of outgoing particles, and quite big value of the slope parameter). The results are presented in Fig. 8 and look quite encouraging.

There are also old IHEP measurements of the antideuteron production in p-Al collisions IHEPpAl69 ; IHEPpAl71 . The data were taken in the forward (or practically forward) direction and we have no information about $`p_t`$ dependence of the antiproton production cross-section (thus we can not take completely into account the anisotropy of angular distributions). Therefore, in Fig. 9, the results of calculations and the data are presented as a function of the total d̄ momentum. Note that, for the Gaussian function, the value of the slope parameter $`Q`$ appears to be very large (of the order of 1.3 GeV/$`c`$) where the deuteron wave function is not known and the results depend strongly on the nuclear potential.

This analysis can be completed by the measurement performed in FNAL FL  where the antideuteron production cross-section was measured for different targets (Be, Ti, W) at intermediate (with respect to the ISR and IHEP experiments) energies ($`\sqrt{s}=23.7`$ GeV/$`c`$). The data were taken at quite high transverse momentum. The theoretical results presented in Fig. 10 are in quite good agreement with the experimental data for the three targets. One of the reasons for this good agreement is a good knowledge of the antiproton production cross-section obtained in this experiment.

It is necessary to note that there are also some experimental results on antideuteron production in nucleus-nucleus collision obtained in the AGS experiment (one point in Si+Al collision Aoki  and two points in Au+Pb AGSAuPb ) and in the NA52 experiment NA52 . Unfortunately here, one has neither a good parameterization of the antiproton production cross-section nor a reliable parameterization for the total inelastic cross-section. The measurements were made in the forward direction and one has no information about the $`p_t`$ dependance of the antiproton production cross-sections (thus it is impossible to take into account correctly the anisotropy of angular distributions which can be very important here). Simple parameterizations (Gaussian and exponential) gave no satisfactory description. In general, the discrepancy between the calculations and the data is of order of a factor 5, which is not very surprising taking into account all these remarks. As an example, in Fig. 11 the calculation of the antideuteron production cross section in Pb+Pb collisions as a function of total momentum in the center-of-mass frame in comparison with the experimental data NA52  are presented. The data cover a very large momentum region and the characteristic momentum in the antiproton production cross-section is very big (of the order of 2 GeV/$`c`$).

## 4 Conclusions

Our main conclusion is that the diagrammatic approach to the coalescence model developed in Kolyb  can be successfully applied to the description of the antideuteron production in proton-proton and proton-nucleus collisions. There are two modifications: firstly, it is necessary to take into account the threshold effect and, secondly, one must include specific consideration of the strong anisotropy of angular distributions of antiproton production. Once these phenomena are taken into account, the model can describe most of existing experimental data on antideuteron production in proton-proton and proton-nucleus collisions.

The successful reproduction of experimental data suggests that good knowledge of the antinucleon production cross-section and of the deuteron wave function allow to describe the antideuteron production cross-section in quite large region of kinematic variables without any additional parameter.

## Acknowledgments

During the preparation of this article, we learnt the death of our friend and colleague V.M. Kolybasov who proposed the model discussed in this article. We would like to dedicate this work to his memory.

The authors would like to thank A. Barrau, M. Buénerd, A.J. Cole, L. Dérome, V.A. Karmanov and P. Salati for very useful and simulating discussions.


format: text

text:
# Scanning Tunneling Spectroscopy on Single Crystal MgB2

## 1 Introduction

Since the first report of superconductivity with a remarkably high $`T_\text{c}=39`$ K in magnesium diboride (MgB2) only two years ago , enormous progress has been made in the studies of this material. From a fundamental point of view, a central issue has been to establish that MgB2 is a two-gap superconductor - a concept which was introduced already in the fifties , and which has now found renewed relevance. Two-gap superconductivity in MgB2 was first predicted theoretically  and is now commonly accepted, confirmed by a large number of experiments .

The two gaps in MgB2 are associated with different parts of the Fermi surface, which is composed of 4 seperate sheets . Two of these are coaxial cylindrical sheets parallel to $`c^{}`$, derived from $`\sigma `$-antibonding states of the boron $`p_{xy}`$ orbitals ($`\sigma `$-band). The other two sheets are derived from $`\pi `$-bonding and antibonding states of the boron $`p_z`$ orbitals ($`\pi `$-band) and are three-dimensional. Calculations of the superconducting gap size for the two different bands yields $`\mathrm{\Delta }_\sigma 7`$ meV, and $`\mathrm{\Delta }_\pi =13`$ meV . These values are respectively above and below what is expected for a BCS $`s`$-wave superconductor, $`\mathrm{\Delta }=1.76k_\text{B}T_\text{c}=5.9`$ meV. Furthermore, and in contrast to many materials or alloys studied earlier, the two bands in MgB2 have roughly equal weight, leading to new and interesting phenomena as we will show in the following. In this respect it has become clear, that local spectroscopic investigations of the mixed state is an ideal way to study the detailed nature of the superconducting state in MgB2 at a microscopic level.

In this paper we summarize the present status of our scanning tunneling spectroscopy measurements, made possible by the recent availability of high quality MgB2 single crystals. The results include spectra obtained by tunneling parallel as well as perpendicular to the crystalline $`c`$-axis, showing how the tunnel direction changes the coupling to the two bands. In addition, vortex imaging was performed with the tunnel current and magnetic field parallel to the $`c`$-axis. The vortices are found to have a number of remarkable properties: An absence of localized states, a very large vortex core size compared to the estimate based on $`H_{\text{c2}}`$, and a strong core overlap. Some of these results were reported earlier .

## 2 Experimental

The scanning tunneling spectroscopy (STS) experiments were performed using a home built scanning tunneling microscope (STM) installed in a 3He, UHV cryostat containing a 14 T magnet . The tunneling measurements were done using electrochemically etched iridium tips, with the bias voltage applied to the sample, i.e. a positive bias corresponds to probing the empty states and negative bias to the occupied states respectively above and below the Fermi level. A tunnel resistance of $`R_\text{t}=0.20.4`$ G$`\mathrm{\Omega }`$ was used, and the differential conductance measured using a standard ac lock-in technique.

Single crystals of MgB2 were grown using a high pressure method described elsewhere , yielding platelike samples with the surface normal parallel to the crystalline $`c`$-axis. The surface area of the crystals are roughly $`0.25\times 0.25`$ mm2, and the thickness varies between a few up to $`100`$ microns. The critical temperature is typically $`T_\text{c}=3738.6`$ K, with a sharp transition, $`\mathrm{\Delta }T_\text{c}=0.40.6`$ K, measured by SQUID magnetometry . Measurements were done both with the tunnel current, $`I_\text{t}`$, parallel and perpendicular to the $`c`$-axis. For $`I_\text{t}c`$, the spectroscopy was done on the surface of an as grown single crystal. In the case of $`I_\text{t}c`$, a relatively thick crystal was cracked to expose a clean surface and immediately mounted and inserted into UHV.

The upper critical field of MgB2 is highly anisotropic , with zero temperature extrapolations of respectively $`H_{\text{c2}}^\text{c}=3.1`$ T and $`H_{\text{c2}}^{\text{ab}}=23`$ T for single crystals , even though the latter value is subject to some uncertainty. Using the Ginzburg-Landau (GL) expression for $`H_{\text{c2}}=\varphi _0/(2\pi \xi ^2)`$, where $`\varphi _0=h/2e`$ is the flux quantum, yields an in-pane coherence length, $`\xi _{\text{GL}}^{\text{ab}}=10`$ nm. Due to the large anisotropy, $`\xi _{\text{GL}}^\text{c}=\xi _{\text{GL}}^{\text{ab}}\times H_{\text{c2}}^\text{c}/H_{\text{c2}}^{\text{ab}}12`$ nm. An estimate of the mean free path, based on the measured residual resistivity  and specific heat , and the calculated Fermi velocity , gives $`l=50100`$ nm, and a first estimate is thus that the samples are in the moderately clean limit . However, the two-band nature of MgB2 complicates the picture since $`H_{\text{c2}}`$ is dominated by the two-dimensional $`\sigma `$-band which gives rise to the large $`H_{\text{c2}}`$-anisotropy, while the transport properties depends on both bands. A detailed analysis of the thermal conductivity by Solugubenko et al. concluded that the mean free path is roughly equal for the two bands with $`l80`$ nm .

## 3 Results

In the following we present the results of our measurements on MgB2 single crystals.

### 3.1 Dependence on the tunnel current direction

The observation of a single or both superconducting gaps depends on the orientation of the sample, since tunneling along different directions changes the coupling to the two bands . Roughly speaking, only bands with a component of the Fermi velocity parallel to the tunneling direction are observed. In particular, this means that coupling to the 2D $`\sigma `$-band is highly suppressed for tunneling parallel to the $`c`$-axis, whereas the 3D $`\pi `$-band will be probed for all tunneling directions.

In Fig. 1(a) we show a superconducting spectrum with $`I_\text{t}c`$.

This is an average of 40 spectra obtained along a 100 nm path, which shows perfect homogeneity. One observes a single gap with coherence peaks at $`\pm 2.9`$ meV associated with the $`\pi `$-band, and additional weak shoulders at $`\pm 6`$ meV as indicated by the arrows. The shoulders are an admixture of the gap in the $`\sigma `$-band.

In Fig. 1(b) we show an average spectrum obtained for $`I_\text{t}c`$. For this direction of tunneling, both superconducting gaps are present. The position of the low-energy coherence peak is unchanged as shown by the dashed lines. The additional peak arising from the $`\sigma `$-band is located at $`\pm 7.3`$ meV. This is the first report of tunneling measurements on single crystals of MgB2 for both $`I_\text{t}c`$ and $`I_\text{t}c`$, which allows a direct correlation between the tunneling direction and the observed gap(s). Our results agree with those of Iavarone et al., who investigated a number of single grains with different, but unknown absolute orientations .

For both tunneling directions the low zero bias conductance indicates a high quality tunnel junction and a low noise level. True vacuum tunneling conditions were confirmed by varying the tunnel resistance, $`R_\text{t}`$, and verifying that the spectra normalized to the conductance outside the superconducting gap collapse onto a single curve.

### 3.2 Temperature dependece for $`I_\text{t}c`$

Tunneling parallel to the $`c`$-axis, we measured the superconducting spectrum for a number of temperatures between 320 mK and $`38.8`$ K. As discussed above this probes mainly the $`\pi `$-band. The results are presented in Fig. 2(a), and shows how the gap is gradually suppressed and finally closes at $`T_\text{c}`$.

We have fitted the measured spectra to the conductance calculated using the BCS expression for the density of states, including a finite quasiparticle lifetime, $`\mathrm{\Gamma }`$ :

$$N(ϵ,\mathrm{\Gamma })=N_0\left|\text{Re}\frac{ϵi\mathrm{\Gamma }}{\sqrt{(ϵi\mathrm{\Gamma })^2\mathrm{\Delta }^2}}\right|.$$
(1)
The conductance is given by

$$\sigma (V)=\frac{dI}{dV}𝑑ϵN(ϵ,\mathrm{\Gamma })\left(\frac{f(ϵ+eV)}{V}\right),$$
(2)
which we furthermore have convoluted by a Gaussian of width $`\omega `$ to take into account experimental broadening. The fits are shown in Fig. 2(a) and are found to be in good agreement with the measurements, except at low temperatures where the height of the coherence peaks is overestimated. The extracted values of the gap are shown in Fig. 2(b), yielding $`\mathrm{\Delta }_\pi (0)=2.22.3`$ meV. We find that $`\mathrm{\Delta }_\pi (T)`$ is in excellent agreement with the temperature dependence of the BCS $`\mathrm{\Delta }(T)`$  when the zero temperature gap is scaled to $`2.3`$ meV. The lifetime parameter $`\mathrm{\Gamma }=0.1`$ meV at low-$`T`$ and increases with temperature. The experimental broadening $`\omega =0.5`$ meV at base-$`T`$ and increases only slightly to $`0.8`$ meV at 10 K. For comparison the ac excitation used for the measurements was $`0.4`$ meV RMS. Above 10 K the effects of $`\mathrm{\Gamma }`$ and $`\omega `$ becomes identical, and the latter is kept fixed.

### 3.3 Vortex imaging in the $`\pi `$-band

We now turn to measurements in an applied magnetic field. These measurements were done with $`I_\text{t}𝑯c`$, and hence selectively probes the $`\pi `$-band. In a type-II superconductor such as MgB2, the magnetic field penetrates into the sample in the form of vortices each carrying one flux quantum, and generally arranged in a periodic array: the vortex lattice. In the core of each vortex, superconductivity is suppressed within a radius roughly given by the coherence length, $`\xi `$. The vortex spacing is determined by the applied field and the flux quantization, and in the case of a hexagonal vortex lattice it is $`d=(2/\sqrt{3\varphi _0/H})^{1/2}`$.

Necessary prerequisites for obtaining vortex lattice images is a high sample homogeneity and a relatively flat sample surface. We verified that the first requirement is fulfilled by measuring perfectly reproducible spectra in zero field at a large number of positions along a several thousand ångstrøm line. A topographic STM image of the sample surface is shown in Fig. 3.

The surface is seen to be very smooth, with an indication of a granular structure with a characteristic size of 5 - 10 nm. No effect of this “granularity” is observed in the spectra.

The magnetic fields were applied at 2 K, and the system allowed to stabilize for at least a few hours. After this time no vortex motion was observed, indicating a fast relaxation and hence low vortex pinning in the crystal. In Fig. 4 we show STS images of vortices induced by three different applied fields.

The images were obtained by measuring the differential conductance at zero bias. Low values of the conductance correpond to superconducting areas, and high values to the vortex cores. Figs. 4(b) and (c) show a well ordered hexagonal vortex lattice, with a lattice constant that corresponds to the applied field within 10%.

The low field in Fig. 4(a) is equivalent to a separation, $`d=220\text{ nm}\xi _{\text{GL}}`$. The vortices can therefore be considered as isolated from each other. Such isolated vortices are expected to contain localized quasiparticle states (Andreev bound states), which should show up as a zero bias conductance (ZBC) peak at the vortex centre , provided that the sample is sufficiently clean to prevent these from being smeared out by scattering. We have measured the evolution of the spectra at a large number of positions along a trace across the vortex core. This is shown in Fig. 5(a).

Contrary to expectations, we find that the normalized ZBC increases to one with no indication of any localized states. Instead, the spectra at the centre of the vortex are absolutely flat, with no excess spectral weight at or close to zero bias. This absence of localized states is striking, considering that $`l=510\times \xi _{\text{GL}}`$. However, as we will show below, the coherence length in the $`\pi `$-band is approximately 50 nm. This is much larger than the GL estimate based on $`H_{\text{c2}}`$, and equal to only 1 – 2 times the mean free path. Nonetheless, for comparison systematic studies of Nb1-xTaxSe2 with $`x=00.2`$, showed that even for $`\xi /l1`$ some excess weight close to zero bias was observed .

The spatial extent of the vortex is shown in Fig. 5(b), where we have plotted the normalized ZBC, $`\sigma ^{}(x,0)`$ from the vortex trace. It is immediately clear that the spatial extension of the vortex core is much larger than the 10 nm estimated from $`H_{\text{c2}}`$. We find that the vortex profile can be fitted by one minus the GL expression for the superconducting order parameter:

$$\sigma ^{}(x,0)=\sigma _0^{}+(1\sigma _0^{})\times (1\mathrm{tanh}x/\xi ),$$
(3)
where $`\sigma _0^{}=0.068`$ is the normalized ZBC measured in zero field. The fit yields a coherence length, $`\xi =\xi _\pi =49.6\pm 0.9`$ nm. Using the GL expression to calculate the upper critical field with this value of the coherence length, yields $`H_{\text{c2}}^\pi =0.13`$ T. This value is in strong contrast to the fact that we have performed vortex lattice imaging up to $`0.5`$ T, and shows that a separate $`H_{\text{c2}}`$ for the $`\pi `$-band does not exist.

The large $`\xi _\pi `$ has however other consequences. In Fig. 6 we have plotted the normalized ZBC between the vortices.

From this it is clear that even modest fields rapidly suppress superconductivity in the region between the vortices, seen both by an increase of the ZBC and by a suppression of the coherence peaks outside the vortex cores. This is consistent with earlier point contact spectroscopy measurements , with the addition that we resolve the local behaviour on a microscopic scale. Such a behaviour is only expected for fields close to $`H_{\text{c2}}`$, and corresponds to a significant core overlap.

## 4 Discussion

The vanishing of $`\mathrm{\Delta }_\pi `$ at the bulk $`T_\text{c}`$ and the sustained superconductivity despite a large vortex core overlap for $`H>0.2`$ T leads to the conclusion that superconductivity in the $`\pi `$-band is induced by the $`\sigma `$-band, either by interband scattering, or Cooper pair tunneling . This means that by itself the $`\pi `$-band would either be non-superconducting or have a very low $`T_\text{c}`$ and upper critical field. Consequently the observed behaviour reflects the state in the $`\sigma `$-band by an interband proximity effect, along the lines of recent theoretical work . This is also consistent with estimates of the coherence lengths, using the BCS expression $`\xi _0=\mathrm{}v_\text{F}/(\pi \mathrm{\Delta }(0))`$ and considering each band separately. Taking the calculated average Fermi velocity in the $`ab`$-plane for the $`\pi `$-band, $`v_\text{F}^\pi =5.35\times 10^5`$ m/s , and the measured gap value $`\mathrm{\Delta }_\pi (0)=2.2`$ meV we get $`\xi _0^\pi =51`$ nm, in excellent agreement with $`\xi _\pi `$ obtained from the vortex profile. A similar analysis for the $`\sigma `$-band, using $`v_\text{F}^\sigma =4.4\times 10^5`$ m/s  and $`\mathrm{\Delta }_\sigma (0)=7.1`$ meV  yields $`\xi _0^\sigma =13`$ nm. This agrees with the coherence length obtained from $`H_{\text{c2}}`$, and reinforces the conclusion that it is mainly the $`\sigma `$-band which is responsible for superconductivity in MgB2, and thus determines the macroscopic parameters $`T_\text{c}`$ and $`H_{\text{c2}}`$. One must however bear in mind that this is a very simple-minded analysis. A detailed theoretical treatment of two weakly coupled superconducting bands in the presence of vortices, and the possibility of having have different coherence lengths, have to our knowledge not yet been presented.

### 4.1 Comparison to specific heat measurements

The vortex core overlap explains the strongly nonlinear field dependence of the electronic specific heat ($`T`$-linear term), $`\gamma `$ . In type-II superconductors core overlap is usually negligible, and in the simplest approximation each vortex creates the same number of quasiparticles at the Fermi surface, hence contributing by the same amount to the specific heat. In that case $`\gamma =\gamma _\text{n}\times H/H_{\text{c2}}`$, where $`\gamma _\text{n}`$ is the electronic specific heat in the normal state. However, in the case of MgB2, with strong core overlap in the $`\pi `$-band, the isolated vortex assumption is violated. Constructing a very simple model for the core overlap by

$$\sigma ^{}(𝒓,0)=\sigma _0^{}+(1\sigma _0^{})\times \left(1\mathrm{\Pi }_i\mathrm{tanh}\frac{|𝒓𝒓_i|}{\xi _\pi }\right),$$
(4)
where $`𝒓_i`$ are the vortex positions for a hexagonal lattice with a density corresponding to the applied field, we can calculate the ZBC at any position in the vortex lattice unit cell. In Fig. 7 we compare the calculated conductance at the midpoint between three vortices with the measured “bulk” ZBC.

This shows a very good agreement, especially taking into account that there are no free parameters in the calculation: $`\xi _\pi `$ is determined from the vortex profile, and $`\sigma _0^{}`$ is given by the finite ZBC in the zero field measurement. The $`\pi `$-band contribution to the specific heat can now be calculated by averaging the normalized ZBC in one vortex lattice unit cell, $`\gamma _\pi =\gamma _\text{n}^\pi \sigma ^{}(𝒓,0)`$ . For this calculation we have set $`\sigma _0^{}`$ equal to zero. On the other hand, the $`\sigma `$-band can be described by the usual linear field dependence $`\gamma _\sigma =\gamma _\text{n}^\sigma \times H/H_{\text{c2}}`$. Adding the terms gives $`\gamma =\gamma _\pi +\gamma _\sigma `$, where $`\gamma _\text{n}^\pi /\gamma _\text{n}^\sigma `$ is the relative weigth of the two bands. The calculated field dependence of $`\gamma `$ is shown in Fig. 7, in perfect agreement with the measured specific heat for polycrystalline MgB2, using $`\gamma _\text{n}^\pi /\gamma _\text{n}^\sigma =0.55/0.45`$. The distribution of weight between the two bands is in agreement with other theoretical and experimental estimates . Furthermore, the successful correlation with specific heat measurements makes us confident that our results reflect bulk properties of MgB2.

## 5 Summary

In summary, we have presented STS data on the $`\pi `$-band in MgB2, including the first vortex imaging in this material. We have demonstrated the absence of localized states in the vortex core, a very large vortex core size and a strong core overlap. The data presents a striking experimental demonstration of the fundamentally different microscopic properties of the two bands in MgB2.

## Acknowledgements

We acknowledge valuable discussions and communication of data prior to publication with F. Bouquet, M. Iavarone, A. Junod and Y. Wang, and thank B. W. Hoogenboom and I. Maggio-Aprile for sharing their experience in STM/STS. This work was supported by Swiss National Science Foundation.


text:
# Spectral Lines from Rotating Neutron Stars

## 1 Introduction

Thermal emission from the surface of a neutron star carries signatures of its strong gravitational field, which become apparent in observations of both its spectral and timing properties. Such measurements, therefore, can be used in principle to infer the masses and radii of neutron stars. Over the past three decades, multiple attempts have been made to constrain the stellar equation of state using observations of the thermal emission from bursting (Lewin, van Paradijs, & Taam 1995), quiescent (e.g., Rutledge et al. 1999), or isolated neutron stars (e.g., Pons et al. 2002; Braje & Romani 2002), the X-ray pulse profiles of rotationally powered pulsars (Page 1995; Pavlov & Zavlin 1997), the high amplitudes of oscillations observed during thermonuclear bursts (Nath, Strohmayer, & Swank 2001), and the frequencies of observed quasi-periodic oscillations (Miller, Lamb, & Psaltis 1998).

Of all the possible methods of measuring the radius of a neutron star or a mass-radius combination, the one that suffers the least from systematic uncertainties and measurement errors is using the gravitational redshift of atomic spectral lines. For a slowly rotating, spherically symmetric neutron star, the gravitational redshift gives directly the stellar compactness (i.e., the ratio $`R_{\mathrm{NS}}/M_{\mathrm{NS}}`$). This method has received a lot of attention recently with the launch of X-ray telescopes with high spectral resolution (such as Chandra and XMM-Newton) and the discovery of thermal emission from nearby, isolated neutron stars (see Becker & Pavlov 2002 for a review). A number of observations of neutron stars have already been carried out, which yielded a potential detection of broad, redshifted absorption lines from the pulsar 1E 1207.4–5209 (Sanwal et al. 2002; Mereghetti et al. 2002). Other observations of isolated neutron stars have typically resulted in featureless X-ray spectra (e.g., Paerels et al. 2001; Drake et al. 2002; Marshall & Schulz 2002).

Apparently, not all targets are created equal. The detection of atomic spectral lines in X-rays requires both heavy metals to be present in the neutron-star atmosphere and the surface layers to have high temperatures for significant thermal emission to be generated. These two requirements can be met most easily in either neutron stars that are young or in ones that are weakly magnetic and accreting steadily from a binary companion. Young neutron stars emit thermally the heat released during their formation. They may also possess heavy-element atmospheres if significant light-element fallback did not occur during the supernova: their short lifetimes and strong magnetic fields render unlikely a significant accumulation of hydrogen rich material from the interstellar medium that could suppress atomic lines. In the case of bursters, the heavy elements in their atmospheres are continually replenished by accretion and the thermonuclear flashes provide large amounts of thermal energy.

Both types of neutron stars that are prime candidates for the detection of spectral lines are fast rotators. (We do not consider here magnetars, for which the presence of ultrastrong magnetic fields introduces large uncertainties in calculating the rest energies of atomic lines). The spin frequencies of known pulsars with ages $`<10^4`$ yr is between $`565`$ Hz (see, e.g., Becker & Pavlov 2002); the inferred spin frequencies of bursters is between $`270620`$ Hz (Strohmayer 2001). These high spin frequencies introduce several relativistic effects such as Doppler boosts, strong self-lensing, frame-dragging and differential gravitational redshift arising from the stellar oblateness. All of them alter the line profiles observed at infinity.

In this Letter, we show the effects of relativistic Doppler boosts and strong gravitational lensing on the width and asymmetry of line profiles originating from the surfaces of rotating neutron stars. We then investigate the systematic uncertainties introduced by these effects in inferring the compactness of neutron stars.

## 2 The Effects of Rotation on the Observed Line Profiles

Neutron stars are the most compact stellar objects and the rotational velocities at their surfaces can reach an appreciable fraction of the speed of light. Therefore, rotational effects on the line profiles originating from the neutron star surface are qualitatively different from the case of a Newtonian slowly spinning star. In particular, the shape of the line profiles observed at infinity is altered in four ways.

First, relativistic Doppler boosts give rise to an asymmetry in the spectral line profiles while broadening them. Second, strong gravitational lensing of surface emission by the neutron star alters the relative contribution of surface elements with different line-of-sight velocities to the line profile. Third, frame dragging in the rotating spacetime of the neutron star affects the photon trajectories and thus the observable surface elements. Finally, the stellar oblateness caused by the rapid spin introduces a difference in the gravitational redshifts of lines that are generated at the rotational equator and at the poles.

In this Letter, we use the numerical methods described in Muno, Özel, & Chakrabarty (2002 and references therein) to calculate spectral line profiles taking into account the first two of the effects discussed above. Investigating the latter two requires the calculation of numerical spacetimes for specific equations of state of neutron-star matter and will be addressed in a forthcoming paper. The results presented here are accurate for spin frequencies smaller than the Keplerian frequency at the surface of a neutron star, i.e., less than a few hundred Hz. It is important to note that for larger spin frequencies, the ellipticity of the neutron stars induced by the rotation can be as large as 0.3, depending on the equation of state (Cook et al. 1994), yielding an additional broadening as large as $`15`$%.

## 3 Results

In the calculations presented here, we show for clarity emission lines; absorption lines are affected in the same way due to the linear character of the equations. Moreover, for numerical reasons, we assume an intrinsic fractional line width of $`0.01`$, which is much smaller than the Doppler width for the cases shown here. Finally, we assume that the observer is on the rotational equator of the star in order to show the maximum rotational effects.

The line profiles measured at infinity coming from the surface of a $`1.4M_{},10`$ km neutron star are shown in Figure 1 for different values of the spin frequency. As the spin frequency is increased, the lines become weaker and broader, and their peak emission shifts towards larger energies. In this figure, the entire neutron star is assumed to be emitting uniformly.

The asymmetry of the line profiles become more prominent when only a fraction of the stellar surface contributes to the line emission. Such a configuration is likely in the case of young neutron stars, which may have lateral composition and temperature gradients owing to their strong magnetic fields. Nonuniform emission is almost certainly relevant also in the case of bursters as indicated by the observations of large amplitude flux oscillations during the thermonuclear bursts (Strohmayer 2001). As Figure 2 shows, in the case of non-uniform emission, phase-integrated line profiles acquire a doubly peaked character, with a brighter blue wing. In this figure, the emitting region, which has an angular radius of $`20^{}`$, is assumed to be at the rotational equator. This configuration yields the largest separation of the two peaks in the line profile because of the largest line-of-sight velocities. Note that due to the strong gravitational lensing, the two maxima of the doubly peaked profile do not correspond to the line emission at rotational phases $`\varphi =\pm \pi /2`$, where $`\varphi =0`$ denotes the phase when the center of the emitting region is aligned with the observer. As a consequence, the separation of the peaks is not proportional to $`2\gamma \mathrm{\Omega }R_{\mathrm{NS}}/c`$, where $`\gamma `$ is the Lorentz factor corresponding to the rotational velocity and $`\mathrm{\Omega }`$ is the angular velocity of the neutron star.

The line-of-sight velocities, and thus the separation of the peaks in the phase-averaged line profiles, are also affected by the actual size of the neutron star, for a given stellar compactness. In Figure 3, we show the line profiles measured at infinity that originate from two neutron stars of the same compactness, $`R/M=2.41G/c^2`$, but of different radii. Clearly, the larger peak separation corresponds to the larger neutron-star radius. This effect allows in principle an independent determination of the mass and radius (modulo emission geometry) of a neutron star of known spin frequency, given the shape and overall redshift of atomic spectral lines.

## 4 Discussion and Conclusions

We studied spectral line profiles from rotating neutron stars taking into account the effects of relativistic Doppler boosts and strong gravitational lensing. We showed that the line profiles are broad, as expected, and also significantly asymmetric. The asymmetry becomes more prominent when the surface emission is non-uniform. Our results have a number of implications for the current searches for gravitationally redshifted line features in the spectra of neutron stars.

First, the large widths and suppressed strengths of the rotationally broadened lines make their detection difficult. This may be able to account for the featureless spectra of a number of isolated neutron stars observed with Chandra and XMM-Newton, such as RXJ 1856–3754 (Braje & Romani 2002; Zavlin & Pavlov 2002). Correspondingly, if narrow line features are detected from rapidly rotating neutron stars, e.g., bursters, they could not have originated from the neutron star surface, unless the emission is restricted to the rotational pole. Searches for lines in such sources should take into account these relativistic effects.

Finally, the asymmetry of the line profiles introduces significant systematic uncertainties in measuring the compactness of a neutron star using gravitational redshifts. As an example, Figure 4 shows that if the peak of the line is used in measuring an apparent redshift, the resulting compactness of the neutron star will be significantly overestimated even when the entire surface is emitting. For the inferred spin frequencies of bursters, in the absence of realistic models, the systematic uncertainties can be as large as $`10\%`$, which are larger than the $`5\%`$ accuracy required to distinguish between the different equations of state (Prakash & Lattimer 2000).

F.Ö. acknowledges support by NASA through Hubble Fellowship grant HF-01156 from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. D.P. acknowledges the support of NSF grant PHY-0070928.


text:
# Soft and hard QCD This talk is taken largely from my new book with Donnachie, Dosch and Nachtmann[1]

## 1 Introduction

There are many approaches to high-energy scattering, among them

$``$Regge theory

$``$dipole models

$``$stochastic vacuum models

$``$saturation models

$``$semiclassical approach

$``$effective field theory

$``$DGLAP

$``$BFKL

They all use different language, but there are many links between them. At present, none of them offers an agreed fundamental explanation for the very striking discovery at HERA, that at high $`Q^2`$ the $`\gamma ^{}p`$ total cross section rises dramatically with increasing energy $`W`$. This is seen in figure 3. At small $`Q^2`$ the rise is compatible with that seen in hadron-hadron total cross sections, $`(W^2)^{ϵ_1}`$ with $`ϵ_10.08`$, but at high $`Q^2`$ the effective power is close to 0.4. Leading-order BFKL predicted this, but unfortunately there is a huge correction in next-to-leading order in $`\alpha _s`$.

## 2 Difficulty with DGLAP

Most fits to the data achieve the rising power from DGLAP evolution, but they do so by making an expansion of the splitting matrix that is mathematically illegal.

The singlet DGLAP equation introduces a two-component quantity

$$𝐮(x,t)=\left(\begin{array}{c}x_f(q_f+\overline{q}_f)\\ xg(x,t)\end{array}\right)t=\mathrm{log}(Q^2/\mathrm{\Lambda }^2)$$
(1)
If one Mellin transforms with respect to $`x`$

$$𝐮(N,Q^2)=_0^1𝑑xx^{N1}𝐮(x,Q^2)$$

$$𝐏(N,\alpha _s(Q^2))=_0^1𝑑zz^N𝐏(z,\alpha _s(Q^2))$$
(2)
the equation becomes very simple:

$$\frac{}{t}𝐮(N,Q^2)=𝐏(N,\alpha _s(Q^2))𝐮(N,Q^2)$$
(3)
The standard approach is to expand the splitting matrix $`𝐏`$ is powers of $`\alpha _s`$, but this is invalid at small $`N`$. Compare

$$\sqrt{N^2+\alpha _s}N=\alpha _s/2N\alpha _s^2/8N^3+\mathrm{}$$
(4)
Here each term in the expansion is singular at $`N=0`$ but the function itself is not: the expansion is illegal when $`N`$ is small. I will discuss later how one might partially overcome this difficulty.

## 3 Dipole model

Figure 3 shows the virtual forward compton amplitude. Each virtal photon couples to a quark-antiquark pair which may be regarded as a colour dipole. This leads to

$$\sigma _{T,L}^{\gamma ^{}p}(x,Q^2)=$$

$$d^2R𝑑z\psi _{T,L}^{}(Q,R,z)\sigma ^{\text{dip}}(x,R)\psi _{T,L}(Q,R,z)$$

(5)
Here, $`T,L`$ denote the polarisation of the photon. $`\psi _{T,L}(Q,R,z)`$ is the wave function at the vertex that couples it to the $`q\overline{q}`$ pair; it depends on the transverse separation $`R`$ of the pair and on the longitudinal momentum fraction $`z`$ of the quark. $`\sigma ^{\text{dip}}(x,R)`$ is the cross section for the interaction of the colour dipole with the target proton.

The literature includes many different choices for $`\sigma ^{\text{dip}}(x,R)`$. Figure 6 shows a few of them. There is no single dipole model!

Nowadays, it is popular to combine the dipole model with the notion of saturation. In hadron-hadron cross sections, such as $`\sigma ^{pp}`$, the Froissart-Martin-Lukaszuk bound is familiar:

$$\sigma ^{pp}(s)\frac{\pi }{m_\pi ^2}\mathrm{log}^2(s/s_0)$$
(6)
It is not a material constraint, because it gives an upper limit of several barns at LHC energies! It is derived from unitarity:

$$\text{Im }a_l(s)|a_l(s)|^2$$
(7)
For $`\gamma p`$ scattering there is no similar inequality because the unitarity relation does not contain an elastic term: to lowest order in $`\alpha _{\text{EM}}`$ only hadronic final states are included. In principle, $`\sigma ^{\gamma p}(W)`$ can get huge at large $`W`$ and $`F_2(x,Q^2)`$ can get huge at small $`x`$ Nevertheless, many people believe that a Froissart-like bound is saturated at an accessible $`W`$ or $`x`$, and they implement this by writing an eikonal-like form

$$\sigma ^{\text{dip}}(x,R)=1\mathrm{exp}\left(R^2/4R_0^2(x)\right)$$
(8)
Now

$$1\mathrm{exp}\left(R^2/4R_0^2(x)\right)=\frac{R^2}{4R_0^2(x)}\frac{R^4}{32R_0^4(x)}+\mathrm{}$$

(9)
and it is natural to identify the first term in this expansion with the most elementary exchange, taken to be two-gluon exchange. Figure 6 shows this for the reaction $`\gamma ^{}pVp`$. At the bottom of the figure there is the proton’s gluon structure function $`xg(x,\mu ^2)`$, and so we have for the dipole cross section

$$\sigma ^{\text{dip}}(x,R)=$$

$$\sigma _0\left\{1\mathrm{exp}\left(\frac{\pi ^2R^2\alpha _s(\mu ^2)xg(x,\mu ^2)}{3\sigma _0}\right)\right\}$$
(10)
with just one parameter $`\sigma _0`$. Because $`xg(x,\mu ^2)`$ obeys DGLAP evolution, this model combines the dipole model with both saturation and DGLAP, and it can give a good fit to experiment. However, there are many ways in which one can successfully fit the DIS data.

## 4 Stochastic vacuum model

The stochastic vacuum model starts from the familiar vacuum gluon condensate

$$0|:g^2F^{\mu \nu }(x)F_{\mu \nu }(y):|0|_{y=x}=M_c^4$$
(11)
with $`M_c`$ a few hundred MeV, and generalises this relation to $`yx`$. This introduces a vacum correlation length. Some rather technical manipulations are needed, for example using the nonabelian Stokes theorem. A particular realisation of the dipole model results, where the soft pomeron is generated from multigluon exchange. The model successfully relates total cross sections to hadron sizes, as is seen in figure 6

## 5 Soft cross sections

Regge theory provides remarkably simple fits to data for all hadron-hadron total cross sections. An example is shown in figure 9. One needs just two powers of $`s`$. One is close to $`1/\sqrt{s}`$ and is identified as resulting from $`\rho ,\omega ,f_2,a_2`$ exchange. The other is close to $`s^{0.1}`$ and its origin is unknown; to give it a name, we say that this term results from soft-pomeron exchange. The extrapolation of the parametrisation shown in figure 9 gives 108 mb at $`\sqrt{s}=20`$ TeV. Although, as I have explained, one can never hope to achieve an energy at which the Froissart bound becomes a relevant constraint, it has often led people to prefer to parametrise the rising component of the cross section with a log$`{}_{}{}^{2}s`$ term rather than a power. It is interesting that the most recent such fit predicts that the cross section at $`\sqrt{s}=20`$ TeV will be some 10 mb greater than given by the simple power fit.

Note the highest-energy points in figure 9: the CDF point is significantly higher than the E710 point. If CDF were to turn out to be correct, this would signal the onset of some new term which would significantly increase the cross section measured at the LHC. The question whether such a hard term is present is of some fundamental importance for the interpretation of the HERA data in figure 3. These data show clearly the presence of a hard term at high $`Q^2`$ and it is generally agreed that it should be understood through pQCD evolution. But does the evolution generate the term, or merely enhance its importance as $`Q^2`$ increases? I am fairly sure that it is the latter that is the case. If this is true, the term should be present in $`\gamma p`$ collisions already at $`Q^2=0`$. While there is room for such an additional term in the data shown in figure 3, the error bars are too large to decide. The LEP data for the $`\gamma \gamma `$ total cross section are similarly unclear. Figure 9 shows the data. The L3 experiment presents two sets of points, corresponding to two different Monte Carlos, which are needed to correct for the fact that the detector’s acceptance is such that only a small fraction of the interactions are visible. The curve represents a sum of the same two powers as the curves in figure 9. The clearest indication that a hard term is indeed present at small $`Q^2`$ is in the ZEUS data for the charm structure function of the proton. As figure 9 shows, already at $`Q^2=1.8`$ GeV2 the rise with increasing $`1/x`$ of $`F_2^c(x,Q^2)`$ is as rapid as that of the complete $`F_2(x,Q^2)`$ at large $`Q^2`$. The same is even true at $`Q^2=0`$. I return to this very important point later.

## 6 Elastic scattering

Regge theory provides a very simple extension to elastic scattering of the total-cross-section fit of figure 9. At high $`s`$, where only the soft-pomeron term survives,

$$\frac{d\sigma ^{pp}}{dt}\frac{d\sigma ^{\overline{p}p}}{dt}\frac{(3\beta _𝐏F_1(t))^4}{4\pi }\left(\frac{s}{s_0}\right)^{2\alpha _𝐏(t)2}$$
(12)
where

$$\alpha _𝐏(t)=1.08+\alpha ^{}t$$
(13)
and $`F_1(t)`$ is the proton’s Dirac elastic form factor. The value $`\alpha ^{}=0.25`$ GeV-2 is fixed by fitting to very accurate ISR data at very small $`t`$. The form (12) then successfully predicts the data at much higher energy. See figure 12.

With no free parameters, we may extend this to $`\pi p`$ elastic scattering. The pion has only two valence quarks, so we replace (12) with

$$\frac{d\sigma }{dt}=\frac{(2\beta _𝐏F_\pi (t))^2(3\beta _𝐏F_1(t))^2}{4\pi }\left(\frac{s}{s_0}\right)^{2\alpha _𝐏(t)2}$$
(14)
Figure 12 shows data for the pion form factor; they fit well to $`F_\pi (t)=1/(1t/m_0^2)`$ with $`m_0^2=0.5`$ GeV2. This leads to the zero-parameter fit shown in figure 12.

Regge theory is remarkably successful!

## 7 Glueballs

Although we do not understand the origin of the soft pomeron, there is a wide feeling that it is just gluon exchange. If that is so, and if its trajectory really is straight, as written in (13), then the value of $`t`$ for which it passes through 2 should be the square of the mass of a $`2^{++}`$ glueball. The WA91 collaboration has a $`2^{++}`$ candidate of exactly the right mass: see figure 13.

Nowadays it is known that Regge trajectories corresponding to ordinary particles are accompanied by daughter trajectories. These are trajectories separated by an integral number of units from the parent trajectory. An example is the $`f`$ family, shown in figure 16. The existence of daughters was predicted from Regge theory at a time when little was known about the meson spectrum. One would expect the pomeron trajectory to have daughters too. The search for glueballs is very important to give more understanding about the pomerons – the hard pomeron is probably associated with glueballs too.

## 8 Odderon

The minimum number of gluons needed to model the pomerons is two, because they represent colourless even-parity exchange. With three gluons, one can model colourless odd-parity exchange, called odderon exchange. There is a clear sign of odderon exchange in $`pp`$ and $`\overline{p}p`$ elastic-scattering data at large $`t`$, but the mystery is that so far odderon exchange has not been identified at $`t=0`$.

Figure 16 shows ISR data for $`pp`$ elastic scattering. There is a very striking dip at $`|t|1.4`$ GeV2. The very last week of running of the CERN ISR showed that $`\overline{p}p`$ elastic scattering is different: the dip is filled in, as is seen in figure 16.

Beyond the dip, the data in the ISR energy range fit very well to perturbative 3-gluon exchange calculated in leading order: see figure 18. That is, they are independent of $`s`$ and vary as $`t^8`$. There are many unanswered questions about this: why does this simple behaviour set in already at such a small $`t`$, why is it not significantly altered by higher-order perturbative corrections, and are the data really energy-independent? It will be interesting to check this at LHC energies. A possibility is that triple-gluon exchange will be replaced at higher energies with triple-hard-pomeron exchange, so that the large-$`t`$ differential cross section actually rises with increasing energy.

## 9 Soft diffraction dissociation

We say that a $`pp`$ scattering event is diffractive if one of the protons loses only an extremely small fraction $`\xi `$ of its momentum. In diffraction dissociation, the other proton breaks up. The mechanism by which this occurs is supposed to be pomeron exchange, as is seen in figure 18. Although the pomeron is not a particle, it is as if it collides with the second proton, and one talks of the pomeron-proton cross section. This cross section should be similar to hadron-hadron-scattering cross sections. In particular, it should rise with energy. But the $`pp`$ diffraction-dissociation data show no sign of this rise. Figure 19 shows data at $`\sqrt{s}=23`$ and 630 GeV and the curves are what is expected to result from the rising pomeron-proton cross section. There is no agreed explanation for this discrepancy, though there have been suggestions that, for some reason, the pomeron flux does not show the expected behaviour with increasing energy.

## 10 Exclusive Higgs production

The exclusive process $`pppHp`$, where both the final-state protons emerge with very high longitudinal momentum, has been discussed extensively over the last ten years or so. This reaction should be generated by double pomeron exchange: see figure 21.

Interest in the reaction has been revived by the suggestion that it might be a good way to discover the Higgs. Higgs searches in hadronic collisions have big background problems, but it is argued that, by measuring the momenta of the final-state protons in figure 21 very accurately, one may determine the missing mass very accurately. So one needs to integrate the background only over a small mass range, so reducing its importance.

The argument now is whether the cross section for the process is large enough to make it visible. In particular, are screening corrections so large as to make the cross section very small? See figure 21. It has been claimed that indeed this is so and that there is a suppression of more than an order of magnitude. However, this claim is based not on a calculation of the screening itself, but on an argument that there is a very large likelihood that the two rapidity gaps in the mechanism of figure 21 will be filled in by the production of extra particles. But if one wants to calculate the amplitude for a given process, it is not relevant what else might happen. If one applied the same argument to pp elastic scattering one would conclude that the cross section should be extremely small, when in fact it is more than a quarter of the total cross section. It is true that in the eikonal model screening corrections are related to the probability of filling in the rapidity gap, but this is special to the eikonal model and there are good reasons not to trust the eikonal model. My own belief is that screening corrections as in figure 21 give a 50% suppression at most. The argument is related to that over whether Froissart-bound considerations have an important effect on how large cross sections are allowed to be. So I think that the cross section for exclusive Higgs production is an order of magnitude bigger than has been claimed recently.

## 11 Hard diffraction

The prediction that there should be a sizeable probability that hard reactions also could lead to a very fast final-state proton was first confirmed in an experiment at the CERN $`\overline{p}p`$ collider. In $`\gamma ^{}p`$ scattering the mechanism is that shown in figure 23. Although the pomeron is not a particle, it is as if the mechanism involves a hard $`\gamma ^{}`$-pomeron collision and so measures the structure function of the pomeron, just as $`\gamma ^{}`$-proton collisions measure the structure function of the proton. This has been studied extensively at HERA, where at small $`x`$ some 10% of the events are found to be diffractive.

The Tevatron experiments have measured the diffractive production of dijets and of the $`W`$. The mechanism of figure 23 suggests that the same pomeron structure function should be involved as in diffractive electroproduction, and that therefore again some 10% of dijet or $`W`$ events should be diffractive. The result that is found is an order of magnitude smaller. Again this has been blamed on the filling in of the rapidity gap by the production of additional particles. Unlike the exclusive Higgs production I have discussed before, these are inclusive processes, for which we have a much less well-defined theoretical formalism, so I do not find this explanation implausible. Nevertheless, I wonder whether things will be different at the much higher energy of the LHC.

## 12 Deep inelastic lepton scattering

I have explained that do not understand how to apply DGLAP evolution at small $`x`$. However, if we combine it with Regge theory and use an important message from the HERA data for the charm structure function $`F_2^c(x,Q^2)`$, it is possible reliably to extract the gluon structure function $`g(x,Q^2)`$ at small $`x`$. It turns out to be larger than nowadays is commonly believed. This is seen in figure 25. The most recent CTEQ and MRST structure functions agree well with each other and with those extracted by the two HERA experiments because they all use similar procedures; however, Donnachie and I believe that the old MRSG structure function is nearer the truth.

When one tries to fit data, it is usually sensible to start with the simplest assumptions and then refine them later. In its simplest form, Regge theory leads to fixed powers of $`x`$ at small $`x`$, and it turns out that two terms are enough:

$$F_2(x,Q^2)f_0(Q^2)x^{ϵ_0}+f_1(Q^2)x^{ϵ_1}$$
(15)
The second term corresponds to soft-pomeron exchange, with $`ϵ_10.08`$ determined from soft reactions. The data need a term that rises more rapidly at small $`x`$; one needs $`ϵ_00.4`$. By fitting the data at each $`Q^2`$, Donnachie and I found that a successful and economical parametrisation of the coefficent functions is provided by

$$f_0(Q^2)=A_0(Q^2)^{1+ϵ_0}/(1+Q^2/Q_0^2)^{1+ϵ_0/2}$$

$$f_1(Q^2)=A_1(Q^2)^{1+ϵ_1}/(1+Q^2/Q_1^2)^{1+ϵ_1}$$
(16)
with $`Q_03\text{ GeV}`$ and $`Q_10.8\text{ GeV}`$. To make the fit, we used real-photon data and DIS data with $`x0.001`$, so that $`Q^2`$ ranges from 0.045 to 35 GeV2. If we then simply multiply the resulting form (15) by $`(1x)^7`$, as is suggested by the dimensional counting rules, it agrees quite well with the HERA data even beyond $`x=0.1`$ and up to $`Q^2=5000`$ GeV2. This is shown in figure 25. Note that this factor $`(1x)^7`$ should not be taken too seriously; it is much too simple.

Data for the charm structure function $`F_2^c(x,Q^2)`$ have the remarkable property that, at all available $`Q^2`$, they fit to just the single hard-pomeron power of $`x`$. Further, to an excellent approximation the coupling of the hard pomeron appears to be flavour blind:

$$F_2^c(x,Q^2)=f_c(Q^2)x^{ϵ_0}$$
(17)
with

$$f_c(Q^2)=\frac{4}{9}/(\frac{4}{9}+\frac{1}{9}+\frac{1}{9}+\frac{4}{9})f_0(Q^2)=0.4f_0(Q^2)$$

(18)
So if we define a charm-production cross section

$$\sigma ^c(W)=\frac{4\pi ^2\alpha _{\text{EM}}}{Q^2}F_2^c(x,Q^2)|_{x=Q^2/(W^2+Q^2)}$$
(19)
it behaves as $`W^{2ϵ_0}`$ at all $`Q^2`$, even down to $`Q^2=0`$: see figure 27. Perturbative QCD directly relates $`F_2^c(x,Q^2)`$ to the gluon structure function, so that at small $`x`$ it too must be dominated by hard-pomeron exchange alone, even at quite small values of $`Q^2`$. This is what causes the rapid rise at small $`x`$ of the DL curve in figure 25.

## 13 DGLAP evolution

I have already explained that the usual procedure introduces spurious singularities into the splitting matrix $`𝐏`$ that appears in the DGLAP equation (3). My own belief is that $`𝐏(N,\alpha _s(Q^2))`$ has no singularities in the complex-$`N`$ plane, or at least no relevant singularities. My reason is that solving (3) would cause a singularity of $`𝐏(N,\alpha _s(Q^2))`$ to induce an essential singularity in $`𝐮(N,Q^2)`$ (that is, a nasty one). The variable $`N`$ is closely related to the orbital angular momentum $`l`$, and I was brought up to believe that matrix elements such as $`𝐮(N,Q^2)`$ do not have essential singularities in the complex $`l`$-plane. This point of view contrasts with that of those who believe that the value of $`ϵ_0`$ is associated with a singularity of $`𝐏(N,\alpha _s(Q^2))`$ and may even be calculated, perhaps by refining the BFKL approach. I think that very probably $`ϵ_0`$ is a nonperturbative quantity that at present cannot be calculated.

A fixed-power behaviour $`x^{ϵ_0}`$ of $`F_2(x,Q^2)`$, such as in (15), corresponds to an $`N`$-plane pole:

$$𝐮(N,Q^2)\frac{𝐟(Q^2)}{Nϵ_0}𝐟(Q^2)=\left(\begin{array}{c}f_0(Q^2)\\ f_g(Q^2)\end{array}\right)$$
(20)
If we insert this into the DGLAP equation (3) and equate the residue of the pole on each side of the equation, we find

$$\frac{}{t}𝐟(Q^2)=𝐏(N=ϵ_0,\alpha _s(Q^2))𝐟(Q^2)$$
(21)
$`ϵ_0`$ is far enough from 0 for the expansion of $`𝐏(N=ϵ_0,\alpha _s(Q^2))`$ to be reasonably safe. So we may easily use the DGLAP equation to calculate the evolution of the hard-pomeron component of $`F_2(x,Q^2)`$. But this is not the case for the soft-pomeron component, because $`ϵ_10.08`$ is too close to 0.

According to figure 25, the various gluon structure functions come together at $`x0.01`$. It is reasonable to assume that for values of $`x`$ larger than this the evolution of the two elements of $`𝐮(x,Q^2)`$ does not use values of $`N`$ close to 0 and therefore the conventional analysis is correct. So we can start at some not-too-large value of $`Q^2`$, 20 GeV2 say. We determine the value of $`f_0(Q^2)`$ there from the phenomenological fit (16) and $`f_g(Q^2)`$ from the MRST gluon structure function $`xg(x,Q^2)`$, which for $`x`$ greater than about 0.01 fits very well to $`x^{ϵ_0}(1x)^5`$. We choose $`\mathrm{\Lambda }`$ such that $`\alpha _s(M_X^2)=0.116`$ and use (21) to calculate the evolution of $`f_0(Q^2)`$ and $`f_g(Q^2)`$ in both directions. The result for $`f_0(Q^2)`$ is the continuous curve in figure 27. The dashed curve is the phenomenological form (16). Provided we adjust $`\mathrm{\Lambda }`$ so that still $`\alpha _s(M_X^2)=0.116`$, LO and NLO evolution give almost identical results.

The agreement between the pQCD calculation and the phenomenological curve is a success not only for the concept of the hard pomeron, but also for pQCD itself. The evolution is from a single value of $`Q^2`$, not the customary global fit, and it introduces far fewer parameters.

Notice that, as $`Q^2`$ increases, the large-$`x`$ behaviour of $`xg(x,Q^2)`$ becomes steadily steeper than $`(1x)^5`$, and so the largest value of $`x`$ for which $`x^{ϵ_0}`$ is a good approximation to the structure function steadily decreases. Figure 28 shows an estimate of this.

We may use the gluon structure function to calculate the charm structure function $`F_2^c(x,Q^2)`$. The result, using just LO photon-gluon fusion with a charm-quark mass $`m_c=1.3`$ GeV, is the solid curves in figure 27. This is an important check on the consistency of the approach. As is seen in figure 9, a steep gluon distribution is needed to fit the data at small $`Q^2`$.

In conclusion, the conventional approach to evolution needs modifying at small $`x`$. It can be corrected if we combine it with Regge theory, but only partly — we can only treat the hard-pomeron part. The resulting gluon distribution is larger at small $`x`$ than has so far been supposed and gives a good description of charm production. I should add that we want good data for the longitudinal structure function, because this gives the most direct window on the gluon distribution.

## 14 Summary

$``$What physics explains the dramatic HERAeffect? $``$Is unitarity a constraint on hard collisions? $``$Do $`pp`$ and/or $`\gamma p`$ total cross sections contain a hard term? $``$Why do we see no odderon at $`t=0`$? $``$How do we understand soft diffraction dissociation? $``$Is diffractive Higgs production large enough to measure? $``$Why does HERA see more hard diffractive events than the Tevatron? $``$The conventional approach to evolution needs modifying at small $`x`$ $``$It can be corrected if we combine it with Regge theory $``$But how do we handle the soft-pomeron part? $``$The gluon distribution is larger at small $`x`$ than has so far been supposed $``$We want good data for the longitudinal structure function


text:
# Lattice QCD with light dynamical quarksCombining the contributions of F. Farchioni and L. Scorzato.

## 1 INTRODUCTION

Predicting the low energy constants of the chiral effective Lagrangian from first principles in QCD is an exciting challenge for lattice gauge theories. This question has attracted a considerable attention in the last few years . However, existing unquenched simulations are typically done in a region where the light quarks are not light enough, in most cases – especially with Wilson-type quarks – with two light quark flavours ($`u`$ and $`d`$) having masses larger than half of the strange quark mass ($`m_{ud}>\frac{1}{2}m_s`$). Since the quark mass dependence in the chiral Lagrangian is explicit, it is not necessary to perform simulations at the physical value of the $`u`$\- and $`d`$-quark masses, which remains a too hard task for the forseeable future. Nevertheless, in order to keep the systematic errors under control, three dynamical quarks should be simulated and, according to , the masses of the light ones should satisfy at least $`m_{ud}\frac{1}{4}m_s`$.

Simulating light dynamical quarks is a difficult task for numerical computations because all known algorithms for QCD have a substantial slowing down towards small quark masses. The present status has been summarized at the Berlin lattice conference . Inspired by that discussion, in a recent paper  we parametrized the cost of an unquenched QCD simulation in two different ways:

$`C=F(r_0m_\pi )^{z_\pi }\left({\displaystyle \frac{L}{a}}\right)^{z_L}\left({\displaystyle \frac{r_0}{a}}\right)^{z_a},`$ (1)
$`C_U=F_U\left({\displaystyle \frac{m_\pi }{m_\rho }}\right)^{z_{\pi \rho }}\left({\displaystyle \frac{L}{a}}\right)^{z_L}\left({\displaystyle \frac{r_0}{a}}\right)^{z_a}.`$ (2)
Here $`r_0`$ is the Sommer scale parameter , $`m_\pi `$ the pion mass, $`L`$ the lattice extension and $`a`$ the lattice spacing. The powers $`z_{\pi ,\pi \rho ,L,a}`$ and the overall constants $`F`$, $`F_U`$ are empirically determined. The value of the constants $`F`$, $`F_U`$ depends on the definition of “cost”. Of course, when the quark masses are very light and the decay $`\rho \pi \pi `$ is possible, the first parametrization is preferable. Typical values of the parameters can be taken from Ukawa’s contribution : $`F_U=5.910^6`$ flop, $`z_{\pi \rho }=6`$, $`z_L=5`$, $`z_a=2`$. Other estimates are in reasonable agreement. However most of the simulations, whereon such estimates are based, were perfomed at quite large quark masses. This is a reason of concern because, as experienced already by several collaborations , at small quark masses new phenomena appear, especially for algorithms based on molecular dynamics integrators.

Given the compelling need of simulations at lighter quarks, we performed algorithmic studies in the regime of very light quark masses (though on very rough lattices)  using the two-step multi-boson algorithm . Here we give an update of the cost analysis. We also include a preliminary study of the volume dependence and a test of a possible improvement of the algorithm.

## 2 COST OF THE SIMULATIONS

### 2.1 Algorithm

Let us first briefly summarize the main features of TSMB. The quark determinant of $`N_f`$ degenerate flavours is represented as

$$\left|det(\stackrel{~}{Q})\right|^{N_f}\frac{1}{detP_{n_1}^{(1)}(\stackrel{~}{Q}^2)detP_{n_2}^{(2)}(\stackrel{~}{Q}^2)}.$$
(3)
$`\stackrel{~}{Q}`$ is the hermitean Wilson-Dirac fermion matrix. The polynomials $`P^{(1)}`$ and $`P^{(2)}`$ satisfy

$$\underset{n_2\mathrm{}}{lim}P_{n_1}^{(1)}(x)P_{n_2}^{(2)}(x)=x^{N_f/2},x[ϵ,\lambda ].$$
(4)
The interval $`[ϵ,\lambda ]`$ covers the spectrum of $`\stackrel{~}{Q}^2`$ on a typical gauge configuration. The first polynomial $`P^{(1)}`$ of order $`n_1`$ is a crude approximation and is realized by the multi-boson representation. The second polynomial $`P^{(2)}`$ of order $`n_2n_1`$ gives a better approximation. It is taken into account in the updates by a global noisy Metropolis correction step. Since for fixed $`n_2`$ (and outside the interval $`[ϵ,\lambda ]`$) the approximation $`P^{(1)}P^{(2)}`$ is not exact, a final correction step is performed by reweighting the gauge configurations which are considered for the evaluation of the expectation values. This is done with a sufficiently high order polynomial $`P_{n_4}^{(2)}`$. For more details about the algorithm in this context we refer to .

### 2.2 Parameters and the definition of cost

The Monte Carlo simulations have been done near the $`N_t=4`$ thermodynamical cross-over line. We tuned the gauge coupling $`\beta `$ and the hopping parameter $`\kappa `$ in order to explore a range of the quark mass going from $`2m_s`$ to $`\frac{1}{5}m_s`$ while keeping $`r_0/a`$ roughly constant between $`r_0/a=1.7`$ and $`r_0/a=1.9`$, which implies $`a0.27`$ fm. Most of the analysis is performed on lattices of size $`8^316`$ implying a physical lattice extension $`L2.2\mathrm{fm}`$. For the study of finite volume effects we used $`12^324`$ and $`16^4`$ lattices.

The cost of a simulation can be measured in units of the floating point operations (flop) which are necessary to produce a new decorrelated configuration (starting from a themalized configuration). The decorrelation quantified by the integrated autocorrelation length $`\tau _{int}`$ depends on the quantity whose autocorrelation is observed. We consider the average plaquette, the lowest eigenvalue of the fermion matrix, the pion correlator at some fixed distance and the pion mass. Instead of the number of flops the integrated autocorrelations can also be quoted by the necessary number of fermion-matrix vector multiplications (MVM’s). We always count mutiplications by the fermion matrix $`\stackrel{~}{Q}`$. In this way the dependence on the computer architecture is reduced. For TSMB we have the following approximate formula for the total amount of MVM’s needed for one update cycle:

$$N\mathrm{\hspace{0.33em}6}(n_1N_\mathrm{\Phi }+N_U)+2(n_2+n_3)N_C+I_GF_G.$$
(5)
Here $`N_\mathrm{\Phi }`$ is the number of local bosonic sweeps per update cycle, $`N_U`$ the number of local gauge sweeps, $`N_C`$ the number global Metropolis accept-reject correction steps, and $`I_G`$ and $`F_G`$ give the number of MVM’s and frequency of the global quasi heatbath. ($`n_3`$ is the order of an auxiliary polynomial $`P^{(3)}`$.) Of course, the number of MVM’s can easily be converted to flops by noting that the number of flops per lattice point (in our code) is 1344. The quantity in (5) is reported for the runs considered here in the columns 7 and 5 of tables 1 and 2 respectively.

### 2.3 Results about the costs

In  power fits have been given for the dependence of the integrated autocorrelations on the quark mass parameter $`M_r(r_0m_\pi )^2`$. For instance, in case of the average plaquette a good fit is

$`\tau _{int}^{plaq}=c_{plaq}M_r^{z_{plaq}},`$
$`c_{plaq}=7.92(68),z_{plaq}=2.02(10),`$ (6)
if the cost is counted in MMVM’s ($`=10^6`$ fermion-matrix vector multiplications). This corresponds to $`z_\pi 4`$ in (1). A reasonable fit of $`\tau _{int}^{plaq}`$ can also be found in the form (2) with $`z_{\pi \rho }6`$ if the heaviest quark mass point is omitted.

Another interesting quantity is the pion correlator and/or the pion mass. For instance, the cost for the pion correlator at timeslice distance $`d=5`$ is shown together with two power fits in figure 1. The fit to all points gives a power similar to $`z_{plaq}`$:

$$c_5=2.68(35),z_5=1.96(17).$$
(7)
However, if the point at $`M_r3`$ is omitted from the fit we obtain a lower power:

$$c_5=2.99(22),z_5=1.15(28).$$
(8)
Fitting the cost for obtaining the pion mass, the result is:

$$c_{m_\pi }=1.99(16),z_{m_\pi }=1.47(16).$$
(9)
(The data and the fit is shown in figure 2). This corresponds to $`z_\pi 3`$ in (1), similarly to the result for the minimal eigenvalue of the squared hermitean fermion matrix:

$$c_{min}=5.36(80),z_{min}=1.48(25).$$
(10)
It is remarkable that the pion mass has a substantially lower cost than, for instance, the average plaquette. In fact, at $`M_r0.5`$ the cost ratio is already almost 10 in favour of the pion mass. This is partly due to the intrinsic fluctuation present in the pion propagator which is originated from the freedom of randomly choosing the position of the source.

### 2.4 Analysis with correction factors

The importance of the final correction step by the reweighting factors depends on the choice of the second polynomial $`P^{(2)}`$. For good enough $`P^{(2)}`$ the reweighting is negligible. In  this has been the case for all but one runs. (The run where the reweighting was not completely negligible has the index (i).) Of course, at very small quark masses $`P_{n_2}^{(2)}`$ becomes good enough only for a very large order $`n_2`$. In such a case it is a question whether it would be more advantageous to keep $`n_2`$ smaller and tolerate a non-trivial reweighting.

The consequences of a poor $`P^{(2)}`$ at small quark masses could, however, be rather unpleasant. An example has been shown in : in run $`(i)`$ at $`M_r0.6`$ one of the four parallel runs had a serious fluctuation with exceptionally small eigenvalues and, as a consequence, very small correction factors in a substantial part of the run (see figure 2 in ). In this part of the run also a few negative eigenvalues of the fermion matrix have been observed. The smallness of the corresponding reweighting factors implies that with a better second polynomial $`P^{(2)}`$ such fluctuations could have been suppressed.

The question is, what is the consequence of such a fluctuation on the integrated autocorrelations and on the related magnitude of statistical errors. This can be investigated, for instance, by considering the linearized deviation of effective pion masses which appears in the linearization method for determining the autocorrelations of functions of primary averages . For an illustration see the case of the effective mass between timeslice distances $`d=3`$ and $`d=7`$ after including the reweighting factors in the above mentionned run of  (figure 3). As this figure shows, the exceptionally small eigenvalues and the corresponding small reweighting factors are accompanied by a large downwards fluctuation of the effective mass $`m_\pi ^{eff}`$. This implies that in this run the average value of $`m_\pi ^{eff}`$ is smaller than in the “normal” runs and at the same time its error is larger. This pulls down the overall mass estimate: the average of all four parallel runs gives for this effective mass $`m_\pi ^{eff}=0.4160(103)`$. For the three other runs the result is: $`m_\pi ^{eff}=0.4298(63)`$. Our conclusion is that it is probably better to suppress the kind of fluctuations in figure 3 by sufficiently improving $`P^{(2)}`$.

## 3 VOLUME DEPENDENCE

The physical volume of the lattices in  is such that strong finite volume effects are not expected. For checking this and for observing the behaviour of our algorithm in increasing volumes new simulations were performed on a $`12^324`$ and a $`16^4`$ lattices (see (f12) and (e16) in table 1). In order to hold all the quantities fixed, except for the volume, run (f12) has the same value of $`\beta `$ and $`\kappa `$ as run (f) and run (e16) the same as run (e). The lattice extensions deduced from $`r_0/a`$ are collected in the last column of table 1.

When going to a larger lattice one has to adjust in TSMB the first polynomial order $`n_1`$. Assuming no large finite volume effects the interval given by $`ϵ`$ and $`\lambda `$ remains unchanged. Therefore it is not needed to change $`n_2`$ and $`n_3`$. This is true for the runs (e) and (e16). For run (f) the lower boundary $`ϵ`$ was not optimal therefore we used a better (i.e. smaller) one in run (f12). Consequently, $`n_2`$ and $`n_3`$ had to be increased there.

A remarkable feature in table 1 is that on the larger lattices the autocorrelation of the pion mass becomes even more favourable compared to the average plaquette. This has the consequence that the cost for an integrated autocorrelation length of the pion mass increases substantially slower than the number of lattice points. As expected only a tiny finite-size effect is observed for the physical quantities measured.

## 4 DETERMINANT BREAKUP

Inspired by  (see also ) we tried to improve autocorrelations by “determinant breakup”. In our case this means, for instance, that simulating two flavours separately should be more effective than two degenerate flavours together. In general, one can write the determinant as

$$\left|det(\stackrel{~}{Q})\right|^{N_f}=\left\{\left|det(\stackrel{~}{Q})\right|^{N_f/N_b}\right\}^{N_b},$$
(11)
with an arbitrary positive integer $`N_b`$ and use the polynomial approximation (3) for the individual factors on the right hand side.

For tests we have chosen run (b) in  where $`M_r3`$ (close to the strange quark mass). Results are shown in table 2. (Note that $`\tau _{int}^{m_\pi }`$ are obtained here from effective masses with the linearization method for autocorrelations and errors and therefore are slightly different from .) Tests at smaller quark masses and/or in larger volumes will be done later.

Choosing some breakup of the flavours we fixed $`n_1`$ by requiring a constant acceptance rate of $`50`$-$`60\%`$. The higher polynomial orders $`n_2`$ and $`n_3`$ were fixed by keeping the relative deviation at the interval ends roughly constant. We found that the autocorrelation measured in update cycles is reduced when breaking up the determinant in more pieces. However the costs per update cycles are increasing. It is remarkable that the cost is rising less than one would expect by just looking at the numbers for $`n_1`$ and $`n_2`$. The reason is that the total number of iterations in the global quasi heatbath $`I_G`$ in (5) are roughly constant, although the number of boson fields $`N_bn_1`$ is increasing. Putting everything together we find that at this simulation point one can gain from determinant breakup almost a factor of two when splitting up two flavours in two pieces.

## 5 EIGENVALUES

The study of the low lying modes of the Wilson-Dirac operator ($`Q`$) is interesting both physically and algorithmically. It is an interesting question how the qualitative features of the spectrum change with increasing volume. As described in , we used the implicitely restarted Arnoldi method  to compute a portion of the eigenvalues of the non-hermitean Wilson-Dirac fermion matrix $`Q`$. Typically we computed the eigenvalues with lower real part of the matrix $`Q`$ and those with smallest modulus of the preconditioned matrix $`\overline{Q}`$ (as defined in ). By use of the analytical relations in  we had access to a reasonable part of the spectrum of $`Q`$ around the origin.

In figure 4 we plot $`O(200)`$ eigenvalues from 10 configurations thermalized at point (f) on the $`8^316`$ lattice. The portion of spectrum to which we have access is that inside the dashed circle and to the left of the dashed vertical line. In figure 5 we plot again $`O(200)`$ eigenvalues from 10 configurations thermalized at point (f12) on $`12^324`$ lattice, which has the same parameters as (f) – except for the volume. As expected, the lowest border of the spectrum does not essentially move, even if the density of eigenvalues is higher (and as a consequence, the portion of spectrum that we can cover is smaller).

## 6 CONCLUSION

In conclusion, the cost increase with the lattice volume is quite acceptable because it is close to the trivial volume factor. In case of the autocorrelation of the pion mass the observed increase turns out to be even smaller.

The determinant breakup is an easy and effective method to speed up the gauge field update.