Production of deuterium, tritium, and $^3$He in central Pb+Pb collisions at 20A, 30A, 40A, 80A, and 158A GeV at the CERN SPS

Production of $d$, $t$, and $^3$He nuclei in central Pb+Pb interactions was studied at five collision energies ($\sqrt{s_{NN}}=$ 6.3, 7.6, 8.8, 12.3, and 17.3 GeV) with the NA49 detector at the CERN SPS. Transverse momentum spectra, rapidity distributions, and particle ratios were measured. Yields are compared to predictions of statistical models. Phase-space distributions of light nuclei are discussed and compared to those of protons in the context of a coalescence approach. The coalescence parameters $B_2$ and $B_3$, as well as coalescence radii for $d$ and $^3$He were determined as a function of transverse mass at all energies.


I. INTRODUCTION
The main goal of the heavy-ion program at the CERN SPS is the experimental investigation of the properties of nuclear matter under extreme conditions. In a head-on inelastic collision of lead nuclei, accelerated to an energy of several tens of GeV per nucleon, a hot and dense fireball of an extraordinary outward pressure gradient is formed. After explosion-like decompression, the fireball * deceased expands well beyond the volume defined by the geometric overlap region of the colliding nuclei resulting in a multiparticle system with strong collective behavior. Experimentally the overall dynamical evolution of the reaction can be probed by measuring particle composition, longitudinal and transverse momentum distributions of different particle species, as well as multi-particle correlations.
The study of light nuclei production is of importance for several reasons. First of all, the mechanism of cluster formation in the interior of the fireball of a heavy-ion collision is not well understood and requires further quantitative investigations. It is likely that a significant fraction of few-nucleon bound states registered near mid-rapidity are produced in a late stage of the reaction when the hadronic matter becomes diluted and most of the newly formed hydrogen and helium isotopes decouple from the source having no subsequent rescatterings. So, light nuclei may serve as probes of the fireball dynamics at the time of the freezeout.
In the simplest coalescence model [1][2][3] the yields of light nuclei are well explained as being determined solely by the distributions of their constituents (protons and neutrons) and an empirical coalescence parameter (B A ) related to the size A of the cluster. Such an approach works very well in proton-induced reactions and in nuclear interactions at low energies where collective flow effects are smaller. However, for expanding nuclear matter, the simple coalescence model needs to be modified. In relativistic heavy-ion collisions the production of nucleon composites depends on the reaction "homogeneity volume" [4], whose characteristics are likely to be sensitive not only to the nucleon phase space distributions at freezeout, but also to the strength of momentum-space correlations induced by collective flow [5]. To get insight into the structure of the source and the characteristics of its density and flow velocity profile, the parameters of the rapidity and transverse momentum distributions of clusters of different sizes need to be obtained over a large phase space region.
Another commonly employed approach to describe particle yields is based on statistical and thermal hadrochemical equilibrium models of particle production [6][7][8][9]. In a conventional thermal model, particle multiplicities are predicted dependent on the bulk thermal parameters of the system -the chemical freezeout temperature T , baryochemical potential µ B and volume V . Though recent versions of statistical thermal models are capable to describe hadron abundances from heavy-ion reactions in the range of collision energies from about 1 to several 10 3 GeV per nucleon [10][11][12], the question was raised [12][13][14] whether this approach is justified when applied to the production of nucleon clusters. The present paper does not discuss this issue, but follows most previous publications on light nuclei production in applying the statistical model to the hadron composition at freezeout and assuming that the yields remain unchanged during the further evolution of the fireball. Since thermal models basically consider particle yields integrated over the full phase space the lack of results on total multiplicities of light nuclei has up to now prevented a straightforward and quantitative test of the applicability of the statistical model approach to light nuclei production in the GeV collision energy range. Thanks to the large acceptance of the NA49 experiment, 3 He and d production can be measured and analyzed in a significant part of the finalstate phase space allowing an extrapolation of the yields to full (4π) phase space. This opens the possibility to examine the statistical approach by these little explored experimental probes.
The study of the production of light nuclei with different proton to neutron ratios in heavy-ion collisions can probe the behavior of the asymmetric dense nuclear matter Equation-of-State (EOS) for a range of densities, temperatures and proton fractions. The dense nuclear matter EOS is of fundamental importance for both nuclear physics and astrophysics. Several experimental observables which potentially reveal information about the density dependence of the symmetry energy associated with the n − p asymmetry ("symmetry energy" term in the nuclear EOS) have been proposed at low and intermediate collision energies: multifragmentation [15][16][17], nucleon directed and elliptic flow [18], the charged pion ratio π − /π + [19], and the isobaric yield ratios of light clusters [20]. Unfortunately, none of these probes is uniquely sensitive to the symmetry energy at all nuclear densities. For example, calculations within an isospin-dependent transport model demonstrated that the n/p ratio is most sensitive to the symmetry energy at subnormal densities (low collision energies), while with increasing collision energy (at supra-normal densities) the sensitivity of the π − /π + ratio to the density dependence of the symmetry energy was found to be stronger [21]. Thus, in order to map out the entire density dependence of the symmetry energy, a combination of several complementary observables in a broad range of collision energies is necessary. NA49 is not capable of detecting and identifying neutrons. Under certain circumstances, however, one might expect that the yield ratio of tritons to 3 He is a measure of the n/p ratio in the fireball, because the freezeout nucleon isospin asymmetry leads to a sizable difference in the production rates for light nuclei of different nucleon composition. This paper presents the results of a study comparing the production rates of A = 3 nuclei in the SPS energy range.
It should be noted that when production of a composite of mass number A is analyzed in the framework of a coalescence approach, the most common assumption is that the yield of neutrons is equal to that of protons. This is only partially true, because in relativistic heavyion collisions relative abundances of nucleon species are expected to change considerably during the dynamical evolution of the reaction due to multiple rescattering effects in dense hadronic matter. Indeed, at AGS energies a mid-rapidity n/p ratio R np = 1.19±0.08 was obtained for central Au+Pb collisions at 11.5A GeV [22], which differs from both the n/p-ratio in the incident nuclei prior to the interaction (≈ 1.5) and R np = 1 used in coalescence studies. The latter assumption introduces an extra systematic error into the results for the coalescence parameters B A which scales as R A−1 np . There is no experimental data on R np for heavy-ion collisions above AGS energies since usually experiments have not been capable of detecting neutrons far from beam rapidity. Thus, new experimental data on the energy dependence of the triton to 3 He asymmetry in the participant region of a central Pb+Pb collision can shed some light on the degree of chemical equilibration attained at SPS energies.
Up to now, light (anti)nuclei production has been studied extensively only in the energy range below √ s N N ≈ 20 GeV at the AGS [23][24][25][26] and SPS [27][28][29][30]. In Ref. [31] the NA49 experiment reported mid-rapidity spectra of deuterons from central Pb+Pb collisions at √ s N N = 8.8, 12.3, and 17.3 GeV. This paper presents results for a wider range of cluster species, which includes also tritons and 3 He nuclei, measured in central Pb+Pb interactions at center of mass energies from 6 to 17 GeV.
The paper is structured as follows. Section II describes the NA49 experiment and the studied data sets. Section III outlines the details of the analysis procedure for light nuclei. The main results of the paper are presented and discussed in Section IV. Section V concludes the paper with a summary of the results.

