text:
# A Comprehensive Statistical Analysis of the Gas Distribution in Lyman-limit and Damped Lyman-𝛼 Absorption Systems

## 1 Introduction

Absorption features in quasar spectra represent a powerful tool to investigate the formation and evolution of gaseous structures in the early Universe (Rauch, 1998, and references therein), the efficiency of processes such as merging (Kauffmann, 1996), gas depletion due to star formation (Wolfe et al., 1995; Storrie-Lombardi and Wolfe, 2000), and ionization by the ultraviolet background radiation field (Weinberg et al., 1997; Savaglio et al., 1997). Also the level of cloud clustering may be investigated analyzing absorption spectra (Cristiani et al., 1997). Lyman-limit absorption systems (hereafter LLS) are cosmological structures which contain enough neutral hydrogen to absorb Lyman continuum photons and produce a break in the QSO continuum flux level. This feature is detectable in moderate-resolution spectra whenever the opacity to Lyman continuum photons satisfies the condition $`\tau _{LL}1`$, equivalent to a neutral hydrogen column density $`N_{\text{H}\text{I}}1.6\times 10^{17}\mathrm{cm}^2`$. The column density of LLS can be determined whenever a residual flux beyond the Lyman break is detectable. Systems with $`N_{\text{H}\text{I}}\stackrel{>}{}5\times 10^{19}\mathrm{cm}^2`$ are also visible in moderate-resolution spectra since they give rise to Damped Lyman-$`\alpha `$ absorption lines with a reference frame equivalent width $`5\text{Å}`$. In this paper we shall refer to these as Damped Lyman-$`\alpha `$ absorption systems (hereafter DLS). The HI column density of absorbers is usually very poorly determined observationally when $`5\times 10^{17}<N_{\text{H}\text{I}}<5\times 10^{19}\mathrm{cm}^2`$ since the Lyman break is saturated, while the Lyman-$`\alpha `$ absorption line is not yet damped. Previous studies of absorbers with poorly determined column densities have been carried out using a coarse binning in $`N_{\text{H}\text{I}}`$. But a coarse binning suffers by statistical problems, related to the fact that different choices of binning may lead to different results (e.g. Stengler-Larrea et al. (1995) for problems related to redshift binning), as well as by some negative consequences on the physical description, because it does not differentiate HI column densities which correspond to different ionization conditions. This is particularly important for $`10^{17}<N_{\text{H}\text{I}}<10^{20}\mathrm{cm}^2`$ since in this range a smoothly varying total gas column density distribution leads to a rapidly varying HI column density distribution due to the rapid change of the $`H`$ ionization fraction (Corbelli et al., 2000, hereafter Paper I).

To tackle this problem, we have collected data on LLS and DLS from all the available literature and perform a comprehensive statistical analysis, based on a Maximum Likelihood procedure, which allows a joined fit of observations with different $`N_{\text{H}\text{I}}`$ sensitivities and uncertainties, with different ionization conditions and at different redshifts.

The outline of this paper is the following: in Section 2 we discuss how the observed HI column density distribution is related to the total hydrogen distribution for face-on absorbers. A set of analytical solutions for the cross section of non-spherical absorbers, randomly oriented in space, is reported in Appendix A. Section 3 contains a description of the database. Section 4 details the statistical algorithm that we use to best fit the HI column density distribution of absorbers. Section 4 also summarizes results relative to our database. In Section 5 we perform realistic simulations of the data, which are used to test the effectiveness and the robustness of the statistical algorithm.

## 2 Column densities relationships

One physically meaningful quantity for investigating the formation and evolution of gas condensations is the distribution of their gaseous mass. However observations measure only the distribution of the neutral hydrogen column density along the line of sight, which may present a rather complex behavior even if the face-on total column density distribution is a power law. In this section we show how to transform one distribution into the other. The relationship between the total and the neutral column density has been investigated in detail in Paper I. Here we focus on the effect of randomly oriented absorbers which are either flat or of ellipsoidal shape.

### 2.1 From the total to the neutral column density distribution

We show first how a face-on total hydrogen column density distribution $`g_{}(N_\text{H})`$ translates into a face-on neutral hydrogen distribution $`f_{}(N_{\text{H}\text{I}})`$. For a generic relationship $`N_\text{H}\left(N_{\text{H}\text{I}}\right)`$ one can write:

$$f_{}\left(N_{\text{H}\text{I}}\right)=g_{}\left(N_\text{H}\right)\frac{\text{d}N_\text{H}}{\text{d}N_{\text{H}\text{I}}}.$$
(1)
If $`g_{}`$ is a pure power-law:

$$g_{}\left(N_\text{H}\right)=KN_\text{H}^\alpha $$
(2)
$`f_{}\left(N_{\text{H}\text{I}}\right)`$ is also a power law when $`N_\text{H}=QN_{\text{H}\text{I}}^\beta `$:

$$f_{}\left(N_{\text{H}\text{I}}\right)=\beta KQ^{(\alpha 1)}N_{\text{H}\text{I}}^\xi ,$$
(3)
with $`\xi =1+\beta (\alpha 1)`$. For $`\beta 1`$ (an ionization fraction monotonically decreasing with $`N_\text{H}`$) and $`\alpha >1`$, $`f_{}\left(N_{\text{H}\text{I}}\right)`$ is always shallower than $`g_{}\left(N_\text{H}\right)`$.

For a generic $`N_\text{H}\left(N_{\text{H}\text{I}}\right)`$ relationship, $`f_{}(N_{\text{H}\text{I}})`$ is not tied to a power law. For example when there is a sharp break in the slope of $`N_\text{H}\left(N_{\text{H}\text{I}}\right)`$ at $`N_b`$, $`f_{}`$ is not continuous at $`N_b`$, being larger on the side where its slope is steeper:

$`f_{}\left(N_{\text{H}\text{I}}\right)=K^{}(`$ $`\left(\xi _11\right)\left(N_{\text{H}\text{I}}/N_b\right)^{\xi _1}\mathrm{\Theta }\left(N_bN_{\text{H}\text{I}}\right)`$ $`+`$ (4)
$`\left(\xi _21\right)\left(N_{\text{H}\text{I}}/N_b\right)^{\xi _2}\mathrm{\Theta }\left(N_{\text{H}\text{I}}N_b\right)`$ $`),`$
where $`\mathrm{\Theta }(x)`$ is the Heaviside step function.

Similarly we can consider multiple breaks: for example if $`N_{\text{H}\text{I}}=N_\text{H}`$ for large $`N_{\text{H}\text{I}}`$ and $`N_{\text{H}\text{I}}=ϵN_\text{H}`$ for small $`N_{\text{H}\text{I}}`$ (as for a constant ionization fraction), $`f_{}`$ follows a power law with index $`\alpha `$ in both regions, but the low-column density distribution is shifted by a factor $`ϵ^{\alpha 1}`$ with respect to the extrapolation of the high-column density one.

### 2.2 Orientation-averaged differential cross section

In this section we consider the effects of the random orientation of non-spherical absorbers on the observed column density distribution. We derive a simple analytic solution for “infinitely thin homogeneous slabs”, namely those whose thickness ($`2h`$) is much smaller than their radial extension ($`2R`$), and whose face-on gas column density is constant. We then test a simple approximation for treating “finite thickness slabs”, using the exact solution for homogeneous ellipsoids, given in Appendix A.

Let $`\sigma (N_{\text{H}\text{I}},\mu )`$ be the differential cross section of an absorber, whose symmetry axis is tilted by an angle $`\theta `$ with respect to the line of sight. Its differential cross section, averaged over orientations, is:

$$\widehat{\sigma }\left(N_{\text{H}\text{I}}\right)=\frac{1}{2}_1^1\sigma (N_{\text{H}\text{I}},\mu )𝑑\mu $$
(5)
(where we have defined $`\mu \mathrm{cos}\theta `$). For infinitely thin slabs of surface area $`A`$ and face-on column density $`N_{\text{H}\text{I}}`$, the averaged cross section reads:

$$\widehat{\sigma }\left(N_{\text{H}\text{I}}\right)=\frac{1}{2}_1^1A|\mu |\delta \left(N_{\text{H}\text{I}}N_{\text{H}\text{I}}/|\mu |\right)𝑑\mu =\frac{AN_{\text{H}\text{I}}^2}{N_{\text{H}\text{I}}^3}\mathrm{\Theta }\left(N_{\text{H}\text{I}}N_{\text{H}\text{I}}\right),$$
(6)
where $`\delta (x)`$ is the Dirac distribution. Therefore the averaged cross section vanishes for $`N_{\text{H}\text{I}}<N_{\text{H}\text{I}}`$, and is proportional to $`N_{\text{H}\text{I}}^3`$ for $`N_{\text{H}\text{I}}>N_{\text{H}\text{I}}`$.

For absorbers of finite thickness the exact behavior of $`\widehat{\sigma }`$ depends on their detailed shape, but the following approximate formula holds:

$$\widehat{\sigma }\left(N_{\text{H}\text{I}}\right)=\frac{AN_{\text{H}\text{I}}^2}{N_{\text{H}\text{I}}^3}\mathrm{\Theta }\left(N_{\text{H}\text{I}}N_{\text{H}\text{I}}\right)\mathrm{\Theta }\left(N_{\text{H}\text{I}}R/hN_{\text{H}\text{I}}\right).$$
(7)
It is shown in Appendix A that the accuracy of Eq. (7) decreases as $`h/R`$ grows larger. When $`hR`$, $`\widehat{\sigma }`$ for ellipsoidal absorbers differs noticeably from Eq. (7); but, once Eq. (7) is convolved with realistic $`f_{}`$ distributions, results keep close to the exact ones.

### 2.3 From the face-on to the line of sight $`N_{\text{H}\text{I}}`$ distribution

The line of sight $`N_{\text{H}\text{I}}`$ distribution, $`f(N_{\text{H}\text{I}})`$, can be written as:

$$f\left(N_{\text{H}\text{I}}\right)=\frac{\widehat{\sigma }\left(N_{\text{H}\text{I}}\right)}{A}f_{}\left(N_{\text{H}\text{I}}\right)𝑑N_{\text{H}\text{I}},$$
(8)
which, using Eq. (7), gives:

$$f\left(N_{\text{H}\text{I}}\right)=\frac{1}{N_{\text{H}\text{I}}^3}_{N_{\text{H}\text{I}}h/R}^{N_{\text{H}\text{I}}}N_{\text{H}\text{I}}^2f_{}\left(N_{\text{H}\text{I}}\right)𝑑N_{\text{H}\text{I}}.$$
(9)
This equation shows that $`f(N_{\text{H}\text{I}})`$ is sensitive to the values of $`f_{}`$ over the whole range $`[N_{\text{H}\text{I}}h/R,N_{\text{H}\text{I}}]`$ (i.e. $`[0,N_{\text{H}\text{I}}]`$ for infinitely thin slabs).

If there is a single break in $`f_{}`$, as given by Eq. (4), the corresponding $`f(N_{\text{H}\text{I}})`$ is:

