text:
# A Radio Telescope Search for Axions

## 1 Introduction

The axion is a well-motivated hypothetical elementary particle. It is the pseudo-Goldstone boson associated with a Peccei-Quinn symmetry invoked to solve the strong-CP problem in QCD pq  weinberg . The axion is a good cold dark matter candidate k&t , it could today dominate potential wells of most astrophysical objects. The axion spontaneously decays into two photons with lifetime k&t

$$\tau _a=6.8\times 10^{24}\left(\frac{0.72}{E/N1.95}\right)^2\left(m_a/\mathrm{eV}\right)^5\mathrm{s},$$
(1)
where $`E/N`$ is a model-dependent parameter. The search described here seeks a spectral line from axion decays in halos of Local Group dwarf galaxies.

There is a model range of perhaps a factor of 10 for the axion-to-two-photon coupling constant $`g_{a\gamma \gamma }`$ kim:1 . In two very different benchmark models, $`g_{a\gamma \gamma }`$ varies from about $`1.35\times 10^{16}(m_a/\mu \mathrm{eV})\mathrm{GeV}^1`$ (the DFSZ model zhit ) to $`3.65\times 10^{16}(m_a/\mu \mathrm{eV})\mathrm{GeV}^1`$ (the KSVZ model kim:2 ). Various terrestrial experiments and astrophysical/cosmological arguments constrain the axion mass to the approximate window $`10^6`$ to $`10^3`$ eV, and $`g_{a\gamma \gamma }`$ to less than $`1\times 10^{10}`$ $`\mathrm{GeV}^1`$ turner . In one such astrophysical argument, leading to the “red giant limit”raffelt , axions alter energy transport in red giants. An axion with $`g_{a\gamma \gamma }`$ greater than $`1\times 10^{10}`$ $`\mathrm{GeV}^1`$ would noticeably alter the abundances in the Hertzsprung-Russell horizontal branch. This limit, while the most sensitive in the upper two decades of the allowed axion mass window, does not reach the benchmark model sensitivities. Further, this limit depends on detailed models of red giant stellar evolution.

By contrast, the main assumption for this radio telescope search is that axions dominate potential wells. This search examines Local Group dwarf galaxies. If axions are present in these galaxies, then photons from their decay will appear as narrow lines in radio telescope power spectra bershady .

## 2 Expected Signal

An estimate of the expected axion signal requires three main input parameters: the spatial distribution of dark matter in a potential well, the particular axion model, and the axion mass. Each halo axion can decay into two photons, each with frequency $`\nu =(1/2)m_ac^2/h`$ in the axion rest frame. The volume power density from axion decay in the halo is $`P_V=(1/2)\rho _ac^2/\tau _a`$ where $`\rho _a`$ is the axion mass density and $`\tau _a`$ is the axion lifetime. Projecting $`P_V`$ onto the sky gives the power per solid angle $`P_S(\mathrm{\Omega })`$.

The observed power from axion decays is then given by convoluting $`P_S(\mathrm{\Omega })`$ with the beam pattern $`B(\mathrm{\Omega })`$ integrated over the surface $`A`$ of the antenna

$$P_{obs}=\frac{1}{4\pi D^2}P_S(\mathrm{\Omega })B(\mathrm{\Omega })𝑑A𝑑\mathrm{\Omega },$$
(2)
where D is the distance from Earth to the galaxy. For the Haystack antenna, $`B(\mathrm{\Omega })`$ is a two-dimensional Gaussian with beamwidth $`1.2\times 7\mathrm{m}\mathrm{m}/36.6\mathrm{m}`$. The factors in equation (2) lead to the relation $`P_{obs}m_a^\alpha `$ with $`3<\alpha <5`$ alpha . The value of the scaling parameter $`\alpha `$ depends on the source size and beamwidth. An axion spectral line would appear with a Doppler-broadened width $`\mathrm{\Delta }\nu /\nu \sigma _{LOS}/c`$ where $`\sigma _{LOS}`$ is the line-of-sight velocity dispersion within the object.

We modeled the halo mass distribution as an isothermal sphere chandra  with radial mass density

$$\rho (r)=\frac{\sigma _{LOS}^2}{4\pi G(r^2+a^2)},$$
(3)
where $`G`$ is the gravitational constant, $`r`$ is the radial distance from the object center, and $`a`$ is the core radius binney&tremaine . We use a cutoff radius of $`10a`$ to avoid divergent total mass; these results are only very weakly dependent on the choice of cutoff radius. There are other models with a steeper mass falloff king , the isothermal sphere model leads to the more conservative limit on coupling constant.

We chose three Local Group dwarf galaxies for study: Pegasus, Leo I, and LGS 3 mateo . Dwarf galaxies were chosen because they have low virial velocities, resulting in a narrow signal line, $`\delta \nu /\nu 10^5`$, easily distinguishable from broader instrument-induced structure. Also, these objects have high mass-to-light ratios, and are likely dominated by dark matter. Further, these three galaxies are not known to be tidally disrupted burkert , and the angular diameter of these galaxies — several acrminutes — is well matched to the antenna beamwidth.

## 3 Observations

Observations were made at Haystack Observatory (Westford, MA) on the 37-meter radio telescope haystack  at frequencies from 35.92 to 44.08 GHz.

Each observation is an “off/on” background-subtracted power spectrum with a 160 MHz bandwidth. Each off/on spectrum was derived from two measurements. First, the telescope was pointed a small angle ahead of the object as it moved across the sky and a 30 second “off source” background spectrum recorded. This off source region of the sky has no known radio sources NED , and represents the same section of the atmosphere and radome which will shortly be tracked with the target source. Second, the telescope was pointed at the object and a 30 second “on source” spectrum recorded. Finally, the off source spectrum subtracted from the on source spectrum yielded the off/on spectrum. This off/on technique, in tracing the same arc of the sky and radome twice, minimizes terrestrial artifacts.

Before making an observation, the antenna temperature was calibrated with a blackbody load placed in front of the receiver feed horn. This calibration was repeated for every change in frequency, change in observing object, and otherwise approximately every 10 minutes. This serves to calibrate the instrumental response as a function of frequency. Typical system temperatures were in the range 100–230 K.

There are two major sources of error. The first is “pointing error”, the angular mismatch between commanded and actual antenna direction. From known sources, such as planets and quasars, the estimated pointing error varied from negligible to 15 arcsec. For this analysis, we used a conservative fixed pointing error of 22 arcsec. A fully realized pointing error results in a signal power degradation of 35%, our dominant uncertainty.

The second source of error is uncertainty in the “efficiency” of the telescope, which represents all sources of signal attenuation. These include the atmosphere, scattering from the radome, and varying gravational deformation of the antenna as the dish moves in elevation. Again from known sources, the estimated efficiency is between 29% and 34%. For this analysis, we conservatively use a constant efficiency of 27%.

## 4 Results

Power spectra were collected at the Haystack Observatory on 21 April 99, 15 October 99, 29-31 December 99, and 21 March 00. The total number of spectra at any one frequency is at least three and sometimes more than ten. Higher numbers of observations were made for frequencies that contained unusual features or candidates. (As explained below, these have all been ruled out as axions.)

Each spectrum has a 160 MHz bandwidth and was derived from a 256-time-lag autocorrelation spectrometer. Thus, each bin has Nyquist width 625 kHz; the expected axion signal has a width of about two bins. The high and low frequency ends of some spectra show bandpass filter skirts, so the first and last quarter of the bins are removed from each spectrum. The remaining 80 MHz of data per spectrum are used for all subsequent analysis.

The baseline of each spectrum is estimated from a least-squares fit to a $`4^{\mathrm{th}}`$-order polynomial. This estimated baseline is then subtracted from the spectrum. In Monte Carlo simulations, the order of the fit polynomial improves the axion signal-to-noise ratio up to $`4^{\mathrm{th}}`$-order, after which the sensitivity changes little.

A typical background-subtracted power spectrum is shown in figure 1. The vertical axis shows power spectral density, the horizontal axis shows frequency. The ends of the spectrum have been removed, leaving the central 128 of the original 256 bins. The curve shows the fitted baseline. In addition, a simulated axion of coupling strength $`g_{a\gamma \gamma }=1.2\times 10^9`$ $`\mathrm{GeV}^1`$ has been superimposed near the middle of the spectrum.

Bins from different spectra at the same frequency are combined into a weighted average combining , forming a “combined spectrum”. The combined spectrum for dwarf galaxy LGS 3 is shown in figure 2. The spectrum shows the relatively large system noise temperature near 39 GHz.

Next, the combined spectrum was normalized using its noise temperature combining . Figure 3 is a histogram of the deviations in the LGS 3 normalized combined spectrum. The line shows the expectation for unit variance Gaussian noise.

We formed a “candidate spectrum” by averaging every two adjacent bins of the normalized combined spectrum. Candidates are defined as those bins in the candidate spectrum exceeding a selected threshold. The filled circles in figure 4 show the number of candidates (left vertical axis) versus threshold for LGS 3. The horizontal axis is threshold in units of the normalized combined spectrum’s standard deviation. The curve shows the expected number of candidates versus threshold for unit variance Gaussian noise.

From figure 4 we selected a threshold of 2.3, resulting in 16 candidates. In Monte Carlo simulations we injected axion peaks with varying power and the expected virial width into the unfitted spectra. We determined that an axion of power of about $`10^{18}`$ W gave a 96% search confidence at the selected threshold; this is the search power sensitivity, $`P_{sens}`$. The injected power varied slightly across the spectrum to maintain a flat 96% search confidence. The open circles in figure 4 show the search confidence for axions injected with $`P_{sens}`$. Like the red giant limit, this power sensitivity is much higher than that required to detect KSVZ axions.

An axion signal should always be present at a given frequency. We returned to Haystack Observatory on 21 Mar 00 to examine the 16 candidates from LGS 3. We determined that in this later data the corresponding power for each candidate was small, this is consistent with the 16 candidates being noise fluctuations. We therefore eliminated all 16 candidates while maintaining overall search confidence at 96%.

$`P_{obs}`$ for all three galaxies are approximetely the same. However, the smaller angular size of LGS 3 yields better coupling constant sensitivity alpha . The search procedure for Leo I and Pegasus is identical to LGS 3, except the power threshold is set high enough to exclude any candidates.

## 5 Coupling constant sensitivity

The relation between measured and predicted (KSVZ) power and coupling constant sensitivity is

$$\frac{P_{sens}}{P_{KSVZ}}=\left(\frac{g_{sens}}{g_{KSVZ}}\right)^\alpha ,$$
(4)
where $`P_{sens}\text{ and }g_{sens}`$ are the axion power and coupling constant sensitivity of this search. $`\alpha `$ is that of section 2 and $`P_{KSVZ}`$ is obtained from equation 2. Figure 5 shows $`\alpha `$ as a function of frequency for LGS 3, Leo I and Pegasus. As frequency increases, the beamwidth of the antenna narrows, the source looks relatively more diffuse, and $`\alpha `$ decreases.

The excluded values of the axion-to-two-photon coupling constant at 96% search confidence for LGS 3, Leo I, and Pegasus are shown in figure 6. The lower horizontal axis is the photon frequency, and the upper horizontal axis is the corresponding axion mass. The vertical axis shows the excluded coupling constant $`g_{a\gamma \gamma }`$. The sensitivity varies with the receiver noise and source. We excluded axion-to-two-photon coupling constants of $`1.0\times 10^9\mathrm{GeV}^1`$ or greater.

## 6 Conclusions

This search ruled out axions of mass 298 to 363 $`\mu \mathrm{eV}`$ with axion-to-two-photon coupling of $`g_{a\gamma \gamma }>1.0\times 10^9`$ $`\mathrm{GeV}^1`$ at 96% confidence. For the first time, a radio telescope was used to search for cold dark matter axions in the allowed mass window.

The sensitivity of this result is comparable to the red giant limit. However, the assumptions are very different; the radio telescope technique does not rely on detailed models of stellar evolution.

The RF-cavity axion search technique sikivie  achieves remarkable sensitivity at the lower axion masses. Radio telescopes are more suited to search for axion of higher mass. Present radio telescope technology allows searches to the very upper end of the allowed axion mass window at sensitivities comparable to this result.

Acknowledgements
We are grateful to the research staff of Haystack Observatory for the use of the facility and helpful discussions. In particular, we thank J. Ball, P. Pratap and P. A. Shute for their guidance. We also thank J. Hewitt for her many helpful suggestions. Support for the use of the Haystack telescope was provided by the National Science Foundation under grant DUE-9952246.


format: text

text:
# Statistical treatment of detection cross-section uncertainties in the analysis of solar neutrino data

