text:
# Some exact results for Boltzmann’s annihilation dynamics

## I introduction

In ballistically controlled reactions, particles with a given initial velocity distribution move freely (ballistic motion) in a $`d`$-dimensional space. When two of them meet, they annihilate and disappear from the system. This apparently simple problem has attracted a lot of attention during the past years EF ; red1 ; red2 ; krap ; pia\_one ; us\_one ; us\_due ; rey\_one ; kafri ; evans ; sire ; chop ; triz ; us\_three  for the following reasons. First, this is one of the few problems of nonequilibrium statistical physics that can be exactly solved in some cases, and second it models some growth and coarsening processes krug\_one .

This field was entered with the pioneering work by Elskens and Frisch EF , where a one-dimensional system with only two possible velocities $`+c`$ or $`c`$ was studied. Using combinatorial analysis, they showed that the density of particles was decreasing according to a power law ($`t^{1/2}`$) in the case of a symmetric initial velocity distribution. The investigation of this one-dimensional problem was generalized by Droz et alus\_one  to the three-velocity case where the initial velocity distribution is given by $`\phi (v;0)=p_+\delta (vc)+p_0\delta (v)+p_{}\delta (v+c)`$ with $`p_+=p_{}`$ (symmetric case) and $`p_++p_0+p_{}=1`$. It turns out that the decay of the particle density depends on the details of the initial velocity distribution. The following analytical results were obtained. For $`p_0<1/4`$, the density $`n(v;t)`$ of particles with velocity $`v=\{0,+c,c\}`$, behaves in the long-time limit as $`n(0;t)t^1,n(\pm c;t)t^{1/2}`$. When $`p_0=1/4`$, $`n(0;t)n(\pm c;t)t^{2/3}`$. Finally, for $`p_0>1/4`$, one finds that $`n(0;t)`$ saturates to a nonzero stationary value, while $`n(\pm c;t)`$ decays faster than a power law. Moreover, it was shown that in one dimension, annihilation dynamics creates strong correlations between the velocities of colliding particles, which excludes a Boltzmann-like approximation. Pairs of nearest neighbor particles have the tendency to align their velocities and propagate in the same direction us\_one .

An analytical investigation of the one-dimensional case with a continuous velocity distribution is much more difficult. A dynamical scaling theory, whose validity was supported by extensive numerical simulations for several velocity distributions, led Rey et alus\_due  to the conjecture that all the continuous velocity distributions $`\phi (v)`$ that are symmetric, regular and, such that $`\phi (0)0`$ are attracted in the long-time regime towards the same Gaussian-like distribution and thus belong to the same universality class.

For higher dimensions, most of the studies are based on an uncontrolled Boltzmann-like description red1 ; red2 ; sire  or numerical simulations chop . However, based on phenomenological mean-field-like arguments, Krapivsky et al. have studied the annihilation kinematics of a bimodal velocity modulus distribution in $`d>2`$ dimensions krap . In the case of a mixture of moving and motionless particles they showed that the stationary particles always persist, while the density of moving particles decays exponentially. This approach contains unknown phenomenological parameters, and thus a complete comparison with the results obtained by numerical simulation is not possible.

In a recent paper, Piasecki et alus\_three  gave an analytical derivation of the hierarchy equations obeyed by the reduced distributions for the annihilation dynamics. In dimension $`d>1`$ for a spatially homogeneous system, and in the limit (the so-called Grad limit) for which the particle diameter $`\sigma 0`$, and the particle density $`n(t)\mathrm{}`$ such that $`n(t)\sigma ^{d1}=\lambda ^1`$, where $`\lambda `$ is the mean free path, the hierarchy reduces to the Boltzmann-like hierarchy. This hierarchy propagates the factorization of the reduced $`k`$-particle distribution in terms of one-particle distribution functions. Thus, if the initial state is factorized, the whole hierarchy reduces to one nonlinear equation for the one-particle Boltzmann distribution. For annihilation kinetics, the ratio of particle diameter to mean free path vanishes in the long-time limit and the situation becomes similar to the Grad limit discussed above for $`\lambda \mathrm{}`$. Thus the long-time limit of the annihilation dynamics (for $`d>1`$) is likely to be adequately described by the nonlinear Boltzmann equation.

A scaling analysis of the nonlinear Boltzmann equation led to analytical expressions for the exponents describing the decay of the particle density and of the root-mean-square velocity in the case of continuous velocity distributions us\_three .

In view of the different behaviors observed in one dimension for discrete or continuous velocity distributions, it is relevant to study the case of distributions with discrete modulus spectrum in dimensions higher than $`1`$. The goal of this paper is to investigate simple examples of this kind in three dimensions for which the non linear Boltzmann equation derived in Ref. us\_three  can be exactly solved. The generalization of this approach to an arbitrary dimension is straightforward, and for the sake of comparison with numerical simulations, we shall also consider the two-dimensional situation in some detail.

The validity of the Boltzmann description in the long-time limit will be confirmed by comparing our analytical predictions with the results obtained by a molecular dynamics simulation.

The paper is organized as follows. In Sect. II we define the model. In Sect. III the three-dimensional Boltzmann equation is solved analytically for a two velocity modulus ($`c_1`$ and $`c_2`$) isotropic distribution. For simplicity we first consider the one velocity model $`c_1=c_2>0`$ that allows to draw interesting comparisons with the same model in one dimension. Then the implicit solution for the particle densities in the general case $`c_1>c_2>0`$ is established. It is shown analytically that in the long-time limit the particle densities decay according to power laws, with exponents depending continuously on the value of the velocity modulus ratio. We also find upper and lower bounds to the particle densities that are compared with the numerical solution of the dynamical equation. The particular case of a mixture of moving ($`c_2>0`$) and motionless ($`c_1=0`$) particles is also investigated. It turns out that the particle densities decay exponentially to zero for the moving particles, and to a nonzero value for the motionless ones. This phenomenology is independent of space dimension, and in Sect. IV, it will be shown explicitly to hold in two dimensions by implementing molecular dynamics simulations. This numerical method has the advantage of being free of the approximations underlying Boltzmann’s dynamics and, therefore, provides an interesting test for the analytical predictions. Section V contains our interpretations and conclusions.

## II The model

We consider a system made of spheres of diameter $`\sigma `$ moving ballistically in three-dimensional space. If two particles touch each other, they annihilate and thus disappear from the system. We consider only two-body collisions. The initial spatial distribution of particles is supposed to be uniform, therefore it remains uniform during the evolution. Indeed, if the state is initially translationally invariant, then the free evolution preserves this symmetry. Annihilation dynamics adds the effect of binary collisions that depends only on the distance between particles, thus preserving the translational invariance. Finally, existing numerical simulations seem to be compatible with this assumption of homogeneity us\_three . We are interested in the time evolution of the number density of particles with a given velocity modulus.

Let $`f_1(𝐯;t)`$ be the distribution function of the density of particles in $`^3`$ with velocity $`𝐯^3`$ at time $`t`$. For spatially homogeneous states, the distribution function has the form

$$f_1(𝐯;t)=n(t)\phi (𝐯;t),$$
(1)
where $`\phi (𝐯;t)`$ is the velocity probability density. In the long-time limit, Piasecki et alus\_three  have shown that the hierarchy satisfied by the reduced distributions approached the Boltzmann hierarchy. If the initial state is factorized, the nonlinear Boltzmann equation provides then the complete description of annihilation dynamics,

$$\frac{}{t}f_1(𝐯_1;t)=\sigma ^2𝑑\widehat{𝝈}\theta (\widehat{𝝈}\widehat{𝐯}_{12})(\widehat{𝝈}\widehat{𝐯}_{12})_^3𝑑𝐯_2|𝐯_{12}|f_1(𝐯_1;t)f_1(𝐯_2;t).$$
(2)
Here $`\theta `$ is the Heaviside function, $`𝐯_{12}=𝐯_1𝐯_2`$ the relative velocity of two particles, $`\widehat{𝐯}_{12}=𝐯_{12}/v_{12}`$ a unit vector, $`v_{12}=|𝐯_{12}|`$, and the integration with respect to $`d\widehat{𝝈}`$ is the angular integration over the solid angle.

We consider spherically symmetric initial conditions $`f_1(v;0)`$, $`v=|𝐯|`$. This symmetry property is propagated by the dynamics. The Boltzmann equation (2) then takes the form

$$\frac{}{t}f_1(v;t)=\frac{2}{3}(\pi \sigma )^2f_1(v;t)_0^{\mathrm{}}𝑑uu^2f_1(u;t)\left[\frac{(u+v)^3|uv|^3}{uv}\right].$$
(3)
Equation (3) is a nonlinear homogeneous integral equation for the distribution function $`f_1(v;t)`$. A simplification arises if the initial velocity distribution has a discrete modulus spectrum. This spectrum is preserved by the annihilation dynamics as no new velocities are created. A simple case is provided by the bimodal distribution

$$\phi (v,0)=\frac{A}{4\pi c_1^2}\delta (vc_1)+\frac{1A}{4\pi c_2^2}\delta (vc_2),$$
(4)
where $`c_2>c_10`$ and $`A`$ denotes the fraction of particles with velocity modulus $`c_1`$.

## III Exact results

Before addressing the general case, we first consider the single-species problem where $`c_2=c_1>0`$.

### III.1 Single-velocity modulus distribution

Setting $`c_2=c_1=c>0`$ in Eq. (4), one obtains from Eq. (1)

$$f_1(v;t)=n(t)\frac{1}{4\pi c^2}\delta (vc).$$
(5)
From the kinetic equation (3), we find

$$\frac{d}{dt}n(t)=\frac{4}{3}\pi \sigma ^2cn(t)^2,$$
(6)
whose solution is

$$n(t)=\frac{n_0}{1+\frac{4}{3}\pi \sigma ^2n_0ct},$$
(7)
where $`n(0)=n_0`$. A striking observation is that, in the limit $`t\mathrm{}`$, the density (7) becomes independent of its initial value $`n_0`$. Note that the same phenomenon is also present for simple diffusion limited annihilation such as $`A+A0`$, when the dimension of the system is larger than $`2`$ reactiondiffusion .

Contrary to the one-dimensional case for which it has been rigorously shown that the density decays proportionally to $`t^{1/2}`$ EF , one sees from Eq. (7) that in three dimensions, Boltzmann’s dynamics is faster as the density decays according to $`t^1`$, which is the mean-field value krap . We note, however, that the same behavior $`n(t)1/t`$ holds in all dimensions within Boltzmann’s kinetic theory (and in fact, more generally within the framework of a scaling analysis of the hierarchy governing the dynamics of ballistic annihilation us\_three ). This discrepancy between Boltmann’s prediction and the exact result in one dimension illustrates the crucial importance of dynamical correlations when $`d=1`$. On the other hand, as suggested in Ref. us\_three  and explicitly shown below by molecular dynamics simulations, the nonlinear Boltzmann equation is relevant for describing the long-time dynamics of ballistic annihilation when $`d2`$. In this case the particles are very diluted and no dynamical correlations can develop during the time evolution, which would violate the molecular chaos hypothesis.

### III.2 Mixture of particles with two nonzero velocity moduli

Consider the case where particles with velocities $`c_1>0`$ and $`c_2>c_1`$ are initially present. Thus $`f_1(v;t)`$ is of the form

$$f_1(v;t)=X(t)\frac{1}{4\pi c_1^2}\delta (vc_1)+Y(t)\frac{1}{4\pi c_2^2}\delta (vc_2),$$
(8)
where $`X(t)`$ and $`Y(t)`$ are, respectively, the densities of particles with velocities $`c_1`$ and $`c_2`$. They add up to the total density $`X(t)+Y(t)=n(t)`$. Upon rescaling the time according to $`\tau =tc_2\pi \sigma ^2/3`$, it follows from Eq. (3) that

