text:
# Testing a hadronic rescattering model for RHIC collisions using the subdivision method

## 1 Introduction

Recent hadronic rescattering model calculations have been shown to agree reasonably well with experimental results obtained in the first physics run for the Relativistic Heavy Ion Collider (RHIC) with $`\sqrt{s}`$ = 130 GeV Au + Au collisions . More specifically, these calculations provide good representations of the data for 1) the particle mass dependence of the $`p_T`$ distribution slope parameters (i.e. radial flow), 2) the $`p_T`$, particle mass, and pseudorapidity dependences of the elliptic flow, and 3) the $`p_T`$ and centrality dependences of the pion Hanbury-Brown-Twiss (HBT) measurements. This is a remarkable accomplishment for a single model.

Many hadrons are initially produced in a relatively small volume in RHIC-type collisions, particularly at early times in the interaction. It is a challenging task for a Monte Carlo calculation of this type to deal accurately with binary collisions between particles in such a large particle density environment without introducing numerical artifacts which may affect the results. One such undesirable artifact that can occur when the interaction range between two particles is much greater than the scattering mean-free-path is non-causality of the collisions . The method of subdivision can be used to minimize these artifacts associated with high particle density.

The goal of this work is to apply the method of subdivision to the present RHIC-energy rescattering calculations to test whether the results from these calculations without subdivision are valid. Calculations will be compared with and without subdivision for the pion observables 1) $`(1/p_T)dN/dp_T`$, 2) elliptic flow vs. $`p_T`$ and $`\eta `$, and 3) HBT vs. $`p_T`$, at a collision impact parameter of 8 fm. Section II will describe the calculational methods used and Section III with give the results of the study.

## 2 Calculational Methods

### 2.1 Hadronic rescattering calculation

A brief description of the rescattering model calculational method is given below. The method used is similar to that used in previous calculations for lower CERN Super Proton Synchrotron (SPS) energies . Rescattering is simulated with a semi-classical Monte Carlo calculation which assumes strong binary collisions between hadrons. The Monte Carlo calculation is carried out in three stages: 1) initialization and hadronization, 2) rescattering and freeze out, and 3) calculation of experimental observables. Relativistic kinematics is used throughout. All calculations are made to simulate RHIC-energy Au+Au collisions in order to compare with the results of the $`\sqrt{s}`$ = 130 GeV RHIC data.

The hadronization model employs simple parameterizations to describe the initial momenta and space-time of the hadrons similar to that used by Herrmann and Bertsch . The initial momenta are assumed to follow a thermal transverse (perpendicular to the beam direction) momentum distribution for all particles,

$$(1/m_T)dN/dm_T=Cm_T/[\mathrm{exp}(m_T/T)\pm 1]$$
(1)
where $`m_T=\sqrt{p_{T}^{}{}_{}{}^{2}+m_{0}^{}{}_{}{}^{2}}`$ is the transverse mass, $`p_T`$ is the transverse momentum, $`m_0`$ is the particle rest mass, $`C`$ is a normalization constant, and $`T`$ is the initial “temperature parameter” of the system, and a gaussian rapidity distribution for mesons,

$$dN/dy=D\mathrm{exp}[(yy_0)^2/(2\sigma _{y}^{}{}_{}{}^{2})]$$
(2)
where $`y=0.5\mathrm{ln}[(E+p_z)/(Ep_z)]`$ is the rapidity, $`E`$ is the particle energy, $`p_z`$ is the longitudinal (along the beam direction) momentum, $`D`$ is a normalization constant, $`y_0`$ is the central rapidity value (mid-rapidity), and $`\sigma _y`$ is the rapidity width. Two rapidity distributions for baryons have been tried: 1) flat and then falling off near beam rapidity and 2) peaked at central rapidity and falling off until beam rapidity. Both baryon distributions give about the same results. The initial space-time of the hadrons for $`b=0`$ fm (i.e. zero impact parameter or central collisions) is parameterized as having cylindrical symmetry with respect to the beam axis. The transverse particle density dependence is assumed to be that of a projected uniform sphere of radius equal to the projectile radius, $`R`$ ($`R=r_0A^{1/3}`$, where $`r_0=1.12`$ fm and $`A`$ is the atomic mass number of the projectile). For $`b>0`$ (non-central collisions) the transverse particle density is that of overlapping projected spheres whose centers are separated by a distance b. The particle multiplicities for $`b>0`$ are scaled from the $`b=0`$ values by the ratio of the overlap volume to the volume of the projectile. The longitudinal particle hadronization position ($`z_{had}`$) and time ($`t_{had}`$) are determined by the relativistic equations ,

$$z_{had}=\tau _{had}\mathrm{sinh}y;t_{had}=\tau _{had}\mathrm{cosh}y$$
(3)
where $`y`$ is the particle rapidity and $`\tau _{had}`$ is the hadronization proper time. Thus, apart from particle multiplicities, the hadronization model has three free parameters to extract from experiment: $`\sigma _y`$, $`T`$ and $`\tau _{had}`$. The hadrons included in the calculation are pions, kaons, nucleons and lambdas ($`\pi `$, K, N, and $`\mathrm{\Lambda }`$), and the $`\rho `$, $`\omega `$, $`\eta `$, $`\eta ^{}`$, $`\varphi `$, $`\mathrm{\Delta }`$, and $`K^{}`$ resonances. For simplicity, the calculation is isospin averaged (e.g. no distinction is made among a $`\pi ^+`$, $`\pi ^0`$, and $`\pi ^{}`$). Resonances are present at hadronization and also can be produced as a result of rescattering. Initial resonance multiplicity fractions are taken from Herrmann and Bertsch , who extracted results from the HELIOS experiment . The initial resonance fractions used in the present calculations are: $`\eta /\pi =0.05`$, $`\rho /\pi =0.1`$, $`\rho /\omega =3`$, $`\varphi /(\rho +\omega )=0.12`$, $`\eta ^{}/\eta =K^{}/\omega =1`$ and, for simplicity, $`\mathrm{\Delta }/N=0`$.

The second stage in the calculation is rescattering which finishes with the freeze out and decay of all particles. Starting from the initial stage ($`t=0`$ fm/c), the positions of all particles are allowed to evolve in time in small time steps ($`dt=0.1`$ fm/c) according to their initial momenta. At each time step each particle is checked to see a) if it decays, and b) if it is sufficiently close to another particle to scatter with it. Isospin-averaged s-wave and p-wave cross sections for meson scattering are obtained from Prakash et al. . The calculation is carried out to 100 fm/c, although most of the rescattering finishes by about 30 fm/c. The rescattering calculation is described in more detail elsewhere .

Calculations are carried out assuming initial parameter values and particle multiplicities for each type of particle. In the last stage of the calculation, the freeze-out and decay momenta and space-times are used to produce observables such as pion, kaon, and nucleon multiplicities and transverse momentum and rapidity distributions. The values of the initial parameters of the calculation and multiplicities are constrained to give observables which agree with available measured hadronic observables. As a cross-check on this, the total kinetic energy from the calculation is determined and compared with the RHIC center of mass energy of $`\sqrt{s}=130`$ GeV to see that they are in reasonable agreement. Particle multiplicities were estimated from the charged hadron multiplicity measurements of the RHIC PHOBOS experiment . Calculations were carried out using isospin-summed events containing at freezeout for central collisions ($`b=0`$ fm) about 5000 pions, 500 kaons, and 650 nucleons ($`\mathrm{\Lambda }`$’s were decayed). The hadronization model parameters used were $`T=300`$ MeV, $`\sigma _y`$=2.4, and $`\tau _{had}`$=1 fm/c. It is interesing to note that the same value of $`\tau _{had}`$ was required in a previous rescattering calculation to successfully describe results from SPS Pb+Pb collisions .

### 2.2 Subdivision method

The method of subdivision is based on the invariance of Monte Carlo particle-scattering calculations for a simultaneous decrease of the scattering cross sections by some factor, $`l`$, and increase of the particle density by $`l`$. As $`l`$ becomes sufficiently large, non-causal artifacts become insignificant. The present rescattering calculation will be tested comparing pion observables for no-subdivision, i.e. $`l=1`$, with subdivision of $`l=5`$. A subdivision of $`l=5`$ was chosen since when applied to other Monte Carlo scattering calculations this value has been shown to produce large changes in observables when non-causal artifacts are present. Since the particle density increase is accomplished by increasing the particle number by a factor $`l`$, the computer CPU time taken per event increases by a factor $`l^2`$, effectively “slowing down” a subdivision study by a factor $`l`$ compared with not using subdivision. For the present study, 1640 events and 124 events were generated for the $`l=1`$ and $`l=5`$ samples, respectively (the statistical value of the $`l=5`$ events being 5 times greater per event). The CPU time taken to generate the $`l=5`$ sample of events was 420 CPU-hours on a 1 GHz PC processor.

## 3 Results

Figures 1-4 show comparisons of $`l=1`$ with $`l=5`$ hadron rescattering model calculations for pion observables. Descriptions of how the observables are extracted from the rescattering calculation are given elsewhere. All calculations are carried out for an impact parameter of 8 fm to simulate a medium non-central collision which should result in significant elliptic flow for the purposes of the present test. Statistical errors are shown either as error bars or are of the order of the marker size when error bars are not shown.

Figure 1 compares the $`p_T`$ distributions for a rapidity cut of $`2<y<2`$. As seen, the $`l=1`$ and $`l=5`$ calculations agree within statistical errors, although there is a slight decreasing trend for $`l=5`$ compared with $`l=1`$ for $`p_T>3`$ GeV/c.

Figures 2 and 3 compare elliptic flow ($`V_2`$). Figure 2 shows $`V_2`$ vs. $`p_T`$ for a rapidity cut of $`2<y<2`$ and Figure 3 shows $`V_2`$ vs. $`\eta `$ for a $`p_T`$ cut of $`0<p_T<3`$ GeV/c. For both cases, the $`l=1`$ and $`l=5`$ calculations are seen to agree within statistical errors, reproducing the “flattening” in $`V_2`$ at large $`p_T`$ and the “peaked” dependence of $`V_2`$ on $`\eta `$.

Finally, Figure 4 compares HBT results vs. $`p_T`$ for a rapidity cut of $`2<y<2`$. Each result is calculated in a $`p_T`$ bin of 0.2 GeV/c. A three-dimensional fit is made to extract the pion source parameters $`R_{Tside}`$, $`R_{Tout}`$, $`R_{long}`$, and $`\lambda `$. As seen in this figure, the $`l=1`$ and $`l=5`$ calculations agree within statistical errors for all parameters at all $`p_T`$ bins, reproducing the decreasing trend of the $`R`$-parameters and increasing trend of the $`\lambda `$-parameter with increasing $`p_T`$.

It can be concluded from Figures 1-4 that non-causal artifacts of the type which can be removed by the subdivision method are not present in the hadronic rescattering code used in this study. It is possible to speculate on why this is the case. There are three main features of the code which might tend to reduce these artifacts: 1) individual particles are allowed to scatter only once per time step, 2) a “scattering time” of two time steps is defined during which particles that have scattered are not allowed to rescatter, and 3) once two particles have scattered with each other, they are not allowed to scatter with each other again in the calculation. The present study was performed for non-central collisions ($`b=8`$ fm), but it is expected that the same conclusion, i.e. that subdivision does not significantly effect the results of the calculation, would also be obtained for central collisions. This is because the particle density is not very different for central collisions in the model than for mid-peripheral collisions since the particle multiplicities are scaled by the overlap volume for $`b>0`$, as mentioned earlier.

In summary, the subdivision method has been used to test the validity of the present rescattering model calculations. It is found that no appreciable changes are seen comparing $`l=1`$ and $`l=5`$ calculations, thus strengthing the confidence in the results presented previously from this rescattering model for RHIC.