## Abstract

We propose a modification to the standard statistical treatment of the detection cross-section uncertainties in the analysis of solar neutrino data. We argue that the uncertainties of the energy-averaged cross sections of the different neutrino fluxes in the same experiment should be treated as correlated. We show that the resulting allowed regions for the neutrino oscillation parameters are significantly larger than the ones obtained with uncorrelated uncertainties.

preprint:     DFTT 23/00  arXiv:hep-ph/0006026

The discrepancy between solar neutrino data  and the predictions of the Standard Solar Model (SSM) , known as the “solar neutrino problem”, represents one of the evidences in favor of neutrino oscillations (see ). The statistical analysis of solar neutrino data allow to obtain information on the values of the neutrino oscillation parameters. The latest solar neutrino data  have been analyzed in the framework of two-neutrino mixing , three-neutrino mixing  and four-neutrino mixing .

One of the crucial ingredients in the analysis of solar neutrino data is the calculation of the statistical covariance matrix $`V`$ of the uncertainties. This matrix determines the $`\chi ^2`$ of the fit through the formula

$$\chi ^2=\underset{j_1,j_2}{}\left(R_{j_1}^{(\mathrm{thr})}R_{j_1}^{(\mathrm{exp})}\right)(V^1)_{j_1j_2}\left(R_{j_2}^{(\mathrm{thr})}R_{j_2}^{(\mathrm{exp})}\right).$$
(1)
Here $`R_j^{(\mathrm{exp})}`$ is the event rate measured in the $`j^{\mathrm{th}}`$ experiment and $`R_j^{(\mathrm{thr})}`$ is the corresponding theoretical event rate, that depends on the neutrino fluxes predicted by the SSM and on the properties of neutrinos that modify the flux along the neutrino propagation. In the case of two-neutrino oscillations the theoretical event rate depends on the mass-squared difference $`\mathrm{\Delta }m^2m_2^2m_1^2`$ and on the mixing angle $`\vartheta `$ (see ).

The covariance matrix $`V`$ includes the uncertainties of the experimental rates $`R_j^{(\mathrm{exp})}`$, and the uncertainties of the theoretical rates $`R_j^{(\mathrm{thr})}`$, quantified by the experimental and theoretical covariance matrices $`V^{(\mathrm{exp})}`$ and $`V^{(\mathrm{thr})}`$, respectively. Since experimental and theoretical errors are independent, the corresponding uncertainties can be added in quadrature:

$$V=V^{(\mathrm{exp})}+V^{(\mathrm{thr})}.$$
(2)
Here we are concerned with the covariance matrix of theoretical uncertainties, $`V^{(\mathrm{thr})}`$, that can be written as the sum of two independent contributions,

$$V^{(\mathrm{thr})}=V^{(\mathrm{cs})}+V^{(\mathrm{fx})},$$
(3)
where $`V^{(\mathrm{cs})}`$ is the covariance matrix of the errors due to the uncertainties of the detection cross sections and $`V^{(\mathrm{fx})}`$ is the covariance matrix of the errors due to the uncertainties of the neutrino fluxes. Very convenient expressions for these matrices have been presented in Refs. . These expressions have been used by many authors for the statistical analysis of solar neutrino data in terms of neutrino oscillations .

Here we would like to examine closely the expression given in  for the covariance matrix of the detection cross section:

$$V_{j_1j_2}^{(\mathrm{cs})}=\delta _{j_1j_2}\underset{i}{}\left(R_{ij_1}^{(\mathrm{thr})}\mathrm{\Delta }\mathrm{ln}C_{ij_1}\right)^2,$$
(4)
where $`R_{ij}`$ is the event rate in the detector $`j`$ due to the neutrino flux produced in the $`i^{\mathrm{th}}`$ thermonuclear reaction in the sun ($`i=pp,pep,\mathrm{He}p,\mathrm{Be},\mathrm{B},\mathrm{N},\mathrm{O},\mathrm{F}`$) and $`C_{ij}`$ is the corresponding energy-averaged cross section. The quantity $`\mathrm{\Delta }\mathrm{ln}C_{ij}=\mathrm{\Delta }C_{ij}/C_{ij}`$ is the relative uncertainty of the cross section $`C_{ij}`$, whose value is given in Ref. .

The expression (4) is based on two assumptions: 1) the errors of the detection cross sections in different experiment are independent; 2) the errors of energy-averaged cross section of the different neutrino fluxes in the same experiment are uncorrelated. Indeed, in Eq. (4) the errors of energy-averaged cross section of the different neutrino fluxes in the same experiment, $`R_{ij}^{(\mathrm{thr})}\mathrm{\Delta }\mathrm{ln}C_{ij}`$, are added in quadrature to the corresponding diagonal element $`V_{jj}^{(\mathrm{cs})}`$.

We think that the first assumption is appropriate, but the second one is questionable. This is due to the fact that the correlations between the uncertainties of each detection cross section at different energies is not known. Moreover, the uncertainties of the energy-averaged cross section of the neutrino fluxes that have an overlap of energy ranges are certainly not uncorrelated.

Since the correlations of the errors of energy-averaged cross section of the different neutrino fluxes in the same experiment are not known, the correct attitude consists in adopting the most conservative approach, assuming a complete correlation. In this case the errors must be added linearly:

$$V_{j_1j_2}^{(\mathrm{cs})}=\delta _{j_1j_2}\left(\underset{i}{}R_{ij_1}^{(\mathrm{thr})}\mathrm{\Delta }\mathrm{ln}C_{ij_1}\right)^2.$$
(5)
This approach has been also recommended in Ref. , that is the standard reference for the detection cross section in Gallium experiments. Let us notice that it reflects a rather realistic possibility: that in which the errors of the cross section in a given experiment have the same sign at all energies.

We have performed a fit of the total rates measured in solar neutrino experiments in terms of oscillations between two active neutrinos ($`\nu _e\nu _\mu `$ or $`\nu _e\nu _\tau `$) in order to study the change of the allowed regions for the neutrino oscillation parameters $`\mathrm{\Delta }m^2`$ and $`\vartheta `$ when the expression (5) for $`V^{(\mathrm{cs})}`$ is used in place of Eq. (4). The results of our fits are presented in Figs. 15. We have used the latest measurements of the event rates in the Homestake  and Super-Kamiokande  experiments, and the weighted average of the rates measured in the two Gallium experiments GALLEX  and SAGE . The values of these rates are given in Table I of Ref. .

Our calculation of the theoretical event rates $`R_{ij}^{(\mathrm{thr})}`$ follows the standard method described in several papers for matter-enhanced MSW transitions  and vacuum oscillations . We calculate the MSW survival probability of $`\nu _e`$’s in the Sun using the standard analytic prescription  and the level-crossing probability appropriate for an exponential density profile . We calculate the regeneration in the Earth using a two-step model of the Earth density profile, that is known to produce results that do not differ appreciably from those obtained with the correct density profile. We have used the tables of neutrino fluxes, solar density and radiochemical detector cross sections available in Bahcall’s web page . We have neglected the matter effects that slightly affect the vacuum oscillation solutions of the solar neutrino problem, as discussed in . In Figs. 15 we have used as abscissa the parameter $`\mathrm{tan}^2\vartheta `$ that allows a better view of the regions at large mixing angles with respect to the usual parameter $`\mathrm{sin}^22\vartheta `$. The parameter $`\mathrm{tan}^2\vartheta `$ has been employed in the past in the framework of three-neutrino mixing  and its use in the framework of two-neutrino mixing has been recently advocated in  because it allows to explore the possible presence of allowed regions with $`\vartheta >\pi /4`$.

The allowed regions in Figs. 15 are calculated through the constraint $`\chi ^2\chi _{\mathrm{min}}^2+\mathrm{\Delta }\chi ^2(\alpha )`$, where $`\chi _{\mathrm{min}}^2`$ is the global minimum of the $`\chi ^2`$, $`\chi _{\mathrm{min}}^2=0.42`$ for $`\mathrm{\Delta }m^2=5.1\times 10^6\mathrm{eV}^2`$ and $`\mathrm{tan}^2\vartheta =1.6\times 10^3`$ in both cases of uncorrelated (Eq. (4)) and correlated (Eq. (5)) cross section uncertainties. We consider confidence levels $`\alpha =0.90,0.95,0.99`$ (the regions inside the solid, short-dashed and long-dashed contours, respectively) for two degrees of freedom, which give $`\mathrm{\Delta }\chi ^2(0.90)=4.61`$, $`\mathrm{\Delta }\chi ^2(0.95)=5.99`$, and $`\mathrm{\Delta }\chi ^2(0.99)=9.21`$.

The linear addition in (5) of the uncertainties of the detection cross section lead to larger values of the diagonal elements of the covariance matrix $`V`$ and to smaller values of the $`\chi ^2`$ in Eq. (1), with respect to the case of uncorrelated uncertainties (Eq. (4)). Since the variation of $`\chi ^2`$ in the space of the neutrino oscillation parameters $`\mathrm{tan}^2\vartheta `$, $`\mathrm{\Delta }m^2`$ near the minimum of the $`\chi ^2`$ is proportional to $`V^1`$ (the contribution of the variation of $`V^1`$ is negligible near the minimum),

$$\delta \chi ^22\underset{j_1,j_2}{}\left(R_{j_1}^{(\mathrm{thr})}R_{j_1}^{(\mathrm{exp})}\right)(V^1)_{j_1j_2}\left(\frac{R_{j_2}^{(\mathrm{thr})}}{(\mathrm{tan}^2\vartheta )}\delta (\mathrm{tan}^2\vartheta )+\frac{R_{j_2}^{(\mathrm{thr})}}{(\mathrm{\Delta }m^2)}\delta (\mathrm{\Delta }m^2)\right),$$
(6)
the slope of the $`\chi ^2`$ is smaller in the case of correlated uncertainties (Eq. (5)) than in the case of uncorrelated uncertainties (Eq. (4)). Therefore, we expect to obtain larger allowed regions in the case of correlated uncertainties, with respect to those obtained with uncorrelated uncertainties.

Figures 1 and 2 show the MSW allowed regions obtained using the expressions (4) and (5) for $`V^{(\mathrm{cs})}`$, respectively. One can see that all the three MSW allowed regions, usually called SMA (for $`\mathrm{\Delta }m^25\times 10^6\mathrm{eV}^2`$ and $`\mathrm{tan}^2\vartheta 10^3`$), LMA (for $`\mathrm{\Delta }m^22\times 10^5\mathrm{eV}^2`$ and $`\mathrm{tan}^2\vartheta 0.3`$) and LOW (for $`\mathrm{\Delta }m^210^7\mathrm{eV}^2`$ and $`\mathrm{tan}^2\vartheta 0.5`$), are larger when the expressions (5) with correlated cross section uncertainties is used. The minimum and maximum values of $`\mathrm{tan}^2\vartheta `$ and $`\mathrm{\Delta }m^2`$ in the 99% CL allowed regions in Fig. 1 (unc) and Fig. 2 (cor) are listed in Table I and the difference between these regions is illustrated in Fig. 3. One can see that the SMA-cor region is slightly larger than the SMA-unc region. The LMA-cor region is appreciably larger than the LMA-unc region, allowing values of $`\mathrm{\Delta }m^2`$ and $`\mathrm{tan}^2\vartheta `$ as high as $`3.5\times 10^4\mathrm{eV}^2`$ and $`0.82`$, in contrast to the upper limits $`2.4\times 10^4\mathrm{eV}^2`$ and $`0.67`$ for the LMA-unc region. The larger effects occur in the LOW region. In the LOW-cor region $`\mathrm{tan}^2\vartheta `$ can be as large as 1, corresponding to maximal mixing, and $`2.4\times 10^8\mathrm{eV}^2\mathrm{\Delta }m^23.5\times 10^7\mathrm{eV}^2`$, whereas in the LOW-unc region $`\mathrm{tan}^2\vartheta 0.82`$ and $`5.6\times 10^8\mathrm{eV}^2\mathrm{\Delta }m^22.6\times 10^7\mathrm{eV}^2`$.

Figures 4 and 5 show the vacuum oscillation allowed regions obtained with uncorrelated (Eq. (4)) and correlated (Eq. (5)) cross section uncertainties, respectively. One can see that the allowed regions obtained with correlated cross section uncertainties are significantly larger that the ones obtained with uncorrelated cross section uncertainties. For example, the 99% CL allowed region at $`\mathrm{\Delta }m^25\times 10^{10}\mathrm{eV}^2`$ is considerably larger with correlated cross section uncertainties (Fig. 5) than with uncorrelated ones (Fig. 4), and the 90% CL allowed region in Fig. 5 at $`\mathrm{\Delta }m^22.7\times 10^{10}\mathrm{eV}^2`$ does not exist in Fig. 4.