$`\dot{X}(\tau )`$ $`=`$ $`4\gamma X(\tau )^2(3+\gamma ^2)X(\tau )Y(\tau ),`$ (9a)
$`\dot{Y}(\tau )`$ $`=`$ $`4Y(\tau )^2(3+\gamma ^2)X(\tau )Y(\tau ),`$ (9b)
where $`0\gamma =c_1/c_2<1`$, and the overdot denotes time derivative with respect to $`\tau `$.

The set of equations (9) is a nonlinear homogeneous system of coupled differential equations with constant coefficients. An implicit solution can be obtained by introducing the function $`V(\tau )`$ defined as $`V(\tau )=Y(\tau )/X(\tau )`$. From Eq. (9) we get

$$\frac{dY}{dX}=\frac{4V^2+(3+\gamma ^2)V}{4\gamma +(3+\gamma ^2)V},$$
(10)
so that

$$\frac{dX}{X}=\frac{4\gamma +(3+\gamma ^2)V}{(1\gamma ^2)V^2+(1\gamma )(3\gamma )V}dV.$$
(11)
Integrating Eq. (11) yields

$$\frac{X_0}{X}=\left(\frac{V_0}{V}\right)^\alpha \left(\frac{V_0+\frac{3\gamma }{1+\gamma }}{V+\frac{3\gamma }{1+\gamma }}\right)^\beta ,$$
(12)
with $`\alpha =4\gamma /[(1\gamma )(3\gamma )]0`$, $`\beta =(3+\gamma ^2)/(1\gamma ^2)\alpha >0`$, $`V(0)=V_0=Y_0/X_0`$, $`X_0=X(0)`$, $`Y_0=Y(0)`$. The special case of $`\gamma =0`$ will be discussed in Sec. III.3, hence from now on we assume that $`\gamma >0`$, so that $`\alpha >0`$. Equations (9) can also be written as

$`{\displaystyle \frac{d}{d\tau }}\left({\displaystyle \frac{1}{X}}\right)`$ $`=`$ $`4\gamma +(3+\gamma ^2){\displaystyle \frac{Y(\tau )}{X(\tau )}},`$ (13a)
$`{\displaystyle \frac{d}{d\tau }}\left({\displaystyle \frac{1}{Y}}\right)`$ $`=`$ $`4+(3+\gamma ^2){\displaystyle \frac{X(\tau )}{Y(\tau )}}.`$ (13b)
Multiplying the right-hand side (RHS) of Eq. (13a) by $`X_0`$ and equating it with the derivative of the RHS of Eq. (12), one obtains upon integration from $`0`$ to $`\tau `$ the relation

$$X_0\tau =_V^{V_0}𝑑u\frac{1}{4\gamma +(3+\gamma ^2)u}\left\{\frac{d}{du}\left[\left(\frac{V_0}{u}\right)^\alpha \left(\frac{V_0+\frac{3\gamma }{1+\gamma }}{u+\frac{3\gamma }{1+\gamma }}\right)^\beta \right]\right\}.$$
(14)
Equation (14) implicitly defines the time dependence of the function $`V(\tau )`$. The procedure to obtain the densities $`X(\tau )`$ and $`Y(\tau )`$ from Eq. (14) is as follows. The integration in Eq. (14) leads to Appel functions, that may be inverted (at least numerically) in order to give $`V(\tau )`$. The insertion of $`V(\tau )`$ in Eq. (12) then gives $`X(\tau )`$. It is then straightforward to obtain $`Y(\tau )`$, having determined $`V(\tau )`$ and $`X(\tau )`$. The structure of the implicit relation (14) permits us to establish interesting analytical results.

First, let us investigate the long time behavior of the particle densities $`X(\tau )`$ and $`Y(\tau )`$. When $`\tau \mathrm{}`$, the LHS of Eq. (14) diverges linearly which implies that $`lim_\tau \mathrm{}V(\tau )=0`$. So, in the long-time limit, the implicit relation (14) leads to the asymptotic formula

$$X_0\tau \stackrel{\tau \mathrm{}}{}\frac{1}{4\gamma }\left(\frac{V_0}{V}\right)^\alpha \left(\frac{1+\gamma }{3\gamma }V_0+1\right)^\beta ,$$
(15)
so that

$$V(\tau )\stackrel{\tau \mathrm{}}{}V_0\left(\frac{1+\gamma }{3\gamma }V_0+1\right)^{\beta /\alpha }\left(4\gamma X_0\tau \right)^{1/\alpha }.$$
(16)
As $`lim_{\tau 0}V(\tau )=0`$, Eq. (13a) takes the asymptotic form $`d/d\tau (1/X)=4\gamma `$. Hence we conclude that

$$X(\tau )\stackrel{\tau \mathrm{}}{}\frac{1}{4\gamma }\tau ^1.$$
(17)
On the other hand, from Eqs. (13b) and (16) we find the long-time relation

$$\frac{d}{d\tau }\left(\frac{1}{Y}\right)\stackrel{\tau \mathrm{}}{}(3+\gamma ^2)\frac{1}{V_0}\left(\frac{1+\gamma }{3\gamma }V_0+1\right)^{\beta /\alpha }\left(4\gamma X_0\tau \right)^{1/\alpha },$$
(18)
which upon integration yields

$$Y(\tau )\stackrel{\tau \mathrm{}}{}\frac{V_0}{4\gamma }(4\gamma X_0)^{1/\alpha }\left(\frac{1+\gamma }{3\gamma }V_0+1\right)^{\beta /\alpha }\tau ^{(3+\gamma ^2)/4\gamma }.$$
(19)
Note that the exponent for the density $`Y(\tau )`$ is a function of the ratio $`\gamma =c_1/c_2`$ and thus is nonuniversal. In the limit $`\gamma 1`$ one recovers the asymptotic behavior of the single-velocity modulus distribution (see Sect. III.1).

Second, we may find analytical upper and lower bounds for $`X(\tau )`$ and $`Y(\tau )`$. Granted that $`\alpha >0`$ and $`\beta >0`$, the integrand of Eq. (14) is a strictly monotonic decreasing positive function of $`u`$, therefore $`V(\tau )<V_0`$ for all $`\tau >0`$. Considering that $`(4\gamma )^1[4\gamma +(3+\gamma ^2)u]^1`$ for $`u0`$, the insertion of Eq. (12) in Eq. (14) provides the inequality $`X_0\tau \left(X_0/X1\right)/4\gamma `$, which leads to an upper bound for $`X(\tau )`$. On the other hand, the inequality $`(4\gamma +(3+\gamma ^2)u)^1[4\gamma +(3+\gamma ^2)V_0]^1`$ yields a lower bound, so that we finally get

$$\frac{X_0}{1+[4\gamma X_0+(3+\gamma ^2)Y_0]\tau }X(\tau )\frac{X_0}{1+4\gamma X_0\tau }.$$
(20)
Note that for times such that

$$4\gamma X_0\tau 1,$$
(21)
the upper bound (20) coincides with the exact asymptotic relation (17). The same kind of analysis as that leading to Eq. (20) yields the upper bound,

$$0Y(\tau )\frac{Y_0}{1+(3+\gamma ^2)X_0\tau }.$$
(22)
The width defined by the difference of the bounds in both cases (20) and (22) is $`𝒪\left(\tau ^1\right)`$. Figures 1 and 2 show the numerical solution for $`X(\tau )`$, $`Y(\tau )`$, the bounds (20) and (22), as well as their asymptotic behaviors (17) and (19) on a logarithmic scale.

The knowledge of the numerical solution (see Figs. 1 and 2) allows to determine the crossover time, separating the early and long-time (power law) regimes.

### III.3 Mixture of moving and motionless particles

We now consider a particular case of Sec. III.2 that we solve exactly in the asymptotic limit $`\tau \mathrm{}`$. The system is now characterized by a certain number of motionless particles (zero velocity, $`c_1=0`$) whereas the rest of the particles have a given nonzero velocity modulus. Thus, setting $`\gamma =0`$ in Eq. (12), inverting the relation in order to find $`V=V(X)`$, then making use of Eq. (13a) with $`\gamma =0`$ leads to

$$\frac{d}{d\tau }\left(\frac{1}{X}\right)=3\left(3+\frac{Y_0}{X_0}\right)\left(\frac{X}{X_0}\right)^{1/3}9.$$
(23)
The integration of Eq. (23) yields

$$3X_0\tau =\frac{1}{3}\left[1\frac{X_0}{X}\right]+\frac{a}{2}\left[1\left(\frac{X_0}{X}\right)^{2/3}\right]+a^2\left[1\left(\frac{X_0}{X}\right)^{1/3}\right]+a^3\mathrm{ln}\left[\frac{a1}{a\left(X_0/X\right)^{1/3}}\right],$$
(24)
with $`a=1+Y_0/3X_0`$. In the asymptotic limit $`\tau \mathrm{}`$, the LHS of Eq. (24) tends to $`+\mathrm{}`$. The density $`X(\tau )`$ cannot thus tend to zero, and must approach a strictly positive value $`X(\mathrm{})=X_{\mathrm{}}>0`$. For $`\tau \mathrm{}`$, all terms on the RHS of Eq. (24) but the logarithmic one approach a finite limit. This implies the asymptotic behavior

$$X_{\mathrm{}}=\underset{\tau \mathrm{}}{lim}X(\tau )=\frac{X_0}{a^3}=X_0\left(\frac{3X_0}{3X_0+Y_0}\right)^3>0.$$
(25)
Replacing $`X`$ by $`X_{\mathrm{}}`$ in all terms of Eq. (24) except the logarithmic one and then inverting the relation $`\tau (X)`$, we find

$$X(\tau )\stackrel{\tau \mathrm{}}{}X_{\mathrm{}}\left[1\epsilon (X_0,Y_0;\tau )\right]^3\stackrel{\tau \mathrm{}}{}X_{\mathrm{}}[1+3\epsilon (X_0,Y_0;\tau )],$$
(26)
where

$$\epsilon (X_0,Y_0;\tau )=\frac{Y_0\mathrm{exp}(\mu /a^3)}{3X_0+Y_0}\mathrm{exp}(3X_{\mathrm{}}\tau ),$$
(27)
and $`\mu =1/3+a/2+a^211a^3/16<0`$. Making use of Eq. (9a) with $`\gamma =0`$, we find

$$Y(\tau )\stackrel{\tau \mathrm{}}{}3X_{\mathrm{}}\epsilon (X_0,Y_0;\tau ).$$
(28)
Hence we have

$$X(\tau )\stackrel{\tau \mathrm{}}{}X_{\mathrm{}}+Y(\tau ).$$
(29)
There is a qualitative difference from the case $`c_1>0`$. As shown in Fig. 3, the density of particles at rest approaches the asymptotic value $`X_{\mathrm{}}>0`$ exponentially fast, while the density of moving particles goes to zero exponentially. Table 1 summarizes the long-time behavior for the different cases.

Note that generalizing our results to any dimension $`d2`$ is straightforward (see below for the case $`d=2`$). The algebraic or exponential decay of the particle densities hold irrespective of $`d`$. In particular, for the general case $`c_1>0`$ the exponent of the density of “slow” particles is independent of $`d`$ so that $`X(\tau )\stackrel{\tau \mathrm{}}{}\tau ^1`$. Finally the relation (29) still holds.

## IV Comparison with numerical simulations

The analytical predictions obtained in the preceding section rely on the validity of the molecular chaos assumption, leading to the Boltzmann equation. It is therefore instructive to compare these predictions to the results of molecular dynamics (MD) simulations, where the exact equations of motion of the particles are integrated (see Ref. us\_three  for more details concerning the method).

### IV.1 Analytical results in two dimensions

MD simulations are most efficiently performed in two dimensions, where the best statistical accuracy can be achieved. We consequently repeat the analysis of Secs. II and III for a two-dimensional system. Introducing the rescaled time $`\tau =2\pi \sigma c_2t`$, one obtains the counterpart of Eqs. (9) in the form

