text:
# Triply deuterated ammonia in NGC 1333

## 1 Introduction

Deuterium-bearing molecules have attracted attention in recent years. Physically, they appear to be good probes of the very cold phases of molecular clouds prior to star formation. Chemically, the isotopic composition of molecules is an important clue to their formation mechanism. There are two ways to make deuterated molecules. First, at temperatures $`\genfrac{}{}{0pt}{}{_<}{^{}}70`$ K, the gas-phase reaction equilibrium $`\mathrm{H}_3^++\mathrm{HD}\mathrm{H}_2\mathrm{D}^++\mathrm{H}_2`$ is shifted in the forward direction. Subsequent deuteron transfer from H2D+ to, e.g., CO and N2, leads to the large observed abundance ratios of DCO+/HCO+ and N2D+/N2H+ of $``$0.1, four orders of magnitude higher than the Galactic D/H ratio (Turner, 2001). The key species of this chemical scheme, H2D+, was recently detected in the Class 0 source NGC 1333 IRAS 4A (Stark et al., 1999). The high densities and low temperatures in this object promote the formation of H2D+ out of H$`{}_{}{}^{+}{}_{3}{}^{}`$ and HD, and also prevent its destruction because the major destroyer of H2D+, CO, is depleted by a factor of $``$100 due to freeze-out onto dust grains.

Alternatively, deuterium-bearing molecules can be formed on dust grains by surface chemistry. Accretion of H and CO onto grains, followed by reaction, is thought to produce solid H2CO and CH3OH. This mechanism favours deuteration because the atomic D/H ratio in the gas phase is much greater than the elemental ratio. The observed abundances of HDCO, CH3OD and D2CO in Orion and IRAS 16293 indicate their synthesis on dust grains (Turner, 1990; Charnley et al., 1997; Ceccarelli et al., 1998). The same mechanism may work for NH3, provided most nitrogen is in atomic form, and support for this route comes from detections of solid NH3 (Lacy et al., 1998; Gibb et al., 2000). However, these detections remain tentative (Dartois & d’Hendecourt, 2001), and observations of N2H+ suggest that in dense clouds, most nitrogen is in molecular form (Womack et al., 1992).

While NH2D has been observed in many sources (Saito et al., 2000; Shah & Wootten, 2001), ND2H has only been detected so far towards the cold, starless cores L134N and L1689N (Roueff et al., 2000; Loinard et al., 2001). The temperatures in these sources of $``$10 K are too low for significant evaporation of even the most volatile ices to occur. Gas-phase reactions can probably account for the observed abundance ratios (Rodgers & Charnley, 2001), but observations of ND3 would present a strong test, as the gas-phase route produces $``$3 times more ND3 than the grain surface route. As a step towards measuring the relative importance of gas-phase and solid-state deuteration, we have observed the $`J_K=1_00_0`$ transition of ND3 towards NGC 1333. Together with similar observations towards Barnard 1 by Lis et al. (2002), this is the first detection of a triply deuterated molecule in interstellar space.

## 2 Observations and Results

The rotational energy levels of symmetric top molecules are labeled by the total angular momentum $`J`$ and its projection on the molecular symmetry axis $`K`$. For NH3 and ND3, inversion motion splits each level further into states which are symmetric (s) and antisymmetric (a) upon reflection in the plane of the H or D atoms. The measured frequencies of the ND3 $`J_K=1_0^a0_0^s`$ transition are 309908.46 ($`F`$=1$``$1), 309909.69 ($`F`$=2$``$1) and 309911.53 ($`F`$=0$``$1) MHz; those of the $`J_K=1_0^s0_0^a`$ transition are 306735.58 ($`F`$=1$``$1), 306736.96 ($`F`$=2$``$1) and 306738.95 ($`F`$=0$``$1) MHz (Helminger & Gordy, 1969). Due to spin statistics, the $`J_K=1_0^a0_0^s`$ transition is stronger than the $`J_K=1_0^s0_0^a`$ transition by a factor of 10.

The Caltech Submillimeter Observatory (CSO) is a 10.4-m single-dish antenna located atop Mauna Kea, Hawaii. At 310 GHz, the CSO has an FWHM beam size of $`25^{\prime \prime }`$, slightly larger than the diffraction limit. Initial observations of ND3 $`J_K=1_0^a0_0^s`$ toward NGC 1333 IRAS 4A ($`\alpha _{2000}=03^h29^m10^s.3`$, $`\delta _{2000}=+31^{}13^{}32^{\prime \prime }`$) were carried out December 5-8, 2001. Additional data were taken January 24, 2002, with the local oscillator frequency shifted to verify the sideband origin of detected features. Weather conditions were average with zenith opacities of $``$0.1 at 225 GHz. An IF frequency of 1594 MHz was used to put the $`J_K=1_0^s0_0^a`$ transition in the image sideband. To subtract atmospheric and instrumental background, positions $`240^{\prime \prime }`$ offset were observed using the chopping secondary mirror. The main beam efficiency at the time of observations, measured through observations of Saturn, was 64%, using a planetary brightness temperature of $`135`$ K at 310 GHz.

Figure 1 (top) shows the combined data from all nights, aligned in signal frequency. This spectrum represents 151 minutes of on-source integration with a mean $`T_{\mathrm{sys}}`$ of 521 K and a source elevation of $`44^{}`$. A feature is detected at 309909.4 MHz, at which frequency neither the JPL (Pickett et al., 1998, spec.jpl.nasa.gov) nor the CDMS (Müller et al., 2001, www.cdms.de) catalog lists any plausible molecular lines other than ND3. The $`307`$ GHz component is not detected to $`T_{\mathrm{mb}}`$$`<14`$ mK (1$`\sigma `$).