In conclusion, we have proposed a modification of the standard statistical treatment of the detection cross-section uncertainties in the analysis of solar neutrino data. We have argued that, because of lack of information, a conservative approach should be used in the treatment of the correlations between the energy-averaged cross sections of the different neutrino fluxes in the same experiment, leading to the assumption of correlated errors and to a linear addition of the uncertainties. We have shown in Figs. 15 that the resulting allowed regions for the neutrino oscillation parameters are significantly larger than those obtained with uncorrelated detection cross sections uncertainties. We think that the appropriate treatment of the detection cross section uncertainties proposed here will be rather important in the future, when more precise solar neutrino data obtained in the Super-Kamiokande, SNO , Borexino , GNO  and other experiments will be published.

###### Acknowledgements.

We would like to thank A. Bottino for useful suggestions and encouragement. C.G. would like to thank P. Krastev, E. Lisi, H. Murayama and S. Petcov for useful discussions on the fit of solar neutrino data.


text:
# Constraints on the physical properties of the damped Ly𝛼 system of Q0000-2619 at 𝑧 = 3.054

## 1 Introduction

Fine structure transitions observed in the absorption spectra of quasars provide unique information on the temperature of the microwave background at the redshift of the absorber, on the intensity of the UV-field and on the density of the absorbing system (Bahcall & Wolf (1968)). The measured excitation temperatures, or upper limits to the excitation temperature, of the fine structure of C i and C ii has been used by several authors to constrain the temperature of the Cosmic Microwave Background radiation (CMB) up to a redshift of 4.38 (Lu et al. (1996)). The most significant published constraints on the CMB temperature at $`z>0`$ are summarized in Fig. 1. The Big Bang cosmological model predicts a simple relationship between the CMB temperature and the redshift $`z`$ (e.g. Peebles (1993)):

$$T_{\mathrm{CMB}}(z)=T_{\mathrm{CMB}}(0)(1+z)$$
(1)
Alternative anisotropic cosmological models (Phillips (1994)) make strong claims for a value of $`T_{\mathrm{CMB}}`$ more than 5 K below the standard prediction at $`z=3`$. So far all the measured excitation temperatures and upper limits are consistent with a Friedmann universe.

The observed excitation temperature (or upper limit) of C i and C ii fine structure has also been used to constrain the density of the absorbing systems (Chaffee et al. (1988); Songaila et al. (1994); Ge et al. (1997)).

Here we discuss the detection of C ii and C ii absorption at $`z=3.0543`$ in the spectrum of Q0000$``$2619 obtained at the NTT. The absorption lines are associated with a damped Ly$`\alpha `$ system with neutral hydrogen column density of $`N`$(H i)$`=1.5\pm 0.5\times 10^{20}\mathrm{cm}^2`$ (Savaglio et al. (1994), SOM hereafter). By measuring the equivalent width of the C ii multiplet absorption lines we derive strict upper limits on the temperature of the CMB at this redshift and constrain the density of the absorbing systems. In the next section the observations are briefly presented. In Sect. 3 the C ii fine structure excitation temperature, which gives directly the upper limit on the CMB temperature at redshift 3.0543, is derived. In Sect. 4 the relative strengths of the different excitation mechanisms are reviewed. The constraints on the density and the UV field in the $`z=3.0543`$ damped Ly$`\alpha `$ system are derived in Sect. 5 and discussed in Sect. 6.

## 2 Observations

In October 1990 echelle observations of Q0000$``$2619 at $`z=4.12`$ were obtained with the ESO Multi Mode Instrument (EMMI) (D’Odorico (1990)) at the ESO NTT telescope. The spectra cover the wavelength range from 4400 Å to 9265 Å with a resolution of 0.2 and 0.3 Å between 4700 and 8450 Å and signal-to-noise ratio S/N$`=1560`$ per resolution element. The data were reduced and analysed by Savaglio et al. (1997) to which we refer the reader for a detailed description of the observations and data reduction procedure. In the spectrum of Q0000$``$2619, SOM have identified nine metal absorption systems; among these two are known damped systems at redshifts 3.054 and 3.390. Eight of the nine systems have redshift greater than 3. We carefully inspected the spectrum looking for absorption from the C i and C ii ground state multiplet. For all systems at redshifts greater than 3, the absorption from the C i ground state multiplet would land redwards of the Lyman-$`\alpha `$ emission. No absorption from C i is detectable in the spectrum.

C ii absorption was detected for both the damped systems at $`z=3.054`$ and $`z=3.390`$ (SOM). For these systems the C ii ($`\lambda 1334`$) absorption line lands bluewards of the Ly$`\alpha `$ emission. The C ii at $`z=3.3913`$ is heavily blended, while the C ii at redshift $`z=3.0543`$ is reasonably clean, despite falling in the Lyman forest at $`\lambda =5410.6`$. At $`\lambda =5415.4`$ we detect a weak absorption line (3.5$`\sigma `$ confidence) consistent with absorption from the excited fine-structure level of C ii at redshift $`z=3.0543`$.

Table 1 summarizes the data on the absorption lines detected for the system at $`z=3.0543`$.

## 3 Excitation temperature

At a redshift of 3.0543, the C ii $`J=1/2`$ ($`\lambda `$1334.53 Å) and $`J=3/2`$ ($`\lambda \lambda `$1335.66,1335.71 Å) absorption lines land respectively at 5410.54 Å and 5415.17,5415.37 Å. Fig. 2 shows the spectrum of Q0000$``$2619 in the vicinity of the C ii multiplet. In this wavelength range the spectral resolution is of 0.2 Å and the S/N per resolution element is about 15.

The multiplet lands in the Lyman forest and in the damping wing of the Ly$`\alpha `$ absorption of the $`z=3.390`$ damped system, at $``$ 5337 Å. The C ii ground state absorption line is slightly blended. We used a multiple Gaussian fit to deblend the C ii line and measured an equivalent width of $`W_\lambda =1.16\pm 0.02`$ Å. The C ii absorption line is detected at $`3.5\sigma `$. For C ii $`W_\lambda =0.077\pm 0.02`$. The Gaussian fit to the line has a FWHM of 0.29 Å which is consistent with the instrument resolution. The two equivalent widths correspond respectively to $`\mathrm{log}(W_\lambda /\lambda )=3.67`$ and $`4.85`$. In both cases a local continuum level corresponding to the damping wing of the Ly$`\alpha `$ absorption at $`z=3.390`$ has been used. The C ii absorption line is well fitted by a Gaussian having a FWHM of $`53\pm 5`$ km s-1, which given the instrumental resolution in this range of 11 km s-1 corresponds to an intrinsic $`b`$ parameter of $`31\pm 3`$ km s-1. In Fig. 3 the theoretical curve of growth of C ii is plotted for three values of $`b`$: 28, 31 and 34 km s-1. The two values of $`\mathrm{log}(W_\lambda /\lambda )`$ are also shown. The corresponding column density values derived for C ii and C ii are summarized in Table 2.

According to the Boltzmann equation, an excitation temperature $`T_{\mathrm{ex}}`$ can be expressed in terms of the column densities $`N_1`$ and $`N_0`$ in the excited and the ground-state level:

$$N_1/N_0=g_1/g_0\mathrm{exp}(\mathrm{\Delta }T_{10}/T_{\mathrm{ex}})$$
(2)
where $`k\mathrm{\Delta }T_{10}`$ is the energy difference between the excited level (1) and the ground level (0). For the fine structure levels $`J=3/2`$ and $`J=1/2`$ of C ii, $`\mathrm{\Delta }T_{10}=91.2`$ K. The weights $`g_J`$ are given by $`2J+1`$. The derived ratio of column densities, $`N_1(J=3/2)/N_0(J=1/2)`$ corresponds to an excitation temperature of 19.6 K (for $`b=28`$ km s-1), 20.5 K (for $`b=31`$ km s-1) and 21.6 K (for $`b=34`$ km s-1). The excitation temperature of 21.6 K provides a strict upper limit on the temperature of the CMB at the absorber redshift of 3.0543; this upper limit would hold even if the C ii absorption were a spurious effect of the Lyman forest. The CMB temperature at this redshift is predicted to be 11.05 K by the Big Bang model.

The excitation temperature provides a strict upper limit to the CMB temperature as other sources may contribute appreciably to the excitation. In fact, if the C ii absorption is not a spurious effect of the line forest, other mechanisms have to be at play to explain an excitation temperature of at least 19.6 K. The possible excitation mechanisms are reviewed in the next section.

## 4 Excitation mechanisms

The higher fine-structure states of a ground state multiplet can be populated by (1) particle collisions, (2) direct excitation by infrared photons and (3) indirect excitation by ultraviolet photons.

The equilibrium between the excitation and de-excitation of the C ii $`J=1/23/2`$ fine structure is described by:

$`N_0({\displaystyle \underset{j}{}}\sigma _{01}vn_j+B_{01}U(\nu _{01})+\mathrm{\Gamma }_{01})=`$ (3)
$`N_1(A_{10}+{\displaystyle \underset{j}{}}\sigma _{10}vn_j+B_{10}U(\nu _{01})+\mathrm{\Gamma }_{10})`$
Here the collision excitation rate is expressed as $`\sigma vn`$, where $`\sigma v`$ is the temperature averaged product between the cross-section and the particle velocity and $`n`$ is the particle density, $`j=`$ H, $`e,p`$. The direct photon excitation rate is expressed as the product between the Einstein probability coefficient for induced transition ($`B`$) and the radiation energy density per frequency interval ($`U(\nu _{01})`$). The UV pumping rate is represented by the coefficient $`\mathrm{\Gamma }`$, which includes the UV energy density term. $`A_{10}`$ is the Einstein probability coefficient for spontaneous transition and for the C ii transition $`J=3/21/2`$, $`A_{10}=2.29\times 10^6`$s-1 (Nussbaumer & Storey (1981)). For convenience we can divide both sides of Eq. (3) by $`A_{10}`$ and rewrite:

$`N_0({\displaystyle \underset{j}{}}q_{01,j}n_j+b_{01}U(\nu _{01})+\gamma _{01})=`$ (4)
$`N_1(1+{\displaystyle \underset{j}{}}q_{10,j}n_j+b_{10}U(\nu _{01})+\gamma _{10})`$
where $`q_{ij,j}`$, $`b_{ij}`$ and $`\gamma _{ij}`$ are respectively the collisional excitation rate, the Einstein probability coefficient for induced transition and the UV pumping rate, all divided by $`A_{10}`$. In order to establish the importance of the various factors in determining the observed population ratio in the fine structure levels of C ii, we need to evaluate the magnitude of each term in this equation.

### 4.1 Collisional excitation and de–excitation

The particles which may be responsible for collisional excitation are essentially electrons, protons and atomic hydrogen, with different contributions dominating in different temperature and ionization regimes.

The expression for the rate of collisional excitation by electrons is given in Bahcall & Wolf (1968) and, by using the detailed computation of the effective collisional strength given in Hayes & Nussbaumer (1984) and Keenan et al. (1986), one obtains typically $`q_{01,e^{}}=0.15`$ $`\mathrm{cm}^3\mathrm{s}^1`$ for $`T_e=100`$ K and $`q_{01,e^{}}=0.051`$ $`\mathrm{cm}^3\mathrm{s}^1`$ for $`T_e=\mathrm{10\hspace{0.17em}000}`$ K.

The excitation rates as a function of electron temperature for neutral hydrogen collisions are also given in Keenan et al. (1986) and for $`T_e=100`$ K $`q_{01,\mathrm{H}}=2.8\times 10^4`$ $`\mathrm{cm}^3\mathrm{s}^1`$, for $`T_e=\mathrm{10\hspace{0.17em}000}`$ K $`q_{01,\mathrm{H}}=1.5\times 10^3`$ $`\mathrm{cm}^3\mathrm{s}^1`$. The electron collision term will dominate the hydrogen collision term for $`n_e0.002n_\mathrm{H}`$ at $`T_e=100`$ K and for $`n_e0.03n_\mathrm{H}`$ at $`T_e=\mathrm{10\hspace{0.17em}000}`$ K. If the plasma is collisionally ionized the electron collision term will be the dominant one for temperature $`T>\mathrm{10\hspace{0.17em}000}`$ K, when the fraction of ionized hydrogen becomes significant (Bahcall & Wolf (1968)). If the absorbing medium is significantly photoionized and $`n_e0.03n_\mathrm{H}`$ the electron collision term will be the dominant term at any kinetic temperature whereas if photoionization is not signficant and $`n_e<0.002n_\mathrm{H}`$ hydrogen collisions will dominate.