$`\dot{X}(\tau )`$ $`=`$ $`4\gamma X(\tau )^2\kappa (\gamma )X(\tau )Y(\tau ),`$ (30a)
$`\dot{Y}(\tau )`$ $`=`$ $`4Y(\tau )^2\kappa (\gamma )X(\tau )Y(\tau ),`$ (30b)
where $`\kappa (\gamma )=_0^\pi 𝑑\phi \sqrt{1+\gamma ^22\gamma \mathrm{cos}\phi }`$. In the limit $`\tau \mathrm{}`$, the solution of the system (30) reads

$`X(\tau )`$ $`\stackrel{\tau \mathrm{}}{}`$ $`\frac{1}{4\gamma }\tau ^1,`$ (31a)
$`Y(\tau )`$ $`\stackrel{\tau \mathrm{}}{}`$ $`{\displaystyle \frac{V_0}{4\gamma }}(4\gamma X_0)^{1/\mu }\left({\displaystyle \frac{4\kappa }{\kappa 4\gamma }}V_0+1\right)^{\nu /\mu }\tau ^{\kappa /4\gamma },`$ (31b)
where $`\mu =4\gamma /(\kappa 4\gamma )`$ and $`\nu =\kappa /(4\kappa )\mu `$. On the other hand, taking the limit $`\gamma 0`$ in Eqs. (30) and solving the corresponding system leads to the long-time behavior

$`X(\tau )`$ $`\stackrel{\tau \mathrm{}}{}`$ $`X_{\mathrm{}}\left[1+\epsilon _2(X_0,Y_0;\tau )\right],`$ (32a)
$`Y(\tau )`$ $`\stackrel{\tau \mathrm{}}{}`$ $`X_{\mathrm{}}\epsilon _2(X_0,Y_0;\tau ),`$ (32b)
where

$$\frac{X_{\mathrm{}}}{X_0}=\frac{1}{(1+V_0\chi )^{1/\chi }}$$
(33)
and $`\chi =4/\pi 1`$. In Eqs. (32), we have

$$\epsilon _2(X_0,Y_0;\tau )=\pi V_0^{(1+V_0\chi )^{1/\chi 1}}\mathrm{exp}\left(JX_{\mathrm{}}/X_0\right)\mathrm{exp}\left(\pi X_{\mathrm{}}\tau \right)$$
(34)
and

$$J=_0^{V_0}𝑑u\mathrm{ln}(u)\left[\frac{d^2}{du^2}\left(\frac{1/\chi +V_0}{1/\chi +u}\right)^{1/\chi }\right].$$
(35)
### IV.2 Molecular dynamics simulations

MD simulations have been implemented with systems of typically $`N=10^5`$ to $`4\times 10^5`$ spheres in two dimensions (discs). Periodic boundary conditions were enforced, and low densities considered, in order to minimize the excluded volume effects discarded at the Boltzmann level (note that these effects are necessarily transient since the density decreases with time).

Figure 4 compares the MD results obtained with $`\gamma =1/10`$ to the predictions of Eqs. (31) (for $`\gamma =1/10`$, the time decay of the “fast” particles is governed by the exponent $`\kappa /4\gamma 7.9`$). Although the large-time behaviors for $`X`$ and $`Y`$ are compatible with those given by Eqs. (31), it may be observed that the corresponding asymptotic regime is difficult to probe, even for large systems. The parametric plot \[or “trajectory” $`Y(X)`$\] shown in the inset is however in agreement with the relation $`YX^{\kappa /4\gamma }`$ deduced from Eq. (31).

We have also performed MD simulations for a mixture of moving and motionless particles ($`\gamma =0`$), where it is expected that the density $`X`$ of particles at rest decreases down to a nonvanishing value $`X_{\mathrm{}}`$. In the situation of an equimolar mixture ($`X_0=Y_0`$), we have $`V_0=1`$ so that according to Eq. (33), $`X_{\mathrm{}}/X_00.414`$. The MD simulations are in agreement with this scenario, and we find $`X_{\mathrm{}}/X_00.408`$ irrespective of the initial conditions for a system with initially $`N=2\times 10^5`$ particles. The results for the time dependence of $`X`$ and $`Y`$ are displayed in Fig. 5. We conclude that the numerical simulations are again in agreement with the prediction of Boltzmann’s kinetic theory.

## V Conclusions

We have shown that for some simple spatially homogeneous systems, characterized by a velocity distribution with a discrete velocity modulus spectrum, it is possible to find the exact solution for the nonlinear integral equation describing the dynamics of ballistic annihilation. These results, obtained at the level of a Boltzmann equation, have been validated by explicit comparison with molecular dynamics simulations in two dimensions.

For a single-velocity modulus distribution, the particle density of the model decays asymptotically as $`n(t)t^1`$, irrespective of space dimension. It was however rigorously shown that in one dimension, the decay is slower, $`n(t)t^{1/2}`$. This difference is a consequence of the fact that in one dimension strong dynamical correlations are created us\_one , which invalidate the approximation underlying Boltzmann’s dynamics. In higher dimensions, the Boltzmann equation becomes exact in the long-time limit.

In the case of a distribution with two different finite nonzero velocity moduli, we found that both particle densities decay for a large time according to a power law. The interesting feature is that the density of the slow particles decays as $`t^1`$, while the density of the fast particles decays more rapidly (e.g., as $`t^{(3+\gamma ^2)/4\gamma }`$ in three dimensions), with a nonuniversal exponent depending continuously on the velocity modulus ratio $`\gamma =c_1/c_2`$. A rough criterion for the crossover time separating the short- and long-time regimes has been given in Eq. (21).

Finally, the case $`c_1=0`$ leads to a particularly interesting behavior. Independently of the initial conditions, the densities of the moving and the motionless particles both decrease exponentially fast; however down to a nonzero value for particles at rest. This behavior is quite different from that observed in the one-dimensional case where the initial value of the density of motionless particles plays an important role in the long-time regime. This difference between one dimension and higher dimensions reflects once again the important role played by the dynamically created correlations for $`d=1`$.

The case with motionless particles can be viewed as a problem of ballistic annihilation of particles with one-velocity moduli moving in a random medium containing immobile traps (the motionless particles) that can capture a moving particle and then disappear. Here again, the situation can be compared to similar problems in diffusion limited annihilation where the presence of traps can modify the long-time dynamics from a power law to an exponential decay OBCM .

It would be interesting to compare the above theoretical predictions with some experimental data. Besides growth and coarsening problems, ballistic annihilation could model other physical systems such as, for example, the fluorescence of laser excited gas atoms with quenching on contact PRM . However, the correspondence between such experimental situations and our model is not yet close enough to allow comparison. We would be highly interested in the knowledge of other physical systems that could be described by the models studied here.

###### Acknowledgements.

This work was partially supported by the Swiss National Science Foundation. M.D. acknowledges the support of the Swiss National Science Foundation and of the CNRS. J.P. acknowledges the hospitality at the Department of Physics of the University of Geneva. We thank P. L. Krapivsky for bringing our attention to Ref. krap .


format: text

text:
# Hadronic structure in 𝜏⁻→𝜋⁻𝜋⁻𝜋⁺𝜈_𝜏 decays

## 1 INTRODUCTION

The $`\tau ^{}(3\pi )^{}\nu _\tau `$ decay has been the subject of much interest over the years. Because of the transformation properties of the weak current under parity and $`G`$-parity, $`\tau `$-lepton decay to an odd number of pions is expected to occur exclusively through the axial vector current, ignoring isospin-violating effects. Since $`J^P,`$ and $`I^G`$ are good quantum numbers of the hadronic current from weak decays, the weak axial current produces systems of pions with $`J^P=0^{}`$ or $`1^+,`$ odd $`G`$-parity and odd numbers of pions.

The $`(3\pi )^{}`$ system in $`\tau ^{}(3\pi )^{}\nu _\tau `$ is expected to be produced dominantly through the poorly understood $`J^P=1^+`$ $`a_1^{}(1260)`$ where $`a_1^{}(1260)\rho ^0\pi ^{}`$ in an $`S`$-wave decay. Although it is possible that the $`(3\pi )^{}`$ system can be produced through the $`J^P=0^{}`$ $`\pi ^{}(1300),`$ this is expected to be suppressed by the PCAC (partially conserved axial current) hypothesis. It has been known for some time that simple models do not provide a satisfactory description. For example, ARGUS  found that the model of Kühn and Santamaria  and that of Isgur, Morningstar, and Reader  did not describe their $`\tau ^{}(3\pi )^{}\nu _\tau `$ data satisfactorily.

Here we report the preliminary results of a model-dependent fit of the substructure in $`\tau ^{}\pi ^{}\pi ^{}\pi ^+\nu _\tau `$ decays. The work is based on the analysis of $`145,000`$ $`\tau ^{}\pi ^{}\pi ^{}\pi ^+\nu _\tau `$ events . In addition to fitting the substructure with various $`a_1^{}(1260)`$ decay modes, fits for the production of the $`(3\pi )^{}`$ system via the $`J^P=0^{}`$ $`\pi ^{}(1300)`$ are included. These results are compared to the substructure in $`\tau ^{}\pi ^0\pi ^0\pi ^{}\nu _\tau `$ previously reported .

Throughout this paper the charge conjugate versions of decays are implied for brevity and clarity.

## 2 EVENT SELECTION

This analysis is based on 4.67 fb<sup>-1</sup> of $`e^+e^{}`$ colliding beam data collected at center-of-mass energy energies $`E_{cm}10.6`$ GeV with the CLEO II detector at the Cornell Electron Storage Ring (CESR). This corresponds to $`4.3\times 10^6`$ $`e^+e^{}\tau ^+\tau ^{}`$ events. $`\tau ^{}\pi ^{}\pi ^{}\pi ^+\nu _\tau `$ events are selected from a topology in which the $`\tau ^+`$ decays into a single charged particle plus a $`\overline{\nu }_\tau `$ and the $`\tau ^{}`$ decays into three charged particles plus a $`\nu _\tau `$ (1 versus 3 charged-particle topology). Selected events are required to have four well-constructed charged tracks coming from the interaction region, and the net charge of the tracks is required to be zero. No more than one track is allowed to be identified as an $`e^{}`$ or $`e^+.`$ The most isolated track is identified, and the remaining three tracks are required to be $`90^{}`$ away from it. Events with photons (isolated electromagnetic showers not associated with a track) are rejected. Events that satisfy the hypothesis $`0.485<M_{\pi ^+\pi ^{}}<0.510`$ GeV are rejected in order to eliminate events containing $`K_s^0\pi ^+\pi ^{}`$ decays. Finally, the 2-dimensional cut shown in Figure 1 is used to reject events with missing tracks and/or energy along the beam pipe. Application of these cuts yields 145,000 $`\tau ^{}\pi ^{}\pi ^{}\pi ^+\nu _\tau `$ events.

The $`\pi ^{}\pi ^{}\pi ^+`$ effective mass spectrum resulting from these cuts is shown in Figure 2. The background, indicated by the shaded histogram, which amounts to about 11% of the events in the spectrum. Monte Carlo studies show that this background is from other $`\tau `$-decays: $`\frac{1}{2}`$ comes from $`\tau `$-events with missing $`\pi ^0`$(’s) and $`\frac{1}{2}`$ is from $`\tau `$-events with $`1`$ charged $`K.`$ Monte Carlo studies indicate that backgrounds from multihadronic $`q\overline{q}`$ events, $`\gamma \gamma `$ events, and QED events are negligible.

## 3 FIT