Data were also obtained at a position $`23^{\prime \prime }`$ East and $`6^{\prime \prime }`$ South of the IRAS 4A position, which is where the DCO+ 3$``$2 emission peaks (Lis et al., in preparation). The bottom panel of Figure 1 shows that the ND3 line is also somewhat stronger than at the IRAS 4A position. The noise level for this spectrum is 23 mK (1$`\sigma `$).

Lines of ND3 exhibit hyperfine structure due to coupling of the 14N nuclear spin with the rotational angular momentum. Splitting due to D, owing to its small quadrupole moment, is only $``$200 kHz and remains unresolved in our data. Thus, we have fitted the observed spectrum with three Gaussian profiles with optical depth ratios of 5:3:1 and assuming that the components have equal widths and excitation temperatures. Such a fit has four parameters: optical depth $`\tau `$, central velocity $`V_0`$, line width $`\mathrm{\Delta }V`$ and intensity $`T_{\mathrm{mb}}`$=$`\tau (T_{\mathrm{ex}}T_{\mathrm{bg}})`$. We used the HFS method inside the CLASS package. Leaving all parameters free gives the values reported in columns 2–5 of Table 1. For the IRAS 4A position, fixing $`V_0`$ and/or $`\mathrm{\Delta }V`$ to values measured in other lines, $`7.0`$ and $`1.2`$ km s-1 (Blake et al., 1995) gives similar results. The fit results indicate that the line shape is consistent with the optically thin hyperfine intensity ratio of 5:3:1, but that the signal to noise is not high enough to constrain the optical depth. Without taking hyperfine broadening into account, $`\mathrm{\Delta }V`$ at the IRAS 4A position would be $`2.6`$ km s-1, much broader than other lines from the cold component of this source (Blake et al., 1995). The line shape thus confirms the assignment of the line to ND3.

From the observed line strengths integrated between 3.5 and 10.0 km s-1, given in column 6 of Table 1, we estimate the ND3 column density assuming an excitation temperature $`T_{\mathrm{ex}}`$=10 K. In the optically thin limit, but outside the Rayleigh-Jeans regime, and using a background temperature of $`T_{\mathrm{bg}}=2.7`$ K, the velocity-integrated optical depth follows from

$$T_{\mathrm{mb}}𝑑V=\frac{h\nu }{k}\left(\frac{1}{e^{h\nu /kT_{\mathrm{ex}}}1}\frac{1}{e^{h\nu /kT_{\mathrm{bg}}}1}\right)\tau 𝑑V$$

so that

$$\tau 𝑑V=\frac{c^3}{8\pi \nu ^3}AN_u(e^{h\nu /kT_{\mathrm{ex}}}1)$$

with

$$N_u=\frac{g_uN}{Q(T_{\mathrm{ex}})}e^{E_{\mathrm{up}}/kT_{\mathrm{ex}}}$$

yields the column density estimates in column 7 of Table 1. Here, $`Q(T_{\mathrm{ex}})`$ is the partition function $`_ig_ie^{E_i/kT_{\mathrm{ex}}}`$, equal to 38.3 for $`T_{\mathrm{ex}}`$=10 K. The Einstein A coefficient of $`2.57\times 10^4`$ s-1 follows from the dipole moment of 1.49 D (di Lonardo & Trombetti, 1981).

The assumed value of $`T_{\mathrm{ex}}`$ represents a kinetic temperature at which chemical fractionation should be efficient. However, at densities well below the critical density of this line, $`10^710^8`$ cm-3, $`T_{\mathrm{ex}}`$ will drop below $`T_{\mathrm{kin}}`$, which changes the column density estimate. As an example, the last column of Table 1 gives the values for $`T_{\mathrm{ex}}`$=5 K.

## 3 Abundance of ND3

To estimate the abundance of ND3 we have used the Monte Carlo radiative transfer program by Hogerheijde & van der Tak (2000, talisker.as.arizona.edu/$``$michiel/ratran.html). Lacking auxiliary data on the DCO+ position, we concentrate on NGC 1333 IRAS 4A, for which we take the temperature and density structure from Stark et al. (1999). Between the outer and inner radii of 3100 and 10 AU, temperatures increase from 13 to 320 K, and densities from $`2\times 10^6`$ to $`4\times 10^{11}`$ cm-3; $`N`$(H2)=$`3.1\times 10^{23}`$ cm-2 in a $`13^{\prime \prime }`$ beam, but strongly depends on beam size due to the $`R^2`$ density distribution. The radiative transfer model for ND3 includes the 30 terms up to $`100`$ cm-1 above ground, including the inversion splitting but not the hyperfine structure. Rate coefficients for de-excitation of NH3 in collisions with H2 from Danby et al. (1988) are used, scaled to the different reduced mass of the ND3-H2 system, and augmented with the terms that are Pauli-forbidden in NH3, and with transitions that would be ortho-para conversions in NH3. Initially, a constant abundance of ND3 (relative to H2) was assumed. The excitation of ND3 as a function of radius is calculated with the Monte Carlo program. The result is integrated over the line of sight and convolved with a $`25^{\prime \prime }`$ beam. The area under the synthetic line profile matches the observed value for ND3/H2$`=3.2\times 10^{12}`$.

As an alternative model, the ND3 abundance was assumed to follow that of H2D+. The major chemical formation path to ND3 starts with the reaction of NH3 with H2D+ and its derivatives DCO+ and N2D+, and proceeds through NH2D and ND2H. As a simple way to model this behaviour, we have assumed a constant H2D+/ND3 ratio. However, this ratio would vary in the case of a varying NH3 abundance, and if the alternative route starting with the reaction of N+ with HD competes, which is slightly endothermic. In the absence of sufficient constraints, we keep H2D+/ND3 constant. Our two assumed ND3 abundance profiles could be tested indirectly by observations of key deuterated molecules such as DCO+ and N2D+.

