format: text

text:
# Constraining the helium abundance with CMB

## I Introduction

Our understanding of the baryon abundance has increased dramatically over the last few years. This improvement comes from two independent paths, namely big bang nucleosynthesis (BBN) and cosmic microwave background radiation (CMB). Absorption features from high-redshift quasars allow to measure precisely the deuterium abundance, D/H. Combined with BBN calculations, this provides a reliable estimate of the baryon to photon ratio, $`\eta `$. An independent determination of the baryon content of the universe from CMB anisotropies comes from the increasingly precise measurements of the acoustic peaks, which bear a characteristic signature of the photon-baryon fluid oscillations. The agreement between these two completely different approaches is both remarkable and impressive (see details below). The time is therefore ripe to proceed and test the agreement between other light elements which are also probed both with BBN and CMB.

The helium abundance has been measured for many years from astrophysical systems. However, the error-bars are seemingly dominated by systematic errors which are hard to assess. Fortunately, the dependence of the helium mass fraction on the CMB anisotropies provides an independent way to measure $`Y_p`$. The aim of this work is to present the first determination of the helium abundance from CMB alone, and to clarify the future potential of this method. The latest CMB data are precise enough to allow taking this further step, and in view of the emerging “baryon tension” between BBN predictions from observations of different light elements olive03  possibly requires taking such a step. The advantage of using CMB anisotropies rather than the traditional astrophysical measurements, is that CMB provide a clear measurement of the primordial helium fraction before it could be changed by any astrophysical process. On the other hand the dependence of the CMB power spectrum on $`Y_p`$ is rather mild, a fact which makes it presently safe to fix the value of the helium mass fraction with zero uncertainity for the purpose of CMB data analysis of other cosmological parameters.

The paper is organized as follows. In section II we review the standard Big Bang Nucleosynthesis scenario. Section III discusses the role of the helium mass fraction for cosmic microwave background anisotropies, the methods used and results. We discuss our forecast for future CMB observations in section III.4, and offer our conclusions in section IV

## II Big Bang Nucleosynthesis

### II.1 The standard scenario

The standard model of big bang nucleosynthesis (SBBN) has only one free parameter, namely the baryon to photon ratio $`\eta _{10}=n_b/n_\gamma 10^{10}`$, which for long has been known to be in the range $`110`$ kolbturner90 . Thus by observing just one primordial light element one can predict the abundances of all the other light elements.

The deuterium to hydrogen abundance, D/H, is observed by Ly-$`\alpha `$ features in several quasar absorption systems at high red-shift, $`D/H=2.78_{0.38}^{+0.44}\times 10^5`$ kirkman2003 , which in SBBN translates into the baryon abundance, $`\eta _{10}=5.9\pm 0.5`$. Using SBBN one thus predicts the helium mass fraction to be in the range $`0.2470<Y_p<0.2487`$. The dispersion in various deuterium observations is, however, still rather large, ranging from $`D/H=1.65\pm 0.35\times 10^5`$ pettini2001  to $`D/H=3.98_{0.67}^{+0.59}\times 10^5`$ kirkman2003 , which most probably indicates underestimated systematic errors.

The observed helium mass fraction comes from the study of extragalactic HII regions in blue compact galaxies. One careful study IzotovThuan  gives the value $`Y_P=0.244\pm 0.002`$; however, also here there is a large scatter in the various observed values, ranging from $`Y_p=0.230\pm 0.003`$ oliveSS1997  over $`Y_p=0.2384\pm 0.0025`$ peimbert2002  and $`Y_p=0.2391\pm 0.0020`$ luri03  to $`Y_p=0.2452\pm 0.0015`$ IT1999 . Besides the large scatter there is also the problem that the helium mass fraction predicted from observation of deuterium combined with SBBN, $`0.2470<Y_p<0.2487`$, is larger than (and seems almost in disagreement with) most of the observed helium abundances, which probably points towards underestimated systematic errors, rather than the need for new physics olive03 ; Barger:2003zg . Figure 1 is a compilation of the above measurements, and offers a direct comparison with the current (large) errors from CMB observations (presented in section III below) and with the potential of future CMB measurements (discussed in section III D).

The observed abundance of primordial $`{}_{}{}^{7}\text{Li}`$ using the Spite plateau is possibly spoiled by various systematic effects Ryan ; Salaris . Therefore it is more appropriate to use the SBBN predictions together with observations to estimate the depletion factor $`f_7=^7\text{Li}_{\text{obs}}/^7\text{Li}_{\text{prim}}`$ instead of using $`{}_{}{}^{7}\text{Li}_{\text{obs}}^{}`$ to infer the value of $`\eta `$ burlesNT01 ; hansen02 ; steig99 ; steig02 .

The numerical predictions of standard BBN (as well as various non-standard scenarios) have reached a high level of accuracy hata95 ; burlesNT01 ; lopezturner1999 ; Emmp1 ; Emmp2 ; Cuoco:2003cu , and the precision of these codes is well beyond the systematic errors discussed above.

### II.2 The role of neutrinos

If the CMB-determined helium mass fraction turns out to be as high as suggested by SBBN calculations together with the CMB observation of $`\mathrm{\Omega }_bh^2`$ (as discussed above), this could indicate a systematic error in the present direct astrophysical helium observations. Alternatively, if the CMB could independently determine the helium value with sufficient precision to confirm the present helium observations, then this would be a smoking gun for new physics. In fact, one could easily imagine non-standard BBN scenarios which would agree with present observations of $`\eta _{10}`$, while having a low helium mass fraction. All that is needed is additional non-equilibrium electron neutrinos produced at the time of neutrino decoupling which would alter the $`np`$ reaction. This could alter the resulting helium mass fraction while leaving the deuterium abundance unchanged. One such possibility would be a heavy sterile neutrino whose decay products include $`\nu _e`$. A sterile neutrino with life-time of $`15`$ sec and with decay channel $`\nu _s\nu _e+\varphi `$ with $`\varphi `$ a light scalar (like a majoron), would leave the deuterium abundance roughly untouched, but can change the helium mass fraction between $`\mathrm{\Delta }Y_p=0.025`$ and $`\mathrm{\Delta }Y_p=0.015`$ if the sterile neutrino mass is in the range $`120`$ MeV dolgov98 . A simpler model would be standard neutrino oscillation between a sterile neutrino and the electron neutrino. The lifetime is about 1 sec when the sterile state has mass about 10 MeV, and the decay channel is $`\nu _s\nu _e+l+\overline{l}`$ (with $`l`$ any light lepton), and such masses and life-times are still unconstrained for large mixing angle dolgov00  (related BBN issues are discussed in refs. dibari ; kirilova ; abaz ; dolvil ). Such possibilities are hard to constrain without an independent measurement of the helium mass fraction.

Another much studied effect of neutrinos is the increased expansion rate of the universe if additional degrees of freedom are present (for BBN), and the degeneracy between the total density in matter and relativistic particles (for CMB). This issue has recently been studied in detail in refs. recentnu ; Barger:2003zg  in view of the new WMAP results, and we need not discuss this further here. We thus fix $`N_\nu =3.04`$ 3neutrinos . Also an electron neutrino chemical potential could potentially alter the BBN predictions kangsteigman , however, with the observed neutrino oscillation parameters the different neutrino chemical potentials would equilibrate before the onset of BBN nuosc , hence virtually excluding this possibility (see however dannygary ).

## III Cosmic Microwave Background

### III.1 Photon recombination and reionization

The recent WMAP data allow one to determine with very high precision the epoch of photon decoupling, $`z_{\mathrm{dec}}`$, i.e. the epoch at which the ionized electron fraction, $`x_e(z)=n_e/n_H`$, has dropped from 1 to its residual value of order $`10^4`$. Here $`n_e`$ denotes the number density of free electrons, while $`n_H`$ is the total number density of H atoms (both ionized and recombined). After this moment, photons are no longer coupled to electrons (last scattering), and they free stream. The redshift of decoupling has been determined to be $`z_{\mathrm{dec}}=1088_2^{+1}`$ spergel03 , which corresponds to a temperature of about 0.25 eV. Helium recombines earlier than hydrogen, roughly in two steps: around redshift $`z=6000`$ HeIII recombines to HeII, while HeII to HeI recombination begins around $`z<2500`$ and finishes just after the start of H recombination (see e.g. lyubarsky83 ; hu95 ; seager99 ; seager00 ).

Denoting by $`n_{He}`$ and $`n_b`$ the number densities per m3 of He atoms and baryons, respectively, the helium mass fraction is defined as $`Y_p=4n_{He}/n_b`$. The baryon number density is related to the baryon energy density today, $`\omega _b`$, by $`n_b=11.3(1+z)^3\omega _b`$ and we have $`n_H=n_b(1Y_p)`$. Usually, the ionization history is described in terms of $`x_e(z)=n_e/(n_b(1Y_p))`$. However, for the purpose of discussing the role of $`Y_p`$, it is more convenient to consider the quantity $`f_e(z)=n_e/n_b`$ instead, the ratio of free electrons to the total number of baryons. For brevity, we will call $`f_e`$ the free electron fraction. Once the baryon number density has been set by fixing $`\omega _b`$, one can think of $`Y_p`$ as an additional parameter which controls the number of free electrons available in the tight coupling regime. The CMB power spectrum depends on the full detailed evolution of the free electron fraction, but we can qualitatively describe the role of helium in four different phases of the ionization/recombination history (see Fig. 2).

* Before HeIII recombination all electrons are free, therefore $`f_e(z>6000)=1Y_p/2`$.
* HeII progressively recombines and just before H recombination begins, $`f_e`$ has dropped to the value $`f_e(z1100)=1Y_p`$.
* After decoupling, a residual fraction of free electrons freezes out, giving $`f_e(30<z<\mathrm{\hspace{0.33em}800})=f_e^{\text{res}}2.710^5\sqrt{\omega _m}/\omega _b`$.
* Reionization of all the H atoms gives $`f_e(z<\mathrm{\hspace{0.33em}20})=1Y_p`$.

During phase (a), the photon-baryons fluid is in the tight coupling regime. However the presence of ionized He increases diffusion damping, therefore having an impact on the damping scale in the acoustic peaks region. When the detailed energy levels structure of HeII is taken into account seager99 , the transition to phase (b) is smoother than in the Saha equation approximation. Therefore the plateau with $`f_e=1Y_p`$ is not visible in Fig. 2. Before H recombination, He atoms remain tightly coupled to H atoms through collisions, with the same dynamical behaviour. In particular, it is the total $`\omega _b`$ which determines the amount of gravitational pressure on the photon-baryons fluid, and which sets the acoustic peak enhancement/suppression. Hence we do not expect the value of $`Y_p`$ to have any influence on the boosting (suppression) of odd (even) peaks. The redshift of decoupling (transition between (b) and (c)) depends mildly on $`Y_p`$ in a correlated way with $`\omega _b`$, since the number density of free electrons in the tight coupling regime (just before H recombination) scales as $`n_e=f_en_b=n_b(1Y_p)`$. Hence an increase in $`\omega _b`$ can be compensated by allowing for a larger helium fraction. An analytical estimate along the same lines as in e.g. kolbturner90  indicates that a $`10\%`$ change in $`Y_p`$ affects $`z_{\text{dec}}`$ by roughly $`0.1\%`$, which corresponds to $`\mathrm{\Delta }z_{\text{dec}}1`$. This is of the same order as the current 1-$`\sigma `$ errors on $`z_{\text{dec}}`$, obtained by fixing $`Y_p=0.24`$.