Although it is desirable to characterize the structure of the $`\pi ^{}\pi ^{}\pi ^+`$ system without model-dependent assumptions, this is very difficult in practice. In a model-independent approach one would assume a general form of the hadronic current that allows the separation and analysis of the complex magnitudes of the axial vector and pseudoscalar components in bins of $`Q^2,`$ $`s_1,`$ and $`s_2,`$ where $`Q^2`$ is the invariant mass-squared of the $`\pi ^{}\pi ^{}\pi ^+`$ system, $`s_1`$ is the invariant mass-squared of the $`\pi ^+\pi ^{}`$ combination that yields the larger value of mass-squared, and $`s_2`$ is the invariant mass-squared of the other $`\pi ^+\pi ^{}`$ combination. There would be no assumptions regarding resonances in the decay. However, in certain chiral limits, scalar effects can be induced in the current from the $`a_1^{}(1260)`$ and non-resonant contributions in the current and thus fake pseudoscalar effects . Secondly, simply binning in three dimensions is problematical: relatively small bins are required to study the $`\tau `$ decay structure and this leads to some bins with low statistics that may be crucial in resolving complicated resonance structure. By using a model-dependent fit, one can make reasonable assumptions about the resonant structure and study the likely existence of each resonance modeled while using the entire dataset in the fitting and without the need for binning.

Using an unbinned maximum likelihood fit, the substructure in $`\tau ^{}\pi ^{}\pi ^{}\pi ^+\nu _\tau `$ decays has been fit with the set of resonances given in Table 1. For comparison purposes, the current Particle Data Group  masses and widths are listed. For the $`f_0(600)`$ and $`f_0(1370)`$ resonance parameters given by the Törnqvist unitarized quark model  have been used. The same resonances and parameters have been used in the analysis of $`\tau ^{}\pi ^0\pi ^0\pi ^{}\nu _\tau `$ events  selected from the same CLEO II set of $`e^+e^{}\tau ^+\tau ^{}`$ events.

## 4 RESULTS

The results of the fit are given in Table 3, where the $`\beta _i`$ are complex coupling constants. The $`\pi ^{}\pi ^{}\pi ^+`$ invariant mass spectrum resulting from the fit is shown as the histogram in Figure 2.

As expected, the $`a_1^{}(1260)\rho ^0\pi ^{}`$ $`S`$-wave decay dominates. There is a significant amount of $`p`$-wave $`a_1(1260)`$ decay into isoscalars \[$`f_0(600),`$ $`f_0(1370)`$\] plus a pion. $`\tau `$-decay via the pseudoscalar $`\pi ^{}(1300)`$ is small, although its inclusion improves the fit. Dalitz plots for sequential $`\pi ^{}\pi ^{}\pi ^+`$ invariant mass regions are shown in Figure 3 for data and the fit.

The isoscalar sub-resonances in the decay of the $`a_1^{}(1260)`$ \[$`f_0(600),`$ $`f_2(1270),`$ and $`f_0(1370)`$\] play a large role in the fits. To demonstrate the need for these isoscalars, the data were fitted without them. Projections of the Dalitz plat without the presence of these isoscalars are shown in Figures 6 and 7. For the projection onto the $`s_1`$ axis (Figure 6), the lack of all three isoscalars results in a poor fit for $`Q^22.25`$ GeV$`^2.`$ For the projection onto the $`s_2`$ axis (Figure 7), the presence of $`f_0(600)`$ is needed for all regions of $`Q^2,`$ especially those of lower $`Q^2`$ values.

These results support the resonance substructure previously reported for the CLEO II $`\tau ^{}\pi ^0\pi ^0\pi ^{}\nu _\tau `$  data, which are based on $`15,000`$ lepton-tagged events. The fit results on these events are given in Table 3 for the reader’s convenience. Both analyses find that axial vector amplitudes with isoscalars, especially $`a_1^{}(1260)f_0(600)\pi ^{},`$ are large.

## 5 SUMMARY

This work is a new measurement of substructure in $`\tau ^{}\pi ^{}\pi ^{}\pi ^+\nu _\tau .`$ A model-dependent unbinned maximum likelihood fit on 145,000 events was performed using an axial vector component \[$`a_1(1260)`$ primary resonance and several subresonances, including isoscalars\] and a pseudoscalar component \[$`\pi (1300)`$ plus several subresonances\]. As expected, this decay is dominated by $`\tau ^{}a_1^{}(1260)\nu _\tau `$ with subsequent $`a_1\rho ^0\pi ^{}`$ $`S`$-wave decay. However, it is also apparent that other $`a_1^{}(1260)`$ decay modes are present. In particular, $`a_1^{}(1260)`$ decay into isoscalars \[$`f_0(600),`$ $`f_2(1270),`$ and $`f_0(1370)`$\] is significant. The $`a_1^{}(1260)f_0(600)\pi ^{}`$ is especially prominent. These results are consistent with the earlier CLEO II results for $`\tau ^{}\pi ^0\pi ^0\pi ^{}\nu _\tau `$ . Both analyses find that axial-vector amplitudes with isoscalars, especially $`a_1^{}(1260)f_0(600)\pi ^{},`$ are large. Regarding non-axial-vector contributions, inclusion of the $`J^P=0^{}`$ $`\pi ^{}(1300)`$ improves the fit to the data in the expected regions, but the need for it is not compelling. It is hoped that the future CLEO-c dataset will clarify some of these issues.

## 6 ACKNOWLEDGEMENTS

This work represents the work of many CLEO collaborators, especially Jason W. Hinson, Jon Urheim, and Alan J. Weinstein. The CESR staff provided us with excellent luminosity and running conditions. This work is supported by the National Science Foundation, the U. S. Department of Energy, and the Natural Sciences and Engineering Research Council of Canada.


text:
# Determination of 𝑆₁₇ from 7Be(𝑑,𝑛)8B reaction: CDCC analyses based on three-body model

## Abstract

The astrophysical factor $`S_{17}`$ for 7Be($`p,\gamma `$)8B reaction is reliably extracted from the transfer reaction 7Be($`d,n`$)8B at $`E=7.5`$ MeV with the asymptotic normalization coefficient method. The transfer reaction is accurately analyzed with CDCC based on the three-body model. This analysis is free from uncertainties of the optical potentials having been crucial in the previous DWBA analyses.

The solar neutrino problem is a central subject in the neutrino physics Bahcall . The major source of the high-energy neutrinos observed by solar neutrino detectors is 8B produced by the 7Be($`p,\gamma `$)8B reaction. The astrophysical factor $`S_{17}`$ for the reaction, however, is one of the most poorly determined reaction rates in the standard solar model; the latest recommendation for the factor $`S_{17}(0)`$ at zero energy, based on recent direct measurements, is $`19_2^{+4}`$ eV b Adelberger , i.e., 10–20 % error exists. This is far from our goal of determining $`S_{17}(0)`$ within 5 % error required in order to determine the neutrino oscillation parameters: the mass difference between $`\nu _e`$ and $`\nu _\mu `$ and their mixing parameter. The main difficulty in the direct measurement comes from ambiguities of determining the effective target thickness of the radioactive 7Be beam. Thus, indirect measurements of $`S_{17}(0)`$ are expected to be essential for determining $`S_{17}(0)`$ accurately.

The transfer reaction 7Be($`d,n`$)8B at low energies is an example of such indirect measurements; once the asymptotic normalization coefficient (ANC) of the overlap function between the 7Be and 8B ground states is determined from the data of the transfer reaction, $`S_{17}(0)`$ can be accurately derived from the ANC, as long as the reaction is peripheral Xu . The reaction has been measured at $`E=7.5`$ MeV and the ANC is extracted with the distorted wave Born approximation (DWBA) Liu . The $`S_{17}(0)`$ obtained with the ANC is $`27.4\pm 4.4`$ eV b, leading to inconsistency with the recommended value. Although the reaction is found with DWBA to be indeed peripheral Gagliardi ; Fernandes , distorting potentials used in DWBA are quite ambiguous Gagliardi ; Fernandes , which prevent one from extracting $`S_{17}(0)`$ accurately. In particular, uncertainties of the $`d`$-7Be optical potential bring about large errors for $`S_{17}`$, typically 30 % in magnitude. The origin of the large ambiguity of distorting potentials is that these are derived from proton and deuteron optical potentials for different targets and/or energies.

In the present paper, we analyze 7Be($`d,n`$)8B at 7.5 MeV with the three-body model, $`p+n+^7`$Be, assuming 7Be to be an inert core. This treatment of the system is justified by Refs. Nunes1 ; Nunes2 , where core excitations of 7Be are shown to be negligible. An advantage of this analysis is that we do not need the ambiguous $`d`$+7Be optical potential in the entrance channel and the ambiguous $`n`$+8B one in the exit channel as shown below.

The three-body dynamics in the entrance channel are explicitly treated by means of continuum-discretized coupled-channels (CDCC) method CDCC ; Piya , the theoretical foundation of which is given in Ref. CDCC-foundation . This theory has been established as a method of solving the three-body system with good accuracy, and extensively applied for various reactions CDCC ; Surrey . Previous CDCC calculation showed that explicit treatment of breakup channels is essential in describing deuteron induced reactions CDCC . The CDCC thus provides a precise description of the wave function in the entrance channel, i.e., $`d`$-7Be system.

The effective Hamiltonian for the entrance channel based on the three-body model contains the optical potential between $`N`$ ($`p`$ or $`n`$) and 7Be. Data of the neutron elastic scattering are available for target 7Li, the mirror nucleus of 7Be, at 4 MeV, approximately half the deuteron energy considered here n+A.dat . First the $`n`$-7Li potential is determined accurately from the data. The potential is then used as an input in CDCC calculation for deuteron elastic scattering on 7Li at 8 MeV, and the numerical result is compared with the experimental data d+A.dat . This is a good test for the neutron optical potential determined above, which is only an input in CDCC calculation for the entrance channel. As for the exit channel, the three-body dynamics are treated with the adiabatic approximation, after testing its accuracy with CDCC calculation. The ANC is then obtained with reliable distorted wave functions in both the entrance and exit channels.

The transition amplitude for the transfer reaction, based on the three-body model ($`p+n+^7`$Be), is

$$T_{\mathrm{fi}}=S_{\mathrm{exp}}^{1/2}<\mathrm{\Psi }_\mathrm{f}^{()}|V_{\mathrm{np}}|\mathrm{\Psi }_\mathrm{i}^{(+)}>.$$
(1)
The three-body wave function $`\mathrm{\Psi }_\mathrm{i}^{(+)}`$ is a solution of the Schrödinger equation $`(H_\mathrm{i}E)\mathrm{\Psi }_\mathrm{i}^{(+)}=0`$ for the three-body Hamiltonian

$$H_\mathrm{i}=K_{\mathrm{np}}+K_{\mathrm{dBe}}+V_{\mathrm{np}}(𝐫_{\mathrm{np}})+U,$$
(2)
with $`U=U_{\mathrm{nBe}}(𝐫_{\mathrm{nBe}})+U_{\mathrm{pBe}}(𝐫_{\mathrm{pBe}})`$. Here $`𝐫_{\mathrm{XY}}`$ is the coordinate of nucleus X relative to nucleus Y. The potential $`V_{\mathrm{np}}`$ is the interaction between $`n`$ and $`p`$, $`U_{\mathrm{pBe}}`$ ($`U_{\mathrm{nBe}}`$) is the proton (neutron) optical potential for the target 7Be at half the deuteron incident energy, and $`K_{\mathrm{np}}`$ and $`K_{\mathrm{dBe}}`$ show kinetic energy operators for two-body systems denoted by the subscripts. The nuclear part of $`U_{\mathrm{pBe}}`$ is assumed to be the same as $`U_{\mathrm{nBe}}`$. The Coulomb part of $`U_{\mathrm{pBe}}`$ is treated approximately by replacing the coordinate $`r_{\mathrm{pBe}}`$ by $`r_{\mathrm{dBe}}`$. Effects of Coulomb breakup of deuteron, not included in the usual treatment above, are found to be quite small in the present system.