The C ii excitation rate for proton collision becomes comparable to the electron contribution only for temperatures $`T_e10^5`$ K (Bahcall & Wolf (1968)), but at these temperatures C ii is completely destroyed by collisional ionization (Sutherland & Dopita (1993)). We will therefore ignore the proton collision contribution.

The collisional de-excitation rate is given by:

$$\sigma _{10,j}v=\frac{1}{2}\sigma _{01,j}v\mathrm{exp}(91.2\mathrm{K}/T_e)$$
(5)
As a first approximation, the collisional de-excitation term on the right hand side of Eq. (4) can be omitted if $`n_e<1`$ cm-3 or $`n_\mathrm{H}<10^3`$ cm-3.

### 4.2 Direct IR photon excitation and de–excitation

The photons responsible for directly populating C ii fine structure excited levels have a frequency of $`\overline{\nu }_{10}=64.0`$ cm-1. Sources of far–infrared photons are the CMB and thermal dust emission. Since the CMB radiation has a Black Body spectrum the direct excitation rate from CMB photons can be expressed as:

$$b_{01}U(\nu _{01})=2\mathrm{exp}(91.2\mathrm{K}/T_{\mathrm{CMB}})$$
(6)
As we have seen, the Big Bang model predicts the CMB temperature to be 11.05 K at this redshift, moreover the upper limit on the CMB temperature at $`z=4.08`$ set by Lu et al. (1996), constrains the CMB temperature empirically. If we assume that the CMB temperature varies monotonically with $`z`$:

$$T_{\mathrm{CMB}}(z)=T_{\mathrm{CMB}}(0)(1+z)^\alpha $$
(7)
then the measure of Lu et al. (1996) gives $`\alpha <1.05`$ between $`z=0`$ and $`z=4.08`$, that is $`T_{\mathrm{CMB}}<11.85`$ K at $`z=3.0543`$. This implies that $`b_{01}U(\nu _{01})`$ must be lower than $`9.1\times 10^4`$ s-1, i.e. only a minor contribution to the derived population ratio in our system, where $`N_1/N_0`$ for C ii is between 0.02 and 0.03 (Table 2). As a first approximation therefore this term can be neglected and the same is true for $`b_{10}U(\nu _{01})=\frac{1}{2}b_{01}U(\nu _{01})`$.

An infrared photon flux with an intensity comparable to the one measured in the Galactic plane (Bennett et al. (1992), Kogut et al. (1996)) would correspond to an excitation rate at least 2 orders of magnitude smaller than the one for a CMB photon flux with $`T_{\mathrm{CMB}}=10`$ K and can be ignored.

### 4.3 Indirect UV photon excitation and de–excitation

The other important type of photon excitation is UV photon pumping. After the absorption of a photon an atom will usually cascade back through a variety of states, sometimes reaching levels that could not be populated by direct radiative upward transition from the ground state. If $`m`$ represents all the quantum numbers for one of the upper levels, reached by photon absorption, the transition rate from level 0 to level 1, is given by (Spitzer (1978)):

$$\mathrm{\Gamma }_{01}=\underset{m}{}B_{0m}U(\nu _{0m})ϵ_{m1}$$
(8)
where $`ϵ_{m1}`$ is the fraction of downward transitions from level $`m`$ that populate level 1, when the atom first reaches the group of lower levels. For transitions within a multiplet the values of $`ϵ_{mj}`$ are tabulated (e.g. Allen (1963)). To evaluate $`\mathrm{\Gamma }_{01}`$ we considered all the direct upward transitions from C ii ground state $`2p^2P^0`$, longwards of 900 Å: $`{}_{}{}^{2}P_{}^{0}`$ $`{}_{}{}^{2}P`$ ($`\lambda \lambda 903,904`$), $`{}_{}{}^{2}P_{}^{0}`$ $`{}_{}{}^{2}S`$ ($`\lambda \lambda 1037,1036`$), $`{}_{}{}^{2}P_{}^{0}`$ $`{}_{}{}^{2}D`$ ($`\lambda \lambda 1335,1334`$). For the UV field we adopted the Milky Way spectral energy distribution (SED) given by Black (1987) with a Milky Way intensity at 912 Å of $`4.7\times 10^{19}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1\mathrm{Hz}^1`$ (Mathis et al. (1983)), obtaining:

$$\mathrm{\Gamma }_{01}=5.3\times 10^{10}\mathrm{s}^1,$$
(9)
that is $`\gamma _{01}=2.3\times 10^4`$ 111Keenan et al. (1986) obtained $`\mathrm{\Gamma }_{01}=2.4\times 10^{10}\mathrm{s}^1`$ by using the UV intensity field and SED given in Gondhalekar et al. (1980). Since the UV pumping de-excitation rate $`\gamma _{10}`$ is of the same order of magnitude as $`\gamma _{01}`$, the UV de-excitation term can also be omitted in Eq. (4), for any likely UV flux intensity and SED.

The final shape of the balance equation is:

$$N_1/N_0=\underset{j}{}q_{01,j}n_j+\gamma _{01}$$
(10)
with the collisional term being dominated by the electrons or hydrogen atoms contribution according to the absorber ionization degree. The UV field can also be expressed in terms of the hydrogen density through the ionization parameter, $`U=\varphi (\mathrm{H})/n_\mathrm{H}c`$, where $`\varphi (\mathrm{H})`$ is the surface flux of hydrogen-ionizing photons. If we assume a UV flux having the Milky Way SED given in Black (1987), $`\gamma _{01}=0.7n_\mathrm{H}U`$ and:

$$T_{\mathrm{ex}}=\frac{91.2\mathrm{K}}{\mathrm{log}[0.5(_jq_{01,j}f_jn_\mathrm{H}+0.7n_\mathrm{H}U)]}$$
(11)
where $`f_j`$ is the fraction of particle $`j`$ with respect to the hydrogen density.

## 5 Results

From the ratio C ii/C iv observed in our data at $`z=3.0543`$ SOM derived $`\mathrm{log}U3.2`$ and assuming solar abundance ratios a consistent fit to the data was obtained with $`\mathrm{log}U=2.8`$ and $`Z0.001Z_{}`$ (SOM). As shown in Table 2, the derived excitation temperature for the levels $`J=3/2`$ and $`J=1/2`$ of C ii for the damped Ly$`\alpha `$ absorber toward Q0000$``$2619 at $`z=3.0543`$ is between 19.6 and 21.6 K. We can use these values and expression 11 to constrain the density and UV flux at the absorber. However, in order to evaluate Eq. (11), we need to make an assumption about the ionization degree of the gas in the absorber. Since the efficiency of the electron and hydrogen collisional is very different, the collisional excitation term depends critically on the gas ionization degree. To illustrate this we will consider two limiting cases for the absorbing cloud: a quasi-neutral gas, for which $`n_e<0.002n_\mathrm{H}`$, and a 10% ionized plasma.

In Fig. 4 the result of Eq. (11) is plotted, in the case of a quasi-neutral gas, where the collisional excitation is due to collisions with atomic hydrogen. The horizontal continuous line gives the upper limit on the excitation temperature. The points where the curves intersect this line correspond to the density values above which the collisional excitation would be such that the excitation temperature would be higher than the measured value. This is a forbidden region. The density of the absorber must thus be lower than 40 cm-3. This is a strict upper limit. Higher values of the ionization parameters would move the intersection point to lower density values. If the C ii absorbing line is a spurious effect of the Lyman forest, then the C ii equivalent width must be lower than the one we measured and the upper limit on the density must be lower; increasing the cloud kinetic temperature also moves the density-limit to lower values. With $`\mathrm{log}U=2.8`$ the upper limit on the absorber density is 20 cm-3. The shaded area in the figure corresponds to the interval of the derived excitation temperature for the absorber toward Q0000$``$2619. The intersection of the long-dashed curve with the lower bound of this area provides a lower limit of 10 $`\mathrm{cm}^3`$ to the density of the system, if the gas is quasi-neutral ($`n_\mathrm{e}<0.002\mathrm{n}_\mathrm{H}`$).

In Fig. 5, C ii fine structure excitation temperature is plotted as a function of the absorber density (like in Fig. 4) this time for the case of a 10% ionized plasma. Electron collisions dominate collisional excitation in this case and because electron collisions are more efficient than atomic collisions, lower gas densities are required to achieve the same level of C ii fine structure excitation. From Fig. 5 one derives that with ionization parameter $`\mathrm{log}U=3.2`$ the range of possible gas densities is between 1 and 5 $`\mathrm{cm}^3`$. We note however that the exact value of the ionization parameter has little effect on the position of the intersect for the curves in Fig. 5. For any $`\mathrm{log}U2.0`$ (thick–dotted line) we obtain that $`n_\mathrm{H}`$ must be greater than 0.7 $`\mathrm{cm}^3`$, with the same 10% ionization and temperature range. An ionization parameter of $`2.0`$ at a density of 0.7 $`\mathrm{cm}^3`$ correspond to a UV flux at 912 Å of 20 times the Galactic UV field (if the same SED is assumed).

These simple considerations allow us to constrain the absorber density to the range of 0.7 $``$ 40 $`\mathrm{cm}^3`$, for a gas kinetic temperature between 100 and 10 000 K and ionization degree between 0 and 10%. An ionization degree between 0 and 100% would imply possible absorber densities in the range 0.1 $``$ 40 $`\mathrm{cm}^3`$. An ionization parameter $`\mathrm{log}U=3.2`$ and the density range 1.0 $``$ 40 $`\mathrm{cm}^3`$ imply a UV field with H-ionizing photon flux ranging from 2 to 80 times the Galactic H-ionizing flux. The upper value is a strict upper limit on the intensity of the UV flux in the absorbing system.

## 6 Discussion

The observed excitation temperature (or upper limit) of C i and C ii fine structure has been used to constrain the density of the absorbing systems by several authors and the published values are very similar to the upper limit derived here. Songaila et al. (1994) observed absorption from the first fine-structure level of C i in a cloud belonging to a damped Ly$`\alpha `$ system at a redshift of 1.776, towards the quasar Q1331$`+`$170. They measured an excitation temperature of $`7.4\pm 0.8`$ K, while the CMB temperature at this redshift is predicted to be 7.58 K. From their measure they derive a 1$`\sigma `$ upper limit for the cloud density of $`n_\mathrm{H}=7`$ cm-3 for a cloud kinetic temperature of 100 K and of $`n_\mathrm{H}=4`$ cm-3 for a cloud kinetic temperature of 1000 K, given that the CMB is at the predicted temperature. Ge et al. (1997) detected absorption from the ground state and the excited state of C i and C ii in the $`z=1.9731`$ damped Ly$`\alpha `$ system of Q0013$``$004. They measure an excitation temperature of $`11.6\pm 1.0`$ K for C i and of $`16.1\pm 1.4`$ K for C ii. They use the C ii excitation temperature to constrain the density of the cloud to be $`n_\mathrm{H}=21.0\pm 9.6`$ cm-3 if the cloud kinetic temperature is of 100 K and $`n_\mathrm{H}=4.5\pm 2.0`$ cm-3 for a kinetic temperature of 1000 K. With these densities and a photo-ionization parameter $`\mathrm{log}U=3.5`$, their estimate of the H-ionizing photon flux ranges between 3.6 and 17 times the Galactic H-ionizing flux. These values are then combined with the measured C i excitation temperature to derive an upper limit of $`7.9\pm 1.0`$ K or $`10.6\pm 1.0`$ K respectively on the temperature of the CMB at this redshift, with the first value as their best guess, based on photo-ionization modeling of the absorbing cloud.

To constrain the density and UV field of the $`z=3.0543`$ absorber of Q0000$``$2619 we have considered a gas kinetic temperature in the range 100$``$10 000 K. There are not many observational constraints on the gas temperature of damped Ly$`\alpha `$ systems. Cloud photo-ionization models can reproduce the observational data with the gas temperature ranging from 15 to 10 000 K (Chaffee et al. (1988)). However observations of 21 cm absorption line from damped Ly$`\alpha `$ systems at $`z3`$ obtained lower limits on the neutral hydrogen spin temperature of the order of 1000 K (Carilli et al. (1996); Kanekar & Chengalur (1997)). This temperature must not be taken as a measure of the gas mean temperature, but as an indication of the presence of warm gas in the absorbing system. The typical spin temperature measured in the Galactic clouds is $`T_s100`$ K (Braun & Walterbos (1992)). Despite the Galactic cloud spin temperature not being directly comparable with the measure of the spin temperature in a damped Ly$`\alpha `$ system, a $`T_s1000`$ K suggests that for a given total neutral hydrogen column density, the damped Ly$`\alpha `$ system contains a larger percentage of warm phase gas ($`T8000`$ K) than is seen in typical Galactic lines of sight (Carilli et al. (1996)).

