text:
# Retro-MACHOs: 𝜋 in the sky?

## 1 Introduction

In the discovery of MACHOs light has shown its power to reveal dark compact objects. In these events the photon from a distant source suffers a very small angular deflection, small enough to make gravitational lensing the relevant mechanism. The bending power of a black hole is not limited, however, to small angles but reaches to $`\pi `$ and odd multiples of $`\pi `$. Illuminated by a powerful point source of light, the black hole therefore will shine back with a series of concentric rings (we call this retrolensing). A waterdrop, too, likewise illuminated, shines back, but for a different reason: the internal reflections experienced by the photons. That returned light shows to the air traveler flying over a fogbank as a “glory”: a rainbow-like halo surrounding the shadow of the plane on the cloud. Each ray of light that contributes to this sensation has suffered its $`\pi `$ deflection in a different waterdrop, therefore the impression of colored rings that the glory makes is impression only, built up in the last analysis out of multitudes of tiny dots of illumination. No one would be so rash as to expect a detectable backscatter from a single water droplet. It is precisely the search for such backscatter from a single black hole (the putative retro-MACHO) that is, however, the topic of this paper.

How search for retro-MACHOs, for “$`\pi `$ in the sky”, especially when dogged by the negative implications of the familiar phrase? To tell observers to go and look everywhere is to them as dismaying as to be told to look nowhere. Fortunately the successful search for MACHOs provides a helpful model, and the records from that search an immediate place of reference. Recent microlensing events appear to directly confirm that a population of stellar-mass black holes exists in our galaxy (Bennett et al., 2001; Mao et al., 2002; Agol et al., 2002). The Sun, in turn, furnishes a powerful nearby source with whose help one can hope to search a nearby region of the galaxy for black holes not otherwise revealed, taking advantage of times when the Earth interposes itself to spare the registering device from the direct glare of the light source. Although the observation of a retro-MACHO is unlikely, it nonetheless affords one of the few ways to directly image nearby black holes, and if ever observed, would make for an impressive confirmation of Einstein’s theory.

## 2 What is a retro-MACHO?

Consider a point source of light, an observer, and a black hole lens, with all three co-linear and the observer located in between the source and black hole. The observer, looking at the image of the source in the black hole, will see (resulting from symmetry) an (infinite) series of concentric rings, corresponding to bending angles of successively higher odd multiples of $`\pi `$. The outermost ring would be from photons suffering a $`\pi `$ deflection, the next inner ring corresponding to a $`3\pi `$ deflection, and so on. For a point source, and in the limit of geometric optics, the amplification of the rings in infinite. For realistic scenarios with extended sources, however, the amplification is finite. For example, consider an extended spherically symmetric source emitting uniformly across its surface. The image of the source at the observer is an annulus, where the outer and inner radii correspond to the appropriate impact parameters for photons coming from the top and bottom parts of the source, as shown in Fig. 2. Since lensing conserves surface brightness, the brightness of the ring image relative to that of the source is given by the ratio of the areas of the images, as seen at the observer.

The fraction of the sky covered by the ring image is given by:

$$A_I=\frac{\pi b_{o}^{}{}_{}{}^{2}\pi b_{i}^{}{}_{}{}^{2}}{4\pi D_{OL}^{}{}_{}{}^{2}},$$
(1)
where $`b_i`$ and $`b_o`$ are the impact parameters corresponding to the inner and outer radii of the ring, respectively. Note that the impact parameter is not the same as the periastron (closest approach) distance:

$$b=\frac{r}{\sqrt{12M/r}},$$
(2)
where $`b`$ is the impact parameter, $`r`$ is the periastron distance, and $`M`$ is the mass of the black hole. Throughout we measure all quantities in geometric units, with $`G=c=1`$ (e.g., $`1M_{\mathrm{}}=1.5\text{km}`$).

This is to be compared with the sky coverage of the source, as seen directly by the observer:

$$A_S=\frac{\pi R_S^2}{4\pi D_{OS}^2}.$$
(3)
The total amplification, $`\mu `$, of the image is given by the ratio of their sky coverages, yielding

$`\mu `$ $`=`$ $`{\displaystyle \frac{b_o^2b_i^2}{R_{S}^{}{}_{}{}^{2}}}\left({\displaystyle \frac{D_{OS}}{D_{OL}}}\right)^2.`$ (4)
To calculate the brightness of the image we need to know the inner and outer radii of the ring, represented by rays from points $`A`$ and $`B`$ in Fig. 1, respectively. Rays from $`A`$ suffer a deflection angle of $`\pi \alpha `$, while for those from $`B`$ the deflection is $`\pi +\alpha `$, where $`\alpha =R_S/D_{LS}`$. In everything that follows, we restrict our attention to the primary (outermost, brightest) images, although a whole succession of images is present (corresponding to higher odd multiples of $`\pi `$). For the sake of simplicity, we confine ourselves to the case of non-rotating (Schwarzschild) black holes. The calculations for rotating (Kerr) holes are significantly more complex (Rauch & Blandford, 1994), and will be deferred to later work. To calculate the impact parameters that correspond to the necessary deflection angles, we take advantage of an approximate expression in Luminet (1979) (see also Chandrasekhar (1983)) for the case of large deflection:

$$b=b_c+b_de^\mathrm{\Theta },$$
(5)
where $`b_c=3\sqrt{3}M`$, $`b_d=648\sqrt{3}e^\pi M(\sqrt{3}1)^2/(\sqrt{3}+1)^23.4823M`$, and $`\mathrm{\Theta }`$ is the deflection angle. Chandrasekhar also provides exact analytic forms for the solution (in terms of Elliptic functions). We find that the approximate expressions are within $`10\%`$ of the exact values for bend angles satisfying $`\mathrm{\Theta }\pi /4`$, and are good to within $`1/2\%`$ for bend angles with $`\mathrm{\Theta }\pi `$.

Using the results of eq. 5 to calculate $`b_o`$ and $`b_i`$, eq. 4 can be simplified to

$$\mu 3.22M^2\frac{D_{OS}^2}{R_SD_{LS}D_{OL}^2}.$$
(6)
Eq. 6 assumes perfect (colinear) alignment of the source, observer, and lens. The generic case, of course, will involve misalignment, as shown in Fig. 2. We define the misalignment angle, $`\beta `$, to be the angle between the source–observer line and the observer–lens line, as measured at the observer (such that $`\beta =0`$ corresponds to the case of perfect alignment shown in Fig. 1). As $`\beta `$ increases from zero the ring

splits into two arcs (on either side of the lens), centered on the lens–observer–source plane (see the inset of Fig. 2), with the arcs shrinking to small, point-like images as $`\beta `$ approaches $`\pi /2`$. The tangential extent of these arcs, as seen at the observer, is determined by the angular extent of the protrusion of the source out of the plane of Fig. 2: $`\mathrm{\Delta }\mathrm{\Theta }2\mathrm{tan}^1\left(R_S/D_{OS}\mathrm{sin}\beta \right)`$. The total area of each image is thus found to be

$`A_{\mathrm{image}}`$ $`=`$ $`\pi (b_o^2b_i^2)\times {\displaystyle \frac{\mathrm{\Delta }\mathrm{\Theta }}{2\pi }}`$ (7)
$`=`$ $`(b_o^2b_i^2)\mathrm{tan}^1\left({\displaystyle \frac{R_S}{D_{OS}\mathrm{sin}\beta }}\right),`$ (8)
where we have neglected the fact that the image heights taper down to zero at their edges. Here $`b_o`$ and $`b_i`$ represent the outer and inner apparent impact parameters of the images, which correspond to bend angles of $`(\pi \delta )\alpha `$ for the bigger image, and $`(\pi +\delta )\alpha `$ for the smaller image. Note that, since surface brightness is conserved, the bigger image is also the brighter one (and is found further from the center of the lens, on the side opposite that of the source). The magnification of each image is again given by the ratio of the areas of the source and image, namely:

$$\mu =(b_o^2b_i^2)\mathrm{tan}^1\left(\frac{R_S}{D_{OS}\mathrm{sin}\beta }\right)\frac{D_{OS}^2}{\pi R_S^2D_{OL}^2}.$$
(9)
In the limit $`\beta 0`$, and summing up both images, eq. 9 goes over to eq. 4, as expected. Using eq. 5, and taking both $`\alpha 1`$ and $`\beta 1`$, the combined magnification of both images simplifies to

$$\mu 3.22M^2\frac{D_{OS}^2}{R_SD_{LS}D_{OL}^2}2.05\beta M^2\frac{D_{OS}^3}{R_S^2D_{LS}D_{OL}^2}.$$
(10)
## 3 Sources for retro-MACHOs

Every source of light on the sky has the potential to be lensed by a retro-MACHO (in the same manner that

the entire Universe is imaged near the horizon of each and every (unobscured) black hole). However, given the very slight magnification, only the very brightest of sources, with the most felicitous of alignments, will be observable.

We envisage the retro-MACHO to sit at one or another distance and ask for its expected luminosity in order to compare it with familiar objects in the sky. In each case we take the primary source of illumination to be the Sun (output $`I=4\times 10^{33}\text{erg}/\text{sec}`$, apparent visual magnitude $`m=26.8^m`$). In Table 1, we show the observed apparent visual magnitudes of these hypothetical solar retro-MACHOs. In addition to perfect alignment ($`\beta =0`$), it is also useful to consider the case of edge alignment ($`\beta =R_S/D_{OS}=R_{\mathrm{}}/1\text{AU}`$), wherein the observer–lens line grazes the edge of the source. Because of its large angular size (as seen from Earth), the Sun offers the best opportunities for close alignment. In Eq. 9 we find that the Sun’s $`1/2\mathrm{°}`$ angular size dominates over smaller values of $`\beta `$ for other sources.

A number of other potential sources, in addition to the Sun, are worth mentioning. Consider Sirius, the next brightest star in the sky. For a perfectly-aligned retro-MACHO of $`1M_{\mathrm{}}`$, at a distance of $`0.01\text{pc}`$, the retro-image of Sirius would appear at $`34^m`$. A misalignment of $`1\mathrm{°}`$ would put its brightness at $`51^m`$. The reg giant Betelgeuse would have a brightness of $`39^m`$ perfectly aligned, and drop to $`53^m`$ when $`\beta =1\mathrm{°}`$. These are to be compared to the Sun, which has a retro-image at $`31^m`$ for perfect alignment. For $`\beta =1\mathrm{°}`$ this drops by only 3 mag, to $`34^m`$, showing the advantage of the Sun’s large angular size.

As the magnification factors in eqs. 6 and 10 scale as $`M^2`$, one might be tempted to consider the supermassive black hole at the galactic center as a retro-MACHO lens. The Sun, at its closest ($`5\mathrm{°}`$) alignment with the galactic center black hole, has a retro-image brightness of $`50^m`$. This becomes even fainter when the considerable extinction to the galactic center (double the usual value, as the photons travel round trip) is taken into account.

It is to be noted that all of the expressions found in this paper become relevant for black holes at cosmological distances by replacing the physical distances with the appropriate cosmological luminosity distances. One can imagine imaging the Milky Way galaxy in the supermassive black hole of a nearby galaxy.

However, even if Andromeda had an unobscured central supermassive black hole, the image of the Milky Way would be dimmer than $`70^m`$. In more extreme cases, we could imagine a supernovae or quasar imaged (in any frequency band, e.g. radio) in a distant supermassive black hole. Even for perfect alignments, such scenarios produce extremely dim images.

## 4 Properties of retro-MACHOs

A retro-MACHO would appear somewhat like a conventional MACHO. However, rather than magnifying a pre-existing star, a retro-MACHO causes a “star” to appear out of nowhere. A number of sample retro-MACHO lightcurves are shown in Fig. 3. The colors of a retro-MACHO are identical to those of the source, and this will provide an important way to confirm the nature of such events. The structure of a retro-MACHO will be unresolvable (much like the multiple images in microlensing are unresolvable), as the radius of the ring is on the order of the Schwarzschild radius of the lens. Even for very close black holes (say, at the edge of the solar system), the angular extent of a retro-MACHO image remains less than a milliarcsec.