The wave function $`\mathrm{\Psi }_\mathrm{i}^{(+)}`$ is obtained with CDCC, that is, by solving the Schrödinger equation in a model space. In CDCC, deuteron breakup states are classified with linear and angular momenta, $`k`$ and $`\mathrm{}`$, and truncated into $`0kk_{\mathrm{max}}`$ and $`0\mathrm{}\mathrm{}_{\mathrm{max}}`$, respectively. The $`k`$ continuum $`[0,k_{\mathrm{max}}]`$ is further divided into bins with a common width $`\mathrm{\Delta }`$. The total wave function is expanded, in terms of the deuteron ground state and the discretized breakup states, into

$$\mathrm{\Psi }_\mathrm{i}^{(+)}=\underset{i}{}\varphi _i(𝐫_{\mathrm{np}})\chi _i(𝐫_{\mathrm{dBe}}),$$
(3)
where $`\varphi _0`$ is the deuteron ground state and $`\varphi _i`$ is the $`i`$-th discretized breakup state obtained by averaging continuous breakup states in the $`i`$-th bin. The coefficient $`\chi _i`$ represents a center-of-mass motion of $`n`$-$`p`$ pair in the $`i`$-th state. Inserting this form into the three-body Schrödinger equation leads to a set of coupled differential equations,

$$(E_iK_{\mathrm{dBe}})\chi _i=\underset{j}{}F_{ij}\chi _j$$
(4)
with $`F_{ij}=<\varphi _i|U|\varphi _j>_{𝐫_{\mathrm{np}}}`$ and $`E_i=Ee_i`$, where $`e_i`$ is an intrinsic energy of the $`i`$-th $`n`$-$`p`$ state. The coupled equations are soluble, since they have a compact kernel in its integral equation form. The precise formulation of CDCC is shown in Ref. Piya . The present model space is $`k_{\mathrm{max}}=1.7`$ fm-1, $`\mathrm{}=0,2`$ and $`\mathrm{\Delta }=1.7/40`$ fm-1. The CDCC solution converges at these values, as the model space is enlarged.

The exit channel is also treated in the three-body model, that is, the exit channel wave function $`\mathrm{\Psi }_\mathrm{f}^{()}`$ is determined by the three-body Hamiltonian

$$H_\mathrm{f}=K_{\mathrm{pBe}}+K_{\mathrm{nB}}+V_{\mathrm{pBe}}(𝐫_{\mathrm{pBe}})+U_{\mathrm{nBe}}(𝐫_{\mathrm{nBe}}).$$
(5)
In the three-body model, 8B is treated by the two-body ($`p`$+7Be) model with the potential $`V_{\mathrm{pBe}}`$. The spectroscopic factor $`S_{\mathrm{exp}}`$ in Eq. (1) is introduced by taking account of the incompleteness of the model. It should be noted that $`H_\mathrm{f}`$ does not contain $`V_{\mathrm{np}}`$ since the interaction is already treated as a transition operator in Eq. (1). In general, the distorting potential $`U_{\mathrm{nBe}}`$ between an outgoing neutron and 7Be differs from the corresponding one in $`H_i`$, since an outgoing neutron has a different velocity from an incoming deuteron in the ($`d,n`$) reaction.

In principle, $`H_\mathrm{f}`$ allows transitions between the ground and continuum states of the p+7Be system in the exit channel. However, effects of the transitions may be simply estimated with the adiabatic approximation, since the ground state of 8B has a binding energy (0.137 MeV) considerably smaller than an energy (4.18 MeV) of outgoing neutron. Following Johnson and Soper Johnson-Soper , we replace the Hamiltonian of the p+7Be system by the binding energy of 8B. Errors of the adiabatic approximation are estimated with CDCC in the exit channel scattering, 8B($`n,n`$)8B, at $`E=4.18`$ MeV. The breakup effect itself is not so large and errors of the approximation are less than 3 % at forward angles. The Johnson-Soper approximation leads to a simple form $`\mathrm{\Psi }_\mathrm{f}^{()}=\chi _{\mathrm{nB}}^{()}\varphi _{\mathrm{pBe}}`$, where $`\varphi _{\mathrm{pBe}}`$ is the wave function of 8B in its ground state and $`\chi _{\mathrm{nB}}^{()}`$ is the wave function of outgoing neutron distorted by the potential $`U_{\mathrm{nBe}}(r_{\mathrm{nB}}8/7)`$, where the zero-range approximation is made to the transition amplitude (1). We discuss the use of this approximation below. It should be noted that in DWBA $`\chi _{\mathrm{nB}}^{()}`$ is determined by the elastic scattering of neutron from 8B, of which no measurement has been done so far. On the other hand, the three-body model approach can avoid this difficulty, as in the entrance channel.

The transfer reaction 7Be($`d,n`$)8B is calculated with the zero-range approximation with its finite-range correction fr-corr . The integration over $`𝐫_{\mathrm{pBe}}`$ in $`T`$ is made up to a large value 40 fm, since the transferred proton is very weakly bound in 8B. The finite-range correction for the transition amplitude (1) including deuteron ground and breakup channels is straightforward; the resultant correction for the $`i`$-th channel keeps the standard form by regarding $`F_{ii}`$ and $`e_i`$ as the potential and the intrinsic energy of the entrance channel. This prescription is tested by doing finite-range DWBA calculation for the deuteron ground channel which is a main component of the transition amplitude (1). The result of the prescription above agrees with that of the full finite-range calculation within 2 % error at forword angles.

The wave function $`\varphi _{\mathrm{pBe}}`$ is calculated with four types of $`V_{\mathrm{pBe}}`$ Kim ; Tombrello ; Robertson ; Esbensen . We determine the spectroscopic factor $`S_{\mathrm{exp}}`$ comparing the calculated 7Be($`d,n`$)8B cross section with the experimental one, for every type of the four potentials. The astrophysical factor $`S_{17}(0)`$ at zero energy is then obtained from the $`S_{\mathrm{exp}}`$ with

$$S_{17}(0)=\frac{S_{\mathrm{exp}}b^2}{0.026}$$
(6)
in the ANC method Xu , where $`b`$ is defined with the Whittaker function $`W`$ as $`\varphi _{\mathrm{pBe}}(r_{\mathrm{pBe}})=bW(r_{\mathrm{pBe}})`$ at $`r_{\mathrm{pBe}}`$ larger than the range of the nuclear force $`V_{\mathrm{pBe}}`$ between $`p`$ and 7Be. It should be noted that if the reaction is peripheral, ANC should be stable against the change of the 8B internal wave functions, i.e., that of the parameter set of single-particle potentials. Thus, one can estimate the error of the ANC calculation from the deviation of ANC (or $`S_{17}`$) with four different models of 8B above.

Figure 1 shows the result of the optical potential search for neutron scattering at 4.08 MeV from 7Li, the mirror nucleus of 7Be. The resultant potential shows a good agreement with data n+A.dat . The optical potential is then applied for deuteron scattering at 8.0 MeV from 7Li. The CDCC calculation with the potential again gives a good agreement with data d+A.dat  at angles $`\theta <60^{}`$, as shown in Fig. 2. The potential is shown to be reliable especially at the forward angles; we use it in $`H_i`$ and obtain the proper wave function for the entrance channel. As for the exit channel, on the other hand, we need the optical potential of the $`n`$+7Li scattering at $`E=4.18`$ MeV to determine $`H_f`$. The data is available at $`E=4.26`$ MeV that is closest to the proper energy. Figure 3 shows the result of the optical potential search for the scattering. The resulting potential well reproduces the data n+A.dat2  at forward angles.

Figure 4 shows our result for 7Be($`d,n`$)8B cross section using $`V_{\mathrm{pBe}}`$ of Kim et al Kim ., compared with the experimental data Liu . The solid and dashed lines represent the calculated results with $`S_{\mathrm{exp}}=1.0`$ and 0.849, respectively. Parameter sets of the optical potentials at 4.08 and 4.26 MeV and the single-particle potential of Kim et al. used in the calculation are listed in Table 1 together. At forward angles $`\theta <60^{}`$, the calculated cross section well reproduces the data Liu  with the spectroscopic factor $`S_{\mathrm{exp}}=0.849`$, leading to $`S_{17}(0)=21.36`$ eV b. It was found that deuteron breakup states play important roles not only in determining distorting potentials but also in the transfer process. In fact, when the deuteron breakup components is set to zero in $`\mathrm{\Psi }_\mathrm{i}^{(+)}`$, the resultant transfer cross section is reduced by 10 % at the forward angles and, more seriously, the angular distribution cannot be reproduced correctly. The components are obviously not included in the framework of the standard DWBA. We thus conclude that the three-body model approach is inevitable.

We show in Table 2 the list of calculated $`S_{\mathrm{exp}}`$, $`b`$, and $`S_{17}(0)`$ for different $`V_{\mathrm{pBe}}`$. One sees that the calculated values of $`S_{17}(0)`$ are almost consistent, which shows that, as mentioned above, the present reaction is peripheral and the ANC method works well. Taking account of the theoretical errors of ANC, the adiabatic approximation (AD) and the finite-range correction (FRC) and of the systematic error of the experimental data Liu  on the 7Be($`d,n`$)8B cross section, we obtain $`S_{17}(0)=20.96_{0.3}^{+0.4}`$ (ANC) $`\pm 0.63`$ (AD) $`\pm 0.42`$ (FRC) $`\pm 2.7`$ (expt) eV b, in consistent with both the recommended values Adelberger  and the recent result $`S_{17}(0)=22.3\pm 1.2`$ eV b of acccurate direct measurement junghans .

In summary, the present analyses based on the three-body model are free from uncertainties of the optical potentials in both the entrance and exit channels which were the most essential problem in the previous DWBA analyses. The deuteron breakup process in the incident channel is significant for determining $`S_{17}(0)`$ within 5 % error required from the neutrino physics. We then conclude that the three-body model approach is essential and necessary. The present analyses provide a precise value $`S_{17}(0)=20.96_{1.3}^{+1.4}(\mathrm{theor})\pm 2.7(\mathrm{expt})`$ eV b; the theoretical ambiguity of $`S_{17}(0)`$ is 6–7 % slightly beyond the required accuracy. However, we can reduce the theoretical error to $`2\%`$ coming from ANC only, if we do the finite-range calculation with accurate $`\mathrm{\Psi }_\mathrm{i}^{(+)}`$ and $`\mathrm{\Psi }_\mathrm{f}^{()}`$ derived by CDCC. Such full-fledged calculations are highly expected. In the present analyses, however, the experimental error (13 %) is even larger than the theoretical one. It is expected that the peripheral properties, essential for the ANC method, become insufficient, as the incident energy increases Fernandes . Thus, accurate measurements on 7Be($`d,n`$)8B and 7Li($`d,d`$)7Li at about a few tens of MeV and on 7Li($`N,N`$)7Li at the half the corresponding deuteron incident energy are highly expected; the proton and deuteron elastic scattering are necessary to determine the nucleon optical potential accurately.

The authors would like to thank T. Motobayashi, T. Kajino, Y. Watanabe, and M. Kawai for helpful discussions. This work has been supported in part by the Grants-in-Aid for Scientific Research of the Ministry of Education, Science, Sports, and Culture of Japan (Grant Nos. 14540271 and 12047233).


text:
# 1 Introduction

## 1 Introduction

The goal of current investigations in heavy-ion collisions is to establish the existence of the quark-gluon plasma (QGP) phase and to study its properties. Several signals have been proposed as indicators for its creation. However, no unambiguous proof has been established so far. A large part of the problem lies in the difficulty to separate effects of the evolution of the hot and dense medium from characteristic changes in the physics of reactions taking place inside this medium. Therefore, it is mandatory to aim at a consistent description of as many observables as possible within a single model of the medium evolution.

In a recent paper , we have proposed a model for the evolution of a fireball assuming local thermal equilibrium and isentropic expansion. In this approach, the evolution is constrained by two major pieces of information: first, by measured hadronic momentum spectra and Hanbury-Brown Twiss (HBT) interferometry data which reflect the freeze-out state; secondly, by information on the Equation of State (EoS) obtained in lattice simulations and represented in terms of a quasiparticle picture. We have demonstrated that this evolution scenario is in agreement with data on dilepton emission measured by the CERES collaboration . We have also used the same model to describe the emission of hard thermal photons  and found good agreement with the data measured by WA98 .