The H2D+ abundance profile in NGC 1333 IRAS 4A was calculated analytically by Stark et al. (1999), using assumed values for the cosmic-ray ionization rate ($`\zeta =5\times 10^{17}`$ s-1) and the abundances of HD ($`2.8\times 10^5`$) and D ($`2.8\times 10^6`$), and using a CO abundance of $`4\times 10^6`$ estimated from C17O data. Due to the small energy difference between H$`{}_{}{}^{+}{}_{3}{}^{}`$ and H2D+, the H2D+ abundance is strongly peaked toward large radii where temperatures are low. We have re-calculated the H2D+ abundance profile using $`\zeta =2.6\times 10^{17}`$ s-1, the mean of the values implied by observations of H$`{}_{}{}^{+}{}_{3}{}^{}`$ and H13CO+ towards seven massive young stars (van der Tak & van Dishoeck, 2000). This calculation also includes dust radiation which the one by Stark et al. (1999) did not. The new calculations are still consistent with the measured H2D+ line flux, and indicate an H2D+ abundance increasing from $`5\times 10^{19}`$ at a radius of $`10`$ AU to $`5\times 10^{10}`$ at $`R=3100`$ AU. To model the ND3 data, models were run for several values of the H2D+:ND3 ratio, and agreement between observed and calculated line flux was found for H2D+/ND3=46. The ND3 abundance at large radii is then $`1.0\times 10^{11}`$, a factor of 3 higher than that found assuming a constant ND3 abundance.

## 4 Chemistry of ND3

Table 2 summarizes the measured column densities of NH3 isotopomers toward NGC 1333 IRAS 4A. It is seen that $`N`$(NH3)/$`N`$(NH2D) $``$10 and $`N`$(NH3)/$`N`$(ND3) $``$1000; no observations of ND2H exist yet. The available data suggest a trend where with each H$``$D substitution, the column density drops by an order of magnitude. Current models of gas-phase chemistry, on the other hand, predict that ND2H/ND3 $`>`$ NH2D/ND2$`>`$ NH3/NH2(Rodgers & Charnley, 2001), unless deuteron transfer reactions are much more rapid than proton transfers. The same trend is expected in the case of grain surface chemistry. However, the measured column densities are uncertain by $``$30% due to calibration, so their ratios could be off by a factor of two and cannot be used to rule out either mechanism.

Since a straightforward comparison of column densities may be complicated by the differences in beam size of the data in Table 2, we have determined the NH3 abundance toward NGC 1333 IRAS 4A using the approach of § 3. The same temperature and density structure as in § 3 are used, and the original collisional rate coefficients of Danby et al. (1988). Based on the UMIST database (Le Teuff et al., 2000, www.rate99.co.uk), the rates of the major destruction reactions of NH3 do not depend on temperature. Formation of NH3 is mainly by dissociative recombination of NH$`{}_{}{}^{+}{}_{4}{}^{}`$, the rate of which has a $`T^{0.5}`$ dependence, which in the model for NGC 1333 IRAS 4A corresponds to a factor of 5. This factor is dwarfed by the exponential increase in H2D+ so only constant-abundance models have been considered for NH3. The observations of Shah & Wootten (2001) are reproduced for NH3/H2$`=1\times 10^8`$.

Rodgers & Charnley (2001) present a chemical scheme to form deuterated ammonia in the gas phase, which assumes that the branching ratios of dissociative recombination are statistical and that the relevant reaction rates are isotope-independent. Using this scheme, for NH3/NH2D=10 (Table 2), an abundance ratio of NH3/ND3 of $``$10000 is expected. The observed value of $`10^8/10^{11}=1000`$ is inconsistent with this prediction. This disagreement may indicate that a detailed chemical network is needed instead of a statistical treatment. In addition, NH3/NH2D was measured in an arcmin-sized region and may be $`<10`$ within the $`25^{\prime \prime }`$ CSO beam. Based on the observed NH3/ND3 ratio, NH3/NH2D could approach unity on small scales.

In the case of surface chemistry, the observed NH3/ND3 ratio implies an atomic D/H ratio of $``$0.15 in the gas phase (Rodgers & Charnley, 2001). This is significantly higher than the values of $`10^210^3`$ in the chemical models of Roberts & Millar (2000). Based on this discrepancy and on the closer agreement of the observed NH3/ND3 ratio with the gas-phase prediction for NH3/NH2D=10, we tentatively conclude that ion-molecule reactions are presently the preferred formation mechanism of NH3 in NGC 1333. This mechanism can reproduce the observed deuteration levels if deuteron transfer reactions are much faster than proton transfers. In the future, this conclusion should be tested through measurements of the NH3, NH2D, ND2H and ND3 abundances in a larger source sample. Any conclusion drawn from such data will be much stronger if the lines are measured with similar beam sizes. More detailed chemical networks are also needed. Such a project could constrain the relative importance of gas-phase and grain-surface deuteration as a function of environment.

###### Acknowledgements.

The CSO is supported by NSF grant AST 99-80846. HSPM acknowledges support from the Deutsche Forschungsgemeinschaft (DFG) via grant SFB 494.


format: text

text:
# Chemical Freezeout in Heavy Ion Collisions

## I Introduction