After H recombination, the residual ionized electron fraction $`f_e^{\text{res}}`$ does not depend on $`Y_p`$, but is inversely proportional to the total baryon density (phase (c)). As the CMB photons propagate, they are occasionally rescattered by the residual free electrons. The corresponding optical depth, $`\tau ^{\text{res}}`$ is given by

$`\tau ^{\text{res}}`$ $`=`$ $`{\displaystyle _{t_0}^{t_{\text{dec}}}}n_e^{\text{res}}c\sigma _T𝑑t`$
$``$ $`1.8610^6{\displaystyle _0^{z_{\text{dec}}}}{\displaystyle \frac{(1+z)^2}{((1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }/\mathrm{\Omega }_m)^{1/2}}}𝑑z.`$
Performing the integral we can safely neglect the contribution of the cosmological constant at small redshift, since $`z_{\text{dec}}\mathrm{\Omega }_\mathrm{\Lambda }/\mathrm{\Omega }_m`$. Retaining only the leading term, the approximated optical depth from the residual ionization fraction is estimated to be

$$\tau ^{\text{res}}1.2410^6(1+z_{\text{dec}})^{3/2}0.045,$$
(2)
independent of the cosmological parameters and of the helium fraction. Therefore after last scattering we do not expect any significant effect on CMB anisotropies coming from the primordial helium fraction, until the reionization epoch.

Fairly little is known about the exact reionization mechanism and its redshift dependence (for a review see e.g. haiman ). Observation of Gunn-Peterson troughs indicate that the universe was completely ionized after redshift $`z6`$, when the universe seemingly completed the reionization BeckerFan , possibly for the second time cen . If temperature information only is available, CMB anisotropies are sensitive only to the integrated reionized fraction, represented by the optical depth, independent of the specific reionization history. However, specific signatures are imprinted on the E-polarization and ET-cross correlation power spectra by the detailed shape of the reionization history (for a detailed discussion, see Bruscoli02 ; Hu03 ; K03 ; Holder03 ). There are several physically motivated reionization scenarios, which however cannot be clearly distinguished at present Haiman03 ; Hansen:2003yj . In this work we use the most simple model, the sudden reionization scenario: we assume that at the reionization redshift $`z_r`$ all the hydrogen was quickly reionized, thus producing a sharp rise of $`n_e`$ from its residual value to $`n_H`$. More precisely, $`z_r`$ is the redshift at which $`x_e(z_r)=0.5`$. In our treatment we neglect HeII reionization, for which there is evidence at a redshift $`z3`$ (see frieman02  and references therein). This effect is small, since the extra electron released at $`z3`$ would change the reionization optical depth by about only $`1\%`$. We also neglect the increase of the helium fraction due to non-primordial helium production, which has a negligible effect on CMB anisotropies. Those approximations do not affect the results at today’s level of sensitivity of CMB data: for WMAP noise levels, even inclusion of the polarization spectra is not enough to distinguish between a sudden reionization scenario and a more complex reionization history. At the level of Planck a more refined modelling of the reionization mechanism will be necessary Holder03 ; naselsky03 .

In the sudden reionization scenario adopted here, the relation between reionization redshift and reionization optical depth, $`\tau _\text{r}`$, is given by

$`\tau _\text{r}`$ $`=`$ $`{\displaystyle _{t_0}^{t_{\text{reion}}}}n_ec\sigma _T𝑑t`$ (3)
$``$ $`11.3c\sigma _T\omega _b(1Y_p){\displaystyle _0^{z_r}}{\displaystyle \frac{d\eta }{da}}𝑑z,`$
where $`t`$ is physical time, $`\eta `$ is conformal time and $`a`$ the scale factor. Here again, since the number density of reionized electrons scales as $`\omega _b(1Y_p)`$, the redshift of reionization is positively correlated with $`Y_p`$ (for fixed optical depth and baryon density).

As a result of the physical mechanism described above, a $`10\%`$ change in $`Y_p`$ has a net impact on the CMB power spectrum at the percent level. The impact on the CMB temperature and polarization power spectra is highlighted in Fig. 3.

In the temperature panel, we notice that a larger helium fraction slightly suppresses the peaks because of diffusion damping, while it has no impact on large scales. Polarization is induced by the temperature quadrupole component at last scattering. When reionization occurs, there is a generation of polarized power on the scale corresponding to the acoustic horizon size at the reionization redshift. This particular signature is called the “reionization bump”, and is clearly visible in the bottom panel of Fig. 3 in the $`\mathrm{}10`$ region. The position of the bump in multipole space scales as $`\mathrm{}_{\text{bump}}\sqrt{z_r}`$ Zaldarriaga97 . As discussed above, a change in the helium fraction implies a shift of the redshift of reionization for a given (fixed) optical depth, Therefore the value of $`Y_p`$ has an effect on the position of the reionization bump in the polarization power spectrum, but not on its height, which is controlled by the optical depth and is proportional to $`\tau ^2`$. This effect is highlighted in the bottom panel: a 10% change in $`Y_p`$ induces roughly a 10% change in the position of the bump. The subsequent two oscillatory features for $`\mathrm{}<\mathrm{\hspace{0.33em}50}`$ reflect the displacement of further secondary, reionization induced polarization oscillations. However, since the value of polarized power is very low in that region, such secondary oscillations are very hard to detect precisely. In principle, given an accurate knowledge of the reionization history, the effect of $`Y_p`$ on the polarization bump would assist into determining the helium abundance. However, our ignorance of the reionization history prevents us from recovering useful information out of the measured reionization bump. The displacement induced by $`Y_p`$ is in fact degenerate with a partial reionization, or with other, more complex reionization mechanisms (see Holder03 ). Hence constraints on $`Y_p`$ come effectively from the damping tail in the $`\mathrm{}>\mathrm{\hspace{0.33em}400}`$ region of the temperature spectrum, which needs to be measured with very high accuracy.

Other light elements like deuterium and helium-3 are much less abundant, and will therefore have even smaller effect on the CMB power spectrum, at the order of $`10^5`$.

### III.2 Monte Carlo analysis

We use a modified version of the publicly available Markov Chain Monte Carlo package cosmomc cosmomc  as described in lewis02  in order to construct Markov Chains in our 7 dimensional parameter space. We sample over the following set of cosmological parameters: the physical baryon and CDM densities, $`\omega _b\mathrm{\Omega }_bh^2`$ and $`\omega _c\mathrm{\Omega }_ch^2`$, the cosmological constant in units of the critical density, $`\mathrm{\Omega }_\mathrm{\Lambda }`$, the scalar spectral index and the overall normalization of the power spectrum, $`n_s`$ and $`A_s`$ (see section III D below for a more precise definition), the redshift at which the reionization fraction is a half, $`z_r`$, and the primordial helium mass fraction, $`Y_p`$. We restrict our analysis to flat models, therefore the Hubble parameter, $`h=H_0/100\text{ km}\text{ s}^1\text{Mpc}^1`$, is a derived parameter, $`h=[(\omega _c+\omega _b)/(1\mathrm{\Omega }_\mathrm{\Lambda })]^{1/2}`$. We consider purely adiabatic initial conditions, and we do not include gravitational waves. In the CMB analysis, we assume 3 massless neutrino families and no massive neutrinos. We include the WMAP data kogut03 ; hinshaw03  (temperature and polarization) with the routine for computing the likelihood supplied by the WMAP team verde03 . We make use of the CBI pearson02  and of the decorrelated ACBAR acbar ; kuo03  band powers above $`\mathrm{}=800`$ to cover the small angular scale region of the power spectrum.

Since $`Y_p`$ is a rather flat direction in parameter space with present-day data, we find that a much larger number of samples is needed in order to achieve good mixing and convergence of the chains in the full 7D space. We use $`M=4`$ chains, each containing approximately $`N=310^5`$ samples. The mixing diagnostic is done on the same lines as in verde03 , by means of the Gelman and Rubin criterion gelman92 . The burn-in of the chains also takes longer than in the case where $`Y_p`$ is held fixed, and we discard 6000 samples per chain.

### III.3 CMB analysis results

Marginalizing over all other parameters, we find that the helium mass fraction from CMB alone is constrained to be $`Y_p<0.647`$ at 99% c.l. (1 tail limit), and

$$0.160<Y_p<0.501$$
(4)
at 68% c.l. (2 tails). Thus, for the first time the primordial helium mass fraction has been observed using the cosmic microwave background. However, present-day CMB data do not have sufficient resolution to discriminate between the astrophysical helium measurements, $`Y_p0.244`$, and the deuterium guided BBN predictions, $`Y_p0.248`$, which would require percent precision.

In Figure 4 we plot the marginalized and the mean likelihood of the Monte Carlo samples as a function of $`Y_p`$. If the likelihood distribution is Gaussian, then the 2 curves should be indistinguishable. The difference between marginalized and mean likelihood for $`Y_p`$ indicates that the marginalized parameters are skewing the distribution, and therefore that correlations play an important role. Although the mean of the 1D marginalized likelihood is rather high, $`(Y_p)=0.33`$, the mean likelihood peaks in the region indicated by astrophysical measurements, $`Y_p0.25`$. In view of this difference, it is important to understand the role of correlations with other parameters, and we will turn to this issue now.

In Figure 5 we plot joint 68% and 99% confidence contours in the ($`\omega _b,Y_p`$)-space. From the Monte Carlo samples we obtain a small and negative correlation coefficient between the two parameters $`\text{corr}(Y_p,\omega _b)=0.14`$. Baryons and helium appear to be anticorrelated simply because present-day WMAP data do not map the peaks structure to sufficiently high $`\mathrm{}`$. Precise measurements in the small angular scale region should reveal the expected positive correlation between the baryon and helium abundances, which is potentially important in order to correctly combine BBN predictions and CMB measurements of the baryon abundance. We turn to this question in more detail in the next section. In SBBN the baryon fraction and helium fraction are correlated along a different direction (cf. Fig. 5). However, this correlation is very weak, and the SBBN relation gives practically a flat line. Since the two parameters are not independent from the CMB point of view, it is in fact not completely accurate to do the CMB analysis with fixed helium mass fraction of $`Y_p=0.24`$ to get the error-bars on the baryon fraction, and then re-input this baryon fraction (and error-bars) to predict the helium mass fraction from BBN. The most accurate procedure is to analyse the CMB data leaving $`Y_p`$ as a free parameter, thereby obtaining the correct (potentially larger) error-bars on $`\omega _b`$ upon marginalization over $`Y_p`$.