II. EXPERIMENTAL SETUP
The NA49 detector is a large acceptance magnetic spectrometer for the study of hadron production in heavy-ion collisions at the CERN SPS. The detector components are described briefly below and a complete description is given in Ref. [32]. The tracking system consists of four Time-Projection Chambers (TPCs). Two Vertex TPCs (VTPC) are placed inside two superconducting dipole magnets and provide momentum analysis. Downstream of the magnets Main TPCs (MTPC) are positioned on each side of the beam trajectory. These record the tracks of charged particles providing up to 90 measurements of the position and specific energy loss dE/dx of charged particles. A resolution of σ dE/dx ≈ 4% is achieved for the MTPCs allowing identification of charged particles in the relativistic rise region by correlating their dE/dx and momentum. Time-Of-Flight (TOF) detectors, composed each of 891 fast plastic scintillator tiles, are placed behind each MTPC. The TOF walls have a timing resolution of 60 ps and are essential for the identification of deuterons and tritons up to momenta of 12 GeV/c. A zero-degree calorimeter (VCAL) located further downstream is employed to trigger on collision centrality. The trajectory of incident beam ions is measured with three proportional counters (BPD1, 2, 3). A set of scintillation and Cherenkov counters positioned upstream of the target is used for beam definition and provides the start of the timing of the experiment.

A. Data sets
The data used in this analysis were collected in years 1996-2002. The experiment utilized a 208 Pb beam at energies of 20A, 30A, 40A, 80A and 158A GeV impinging on a lead target of 224 mg/cm 2 thickness corresponding to a 1% interaction probability. The interaction trigger selected the 12% most central collisions at 158A GeV, at other beam energies the data were recorded with a 7% central trigger. In order to obtain similar acceptance at all beam energies the strength of the magnetic field in the VTPCs was changed in proportion to the beam energy. An overview of the data sets used in this analysis including collision energy, centrality, and total number of events is given in Table I. Data at 40A and 158A GeV were recorded for two opposite polarities of the magnetic field with approximately equal number of events for each setting.
B. Time-of-flight reconstruction As described in detail in Ref. [30], straight line (MTPC) segments of reconstructed tracks were extrapolated to the TOF walls and matched with TOF hits. Corrections for the position of the hit inside the scintillator (tile) and the amplitude-dependent time-walk effect in the discriminator were applied tile-wise. Values of masssquared were then calculated from the reconstructed momentum p, the flight path to the TOF detector l and the measured time-of-flight t as where c denotes the speed of light.

C. Event and track selection
This section describes the cuts applied to select events and tracks for further analysis. In order to reduce the background from non-target interactions, only events for which the reconstructed primary vertex coordinate along the beam axis is within 1 cm from the nominal target position were retained. The fraction of events remaining after application of this cut varied slightly with bombarding energy and was about 99%.
To ensure optimal momentum resolution, tracks had to be reconstructed in a VTPC and a MTPC and were required to have more than 10 space points in the VTPCs. Short tracks were eliminated by requiring that the track segment in the MTPC was longer than 1.5 m in order to obtain good dE/dx measurements and minimize the effect of track splitting. Additionally, tracks were required to have a good quality trajectory fit. To guarantee precise time measurements and reject tracks depositing too little energy in the tiles because of the edge effect, a cut on the energy deposited in a scintillator was imposed discarding the lowest 10% of the pulse height distribution in a tile. Moreover, if more than one MTPC track candidate was matched to the same scintillator tile, resulting in an ambiguous time-of-flight measurement, these tracks were removed from the analysis.

