text:
# High Energy Neutrino Astronomy: The Experimental Road

## 1 Introduction

Whereas MeV neutrino astronomy has been established by the observation of solar neutrinos and neutrinos from supernova SN1987, neutrinos with energies of GeV to PeV which must accompany the production of high energy cosmic rays still await discovery. Detectors underground have turned out to be too small to detect the feeble fluxes of energetic neutrinos from cosmic accelerators. The high energy frontier of TeV and PeV energy is currently being tackled by much larger, expandable arrays constructed in open water or ice. Detectors tailored to record acoustic, radio, fluorescence or air shower signatures from neutrino interactions at EeV energy (= $`10^9`$ GeV) and above are being designed in parallel . Fig. 1 sketches the energy domains of different techniques.

## 2 Physics Goals

The central goal of high energy neutrino telescopes is to settle the origin of high energy cosmic rays . The directional information of these charged particles - protons, light and heavy nuclei - is lost due to deflection in cosmic magnetic fields (apart from the extreme energies above $`10^{10}`$ GeV where deflection is negligible). Source tracing, i.e. astronomy, is only possible by neutral, stable particles like $`\gamma `$ rays and neutrinos. In contrast to $`\gamma `$ rays which may come from pure electron acceleration, only neutrinos provide incontrovertible evidence of proton acceleration. On top of that, neutrinos do not suffer from absorption by the omnipresent infrared or radio background when propagating through space. The range of TeV $`\gamma `$ rays is only about 100 Mpc, at PeV even only 10 kpc, i.e. about the radius of our Galaxy. Therefore, the topology of the far distant high energy Universe may possibly be investigated only with neutrinos.

The physics goals of high energy neutrino telescopes include:

1. Search for neutrinos from cosmic acceleration processes in galactic sources like micro quasars or supernova remnants (SNR), or extragalactic sources like active galactic nuclei (AGN) or gamma ray bursts (GRB),
2. search for ultra-high energy (UHE) neutrinos from interactions of UHE cosmic rays with the photons of cosmic 3K microwave background (the so called GZK neutrinos ), from topological defects (TD) or from the decay of super-heavy particles,
3. search for neutrinos from the annihilation of Weakly Interacting Massive Particles (WIMPs),
4. search for magnetic monopoles,
5. monitoring our Galaxy for MeV neutrinos from supernova bursts.

Most models related to sources of type a) assume acceleration by shock waves propagating in accretion discs around black holes or along the extended jets emitted perpendicularly to the disk (bottom-up models). Neutrinos are generated in decays of mesons produced by interactions of the accelerated charged particles with ambient matter or with photon gas.

$`p+p(\gamma )p(n)+\pi `$
$`\mu +\nu `$
The neutrino energy spectrum of many models follows an $`E_\nu ^2`$ behaviour, at least over a certain range of energy. Assuming an $`E_\nu ^2`$ form and normalizing the neutrino flux to the measured flux of cosmic rays at highest energies leads to an upper bound of $`dN/dE_\nu 5\times 10^8E_\nu ^2`$ GeV-1 cm-2 s-1 sr-1 to the diffuse neutrino flux (i.e. the flux integrated over all possible sources) . Reasonably weakened assumptions loosen this bound by more than one order of magnitude to $`10^6E_\nu ^2`$ GeV-1 cm-2 s-1 sr-1 (see also Fig.14). The so-called top-bottom scenarios of type b) are suggestive for the explanation of highest energy cosmic rays. In this scenario, high energy particles would be “born” with high energies, and not accelerated from low to high energies, as in the standard bottom-up scenarios.

We will focus to a) and b) in the following and refer to  and references therein for more information on c) \- e).

## 3 Cherenkov telescopes under water and ice

Optical underwater/ice neutrino detectors consist of a lattice of photomultipliers (PMs) housed in transparent pressure spheres which are spread over a large open volume in the ocean, in lakes or in ice. In most designs the spheres are attached to strings which - in the case of water detectors - are moored at the ground and held vertically by buoys. The typical spacing along a string is 10-20 meters, and between strings 30-100 meters. The spacing is incomparably large compared to Super-Kamiokande. This allows to cover large volumes but makes the detector practically blind with respect to phenomena below 10 GeV.

The PMs record arrival time and amplitude of Cherenkov light emitted by muons or particle cascades. The accuracy in time is a few nanoseconds. Fig. 2 sketches the two basic detection modes.

In the muon mode, high energy neutrinos are inferred from the Cherenkov cone accompanying muons which enter the detector from below. Such upward moving muons can have been produced only in interactions of muon neutrinos having crossed the earth. The effective volume considerably exceeds the actual detector volume due to the large range of muons (about 1 km at 300 GeV and 24 km at 1 PeV). Muons which have been generated in the earth atmosphere above the detector and punch through the water or ice down to the detector, outnumber neutrino-induced upward moving muons by several orders of magnitude and have to be removed by careful up/down assignment. At energies above a few hundred TeV, where the earth is going to become opaque even to neutrinos, muons arrive only from directions close to the horizon, at EeV energies even only from the upper hemisphere. Most of these muons can be distinguished from down going atmospheric muons due to their higher energy deposition.

Apart from elongated tracks, cascades can be detected. Their length increases only like the logarithm of the cascade energy. With typically 5-10 meters length, and a diameter of the order of 10 cm, cascades may be considered as quasi point-like compared to the spacing of photomultipliers in Cherenkov telescopes. The effective volume for cascade detection is close to the geometrical volume. While for present telescopes it therefore is much smaller than that for muon detection, for kilometer-scale detectors and not too large energies it can reach the same order of magnitude like the latter.

Underwater/ice telescopes are optimized for the detection of muon tracks and for energies of a TeV or above, by the following reasons:

1. The flux of neutrinos from cosmic accelerators is expected to be harder than that of atmospheric neutrinos above 1 TeV, yielding a better signal-to-background ratio at higher energies.
2. Neutrino cross section and muon range increase with energy. The larger the muon range, the larger is the effective detection volume.
3. The mean angle between muon and neutrino decreases with energy like $`E^{0.5}`$, with a pointing accuracy of about one degree at 1 TeV.
4. Mainly due to pair production and bremsstrahlung, the energy loss of muons increases with energy. Above 1 TeV, this allows to estimate the muon energy from the larger light emission along the track.

The development in this field was stimulated by the DUMAND project close to Hawaii which was cancelled in 1995. The breakthrough came from the other pioneering experiment located at a depth of 1100 m in the Siberian Lake Baikal. The Baikal collaboration not only was the first to deploy three strings (as necessary for full spatial reconstruction ), but also reported the first atmospheric neutrinos detected underwater (, see Fig. 4, left). At present, NT-200 is taking data, an array comprising 192 mushroom shaped 15”-PMs at eight strings. A moderate upgrade (NT200+) is planned for 2003/04 (see fig.3, which shows at the top the small, compact NT-200 array plus three sparsely instrumented distant strings forming together NT200+). NT200+ will allow a significantly improved cascade reconstruction within the volume framed by the new strings.

With respect to its size, NT-200 has been surpassed by the AMANDA detector . Rather than water, AMANDA uses the 3 km thick ice layer at the geographical South Pole. Holes are drilled with hot water, and strings with PMs are frozen into the ice. With 677 PMs at 19 strings, most at depths between 1500-2000 m, the present AMANDA-II array reaches an area of a few 104 m2 for 1 TeV muons. Although still far below the square kilometer size suggested by most theoretical models, AMANDA-II may be the first detector with a realistic discovery potential with respect to extraterrestrial high energy neutrinos. Limits obtained from the analysis of data taken with the smaller ten-string detector AMANDA-B10 in 1997 are similar to or below those limits which have been obtained by underground detectors over more than a decade of data taking. The limit on the diffuse flux from unresolved sources with an assumed $`E^2`$ spectrum is $`0.810^6E_\nu ^2`$ GeV-1 cm-2 s-1 sr-1 , below loosest theoretical bounds , slightly below the corresponding Baikal limit and nearly an order of magnitude below limits from underground experiments. AMANDA limits on point sources on the Northern sky  complement the limits obtained from detectors on the Northern hemisphere for the Southern sky (see Fig. 15). The sensitivity of AMANDA-B10 has been verified by samples of events which are dominated by atmospheric neutrinos . Fig. 4 (right) shows a neutrino event taken with AMANDA-B10, Fig. 5 the sky map of first neutrino candidates taken in 1997.

Based on the experience from AMANDA, a cubic kilometer detector, ICECUBE , is going to be deployed at the South Pole. It will consist of 4800 PMs at 80 vertical strings, with 125 m inter-string-distances and a 16 m spacing between the PMs along a string (see fig.6). The 8-inch AMANDA PMs will be replaced by 10-inch PMs. As for the recently upgraded Amanda read-out, full transient waveforms will be recorded from every PM.

Two projects for large neutrino telescopes are under construction in the Mediterranean - ANTARES  and NESTOR  \- see fig.3. Both have assessed the relevant physical and optical parameters of their sites, developed deployment methods, performed a series of operations with a few PMs and layed underwater cables to the future locations of the detectors. ANTARES and NESTOR envision different deployment schemes and array designs. The NESTOR group plans to deploy a tower of several floors, each carrying 12 PMs at 16 m long arms. Later, a ring consisting of 72 former DUMAND PMs is planned. The ANTARES detector will consist of 12 strings, each equipped with 30 triplets of PMTs. This detector will have an area of about 2 $``$ 104 m2 for 1 TeV muons - similar to AMANDA-II - and is planned to be fully deployed by the end of 2004. In addition to these two advanced projects, there is an Italian initiative, NEMO, which finished site investigations at a location 80 km from Sicily and is now in the phase of prototype studies for a cubic kilometer detector . At the same time, also ANTARES, BAIKAL and NESTOR envisage larger arrays, possibly of cubic kilometer size.