In view of the emerging baryon tension between CMB and BBN, it is important to check whether allowing helium as a free parameter can significantly change the CMB determination of the baryon density or its error. In order to evaluate in detail the impact of $`Y_p`$ on the error-bars for $`\omega _b`$, we consider the following 3 cases.

* The usual case, when the helium fraction for the CMB analysis is assumed to be known a priori and is fixed to the canonical value $`Y_p=0.24`$.
* A case with a weak astrophysical Gaussian prior on the helium fraction, which we take to be $`Y_p=0.24\pm 0.01`$. As discussed above, the error-bars of the astrophysical measurements are typically a factor 5 tighter than this, but our prior is chosen to encompass the systematic spread between the different observations.
* The case in which we assume a uniform prior for $`Y_p`$ in the range $`0Y_p1`$, i.e. $`Y_p`$ is considered as a totally free parameter.

We do not find any significant change in the error-bars for $`\omega _b`$ in the 3 different cases. The confidence intervals on $`\omega _b`$ alone are determined to be (case (c)) $`0.0221<\omega _b<0.0245`$ at 68% c.l. ($`0.0204<\omega _b<0.0276`$ at 99 % c.l.). The standard deviation of $`\omega _b`$ as estimated from the Monte Carlo samples is found to be $`\widehat{\sigma }_b=1.310^3`$. This is in complete agreement with the error-bars on $`\omega _b`$ obtained by the WMAP team for the standard $`\mathrm{\Lambda }`$CDM case spergel03 . We conclude that at the level of precision of present-day CMB data, it is still safe to treat the baryon abundance and the helium mass fraction as independent parameters. This result is non-trivial, since the fact that the damping tail is not yet precisely measured above the second peak would a priori suggest that degeneracies between $`Y_p`$, $`\omega _b`$ and $`n_s`$ could potentially play a role once the assumption of zero uncertainity on $`Y_p`$ is relaxed. The impact of $`Y_p`$ is small enough, and the error-bars on $`\omega _b`$ large enough that a uniform prior on $`Y_p`$ can still be accomodated within the uncertainity in the baryon abundance obtained for case (a). However, the $`Y_p\omega _b`$ correlation will have to be taken into account to correctly analyse future CMB data, with a quality such as Planck. We discuss this potential in the next section.

We observe the expected correlation between the redshift of reionization and the helium fraction (Fig. 6), which is discussed above. The correlation coefficient between the two parameters is found to be rather large and positive, $`\text{corr}(Y_p,z_r)=0.40`$. This correlation produces a noticeable change in the marginalized 1D-likelihood distribution for $`z_r`$ as we go from case (a) to case (c). Marginalization over the additional degree of freedom given by $`Y_p`$ broadens considerably the error-bars on $`z_r`$. In fact, the 68% confidence interval for $`z_r`$ increases by roughly 20% (and shifts to somewhat higher values), from $`10.220.9`$ (case (a)) to $`10.623.3`$ (case (c)). Case (b) exhibits similar error-bars as case (a). On the other hand, the determination of the reionization optical depth is not affected by the inclusion of helium as a free parameter, giving in all cases $`0.08<\tau _\text{r}<0.23`$. Correspondingly, the correlation is less significant, $`\text{corr}(Y_p,\tau _\text{r})=0.11`$. We therefore conclude that the differences in the determination of $`z_r`$ are due only to the variation of the amount of electrons available for reionization as $`Y_p`$ is changed.

Leaving $`Y_p`$ as a free parameter also has an impact on the relation between $`\omega _b`$ and the scalar spectral index, $`n_s`$. The extra power suppression on small scales which is produced by a larger $`Y_p`$ can be compensated by a blue spectral index (see Fig. 7).

### III.4 Potential of future CMB observations

In order to estimate the precision with which future satellite CMB measurements will be able to constrain the helium mass fraction we perform a Fisher matrix analysis (FMA). This technique approximates the likelihood function with a Gaussian distribution around a fiducial model, which is assumed to be the best fit model. The Fisher information matrix $`F`$ gives the second order expansion of the likelihood around its peak, and it is computed from the derivatives of the power spectrum with respect to the cosmological parameters. The expected performance of the experiment can be modelled with a noise contribution to the likelihood function, which is described in terms of a few experimental parameters. The covariance matrix $`C`$ is then given by the inverse of the Fisher matrix, $`C=F^1`$. It is then straightforward to evaluate the expected 1-$`\sigma `$ error on parameter $`i`$, which is given by $`\sqrt{c_{ii}}`$ (all other marginalized). The main advantage of the FMA is that it gives reliable and accurate predictions (including information on the expected degeneracies) with minimal computational effort. For further details on the Fisher Matrix formalism, see e.g. knox95 ; bond97 ; zaldarriaga97 ; tegmark96 ; EB99 ; bowen ; alpha1 ; alpha2 .

#### III.4.1 Parameters set

In order to obtain a reliable prediction, it is extremely important to choose a parameter set wrt which the dependence of the CMB power spectrum is as linear and uncorrelated as possible. This issue has been discussed exhaustively in Ref. kosowsky02 , where the authors introduce a set of “physical parameters” which satisfies the above requirements. In the present work we retain most of the physical parameters defined in ref. kosowsky02 : the ratio between the sound horizon at decoupling and the angular diameter distance $`𝒜`$, the baryon density $`=\mathrm{\Omega }_bh^2`$, the energy density in the cosmological constant $`𝒱=\mathrm{\Omega }_\mathrm{\Lambda }h^2`$, the matter-radiation density at decoupling $``$ and $``$, which is mainly a function of the matter and radiation content. We adopt a slightly different choice for the physical parameter describing reionization. For adiabatic perturbations, the initial power spectrum of the gauge invariant curvature perturbation $`\zeta `$ is written as

$$P_\zeta (k)=A_s\left(\frac{k}{k_0}\right)^{n_s1}$$

(and we do not take running into account). The quantity $`\zeta `$ corresponds to the intrinsic curvature perturbation on comoving hypersurfaces, and at the end of inflation is related to the gravitational potential perturbation, $`\mathrm{\Psi }`$, by $`\zeta =\frac{3}{2}\mathrm{\Psi }`$ (see e.g. durrer  for more details). We take the pivot-scale $`k_0`$ to be $`k_0=0.05\text{ Mpc}^1`$. If $`\tau _\text{r}`$ denotes the optical depth to reionization, then defining $`𝒯=A_s\mathrm{exp}(2\tau )`$ is a good way to take into account the degeneracy between the optical depth and normalization. Our parameter set contains then the six above physical parameters ($`𝒜,,𝒱,,,𝒯`$), the power spectrum normalization $`A_s`$, the scalar spectral index $`n_s`$ and the helium fraction $`Y_p`$.

The choice of the physical parameter set makes it easy to implement in the FMA interesting theoretical priors. For instance, we are interested in imposing flatness in our forecast, in order to be able to directly compare present-day accuracy on $`Y_p`$ with the potential of Planck and CVL. The prior on the curvature of the universe is imposed in the FMA by fixing the value of the parameter $`𝒜`$ to the one of the fiducial model. In fact, the parameter $`𝒜`$ is a generalization of the shift parameter, which describes the sideway shift of the acoustic peak structure of the CMB power spectrum as a function of the geometry of the universe and its content in matter, radiation and cosmological constant. Although imposing $`𝒜=`$ const is not the same as having curvature=constant over the full range of cosmological parameters, for the purpose of evaluating derivatives the two conditions reduce to the same. The fact that our fiducial model is actually slightly open (see below), does not make any substantial difference in the results, apart from reducing the numerical inaccuracies which would arise had we computed the derivatives around an exactly flat model. We can also easily impose a prior knowledge of the helium fraction, by fixing the value of $`Y_p`$, as it is usually done in present CMB analysis, and investigate how this modifies the expected error on the the baryon density.

#### III.4.2 Accuracy issues

We numerically compute double sided derivative of the power spectrum around the fiducial model with cosmological parameters given in Table 1. We find it necessary to increase the accuracy of CAMB by a factor of 3 in each of the “accuracy boost” values. As a fiducial model, we use the best fit model to the WMAP data for the standard $`\mathrm{\Lambda }`$CDM scenario, as given in Table 1 of ref. spergel03 . However, in order to avoid numerical inaccuracies which arise when differentiating around a flat model, we reduce slightly the value of $`\mathrm{\Omega }_\mathrm{\Lambda }`$ by imposing an open universe, $`\mathrm{\Omega }_{\text{tot}}=0.99`$. We perform the FMA for the expected capabilities of Planck’s High Frequency Instrument (HFI) and for an ideal CMB measurement which would be cosmic variance limited (CVL) both in temperature and in E-polarization (and we do not consider the B-polarization spectrum), and therefore represents the best possible parameter measurement from CMB anisotropies alone. The complicated issues coming from foreground removals, point source subtractions, etc. are assumed to be already (roughly) taken into account in the experimental parameters for the experiment. Those are the effective percentual sky coverage $`f_{\text{sky}}`$, the number of channels, the sensitivity of each channel $`\sigma _c^{T,E}`$ for temperature (T) and E-polarization (E) in $`\mu `$K and the angular resolution $`\theta _c^{T,E}`$ (in arcmin). For Planck HFI, we take the 3 channels with frequencies 100, 143 and 217 GHz, with respectively $`\sigma _{c=1,2,3}^T=5.4,6.0,13.1`$ and $`\sigma _{c=2,3}^E=11.4,26.7`$ and we have $`f_{\text{sky}}=0.85`$ PlanckHP  Since the CVL is an ideal experiment, we put its noise to zero and assume perfect foregrounds removal, so that $`f_{\text{sky}}=1`$. In order to test the accuracy of our predictions and compare present-day results with the forecasts, we also perform an FMA with WMAP first year parameters, obtaining excellent agreement between the FMA results and the error-bars from actual data. For the purpose of comparison, we include forecasts for the full WMAP 4 years mission, which will also measure E-polarization and reduce present-day errors on the temperature spectrum by a factor of 2. We limit the range of multipoles to $`\mathrm{}<2000`$, because at smaller angular scales non-primary anisotropies begin to dominate (Sunyaev-Zeldovich effect). The authors of ref. precision03  discuss the issue of numerical precision of 3 different CMB codes and conclude that they are accurate to within $`0.1\%`$. While this is encouraging, it is not of direct relevance to this work, since what matters in the computation of derivatives is not much the absolute precision of the spectra, but rather their relative accuracy.

#### III.4.3 Forecasts and discussion

Table 2 summarizes our forecasts for the future measurements and compares them with the results obtained from WMAP actual data.