Currently, by colliding heavy ions at the Super Proton Synchrotron (SPS) and the Relativistic Heavy Ion Collider (RHIC), experimentalists QuarkMatter99 ; QM2001  have endeavored to create a deconfined state of quarks and gluons – the Quark Gluon Plasma (QGP). Some degree of thermalization in the heavy ion reaction is a prerequisite for QGP formation. It is an experimental fact that the ratios of the twenty or so hadron species produced in the heavy ion collision are close to the thermal ratios expected of an ideal gas of hadrons at a temperature of $`T160180\text{MeV}`$ Thermal . It is exciting that this temperature is close to the transition temperature to the QGP, $`T_c160\text{MeV}`$. However, after accounting for finite size corrections, the hadron ratios in $`p\overline{p}`$ and $`e^+e^{}`$ reactions are reproduced by the same thermal models used to describe heavy ion data. It seems that statistical mechanics provides a universal description of hadronization. Unlike $`p\overline{p}`$ and $`e^+e^{}`$ reactions, the $`5000`$ hadrons produced in the heavy ion collision rescatter after hadronization. Ideally, this ensemble of hadrons may be considered a hot hadronic gas; hydrodynamics then describes the subsequent evolution. With an understanding of the final hadronic expansion, the global properties of the collision can be quantified.

We first review the notion of chemical freezeout in the hadron gas Bebie-Chemical1 ; Ed-Chemical ; Leutwyler-Chemical1 ; Gavin-Chemical . In the hadron gas, the principal hadronic reactions, e.g. $`\pi \pi \rho \pi \pi `$, $`\pi KK^{}\pi K`$, $`\pi N\mathrm{\Delta }\pi N`$, do not change the yields of pions, kaons, and nucleons. Nevertheless, these reactions provide a mechanism for thermalization during the hadronic evolution. Thermal equilibration times at a temperature of $`T160\text{MeV}`$ are typically $`2\text{fm/c}`$ Raju-Nonequilibrium . In the strong interactions of light hadrons there are only three conserved currents: baryon number ($`J_B^\mu `$), strangeness ($`J_S^\mu `$), and isospin ($`J_I^\mu `$). Other hadronic reactions (for instance $`\pi NN^{}(1530)\mathrm{\Delta }\pi N\pi \pi `$) do change the total yield of pions, kaons and nucleons. Because of these reactions the system approaches chemical equilibrium. The yield of pions, kaons, and nucleons changes via such reactions until the Gibbs free energy reaches a minimum. However, the time scale of chemical equilibration is much longer than that of thermal equilibration; chemical equilibration times are typically $`\tau _{ch}200\text{fm/c}`$ Song-Chemical ; Pratt-Chemical . Therefore in the hadron gas, there are two disparate time scales, $`\tau _{th}`$ and $`\tau _{ch}`$.

In heavy ion collisions, the lifetime of the hadronic stage is approximately $`10\text{fm/c}`$, which is very short compared to $`\tau _{ch}`$ , but longer than $`\tau _{th}`$. Therefore, on the time scale of the collision, although the chemical composition is fixed at the time of hadronization, the system continues to evolve kinetically for some time until the particles breakup. The stages of the collision have been described with the following picture. First, there is chemical freezeout (hadronization). Then the hadrons evolve as a hadronic fluid until thermal freezeout (breakup). In this picture, the particles develop chemical potentials during the hadronic evolution since the total number of particles is fixed and the temperature decreases.

There is some evidence for this chemical/thermal freezeout picture. First, hadronic cascades indicate that following hadronization pions and nucleons rescatter for $`10`$ fm/c. Consequently, the pion $`p_T`$ decreases while the nucleon $`p_T`$ increases. Thus, microscopic calculations indicate that pions cool and increase the transverse flow of the nucleons HydroUrqmd . Second, Rapp and Shuryak have argued that a pion chemical potential can explain the anti-baryon yield at the SPS, since $`p`$ and $`\overline{p}`$ can be produced via the forward and backward reactions of $`\overline{p}p5\pi `$ Rapp-Chemical . Finally, the most compelling experimental evidence for the pion chemical potential comes from a combined analysis of spectra and HBT correlations Ferenc-UniversalF . This study indicates that the hadron density at freezeout is substantially lower than the density of a hadron gas at $`T160\text{MeV}`$. Furthermore, the pion phase space density, which can be extracted from from two pion correlations, is overpopulated. An overpopulated phase space is expected from a Bose gas with a chemical potential.

Hydrodynamics has been used extensively to model heavy ion reactions. However, these hydrodynamic calculations assume thermal $`and`$ chemical equilibrium. While thermal equilibration is at least plausible, chemical equilibration is certainly impossible. Typically a hydrodynamic simulation is run until a universal freezeout temperature $`T_f120140\text{MeV}`$, which is adjusted to match the pion and nucleon spectra. Subsequently, even if the spectral shape is correct, the yields of $`\overline{p},\mathrm{\Lambda },K`$ for example, are typically wrong by factors of two. To account for this discrepancy, a comparison to the data is then made by rescaling the particle yields (by hand) to their value at $`T160\text{MeV}`$. Energy-momentum and number conservation are violated in this inconsistent procedure.

One approach to the problem of chemical freezeout is to stop the hydrodynamic evolution at $`T_c`$ and then continue the evolution with a hadronic cascade HydroUrqmd ; Htoh . This approach incorporates chemical freezeout as chemical equilibration times are encoded into the hadronic cross sections in the cascade. Although hydro+cascade provides a comprehensive description of the heavy ion data and has been used extensively, the final spectra are sensitive to exactly where the switch to the cascade is made HydroUrqmd . Ideally, the switching temperature, $`T_{switch}`$, should be varied, with final results independent of this artificial parameter. However, the cascade conserves $`\pi ,K,\mathrm{\Lambda }\mathrm{}`$ number in the dominant reactions; therefore to achieve a kind of dual description chemical freezeout must be incorporated into the hydrodynamic evolution.

This incorporation may be achieved by including additional conservation laws for $`T<T_c`$. For each conserved species $`\pi ,K,\mathrm{\Lambda }\mathrm{}`$, we have