There have been longstanding discussions about the best location for a future large neutrino telescope. What concerns geographic location, one detector on each hemisphere are necessary for full sky coverage. With respect to optical properties, water detectors in oceans seem to be favoured: although the absorption length of Antarctic ice at Amanda depths is nearly twice as long as in oceans (and about four times that of Baikal), ice is characterized by strong light scattering, and its optical parameters vary with depth. Light scattering leads to a considerable delay of Cherenkov photons. On the other hand ice does not suffer from the high potassium content of ocean water or from bioluminescence. These external light sources result in counting rates ranging from several tens of kHz to a few hundred kHz, compared to less than 500 Hz pure PM dark count rate in ice. Depth arguments favour oceans. The seabed at the NESTOR site is deepest (4 km), closely followed by NEMO (3.5 km). With only 2.5 km, the ANTARES ground is at about the same depth as the lowest AMANDA modules. Note, however, that the depth argument lost some of its initial strength after BAIKAL and AMANDA had developed reconstruction methods which effectively reject even the high background at shallow depths. Actually, the main advantage of great depths is the possibility to look higher above horizon, i.e. an increased angular acceptance. For water, a detector at greater depth suffers less from sedimentation of biomatter and from increased noise rates due to bio-luminescence. What counts most, at the end, are basic technical questions like deployment, or the reliability of the single components as well as of the whole system. Systems with a non-hierarchical structure like AMANDA (where each PM has its own 2 km cable to surface) will suffer less from single point failures than water detectors do. In the case of water, longer distances between detector and shore station have to be bridged. Consequently, not every PM can get its own cable to shore, resulting in a hierarchical system architecture. This drawback of water detectors may be balanced by the fact that they allow retrieval and replacements of failed components, as the BAIKAL group has demonstrated over many years.

Most likely, the present efforts will converge to two cubic kilometer detectors for very high energy neutrinos, ICECUBE at the South Pole and one in the Mediterranean.

## 4 Acoustic detection

Acoustic particle detection was proposed first in the fifties  and experimentally proven two decades later . A high energy particle cascade deposits energy into the medium via ionization losses, which is immediately converted into heat. The effect is a fast expansion, generating a bipolar acoustic pulse with a width of a few ten microseconds in water or ice (see fig.7). Transverse to the pencil-like cascade (diameter about 10 cm) the radiation propagates within a disk of about 10 m thickness (the length of the cascade) into the medium. The signal power peaks at 20 kHz where the attenuation length of sea water is a few kilometers, compared to a few tens of meters for light. Given a large initial signal, huge detection volumes can be achieved. Provided efficient noise rejection, acoustic detection might be competitive with optical detection at multi-PeV energies .

Present initiatives (see ) envisage combinations of acoustic arrays with optical Cherenkov detectors (NESTOR, ANTARES, ICECUBE) or the use of existing sonar arrays for submarine detection close to Kamchatka and in the Black Sea . Most advanced is AUTEC, a project using a very large hydrophone array of the US Navy, close to the Bahamas . The existing array of 52 hydrophones spans an area of 250 km2 and has good sensitivity between 1-50 kHz. It is expected to trigger on events above 100 EeV with a tolerable false alarm rate.

## 5 Radio detection

Electromagnetic showers generated by high energy electron neutrino interactions emit coherent Cherenkov radiation. Radio Cherenkov emission was predicted in 1962  and confirmed by recent measurements at SLAC and ANL . Electrons are swept into the developing shower, which acquires a negative net charge from the added shell electrons . This charge propagates like a relativistic pancake of 1 cm thickness and 10 cm diameter. Each particle emits Cherenkov radiation, with the total signal being the resultant of the overlapping Cherenkov cones. For wavelengths larger than the cascade diameter, coherence is observed and the signal rises proportional to $`E^2`$, making the method attractive for high energy cascades. The bipolar radio pulse has a width of 1-2 ns. In ice as well as in salt domes, attenuation lengths of several kilometers can be obtained, depending on the frequency band, the temperature of the ice, and the salt quality. Thus, for energies above a few tens of PeV, radio detection in ice or salt might be competitive or superior to optical detection .

A prototype Cherenkov radio detector called RICE is operating at the geographical South Pole . Twenty receivers and emitters are buried at depths between 120 and 300 m (fig.8). From the non-observation of very large pulses, a limit of about $`10^4E_\nu ^2`$ GeV-1 cm-2 s-1 sr-1 has been derived for energies above 100 PeV.

SALSA, a R&D project study for radio detection in natural salt domes, promises to get a limit about three orders of magnitude better . ANITA (ANtarctic Impulsive Transient Array) is an array of radio antennas planned to be flown at a balloon on an Antarctic circumpolar path in 2006 . From 35 km altitude it may record the radio pulses from neutrino interactions in the thick ice cover and monitor a really huge volume (see fig.9). The expected sensitivity from a 30 day flight is about $`10^7E_\nu ^2`$ GeV-1 cm-2 s-1 sr-1 at 10 EeV.

Most exotic is the Goldstone Lunar Ultrahigh Energy Neutrino Experiment, GLUE (see fig.10). It has searched for radio emission from extremely-high energy cascades induced by neutrinos or cosmic rays skimming the moon surface . Using two NASA antennas, an upper limit of $`10^4E_\nu ^2`$ GeV-1 cm-2 s-1 sr-1 at 100 EeV has been obtained.

## 6 Detection of neutrino energies via air showers

At supra-EeV energies, large extensive air shower arrays like the AUGER detector in Argentina  or the telescope array  may seek for horizontal air showers due to neutrino interactions deep in the atmosphere (showers induced by charged cosmic rays start on top of the atmosphere). Figure 11 explains the principle. AUGER consists of an array of water tanks going to span an area of more than 3000 km2 and will record the Cherenkov light of air-shower particles crossing the tanks. It is combined with telescopes looking for the atmospheric fluorescence light from air showers. The optimum sensitivity window for this method is at 1-100 EeV, the effective detector mass is between 1 and 20 Giga-tons, and the estimated sensitivity is of the order of $`10^8E_\nu ^2`$ GeV-1 cm-2 s-1 sr-1. An even better sensitivity might be obtained for tau neutrinos, $`\nu _\tau `$, scratching the Earth and interacting close to the array. The charged $`\tau `$ lepton produced in the interaction can escape the rock around the array, in contrast to electrons, and in contrast to muons it decays after a short path into hadrons. If this decay happens above the array or in the field of view of the fluorescence telescopes, the decay cascade can be recorded. Provided the experimental pattern allows clear identification, the acceptance for this kind of signals can be large. For the optimal energy scale of 1 EeV, the sensitivity might reach $`10^8E_\nu ^2`$ GeV-1 cm-2 s-1 sr-1. A variation of this idea is to search for tau lepton cascades which are produced by horizontal PeV neutrinos hitting a mountain and then decay in a valley between target mountain and an “observer” mountain – see fig.12 .

Already eight years ago, the Fly’s Eye collaboration , and more recently, the Japanese AGASA collaboration  have practiced the search mode of horizontal air showers. AGASA derived an upper limit of the order of 10-5 in the units given above \- only just one order of magnitude above some predictions for AGN jets and for topological defects.

Heading to higher energies leads to space based detectors monitoring larger volumes than visible from any point on the Earth surface. The projects EUSO  and OWL  foresee to launch large mirrors with optical detectors to 500 km height. The mirrors would look down upon the atmosphere and search for nitrogen fluorescence signals due to neutrino interactions (see fig.13). The monitored mass would be up to 10 Tera-tons, with an energy threshold of about 1010 GeV.

Finally, I mention the idea to detect the radio emission from cosmic ray and neutrino induced air showers with low frequency radio telescopes, as dicussed in .

## 7 Scenario for the next decade

The next ten years promise to be a particularly exciting decade for high energy neutrino astrophysics. Figures 14 and 15 sketch possible scenarios to move the frontier towards unprecedented sensitivities.

Like every estimate of time scales and expected physics performance, this scenario should be taken with caution. Notorious time delays in the realization of projects on the one hand, and possible new approaches on the other, will likely modify the evolution. In addition, new techniques have to be fully understood – their sensitivity to signals as well as the backgrounds! With AMANDA and BAIKAL, the optical underwater/ice technique has mastered this phase: downgoing muons and atmospheric neutrinos provided an invaluable calibration source. The energy range of acoustic or radio techniques as described in sections 4 and 5, however, is beyond the range covered by atmospheric neutrinos or normal cosmic rays. Calibration of the detectors in the absence of a surefire signal will be a challenge.

Figure 14 addresses the sensitivity to diffuse fluxes, i.e. integrating over the full angular acceptance of the detectors. The scale is set by the known flux of atmospheric neutrinos, by the bounds derived from observed fluxes of charged cosmic rays (W&B  and, with less stringent assumptions, the lower MPR curve ), by gamma rays (horizontal MPR line which assumes that cosmic rays are mostly confined in the cosmic source region and only gammas and neutrinos escape), and by specific model predictions . The figure shows two of the latter, one for the predicted flux of GZK neutrinos at ultra-high energies (see item b) above), the other for a model on neutrinos from Active Galactic Nuclei, peaking at TeV-PeV energies (Stecker and Salomon, SS ). The majority of the limits shown are published as “differential” limits, defining the flux sensitivity as the neutrino flux which gives at least one observed event per decade of energy per year (assuming negligible background). The limits published for AMANDA assume an $`E^2`$ flux. They come from two separate analyses, the one studying upward muons tracks , the other downward tracks of very high energy which are unlikely being due to muons generated in the atmosphere . Both analyses, however, properly account for the background of these atmospheric muons. The lines marked 1 and 2 extend over the energy range containing 90% of the events expected from an $`E^2`$ flux. For better illustration of the progress in time and over all the energy range, the AMANDA limits as well as the limits expected for ICECUBE have been translated to limits differentially per energy decade.

The present frontier is defined by TeV-PeV limits obtained by AMANDA and BAIKAL, and by PeV-EeV limits from the South Pole radio array RICE and the Japanese AGASA air shower array. Note that the Baikal/Amanda limits just reached a level sufficient to test (and actually to exclude) the model shown. The progress over the next 2 years will come from AMANDA and BAIKAL. After that, the Mediterranean telescopes – ANTARES and NESTOR – will start to contribute, flanked by AUGER and the ANITA balloon mission at high energies. This could result in an improvement of about two orders of magnitude over the full relevant energy range. Actually, five years from now a large variety of models might have been tested, including several predictions for neutrinos at GZK energies . Finally, in ten years from now, the TeV-PeV sensitivity will be defined by the cubic kilometer arrays at the South Pole and in the Mediterranean (marked as “km3”). At the high energy frontier, a SALSA-like experiment, and still higher, satellite detectors, might push the limit down by about three orders of magnitude compared to 2002.