As in the MACHO case, the timescale of retro-MACHO events is determined by the relative motion of the source, lens, and observer. However, rather than the Einstein angle of the lens being the relevant angular scale, in the retrolensing case we are interested in the angular extent of the source measured at the observer. We define the timescale of an event to be given by the length of time it takes to go from perfect alignment to edge alignment:

$$t=\frac{D_{OL}}{v}\frac{R_S}{D_{LS}},$$
(11)
where $`v`$ is the transverse velocity of the observer with respect to the source–lens line-of-sight. The Sun’s angular extent is $`1/2\mathrm{°}`$, which means that a solar retro-MACHO lasts about $`1/2`$ day, assuming the source–lens motion is stationary as compared to the motion of the Earth about the Sun. For lens motion at galactic velocities ($`200\text{km}/\text{sec}`$), the timescale for the lens to move an appreciable angle on the sky is considerably longer (unless the lens is within the solar system). For the case of distant stars being retro-lensed by distant black holes, however, the timescales can shorten because of the extremely small angles needed for edge alignment. The retro-image of Sirius in a black hole about a parsec away can be expected to last on the order of hours, with the precise value depending on the relative projected velocities.

Another point of consideration is the light travel-time of a given retro-MACHO event. For example, for a solar retro-MACHO at a distance of a parsec, the event will appear on the sky corresponding to a position directly opposite to where the Sun was located $`2\text{pc}/c6\text{years}`$ ago! For a retro-MACHO to occur within $`1\mathrm{°}`$ of the current antipodal position of the Sun, the black hole must be within about a day’s round-trip light travel time, putting the black hole lens at the distance of Pluto.

In addition, solar retro-MACHO events will repeat annually, as the Earth returns to the appropriate alignment spot in the ecliptic. This annual repetition may become one of the most valuable signatures of retro-MACHOs, and can be used both to identify candidates in sky surveys, and to confirm the nature of the events.

## 5 Observing strategies

From Table 1, it is apparent that retro-MACHOs are exceedingly dim, even in the best of cicumstances. The most likely observed retro-MACHO would be our Sun imaged in a nearby stellar-mass black hole. Simple theoretical estimates give a value for the local number density of stellar-mass black holes of $`N_{\text{BH}}8\times 10^4/\text{pc}^3`$ (Shapiro & Teukolsky, 1983). Based on the microlensing event rates towards the galactic bulge of the MACHO and OGLE teams, we find $`N_{\text{BH}}2\times 10^4/\text{pc}^3`$ (Agol & Kamionkowski, 2001; Alcock et al., 2001). Both these rates include black holes in the mass range $`5M_{\mathrm{}}M_{\text{BH}}15M_{\mathrm{}}`$. Although it is unlikely for there to be a stellar-mass black hole within a parsec, it is reasonable to expect one within $`10\text{pc}`$ of the Sun. Such a black hole, if perfectly aligned with the Sun (i.e. in the ecliptic plane), would produce a $`49^m`$ retro-image of the Sun. Even at a distance of $`1\text{pc}`$, the image only brightens to $`41^m`$. The earliest population of stars are thought to have masses in the 300–1,000$`M_{\mathrm{}}`$ range, and probably collapse to form black holes in a similar mass range (see, e.g., Fryer, Holz, & Hughes (2001)). The retro-MACHOs from these holes would be visible at $`32^m`$ to a distance of over a parsec. It is apparent, however, that even in the best of circumstances only a chance encounter with an unusually nearby or massive black hole will lead to an observable retro-MACHO.

Microlensing surveys have produced a number of candidate black hole lenses. These arise from long-timescale events, where the parallax effect of the Earth’s revolution about the Sun is strong enough to break some of the degeneracies in the microlensing fits. For example, the best fit to MACHO galactic bulge event 96-BLG-5 is a black hole of mass $`6M_{\mathrm{}}`$ at a distance of $`1\text{kpc}`$ (Bennett et al. 2001; see also Agol et al. 2002). Although this is too distant to be observed, it is conceivable that a future microlensing event may provide us with a definite retrolensing candidate. In these cases we know exactly where and when to look, and what to expect.

Although the likelihood is small, a close approach of a stellar-mass black hole would have a potentially catastrophic impact on our solar system (both through disruption of the Oort cloud, and through direct effects on the orbital stability of the planets). Retrolensing offers perhaps the only direct observational signature of such an incoming “rogue” black hole. The Sun can be considered a searchlight, sweeping out a $`1/2\mathrm{°}`$-wide swath of the sky centered on the ecliptic plane in search for black holes. For a magnitude limit of $`\overline{m}`$, this reveals all black holes of mass $`M`$ or greater out to a limiting distance of $`0.02\text{pc}\times \sqrt[3]{10^{(\overline{m}30)/2.5}(M/10M_{\mathrm{}})^2}`$. Through repeated, deep observations of the direction anti-podal to the Sun, retrolensing can be utilized to construct a black hole early warning system.111 It is to be noted that a very powerful ($`>\text{ Gigawatt}`$) laser can be used to actively probe for retro-MACHOs. One could scan the entire sky systematically, and monitor for returning glints of light (in the same frequency band as that of the laser, and with time delays a direct measure of distance to the black holes). This would constitute a true black hole early warning system. Consider Nemesis, a hypothetical solar companion (Hills, 1985). If we take Nemesis to be a $`1M_{\mathrm{}}`$ black hole at the outer edge of the Oort cloud ($`D_{OL}=9\times 10^4\text{AU}`$), and lying near the ecliptic, it would be directly imaged as a retro-MACHO at $`43^m`$. However, such a Nemesis would be revealed ($`m<30^m`$) if it approached to a distance of $`D_{OL}10^3\text{AU}`$. It is to be noted that microlensing, although quite proficient at detecting compact objects in the direction of the galactic bulge at distances of kiloparsecs, is not particularly efficient at detecting nearby black holes. The Einstein angle of a point mass lens is given by $`\theta _E=\sqrt{4MD_{LS}/(D_{OL}D_{OS})}`$, which for distant stars lensed by nearby masses simplifies to $`\theta _E\sqrt{4M/D_{OL}}`$. For a black hole at $`D_{OL}=10^3\text{AU}`$, a star would have to be aligned to within $`1.5\text{arcsec}`$ to be microlensed. For a solar retro-MACHO, on the other hand, edge-alignment is automatically satisfied for any black hole lying within $`1/4\mathrm{°}`$ of the ecliptic. Retrolensing is thus better than microlensing for detecting approaching nearby black holes.

We are entering an age where the entire sky can be expected to be imaged deeply, and all transient objects noted (Blandford, 2001). Within any such large, repetitive sample, a retro-MACHO would look like a MACHO event with zero baseline flux (i.e. no background star). Upon detecting such an event, one would try to identify a possible light source at the diametrically opposite position on the sky (“through” the Earth). Does the source have identical colors to the transient image? If so, retro-MACHO! For example, in the case of the Sun, we would look for images with precisely solar spectra. It may be of interest to perform a thorough search of extant microlensing survey databases for events that do not correspond to amplification of background stars. Presuming a non-detection of such events, and taking the depth of current surveys to be $`22^m`$, we can directly rule out unobscured black holes with mass greater than $`10M_{\mathrm{}}`$ in the direction of the galactic bulge ($`\beta 5\mathrm{°}`$) at distances less than $`150\text{AU}`$. These distance limits improve dramatically with deeper surveys.

Should a retro-MACHO ever be observed, it would afford an unprecedented opportunity to study the strong-field aspects of general relativity. These objects probe very near the horizon of the holes—in the perfect alignment case, the image represents photons that have passed within $`3.5m=1.75r_s`$ of the singularity (where $`r_s`$ is the Schwarzschild radius), as deep into the strong field regime as we are ever likely to probe via light. A retro-MACHO observation would be an impressive confirmation of the theory of general relativity, and would leave little doubt as to the existence of Einstein black holes.

It is a pleasure to acknowledge discussions with Eric Agol, Doug Eardley, Warner Miller, Bohdan Paczyński, James Peebles, Sterl Phinney, Bob Rutledge, and Johndale Solem in the course of this work. We thank Emily Bennett for providing invaluable assistance in the preparation of the manuscript.


format: text

text:
# Microscopic structure of liquid hydrogen: a neutron diffraction experiment

## I Acknowledgments

###### Acknowledgements.

The skillful technical assistance of Mr. Pierre Palleau (ILL, Grenoble, France), during the setup of the experiment, is gratefully acknowledged.


text:
# Radiative accelerations on Ne in the atmospheres of late B stars

## 1 Introduction

Radiatively driven diffusion processes are generally accepted to be responsible for chemical peculiarities observed among upper main sequence CP stars (Michaud, 1970). The upward radiative acceleration and downward gravity acting on different chemical species compete and drive their microscopic diffusion, resulting in the stratification of these elements in the stellar atmospheres or envelopes of these stars. If these atmospheres are sufficiently stable, chemical inhomogeneities can persevere, and are not wiped out by various mixing processes such as convective zones, stellar winds or rotationally induced mixing. In the particular case of HgMn star atmospheres, temperatures are high enough for hydrogen to be largely ionized, thus removing its convection zone. Also, if the He ionization convection zone diminishes due to He settling in late B stars rotating slower than c. 75 km s-1 (Michaud, 1982; Charbonneau & Michaud, 1988), diffusion can operate very effectively in the atmosphere. This is the case for HgMn stars as they are also known as slow rotators. Such diffusion processes are generally time dependent problems, but the characteristic time-scales in the atmosphere are much shorter than in the envelope ($`10^3`$ yr) due to the considerable lower density, smaller collision rates and higher diffusion velocities (c. 1 cm s-1). This makes the time dependent calculations extremely difficult. So far, diffusion can be followed in time only in the envelope, as shown in the case of neon by (1997, 1999); Landstreet, Dolez & Vauclair (1998) and LeBlanc, Michaud & Richer (LeBlanc et al.2000). These calculations suggest that radiative accelerations on Ne are much lower than gravity below the atmosphere and Ne should sink and be observed as underabundant unless a weak stellar wind transports upwards the Ne accumulated in deep layers.

Nevertheless, observations are restricted to the atmospheres. This important stellar region bridges the effects in the envelopes with those outside the stars, such as stellar winds, and one clearly needs also to know the radiative accelerations here in order to compare the theory with observations. The case of Ne is worth studying for two reasons: it is an important element having very high standard abundance (Anders & Grevesse, 1989; Grevesse, Noels & Sauval, Grevesse et al.1996; Dworetsky & Budaj, 2000) of about A(Ne/H)=$`1.23\times 10^4`$, and because one can expect similar effects as in the case of He. Ne i (along with He i) has a very high ionization potential and first excited level well above the ground state. Ne i is thus the dominant ionization state in the atmospheres of these stars and all Ne i resonance lines are in the Lyman continuum. This part of the spectrum is very sensitive to non-LTE (NLTE) effects and the temperature structure of the atmosphere. It is known that there is a big difference in the temperature behaviour between LTE and NLTE model atmospheres at small optical depths and that temperature can rise outwards in the NLTE models. Ne i is also known to have strong NLTE effects on lines originating from excited levels in the visible part of the spectrum even in late B-type main sequence stars (Auer & Mihalas, 1973; Sigut, 1999; Dworetsky & Budaj, 2000). Thus, the main complication is clear in doing the task. One cannot rely upon the diffusion approximation when calculating energy flux as this is a non-local quantity in the atmosphere, in contrast to envelope calculations. One needs to solve the radiative transfer as, e.g., in Hui-Bon-Hoa, LeBlanc & Hauschildt (Hui-Bon-Hoa et al.2000), Hui-Bon-Hoa, Alecian & Artru (Hui-Bon-Hoa et al.1996) or Hui-Bon-Hoa et al. (2002), including NLTE effects if possible.