D. Identification of light nuclei
The identification of light nuclei (d, t, 3 He) was based on momentum, dE/dx, and time-of-flight measurements. Deuteron and triton candidates were required to have a TOF hit matched to the MTPC track. Identification of deuterons and tritons was performed in momentum bins of 2 GeV/c width by selecting particles with measured values of dE/dx and m 2 within three standard deviations of the expected values (see Fig. 1(a)). The background contamination in both the deuteron and triton samples was estimated by analyzing the projection of the dE/dx versus m 2 histogram onto the m 2 -axis with an upper limit dE/dx cut applied: dE/dx < ( dE/dx t +3σ t ), where dE/dx t is the predicted value for tritons and σ t is the dE/dx resolution. The obtained distribution (see example in Fig. 1(b)) consisting of two signal peaks for d and t plus some background was then fitted to a sum of two Gaussians superimposed on an exponential plus firstorder polynomial function. The raw yield was calculated by counting the entries in mass windows of 2.5<m 2 <4.5 and 6.9<m 2 <7.9 GeV 2 /c 4 for d and t, respectively. The percentage of counts outside of the mass window was estimated from the Gaussian signal shape; the background contribution, not exceeding 10% within the studied momentum range, was subtracted from the data.
The identification of 3 He candidates can rely completely on the specific energy loss measurement in the MTPC gas on account of their double charge. Since matching to TOF is not required the kinematic acceptance is much larger than for deuterons and tritons. Due to overlap of the dE/dx bands for 3 He and 4 He at momenta above 10 GeV/c, the latter species contaminates the 3 He selection band. Based on the scaling behavior of light nuclei yields with increasing mass number A (see Section IVB, Eq. 5), the 4 He contamination in the selected candidates was, however, estimated to be below 3% at all collision energies and was neglected. Helium nuclei were selected by a 3σ cut around the predicted dE/dx position as indicated in the example shown in Fig. 2(a). In order to estimate the background of misidentified particles in the 3 He samples the distributions were projected onto the dE/dx axis in bins of momentum. The projected distributions were then fitted by a Gaussian for the signal plus a sum of a first-order polynomial and an exponential function for the background (see an example in Fig. 2(b)). The background contribution, varying with beam energy between 2% at 20A GeVand 20% at 158A GeV, was subtracted from the data.

E. Corrections
After selecting the d, t, and 3 He samples, all light nuclei candidates are binned in rapidity y, transverse momentum p t , and transverse mass m t −m defined as: where E and p z are the energy and longitudinal momentum component in the center-of-mass system, p t , p x , p y p/Z (GeV/c) 1 10 dE/dx (arb. units) the transverse momentum components and m the rest mass of the candidate. The raw yields of clusters need to be corrected for geometrical acceptance, detector efficiency, and for the losses due to the applied cuts and PID selection criteria. The correction for the limited geometrical acceptance was obtained from Monte Carlo (MC) simulations. In order to have similar statistics of simulated tracks in different phase space bins, MC tracks were taken from a flat distribution over momentum p, polar angle θ, and azimuthal angle φ. The particles were propagated through the detector setup with the GEANT3 package to determine the measurable track length and the potential number of space points. Tracks were then checked to satisfy all the geometrical fiducial cuts and number of space point cuts imposed on real data. Ionization energy loss and multiple scattering were taken into account by the GEANT3 tracking. Dead or inefficient TOF tiles not used in data analysis were removed from simulations. An acceptance map was generated for each cluster species and at each beam energy (i.e. magnetic field setting). Figure 3 shows the NA49 phase space coverage for d, t and 3 He in terms of transverse mass and rapidity at the lowest and the highest beam energy.
The tracking efficiency was studied by embedding simulated tracks into real data; it was found to be close to 100%. All the corrections due to quality cuts were evaluated from the data. The inefficiency due to multiple tracks hitting the same TOF tile varies from 6% (20A GeV) to 11% (158A GeV); the corresponding corrections were obtained by counting rejected tracks. The correction for the TOF hit pulse height cut is typically 10% and is weakly dependent on the beam momentum and particle species. Corrections for absorption along the trajectory through the detector material were, however, not taken into account since nucleus-nucleus interactions are not well described in GEANT3. Nevertheless, a rough estimate of the losses due to inelastic interactions of d, t, and 3 He in the material can be made based on the simulated absorption of protons. The latter was obtained by switching on and off nuclear interactions in GEANT3 and then comparing the fraction of protons that survived when passing through the detector reaching the TOF wall. Because the momenta of protons registered in the TOF wall is above 1 GeV/c, the absorption loss for protons having a TOF hit was found to be weakly dependent on rapidity and p t and did not exceed 3.5%.
Making an assumption on how the inelastic cross-section scales with atomic mass number A [33], the upper limit of the absorption loss was estimated to be of 5% and 7% for A=2 and A=3 nuclei, respectively. The estimated uncertainty associated with the absorption losses enters as a contribution to the overall systematic error.
Since the secondary nuclei knocked out by hadronic interactions in the material have momenta considerably lower than the low-p cut-off of the TOF acceptance of 1 GeV/c, their contribution in the analyzed data samples is negligible.

F. Systematic uncertainties
The main sources of the overall systematic uncertainty originate from extraction of the raw yields (including background subtraction) and from efficiency correction factors. Each particular contribution to the systematic uncertainty associated with the extraction procedure as well as with corrections due to the applied PID and quality cuts (described in the previous section) was estimated by varying one-by-one the respective selection criteria and repeating the analysis procedure. Thus, the uncertainty of background subtraction was estimated by varying the value of the dE/dx PID cut and changing the fit range and functional shapes for the background: the resulting raw yields were found to be consistent within 3% for d and 3 He and 5% for t, respectively. The error related to the pulse height cut is estimated at 1-2% (depending on the beam momentum). An estimate for the systematic uncertainty due to the TOF multi-hit cut was obtained from the spread of the difference of individual tile-wise multi-hit corrections with respect to the one averaged over all the scintillators. The spread was found to vary within 2% over the considered p t range and slightly depended on the collision energy. The various contributions, including uncertainties associated with the absorption losses, were added quadratically resulting in a total systematic uncertainty of 6% for deuterons and 9% for tritons and 3 He.
In order to investigate the reliability of these estimates, the datasets at 40A and 158A GeV were divided into two parts that were taken with opposite directions of the magnetic field in the VTPCs (named as ST D+ and ST D− setting). At opposite polarities of the field the produced charged particles are registered in the opposite halves (relative to the beam line) of the detector. The 3 He analysis was then performed in each data subset separately and the difference in the fiducial yields averaged over all rapidity bins was of order 8% and 10% at 40A and 158A GeV, respectively. As the final yields are the average for the ST D+ and ST D− data sets the results agree within the uncertainties estimated above.