The constraints on the density of the damped system can be used to estimate its size. In the simplistic hypothesis that the absorber is homogeneous a neutral hydrogen column density of $`1.5\pm 0.5\times 10^{20}\mathrm{cm}^2`$ and a density range of 0.7 $``$ 40 $`\mathrm{cm}^3`$ imply that the size of the system along the line of sight is between 1$``$100 pc. If the filling factor is significantly lower than 1, this estimate is a lower limit. It has been proposed that the DLAs with neutral hydrogen column density $`N`$(H i)$`2\times 10^{20}\mathrm{cm}^2`$, are large and massive galactic disks (e.g., Prochaska & Wolfe (1997)) with typical sizes of a few kpc. The constraints we derived for this damped system at $`z=3.0543`$ indicate that the size of this absorber is of the order of the size of giant hydrogen clouds in our galaxy or of a galactic disk seen face-on. If this is not a galactic disk, at these redshift, systems of this scale could be protogalactic clumps, that is the building blocks of the various type of galaxies (ellipticals, spirals, etc.) that are observed at present epoch (Khersonsky & Turnshek (1996), Haehnelt et al. (1998)).

## 7 Conclusion

In the damped Ly$`\alpha `$ system at $`z=3.0543`$ of Q0000$``$2619 studied by SOM C ii absorption was detected. In this paper we report the detection at 3.5$`\sigma `$ of absorption from the excited fine structure level of C ii at $`z=3.0543`$. From the measure of the equivalent width of the two lines we derived an upper limit of 21.6 K on the fine structure excitation temperature. This value provides a strict upper limit on the temperature of the CMB at $`z=3.0543`$, which at this redshift is predicted to be 11.05 K. We then used the derived relative populations of the fine structure levels of C ii to set constraints on the absorber density and on the UV field in the absorbing cloud. Assuming an ionization degree ranging from 0 to 10%, a gas kinetic temperature between 100 and 10 000 K and a UV field with a Milky Way spectrum, the density of the absorber is constrained to be between 0.7 and 40 cm-3 and the H-ionizing flux between $`1\times 10^7`$ and $`8\times 10^8\mathrm{cm}^2\mathrm{s}^1`$, that is between 1 and 80 times the intensity of the Galactic UV field. The upper limits hold even if the detected absorption from C ii at $`z=3.0543`$ is seriously contaminated by Lyman forest absorption. If the damped Ly$`\alpha `$ system is assumed to be homogeneous, a density value between 0.7 and 40 cm-3 constrains the size of this absorber at $`z3`$ to be between 1 and 100 pc in the direction of the line of sight.

###### Acknowledgements.

We would like to thank S. Savaglio for providing the NTT spectra of Q0000$``$2619 and for her helpful comments to the manuscript, and P. Jakobsen for many illuminating discussions.


text:
# R-modes in Neutron Stars with Crusts: Turbulent Saturation, Spin-down, and Crust Melting

## 1. Introduction

The possibility of gravitational radiation from rapidly-rotating neutron stars has recently come under intense scrutiny. This interest has been fueled partly by the need to explain both the limited range of rotation periods in low-mass X-ray binary systems (LMXBs) with various accretion rates (Bildsten 1998; Andersson 1998) and the observed upper limit in spin rate of young neutron stars, and partly by the possibility that the gravitational radiation from LMXBs or newborn neutron stars may be detectable by gravity observatories such as LIGO (Owen et al. 1998; Bildsten 1998; Andersson, Kokkotas, & Stergioulas 1999). For neutron stars with crusts, there are two potential sources of gravitational radiation: a mass quadrupole originating in the solid crust of the star (Bildsten 1998), or a current quadrupole from Rossby waves, which become overstable under the influence of their own gravitational radiation (Andersson 1998; Friedman & Morsink 1998).

The outer layer of a neutron star solidifies if the typical energy of Coloumb interactions between nuclei, $`Z^2e^2/r`$, exceeds the thermal energy, $`k_bT`$, by a critical factor $`172`$ (Farouki & Hamaguchi 1993). Here, $`Z`$ and $`r`$ are the mean charge and the mean spacing of nuclei, respectively. The crystallization (also the melting) temperature is therefore

$$T_{\mathrm{melt}}=6.0\times 10^9Z_{20}^2A_{100}^{1/3}\mathrm{K},$$
(1)
where we have scaled $`Z`$ and $`A`$ (mean atomic weight) with values suitable for the bottom of the LMXB crust (Haensel & Zdunik 1990). Young neutron stars (Negele & Vautherin 1973), and LMXBs whose crust have been sufficiently heated, may have higher values of $`Z`$ and therefore higher $`T_{\mathrm{melt}}`$ Nucleation of a crystalline crust begins at around $`T_{\mathrm{melt}}`$, and crystallization accelerates exponentially as the temperature is decreased further (de Blasio 1995). All but the hottest neutron stars are expected to possess a solid crust of thickness $`1\mathrm{km}`$. Our investigation continues along the lines of Bildsten & Ushomirsky (2000; hereafter, BU) to explore the implications of such a solid crust for r-modes. We also consider the fate of the crust during r-mode instability.

We shall focus on the fastest growing r-mode with spherical indexes $`\mathrm{}=m=2`$, whose rotating frame frequency is $`\omega =2\mathrm{\Omega }_s/3`$. Here $`\mathrm{\Omega }_s`$ is the spin frequency of the neutron star; we denote its angular counterpart as $`\nu _s\mathrm{\Omega }_s/2\pi `$. In the fluid core, the r-mode is largely horizontal, and the modulus of its displacement vector is (Owen et al. 1998)

$$|\xi (r,\theta )|=\sqrt{\frac{45}{64\pi }}\alpha R\left(\frac{r}{R}\right)^2\left[1\mathrm{cos}^4\theta \right]^{1/2},$$
(2)
where $`\theta `$ is the co-latitude, $`r`$ the distance from the center and $`R`$ the stellar radius. The amplitude of the velocity perturbation is related to the displacement by $`|v|=\omega |\xi |`$. If the neutron star crust is perfectly rigid, the r-mode produces a periodic rubbing at the fluid-solid boundary with velocity $`|\mathrm{\Delta }v||v(R)|`$. However, Levin & Ushomirsky (2000) have argued that the crust of a neutron star is not perfectly rigid, and it will participate in the lateral motion of the mode. To account for this, we write $`|\mathrm{\Delta }v|=\eta |v|`$, with $`\eta =1`$ in the limit where the shear force in the crust much exceeds the local Coriolis force.

BU considered energy dissipation in a thin viscous boundary layer at the above core-crust interface. Dissipation due to “molecular” viscosity is more intense in this laminar boundary layer than in the bulk flow caused by the mode. Therefore, the presence of the boundary layer increases the stellar rotation rate above which the mode becomes unstable, relative to the case when the crust is absent (Owen et al. 1998). The mode is pumped by gravitational radiation and damped by viscosity at rates that both scale as $`\alpha ^2`$. So a laminar boundary layer cannot provide a saturation mechanism for an unstable r-mode.

In contrast, a turbulent boundary layer provides nonlinear dissipation that leads to r-mode saturation. A turbulent boundary layer occurs once the mode amplitude grows above a critical value. It removes kinetic energy from the mode with a rate that is cubic in $`\alpha `$. As $`\alpha `$ grows, turbulent dissipation rapidly catches up with the amplification of the mode energy due to gravitational back-reaction. We find this mechanism typically saturates the mode amplitude $`\alpha `$ at values much less than unity, i.e., well below what has been assumed in the literature (Levin 1999; Owen et al. 1998)

The kinetic energy of the r-mode is converted into heat in the thin turbulent boundary layer.111The heat input due to the viscous boundary layer is insignificant except when the r-mode is marginally unstable. If conduction and local neutrino emission cannot carry the heat away sufficiently quickly, the local temperature may increase above the melting temperature (eq. ) and the crust will begin to melt. In §3, we calculate the thermal profile induced by boundary-layer heating at the core-crust interface.

Finally, in §4, we consider the implications of turbulent saturation and local heating on the thermal and spin evolution of LMXBs as well as new-born neutron stars. We conclude in §5.

Throughout this paper we employ as a fiducial model an $`n=1`$ polytrope of mass $`M=1.4M_{}`$ and radius $`R=12.53\mathrm{km}`$ (as in Lindblom, Mendell, & Owen 1999). The bottom of the crust is assumed to have density $`\rho =1.5\times 10^{14}g\mathrm{cm}^3`$ (Brown 2000).

## 2. Viscous and turbulent boundary layers

Going outward from the fluid core to the solid crust, the shear modulus of the material changes from zero to a large (but finite) value. This implies a discontinuity $`\mathrm{\Delta }v`$ in the horizontal velocity of the r-mode across the core-crust boundary. The r-mode periodically rubs the fluid core against the solid crust with this velocity. A viscous boundary layer develops for small velocity jump, and a turbulent boundary layer sets in when the jump increases above some critical value.

The effect of a viscous boundary layer on r-modes was first analyzed by BU. We reiterate some of their results here. Ignoring the Coriolis force, the thickness of this layer is given by the diffusion length during one mode period (Stokes 1851),

$$\delta =\sqrt{2\nu /\omega },$$
(3)
where $`\nu `$ is the molecular viscosity at the base of the crust, $`\nu 1.8\times 10^4T_8^2\mathrm{cm}^2\mathrm{s}^1`$ for the chosen density (Flowers & Itoh 1979; Cutler & Lindblom 1987), and $`T_8`$ stands for $`T/10^8\mathrm{K}`$. The rate of energy dissipation per unit area is equal to the mean relative kinetic energy contained within this layer (half of $`\rho |\mathrm{\Delta }v|^2/2`$ times $`\delta `$) multiplied by the mode frequency:

$$\frac{d\dot{E}_{\mathrm{vbl}}}{dA}=\frac{1}{4}\rho |\mathrm{\Delta }v|^2\delta \omega .$$
(4)
This yields a global energy dissipation rate of

$$\dot{E}_{\mathrm{vbl}}0.03\mathrm{s}^1\mathrm{E}\eta ^2\mathrm{T}_8^1\nu _{\mathrm{kHz}}^{1/2},$$
(5)
where $`\nu _{\mathrm{kHz}}\nu _s/1\mathrm{kHz}`$, and $`E`$ refers to the energy in the mode (Owen et al. 1998),

$$E=\frac{1}{2}\alpha ^2\mathrm{\Omega }_s^2MR^2\stackrel{~}{J},$$
(6)
with $`\stackrel{~}{J}=0.016`$.222This expression for $`E`$ assumes that the r-mode satisfies equation (2) in the whole star; strictly speaking, an impenetrable crust reduces the mode energy by a factor of $`2`$. We ignore this correction throughout the paper. BU found that the viscous boundary layer stabilizes r-modes against gravitational back-reaction excitation at low spin frequencies and low stellar temperatures.

When an unstable r-mode acquires sufficiently large amplitude, the laminar boundary layer will make a transition to turbulence. This requires the local Reynolds number to exceed a critical value

$$\mathrm{}\frac{|\mathrm{\Delta }v||\mathrm{\Delta }\xi |}{\nu }\mathrm{}_{\mathrm{crit}},$$
(7)
where the amplitude of the relative displacement $`|\mathrm{\Delta }\xi |=|\mathrm{\Delta }v|/\omega `$ and the critical Reynolds number is determined experimentally to be $`\mathrm{}_{\mathrm{crit}}2\times 10^5`$ (Jensen, Sumer, & Fredsoe 1989).333Note that if we define a Reynolds number using the width of the turbulent layer (as is done in a pipe flow), as opposed to using $`\mathrm{\Delta }\xi `$, we find $`\mathrm{}_{\mathrm{crit}}500`$ at the transition to turbulence (Landau & Lifshitz 1959).

Combining the Reynolds number criterion (7) with the definition of $`|\mathrm{\Delta }\xi |`$, we find that the onset of turbulence at the equator requires a mode amplitude $`\alpha `$ greater than a critical value $`\alpha _{\mathrm{crit}}`$:

$$\alpha _{\mathrm{crit}}=1.6\times 10^3\eta ^1T_8^1\nu _{\mathrm{kHz}}^{1/2},$$
(8)
which corresponds to a thickness of the boundary layer $`d=31\mathrm{cm}T_8^1\nu _{\mathrm{kHz}}^1`$ (see Appendix A). At other latitudes, turbulence sets in at a larger mode amplitude.