Unfortunately, observations of Ne in HgMn stars are very scarce and apart from recent LTE Ne abundances in $`\kappa `$ Cnc, HR7245 and HD29647 by Adelman & Pintado (2000) and Adelman et al. (2001) the only study of Ne was undertaken recently by Dworetsky & Budaj (2000). Their NLTE analysis of 20 HgMn and 5 normal stars revealed that Ne is clearly underabundant in HgMn stars and the deficit reaches its peak values in the middle of the $`T_{\mathrm{eff}}`$ domain of HgMn stars.

In this paper we deal with the radiative accelerations of neon in the atmosphere. This might enable us to model the vertical element stratification in the future, e.g. considering the atmosphere as being in a sequence of stationary states which are to be tailored to the much more slowly varying envelope as proposed by Budaj et al. (1993). Nevertheless, some conclusions can be drawn already from the radiative accelerations alone. In the following, cgs units and LTE approximations (in a NLTE atmosphere model) are used if not specified otherwise.

## 2 Radiative acceleration

### 2.1 Theory

The physics and calculations of radiative accelerations in the envelopes are clearly explained in (1997), recent improvements are described in Gonzales et al. (1995) and some theoretical aspects are discussed in (1995). Our approach to the problem of calculating accelerations in the atmospheres is described here. Radiative acceleration on the element is due to the momentum transferred to the element ions from photons absorbed via bound-bound, bound-free and free-free transitions. Bound-bound transitions are usually the most important. An expression for the radiative acceleration on the ion $`i`$ through its spectral lines can be derived from the radiative transfer equation. For a particular spectral line in the case of complete redistribution we have:

$$dI_\nu =[n_{u,i}A_{ul}(n_{l,i}B_{lu}n_{u,i}B_{ul})I_\nu ]h\nu \phi _{lu}(\nu )ds$$
(1)
where the three terms describe spontaneous isotropic emission, absorption, and stimulated emission, respectively. Note that stimulated emission is in the direction of the photon absorbed. Here, and in the following, $`I_\nu `$ is monochromatic intensity, $`F_\nu `$ is monochromatic energy flux, $`h\nu `$ is the energy of the transition from the lower level $`l`$ to the upper level $`u`$, $`\phi _{lu}(\nu ,\tau )`$ is the normalized Voigt profile, $`\tau `$ is the Rosseland optical depth, $`s`$ is the distance along the beam, $`n_{l,i}`$ is the NLTE population of the $`l`$-th state of $`i`$-th ion, $`n_i`$ is the total NLTE $`i`$-th ion population, $`b_{l,i},b_{u,i}`$ are b-factors defined as a ratio of NLTE to LTE level population ($`b_{l,i}n_{l,i}/n_{l,i}^{}`$), $`m`$ is the mass of the $`i`$-th ion, $`c`$ is the speed of light, $`k`$ is the Boltzmann constant, and $`T`$ is the temperature. $`A_{ul},B_{ul},B_{lu}`$ are the Einstein coefficients where $`B_{lu}`$ is related to the oscillator strength $`f_{lu}`$ and $`B_{ul}`$ by:

$$B_{lu}=\pi e^2f_{lu}(m_\mathrm{e}ch\nu )^1,g_{l,i}B_{lu}=g_{u,i}B_{ul}$$
(2)
where $`e,m_\mathrm{e}`$ are the electron charge and mass respectively and $`g_{l,i},g_{u,i}`$ are statistical weights of the lower and upper levels. Introducing the vertical geometrical depth $`dr=\mathrm{cos}\theta ds`$, multiplying Eq. 1 by $`\mathrm{cos}\theta d\omega `$ and integrating through all solid angles $`d\omega `$ we get for the radiation pressure gradient $`dp_\nu ^R/dr`$:

$$c\frac{dp_\nu ^R}{dr}=\left(1\frac{n_{u,i}}{n_{l,i}}\frac{g_{l,i}}{g_{u,i}}\right)n_{l,i}B_{lu}h\nu F_\nu \phi _{lu}(\nu )$$
(3)
If $`p_{lu}`$ is the momentum removed from the radiation field through one spectral line in unit volume and $`a_{lu}`$, $`a_i`$ are the corresponding accelerations acquired by the ion $`i`$ then we get:

$$a_i(\tau )=\underset{lu}{}a_{lu}=\frac{1}{mn_i}\underset{lu}{}\frac{dp_{lu}}{dt}\mathrm{or}$$
(4)
$$a_i(\tau )=\frac{1}{mn_i}\underset{lu}{}_0^{\mathrm{}}\frac{dp_\nu ^R}{dr}𝑑\nu $$
(5)
The sum runs through all $`lu`$ transitions in the ion. Applying the Boltzmann formula to LTE populations, in this general NLTE case, we finally have:

$$a_i(\tau )=\underset{lu}{}\frac{h\nu }{mc}\left(1\frac{b_{u,i}(\tau )}{b_{l,i}(\tau )}e^{h\nu /[kT(\tau )]}\right)B_{lu}\frac{n_{l,i}(\tau )}{n_i(\tau )}$$

$$_0^{\mathrm{}}F_\nu (\tau )\phi _{lu}(\nu ,\tau )𝑑\nu .$$
(6)
Note that, the particle flux of the element associated with such definition of radiative acceleration would be: 111One can define the radiative acceleration in some other way as e.g. (1997), corresponding to different expressions for the particle flux.

$$J=\underset{i}{}n_iv_i=\underset{i}{}n_iD_i\frac{m}{kT}a_i$$
(7)
where $`v_i,D_i`$ are the $`i`$-th ion diffusion velocity and diffusion coefficient, respectively.

The above expression for radiative acceleration includes the correction for stimulated emission which is often omitted by other authors, as pointed out by Budaj (1994) and (1997), but does not include the redistribution effect pointed out by (1976) and generalized to include gravity by Alecian & Vauclair (1983). It has been substantially revisited by Gonzales et al. (1995). The main idea of the redistribution effect is that, after the photon absorption, ion $`i`$ is in the excited state $`u`$ and has non-negligible probability $`(1r_{u,i})`$ of being ionized (mainly by collisions with fast electrons) before losing its momentum (mainly by collisions with protons). The probability of its remaining in the state $`i`$ is $`r_{u,i}`$ as the probability of recombination or further ionization is negligible. Momentum absorbed during the transition in the state $`i`$, $`p_{lu}`$, should then be redistributed to the two mass reservoirs $`mn_{i+1}`$ and $`mn_i`$ of $`i+1`$ and $`i`$ ions in the proportions $`(1r_{u,i})`$ and $`r_{u,i}`$, respectively. Consequently, one gets: 222Note that these formulae are valid for our radiative accelerations, i.e. as defined by Eq. 4. Other formulae for the redistribution effect without $`\frac{n_i}{n_{i+1}}`$ factor would also hold but if combined with a slightly different definition of radiative acceleration, e.g. that of (1997).

$$a_i=\frac{1}{mn_i}\underset{lu}{}\frac{dp_{lu}}{dt}r_{u,i}=\underset{lu}{}a_{lu}r_{u,i}$$
(8)
$$a_{i+1}=\frac{1}{mn_{i+1}}\underset{lu}{}\frac{dp_{lu}}{dt}(1r_{u,i})=\underset{lu}{}a_{lu}(1r_{u,i})\frac{n_i}{n_{i+1}}$$
(9)
The redistribution function, $`r_{u,i}`$, depends on the ionization rate from the particular level, $`\beta _{u,i}`$, and the collision rate of the ion, $`\beta _i`$:

$$r_{u,i}=\frac{\beta _i}{\beta _{u,i}+\beta _i}=\frac{1}{\frac{\beta _{u,i}}{\beta _i}+1}$$
(10)
It is often useful to picture the radiative accelerations on many ions being applied via some effective radiative acceleration on the element as a whole. Consequently, to compare the results with other authors and to display them, we also calculated the weighted mean radiative acceleration (although the results are stored for each ion and abundance separately). We adopted the following expression suggested by Gonzales et al. (1995):

$$a(\tau )=\frac{_in_iD_ia_i}{_in_iD_i}$$
(11)
$`D_i`$ is obtained from $`D_{i,\mathrm{HI}},D_{i,\mathrm{HII}}`$, which are diffusion coefficients which describe the diffusion of the $`i`$-th ion in neutral and ionized hydrogen, respectively:

$$\frac{1}{D_i}=\frac{1}{D_{i,\mathrm{HI}}}+\frac{1}{D_{i,\mathrm{HII}}}$$
(12)
For the diffusion of the $`i`$-th ion (excluding neutral) in protons we used (Aller & Chapman, 1960):

$$D_{i,\mathrm{HII}}=\frac{1.947\times 10^9T^{5/2}}{n_{\mathrm{HII}}M_2Z_i^2A_1(2)}$$
(13)
where $`n_{\mathrm{HII}}`$ is the proton number density, $`Z_i`$ is the ionization degree (0 = neutral, …) and:

$$A_1(2)=\mathrm{ln}(1+x_\mathrm{D}^2);M_2=\frac{A}{1+A}$$

where $`A`$ is atomic weight in atomic mass units and:

$$x_\mathrm{D}=\frac{4d_\mathrm{D}kT}{Z_ie^2}\mathrm{and}d_\mathrm{D}=\sqrt{\frac{kT}{4\pi n_\mathrm{e}e^2}},$$

where $`d_\mathrm{D}`$ is Debye length and $`n_\mathrm{e}`$ is electron number density. For the diffusion of neutrals in protons we used (Michaud, Martel & Ratel, Michaud et al.1978):

$$D_{0,\mathrm{HII}}=3.3\times 10^4\frac{T}{n_{\mathrm{HII}}}\sqrt{\frac{1+A}{\alpha A}}$$
(14)
where $`\alpha `$ is the polarizability of neutral species ($`\alpha =0.395,0.667\times 10^{24}`$ cm3 for Ne i and H i respectively, from Teachout & Pack (1971)). The same formula was used for the diffusion of ions in neutral hydrogen with $`\alpha \alpha (\mathrm{H}\mathrm{I})`$ and $`n_{\mathrm{HII}}n_{\mathrm{HI}}`$, where $`n_{\mathrm{HI}}`$ is the neutral hydrogen number density. For the diffusion of neutrals in neutral hydrogen we used (Landstreet et al., 1998):

$$D_{0,\mathrm{HI}}=\frac{5.44\times 10^3\sqrt{T}}{n_{\mathrm{HI}}(\delta _{\mathrm{HI}}+\delta _i)^2}\sqrt{\frac{1+A}{A}}$$
(15)
where $`\delta _{\mathrm{HI}},\delta _i`$ are H i and diffusing atom diameters, respectively. The following values calculated using the method of Clementi, Raimondi & Reinhardt (Clementi et al.1963)333http://www.webelements.com/webelements.html were adopted: $`\delta _{\mathrm{HI}}=\mathrm{1.06\hspace{0.17em}10}^8\mathrm{cm},\delta _{\mathrm{NeI}}=\mathrm{0.76\hspace{0.17em}10}^8\mathrm{cm}`$.

### 2.2 Line selection and numerical calculations

As one goes from the envelope higher into the atmosphere, the electron density and associated Stark broadening drop sharply and the spectral lines become much narrower. One needs to consider fine structure to integrate properly the contributions from bound-bound transitions. Atomic data for the Ne i-iv lines were extracted from the Kurucz line list (Kurucz, 1990, CDROM 23). However, one usually does not need to consider all the lines of the element and it is of high practical importance to select just the most important ones. Lines which are important for a particular atmosphere model were selected in the following way.