$`_\mu J_\pi ^\mu =_\mu J_K^\mu =_\mu J_\mathrm{\Lambda }^\mu =\mathrm{}=0.`$ (1)
The Equation of State(EOS) is modified and now depends on $`n_\pi ,n_K,n_\mathrm{\Lambda }\mathrm{}`$, in addition to $`ϵ`$ and $`n_B`$. Furthermore, the relationship between energy density and temperature is dramatically different. In this work, we consistently incorporate chemical freezeout into the hydrodynamic evolution of the hydro+cascade model, H2H Htoh . Similar recent efforts have been presented Grassi-Mu ; Hirano-Mu . H2H Htoh  and other hydrodynamic models Huovinen ; Hirano-Flow  have been compared extensively to available data both at the SPS and RHIC. The purpose of this work is to illustrate how chemical freezeout changes the results of these works.

Briefly, in the H2H model the initial stages of the collision are modeled with (2+1) ideal fluid dynamics, assuming Bjorken scaling in the longitudinal direction Bjorken-83 ; Ollitrault-Elliptic . The initial entropy is distributed in the transverse plane according to a Glauber model at a time, $`\tau _o=1`$ fm/c. At particular temperature, $`T_{switch}`$ (which is less than the critical temperature $`T_c=165\text{MeV}`$), the fluid is converted into hadrons via the Cooper-Frye formula CooperFrye . The Cooper-Frye formula is appended with a theta function rejecting backward moving hadrons. Initial conditions for the hydrodynamic evolution have been chosen to match the model charged particle multiplicity and net baryon number to the experimental values at the SPS (PbPb with $`\sqrt{s}=17`$ GeV A) and RHIC (AuAu with $`\sqrt{s}=130`$ GeV A) Htoh .

The calculation presented here are for the SPS (PbPb with $`\sqrt{s}=17`$ GeV A) at an impact parameter of b=6 fm. (A non-central impact parameter is taken in order to study elliptic flow. Final flow velocities and freezeout conditions do not change rapidly with impact parameter Htoh ; Kolb-Centrality .) The initial conditions (see Ref. Htoh  for more detail) are characterized by an entropy per baryon ratio, $`s/n_B=42`$, and an average initial energy density, $`<ϵ>|_{\tau _o}=5.4`$ $`\text{GeV/fm}^3`$. The EOS has an $`0.8`$ $`\text{GeV/fm}^3`$ latent heat and is identical (above $`T_c`$) to the equation of state LH8 which was used previously Htoh . As this work concentrates on the freezeout stages of the collision, the differences between the presented results and the analogous calculations at RHIC are small Teaney-Thesis . The slightly larger baryon density at the SPS does not significantly alter the freezeout dynamics which is already meson dominated at the SPS Htoh ; HydroUrqmd

Sect. II discusses the EOS in which all the hadronic species are conserved. The hydrodynamic equations are then solved with and without the additional conservation laws and the solutions contrasted. The principal result is that the particle yields of a hydrodynamic calculation can be consistently modified (increased), provided the freezeout temperature is also modified (decreased), to keep the energy density the same. Returning to the hydro+cascade approach, in Sect. III, we vary the switching temperature to the cascade, $`T_{switch}`$. For $`115\text{MeV}<T_{switch}<160\text{MeV}`$, the model spectra are insensitive to variations. Elliptic flow remains somewhat sensitive to the switching temperature, although the sensitivity is reduced when chemical freezeout is incorporated.

With these basic results, the freezeout conditions of the cascade RQMD are analyzed in Sect. III. At freezeout, the density of pions in the cascade is equal to the density of an ideal hadron gas with $`T_f100\text{MeV}`$ and $`\mu _\pi 80\text{MeV}`$. This freezeout density is universal and is independent of impact parameter, collision energy and the switching temperature, $`T_{switch}`$ Htoh ; HydroUrqmd . In evaluating the extent to which the cascade reproduces the hydrodynamics in this low temperature region, we find that the hydrodynamics cools much more quickly than the cascade. This rapid cooling impacts both spectra and $`v_2(p_T)`$.

## II Implementing Chemical Freezeout in the EOS

In this section, we construct an EOS which incorporates chemical freezeout. (For similar constructions see Bebie-Chemical1 ; Pratt-Chemical ; Grassi-Mu ; RappWambach .) Above $`T_c`$, the EOS is identical to LH8 which was described in Htoh . Only the hadronic portion of the EOS is modified.

For the purposes of this work, we consider expanding hadronic gas with only the lowest hadron multiplets. Specifically, we consider only the $`0^{}`$ and $`1^{}`$ meson octets, $`\eta ^{}`$ and $`\varphi `$, the $`\frac{1}{2}^+`$ baryon and anti-baryon octets and the $`\frac{3}{2}^+`$ baryon and anti-baryon decuplets. The strong resonant reactions in this ensemble excite mesons and baryons from the lower hadronic multiplet to the higher multiplet and are given (up to isospin and baryon/anti-baryon symmetry) by

$`\pi \pi \rho \pi \pi `$ (2)
$`\pi K^+K^+\pi K^+`$
$`\pi K^{}K^{}\pi K^{}`$
$`\pi N\mathrm{\Delta }\pi N`$
$`\pi Y\mathrm{\Sigma }^{}\pi Y`$
$`\pi \mathrm{\Xi }\mathrm{\Xi }^{}\pi \mathrm{\Xi },`$
where Y denotes the hyperons $`\mathrm{\Lambda },\mathrm{\Sigma }`$. These are certainly the most important reactions for the late stage of the heavy ion collision. If only these reactions and elastic collisions are included, then the system has 16 conserved currents: baryon number($`J_B^\mu `$), strangeness($`J_S^\mu `$), isospin ($`J_I^\mu `$) and 13 other conserved numbers, $`J_{H_i}^\mu `$, where $`H_i`$ runs over the hadron species,