A. Transverse momentum spectra and yields
The invariant p t spectra of identified 3 He nuclei at five collision energies in different rapidity intervals are shown in Fig. 4. The interval sizes range from 0.3 at 20A-40A GeV to 0.4 at 80A and 158A GeV. The experimental distributions were scaled down successively for clarity of presentation. The same scaling factors were used for the rapidity slices located symmetrically with respect to the center of mass. First a test was performed whether the yields in the rapidity bins symmetric relative to y = 0 are consistent. The test procedure includes the following steps. The distribution at a particular forward rapidity bin (source) was first fitted with an appropriate function: a sum of two exponential functions was employed. Then the fit parameters defining the shape of the spectrum were fixed and the "mirrored" spectrum (target) at corresponding backward rapidity was fitted with only the normalization parameter allowed to vary. Finally, χ 2 per point for the target spectrum with respect to the shape of the source distribution was determined within the common p t acceptance range. The procedure was then repeated switching the source and target spectrum. The result was that for all colliding energies and rapidity intervals the χ 2 /NDF values ranged from 0.4 to 1.9, with rapidity averaged values χ 2 /NDF of 1.3, 1.2, 1.6, 1.1, and 1.2 at 20A, 30A, 40A, 80A and 158A GeV, respectively. Such a low value of the χ 2 per degree of freedom indicates that the results from forward rapidities are consistent within their uncertainties with those at backward rapidities.
The NA49 acceptance for deuterons is sufficient to examine p t spectra of d in three rapidity intervals. Figure 5 presents the corresponding invariant p t spectra of deuterons at five collision energies. The intervals are specified in the legends of the figure.
The statistics for tritons is, however, too low to obtain meaningful p t spectra in several rapidity bins. Thus, for t the results for each p t bin were integrated over the TOF rapidity acceptance. Figure 6 shows the invariant yield of tritons versus p t at all bombarding energies.
In order to obtain dN/dy the measured p t distributions need to be extrapolated into unmeasured p t regions exploiting information on the spectral shape. For this, 3 He spectra were first tested with an exponential function: where dN/dy and T are two fit parameters, m t = p 2 t + m 2 is the transverse mass and m is the 3 He rest mass. Such a functional form reproduces most meson spectra from heavy-ion collisions quite well [34,35]. However, the description of the shapes of the mid-rapidity spectra of light nuclei by Eq. 2 is not satisfactory. As can be seen in Fig. 7, single-exponential fits (plotted with dotted lines) overestimate mid-rapidity p t spectra of 3 He at low and high transverse momenta. The degree of deviation is indicated by a typical χ 2 /N DF of about 7. A sum of two exponentials describes the midrapidity spectra much better (see dashed lines in Fig. 7). Thus such a parameterization was used for extrapolation of the p t spectra of light nuclei with respect to midrapidity. The observed difference between the two-and single-exponential fits, however, diminishes towards forward rapidities (see the results for blue down-pointing triangles and pink stars in Fig. 7). Since the singleexponential function of Eq. 2 produces much more stable fit results for spectra with limited p t coverage at low transverse momenta it was used for the extrapolation of spectra at very forward rapidities down to p t = 0. The extrapolation amounts to 3-7% of the dN/dy value for the 3 He spectra near mid-rapidity and increases to almost 40% for the results at the most forward rapidity.
For the case of deuterons the spectra from two adjacent rapidity bins were combined to obtain the parameters defining the spectral shape from a fit to a sum of two exponentials. These parameters were then fixed for the extrapolation of each spectrum from the combination to the unmeasured region. The extrapolation for the not-covered p t region is less than 10% at 158A GeV, the amount of extrapolation at lower energies varies from 5 to 25% at mid-rapidity and increases up to 70% for the most backward rapidity bin at 20A GeV. The systematic uncertainty of dN/dy arises from the uncertainty of spectra normalization and of the extrapolation procedure. Regarding the first contribution, the overall systematic uncertainty for the yields of clusters was estimated to 6-9% (see Section III F).
The systematic uncertainty of dN/dy for the data with a limited p t coverage is largely determined by the extrapolation procedure. The uncertainty associated with the extrapolation was estimated by using different functions: single exponential, sum of two exponentials and Boltzmann form. It was found that the difference in the results for the extrapolation using different fit functions for t is about 7%. For 3 He the method gives a typical uncertainty of 1-3% at mid-rapidity and approximately 10% for the most forward rapidity bin. For deuterons this systematic uncertainty varies from 5% to 15%.
The results on dN/dy are tabulated in Table II for  3 He and in Table III  The NA44 experiment studied the production of deuterons and tritons in central Pb+Pb interactions at the top SPS energy. Their experimental data [28] on d in the 10% and t in the 20% most central Pb+Pb collisions are plotted along with the present measurements in Fig. 9(e) and Fig. 8(e), respectively. Although centrality selections differ slightly and the yield of tritons is somewhat higher than that of 3 He at SPS energies (see Section IV D), the agreement between the two experiments can be considered reasonable.
In order to extrapolate the integral of dN/dy to full phase space two different parameterizations of the rapidity spectra were employed which provide lower and upper limits of the integral. For the lower limit the same parameterization was used as in Ref. [30]. There it was found that a sum of three Gaussians (one centered at midrapidity and two others displaced symmetrically relative to y = 0) describes the rapidity spectrum of deuterons from mid-central Pb+Pb collisions at 158A GeV quite well. This picture is based on the assumption that the observed particle emission pattern requires (at least) three sources: two located close to the (quasi)projectile/target rapidity and one at mid-rapidity. Fits using this function ('Fit A') are indicated in Figs. 8 and 9 by dot-dashed red lines. Extrapolations obtaining the upper limit ('Fit B') are based on the assumption that the longitudinal freezeout distribution spans the entire rapidity range with a broad minimum at mid-rapidity. This behavior was parameterized by a parabolic function with a sharp drop at  ±y beam (see black dashed histograms in Figs. 8 and 9). The total yield was then obtained by summing the measured values with the integral of the corresponding extrapolation function over the unmeasured region. The extrapolation accounts from 30% to 63% and from 20% to 85% of the 4π yield for 3 He and d, respectively, depending on the collision energy and type of extrapolation The resulting estimates for the total yields (multiplicities) of 3 He and deuterons are tabulated in Table IV for the two extrapolation functions discussed above. The average between these two estimates is plotted versus √ s N N in Figs. 10 and 11 for 3 He and d, respectively. The plotted overall uncertainty for the mean is a combination of the squares of the data point uncertainties and the extrapolation uncertainty caused by the lack of knowledge about the true shapes of the dN/dy distributions near the beam(target) rapidity. The latter uncertainty was estimated as half of the difference between the 'Fit A' and 'Fit B' extrapolations over the uncovered portion of the rapidity spectra. As can be seen from Figs. 10 and 11, cluster multiplicities decrease very fast as collision energy increases. These results may indicate a decrease of the average nucleon phase space density which determines the number of pn and pnp combinations for potential coalescence into d and 3 He, respectively.
In the framework of a statistical thermal model the abundance N C of a nucleon cluster of mass m, degeneracy factor g, charge q, and baryon number B is given by where V , T , µ B , µ q , and K 2 are the source volume, temperature, baryochemical potential, charge potential, and Bessel function of the second kind. Such models have been able to reproduce the multiplicities of differ- ent types of particles in elementary and heavy-ion interactions. There are several parameterizations for the thermal fireball parameters T , µ B , and V (or equivalently the fireball radius R) over a wide range of nuclear collision energies from AGS to LHC [8,[36][37][38]. The overall average of these predictions at SPS energies is given in Table V. The listed uncertainty is taken as half of the difference between the highest and lowest value for the fireball parameters provided by the various parameterizations. The µ q /T values were obtained from NA49 measurements of the π + /π − ratio [35] as Using these fireball parameters the mean multiplicities of d and 3 He were computed at all five collision energies according to Eq. 3. The results are plotted in Figs. 10 and 11 with blue circles and are listed in the last two columns of Table V. The overall uncertainties of the total yields were estimated by standard error propagation.
As can be seen, thermal model calculations are capable of reproducing the energy dependence of the cluster multiplicities not only qualitatively but also quantitatively. The deviation of the calculations from the measured abundances does not exceed 2 standard deviations (see insets in Figs. 10 and 11). It seems that there might be a systematic underprediction for the yield of d, which however cannot be claimed to be significant due to the correlated systematic uncertainty of the extrapolation to full phase space.  To inspect how the shape of the rapidity spectra for light nuclei varies with collision energy and atomic mass number, cluster yields are plotted in Fig. 12 as a function of the normalized rapidity y/y beam . All the rapidity distributions are concave. In order to quantify the changes the data were fitted with a parabola a+ b (y/y beam ) 2 (the fits are shown by dashed lines). The ratio of the fit parameters b/a (relative concavity) for 3 He and d is plotted in Fig. 13 as a function of √ s N N . One observes that the relative concavity of the rapidity distributions for light nuclei tends to increase with increasing beam momentum and cluster mass. Such a behavior of the invariant yields versus rapidity was earlier observed at AGS energies [24], where the relative concavity of the yield of clusters with atomic mass number from A = 2 to 4 progressively increases with A in the range of transverse momentum 0.1 < p t /A < 0.2 GeV/c.
Assuming that coalescence is the dominant process of cluster formation close to mid-rapidity, one expects the relative concavity of the rapidity spectra for 3 He to increase by the power of 3/2 relative to that for d (2/3 in case of the reverse order). In Fig. 13, the shaded area shows the b/a ratio for the 3 He spectra to the power of  Typically cluster production yields change drastically with the atomic mass number A and can be characterized by a parameter P , the 'penalty factor' for adding an extra nucleon to the system. Figure 14 presents the A dependence of the mid-rapidity yield dN/dy for p, d, and 3 He. The NA49 measurements for protons are taken from Refs. [39,40]. The data points for d and 3 He are  the numerical values of the parabolic fits to the rapidity spectra in Figs. 8 and 9 at mid-rapidity (y = 0). In a statistical approach the particle production rate is proportional to its spin degeneracy factor (2J+1), so it is reasonable to divide the deuteron rates by the factor 3/2 as it was done in Ref. [43]. The penalty factor P was then obtained from a fit to the atomic mass number dependence of dN/dy at mid-rapidity with an exponential function of the form: The fit results are drawn in Fig. 14 as dashed lines and the fitted values of the parameter P are shown in Fig. 15 as a function of √ s N N .
The same analysis can be performed on the total yields of nucleon clusters. For protons the rapidity spectra at 40A and 158A GeV from [40] were extrapolated to the full phase space employing a parametrization by the sum of three Gaussian distributions (as described above). At other energies, however, proton measurements over a phase space region sufficient for extrapolation to 4π are not available. Thus at these energies the penalty factors and their uncertainties had to be calculated from the integrated yields for 3 He and d only. The results are plotted in Fig. 15 as open symbols supplementing the data points obtained at AGS, SPS and LHC energies de- The excitation function for the penalty factor rises rapidly at small collision energies (the slope for the points based on mid-rapidity dN/dy is greater than for the 4π multiplicity data) and appears to level off at higher energies. Such a saturation behavior can be explained within thermal statistical models where the penalty factor P for cluster yields is determined by the Boltzmann factor exp[(m − µ B )/T ] (see e.g. Ref. [42]), with µ B , T , and m being the baryochemical potential, freezeout temperature, and nucleon mass, respectively. Employing the parameterizations for the energy dependence of T and µ B established in Refs. [8,[36][37][38], the Boltzmann factor was computed over the region of collision energies from √ s N N = 4 GeV to 3 TeV. The calculated excitation functions are drawn in Fig. 15 with lines of different types. As can be seen, thermal model predictions are in qualitative agreement with the measured penalty factors.
For a complete picture, one should bear in mind that there exist more data on the penalty factor for nucleon clusters detected in more restricted phase space regions. For example, analyzing the yields of light nuclei at p t /A < 300 MeV/c, a penalty factor of 48 ± 3 was found by the E864 experiment in the 10% most central Au+Pb interactions at beam energy of 11.5A GeV [43]. A value of about 223 ± 38 was obtained near zero p t by the NA52 experiment from minimum bias Pb+Pb collisions at 158A GeV [44]. Reference [45] reports a penalty factor of about 625 (at p t ∼ 0.8 GeV/c per nucleon) that was deduced using measurements by the STAR experiment in the 12% most central Au+Au collisions at √ s N N = 200 GeV [46]. All the above mentioned values of P were obtained for small regions of the final state phase space, thus they cannot be directly compared to those extracted from the integrated data. This is simply due to the strong radial flow of baryons in heavy-ion collisions resulting in a non-trivial pattern of space-momentum correlations at freezeout. This affects the probability of cluster formation differently in different phase space cells. In order to illustrate this, light nuclei production yields were studied more differentially in bins of p t /A. The results are presented in Fig. 16. As an example, panel a) shows invariant p t spectra of p, d and 3 He near midrapidity from central Pb+Pb collisions at 30A GeV (data for protons were taken from Ref. [39]). The distributions were fitted to a sum of two exponentials (shown by lines in Fig. 16)) with proper error estimates assigned to each data point of the p t spectrum. The invariant yields (with uncertainties) for p, d, and 3 He were calculated from the integrals of the fit functions at several values of p t per nucleon in the range p t /A from 0 to 1.0 GeV/c . Each triple is shown in Fig. 16(b) and was fitted to Eq. 5. The resulting values of the penalty factor P for each energy are plotted in Fig. 16(c) as a function of p t /A. An interesting regularity is observed for the p t dependence of the penalty factor: 1) the dependence on p t /A is similar for all energies, but the magnitude increases with increasing collision energy; 2) the penalties vary slowly at low p t /A and begin to rise faster above p t /A ∼ 0.5 GeV/c. The excitation functions of the penalty factor for zero p t and p t /A = 0.8 GeV/c are plotted in Fig. 16(d). As can be seen, the measurements at AGS and RHIC energies are consistent with the trend shown by the data from the SPS.
In the framework of both the thermal and coalescence approaches the penalty factor is related to the average phase space density of single nucleons f N (x, p) . becomes weaker, while the collective transverse motion gets stronger, thus explaining the observed trend of the penalty factor to increase with √ s N N . It was also found experimentally that the ratio of deuterons to protons at mid-rapidity is nearly constant over the whole collision centrality range in Pb+Pb interactions at the top SPS energy [31]. This finding (taking into account that the d/p ratio can be related to f N ) implies a small variation of the average baryon phase space density with collision centrality and thus offers an explanation for the good agreement between the NA49 measurement of the penalty factor near zero p t from central Pb+Pb collisions at 158A GeV and the one obtained by the NA52 Collaboration from a minimum bias data set (see Fig. 16(d)).