For practical purposes it is convenient to ascribe three characteristic Rosseland optical depths $`\tau _1,\tau _2,\tau _3`$ or characteristic temperatures $`T_1(\tau _1)T_2(\tau _2)T_3(\tau _3)`$ to each ion $`i`$ of the element of interest. The range $`T_1`$ to $`T_3`$ spans the region where the ion $`i`$ is sufficiently populated and its contribution to the total radiative acceleration of the element is thus not negligible, while $`T_2`$ is the approximate temperature at which the ion dominates the element and the $`i`$-th ion opacity is maximal. For each line $`lu`$ of the ion $`i`$ two parameters $`RA_{i,lu,1}`$ and $`RA_{i,lu,3}`$ were calculated which correspond to two characteristic temperatures $`T_1,T_3`$:

$$RA_{i,lu,1,3}=8.853\times 10^{13}m^1\left[1e^{h\nu /(kT_{1,3})}\right]$$

$$f_{lu}g_{l,i}e^{E_{l,i}/(kT_{1,3})}F_\nu (\tau _{1,3})$$
(16)
where $`E_{l,i}`$ is the excitation potential of the lower level and:

$$F_\nu (\tau _{1,3})=\pi B_\nu (T_{\mathrm{eff}})\mathrm{if}T_{1,3}<T_{\mathrm{eff}},\mathrm{or}$$
(17)
$$F_\nu (\tau _{1,3})=\frac{4}{3}\frac{\pi }{\chi _{c,\nu }}\frac{B_\nu (T)}{T}\frac{T}{r}\mathrm{if}T_{1,3}>T_{\mathrm{eff}},$$
(18)
where $`B_\nu `$ is the Planck function, $`T_{\mathrm{eff}}`$ is the effective temperature, and $`\chi _{\mathrm{c},\nu }`$ the approximate continuous opacity after Borsenberger, Michaud & Praderie (Borsenberger et al.1979).

The parameter $`RA_{i,lu,1,3}`$ is essentially a rough estimate of the expected radiative acceleration through weak lines, as Eq.18 takes only opacity in the continuum into account and partition functions are set to 1.444The first simplification does not have any serious effect on this selection of the most important lines of a particular ion as it means we overestimate the expected accelerations through strong lines, so we cannot omit them. The second simplification also does not have any effect on line selection of an ion as $`RA_{i,lu,1,3}`$ or radiative accelerations for all the lines of a particular ion are scaled by the same value of ion population which is proportional to its partition function. Again, it also means we cannot omit important lines. Then we calculate a sum $`RA_{i,1,3}=_{lu}RA_{i,lu,1,3}`$ through all the lines of the $`i`$-th ion and skip the lines which do not contribute significantly to the sum. We ended up with about 103 Ne lines for a particular model atmosphere.

To evaluate numerically the integral in Eq. 6 it is necessary to estimate the integration step and integration boundaries. Thermal Doppler broadening at the coolest layer and classical radiative damping put constraints on the width of any spectral line and we use minimum frequency spacings equal to

$$\mathrm{\Delta }\nu _{\mathrm{step}}(\mathrm{\Delta }\nu _\mathrm{D}(\tau _1)+\mathrm{\Gamma }_\mathrm{R}/4\pi )/2$$
(19)
where: $`\mathrm{\Delta }\nu _\mathrm{D}=\nu c^1\sqrt{2kT_1/m}`$ and $`\mathrm{\Gamma }_\mathrm{R}=2.47\times 10^{22}\nu ^2`$. Integration boundaries are associated with our line width, which is generally the broadest at the place of highest electron density due to Stark broadening and we use:

$$\mathrm{\Delta }\nu _\mathrm{w}=k_1[\mathrm{\Gamma }(\tau _3)/4\pi +\mathrm{\Delta }\nu _\mathrm{D}(\tau _3)]$$
(20)
where

$$\mathrm{\Gamma }(\tau _3)=\mathrm{\Gamma }_\mathrm{R}+\mathrm{\Gamma }_\mathrm{S}(\tau _3)+\mathrm{\Gamma }_\mathrm{W}(\tau _3)$$
(21)
is the sum of the radiative, Stark and Van der Waals damping parameters. Here, $`k_1`$ is an adjustable constant for which numerical tests revealed that $`k_123`$ is high enough. This is, however, just the case of a weak line. If the opacity in the center of our line is much greater than the opacity in the continuum one must usually integrate much further until the line opacity drops significantly below continuum. The reason why is apparent if one expresses the flux in Eq.6 following the diffusion approximation and neglects continuum opacity; it is inversely proportional to the Voigt function. Consequently, the integrand in Eq.6 is a constant function of frequency. Only when the integration is carried far enough into the line wing, where continuum opacity again dominates, does flux approach a constant function of frequency and the integrand sharply fall. Finally, the integration boundaries should be the maximum value of the mentioned effects:

$$\mathrm{\Delta }\nu _{\mathrm{width}}=\pm \mathrm{MAX}[\mathrm{\Delta }\nu _\mathrm{w}(\tau _3);\mathrm{\Delta }\nu _\mathrm{s}(\tau _2)]$$
(22)
where $`\mathrm{\Delta }\nu _s(\tau _2)`$ is defined for strong lines only and is $`\mathrm{\Delta }\nu `$ for which the line opacity drops significantly below the continuum opacity at the depth where the involved ion population reaches its peak values. Typically, anything from 40 to 103 frequency points were necessary to calculate the radiative acceleration through a Ne line.

The models adopted here were NLTE model atmospheres calculated with tlusty195 (Hubeny, 1988; Hubeny & Lanz, 1992, 1995) for standard solar composition and zero microturbulence. We treated H i, He i, He ii, C i, C ii, C iii as explicit ions, which means that their level populations were solved in NLTE under the assumption of statistical equilibrium and their bound-bound, bound-free, and free-free opacities were taken into account. Carbon was included as it is an important opacity source at cooler effective temperatures (Hubeny, 1981). In Fig. 1 we plot the NLTE model against the Kurucz (1993) LTE line blanketed model for the same effective temperature and gravity. The two vertical bars point to the area where the continua around Ne i $`\lambda `$736 and Ne i $`\lambda `$6402 lines are formed. The first is a strong resonance line in the Lyman continuum and the second is in the Paschen continuum in the visible. Abundance analyses of Ne have been performed using $`\lambda `$6402 (Dworetsky & Budaj, 2000). The Lyman continuum is very opaque and models which go very high in the atmosphere are needed to solve properly the radiative transfer at these wavelengths.

$`F_\nu (\tau )`$ in Eq. 6 was calculated using the synspec36 and synspec42 (Hubeny, Lanz & Jeffery, 1994) codes to solve the radiative transfer modified to get a nonstandard output of flux at each frequency and depth (I. Hubeny, priv. comm.), and to take into account the elements with atomic number $`>30`$ (Krtička, 1998). This enabled us to take into account the bound-free and free-free opacity from explicit levels specified in the atmosphere model and the line blending using the full line list of Kurucz (1990). All the Kurucz lines with line-to-continuum opacity ratio greater than $`10^8`$ at $`\tau =\mathrm{MAX}[1,\tau _2]`$ were considered. All the resonance lines of Ne and other elements were calculated in approximate NLTE (allowing the source function to deviate from $`B_\nu `$ following the ‘second-order escape probability method’ of (1984), see Eq. 25).

Radiative accelerations were then calculated following Eq. 6 assuming $`b_{u,i}=b_{l,i}=1`$ and LTE level populations. Partition functions were taken from the uclsyn code (Smith & Dworetsky, 1988; Smith, 1992). Radiative accelerations were calculated for Ne i-iv and four abundances: $`50\times A_{\mathrm{}},A_{\mathrm{}},10^2\times A_{\mathrm{}},10^5A_{\mathrm{}}`$ where $`A_{\mathrm{}}=1.23\times 10^4`$ and for the following model atmosphere parameters: $`T_{\mathrm{eff}}=\mathrm{11\hspace{0.17em}000},\mathrm{12\hspace{0.17em}000},\mathrm{13\hspace{0.17em}000},\mathrm{14\hspace{0.17em}000},\mathrm{15\hspace{0.17em}000}`$ K and $`\mathrm{log}g=4`$, plus models with $`T_{\mathrm{eff}}=\mathrm{12\hspace{0.17em}000}`$ K and $`\mathrm{log}g=3.5,4.5`$. These are available in digital form in Budaj, Dworetsky & Smalley (2002) including a fortran77 code containing the partition function routines.

## 3 Results and discussion

### 3.1 Radiative accelerations on different Ne ions

The LTE populations of various Ne ions throughout the atmosphere are depicted in Fig. 2 for a representative atmosphere model with $`T_{\mathrm{eff}}=\mathrm{12\hspace{0.17em}000}`$ K, $`\mathrm{log}g=4`$. Ne i is definitely dominant in the atmosphere. Ne iv becomes populated at the base of our models while Ne ii may dominate at very low densities and electron concentrations for $`dm<10^4`$ g cm-2. But above $`dm<10^3`$g cm-2 strong departures of Ne ionization and excitation equilibrium from LTE occur (Sigut, 1999) and our accelerations on different Ne ions should be considered as approximate above this level (see Fig.9 for an illustration of expected deviations from fully consistent NLTE case).

Accelerations in the representative model are plotted in Fig. 3. One can see that none of them exceeds the gravitational acceleration except at the very deepest layers. The acceleration on Ne i is strongest in the atmosphere where this ion is also most populated. The contribution of the neutral species like Ne i to the total effective radiative acceleration of the element is both favoured by weighting with their diffusion coefficient (see Eq. 11) which is approximately two orders of magnitude larger than that for the charged species and reduced by the redistribution effect. Nevertheless, acceleration of Ne i will be the most important contribution to the total acceleration not only in the atmosphere but also at rather great depth up to $`dm10`$ g cm-2. The low radiative accelerations found on Ne i are a consequence of its atomic structure because all resonance transitions, which are usually the most important, are in the Lyman continuum where Ne i atoms see little photon flux, while the rest of the lines originate from highly excited and weakly populated states. Nor can the redistribution effect discussed earlier increase the Ne I acceleration (see Eq. 8 and Fig. 7); just the opposite occurs. Radiative acceleration gained in the Ne i state is redistributed between Ne i and Ne ii, which reduces the Ne i acceleration and adds it to the Ne ii state where, however, the acceleration is not so effective due to the much lower diffusion coefficient of ionized species; see Eqs. 7, 11.

We would like to stress that radiative acceleration (Eq. 6) on a particular ion is very sensitive to the partition functions. These are generally precise in the region where the ion dominates the element but have large uncertainties outside the region. Displaying effective accelerations overcomes the problem as partition functions drop out in Eq. 11 outside the region. Thus, if the reader wishes to use our radiative accelerations stored separately for each ion, the same partition functions should be used also.

### 3.2 Dependence on effective temperature and surface gravity

Now we can proceed further and explore how the situation changes with the effective temperature of the star. It is clear from Fig. 4 that the total acceleration rapidly increases with the effective temperature. However, within our range of interest it never exceeds the gravity in the line forming region. Nevertheless, this suggests that Ne could accumulate at deeper layers in cooler stars than in hotter ones, the latter thus being more vulnerable to departures from the ideal stable atmosphere. Various mixing processes or radiatively driven stellar winds whose intensity also increases with the effective temperature (Babel, 1995; Krtička & Kubát, 2001) might more easily enrich Ne in the atmospheres of much hotter stars. This resembles the case of He which has atomic structure similar to Ne. The abundance of He progressively increases with effective temperature from He-weak towards He-rich stars. This suggests that it might be interesting to search for Ne-rich analogues of He-rich stars among early B stars. On the other hand, for $`T_{\mathrm{eff}}<\mathrm{11\hspace{0.17em}000}`$ K, hydrogen becomes partially ionized, which triggers ambipolar diffusion and hydrogen superficial convective zones. This may also transport some Ne to the atmosphere from below. Our observations (Dworetsky & Budaj, 2000) indicate that Ne underabundances are most pronounced in the middle of the HgMn domain and tend to be less extreme towards the cool and hot ends of the HgMn temperature region. This pattern could therefore be qualitatively expected. Observing Ne below $`T_{\mathrm{eff}}=\mathrm{11\hspace{0.17em}000}`$ K, however, becomes very difficult.