$`f\left(N_{\text{H}\text{I}}\right)={\displaystyle \frac{K}{N_{\text{H}\text{I}}^3}}(`$ $`{\displaystyle \frac{\xi _11}{3\xi _1}}N_b^{\xi _1}(\mathrm{min}(N_{\text{H}\text{I}},N_b)^{3\xi _1}\mathrm{min}(N_{\text{H}\text{I}}h/R,N_b)^{3\xi _1})+`$ (10)
$`{\displaystyle \frac{\xi _21}{3\xi _2}}N_b^{\xi _2}(\mathrm{max}(N_{\text{H}\text{I}},N_b)^{3\xi _2}\mathrm{max}(N_{\text{H}\text{I}}h/R,N_b)^{3\xi _2})).`$
For $`N_{\text{H}\text{I}}<N_b`$ and for $`N_bR/h<N_{\text{H}\text{I}}`$, $`f`$ follows power laws with indices $`\xi _1`$ and $`\xi _2`$ respectively. In the intermediate $`N_{\text{H}\text{I}}`$ range the behavior of $`f`$ is more complex: if $`f_{}`$ is steeper at $`N_{\text{H}\text{I}}<N_b`$, its slope approaches $`3`$ for $`N_b<N_{\text{H}\text{I}}N_bR/h`$ while it bends to connect to the power-law regime at $`N_{\text{H}\text{I}}>N_bR/h`$; in the opposite case $`f`$ suddenly rises for $`N_{\text{H}\text{I}}>N_b`$, and then smoothly connects with the high column density power-law branch.

Fig. 1 shows qualitatively the features present in $`f(N_{\text{H}\text{I}})`$ for a realistic $`N_\text{H}(N_{\text{H}\text{I}})`$ relationship (see Paper I): $`N_\text{H}N_{\text{H}\text{I}}^\beta `$ for $`N_{b1}N_{\text{H}\text{I}}N_{b2}`$, $`N_\text{H}N_{\text{H}\text{I}}`$ elsewhere. We choose $`\alpha =1.5`$, $`h/R=0.1`$, $`N_{b1}=10^{17}\mathrm{cm}^2`$, $`N_{b2}=10^{20}\mathrm{cm}^2`$, and $`\beta =0.001`$ to emphasize the various regimes described above. Note that $`f(N_{\text{H}\text{I}})`$ at low column densities is lower by a factor $`\left(N_{b1}/N_{b2}\right)^{(1\beta )(\alpha 1)}`$ with respect to the extrapolation of the distribution at high column densities (dashed line). At intermediate column densities a dip in the distribution appears, approaching a $`3`$ slope towards the low column density side.

## 3 DLS&LLS surveys and our data sample

In this section we show which information should be extracted from the large number of DLS&LLS surveys available in order to build up a homogeneous database for a global statistical approach.

In a LLS the value of $`N_{\text{H}\text{I}}`$ can be determined from the ratio between $`I`$, the residual Lyman continuum flux on the blue side of the break, and $`I_0`$, the unabsorbed continuum flux:

$$N_{\text{H}\text{I}}=1.6\times 10^{17}\tau _{LL}\mathrm{cm}^2=1.6\times 10^{17}\mathrm{ln}\left(I_0/I\right)\mathrm{cm}^2.$$
(11)
If the residual flux $`I`$ is too small to be measured, only a lower limit to $`N_{\text{H}\text{I}}`$ can be derived: this usually happens for $`\tau _{LL}\stackrel{>}{}3`$. Since there are measurements sensitive to Lyman continuum optical depths as small as 0.4, we include in our database all LLS surveys which are sensitive to $`N_{\text{H}\text{I}}10^{16.81}\mathrm{cm}^2`$.

We include in the compilation all available searches for absorption lines with a rest frame equivalent width $`W5\text{Å}`$, for which $`N_{\text{H}\text{I}}`$ can be estimated from:

$$N_{\text{H}\text{I}}=1.88\times 10^{18}W^2\mathrm{cm}^2.$$
(12)
A complication comes from the fact that spurious lines with large equivalent width may result from a blending of weaker lines with metal line systems or Lyman-$`\alpha `$ forest lines. Only some surveys have sufficient spectral resolution to test, by fitting a Voigt profile, whether lines with $`W5\text{Å}`$ are really damped. In our database we include both the column density value determined directly from Eq. (12) and that derived from the Voigt profile fit, whenever this is available. By comparing these two values we find that there is a bias in the $`N_{\text{H}\text{I}}`$ values estimated from the equivalent width, and derive the following statistical correction:

$$\mathrm{lg}_{10}N_{\text{H}\text{I},corr}=6.261+0.705\mathrm{lg}_{10}N_{\text{H}\text{I},W}.$$
(13)
We apply this correction to all absorption lines for which only the $`W`$ value is available. If only a lower limit to $`\tau _{LL}`$ can be derived from the Lyman continuum absorption, searches for the corresponding damped Lyman-$`\alpha `$ line can be used to establish a value or an upper limit to the HI column density. If the corresponding Lyman-$`\alpha `$ absorption line has not been searched for, we take an upper limit of $`5\times 10^{21}\mathrm{cm}^2`$, which is the highest measured value of $`N_{\text{H}\text{I}}`$ in our database.

Unfortunately upper limits to saturated LLS leave large uncertainties in the column density range $`5\times 10^{17}\stackrel{<}{}N_{\text{H}\text{I}}\stackrel{<}{}5\times 10^{19}\mathrm{cm}^2`$, where theoretical predictions are very sensitive to the detailed shape of the $`N_\text{H}(N_{\text{H}\text{I}})`$ relationship. Hence it is essential to find a correct statistical approach for evaluating the $`N_{\text{H}\text{I}}`$ distribution in the presence of large uncertainties to determine the ionization conditions of the gas as well as the total gas column density distribution in LLS and DLS.

We have collected data from the following references: Tytler (1982); Bechtold et al. (1984); Wolfe et al. (1986); Tytler (1987); Lanzetta (1988); Sargent, Steidel and Boksenberg (1989); Turnshek et al. (1989); Lanzetta (1991); Lanzetta et al. (1991); Courvoisier and Paltani (1992); Bahcall et al. (1993); Storrie-Lombardi et al. (1994); Lanzetta et al. (1995); Stengler-Larrea et al. (1995); Wolfe et al. (1995); Storrie-Lombardi et al. (1996); Jannuzi et al. (1998); Storrie-Lombardi and Wolfe (2000). In order to ensure an accurate estimate of the coverage of the survey we have excluded from the compilation the full spectrum or part of it whenever the continuum was visually found to be very noisy or non-uniform (e.g. some spectrum edges where the flux is below the 2-$`\sigma `$ level in Lanzetta et al. (1995) or the objects 1130-106Y, PKS1206+459, MC1215+113 in Bahcall et al. (1993)). The final number of directions i.e. of QSOs included is 661. We have decided to split each direction in parts (hereafter called “paths”) such that each path is defined by homogeneous search parameters and contains at most one detection. In other words one direction may correspond to more than one path if different sensitivities are used, or if more than one absorber is present along the line of sight. Our database is available in electronic form upon request to the authors. Parameters contained in the data table can be divided into two groups: those describing the search, and those related to the detection.

The main parameters associated with the search coverage are: (i) the redshift path limits, i.e. the upper and lower redshift for one direction in space observed with a given sensitivity; (ii) the sensitivity, specified by the threshold above which $`N_{\text{H}\text{I}}`$ was detectable.

For each direction the upper redshift was set to $`5000\mathrm{km}\mathrm{s}^1`$ less than the QSO redshift, whenever the spectral coverage extends beyond it. We have usually taken as the lowest redshift that given by the observers. This does not coincide with the spectral coverage boundary in the presence of “shadowing effect”. The “shadowing effect” intervenes whenever the continuum flux, absorbed by a LLS at $`z_{}`$, gets below detectability: in this case the lower limit of the path has been set to $`z_{}`$. If no DLS searches have been performed no further paths are appended along that direction; otherwise, due to the different frequency of the Lyman-$`\alpha `$ line respect to the H ionization threshold, further paths are considered between $`z_{}`$ and the Lyman-$`\alpha `$ shadowing redshift ($`z=0.75z_{}0.25`$). The lowest redshift of the LLS search in the absence of shadowing is reported in the database as well, since in Section 5 it will be used for simulating data samples. The sensitivity for a Lyman-limit search depends on the minimum opacity detectable in the survey: typically $`\tau _{LL}=0.4`$, 1.0 or 1.5, corresponding to column density thresholds $`\mathrm{lg}_{10}N_{th}=`$16.81, 17.20 and 17.38 respectively. We have set to 0.4 the minimum detectable $`\tau _{LL}`$ for observations listed in Bahcall et al. (1993) and Tytler (1987); we have instead set to 1.0 the threshold for all the observations listed in Lanzetta (1991) and Sargent, Steidel and Boksenberg (1989) which were revised by Stengler-Larrea et al. (1995). For absorption line searches with rest frame equivalent width $`W5\text{Å}`$ the sensitivity of the path was set to $`\mathrm{lg}_{10}N_{\text{H}\text{I}}=20.13`$ (since we have applied the correction given in Eq. (13)). If both LLS and DLS searches were performed on the same path the threshold of the path was set to the lowest of the two, except in the presence of shadowing effects.

The parameters associated with the detections are: (i) the redshift of the detected absorber; (ii) the decimal logarithm of the estimated HI column density of the absorber along the line of sight ($`\mathrm{lg}_{10}N_{\text{H}\text{I}}`$); (iii) the minimum and the maximum allowed values for $`\mathrm{lg}_{10}N_{\text{H}\text{I}}`$; (iv) the Voigt correction factors.

When $`\mathrm{lg}_{10}N_{\text{H}\text{I}}`$ could not be determined from $`\tau _{LL}`$ or $`W`$ measures, but only limits for it are available, we have reported as $`\mathrm{lg}_{10}N_{\text{H}\text{I}}`$ the arithmetic mean between the minimum and maximum allowed values. When instead the value of $`\mathrm{lg}_{10}N_{\text{H}\text{I}}`$ could be directly determined we have estimated the minimum and maximum allowed values by assuming a standard spread. In the case of non-saturated LLS we have used a spread of $`\pm 0.2`$ in decimal logarithm, except for the cases listed in Jannuzi et al. (1998) for which we took the uncertainties quoted by the authors. In the case of DLS we have used again a spread of $`\pm 0.2`$ when the line profile has been resolved and the Voigt profile has been fitted; otherwise we have used $`\pm 0.4`$. The uncertainties in $`\mathrm{lg}_{10}N_{\text{H}\text{I}}`$ that we quote are usually larger than those that can be directly derived from uncertainties on either $`\tau _{LL}`$ or $`W`$. This is because there are other sources of errors, such as the determination of the level of the continuum (which affects both $`W`$ and $`\tau _{LL}`$ estimates), or the possible blending with other small absorption lines. The Voigt correction factors are the differences between the $`\mathrm{lg}_{10}N_{\text{H}\text{I}}`$ values derived from the Voigt-profile fit and those from $`W`$ (Eq. (12)).

Two further columns contain flags, which are different from 0 whenever LLS and DLS searches respectively have been performed on a given path. A value of $`1`$ implies that no systems were found, a value of 1 means that the column density estimate is quite accurate, 2 that it is less accurate and 3 that it is very uncertain, as for DLS declared “non damped”. In this last case the column density value for the DLS has been set to zero, or used as an upper limit if saturated LLS are detected at the same redshift.

The coverage of the whole sample, with and without shadowing, for various search sensitivities is shown in Fig. 2, where the inhomogeneity of the actual coverage is evident.

## 4 Statistical Analysis of the Data

Various authors have analyzed Lyman-$`\alpha `$ absorbers data to determine the $`N_{\text{H}\text{I}}`$ distribution and/or its evolution with time. Nevertheless no statistical analysis of the $`N_{\text{H}\text{I}}`$ distribution in LLS and DLS has been attempted so far using all available data and without a pre-defined coarse binning. In the following we show how to perform a global unbinned statistical analysis of heterogeneous data.

### 4.1 Maximum Likelihood Analysis

Our approach for estimating the $`N_{\text{H}\text{I}}`$ distribution is similar to the Maximum Likelihood analysis used by Storrie-Lombardi et al. (1996), with the following two improvements: (i) we use a generic (not just a pure power-law) distribution for $`N_{\text{H}\text{I}}`$; (ii) we take into account uncertainties in the $`N_{\text{H}\text{I}}`$ measurements. In order to simplify the notation we shall use the symbol $`N`$ instead of $`N_{\text{H}\text{I}}`$ in this section as well as in Section 4.2.