We notice that when the WMAP full 4 years data will be available (including E-polarization), the error on the baryon density is expected to decrease by a factor of 2 to $`2.86\%`$, compared to today’s $`5.04\%`$ (assuming flatness). Nevertheless, inclusion of $`Y_p`$ as a free parameter will still have no effect on the determination of $`\omega _b`$ for WMAP, i.e. $`Y_p`$ will remain an essentially flat direction when marginalized over. While the determination of the helium fraction will improve, the FMA cannot reliably assess quantitatively how much, since for such large errors the likelihood distribution is not Gaussian and the quadratic approximation breaks down. In the table we therefore give the FMA estimation as an indication, with the caveat that the Fisher approximation is likely to be inaccurate for the real errors on $`Y_p`$ from WMAP’s 4 years data.

It is interesting that for Planck the effect of the helium fraction can no longer be neglected. Inclusion of the helium fraction increases the error on $`\omega _b`$ by roughly $`80\%`$, from $`0.70\%`$ to $`1.26\%`$. The correlation between the two parameters will have to be taken into account, as is evident from Fig. 8. The expected correlation coefficient is $`\text{corr}(Y_p,\omega _b)=0.84(0.91)`$ for Planck (for CVL). The expected $`1\sigma `$ error on $`Y_p`$ is about $`5\%`$ for Planck, or $`\mathrm{\Delta }Y_p0.01`$. This is of the same order as the spread in current astrophysical measurements. We conclude that in Planck-accuracy data analysis it will be necessary to include the uncertainty in the determination of the helium mass fraction, at least in the form of a Gaussian prior over $`Y_p`$ of the type we used in the CMB data analysis presented above.

Finally, measuring CMB temperature and polarization with cosmic variance accuracy would allow to constrain $`Y_p`$ to within $`1.5\%`$, or $`\mathrm{\Delta }Y_p0.0036`$ (assuming flatness). Such an ideal measurement would be able to discriminate between the BBN-guided, deuterium based helium value and the current lowest direct helium observations (cf. Fig. 1).

Our forecasts for the uncertainty in the Helium mass fraction from future observations are in excellent agreement with the findings of Ref. knox03 . There, the standard deviation on $`Y_p`$ for Planck is estimated to be $`\mathrm{\Delta }Y_p=0.012`$. The authors of Ref. knox03  also consider an experiment (CMBPol) with characteristics similar to our CVL, for which they forecast $`\mathrm{\Delta }Y_p=0.0039`$, again in close agreement with our result. An earlier work eisenstein  found for Planck (temperature and polarization) $`\mathrm{\Delta }Y_p=0.013`$, in satisfactory concordance with our result. It should be noticed that the forecast reported for MAP in Table 2 of Ref. eisenstein , namely $`\mathrm{\Delta }Y_p=0.02`$, is nothing but the Gaussian prior $`Y_p=0.24\pm 0.02`$ which was assumed in their analysis.

The main source of improvement for the determination of $`Y_p`$ will be the better sampling of the temperature damping tail provided by Planck and the CVL. Polarization measurements have mainly the effect of reducing the errors on other parameters. In fact, we have checked that excluding from our FMA the $`2\mathrm{}50`$ region of the E-polarization and ET-correlation spectra changes the forecast precision on $`Y_p`$ less than about 10-15% for Planck and less than a few percent for CVL. This supports the conclusion that the low-$`\mathrm{}`$ reionization bump is not very useful in measuring the helium abundance, because of the degeneracy with $`z_r`$.

## IV Conclusions

We have analysed the ability of CMB observations to determine the helium mass fraction, $`Y_p`$. We find that present data only allow a marginal detection, $`0.160<Y_p<0.501`$ at 68% c.l.. This determination is completely independent from the usual astrophysical observations and uses CMB data only. We discuss degeneracies between $`Y_p`$ and other cosmological parameters, most notably the baryon abundance, the redshift and optical depth of reionization and the spectral index. We conclude that present-day CMB data accuracy does not require the inclusion of $`Y_p`$ as a free parameter. We find that Planck will determine the helium mass fraction within $`5\%`$ (or $`\mathrm{\Delta }Y_p0.01`$), which however will only allow a marginal discrimination between different astrophysical measurements. Nevertheless, we point out that the uncertainty of the helium fraction will have to be taken into account in order to correctly estimate the errors on the baryon density from Planck. To determine if the emerging baryon tension (from BBN) is related to underestimated systematic error-bars or whether it is an indication of new physics, CMB observation will have to be pushed very close to the cosmic variance limit in both temperature and polarization.

## Acknowledgement

It is a pleasure to thank Ruth Durrer, Lloyd Knox, Samuel Leach, Anthony Lewis, Christophe Ringeval, Dominik Schwarz and Gary Steigman for useful comments and discussions. We thank the anonymous referee for many useful suggestions. This work was performed on the SUN Enterprise 10000 Supercomputer owned and operated by the University of Geneva. R.T. is partially supported by the Swiss National Science Fundation, the Schmidheiny Fundation and the European Network CMBNET. S.H. thanks the Tomalla foundation for support.


text:
# Substructures in Cold Dark Matter Haloes

## 1 Introduction

The formation and evolution of structure in the Universe is a topic of fundamental interest. In the last decades, the Cold Dark Matter (CDM) model has been extensively studied and has had considerable success in reproducing observational results, both on galactic and on cluster scales. In fact, the CDM model with the ‘concordance’ set of cosmological parameters ($`\mathrm{\Lambda }`$CDM) has been so successful that it can now be considered a standard paradigm for the formation of structure in the Universe. According to this model, the dominant force that drives structure formation is gravity, and large systems like galaxy clusters are formed via hierarchical merging of smaller structures.

Numerical simulations of gravitational clustering of dark matter are an indispensable tool for investigating the non-linear growth of structures in its full geometrical complexity. Until recently, dissipationless simulations suffered from the so-called *overmerging* problem, i.e. substructures disrupt very quickly within dense environments (Katz & White, 1993). However, both analytic work (Moore, Katz & Lake 1996) and high resolution simulations (Tormen, Bouchet & White 1997; Ghigna et al., 1998; Klypin et al., 1999b; Ghigna et al., 2000) have demonstrated that the cores of dark matter haloes that fall into a cluster can actually survive as self-gravitating objects orbiting in the smooth dark matter background of the cluster, provided high enough force and mass resolution are used. Recent high-resolution simulations (Ghigna et al., 2000; Springel et al., 2001b) have also shown that the abundance of these substructures is in agreement with the observed abundance of galaxies in clusters, suggesting a natural one-to-one identification of luminous cluster galaxies and dark matter substructures.

Another interesting claim is that the shape of the substructure mass function is independent of the mass of the parent halo (Moore et al., 1999b). It is not obvious that this should be the case, because in CDM cosmologies, the initial conditions do depend on scale and galaxies form several billion years before clusters. As Moore et al., (1999b) pointed out, the logarithmic slope of the power spectrum asymptotically approaches $`3`$ on small scales, so clumps of widely different (yet sufficiently small) mass tend to collapse simultaneously, and as a result the timescale between the collapse of the first substructures and their incorporation into larger haloes becomes shorter. One might therefore expect that substructures were more easily disrupted in low mass haloes, producing a substructure mass function which depends on the mass of the parent halo.

Observationally, the predicted abundance of substructures in clusters is one of the major successes of the CDM model (Springel et al., 2001b), but on galactic scales, it appears that simulations predict more substructures than are visible by almost two orders of magnitude (Kauffmann, White & Guiderdoni 1993; Moore et al., 1999a; Klypin et al., 1999b; Tasitsiomi, 2002). This is commonly referred to as the ‘dwarf galaxy crisis’ of CDM. There have been suggestions that the solution to this problem lies in processes such as heating by a photo-ionising background that suppresses star formation in small haloes at early times (Efstathiou, 1992; Bullock, Kravtsov & Weinberg 2000; Somerville, 2002; Benson et al., 2002). Alternatively, it has been suggested that the nature of dark matter may be different than assumed in the canonical $`\mathrm{\Lambda }`$CDM model, for example by being warm or self-interacting, both of which could selectively eliminate small-scale structure. However, self-interactions appear to be relatively ineffective in reducing the number of subhalos, unless the assumed cross-section is unreasonably large (Colin et al., 2002).

A less drastic resolution was suggested by Stoehr et al., (2002) who noted that it might be possible to identify the observed Galactic satellites with the few most massive subhalos and that the rest contain no stars. Direct evidence for the large population of dark satellites predicted by CDM models comes from the anomalous flux ratios of multiply imaged quasars (Mao & Schneider, 1998; Chiba, 2002; Dalal & Kochanek, 2002).

So far, a detailed numerical analysis of substructures has only been carried out in high resolution re-simulations of a few individual haloes (Moore et al., 1999b; Ghigna et al., 2000; Springel et al., 2001b). In this paper, we carry out a systematic analysis of substructures as a function of the mass of the parent halo and as a function of environment. We study the mass functions of subhalos, their radial distributions and their merging and mass accretion histories. These quantities are of fundamental interest for galaxy formation, because dark matter haloes and substructures represent the birth places of luminous galaxies. Their accretion and merging histories regulate the rate at which baryons can cool, determining in this way the rate at which stars form in galaxies as a function of cosmic time.

The layout of the paper is as follows: in Section 2 we describe the simulations that we use; in Section 3 we give a brief description of the algorithm used to find substructures in haloes; in Section 4 we present the subhalo mass function, in Section 5 we analyse the mass distribution of the largest substructures; in Section 6 we study the radial distribution of substructures; in Section 7 we discuss the merging and mass accretion histories of substructures, both as a function of mass and as a function of environment. A summary and a discussion of the results obtained are presented in Section 8.

## 2 Cluster simulations

In this study, we use collisionless simulations of clusters generated using the ‘zoom’ technique (Tormen et al. 1997; Katz & White, 1993). First, a cosmological simulation of a large region is used to select a suitable target cluster. The particles in the final cluster and its surroundings are then traced back to their initial Lagrangian region and are replaced by a larger number of lower mass particles. These are perturbed using the same fluctuation distribution as in the parent simulation, but now extended to smaller scales to account for the increase in resolution. This resampling of the initial conditions of the Lagrangian region of the cluster thus allows a localised increase in mass and force resolution. Outside the *high-resolution* region, particles of variable mass, increasing with distance, are used so that the computational effort is concentrated on the cluster of interest, while still maintaining a faithful representation of the large-scale density and velocity field of the parent simulation.

In this paper, we study a set of $`11`$ high-resolution re-simulations of galaxy clusters ($`5`$ of mass $`10^{14}h^1\mathrm{M}_{}`$ and $`6`$ of mass $`10^{15}h^1\mathrm{M}_{}`$), and a high resolution re-simulation of a ‘typical’ region of the Universe. The simulations were carried out with the parallel tree-code GADGET (Springel, Yoshida & White 2001).