Turbulence implies that eddies, rather than viscosity, advect momentum toward the crust. The crust feels a drag force very nearly in phase with $`\mathrm{\Delta }v`$ with a magnitude (eq. \[A6\])

$$\mathrm{Drag}C_D\rho |\mathrm{\Delta }v|^2/2,$$
(9)
where the drag coefficient $`C_D`$ is measured to be $`C_D5\times 10^3`$ (Jensen et al. 1989) at the onset of turbulence and varies logarithmically with the Reynolds number beyond the critical point, as is discussed in Appendix A. Kinetic energy is removed from the r-mode with a time-averaged rate per unit area of

$$\frac{d\dot{E}_{\mathrm{turb}}}{dA}=\frac{1}{2}|\mathrm{\Delta }v||\mathrm{Drag}|,$$
(10)
where time-averaging gives the factor $`1/2`$ as the drag is in phase with $`\mathrm{\Delta }v`$. Integrating over the sphere assuming a constant $`C_D`$, we find the total dissipation rate

$$\dot{E}_{\mathrm{turb}}\frac{3\pi }{4}\rho R^2C_D|\mathrm{\Delta }v_1|^3,$$
(11)
where $`|\mathrm{\Delta }v_1|=(5/16\pi )^{1/2}\alpha \mathrm{\Omega }_sR`$ is the velocity jump at the equator. In general, we adopt for $`C_D`$ the value at the equator, where the Reynolds number is the largest.

Turbulent dissipation increases with $`\alpha `$ more rapidly than does the pumping of the mode by gravitational back-reaction, and this leads to saturation. Recall that the mode gains energy due to gravitational wave radiation at the rate (cf. Owen et al. 1998)

$$\dot{E}_{\mathrm{gr}}=0.11\mathrm{s}^1\mathrm{E}\nu _{\mathrm{kHz}}^6.$$
(12)
Balancing this energy gain with the nonlinear energy drain from turbulence, we find a saturation amplitude

$$\alpha _{\mathrm{sat}}3.5\times 10^3\eta ^3\nu _{\mathrm{kHz}}^5\left(\frac{5\times 10^3}{C_D}\right).$$
(13)
Note that $`\alpha _{\mathrm{sat}}`$ is independent of temperature up to the logarithmic dependence of $`C_D`$ on $`\mathrm{}`$. Typically $`\alpha _{\mathrm{sat}}\alpha _{\mathrm{crit}}`$ for an unstable r-mode, in which case the r-mode grows rapidly until it reaches a steady state with $`\alpha =\alpha _{\mathrm{sat}}`$. The scaling $`\alpha _{\mathrm{sat}}\eta ^3\nu _{\mathrm{kHz}}^5`$ seems to imply that the saturation amplitude will increase enormously (at a given spin frequency) if the crust’s rigidity parameter $`\eta `$ is reduced from unity to $`0.1`$. Newborn, rapidly spinning neutron stars will therefore spin down at a rate that is very sensitive to $`\eta `$ (and may be saturated by other means if eq.  predicts $`\alpha _{\mathrm{sat}}>1`$). In contrast, the spins of LMXBs are confined within a narrow frequency range that varies with $`\eta `$ (as we will show in §4). The spin frequency at which they first become unstable scales as $`\nu _s\eta ^{4/11}`$ (equating eqs. with ). Taking this change of spin frequency into account, we find that the maximum r-mode amplitude in LMXBs varies as $`\eta ^{1.2}`$. Low saturation amplitudes are expected in LMXBs even for $`\eta 0.1`$.

We return now to comment on a few issues relevant to turbulence onset. The first is related to the local stratification. We find the Richardson number

$$𝒥\frac{N^2}{\left|dv/dr\right|^2}C_DN^2/\mathrm{\Omega }_s^2\frac{1}{4}$$
(14)
at $`\alpha =\alpha _{\mathrm{sat}}`$ for a Brunt-Väisälä frequency of $`N500\mathrm{rad}\mathrm{s}^1`$ (Reisenegger & Goldreich 1992). This indicates that the stable stratification at the core-crust boundary does not prevent the onset of turbulence.

Notice that the Reynolds number (eq. ) can also be written as $`\mathrm{}=4(|\mathrm{\Delta }\xi |/\delta )^22C_D(d/\delta )^2`$, where the thickness of the turbulent layer $`d`$ is related to $`|\mathrm{\Delta }\xi |`$ as $`d(C_D/2)^{1/2}|\mathrm{\Delta }\xi |`$ (see Appendix A). The condition $`\mathrm{}\mathrm{}_{\mathrm{crit}}`$ implies $`d\delta `$. Viscosity is unimportant in most of the turbulent region. Moreover, the turn-over time of the energy-bearing (also the largest) eddies in the turbulent region is much shorter than the mode period, $`|\mathrm{\Delta }v|/d\omega `$. This allows for a well-developed turbulent cascade.

The interaction between shear turbulence and stellar rotation, the possibility of equipartition magnetic fields in the turbulent layer, and the possible superfluid nature of the material, are potentially important questions that lie beyond the scope of this paper.

## 3. Temperature Profile from Boundary Layer Heating

Turbulence converts the kinetic energy of the r-mode into heat and deposits it in a thin region of width $`d100\mathrm{cm}\eta ^2\nu _{\mathrm{kHz}}^5`$ below the crust. If thermal conduction and local neutrino emission cannot carry away the heat efficiently, the temperature in this layer may be driven up to past the melting temperature of the crust (eq. ). As our saturation mechanism (§2) depends on the presence of such a solid crust, we wish to establish the temperature profile of the boundary layer. We assume a steady state solution because the conduction time over the regions of interest to us is much shorter than the evolution time.

The conduction of heat away from the boundary layer requires a temperature profile that declines in each direction. This, in turn, implies that the rate of local neutrino cooling decreases with distance away from the boundary layer. We assume that the scale length of the temperature perturbation is much shorter than the density scale height, and therefore much smaller than the stellar radius as well. This assumption, whose validity must be checked a posteriori, allows us to consider only planar heat conduction in a region of constant density.

The following toy model gives insight into the full solution. Let us consider a piecewise linear temperature profile, one that reaches a temperature $`T_{\mathrm{bl}}`$ within the boundary layer (of width $`d`$), and falls to zero over lengths $`l_{\mathrm{core}}`$ and $`l_{\mathrm{crust}}`$ to the inside and outside, respectively. The two lengths differ because the thermal conductivity takes different values in the two directions.

In a state of equilibrium, any energy not radiated as neutrinos within the boundary layer must be carried away by conductive fluxes:

$$\frac{d\dot{E}_{\mathrm{turb}}}{dA}=dϵ_0T_{\mathrm{bl}}^8+F_{\mathrm{core}}+F_{\mathrm{crust}},$$
(15)
where neutrino emissivity due to modified Urca reactions is $`ϵ_\nu ϵ_0T^8\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1`$ with $`ϵ_0=7.4\times 10^{52}`$ (Friman & Maxwell 1979). The heat flux leaving the boundary layer in each direction is

$$F_{\mathrm{core}}=\kappa _{\mathrm{core}}\frac{T_{\mathrm{bl}}}{l_{\mathrm{core}}}\mathrm{and}F_{\mathrm{crust}}=\kappa _{\mathrm{crust}}\frac{T_{\mathrm{bl}}}{l_{\mathrm{crust}}},$$
(16)
respectively. The thermal conductivity of the crystalline crust is dominated by electron-phonon scattering and scales with the local temperature inversely, $`\kappa _{\mathrm{crust}}=\overline{\kappa }_{\mathrm{crust}}T_8^110^{20}T_8^1\mathrm{erg}\mathrm{cm}^1\mathrm{s}^1\mathrm{K}^1`$, while that of the liquid core is dominated by electron-proton and electron-neutron scattering with the same temperature dependence, $`\kappa _{\mathrm{core}}=\overline{\kappa }_{\mathrm{core}}T_8^110^{23}T_8^1\mathrm{erg}\mathrm{cm}^1\mathrm{s}^1\mathrm{K}^1`$ (Flowers & Itoh 1979).

Another consequence of thermal equilibrium is that the heat conducted into the crust or core must be radiated there in neutrinos. With our approximation of a linear temperature profile, this implies $`Fϵ_0lT_{\mathrm{bl}}^8/9`$ in each direction with the conductive lengths $`l`$ defined in equation (16). A more accurate calculation in which no temperature profile is assumed a priori (Appendix B) gives $`F=ϵ_0lT_{\mathrm{bl}}^8/4`$. Combining this with equation (16), we find the conductive lengths to be

$`l_{\mathrm{core}}`$ $``$ $`2.4\times 10^5\mathrm{cm}\left(10^9\mathrm{K}/T_{\mathrm{bl}}\right)^4,`$
$`l_{\mathrm{crust}}`$ $``$ $`7.4\times 10^3\mathrm{cm}\left(10^9\mathrm{K}/T_{\mathrm{bl}}\right)^4.`$ (17)
As long as $`T_{\mathrm{bl}}5\times 10^8\mathrm{K}`$, $`l_{\mathrm{crust}}<l_{\mathrm{core}}R`$ and our assumption of a thin heated region is justified. The heating becomes increasingly local as the boundary layer rises in temperature: when more heat is deposited in the boundary layer, a larger temperature gradient is required to conduct the heat outward leading to a smaller conduction length. The assumption of an isothermal star is thus not justified in the presence of any significant localized heating. In §4.1 and §4.2 we shall compare the implications of the isothermal approximation with those of localized temperature profile for neutron star spin evolution.

Substituting equations (16) and (17) into equation (15), adopting $`\alpha =\alpha _{\mathrm{sat}}`$ and evaluating $`d\dot{E}_{\mathrm{turb}}/dA`$ at the equator, we find a quadratic equation for $`T_{\mathrm{bl}}^4`$, whose solution gives

$$T_{\mathrm{bl}}=1.5\times 10^{10}\mathrm{K}\eta ^{1/2}\nu _{\mathrm{kHz}}^{13/8}\left(\sqrt{1+\mathrm{b}^2}\mathrm{b}\right)^{1/4},$$
(18)
where $`b=7.4\times 10^3\eta ^4\nu _{\mathrm{kHz}}^{23/2}`$ and $`C_D`$ is taken to be $`0.005`$. Neutrino losses from the boundary layer itself dominate over the conductive fluxes for $`b1`$, or $`\nu _s650\mathrm{Hz}\eta ^{8/23}`$. For a given mode amplitude, the boundary layer temperature decreases towards the poles.

We find that the flux going to the core is roughly $`(\kappa _{\mathrm{core}}/\kappa _{\mathrm{crust}})^{1/2}30`$ larger than that conducted to the crust. However, if the above equilibrium temperature exceeds $`T_{\mathrm{melt}}`$, the crust will begin to melt. Once melting begins, we expect that the energy input from the boundary layer will be used almost entirely to melt the crust. The temperature in this layer will be kept very close to $`T_{\mathrm{melt}}`$ by the continual inflow of cool, freshly molten crust, as the boundary layer moves upward following the bottom of the receding crust. The time it takes to completely melt the crust can be estimated to be a few days (the ratio $`10^{48}\mathrm{erg}/\dot{\mathrm{E}}_{\mathrm{gr}}`$ where $`10^{48}\mathrm{erg}`$ is roughly the chemical binding energy of the crust).

## 4. Implications for Neutron Star Spin and Temperature Evolution

We now consider the implications of turbulent saturation for the evolution of young neutron stars and LMXBs in the $`(\mathrm{\Omega }_s,T)`$ plane. In this section we assume the crust is present and is not melted, an assumption we examine in §4.2. Once the neutron star is pushed into the instability region, either by cooling or by accretion, its r-mode quickly grows to its saturated value as given by equation (13). Inserting it into equation (C2) we find that the time to spin down to a frequency $`\nu _{\mathrm{kHz}}`$ is fairly independent of the initial spin frequency of the star, and is given by

$$t=5.2\times 10^5\mathrm{s}\eta ^6\nu _{\mathrm{kHz}}^{16}.$$
(19)
The time to spin the star down to, say, $`\nu _s=500\mathrm{Hz}`$ is roughly $`10^3\eta ^6\mathrm{yr}`$. As $`\alpha _{\mathrm{sat}}`$ decreases sharply with decreasing frequency, most of the spin-down time is spent near the lowest frequencies.