No explicit statement about the hadrochemical composition of the fireball is made in this model. Instead, we use phenomenological arguments to include the effects of enhanced pion phase space density into the EoS of hot hadronic matter. On the other hand, statistical models are extremely successful in describing the measured ratios of different hadron species for a range of collision energies from SIS to RHIC (see e.g. ). It is the aim of this paper to show that our fireball evolution model is consistent with the idea of statistical hadronization.

The paper is organized as follows: First, we outline the version of statistical hadronization used in our model. After comments on technical details, we present results and compare to data both for SPS and RHIC conditions. We then investigate the sensitivity of these results to the model parameters, specifically we study the role of possible in-medium modifications of particle properties. We conclude by sketching a possible scenario in which in-medium modifiactions are possible while statistical hadronization is still in agreement with the data.

## 2 The statistical hadronization model

The fireball evolution model outlined in  contains a phase transition from partonic to hadronic matter. If one is interested in the composition of the hadronic system, one has to introduce a model of hadronization. In the present paper, we assume that the particle content of the fireball at the critical temperature $`T_C`$ can be found by considering a system of (non-interacting) hadrons in chemical equilibrium, described by the grand canonical ensemble, which subsequently undergoes decay processes. This is based on ideas presented in , where it was shown that a huge number of hadron ratios can be explained by a fit of only two model parameters, temperature $`T`$ and baryochemical potential $`\mu _B`$.

Clearly, the notion of non-interacting hadrons close to the phase boundary is not suitable for the description of the phase transition thermodynamics on a fundamental level, but should rather be regarded as an effective prescription. The underlying assumption is that the level-spacings and degeneracies of the QCD excitation spectrum along with the energy scale provided by the temperature already contain enough information to determine the relative contributions of different hadron species and that details of the interactions can either be absorbed into hadronic quasiparticle properties such as temperature dependent masses and widths or average out.

The calculation of hadronic abundancies in this framework is a two step process: First, we have to specify the properties of the medium (the QGP) at the transition in terms of the densities of quarks, antiquarks and gluons (This is most conveniently done by specifying the temperature $`T_C`$ and baryo- and strange chemical potentials $`\mu _B,\mu _S`$). These parameters can be fit to the observed hadron ratios (as done in ), calculated within a dynamical model (as done in the present paper for SPS) or inferred from experiment (as done in the present paper for RHIC).

Once the densities in the QGP have been specified, the hypothesis of statistical hadronization is used to map those into the hadron ratios observable in experiment. In detail, we use the following procedure to calculate ratios for SPS:

Employing the grand canonical ensemble, we expect the density for each particle species $`n_i`$ to be given by

$$n_i=\frac{d_i}{2\pi ^2}_0^{\mathrm{}}\frac{p^2dp}{\mathrm{exp}\{[E_i(p)\mu _i]/T_C\}\pm 1}.$$
(1)
Here, $`d_i`$ denotes the degeneracy factor of particle species $`i`$ (spin, isospin, particle / antiparticle), the +(-) sign is used for fermions (bosons) and $`E_i(p)=\sqrt{m_i^2+p^2}`$. $`m_i`$ stands for the particle’s vacuum mass. We use a value of 170 MeV for the critical temperature $`T_C`$ as determined in lattice calculations for two light and one heavy flavours  and also used in the fireball evolution model. The chemical potential $`\mu _i`$ takes care of conserved baryon number $`B_i`$ and strangeness $`S_i`$ for each species:

$$\mu _i=\mu _BB_i\mu _SS_i$$
(2)
We neglect a (small) contribution $`\mu _{I_3}I_i^3`$ coming from the isospin asymmetry in the colliding nuclei.

The baryochemical potential $`\mu _B`$ is then fixed by the requirement that the net number of baryons inside the thermalized region is equal to the number of collision participants $`N_{part}`$ if the total volume $`V`$ is known:

$$V\underset{i}{}n_iB_i=N_{part}.$$
(3)
Similarly, strangeness conservation demands

$$V\underset{i}{}n_iS_i=0.$$
(4)
Thus, the only parameter of the model is the fireball volume $`V`$. At the phase transition, however, this volume can be determined from the EoS of the QGP evaluated at $`T=T_C`$ using the entropy density $`s(T)`$ as

$$V(T_C)=S_{tot}/s(T_C),$$
(5)
if the total entropy content $`S_{tot}`$ of the fireball is known. This quantity, however, can be obtained from measuring charged particle multiplicities $`N^+`$ and $`N^{}`$ in suitable rapidity bins and calculating

$$D_Q=\frac{N^+N^{}}{N^++N^{}}.$$
(6)
The quantity $`D_Q`$ stands for the inverse of the specific entropy per net baryon $`S/B`$, and the product $`D_Q(S/B)`$ roughly measures the entropy per pion . For SPS collisions at 160 AGeV, we find an entropy per net baryon $`S/B=26`$ for central collisons.

At RHIC, hadron ratios are measured at midrapidity. It is known that the distribution of baryon number across rapidity is very inhomogeneous at collider energies, hence a prediction of the quark and gluon densities within a globally averaged framework as outlined for the SPS case using Eqs. (3) and (5) is bound to fail. For the sake of simplicity, we will not attempt to calculate $`\mu _B`$ in a (more complicated) dynamical model but rather infer it from the experimentally accessible specific entropy end entropy, keeping in mind that our results for RHIC do not test a dynamical model but only the statistical hadronization hypothesis (in principle, however, the baryon distribution in rapidity is calculable, e.g. ).

For RHIC 6% central Au-Au collisions at 130 AGeV, the specific entropy $`S/B=220`$ at midrapidity is substantially higher due to the larger particle multiplicity and the smaller net baryon content in the central region. The transition volume is then estimated by using the EoS to determine the volume at which the measured entropy in the midrapidity slice correspond to the temperature $`T_C`$ if a backward extrapolation of the expansion in time is made. Note, however, that the transition volume is only relevant for excluded volume corrections and can be neglected for a first order approximation since we are interested in hadron ratios and not total abundancies.

Thus, all ingredients entering Eq. (1) are determined (without additional free parameters) and we can evaluate the expression for a suitable choice of hadrons and resonances.

We include all mesons and mesonic resonances up to masses of 1.5 GeV and all baryons and baryonic resonances up to masses of 2 GeV. This amounts to 30 (strange and nonstrange) mesonic states and 36 (nonstrange to multistrange) baryonic states. In order to compare to experimental results, we calculate their decay into particles which are long-lived as compared to the fireball, such as $`\pi ,K,\eta ,N,\mathrm{\Lambda },\mathrm{\Sigma }`$ and $`\mathrm{\Omega }`$.

In order to account for interactions between particles, which at small distances become repulsive, we assume a hard core radius $`R_C`$ of 0.3 fm for all particles and resonances. The corresponding excluded volume $`V_{ex}=_iN_iV_{ex}^i`$, with $`V_i^{ex}=\frac{4\pi }{3}R_C^3`$ for all species, is subtracted from the volume obtained in Eq. (5), which in turn affects the total number of produced particles. As the excluded volume itself depends on the total number of particles, we iterate the correction for a self-consistent result. The hard core radius of 0.3 fm for protons is determined by comparison with p-p collisions . In the absence of such information for the other mesons and baryons, we assume its universality.

All data on particles is taken from . For many higher-lying states, the properties as well as the decay channels are poorly known. In these cases, we proceed as follows: If a quantity (e.g. masses and widths) is given only within a certain range, the arithmetic mean of this range is used in the model. Decay channels which are reported to be ’seen’ are assumed to receive equal contributions from the branching ratio which is left after all known channels have been accounted for. Branching ratios less than 1% have been neglected. Decay chains (such as $`a_2\rho \pi \pi \pi \pi `$) have been followed through. For resonances with large width, we integrate Eq. (1) over the mass range of the resonance using a Breit-Wigner distribution.

## 3 Some remarks on thermodynamics

As mentioned earlier, the notion of a non-interacting hadronic resonance gas cannot be expected to describe the full thermodynamics of the system close to the phase boundary. Therefore, one cannot expect all thermodynamical parameters to match smoothly coming from above and below the phase transition. An ensemble of free hadrons should therefore not be used to calculate thermodynamical properties of the hot hadronic medium, as there is strong evidence for re-scattering processes in the hadronic phase from dilepton emission (see e.g. our result in ) and these interactions certainly modify the thermodynamics.

This becomes apparent once we compare thermodynamic parameters directly. Coming from above $`T_C`$, we find from the quasiparticle model of the QGP

$$p/T_C^4=0.3,ϵ/T_C^4=7.0\text{and}s/T_C^3=7.3.$$

Calculating the corresponding hadronic resonance gas values at $`T_C`$ (using the ensemble of particles and resonances described above), one obtains

$$p/T_C^4=1.0,ϵ/T_C^4=6.3\text{and}s/T_C^3=7.3.$$

The crucial observation is that the entropy density matches smoothly (one should keep in mind, however, that there are errors given by the extraction of $`T_C`$ from lattice data, the extrapolation of lattice data to physical quark masses and the poor knowledge of high-lying resonances). This implies that we have not introduced a gross violation of the number of available degrees of freedom (which is important, as the goal of this study is to obtain particle numbers out of an *isentropic* expansion scenario with a crossover or second order phase transition). If the entropy density did not match at $`T_C`$, we would have introduced a first order phase transition which would certainly lead to a breakdown of the simple, ’instantaneous’ hadronization scenario used here.

On the other hand, in the pressure and energy density we find discontinuities at $`T_C`$. The difference of the energy density in both phases can be regarded as acceptable given the uncertainties, but in the case of the pressure the discrepancy is large. It is, however, certainly conceiveable that attractive interactions in the system lower the pressure. Thermodynamical consistency then fixes the energy density via $`ϵ+p=sT`$, achieving continuity in both quantities. A more fundamental approach in order to achieve continuity of all parameters across the phase transition, however, faces other difficulties:

* Any interacting model of hadrons will have to introduce additional coupling constants, which, at least for the high-mass resonance states, are poorly known.
* Incorporating all *known* interactions is of limited help: If one considers the free ensemble at $`T_C`$, one finds e.g. that the sum of all excited $`\mathrm{\Delta }`$ states is larger than the number of $`\mathrm{\Delta }(1232)`$. This is also true for the nucleon and nuclear resonances. Hence, resonance states are no small contribution and cannot be neglected. If one is interested in the number of pions resulting after resonance decays, the role of heavy resonances is even further enhanced due to very pion-rich decay channels of these heavy states.
* Accepting the limited knowledge of heavy-resonance properties, one might think of using some effective *ad hoc* prescription of restoring thermodynamical consistency. However, there is no unique way of doing so, therefore the model looses all predictive power.

In contrast, statistical hadronization is a simple prescription and completely determined by at most two parameters, which in the present context are given by the fireball evolution model. Therefore, it appears useful to accept the limitations of the present model for the benefit of a parameter-free calculation.

## 4 Results

The resulting hadron ratios are shown in Fig. 1 for the case of 158 AGeV central Pb-Pb collisions at SPS, compared with the experimentally measured values . One observes that the overall agreement with data is satisfactory with few exceptions, though not quite as good as in e.g.  where a two parameter fit to the data was performed. Note that not only the ratios of hadron yields agree to experiment in the present approach but also the absolute numbers, as the baryochemical potential $`\mu _B`$ is explicitly linked to the (known) number of participants.

The calculation yields a baryochemical potential $`\mu _B=265`$ MeV and a strange chemical potential $`\mu _S=59`$ MeV. Note that these quantities depend on the number of resonances included into the calculation; therefore, a direct comparison with the values obtained in other types of models is not necessarily meaningful.