$$\pi ,\overline{N},Y,\overline{Y},\mathrm{\Xi },\overline{\mathrm{\Xi }},K^{},\eta ,\eta ^{},\omega ,\varphi ,\mathrm{\Omega },\overline{\mathrm{\Omega }}.$$
(3)
Nucleon number and $`K^+`$ number are not included in the list of hadrons, but their conservation follows from the conserved hadron currents already specified and from baryon and strangeness conservation. Perhaps some of the stranger species ($`\mathrm{\Omega },\varphi `$) should not be included in the list of hadrons since they decouple early in the hadronic evolution. In addition, $`\omega `$ does not neatly fit into the thermal/chemical freezeout picture, as its lifetime is comparable to the collision lifetime. Nevertheless, these particles are included for theoretical consistency and we have found that they do not affect the bulk properties.

With the assumption of equilibrium the hydrodynamic equations become

$`_\mu T^{\mu \nu }=0`$ (4)
$`_\mu (n_iu^\mu )=0,`$ (5)
where $`T^{\mu \nu }=(e+p)u^\mu u^\nu pg^{\mu \nu }`$. This set of equations, together with the assumption of smooth flow, imply entropy conservation,

$`_\mu (su^\mu )=0.`$ (6)
An EOS – a relation between the pressure and the energy density and number densities associated with the 16 conserved currents – is required to complete the system.

The hadronic EOS is taken as a sum of independent ideal gases over the hadrons considered. This approximation is based upon the fact that the thermodynamics of pions interacting via $`\pi \pi \rho \pi \pi `$ is nearly equivalent to the thermodynamics of an ideal gas of $`\pi `$-s and $`\rho `$-s Raju-ResonanceGas . In all the ideal gas formulas, the 16 chemical potentials all appear in the combination,

$`\mu _BB+\mu _SS+\mu _II+{\displaystyle \underset{i}{}}\mu _{H_i}H_i,`$ (7)
where $`H_i`$ runs over the hadron list given in Eq. 3. For instance, $`\pi ^{}`$ has $`I=1/2,H_\pi =1`$, and all other quantum numbers zero. On the other hand, $`K^{}`$ has $`B=0`$, $`S=1`$, $`I=1/2`$, $`H_\pi =1`$, $`H_K^{}=1`$ and all other quantum numbers zero. The EOS is a relationship between the pressure $`p`$, and the energy density $`e`$, and the number densities $`\{n_i\}`$, and is therefore a function of 1+16=17 variables.

In the heavy ion collision, only a small fraction of the total 17 variable phase space is sampled. Indeed, if at some time $`n_i/s`$ is constant everywhere in space, then Eq. 6 and Eq. 4 imply that $`n_i/s`$ is a constant for all later times Greiner ; LL-HydroPg2 . (This follows directly since $`_\mu \mathrm{log}(s)=_\mu \mathrm{log}(n_i)=_\mu u^\mu `$ .) Thus with the assumption of entropy conservation, it is sufficient to know the EOS along the trajectory where $`n_i/s`$ is constant. We now specify $`n_i/s`$ and construct the hadronic EOS along this trajectory. This unique trajectory in the 17 variable phase space is called the adiabatic path below.

The ratio $`s/n_B=42`$ is fixed by the initial conditions and is chosen to reproduce the experimental proton to pion ratio Htoh . The net strangeness is zero, $`n_s=0`$. For simplicity, the net isospin is taken to be zero, $`n_I=0`$. As discussed in the introduction, at hadronization the system is born into chemical equilibrium. Therefore at $`T_c`$, $`n_H/s`$ is taken to be its chemically equilibrated value. Now, with $`n_i/s`$ a specified constant on the $`T=T_c`$ hyper-surface, as the system expands and cools, $`n_i/s`$ remains constant. The particles develop chemical potentials to ensure this constancy. At a computational level, the procedure for constructing the EOS along the adiabatic trajectory is the following. First, at $`T=T_c`$ adjust $`\mu _B`$,$`\mu _S`$, and $`\mu _I`$, until $`n_S=n_I=0`$ and $`s/n_B=42`$. Then, calculate $`n_H/s`$ at $`T_c`$. Then in small increments lower the temperature, adjust all the chemical potentials to leave $`n_i/s`$ constant, and tabulate all the thermodynamic quantities (i.e, pressure, entropy, energy density) along the way.

Consider the results of this procedure. First, a few of the chemical potentials are shown as a function of the switching temperature by the thick lines in Fig. 1(a).

We can derive an approximate formula for the chemical potentials as a function of $`T`$, for a collection of non-relativistic ideal gases. For non-degenerate and non-relativistic ideal gases the partial pressure $`p_i`$, and partial energy density $`e_i`$ of the i-th species are given by,

$`p_i`$ $`=`$ $`n_iT`$ (8)
$`e_i`$ $`=`$ $`n_im_i+(\text{Const})n_iT.`$ (9)
Since $`Ts=_iϵ_i+p_i\mu _in_i`$ ,