The radiative accelerations also depend on the surface gravity (Fig. 5). Lowering the surface gravity of the model shifts the radiative acceleration behaviour to deeper layers and reduces the radiative acceleration at the bottom of the atmosphere. However, the gravity to radiative acceleration ratio does not change very much and in this context effective temperature is a more important atmospheric parameter than $`\mathrm{log}g`$. Apparently, within the main sequence (approximately $`3.5\mathrm{log}g4.5`$), radiative acceleration of Ne never overcomes gravity.

### 3.3 Effects of homogeneous Ne abundance

Radiative acceleration on the element also depends on its abundance. This is a simple consequence of the fact that the energy flux in the integral in Eq. 6 depends on the element abundance. This is, however, only important in the case of strong lines where the line opacity is comparable to or larger than opacity in the continuum. One can then expect that radiative acceleration will asymptotically increase as element abundance decreases, until it reaches some maximum value corresponding to the sufficiently low abundance when all the element lines disappear and radiative acceleration will no longer be sensitive to the abundance (for more detail see Alecian & LeBlanc, 2000). The situation is illustrated in Fig. 6. Surprisingly, the accelerations are very little dependent on abundance in the photospheres above $`dm0.1`$ g cm-2 and do not exceed gravity for any abundance except at the extreme base of our models. Here a Ne deficit ranging from -0.5 dex for the hottest to -2 dex for the coolest model could be supported by radiation (i.e. radiative acceleration roughly equals gravity for this abundance and depth, $`dm20`$ g cm-2). This weak abundance sensitivity results from the fact that resonance Ne i lines are not very strong because they are in the Lyman continuum where the continuous opacity is very large, effectively reducing line to opacity ratio. We have found that contributions to the radiative acceleration of lines originating from excited Ne i levels is also important, but these weak lines do not introduce strong dependence on the abundance because their line-to-continuum opacity ratio is low due to the high excitation energies and low populations of the levels from which they originate. One can see from the figure that there is no abundance for which the radiative acceleration could balance the gravity in the line forming region above $`dm0.1`$ g cm-2. The implication is that in this region, assuming a stable atmosphere devoid of any motions, Ne should sink and should be almost completely depleted and only its strong concentration gradient could balance the downward flux of Ne atoms. The Ne atom is at least an order of magnitude heavier than the mean ‘molecular’ mass of the gas, so the characteristic scale of the Ne abundance gradient would be an order of magnitude less than the pressure scale height, i.e. $`\mathrm{log}dm0.1`$. Consequently, the fact that Ne is detected in most HgMn stars (Dworetsky & Budaj, 2000) suggests the presence of some kind of instability competing with diffusion and suppressing the efficiency of diffusion processes. An example of such an instability is the weak stellar wind suggested by Landstreet et al. (1998).

One may also notice a sudden drop in the radiative accelerations at approximately $`10^4dm10^3`$ g cm-2. This is not a numerical artifact but a very interesting effect when radiative acceleration may acquire negative values; this will be discussed in more detail below.

### 3.4 Redistribution effect

To take the redistribution effect into account properly one would need to calculate probabilities of all possible radiative and collisional transitions from the upper level to continuum including multiple cascade transitions as suggested by Gonzales et al. (1995). That is a very difficult task and the above authors developed a method in which they set $`r_{u,i}=0`$ for $`n3`$ and $`r_{u,i}=1`$ for lower states. Thus their $`r_{u,i}`$ is a simple step function of level energy and does not depend on e.g. temperature. To cope with the problem and estimate the influence of the redistribution effect we propose a slightly different and more general approach. The idea is simple. As mentioned in Section 2.1, collisions with electrons may change the ionization state before the momentum absorbed by the atom is lost in collisions with protons. Consequently, the redistribution function from a particular excited level depends mainly on the electron to proton collisional rates ratio. Higher atomic levels with energies corresponding to the electron kinetic energy $`kT`$ below the ionization threshold are most strongly coupled to the continuum via collisional excitation and ionization. Probabilities of such processes are strongly temperature dependent via the Boltzmann factor $`e^{\mathrm{\Delta }E/kT}`$ (, 1962; Van Regemorter, 1962) and still deeper atomic levels are affected as one goes deeper into the star. Electron to proton number density ratio almost does not change with depth in such hot stars and the Boltzmann factor thus embraces the essence of the electron to proton collision rates ratio and one may write $`\beta _{u,i}/\beta _i=Ce^{\mathrm{\Delta }E/BkT}`$. Consequently, based on Eq. 10 we suggest the following expression for the $`r_{u,i}`$:

$$r_{u,i}=\frac{1}{Ce^U+1}\mathrm{with}U=(I_iE_{u,i})\left(\frac{1}{A}+\frac{1}{BkT}\right)$$
(23)
where $`I_i`$ is the $`i`$-th ion ionization potential and $`A,B,C`$ are adjustable constants. It is beyond the scope of the present paper to calibrate precise values for the above mentioned constants. $`B`$ is most important as it parameterises the temperature dependence and is of the order of 1. The $`1/A`$ term was introduced for those who prefer to use a Gonzales et al. (1995) type of redistribution but do not like its step function shape and $`A`$ is of the order of $`kT`$, but we use $`A\mathrm{}`$. $`C`$ parameterises the ratio of the total electron to proton collision rates. If $`E_{u,i}I_i`$ the level will register all electron impacts and $`\beta _{u,i}/\beta _i=C`$. Because electrons are $`\sqrt{m_p/m_e}`$ faster than protons, their collisions will be more frequent by approximately the same factor, and C is thus of the order of 40. The method of Gonzales et al. (1995) is equivalent to the approximation of Eq.23 by a temperature independent rectangle and could be approached as a special case when $`B\mathrm{}`$ and $`A0`$ or $`A\mathrm{}`$ depending on whether $`E_{u,i}>I_iE_{crit}`$ or $`E_{u,i}<I_iE_{crit}`$, respectively, where $`E_{crit}`$ is some threshold energy. Note that $`A0`$ means $`r_{u,i}1`$ i.e. no redistribution at all and $`A\mathrm{}`$ means $`r_{u,i}1/(C+1)`$ a small value, i.e. total redistribution is approached but never fully realised. Assuming that Gonzales et al. (1995) took the lowest $`n=3`$ level of C ii as a cutoff in their calculation for $`T=31000`$K which corresponds to $`E_{crit}=80000`$ cm-1 below continuum, we can simulate their method (using our Eq. 89 and forcing $`r_{u,i}=1`$ or $`r_{u,i}=0`$) and calibrate our method so that it would give $`r_{u,i}=0.5`$ at the same temperature. Assuming $`C=40,A\mathrm{}`$ we get from Eq. 23 that $`B=1.0`$ – a reasonable value indeed. Results of both methods are illustrated in Fig. 7 and one can see that while both methods may give similar results for $`T=31000`$K the Gonzales et al. (1995) method strongly overestimates the redistribution effect for lower temperatures. A simple test of the redistribution effect treatment is recommended by calculating radiative accelerations without weighting them by diffusion coefficients (setting $`D_i=1`$ for example in Eq. 11). This is because such accelerations are an invariant of the redistribution process.

### 3.5 Comparison with other work

We calculated an atmosphere model for $`T_{\mathrm{eff}}=\mathrm{11\hspace{0.17em}530}`$ K, $`\mathrm{log}g=4.43`$, corresponding to one of the envelope models of Landstreet et al. (1998). Then we calculated the radiative accelerations, but in the case of Ne i we omitted all except the resonance lines, to simulate the assumption adopted by Landstreet et al.. This assumption is not really justified as the omitted lines contribute significantly to the radiative acceleration. Inclusion of Ne i lines originating from excited levels, which absorb at the wavelengths where a star radiates much more than in the Lyman continuum, may increase the acceleration on Ne i by more than 1.5 dex. However, this does not have any serious consequences on the conclusions of Landstreet et al., as the radiative accelerations are so small they will still remain below the gravity. The comparison is shown in Fig. 8. In this case we also considered $`D_i=1`$ in Eq. 11 to follow their method. The agreement is quite good and small departures at the base of the atmosphere are mainly caused by slightly different temperature behaviour between our model atmosphere and their envelope model. The curious gap in our curve at $`10^4dm10^3`$ g cm-2 appears because the acceleration is negative here and cannot be plotted on a logarithmic scale.

### 3.6 NLTE effects

The aforementioned negative acceleration is a very interesting NLTE effect which may lower still further the radiative accelerations on some ions in this part of the atmosphere. It is a consequence of the temperature inversion in the model atmosphere and the atomic structure of Ne i. As shown in Fig. 1, temperatures in the outer layers rise in a NLTE model for $`10^4dm10^3`$ g cm-2. At this depth, and in these stars, the atmosphere is still optically thick in the Lyman continuum. Thus, according to the diffusion approximation, the star radiates into itself towards cooler regions at these wavelengths. If an element happens to have its important transitions in the Lyman continuum it may be pushed into the star by the radiation and its acceleration is thus negative. Similar effects can be expected also in the cores of very strong lines not necessarily in the Lyman continuum. To explore NLTE effects on the results we also calculated radiative accelerations on Ne i for three different assumptions in a representative NLTE model atmosphere with $`T_{\mathrm{eff}}=\mathrm{15\hspace{0.17em}000}\mathrm{K},\mathrm{log}g=4`$ and zero microturbulence:

1. LTE case: neon populations in Eq. 6 are in LTE and the neon lines’ source functions for flux calculations is equal to the Planck function

$$S_\nu =B_\nu ;$$
(24)
2. approximate NLTE: the same, but the line source function for resonance transitions is allowed to deviate from the Planck function following second order escape probability methods

$$S_\nu =\sqrt{\frac{ϵ}{ϵ+(1ϵ)K_2(\tau )}}B_\nu $$
(25)
where $`ϵ`$ is the photon destruction probability and $`K_2(\tau )`$ is the kernel function (see Hubeny et al., 1994, the synspec manual, for more detail). Note that this option was used in the calculations described in Section 2.2 above;
3. fully consistent NLTE case: with NLTE model atmosphere, NLTE populations and NLTE source functions

$$S_\nu =\frac{2h\nu ^3}{c^2}\left(\frac{g_{u,i}n_{l,i}}{g_{l,i}n_{u,i}}1\right)^1.$$
(26)
In the full NLTE case the NLTE model atmosphere and Ne i level populations were calculated with tlusty195. As an input for tlusty195 we considered Hubeny’s H i atom model with 9 explicit levels and continuum, the 31 level Ne i atom model taken from Dworetsky & Budaj (2000), and a simple 4 level atom of Ne ii with continuum (Table 1). Generally, the data for Ne bound-bound and bound-free transitions were taken from the TOPbase (Cunto et al., 1993; Hibbert & Scott, 1994) but Ne i oscillator strengths were from (1998). The modion IDL interface written by Varosi et al. (1995) was particularly useful in constructing Ne i and Ne ii atom models. Populations of those very high Ne i levels not included in our Ne i atom model were considered in LTE relative to the ground state of Ne ii, of which the population was calculated in NLTE, when calculating the energy flux and acceleration. All Ne i level populations were set to LTE below $`dm>0.3`$ g cm-2 when calculating the radiative accelerations.