Let $`f(N,z)`$ be the distribution of column density $`N`$ at redshift $`z`$. The average number of detections along a given direction is $`f(N,z)\delta v`$, where $`\delta v=\delta N\delta z`$ is the area of an infinitesimal cell in the $`N`$$`z`$ plane. For infinitesimal cells the probability of having more than one detection in the same cell is negligible. The probabilities of having zero or one detection in a given cell are respectively:

$$P_0(N,z)=\mathrm{exp}\left(f\delta v\right);P_1(N,z)=(f\delta v)\mathrm{exp}\left(f\delta v\right).$$
(14)
The likelihood for a given data set is the product of the probabilities over “all” cells in the $`N`$$`z`$ plane and over all directions in space. This product can be separated into two factors, one relative to the “empty” cells, and the other to the “full” cells:

$$_o=\underset{i=\text{all}}{}P_i=\underset{i=\text{empty}}{}P_{0,i}\underset{j=\text{full}}{}P_{1,j}=\underset{i=\text{all}}{}\mathrm{exp}\left(f_i\delta v\right)\underset{j=\text{full}}{}\left(f_j\delta v\right).$$
(15)
We know the total number of full cells (i.e. the number of detections, hereafter $`p`$) but we might not be able to position the detections into specific cells in the $`N`$$`z`$ plane due to the large uncertainties of $`N`$. For instance, if a detection has $`N_j`$ in the range $`[N_{min,j},N_{max,j}]`$, the absorber is in one of the cells of this $`N`$ range, but we don’t know in which one: the probability for this event is the sum of probabilities relative to the cells (the same considerations would apply to errors in $`z`$ if redshifts are not well determined). In the limit of small cell sizes, the sum is then approximated by an integral of $`f`$ over the allowed range of $`N`$, and the logarithmic value of $`_o`$ reads:

$$\mathrm{ln}_o=w(N,z)f(N,z)𝑑N𝑑z+\underset{j=1}{\overset{p}{}}\mathrm{ln}\left(_{N_{min,j}}^{N_{max,j}}f(N,z_j)𝑑N\delta z\right),$$
(16)
where $`w(N,z)`$ is the number of different paths at a given $`z`$ with a column density threshold smaller than $`N`$. The first term in Eq. (16) is related only to the parameters of the search, and can be evaluated using the normalization condition for $`f`$:

$$w(N,z)f(N,z)𝑑N𝑑z=p.$$
(17)
We can rewrite Eq. (16) using the above normalization condition:

$$\mathrm{ln}_o=p\left(\mathrm{ln}(\delta z)1\right)+\underset{j=1}{\overset{p}{}}\mathrm{ln}\left(_{N_{min,j}}^{N_{max,j}}f(N,z_j)𝑑N\right).$$
(18)
If column densities have been measured with a rather small uncertainty for a number $`p_1`$ of absorbers, Eq. (18) can be approximated as:

$$\mathrm{ln}_o=p\left(\mathrm{ln}(\delta z)1\right)+\underset{j=1}{\overset{p_1}{}}\mathrm{ln}(\delta N_j)+\underset{j=1}{\overset{p_1}{}}\mathrm{ln}\left(f(N_j,z_j)\right)+\underset{j=p_1+1}{\overset{p}{}}\mathrm{ln}\left(_{N_{min,j}}^{N_{max,j}}f(N,z_j)𝑑N\right),$$
(19)
which allows a faster numerical computation. The terms with $`\delta N_j`$ and $`\delta z`$ are constant and can be neglected while searching for the maximum value of $`_o`$. The limiting case $`p_1=p`$ reduces Eq. (19) to the formula in Storrie-Lombardi et al. (1996).

### 4.2 The fit goodness and the true column density distribution

The Maximum Likelihood analysis does not contain any statistical test on the goodness of the “best” solution for $`f`$, since the absolute value of $`_o`$ is undetermined by a constant factor. Such a test can be done for example by comparing the theoretical versus the observational cumulative function. For large uncertainties of the column density values, this comparison requires a reasonable guess at the “true” value of $`N`$ for each individual detection. A “reasonable” guess does not necessarily mean that each $`N`$ value should reproduce the true one but that their overall distribution cannot be distinguished from the original one. We will show how to find this “reasonable” guess in the presence of large measurement uncertainties in $`N`$ in order to avoid the use of a coarse binning.

For a given absorber $`j`$ the probability of having a true column density smaller than $`N`$ in the allowed range $`[N_{min,j},N_{max,j}]`$ is:

$$R(N,j)=_{N_{min,j}}^Nf(N^{},z_j)𝑑N^{}/_{N_{min,j}}^{N_{max,j}}f(N^{},z_j)𝑑N^{}.$$
(20)
$`R_j`$ behave as random quantities and should have a uniform distribution in the interval $`[0,1]`$, independently of the shape of the distribution function $`f`$.

The probability of having a detection with column density smaller than $`N`$ in the overall $`N`$ range, is:

$$C(N,j)=_{N_{th,min}}^Nw(N^{},z_j)f(N^{},z_j)𝑑N^{}/_{N_{th,min}}^{\mathrm{}}w(N^{},z_j)f(N^{},z_j)𝑑N^{}.$$
(21)
The distribution of $`C_j`$ should also be uniform in the interval $`[0,1]`$, and is equivalent to the $`N`$ cumulative distribution for a narrow range of redshift.

For a sample of $`x_j`$ (sorted for increasing $`x_j`$ values, with $`0x_j1`$ and $`1jp`$) the cumulative distribution can be written as:

$$S(x)=\{\begin{array}{cc}0,\hfill & x<x_1\hfill \\ j/p,\hfill & x_jx<x_{j+1}\hfill \\ 1,\hfill & x_px.\hfill \end{array}$$
(22)
The conventional way to test whether a given a distribution is compatible with the uniform one is to use the Kolmogorov-Smirnov test on the maximum deviation between the two cumulative functions:

$$D=\underset{0x1}{\mathrm{max}}\left|S(x)x\right|;$$
(23)
the significance of $`D`$ being given by the function $`Q_{KS}(\sqrt{p}D)`$ (Kendall and Stuart, 1967).

However if the maximum discrepancy between the two distributions originates from some bias, the Kolmogorov-Smirnov test becomes completely insensitive to the overall match. Such effect is reduced if instead of $`D`$ we use the quantity $`U`$:

$$U^2=p_0^1\left(S(x)x_0^1\left(S(x^{})x^{}\right)𝑑x^{}\right)^2𝑑x.$$
(24)
$`U`$ varies with $`x_j`$ more smoothly than $`D`$ does, and its significance is $`Q_{KS}(\pi pU)`$ (Kendall and Stuart, 1967). We will use $`U`$ for testing the uniformity of $`R_j`$ and $`C_j`$ since it is practically more effective and robust than $`D`$ to ensure a good matching over the whole $`N`$ range.

For any given distribution $`f`$ we compute a “best guess” for the “true” column densities by maximizing the product of the significances relative to $`U_R`$ and to $`U_C`$ (hereafter referred to as $`_R`$ and $`_C`$, respectively). There are always infinite guesses which allow a uniform distribution for $`R_j`$. However $`R_j`$ and $`C_j`$ are mutually dependent and their relationship depends on $`f`$ (see Eqs. (20) and (21)). For an incorrect choice of $`f`$ or of the estimated ranges $`[N_{min,j},N_{max,j}]`$ it may then become impossible to derive uniform distributions for $`R_j`$ and $`C_j`$ at the same time.

### 4.3 Numerical results

Given an original distribution function for the face-on total column density of absorbers of the form:

$$g_{}(N_\text{H},z)=K(1+z)^\gamma N_\text{H}^\alpha $$
(25)
and a functional relationship between $`N_\text{H}`$ and $`N_{\text{H}\text{I}}`$, a routine derives the Maximum Likelihood values of $`\alpha `$ and $`\gamma `$ for an average absorber cross section $`\widehat{\sigma }`$ (Eq. (7)). $`K`$ is fixed by the normalization condition to the total number of detections (Eq. (17)); $`\gamma `$ accounts both for cosmological effects, and for the physical evolution of the absorbers. Here we assume no evolution in the power-law index $`\alpha `$ in a given $`z`$ interval and a pure power-law dependence for $`g_{}`$, but it would be straightforward to implement our routine if necessary to change these assumptions. Physically meaningful $`N_{\text{H}\text{I}}`$$`N_\text{H}`$ profiles in various redshift ranges have been computed by simulating the ionization structure of gaseous slabs (see Paper I). To identify each profile we use the parameter $`X`$ defined as:

$$X=\mathrm{lg}_{10}\left(N_\text{H}/N_{\text{H}\text{I}}\right)\text{for}N_{\text{H}\text{I}}=1.6\times 10^{17}\mathrm{cm}^2.$$
(26)
(the corresponding $`N_{\text{H}\text{I}}`$ depends on $`\widehat{\sigma }`$). An increasing value of $`X`$ corresponds to a neutral-to-ionized transition getting sharper and occurring at higher $`N_\text{H}`$. After determining the best $`\alpha `$ for fixed values of the parameters $`X`$, $`\gamma `$ and ($`h/R`$) we can determine the best values for these 3 parameters by comparing the relative Maximum Likelihood solutions. We use the Likelihood given in Eq. (19) with $`p_1`$ as the number of detections with an $`N_{\text{H}\text{I}}`$ uncertainty, $`\mathrm{lg}_{10}N_{max}\mathrm{lg}_{10}N_{min}\mathrm{\Delta }\left(\mathrm{lg}_{10}N_{\text{H}\text{I}}\right)`$. Values for $`R_j`$ and $`C_j`$ are computed for all the remaining $`(pp_1)`$ detections. Such a computation is not relevant for the Likelihood maximization itself, but it is useful in order to compare the theoretical and the observational cumulative. Since it is not always possible to obtain uniformly distributed $`R_j`$ and $`C_j`$ for any set of model parameters, the uniformity level of $`R_j`$ and $`C_j`$ distributions is a further test on the quality of the fit to $`f(N_{\text{H}\text{I}})`$ distribution. For this reason the routine maximizes the product of $`_R`$ and $`_C`$.

The main results of our statistical fit procedure to data in the redshift range \[1.75,3.25\], where most of the observations are, have been discussed in Paper I, together with the relative cosmological implications. The Maximum Likelihood solution, $`\alpha =2.70`$ and $`X=2.82`$, is obtained for $`\gamma =1.0`$ and $`h/R=0.2`$ but $`\alpha `$ and $`X`$ are only weakly sensitive to reasonable $`\gamma `$ and $`h/R`$ variations. We now investigate the possible dependence of $`\alpha `$ and of the cumulative function on the limiting value of $`\mathrm{\Delta }(\mathrm{lg}_{10}N_{\text{H}\text{I}})`$ for the treatment of errors, keeping $`X`$ fixed to 2.82.

We consider 4 different cases: (a) no error treatment, i.e. for all absorbers we use the estimated HI column density (either the measured value or the middle point in the allowed $`\mathrm{lg}_{10}N_{\text{H}\text{I}}`$ range); (b) and (c) error treatment has been applied to all the detections with $`\mathrm{\Delta }(\mathrm{lg}_{10}N_{\text{H}\text{I}})`$ equal to 0.9 and 0.5 respectively; (d) error treatment has been applied to all detections ($`\mathrm{\Delta }(\mathrm{lg}_{10}N_{\text{H}\text{I}})=0`$). Case (c) has been used for the results derived in Paper I. Fig. 3 shows the resulting value of $`\alpha `$ for the 4 cases and the comparison between the theoretical and observational cumulatives. Notice that, when the treatment of errors is operative, there is a good match between the two cumulatives. The “observational” cumulatives in these cases (dotted histograms, in Fig. 3) depend on the assumed distribution, because the $`N_{\text{H}\text{I}}`$ positions vary accordingly.

The logarithm of the Likelihood $`_o`$, normalized to its maximum value, is shown in Table 1 as we vary $`\alpha `$ for all the 4 cases. We can see that the maximum of the Likelihood is very well determined, without the presence of secondary maxima, and with 1-$`\sigma `$ uncertainties of the order of 0.1. In Table 1 we also give the significances $`_R`$ and $`_C`$ of the uniformity of the $`R_j`$ and $`C_j`$ distributions. $`R_j`$ and $`C_j`$ can be uniformly distributed for a very wide range of $`\alpha `$ when the number of detections with variable $`N_{\text{H}\text{I}}`$ is large; however only near the Maximum Likelihood value of $`\alpha `$ $`R_j`$ and $`C_j`$ are both uniformly distributed. The maximum of $`_C`$ is better defined than that of $`_R`$, and its position agrees more closely with that of the maximum value of $`_o`$. The uniformity of the $`C_j`$ distribution represents a valid test of the goodness of the Maximum Likelihood solution and of the fitting routine.

In Paper I we have shown that, in the $`\alpha `$-$`X`$ plane, the constant-probability contours are elongated, giving a global uncertainty on $`\alpha `$ of about $`\pm 0.3`$ when $`X`$ uncertainties are taken into account. Fig. 4 shows the probability contours in the $`\alpha `$-$`X`$ plane for 3 redshift intervals which together cover all the redshifts in our compilation, each containing the same number of detections ($`90`$). We find that at lower redshifts (Subsample A) the probability contours are shifted towards higher values of $`\alpha `$ and $`X`$; at higher redshifts (Subsample C) instead the allowed $`\alpha `$ values do not show significant changes, while $`X`$ decreases. However, since these variations are of the same order of parameter uncertainties, they could just be statistical fluctuations. Further insight will come out from the simulated data analysis in the next session.

An issue that will be discussed more in detail in a future paper is whether line blending is present in the data and could affect the above results. The probability of blending of Lyman-$`\alpha `$ lines is relatively low for randomly distributed absorbers, but clustering or blending with some metal lines might imply a somewhat higher probability and affect lines with $`5W<10\text{Å}`$, sometimes declared as non damped. To address this point we have performed an analysis of the data without considering low-resolution searches for lines with $`5W<10\text{Å}`$: the $`\alpha `$ and $`X`$ values that we find are consistent with those relative to the whole data sample. Therefore we tend to believe that our data sample closely reflects the real HI distribution.

## 5 Simulations

Possible biases might affect the data analysis, either related to our assumptions about the physical properties of the absorbers, or to our interpretation of the observations. In this section we show how to use simulated data for investigating the latter type of biases. These might derive from: (i) an uneven coverage of $`N_{\text{H}\text{I}}`$ and $`z`$ since the database contains results from surveys with different sensitivities; (ii) non uniform uncertainties in $`N_{\text{H}\text{I}}`$; (iii) incorrect estimates of detection thresholds or of $`N_{\text{H}\text{I}}`$ uncertainties. Simulated data can also be used to test the statistical procedure described in the previous section, and the confidence levels properties for the best fitting parameters.

### 5.1 The numerical procedure

Simulations supply, in a controlled way, a large number of samples to which we can apply our statistical analysis. We have devised a routine to simulate data in a highly realistic way: starting from an “a priori” distribution function for the absorbers total column densities, it generates a data set compatible with a given column density and redshift coverage. In the present work we use a coverage consistent with the observations listed in our database. The basic ingredients of the simulation are: 1. an original distribution function of $`N_\text{H}`$ for absorbers of a given geometrical shape; 2. an $`N_\text{H}`$$`N_{\text{H}\text{I}}`$ relationship for each redshift bin; 3. a $`N_{\text{H}\text{I}}`$$`z`$ coverage for the simulated survey. After the total face-on column density has been transformed into a line of sight HI value, we introduce measurement errors similar to those present in the actual survey compilation. These errors are simulated by deriving first the random quantities $`R_j`$, uniformly distributed, and then locating the uncertainties $`[N_{\text{H}\text{I},min,j},N_{\text{H}\text{I},max,j}]`$ such that the position of “true” values $`N_{\text{H}\text{I}}`$ corresponds to $`R_j`$. When $`(\mathrm{lg}_{10}N_{\text{H}\text{I},min,j}+\mathrm{lg}_{10}N_{\text{H}\text{I},max,j})/2`$, is above $`10^{17.68}\mathrm{cm}^2`$ the absorption is considered as saturated. In this case, if there has been a search for the corresponding damped Lyman-$`\alpha `$ but the line has not been found, the uncertainty range is set to $`[10^{17.68},10^{20.13}]\mathrm{cm}^2`$; if instead the damped Lyman-$`\alpha `$ line has not been searched for, the uncertainty range is set to $`[10^{17.68},10^{21.7}]\mathrm{cm}^2`$.

For saturated absorption we simulate also the shadowing effect. Shadowed absorbers are eliminated from the sample and the $`N_{\text{H}\text{I}}`$$`z`$ coverage of the survey is modified accordingly for LLS and DLS searches: this “after-shadowing” coverage is what should be compared with the real data coverage.

### 5.2 Simulation results

We first simulate samples of data consistent with the Maximum Likelihood solution $`\alpha =2.683`$, $`X=2.82`$, obtained for real data in the $`z`$ range \[1.75,3.25\] for $`p_1=0`$. Results of the statistical analysis on one of our simulated data samples are shown in Fig. 5, where the efficiency of our method in recovering the original distribution is clear. Fig. 5 is the analogous of Fig. 3 for real data, the only difference being that now in (a) the “true” values of $`N_{\text{H}\text{I}}`$ for simulated data are used, without any degradation by measurement processes and are compared directly with the model. The slight disagreement between the two curves in (a) is indicative of the magnitude of statistical fluctuations in the simulation.

$`_o`$, $`_R`$ and $`_C`$ for the simulated sample show a behavior very similar to that for real data (Table 1). An increase in the number of detections on which the treatment of errors is applied improves the match between the “measured” distributions and the original ones. Our algorithm is very efficient in recovering the original $`R_j`$ distribution as well but there is no way to recover the original $`R_j`$ values of individual detections. Detections with the same range of column densities are in all respects identical among themselves, and any permutation between their respective $`R_j`$ values cannot be traced (i.e. there is no correlation between the “true” and derived $`R_j`$ values).

To investigate if the selection of an error threshold affects the determination of $`\alpha `$ we have produced 5000 different simulated samples with $`\alpha =2.683`$. On each sample we apply the analysis with the 4 different error thresholds as in Fig. 5. The results, reported in Table 2, show that the best fit $`\alpha `$ values are only slightly smaller than the original ones, and the dispersion for the best $`\alpha `$ values found ($`\sigma (\alpha )`$) is consistent with the average estimated uncertainty ($`<\sigma >`$).

We can also use the simulated data to check the uncertainties in the best-fitting parameters ($`\alpha `$,$`X`$) by producing a confidence level map for the simulated data, for 3 $`z`$-selected subsamples, analogously to what has been done in Section 4.3 for the real data (Fig. 4). Fig. 6 shows for one simulated sample that the confidence levels present an elongation similar to the levels in Fig. 4. This means that in the real data there is no evidence of deviations of the gas distribution from our model assumptions. The elongation depends on a partial degeneracy of models in the $`\alpha `$-$`X`$ plane: $`\alpha `$ should increase as $`X`$ increases to give the same slope for the distribution of low HI column density LLS, a drawback that will be difficult to cure by a moderate increase of the numbers of known absorbers.

By looking at the three subsamples (Fig. 6) it is clear that statistical fluctuations are present mostly in the direction of the elongation, consistently with the estimated uncertainties. With the present sample coverage the real data are still consistent with no evolution of $`\alpha `$ and $`X`$, but some evolution cannot be excluded, especially for the high redshift sample. For this sample in fact the shift of the probability contours is not in the direction of elongation but mostly towards lower values of $`X`$. This would imply that either the ionizing flux is lower or that the absorbing gas has higher volume densities than at lower redshifts. Our analysis shows that in order to have a definite answer on the evolution of $`\alpha `$ and $`X`$ closer contour levels are needed. This can be achieved by determining observationally the column density of saturated LLS (through lines in other spectral regions or through the recover of the quasar flux at shorter wavelengths) or by increasing the number of known high $`N_{\text{H}\text{I}}`$ DLS.

## Appendix A Averaged cross section for homogeneous ellipsoids

In this Appendix we compute analytically the differential cross section for absorbers which can be modelled as axisymmetric ellipsoids of constant density. We use a cartesian coordinate system $`(x,y,z)`$ in which $`z`$ is parallel to the line of sight. The ellipsoid has a semi-minor axis $`h`$ and semi-major axes $`R`$ along the orthogonal directions. If its axis of symmetry is tilted along the $`y`$ axis by an angle $`\theta `$ with respect to the line of sight, the surface of the ellipsoid is:

$$\left(\frac{\mathrm{sin}^2\theta }{h^2}+\frac{\mathrm{cos}^2\theta }{R^2}\right)x^2+2\mathrm{sin}\theta \mathrm{cos}\theta \left(\frac{1}{h^2}\frac{1}{R^2}\right)xz+\left(\frac{\mathrm{cos}^2\theta }{h^2}+\frac{\mathrm{sin}^2\theta }{R^2}\right)z^2=1\frac{y^2}{R^2}.$$
(A1)
The spatial extent of the ellipsoid in the $`z`$ direction is $`2\zeta `$. Defining:

$$\xi ^2\mathrm{cos}^2\theta +\mathrm{sin}^2\theta \frac{h^2}{R^2},$$
(A2)
we can write:

$$\zeta ^2=\frac{h^2}{\xi ^2}\left(1\frac{x^2}{\xi ^2R^2}\frac{y^2}{R^2}\right).$$
(A3)
The column density along the line of sight is $`N_{\text{H}\text{I}}=2n\zeta `$. Since $`\zeta h/\xi `$ the $`N_{\text{H}\text{I}}`$ distribution for a fixed orientation presents an upper cut of $`N_{\text{H}\text{I}}/\xi `$ (where $`N_{\text{H}\text{I}}=2nh`$ is the maximum face-on column density).

The ellipsoid cross section for column densities larger than $`N_{\text{H}\text{I}}`$ is:

$$\mathrm{\Sigma }(N_{\text{H}\text{I}})=\pi R^2\xi \left(1\left(\frac{\xi N_{\text{H}\text{I}}}{N_{\text{H}\text{I}}}\right)^2\right)$$
(A4)
which gives the differential cross section for one fixed orientation in space:

$$\sigma (N_{\text{H}\text{I}},\mu )=\pi R^2\frac{\xi ^3N_{\text{H}\text{I}}}{2n^2h^2},$$
(A5)
($`\mu \mathrm{cos}\theta `$) and an orientation-averaged cross section:

$$\widehat{\sigma }\left(N_{\text{H}\text{I}}\right)=\pi R^2_0^1\frac{\xi ^3N_{\text{H}\text{I}}}{2n^2h^2}\mathrm{\Theta }\left(N_{\text{H}\text{I}}\xi N_{\text{H}\text{I}}\right)𝑑\mu .$$
(A6)
Introducing the variables:

$$w=\frac{R\xi }{h}=\sqrt{1+\left(\frac{R^2}{h^2}1\right)\mu ^2};W=\frac{R}{h}\mathrm{min}(1,\frac{N_{\text{H}\text{I}}}{N_{\text{H}\text{I}}}),$$
(A7)
the average differential cross section for randomly oriented absorbers is:

$`\widehat{\sigma }(N_{\text{H}\text{I}})`$ $`=`$ $`\pi R^2{\displaystyle \frac{2h^3N_{\text{H}\text{I}}}{N_{\text{H}\text{I}}^2R^3\sqrt{R^2/h^21}}}{\displaystyle _1^W}{\displaystyle \frac{w^4}{\sqrt{w^21}}}𝑑w`$ (A8)
$`=`$ $`\pi R^2{\displaystyle \frac{3h^4N_{\text{H}\text{I}}\left(W\left(1+{\displaystyle \frac{2}{3}}W^2\right)\sqrt{W^21}+\mathrm{lg}\left(W+\sqrt{W^21}\right)\right)}{N_{\text{H}\text{I}}^2R^4\sqrt{1h^2/R^2}}}.`$
For $`N_{\text{H}\text{I}}<2nh`$ the $`\widehat{\sigma }(N_{\text{H}\text{I}})`$ increases linearly with $`N_{\text{H}\text{I}}`$ while it vanishes for $`N_{\text{H}\text{I}}>2nR`$. In the intermediate range the behavior of $`\widehat{\sigma }`$ is more complex, but for $`N_{\text{H}\text{I}}2nR`$, the following approximation holds:

$$\widehat{\sigma }\left(N_{\text{H}\text{I}}\right)=\frac{AN_{\text{H}\text{I}}^2}{N_{\text{H}\text{I}}^3}\left(1+𝒪\left(\left(\frac{hN_{\text{H}\text{I}}}{RN_{\text{H}\text{I}}}\right)^2\right)\right),$$
(A9)
where $`A=\pi R^2/2\sqrt{1h^2/R^2}`$. Notice the similarity between Eq. (A9) and Eq. (6) for infinitely thin slabs. In Fig. 7 we compare the exact solution for the ellipsoids of various axial ratios $`h/R`$ to the approximate formula we use in our routines (Eq. (7)).

Although $`\widehat{\sigma }`$ depends on the detailed shape of the absorbers, the approximation by a $`N_{\text{H}\text{I}}^3`$ power law with an upper cutoff is appropriate for ellipsoids as well as for other geometrical shapes, especially for thickness $`h/R1`$. The average value of the projected column density is instead more geometry dependent. For ellipsoidal absorbers this is:

$`\widehat{N}_{\text{H}\text{I}}`$ $`=`$ $`{\displaystyle \frac{_0^{N_{\text{H}\text{I}}R/h}N_{\text{H}\text{I}}\widehat{\sigma }(N_{\text{H}\text{I}})𝑑N_{\text{H}\text{I}}}{_0^{N_{\text{H}\text{I}}R/h}\widehat{\sigma }(N_{\text{H}\text{I}})𝑑N_{\text{H}\text{I}}}}`$ (A10)
$`=`$ $`{\displaystyle \frac{4N_{\text{H}\text{I}}R\sqrt{R^2h^2}}{3\left(R\sqrt{R^2h^2}+h^2\mathrm{ln}\left(\left(R+\sqrt{R^2h^2}\right)/h\right)\right)}}.`$
For $`hR`$, $`\widehat{N}_{\text{H}\text{I}}`$ approaches $`4N_{\text{H}\text{I}}/3`$. This result differs from that for a slab, for which:

$$\widehat{N}_{\text{H}\text{I}}=2N_{\text{H}\text{I}}\left(1h/R\right)/\left(1h^2/R^2\right),$$
(A11)
whose limiting value is $`\widehat{N}_{\text{H}\text{I}}=2N_{\text{H}\text{I}}`$ for infinitesimal thickness.


format: text

text:
# 1 Introduction

## 1 Introduction

A central task of modern cosmology is to determine what the universe made of. A number of observations suggest that the bulk of the matter in the universe is not luminous . Direct searches for baryonic matter in the form of massive compact halo objects (MACHOs) cannot account for the matter that seems necessary to explain these observations .

An alternative explanation is that weakly interacting massive particles (WIMPs) exist copiously in the halo of our galaxy but only rarely interact with ordinary matter . Candidate WIMPs from particle theory include the axion and the lightest supersymmetric particle (LSP) in supersymmetric theories with R-parity conservation.

Numerous experiments have been set up in attempts to directly detect WIMPs . The two which are sensitive to the smallest spin independent cross sections are the CDMS Ge experiment  and the DAMA NaI experiment . Recently, DAMA reported the presence of a signal consistent with a WIMP at better than $`4\sigma `$. When interpreted as a standard WIMP with spin independent interactions, CDMS rules out nearly all of the DAMA $`3\sigma `$ preferred region at $`90\%`$ confidence and all of it at $`84\%`$ confidence. Attempts to reconcile these experiments using spin dependent interactions have been shown to be in gross conflict with indirect detection experiments and previous direct searches .

In this paper, we will show that a simple modification to the properties of the dark matter particle can change the kinematics of the scattering sufficiently to reconcile the two experiments. In particular, we explore the possibility of inelastic dark matter: relic particles that cannot scatter elastically off of nuclei. The outline of the paper is as follows: we begin by comparing the details of the two experiments and give a naive argument why inelastic dark matter can reconcile them. In section two we explicitly calculate the event rate at CDMS and DAMA taking into account the inelasticity of the scattering. In section three we use this calculation to study what differences can arise relative to the elastic case and to examine whether there are regions of parameter space that give a signal at DAMA but a null result at CDMS. In section four we discuss how inelastic dark matter could arise from a massive complex scalar split into two approximately degenerate real scalars, or from a Dirac fermion split into two approximately degenerate Majorana fermions. We also present a specific model, featuring a real component of the sneutrino as the dark matter, in which the mass splitting required to reconcile DAMA and CDMS arises naturally. In section five we discuss direct detection possibilities at future experiments.

### 1.1 CDMS and DAMA

If we are to understand the DAMA signal as evidence of dark matter, but simultaneously accept the null result of CDMS, we must reconsider some basic element of the WIMP hypothesis. Before we address such a modification, we should understand the differences between the DAMA and CDMS experiments.

The DAMA experiment utilizes a set of NaI crystals at the Gran Sasso National Laboratory of INFN to search for WIMPs. The basic premise of the experiment is that if WIMPs are present in the galaxy, as the galaxy rotates we feel a “wind” of WIMPs which will scatter elastically off of the target nuclei. As the Earth moves in its orbit about the sun, the flux and velocity distribution (as seen by a terrestrial observer) vary. Rather than attempt to directly discriminate signal events against background, the DAMA experiment seeks to measure this modulation. There are two basic controls to this experiment. First, the signal phase must coincide with the Earth’s motion in the solar system, which moves maximally with the galactic rotation on June 2, and maximally against on December 2. The second requirement is that the signal must lie dominantly in the lowest energy bins - a characteristic signal of WIMP scattering.

In contrast, CDMS uses a smaller Ge target, but has excellent background rejection capable of distinguishing nuclear recoils from electron scatterings for scattering energies greater than $`10\mathrm{k}\mathrm{e}\mathrm{V}`$. As a consequence their limits are comparable to those that would have been expected from a null DAMA result.

Exclusion plots are typically given in the $`m_\chi \sigma _n`$ plane, where $`m_\chi `$ is the mass of the candidate and $`\sigma _n`$ is the scattering cross section per nucleon. Implicit is the assumption that there are no great modifications in the scattering process between the two experiments.

However, if the dark matter cannot scatter elastically, then kinematical effects substantially distinguish the experiments. Consider two states, $`\chi _{}`$ and $`\chi _+`$, with $`\chi _+`$ only slightly heavier than $`\chi _{}`$, such that $`\chi _{}`$ can only scatter by transitioning to $`\chi _+`$. It is a simple kinematical constraint that $`\chi _{}`$ can only scatter inelastically off of a nucleus with mass $`m_N`$ if

$$\delta <\frac{\beta ^2m_\chi m_N}{2(m_\chi +m_N)},$$
(1)
where $`\delta `$ is the mass splitting between $`\chi _{}`$ and $`\chi _+`$. The possibility of evading direct detection by having a large enough splitting $`\delta `$ was pointed out in . Here we focus on the fact the constraint of equation (1) becomes increasingly severe as $`m_N`$ is decreased. Since iodine has an atomic number of 127, while germanium has an atomic number of 73, we have the prospect of a situation where particles will scatter at DAMA but not at CDMS. For $`\beta c220\mathrm{k}\mathrm{m}/\mathrm{s}`$ (a typical dark matter particle velocity), and $`m_\chi =100\mathrm{G}\mathrm{e}\mathrm{V}`$, the limits are $`11\mathrm{k}\mathrm{e}\mathrm{V}`$ for CDMS and $`15\mathrm{k}\mathrm{e}\mathrm{V}`$ for DAMA. If the mass splitting $`\delta `$ were $`13\mathrm{k}\mathrm{e}\mathrm{V}`$, such a particle would be visible to DAMA but not CDMS.

Of course, in the halo of the galaxy there is a distribution of velocities, so the calculation is not as simple as we have just illustrated. In the full calculation, we will find that the values of $`\delta `$ relevant for reconciling the experiments are somewhat larger than 15 keV, and that the window for $`\delta `$ has a size $`50100`$ keV rather than $`5`$ keV.

## 2 Direct Detection Rates

In this section we review the standard calculation of event rates at direct detection experiments . The differential rate per unit detector mass is given by

$$\frac{dR}{dE_R}=N_T\frac{\rho _\chi }{m_\chi }_{v_{min}}𝑑vvf(v)\frac{d\sigma }{dE_R}.$$
(2)
Here $`E_R`$ is the recoil energy of the target nucleus, $`N_T`$ is the number of target nuclei per unit mass, $`\rho _\chi `$ is the local density of dark matter particles of mass $`m_\chi `$, $`\frac{d\sigma }{dE_R}`$ is the differential cross section for relic-nucleus scattering, and $`v`$ and $`f(v)`$ are the relic speed and speed distribution function in the detector rest frame. We take $`\rho _\chi =.3`$ GeV/cm3.

Because we are interested in spin-independent scattering, the differential cross section may be written

$$\frac{d\sigma }{dE_R}=\frac{m_N}{2v^2}\frac{\sigma _n}{\mu _n^2}\frac{\left(f_pZ+f_n(AZ)\right)^2}{f_n^2}F^2(E_R),$$
(3)
where $`m_N`$ is the nucleus mass, $`\mu _n`$ is the reduced mass of the relic-nucleon system, $`f_n`$ and $`f_p`$ are the relative coupling strengths to neutrons and protons, and $`\sigma _n`$ is the relic-neutron cross section at zero momentum transfer, in the elastic $`(\delta =0)`$ limit. We use the Helm form factor

$$F^2(E_r)=\left(\frac{3j_1(qr_0)}{qr_0}\right)^2e^{s^2q^2},$$
(4)
with $`q=\sqrt{2m_NE_R}`$, $`s=1`$ fm, $`r_0=\sqrt{r^25s^2}`$, and $`r=1.2A^3`$.

We assume a standard Maxwell-Boltzman distribution for the relic velocities in the galactic rest frame, with a root-mean-squared velocity $`v_{rms}=\sqrt{\frac{3}{2}}v_0`$, where we take $`v_0=220`$ km/s to be the rotational speed of the local standard of rest (LSR). In our calculation we take the escape velocity to infinity for simplicity, when one really should take $`v_{esc}650`$ km/s. By doing so we overestimate the signal for large values of the mass splitting $`\delta `$. For a 100 GeV relic, this is a 10% effect at CDMS for $`\delta =100`$ keV and a factor of two effect for $`\delta =150`$ keV. Because iodine is heavier than germanium, the effect is far milder at DAMA, roughly 10% at $`\delta =150`$ keV.

The earth’s speed relative to the galactic rest frame is

$$v_e=v_{}+v_{orb}\mathrm{cos}\gamma \mathrm{cos}\left(\omega (tt_0)\right).$$
(5)
Here $`v_{}=v_0+12`$ km/s, $`v_{orb}=30`$ km/s, $`\omega =2\pi `$/year, $`t_0`$ June 2nd, and $`\mathrm{cos}\gamma =.51`$. Defining the dimensionless variables $`\eta =v_e/v_0`$ and $`x_{min}=v_{min}/v_0`$, performing the velocity integration in (2), and applying (3), one obtains

$$\frac{dR}{dE_R}=\frac{N_Tm_N\rho _\chi }{4v_0m_\chi }F^2(E_R)\frac{\sigma _n}{\mu _n^2}\frac{\left(f_pZ+f_n(AZ)\right)^2}{f_n^2}\left(\frac{\mathrm{erf}(x_{min}+\eta )\mathrm{erf}(x_{min}\eta )}{\eta }\right).$$
(6)
For the DAMA detector, one should take into account there being two species of target nuclei with different quenching factors.

Often one considers the case where $`f_n=f_p`$ (so that the rate is proportional to $`A^2`$), and presents results in the $`m_\chi `$ \- $`\sigma _n`$ plane. Below we will be particularly interested in models in which the scattering is dominated by vector interactions arising from $`Z`$ boson exchange, giving $`f_n/f_p=(14\mathrm{sin}^2\theta _W).08`$ (and yielding a rate that is instead nearly proportional to $`(AZ)^2`$). In all of our calculations we take this value for $`f_n/f_p`$.

The differential rate of equation (6) depends on the mass splitting parameter $`\delta `$ through $`x_{min}`$, which is given by

$$x_{min}=\frac{1}{v_0}\sqrt{\frac{1}{2m_NE_R}}\left(\frac{m_NE_R}{\mu }+\delta \right),$$
(7)
where $`\mu `$ is the reduced mass of the relic-nucleus system. A non-zero $`\delta `$ increases the minimum relic speed required to produce a given nuclear recoil energy. In the following section we explore potential consequences for direct detection signals arising due to this modification.

## 3 Signals at CDMS and DAMA

Before we study whether there are regions of parameter space that are consistent with both DAMA and CDMS, it is worthwhile to investigate the differences arising when compared with the elastic case. We have seen that for a given velocity of dark matter particle, it might be that only DAMA is be able to detect the particle, and not CDMS. Given the distribution of relic velocities, we can now determine what effect the inelasticity has on the full signal integrated over all velocities. The simplest quantity to consider is the level at which the signal111For our purposes here, we will consider the signal to be the events falling in the $`10\mathrm{k}\mathrm{e}\mathrm{V}100\mathrm{k}\mathrm{e}\mathrm{V}`$ bins for CDMS and $`2\mathrm{k}\mathrm{e}\mathrm{V}10\mathrm{k}\mathrm{e}\mathrm{V}`$ for DAMA. is suppressed when compared with the elastic case. We plot these suppressions for CDMS and DAMA in figure 1.

We can easily see that our basic intuition is borne out. The greater the splitting between $`\chi _{}`$ and $`\chi _+`$, the greater the suppression for CDMS compared to that of DAMA. Since the CDMS excluded region only just covers the DAMA preferred region, even a factor of a few can dramatically improve the consistency of the experiments.

However, the relative suppression is not the only relevant quantity because DAMA is not sensitive to the total flux, but rather to the modulation of the flux. Because of the inelasticity, DAMA only sees those particles on the high tail of the Maxwellian distribution. Consequently, a small modulation in the average velocity can lead to much higher modulation for a given signal when compared with the elastic case. This effect is demonstrated in figure 2. The combination of these two effects results in DAMA having significant regions of sensitivity that are inaccessible to the existing Ge experiments.

The DAMA signal can be decomposed into background, unmodulated signal and modulated signal as

$$\mu _k=b_k+S_{0,k}+S_{m,k}\mathrm{cos}(\omega t),$$
(8)
where $`k`$ indexes the energy bin of each piece of the total measured events $`\mu `$. Recently, the DAMA collaboration published its best fit values of $`S_{0,k}`$ and $`S_{m,k}`$ for the energy bins $`23\mathrm{k}\mathrm{e}\mathrm{V}`$, $`34\mathrm{k}\mathrm{e}\mathrm{V}`$, $`45\mathrm{k}\mathrm{e}\mathrm{V}`$ and $`56\mathrm{k}\mathrm{e}\mathrm{V}`$ . It is tempting to fit the inelastic scattering case to these values, but to do so would be misleading. These best fit values are derived assuming the energy spectrum and relative size of the modulated piece to be given by the known relations for an elastically scattering WIMP. We have already seen in figure 2 that the standard WIMP and inelastic WIMP cases can lead to very different predictions for the relative size of the modulated piece, so any fit to the published best-fit values would not be rigorous.

Another, potentially more significant reason that we cannot use the standard WIMP $`S_{m,k}`$ values comes from changes in the energy spectrum of the events. Because the scattering is inelastic, the total number of events may not rise exponentially at low energy. In other cases, the spectra will be nearly identical. As examples we compare in figure 3 the expected WIMP spectrum of the modulation signal to the spectrum in the inelastic WIMP scenario for two values of $`\delta `$. The potential differences revealed in figure 3 make it possible to fit only to the model independent data recently published. We will discuss the details of this fit shortly.

These spectrum differences carry over to germanium experiments. As we show in figure 4, the changes can again be significant, and can again alter the interpretation of the experimental data. For instance, in the elastic case one expects an exponential rise in the number of relic scattering events for lower energies. Were CDMS to see many events in the $`4060\mathrm{k}\mathrm{e}\mathrm{V}`$ bin, but essentially an absence of events below $`40\mathrm{k}\mathrm{e}\mathrm{V}`$, this would be inconsistent with an elastic dark matter signal, but not with an inelastic dark matter signal. Again, we do not perform a rigorous fit to the CDMS data as this would require an ability to accurately simulate the correlation of multiple scatterers with single scatterers, which we lack.

For the purposes of generating allowed regions, we will thus use the following limits: for DAMA, we will use the published model independent modulation in the $`26\mathrm{k}\mathrm{e}\mathrm{V}`$ bins of $`0.088\pm 0.02\mathrm{counts}/\mathrm{day}/\mathrm{kg}`$  and consider the three sigma region to be allowed. DAMA claims not to have modulation in the higher energy bins. Although the measured modulation for energies above $`6\mathrm{k}\mathrm{e}\mathrm{V}`$ is not published, we will take an upper limit of $`0.003\mathrm{counts}/\mathrm{day}/\mathrm{kg}`$, which we consider quite reasonable given the errors on the best fit values for the higher energy bins. For CDMS, we will require a predicted mean of fewer than six total events, consistent with the published limits .

DAMA has also reported null results arising from a pulse shape analysis (PSA) of a portion of their NaI data  and of data from an experiment with Xe $`(A=129)`$ . Using the pulse shape, they can discriminate signal from background, and place a limit on the total number of events. Both of these studies affect the elastic WIMP preferred region for DAMA. Extracting rate limits from tables and plots of  and , we find the Xe studies have the dominant impact on our allowed regions. For the Xe experiment we require the signal to be less than $`0.7`$ $`\mathrm{counts}/\mathrm{day}/\mathrm{kg}`$ for the $`1315\mathrm{k}\mathrm{e}\mathrm{V}`$ bin, $`0.25`$ $`\mathrm{counts}/\mathrm{day}/\mathrm{kg}`$ for the $`1520\mathrm{k}\mathrm{e}\mathrm{V}`$ bin, $`0.15`$ $`\mathrm{counts}/\mathrm{day}/\mathrm{kg}`$ for the $`2025\mathrm{k}\mathrm{e}\mathrm{V}`$ bin, and $`0.075`$ $`\mathrm{counts}/\mathrm{day}/\mathrm{kg}`$ for the $`2530\mathrm{k}\mathrm{e}\mathrm{V}`$ bin, consistent with published limits .

We show the allowed regions subject to these constraints for various values of $`m_\chi `$ in figure 5. As expected, there are broad regions that fit the DAMA data and which are not excluded by CDMS. It is important to note that our qualitative results are not very sensitive to the details of the criteria used to determine what signals are consistent with experiment. The general features of figure 5 remain essentially intact even if we are more conservative in our estimates of the allowed counts at CDMS, or of the accuracy of the measured modulation at DAMA.

As an explicit example, let us consider the point $`m_\chi =70\mathrm{G}\mathrm{e}\mathrm{V}`$, $`\delta =105\mathrm{k}\mathrm{e}\mathrm{V}`$ and $`\sigma _n=5\times 10^{40}\mathrm{cm}^2`$. Here the modulation is quite consistent with the DAMA best fit point, but CDMS has only an expected signal of $`0.5`$ events, and the Xe pulse shape analysis constraints are evaded. A comparison between the inelastic point and the DAMA best fit values for the elastic case is given in table 1.

### 3.1 Cosmological Uncertainties

Unlike the ordinary WIMP scenario, the only inelastic WIMPs that scatter in existing experiments are those on the high end of the Maxwell-Boltzman velocity distribution. As such, there is greater uncertainty in the precise values of $`\sigma _n`$ that fit the data than for an ordinary WIMP.

In particular, there is significant uncertainty in the dispersion velocity $`v_{rms}`$ and in the local halo velocity $`v_0`$. Although these uncertainties are $`O(10\%)`$, the effects can be amplified because of the presence of the exponential in the distribution. We have investigated these effects and found that the preferred cross sections can can shift by as much as a factor of three for $`m_\chi =100\mathrm{G}\mathrm{e}\mathrm{V}`$ and a factor of seven for $`m_\chi =50\mathrm{G}\mathrm{e}\mathrm{V}`$. Likewise the local density $`\rho _\chi `$ is uncertain to a factor of approximately two, and moreover, the presence of substructure in the halo of the galaxy can lead to amplifications of the local density relative to the average halo density by a factor of three or more .

Finally, we must restate that we have assumed a Maxwell-Boltzmann distribution, which arises in the isothermal sphere model of dark matter. Changes to the velocity profile of the dark matter can have significant effects on the modulation for standard WIMPs , and inelastic dark matter potentially is even more sensitive to these changes. Such uncertainties are difficult to quanitify and we do not discuss them further.

Altogether these uncertainties can amount to a change in the preferred $`\sigma _n`$ values, but we should emphasize that the sizes of the regions of parameter space that yield consistency between DAMA and CDMS do not change dramatically.

## 4 Models of Inelastic Dark Matter

Up to this point, we have considered inelastic dark matter as an interesting phenomenological possibility, but have not addressed how such particles might arise in a reasonable model. One possibility is that the relic particle is a real scalar, so that its vector coupling to nuclei is forbidden by Bose symmetry. Consider a complex scalar $`\varphi =\frac{1}{\sqrt{2}}(a+ib)`$ coupled to an Abelian gauge field $`A_\mu `$. Its vector interaction comes from

$$|D_\mu \varphi |^2gA_\mu (a^\mu bb^\mu a).$$
(9)
That is, the real scalars $`a`$ and $`b`$ couple to each other, but neither couples to itself.

These real scalars are degenerate if the only mass term is $`m^2|\varphi |^2`$, but introducing a small additional mass term $`\mathrm{\Delta }^2\varphi ^2+h.c.`$, splits this degeneracy222Of course, $`\mathrm{\Delta }`$ violates gauge invariance, and can only arise once the gauge symmetry of the theory has been broken.. If $`m`$ is roughly $`100`$ GeV, and we want a splitting $`100`$ keV, then we require $`\mathrm{\Delta }^2(100\mathrm{MeV})^2`$. In the model of section 4.1, which features a real component of a sneutrino as the dark matter, this scale for $`\mathrm{\Delta }`$ arises naturally.

Before discussing this model, we note that the inelastic dark matter could instead be fermionic. Consider a Dirac fermion $`\psi =\left(\begin{array}{cc}\eta & \overline{\xi }\end{array}\right)`$ that has vector and axial-vector couplings to quarks:

$$\overline{\psi }\gamma _\mu (g_V^{}+g_A^{}\gamma _5)\psi \overline{q}\gamma ^\mu (g_V+g_A\gamma _5)q.$$
(10)
Assuming for simplicity that the various $`g`$’s are of comparable size, the largest contribution to the low-energy scattering of $`\psi `$ off of nuclei will come from the vector-vector piece, which will yield an amplitude that scales roughly as the number of nucleons. The axial-axial piece yields a smaller spin-dependent contribution that lacks this enhancement, while the vector-axial pieces vanish in the extreme non-relativistic limit.

Now suppose that in addition to a Dirac mass $`100`$ GeV for $`\psi `$, the Lagrangian also contains a very small Majorana mass term $`\frac{\delta }{2}(\eta \eta +\overline{\eta }\overline{\eta })`$, with $`\delta 100`$ keV. Then the Majorana fermion mass eigenstates are

$`\chi _1{\displaystyle \frac{i}{\sqrt{2}}}(\eta \xi )m_1=m\delta `$ (11)
$`\chi _2{\displaystyle \frac{1}{\sqrt{2}}}(\eta +\xi )m_2=m+\delta .`$ (12)
The vector current essentially couples $`\chi _1`$ to $`\chi _2`$, with only a small additional piece $`\delta /m`$ coupling each mass eigenstate to itself:

$$\overline{\psi }\gamma _\mu \psi i(\overline{\chi }_1\overline{\sigma }_\mu \chi _2\overline{\chi }_2\overline{\sigma }_\mu \chi _1)+\frac{\delta }{2m}(\overline{\chi }_2\overline{\sigma }_\mu \chi _2\overline{\chi }_1\overline{\sigma }_\mu \chi _1).$$
(13)
Because $`\delta /m10^6`$, we ignore the second term, and find that the only way for $`\chi _1`$ to scatter coherently off of nuclei is to make a transition into the heavier $`\chi _2`$ state. This inelastic process will dominate relative to the elastic, spin dependent scattering provided that the coherence enhancement, which gives a factor $`A^2510^3`$ in the cross section, overcomes than the suppression due to the inelasticity. In this case, the rate can depend sensitively on the mass of the target nucleus, as desired.

### 4.1 Sneutrino Dark Matter

Interestingly enough, a suitable candidate for inelastic dark matter has already been discussed in the literature. In supersymmetric theories with lepton number violation, the LSP can be a real component of the sneutrino : a Lagrangian term $`\mathrm{\Delta }^2\stackrel{~}{\nu }\stackrel{~}{\nu }+h.c.`$ lifts the degeneracy between the sneutrino’s odd and even CP eigenstates $`\stackrel{~}{\nu }_{}`$ and $`\stackrel{~}{\nu }_+`$. This splitting prevents elastic scattering of the lightest state, $`\stackrel{~}{\nu }_{}`$, off of nuclei through $`Z`$ exchange333There are contributions that will induce an elastic scattering, for instance from Higgs exchange, but these are all small and can be ignored for our purposes here.. There is still the challenge of achieving a cosmologically interesting relic abundance - for an ordinary $`100`$ GeV sneutrino $`\mathrm{\Omega }_{\stackrel{~}{\nu }}`$ comes out too small. In , this problem was resolved by taking the mass splitting between $`\stackrel{~}{\nu }_{}`$ and $`\stackrel{~}{\nu }_+`$ to be large enough to prevent coannihilation via an s-channel $`Z`$ in the early universe, $`\delta >5`$ GeV, leading to a radiatively generated neutrino neutrino mass $`m_\nu >`$ 5 MeV. Different approaches were taken in the models of . These models feature standard model singlet scalars $`\stackrel{~}{n}`$ that are kept light by a global symmetry  in analogue to the Giudice-Masiero solution to the $`\mu `$ problem , or by a gauged B-L symmetry . The singlet states mix with ordinary sneutrinos through weak scale $`A`$ terms, so that the gauge interactions of the mass eigenstates are suppressed by mixing angles. This suppression allows for an interesting relic abundance even for values of $`\delta `$ too small to prevent coannihilation between $`\stackrel{~}{\nu }_{}`$ and $`\stackrel{~}{\nu }_+`$.

For concreteness we will specialize to the model of . The global symmetry that prevents a tree level mass for the singlet $`\stackrel{~}{n}`$ states is broken by the vev of a spurion $`X`$ that also breaks supersymmetry. We assume that the $`A`$ and $`F`$ components of $`X`$ both have intermediate scale vevs: $`A_X\sqrt{F_X}m_I\sqrt{vM_{Pl}}`$. The spurion couples to the neutrino and singlet superfields according to

$$\frac{1}{M_{Pl}}\left[XLNH_u\right]_F+\frac{1}{M_{Pl}}[X^{}NN\left(1+\frac{X^{}X}{M_{Pl}^2}+\mathrm{}\right)]_D+h.c.$$
(14)
The operators of (14) can be justified by ordinary R parity (under which $`N`$ is odd and $`X`$ is even), together with an R symmetry where $`N`$ has R charge 2/3, $`X`$ has charge 4/3, and $`L`$ and $`H_u`$ have R charge 0. As discussed in , at tree level (14) yields a neutrino mass matrix whose light eigenvalue is $`v^2/M_{Pl}`$. However, (14) also contains

$$A\stackrel{~}{l}\stackrel{~}{n}h_u\mathrm{\Delta }^2(\stackrel{~}{n}\stackrel{~}{n}+h.c.),$$
(15)
with $`A`$ roughly weak scale and $`\mathrm{\Delta }^2m_I^5/M_{Pl}^3`$. These interactions radiatively induce a Majorana mass for the left-handed neutrino

$$m_\nu \frac{g^2}{384\pi ^2}\frac{v^{3/2}}{M_{Pl}^{1/2}}$$
(16)
that is larger than that obtained from the tree-level seesaw, and moreover, in roughly the correct range for explaining the atmospheric neutrino anomaly.

For our present purpose, however, the impact of (15) on the scalar masses is what matters most. Neglecting the small lepton number violating mass parameter $`\mathrm{\Delta }`$, the sneutrino mass-squared matrix is

$$\left(\begin{array}{cc}\stackrel{~}{\nu }^{}& n\end{array}\right)\left(\begin{array}{cc}m_L^2& \frac{1}{\sqrt{2}}Av\mathrm{sin}\beta \\ \frac{1}{\sqrt{2}}Av\mathrm{sin}\beta & m_R^2\end{array}\right)\left(\begin{array}{c}\stackrel{~}{\nu }\\ n^{}\end{array}\right).$$
(17)
The $`A`$ term coupling induces a mixing between $`\stackrel{~}{\nu }`$ and $`\stackrel{~}{n}`$, yielding a lighter mass eigenstate

$$\stackrel{~}{\nu }_1=\stackrel{~}{\nu }\mathrm{sin}\theta +\stackrel{~}{n}^{}\mathrm{cos}\theta .$$
(18)
The coupling of $`\stackrel{~}{\nu }_1`$ to the $`Z`$ boson is thus suppressed by $`\mathrm{sin}^2\theta `$. The lepton number violating parameter $`\mathrm{\Delta }^2`$ lifts the degeneracy between the CP-even and CP-odd components of $`\stackrel{~}{\nu }_1`$, leading to a small mass difference

$$\delta 2\mathrm{cos}^2\theta \frac{\mathrm{\Delta }^2}{m_1}.$$
(19)
For this splitting to resolve the conflict between CDMS and DAMA, one needs $`\delta 50100`$ keV, roughly444The lifetime of $`\stackrel{~}{\nu }_+`$ is $`\tau \left(\frac{1/3}{\mathrm{sin}\theta }\right)^4\left(\frac{100\mathrm{keV}}{\delta }\right)^5(410^2\mathrm{years})`$, so for the mass splittings and mixing angles of interest, it is safe to assume that only $`\stackrel{~}{\nu }_{}`$ is present today. Photons can be produced in these decays, but the decays take place before recombination for the parameters of interest, and the photons are soft enough to render negligible the effect on the CMBR spectrum.. For a 100 GeV sneutrino, this implies $`[X^{}XX^{}]_Dm_I^5(310^{10}\mathrm{GeV})^5`$, corresponding to a reasonable value for the intermediate scale.

To explore the feasibility of this scenario, we apply the same criteria used in section 3 to establish consistency with CDMS and DAMA for $`\delta =50`$ and 100 keV, and display the allowed regions in the ($`m_{\stackrel{~}{\nu }}`$, $`\mathrm{sin}\theta `$) plane. Note that because the scattering off of nuclei is suppressed both by the inelasticity of the reaction and by a $`\mathrm{sin}^4\theta `$ factor, the ability to obtain a large enough signal at DAMA depends crucially on the fact that ordinary sneutrinos give a signal roughly three orders of magnitude above present bounds. We also calculate the relic abundance as a function of $`m_{\stackrel{~}{\nu }}`$ and $`\mathrm{sin}\theta `$ using standard methods. The results shown in figures 6a)-d) indicate that there are indeed regions of parameter space featuring interesting relic abundances and acceptable direct detection rates. In the early universe, the efficiency of annihilation processes that occur via s-channel Higgs exchange, such as $`\stackrel{~}{\nu }_{}\stackrel{~}{\nu }_{}b\overline{b},ZZ,W^+W^{}`$, is sensitive to the size of the trilinear scalar coupling $`A`$, leading to the dependence of the relic abundance on $`A`$ evident in the figures.