format: text

text:
# Inequalities for solutions to some nonlinear equations key words: inequalities, nonlinear equations, dynamical systems method, evolution equations AMS subject classification: 37C35, 37L05, 37N30, 47H15, 58C15, 65J15,

## 1 Introduction

In this paper Theorem 1 gives a method for proving existence of the solution to nonlinear equation $`F(u)=0`$ and for computing it. Our method (dynamical systems method: DSM) consists of solving a suitable Cauchy problem which has a global solution $`u(t)`$ such that $`y:=u(\mathrm{})`$ does exist and $`F(y)=0`$. Theorem 1 generalizes Theorem 1 in  and its proof is based on the idea in . The result of Example 1 in Section 3 below was obtained in . The result of Example 5 was essentially obtained in , where the existence of the solution was assumed, while in Theorem 1 below the existence of the solution is proved. In Example 6 some of the results from  are obtained. In Remark 2.4 in  some of the examples we discuss in Section 3 were mentioned. In  DSM was developed for solving ill-posed operator equations with not necessarily monotone operators and for constructing convergent iterative methods for their solution.

Let $`F`$ be a nonlinear Frechet differentiable map in a real Hilbert space. Consider the equation

$$F(u)=0.$$
$`(1.1)`$
Let $`\mathrm{\Phi }(t,u)`$ be a map, continuous in the norm of $`H`$ and Lipschitz with respect to $`u`$ in the ball $`B:=\{u:uu_0R,uH\}`$. Weaker conditions, which guarantee local existence and uniqueness of the solution to (1.6) below, would suffice. Assume that:

$$(F^{}(u)\mathrm{\Phi }(t,u),F(u))g_1(t)F(u)^2uB,$$
$`(1.2)`$
and

$$\mathrm{\Phi }(t,u)g_2(t)F(u)uB,$$
$`(1.3)`$
where $`g_j,j=1,2,`$ are positive functions on $`R_+:=[0,\mathrm{})`$, $`g_2`$ is continuous, $`g_1L_{loc}^1(R_+)`$,

$$_0^{\mathrm{}}g_1𝑑t=+\mathrm{},$$
$`(1.4)`$
and

$$G(t):=g_2(t)\mathrm{exp}(_0^tg_1𝑑s)L^1(R_+).$$
$`(1.5)`$
Remark:  Sometimes the assumption (1.3) can be used in the following modified form:

$$\mathrm{\Phi }(t,u)g_2(t)F(u)^buB,$$
$`(1.3^{})`$
where $`b>0`$ is a constant. The statements and proofs of Theorems 1-3 in Sections 1 and 2 can be easily adjusted to this assumption.

Consider the following Cauchy problem:

$$\dot{u}=\mathrm{\Phi }(t,u),u(0)=u_0,\dot{u}:=\frac{du}{dt}.$$
$`(1.6)`$
Finally, assume that

$$F(u_0)_0^{\mathrm{}}G(t)𝑑tR.$$
$`(1.7)`$
The above assumptions on $`F,\mathrm{\Phi },g_1`$ and $`g_2`$ hold throughout and are not repeated in the statement of Theorem 1 below. By global solution to (1.6) we mean a solution defined for all $`t0`$.

Theorem 1. Under the above assumptions problem (1.6) has a global solution $`u(t)`$, there exists the strong limit $`y:=lim_t\mathrm{}u(t)=u(\mathrm{})`$, $`F(y)=0`$, $`u(t)B`$ for all $`t0`$, and the following inequalities hold:

$$u(t)yF(u_0)_t^{\mathrm{}}G(x)𝑑x,$$
$`(1.8)`$
and

$$F(u(t))F(u_0)\mathrm{exp}(_0^tg_1(x)𝑑x).$$
$`(1.9)`$
Theorem 1 generalizes Theorem 1 in .

In Section 2 proof of Theorem 1 is given, and two other theorems are proved, and in Section 3 examples of applications are presented. In Section 4 a linear equation and in Section 5 a nonlinear equation are discussed.

## 2 Proof of Theorem 1 and additional results.

The assumptions about $`\mathrm{\Phi }`$ imply local existence and uniqueness of the solution $`u(t)`$ to (1.6). To prove global existence of $`u`$, it is sufficient to prove a uniform with respect to $`t`$ bound on $`u(t)`$. Indeed, if the maximal interval of the existence of $`u(t)`$ is finite, say $`[0,T)`$, and $`\mathrm{\Phi }(t,u)`$ is locally Lipschitz with respect to $`u`$, then $`u(t)\mathrm{}`$ as $`tT`$. Let $`g(t):=F(u(t))`$. Since $`H`$ is real, one uses (1.6) and (1.2) to get $`g\dot{g}=(F^{}(u)\dot{u},F)g_1(t)g^2`$, so $`\dot{g}g_1(t)g`$, and integrating one gets (1.9), because $`g(0)=F(u_0)`$. Using (1.3), (1.6) and (1.9), one gets:

$$u(t)u(s)g(0)_s^tG(x)𝑑x,G(x):=g_2(x)\mathrm{exp}(_0^xg_1(z)𝑑z).$$
$`(2.1)`$
Because $`GL^1(R_+)`$, it follows from (2.1) that the limit $`y:=lim_t\mathrm{}u(t)=u(\mathrm{})`$ exists, and $`yB`$ by (1.7). From (1.9) and the continuity of $`F`$ one gets $`F(y)=0`$, so $`y`$ solves (1.1). Taking $`t\mathrm{}`$ and setting $`s=t`$ in (2.1) yields estimate (1.8). The inclusion $`u(t)B`$ for all $`t0`$ follows from (2.1) and (1.7). Theorem 1 is proved. $`\mathrm{}`$

If condition (1.2) is replaced by

$$(F^{}\mathrm{\Phi },F)g_1(t)F^a,0<a<2,$$
$`(2.2)`$
then the proof of Theorem 1 yields $`g^{1a}\dot{g}g_1(t)`$, so

$$0g(t)[g^{2a}(0)(2a)_0^tg_1(s)𝑑s]^{\frac{1}{2a}}.$$
$`(2.3)`$
If (1.4) holds, then (2.3) implies $`g(t)=0`$ for all $`tT`$, where $`T`$ is defined by the equation:

$$g^{2a}(0)(2a)_0^Tg_1(s)𝑑s=0.$$
$`(2.4)`$
Thus $`F(u(t))=0`$ for $`tT`$. So, by (1.3), $`\mathrm{\Phi }=0`$ for $`tT`$. Thus, by (1.6), $`u(t)=u(T)`$ for $`tT`$. Therefore $`y:=u(T)`$ solves equation (1.1), $`F(y)=0`$, and $`u(T)u(0)F(u_0)_0^Tg_2𝑑s`$. If $`F(u_0)_0^Tg_2𝑑sR`$, then $`u(t)B`$ for all $`t0`$. We have proved:

Theorem 2. If (1.2) is replaced by (2.2), (1.4) holds, and $`F(u_0)_0^Tg_2𝑑sR`$, where $`T`$ is defined by (2.4), then equation (1.1) has a solution in $`B=\{u:uu_0R\}`$, the solution $`u(t)`$ to (1.6) exists for all $`t>0`$, $`u(t)B`$, $`u(t):=y`$ for $`tT`$, and $`F(y)=0`$, $`yB`$.

If (2.2) holds with $`a>2`$, and (1.4) holds, then a similar calculation yields:

$$0g(t)[g^{(a2)}(0)+(a2)_0^tg_1(s)𝑑s]^{\frac{1}{2a}}:=h(t)0t\mathrm{},$$
$`(2.5)`$
because of (1.4). Assume that

$$_0^{\mathrm{}}g_2(s)h(s)𝑑sR.$$
$`(2.6)`$
Then (1.3) and (1.6) yield $`u(t)u(0)R`$, $`u(t)u(\mathrm{})_t^{\mathrm{}}g_2(s)h(s)𝑑s0`$ as $`t\mathrm{}`$. Therefore an analog of Theorem 1 is obtained:

Theorem 3. If (2.2) holds with $`a>2`$, (1.4) and (2.6) hold, then the solution $`u(t)`$ to (1.6) exists for all $`t>0`$, $`u(t)B`$, there exists $`u(\mathrm{}):=y`$ , and $`F(y)=0`$, $`yB`$.

## 3 Applications.

If $`g_j=c_j,j=1,2,`$ and $`c_j>0`$ are constants, then (1.4) and (1.5) hold, $`_0^{\mathrm{}}G𝑑x=c_2/c_1`$, so condition (1.7) is:

$$\frac{c_2}{c_1}F(u_0)R,c_j=c_j(R),j=1,2.$$
$`(3.1)`$
Let us give examples of applications of Theorem 1 using its simplified version with $`g_j=c_j>0,j=1,2`$.

Example 1. Continuous Newton-type method.

Let $`\mathrm{\Phi }=[F^{}(u)]^1F(u)`$, and assume

$$[F^{}(u)]^1m_1=m_1(R),uB.$$
$`(3.2)`$
Assumption (3.2) holds in all the examples below. It implies ”well-posedness” of equation (1.1). The case of ill-posed equation (1.1), when (3.2) is not valid, was studied in  for linear selfadjoint nonnegative-definite operator, in  for nonlinear monotone hemicontinuous operator, and in  for some other cases. Then $`c_1=1`$, $`c_2=m_1`$, $`\mathrm{\Phi }`$ is locally Lipschitz if, for example, one assumes

$$F^{\prime \prime }(u)M_2,uB,$$
$`(3.3)`$
where $`M_2=M_2(R)`$ is a positive constant, and (3.1) is:

$$m_1(R)F(u_0)R.$$
$`(3.4)`$
In the examples below condition (3.3) is assumed and not repeated. As we have mentioned in Section 1, this condition can be weakened.

Conclusion 1: By Theorem 1, inequality (3.4) implies existence of a solution $`y`$ to (1.1) in $`B`$, global existence and uniqueness of the solution $`u(t)`$ to (1.6), convergence of $`u(t)`$ to $`y`$ as $`t\mathrm{}`$, and the error estimate (1.9). Condition (3.4) is always satisfied if equation (1.1) has a solution $`y`$ and if $`u_0`$ is chosen sufficiently close to $`y`$.

Example 2. Continuous simple iterations method.

Let $`\mathrm{\Phi }=F`$, and assume $`F^{}(u)c_1(R)>0`$ for all $`uB`$. Then $`c_2=1`$, $`c_1=c_1(R)`$, and (3.1) is:

$$[c_1(R)]^1F(u_0)R.$$
$`(3.5)`$
If this inequality holds, then Conclusion 1 holds with (3.5) replacing (3.4).

Example 3. Continuous gradient method.

Let $`\mathrm{\Phi }=[F^{}]^{}F`$ and assume (3.2). Then $`c_2=M_1(R)`$, because $`[F^{}(u)]^{}=F^{}(u)M_1(R)`$, and $`(F^{}\mathrm{\Phi },F)=[F^{}(u)]^{}F^2m_1^2F^2`$, so $`c_1=m_1^2`$, where $`m_1`$ is the constant from (3.2). Here we have used the estimates $`f=A^1AfA^1Af,AfA^1^1f`$, with $`A:=F^{}(u)`$ and $`f=F(u)`$. Estimate (3.1) is:

$$M_1m_1^2F(u_0)R.$$
$`(3.6)`$
If this inequality holds, then Conclusion 1 holds with (3.6) replacing (3.4).

Example 4. Continuous Gauss-Newton method.

Let $`\mathrm{\Phi }=([F^{}]^{}F^{})^1[F^{}]^{}F`$, and assume (3.2). Then $`c_1=1`$, $`c_2=m_1^2M_1`$, and (3.1) is:

$$M_1m_1^2F(u_0)R.$$
$`(3.7)`$
If this inequality holds, then Conclusion 1 holds with (3.7) replacing (3.4).

Example 5. Continuous modified Newton method.

$`\mathrm{\Phi }=[F^{}(u_0)]^1F(u)`$, and assume $`[F(u_0)]^1m_0`$. Then $`c_2=m_0`$. To find $`c_1`$, let us note that:

$$(F^{}\mathrm{\Phi },F)=F(u)^2((F^{}(u)F^{}(u_0))[F(u_0)]^1F,F)0.5F(u)^2,$$

provided that

$$|((F^{}(u)F^{}(u_0))[F(u_0)]^1F,F)|M_2Rm_0F(u)^20.5F(u)^2,$$

that is, $`M_2Rm_00.5`$. If $`R=(2M_2m_0)^1`$, then the last inequality becomes an equality. Choosing such $`R`$, one has $`c_2=m_0`$, $`c_1=0.5`$, and (3.1) is: $`2m_0F(u_0)(2M_2m_0)^1`$, that is,

$$4m_0^2M_2F(u_0)1.$$
$`(3.8)`$
If this inequality holds, then Conclusion 1 holds with (3.8) replacing (3.4).

Example 6. Descent methods.

Let $`\mathrm{\Phi }=\frac{f}{(f^{},h)}h`$, where $`f=f(u(t))`$ is a differentiable functional $`f:H[0,\mathrm{})`$, and $`h`$ is an element of $`H`$. One has $`\dot{f}=(f^{},\dot{u})=f`$. Thus $`f=f_0e^t`$, where $`f_0:=f(u_0)`$. Assume $`\mathrm{\Phi }c_2|f|^b,b>0`$. Then $`\dot{u}c_2|f_0|^be^{bt}`$. Therefore $`u(\mathrm{})`$ does exist, $`f(u(\mathrm{}))=0`$, and $`u(\mathrm{})u(t)ce^{bt}`$, $`c=const>0`$.

If $`h=f^{}`$, and $`f=F(u)^2`$, then $`f^{}(u)=2[F^{}]^{}(u)F(u)`$, $`\mathrm{\Phi }=\frac{f}{f^{}^2}f^{}`$, and (1.6) is a descent method. For this $`\mathrm{\Phi }`$ one has $`c_1=\frac{1}{2}`$, and $`c_2=\frac{m_1}{2}`$, where $`m_1`$ is defined in (3.2). Condition (3.1) in this example is condition (3.4).

If (3.4) holds, then Conclusion 1 holds.

We have obtained some results from . Our approach is more general than the one in , since the choices of $`f`$ and $`h`$ do not allow one, e.g., to obtain $`\mathrm{\Phi }`$ used in Example 5.

## 4 Remark about linear equations.

The following result was proved in : if equation $`Ay=f`$ in a Hilbert space has a solution $`y`$, and $`A0`$ is a linear selfadjoint operator, then the global solution $`u`$ to the regularized Cauchy problem

$$\dot{u}=Au\alpha (t)u+f,u(0)=u_0,$$
$`(4.1)`$
has a limit $`lim_t\mathrm{}u(t):=u(\mathrm{})`$, and $`A(u(\mathrm{}))=f`$. In  $`u_0H`$ is arbitrary, $`\alpha >0`$ is a continuously differentiable, monotonically decaying to zero as $`t\mathrm{}`$, function on $`R_+`$, $`_0^{\mathrm{}}\alpha 𝑑t=+\mathrm{},`$ and $`\alpha ^2\dot{\alpha }0`$ as $`t\mathrm{}`$.

Unfortunately the author of  was not aware of the paper . Some of the proofs in  are close to these in . The author thanks Dr. Ya. Alber for pointing out references  and .

If $`\alpha >0`$, $`\dot{\alpha }0,`$ and $`\alpha ^2|\dot{\alpha }|c,`$ where $`c=const`$, then $`\alpha ^1(t)\alpha ^1(0)ct`$, so $`\alpha (t)[ct+\alpha ^1(0)]^1`$ and consequently $`_0^{\mathrm{}}\alpha 𝑑t=+\mathrm{}`$ (cf ). Therefore the condition $`_0^{\mathrm{}}\alpha 𝑑t=+\mathrm{}`$ in  can be dropped.

In this Section we give a new derivation of the result in  under weaker assumptions about $`\alpha `$, and show that the regularization in (4.1) is not necessary.

First, let us prove that the regularization in (4.1) is not necessary: the result holds with $`\alpha =0`$. Below $``$ denotes strong convergence in $`H`$.

The solution to (4.1) with $`\alpha =0`$ is $`u(t)=U(t)u_0+_0^tU(ts)f𝑑s`$, where $`U(t):=\mathrm{exp}(tA)`$. If $`E_\lambda `$ is the resolution of the identity of the selfadjoint operator $`A`$, then $`U(t)u_0=_0^{\mathrm{}}e^{t\lambda }𝑑E_\lambda u_0Pu_0`$ as $`t\mathrm{}`$, where $`P`$ is the operator of the orthogonal projection on $`N`$, and $`N`$ is the null-space of $`A`$. Also $`_0^tU(ts)f𝑑s=_0^{\mathrm{}}(1e^{t\lambda })𝑑E_\lambda yyPy`$ as $`t\mathrm{}`$, by the dominated convergence theorem. Thus, $`u(\mathrm{})=yPy+Pu_0`$ and $`A(u(\mathrm{}))=f`$. $`\mathrm{}`$

Consider now the case $`0<\alpha 0`$ as $`t\mathrm{}`$. If $`h(t):=\mathrm{exp}(_0^t\alpha (s)𝑑s)`$, and $`u`$ solves (4.1), then

$$u(t)=h^1(t)U(t)u_0+h^1(t)_0^{\mathrm{}}\mathrm{exp}(t\lambda )_0^te^{s\lambda }h(s)𝑑s\lambda 𝑑E_\lambda y.$$
$`(4.2)`$
Using L’Hospital’s rule one checks that

$$\underset{t\mathrm{}}{lim}\frac{\lambda _0^te^{s\lambda }h(s)𝑑s}{e^{t\lambda }h(t)}=\underset{t\mathrm{}}{lim}\frac{\lambda e^{t\lambda }h(t)}{\lambda e^{t\lambda }h(t)+e^{t\lambda }h(t)\alpha (t)}=1\lambda >0.$$
$`(4.3)`$
From (4.2), (4.3), and the dominated convergence theorem, one gets $`u(\mathrm{})=yPy`$. The first term on the right-hand side of (4.2) tends to zero as $`t\mathrm{}`$ (even if $`Pu_00`$), if $`h(\mathrm{})=\mathrm{}`$. To apply the dominated convergence theorem, one checks that $`\frac{\lambda _0^te^{(ts)\lambda }h(s)𝑑s}{h(t)}=\frac{\lambda _0^te^{s\lambda }h(ts)𝑑s}{h(t)}1`$ for all $`t>0`$ and all $`\lambda >0`$, where the inequality $`0<h(ts)h(t)`$, valid for $`s0`$, was used. $`\mathrm{}`$

Our derivation uses less restrictive assumptions on $`\alpha `$ than in  and : we do not assume differentiability of $`\alpha `$, and the property $`lim_t\mathrm{}\alpha ^2\alpha ^{}=0`$. The property $`_0^{\mathrm{}}\alpha 𝑑t=+\mathrm{}`$, which is equivalent to $`h(\mathrm{})=\mathrm{}`$, was used above only to prove that $`lim_t\mathrm{}h^1(t)U(t)u_0=0`$. If $`_0^{\mathrm{}}\alpha 𝑑t:=q<\mathrm{}`$, then $`u(\mathrm{})=yPy+e^qPu_0`$, and $`Au(\mathrm{})=f`$, so that the basic conclusion holds without the assumption $`h(\mathrm{})=\mathrm{}`$.

Finally, let us prove a typical for ill-posed problems result: the rate of convergence $`u(t)y`$ can be as slow as one wishes, it is not uniform with respect to $`f`$. Assume $`\alpha =0`$, but the proof is essentially the same for $`0<\alpha 0`$ as $`t\mathrm{}`$. Assume that $`A>0`$ is compact, and $`A\phi _j=\lambda _j\phi _j`$, $`(\phi _j,\phi _m)=\delta _{jm}`$. Then (4.2) with $`y=y_m:=\phi _m`$ and $`u_0=0`$ yields $`u(t)=\phi _m(1e^{t\lambda _m})`$. Thus $`u(\mathrm{})=y`$, but for any fixed $`T>0`$, however large, one can find $`m`$ such that $`u(T)y_m>0.5`$, that is, convergence is not uniform with respect to $`f`$.

## 5 Remark about nonlinear equations.

In this Section we give a short and simple proof of the basic result in , and close a gap in the proof in , where it is not explained why one can apply the L’Hospital rule the second time.

The assumptions in  are: the operator $`A`$ is monotone (possibly nonlinear), defined on all of $`H`$, hemicontinuous, problem (4.1) has a unique global solution, equation $`A(y)=f`$ has a solution, $`\alpha (t)>0`$ decays monotonically to zero, $`lim_t\mathrm{}\dot{\alpha }\alpha ^2=0`$, and $`\alpha `$ is convex.

We refer below to these assumptions as A3). If A3) hold, the basic result, proved in , is the existence of $`u(\mathrm{}):=y`$, and the relation $`A(y)=f`$. In , p. 184, under the additional assumption (1.24) from , the global existence of the solution to (4.1) is proved. Actually, the assumption about global existence of the solution to (4.1) can be dropped altogether: in , p.99, it is proved that A3) (and even weaker assumptions) imply that problem (4.1) has a unique global solution.

Let us give a short proof of the basic result from . It is well known that A3) imply that the problem $`A(v_\alpha )+\alpha v_\alpha f=0,`$ for any fixed number $`\alpha >0`$, has a unique solution, there exists $`lim_{\alpha 0}v_\alpha :=y`$, $`A(y)=f,`$ and $`yz`$, for any $`z\{z:A(z)=f\}`$. Thus, for any small $`\delta >0,`$ one can find $`\alpha _\delta `$ such that $`v_\alpha y<\delta /2`$ for all $`\alpha >\alpha _\delta `$, $`lim_{\delta 0}\alpha _\delta =0`$. Let $`w:=uv_\alpha `$, where $`u`$ solves (4.1) and $`v_\alpha `$ does not depend on $`t`$. Then $`\dot{w}=[A(u)A(v_\alpha )+\alpha (t)(uv_\alpha )+(\alpha (t)\alpha )v_\alpha ]`$. Multiply this by $`w`$, use the monotonicity of $`A`$, and let $`w:=g`$. Then $`g\dot{g}\alpha (t)g^2+c|\alpha (t)\alpha |g`$, $`c=y`$. Indeed, multiply $`A(v_\alpha )+\alpha v_\alpha A(y)=0`$ by $`v_\alpha y`$ and use monotonicity of $`A`$ to get $`\alpha (v_\alpha ,v_\alpha y)0`$. Thus $`v_\alpha y`$, so $`c=y`$.

Since $`\alpha (t)`$ is convex, one has $`|\alpha (t)\alpha ||\dot{\alpha }(t)|(t_\alpha t)`$, where $`t_\alpha t`$ is defined by the equation $`\alpha =\alpha (t_\alpha )`$, $`lim_{\alpha 0}t_\alpha =\mathrm{}`$. Thus, $`g\dot{g}\alpha (t)g^2+c|\dot{\alpha }(t)|(t_\alpha t)g`$, and, taking $`u(0)=v_\alpha `$, one gets

$$g(t_\alpha )ce^{_0^{t_\alpha }\alpha (x)𝑑x}_0^{t_\alpha }e^{_0^s\alpha (x)𝑑x}|\dot{\alpha }(s)|(t_\alpha s)𝑑s.$$
$`(5.1)`$
We prove below that