For central collisions at RHIC at 130 AGeV beam energy we present the results in Fig. 2 and compare to experimental data measured around midrapidity . As the distribution of net baryon number in rapidity is very inhomogeneous at RHIC, the naive application of Eq. (3) fails, as we obtain a low entropy per baryon, $`S/B=75`$. At midrapidity, however, this ratio is closer to $`S/B=220`$. Using this value to calculate the ratios in a suitable interval around midrapidity, we find much better agreement to the data (see Fig. 2). Naturally, Eqs. (3) and (5) cannot be used in this case. Instead, one has to consider the volume of the rapidity slice in question and the number of baryons found inside this interval (which can either be predicted by a model of measured in experiment). Starting with this number and $`S/B`$, all ratios of produced hadrons within this rapidity slice are then predicted. Clearly, a more detailed fireball evolution model for the RHIC scenario is desirable to increase the predictive power.

## 5 In-medium modifications

In the last section, we have demonstrated that the idea of statistical hadronization combined with our fireball evolution model is able to achieve good agreement in comparison with data, provided that one uses the vacuum masses and widths of all resonances in Eq. (1). On the other hand, we have used the same fireball evolution to calculate dilepton emission and here we found that a significant broadening of the $`\rho `$ meson was responsible for the observed enhancement in the invariant mass region below 700 MeV . The underlying thermal field theory calculations of the modifications of the vector meson properties at finite temperature  and density  indicate not only a modified $`\rho `$ but also broadening and mass shift of the $`\omega `$ and broadening of the $`\varphi `$ due to the interaction with the medium.

All these calculations (as any other perturbative expansion) become unreliable in the vicinity of the phase transition, therefore these results cannot strictly be taken over to the current calculations where we require these quantities close to $`T_C`$. However, we may take them as a hint that two possible modifications of particle properties in a hot and dense medium may take place:

*Mass shifts* of particles in the medium are commonly related to the restoration of chiral symmetry. In e.g. , by investigating scale invariance of an effective Lagrangian, the in-medium scaling laws

$$m_\rho ^{}/m_\rho m_\omega ^{}/m_\omega m_N^{}/m_N\left(\overline{q}q^{}/\overline{q}q\right)^{1/3}$$
(7)
were established (Brown-Rho scaling). In this equation, asterisks denote quantities at finite density, $`m_\rho ,m_\omega `$ and $`m_N`$ the masses of the $`\rho `$ and $`\omega `$ meson, $`m_N`$ is the nucleon mass and $`\overline{q}q`$ stands for the chiral condensate.

*Decay widths* of particles are in general increased in a medium due to the presence of new interaction channels. In  e.g., it was shown that the $`\omega `$ meson resonance experiences strong broadening in a hot environment due to the presence of the scattering process $`\omega \pi \pi \pi `$. Even in the absence of such effects, the decay width at finite temperature is enhanced if the decay products are bosons due to the presence of bosons of the same type in the heat bath (Bose-enhancement). As pions are the most abundant species in a thermal environment and most decays involve one or more pions, this effect should influence almost all resonances in a hot medium. The modification of decay widths is only relevant in Eq. (1) if the in-medium width is a sizable fraction of the particle mass. In this case, a large contribution to the particle yield comes from masses lower than the peak mass in the Breit-Wigner distribution, which are exponentially enhanced. This enhancement more than counterbalances the suppression of contributions shifted into the higher mass region. Therefore, an increase in the decay width of broad resonances acts similarly as a mass reduction.

There is a third possibility of in-medium particle properties which has no direct influence on the particle spectral function (and is therefore not visible in the dilepton data): The binding potential between the constituents of hadrons could be partially screened by thermal fluctuations, leading to an increase of the hadronic core radius.

Such *increased radii* would lead to an enhanced excluded volume correction. This effect hardly influences the results of , as the fireball volume is implicitly determined by matching the fitted baryochemical potential to the number of participants, but in the present approach we can expect to observe the influence of increased hard core radii, as the volume is kept fixed.

In a first step, we examine the effects of in-medium mass shifts in a qualitative way by tentatively multiplying the vacuum masses of all hadrons with the exception of the pion by a constant $`c`$. The result is shown in Fig. 3.

The result shows that even a moderate mass reduction of 10% in the medium is not in line with the observed hadron ratios. In particular, particle-antiparticle ratios are strongly affected. In the case of $`\overline{p}/p`$, one might argue that the relevant inelastic (annihilation) cross section is not small as compared to the elastic one and therefore this ratio cannot really be fixed at $`T_C`$, but must be adjusted dynamically in the subsequent evolution. Indeed, it was shown in  that this is possible if one takes the statistical hadronization prediction as an initial condition for rate equations. It is unclear if this is still possible for different initial conditions, but even if this is the case, this is not an option for the multistrange particle/antiparticle ratios.

The overall behaviour of the result can be qualitatively understood as follows: The reduced nucleon mass implies a lower value of the baryochemical potential in order to produce the observed number of participants, this in turn affects single and double strange particles and implies changes in $`\mu _S`$ via Eq. (4). Therefore, ratios of particles and antiparticles with 2 or 3 non-strange valence quarks, such as $`\overline{p}/p`$ or $`\overline{\mathrm{\Lambda }}/\mathrm{\Lambda }`$ are most affected. Multistrange particle/antiparticle pairs follow the trend, though in a way less pronounced.

In a second run, we investigate the effect of thermal broadening of resonances, increasing all widths by a constant multiplicative factor $`c_\mathrm{\Gamma }`$. The resulting hadron ratios are shown in Fig. 4.

One observes the same qualitative behaviour as for a mass reduction, as we have argued before. The effects of the increased width are, however, less dramatic. An increase by 20% in the width of all resonances is still in line with all data except $`\overline{p}/p`$ and even an increase of 50% is still acceptable for most of the ratios. This is reassuring, as there is almost certainly thermal broadening of resonances in a hot medium.

As the behaviour for both broadening of resonances and mass reduction (which are both in line with the dilepton data) is qualitatively the same for the hadron ratios, the effects of reduced masses cannot be compensated by introducing additional broadening. So if any of the effects fails in the description of the data, a combination of both will also fail.

In the third run, we explore the effect of different choices for the core radius $`R_C`$ on the ratios. In order to investigate the sensitivity of our results to the initial choice of $`R_C`$, we do not only consider thermally increased radii but also a reduced initial choice. In Fig. 5, we show the model predictions for a reduction of $`R_C`$ by 25% and for an increase of the same amount.

We observe that the overall sensitivity of the resulting hadron ratios to the hadronic core radius is rather weak. On the other hand, an enhanced excluded volume correction (dotted line in Fig. 5) acts in a rather peculiar way: In order to arrive at the same number of participants, the baryochemical potential has to *increase*. This effect is opposite to the behaviour for in-medium mass reductions or increase of decay width. Therefore, one might expect that the net effect of both thermal broadening of resonances *and* increased core radius due to a thermally screened binding potential is moderate and partial compensation may occur.

In order to test this conjecture, we study a scenario in which the width of resonances has been increased by 50% and simultaneously the core radius has been set to 0.45 fm instead of 0.3 fm. The result is shown in Fig. 6.

One observes that indeed the expected compensation occurs and the model prediction approaches the data points. No effort has been made to obtain a best description of the data using $`c_\mathrm{\Gamma }`$ and $`R_{hc}`$ as fit parameters.

## 6 Conclusions

Assuming a thermalized system and statistical hadronization at the phase boundary with subsequent resonance decays, we have demonstrated that our recently proposed fireball evolution scenario  leads to a reasonable description of the data on hadron ratios, if the known vacuum properties of particles and resonances are used. This provides a valuable consistency check of three pieces of information: 1) The entropy density at $`T_C=170`$ MeV can be calculated using the hadronic ensemble. The entropy density of the QGP has to agree with that number in order to produce a second order transition or a crossover, as assumed in the fireball evolution scenario. 2) The quasiparticle model predicts the entropy density at $`T_C`$, but without specifying an absolute value for $`T_C`$. 3) The essential information of the absolute value of $`T_C`$ is obtained in lattice calculations , and using this value, we find a consistent scenario. In that sense, hadron ratios provide information on the EoS of the quark-gluon plasma.

On the other hand, in-medium modifications of hadron properties are mandatory if one tries to explain the dilepton invariant mass spectrum measured by the CERES collaboration . Specifically, the $`\rho `$ channel requires strong broadening. Thermal effective field theory calculations indicate the presence of such effects for other particles as well .

In a schematic investigation, we have demonstrated that such strong broadening of resonances or mass reductions are not in line with the measured data, and no combination of these two effects can be. Such changes in particle properties are also inconsistent with the phase transition scenario used here, as they destroy the match of the entropy density and require in general a different choice of $`T_C`$. This is also apparent from the fact that the *absolute* numbers of particle production do not agree with measurements.

However, if one also consideres the effect of a screened binding potential in a thermal environment, one can make the assumption that the core radius of hadrons grows. This effect can then compensate the effects of both increased decay widths and mass reductions to some degree, as demonstrated in Fig. 6. This result is reassuring, as it allows to reconcile the in-medium modifications observed in the dilepton data with the statistical model description of measured hadron ratios in a consistent model framework.

## Acknowledgements

I would like to thank W. Weise, A. Polleri, R. A. Schneider, P. Braun-Munzinger and J. Stachel for interesting discussions and helpful comments.


text:
# Relativistic Turbulence: A Long Way from Preheating to Equilibrium

## Introduction.

The dynamics of equilibration and thermalization of field theories is of interest for various reasons. In high-energy physics understanding of these processes is crucial for applications to heavy ion collisions and to reheating of the early universe after inflation. Inflation solves the flatness and the horizon problems of the standard big bang cosmology and provides a calculable mechanism by which initial density perturbations were generated . At the end of inflation the Universe was in a vacuum-like state. In the process of decay of this state and subsequent thermalization (reheating) the matter content of the universe is created. It was realized recently that the initial stage of reheating, dubbed preheating , is a fast, explosive process. This initial stage by now is well understood . Strong and fast amplification of fluctuation fields at low momenta may lead to various interesting physical effects, like non-thermal phase transitions , peculiar baryogenesis , generation of high-frequency gravitational waves , etc.

Understanding of the subsequent stages of reheating and thermalization processes and calculation of the final equilibrium temperature is important for various applications, most notably baryogenesis and the problem of over-abundant gravitino production in supergravity models . Thermalization of field theories was discussed already, see e.g. Refs. . However, at present the process of thermalization after preheating is still far away from being well understood and developed. The problem is that at the preheating stage the occupation numbers are very large, of order of the inverse coupling constant. In addition, in many models the zero mode does not decay completely. Therefore, a simple kinetic approach is not applicable.

Fortunately, the description in terms of classical field theory is valid in this situation , and the process of preheating, as well as subsequent thermalization, can be studied on a lattice. In this paper we adopt this approach. Our goal is to integrate the system on a lattice sufficiently accurately and sufficiently far in time to be able to see generic features, and possibly to the stage, at which the kinetic description becomes a good approximation scheme. Lattice studies of thermalization, similar to ours, were done in Ref. . Several generic rules of thermalization were formulated, like the early equipartition of energy between coupled fields. However, the problem is very complicated and there are other unanswered important questions like what is the final thermalization temperature, at what stage the kinetic description becomes valid, what is the functional form of particle distributions during the thermalization stage, etc.

For our study we use a higher accuracy, improved version of the LATTICEEASY code . We show that the distribution functions follow a self-similar evolution related to the turbulent transport of wave energy. This property enables us to estimate the physical reheating temperature, which turns out to be very low. The concept should be rather model independent since typical ranges of particle momenta at preheating and in thermal equilibrium are widely separated. However, in this letter we will restrict our numerical integration (but not the discussion) to the “minimal” inflationary model, the massless $`\lambda \mathrm{\Phi }^4`$-theory.

## The Model.

With conformal coupling to gravity and after a rescaling of the field, $`\phi \mathrm{\Phi }a`$, where $`a(t)`$ is the cosmological scale factor, the equation of motion in comoving coordinates describes a $`\phi ^4`$-theory in Minkowski space-time,