Just as one specific illustrative example, consider the parameters $`m_{\stackrel{~}{\nu }}=70`$ GeV, $`\delta =70`$ keV, and $`\mathrm{sin}^4\theta =1/70`$. This choice of $`\mathrm{sin}\theta `$ leads to an interesting relic abundance for a broad range of SUSY parameters. For this choice of $`m_{\stackrel{~}{\nu }}`$, $`\delta `$, and $`\mathrm{sin}\theta `$, we calculate a mean of less than 2 events at CDMS and satisfy the constraints from the Xe pulse shape analysis. Moreover, as shown in table 2,

the values of $`S_{m,k}`$ we obtain for DAMA in the 2-6 keV energy range are nearly identical to those we obtain in the elastic case using DAMA’s best fit point $`\sigma _n=7.210^6`$ pb and $`m_\chi =52`$ GeV.

### 4.2 Indirect detection

As dark matter passes through the sun, it can scatter off of nuclei and be captured in the sun’s potential well . After a significant amount of dark matter has been captured, it can annihilate into other particles. If muon neutrinos are produced, those that reach the Earth can produce high energy muons through charged-current interactions. A number of experiments have attempted to direct WIMP matter indirectly by looking for these upward-going muons, leading to a current limit on their flux of $`10^{14}`$ cm-2 s-1 .