C. Analysis of transverse mass spectra
This section discusses systematic dependences of transverse mass spectra of clusters on the collision energy, rapidity, and particle mass. Commonly, m t distributions are examined either individually in terms of the characteristic inverse slope parameter T ef f (effective temperature) or simultaneously in the framework of a blast-wave (BW) model.
The first method was applied to extract both T ef f and the mean transverse kinetic energy m t − m. Figure 17 shows the fully corrected m t spectra of 3 He nuclei in rapidity slices at five bombarding energies. The rapidity  Red circles represent the NA49 data, the AGS measurement (blue triangle) is from [24], and the green star indicates the result from the ALICE experiment [41]. The thermal statistical model estimates are from [8,[36][37][38] (see text for detail).
binning, indicated in the figure, was the same as in the previous section and data from the bins symmetric relative to mid-rapidity y = 0 were combined in order to decrease statistical fluctuations. The average transverse energy m t − m was deduced from the measured data points combined with the integral over the unmeasured m t -range. The latter was computed by fitting the spectra with a sum of two exponential functions. The results for m t − m of 3 He are presented in Fig. 18 as a function of normalized rapidity y/y beam . As can be seen, the rapidity dependence at all energies follows a belllike shape. Thus, Gaussian fits were applied keeping the position of the maximum fixed at mid-rapidity y = 0. The two other parameters of the Gaussians ( m t − m at mid-rapidity and width σ y )) are plotted in Fig. 19 as a function of center-of-mass energy √ s N N . The drawn uncertainties are the fit errors. The overall systematic uncertainty for m t − m at mid-rapidity was estimated to amount to less than 5%. As can be seen from Fig. 19(a), where the present results are shown along with the measurements at AGS [24,26] and RHIC [47] energies, the mean transverse mass for 3 He as a function of the collision energy qualitatively follows the trend observed for hadrons in central A+A collisions [35]: rising at low energies and leveling off in the SPS energy region. The shape of the transverse mass spectra is mainly determined by the parameters characterizing the source (temperature, pressure and collective velocity profile). Consequently, the results suggest that at SPS energies the variations in these basic fireball parameters are small.
For d and t the acceptance (defined by requiring TOF information) is more restricted in rapidity and transverse mass (see Fig. 3). Therefore, in order to get coverage in m t −m sufficient for examination of the spectra shapes, the rapidity interval was enlarged to ∆y ≈ 0.6 and ∆y > 1.0 in case of d and t, respectively. Figure 20 presents the fully corrected transverse mass spectra for d and t in these rapidity intervals from central Pb+Pb collisions at five bombarding energies.
The numerical values of m t −m for d and 3 He at midrapidity are given in Table VI. For 3 He the numbers for m t − m along with their uncertainties are taken from the Gaussian fits in Fig. 18. The values for deuterons represent extrapolations to mid-rapidity from the measured rapidity range (the average rapidity y is -0.3, -0.35, -0.4, -0.7, and -0.7 at 20A, 30A, 40A, 80A, and 158A GeV, respectively) under the explicit assumption of a Gaussian rapidity dependence of m t − m with the width parameter σ y taken from the results for 3 He.
Since transverse mass spectra of clusters are not well described by an exponential function (see Fig. 17) they cannot be characterized by a single slope parameter in most cases. The estimates for T ef f were therefore obtained by fitting the spectra with exponential functions excluding the region of m t − m < 0.5 GeV. The results are listed in Table VI. For the case of tritons, however, extraction of the slope parameter of the spectra becomes problematic. At low m t the yields are measured far away from the central rapidity, while at larger m t the acceptance for tritons is near mid-rapidity. Therefore the shape of the triton spectra is strongly modified (becoming steeper) due to the rapidity dependence of cluster yields (see Fig. 12), thus making a reliable estimate of m t − m and its uncertainty impossible.
In central heavy-ion collisions the pressure gradient in the system generates strong transverse radial flow. Particles inside a collective velocity field acquire additional momentum proportional to the particle's mass. This implies that the average transverse kinetic energy E t kin depends on both the strength of radial flow and random thermal motion as  [24], the result from the STAR Collaboration was reported in [45].
where γ = 1/ 1 − β 2 and β is the average radial collective velocity and T the temperature. Figure 21 shows the NA49 results for the mid-rapidity value of m t − m for protons and light nuclei from central Pb+Pb collisions at 20A-158A GeV. The data points for protons were taken from Refs. [39,40], the values for d and 3 He were obtained in this study. Evidently, m t − m rises approximately linearly with mass at all collision energies. These results may look surprising because it seems unlikely that objects of a few MeV binding energy per nucleon are participating in multiple thermalization collisions which generate the common velocity field inside fireballs of about 120-140 MeV temperature. On other hand, it was demonstrated in Ref. [48] that in the frame-work of the coalescence approach the choice of a suitable parametrization for the spatial dependence of the single nucleon density can reproduce the observed mass dependence of the inverse slope parameter T ef f (or m t −m) of composites. For example, an interplay between a linear collective flow profile and a uniform density distribution gives an effective temperature rising linearly with mass.
In order to separate the contributions from random thermal and radial collective motion the data on m t −m at each collision energy in Fig. 21 were tested against Eq. 6 with two fit parameters: T and β . However, as noted in [49], the extrapolation of linear fits to zero mass (i.e. the temperature parameter T ) cannot be directly related to the source temperature since the apparent tem- perature in expanding fireballs is blue shifted as Thus, in order to obtain the true temperature, the first fit parameter was corrected by the blue-shift factor according to Eq. 7. The average transverse velocity β and source temperature at the kinetic freezeout extracted from these fits are given in Table VII and plotted in Fig. 22 with green circles.
The discussed source parameters T and β can also be estimated in the framework of a hydrodynamically inspired blast-wave (BW) model [49] by fitting the transverse mass spectra of particles of different masses simultaneously to the function where C i is the normalization for particle of type i and T is the freezeout temperature. The parameter ρ is defined as ρ = tanh −1 (β t ξ n ), where β t is the surface velocity and ξ = r/R with R the fireball radius. Furthermore, a boxlike spatial density distribution (f (ξ) = 1) and a linear velocity profile (n = 1) were assumed and thus β = 2 3 β t . As an example, the results of a BW-analysis of the NA49 experimental data on charged π and K mesons as well as protons and antiprotons [35,39,50] from central Pb+Pb collisions at 40A GeV are shown in Fig. 23. BW fits are drawn by solid curves and fit parameters (T, β t ) are listed in the inset. The systematic uncertainties of the fit parameters were estimated by varying the lower bound of the fitting interval for some species and by excluding different particles from the analysis. These uncertainties do not exceed 3-4% in most cases. The results of the BW fits are tabulated in Table VII and plotted in Fig. 22 by blue squares. In addition, recent data from the STAR Beam Energy Scan (BES) program [51] are shown by red stars. As the results of different analyses consistently indicate, the freezeout kinetic parameters (T kin , β ) do not vary significantly within the energy  Ratios of the yields of nuclear clusters with the same A but different nucleon content (such as the ratio t/ 3 He) can serve as an indicator of the isospin asymmetry in the source. The initial n/p ratio of 1.54 in lead nuclei can vary dramatically in the course of Pb+Pb reactions. During the hadron phase, multiple nucleon-nucleon and pion- nucleon inelastic collisions inside the interaction zone change this ratio. The value of n/p at freezeout can be deduced from comparing the yield of tritons (a composite of two neutrons and one proton) to that of 3 He clusters (two protons and one neutron) because the yield of each species is proportional to different combinations of the phase space densities of the isospin partners. For extracting information on the n/p ratio the shapes of the transverse momentum distributions for t and 3 He are studied first. Figure 24 presents the ratio of yields of t to 3 He as function of p t . For this particular study, the data for 3 He at each beam energy were averaged over the rapidity range of the measurements for tritons. Be-  cause of the TOF acceptance (see Fig. 3) the average rapidity for tritons depends on p t for small transverse momenta ("banana"-like acceptance) and the dependence is stronger at low beam energies. To avoid extra complications due to the change of yields with rapidity, the ratio was computed above p t ≈ 0.5 GeV/c and 0.3 GeV/c at 20A and 30A-40A GeV, respectively. The uncertainties shown in Fig. 24 are mainly associated with the triton statistics and within these uncertainties there is no evident trend with p t in the ratio. For each beam energy the dependence of the t/ 3 He ratio was fitted to a constant indicated by dashed lines in Fig. 24. The ratio of triton to 3 He yields averaged over the transverse momentum inter-   Fig. 25 (red dots) as a function of the center-of-mass energy. The decreasing trend with √ s N N suggests that a complete isospin equilibration may eventually be achieved at an energy above the SPS range. The data points from the E864 [22,24] and E878 [52] experiments give an impression of how close the t/ 3 He and n/p ratios are at AGS energies.
It is also expected that in heavy-ion collisions the n/p ratio and the π − /π + ratio should resemble each other since all these species are involved in the process of dynamical evolution of the overall isospin balance. Figure 25 also shows the NA49 data on the π − /π + ratio at mid-rapidity [35,50] (green stars) together with the measurement at lower energies from the E895 experiment [53]. The measurements indicate that, indeed, both ratios remain coupled over the AGS and SPS energy ranges.