Most likely, the first signal with clear signature will be a point source, possibly a transient signal which is easiest to identify. Figure 15 sketches a possible road until 2012. Best present limits are from MACRO, Super-Kamiokande (Southern sky) and AMANDA-B10 (Northern sky). This picture will not change until the Mediterranean detectors come into operation. AMANDA and ANTARES/NESTOR have the first realistic chance to discover an extraterrestrial neutrino source. The ultimate sensitivity for the TeV-PeV range is likely reached by the cubic kilometer arrays. This scale is set by many model predictions for neutrinos from cosmic accelerators or from dark matter decay. However, irrespective of any specific model prediction, these gigantic detectors, hundred times larger than AMANDA and thousand times larger than underground detectors, will hopefully keep the promise for any detector opening a new window to the Universe: to detect unexpected phenomena.

Acknowledgment: I thank S. Barwick, P. Gorham, M. Kowalski and J. Learned for stimulating and helpful discussions. Furthermore, I acknowledge useful information and suggestions obtained from J. Brunner, A. Ringwald, D. Semikoz, A. Tsirigotis and Sh. Yoshida. This paper is a written version of talks given at the 18th European Cosmic Ray Conference, Moscow 2002, and the 8th Topical Seminar on Innovative Particle and Radiation Detectors, Siena 2002. I thank L. Kuzmichev and F. Navarria, respectively, for their kind invitations and their support.


format: text

text:
# The Spin-Kinetic Temperature Coupling and the Heating Rate due to Lyman Alpha Scattering before Reionization: Predictions for 21cm Emission and Absorption

## 1 Introduction

After recombination of the primordial plasma at $`z1000`$, the Universe entered an age of darkness. The primordial density fluctuations continued to grow, until non-linear gravitational collapse led to the formation of the first dark matter halos in which the gas could collapse, and then radiatively cool. The dark ages ended at the time when the first stars formed in these halos, with a characteristic mass of $`10^6M_{\mathrm{}}`$ (Couchman & Rees, 1986; Tegmark et al., 1997; Abel et al., 2000; Bromm et al., 2002). These first stars, as well as any black holes or neutron stars formed subsequently which could accrete surrounding matter, produced photons capable of ionizing hydrogen (with wavelength $`\lambda <91.2\mathrm{nm}`$). As soon as the ionizing radiation produced in the sites of star formation was able to escape into the low-density intergalactic medium, the process of reionization started. During reionization, the ionized volume fraction of the universe increased gradually, until all the low-density regions were ionized. Recent observations of the Gunn-Peterson trough (Gunn & Peterson, 1965) at $`z6`$ (Becker et al., 2001; Fan et al., 2002, 2003), indicate that the end of reionization probably occurred at $`z6`$, although the epoch at which reionization began is still unknown (e.g., Gnedin, 2000; Miralda-Escudé, 2003).

The redshifted 21 cm line opens up a promising window for observing this period of cosmic history. This line is produced by the transition between the singlet and triplet hyperfine levels of the hydrogen atom at the electronic ground state. The line can be observed in absorption or emission against the Cosmic Microwave Background (CMB). The brightness temperature is given by

$$\delta T=\frac{T_sT_{CMB}}{1+z}(1e^\tau )(0.025\mathrm{K})\left(\frac{\mathrm{\Omega }_bh_0}{0.03}\right)\left(\frac{0.3}{\mathrm{\Omega }_{m0}}\right)^{1/2}\left(\frac{1+z}{10}\right)^{1/2}\frac{\rho _{HI}}{\overline{\rho }_H}\frac{T_sT_{CMB}}{T_s}.$$
(1)
where $`T_s`$ and $`T_{CMB}`$ are the hydrogen spin temperature and the CMB temperature, $`\tau `$ is the optical depth, $`\rho _{HI}`$ is the density of atomic hydrogen and $`\overline{\rho }_H`$ is the mean hydrogen density of the universe (we have used the approximations $`\tau 1`$ and $`\mathrm{\Omega }_{m0}(1+z)^31`$). Therefore, the detailed spectrum of the CMB at the long wavelengths of the redshifted 21cm, when observed with high angular resolution, contains precious information on the patchwork of density, spin temperature, and ionization inhomogeneities during the era when the intergalactic medium was still atomic (Hogan & Rees, 1979; Scott & Rees, 1990; Kumar, Padmanabhan, & Subramanian, 1995; Madau, Meiksin, & Rees, 1997; Shaver et al, 1999; Tozzi et al., 2000; Iliev et al., 2002).

The hydrogen spin temperature plays a key role in determining the amplitude of the 21cm emission or absorption signal against the CMB. At the same time, the evolution of the spin temperature depends on the gas kinetic temperature, and the way that the hyperfine structure level populations couple to the CMB and the kinetic temperatures. The adiabatic cooling of the hydrogen due to the expansion of the universe causes the gas temperature to decline as $`(1+z)^2`$, while the CMB temperature is decreasing only as $`1+z`$, so the gas cools below the CMB. At redshifts $`z100`$, the residual ionization left over from the epoch of recombination (Peebles, 1968) keeps the gas kinetic temperature at the CMB temperature owing to the heating produced by electron scattering of the CMB, but below this redshift this heating source becomes small and the gas temperature drops below the CMB. If the spin temperature can be coupled to the kinetic temperature at this stage, the atomic medium would be observed in absorption.

The coupling of the spin and kinetic temperatures can be achieved by atomic collisions or by scattering of Ly$`\alpha `$ photons. The coupling due to atomic collisions is small except at very high redshift, where observations of the redshifted 21cm radiation becomes more difficult, or for unrealistically high values of $`\mathrm{\Omega }_b`$ (Scott & Rees, 1990). A stronger coupling can be induced by Ly$`\alpha `$ photons, which should be produced as soon as the first stars are formed. The Ly$`\alpha `$ photons couple the spin and kinetic temperatures by being repeatedly scattered in the gas. The resonance scattering of a photon consists of a transition from $`n=1`$ to $`n=2`$, followed by the opposite downward transition. This scattering process can change the population of the hyperfine structure levels (Wouthuysen, 1952; Field, 1958). The probabilities for a transition back to the ground or excited hyperfine structure states depend on the slope of the radiation spectrum near the Ly$`\alpha `$ line center, or the “color temperature”. When the radiation spectrum reaches a steady state, this color temperature near the Ly$`\alpha `$ line center is equal to the gas kinetic temperature, so the Ly$`\alpha `$ scatterings introduce a thermal coupling between the spin and kinetic temperature. The result of these couplings is that the spin temperature is given by (Field, 1958, 1959)

$$T_S=\frac{T_{CMB}+y_\alpha T_k+y_cT_k}{1+y_\alpha +y_c},$$
(2)
where

$$y_\alpha \frac{P_{10}T_{}}{A_{10}T_k};y_c\frac{C_{10}T_{}}{A_{10}T_k}.$$
(3)
Here, $`T_k`$ is the kinetic temperature, $`T_{}=0.0682`$K is the hyperfine energy splitting, $`A_{10}=2.87\times 10^{15}\mathrm{s}^1`$ is the spontaneous emission coefficient of the 21cm line, and $`C_{10}=\kappa _{10}n_H`$ is the collisional de-excitation rate of the excited hyperfine level. The value of $`\kappa _{10}`$ ranges from $`2\times 10^{14}`$ to $`2.5\times 10^{10}\mathrm{cm}^3\mathrm{s}^1`$ for $`T`$ in the range $`11000`$ K and was tabulated in Allison & Dalgarno (1969).

The indirect de-excitation rate $`P_{10}`$ of the hyperfine structure levels is related to the total Ly$`\alpha `$ scattering rate $`P_\alpha `$ by $`P_{10}=4P_\alpha /27`$ (Field, 1958), which is given by $`P_\alpha =𝑑\nu cn_\nu \sigma (\nu )`$, where $`n_\nu `$ is the number density of photons per unit frequency, and $`\sigma (\nu )`$ is the cross section for Ly$`\alpha `$ scattering (Madau, Meiksin, & Rees, 1997, hereafter MMR).

Besides coupling the spin temperature and kinetic temperature of the gas, the scattering of Ly$`\alpha `$ photons may also heat up the gas. In this paper, we present the scattering rate and heating rate due to Ly$`\alpha `$ photons in a uniform atomic intergalactic medium as a function of the two parameters they depend on: the gas temperature and the Ly$`\alpha `$ scattering optical depth of the medium (or Gunn-Peterson optical depth).

There are two different ways in which photons propagating through the atomic medium may reach the Ly$`\alpha `$ line, and start being repeatedly scattered by hydrogen: (a) Photons that are emitted between the Ly$`\alpha `$ and Ly$`\beta `$ wavelengths by sources such as massive stars will be redshifted as they penetrate the atomic medium until they reach the Ly$`\alpha `$ wavelength. These photons enter the Ly$`\alpha `$ line from its blue wing, and will be called continuum photons in this paper. (b) Photons that are emitted at wavelengths shorter than Ly$`\beta `$, but still longer than the Lyman limit so that they can penetrate deep into the atomic medium, will redshift until they reach the Ly$`\beta `$ or a higher Lyman series line, and will then be converted into a Ly$`\alpha `$ photon when the hydrogen atom they excite decays first to the 2p state before reaching the ground state. These photons will be injected near the Ly$`\alpha `$ line center with a distribution determined by the Voigt function, and will be called injected photons. These two types of photons result in different scattering and heating rates, so results will be presented separately for both of them.

MMR estimated the heating rate due to Ly$`\alpha `$ photons by assuming that the average relative change in a Ly$`\alpha `$ photon energy in each scattering is the same as when the photon is scattered by an atom at rest: $`\mathrm{\Delta }E/E=h\nu _\alpha /m_Hc^2`$. Hence, they assumed that the total energy transfer rate is

$$\dot{E}=\frac{\mathrm{\Delta }E}{E}h\nu _\alpha P_\alpha .$$
(4)
At the thermalization rate $`P_{th}=27A_{10}T_{CMB}/4T_{}`$ (which is the rate $`P_\alpha `$ required to bring the spin temperature close to the kinetic temperature), the heating rate assumed by MMR is $`\dot{E}/k_B220\mathrm{K}\mathrm{Gyr}^1(1+z)/7`$, which would heat up the gas kinetic temperature above the CMB temperature in a fraction of a Hubble time at $`z6`$.