Even within the model of section 4.1, the expected flux of upward-going muons is quite uncertain, for a number of reasons. First, the capture rate in the sun is sensitive to the parameters $`m_{\stackrel{~}{\nu }}`$, $`\mathrm{sin}\theta `$ and $`\delta `$. Second, if the captured sneutrinos annihilate directly into neutrinos, the flavor of the lightest sneutrino determines what flavor of neutrino is produced, and details of the neutrino masses and mixings impact the flavor of the neutrino detected at the Earth. Third, cosmological uncertainties mentioned in section 3.1 can change the preferred region of $`\sigma _n`$ and thus the capture rate. Finally, relatively minor extensions to the model of section 4.1 can also complicate matters. As a consequence of these various sources of uncertainty, indirect techniques do not rule out inelastic sneutrino dark matter. However, they do impose strong constraints, as broad regions of parameter space lead to signals above experimental bounds. Moreover, indirect detection experiments offer the strong possibility of detection if the bound on the muon flux improves considerably .

One might expect that the same inelasticity that suppresses the signal at CDMS should be even more effective in suppressing the capture rate by the sun, which is mainly composed of relatively light nuclei. In fact, this is typically not the case. Because particles passing through the sun are unusually energetic (the escape velocity at the surface of the sun is much larger than the average velocity of a halo particle), the inelasticity is less relevant in the sun than at direct detection experiments.