E. Coalescence
In a coalescence approach [1,2] the invariant yield N A of clusters with charge Z and atomic mass number A is related to the product of the yields of protons N pr and  neutrons N n through the coefficient B A , the so-called coalescence parameter: where p = P A /A. Assuming that the ratio of neutrons to protons is unity, B A is then calculated by dividing the cluster yield at a given momentum P A by the A-th power of the proton yields at P A /A. Results of such a combined analysis of clusters from this study and the proton spectra measured in Ref. [39] are presented in  Figs. 26 and 27, which show B 2 and B 3 as a function of transverse mass at five collision energies. It should be noted, that in a coalescence analysis the data used for clusters and protons need to be measured in the same rapidity interval since there is in general a non-negligible rapidity dependence of the particle yields at a given m t . The available NA49 spectrometer acceptance, however, allows a common m t coverage only in the region of cluster m t −m > 0.25 GeV.
It is seen that for all collision energies the coalescence parameters are rising with transverse mass in accordance with the expectation that strong position-momentum correlations are present in the expanding source leading to a higher coalescence probability at larger values of m t [5]. When calculating the systematic uncertainty of the presented values of B 2,3 , an uncertainty associated with the proton yields needs to be included. This was estimated by comparing the NA49 results on proton yields obtained with two different analysis methods. The dN/dy values for protons from an analysis using dE/dx measurements reported in Ref. [40] differ from those based on the combined dE/dx+TOF analysis published in Ref. [39] by 5% and 6% at 40A GeV and 158A GeV, respectively. Based on these differences a systematic uncertainty of 6% was assigned to the proton yields and was further assumed not to vary with energy. Standard error propagation then led to an estimated uncertainty of 12% and 18% for B 2 and B 3 , respectively.
Published results on coalescence factors in heavy-ion experiments have been measured in different phase space regions since the experiments differed in their rapidity and p t coverage. In order to compare the present measurements for B 2 and B 3 with previously obtained results, the dependences on m t −m shown in Figs. 26 and 27 were extrapolated down to m t − m = 0 (p t = 0). For this purpose, a functional form of a 0 exp[a 1 (m t − m)] was fitted to the results obtained at each energy and is plotted by dashed lines. The fit parameter a 0 equals the coalescence parameter at p t = 0, while the value of the parameter a 1 depends on the difference between the slope parameters of the spectra for clusters and protons (i.e. a 1 ≈ (1/T prot − 1/T A )). The results on B A at p t = 0 are listed in Table VIII and plotted in Fig. 28.
As was pointed out in the introduction, the lack of knowledge about the production of neutrons in heavyion reactions introduces a bias in the determination of the coalescence parameters when employing only the proton yield. Using the results on the t/ 3 He ratio from Section IV D one can correct the values for B 2,3 obtained under the assumption of equal yields for nucleons of both types. The results in parentheses in Table VIII are the coalescence parameters for d and 3 He corrected by the ratio R np ≈ t/ 3 He.  Within the SPS energy range the variation of the coalescence parameter is less than 40% and 60% for B 2 and B 3 , respectively. Figure 28 compares the results for B 2 and B 3 at p t = 0 (not corrected by R np ) obtained here to experimental data from the Bevalac [23], AGS [24,25], SPS [27,28], and RHIC [47,54,55]. One concludes from this compilation that the coalescence parameters decrease only slowly with √ s N N over a broad range of collision energies.
In the framework of thermal models of cluster production [6,7] the coalescence parameter is a measure of the source size: B A ≈ (1/V ) A−1 . Thus, the observed energy dependence of B A implies that the transverse size of the emitting source does not change much in this energy domain. This behavior is consistent with that found in two-pion interferometry (HBT) measurements [56]. In Ref. [5], calculations implementing collective expansion of the reaction zone within the density matrix formalism demonstrated a close relation of the HBT radii to those obtained from the coalescence analysis. Using the prescription given in Ref. [5] the coalescence radii (R coal ) for deuterons and 3 He were calculated at all collisions energies. The results are shown in Fig. 29 along with the data from the AGS and RHIC. One observes that the values of R coal for d and 3 He agree within their uncertainties and increase gradually with the collision energy. The latter may indicate a small increase of the freezeout volume in this energy domain.

V. SUMMARY
This paper presents results on the production of d, t and 3 He nuclei in central Pb+Pb collisions at 20A-158A GeV recorded with the NA49 detector at the CERN SPS. The results for 3 He cover a wide range of rapidity and transverse momentum, while the measurements for d and t were possible only in regions closer to mid-rapidity and more restricted in transverse momentum. Cluster yields were determined and exhibit a concave shape as function of rapidity with an increase of the degree of concavity for heavier systems. The yields of d and 3 He integrated over the full phase space agree with thermal model predictions at all collision energies. The transverse mass spectra of clusters were measured and the average values <m t >−m were found to increase linearly with the mass. This behavior favors a combination of a box density profile with a linear velocity profile in the source of the clusters. The evolution of the isospin asymmetry in the fireball was studied using the triton to 3 He ratio. It was found to change gradually with collision energy following the trend observed in the ratio of π − to π + yields. The coalescence parameters B 2,3 were derived showing a weak energy dependence. This observation suggests only a small increase of the freezeout volume from AGS to RHIC energies.