$$\underset{t\mathrm{}}{lim}\alpha (t)e^{_0^t\alpha (s)𝑑s}=\mathrm{}.$$
$`(5.2)`$
This allows one to apply twice L’Hospital’s rule to the right-hand side of (5.1), and to get: $`lim_{\alpha 0}g(t_\alpha )=lim_{t_\alpha \mathrm{}}\frac{\dot{\alpha }(t_\alpha )}{\dot{\alpha }(t_\alpha )+\alpha ^2(t_\alpha )}=0`$. Now, $`u(t_\alpha )yu(t_\alpha )v_\alpha +v_\alpha y`$, $`v_\alpha y\delta /2`$, and $`u(t_\alpha )v_\alpha \delta /2`$, for sufficiently large $`t_\alpha `$. Since $`\delta >0`$ is arbitrarily small, it follows that $`lim_t\mathrm{}u(t)y=0`$.

Let us prove (5.2). From our assumptions about $`\alpha `$, it follows that for all sufficiently large $`t`$, one has $`\dot{\alpha }\alpha ^2c`$, where $`0<c<1`$, so $`\alpha (t)(c_1+t)^1b`$, where $`b:=c^1>1`$, $`c_1>0`$ is a constant, and $`e^{_0^t\alpha (s)𝑑s}(c_1+t)^b`$. Thus, (5.2) holds. $`\mathrm{}`$

If one assumes additionally that $`A`$ is Frechet differentiable, then the proof is shorter. Namely, let $`h(t):=A(u(t))+\alpha (t)u(t)f:=\psi .`$ Then $`h\dot{h}=((A^{}(u(t))+\alpha (t))\psi ,\psi )\alpha (t)h^2`$, because $`A^{}0`$ due to the monotonicity of $`A`$. Thus $`h(t)\varphi (t)`$, where $`\varphi (t):=h(0)e^{_0^t\alpha 𝑑s}`$. As we proved in Section 4, the assumptions on $`\alpha (t)`$ imply $`\alpha (t)(c_1t+c_2)^1`$, where $`c_1`$ and $`c_2`$ are positive constants, and $`c_1`$ can be chosen so that $`0<c_1<1`$, due to the assumption $`lim_t\mathrm{}\dot{\alpha }\alpha ^2=0`$. Therefore $`_0^{\mathrm{}}\varphi (t)𝑑t<\mathrm{}`$. From (4.1) one gets: $`\dot{u}\varphi (t)`$. Because $`_0^{\mathrm{}}\varphi (t)𝑑t<\mathrm{}`$, it follows that $`u(\mathrm{}):=y`$ exists, and $`u(\mathrm{})u(t)_t^{\mathrm{}}\varphi (s)𝑑s`$. Finally, $`A(y)=f`$ because $`h(\mathrm{})=0=A(y)f`$. Any choice of $`\alpha `$, for which $`_0^{\mathrm{}}\varphi (t)𝑑t<\mathrm{}`$, is sufficient for the above argument. $`\mathrm{}`$


text:
# The initial Helium content of Galactic Globular Cluster stars from the 𝑅-parameter: comparison with the CMB constraint

## 1 Introduction

Galactic Globular Cluster (GGC) stars are the oldest objects in the Galaxy, and their initial He abundance ($`Y_{GGC}`$, where $`Y`$ denotes the mass fraction of He) is supposed to be approximately equal to the primordial He abundance ($`Y_p`$) produced during the Big Bang Nucleosynthesis (BBN).

The value of $`Y_p`$ derived from spectroscopy of low-metallicity, extragalactic H II regions, appears to be still subject to systematic uncertainties (see, e.g., the discussion in Bono et a. 2002 and references therein); as an example, Olive, Steigman & Skillman (1997) have determined $`Y_p`$=0.234$`\pm `$0.002, while Izotov & Thuan (1998) obtained $`Y_p`$=0.244$`\pm `$0.002, consistent with earlier findings by Kunth & Sargent (1983).

On the other hand, recent determinations of the cosmological baryonic matter density ($`\mathrm{\Omega }_b`$) from the cosmic microwave background (CMB) power spectrum obtained by the BOOMERANG, DASI and MAXIMA experiments provide consistently (e.g. Pryke et al. 2002, Ödman et al. 2002; Sievers et al. 2002) a value $`\mathrm{\Omega }_bh^2`$=0.022$`\pm `$0.003 ($`h`$ is the Hubble constant in units of 100 $`KmMpc^1s^1`$). This baryon density coupled with standard BBN calculations (Burles, Nollett & Turner 2001) provides $`Y_p`$=0.248$`\pm `$0.001; Such an independent determination of $`Y_p`$ is close to the spectroscopic determination by Izotov & Thuan (1998).

As for the value of $`Y_{GGC}`$, empirical estimates are necessarily indirect, since He-lines are not detectable in GGC star spectra, apart from the case of hot Horizontal Branch (HB) objects, whose atmospheres are however affected by gravitational settling and radiative levitation, which strongly alter the initial chemical stratification (see, e.g., Michaud, Vauclair & Vauclair 1983; Moehler et al. 1999). $`Y_{GGC}`$ estimates make use of results from stellar evolution, taking advantage of the fact that the evolution of low mass Population II stars is affected by their initial He-content. More in detail, the so-called $`R`$-parameter (Iben 1968), defined as the number ratio of HB stars to Red Giant (RGB) stars brighter than the HB level ($`R=N_{HB}/N_{RGB}`$), can be employed in order to determine $`Y_{GGC}`$ and therefore $`Y_p`$. The basic idea behind the use of this parameter as He indicator is that, at a given metallicity, a higher initial He-content implies a brighter HB and, in turn, a lower value of $`N_{RGB}`$ ($`N_{HB}`$ is only slightly affected), with the consequent increase of $`R`$. Other parameters derived from stellar evolution can also be employed (see, e.g., the discussion in Sandquist 2000 and Zoccali et al. 2000), but they are better suited to determine relative He-abundances than absolute ones.

Following earlier analyses by Buzzoni et al. (1983) and Caputo, Martinez Roger & Paez (1987), Sandquist (2000 – S00) and Zoccali et al. (2000 - Z00) have recently estimated $`Y_{GGC}`$ by measuring the value of $`R`$ in large samples of GGCs (43 objects taken from various sources, in case of the S00 paper; 26 objects with homogeneous $`HST`$ photometry for the Z00 paper), and compared the observational values with results from stellar evolution models. In both cases a value $`Y_{GGC}`$0.20 was found, in severe disagreement with the CMB result and spectroscopic determinations for H II regions. This large discrepancy between CMB and $`R`$-parameter results casts doubt on the ability of stellar models to accurately predict the evolutionary times of these crucial phases of stellar evolution. Reasons for this low He-content inferred from the $`R`$ parameter have been ascribed to the uncertainty of the 12C$`(\alpha ,\gamma )^{16}`$O reaction rate, which is relevant during the late stages of central He-burning, but may in principle also be due, for example, to an improper treatment of the mixing in the central convective region of HB stars, which affect the HB evolutionary timescale, and hence $`N_{HB}`$ (see, e.g., Z00).

In this paper we present a new determination of $`Y_{GGC}`$ from the $`R`$-parameter, employing in the model computation the most recent determination of the 12C$`(\alpha ,\gamma )^{16}`$O reaction rate (Kunz et al. 2002) together with an improved equation of state (EOS) for the stellar matter (Irwin et al. 2002). We show that with these two improved physical inputs the discrepancy between CMB and $`R`$-parameter results almost completely disappears. This gives strong support to the accuracy of present HB stellar models and the adequacy of our convective core mixing treatment during the HB phase. In §2 we briefly discuss the observational data and the theoretical stellar models, while in §3 and §4 the results are presented and discussed. Conclusions follow in §5.

## 2 Observational data and theoretical models

We have determined $`Y_{GGC}`$ making use of our theoretical models and the observational databases presented by Z00 and S00. Z00 provide empirical $`R`$ values for 26 GGCs spanning all the relevant metallicity range. The number of RGB stars is computed starting from the level of the observed V magnitude of the Zero Age HB (ZAHB – lower envelope of the observed HB star distribution). This means that to compare the models with Z00 data one has first to transform theoretical bolometric luminosities to $`V`$ magnitudes, and then determine the theoretical $`R`$ values following the definition of $`R`$ used by Z00.

S00 provides $`R`$ values for 43 GGCs, but employs a slightly different definition of $`R`$, the same as Buzzoni et al. (1983): the reference level for the RGB counts is the bolometric luminosity corresponding to the average $`V`$-magnitude of HB stars. In order to determine the level of the RGB corresponding to this bolometric luminosity S00 have applied to the observational data a relationship for the difference in bolometric corrections between HB and RGB stars (which is different from the one we applied to our models when using the Z00 definition for the $`R`$-parameter).

We have compared our theoretical results with these two databases taking into account these two different definitions of $`R`$. The derived He-abundances will thus reflect the theoretical uncertainties related to the the different bolometric corrections employed in the two methods, and the different observational samples.

In order to take into account current empirical uncertainties on the GGC metallicity scale we have used for the individual clusters \[Fe/H\] values given by Rutledge, Hesser & Stetson. (1997) on both the Carretta & Gratton (1997 – CG97) and Zinn & West (1984 – ZW84) scales (the internal accuracy of these \[Fe/H\] values is of the order of 0.10 dex). For clusters not listed by Rutledge et al. (1997) we have used the original ZW84 values transformed then to the CG97 scale using the conversion formula given by CG97.

We have determined the existence of possible He-abundance vs \[Fe/H\] correlations by computing the corresponding correlation coefficient and evaluating its significance. In case of no correlation, we have determined $`Y_{GGC}`$ by means of a weighted mean of the individual values, with weights inversely proportional to the square of the individual errors. We notice that in Z00 the value of $`Y_{GGC}`$ was determined by simply considering the constant value of the He-abundance best fitting the individual datapoints, without taking into account individual errors. The existence of a significant spread in the individual cluster He-abundance has been studied by means of the statistical F-test.

We computed new theoretical models and isochrones for Y=0.23, 0.245, 0.26, using – like in Z00 analysis – the same code, input physics and bolometric corrections as in Cassisi & Salaris (1997), with the following modifications:

* We have updated the energy loss rates from plasma-neutrino processes using the most recent and accurate results provided by Haft, Raffelt & Weiss (1994).
* We have updated the nuclear reaction rates using the NACRE database (Angulo et al. 1999), with the exception of the 12C$`(\alpha ,\gamma )^{16}`$O reaction. For this reaction we employ the more accurate recent determination by Kunz et al. (2002), based on $`\gamma `$ angular distribution measurements of 12C$`(\alpha ,\gamma )^{16}`$O and a consistent $`𝐑`$-matrix analysis of the process. The claimed relative uncertainty of this new rate is half of the uncertainty quoted in previous determinations.
* An improved EOS whose code is made publicly available at ftp://astroftp.phys.uvic.ca /pub/irwin/eos/code/eos\_demo\_fortran.tar.gz under the GNU General Public License (GPL), has been used. A full description of this EOS is in preparation (Irwin et al. 2002) so we will only summarize its principal characteristics here. The EOS is calculated using an equilibrium-constant approach to minimize the Helmholtz free-energy. For realistic abundance mixtures, this approach greatly reduces the number of linear equations that must be solved per iteration so that the solution can be rapidly obtained. This speed makes it practical to call the EOS directly from the stellar-interior code without introducing the errors associated with interpolating EOS tables (Cassisi & Irwin 2002, see also Dorman, Irwin, & Pedersen 1991). The equilibrium-constant approach gives numerical solutions of high quality with thermodynamic consistency which is typically better than 1 part in $`10^{11}`$. The “EOS1” mode of the free-energy model that is used for the present calculations includes the following: arbitrarily relativistic and degenerate free electrons (Eggleton, Faulkner, & Flannery 1973); a pressure-ionization occupation probability similar to that of Mihalas, Dappen, & Hummer (1988); a Planck-Larkin occupation probability (Rogers 1986); the exchange effect for arbitrarily relativistic and degenerate electrons (Kovetz, Lamb, & Van Horn 1972); and the Coulomb effect. The Coulomb effect is treated with the Debye-Hückel approximation in the weak coupling limit and an approximation (Pols et al. 1995) of the multicomponent combination of the one-component plasma result (DeWitt, Slattery, & Chabrier 1996) in the strong-coupling limit. A spline fit is used to interpolate between the weak and strong coupling limits. The size of the intermediate coupling region and the size of the interaction radii that characterize the pressure-ionization occupation probability are adjusted to fit the OPAL EOS tables distributed at ftp://www-phys.llnl.gov/pub/opal/eos/.
* We have explicitly taken into account the $`\alpha `$-enhanced chemical composition typical of Population II stars, using the same initial metal mixture employed by Salaris & Weiss (1998), and their same opacity tables; the heavy element distribution has an average $`\alpha `$-enhancement equal to 0.4 dex. This is potentially important for the upper metallicity end of the GGCs, since in that regime the well known equivalence between low-mass scaled-solar and $`\alpha `$-enhanced models with the same total metallicity (Salaris, Chieffi & Straniero 1993) is no longer valid (Salaris & Weiss 1998, VandenBerg et al. 2000). We just notice, in passing, that an average enhancement of 0.4 dex is in full agreement with abundance data from Halo field stars (e.g., the discussion in Salaris & Weiss 1998), while the GGCs data compiled by Carney (1996) seem to point out to an $`\alpha `$-enhancement slightly lower, of about 0.3 dex. Such a small difference – if real – does not introduce any serious bias in our final $`Y_{GGC}`$ estimates, because such a small difference in the $`\alpha `$-enhancement between clusters and models can be fully compensated (also at high metallicities) by a small rescaling of the Z-\[Fe/H\] relationship of the theoretical models, which introduces a systematic effect of less than 0.001 on the individual cluster He-abundance estimates.
* We have accounted for the different evolutionary times characterizing the red and blue parts of the HB. Ordinarily, the theoretical values of $`R`$ are computed – as in Z00 – by considering the HB evolutionary time of a star populating the middle of the RR Lyrae instability strip (log$`(T_{eff})`$=3.85). This is strictly adequate only for those clusters with an HB populated at the RR Lyrae instability strip and redward (increasing total stellar mass), since the HB evolutionary timescales are basically unchanged when moving from the instability strip towards the red (see the discussion in Z00). However, stars populating the HB blueward of the instability strip do show different evolutionary times, which increase for decreasing total stellar mass. At the bluest end of a typical blue HB the increase of the HB evolutionary time with respect to the RR Lyrae strip counterpart can amount to about 20 % (Z00). We will discuss in §3.3 the implications for the derived He-abundance in GGCs with a blue HB.

## 3 The value of $`Y_{GGC}`$ from the $`R`$-parameter

In this section we present separately our determination of $`Y_{GGC}`$ using the Z00 and S00 samples.

### 3.1 The Z00 sample

Fig. 1 displays the run of the empirical data by Z00 together with theoretical predictions for ages of 11 and 13 Gyr and $`Y`$=0.245 (solid lines), as a function of \[Fe/H\]. At fixed age and $`Y`$ the theoretical value of $`R`$ is very slowly decreasing up to \[Fe/H\]$``$1.15. Between \[Fe/H\]$``$1.15 and \[Fe/H\]$``$0.85 $`R`$ increases steeply; this increase is due to the fact that the RGB bump, previously located at brightnesses larger than the ZAHB, moves below the ZAHB level with increasing metallicity, thus causing an abrupt decrease of the number of RGB stars brighter than the ZAHB (see, e.g., the discussion in Z00). At higher \[Fe/H\] values $`R`$ is again only very mildly decreasing with increasing \[Fe/H\]. It is also interesting to notice how the dependence of $`R`$ on age is restricted to the interval ranging from \[Fe/H\]$``$1.15 to \[Fe/H\]$``$0.85, which is exactly the metallicity range where the RGB bump crosses the ZAHB level. This is easily explained by the fact that the RGB bump luminosity does depend on the stellar age (e.g., Cassisi & Salaris 1997) while the ZAHB level is basically unaffected for ages typical of GGCs; in general higher ages shift the RGB bump location towards lower luminosities.

In the same Fig. 1 we also display the theoretical $`R`$ values for an age of 13 Gyr and $`Y_{GGC}`$=0.23, to show the sensitivity of $`R`$ to the model initial He-content. The average value of the derivative $`\delta R`$/$`\delta Y`$ is $``$10.

We have estimated $`Y_{GGC}`$ and the associated 1$`\sigma `$ dispersion of the individual abundances (the latter will be thoroughly discussed at the end of this section), by first assuming an average age of 13 Gyr for the clusters (see, e.g., the analyses by VandenBerg 2000; Salaris & Weiss 1998, 2002); the error on the individual cluster $`Y`$ values has been obtained from the quoted errors on the value of $`R`$. When considering all 26 clusters together with the CG97 \[Fe/H\] scale, we obtained a weighted mean $`Y_{GGC}`$=0.240$`\pm `$0.003 (1$`\sigma `$ error). However, we found a clear correlation between $`Y_{GGC}`$ and \[Fe/H\] in the sense that the mean $`Y_{GGC}`$ obtained for clusters with \[Fe/H\] between $``$1.15 and $``$0.85 (the metallicity range influenced by the assumed cluster age) is $`Y_{GGC}`$=0.231$`\pm `$0.005, while for \[Fe/H\]$`<`$1.15 and \[Fe/H\]$`>`$0.85 we found, respectively, $`Y_{GGC}`$=0.247$`\pm `$0.005 and $`Y_{GGC}`$=0.244$`\pm `$0.003 (no correlation of the individual Y estimates with \[Fe/H\] has been found in these latter two metallicity ranges). It is evident that the mean values of Y determined for \[Fe/H\]$`<`$1.15 and \[Fe/H\]$`>`$0.85 are in good agreement while a substantial lower value is obtained for the clusters whose $`R`$ parameter is affected by the precise value of the age. We have therefore rederived the He-content with a different assumption about the cluster ages. Rosenberg et al. (1999) and Salaris & Weiss (2002) have shown how clusters with \[Fe/H\] larger than $``$1.2 (on the CG97 metallicity scale) display a large age spread and are on average younger by $``$2 Gyr than the more metal poor clusters. We have therefore recomputed the values of Y assuming an age of 13 Gyr for the clusters with \[Fe/H\]$`<`$1.2 and 11 Gyr for more metal rich ones. As expected, the mean Y values for \[Fe/H\]$`<`$1.15 and \[Fe/H\]$`>`$0.85 are unchanged, but this time, in the \[Fe/H\] range between $``$1.15 and $``$0.85, we obtain a mean $`Y_{GGC}`$=0.239$`\pm `$0.004 which, within the 1$`\sigma `$ error bar, is in better agreement with the results at higher and lower metallicities. It is therefore important to notice that the precise individual cluster ages do matter when determining an accurate $`Y_{GGC}`$ value for clusters in this \[Fe/H\] range.

We have repeated the previous analysis by employing the ZW84 \[Fe/H\] scale. Adopting an age of 13 Gyr for all clusters, we derive a mean $`Y_{GGC}`$=0.242$`\pm `$0.003 for the whole cluster sample, and we do not find any correlation between Y and \[Fe/H\]. However, Salaris & Weiss (2002) have shown that, when considering the ZW84 metallicity scale, clusters with \[Fe/H\]$`>`$1.6 show a large age spread and are on average younger than the more metal poor ones (see also VandenBerg 2000). We therefore repeated the previous calculation considering an age of 13 Gyr when \[Fe/H\]$`<`$1.6 and 11 Gyr at higher \[Fe/H\]; we obtain a mean $`Y_{GGC}`$=0.243$`\pm `$0.003, consistent with the value determined for a constant age of 13 Gyr. This result comes from the fact that, when using the ZW84 metallicities, there are no clusters populating the \[Fe/H\] range which is strongly affected by age.

Another important question to be addressed is the significance of the dispersion of the individual cluster values around the mean $`Y_{GGC}`$. In particular, it is important to know if the observed 1$`\sigma `$ dispersion, of the order of 0.02, is entirely due to the error on the individual cluster estimates. To address this point we have applied the statistical F-test (see, e.g., an application to the case of GGC ages by Chaboyer et al. 1996; Salaris & Weiss 1997, 2002) to our sample of He determinations. In case of the CG97 \[Fe/H\] scale we have restricted the analysis to the clusters within the metallicity range unaffected by age, so that an age spread will not affect the observed He-abundance dispersion. For each individual cluster we have calculated a set of synthetic He-abundances by randomly generating – using a Monte Carlo procedure – 10000 abundance values, according to a Gaussian distribution with mean value equal to the observed mean $`Y_{GGC}`$, and $`\sigma `$ equal to the individual He-abundance error. This is repeated for all clusters in the selected sample and the 10000 values for each individual clusters are joined to produce an “expected” $`Y_{GGC}`$ distribution for the entire cluster sample, on the assumption that the detected He-abundance spread is not intrinsic, but due just to the individual error bars. The F-test has been then applied in order to determine if this “expected” distribution has the same variance as the observed one. We state that a $`Y_{GGC}`$ range does exist if the probability that the two distributions have different variance is larger than 95%. In case this condition is verified, the size of the true $`Y_{GGC}`$ range ($`\sigma _Y`$) can be estimated according to $`\sigma _Y`$=$`(\sigma _{obs}^2\sigma _{exp}^2)^{0.5}`$, where $`\sigma _{obs}`$ and $`\sigma _{exp}`$ are, respectively, the 1$`\sigma `$ dispersion of the actual data and of the “expected” distribution.

The result of this test applied to the Z00 sample with our two choices of the \[Fe/H\] scale indicates that the observed dispersion around the mean $`Y_{GGC}`$ is entirely due to the formal errors (the probability that the observed and the synthetic distributions have different variance is below 70% in both cases) on the individual determinations; therefore no statistically significant spread in the individual He-abundances is found.

### 3.2 The S00 sample

Fig. 2 displays the run of the empirical data by S00 together with theoretical predictions for ages of 11 and 13 Gyr and $`Y`$=0.245 (solid lines), as a function of \[Fe/H\]. At fixed age and $`Y`$ the theoretical value of $`R`$ is very slowly decreasing up to \[Fe/H\]$``$0.85. At higher metallicities the value of $`R`$ increases, due again to the fact that the bolometric luminosity of the RGB bump crosses the reference HB bolometric luminosity. The shift to higher metallicities of this crossing region with respect to the $`R`$ definition previously used, arises from the fact that the bolometric luminosity of the RGB reference level corresponds to $`V`$ magnitudes fainter than the $`V`$ magnitude level of the HB. This implies that RGB bump stars are included into the determination of $`R`$ up to higher metallicities than the case of Z00 definition. This high metallicity region is also the only one affected by age (see discussion in §3.1); therefore the precise choice of the GGC ages does not affect at all the results when using the S00 definition of the $`R`$-parameter, since only very few clusters show these high values of \[Fe/H\] (see Fig. 2), and only on the ZW84 scale.