In the present investigation we calculate the heating rate of the continuum and injected Ly$`\alpha `$ photons by computing the background spectrum near the Ly$`\alpha `$ line using the Fokker-Planck method introduced by Rybicki & Dell’antonio (1994) (hereafter RD), and using energy conservation. We find that the heating rate due to Ly$`\alpha `$ scattering is much smaller than the estimate of eq. (4), and as a consequence Ly$`\alpha `$ photons do not constitute a significant heating source of the gas. Calculating the heating rate caused by Ly$`\alpha `$ photons as if every atom were at rest, without taking into account the thermal motions, is incorrect, and as we shall see the heating rate per unit volume is not proportional to the atomic density. We explain in the Appendix the reason why the argument used by MMR in their Appendix B to justify using eq. (4) for the heating rate is invalid. As pointed out by MMR, X-rays are important in pre-heating the atomic medium before it is reionized. In fact, our work leads to the conclusion that X-rays are the only important source of heating before the atomic medium is reached by an ionization front.

This paper is organized as follows: in §2 we describe our method of calculation. We present our result on scattering and heating rates in §3, explaining the physical reason why gas is not heated by Ly$`\alpha `$ photons as if every atom were at rest. We illustrate the consequences for future 21cm observations in §4, with an example for the evolution of the spin temperature with redshift, including also a model for X-ray heating. Our conclusions are summarized in §5. We adopt a flat $`\mathrm{\Lambda }`$CDM model with the following cosmological parameters: $`h_0=0.7`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_bh_0^2=0.021`$, $`Y_{\mathrm{He}}=0.24`$.

The code for solving the line profile and calculating the scattering and heating rates is available from the authors111Our code, named “Lyman Alpha scattering and Spin Temperature”, or LAST, can be downloaded at http://theory.kitp.ucsb.edu/$$xuelei/LAST..

## 2 Method of Calculation

Our goal in this work is to compute the Ly$`\alpha `$ scattering rate and the heating rate of hydrogen induced by this scattering, from which one can then compute the evolution of the spin temperature. The scattering and heating rates can be calculated once the background spectrum near the Ly$`\alpha `$ line is known. Let us consider the background of photons produced by constant, flat spectrum sources in an expanding universe (in practice the intrinsic source spectrum does not matter because Ly$`\alpha `$ scattering is important only on a very narrow range of frequencies around Ly$`\alpha `$ , in which the intrinsic source spectrum can be considered to be constant). In the absence of scattering, the photons would simply redshift and would produce a photon number density per unit frequency $`n(\nu )=n_0`$ that is independent of frequency near Ly$`\alpha `$ . The presence of scatterings by hydrogen introduces a feature in the spectrum around the Ly$`\alpha `$ line, which has a constant shape $`n(\nu )`$ once a steady-state is reached. Given this spectrum, the rate at which Ly$`\alpha `$ photons are being scattered is given by the integral

$$P_\alpha =c𝑑\nu n(\nu )\sigma (\nu ).$$
(5)
Here, $`\sigma (\nu )`$ is the Ly$`\alpha `$ scattering cross section as a function of frequency. The cross section has a thermal width $`\mathrm{\Delta }\nu _D`$ determined by the temperature $`T`$ of hydrogen atoms of mass $`m_H`$,

$$\mathrm{\Delta }\nu _D=\frac{\sqrt{2k_BT/m_H}}{c}\nu _\alpha ,$$
(6)
where $`\nu _\alpha `$ is the Ly$`\alpha `$ line central frequency and $`k_B`$ is the Boltzmann constant. The cross section at the line center is

$$\sigma _0=\frac{\pi e^2}{m_ec}f_{12}[\mathrm{\Delta }\nu _D]^1,$$
(7)
where $`f_{12}=0.416`$ is the oscillator strength for the Ly$`\alpha `$ transition. The dependence of the cross section on frequency is given by the normalized Voigt function,

$$\varphi (x)=\frac{a}{\pi ^{3/2}}_{\mathrm{}}^{\mathrm{}}𝑑y\frac{e^{y^2}}{(xy)^2+a^2},$$
(8)
where $`x=(\nu \nu _\alpha )/\mathrm{\Delta }\nu _D`$, the Voigt parameter is $`a=A_{21}/(8\pi \mathrm{\Delta }\nu _D)`$, and $`A_{21}=6.25\times 10^8`$ Hz is the Einstein spontaneous emission coefficient. We will also use the photon intensity $`J(\nu )=cn(\nu )/(4\pi )`$. With these definitions, the scattering rate is

$$P_\alpha =4\pi \sigma _0\mathrm{\Delta }\nu _DJ(x)\varphi (x)𝑑x.$$
(9)
The heating rate of the gas induced by Ly$`\alpha `$ scattering can be simply obtained by requiring conservation of energy. The scatterings cause the cosmic background to have a narrow spectral feature around the Ly$`\alpha `$ wavelength. Were the scatterings to cease at a certain instant, this feature would simply be redshifted to longer wavelengths. Instead, when a steady-state spectrum is reached, the feature is fixed at a constant wavelength. The reduction in energy per unit volume of the background due to this feature is $`\mathrm{\Delta }ϵ=𝑑\nu [n_0n(\nu )]h\nu 𝑑\nu `$ (for the purpose of this discussion we ignore the photons injected at the Ly$`\alpha `$ wavelength by recombinations, so that in the absence of scatterings, $`n(\nu )=n_0`$; injected photons will be fully considered later). Over an interval of time $`dt`$, the energy per unit volume that the background photons must lose to keep the spectral feature at a fixed wavelength is $`\mathrm{\Delta }ϵ(d\nu /\nu )=\mathrm{\Delta }ϵHdt`$, where $`d\nu `$ is the change in frequency due to redshift over the interval $`dt`$, and $`H`$ is the Hubble constant at time $`t`$. The energy lost by the photons must be gained by the hydrogen atoms, so the gas heating rate per unit volume is

$$\mathrm{\Gamma }=H\mathrm{\Delta }ϵ=H\left[n_0n(\nu )\right]h\nu 𝑑\nu =\frac{4\pi Hh\mathrm{\Delta }\nu _D}{c}\left[J_0J(\nu )\right]\nu 𝑑\nu .$$
(10)
We use these equations in §3 to compute the scattering and heating rates once the line spectrum has been calculated. An alternative derivation of the heating rate in equation (10) is given in the Appendix.

The spectrum of the cosmic background near Ly$`\alpha `$ is determined by the radiative transfer equation for resonant scattering for a homogeneous and isotropic expanding universe, which is (RD)

$$\frac{J(x,\tau )}{\tau }=\varphi (x)J(x,\tau )+\gamma _S\frac{J(x,\tau )}{x}+R_{II}(x,x^{})J(x^{},\tau )𝑑x^{}+C(\tau )\psi (x),$$
(11)
where $`\tau =(cn_1\sigma _0)t`$ is the mean free time between scatterings for a photon at the line center, $`t`$ is the physical time, and $`n_1n_H`$ is the number density of hydrogen atoms in the ground state. The Sobolev parameter is $`\gamma _S=\tau _{GP}^1`$, where $`\tau _{GP}=(\pi e^2n_1f_{12})/(Hm_e\nu _\alpha )`$ is the Gunn-Peterson optical depth. The first term in the right-hand-side of equation (11) describes the removal of photons that are scattered, and the second term accounts for the redshift. The third term represents the addition of the photons redistributed over frequency after scattering (Rybicki & Hummer, 1992), and the fourth term is due to the injection of new Ly$`\alpha `$ photons with a spectral profile $`\psi (x)`$.

In the Fokker-Planck approximation, the scattering problem is treated as diffusion in the approximation that the width of the spectral feature is large compared to the frequency change in typical scatterings. The redistribution term then becomes (RD)

$$R_{II}(x,x^{})J(x^{},t)𝑑x^{}=\varphi (x)J(x)+\frac{1}{2}\frac{}{x}\left[\varphi (x)\frac{J(x)}{x}+2\eta \varphi (x)J(x)\right],$$
(12)
where the recoil term is included, with $`\eta =h\nu _\alpha /(c\sqrt{2m_Hk_BT})`$. This approximation conserves photon number. The equation is then reduced to \[see eq. (34) of RD\]

$$\frac{J}{\tau }=\frac{1}{2}\frac{}{x}\left[\varphi \frac{J}{x}+2\eta \varphi J\right]+\gamma _S\frac{J}{x}+C(t)\psi (x).$$
(13)
When the equilibrium state is reached, $`J/\tau =0.`$

Since this equation is linear, we can write its solution as $`J=J_c+J_i`$, where $`J_c`$ is the spectrum of the continuum photons, and $`J_i`$ is the spectrum of the recombination photons injected near the line center.

For the continuum photons, there is no injection and so $`C(t)=0`$. Equation (13) is reduced to

$$\varphi (x)\frac{d}{dx}J_c(x)+2\left[\eta \varphi (x)+\gamma _S\right]J_c(x)=A,$$
(14)
where $`A`$ is a constant of integration determined by the boundary condition: as $`x\mathrm{}`$, $`\varphi (x)0`$, $`J_cJ_{c0}`$, where $`J_{c0}=(4\pi /c)n_{c0}`$ is the unperturbed intensity. Therefore, $`A=2\gamma _SJ_0`$. The equation is then solved as an initial value problem of an ordinary differential equation with the condition $`J_c(x)J_{c0}`$ for $`x\mathrm{}`$, required by photon number conservation.

In Fig. 1, we plot the solution for a typical case: $`T=2.65`$K and $`\gamma _S=1.44\times 10^6`$, corresponding to the physical conditions of unheated gas at $`z=10`$ in a flat $`\mathrm{\Lambda }`$CDM model. An absorption feature is produced near the center of the line.

For the injected photons, we assume a constant source $`C(t)=C`$. This should generally be an accurate approximation because the timescale for the evolution of the production rate of Ly$`\alpha `$ photons from sources is likely be of order the Hubble time, which is much longer than the relaxation time, which is of order the time to redshift across the width of the Ly$`\alpha `$ spectral profile of the cosmic background. We also inject photons with a single frequency at the Ly$`\alpha `$ line center, i.e., $`\psi (x)=\delta (x)`$. Because photons change randomly in frequency over a thermal width at each scattering, and the number of scatterings each photon undergoes is very large, injecting the photons with a thermal width instead does not make any difference to the final profile \[we have verified this by obtaining solutions with $`\psi (x)`$ given by a Voigt profile, which did not result in any significant change to the resulting spectral profile\]. The solution is

$$\varphi (x)\frac{d}{dx}J_i(x)+2\left[\eta \varphi (x)+\gamma _S\right]J_i(x)=2C\mathrm{\Theta }(x)+B,$$
(15)
where $`\mathrm{\Theta }(x)`$ is the Heaviside step function. The value of the integration constant $`B`$ can again be determined from boundary conditions: as $`x+\mathrm{}`$, we must now have $`J_i(x)0`$, from which we obtain $`B=2C`$. The profile of the line can again be obtained by solving eq. (15), with the initial condition $`J_i(x=\mathrm{})=J_{i0}`$, where $`J_{i0}`$ is the photon intensity due to the injected photons. The result is shown in Fig. 2, for two different temperatures; we see that some of the photons injected at the line center are scattered to the right (blue) side.

We note that the solution for the spectral profile of continuum and injected photons has a characteristic width in the variable $`x=\nu /\mathrm{\Delta }\nu _D`$ that is substantially larger than unity, implying that the Fokker-Planck approximation we have used should be accurate.

## 3 Results

The heating and scattering rates are of course proportional to the background intensity. During the epoch of reionization, at least one ionizing photon has to be emitted for every baryon in the universe, and more if the typical baryon can recombine many times. The number of Ly$`\alpha `$ photons emitted is of the same order for sources without a large Lyman break, and can be much larger if many ionizing photons are internally absorbed. We use the intensity corresponding to one photon per hydrogen atom in the universe as a fiducial value:

$$\stackrel{~}{J}_0(z)=\frac{n_Hc}{4\pi \nu _\alpha }=1.65\times 10^{13}(1+z)^3\left(\frac{\mathrm{\Omega }_bh_0^2}{0.02}\right)\mathrm{cm}^2\mathrm{s}^1\mathrm{sr}^1\mathrm{Hz}^1,$$
(16)
where $`n_H`$ is the hydrogen number density. For reference to other work, we quote also the value of the energy intensity of this unit: $`\stackrel{~}{J}_0h\nu _\alpha =2.69\times 10^{24}(1+z)^3(\mathrm{\Omega }_bh_0^2/0.02)\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1\mathrm{sr}^1\mathrm{Hz}^1`$. The heating rate in equation (10) can be reexpressed as a heating rate per hydrogen atom and per Hubble time. For the continuum photons, it is:

$$\frac{\mathrm{\Gamma }_c}{Hn_Hk_B}=\frac{h\mathrm{\Delta }\nu _D}{k_B}\frac{J_{c0}}{\stackrel{~}{J}_0}𝑑x\left[1\frac{J_c(x)}{J_{c0}}\right]\sqrt{2TT_0}\frac{J_{c0}}{\stackrel{~}{J}_0}I_c,$$
(17)
where we have defined $`I_c=𝑑x[1J_c(x)/J_{c0}]`$, and $`T_0=(h\nu _\alpha )^2/(k_Bm_Hc^2)=1.29\times 10^3`$ K. For the example shown in Figure 1, the integral is the shaded area, $`I_c=16.7`$. It is clear from equation (17) that if the number of photons in the background is not greater than one per hydrogen atom, the heating rate due to Ly$`\alpha `$ scattering is not very large because $`T_0T`$. In the example of $`T=2.65`$ K at $`z=10`$, equation (17) yields $`\mathrm{\Gamma }_c/(Hn_Hk_B)=1.4\mathrm{K}(J_{c0}/\stackrel{~}{J}_0)`$.

Similarly, for the injected photons with a $`\delta `$-function profile, the heating rate is

$$\frac{\mathrm{\Gamma }_i}{Hn_Hk_B}=\frac{h\mathrm{\Delta }\nu _D}{k_B}\frac{J_{i0}}{\stackrel{~}{J}_0}\left\{_{\mathrm{}}^0𝑑x\left[1\frac{J_i(x)}{J_{i0}}\right]_0^+\mathrm{}𝑑x\frac{J_i(x)}{J_{i0}}\right\}\sqrt{2TT_0}\frac{J_{i0}}{\stackrel{~}{J}_0}I_i.$$
(18)
We note that, for the more exact case in which photons are injected with a Voigt profile in frequency, the heating rate is obtained by subtracting from $`J_i(x)`$ the steady-state solution of the background with no scattering (shown in Fig. 3 of RD). However, one can easily prove that the result is exactly the same as what is obtained by subtracting instead the step function spectrum that results from a $`\delta `$ function of the injected photons.

In Fig. 2, we show the spectrum $`J_i`$ for two different temperatures: $`T_k=2.65`$K and $`T=1`$K, for $`\gamma _S=1.44\times 10^6`$. In the $`T_k=2.65`$K case, we have $`\mathrm{\Gamma }_i/(Hn_Hk_B)=0.104\mathrm{K}(J_{i0}/\stackrel{~}{J}_0)`$. This negative heating rate means that the scattering of these injected photons actually cools the gas, because some of the photons are scattered to the blue side of the spectrum and gain energy from the hydrogen gas. At lower gas temperatures, as in the $`T_k=1`$K case shown in the figure, $`\mathrm{\Gamma }_i/(Hn_Hk_B)=0.40\mathrm{K}(J_{i0}/\stackrel{~}{J}_0)`$, the photons heat the gas. However, at $`z=10`$ the temperature of adiabatically expanding gas is 2.65 K, and there is no known mechanism which could cool the gas below this temperature, so the net effect of the injected photons must be to cool the gas. In any case, both the heating rate due to continuum photons and the cooling rate due to injected photons are very small.

To calculate the spin temperature from equation (2), we also need to know the scattering rate of the Ly$`\alpha `$ photons in equation (9), which can be reexpressed as

$$P_\alpha =H\tau _{GP}\frac{J(x)}{\stackrel{~}{J}_0}\varphi (x)𝑑x.$$
(19)
We define

$$S_c\frac{J_c(x)}{J_{c0}}\varphi (x)𝑑x,S_i\frac{J_i(x)}{J_{i0}}\varphi (x)𝑑x.$$
(20)
We plot the integrals $`I_c`$, $`I_i`$, $`S_c`$, and $`S_i`$ as a function of $`\gamma _S=\tau _{GP}^1`$, for several different temperatures, in Figures 34. Our code calculates the values of these integrals for any given gas temperature and optical depth.

The intensity corresponding to the thermalization rate required to bring the spin temperature down to the kinetic temperature (see the introduction) is

$$\frac{J_0}{\overline{J}_0}=\frac{P_\alpha }{H\tau _{GP}S}=\frac{27A_{10}T_{CMB}\gamma _S}{4T_{}HS},$$
(21)
where the subscripts $`c`$ or $`i`$ can be applied. For the case $`\gamma _S=1.44\times 10^6`$ and $`T=2.65`$ K at $`z=10`$, $`S_i`$ or $`S_c`$ are close to one and we find $`J_0/\overline{J}_00.3`$. Thus, from our earlier result that the heating rate in this case is $`\mathrm{\Gamma }_c/(Hn_Hk_B)=1.4\mathrm{K}(J_{c0}/\stackrel{~}{J}_0)`$ for the case of continuum photons, we find that over the time that continuum Ly$`\alpha `$ photons bring the spin temperature down to the gas temperature, they should heat the gas by only $`0.5`$ K, and injected photons should cool the gas by $`0.03`$ K.

The reason why the heating rate is much less than the value obtained by assuming that each atom loses the same average energy at each scattering as if it were at rest is that the intensity of the background spectrum drops steeply at the Ly$`\alpha `$ line center, as shown in Figures 1 and 2. The vast majority of the scatterings occur within a few thermal widths of the line center, where the color temperature of the spectrum (reflecting its slope) is almost equal to the kinetic temperature. When these two temperatures are the same, there obviously cannot be any heat exchange between photons and atoms because they are in thermodynamic equilibrium. The heating from continuum photons occurs because the color temperature is slightly larger than the kinetic temperature. Injected photons have a temperature slightly lower than the kinetic temperature, so they cool the gas. Cooling is caused by photons that are on the red side of the line center because these photons are more likely to be scattered by atoms moving in the direction opposite to the photon (for which the Doppler effect can shift the photon frequency to the line center), and these atoms will be slowed down by the scattering. When the color temperature is below the kinetic temperature, the cooling effect due to the preferential scattering of photons moving in the opposite direction than the atoms is greater than the heating effect due to recoil from photons moving in a random direction.

The heating or cooling rates per unit volume are not proportional to the atom density. If the atom density is increased, not only the number of scatterings increases, but also the color temperature is brought closer to the kinetic temperature because of the increased rate of interaction between atoms and photons, and the average energy exchanged between them at each scattering is reduced. The heating rate changes only with the width of the background spectral feature that determines the integrals $`I_c`$ and $`I_i`$, according to the energy conservation argument used in §2.

In the calculation of the spin temperature evolution from equation (2), we assume that the color temperature of the Ly$`\alpha `$ photons is equal to the gas kinetic temperature. This is true at the line center in the limit of large optical depth (Field 1958, 1959), but the slope of the photon distribution varies outside the line center. In the more general case, equations (2)-(3) can be generalized as

$$y_\alpha T_k𝑑xy_\alpha (x)T_\alpha (x)=\frac{T_{}}{A_{10}}\frac{4}{27}𝑑P_\alpha (x)$$
(22)
$$y_\alpha 𝑑xy_\alpha (x)=\frac{T_{}}{A_{10}}\frac{4}{27}𝑑P_\alpha (x)/T_\alpha (x)$$
(23)
where the color temperature $`T_\alpha `$ at $`x`$ is defined by

$$\frac{J^{}(x)}{J(x)}\frac{h\nu }{k_BT_\alpha (x)}\frac{h\nu _0}{k_BT_\alpha (x)},$$
(24)
and

$$dP_\alpha (x)=4\pi \sigma _0\mathrm{\Delta }\nu _DJ(x)\varphi (x)dx$$
(25)
Because the vast majority of scatterings occur near the line center, where the color temperature and kinetic temperature are very nearly the same, this generalization does not significantly change the result for the spin temperature. We have found by a numerical test that this correction affects the spin temperature by only a small fraction of a percent.

## 4 Simple Models for the Thermal Evolution

We will now consider a simple model where the emissivity of Ly$`\alpha `$ photons turns on at a redshift $`z_i`$ and increases linearly with redshift thereafter. We will also consider the effect of X-ray emission from the same star formation regions, which dominates the heating rate of the gas. This will illustrate a plausible thermal history of the gas and the expected strength of the absorption or emission signal in the redshifted 21cm line.

The kinetic temperature of the gas in the expanding Universe evolves as

$$(1+z)\frac{dT_k}{dz}=2T_k\frac{2}{3}\frac{(\mathrm{\Gamma }_{tot}\mathrm{\Lambda }_{tot})}{(Hnk_B)},$$
(26)
where $`\mathrm{\Gamma }_{tot}`$ and $`\mathrm{\Lambda }_{tot}`$ are the total heating and cooling rates, respectively, and $`n`$ is the total gas number density.

In the absence of heating and cooling, the gas temperature decreases adiabatically with $`T(1+z)^2`$. We use the code RECFAST 222http://www-cfa.harvard.edu/$$sasselov/rec (Seager, Sasselov, & Scott, 1999; Seager, Sasselov & Scott, 2000) to calculate the temperature evolution of the gas before the first stars and quasars turn on. We note that at $`z40`$ the spin temperature drops below the CMB temperature due to the collisional coupling (in Fig. 6, this drop is seen at the highest redshift end). The corresponding absorption is, however, difficult to observe at the very long wavelengths corresponding to this redshift. At lower redshifts, the spin temperature is practically equal to the CMB temperature until the first Ly$`\alpha `$ photons are produced.

### 4.1 Ly$`\alpha `$ emissivity

If the comoving photon emissivity (defined as the number of photons emitted per unit comoving volume, time and frequency) at the Ly$`\alpha `$ line from stars at redshift $`z`$ is $`ϵ(\nu _\alpha ,z)`$, then the comoving Ly$`\alpha `$ photon intensity at redshift $`z`$ is given by

$$J_{c0}(z)=\frac{1}{4\pi }(1+z)^{\alpha _S}_z^{z_{max}}\frac{c}{H(z^{})}ϵ(\nu _\alpha ,z^{})(1+z^{})^{(1+\alpha _S)}𝑑z^{},$$
(27)
where the emissivity spectrum is assumed to be $`ϵ(\nu )\nu ^{\alpha _S1}`$, and $`z_{max}`$ is determined by the condition that only the photons emitted at frequencies lower than the Ly$`\beta `$ line will not interact with the intergalactic medium:

$$(1+z_{max})=\frac{\nu _\beta }{\nu _\alpha }(1+z)=\frac{32}{27}(1+z)$$
(28)
To discuss the plausible range of values of the comoving emissivity, it is convenient to introduce the unit

$$\stackrel{~}{ϵ}\frac{n_HH}{\nu _\alpha }=\frac{4\pi }{c}H\stackrel{~}{J}_0.$$
(29)
Thus, if $`ϵ=\stackrel{~}{ϵ}`$, one photon is being emitted per hydrogen atom and per Hubble time, per unit $`\mathrm{log}(\nu )`$.

During the epoch of reionization, at least one photon must be emitted for every baryon in the universe. For a homogeneous medium, the recombination rate $`\alpha n_e^2`$ is equal to the Hubble rate, $`H`$, at $`z6`$; hence, if reionization takes place during the interval $`20z6`$, a few photons should be emitted for every baryon to ionize them and to compensate for recombination. The number of recombinations can be increased due to the clumping factor of the ionized intergalactic gas (see, however, Miralda-Escudé, Haehnelt & Rees, 2000 for a discussion of why the clumping factor of ionized gas is probably not much larger than unity during the epoch of reionization). If these ionizing photons are produced by star-forming regions, then the photon emissivity per unit $`\mathrm{log}(\nu )`$ between Ly$`\alpha `$ and the Lyman limit should be 3 to 5 times larger than the emissivity of ionizing photons, owing to the average Lyman break of the atmospheres of massive stars with a Salpeter Initial Mass Function (e.g., Bruzual & Charlot, 1993). The photons between Ly$`\alpha `$ and the Lyman limit may have a larger escape fraction from the star formation regions than the ionizing photons (if some of the absorption is due to local hydrogen with low dust content), in which case the ratio of the emissivity between Ly$`\alpha `$ and Lyman limit to the ionizing emissivity would be further increased. This suggests that the emissivity near Ly$`\alpha `$ should be at least $`ϵ10\stackrel{~}{ϵ}`$.

The total number of ionizing photons emitted is also related to the metallicity produced by massive stars. To produce one ionizing photon per baryon, the mean metallicity of the universe must be raised by $`10^5`$, or $`5\times 10^4Z_{}`$ (e.g., Madau & Shull, 1996). If the typical metallicity in the Ly$`\alpha `$ forest and damped Ly$`\alpha `$ systems at $`z3`$, which is $`Z10^2Z_{}`$, is not to be exceeded, one cannot produce more than $`20`$ ionizing photons per baryon, implying that $`ϵ100\stackrel{~}{ϵ}`$ at Ly$`\alpha `$ .

As a simple model of the Ly$`\alpha `$ background evolution, we assume that the Ly$`\alpha `$ emissivity increases linearly from 0 at some initial redshift $`z_i`$,

$$ϵ(z)=\{\begin{array}{cc}0;\hfill & z>z_i\hfill \\ ϵ_r\frac{z_iz}{z_iz_r};\hfill & z_r<z<z_i\hfill \end{array}$$
(30)
We will use two values for the peak emissivity at Ly$`\alpha `$ , $`ϵ_r`$, reached at $`z_r`$: $`ϵ_r=10\stackrel{~}{ϵ}`$ and $`ϵ_r=100\stackrel{~}{ϵ}`$, based on the previous discussion. Assuming also a flat spectrum, $`\alpha _S=0`$, the Ly$`\alpha `$ intensity due to the continuum photons is then

$$\frac{J_{c0}(z)}{\stackrel{~}{J}_0}=\frac{ϵ_r}{\stackrel{~}{ϵ}}g(1+z),$$
(31)
where the function $`g`$ is

$$g(x)=\frac{x^{3/2}}{z_iz_r}\left[\frac{2}{3}(1+z_i)\left(x^{3/2}x_<^{3/2}\right)2\left(x^{1/2}x_<^{1/2}\right)\right],z<z_i$$
(32)
and

$$x_<=\mathrm{min}[(1+z_i),\left(\frac{32}{27}x1\right)].$$
(33)
The resulting proper Ly$`\alpha `$ intensity is plotted as the function $`(1+z)^3J_{c0}/\stackrel{~}{J}_0`$ in Fig. 5, for $`z_r=6`$ and $`ϵ_r/\stackrel{~}{ϵ}=10`$, for two cases of initial redshift: $`z_i=20`$ and $`z_i=12`$. The proper Ly$`\alpha `$ intensity increases rapidly after $`z_i`$ and then decreases as $`z`$ decreases due to the $`(1+z)^3`$ factor.

We will assume that $`J_{i0}=J_{c0}`$. The injected photons have a comparable intensity to the continuum ones, because the photons emitted between Ly$`\beta `$ and the Lyman limit are similar to those emitted between Ly$`\alpha `$ and Ly$`\beta `$, and a large fraction of the Lyman continuum photons will eventually result in an injected Ly$`\alpha `$ photon when the atom that they ionize recombines again.

### 4.2 Heating by soft X-rays

Since the scattering of Ly$`\alpha `$ photons produces a negligible heating rate for the gas, the main heating mechanism should be soft X-ray photons. The soft X-rays can be emitted by a number of sources associated with star-forming regions, such as X-ray binaries, supernova remnants and QSOs.

X-ray photons have a photoionization cross section that is much smaller than that of the ionizing photons near the hydrogen ionization potential, so they can penetrate deep into the neutral regions of the intergalactic medium. There, they photoionize hydrogen or helium, and the high-energy electron that is produced can then transfer a fraction of its energy to the gas by Coulomb collisions with other electrons and protons. Once the energy has been thermalized among the charged particles at a temperature low enough to avoid losses by collisional excitation of hydrogen, the thermal energy can be transferred to the hydrogen atoms in slow collisions involving dipole interaction. The fraction of the X-ray energy that can be converted to the gas thermal energy (rather than being used in collisional ionizations and excitations of hydrogen or helium) depends on the fractional ionization of the gas. Using the RECFAST code, we found that the initial ionized fraction from the recombination era (before more ionizations are caused by the X-rays) is $`2\times 10^4`$, which implies that a fraction of 14% of the X-ray energy is used to heat the gas (Shull & Van Steenberg, 1985). As the fraction of free electrons increases, this fraction also increases, but for simplicity we shall assume it a constant.

We parameterize the emissivity in X-rays in terms of the fraction of the energy that is emitted in X-rays compared to the energy emitted at Ly$`\alpha `$ per unit $`\mathrm{log}\nu `$, which we designate as $`\alpha _x`$. The X-ray heating rate is then

$$\mathrm{\Gamma }_x(z)=1.4\times 10^3\left(\frac{\alpha _x}{0.01}\right)hϵ(z),$$
(34)
where $`ϵ(z)`$ is given in Eq. (30), and $`h`$ is the Planck constant. At $`z=6`$, the heating rate due to X-rays in this model is

$$\frac{\mathrm{\Gamma }_x}{n_HHk_B}=1.7\times 10^3\left(\frac{\alpha _x}{0.01}\right)\frac{ϵ_r}{10\stackrel{~}{ϵ}}\mathrm{K}.$$
(35)
Therefore, by emitting only $``$ 1% of the energy as X-rays, the sources can rapidly heat the intergalactic medium to temperatures above the CMB.

### 4.3 Results for thermal evolution and the 21 cm Antenna Temperature

The evolution of the gas kinetic and spin temperatures in the absence of X-ray heating is plotted in Fig. 6, for two values of the initial redshift ($`z_i=20`$ and $`z_i=12`$), and two values of the emissivity ($`ϵ_r/\stackrel{~}{ϵ}=10`$ and $`100`$).

In the absence of X-rays, the heating due to Ly$`\alpha `$ scattering is very low and the kinetic temperature of the gas deviates only slightly from the adiabatic curve at late times. For our assumed emissivity of 10 to 100 photons per baryon per Hubble time, the spin temperature drops to values very close to the kinetic temperature, producing a large absorption signal on the CMB of up to $`200`$ mK, as seen in Figure 7.

Figure 8 shows the thermal evolution and the resulting antenna temperature when we include X-ray heating with $`\alpha _x=0.01`$, for the case $`z_i=12`$ and $`ϵ_r/\stackrel{~}{ϵ}=10`$. With only 1% of the energy emitted as X-rays, the kinetic temperature of the gas increases rapidly when the first sources turn on, leaving only a brief period of absorption in the beginning of small amplitude, which then turns to emission, with an amplitude of less than 20 mK.

## 5 Conclusion

We have calculated the Ly$`\alpha `$ photon scattering and heating rates in the atomic intergalactic medium, including both redshifted continuum photons emitted between Ly$`\alpha `$ and Ly$`\beta `$, and injected photons produced by recombinations or excitations by Ly$`\beta `$ and higher Lyman series photons. The heating rate of the gas due to the scattering of these Ly$`\alpha `$ photons is very small and essentially negligible. In the absence of other heating sources, the gas would continue to cool down almost adiabatically while the spin temperature would rapidly drop to nearly the kinetic temperature given a reasonable emission of photons during the reionization epoch. This would produce a 21cm absorption against the CMB as large as $``$ 200 mK, at the average hydrogen density of the universe.

Of course, if the absorption signal can be observed at high angular and frequency resolution, the antenna temperature should fluctuate around its mean value owing to density fluctuations, as well as spin temperature fluctuations that would be produced by spatial variations in the history of the photon emissivity between Ly$`\alpha `$ and the Lyman limit (e.g.Tozzi et al., 2000; Ciardi & Madau, 2003). At the same time, during the epoch of reionization, there should be a gradual increase of the fraction of the volume that is ionized (e.g., Miralda-Escudé, 2003), and the average antenna temperature should then be multiplied by the fraction of baryons that remain in atomic form. Because this fraction of baryons should decrease gradually, any 21cm absorption or emission should also vary smoothly with redshift.

Unfortunately, the large absorption signal predicted in the absence of X-ray heating is unlikely to be present. It is very likely that significant X-ray emission is produced within only a few million years of the formation of the first stars in the universe: some of the stars may be in binaries which will evolve into high-mass X-ray binaries, and supernova explosions can also emit soft X-rays (see e.g. Oh, 2001; Venkatesan, Giroux, & Shull, 2001). Photons with energies as low as $`0.2`$ keV can propagate as far as 3 comoving Mpc over the atomic intergalactic medium at $`z10`$ before they are absorbed by hydrogen or helium.

There may still be some ways in which certain regions of the intergalactic medium may be illuminated by Ly$`\alpha `$ photons from the first stars, without receiving soft X-rays as well. For example, if stars form in a region with a high column density of hydrogen and helium with very little dust, and the hydrogen and/or helium is not fully ionized for a long time, the photons between Ly$`\alpha `$ and the Lyman limit from stars would escape, while the ionizing and soft X-ray photons might all be locally absorbed. In particular, there is the possibility that singly ionized helium remains self-shielded as massive stars ionize the gas around them. We note from Figure 6 that the time required to bring down the spin temperature can be very short ($`10^8`$ years or less) for typical values of the Ly$`\alpha `$ background intensity, and it is possible that a complete absorption of soft X-rays from the first star-forming regions would persist long enough to allow for a strong 21cm absorption signal in local regions of the intergalactic medium, when they are illuminated by the first light from stars forming in their vicinity. Nevertheless, it seems likely that any decrease of the spin temperature will be short-lived, and that X-ray heating will cause the 21cm signal to switch to emission a short time after the first stars are formed.

Future observations of the 21cm signal in absorption and emission at high angular and frequency resolution have an enormous potential for revealing the details of the reionization epoch. Fluctuations in the antenna temperature arise from both gas density and spin temperature variations. The density fluctuations will inform us of the state of gravitational evolution of primordial fluctuations at the highest attainable redshifts after recombination, while the spin temperature fluctuations will reveal the spatial distribution of the first sources of ultraviolet and X-ray emission in the universe.

This work was supported in part by NSF grants NSF-0098515 and PHY99-07949.

## Appendix

MMR assumed that the heating rate per unit volume of the hydrogen gas due to Ly$`\alpha `$ scattering is equal to the recoil energy given to a hydrogen atom at rest, $`(h\nu _\alpha )^2/(m_Hc^2)`$, times the scattering rate:

$$\mathrm{\Gamma }_{MMR}=\frac{(h\nu _\alpha )^2}{m_Hc^2}n_HP_\alpha =\frac{(h\nu _\alpha )^2}{m_Hc^2}\mathrm{\hspace{0.17em}4}\pi n_H\sigma _0\mathrm{\Delta }\nu _DJ(x)\varphi (x)𝑑x.$$
(1)
This heating rate was justified from equation (13) for the time evolution of the photon intensity $`J`$. MMR assumed that the rate of heat transfer to the atoms is determined by the change in the spectrum due to the recoil effect, which is the second term in the right-hand-side of equation (13):

$$\left(\frac{J}{\tau }\right)_{recoil}=\frac{(\eta \varphi J)}{x}.$$
(2)
They then computed the heating rate by integrating the variation in photon energy density per unit physical time over frequency:

$$\mathrm{\Gamma }_{MMR}=4\pi h\mathrm{\Delta }\nu _Dn_H\sigma _0𝑑xx\frac{(\eta \varphi J)}{x}=(Hn_Hk_B)T_0\tau _{GP}\frac{J_{c0}}{\stackrel{~}{J}_0}S_c,$$
(3)
where the integral $`S_c`$ is given in equation (20), and $`\tau _{GP}`$ is the Gunn-Peterson optical depth (all symbols are defined in §2 and 3; the last equality involves an integration by parts).

However, the first term in equation (13) must also be included to compute the heating rate of the gas. In fact, the heating rate of the gas is not just determined by the recoil, but also by the fact that the background intensity is greater on the red wing of the Ly$`\alpha `$ line than on the blue wing, and scatterings have an average tendency to bring photons back into the line center. The energy change of the photons due to the first term in equation (13) can only go into the gas (whereas that of the third term is due to redshift and goes into the expansion of the universe, and that of the fourth term is provided by a source). Thus, the correct heating rate is

$$\mathrm{\Gamma }=4\pi h\mathrm{\Delta }\nu _Dn_H\sigma _0𝑑xx\frac{}{x}\left(\frac{\varphi }{2}\frac{J}{x}+\eta \varphi J\right).$$
(4)
Under steady-state conditions (and no source term; here, we consider continuum photons only, to simplify the discussion), we can replace the quantity inside parenthesis by $`\tau _{GP}^1J/x`$, and integrating by parts, we obtain

$$\mathrm{\Gamma }=(Hn_Hk_B)\sqrt{T_0T}\frac{J_{c0}}{\stackrel{~}{J}_0}I_c,$$
(5)
which is the expression in equation (17). Note that, because $`\stackrel{~}{J}_0=(n_Hc)/(4\pi \nu _\alpha )`$, the heating rate per unit volume is independent of the atomic density, except for the weak dependence of $`I_c`$ on $`\tau _{GP}`$, whereas MMR assumed that the heating rate per atom is independent of $`n_H`$.


text:
# Experimental Lagrangian Acceleration Probability Density Function Measurement

## 1 Introduction

Recent developments of new experimental techniques have allowed significant progress in the Lagrangian description of turbulence . Although this aspect has been investigated theoretically since the beginning of the 20th century, the extreme difficulty of recording particle trajectories in high Reynolds number turbulence restricted experimental investigations. The availability of very high speed imaging sensors gives access to the Lagrangian acceleration in fully developed turbulence. It has been shown that the probability density function (PDF) of a Lagrangian acceleration component is strongly non-Gaussian, i.e., the PDF exhibits extremely large tails . This has been also observed in direct numerical simulations of fully developed turbulence

These observations resemble features obtained in the field of non-extensive statistical mechanics . This opens two promising perspectives: first toward a statistical theory of turbulence and second toward a dynamical modeling of the motion of fluid particles. The latter comes from the theoretical link first established by C. Beck  between the Tsallis entropy and a stochastic equation for the Lagrangian acceleration. These developments led to a series of models for the PDF of Lagrangian accelerations. The multitude of models needs to be tested by an accurate comparison with experimental data, especially in the low probability tails of the acceleration PDF.

Here we present new experimental results that have allowed us to improve the estimation of the acceleration PDF. First we describe the modifications that have been implemented on the previous experiment, then we show the acceleration PDF and compare it to the different classes of models. Finally we show and discuss some PDF’s of acceleration conditioned on velocity which are improvements over results presented previously in Sawford et al. .

## 2 The experiment

The experimental set-up has been described in detail in previous publications . The reader should refer to these for detailed information on the experiment. Let us recall briefly the main features. The flow under consideration is of the Von Kármán kind: the water is driven by two counter-rotating disks, 20 cm in diameter, 33 cm apart, mounted in a cylindrical container. The rotation frequency of the disks is regulated at 3.5 Hz for the experiment under consideration here. The integral scale of the flow has been shown to be 7 cm and the energy dissipation rate to be 1.14 m2s-3 . The transverse and axial spatial coordinates of small tracer particles (diameter 46 $`\mu `$m, density $`\rho =1.06g/cm^3`$) are recorded by two silicon strip detectors. Each detector provides a 1D recording of the particle trajectories at a sampling rate up to 70 kHz. The root-mean-square velocity is 0.43 ms-1, so that the Taylor based Reynolds number is 690. The Kolmogorov time is $`\tau _\eta =0.93`$ ms, the sampling frequency is $`F_s=69.9`$ kHz so that $`\tau _\eta `$ corresponds to 65 samples. We study here only one transverse component of the acceleration. Even though the velocity components have been shown to be anisotropic, the acceleration components are close to isotropic at that Reynolds number .

Compared to the experiment reported earlier, we are now able to record data three time faster (in terms of global efficiency, i.e., including transfer rate from the digitizer to the computer, repetition rate, etc.) because of a hardware upgrade of the data-acquisition computers. This allowed us to obtain more data in a reasonable time. The laser used for illumination has also been upgraded to a frequency doubled pulsed YAG laser with mean output power as high as 35 watts, a typical pulse duration of 200 ns and a pulse repetition rate of up to 70 kHz. This pulsed laser improves the efficiency of the silicon strip detectors due to the short duration of the pulses.

Another difference lies in the signal processing. As stressed in Voth et al. , one has to low-pass filter the trajectories to filter the noise in order to have access to the acceleration. In our previous work , the filtering and differentiation were achieved by fitting a parabola to the data, in a time window whose duration had to be appropriately chosen. Here, this procedure has been modified to improve the performance of the data processing algorithm. We now use a differentiating and filtering kernel $`k(\tau )`$ so that the acceleration is obtained by convolution $`a=kx`$, where $``$ is the convolution product and $`x`$ is the spatial coordinate of the particle. We have chosen a Gaussian kernel

$$g(\tau )=\frac{1}{\sqrt{\pi w^2}}\mathrm{exp}\left(\frac{\tau ^2}{w^2}\right).$$
(1)
To achieve both a low pass filter and the second order differentiation of trajectories, one simply has to differentiate twice the Gaussian kernel

$$k(\tau )=\frac{2}{\sqrt{\pi }w^3}\left(\frac{2\tau ^2}{w^2}1\right)\mathrm{exp}\left(\frac{\tau ^2}{w^2}\right).$$
(2)
In practice, however, time is discrete and one has to use a finite time window, therefore the kernel is truncated. This requires a modification of the normalization. The following kernel was used in the data analysis (for discrete times):

$`k(\tau )`$ $`=`$ $`A\left({\displaystyle \frac{2\tau ^2}{w^2}}1\right)\mathrm{exp}\left({\displaystyle \frac{\tau ^2}{w^2}}\right)+B\text{ if }|\tau |L`$ (3)
$`k(\tau )`$ $`=`$ $`0\text{ otherwise}`$ (4)
where $`A`$ and $`B`$ are normalization constants chosen to ensure that $`k1=0`$ and $`k\tau ^2=2`$, i.e., the second derivative of a constant is 0 and the second derivative of a parabola is 2. Note that, if $`L`$ equals one time sample $`1/F_s`$, one recovers the simple 3 points discrete second derivative (independently of the value of $`w`$). If $`Lw1/F_s`$ one converges to the continuous kernel $`g`$ when $`F_s`$ goes to infinity. The results obtained from the filtering operation are found to be in agreement with those obtained by the parabolic fitting procedure. Note that due to the truncation, the actual width of the filter is a slightly lower that $`w`$.

## 3 Probability density function of Lagrangian acceleration

### 3.1 Experimental estimation of the PDF

The results presented here required continuous acquisition of data for two weeks. The data set consists of $`\mathrm{1.7\hspace{0.33em}10}^8`$ points, which include the new measurements and the data reported in the previous publications . The experimental probability density function (PDF) is presented in Fig. 1 as dots. The new data allows the estimation of the PDF to reach probabilities more than one order of magnitude lower than the results in , i.e., less than $`10^7`$. Fig. 1 also shows the stretched exponential fit:

$$P(a)=C\mathrm{exp}\left(a^2/\left(1+\left|\frac{a\beta }{\sigma }\right|^\gamma \right)\sigma ^2\right)$$
(5)
with $`\beta =0.513\pm 0.003`$, $`\gamma =1.600\pm 0.003`$, $`\sigma =0.563\pm 0.02`$ and $`C`$ is a normalization constant ($`C=0.733`$). This expression is the same as the one proposed earlier , with slightly different values due to improved accuracy of the fit and perhaps also a difference in the Reynolds number. The initial values were $`\beta =0.539`$, $`\gamma =1.588`$, $`\sigma =0.508`$ at $`R_\lambda =970`$. This stretched exponential fit is shown to follow the experimental data very well.

In Fig. 2, an estimation of the statistical dispersion of the PDF is shown. It has been obtained by splitting the whole data set into 6 subsets. The relative rms uncertainty is displayed. One observes a minimum of 3% for small values of the acceleration, $`|a|<10a^2^{1/2}`$, for which the number of samples is maximum. The uncertainty then increases as $`|a|`$ increases. It remains lower than 40% for $`|a|<40a^2^{1/2}`$.

One consequence of our better estimation of the acceleration PDF is that we can conclude that the fourth moment converges. Fig. 3 displays the contribution to the fourth moment: $`a^4P(a)`$. The first observation is that this function decreases as $`|a|`$ increases. From our fit, the decay is consistent with a stretched exponential:

$$P(a)\mathrm{exp}(|a|^{0.4}).$$
(6)
This means that the integral $`a^4P(a)`$ converges. The data gives a flatness of

$$F=\frac{a^4}{a^2^2}=55\pm 8.$$
(7)
and the stretched exponential fit gives a flatness value equal to $`F=55.2`$. The possible influence of the particle response time on this value is expected to be low as the characteristic time of the particle is $`0.12\tau _\eta `$. The filtering procedure has a limited impact because of the low value of the filter width. Changing the width from $`0.23\tau _\eta `$ to $`0.31\tau _\eta `$ decreases the variance by 10$`\%`$ and the flatness by 15$`\%`$. The shape of the PDF is unchanged for $`a/a^2^{1/2}<25`$. Thus we expect the flatness to be slightly underestimated but do not expect any qualitative change in the shape of the PDF.

### 3.2 Comparison with models inspired by non-extensive statistical mechanics

Among the large class of superstatistics , four specific examples have been proposed to model the observed acceleration PDF. The first one, proposed by C. Beck  (and referred to as Beck-$`\chi ^2`$ here after), is based on the $`\chi ^2`$-distribution. One very interesting feature of this analysis is the possibility of making a connection between a statistical approach and a dynamical approach. It allows us to build a stochastic equation for the Lagrangian time dynamics of acceleration. In the same spirit, A. M. Reynolds  and then C. Beck  proposed a second model (referred as Reynolds and Beck-Log resp.) based on log-normal statistics. The approach of Reynolds is much more ambitious. He tries to build a full stochastic equation that models the acceleration. This model takes into account not only the small time scales dynamics but also the large scales. He takes all numerical parameters from experiments or direct numerical simulations reported in the literature. Therefore there are no free parameters in his model, once the Reynolds number is chosen. Thus the result of the model is not a fit to the experimental data.

Another more classical approach has been proposed by T. Arimitsu and N. Arimitsu  (referred as Arimitsu). They use the properties of scale invariance of the Navier-Stokes equation to build a multifractal description of the acceleration PDF.

These four models are displayed in figures 4 and 5 together with the experimental data. The first obvious observation is that the $`\chi ^2`$-distribution model (with $`q=3/2`$) does not fit the data. The decay of the tails displays qualitatively incorrect behavior. Indeed, for any value of the parameter $`q`$, the model predicts a power-law decay which is not observed. For $`q`$ close to 3/2, which gives the best fit, the fourth moment diverges which is not the case for the experimental data.

The modification of the distribution to log-normal, introduced later by C. Beck, makes the fit qualitatively much better. In particular the flatness, for $`s^2=3.0`$ takes a finite value $`F=60.3`$. The Reynolds model does not agree as well with the data as the others. One reason is that there are no free parameters in this model. Second, we low-pass filter the acceleration at the scale $`0.23\tau _\eta `$. This scale is small but not zero. If the temporal position signals from the Reynolds model are low-pass filtered at the same scale then the overlap with the experimental data is highly improved . The multifractal model by Arimitsu also provides a good agreement with the experimental results. Concerning the region of low acceleration region, the Arimitsu model captures the correct behavior, contrary to the other models. Reynolds’ and Beck’s log-normal models overestimate the zero acceleration probability whereas the Tsallis statistics model underestimates it.

For all these models the departure from the experimental curve for large accelerations is best viewed in Fig. 5. Most of the proposed models do not peak at the same value of acceleration as the experimental curve. Nevertheless most of them are very close to the experiment.

### 3.3 Conditional probabilities

Conditional acceleration statistics are important in many second order stochastic models . A very large number of statistics is required to obtain acceleration PDF’s conditional on large velocities since large velocities are rare events. The new large data set allowed us to calculate distributions conditional on larger velocities than those shown in Sawford et. al. . Figure 6 shows PDF’s of the transverse component of acceleration conditional on the transverse component of velocity. The conditional probability distributions are highly non-Gaussian. They have almost the same shape as the unconditional distribution. The conditional acceleration variance increases with increasing velocity which seems to be due mainly to a widening of the tip of the distribution as can be see in in Figure 6. For homogeneous isotropic turbulence the conditional mean acceleration should be zero. Departures from zero reflect the anisotropy of our flow although DNS of homogeneous isotropic turbulence has also shown slight departures from zero . Sawford et. al. present a scaling argument that indicates that the conditional acceleration component variance $`A_i^2|u_iu_i^6`$ to leading order. A fit to our data in Fig. 6c. shows good agreement with this. The dependence of acceleration distributions on the velocity violates local homogeneity  which is an assumption of Kolmogorov 41 theory.

## 4 Conclusion

Our new set of data helps to discriminate among different classes of models. In particular, the model based on the Tsallis entropy is shown to disagree with the experimental observations: the acceleration PDF does not display any power law decrease and is well fitted by a stretched exponential decay. One consequence is that the fourth moment of the acceleration component is finite and we observe a value of the flatness close to 55. The models based on log-normal statistics proposed by Beck and Reynolds or multifractal analysis (Arimitsu) are indeed close to our experimental observations. Since the only difference between Beck’s two models is the underlying statistics (Tsallis or log-normal), the hypothesis of Tsallis statistics must be discarded in order to reproduce the observed behavior of the acceleration PDF. The use of a log-normal statistics and the physical link of the $`\beta `$ parameter in Beck’s models to the dissipation rate are in agreement with the refined 1962 Kolmogorov-Obukhov theory of turbulence which is known to describe well the features of intermittency.

We also show that the joint PDF of acceleration and velocity is not gaussian. The acceleration variance conditioned on velocity is also not constant. These results have important consequences on the structure of the stochastic equation required to model dispersion by turbulent flows.

This work has been supported by NSF grant #9988755. We thank Arthur La Porta for valuable discussions and writing substantial parts of the data analysis code.