In what follows, we have followed  in calculating solar capture rates and the induced muon flux, but have modified the approach to approximate the suppression of the capture rate due to the inelasticity (see Appendix). This suppression depends on the mass of the nucleus. For example, for $`\delta =100`$ keV and $`m_{\stackrel{~}{\nu }}=100`$ GeV, we find a factor $`20`$ suppression for scattering off of oxygen in the sun, and a factor $`2`$ suppression for scattering off of iron.

We will separately consider first, the case in which the sneutrinos cannot annihilate into W’s, and second, the case in which they can. For $`m_{\stackrel{~}{\nu }}<m_W`$, sneutrinos in the sun typically annihilate dominantly to neutrinos via t-channel neutralino exchange. If we neglect cosmological uncertainties, we find that for values of $`m_{\stackrel{~}{\nu }}`$, $`\delta `$ and $`\mathrm{sin}\theta `$ that lead to interesting relic abundances and consistency with CDMS and DAMA, the flux of neutrinos produced is quite large. If these are all muon flavor, we would expect a flux of upward-going muons of at least $`6\times 10^{13}\mathrm{cm}^1\mathrm{s}^1`$, in conflict with experimental results. On the other hand, if the sneutrinos annihilate into electron neutrinos that do not oscillate into muon neutrinos, bounds from direct detection are evaded entirely. If this is the case, only direct detection experiments will be able to yield a positive signal.