The above conclusion on the evolution of the spin frequency may need modification if the crust can be melted by the r-mode heating. This prompts us to study the temperature evolution of the star. This we do by first adopting the isothermal approximation for the whole star. It applies in the limit of infinite thermal conductivity or large-scale heating. We then apply our results from §3 to the realistic case in which the conductivity is finite and r-mode heating is local.

Evolution in the unrealistic isothermal case can be easily integrated using equations in Appendix C. It is useful for comparison with previous work. More importantly, insights gained from the isothermal case can be applied to the realistic case.

### 4.1. Isothermal Approximation

After the r-mode reaches its saturation amplitude, it provides the neutron star with a source of heat from turbulent dissipation. The neutron star heats up quickly until it reaches a state in which the r-mode heating is balanced by the neutrino losses. With time the temperature will decrease slowly because heating by r-mode decreases with the spin frequency. We set the r-mode heating term equal to the neutrino cooling term in equation (C3) to find the following scaling between the temperature and the spin rate of the star,

$$T_8=47\eta ^{3/4}\nu _{\mathrm{kHz}}^{9/4},$$
(20)
where we take $`\alpha =\alpha _{\mathrm{sat}}`$ and ignore the logarithmic dependence of $`C_D`$ on temperature by using $`C_D=0.005`$. This relation gives us the thermal equilibrium spin-down line. It intersects the instability curve at the terminal frequency $`\nu _T520\mathrm{Hz}\eta ^{0.35}`$. All neutron stars entering the instability region with an initial spin frequency $`\nu _s\nu _T`$ will converge onto the equilibrium line and exit at $`\nu _T`$. The ones that enter at $`\nu _s<\nu _T`$ find their spin rates hardly affected by r-mode instability. This is because their spin down time (eq. ) much exceeds their neutrino cooling time.

Figures 1 and 2 exhibit the results of our evolutionary calculations with a variety of mass accretion rates in LMXBs and a range of initial spin rates for young neutron stars. The evolution equations we adopted for the integrations are listed in Appendix C. They are similar to those derived in Owen et al. (1998) apart from the additional terms that arise from boundary layer dissipation. We consider only $`\eta =1`$ here. We find that, up to logarithmic dependence of $`C_D`$ on the Reynolds number, the numerical results confirm equations (19) and (20).

Note that the presence of an equilibrium spin-down line is a universal feature independent of the saturation mechanism. If r-modes are saturated to constant values, for instance, we find the equilibrium temperature $`T\nu _s`$.

Levin (1999) was the first to point out that LMXBs undergo limit cycles in which the angular momentum accreted over the whole cycle is radiated by gravitational waves when the r-mode is unstable. The r-mode phase lasts a time that depends on the saturation mechanism and is relatively insensitive to the initial spin frequency at which the r-mode is destabilized. Equation (19) yields the spin-down time at the terminal frequency $`\nu _T`$ to be $`10^3\eta ^{10/31}\mathrm{yr}`$. Compared to an accretion spin-up time of $`10^7\mathrm{yr}`$ at the Eddington accretion rate, we find a maximum duty cycle of $`10^4`$ during which time the gravitational wave emitted by the r-mode may be potentially detectable, with the duty cycle decreasing at lower accretion rate. This duty cycle is much higher than that found by previous investigators (e.g., Levin 1999) mostly due to our much reduced saturation amplitude.

### 4.2. Realistic Evolution Accounting for Local Heating

For simplicity, we consider only the evolution of LMXBs in this section. Numerical integration for the evolution tracks in the realistic case is much more difficult than in the isothermal case. Fortunately, they are not necessary. As is evident from Figure 1, given initial conditions, the equilibrium line and the instability curve completely determine the evolution tracks.

The r-mode instability curve in the low temperature range depends only on the boundary layer temperature because damping from the viscous boundary layer dominates over other damping mechanisms. The equilibrium spin-down line is also given in terms of the boundary layer temperature (eq. ). So we can simplify the problem in hand by considering only the temperature at the boundary layer.

There are two main differences between the realistic and the isothermal case: a lower terminal frequency and an enhanced heating in the boundary layer.

The new equilibrium spin-down line (eq. ) intersects the instability curve at a terminal frequency

$$\nu _T490\eta ^{0.35}\mathrm{Hz}$$
(21)
that is lower than in the isothermal case. It takes slightly longer to spin the star down to this lower value of $`\nu _T`$.

The r-mode heats the boundary layer to a higher temperature than under the isothermal approximation. This results from the fact that, unlike in the isothermal case, there is now only a shell of thickness $`(l_{\mathrm{core}}+l_{\mathrm{crust}})/4+d`$ which receives heat and cools by neutrinos (§3).

We present the LMXB limit cycles for the two cases $`\eta =1`$ and $`\eta =0.1`$ in Figure 3. One observational consequence suggests itself in the figure: all LMXBs are expected to have spin frequencies falling in the narrow range of $`490700\mathrm{Hz}`$ when $`\eta =1`$, and $`220300\mathrm{Hz}`$ when $`\eta =0.1`$, independent of their accretion history. This continues the idea proposed by Bildsten (1998) and Andersson et al. (1999) that r-modes may be instrumental in halting the LMXB spin-up. The resemblance between the above frequency range for $`\eta =0.1`$ and the observed LMXB spin rates (van der Klis 2000) is intriguing. However, it may be difficult to explain the fastest millisecond pulsars (period of $`1.5`$ ms) using the same value of $`\eta `$.

The above results may need modification if the crust melts during the r-mode evolution. Melting the crust is more likely than is suggested by the isothermal approximation. As Figure 3 demonstrates, the crust is most likely to melt if $`\eta `$ is large and if the accretion rate is low. We define a melting frequency, $`\nu _m`$, at which the equilibrium line intersects the melting line $`T=T_{\mathrm{melt}}`$. We find for $`T_{\mathrm{melt}}=6\times 10^9\mathrm{K}`$, $`\nu _m=660\mathrm{Hz}`$ for $`\eta =1`$ and $`315\mathrm{Hz}`$ for $`\eta =0.1`$. We speculate below on some possible evolutionary consequences when crust melting is taken into account.

### 4.3. Possible Outcomes Including Crust Melting and Forming

The rotational energy ($`10^{51}\nu _{\mathrm{kHz}}^2\mathrm{erg}`$) far exceeds the chemical binding energy of the crust ($`10^{48}\mathrm{erg}`$) for spin rates of interest: crust melting is not inhibited by lack of energy. Moreover, as discussed in §3, the temperature at the bottom of the crust may be driven above the melting temperature by turbulent heating; this causes melting. For our fiducial parameters, this occurs for slowly accreting LMXBs if $`\eta =1`$ (Fig. 3).

After the crust is melted, some other nonlinear dissipation must cause the r-mode to saturate. It is beyond the scope of this paper to predict the evolution when the crust is molten. However, we can be confident that the star will cool and spin down sufficiently to form a crust again. Once this occurs, the subsequent evolution is again within our jurisdiction, and it depends on the stellar spin frequency ($`\nu _e`$) at the moment the crust forms. We discuss five representative evolution scenarios that begin with different values of $`\nu _e`$, labeled A through E in Figure 4.

In case A, the stellar temperature drops below the melting temperature while the star is still spinning so fast ($`\nu _e>\nu _m`$) that the boundary layer temperature would exceed $`T_{\mathrm{melt}}`$ if a crust did form. Therefore, the material that could form a crust is kept hot (by some means) with $`TT_{\mathrm{melt}}`$. This implies that the formation of an appreciable crust is delayed until the rotation has slowed to $`\nu _s=\nu _m`$, where $`\nu _m`$ is the frequency at which the equilibrium line hits the melting line. After this point the star will spin down with a saturated r-mode (eq. ) and cool in a state of thermal equilibrium (eq. ); this behavior is represented by the solid line in Figure 4. The star leaves the instability region with a specific (terminal) spin frequency when the mode is no longer unstable.

Case B represents the outcome if the initial spin frequency $`\nu _e`$ is low enough to allow crust formation, yet higher than the terminal spin frequency. A crust forms and the star evolves towards the equilibrium spin-down line, with the the cooling of the boundary layer regulated by that of the whole star. The subsequent evolution is identical to case A.

Case C enters with $`\nu _s<\nu _T`$. Its spin-down time is longer than the cooling time and it exits the r-mode instability curve without being significantly spun down. Case D has a similar outcome: r-mode instability is quenched entirely once the crust forms.

Case E hits the instability curve to the right of the melting line. It evolves along the critical stability curve until the melting line, after which the r-mode becomes stable. This is the only scenario in which we can predict the evolution of an r-mode in the absence of a crust.

From these considerations we conclude that the final rotation rates of young neutron stars, and of those LMXBs capable of melting their crusts, depend on their spin evolution while the crust is absent. In the case that $`\eta =1`$ and $`T_{\mathrm{melt}}=6\times 10^9\mathrm{K}`$, the final spin frequencies should lie in the range of $`50`$ to $`480\mathrm{Hz}`$, where the upper limit is the terminal frequency $`\nu _T`$. The chosen value of $`\eta `$ should be considered uncertain, but recent calculations (Levin & Ushomirsky 2000) point to $`\eta 0.1`$. Also, a minimum value of $`\eta 0.003`$ can be derived from equation (21) by assuming that $`\nu _T`$ is larger than the frequency of the fastest known young neutron star ($`16\mathrm{m}\mathrm{s}`$; Marshall et al. 1998). Moreover, young neutron stars, or LMXBs that have been heated above the temperature necessary for nuclear burning, are more likely to be in a state of chemical equilibrium at the base of their crusts; this implies higher values of $`Z`$, thus higher melting temperatures. However, the above discussion demonstrates that this would only raise the lower bound for the possible range of final rotation rates without changing the upper bound ($`\nu _T`$).

## 5. Conclusions

R-modes in the fluid core of a neutron star generate velocity shear below the solid crust. The primary conclusion of this paper is that this velocity shear can drive turbulence which in turns limits the amplitude that an overstable r-mode can obtain. This saturation amplitude is easily calculable. It rises steeply with mode frequency and depends on the temperature at the turbulent boundary layer only logarithmically. For most spin frequencies of interest, we find this amplitude falls well below unity.

We also find that heating due to turbulent dissipation in the boundary layer leads to a localized temperature enhancement relative to the surrounding material, and the heat is removed by neutrinos from both the boundary layer itself and from the nearby regions to which the heat is conducted. It is erroneous to assume that the whole star remains isothermal in the presence of a localized heat source. The boundary layer temperature may exceed the melting temperature of the crust.

When studying the spin and temperature evolution of neutron stars undergoing r-mode instability, we find the following remarkable behaviors: because of the strong dependence of neutrino cooling on the temperature, the temperature of the boundary layer rapidly reaches an equilibrium value that depends on the mode amplitude and spin frequency, and hence, for a saturated mode, only on the stellar spin. This defines an equilibrium spin-down line along which all neutron stars evolve. We define the spin frequency at which this line intersects the instability line as $`\nu _T`$, the terminal frequency. This frequency is a significant feature of the spin evolution. We expect all neutron stars with crusts to exit the r-mode instability range with $`\nu _s=\nu _T`$, as long as their initial spin frequencies are higher than $`\nu _T`$. Neutron stars born with $`\nu _s<\nu _T`$ will never experience significant spin-down. If the crust can be melted and reformed, the terminal frequency defines the upper limit to the final spin rates of these neutron stars. We also establish a spin-down law for which the time required to spin the star down to $`\nu _T`$ is roughly independent of the initial spin rate. These features are expected to persist for other saturation mechanisms.

This work has relied on several assumptions and approximations which should be noted. We have ignored any dissipation that may occur within the solid crust itself, the possibility that the crust may be melted by mechanical strain or heat from nuclear reactions, and additionally any nonlinear dissipation caused by coupling between modes in the bulk of the fluid core (Schenk et al. 2000). We have ignored the possibility that gravitational radiation may be dominated by a permanent mass quadrupole within the crust. We have not considered the effects of superfluidity in the material below the crust. The detailed value of the boundary layer temperature depends on the exact value of the thermal conductivity as well as on the parameter $`\eta `$ that describes the crust motion; these we consider relatively uncertain. Moreover, melting of the crust is sensitive to the composition at its base, and this may differ between young and accreting neutron stars. Finally, our discussions in §4.3 do not include the possibility that chemical composition in the crust may change when melting occurs.

We wish to thank Peter Goldreich, Norm Murray, Greg Ushomirsky and Marten van Kerkwijk for comments and encouragements. Each of the authors received support from an NSERC fellowship. CDM was also supported by a visiting postdoctoral fellowship at Caltech; he wishes to thank Sterl Phinney and Roger Blandford for their gracious hospitality during this visit. PA wishes to thank Katrin Schenk, Ira Wasserman, Eanna Flanagan, and Saul Teukolsky for useful conversations.