The differences between the three assumptions are exhibited in Fig. 9 and become crucial for $`dm<10^3`$ g cm-2, where the radiative accelerations in LTE and NLTE may differ by more than 3 dex. In the case of LTE we observe negative acceleration in the region of the temperature inversion. This was explained above. The approximate NLTE case is parallel to the LTE acceleration but about 1.5 dex smaller. This is caused by the difference in computed energy fluxes. While flux in the approximate NLTE case is real - lines are in absorption - and surprisingly close to the flux in the full NLTE case, the resonance Ne i lines are in emission in LTE as soon as the atmosphere becomes optically thin in the Lyman continuum and acceleration soars immediately. The difference between approximate NLTE and full NLTE originates mainly in the differences in level populations.

## 4 Summary

We have calculated radiative accelerations on Ne i-iv ions in the atmospheres of late B main sequence stars. We take into account the fine structure and calculations include the effects of line blending using the whole Kurucz line list, bound-free and free-free opacity of explicit elements (H, He, C) as well as some NLTE effects. We explored how the acceleration changes with respect to effective temperature, surface gravity and homogeneous Ne abundance and found that it is much smaller than the gravity in and above the observed line forming region. Only at the base of our models ($`dm20`$ g cm-2) could a Ne deficit of about 0.5-2.0 dex be supported by the radiation, depending on the effective temperature. This implies that, in the stable atmospheres of late B stars, Ne should be almost completely depleted from the photosphere. This is qualitatively in agreement with the observations of Dworetsky & Budaj (2000) but the fact that we detected Ne in some HgMn stars suggests the presence of some weak transport or mixing mechanism contaminating observed layers by Ne from the reservoir underneath. Possibly, this mechanism is associated with the partial hydrogen ionization at the cool end and a very weak stellar wind at the hot end of the HgMn $`T_{\mathrm{eff}}`$ span.

Finally, we demonstrated by a full NLTE calculation the importance of such effects on radiative accelerations and slightly improved a current treatment of the redistribution effect which should be used in future calculations.

## Acknowledgments

We thank I. Hubeny and J. Krtička for patient discussions and help with tlusty and synspec codes, B. Smalley for his work on fitting partition functions for uclsyn, J.D. Landstreet and G. Alecian for discussing details about their calculations and an anonymous referee for many useful comments. J.B. gratefully acknowledges the support of the Royal Society/NATO fellowship scheme (Ref. 98B) and partial support by the VEGA grant No. 7107 from the Slovak Academy of Sciences and APVT-51-000802 project.


text:
# Nonlinear clustering during the cosmic Dark Ages and its effect on the 21-cm background from minihalos

## 1 Introduction

Despite striking progress over the past decade, cosmologists have still failed to see into the “Dark Ages” of cosmic time. Although the details of this epoch between recombination at redshift $`z10^3`$ and reionization at $`z6`$ are crucial for understanding issues ranging from early structure formation to the process of reionization itself, no direct observations of any kind have been made in this redshift range. Recently, we made the first proposal for such direct observations (Iliev et al. 2002, from hereafter Paper I) based on collisional excitation of the hydrogen 21-cm line in the warm, dense, neutral gas in virialized minihalos (halos with virial temperature $`T_{\mathrm{vir}}10^4\mathrm{K}`$), the first baryonic structures to emerge in the standard CDM universe. We showed that collisional excitation is sufficient to increase the spin temperature of hydrogen atoms inside minihalos above the temperature the cosmic microwave background (CMB). Minihalos should generally appear in emission with respect to the CMB, producing a background “21-cm forest” of redshifted emission lines, well-separated in frequency. Unlike all previous works (Hogan & Rees, 1979; Scott & Rees, 1990; Subramanian & Padmanabhan, 1993; Madau, Meiksen & Rees, 1997; Shaver et al., 1999; Tozzi et al., 2000), this mechanism does not require sources of Ly$`\alpha `$ radiation for “pumping” the 21-cm line and decoupling it from the CMB.

In Paper I we calculated the 21-cm emission properties of individual minihalos in detail, along with the total background and its large-scale fluctuations due to clustering. Although individual lines and the overall background are too weak to be readily detected, the background fluctuations should be measurable on $`10^{}100^{}`$ scales with the currently-planned radio arrays LOFAR and SKA. We demonstrated that such observations can be used to probe the details of reionization as well as measure the power spectrum of density fluctuations at far smaller scales than have been constrained previously.

In the current paper we extend these results to smaller scales. Our previous calculations showed that the fluctuation signal increases as the beam size and frequency bandwidth decrease, which corresponds to sampling the 21-cm emission from minihalos within smaller volumes. As these volumes correspond in turn to length scales that are more nonlinear, nonlinear effects have the greatest impact on the angular scales and frequency bandwidths at which the signal is the strongest.

This issue is of particular importance as our previous investigations relied on the standard, simplified nonlinear bias of Mo & White (1996), which breaks down at many of the relevant redshifts and length scales. For example, an rms density fluctuation in a Cold Dark Matter universe with a cosmological constant ($`\mathrm{\Lambda }`$CDM) at $`z=6(8),`$ is strongly nonlinear ($`\sigma (M)1`$) for $`M10^9M_{}`$ ($`6.2\times 10^7M_{}`$), which roughly corresponds to the region sampled by beam sizes 9” (3”), for frequency bandwidths 12 kHz (3 kHz), respectively. Nonlinear effects can influence the predicted background fluctuations on even larger scales, up to few comoving Mpc (corresponding to few hundred kHz frequency bandwidths and 1-10 arc min beams), even if these scales are still not strongly nonlinear. Additionally, rare halos are always more strongly clustered than the underlying density distribution (i.e. the bias is $`>1`$), again bringing nonlinear issues to the fore.

In order to address these important issues we have carried out a series of high-redshift N-body simulations of small scale structure formation, which we use to construct 21-cm line radio maps that illustrate the expected fluctuations in the emission. As no current simulations are able to span the full dynamic range relevant to 21-cm emission, however, we extend our results by developing a new formalism for calculating the nonlinear Eulerian bias of halos, which is based on the Lagrangian bias formalism of Scannapieco & Barkana (2002, hereafter SB02). We describe this approach in detail in this paper, verify it by comparing it with the results of our N-body simulations, and apply it to to calculate improved predictions for the fluctuations in the 21-cm emission. These are then compared to the predictions given in Paper I, quantifying the impact of nonlinear effects on minihalo emission.

The structure of this work is as follows. In § 2 we describe a set of simulations of minihalo emission at high redshift, and use these to construct simulated maps at small angular scales. In § 3 we present our improved calculation of the Eulerian bias and verify it by comparison with the numerical simulations presented in § 2. In § 4 we modify the calculation of the radiation background from minihalos to incorporate the contribution due to nonlinear clustering of sources, according to the formalism described in § 3, and describe the results of our nonlinear formalism. Conclusions are given in § 5.

## 2 Numerical simulations and simulated 21-cm radio maps

We simulated the formation of minihalos in three cubic computational volumes of comoving size 1 Mpc, 0.5 Mpc, and 0.25 Mpc, respectively. We used a standard Particle-Particle/Particle-Mesh (P3M) algorithm (Hockney & Eastwood 1981), with $`128^3`$ equal-mass particles, a $`256^3`$ PM grid, and a softening length of 0.3 grid spacing. Here and throughout this paper we consider a flat $`\mathrm{\Lambda }`$CDM model with density parameter $`\mathrm{\Omega }_0=0.3`$, cosmological constant $`\lambda _0=0.7`$, Hubble constant $`H_0=70\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, baryon density parameter $`\mathrm{\Omega }_b=0.02h^2`$ (where $`h=H_0/100\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$), and no tilt. The initial conditions where generated using the transfer function of Bardeen et al. (1986) with the normalization of Bunn & White (1997) (for details, see Martel & Matzner 2000, §2). All simulations started at redshift $`z=50`$ and terminated at redshift $`z=9`$. To identify minihalos, we used a standard friends-of-friends algorithm with a linking length equal to 0.25 times the mean particle spacing. We rejected halos composed of 20 particles or less. Table 1 lists the comoving size $`L_{\mathrm{box}}`$ of the box, the comoving value of the softening length $`\eta `$, the total mass $`M_{\mathrm{tot}}`$ inside the box, and the particle mass $`M_{\mathrm{part}}`$.

Once the halos are identified, we can use their properties and distribution to compute the corresponding radio maps. We assign to each halo a 21-cm flux based on its mass and redshift of formation by modeling the halos as Truncated Isothermal Spheres (Shapiro, Iliev & Raga, 1999; Iliev & Shapiro, 2001). The 21-cm fluxes from individual minihalos are obtained by solving the radiative transfer equation self-consistently through each halo to obtain the line-integrated flux, as described in Paper I. Let us consider a cylindrical volume with radius $`R=\mathrm{\Delta }\theta _{\mathrm{beam}}(1+z)D_A(z)/2`$, and length $`L_{\mathrm{box}}`$ , where $`D_A(z)`$ is the angular diameter distance, $`\mathrm{\Delta }\theta _{\mathrm{beam}}`$ is the beam size of the observation. Then the beam-averaged differential antenna temperature $`\overline{\delta T}_b`$ is calculated using equation 6 of Paper I, which can be discretized as follows:

$$\overline{\delta T}_b(𝐱)=\frac{c(1+z)^4}{2\nu _0H(z)}\frac{1}{\pi R^2L_{\mathrm{box}}}\underset{i}{}(\mathrm{\Delta }\nu _{\mathrm{eff}})_i(\delta T_{b,\nu _0})_iA_ie^{(𝐱𝐱_i)^2/2R^2},$$
(1)
where $`H(z)`$ is the Hubble constant at redshift $`z`$, $`\nu _0`$ is the rest-frame line frequency, $`\pi R^2L_{\mathrm{box}}`$ is the comoving volume of the beam $`(\mathrm{\Delta }\nu _{\mathrm{eff}})_i`$, $`(\delta T_{b,\nu _0})_i`$, and $`A_i`$ are the effective line width, the line-center differential brightness temperature, and the geometric cross-section of halo $`i`$, respectively, and we use Gaussian filter in order to ensure that the contribution of a particular minihalo to a particular pixel varies smoothly with the location of the minihalo.

The resulting maps are shown in Figure 1. Due to the small box sizes that are required to resolve the minihalos, the beams used to produce the maps are also very small, ranging from 4” to 0.25”. The fluxes from such small beam sizes are well below the sensitivity limits of the currently planned radio arrays LOFAR and SKA. Additionally, the larger-box (1 Mpc and 0.5 Mpc) simulations do not have sufficient mass resolution to resolve the smallest halos, while the larger-mass minihalos are not present in the smaller box simulations. Therefore the flux levels are somewhat underestimated in all cases, and should be considered only as illustrative of fluctuations on very small scales.

## 3 Nonlinear clustering regime and comparison to simulations

### 3.1 An analytic approach to nonlinear clustering

As our simulations are unable to span a sufficient range of scales to both resolve all minihalos and reproduce the physical scales accessible to observations, we adopt instead an alternate technique to estimate the full nonlinear signal. We apply an approximate formalism that is an extension of the one developed in SB02. A more complete discussion of this method, refinement of the basic formalism, and detailed comparisons with simulations will be presented in a separate paper. In this section we summarize the basic features of this approach and show that it is sufficient for our purposes in this investigation. For an alternative approach see Catelan, Matarese, & Porciani (1998).

In SB02 the authors extended the Press-Schechter formalism as reinterpreted by Bond et al. (1991) to calculate the number density of collapsed objects in two regions of space initially separated by a fixed comoving distance. In this approach virialized halos are associated with linear density peaks that fall above a critical value usually taken to be $`\delta _c=1.686.`$ Here the linear overdensity field is defined as $`\delta (𝐱)D^1(z)\rho (𝐱,z)/\overline{\rho }(z)1`$, where $`\rho (𝐱)`$ is the linear mass density as a function of position and and $`\overline{\rho }`$ is the mean mass density. Note that here the linear overdensity field $`\delta (𝐱)`$ is described in terms of its value extrapolated to $`z=0`$, and its evolution with time is subsumed by the “linear growth factor” $`D(z)=\delta (0)/\delta (z).`$ Other choices for $`\delta _c`$ and its evolution have been explored in the literature, by fitting to simulations (Sheth, Mo, & Tormen, 2001; Jenkins et al., 2001) or incorporating the Zel’dovich or adhesion approximations (Lee & Shandarin, 1999; Menci, 2001). While these methods improve the accuracy of the Press-Schechter technique, they provide no information about nonlinear clustering. Here we strive to develop a formalism that is applicable to nonlinear clustering regardless of the recipe chosen for $`\delta _c`$.