$`1`$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{n_i}{s}}\left({\displaystyle \frac{m_i+\mu _i}{T}}+\text{Const}\right).`$ (10)
$`n_i/s`$ is constant below for $`T<T_c`$. To keep all the terms in parenthesis constant below $`T_c`$, we require

$`\mu _i=m_i\left(1{\displaystyle \frac{T}{T_c}}\right)+{\displaystyle \frac{T\mu _c}{T_c}},`$ (11)
where $`\mu _c`$ is the chemical potential at the critical temperature. The approximation is shown by the dashed lines and works well for all particles except pions. Being the most important, pions are shown separately in Fig. 1(b). Thus, a pion chemical potential of nearly $`80\text{MeV}`$ may be acquired. Similar values have been found previously RappWambach ; Grassi-Mu . See Pratt-Chemical  for a discussion of how inelastic reactions can reduce these values.

Now we return to the EOS. The energy density as a function of temperature, with and without the chemical potentials, is shown in Fig. 2. The principal observation is

clear: without chemical freezeout the energy density drops very rapidly as a function of temperature, since the total number of particles drops. For later discussion, the temperature is tabulated as a function of energy density with and without chemical potentials, in Table 1.

However, to find the hydrodynamic solution the most important quantity is not the relationship between energy density and temperature, but the relation between the energy density and pressure – the EOS. Fig. 3 shows the pressure,

the squared sound speed $`c_s^2`$, and the entropy density versus the energy density along the adiabatic path. Along the adiabatic path, these quantities are all related because

$`\left({\displaystyle \frac{dp}{de}}\right)_{\left\{n_i/s\right\}}`$ $``$ $`c_s^2`$ (12)
$`\left({\displaystyle \frac{ds}{de}}\right)_{\left\{n_i/s\right\}}`$ $`=`$ $`{\displaystyle \frac{s}{p+e}}.`$ (13)
(The second relation follows by noting that $`\left(\frac{ds}{de}\right)_{n_i/s}=\left(\frac{ds}{de}\right)_{n_i}+\left(\frac{ds}{dn_i}\right)_e\frac{n_i}{s}\left(\frac{ds}{de}\right)_{n_i/s}`$, solving for $`\left(\frac{ds}{de}\right)_{n_i/s}`$ and using thermodynamic identities). Unlike energy density versus temperature, pressure is scarcely modified by the chemical potentials. Differentiating the pressure, we see that $`c_s^2`$ is reduced somewhat at low temperatures when chemical freezeout is incorporated. Integrating to find the entropy, we see that the entropy as a function of energy is nearly identical with and without chemical freezeout.

From the point of view of dynamics, this means that the hydrodynamic solutions, with and without chemical freezeout, are nearly the same when expressed as a function of energy density. The flow velocities on a freezeout surface of constant energy density are independent of whether or not chemical potentials are included. However, the temperature on that freezeout surface depends dramatically on chemical freezeout, as can be seen from Table 1. To make this point clear, the hydrodynamic solutions are illustrated in Fig. 4 with and without chemical

potentials. The solutions are similar; the small differences can be traced to small differences in the speed of sound.

Often in hydrodynamic simulations not incorporating chemical freezeout, the freezeout surface $`T130\text{MeV}`$ is taken Huovinen ; Hirano-Flow . Using Table 1, this roughly corresponds to a temperature of $`T100\text{MeV}`$. Although this temperature seems low, it is not out of keeping with the phenomenological freezeout temperatures extracted from thermal fits to radial SSH-FlowProfile  and recent elliptic flow data STAR-EllipticParticle . The extent to which this low freezeout temperature is seen in a hadronic cascade is addressed in the next section.

## III RQMD and Chemical Freezeout

Returning to hydro+cascade, we can study the extent to which chemically hydrodynamics reproduces the dynamic response of the cascade. An important observation is that with the assumption of chemical equilibrium, the hydrodynamics is independent of the dominant reactions in the cascade. With chemical freezeout, reactions such as $`\pi N\mathrm{\Delta }\pi N`$, are encoded into the conservation laws and ultimately into the chemical potentials which affect the spectra. Baryons are accelerated only by reducing the kinetic energy of the pions. In a hydrodynamic language, the flow velocity is increased only by rapidly decreasing the temperature.

One of the problems with the hydro+cascade approach HydroUrqmd  is the sensitivity to the switching temperature. To make a smooth transition from hydrodynamics into the cascade, the same conservation laws have to be implemented in both approaches. In Fig. 5, we adjust the switching

temperature from $`T_{switch}=160\text{MeV}`$ to $`T_{switch}=117\text{MeV}`$ , and compare spectra in which chemical freezeout is incorporated, $`\mu _\pi >0`$, to spectra in which chemical freezeout is not incorporated, $`\mu _\pi =0`$. Ideally, the spectra should be insensitive to $`T_{switch}`$.

First, notice that if chemical freezeout is not incorporated the yields of $`K`$ and $`\overline{p}`$ are reduced by factors of $`e^{\frac{M_K}{T}\frac{M_K}{T_c}}`$ and $`e^{\frac{M_N}{T}\frac{M_N}{T_c}}`$, or 3 and 7 for temperatures of $`T120\text{MeV}`$. Furthermore, the final flow is too strong, since both the the mass energy and the kinetic energy are converted into flow. Second, notice that once chemical freezeout is incorporated, the sensitivity to the switching temperature is small. This is because the freezeout energy density is higher; therefore it is “as if” the freezeout temperature were higher. Close inspection of the pion spectrum shows that the final spectrum with chemical freezeout is too cool relative to the Hydro+cascade model. This rapid cooling will be addressed below.

A naive explanation for the insensitivity to $`T_{switch}`$ is that the decrease in temperature is compensated by flow. Fig. 6 studies this explanation in detail.

Hadronic rescattering increases the mean $`M_T`$ of the nucleons while slightly decreasing the mean $`M_T`$ of the pions. Although these general trends are reproduced by the hydrodynamics, with the hydrodynamics the pion cooling is too rapid and the nucleon acceleration is slightly too large. These facts may be gleaned from a very close inspection of the four curves in Fig. 6. However, the differences are all small and it is difficult to tell the difference between chemically frozen hydrodynamics and free streaming from spectra alone.

Instead of spectra, we study the elliptic flow when chemical potentials are included. Fig. 7

shows the sensitivity of elliptic flow to $`T_{switch}`$. The elliptic flow remains fairly sensitive to the switching surface and increases by 25% as the switching temperature is lowered. The hydrodynamic contribution increases while the RQMD contribution decreases as $`T_{switch}`$ is lowered. The fact that the Hydro+RQMD curve is flatter than the Hydro Only curve indicates that hydro+cascade model is partially successful.

Additional information is gained with Fig. 8, which examines

the elliptic flow spectrum, $`v_2(p_T)`$. The rapid cooling of the hydrodynamic response with chemical freezeout (T=117 $`\mu _\pi >0`$, Hydro+RQMD) is seen as an increase in elliptic flow at modest $`p_T`$, relative to the normal curves (T=160 Hydro+RQMD). By comparing the $`T=160\text{MeV}`$ Hydro Only and Hydro+RQMD points, we see that in the cascade the pions cool only slightly. Nevertheless, the general trend is the same for both the hydro with chemical freezeout and the cascade with chemical freezeout.

We now turn to the distribution of pions at freezeout. The density of pions at freezeout is approximately constant as a function of impact parameter and beam energy Htoh . Now we interpret this density in thermal terms. In the cascade, the density of pions at freezeout is the same as the density of pions in a hadron gas with $`T110\text{MeV}`$ and $`\mu _\pi 70\text{MeV}`$. Notice Fig. 9, which shows the mean emission time

as the switching temperature is lowered. As with the spectra, the mean emission time after cascading is insensitive to the switching temperature. The system of hadrons introduced into the cascade expands until the density reaches a certain value and subsequently breaks up. Since the number of hadrons is conserved during the hydrodynamic evolution, the final breakup does not depend on where the hydrodynamics is stopped and where the cascade begins.

We now try to measure the RQMD pion density at freezeout. The cascade does not have a freezeout surface, as is normally assumed in an idealized hydro picture. Rather, particles are emitted per unit time and volume, i.e. $`\frac{dN}{\tau d\eta dxdyd\tau }`$. We define the freezeout density in the transverse plane as

$`n_\pi ^{RQMD}`$ $``$ $`{\displaystyle \frac{1}{\pi R_o^2}}{\displaystyle _{\sqrt{x^2+y^2}<R_o}}𝑑x𝑑y𝑑\tau {\displaystyle \frac{dN}{\tau d\eta dxdyd\tau }},`$ (14)
where $`R_o`$ is taken to be 3 fm. This density, $`n_\pi ^{RQMD}`$, with and without hadronic rescattering is shown in Fig. 10.

Of course, without rescattering, the density of pions should reflect the density of pions on the freezeout surface. With the EOS used in this work, this number density is simply the zero-th component of $`J_{H_\pi }^\mu `$ or

$`J_{H_\pi }^0=n_\pi +2n_\rho +n_\mathrm{\Delta }\mathrm{}.`$ (15)
This thermal density of pions is also shown in Fig. 10 and gives a reasonable description of the Hydro Only points. When rescattering is included, the density decreases until a constant value $`0.1\text{fm}^3`$, is reached. Comparing the thermal model curves to the Hydro+RQMD points, we see that the freezeout densities are equal for $`T100110\text{MeV}`$ and $`\mu _\pi 80\text{MeV}`$. Although this temperature is low, it is not out of keeping with the numbers extracted from hydrodynamic fits to the data STAR-EllipticParticle ; SSH-FlowProfile .

## IV Summary of Chemical Freezeout

We have constructed an EOS consistent with chemical freezeout. The relationship between pressure and energy density is approximately the same as the standard EOS. However, these two EOS exhibit dramatically different relationships between energy density and temperature. At a practical level, this means that the results from a hydrodynamic calculation without chemical freezeout may be consistently converted to a calculation with chemical freezeout by simultaneously adjusting the yields and the temperature to keep the energy density constant. Table 1 gives the required conversion factors.

We used this chemically consistent EOS to study the hydrodynamic + hadronic cascade model, H2H Htoh . Since, chemical freezeout is incorporated via the hadronic cross sections into the cascade, the same conservation laws are implemented in both the hydrodynamics and the cascade. With this congruence, we examined the model sensitivity to the switching temperature from hydro to cascade, $`T_{switch}`$. The spectra are insensitive to the switching temperature for $`120\text{MeV}<T_{switch}<160\text{MeV}`$. On the other hand, elliptic flow at the SPS remains mildly sensitive to $`T_{switch}`$ even when chemical freezeout is incorporated. The insensitivity of the model’s results to $`T_{switch}`$ partially validates the hydro+cascade approach.

An important feature of a chemically frozen EOS is that the principal reactions in the cascade, e.g. $`\pi N\mathrm{\Delta }\pi N`$, are encoded into the hydrodynamics through the chemical potentials. Thus, in the cascade, the nucleons are accelerated only by reducing the kinetic energy of the pions. With a chemically frozen EOS, this qualitative feature is reproduced by the hydrodynamics. However, as might be expected, the cooling and attendant acceleration are larger in the hydrodynamics than in the cascade. When the cascade is replaced with the hydrodynamics, the pion spectrum is slightly steeper while the nucleon spectrum is slightly flatter. Pion cooling may also be seen with $`v_2(p_T)`$. Here again, the small changes in $`v_2(p_T)`$ due to hadronic rescattering are strongly magnified by the hydrodynamics with chemical freezeout.

By incorporating chemical freezeout, many of the qualitative features of the cascade’s evolution are reproduced with the hydrodynamics (see Pasi-Photons  for further investigations of this point). Thus, a thermal interpretation of the cascade’s response is partially justified. Within RQMD, we find the density of pions at freezeout to be $`n_\pi 0.12\text{fm}^3`$, which corresponds to $`T110\text{MeV}`$ and $`\mu _\pi 70\text{MeV}`$ in thermal terms. Precisely these parameters have been extracted from thermal fits to hadronic data. Therefore, an understanding of the properties of a chemically frozen hadronic gas helps to bridge the cascade and thermal descriptions of the final stages of the heavy ion collision.

## V Acknowledgments

I thank Edward Shuryak for continued support.