By assuming t=13 Gyr for all GGC we find again a $`Y_{GGC}`$ distribution uncorrelated with \[Fe/H\]. A mean value $`Y_{GGC}`$=0.246$`\pm `$0.005 is obtained when considering the CG97 \[Fe/H\] scale, while $`Y_{GGC}`$=0.241$`\pm `$0.004 is derived when the ZW84 metallicities are employed. As for the dispersion of the $`Y_{GGC}`$ values around these means, we obtain in both cases $`\sigma _Y`$=0.04; we have applied the F-test also in this case and derived that the dispersion can’t be completely explained by the formal errors on the individual determinations (the probability that the variance of the He-abundance distribution in the observed sample and in the synthetic one are different is larger than 99%), and has an intrinsic component equal to $`\sigma _Y`$=0.03 (analogous conclusions were reached by S00).

### 3.3 Clusters with a blue HB

All the $`Y_{GGC}`$ values given before have been obtained by considering the evolutionary time of HB stars in the instability strip when computing the theoretical value of $`R`$; this is also what has been done by Z00 and S00.

While this assumption is well founded in case of HBs populated at the strip and redward, it is less adequate in case of very blue HBs (see the discussion in Caputo et al. 1987 and Z00); this is particularly true when the location of the average mass populating the observed HB corresponds to $`V`$ about 0.5-1.0 mag fainter than the instability strip level. This is due to the fact, as discussed in previous section, that the HB evolutionary times increase when one greatly reduces the total stellar mass with respect to the values attained at the instability strip. To correct for this possible systematic uncertainty caused by our assumption we have applied the following procedure.

For both the Z00 and S00 samples we have identified those clusters whose HB is mainly populated at the blue side of the instability strip; among these clusters, through comparisons with our HB models, we have identified the objects whose average HB mass is located more than 0.7 mag below the RR Lyrae level. For these clusters, we have recomputed the theoretical $`R`$-values by taking as representative of the HB evolutionary lifetime the corresponding value for the average mass. There are only 6 clusters in the Z00 sample, and 8 clusters in the S00 sample that satisfy this condition.

When applying these corrected evolutionary times to the blue HB clusters we find that the $`Y_{GGC}`$ values obtained in the previous analysis are reduced by only 0.001-0.002. The size and significance of the $`Y_{GGC}`$ spread, and the behaviour with the respect to \[Fe/H\], are unchanged with respect to the previous results. In Table 1 we summarize the $`Y_{GGC}`$ results, with and without the correction for the blue HB clusters. In case of the Z00 sample and the CG97 \[Fe/H\] scale we display the results for the metallicity range that is insensitive to the choice of the cluster ages. Fig. 3 displays the distribution of the individual GGC He-abundance, for both samples and both choices of the \[Fe/H\] scale, taking into account the correction for the blue HB GGCs. The Z00 sample clearly has a significantly narrower abundance range than the S00 sample.

## 4 Discussion

In the previous section we found that the $`R`$ parameter provides values of $`Y_{GGC}`$ between $``$0.240 and $``$0.245, independent of \[Fe/H\]; the exact values are summarized in Table 1, together with the size of the intrinsic spread $`\sigma _Y`$ of the individual cluster He-abundances.

The mean values of $`Y_{GGC}`$ deduced from the Z00 and S00 sample are in excellent agreement within the associated 1$`\sigma `$ error, in spite of the - in principle - different bolometric corrections applied to the data analysis and the different photometric samples employed. It is however important to mention the fact that the Z00 data do not provide any indication of a statistically significant spread of $`Y_{GGC}`$, while the opposite is true for the S00 data. One possible reason for this occurrence may be the inhomogeneity of the S00 sample, which is made of photometries taken with very different instruments and detectors (photographic, photoelectric and CCD photometries), reduced with different procedures in the course of the last 25 years, and with possibly different methods to correct for incompleteness, as opposed to the homogeneously observed, reduced and analyzed $`HST`$ sample by Z00.

Another possibility to explain this He-abundance spread is related to the existence of population gradients within the observed clusters, coupled with the fact that the $`HST`$ data employed by Z00 mainly sample regions of the clusters’ cores, whereas the ground based photometries adopted by S00 sample more external regions located at various distances from the cores.

We have also performed another test, by comparing the individual He-abundance for 13 clusters in common between the Z00 and S00 samples. In Fig. 4 we display the abundances for these 13 clusters derived from the Z00 and S00 data (we have chosen to use in this figure the ZW84 \[Fe/H\] scale, but this choice does not affect the result of the comparison), considering the corrections for the blue HB clusters. The Z00 data provide a mean $`Y_{GGC}`$=0.237$`\pm `$0.004, in very good agreement with the result from the whole sample displayed in Table 1 ($`Y_{GGC}`$=0.240$`\pm `$0.003); the dispersion around the mean is again due (as for the whole sample) only to the error on the individual determinations. In case of the S00 data for the same 13 clusters, the mean $`Y_{GGC}`$=0.224$`\pm `$0.006 is smaller than for the Z00 data, and also significantly smaller than the mean value for the whole sample ($`Y_{GGC}`$=0.240$`\pm `$0.004); the dispersion around the mean $`Y_{GGC}`$ is larger than in the Z00 case. Therefore, whereas a somewhat random selected sizable cluster subsample (the 13 common clusters span all the relevant \[Fe/H\] range as well as show both red and blue HBs) show the same properties of the whole sample in case of the Z00 data, the opposite is true for the S00 data. This may lend some support to the idea that the significant dispersion of $`Y_{GGC}`$ for the whole S00 sample is due to some inhomogeneity intrinsic to the data used for determining the observed $`R`$ values. On the other hand, when the 4 most metal rich clusters (\[Fe/H\]$`>`$1) are excluded from the comparison shown in Fig. 4, the dispersion of the S00 data becomes comparable with the Z00 one, while the mean He-abundance is similar to the value for the whole S00 sample.

This seems to point to some metallicity-related effect, which however does not explain the dispersion for the whole S00 sample. In fact, if we apply the F-test discussed in Sections 3.1 and 3.2 to the S00 sample without the clusters with \[Fe/H\]$`>`$1, we still obtain a statistically significant dispersion of the individual He-abundances.

In spite of this difference regarding the spread in the cluster He-abundance for the Z00 and S00 samples, our results clearly indicate a mean value of $`Y_{GGC}`$ which is not in significant contradiction with the CMB constraint. This is very different from the conclusions reached by Z00 and S00 analyses, which derived an unrealistically low He-abundance, namely $`Y_{GGC}`$0.20, completely inconsistent with the CMB constraint. When we redetermine $`Y_{GGC}`$ by using the same observational data and theoretical scenario adopted by Z00, but using the same weighted average method employed in our analysis and considering the metallicity range unaffected by the selected cluster age, we obtain $`Y_{GGC}`$=0.21, still largely incompatible with the CMB constraint.

The new 12C$`(\alpha ,\gamma )^{16}`$O reaction rate and the new EOS are the two physical ingredients that have strongly modified the theoretical $`R`$ values with respect to the results by Z00, whose employed stellar models we have updated for this work. In particular, the recent estimate of the 12C$`(\alpha ,\gamma )^{16}`$O reaction rate (Kunz et al. 2002) has reduced the HB evolutionary times (at a fixed core mass and envelope composition) by $`78`$%, and the new EOS has further reduced the HB evolutionary times by $``$ 10%. On the other hand, the new EOS also slightly reduces the value of the He-core mass at the He-flash for a given age, which has the effect of increasing by $``$2% the HB evolutionary time; there is also a further increase by $``$4 % for the value of $`N_{RGB}`$ because a larger portion of the RGB is considered in the evaluation of the $`R`$-parameter. These effects cause a total reduction of $`R`$ by $``$20 % which, for a typical average observed value of $`R`$ (i.e., with the Z00 definition of $`R`$) equal to $`1.41.5`$, corresponds to an increase of the estimated $`Y_{GGC}`$ by about 0.03.

The 12C$`(\alpha ,\gamma )^{16}`$O reaction rate by Kunz et al. (2002) has a relative uncertainty of about $`\pm `$30%, which translates into a systematic uncertainty of about $`\pm `$0.008 around the values obtained in our analysis. It is also very interesting to notice at this point that Metcalfe & Handler (2002) find, from asteroseismology data for two local white dwarfs, central Oxygen abundances consistent with value obtained by using the 12C$`(\alpha ,\gamma )^{16}`$O rate by Kunz et al. (2002) during the progenitor He-burning phase.

We believe the EOS calculations for the current set of stellar models do not contribute significant errors to the final results. The new EOS has been adjusted to fit the tabulated OPAL results. The quality of the fit is quite good. For example, the residuals for solar conditions are less than 0.06% in the pressure, and this good agreement should also extend to all but the highest density portions of evolved models where there are some EOS uncertainties in the treatment of the Coulomb and electron exchange effects. However, the large variety of effects on the models caused by these non-ideal effects largely cancel each other so the calculated $`R`$ values are insensitive to these uncertainties.

The Coulomb effect arises because the attractive Coulomb force between ions and electrons tend to correlate the two kinds of particles. The exchange effect arises because the total eigenfunction of electrons, which is antisymmetric with respect to electron exchange, anti-correlates the electrons with each other. For fixed density and temperature, both the Coulomb correlation and the exchange anti-correlation reduce the amount of pressure required to confine the gas to its volume and also reduce the adiabatic gradient.

To determine how Coulomb and exchange effects alter the theoretical R values, we did some test stellar-evolution calculations with and without the Coulomb or exchange effects for main sequence, RGB and HB phases. The calculated $`R`$ value is equal to $`t_{HB}/t_{RGB}`$, where $`t_{HB}`$ is the duration of the HB evolutionary phase, and $`t_{RGB}`$ is the duration of that part of the RG phase whose luminosity is greater than the luminosity of the HB. By analyzing the various numerical experiments, we found that although the evolutionary rate along the RGB is slightly affected by Coulomb and exchange effects, the quantity $`t_{RGB}`$ is not significantly changed as a consequence of the variation of the HB luminosity level which compensates the change in the RGB evolutionary rate. On the contrary, we found the following results for $`t_{HB}`$:

* Coulomb and exchange effects for the precursor phase (i.e. the MS and the RGB) decrease the core mass by a small amount, but this reduction in fuel is more than compensated by the accompanying decrease in the helium burning luminosity so the total precursor effect for Coulomb and exchange is a 4% increase in $`t_{HB}`$. The exchange effect accounts for about one sixth of this total.
* The Coulomb effect (which for weak coupling and full ionization is proportional to the cube of the atomic number of the element) is considerably enhanced for later stages of the HB phase because He-burning in the core substantially increases the abundance of C and O.
* Coulomb and exchange effects for the HB phase increase the convective core mass by roughly 5%, but that is more than compensated by a helium burning luminosity increase of roughly twice as much. Thus, the total HB effect for Coulomb and exchange is a 6% decrease in $`t_{HB}`$. The exchange effect accounts for about one third of this total.
* When one combines the opposite precursor and HB effects together, a further cancellation occurs so the total effect for Coulomb and exchange is only a 2% decrease in $`t_{HB}`$.
* An alternative spline fit to the Coulomb effect (see the EOS description in §2) with a substantially enlarged range of intermediate coupling, changed the Coulomb results by about 10 per cent of their size. This translates to a 0.2% uncertainty in $`t_{HB}`$ and calculated $`R`$, and a negligible uncertainty in the derived $`Y_{GGC}`$ value.

Models of stellar interiors are sensitive to EOS interpolation errors (Dorman, Irwin, & Pedersen 1991) so the most reliable calculational procedure is to eliminate EOS interpolation errors by calling the EOS code directly from the stellar interior code. The present EOS is fast enough so that such direct use is practical on workstation type computers, but of course still substantially slower than calculations done with interpolated EOS tables. Thus, in the interests of reducing the required computer time for the computations we interpolated tables of EOS results that were tabulated with the present EOS for the required ranges of pressure, temperature, $`Y`$, $`X_C`$, $`X_N`$, and $`X_O`$ for a fixed non-CNO metal abundance mix. The adopted grid spacings are small enough in all coordinates so that the resulting $`t_{HB}`$ values gave excellent agreement with one test calculation using direct EOS results.