Finally, we note that there are specific parameter choices for which the dominant annihilation of sneutrinos in the sun is through s-channel Higgs to $`b\overline{b}`$. We find that this allows the flux of upward-going muons to be as small as a factor of $`2`$ above current limits for parameters that yield an acceptable abundance and acceptable direct detection signals.

With cosmological uncertainties included, more scenarios are allowed. Relatively small variations in $`v_{rms}`$ can accomodate factors of two, such as if the dominant annihilation is into $`b\overline{b}`$. Direct annihilation into muon neutrinos would require a more specious conspiracy of errors. For instance, if the solar system were presently in an anomalously high density region of the galaxy arising from substructure, and if $`v_{rms}`$ were $`3\sigma `$ above the value we have used, experiments could accomodate as much as one-third of the neutrinos produced being muon flavored. This seems quite unlikely, but is, at least in principle, still allowed.

However, for these lighter sneutrinos, if indirect detection experiments improve by an order of magnitude, they will be able to probe almost all of the parameter space, even accounting for a broad class of cosmological uncertainties, and situations where there is annihilation to $`b\overline{b}`$.

For heavier sneutrinos ($`m_{\stackrel{~}{\nu }}>m_W`$), the dominant annihilation processes in the sun can easily be s-channel Higgs exchange to W’s and Z’s. In this case we find that it is possible to reduce the expected signal at direct-detection experiments to a factor of $`3`$ above current limits for parameters consistent with DAMA, CDMS, and $`\mathrm{\Omega }_{\stackrel{~}{\nu }}h^2.1`$. These heavier sneutrinos are less affected by cosmological uncertainties, but these uncertainties still make it impossible to rule out this scenario. Future improvements in indirect detection could rule out this region of parameter space, especially as the experimental signal is less sensitive to the flavor of the lightest sneutrino than for the case of $`m_{\stackrel{~}{\nu }}<m_W`$.

Of course, this discussion applies only to the model of section 4.1, and one can consider modifications to the model that suppress the indirect detection signal. The premise of the model is that light standard model singlets are natural. Given this, if we add to that model another standard model singlet $`\eta `$, with the same R-charge as the right handed neutrino $`N`$, but opposite R-parity, we expect a superpotential interaction $`\eta NN`$. Then through t-channel $`\eta `$ exchange, $`\stackrel{~}{\nu }`$’s can annihilate to right-handed neutrinos. If these decay dominantly into muons or electrons (rather than tau’s) and off-shell W’s, we find that it is possible to bring the flux of upward going muons induced by the W decay products down to current limits555In such a scenario, the relic abundance is modified, but it is still possible to have $`\mathrm{\Omega }_{\stackrel{~}{\nu }}h^20.1`$.. Future indirect experiments would still likely be able to see the decay products of these right handed neutrinos.

This is just one example of a modification to the model which diminishes the signal, and there may be others, but such uncertainty is difficult to quantify. While indirect experiments offer a good opportunity to test specific models and regions of parameter space, there is an excellent likelihood that upcoming direct detection experiments will be able to determine whether inelastic dark matter is the resolution of the conflict between DAMA and CDMS.

## 5 Future Experiments

In the inelastic dark matter scenario, the boundaries of the DAMA preferred region are not far from the current limits from CDMS. Planned experiments should be able to cover the existing DAMA region. Most important are planned improvements to germanium experiments, and the CRESST experiment, which will use the heavy element tungsten.

CDMS will soon be moving to the Soudan mine, and should be able to improve its limits by at least two orders of magnitude . The GENIUS Ge experiment  should go well below that, likely allowing both to test much of the preferred regions discussed in section 3.

There is a caveat in this statement: in generating the plots of figure 5, we neglected to include the effect of a finite galactic escape velocity. This was a harmless simplification for our purposes there, because the effects at DAMA due to the finite galactic escape velocity are relatively minor. The effects can be much larger at CDMS.

Recall that the requirement for scattering is

$$\delta <\frac{\beta ^2}{2}\frac{m_Nm_\chi }{m_N+m_\chi }.$$
(20)
This constraint is particularly stringent for light candidates. For instance, with $`m_\chi =50\mathrm{G}\mathrm{e}\mathrm{V}`$ and $`v_{esc}=650\mathrm{k}\mathrm{m}/\mathrm{s}`$ equation (20) tells us that only for $`\delta <122\mathrm{k}\mathrm{e}\mathrm{V}`$ can one hope to obtain any signal at all at a germanium detector (recall that the highest velocity of particles incident on the earth is $`v_{esc}+v_{}`$). Thus, the higher $`\delta `$ regions may not be testable at CDMS.

For heavier candidates, the finite galactic escape velocity is not especially important, even at CDMS. With a galactic escape velocity of $`650\mathrm{k}\mathrm{m}/\mathrm{s}`$, and $`m_\chi =100\mathrm{G}\mathrm{e}\mathrm{V}`$, the cutoff for $`\delta `$ is $`172\mathrm{k}\mathrm{e}\mathrm{V}`$. On the other hand, the galactic escape velocity is not a particularly well known quantity, and if instead we take $`v_{esc}=450\mathrm{k}\mathrm{m}/\mathrm{s}`$, the cutoff for $`\delta `$ is only $`102\mathrm{k}\mathrm{e}\mathrm{V}`$.

These uncertainties make the CRESST experiment , using tungsten, especially significant. Because tungsten ($`A=183`$) is heavier than iodine ($`A=127`$), given adequate exposure time, CRESST should cover the DAMA preferred region, irrespective of cosmological uncertainties.

A very real possibility is that both germanium and tungsten experiments will have signals, which, when interepreted as elastic scatterings, would be inconsistent with one another. The most striking possibility of all is a spectrum deformation at the germanium detectors, as discussed in section 3. If CDMS were to see an excess of events in the $`3070\mathrm{k}\mathrm{e}\mathrm{V}`$ region, but no excess below $`30\mathrm{k}\mathrm{e}\mathrm{V}`$, it would be a compelling signature of this scenario.

## 6 Conclusions

If in fact the majority of the matter of the universe is non-baryonic, the attempt to determine its nature is one of the most exciting endeavors of modern cosmology. Existing dark matter searches have already begun to probe interesting regions of parameter space for candidate particles such as neutralinos and axions.

The positive result from the DAMA experiment is difficult to a understand in terms of these candidates, as it is in seeming conflict with constraints arising from the CDMS experiment. We have seen that this conflict vanishes if we allow for the possibility that the dark matter particle can only scatter inelastically.

We have shown that the sneutrino, when mixed with a singlet scalar with weak lepton number violation, is a viable candidate for inelastic dark matter. The regions of parameter space which give an interesting relic sneutrino abundance overlap with the regions which give a positive DAMA signal. Indirect detection experiments tightly constrain models of sneutrino dark matter, but do not rule them out.

Even absent a particular model, we find it interesting that such a simple modification of the dark matter’s properties can give remarkably different predictions, including the suppression of a signal at CDMS. We consider these results sufficiently interesting as to warrant an analysis of the full DAMA data set should the raw data become available.

Acknowledgments

We thank L. Hall, W. Haxton, C. Hogan, and T. Quinn for useful discussions, and H. Murayama and A. Nelson for reviewing the paper and providing valuable comments.

## Appendix A Appendix

Here we describe how we approximate the suppression of the rate of relic capture by the sun due to the inelasticity of the scattering. Ignoring nuclear form factors, the scattering probability for a given relative velocity $`w`$ is equally distributed between the minimum and maximum nuclear recoil energies $`\mathrm{\Delta }E_{min}`$ and $`\mathrm{\Delta }E_{max}`$ (these parameters depend on the nucleus mass, the relic mass $`m`$, $`w`$, and $`\delta `$). Ordinarily, the low-energy scattering cross section is independent of $`w`$, but in the inelastic case there is an additional phase space factor $`\sqrt{12\delta /(\mu w^2)}`$, where $`\mu `$ is the reduced mass. Capture only occurs when $`\mathrm{\Delta }E>\mathrm{\Delta }E_{capture}1/2\left(mw^2(m+\delta )v_{esc}^2(r)\right)\delta `$ holds. Here $`v_{esc}(r)`$ is the (position-dependent) escape velocity, which we approximate as

$$v_{esc}(r)=v_c^2\frac{M(r)}{M_{}}(v_c^2v_s^2),$$
(21)
where $`v_c=1354`$ km/s, $`v_s=795`$ km/s, and $`M(r)`$ is the mass contained within the radius $`r`$. The capture rate off of a given species of nuclei is then proportional to

$$_0^R_{}𝑑rr^2\rho (r)_{v_{esc}}^{\mathrm{}}𝑑ww^3e^{\frac{(w^2v_{esc}^2)}{v_0^2}}\sqrt{12\delta /(\mu w^2)}\left(\frac{\mathrm{\Delta }E_{max}\mathrm{\Delta }E_{capture}}{\mathrm{\Delta }E_{max}\mathrm{\Delta }E_{min}}\right),$$
(22)
where $`\rho (r)`$ is the mass density of the species and $`v_0`$ is the rotational speed of the local standard of rest. We calculate this factor (which does not account for form factor suppressions) in the elastic ($`\delta =0`$) and inelastic cases to estimate the suppression coming from the inelasticity. We then obtain capture rates by multiplying this suppression with the rate obtained for the elastic case using the formulae of  (which do include form factor suppressions).