One of our massive clusters was taken from the ‘S-series’ studied by Springel et al., (2001b), where the parent simulation employed was the GIF-$`\mathrm{\Lambda }`$CDM model carried out by the Virgo Consortium (Kauffmann et al., 1999). This parent simulation followed $`256^3`$ particles of mass $`1.4\times 10^{10}h^1\mathrm{M}_{}`$ within a comoving box of size $`141.3h^1`$Mpc on a side. The other cluster re-simulations and the simulation of the field region were selected from the VLS simulation carried out by the Virgo Consortium (Jenkins et al., 2001; Yoshida, Sheth & Diaferio 2001). The simulation was performed using a parallel P3M code (Macfarland et al., 1998) and followed $`512^3`$ particles with a particle mass of $`7\times 10^{10}h^1\mathrm{M}_{}`$ in a comoving box of size $`479h^1`$Mpc on a side. In all cases, the parent simulation and the re-simulations were characterised by the following cosmological parameters: $`\mathrm{\Omega }_0=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, spectral shape $`\mathrm{\Gamma }=0.21`$, $`h=0.7`$ (we adopt the convention $`H_0=100h\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$) and normalisation $`\sigma _8=0.9`$.

In Table 1, we summarise some important numerical parameters of the simulations used. We will refer to our five high-mass clusters of mass resolution $`2\times 10^9h^1\mathrm{M}_{}`$ as type ‘B1’, and to the low-mass clusters as type ‘B2’. These simulations were carried out by Barbara Lanzoni as part of her PhD thesis and were previously used in Lanzoni et al., (2003). The ‘S2’ simulation is taken from the ‘S-series’ of Springel et al., (2001b). The field region ‘M3’ was adopted from the ‘M-series’ studied by Stoehr, (2003).

## 3 Identification of dark matter substructures