Once $`\delta _c`$ determined, the presence of a halo can be associated with a random walk procedure. At a given redshift we consider the smoothed density in a region around a point in space. We begin by averaging over a large mass scale $`M`$, or equivalently, by including only small comoving wavenumbers $`k`$. We then lower $`M,`$ adding $`k`$ modes, and adjusting $`\delta (M)`$ accordingly. This amounts to a random walk in which the “time variable” is the variance associated with the filter mass and the “spatial variable” is the overdensity itself. This walk continues until the averaged overdensity crosses $`\delta _cD^1(z)`$ and we assume that the point belongs to a halo with a mass $`M`$ corresponding to this filter scale.

The problem of halo collapse at two positions in space can be similarly associated with two random walks in the presence of a barrier of height $`\delta _cD^1(z)`$. These random walks are extremely correlated at large scales, and become less and less so as the modes of smaller and smaller scales are added. In SB02 this process was approximated by a completely correlated random walk down to an overall variance $`S`$ = $`\xi `$, where $`\xi `$ is the correlation function that accounts for the excess probability of finding a second halo at a fixed distance from a given halo. This correlated random walk was then followed by uncorrelated walks down to the final variances associated with the masses of the objects, $`S=\sigma ^2(M_1)`$ and $`\sigma ^2(M_2).`$

This procedure lends itself directly to calculating the overall increase in number density due to infall. The fact that the overdensity on large scales will contract the distances between halos provides an additional contribution at the length scale $`x`$ such that $`\sigma ^2(x)=\xi `$, the scale down to which the evolution of the two points is completely correlated. Let $`\delta _0`$ be the linear overdensity at redshift zero at that scale. The probability distribution of $`\delta _0`$, $`Q_0`$ is then simply given by a Gaussian, with a barrier imposed at $`\nu (z)\delta _cD^1(z)`$:

$$Q_0(\nu ,\delta _0,\xi )=[G(\delta _0,\xi )G(2\nu \delta _0,\xi )]\theta (\nu \delta _0),$$
(2)
where $`G(\delta ,\xi )(2\pi \xi )^{1/2}e^{\delta ^2/2\xi }`$, and $`\theta `$ is the Heaviside step function.

We now consider $`f_2(\nu (z),\sigma _1^2,\sigma _2^2,\delta _0,\xi )d\sigma _1^2d\sigma _2^2`$, the joint probability of having point $`A`$ in a halo with mass corresponding to the range $`\sigma _1^2`$ to $`\sigma _1^2+d\sigma _1^2`$ and point $`B`$ in a halo with mass corresponding to the range $`\sigma _2^2`$ to $`\sigma _2^2+d\sigma _2^2`$, whose random walks pass though a value of $`\delta _0`$ at the $`\sigma ^2(M)=\xi `$ scale. In this case we obtain an expression analogous to equation (37) of SB02

$`f_2(\nu (z),\delta _0,\sigma _1^2,\sigma _2^2,\xi (r))`$ $`=`$ $`{\displaystyle \frac{}{\sigma _1^2}}{\displaystyle \frac{}{\sigma _2^2}}\left[2{\displaystyle _{\mathrm{}}^{\nu (z)}}𝑑\delta _1{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\delta _1\right]\left[2{\displaystyle _{\mathrm{}}^{\nu (z)}}𝑑\delta _2{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\delta _2\right]Q_{12}(\nu ,\delta _1,\delta _2,\sigma _1^2,\sigma _2^2,\xi (r)),`$ (3)
where now $`Q_{12}G(\delta _1\delta _0,\sigma _1^2\xi )G(\delta _2\delta _0,\sigma _2^2\xi )`$, and the reader is referred to SB02 for a more detailed derivation of $`f_2`$ from the underlying probability distribution.

Averaging equation (3) over the probability distribution for $`\delta _0`$ as given in equation (2) and carrying out the partial derivatives we obtain the Eulerian “bivariate” number density, the joint probability of having point $`A`$ lie in a halo in the mass range $`M_1`$ to $`M_1+dM_1`$ and point $`B`$ at a comoving distance $`r`$ lie in a halo in the mass range $`M_2`$ to $`M_2+dM_2`$ at a redshift $`z`$:

$$\frac{d^2n_{12,E}^2}{dM_1dM_2}(r,z)=\frac{\overline{\rho }}{M_1}\left|\frac{d\sigma _1^2}{dM_1}\right|\frac{\overline{\rho }}{M_2}\left|\frac{d\sigma _2^2}{dM_2}\right|\left(1\frac{\xi }{\sigma _1^2}\right)\left(1\frac{\xi }{\sigma _2^2}\right)f_{2,E}(\nu (z),\sigma _1^2,\sigma _2^2,\xi (r)),$$
(4)
where

$$f_{2,E}(\nu ,\sigma _1^2,\sigma _2^2,\xi )_{\mathrm{}}^\nu 𝑑\delta _0Q_0(\delta _0,\xi )g(D^1\delta _0)f(\nu \delta _0,\xi \sigma _1^2)g(D^1\delta _0)f(\nu \delta _0,\xi \sigma _2^2),$$
(5)
and $`f(\nu ,\sigma ^2)(2\pi )^{1/2}(\nu /\sigma ^3)e^{\nu ^2/2\sigma ^2}`$. Here $`g(D^1\delta _0)`$ is $`\rho /\overline{\rho }`$ at the scale at which the points $`A`$ and $`B`$ are completely correlated. This can be thought of as the contraction of a large spherical perturbation that surrounds both points and contains a total mass $`M`$ of material, where $`\sigma ^2(M)=\xi `$. This contribution is then just $`(1+D^1\delta _0)`$ in the linear regime, and for our purposes here we assume $`g(D^1\delta _0)=(1+D^1\delta _0)`$ for all values as this reproduces well the results from simulations.

Finally we define the Eulerian bias as

$$b_E^2\xi D^2\frac{d^2n_{12,E}^2}{dM_1dM_2}\left(\frac{dn_{1,E}}{dM_1}\frac{dn_{2,E}}{dM_2}\right)^11$$
(6)
where we divide by $`dn_E/dM`$, a single point number density that accounts for the self correlations between the collapsed peaks and the overdense sphere, which are over-counted in equation (4). To compute this probability we again carry out an average over the probability distribution (2), but in this case we consider $`f_1(\nu ,\delta _0,\sigma _1^2)d\sigma _1^2=f(\nu \delta _0,\sigma ^2\xi )d\sigma _1^2,`$ the probability of having a single point in a halo with a mass corresponding to the range $`\sigma _1^2`$ to $`\sigma _1^2+d\sigma _1^2`$ whose random walk passes through $`\delta _0`$ at the $`\sigma ^2(M)=\xi `$ scale. This gives $`dn_E/dM=(\overline{\rho }/M)|d\sigma ^2/dM|(1+\xi /\sigma ^2)f_E(\nu ,\sigma ^2)`$ where

$$f_E(\nu ,\sigma ^2,\xi )_{\mathrm{}}^\nu 𝑑\delta _0Q_0(\delta _0,\nu ,\xi )g(\delta _0D^1)f(\nu \delta _0,\xi \sigma ^2).$$
(7)
Combining these expressions yields our final expression for the bias

$$b_E^2\xi D^2=\frac{f_{2,E}(\nu ,\sigma _1^2,\sigma _2^2,\xi )}{f_E(\nu ,\sigma _1^2,\xi )f_E(\nu _2,\sigma _2^2,\xi )}1.$$
(8)
Note that at large distances we can work to order $`\delta _0^2`$ to determine the asymptotic limit as $`\xi /\sigma ^20.`$ In this limit

$$g(\delta _0D^1)f(\nu \delta _0,\sigma ^2\xi )\left[1D^1\delta _0\left(\frac{\nu ^2}{\delta _c\sigma ^2}\frac{1}{\delta _c}+1\right)+𝒪(\delta _0^2)\right]f(\nu ,\sigma ^2\xi ),$$
(9)
and the only surviving term in equation (8) is the cross term between terms of order $`\delta _0`$. All other terms cancel out between the numerator and the denominator, giving

$$b_E^2=\left(1+\frac{\nu /\sigma _1^21}{\delta _c}\right)\left(1+\frac{\nu /\sigma _2^21}{\delta _c}\right),$$
(10)
Thus our formalism reproduces the usual bias as in Mo & White (1996) at large distance, as was used in Paper I.

Finally, we define the flux-averaged correlation between objects as

$$\overline{b}_E^2\xi D^2=\frac{𝑑M_1𝑑M_2F(M_1)(dn/dM_1)F(M_2)(dn/dM_2)b_E^2\xi D^2}{\left[𝑑MF(M)(dn/dM)\right]^2},$$
(11)
where $`F(M)`$ is the line-integrated flux of a minihalo of mass $`M`$.

### 3.2 Nonlinear properties and comparisons with simulations

In Fig. 2 we plot $`b_E^2\xi D^2`$ as calculated from equation (8) and the standard expression, equation (10), using the fit to the density fluctuation power spectrum given by Eisenstein & Hu (1999). In each panel we show the behavior of $`10^7M_{}`$ and $`10^6M_{}`$ objects at four different redshifts, which span the detectable range. At large distances, the clustering expressions approach each other asymptotically, as expected from equation (10). On the other hand, the solid lines shoot up dramatically at the smallest distances. This is because the separation between the halos is comparable to their radii, and thus the likelihood of finding a second collapsed halo at the same point becomes infinite as $`r0.`$ At intermediate distances the behavior is more complex. For relatively rare objects, corresponding to high redshifts, the nonlinear values consistently exceed the standard ones. For example, at a distance of 0.1 comoving Mpc, the nonlinear $`10^7M_{}`$ value is almost twice that of the standard result at $`z=20`$, while it is only about 70% of the standard result at $`z=6`$.

To compare this behavior to numerical results, we calculated the correlations between halos for each of our simulations. In this case $`b_E^2\xi D^2`$ is directly comparable to the halo correlation function $`\xi _{i,j}(r)`$, the excess probability of finding a particle within a mass bin $`m_i`$ and a particle within a mass bin $`m_j`$ separated by a distance $`r`$. This expression is symmetric, so that $`\xi _{i,j}(r)=\xi _{j,i}(r)`$. We define $`N_{i,j}(r)`$ as the mean number of particles in a mass bin $`m_j`$, located within a radius $`r`$ of a particle in a mass bin $`m_i`$, averaged over all particles in a mass bin $`M_i`$. This is is obtained by dividing the total number of pairs with separation $`r^{}<r`$, by the number $`N_i`$ of $`i`$ particles,

$$N_{i,j}(r)=\frac{4\pi N_j}{3V_{\mathrm{box}}}\left[r^3+3_0^r\xi _{i,j}(r^{})r^2𝑑r^{}\right],$$
(12)
where $`N_j`$ is the number of particles of mass $`m_j`$ and $`V_{\mathrm{box}}`$ is the volume. Finally, we differentiate equation (12) to get $`\xi _{i,j}(r)`$.

Using these definitions we examined the distribution of halos in each of our simulations at $`z=15`$ and $`z=9`$. In each case we divided the halos into mass bins containing roughly the same number ($`600`$ ) of objects and compared their correlations to those given by equation (8). In the left panels of Fig. 3 we compare the mass averaged correlation between halos,