$$\mathrm{}\phi +\lambda \phi ^3=0.$$
(1)
At the end of inflation the field is homogeneous, $`\phi =\phi _0(t)`$. Later on fluctuations develop, but the homogeneous component of the field, which corresponds to the zero momentum in the Fourier decomposition, may be dynamically important and is referred to as the “zero-mode.” In such situations it is convenient to make a further rescaling of the field, $`\varphi \phi /\phi _0(t_0)`$, and of the space-time coordinates, $`x^\mu \sqrt{\lambda }\phi _0(t_0)x^\mu `$, which transforms the equation of motion (1) into dimensionless and parameter free form,

$$\mathrm{}\varphi +\varphi ^3=0.$$
(2)
Here $`t_0`$ corresponds to the initial moment of time (end of inflation), and in what follows we denote dimensionless time as $`\tau `$. With this rescaling the initial condition for the zero-mode oscillations is $`\varphi _0(\tau _0)=1`$. All model dependence on the coupling constant $`\lambda `$ and on the initial amplitude of the field oscillations now is encoded in the initial conditions for the small (vacuum) fluctuations of the field with non-zero momenta . The physical normalization of the inflationary model corresponds to a dimensionful initial amplitude of $`\phi _0(t_0)0.3M_{\mathrm{Pl}}`$ and a coupling constant $`\lambda 10^{13}`$ . The re-parametrization property of the system allows to chose a larger value of $`\lambda `$ for numerical simulations. We have used $`\lambda =10^8`$.

## Numerical Procedure and Results.

We use a 3-D cubic lattice with periodic boundary conditions. The finite-differences scheme that was used is 2nd order in time and 4-th order in space. The results displayed here are taken from a simulation with $`256^3`$ sites and a physical box size $`L=14\pi `$. With this box size the infrared modes which belong to the resonance band are still well represented, while the ultraviolet lattice cut-off is sufficiently far away from the occupied modes, such that the particle spectra are not distorted even at late times. We have studied the dependence of our results on the lattice- and the box size to avoid lattice artifacts. Various quantities are measured and monitored both in configuration space (zero mode, $`\varphi _0\varphi `$, and the variance, $`var(\varphi )\varphi ^2\varphi _0^2`$) and in the Fourier space. Using fourier transformed fields we first define the wave amplitudes (which correspond to annihilation operators in the quantum problem),

$$a(\stackrel{}{k})\frac{\omega _k\varphi _\stackrel{}{k}+i\dot{\varphi }_\stackrel{}{k}}{(2\pi )^{3/2}\sqrt{2\omega _k}}.$$
(3)
The effective frequency $`\omega _k\sqrt{k_\mathrm{D}^2+m_{eff}^2}`$ is determined by the effective mass $`m_{eff}^2=3\lambda \varphi ^2`$ and the inverse $`k_\mathrm{D}^2`$ of the lattice Laplacian. In our numerical scheme the latter is given by

$$k_\mathrm{D}^2=b^2\underset{i\{1,2,3\}}{}\left(\frac{5}{2}\frac{8}{3}\mathrm{cos}(bk_i)+\frac{1}{6}\mathrm{cos}(2bk_i)\right).$$
(4)
Here $`b=2\pi /L=1/7`$ is the lattice constant. Making use of $`a_k`$, we calculate various correlators, $`n(k)a^{}a`$, $`\sigma (k)aa`$, $`a^{}a^{}aa`$, etc. The first one, which corresponds to the particle occupation numbers, is of prime interest.

We begin the discussion of our numerical results with the evolution of the zero-mode and the variance of the field, which are shown in Fig. 1. Initially we see an exponentially fast transfer of the zero-mode energy into fluctuations during preheating (up to $`\tau 300`$). It is followed by a long and slow relaxation phase. In this regime ($`\tau >1500`$) the amplitude of the zero mode oscillations decreases according to $`\tau ^{1/3}`$, the variance of the field (averaged over high-frequency oscillations) drops according to $`\tau ^{2/5}`$. This is consistent with previous results . In addition we find that in this regime the zero-mode is in a non-trivial dynamical equilibrium with the bath of highly occupied modes: when the zero-mode is artificially removed, it is recreated on a short time-scale (Bose condensation).

At early times the distribution functions of particles over momenta, see Fig. 2, have peaky structure. The first peak which corresponds to the parametric resonance is initially at the theoretically predicted value of $`k1.27`$ . Later ($`\tau >1500`$) the spectra become smooth and at small $`k`$ approach a power-law, $`n_kk^s`$, where $`s`$ fluctuates in the range of $`1.51.7`$, depending on time and the range of $`k`$ where it is fitted. This power law clearly differs from the the classical thermal equilibrium, $`n_k\omega _k^1`$. It is followed by the exponential cut-off, whose position monotenously shifts with time towards higher $`k`$. Pumping of energy from the zero mode stays effective all the times (note a small bump in the particle distributions in Fig. 2 at $`k1`$). It corresponds to the annihilation of four condensed particles into two quanta. Rescattering of two particles into two particles is also effective. One of the two can belong to the zero-mode condensate either in initial or in the final two particle state. We also can see in Fig. 2 that in the power-law region $`n_k`$ is a function of $`k/\overline{\varphi }_0`$ only, where $`\overline{\varphi }_0(\tau )`$ represents the amplitude of the zero-mode at time $`\tau `$. This effect can be related to the above described dynamical equilibrium between zero-mode and the bath of particles. Indeed, going from Eq. (1) to Eq. (2) we can rescale by the current amplitude of the zero-mode.

The picture presented in Fig. 2 at late times resembles stationary Kolmogorov turbulence. It appears as such due to rescaling of momenta by the amplitude of the zero-mode, but in fact in the present model the turbulence can not be stationary because the amplitude of the zero-mode (i.e. the strength of the source of turbulence) decreases. Further examination of Fig 2 suggests that the evolution of particle spectra may be self-similar. We have tried therefore the following anzatz

$$n(k,\tau )=\tau ^qn_0(k\tau ^p).$$
(5)
Spectra rescaled at several moments of time by the relation inverse to Eq. (5) are shown in Fig. 3. We have found that the evolution is indeed self-similar with $`q3.5p`$ and $`p1/5`$.

## Discussion.

Here we discuss the question whether a simple kinetic theory gives predictions for turbulence and self-similarity exponents in agreement with the lattice calculations. Our lattice study of higher order correlators, like $`a^{}a^{}aa`$, shows that the field distribution is very close to Gaussian, see also . This facilitates the use of the kinetic approach. On the other hand we have found that the magnitude of $`\sigma (k)`$ is of order of a few percent compared to $`n(k)`$, and it is even larger in the region of resonant momenta. This means that the strict kinetic approach should include $`\sigma (k)`$. Nevertheless we neglect these effects and write a kinetic equation in a simple form $`\dot{n_k}=I_k`$, where the collision integral for a $`m`$-particle interaction is given by

$$I_k=𝑑\mathrm{\Omega }_kU_kF[n].$$
(6)
In $`d`$ spatial dimensions the integration measure $`d\mathrm{\Omega }_k`$ is given by $`m1`$ integrations over $`d`$-dimensional Fourier space. We include in it the energy-momentum conservation $`\delta `$-functions. But we do not include there the relativistic $`1/\omega (k_i)`$ “on-shell” factors, which instead appear in the “matrix element” of the corresponding process, $`U_k`$. This will make the discussion of relativistic and non-relativistic cases uniform. The function $`F[n]`$ is a sum of products of the type $`n_{k_j}^1_{i=1}^mn_{k_i}`$, where $`j\{1,\mathrm{},m\}`$ with appropriate signs and permutations of indices for incoming and outgoing particles. All dynamical aspects of turbulence follow from the scaling properties of the system . Let $`\omega _k,n_k`$ and $`U_k`$ have defined weights under a $`\xi `$-rescaling of Fourier-space,

$`\omega (\xi k_i)`$ $`=`$ $`\xi ^\alpha \omega (k_i),`$ (7)
$`U(\xi k_1,\mathrm{},\xi k_m)`$ $`=`$ $`\xi ^\beta U(k_1,\mathrm{},k_m),`$ (8)
$`n(\xi k_i)`$ $`=`$ $`\xi ^\gamma n(k_i).`$ (9)
The weight of the full collision integral under this re-parametrization is

$$I_{\xi k}=\xi ^{d(m2)\alpha +\beta +(m1)\gamma }I_k.$$
(10)
It follows that the stationary turbulence with constant energy flux over momentum space is characterised by a power-law distribution function, $`n_kk^s`$, where $`s=d+\beta /(m1)`$. The scaling properties also give the exponents of the self-similar distribution, Eq. (5). Assuming energy conservation in particles and with $`\xi =\tau ^p`$ we find $`q=4p`$ and $`p=1/((m1)\alpha \beta )`$. For stationary turbulence we find that $`p`$ should be $`(m1)`$ times larger.

For a $`\lambda \varphi ^4`$-theory in three spatial dimensions and four-particle interaction we have $`m=4`$, $`\beta =4\alpha `$ and $`\alpha =1`$. In this case $`s=5/3`$ and $`p=1/7`$. For three-particle interaction (the fourth particle belongs to the condensate in this case and the matrix element contains an additional factor of $`\overline{\varphi }_0^2`$) we find $`s=3/2`$ and a smaller value for $`p`$ compared to the previous case. We can not distinguish between $`5/3`$ and $`3/2`$ for s in our numerical integrations, $`s`$ rather fluctuates between these two numbers, while $`1/7`$ for $`p`$ gives a fit to the data not as good as displayed in Fig. 3. However, during the integration time the energy in particles is neither conserving, nor there is a stationary source of energy. Namely, starting from the time at which the solution becomes self-similar, $`\tau 3000`$, to the end of our integration, the energy influx from the zero mode to particles is about 20 %. Correcting for this energy influx we find $`q3.5p`$ and $`p1/6`$. This should be considered as satisfactory agreement given the simplifications which were made.

## Equilibration time and temperature.

At late times the influence of the zero-mode should become negligible, but we still may expect the self-similar character of the evolution. Solution Eq. (5) with $`p=1/7`$ should be valid in this case. This allows us to find the time needed to reach equilibrium. Indeed, the classical evolution will continue until the occupation numbers in the region of the peak in Fig. 3 will become of order one. At this time quantum effects become important and the distribution relaxes to thermal. Values of momenta were this happens are $`k_{\mathrm{max}}\lambda ^{1/4}\phi _0(\tau _0)`$. On the other hand the initial distribution is centred around $`k_0\lambda ^{1/2}\phi _0(\tau _0)`$ and moves to ultraviolet according to Eq. (5) as $`k_0\tau ^p`$. It follows that the time to reach equilibrium is $`\tau \lambda ^{7/4}10^{23}`$, where in the second equality we assumed the normalization to the inflationary model. For the reheating temperature we find, rotating back from the conformal reference frame, $`T_Rk_{max}/a(\tau )\lambda ^2\phi _0(t_0)10^{26}M_{\mathrm{Pl}}100`$ eV, where for the conformal scale factor we have used $`a(\tau )=\tau `$.

## Conclusions.

Reheating after preheating appears to be a rather slow process. Although the “effective temperature” measured at low momentum modes during preheating may be high, in the model we have considered the resulting true temperature is parametrically the same as what could have been obtained in “naive” perturbation theory. Namely, equating the rate of scattering in thermal equilibrium to the Hubble expansion rate one obtains $`T\lambda ^2M_{\mathrm{Pl}}`$ in this model. We anticipate this result should be applicable to more realistic models of inflation. Note that realistic models involve many fields and interactions and larger coupling constants will determine the true temperature.

## Acknowledgements.

We thank A. Riotto, C. Schmid and D. Semikoz for many useful discussions during various stages of this project. R.M. thanks the Tomalla Foundation for financial support.