The identification of substructures in dark matter haloes is a difficult technical problem and many different algorithms have been developed to accomplish this task, for example the hierarchical friends-of-friends algorithm (HFOF) (Gottlöber, Klypin & Kravtsov 1998), the bound density maximum algorithm (BDM) (Klypin et al., 1999a), and SKID (see http ref: http://www-hpcc.astro.washington.edu/tools). Each of these algorithms has its own advantages and weaknesses, so that arguably none of them is completely satisfactory yet. In this work, we use the algorithm SUBFIND proposed by Springel et al., (2001b), which combines ideas used in other group finding techniques with a topological approach for finding substructure candidates. SUBFIND can handle haloes of arbitrary shape, does not require an iterative procedure for finding subhalo candidates, and is capable of detecting arbitrary levels of ‘subhalos within subhalos’. In this section, we briefly summarise how the method works.

In a first step, a standard friends-of-friends (FOF) algorithm is used to identify virialized parent haloes. The FOF algorithm links together all particle pairs with separation less than a linking length $`b`$. We adopt the standard value $`b=0.2`$ in units of the mean particle separation, which selects groups of particles with overdensities close to the value predicted by the spherical collapse model for the virialized regions of haloes. The next step is to compute an estimate of the local density at the position of each particle in the group. To this end, we employ an adaptive kernel interpolation method similar to the one used in smoothed particle hydrodynamics. In the resulting density field, we define as substructure candidates locally overdense regions which are enclosed by isodensity contours that traverse a saddle point. Our method for finding these regions can be visualised as follows: we reconstruct the density field by considering particles in order of decreasing density, thus working our way from high to low density. This corresponds to gradually lowering a global threshold in the density field sampled by the simulation particles. Isolated overdense regions grow slowly in size during this process. When two such separate regions coalesce to form a single region, their density contours join at a saddle point. Each time such an event occurs, we have found two substructure candidates.

After the regions containing substructure candidates have been identified, we apply an unbinding procedure where we iteratively reject all particles with positive total energy in order to eliminate ‘background’ particles that do not belong to the subhalo. For the purposes of this study, we consider all substructures that survive this unbinding procedure, and still have at least $`10`$ self-bound particles, to be genuine subhalos.

In summary, the algorithm SUBFIND decomposes a given particle group into a set of disjoint and self-bound substructures, each of which is identified as a locally overdense region in the density field of the background halo. Note that SUBFIND classifies all the particles inside a FOF group either as belonging to a bound substructure or as being unbound. The self-bound part of the FOF background halo itself will then also appear in the substructure list. We will exclude it when referring to subhalos or substructures in the following analysis.

## 4 The subhalo mass function

The sample of parent haloes used for studying the mass function analysis consists of $`6`$, $`5`$, $`34`$ and $`100`$ haloes in the mass ranges $`8.68\times 10^{14}`$$`1.79\times 10^{15}`$ (from simulations B1 and S2), $`6.99\times 10^{13}`$$`1.27\times 10^{14}`$ (from simulations B2), $`7.0\times 10^{12}`$$`2.0\times 10^{13}`$ (from simulation M3) and $`7.0\times 10^{11}`$$`2.0\times 10^{12}h^1\mathrm{M}_{}`$ (from simulation M3).

The resulting subhalo mass functions are shown in Fig. 1. In the first four panels, we plot the differential mass functions for parent haloes of different mass. The histograms are computed by stacking all the haloes in the given range of mass and the error bars represent Poisson errors. The solid line in each of the panels is a power-law fit to the measured differential mass function; the fit is performed using the least absolute deviation method over the range of mass shown by the line. In all the cases the slope of this unrestricted fit is close to $`1`$ (it is equal to $`0.98`$ for the top left panel, $`0.97`$ for the top right panel, $`1.11`$ for the middle left panel and $`1.13`$ for the middle right panel). However, we note that the lowest mass bins, which have the smallest statistical errors, are best fit with a slightly shallower slope: if we restrict the fit to the $`4`$ lowest mass bins the slope is $`0.94`$ for the top left panel and $`0.85`$ for the middle left panel. These are closer to the value $`0.8`$, measured by Helmi, White & Springel (2002) for a single cluster simulation of extremely high-resolution. Also note that a slope shallower than $`1`$ at the low-mass end implies that the integrated mass in substructures remains bounded and is dominated by the most massive subhalos. It is likely that our subhalo mass functions are steepened somewhat by a cut-off in abundance for very massive substructures.

The bottom left panel of Fig. 1 shows the cumulative mass function for all the haloes used in the sample. To compare the different subhalo mass functions, we have rescaled the subhalo mass by dividing by the virial mass of the parent halo. Each line represents the average cumulative mass function over all the haloes in each mass bin. Note that in this paper we define the ‘virial mass’ as $`M_{200}`$, the mass within a sphere of density $`200`$ times the critical value at redshift zero. The lines end at different places because of the differing mass resolution of the simulations (see Table 1). We find that all four cumulative mass functions agree within the statistical errors. Finally, in the bottom right panel, we show the differential subhalo mass functions in units of rescaled mass.

We note that the ‘universality’ of the subhalo mass function seen here appears to be quite robust with respect to numerical resolution. In Fig. 2 we compare the average cumulative mass functions for haloes with mass $`10^{14}h^1\mathrm{M}_{}`$ from simulations B2 and M3. We here averaged $`5`$ haloes for simulations B2, and $`4`$ for simulation M3, to reduce the object-to-object scatter that is unavoidable for small numbers of subhalos. Despite an order of magnitude difference in numerical resolution, the agreement between the simulations is good. We are able to resolve the right number of objects in the low-resolution simulation above its resolution limit (shown as a vertical dotted line in the figure). A similar result was obtained by Springel et al., (2001b, see their Fig. 5) using a set of $`4`$ re–simulations of the same cluster with systematically increasing resolution, thereby allowing a direct study of numerical convergence. This showed in particular that the S2 simulation used here has well converged to the properties of a much higher resolution simulation above its own resolution limit, as used here. Further support for our results was also found by Stoehr et al., (2003, see their Fig. 3). They compared the S–series simulations from Springel et al., (2001b) with an extremely well resolved re–simulation of a Milky–Way sized halo. This latter simulation used an updated version of the simulation code and more conservative integration parameters than used here (following Power et al., 2003), suggesting that the subhalo mass function is a relatively robust quantity and that the simulations we discuss here are adequate for our purposes.

As a further check of the robustness of our results we also compare the internal structure of subhalos drawn from our different simulations. Fig. 3 shows the correlation between the substructure mass and the third power of the maximum circular velocity, $`V_{\mathrm{max}}`$, measured directly from the circular velocity curve of each subhalo. Different symbols are used for substructures drawn from different simulations. Note that for the range of masses shown in the plot, substructures drawn from simulation M3 contain at least $`60`$ particles. While the scatter is clearly large for haloes with such a low number of particles, the good general agreement between the runs suggests that the smallest substructures in our lower resolution simulations have an internal structure that is still reliably resolved, at least in a statistical fashion.

Our results confirm the conclusion drawn by Moore et al., (1999b): the mass function of substructures appears to be almost independent of the mass of the parent halo. While our results are consistent with such a ‘scale–free’ subhalo mass function, the halo–to–halo scatter in our simulation set is quite large, preventing us from putting tight constraints on the accuracy with which the ‘scale–free’ subhalo mass function is preserved when haloes of different mass are considered. And hence there is still room for weak trends with mass. A clear detection of these would require simulations with larger dynamic range, and larger samples of simulated haloes for each mass bin.

As we discuss in more detail in Sec. 7.2, our findings suggest that the destruction of satellites due to the physical processes of dynamical friction and tidal stripping on one hand, and the accretion of new satellites on the other hand, cancel out in such a way that the subhalo mass function does not depend or at best very weakly depends on the mass of the parent halo. The reason for the invariance of the subhalo mass function may lie in the physical nature of this dynamical balance, which may be insensitive to the slightly broken scale-invariance of dark haloes themselves. This shows up as a mass-dependence of halo concentrations, for example. Some fully analytic models for the subhalo abundance have been developed (e.g. Sheth, 2003), but they are presently not able to account for mass-loss and dynamical friction self-consistently, and so provide little guidance to answer this interesting question. A full understanding of the apparent ‘conspiracy’ that establishes an almost mass-invariant subhalo mass function will therefore require further simulations.

## 5 The most massive substructures

In this section we investigate whether the properties of the largest substructures depend on the mass of the parent halo. This is interesting since the largest substructures mark the sites where one expects to find the brightest galaxies.

In the following, $`M_1`$ refers to the mass of the most massive subhalo and $`M_2`$ to the mass of the second most massive subhalo within the virial radius of a given object of virial mass $`M_{200}`$. Note that we have excluded from our analysis the subhalo associated with the FOF group itself. In a semi-analytic scheme, it is this ‘subhalo’ that would host the brightest cluster galaxy (BCG). In Sec. 7.2 we will show that, once accreted onto a massive halo, substructures suffer significant stripping, an effect that is more important for substructures accreted at higher redshift. It is then likely that the largest substructures we find within the virial radius at the present time were accreted at relatively low redshift.

In Fig. 4, we plot $`M_1/M_{200}`$ as a function of $`M_{200}`$ for $`434`$ haloes drawn from all the simulations listed in Table 1. This sample includes not only the central clusters in our re-simulations, but also the other haloes found in the high-resolution regions around the re-simulated objects down to a mass limit of $`10^{13}h^1\mathrm{M}_{}`$. We took care however to exclude *contaminated* haloes that contained low resolution particles. In simulation M$`3`$, we selected only haloes with a mass larger than $`10^{12}h^1\mathrm{M}_{}`$. The small symbols in Fig. 4 indicate the value of $`M_1/M_{200}`$ measured for each individual halo, while the filled circles represent the median of the distribution. We have taken bins in $`M_{200}`$ such that there are an equal number ($`143`$) of haloes in each bin, except for the last six points, which we treated as a separate bin, corresponding to the central cluster haloes in simulations B$`1`$ and S$`2`$. The error bars mark the $`20`$ and $`80`$ per cent percentiles of the distribution.

The results in Fig. 4 suggest that $`M_1/M_{200}`$ depends very little on the mass of the parent halo. Interestingly, $`M_1/M_{200}`$ appears to exhibit less scatter for the most massive haloes, but the number of simulated clusters we have in this high mass regime is rather small, so it is unclear whether this effect is statistically significant.

The results in Fig. 4 imply that the median value of $`M_1`$ increases in proportion to the mass of the cluster, suggesting that second ranked galaxies will be more luminous in more massive haloes. Note also that the mass of the largest substructure within the virial radius is typically only a few per cent of the virial mass.

In Fig. 5, we show the ratio $`M_2/M_1`$ as a function of the mass of the parent halo. Once again there is rather little dependence on $`M_{200}`$ with a possible decrease in the scatter for more massive haloes. Note that the median value of $`M_2/M_1`$ is around $`0.5`$. If the stellar masses of the second and third brightest galaxies in a cluster scale simply with the masses of their dark subhalos, they should have K-band luminosities that are equal to within $`0.5`$ mag.

## 6 The radial distribution of subhalos

The large sample of subhalos in our simulations allows us to study their radial distribution and to investigate how it depends on the mass of the parent halo. In Fig. 6, we plot the number density of substructures as a function of the normalised distance $`R/R_{200}`$ from the centre of the halo, defined here as the position of the most bound particle in the halo. We show averaged results for haloes with masses $`10^{15}h^1\mathrm{M}_{}`$, $`10^{14}h^1\mathrm{M}_{}`$, and $`10^{13}h^1\mathrm{M}_{}`$, and we limit the analysis to subhalos with masses greater than a fixed fraction ($`2\times 10^4`$) of the virial mass of the parent halo. This fraction is chosen because it lies just above the mass limit where it is possible to identify substructures in all of our simulations. As Fig. 1 shows, there are typically $`50`$ subhalos per parent halo with $`M_{\mathrm{sub}}/M_{200}>2\times 10^4`$, so by stacking a large sample of haloes, it is possible to calculate an average density profile that has rather little noise. Note that the densities plotted in Fig. 6 have been normalised to the mean density inside the virial radius. The solid lines with symbols show results for the $`3`$ different parent halo mass ranges defined above. For comparison, we have also plotted the dark matter radial profiles as dashed lines. Note that a small shift in the abscissa has been added to make the plot more readable.

We find that the subhalo profiles are ‘anti-biased’ relative to the dark matter in the inner regions of the haloes. This agrees with the results of Ghigna et al., (2000). Surprisingly we also find that the radial number density profiles are steeper in low mass haloes than in high mass haloes, a finding that reserves further investigation.

We now use our highest resolution cluster simulation to investigate whether subhalos of different mass have different radial profiles. In Fig. 7, we show the cumulative fraction of substructures as a function of $`R/R_{200}`$ for subhalos with $`M_{\mathrm{sub}}>0.01M_{200}`$ (solid line) and $`M_{\mathrm{sub}}0.01M_{200}`$ (dashed line).

As Fig. 1 already made clear, there are many more substructures with $`M_{\mathrm{sub}}0.01M_{200}`$ than with $`M_{\mathrm{sub}}>0.01M_{200}`$ ($`9749`$ versus $`96`$). Fig. 7 now shows that more massive substructures are preferentially located in the external regions of their parent haloes. This can be understood as a consequence of tidal truncation and stripping effects that quickly decrease the mass of subhalos as they fall into the cluster and reach the dense inner cores of the parent haloes (see Section 7.2 for a more quantitative analysis of mass-loss due to stripping).

Also note that this finding can be naturally explained as a consequence of the orbital decay experienced by substructures. As shown in Tormen, Diaferio & Syer (1998), the orbital decay is consistent with expectations based on the combined effects of dynamical friction and mass-loss. As a result, massive substructures are driven to the centre more rapidly than less massive ones: Tormen et al. show that the orbital decay occurs in less than a Hubble time if the initial mass of substructures is larger than $`1`$ per cent of the mass of the main cluster, while the substructures can retain their identity for a significant fraction of the Hubble time if their mass is smaller than $`5`$ per cent of the main cluster mass. Once driven to the centre, massive substructures are destroyed and no longer distinguishable from the central halo; this naturally explains the mass segregation that we see in our simulations.

In Fig. 8, we show the cumulative fraction of the total mass of the parent halo that is in substructures as a function of normalised distance from the halo centre. The different lines and symbols have the same meaning as in Fig. 6 and represent median relations for all the haloes in each mass bin. The error bars mark the $`20`$ and $`80`$ per cent percentiles of the distribution. The mass fraction in substructures rises from $`1`$ per cent at a radius $`0.3R_{200}`$ to $`6`$ per cent at $`rR_{200}`$. Note that the total mass fraction is dominated by the small number of most massive subhalos, and is hence a rather noisy quantity that shows large variations from system to system. The mass fraction may also reach values above $`10`$$`15`$ per cent, but then the underlying FOF parent halo is typically quite aspherical, with the most massive subhalo often lying outside the formal radius $`R_{200}`$ of the parent halo.

## 7 Subhalo histories

So far we have analysed the properties of subhalos only at redshift $`z=0`$. In this section, we turn to the time evolution of the masses of subhalos and their merging histories.

In order to carry out this analysis, we have measured *merger trees* for each subhalo from the simulations. These trees allow us to specify all the progenitors (or the descendants) of a substructure at each epoch. To build the merger trees, we use a slightly modified version of the code described in Springel et al., (2001b). We briefly review the main features of the relevant algorithms in the following section.

### 7.1 Constructing merging trees

Following Springel et al., (2001b) we define a subhalo $`S_\mathrm{B}`$ at redshift $`z_\mathrm{B}`$ to be the *progenitor* of a subhalo $`S_\mathrm{A}`$ at redshift $`z_\mathrm{A}`$ (with $`z_\mathrm{A}<z_\mathrm{B}`$) if more than one half of the $`N_{\mathrm{link}}`$ most bound particles belonging to $`S_\mathrm{B}`$ end up in $`S_\mathrm{A}`$. Springel et al., (2001b) adopted $`N_{\mathrm{link}}=10`$. However, we obtained considerably better results with a value of $`N_{\mathrm{link}}`$ that varies between $`10`$ for the less massive substructures to $`100`$ for more massive ones. In this way, occasional failures of the code when building the merger trees were more robustly avoided. Particularly if substructures undergo major mergers, the code identified in some cases the wrong progenitor for $`N_{\mathrm{link}}=10`$, or lost track of an entire subhalo. We have also updated the code so that volatile links to ‘evanescent’ substructures (i.e. objects close to the resolution limit that occasionally appear and then disappear) are avoided.

With these choices, we manage to trace $`85`$$`90`$ per cent of all the substructures in our simulations back to the point when they were first accreted. This fraction goes up to $`87`$$`93`$ per cent if we only consider substructures with more than $`100`$ particles.

### 7.2 Mass accretion history

We now use our merging trees to study the *mass accretion histories* of the subhalos in our simulations. Van den Bosch, (2002) has carried out a similar analysis for dark matter haloes and has proposed an analytic expression for the mass accretion function based on the extended Press-Schechter formalism (Press & Schechter, 1974; Bond et al., 1991; Bower, 1991). This function was found to be in excellent agreement with the results of high-resolution $`N`$-body simulations. Our aim is to study the mass accretion function for subhalos and study whether there is any dependence on mass or on environment.

We have selected subhalos at redshift $`z=0`$ in two different mass ranges ($`10^{11}h^1\mathrm{M}_{}`$ and $`10^{12}h^1\mathrm{M}_{}`$). In order to test for the effects of environment, we selected on one hand subhalos that reside within the virial radius of the massive clusters that formed in simulations B$`1`$ and B$`2`$ by the present day, and on the other hand subhalos located within the smaller haloes found in simulation M$`3`$. In the following, we will refer to the substructures in the clusters as *cluster subhalos* and to the substructures inside the smaller haloes of $`M`$3 as *field subhalos*. Note that since we have excluded from our analysis the main subhalo associated with the FOF group and since on average the most massive substructure has a mass a few per cent of $`M_{200}`$ (see Fig. 4) we end up with very few substructures selected from M$`3`$. In particular we have only $`5`$ substructures with a mass $`10^{12}h^1\mathrm{M}_{}`$ and $`38`$ substructures with a mass $`10^{11}h^1\mathrm{M}_{}`$. The corresponding numbers for the substructures selected from simulations B$`1`$ and B$`2`$ are $`62`$ and $`338`$.

For each of these samples, we build the mass accretion histories as follows: we start from a particular subhalo at redshift $`z=0`$ and construct its merger tree as described in the previous section. At each redshift we track the history of the selected subhalo by linking it with its most massive progenitor.

In Fig. 9 we show a typical example of a mass accretion history for a subhalo with mass $`2\times 10^{11}h^1\mathrm{M}_{}`$ at redshift zero. In the lower panel, we show the mass accretion history of the subhalo and in the upper panel, the corresponding mass of the halo in which the subhalo resides at each redshift.

In this example, the subhalo was accreted onto a larger halo at redshift $`1`$ (shown as a dotted line on the plot). For times prior to this event the substructure was a main subhalo, i.e. the subhalo corresponding to a FOF group, and its mass grew monotonically in time. From now on, we will refer to this event as the *accretion time* ($`t_{\mathrm{accr}}`$) of the subhalo. A few snapshots later, at redshift $`0.8`$, the substructure and its host halo were accreted onto the main progenitor of the cluster (shown as a solid line on the plot).

After the subhalo is accreted, it suffers significant tidal stripping and decreases in mass. In this particular example, the final mass of the subhalo is $`40`$ per cent of the value at its accretion time.

We find that for $`60`$ per cent of the subhalos in the $`10^{11}`$ mass bin and $`80`$ per cent of the subhalos in the $`10^{12}`$ bin the accretion event corresponds to the accretion of the substructure onto the main progenitor of the cluster itself. For most of the rest, the time elapsed between these two events is fairly short. The results we will show later are essentially unchanged if we adopt as definition of the accretion time, the accretion of the substructure onto the main progenitor of the cluster itself.

Fig. 10 shows the distribution of the accretion redshifts for the cluster substructures in our sample. Interestingly, we find that a large fraction of the substructures are accreted at redshift $`z<1`$. As noted above, for most of these substructures this accretion event corresponds to the accretion onto the main cluster itself. Our results hint that substructures are constantly erased in the cluster, being replenished by newly infalling haloes.

In Fig. 11, we plot the distribution of $`M(t=t_0)/M(t=t_{\mathrm{accr}})`$, i.e. the ratio of the mass of the subhalo at the present day to the mass it had when it was first accreted. The histograms show that this ratio has a quite broad distribution, varying from a value of $`1`$ to $`0.1`$. We note that most of the subhalos that have lost only small amounts of mass have been accreted very recently.

This is more clearly shown in Fig. 12 where we plot the average mass accretion function for the cluster subhalos in the two mass bins considered. Three different accretion redshift intervals are considered and in all cases the subhalo masses are normalised to the mass of the subhalo at $`t_{\mathrm{accr}}`$. The thick solid line shows the mean relation for subhalos with mass $`10^{11}h^1\mathrm{M}_{}`$, while the thin line shows the relation for $`M_{\mathrm{sub}}10^{12}h^1\mathrm{M}_{}`$. The mass accretion function monotonically increases for times prior to the accretion event. Once the substructures are accreted, tidal stripping is effective and operates on short time scales. The longer the substructure spends in a more massive halo, the larger is the destructive effect of tidal stripping. Substructures remaining at $`z=0`$ that were accreted at redshift larger than $`1`$ (panel c) have been typically stripped of $`80`$ per cent of their mass. There is also a slight indication that stripping is more effective for more massive substructures: panels (a) and (b) show that more massive substructures have been stripped significantly more than less massive substructures accreted at the same redshift. This effect does not appear in panel (c) but note that we have very few objects accreted in this redshift bin in our sample of more massive subhalos.

In Fig. 13 we compare the mass accretion histories of field and cluster subhalos. We limit the analysis to substructures with mass $`10^{11}h^1\mathrm{M}_{}`$. Again the mass accretion function is normalised to the mass of the subhalo at the accretion time. We find that field subhalos and cluster subhalos have remarkably similar histories suggesting that the efficiency of the tidal stripping is largely independent of the mass of the parent halo.

### 7.3 Merging histories

In hierarchical models of galaxy formation, galaxies and their associated dark matter haloes form hierarchically through merger and accretion processes. In this context, the term *merger* is usually used to refer to an interaction between two objects of similar mass, while the term *accretion* is used to describe the infall of small objects onto much more massive haloes.

Observational results and numerical simulations confirm that interactions (such as tidal truncation or collisions) play an important role in the evolution of galaxies. There is, for example, solid evidence (Schweizer & Seitzer, 1992; Whitmore et al., 1997; Barnes, 1999) that at least some elliptical galaxies are the result of mergers between disk galaxies of similar mass. Mergers may also have a strong effect on the baryonic component of galaxies; they can trigger bar-instabilities in stellar disks and cause an inflow of gas into the galaxies centres, fuelling AGN activity or starbursts.

Following the merger tree of the substructures in our sample we can analyse in detail the merger history, distinguishing between mergers and accretion events. To build the merger history we proceed as follows: we consider all the substructures within the virial radius at redshift zero and follow their merger trees back in time, checking as substructures are accreted onto the main progenitor. We count as mergers all accretion events involving haloes with mass larger than $`2\times 10^{10}h^1\mathrm{M}_{}`$ and mass ratio smaller than $`5:\mathrm{\hspace{0.17em}1}`$. Note that the lower limit on the mass of the merging haloes correspond to the resolution limit of our cluster simulations (see Table 1).

In order to have enough information without running into numerical resolution effects, we limit the present analysis to subhalos with mass $`10^{12}h^1\mathrm{M}_{}`$. Fig. 14 shows the mean number of mergers per subhalo, identified at redshift $`z=0`$, occurring after the redshift plotted in the abscissa. As in the previous section, the sample is split into three different accretion redshift intervals.

Merger events are less frequent once a halo is accreted onto a more massive structure. This is because the merging efficiency is higher in environments where the relative velocities of subhalos are similar to their internal virial velocities. Once haloes are accreted by a massive halo, merging is suppressed by the large velocity dispersion they aquire. This effect is clearly visible in Fig. 14 in the change in slope of the curves near the accretion redshift. Note however that when one goes to significantly higher redshift the difference between the three curves vanishes.

In Fig. 15 we show the fraction of subhalos that have had at least one merger after the redshift plotted in the abscissa. Note that the merger events we are considering are characterised by similar masses and will most likely influence the morphology of the main galaxy, leading to the formation of a bulge component. The final morphology of the galaxies residing in these substructures will depend on the time between the last major merger and accretion onto the cluster: the longer this time, the larger is the likelihood that the galaxy can grow a new disk.

The results in Fig. 15 show that $`80`$ per cent of the $`10^{12}h^1\mathrm{M}_{}`$ haloes in our $`z=0`$ sample have had at least one major merger at redshift below $`2`$; this fraction decreases to $`40`$ per cent for redshift $`<1`$. Surprisingly the fraction is almost independent of the accretion time. These results suggest that a large fraction of subhalos in our cluster sample will host early-type galaxies with a significant bulge component.

## 8 Summary and discussion

We have used a set of high resolution numerical simulations in a $`\mathrm{\Lambda }`$CDM Universe to study dark matter halo substructures. Such dark matter substructures mark the sites where luminous galaxies are expected to be found, so the analysis of their mass functions, radial distributions, merging and mass accretion histories should help us to better understand the properties of the galaxies that form in hierarchical galaxy formation models. Comparison with observational data should then suggest the physics that needs to be included in viable models of galaxy formation and evolution.

In agreement with previous work (Moore et al., 1999b), we find that the shape of the subhalo mass function is almost independent of the mass of the parent halo, with galactic haloes being essentially scaled versions of cluster haloes.

We find that the average mass of the largest substructure within the virial radius (excluding the BCG) scales linearly with the mass of the parent halo. If the stellar masses scale linearly with the dark matter mass of the parent substructure, the second ranked galaxies should have K-band luminosities that increase roughly linearly with the mass of the main halo and are equal to those of the third ranked galaxies to within $`0.5`$ mag.

Note that the assumption that the stellar mass scales linearly with the dark matter mass of the parent substructure cannot generally be true since stars are typically much more concentrated than dark matter. The relation between stellar mass and the mass of parent substructure is then quite complex and should be followed considering the details of star formation and feedback process as is done for example in Springel et al., (2001b).

We have also used the simulations to study the radial distribution of substructures. In agreement with previous work (Ghigna et al., 2000), we find that the radial profile of substructure number density is ‘anti-biased’ relative to the dark matter profile in the inner regions of haloes. The most massive substructures reside preferentially in the outer regions of haloes. This is, at least in part, because substructures undergo substantial tidal stripping in the dense inner regions of haloes.

We have studied the evolution with time of this stripping process and find that once a subhalo is accreted by a larger system, tidal stripping is highly effective; the longer a substructure spends in a more massive halo, the larger is the destructive effect. This suggests that substructures are constantly erased in clusters, being replenished by newly infalling galaxies.

Interestingly, a significant fraction of the substructures found in present-day clusters were first accreted at redshifts $`z<1`$, implying that tidal truncation of the dark haloes of cluster galaxies happened relatively recently, and that at $`z>1`$ many present-day ‘passive’ cluster members should still have been central galaxies of their own extended dark haloes. If gas was able to cool in these haloes, the galaxies may have been considerably more active at $`z1`$ than at present. Substructures in smaller haloes have histories remarkably similar to those in cluster haloes, suggesting that the efficiency of tidal stripping is largely independent of the mass of the main halo.

Once a substructure is accreted onto a cluster, its merging probability decreases because of the large velocity dispersion of the system. Observational data are currently being collected on merger rates in different environments (van Dokkum et al., 1999; Patton et al., 1997). In future work, we plan to carry out a more direct comparison with observational data of this kind, and so to test the hierarchical paradigm for galaxy formation in new ways.

## Acknowledgements

The simulations presented in this paper were carried out on the T3E supercomputer at the Computing Center of the Max-Planck-Society in Garching, Germany and on the IBM SP2 at CINES in Montpellier, France.
We thank Antonaldo Diaferio, Gao Liang and the referee, Fabio Governato, for useful comments and discussions that significantly improved the presentation.
G. D. L. thanks the Alexander von Humboldt Foundation, the Federal Ministry of Education and Research, and the Programme for Investment in the Future (ZIP) of the German Government for financial support.


text:
# Explaining UXOR variability with self-shadowed disks

## 1 Introduction

UX Orionis objects (UXORs, see e.g. reviews by Thé 1994 and Waters & Waelkens 1998) are mostly intermediate mass pre-main-sequence stars displaying a peculiar kind of photometric variability: their V-band light curves are characterized by sudden drops in brightness of up to 3 magnitudes with durations of days to many weeks. These events are separated by relatively long periods of persistence. During these so-called ‘Algol-type minima’ (the name refers to the similarity of the rapid decrease in brightness to those exhibited by eclipsing binaries such as Algol) their spectrum is reddened, and is accompanied by a strong increase in polarization (Grinin et al. 1991). In very deep minima the star often becomes bluer again (e.g. Bibo & Thé 1991). The origin of these brightness drops and the bluing effect has been debated for a long time. The currently favored view is that variations of the column density of dust in the line of sight to the star is to be held responsible. This idea was first put forward by Wenzel (1969) and has since been worked out in more detail by many authors (e.g. Grinin 1988; Voshchinnikov 1990; Grinin et al. 1998). If these obscuration events are due to localized dust clumps of filaments passing through the line of sight, then all phenomenological properties are reproduced: the initial reddening due to dust absorption, and the increase of polarization and later bluing due to the unobscured scattered radiation. While this picture is attractive, it does not explain the nature and origin of these dust clumps or filaments.

One of the leading theories is that Herbig Ae/Be stars are surrounded by many large proto-cometary-clouds or cometary bodies (Grady et al. 2000 and references therein). When one of these objects happens to cross the line of sight towards the star, then an absorption event is expected with precisely the characteristics seen in UXOR stars. The absorption initially reddens the stellar light, but for very large extinction the light scattered off other dust clouds starts to dominate, restoring the color again. One of the main problems of this model is making it consistent with the observed infrared SED. Most of these UXORs have rather strong far-IR excess, with $`L_{\mathrm{dust}}/L_{}0.2\mathrm{}0.3`$. This ratio is an indication of the covering fraction of the dust with respect to the star, so that for UXORs at least 20% of the sky is covered with dust grains. If this all is due to comets, their orbits must be almost isotropically distributed, and there must be a very large number of cometary bodies orbiting around the star. Moreover, the chemical composition of gaseous material moving in and out of the line of sight does not appear to resemble that expected in comets (Natta et al. 2000).

Currently a more favored model for the obscuring clouds/filaments seems to be that of a protoplanetary disk seen nearly edge-on (Grinin et al. 1991, 1996; Bertout 2000; Natta & Whitney 2000; Kozlova et al. 2000). If the inclination angle of the disk is so large that the line of sight towards the star skims the surface of the disk, then it is conceivable that hydrodynamic fluctuations of the surface of the disk could cause dust filaments to pass through the line of sight. Since protoplanetary disks are usually assumed to have a flaring geometry (Kenyon & Hartmann 1987) dust filaments held responsible for these obscuration events must reside in the outer regions of the disk (Thé 1994). A problem with this picture is that one does not expect obscuration events of this kind to be on a time scale of days or weeks. Turbulent eddies and filaments in the disk at 100 AU will tend to move across the star on a time scale of multiple years. One would have to invoke turbulent eddies with a size of about 0.01 times the pressure scale height at 100 AU, and which are very compact (hundreds of times the typical densities at those radii) in order to explain these phenomena with the outer edge of the disk.

In this Letter we propose a new version of the nearly-edge-on disk hypothesis, which does not suffer from the problems mentioned above, and which in fact arises quite naturally from self-consistent models of protoplanetary disks around Herbig Ae/Be stars.

## 2 The model

The idea we wish to present in this Letter is inspired by recent models of passive disks around Herbig Ae/Be stars (Dullemond, Dominik & Natta 2001, henceforth DDN01; Dullemond 2002). In these models, dust evaporation by the radiation of the central star has removed the dust inward of the dust evaporation radius. The dusty part of the disk therefore has an inner rim (at around 0.5 to 1 AU from the central star) which is irradiated frontally by the central star, and has a puffed-up geometry (Natta et al. 2001; DDN01). This puffing-up is a result of the fact that the rim is much hotter than the rest of the disk, since the latter is only irradiated under a shallow angle (see e.g. Calvet et al. 1991; Chiang & Goldreich 1997). Bell et al. (1997) showed that viscous dissipation by accretion could also cause the inner regions of the disk to puff up. The puffed-up geometry of the hot inner rim makes it a good candidate for the origin of the obscuring clouds. The time scales are right: a Kepler time scale is a few months, so that it is very well possible that turbulent filaments pass through the line of sight in a matter of days to weeks. Also the densities are right: the inner rim is expected to have a very high density, allowing that even relatively tenuous hydrodynamic perturbations can have optical depths much larger than unity. The fact that this puffed-up inner rim is a good candidate source for obscuring clouds was already suggested by Natta et al. (2001). But they admit that for inclinations necessary for the inner rim to marginally obscure the star, the outer flaring disk must already strongly obscure both the star and the inner rim. The star would then not have been classified as a Herbig Ae/Be star. A pictographical representation of the idea and the corresponding problem with the flaring outer disk is shown in the left panel of Fig. 1.

A solution to this flaring disk problem may lie in recent 2-D models of Herbig Ae/Be star disks (Dullemond 2002, henceforth D02; Dullemond & Dominik in prep., henceforth DD03). In these models the structure and SED of the disk is computed self-consistently by coupling 2-D continuum radiative transfer to the equations of vertical hydrostatic equilibrium. It was shown that, in addition to the flaring disk + inner rim geometry, there can also exist disks that have sunk entirely into the shadow of their own puffed-up inner rim. These “self-shadowed disks” have a similar inner rim structure as the previous models, but they don’t have the flaring outer regions. For these disks one can therefore easily find inclinations in which the line of sight skims the inner rim of the disk, without passing through a flaring outer disk. In this case the idea, that hydrodynamic turbulent filaments from the puffed-up inner rim are the root cause of UXOR variability, may in fact work. An illustration is shown in the right panel of Fig. 1.

Interestingly, this hypothesis makes a prediction for the shape of the SEDs of UXORs. Self-shadowed disks were shown in D02/DD03 to have relatively weak far-IR excess. They may be classified as “Group II” sources, in the classification of Meeus et al. (2001) 111Note that both Meeus et al. Group I and Group II sources are optically visible and hence belong to the Lada (1987) Class II YSOs.. Flaring disks, on the other hand, have a relatively strong far-IR excess, and belong to the “Group I” sources. The hypothesis put forward in this Letter therefore predicts that UXORs belong to Group II. There may be marginal cases in which the flaring is present, but weak enough to allow the observer to look through the flaring part without too strong extinction. But these should be a minority. The majority of UXOR sources should belong to Group II.

Incidentally, we note that emission from polyaromatic hydrocarbons (PAHs), widely observed in Herbig Ae/Be stars, is believed to arise mainly in the part of the surface layer of the disk that is exposed to direct stellar radiation. Therefore, our 2-D models predict that Group II sources should have only weak or non-existent emission of PAHs, while Group I sources may have strong emission of this kind (DD03). Our proposed model for UXOR variability therefore also predicts UXORs to have weak to non-existent PAH emission.

## 3 SEDs of UXORs and other Herbig Ae/Be stars

Guided by the hypothesis put forward in the previous section, we have investigated the correlation between the shape of the infrared energy distribution and the occurence of UXOR phenomena in HAeBes. Following van Boekel et al. (2003), we characterize the infrared energy distribution of HAeBes by two quantities: the ratio of $`L_{\mathrm{NIR}}`$ (the integrated luminosity as derived from broad-band $`J`$, $`H`$, $`K`$, $`L`$ and $`M`$ photometry from the literature) to $`L_{\mathrm{IR}}`$ (the corresponding quantity derived from IRAS 12, 25 and 60 $`\mu `$m photometry), and the (non color-corrected) IRAS - color. Since Meeus et al. Group I sources show an energy distribution closer to a power-law than the more “double-peaked” energy-distributions of their Group II sources, these two groups will naturally separate in a diagram of $`L_{\mathrm{NIR}}`$/$`L_{\mathrm{IR}}`$ versus the IRAS - color. This will be discussed in detail in van Boekel et al. in prep. In Fig. 2 the diagram is shown for a sample of 86 HAeBes.

Using $`JHKLM`$ and IRAS photometry from literature, we have computed $`L_{\mathrm{NIR}}`$/$`L_{\mathrm{IR}}`$ and IRAS - for all probable HAeBes of spectral types B, A, and F from the catalogue of Herbig Ae/Be stars by Thé et al. (1994) for which sufficient photometry was available to compute these quantities. Empirically, we find that the line $`L_{\mathrm{NIR}}`$/$`L_{\mathrm{IR}}>([12][60])+0.9`$ provides the best separation between the sources known from visual inspection of their energy distribution to belong to Group I and those belonging to Group II.

Using this relation, we find that our sample of 86 HAeBes contains 47 Group I sources and 39 Group II sources. If we define UX Orionis stars as stars of spectral type B9 or later (earlier-type stars do not show the UXOR phenomen; Bibo & Thé 1991; van den Ancker et al. 1998) showing optical variations larger than 1 magnitude on time-scales of days to weeks, 18 of the sources in our sample can be classified as UXOR. Of those, 14 are located in the part of the $`L_{\mathrm{NIR}}`$/$`L_{\mathrm{IR}}`$ versus IRAS - diagram occupied by Meeus et al. Group II sources, whereas the other four are close to the line separating groups I and II. If we limit ourselves to those sources for which we know that they are reasonably ‘isolated’ (i.e. they are certain to be the dominant IR source in the IRAS beam), then there are ten UXORs in Group II and one in Group I.

## 4 Discussion and conclusion

On the basis of the findings of Section 3 we conclude that the 2-D disk-model prediction, that UXOR-type phenomena should only occur in self-shadowed disks, i.e. the HAeBes with relatively modest far-infrared excesses, is consistent with the observational data present in literature.

We also note that a new study (Acke et al., in prep.) of all 48 Herbig Ae/Be stars observed spectroscopically by the Infrared Space Observatory finds that, while Group I sources often show strong PAH emission, these spectral features are absent or weak in Group II sources. Moreover, they find that the vast majority of the UXORs do not have strong PAH features in their infrared spectra. This is consistent with the model predictions.

An important issue that we have not addressed so far is the percentage of Group II sources that we expect on the basis of our model to display UXOR-type variability. This question is difficult to answer, since this requires knowledge of the turbulent behavior of the disk. Typically the surface height of the inner rim of such a disk is of the order of $`H_{\mathrm{rim}}0.2R_{\mathrm{rim}}`$. If hydrodynamic fluctuations are of the order of $`\delta H_{\mathrm{rim}}0.1R_{\mathrm{rim}}`$, then one expects about 13% of Group II stars to display UXOR variability. From Fig. 2 it seems that about 33% of the Group II sources in the catalog of Thé et al. 1994) have UXOR variability. However, we note that the historical selection criteria for Herbig Ae/Be stars clearly favour strongly variable stars, such as UXORs. Therefore it could well be that the real fraction of UXORs is smaller than the 33% in our current sample. Clearly a study of the UXOR fraction in an unbiased sample of Herbig Ae/Be stars could place more stringent contraints on the magnitude of the hydrodynamic fluctuations required to explain the fraction of UXORs.