As an additional test for the adequacy of our models and therefore of our inferred $`Y_{GGC}`$, we have also considered the so-called $`R_2`$ parameter, defined as the number ratio of Asymptotic Giant Branch (AGB) to HB stars (Caputo et al. 1989). This parameter is strongly sensitive to the extension of the convective cores during the HB phase, while it is fairly insensitive to the initial metal and He-abundance of the models, and the precise value of the age. A test for our treatment of the convection in the HB stellar cores is of fundamental importance, since the extension of the convective core strongly affects the evolutionary time along the HB phase. An underestimate of the size of the HB convective cores would cause an underestimate of the HB evolutionary times, with a consequent spurious increase of $`Y_{GGC}`$. To compare theory with observations we have used the database by S00, which also provides the number of AGB stars for each cluster. These empirical data confirm the negligible effect of \[Fe/H\] on $`R_2`$; the mean value for the 43 clusters by S00 is $`R_2`$=0.14$`\pm `$0.05.

In our models we have treated the convective mixing during central He-burning by including semiconvection, following the prescriptions by Castellani et al. (1985). We have suppressed the onset of the breathing pulses during the latest phases of central He-burning by imposing that the allowed extension of the convective core does not lead to an increase of the central He abundance from one model to the next one (Caputo et al. 1989). All our models have reached the thermal pulses phase along the AGB; from this moment on the evolution is so fast that neglecting the computation of the thermal pulses does not affect the theoretical value of $`R_2`$. Our computations provide $`R_2=`$0.12, in good agreement with the empirical result; this confirms the adequacy of our mixing treatment in the HB stellar cores. If breathing pulses are not inhibited, HB evolutionary times are longer, due to the ingestion of fresh He into the central convective region following the onset of the pulses. We obtain in this case $`R_2`$0.08, in disagreement with observations (a similar conclusion was reached by Caputo et al. 1989 by comparing their models with data about the GGC M5).

We have also experimented with an alternative procedure to inhibit the onset of the breathing pulses; following the suggestions by Dorman & Rood (1993) we have set to zero the gravitational term in the energy generation equation for the central stellar regions during the later stage of core He-burning. In this way the breathing pulses are also effectively inhibited (see the detailed discussion by Dorman & Rood 1993), and we obtained a decrease of both AGB and HB evolutionary time by about 2% with respect to the procedure followed in our reference models; this leaves the value of $`R_2`$ unchanged ($`R_2`$0.12), and causes a systematic increase of $`Y_{GGC}`$ by $``$0.003.

The error on the individual $`Y_{GGC}`$ values displayed in Table 1 comes basically from the random error on the individual He-abundance determinations (i.e., from the random error on the individual $`R`$-parameter estimates). In order to give a best estimate for $`Y_{GGC}`$ including also the sources of systematic error described before (associated to uncertainties in the theoretical models) and the effect of the still uncertain \[Fe/H\] scale, we have resorted to a Monte Carlo technique briefly explained in the following. We have considered as reference values for $`Y_{GGC}`$, the ones determined by adopting the CG97 \[Fe/H\] scale, taking into account the corrections for the blue HB clusters (lines 5 and 7 of Table 1 for, respectively, the Z00 and S00 sample); we notice that in case of the Z00 data we consider the subsample unaffected by the precise choice of the clusters’ age. Starting from each of these two reference $`Y_{GGC}`$ we have generated a set of 10000 synthetic He-abundance values, by applying (through a Monte Carlo simulation) to each generated abundance value a set of random and systematic errors, according to a given probability distribution. In particular, random errors have been modeled according to a Gaussian distribution with mean value equal to the reference one, and $`\sigma `$ equal to the corresponding random error on $`Y_{GGC}`$ provided in Table 1. The systematic uncertainties due to the choice of the \[Fe/H\] scale (which causes a decrease of $`Y_{GGC}`$ by 0.003 with respect to the reference value), 12C$`(\alpha ,\gamma )^{16}`$O reaction rate (variation by $`\pm `$0.008), and breathing pulses suppression technique (increase by 0.003) have been modeled using an uniform distribution spanning the appropriate range.

The mean values for the two final synthetic distributions of He-abundances are $`Y_{GGC}`$=0.243$`\pm `$0.006 in case of the Z00 sample, and $`Y_{GGC}`$=0.244$`\pm `$0.006 for the S00 sample. These values are, as expected, in very good reciprocal agreement and moreover compare well with the primordial He-abundance $`Y_p`$=0.248$`\pm `$0.001 inferred from the CMB in conjunction with primordial nucleosynthesis computations.

Another important result of our analysis is the fact that there is no statistically significant increase of $`Y_{GGC}`$ with \[Fe/H\], at least within the analyzed GGC samples. This bears considerable interests for studies about Galactic chemical evolution. As a test for the reliability of this result we have performed the following numerical experiment. We have considered the clusters of the Z00 sample and the ZW84 metallicity scale (we obtain an analogous result when using the CG97 scale); for each cluster we have considered a reference $`R`$ value obtained from our theoretical models, using a primordial He-mass fraction of 0.245 and assuming a value for the chemical enrichment ratio $`\mathrm{\Delta }Y/\mathrm{\Delta }Z`$. We have then generated, using a Monte Carlo procedure, 10000 synthetic He-abundances for each individual cluster and a given choice of $`\mathrm{\Delta }Y/\mathrm{\Delta }Z`$, using a Gaussian distribution with mean value equal to the reference theoretical $`R`$ value and $`\sigma `$ equal to the actual random error on observed $`R`$ value.

For each of these synthetic samples we have then tried to recover the input $`\mathrm{\Delta }Y/\mathrm{\Delta }Z`$ value; we concluded from this analysis that ratios $`\mathrm{\Delta }Y/\mathrm{\Delta }Z>`$1 should have been unambiguously detected even taking into account the actual observational errors on the determination of $`R`$.

## 5 Summary

Following recent precise determinations of the primordial He-abundance coming from CMB analyses and primordial nucleosynthesis computations, we have rederived the initial He-abundance for stars in two samples of GGCs (Z00 and S00), using the $`R`$-parameter as abundance indicator. We have employed theoretical models computed adopting new and more accurate determinations of the EOS for the stellar matter and of the crucial 12C$`(\alpha ,\gamma )^{16}`$O reaction rate. Our models include semiconvection, while the breathing pulses are suppressed, in agreement with the observational constraints coming from the measurements of the $`R_2`$ parameter in the S00 sample.

By taking into account the uncertainties in the observed individual $`R`$ value, as well as the uncertainties in the GGC metallicity scale, the 12C$`(\alpha ,\gamma )^{16}`$O reaction rate and the method for the breathing pulses suppression, we obtain $`Y_{GGC}`$=0.243$`\pm `$0.006 in case of the Z00 sample, and $`Y_{GGC}`$=0.244$`\pm `$0.006 in case of the S00 sample. These abundances are in good reciprocal agreement and fully consistent with $`Y_p`$=0.248$`\pm `$0.001 recently determined from CMB analyses and primordial nucleosynthesis computations. Within the S00 sample we find a statistically significant spread of the individual He-abundances. This spread in the He-abundances is not found in the Z00 sample, and we argue that it is due to the inhomogeneity of the observational database used by S00, as opposed to the homogeneously observed and reduced photometry employed by Z00.

It is important to remark that none of the two samples show any statistically significant increase of $`Y_{GGC}`$ with the cluster \[Fe/H\], a fact that is relevant in the context of the chemical evolution of the Galaxy.

S.C. and M.S. gratefully acknowledge the hospitality of the Max-Planck-Institut für Astrophysik, where a large part of this work has been carried out. A.W.I. gratefully acknowledges partial financial support from an operating grant to Don A. VandenBerg from the Natural Sciences and Engineering Research Council of Canada. S.C. has been supported by MURST (Cofin2002). We wish thank the referee, E. Sandquist, for useful remarks which helped to improve the presentation of the paper.


text:
# Time-dependent wave-packet approach for fusion reactions of halo nuclei

## 1 INTRODUCTION

Since the discovery of the halo nuclei, much effort has been devoted to develop reaction theories appropriate for weakly bound projectiles. The eikonal approximation has been found to be useful for reactions at medium and high incident energies. However, even a basic picture has not yet been established for the reactions at low incident energies. For example, there have been controversial arguments on whether the fusion probability is enhanced or suppressed by the presence of the halo nucleon.

We have been analyzing the low energy reactions of halo nuclei in a three-body direct reaction model. A time-dependent wave packet approach has been developed for this purpose. We first analyzed a one-dimensional three-body model. We then made analysis of the fusion reaction of 11Be on a medium mass target, 40Ca, in a three-dimensional model. The calculation was restricted to the cases of total angular momentum being equal to zero (head-on collision). In these analyses, we have found that the fusion probability is hindered by the presence of the halo neutron. This conclusion is opposite to many other theoretical works which have claimed the increase of the fusion probability.

We here report our extended analyses with the wave packet method for the reaction of 11Be on 208Pb. For this reaction, the analyses have been made recently with the three-body model by other groups. The measurements are also available for systems close to this.

The Hamiltonian of the three-body model is time-independent. Still, there are several reasons to employ the time-dependent approach for the static problem. The time evolution of the wave packet gives us intuitive pictures for the reaction mechanism. It is not necessary to prepare complicated scattering boundary conditions for the three-body final states. The reaction probabilities for a certain energy region can be obtained at once from a single wave-packet solution. The trade-off for these advantages is a heavy computational cost. The method has also been developed and applied in the field of chemical reactions.

## 2 THE THREE-BODY MODEL

We describe a reaction of a single-neutron halo nucleus by the three-body model consisting of the halo neutron (n), the core nucleus (C), and the target nucleus (T). The projectile (P) is composed of the halo neutron and the core nucleus, P=n+C. We consider the reaction of 11Be on 208Pb. Namely, the core is the 10Be and the target is the 208Pb. The time-dependent Schrödinger equation is expressed as

$$i\mathrm{}\frac{}{t}\psi (𝐑,𝐫,t)=\left[\frac{\mathrm{}^2}{2\mu }_𝐑^2\frac{\mathrm{}^2}{2m}_𝐫^2+V_{nC}(r)+V_{CT}(R_{CT})+V_{nT}(r_{nT})\right]\psi (𝐑,𝐫,t),$$
(1)
where we denote the relative n-C coordinate as $`𝐫`$ and the relative P-T coordinate as $`𝐑`$. The reduced masses of n-C and P-T motions are $`m`$ and $`\mu `$, respectively. The n-C potential, $`V_{nC}(r)`$, is real. The n-T potential, $`V_{nT}(r_{nT})`$, is also chosen to be real. The C-T potential, $`V_{CT}(R_{CT})`$, is complex. The real part of $`V_{CT}`$ includes the Coulomb and nuclear potentials. The imaginary part of $`V_{CT}`$ describes the fusion process as a loss of the flux in the reaction between the core and the target nuclei. All the nuclear potentials are taken to be of Woods-Saxon shape. For the n-C potential, $`r_{nC}=1.3`$ fm and $`a_{nC}=0.7`$ fm. The strength of the potential is set to give the $`2s`$ orbital energy at $`0.6`$ MeV, close to the neutron separation energy of the 11Be nucleus. For the n-T potential, $`r_{nT}=1.27`$ fm and $`a_{nT}=0.67`$ fm. The strength will be varied. For the real part of the C-T potential, we employ $`r_{CT}^R=1.192`$ fm, $`a_{CT}^R=0.63`$ fm, and $`V_{CT}=46.764`$ MeV. The Coulomb potential is that of the uniformly charged sphere with the radius parameter, $`r_{CT}^C=1.26`$ fm. For the imaginary part of the C-T potential, we employ $`r_{CT}^I=0.8`$ fm, $`a_{CT}^I=0.4`$ fm, and $`W_{CT}=50`$ MeV.