$$\overline{\xi }_{i,j}=\left(_{m_i}\frac{dn}{dM_i}_{m_j}\frac{dn}{dM_j}b_E^2\xi D^2\right)\left(_{m_i}\frac{dn}{dM_i}_{m_j}\frac{dn}{dM_j}\right)^1,$$
(13)
to the numerical value, $`\xi _{i,j}`$, for objects in three mass bins in our 1 $`\mathrm{Mpc}^3`$ simulation at $`z=15`$ (Fig. 3). For all mass ranges, there is a good match between the numerical $`\xi _{i,i}`$ and our analytic expression, with both values equaling or exceeding the linear predictions at most distances. The only case in which the predictions of the simulations fall below equation (10), is at large distances, where the separations between halos begin to approach the simulation box size, and some damping is to be expected. Comparisons with our smaller box simulations yielded comparable results, but with larger damping at these distances, indicating that the analytic and simulated correlations are a good match at all reliable separations. On the other hand, the cross-correlation functions are a poorer match, particularly if the mass bins are very different. For all correlations however, our nonlinear estimate falls between the standard and simulated results, indicating that the correlations between halos exceed those predicted by equation (10), and that our approach gives a conservative estimate of this excess.

In the right panels of Fig. 3, we repeat this comparison with our 1 $`\mathrm{Mpc}^3`$ simulation at $`z=9`$. In this case, we divided our halos into five mass bins, three of which are shown in this figure. In this case there is good agreement between the correlations and cross-correlations obtained from the simulations and our nonlinear estimates, for all mass scales considered. In fact the analytical estimates are consistent with the results of the simulations for all distances that are small relative to the box, apart from the dip at small distances in the $`m_1m_1`$ case, which is most likely due to small-number statistics. Note that unlike the higher redshift case, our nonlinear estimates now fall below the standard analytic approach at larger distances, as we saw in the lower redshift cases in Fig. 2.

## 4 Fluctuations of the 21-cm emission from minihalos: effects of the nonlinear bias

The amplitude of $`q`$-$`\sigma `$ (i.e. $`q`$ times the rms value) angular fluctuations in the differential brightness temperature $`\delta T_b`$ (or, equivalently, of the flux) in the linear regime are given by

$$\delta T_b^2^{1/2}=qb(z)\sigma _p\overline{\delta T_b}$$
(14)
(Paper I), where $`b(z)`$ is the mean flux-weighted bias, and

$`\sigma _p^2={\displaystyle \frac{8D^2(z)}{\pi ^2L^2R^2}}{\displaystyle _0^{\mathrm{}}}𝑑k{\displaystyle _0^1}𝑑x{\displaystyle \frac{\mathrm{sin}^2(kLx/2)J_1^2[kR(1x^2)^{1/2}]}{x^2(1x^2)}}(1+fx^2)^2{\displaystyle \frac{P(k)}{k^2}},`$ (15)
where $`P(k)`$ is the linear power spectrum at $`z=0`$, and the factor $`(1+fx^2)^2`$, where $`f[\mathrm{\Omega }(z)]^{0.6}`$ is the correction to the cylinder length for the departure from Hubble expansion due to peculiar velocities (Kaiser, 1987). In order to apply the formalism developed in §3, we use that the power spectrum is the Fourier transform of the correlation function $`\xi (r)`$, obtaining

$$\sigma _p^2=\frac{32D^2(z)}{\pi L^2R^2}_0^{\mathrm{}}𝑑rr^2\xi (r)f(r,R,L),$$
(16)
where

$$f(r,R,L)_0^1𝑑x(1+fx^2)^2_0^{\mathrm{}}𝑑k\frac{\mathrm{sin}(kr)}{kr}\frac{\mathrm{sin}^2(kLx/2)}{x^2}\frac{J_1^2[kR(1x^2)^{1/2}]}{k^2(1x^2)}.$$
(17)
Using equations (14) and (16) the mean squared angular fluctuations become

$$\delta T_b^2=\left(\frac{32}{\pi L^2R^2}\right)^2_0^{\mathrm{}}𝑑rr^2\left[D^2b^2(z)\xi (r)(\overline{\delta T_b})^2\right]f(r,R,L).$$
(18)
We apply the formalism developed in §3 by simply replacing $`\left[D^2b^2(z)\xi (r)(\overline{\delta T_b})^2\right]`$ with the corresponding quantity $`D^2\overline{b_E^2}\xi (r)(\overline{\delta T_b})^2`$ for the mean Eulerian flux-weighted bias calculated in §3. As we showed, at large distances the two expressions are equivalent, while at smaller distances the new formalism better reproduces the numerical results.

We start with deriving the optimal frequency bandwidth for observing the 21-cm emission from minihalos. As long as $`\mathrm{\Delta }\nu _{\mathrm{obs}}`$ and the beam size $`\mathrm{\Delta }\theta _{\mathrm{beam}}`$ are large enough to provide a fair sample of the halo distribution, the mean 21-cm flux is independent of $`\mathrm{\Delta }\nu _{\mathrm{obs}}`$, while the fluctuations grow as $`\mathrm{\Delta }\nu _{\mathrm{obs}}`$ decreases. However the sensitivity of the radio array deteriorates for smaller bandwidths as $`\mathrm{\Delta }\nu _{\mathrm{obs}}^{1/2}`$ (e.g. Shaver et al. 1999; Tozzi et al. 2000). The results for $`3\sigma `$ fluctuations at redshifts $`z=7`$ and $`z=8.5`$ vs. observed frequency bandwidth $`\mathrm{\Delta }\nu _{\mathrm{obs}}`$ for beam sizes $`\mathrm{\Delta }\theta _{\mathrm{beam}}=9^{}`$ and 25’ are shown in Fig. 4, along with the results obtained using the linear bias in order to facilitate comparison. As the bandwidth increases, the integration time required to detect the signal decreases, reaching minimum at about $`\mathrm{\Delta }\nu _{\mathrm{obs}}=2`$ MHz and remaining roughly unchanged thereafter. The fluctuations decrease with increasing bandwidth, however,hence the optimal bandwidth is about $`\mathrm{\Delta }\nu _{\mathrm{obs}}2`$ MHz. We see that at large scales the current results are largely indistinguishable from the linear results, as expected, since at such large scales the bias calculated in § 3 reduces to the linear expression in equation (10). For small values of $`\mathrm{\Delta }\nu _{\mathrm{obs}}`$ the differences between the two results grow to $`10\%`$. However, in the observable range the two results hardly differ.

Illustrative results for $`3\sigma `$ fluctuations at redshifts $`z=8.5`$ and $`z=20`$ vs. the beam size $`\mathrm{\Delta }\theta _{\mathrm{beam}}`$ for fixed observed frequency bandwidth $`\mathrm{\Delta }\nu _{\mathrm{obs}}=2`$ MHz are shown in Fig. 5. The approximate sensitivities for the LOFAR and SKA radio arrays are also shown where appropriate (see Paper I for details), but from here on we concentrate on the differences between the results from the two approaches and the robustness of our original predictions rather than observability, which was discussed in Paper I. At redshifts $`z15`$ using the linear bias gives an overestimate, although a small one compared to the various uncertainties of the calculation, while at higher redshifts using the linear bias gives a small underestimate of the fluctuations.

In Fig. 6 (left panels) we plot the $`3\sigma `$ fluctuations at redshifts $`z=8.5`$ and $`z=20`$ vs. the observer frequency bandwidth $`\mathrm{\Delta }\nu _{\mathrm{obs}}`$ for fixed beam size $`\mathrm{\Delta }\theta _{\mathrm{beam}}=1^{}`$. We have chosen such a small beam size in order to investigate the upper limit of the deviation of the results due to the nonlinear clustering of halos. For small values of both $`\mathrm{\Delta }\nu _{\mathrm{obs}}`$ and $`\mathrm{\Delta }\theta _{beam}`$ the differences are as large as 20-30%. In Fig. 6 (right panels) we plot the $`3\sigma `$ fluctuations for beam sizes $`\mathrm{\Delta }\theta _{\mathrm{beam}}=9^{}`$ and $`25^{}`$ vs. redshift $`z`$ (i. e. the spectrum of fluctuations) for observer frequency bandwidth $`\mathrm{\Delta }\nu _{\mathrm{obs}}=2`$ MHz. Again, the same patterns emerge, where the linear bias approach overestimates the fluctuations at the lower end of the range of considered redshifts by up to 10 %, while underestimating the fluctuations at the higher redshifts, although by a smaller fraction.

## 5 Conclusions

We have developed a new analytical method for estimating the nonlinear bias of halos and used it to calculate the fluctuations of the 21-cm emission from the clustering of high-$`z`$ minihalos. This method is likely to be useful in tackling a much wider range of problems that lie beyond the capabilities of current numerical simulations, and will be refined and further verified in a future publication. The minihalo bias predicted by our method at large scales reproduces the linear bias, as derived in Mo & White (1996), and is much larger at small scales, confirming naive expectations. At intermediate scales, however, its behavior is more complex and both mass- and redshift-dependent. For very rare halos at high-$`z`$ (roughly $`z>15`$), the standard linear bias is lower than our nonlinear prediction at all length scales. At the lower end of the redshift range we consider (down to $`z=6`$) the linear bias is higher at intermediate scales (few hundred kpc to few Mpc comoving) and lower at smaller scales.

We have also compared the predictions of our new method with the results of N-body numerical simulations, which we used both to verify our approach and to produce sample radio maps. Due to the limited dynamical range of our simulations, however, these maps are only illustrative of the behavior of the fluctuations on very small scales, below the sensitivity limits of LOFAR and SKA. On the scales at which the simulations are reliable, we find excellent agreement between these results and our analytical approach. The most significant discrepancies occur in the calculation of cross-correlation functions of very different mass bins at higher redshifts. However, these departures are relatively modest, and in all cases our method reproduces the simulation results significantly better than the linear bias.

Despite these differences, we find our original linear bias predictions for the fluctuations of the 21-cm emission from minihalos to be robust at the scales and frequencies corresponding to observable signals. The prediction using the flux-weighted nonlinear bias never departs from the linear prediction by more than few percent in that range, well within the other uncertainties of the calculation. This robustness is partly accidental, however, and is due to the nonlinear bias varying above and below the linear one depending on the length scale, leading to partial cancellation of the differences when the correlation function is integrated over the length scales. For small observational bandwidths, $`\mathrm{\Delta }\nu _{\mathrm{obs}}`$, there is less cancellation and the differences in the two predictions grow at all beam sizes, $`\mathrm{\Delta }\theta _{\mathrm{beam}}`$. Similarly, if $`\mathrm{\Delta }\theta _{\mathrm{beam}}`$ is small, there is little cancellation, and the discrepancies are larger at all bandwidths. In these cases the linear bias gives an overestimate of the 21-cm emission fluctuations from minihalos at the low end of the redshift range, and an underestimate at high-$`z`$, even at large values of $`\mathrm{\Delta }\nu _{\mathrm{obs}}`$. However, when both the beam size and the frequency bandwidth are large the differences in the two approaches at small scales are diluted, and the resulting 21-cm emission fluctuations become identical. Finally, we predict that the best observational frequency bandwidth for improving the chances for detection is $`\mathrm{\Delta }\nu _{\mathrm{obs}}2`$ MHz. Thus it may be through a such frequency window that astronomers get their first glimpses of the cosmological Dark Ages.

## Acknowledgments

We would like to thank Andrea Ferrara and Rennan Barkana for useful comments and discussions. This work was supported in part by the Research and Training Network “The Physics of the Intergalactic Medium” set up by the European Community under the contract HPRN-CT2000-00126 RG29185, NASA grants NAG5-10825 to PRS and NAG5-10826 to HM, and Texas Advanced Research Program 3658-0624-1999 to PRS and HM. ES has been supported in part by an NSF MPS-DRF fellowship.