We stress that so far our model is not dependent on details of the dust properties. However, for an exact prediction of the degree of polarization and the detailed shape of color-magnitude diagrams, the scattering properties of the dust grains, in addition to the geometry of the disk, play an important role. In typical 2-D axisymmetric models of the structure of disks of Herbig Ae stars described in D02/DD03 we find that the disk covering fractions are within the range required by Natta & Whitney (2000) to reproduce the color magnitude diagrams and polarization behavior of UXORs. Therefore we expect our model to be able to reproduce this aspect of the observational characteristics of UXORs as well.

Finally, for the model to work, the disk must have a puffed-up inner rim. For Herbig Ae/Be stars this is a natural consequence of the 2-D self-consistent model. For T Tauri stars there might exist such a rim as well, but it is far more difficult to find parameters for which the entire disk lies in the shadow of the inner rim. On the basis of these arguments it is therefore to be expected that T Tauri stars are only rarely seen to be UXOR type stars. There appear to be some indications that this is indeed the case: Herbst et al. (1994) report that the dominant sources of photometric variations seen in T Tauri stars are rotational modulation due to cool spots on the stellar surface, and changes in the excess or veiling continuum. UXOR variations, although present in their sample, appear to be rare, and are limited to stars with spectral type earlier than K0.

We conclude that the explanation for UXOR variability presented in this paper seems to work better than explanations published in the literature so far. It is a natural consequence of self-consistent 2-D models of disks around Herbig Ae/Be stars, and it makes a number of predictions, which seem to be confirmed by observations.

CPD wishes to thank A. Natta and C. Dominik for discussions and suggestions.