In this three-body model, we define the fusion process of the three-body reaction as a loss of the flux caused by the absorptive C-T potential. Namely, the fusion is supposed to occur between the core and the target nuclei when they come closer beyond the Coulomb barrier. In this definition, the fusion probability includes both the complete and incomplete fusions, since the lost flux does not distinguish final states of the neutron.

In practical calculations, we make a partial wave expansion. There are two options for this; the body-fixed representation and the space-fixed one. For computations at finite total angular momenta, the former is superior to the latter. Here we only consider the reactions with zero total angular momentum, for which both representations merge into the same expression,

$$\psi ^{J=0}(𝐑,𝐫,t)=\underset{l=0}{\overset{l_{max}}{}}\frac{u_l(R,r,t)}{Rr}\frac{\sqrt{2l+1}}{4\pi }P_l(\mathrm{cos}\theta ),$$
(2)
where $`\theta `$ represents the angle between the two vectors $`𝐑`$ and $`𝐫`$. The $`l`$ is the relative angular momentum of the n-C motion, which is now identical to that of the P-T motion in the case of $`J=0`$. In numerical calculations, one has to restrict the summation within a finite range, $`ll_{max}`$.

As the initial wave packet, we employ a Gaussian packet for the P-T motion,

$$u_0(R,r,t_0)=\mathrm{exp}[\gamma (RR_0)^2iK_0R]\varphi _0(r),$$
(3)
where $`\varphi _0(r)`$ is the $`2s`$ orbital of the n-C motion and $`u_l=0`$ for $`l0`$ at $`t=t_0`$. The $`R_0`$ specifies the P-T separation at the center of the wave packet and $`K_0`$ indicates the P-T wave number on average at $`t=t_0`$. This wave packet state has the average total energy approximately given by $`E_0=\mathrm{}^2K_0^2/2\mu +V_{CT}(R_0)ϵ_0`$, where $`ϵ_0`$ is the binding energy of the orbital $`\varphi _0`$. The $`R_0`$ is set to be large so that $`V_{CT}`$ and $`V_{nT}`$ vanish except for the C-T Coulomb potential. We set $`R_0=25`$ fm and $`\gamma =0.05`$ fm-2.

To solve the time-dependent Schrödinger equation for the radial wave function $`u_l(R,r,t)`$, we discretize the radial variables. We treat radial region $`0<R<50`$ fm and $`0<r<60`$ fm with radial step $`\mathrm{\Delta }R=0.5`$ fm and $`\mathrm{\Delta }r=1`$ fm, respectively. The discrete variables representation is used for the second-order differential operator. It is important to employ very high partial waves $`l`$ to describe the transfer reactions in the Jacobi coordinate of the incident channel. We take up to $`l_{max}=70`$. The wave function $`u_l(R,r,t)`$ is thus described on the grid points of about 420,000. For the time evolution, we employ the Taylor expansion method. Previously we employed the Crank-Nicholson formula. The Taylor expansion method requires usage of the smaller time step and gives more accurate results. The time evolution is continued until most parts of the wave packet leave the interaction region after scatterings and the potentials from the target again becomes negligible except for the Coulomb.

## 3 WAVE PACKET DYNAMICS

The wave packet is a superposition of the solutions for various total energies around $`E_0`$. It provides us an intuitive picture for the reaction dynamics around this energy. In Fig. 1, we show a time evolution of the wave packet. The $`K_0`$ value is so chosen that the average total energy $`E_0`$ is about 38 MeV, which is close to the barrier top energy. The left-panels show the density distribution integrated over the angle $`\theta `$, $`\rho (R,r,t)=\frac{1}{R^2r^2}_l|u_l(R,r,t)|^2`$. Here the abscissa is $`r`$ and the ordinate is $`R`$. The right-panels show the density distribution integrated over $`R`$, $`\rho (r,\theta ,t)=\frac{1}{r^2}𝑑R|_lu_l(R,r,t)\left[\frac{2l+1}{2}\right]^{1/2}P_l(\mathrm{cos}\theta )|^2`$. Here the abscissa is $`x=r\mathrm{cos}\theta `$ and the ordinate is $`y=r\mathrm{sin}\theta `$. The center corresponds to the center-of-mass of the projectile, and the $`x`$ direction is parallel to the $`𝐑`$ vector. In the figure, the target nucleus approaches to the center from the right, then returns back. The time evolves from the top to the bottom.

Let us first look at the five panels in the left. The top-left panel shows the initial wave packet. There is a nodal structure of the 2s orbital in the $`r`$ direction. The distribution extends to large $`r`$ reflecting the halo structure. In the $`R`$ direction, the wave packet shows the Gaussian distribution around $`R_0`$. As the time evolves, the wave packet moves towards smaller $`R`$. The third (middle) panel roughly corresponds to the classical turning point. Some components of the wave packet pass through the C-T Coulomb barrier, while some are reflected by the barrier. The components going beyond the barrier are absorbed by the imaginary C-T potential and disappear. After the wave packet reaches the closest point, a flow to larger $`r`$ is seen. This corresponds to the breakup process. Performing the calculation with the n-T potential being switched off, we have confirmed that this is induced by the C-T Coulomb force. One can also see the wave packet component in the diagonal $`(Rr)`$ direction. This is a transfer component. By taking account of high partial waves ($`l_{max}=70`$), one can describe the transfer processes in the Jacobi coordinate system of the incident channel. A nodal structure near the right end in the bottom-left panel is due to a spurious reflection at the boundary in the $`r`$ coordinate.

The five panels in the right-hand side are useful to see behaviors of the halo neutron. The initial distribution is spherical since we assume $`2s`$ orbital in 11Be. Because of the acceleration of the 10Be core induced by the target Coulomb field, the halo neutron appears to be accelerated to the right (towards the target). Note that the center in these figures corresponds to the center of mass of the projectile. This acceleration induces the Coulomb breakup. As the target nucleus gets closer, one can also see the formation of molecule-like orbital in the third (middle) panel. In the bottom panel, the breakup component shows the fan-shaped, multiple waves. This is again caused by the spurious reflections at the boundary in $`r`$ coordinate. The transfer component is also seen as an oval distribution which locates its center at the average position of the target nucleus.

## 4 FUSION PROBABILITY

To obtain the fusion probability for a fixed incident energy from the wave packet solution, we apply the energy projection procedure. For a wave packet state $`\psi _a`$, we define the energy distribution function $`W_a(E)`$ as

$$W_a(E)=\psi _a|\delta (EH)|\psi _a,$$
(4)
where $`H`$ is the three-body Hamiltonian in Eq. (1). We evaluate Eq. (4) for the initial ($`a=i`$) and the final wave packet ($`a=f`$). The fusion probability is then obtained as a fraction of the energy distribution that was absorbed by the C-T absorptive potential.

$$P_{\mathrm{fusion}}(E)=\frac{W_i(E)W_f(E)}{W_i(E)}.$$
(5)
In calculating the energy distribution function $`W_a(E)`$, one may also employ the time evolution of the wave packet. Fourier transforming the delta function, one obtains

$$W_a(E)=\frac{1}{\pi \mathrm{}}\mathrm{Re}_0^{\mathrm{}}𝑑te^{iEt/\mathrm{}}\psi _a\left(\frac{t}{2}\right)|\psi _a\left(\frac{t}{2}\right),$$
(6)
where $`\psi _a(t)`$ is the solution of the time-dependent Schrödinger equation (1) with the initial condition $`\psi _a(0)=\psi _a`$.

Let us elucidate a basic picture of the reaction dynamics from the calculated fusion probabilities. Figure 3 shows the fusion probability of central collision (zero total angular momentum) as a function of the P-T relative incident energy. The fusion probabilities of three-body calculations are compared with that of the two-body 10Be-208Pb reaction without a halo neutron. The three-body calculations have been carried out with different n-T potentials: (i) $`V_{nT}=0`$ (dotted curve), (ii) $`V_{nT}=48`$ MeV (dot-dashed), and (iii) $`V_{nT}=42`$ MeV (dot-dot-dashed). The Fig. 1 corresponds to the case (iii), where the substantial transfer component is seen. There is small transfer component in the case of (ii).

As is seen from the figure, the fusion probabilities in the three-body calculation are always smaller than that of the two-body calculation. In other words, the fusion probability is suppressed by the presence of the halo neutron. Furthermore, this suppression is caused principally by the C-T Coulomb interaction. The n-T potential plays a minor role in this fusion reaction.

We could imagine two possible mechanisms that are responsible for the fusion suppression. The first one is a spectator role of the halo neutron. As is seen in Fig. 1, the halo neutron is emitted to the forwared direction. This indicates that the halo neutron proceeds straight keeping the incident velocity, while the core nucleus is reflected by the target Coulomb field. In this spectator picture, the relative energy between the core and the target nuclei is effectively smaller than the projectile-target relative energy. This mechanism, which we discussed previously, is expected to shift the energy dependence of the fusion probability by a constant amount. The calculated fusion probability, however, does not show such a simple shift in energy but shows the stronger suppression above the barrier. The second mechanism that could explain this property is the excitation of the n-C motion. The Coulomb interaction between the core and target nuclei induces the internal excitation of the projectile. The excitation energy will distribute in a certain energy region with a fluctuation. The energy conservation requires that the P-T relative energy decreases to cancel the projectile excitation energy. This results in the decrease of the effective incident energy with the fluctuation. The suppression of the fusion probability persists even high above the barrier, if there is a component of high excitation energy in the n-C motion. Of course the above two mechanisms are not independent but mutually related. More analyses are necessary to obtain clear understanding.

Recently, three-body calculations for the fusion reaction of the same system have been reported. These studies, however, report the conclusion opposite to ours, namely the strong enhancement of the fusion probability around the barrier energy. Let us next consider a reason for this discrepancy. In these three-body calculations, the continuum states are discretized in the momentum space. The partial wave expansion of the n-C motion is truncated at very low angular momentum, a few units in the $`\mathrm{}`$. We believe that this truncation is the origin of the discrepancy, causing even an opposite conclusion. Figure 3 shows the fusion probability calculated with $`l_{max}=70`$ (dashed curve) and that with $`l_{max}=2`$ (dotted curve). As is seen clearly, the calculation taking only low partial waves gives strong enhancement of the fusion probability at low incident energy. However, increasing the number of partial waves, this enhancement completely disappears and the fusion probability is actually suppressed in the entire energy region. Therefore, it is highly probable that the fusion enhancement obtained in Refs. is a false account derived from the calculations which did not converge with respect to the partial waves.

We have found that the partial wave expansion up to small $`l`$-value gives a qualitatively correct answer when one includes only Coulomb breakup mechanism switching-off the n-T nuclear potential. Therefore, it is the n-T nuclear potential that does not allow the partial wave truncation. A similar behavior is reported in Ref.  which points out that the inclusion of the high partial waves is necessary to obtain converged breakup cross section at low incident energies.

## 5 OUTLOOK

To calculate the fusion cross section, we need to make calculations for reactions at finite total angular momenta. This is under progress in the body-fixed frame descrption, where the computational cost scales linearly to the total angular momentum.

For the three-body direct reactions, coupled-channel frameworks discretizing the continuum channels in momentum space have been developed and successfully applied. The present real-space approach for the three-body problem is more straightforward but computationally expensive. As the rapid progresses of the computational resources, we expect the real-space approach will be more powerful in the future. We would like also to mention that we have recently developed an alternative real-space approach for the breakup reaction employing the absorbing boundary condition.