## Appendix A Turbulent Boundary Layer

In this appendix, we give a short description of the turbulent boundary layer caused by shear flows. We provide scalings for the thickness of the layer, the velocity profile inside the layer and the drag force acting on the wall. These scalings are confirmed by numerical experiments and we use them in the main text to derive the energy dissipation rate due to shear turbulence.

Let $`\mathrm{\Delta }v`$ be the oscillatory rubbing velocity between the neutron star’s fluid core and its solid crust. The horizontal displacement is $`\mathrm{\Delta }\xi =\mathrm{\Delta }v/\omega `$. When the boundary layer becomes turbulent, we can describe it with results from a similar situation: the turbulent boundary layer generated by a steady flow along a semi-infinite plane. Our discussion therefore follows closely that of Landau & Lifshitz (1959, §42-44). In this case, $`\mathrm{\Delta }\xi `$ can be understood as the distance away from the edge of the plane, and $`\mathrm{\Delta }v`$ as the velocity of the shear flow. Although the steady solution is expected to reproduce the correct scalings for the boundary layer structure, it is time-independent and therefore cannot predict the phase shift between the shear velocity and the drag force. We obtain this phase shift from the experimental data of Jensen et al. (1989).

The thickness of the turbulent boundary layer ($`d`$) grows linearly with the distance away from the edge of the plane ($`\mathrm{\Delta }\xi `$). We write

$$d=\lambda \mathrm{\Delta }\xi ,$$
(A1)
where $`\lambda `$ is a small dimensionless number that describes the expansion of the turbulent layer relative to the horizontal displacement away from the edge. If we define $`v_{}`$ to be the root mean squared fluctuating velocity in the turbulent layer, we find $`\lambda =v_{}/\mathrm{\Delta }v`$

Within the boundary layer, the time-averaged velocity $`u`$ at distance $`y`$ away from the plane obeys a logarithmic profile,

$$u(y)=\frac{v_{}}{c_1}\mathrm{ln}\left(\frac{c_2v_{}y}{\nu }\right),$$
(A2)
where $`c_10.4`$ and $`c_27.69`$ are experimentally determined constants. At $`y=d`$, we should recover the flow velocity, so $`\mathrm{\Delta }v=v_{}c_1^1\mathrm{ln}(c_2v_{}d/\nu )`$. Combining this expression with equation (A1), we find the following solution for $`\lambda `$,

$$\frac{c_1}{\lambda }=\mathrm{ln}\left(c_2\lambda ^2\mathrm{}\right),$$
(A3)
where we define an experimentally controlled Reynolds number for the flow by

$$\mathrm{}\frac{\mathrm{\Delta }v\mathrm{\Delta }\xi }{\nu }.$$
(A4)
A rough scaling for $`\lambda `$ can be found in the limit of very large Reynolds number where the solution to equation (A3) is

$$\lambda \frac{c_1}{\mathrm{ln}(c_2\mathrm{})}.$$
(A5)
The fact that $`\lambda `$ is a slowly varying function of the Reynolds number helps simplifying the calculation for energy dissipation rate.

The drag force per unit area on the surface due to the fluctuating velocity field is

$$\mathrm{Drag}=\rho v_{}^2=\lambda ^2\rho \mathrm{\Delta }v^2\frac{1}{2}C_D\rho \mathrm{\Delta }v^2,$$
(A6)
which identifies the drag coefficient in the turbulent regime to be $`C_D=2\lambda ^22[c_1/\mathrm{ln}(c_2\mathrm{})]^2`$, also a slowly decreasing function of the Reynolds number. Jensen et al. (1989) show that in the fully turbulent regime, the drag force is nearly in phase with the shear velocity.

The scalings $`C_D=2\mathrm{}^{1/2}`$ in the laminar regime (Landau & Lifshitz 1959) and $`C_D=2\lambda ^2`$ in the turbulent regime compare very well with the experimental determinations of $`C_D`$ by Jensen et al. (1989) over a large range of Reynolds number. Furthermore, if we estimate the critical Reynolds number for the onset of turbulence by equaling the two $`C_D`$ expressions, we find $`\mathrm{}_{\mathrm{crit}}=1.6\times 10^5`$, in good agreement with the experiment. At $`\mathrm{}=\mathrm{}_{crit}`$, $`\lambda =0.05`$ and $`C_D=0.005`$. Beyond the onset of turbulence, $`C_D`$ decreases logarithmically with the Reynolds number.

## Appendix B Equilibrium Temperature Profile Around a Localized Source of Heat

We assume thermal equilibrium has been reached. On either side of a localized heat source, the heat conduction and energy conservation equations read,

$`\stackrel{}{F}`$ $`=`$ $`\kappa \stackrel{}{}T,`$
$`\stackrel{}{F}`$ $`=`$ $`ϵ_\nu ,`$ (B1)
where $`\kappa `$ is the thermal conductivity and scales inversely with $`T`$ in the neutron star interior, while $`ϵ_\nu `$ is the neutrino emissivity with $`ϵ_\nu T^8`$ for modified Urca reactions. We define two constants $`\kappa _0=\kappa T`$ and $`ϵ_0=ϵ_\nu /T^8`$.

Under the planar approximation, we manipulate the above equations to find

$$\frac{d\mathrm{ln}T}{dr}\frac{d^2\mathrm{ln}T}{dr^2}=\frac{ϵ_0}{\kappa _0}T^8\frac{d\mathrm{ln}T}{dr}.$$
(B2)
where $`r`$ is the radius. Integrating both sides, we find

$$\left(\frac{d\mathrm{ln}T}{dr}\right)^2|_b^a=\frac{ϵ_0}{4\kappa _0}T^8|_b^a,$$
(B3)
where $`a`$ and $`b`$ are two arbitrary points within the integration range. Here, we take point $`a`$ to be at the boundary layer, which is the source of heat, and point $`b`$ to be at many conduction lengths away from the boundary layer so that the flux and temperature at point $`b`$ may be regarded as negligible. The flux and temperature at the boundary are then related by

$$F=\frac{1}{2}(ϵ_0\kappa _0)^{1/2}T_{\mathrm{bl}}^4.$$
(B4)
If we define a conduction length $`l=\kappa T_{\mathrm{bl}}/F=\kappa _0/F`$, we find $`F=ϵ_0lT_{\mathrm{bl}}^8/4`$, a result we adopt in §3.

## Appendix C Evolution Equations

We have modified the phenomenological model of Owen et al. (1998) to calculate the effect of the r-mode on the spin and temperature of the neutron star in the isothermal approximation. The evolution equations for mode amplitude $`\alpha `$, spin frequency $`\mathrm{\Omega }_s`$, and stellar temperature $`T`$ are

$`\dot{\alpha }`$ $`=`$ $`{\displaystyle \frac{\alpha }{\tau _{\mathrm{gr}}}}{\displaystyle \frac{\alpha }{\tau _v}}{\displaystyle \frac{1\alpha ^2Q}{1+\alpha ^2Q}},`$ (C1)
$`\dot{\mathrm{\Omega }_s}`$ $`=`$ $`{\displaystyle \frac{2\mathrm{\Omega }_s}{\tau _v}}{\displaystyle \frac{\alpha ^2Q}{1+\alpha ^2Q}}+{\displaystyle \frac{1}{\stackrel{~}{I}}}\left({\displaystyle \frac{4}{3}}\right)^{1/2}{\displaystyle \frac{\dot{M}}{M}}\mathrm{\Omega }_d,`$ (C2)
$`C_TT_9\dot{T}_9`$ $`=`$ $`L_{\nu ,9}T_9^8+ϵ_a(\dot{M}/M)(GM^2/2R)+2E/\tau _{\mathrm{heat}},`$ (C3)
where $`Q=0.094`$, $`\stackrel{~}{I}=0.261`$, $`\dot{M}`$ is the mass accretion rate onto the surface of the neutron star (equals $`0`$ in the young neutron star case), $`C_T=1.4\times 10^{48}\mathrm{erg}\mathrm{s}^1`$ is twice the specific heat at a temperature of $`10^9K`$, and $`L_{\nu ,9}=7.4\times 10^{39}\mathrm{erg}\mathrm{s}^1`$ is the neutrino cooling rate at $`T=10^9K`$. The terms involving $`\dot{M}`$ describe the torque and heating by accretion (Levin 1999), and $`ϵ_a`$ is the fraction of accreted energy that goes into heating the star, we take it to be unity. The mode energy $`E=0.5\alpha ^2\mathrm{\Omega }_s^2MR^2\stackrel{~}{J}`$ with $`\stackrel{~}{J}0.016`$. The initial LMXB temperature is determined self-consistently using equation (C3) setting $`\alpha =0`$ and $`\dot{T}=0`$. It varies by a factor of a few over three orders of magnitude in accretion rates. The various time-scales in the above equations are listed below.

Following Lindblom et al. (1999), we write the energy input rate due to gravitational wave back-reaction to be

$$\frac{1}{\tau _{\mathrm{gr}}}\frac{\dot{E}_{\mathrm{gr}}}{2E}=\frac{1}{\stackrel{~}{\tau }_{\mathrm{gr}}}\nu _{\mathrm{kHz}}^6,$$
(C4)
with $`\stackrel{~}{\tau }_{\mathrm{gr}}=18.7\mathrm{s}`$ and $`\nu _{\mathrm{kHz}}=\nu _s/1\mathrm{kHz}`$. The negative sign indicates driving. Similarly, the bulk viscosity and molecular shear viscosity rates are

$$\frac{1}{\tau _{\mathrm{bv}}}=\frac{1}{\stackrel{~}{\tau }_{\mathrm{bv}}}\nu _{\mathrm{kHz}}^2T_8^6,\mathrm{and}\frac{1}{\tau _{\mathrm{sv}}}=\frac{1}{\stackrel{~}{\tau }_{\mathrm{sv}}}T_8^2,$$
(C5)
with $`\stackrel{~}{\tau }_{\mathrm{bv}}=3.5\times 10^{17}\mathrm{s}`$ and $`\stackrel{~}{\tau }_{sv}=2.5\times 10^6\mathrm{s}`$. We will continue to use these estimates in the presence of a crust even though the r-mode occupies a smaller volume.

BU and this article have pointed to the importance of the boundary layers at the core-crust interface. We model the dissipation rate in this layer as a function that is continuous when the layer changes from laminar (viscous) to turbulent,

$$\frac{1}{\tau _{\mathrm{bl}}}=\mathrm{\Theta }(\alpha _{\mathrm{eq}}\alpha )\frac{1}{\tau _{\mathrm{vbl}}}+\mathrm{\Theta }(\alpha \alpha _{\mathrm{eq}})\frac{1}{\tau _{\mathrm{tbl}}},$$
(C6)
with $`\mathrm{\Theta }(x)=1`$ when $`x>0`$ and $`0`$ otherwise. The viscous and the turbulent damping rates are, respectively,

$$\frac{1}{\tau _{\mathrm{vbl}}}=\frac{1}{\stackrel{~}{\tau }_{\mathrm{vbl}}}\eta ^2\nu _{\mathrm{kHz}}^{1/2}T_8^1\mathrm{and}\frac{1}{\tau _{\mathrm{tbl}}}=\frac{1}{\stackrel{~}{\tau }_{\mathrm{tbl}}}\eta ^3C_D\alpha \nu _{\mathrm{kHz}},$$
(C7)
where the fiducial time-scales $`\stackrel{~}{\tau }_{\mathrm{vbl}}=65\mathrm{s}`$ and $`\stackrel{~}{\tau }_{\mathrm{tbl}}=3.3\times 10^4s`$. We equate $`\tau _{\mathrm{vbl}}`$ with $`\tau _{\mathrm{tbl}}`$ to obtain the swith-over amplitude $`\alpha _{\mathrm{eq}}=1.0\times 10^3\eta ^1\nu _{\mathrm{kHz}}^{1/2}T_8^1`$. In our calculations, we evaluate $`C_D`$ using the Reynolds number at the equator.

The total rate of viscous damping for the r-mode as well as the rate of heating due to the r-mode (ignoring bulk viscosity heating) are given by

$$\frac{1}{\tau _v}=\frac{1}{\tau _{\mathrm{bv}}}+\frac{1}{\tau _{\mathrm{sv}}}+\frac{1}{\tau _{\mathrm{bl}}}\mathrm{and}\frac{1}{\tau _{\mathrm{heat}}}\frac{1}{\tau _{\mathrm{sv}}}+\frac{1}{\tau _{\mathrm{bl}}}.$$
(C8)