text:
# Effective Screening due to Minihalos During the Epoch of Reionization

## 1 Introduction

The spectra of the quasars SDSS 1030+0524 at $`z=6.28`$ (Becker et al., 2002) and SDSS 1044-0125 at $`z=5.73`$ (Djorgovski et al., 2002) feature broad regions of high optical depth, possibly implying the existence of a Gunn-Peterson trough of complete resonant Ly$`\alpha `$ absorption by a neutral intergalactic medium (Gunn & Peterson, 1965). However, the discovery of a Ly$`\alpha `$ emission line from a galaxy at $`z=6.56`$ (Hu et al., 2002) implies that the gas in the vicinity of that galaxy is mostly ionized and that the neutral fraction is highly non-uniform. The appearance of a Gunn-Peterson trough is generally attributed to absorption by a smooth intergalactic medium (hereafter IGM) with a significant neutral fraction, and it indicates that observations are approaching the reionization era.

The reionization of hydrogen is expected to involve several distinct stages. The initial, “pre-overlap” stage \[using the terminology of Gnedin (2000)\] consists of individual ionizing sources turning on and ionizing their surroundings. The radiation from the first galaxies must make its way through the surrounding gas inside the host halo, then through the high-density region which typically surrounds each halo. Once they emerge, the ionization fronts propagate more easily into the low-density voids, leaving behind pockets of neutral, high-density gas. During this period the IGM is a two-phase medium characterized by highly-ionized regions separated from neutral regions by sharp ionization fronts. The main, relatively rapid “overlap” phase of reionization begins when neighboring H II regions start to overlap. Whenever two ionized bubbles are joined, each point inside their common boundary becomes exposed to ionizing photons from both sources. Therefore, the ionizing intensity inside H II regions rises rapidly, allowing those regions to expand into high-density gas. This process leads to a state in which the low-density IGM is highly ionized and ionizing radiation reaches everywhere except for gas located inside self-shielded, high-density clouds. The ionizing intensity continues to grow during the “post-overlap” phase, as an increasing number of ionizing sources becomes visible to every point in the IGM.

The recent observations of strong Ly$`\alpha `$ absorption at $`z6`$ indicate a significant neutral fraction in the IGM (e.g., Fan et al., 2002; Lidz et al., 2002). However, they do not necessarily correspond to the overlap stage of reionization, since even after the end of overlap, variations in the gas density and the ionizing intensity in the IGM may be able to produce such large regions of absorption (Barkana, 2002). When the voids become transparent to ionizing photons, the naive expectation is for the ionizing intensity to increase quickly as multiple ionizing sources begin to reach every point in the IGM. However, such a sudden, homogeneous end to the Gunn-Peterson trough may conflict with the observations, which imply that along one line of sight, overlap occurs at $`z>5.9`$ (Becker et al., 2002), while a relatively large region of strong absorption is found along a second line of sight at $`z=5.3`$ (Djorgovski et al., 2002); furthermore, the recent $`z=6.56`$ galaxy may indicate the existence of a transparent region of the IGM at this much higher redshift (Hu et al., 2002). Thus, the ionizing intensity during post-overlap must have been inhomogeneous and must have increased gradually, with both effects possibly caused by the shadows due to the remaining neutral gas. Much of this gas may have been located inside the dense halos of collapsed objects.

In the popular $`\mathrm{\Lambda }`$CDM cosmology, most of the collapsed gaseous halos prior to reionization had a virial temperature below $`10^4`$K, at which cooling due to atomic transitions is heavily suppressed \[Barkana & Loeb (2001), and references therein\]. The gaseous low-mass halos (hereafter “minihalos”) thus remain intact until photoionization heating by the cosmic UV background photo-evaporates much of this gas back into the IGM during the overlap phase (Barkana & Loeb, 1999). Haiman, Abel, & Madau (2000) noted that these halos may consume a significant number of ionizing photons if they are assumed to be completely photo-ionized throughout the time it takes them to expand back into the IGM. However, Cen (2001) pointed out that as an ionizing front approaches a dense gas halo, it slows down and is eventually halted outside the virial radius, with a shock front entering the gas halo instead. Regardless of the precise dynamics, the process of photoevaporation takes an amount of time which is set by the characteristic sound speed of $`10\mathrm{km}\mathrm{s}^1`$ (corresponding to $`T10^4`$ K). Whether minihalos play a substantial global role during this timespan depends on their ability to screen ionizing sources.

The limited mass resolution of current numerical simulations of reionization (Gnedin, 2000; Razoumov et al., 2002) does not allow them to resolve these minihalos and include their cumulative screening effect; Haiman, Abel, & Madau (2000) made a crude estimate of the covering factor of minihalos and found that it could be significant. In this paper we gauge the importance of absorption by minihalos with a detailed semi-analytic calculation. In §2 we describe our model for the basic halo properties and the statistics of their screening. In §3 we describe our numerical results. Finally, §4 discusses the main conclusions of this work. Throughout the paper we use the $`\mathrm{\Lambda }`$CDM cosmological parameters of $`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, and $`\mathrm{\Omega }_b=0.05`$, for the density parameters of matter, cosmological constant, and baryons, respectively. We also assume a Hubble constant $`h=0.7`$ in units of $`100\text{ km s}^1\text{Mpc}^1`$, and a primordial scale-invariant ($`n=1`$) power spectrum with a value $`\sigma _8=0.8`$ for the root-mean-square amplitude of mass fluctuations in spheres of radius $`8h^1`$ Mpc.

## 2 Theoretical Model

### 2.1 Basic Halo Properties

We begin with a description of our model for the properties of minihalos and of ionizing sources in the context of hierarchical galaxy formation in $`\mathrm{\Lambda }`$CDM. The input for our opacity calculation involves the distribution of halo masses for halos that contain uncooled gas and for halos which host galaxies, as well as knowledge of the density distribution inside each minihalo.

We assume that the abundance of halos is given by the model of Press & Schechter (1974), with the modifications which Sheth & Tormen (1998) and Jenkins et al. (2001) used in order to fit more accurately the mass function seen in numerical simulations. Galaxies form in halos in which gas can accumulate and cool. At high redshift, gas can cool efficiently in halos down to a virial temperature of $`10^4`$ K or a circular velocity of $`V_c16.5\mathrm{km}\mathrm{s}^1`$ with atomic cooling, which we assume to be the dominant cooling mechanism (we consider the effect of H2 cooling in §3). Before reionization the IGM is cold and neutral, and these cooling requirements set the minimum mass for halos which can host galaxies. During reionization, the ionized IGM is heated to a temperature $`T_{\mathrm{IGM}}1`$$`2\times 10^4`$ K, and gas infall into halos is suppressed by the increased pressure force. However, this pressure suppression is not expected to cause an immediate suppression of the cosmic star formation rate, since even after fresh gas infall is halted the gas already in galaxies continues to produce stars, and mergers among already-formed gas disks also trigger star formation. We neglect the suppression effect in this paper, since we consider only the era of reionization and its immediate aftermath.

We use the term ’minihalos’ to denote halos in which gas accretes but is unable to cool. These halos contain virialized dark matter and gas, where the gas does not collapse onto a disk and does not participate in star formation. The maximum total mass of a minihalo is $`1.6\times 10^8M_{}`$ at $`z=6`$, corresponding to the atomic cooling threshold (e.g. Barkana & Loeb, 2001). This implies that the gas in minihalos is almost entirely neutral prior to reionization. The minimum mass is determined by the IGM temperature, which prevents gas from collecting inside the smallest halos. The Jeans mass is the minimum mass-scale of a perturbation for which gravity overcomes pressure, in linear perturbation theory. The Jeans mass is related only to the evolution of perturbations at a given time, while a collapsing halo is the end result of the long evolution history of a perturbation. When the Jeans mass varies during this time, the overall suppression of the growth of perturbations depends on a time-averaged Jeans mass. Gnedin & Hui (1998) showed that the correct time-averaged mass is the filtering mass $`M_F`$, which we adopt as the minimum mass of minihalos. At $`z=6`$, for example, we find $`M_F=4.6\times 10^4M_{}`$ in areas that have not yet been reionized and heated.

If gas collects inside a halo and also cools, it can collapse to high densities and form stars (or possibly also a quasar). The ability of stars to form is determined by gas accretion which, in a hierarchical model of structure formation, is driven by mergers of dark matter halos. In order to determine the lifetime of a typical source, we first determine the time spent by gas (i.e. the gas ‘age’) within a given halo using the average rate of mergers which built up the halo. Based on the extended Press-Schechter formalism (Lacey & Cole, 1993), for a halo of mass $`M`$ at redshift $`z`$, the fraction of the halo mass which by some higher redshift $`z_2`$ had already accumulated in halos with galaxies is

$$F_M(z,z_2)=\mathrm{erfc}\left(\frac{1.69/D(z_2)1.69/D(z)}{\sqrt{2(S(M_{\mathrm{min}}(z_2))S(M))}}\right),$$
(1)
where $`D(z)`$ is the linear growth factor at redshift $`z`$, $`S(M)`$ is the variance on mass scale $`M`$ (defined using the linearly-extrapolated power spectrum at $`z=0`$), and $`M_{\mathrm{min}}(z_2)`$ is the minimum halo mass for hosting a galaxy at $`z_2`$ (determined by cooling as discussed above).

We estimate the total age of gas in the halo as the time since redshift $`z_2`$ where $`F_M(z,z_2)=0.2`$, so that most ($`80\%`$) of the gas in the halo has fallen into galaxies since then. Low-mass halos form out of gas which has recently cooled for the first time, while high-mass halos form out of gas which has already spent previous time inside small galaxies. We emphasize that the age we have defined here is not the formation age of the halo itself, but rather it is an estimate for the total period during which the gas which is incorporated in the halo has participated in star formation. However, the rate of gas infall is not constant, and even within the galaxy itself the gas likely does not form stars at a uniform rate. Indeed, galaxies could go through repeated cycles of a star formation burst followed by feedback squelching, followed by another cycle of cooling, fragmentation and star formation. The details involve complex astrophysics, which are not understood even at low redshift, so we account for the general possibility of bursting sources by adding a parameter $`\zeta 1`$, the duty cycle. When $`\zeta <1`$, compared to $`\zeta =1`$ there are fewer active sources at any given time (by a factor of $`\zeta `$) but each has a larger star formation rate (by a factor of $`1/\zeta `$). In addition, the star formation rate is proportional to the efficiency parameter $`\eta `$, which is the fraction of gas in galaxies which turns into stars.

Specifically, the ionizing source in a halo of mass $`M`$ emits, over its entire lifetime (which at $`z=6`$ is typically $`t_S1\times 10^8(\zeta /0.25)`$ yr), the following number of photons per unit solid angle:

$$\frac{dP}{d\mathrm{\Omega }}=6.4\times 10^{64}\left(\frac{M}{10^8M_{}}\right)\left(\frac{N_{\mathrm{ion}}}{40}\right)\left(\frac{\mathrm{\Omega }_b/\mathrm{\Omega }_m}{0.17}\right).$$
(2)
In this equation, $`N_{\mathrm{ion}}`$ is the overall number of ionizing photons produced per baryon in galactic disks; it includes 4000 ionizing photons produced per baryon in stars \[for a Scalo (1998) stellar IMF and a metallicity equal to $`1/20`$ of the solar value111This low metallicity is expected to roughly characterize $`z6`$ galaxies, which form out of an IGM which may have been enriched by an even earlier generation of galaxies; see, e.g., Barkana & Loeb (2001) for further discussion.\], times an efficiency $`\eta =10\%`$, times an escape fraction $`f_{\mathrm{esc}}^{\mathrm{disk}}=10\%`$ of ionizing photons from the disk itself \[see Barkana & Loeb (2001) and references therein\]. Our choice of typical parameters $`\zeta =0.25`$ and $`\eta =10\%`$ is loosely based on numerical simulations and on the low redshift cosmic star formation rate, respectively, but the exact parameters valid during reionization can only be determined from future observations.

### 2.2 Screening due to Minihalos

At any given redshift we consider a single halo emitting ionizing photons, surrounded by a distribution of mini-halos of various masses at various distances $`R`$ from the center of the source halo. Except where indicated otherwise, we neglect the ionizing radiation from other sources when calculating the transmission of the radiation from a source through the minihalos. This assumption is justified as long as we consider the early stage of reionization when the ionizing background has not yet been established, and as long as we consider a small enough distance that the ionizing intensity is dominated by the central ionizing source.

In order for the radiation to reach a large distance $`R`$, it must overcome two obstacles caused by the minihalos: (i) it must ionize the neutral gas within the minihalos, and (ii) it must overcome the recombination rate of this gas in order to keep it ionized. Our statistical treatment of both quantities is similar, so in general we label the quantity we are considering by $`\nu `$. Thus, for example, $`\nu _1`$ is the sum of the number of hydrogen atoms per unit solid angle encountered by a line of sight going out to radius $`R`$. This sum is then compared to the total photon emission per solid angle $`dP/d\mathrm{\Omega }`$ of the source, and the photons escape out to radius $`R`$ along those directions for which $`dP/d\mathrm{\Omega }>\nu _1`$. Separately we also calculate $`\nu _2`$, the sum of the number of recombinations per second of hydrogen atoms (assuming the gas starts out ionized) per unit solid angle along a line of sight out to radius $`R`$. Photons can only escape along those directions for which $`d^2P/[d\mathrm{\Omega }dt]>\nu _2`$, where the (average) emission rate of photons per solid angle is

$$\frac{d^2P}{d\mathrm{\Omega }dt}=\frac{1}{t_S}\frac{dP}{d\mathrm{\Omega }},$$
(3)
in terms of the source lifetime $`t_S`$.

Both $`\nu _1`$ and $`\nu _2`$ are determined by the density distribution of gas in the minihalos. Numerical simulations of hierarchical halo formation indicate a roughly universal spherically-averaged density profile for dark matter in the resulting halos (Navarro, Frenk, & White 1997, hereafter NFW), though with considerable scatter among different halos. Such high-resolution simulations have not been done with gas. However, in one-dimensional simulations of dark matter and adiabatic gas, we find \[based on the methods of Thoul & Weinberg (1996) and Haiman, Thoul, & Loeb (1996)\] similar final density profiles for the gas and dark matter, where both components form a cusp and the gas profile is only slightly shallower than the dark matter profile. Note that this result is found even though the physics of the two components is rather different, with the dark matter undergoing violent relaxation due to shell crossing, while the gas is halted at a virialization shock which moves outward. For simplicity, we therefore assume that even in the three-dimensional case the virialized gas profile is similar to that of the dark matter, and so we adopt the radial NFW profile for the gas density in minihalos,

$$\rho (r)=\frac{\rho _0}{c_\mathrm{N}x(1+c_\mathrm{N}x)^2},$$
(4)
where $`x=r/r_{\mathrm{vir}}`$, $`c_\mathrm{N}`$ is the concentration parameter, and the virial radius is

$$r_{\mathrm{vir}}=2.1\left(\frac{M}{10^8M_{\mathrm{}}}\right)^{\frac{1}{3}}\left(\frac{\mathrm{\Omega }_m}{0.3}\right)^{\frac{1}{3}}\left(\frac{1+z}{7}\right)^1\left(\frac{h}{0.7}\right)^{\frac{2}{3}}\mathrm{kpc},$$
(5)
where in this equation and in what follows we assume a sufficiently high redshift for the Einstein-de Sitter limit \[$`\mathrm{\Omega }_m(z)1`$\]. The value of the concentration parameter for minihalos at high redshift is unknown. We extrapolate from simulations of Bullock et al. (2000), which indicate

$$c_\mathrm{N}4.5\left(\frac{M}{10^8M_{\mathrm{}}}\right)^{0.1}\left(\frac{1+z}{7}\right)^1.$$
(6)
For the population of minihalos we adopt for simplicity a single value of $`c_\mathrm{N}`$, the value expected for a minihalo of mass equal to the mass-weighted mean minihalo mass.

Consider a line of sight which passes through a minihalo of mass $`M`$, with an impact parameter of $`\alpha `$ (in units of $`r_{\mathrm{vir}}`$) relative to the center of the minihalo. The total hydrogen column density is then

$$N_{\mathrm{HI}}=5.6\times 10^{19}\frac{c_\mathrm{N}^2}{g(c_\mathrm{N})}f_1(\alpha ,c_\mathrm{N},x_m)\left(\frac{M}{10^8M_{}}\right)^{\frac{1}{3}}\left(\frac{1+z}{7}\right)^2\left(\frac{\mathrm{\Omega }_b}{0.05}\right)\left(\frac{\mathrm{\Omega }_m}{0.3}\right)^{\frac{1}{3}}\left(\frac{h}{0.7}\right)^{\frac{4}{3}}\mathrm{cm}^2,$$
(7)
where we include gas out to a maximum radius $`x_m`$ (in units of $`r_{\mathrm{vir}}`$), and where

$$g(c_\mathrm{N})=\mathrm{log}(1+c_\mathrm{N})\frac{c_\mathrm{N}}{1+c_\mathrm{N}}.$$
(8)
Defining

$$\omega _1=1\alpha ^2c_\mathrm{N}^2$$
(9)
and

$$\omega _2=1+c_\mathrm{N}x_m,$$
(10)
we derive the function

$$f_1(\alpha ,c_\mathrm{N},x_m)=\frac{c_\mathrm{N}}{\omega _2}\frac{\sqrt{x_m^2\alpha ^2}}{\omega _1}+\frac{1}{|\omega _1|^{\frac{3}{2}}}\times \{\begin{array}{cc}\mathrm{log}\left|\frac{\alpha \omega _2}{\alpha ^2c_\mathrm{N}+x_m\sqrt{\omega _1(x_m^2\alpha ^2)}}\right|\hfill & \mathrm{if}\alpha <\frac{1}{c_\mathrm{N}}\hfill \\ \mathrm{sin}^1\left[\frac{\alpha ^2c_\mathrm{N}+x_m}{\alpha \omega _2}\right]\frac{\pi }{2}\hfill & \text{otherwise.}\hfill \end{array}$$
(11)
Similarly, the total recombination rate per unit area, along the same line of sight, is

$$\mathrm{\Gamma }_{\mathrm{rec}}=5.4\times 10^4\frac{c_\mathrm{N}^4}{g^2(c_\mathrm{N})}f_2(\alpha ,c_\mathrm{N},x_m)\left(\frac{M}{10^8M_{}}\right)^{\frac{1}{3}}\left(\frac{1+z}{7}\right)^5\left(\frac{\mathrm{\Omega }_b}{0.05}\right)^2\left(\frac{\mathrm{\Omega }_m}{0.3}\right)^{\frac{1}{3}}\left(\frac{h}{0.7}\right)^{\frac{10}{3}}\mathrm{cm}^2\mathrm{s}^1.$$
(12)
Using $`\omega _1`$ and $`\omega _2`$ from above and defining

$$\omega _3=1\frac{2}{\omega _2}+\frac{\omega _1}{\omega _2^2},$$
(13)
we find

$`f_2(\alpha ,c_\mathrm{N},x_m)`$ $`=`$ $`{\displaystyle \frac{\sqrt{\omega _3}}{\omega _1}}\left(1+{\displaystyle \frac{1}{2\omega _2}}+{\displaystyle \frac{3}{2\omega _1}}\right)+{\displaystyle \frac{\frac{\pi }{2}\mathrm{sin}^1\left(\frac{\alpha }{x_m}\right)}{\alpha c_\mathrm{N}}}+\left(1+{\displaystyle \frac{1}{2\omega _1}}+{\displaystyle \frac{3}{2\omega _1^2}}\right)|\omega _1|^{\frac{1}{2}}`$ (16)
$`\times \{\begin{array}{cc}\mathrm{log}\left|\frac{1\frac{\omega _1}{\omega _2}\sqrt{\omega _1\omega _3}}{\alpha c_\mathrm{N}}\right|\hfill & \mathrm{if}\alpha <\frac{1}{c_\mathrm{N}}\hfill \\ \mathrm{sin}^1\left[\frac{\alpha c_\mathrm{N}+x_m/\alpha }{\omega _2}\right]\frac{\pi }{2}\hfill & \text{otherwise.}\hfill \end{array}`$
### 2.3 Statistics of Screening

The quantities $`\nu _1`$ and $`\nu _2`$ introduced in the previous subsection are sums over the contributions from all minihalos encountered along a line of sight, where a single cross-section of a minihalo at proper distance $`R`$ contributes $`N_{\mathrm{HI}}R^2`$ to $`\nu _1`$ and $`\mathrm{\Gamma }_{\mathrm{rec}}R^2`$ to $`\nu _2`$. We wish to calculate the optical depth $`d\tau `$ for getting a value in the range $`\nu _1`$ to $`\nu _1+d\nu _1`$ by intersecting a minihalo along a line of sight. This involves an integral along the line of sight of the minihalo abundance times cross-sectional area. We first note that the cross-sectional area of a single minihalo, for a given fractional impact parameter $`\alpha `$, is $`\pi \alpha ^2r_{\mathrm{vir}}^2`$. Thus, the contribution to a given value of $`\nu _1`$ is $`2\pi \alpha d\alpha r_{\mathrm{vir}}^2`$; determining the range $`d\alpha `$ that contributes to a given range $`d\nu _1`$, we arrive at the total optical depth:

$$\nu _1\frac{d\tau }{d\nu _1}=_0^{R_{\mathrm{max}}}𝑑R_{M_{\mathrm{min}}}^{M_{\mathrm{max}}}𝑑n(M,z)b(R,M,z)\mathrm{\hspace{0.17em}2}\pi \alpha r_{\mathrm{vir}}^2(M,z)\frac{f_1(\alpha ,c_\mathrm{N},x_m)}{\left|\frac{}{\alpha }f_1(\alpha ,c_\mathrm{N},x_m)\right|},$$
(17)
where $`dn`$ is the halo abundance, $`b(R,M,z)`$ is the bias of minihalos, and everywhere in the integrand $`\alpha `$ is determined by inverting the equation $`N_{\mathrm{HI}}(M,z,\alpha )R^2=\nu _1`$. We get a similar formula for $`\nu _2`$ as well.

A crucial component of the minihalo abundance is the bias $`b(R,M,z)`$, which accounts for the biased distribution of minihalos near the source halo. We discuss here the calculation of bias in some detail. First, we set the bias equal to zero within the virial radius of the source halo, since minihalos cannot exist inside the source halo. Outside the virial radius, the bias can be divided into two components. One can be called “Lagrangian bias”, or “initial condition bias”, and results from the fact that the initial perturbations in nearby points are correlated. Although various analytical approaches to Lagrangian correlations have been suggested, we use the natural and general approach of Scannapieco & Barkana (2002), who worked directly within the standard theoretical framework of excursion sets (Bond et al., 1991; Lacey & Cole, 1993). Their results yield the bias of two halos at a given distance as a function of the two halo masses and the two formation redshifts (although in this application we require only the equal redshift case). Most importantly for the current application, this approach to bias explicitly accounts for the existence of the source halo. For example, there is a strong correlation between nearby rare peaks in the initial density field, which tends to give a positive bias among massive halos at high redshift; on the other hand, all halos are anti-biased at initial positions that are very close to a given halo $`A`$, since much of the mass there is likely to belong to the halo $`A`$ itself when it finally forms. Because of this effect, minihalos are anti-biased at initial positions that are close to a source halo, while source halos themselves are positively correlated at the same distances (though both are anti-biased at the very smallest distances).

The second part of the bias can be called “Eulerian bias” or “infall bias”, and is a result of the high-density environment around the massive source halo when it forms, due to infall of mass from a large initial region in the direction of the forming halo. If, for example, minihalos contain 25% of the mass density in an unbiased region, then in the absence of Lagrangian bias, minihalos would contain 25% of the mass density also near the source halo; thus, the absolute number density of minihalos would increase in proportion to the general mass overdensity near the source halo. In general, the Lagrangian bias must be included as an additional factor beyond this. In the spherical approximation, the Eulerian bias depends only on the final density profile near the virialized source halo; assuming there is no shell crossing (a reasonable approximation outside the virial radius), we can map each mass shell in the final density profile to its initial “Lagrangian” position, and we then use this initial position in order to calculate the Lagrangian bias factor. There are a number of possible approaches to estimate the density profile outside the virial radius of the source halo, but approaches based on spherical collapse, for example, do not match the profile from numerical simulations near the virial radius. However, all approaches give three-dimensional density profiles with power-law slope between $`R^2`$ and $`R^3`$, and in any case, realistically we expect a large scatter in the profiles around halos. For simplicity, here we extrapolate the NFW profile of the source halo to larger radii, an approach which also gives a fall-off between $`R^2`$ and $`R^3`$ in the relative overdensity (i.e., we assume that at large $`R`$ the density approaches its cosmic mean value).

In order to calculate the fraction of ionizing photons from the source that escape out to radius $`R`$, we require the probability distributions of $`\nu _1`$ and $`\nu _2`$. However, the above optical depth $`d\tau `$ does not directly yield the probability; it only accounts for the possibility of intersecting a single minihalo, while there can be multiple intersections along a single line of sight. In addition, the total optical depth can be greater than unity. We determine the actual probability distribution with a Monte Carlo approach. Specifically, we divide the total optical depth to $`100`$ bins, each of which contains a $`\mathrm{\Delta }\tau 1`$. In each Monte Carlo trial, we obtain a total $`\nu _1`$ (for example) by summing over the bins; a given bin makes a contribution with probability $`\mathrm{\Delta }\tau `$, and the contributed value is determined by the distribution $`d\tau /d\nu _1`$ within the bin. We determine the probability distribution using $`100,000`$ trials (for each set of parameters). Note that our calculation assumes that the biasing of the minihalo abundance is determined by the correlation with the nearby source halo, and we neglect any additional correlations among the minihalos themselves.

Once we have obtained the probability distributions of $`\nu _1`$ and $`\nu _2`$, we can calculate the probability distribution of the escape fraction of photons. We often find it convenient to summarize this distribution with a single quantity, the mean escape fraction of photons, averaged over the probability distribution (which corresponds to averaging over solid angle). Including only the constraint of initial ionization (i.e., $`\nu _1`$), the escape fraction is

$$f_{\mathrm{esc}}^{(1)}=_{\nu _1=0}^{dP/d\mathrm{\Omega }}\left(1\frac{\nu _1}{dP/d\mathrm{\Omega }}\right)\frac{dp}{d\nu _1}𝑑\nu _1,$$
(18)
where $`dp/d\nu _1`$ denotes the probability distribution of $`\nu _1`$. We similarly get $`f_{\mathrm{esc}}^{(2)}`$, with $`\nu _1`$ and $`dP/d\mathrm{\Omega }`$ replaced by $`\nu _2`$ and $`d^2P/[d\mathrm{\Omega }dt]`$, respectively. Combining the two constraints exactly would be difficult given the fact that they are correlated, since a given minihalo contributes to both $`\nu _1`$ and $`\nu _2`$; in practice, we estimate the final escape fraction by taking the lower of the two, which is usually $`f_{\mathrm{esc}}^{(2)}`$.

### 2.4 Characteristic Distances

When we consider the escape of photons out to various distances from a source, some values of distance are particularly important. We first want to determine a fiducial distance $`R`$ that represents the beginning of the overlap stage of reionization. The overlap stage starts once neighboring H II regions begin to merge, thus increasing the ionizing intensity within these regions because of the access to multiple ionizing sources. The rising intensity allows the H II regions to expand faster and into denser regions, leading to additional mergers, with the process accelerating until the entire IGM becomes highly ionized except for gas located inside self-shielded, high-density clouds.

One guess for the distance $`R`$ which characterizes the beginning of overlap is simply the typical distance to the nearest other ionizing source. However, this value of $`R`$ does not properly represent the beginning of overlap; first, it may be too large since the photons need not reach the nearby source itself but only the H II region produced by that source; and second, this $`R`$ may be too small, since if the nearby source has a mass much smaller than that of the first source then it will make a negligible addition to the ionizing intensity and this will not advance the process of overlap. We thus choose a distance $`R`$ which more accurately reflects these considerations. First we find that distance where we would typically expect the total mass in additional ionizing sources to be equal to the mass $`M`$ of the first source222We avoid the case where the average mass $`M`$ corresponds to a halo of mass $`M_2M`$, which is present only rarely (with probability $`M/M_2`$); thus, we choose the distance where we expect to either find a single halo with $`M_2>M`$, or a total mass $`M`$ in source halos for which $`M_2<M`$ individually.. Note that we only count active sources (when the duty cycle $`\zeta <1`$). Now, since the other sources also send photons toward the first source, we take only a fraction of the distance (where just the portion beyond the virial radius is reduced). The precise fraction depends on the actual masses of the additional sources; we adopt a fraction of $`1/2`$, which is exact for the symmetric case of having a second source halo of mass equal to $`M`$. This final choice of proper distance $`R`$, which represents the beginning of overlap, is referred to hereafter as half the equal-mass distance.

A second important distance is related to the Gunn-Peterson trough. The essential feature that defines Gunn-Peterson absorption is the presence of long stretches of continuous, strong Ly$`\alpha `$ absorption. Such absorption has been revealed in recent observations at $`z6`$ (Becker et al., 2002; Djorgovski et al., 2002). Neutral hydrogen absorbs photons near the Ly$`\alpha `$ resonance with a very high opacity as given by Gunn & Peterson (1965). As a result, the IGM itself produces strong absorption unless it is highly ionized. The IGM is highly ionized only if the ionizing intensity is sufficiently high. For a given distribution of ionizing sources, each point in the IGM must see sources out to some distance $`R`$ in order to reach the necessary total ionizing intensity. This proper distance $`R`$ is given by \[see equation (18) of Barkana (2002)\]

$$R_{\mathrm{GP}}=3.4\left(\frac{\tau _{\mathrm{Ly}\alpha }}{2.5}\right)^1\left(\frac{1+z}{7}\right)^{\frac{3}{2}}\left(\frac{\mathrm{\Omega }_m}{0.3}\right)^{\frac{1}{2}}\left(\frac{\mathrm{\Omega }_bh}{0.035}\right)\left(\frac{t_S}{5\times 10^8\text{ yr}}\right)\left(\frac{N_{\mathrm{ion}}}{40}\frac{F_{\mathrm{col}}}{0.1}\right)^1\left(\frac{\mathrm{\Delta }_\tau }{0.4}\right)^2\mathrm{Mpc}.$$
(19)
We have defined here the end of the Gunn-Peterson trough as the time when the optical depth drops below $`\tau _{\mathrm{Ly}\alpha }=2.5`$. The total gas fraction inside source halos is $`F_{\mathrm{col}}`$. For the purposes of this estimate we have also used the approximation that all sources have the same lifetime $`t_S`$, where everywhere in this equation we have assumed $`\zeta =1`$ since the value of $`\zeta `$ does not affect the total ionizing intensity. The density in voids which typically dominates the Ly$`\alpha `$ transmission is $`\mathrm{\Delta }_\tau =0.4`$ in units of the cosmic mean gas density. We have neglected source clustering, but such clustering only makes voids relatively empty of sources, increasing the distance $`R`$ required to prevent transmission. For additional discussion of equation (19), see Barkana (2002).

When we investigate the end of the Gunn-Peterson trough, we must also consider the effect of the ionizing background in partially ionizing the minihalos. This is required because the end of the Gunn-Peterson trough occurs relatively late, after an ionizing background may have already been established; furthermore, it is difficult for the central ionizing source to dominate the ionizing intensity at the rather large distance $`R_{\mathrm{GP}}`$. We estimate the effect of the ionizing background as follows. First, we find the distance $`R_{0.5}`$ at which $`f_{\mathrm{esc}}=0.5`$ for a typical source (i.e., a source with mass equal to the mass-weighted mean mass of all sources). We assume that the ionizing background at a random point in the IGM is given by the average integrated intensity due to all sources out to the distance $`R_{0.5}`$, times one half (since $`f_{\mathrm{esc}}=0.5`$ at this distance). With this background intensity, ionization equilibrium yields the neutral fraction at each radius of a given minihalo. For simplicity, we find the neutral fraction at the virial radius of a typical minihalo, and we then assume this neutral fraction to apply throughout every minihalo. This overestimates the effect of the ionizing background, since the central portions of minihalos are in reality denser and more highly neutral; however, as shown in §3, even with this overestimate the overall effect of the partial ionization on screening is small. Once partial ionization is included, $`f_{\mathrm{esc}}`$ increases at a given $`R`$, and thus $`R_{0.5}`$ also increases; we iterate in order to find the self-consistent value of $`R_{0.5}`$, and the corresponding value of the neutral fraction in minihalos. We then use this neutral fraction when we compute screening at $`R=R_{\mathrm{GP}}`$.

### 2.5 Direct Gunn-Peterson Absorption

In addition to their screening of ionizing photons, minihalos also absorb Ly$`\alpha `$ photons and may therefore contribute to the Gunn-Peterson trough. Direct absorption near the Ly$`\alpha `$ resonance is particularly strong (Gunn & Peterson 1965; see also the related discussion in the previous subsection). However, minihalos contain $`20\%`$ of the gas in the universe at an average overdensity of $`200`$; since for any population of absorbers the covered fraction of a line of sight equals the volume fraction of the population333More precisely, the covering factor of a line of sight equals the volume filling factor. Consider, for simplicity, a number density $`n`$ of halos which all have the same shape and inclination to the line of sight. A single halo has a cross-sectional area $`A=_A𝑑x𝑑y`$, where we use two-dimensional $`(x,y)`$ coordinates. Denote the total depth of the halo at each point $`D(x,y)`$. On a random line of sight, the halos cover a fraction $`=n_A𝑑x𝑑yD(x,y)`$. This integral equals the volume of a single halo. This proof is easily generalized to a population of halos with various shapes., minihalos directly cover only $`0.1\%`$ of a random line of sight to a quasar. More interesting, therefore, is the possibility of absorption due to the damping wings of the Ly$`\alpha `$ line (Miralda-Escudé, 1998).

To estimate the probability that the population of minihalos produces strong absorption, we first consider the optical depth due to a single minihalo; we do this as a function of distance $`R`$ from the minihalo, where the distance represents different wavelengths $`\lambda `$ by the correspondence of $`\lambda `$ and $`z`$ via a pure Hubble velocity flow \[see the similar approach, e.g., in Loeb & Rybicki (1999)\]. We first note a number of simplifications: the size of each halo (typically a few kpc) is much smaller than the region of significant absorption due to the halo’s damping wings (typically a fraction of a Mpc); the latter region, in turn, is much smaller than the horizon at the reionization redshift. Under these conditions, the damping wing produces an optical depth $`\tau _{\mathrm{Ly}\alpha }`$ which is simply proportional to the total column density of H I divided by the distance squared from the minihalo. This results from the Lorentzian damping profile of the line, and it corresponds to the well-known result that when the damping wings dominate, the equivalent width of an absorption line is proportional to the square root of the column density (note that at small distances from the absorber, distance is proportional to change of wavelength $`\mathrm{\Delta }\lambda `$). Specifically, at a distance $`R`$ from a column density $`N_{\mathrm{HI}}`$,

$$\tau _{\mathrm{Ly}\alpha }=0.51\left(\frac{N_{\mathrm{HI}}}{10^{20}\mathrm{cm}^2}\right)\left(\frac{R}{1\mathrm{Mpc}}\right)^2\left(\frac{1+z}{7}\right)^3\left(\frac{\mathrm{\Omega }_m}{0.3}\right)^1\left(\frac{h}{0.7}\right)^2.$$
(20)
Including the damping wings of the entire population of minihalos, we calculate the total optical depth (or covering factor) $`C`$ along the line of sight, for getting a Ly$`\alpha `$ optical depth above some value $`\tau _{\mathrm{Ly}\alpha }`$. The cross-sectional area contributed by the impact parameter range $`d\alpha `$ is $`2\pi \alpha d\alpha r_{\mathrm{vir}}^2`$; for a given $`\alpha `$, if we determine the maximum distance $`R`$ from the minihalo for which the optical depth is greater than $`\tau _{\mathrm{Ly}\alpha }`$, then the minihalo covers with this optical depth a total line-of-sight distance $`2R`$ (including the blue and red damping wings). Thus, the covering factor equals

$$C(>\tau _{\mathrm{Ly}\alpha })=_0^{x_{\mathrm{max}}}𝑑\alpha _{M_{\mathrm{min}}}^{M_{\mathrm{max}}}𝑑n(M,z)4\pi \alpha r_{\mathrm{vir}}^2(M,z)R,$$
(21)
where $`R`$ in the integrand is determined from $`\tau _{\mathrm{Ly}\alpha }`$ by inverting equation (20) for each $`M`$ and $`\alpha `$, using equation (7) for $`N_{\mathrm{HI}}`$. As noted in §2.2, we also adopt a single concentration parameter $`c_\mathrm{N}`$ for all minihalos; this gives them similar density profiles, which simplifies the calculation (i.e., the double integral in equation (21) can then be factored into two single-variable integrals). Note that the covering factor is inversely proportional to $`\sqrt{\tau _{\mathrm{Ly}\alpha }}`$.

As in §2.3, here too the optical depth includes only the contribution from one minihalo at a time, and we determine the full probability distribution of values of $`\tau _{\mathrm{Ly}\alpha }`$ with Monte Carlo trials. Note that in this calculation, we do not include clustering among the minihalos themselves. Clustering will change the probability distribution somewhat, increasing the optical depth in regions of clustered halos while decreasing the optical depth in relatively empty regions. However, the damping wings effectively smooth the absorption on the scale of a few comoving Mpc, reducing the effect of clustering. There is a similar smoothing of the effect of peculiar velocities, which we likewise neglect.

### 2.6 Photoevaporation of Minihalos

As discussed above, the minihalos can affect the evolution of the IGM in a number of ways; the gas, however, does not remain in minihalos forever. Barkana & Loeb (1999) showed that photoionization heating by the cosmic UV background can evaporate much of this gas back into the IGM during the overlap era of reionization. This work also showed that this process affects a broad range of halo masses, with only a small gas fraction evaporating out of $`10^8M_{\mathrm{}}`$ halos, but with halos below $`10^6M_{\mathrm{}}`$ losing their entire gas content because of their shallow gravitational potential wells. Even after the initial loss of unbound gas, the bound gas in the minihalo expands and may further evaporate; however, this stage proceeds more slowly.

In this paper we consider two main stages in this process. In the first stage, before reionization has begun, the minihalos are fully intact. In the second stage, only the bound gas remains, a situation which we approximate by keeping only the minihalos in which $`>50\%`$ of the gas is initially bound after UV heating; this gives a minimum circular velocity $`V_c11\mathrm{km}\mathrm{s}^1`$ which is independent of redshift (Barkana & Loeb, 1999). In the transition from the first to the second stage, the typical mass of the minihalos that photoevaporate (i.e., their mass-weighted mean) is $`7\times 10^6M_{\mathrm{}}`$. Such a halo has a virial radius of 0.9 kpc, and once it is heated it should evaporate with a $`10\mathrm{km}\mathrm{s}^1`$ outflow in $`90`$ million years, which is equivalent to the time between $`z=6.5`$ and $`z=6`$. Additional evolution during the second stage is even slower; the expansion of the remaining minihalos (typical $`M9\times 10^7M_{\mathrm{}}`$) takes place on a timescale of $`200`$ million years (equivalent to the redshift interval between $`z=7.2`$ and $`z=6`$), and during this time each minihalo is expected to accrete additional dark matter or merge with other minihalos. In any case, it is clear that if minihalos obstruct reionization then their relatively slow photoevaporation can stretch the reionization era over a significant redshift range, as we discuss further below.

The two-stage process that we have defined is only an approximation and it neglects some complications. For example, when an ionizing source begins to ionize the minihalos around it, those minihalos can immediately begin to photo-evaporate; however, as just discussed, most of the gas in minihalos takes a substantial amount of time to photoevaporate, even in the first stage which involves the unbound gas. Our calculations apply during this period.

## 3 Results

Our main goal is to determine quantitatively the degree of effective screening of ionizing sources due to the minihalo population. Figure 1 shows the escape fraction $`f_{\mathrm{esc}}`$ of ionizing photons as a function of the physical distance $`R`$ from the source. We pick a representative of the population of ionizing sources, namely a halo of mass $`M=2.2\times 10^{10}M_{}`$, which is the mass-weighted mean source halo mass (assuming that the minimum halo mass for sources is $`1.6\times 10^8M_{}`$ as determined by atomic cooling). As shown in the figure, the typical distance from this halo to the nearest other ionizing source is $`R=24`$ kpc, giving a high $`f_{\mathrm{esc}}=92\%`$. However, as discussed in §2.4, the distance which best represents the beginning of overlap is half the equal-mass distance. For the case shown in the figure, this distance is $`R=73`$ kpc, giving a low $`f_{\mathrm{esc}}=18\%`$.

Figure 1 also illustrates the importance of a full calculation of the biased distribution of minihalos and of other source halos around the first source halo. If we did not include bias, we would calculate half the equal-mass distance to be $`R=0.13`$ Mpc, and $`f_{\mathrm{esc}}=15\%`$. Including only the infall bias (see §2.3) increases the abundance of both minihalos and source halos444When the second halo mass $`M_2>M`$, infall bias is computed using infall onto $`M_2`$. Also, regardless of the choice of bias, $`M`$ is not allowed to fall within the virial radius of $`M_2`$.. However, $`R`$ is reduced only slightly to 0.12 Mpc, since at this large distance from the virial radius, the total enclosed mass is essentially unchanged by infall; however, $`f_{\mathrm{esc}}=0.03\%`$ is still strongly reduced at these large distances because of the screening due to the large concentration of minihalos very close to the source. To better understand this, consider moving a given minihalo twice closer to the source; the resulting halo then contributes 4 times the covering factor, but produces values of $`\nu _1`$ and $`\nu _2`$ that are reduced by $`1/4`$. We indeed find that infall bias increases the covering factor most significantly at the lowest values of $`\nu _1`$ and $`\nu _2`$. Finally, when Lagrangian bias is included, $`R`$ is further reduced due to the strong positive correlation among source halos; $`f_{\mathrm{esc}}`$, however, actually increases (at a given $`R`$) because of a slight anti-bias of minihalos around the source halo (see §2.3).

Figure 2 shows the full probability distribution of $`\nu _2`$, from which $`f_{\mathrm{esc}}`$ (shown in Figure 1) is derived using equation (18). The figure compares the results of the Monte Carlo calculation, which includes the effect of screening due to multiple minihalos, to the naive probability $`P(>\nu _2)=1e^\tau `$, given an optical depth $`\tau `$ for having a recombination rate per unit solid angle greater than $`\nu _2`$. Clearly, the effect of multiple screening changes completely the probability distribution at high values of $`\nu _2`$; however, whenever the probability of screening the source is significant, it is dominated by the contribution of relatively low values of $`\nu _2`$, and the mean $`f_{\mathrm{esc}}`$ is only moderately affected by multiple scattering. For example, for the three cases shown in the figure, the naive probability gives $`f_{\mathrm{esc}}=5.8\%`$, $`21\%`$, and $`93\%`$ (in order of increasing $`R`$), while the exact probability distribution yields $`f_{\mathrm{esc}}=4.2\%`$, $`18\%`$, and $`92\%`$, respectively. Note also that our inclusion of partial blocking as in equation (18) is crucial for getting the correct value of $`f_{\mathrm{esc}}`$. For example, using the total covering factor of all values of $`\nu _2`$, the naive probability of having no blocking at all is $`2.0\%`$, $`7.7\%`$, and $`64\%`$, respectively; this probability is not an accurate estimate for $`f_{\mathrm{esc}}`$.

Although Figure 1 illustrates the results for a single representative source halo mass, a wide range of halo masses may in fact contribute to reionization. Figure 3 illustrates how the effectiveness of screening varies with halo mass. The variation with halo mass is complicated by several competing effects. On the one hand, larger halos produce a higher flux of ionizing photons. On the other hand, larger halos must in general send their photons farther before they can effectively participate in overlap, i.e., there is an increase with halo mass in the equal-mass distance; however, the rate of increase declines at the high-mass end because of the powerful clustering of extremely massive halos. The end result is that for the case of $`z=6`$, minihalos screen all source halos up to $`2\times 10^{12}M_{}`$ roughly equally, with $`13\%<f_{\mathrm{esc}}<27\%`$; although halos with even higher masses are not strongly screened, most ionizing photons are produced by halos for which $`f_{\mathrm{esc}}`$ is low.

Figure 3 also shows the separate escape fractions which result from the obstacles (see §2.2) caused by initial ionizations ($`\nu _1`$) and by recombinations ($`\nu _2`$). In each case, recombinations give the stronger constraint, hence we have generally used recombinations to determine escape fractions. The figure shows as well the dependence of $`f_{\mathrm{esc}}`$ on redshift, where we make the individual sources more efficient at high redshift, so that the total production rate of ionizing photons remains fixed, and thus reionization can occur at the higher redshift. At high redshift, the total mass fraction in minihalos is smaller, but the biasing around the source is stronger both for minihalos and for other sources. The final result is that for a typical source (whose mass declines with redshift), $`f_{\mathrm{esc}}`$ increases with redshift, from $`18\%`$ at $`z=6`$ to $`28\%`$ at $`z=10`$ and $`44\%`$ at $`z=15`$.

We have thus far adopted the minimum source mass of $`1.6\times 10^8M_{}`$ given by atomic cooling, and have assumed that all halos above this mass produce ionizing photons with the same overall efficiency. However, internal stellar feedback effects may be especially disruptive in source halos with relatively low masses and low binding energies. Therefore, we consider the possibility of requiring a higher minimum mass for halos that host efficient ionizing sources. At the same time we make the individual sources more efficient, so that the reionization redshift remains approximately fixed. Figure 4 shows that the resulting $`f_{\mathrm{esc}}`$, for the characteristic source halo mass, increases at a given $`R`$ as the minimum source halo mass is increased. The ionizing photons, however, must reach a more distant $`R`$ before overlap begins, since the abundance of ionizing sources is reduced. The final result is that the screening depends only mildly on the minimum source halo mass. For the four cases between $`1.6\times 10^8M_{}`$ and $`1.6\times 10^{11}M_{}`$ shown in Figure 4, $`f_{\mathrm{esc}}=18\%`$, $`14\%`$, $`11\%`$, and $`9.2\%`$, respectively, i.e., $`f_{\mathrm{esc}}`$ is reduced by only a factor of two when the minimum halo mass is increased by three orders of magnitude.

Clearly, screening due to fully intact minihalos is a substantial barrier and can delay the beginning of overlap. Early on, when the dense gas surrounding each minihalo is still neutral, the screening is even more effective. For example, when the minimum source mass is given by atomic cooling, $`f_{\mathrm{esc}}=18\%`$ if we include the gas out to the virial radius of each minihalo; if we extrapolate the NFW profile out to twice the virial radius, around each minihalo, the result is $`f_{\mathrm{esc}}=4.1\%`$. Thus, the surroundings of the minihalos will be ionized before overlap can make substantial progress, and at that point the minihalos themselves will begin to photoevaporate. Figure 4 also shows the results for the later stage discussed in §2.6, when only the bound gas remains in the minihalos. Once the unbound gas has photoevaporated, the effect of screening no longer substantially impedes the start of overlap. In particular, the four cases shown in Figure 4 give $`f_{\mathrm{esc}}=84\%`$, $`82\%`$, $`78\%`$, and $`76\%`$, respectively, in order of increasing minimum halo mass.

The results also depend on additional parameters. For example, our standard case of $`\zeta =0.25`$ and $`\eta =10\%`$ gives $`f_{\mathrm{esc}}=18\%`$ for intact minihalos and $`84\%`$ for bound minihalos only. Increasing $`\zeta `$ to 0.5 gives $`38\%`$ and $`90\%`$ for these two cases, respectively. Note that increasing $`\zeta `$ at a fixed efficiency makes each source dimmer but increases the number density of sources, thus reducing the distance that photons are required to reach in order for the H II bubbles of nearby sources to overlap. Reducing $`\zeta `$ to 0.25 gives $`7.1\%`$ and $`77\%`$. The results are less sensitive to changing $`\eta `$; setting $`\eta =20\%`$ gives $`24\%`$ and $`86\%`$, while setting $`\eta =5\%`$ gives $`13\%`$ and $`83\%`$, respectively. Thus, the detailed numbers depend on the uncertain source parameters, but the qualitative conclusion remains the same: intact minihalos substantially impede the overlap stage, while bound minihalos do not.

Although overlap can proceed once only bound minihalos remain, observationally the universe remains opaque to Ly$`\alpha `$ radiation. In other words, the ionizing intensity in the IGM is initially low, and the small remaining H I fraction continues to produce a Gunn-Peterson trough. The trough disappears only when sources can send their photons out to the rather large distance $`R_{\mathrm{GP}}`$ (see §2.4). Even the bound minihalos on their own produce (for our standard parameters, and a typical source) a covering factor of 71 for total blocking, which implies absolutely total screening for $`R=R_{\mathrm{GP}}`$. Note that we include in this calculation the effect of the ionizing background in partially ionizing the minihalos (see §2.4). The effect of this on the overall screening is small since in the current context, when a line of sight from a source encounters a minihalo the typical result is complete blocking even though the neutral fraction at the virial radius of the minihalo is $`1\%`$; this results from the fact that we are considering only bound minihalos (which are relatively massive, and thus have high column densities), and we are considering large distances from the source (so that the column density or recombination rate per unit solid angle is high).

Thus, the screening effect can substantially delay the end of the Gunn-Peterson trough, and this remains true even if the gas in some minihalos cools due to H2 molecules and only a small fraction of the minihalos are left intact. The history of molecular cooling in minihalos prior to reionization involves complicated effects of radiative feedback, effects which are now beginning to be explored with a new generation of numerical simulation codes. Although only redshifts $`10`$ can be computed and the results must be considered tentative, Ricotti, Gnedin, & Shull (2002) found that molecular cooling does occur in some minihalos, but the total ionizing intensity is limited by a strong negative feedback of direct ionization of H2 molecules; in particular, they argued that minihalos cannot contribute substantially to overlap compared to more massive ionizing sources. Future simulations may determine the precise fraction of minihalos in which the gas remains uncooled. We find that even if only $`5\%`$ of the minihalos are intact, the bound ones produce essentially complete blocking (i.e., $`f_{\mathrm{esc}}=2.5\%`$) at $`R=R_{\mathrm{GP}}`$; with $`2\%`$ of the minihalos the escape fraction is still small ($`23\%`$), while $`1\%`$ of the bound minihalos block around half the photons (i.e., $`f_{\mathrm{esc}}=48\%`$).

Finally, we consider the absorption of Ly$`\alpha `$ photons by minihalos (§2.5). Since the covering factor (or optical depth) $`C`$ is a power law in $`\tau _{\mathrm{Ly}\alpha }`$, the conversion via Monte Carlo to a cumulative probability $`P(>\tau _{\mathrm{Ly}\alpha })`$ involves a universal function, where $`P`$ is a function of $`C`$ only and does not depend separately on $`\tau _{\mathrm{Ly}\alpha }`$. This universal function is plotted in Figure 5, where this function is shown to be significantly different from the naive probability of $`1e^C`$. The probability is $`50\%`$ when $`C=0.53`$, $`68\%`$ when $`C=0.78`$, $`95\%`$ when $`C=1.5`$, and $`99\%`$ when $`C=2.0`$.

For the covering factor itself we obtain with fully intact minihalos, $`C(>\tau _{\mathrm{Ly}\alpha })=0.97/\sqrt{\tau _{\mathrm{Ly}\alpha }}`$. Therefore, $`P(>1)=78\%`$ of the line of sight is covered with $`\tau _{\mathrm{Ly}\alpha }>1`$ by minihalo damping wings; the number is $`P(>5)=42\%`$ with $`\tau _{\mathrm{Ly}\alpha }>5`$. In the earliest stage, when we can extrapolate the NFW profile out to twice the virial radius around each minihalo, the corresponding results are $`C(>\tau _{\mathrm{Ly}\alpha })=2.3/\sqrt{\tau _{\mathrm{Ly}\alpha }}`$, with $`P(>1)=99.6\%`$ and $`P(>5)=81\%`$. On the other hand, in the late stage when only the bound gas remains inside minihalos, $`C(>\tau _{\mathrm{Ly}\alpha })=0.083/\sqrt{\tau _{\mathrm{Ly}\alpha }}`$, with $`P(>1)=8.3\%`$ and $`P(>5)=3.7\%`$. The low optical depth produced by the bound gas results from the fact that, for a given total mass, small halos are more effective Ly$`\alpha `$ absorbers than large halos. This results from equation (20): a reduction in column density by a factor of 2 reduces the distance covered (by a given $`\tau _{\mathrm{Ly}\alpha }`$) by only a factor of $`\sqrt{2}`$. We conclude that as long as the minihalos are intact, they can produce significant Ly$`\alpha `$ absorption, although they are unlikely to themselves cause full absorption (i.e., $`\tau _{\mathrm{Ly}\alpha }>5`$ with a probability near unity). After the unbound gas photoevaporates, however, the minihalos can make only a minor contribution to Ly$`\alpha `$ absorption.

## 4 Conclusions

In this paper we have shown that the population of gaseous minihalos can produce substantial screening of ionizing photons during reionization. After accounting for the biased distribution of minihalos around sources, we have found that a fully intact population can delay the start of overlap; for a typical source at redshift $`z=6`$, advancing overlap requires getting the ionizing photons out to a proper distance $`R=73`$ kpc, for which the minihalos limit the escape fraction to a low value of $`f_{\mathrm{esc}}=18\%`$. The escape fraction out to the appropriate distance is roughly the same for source halos over a wide range of mass, with only rare $`M>10^{13}M_{}`$ halos having a much larger escape fraction. Once the smallest minihalos photo-evaporate, gas remains in the bound minihalos for a substantially longer period. The bound minihalos are unable to delay the start of overlap, producing a value of $`f_{\mathrm{esc}}=84\%`$ (for our standard parameters at $`z=6`$).

To obtain the escape fraction, we have included the effect of screening due to multiple minihalos (using a Monte Carlo approach), and also the effect of partial blocking \[see equation (18)\]. We typically adopted the minimum source mass of $`1.6\times 10^8M_{}`$ given by atomic cooling; however, if stellar feedback favors higher-mass halos, we have shown that $`f_{\mathrm{esc}}`$ is insensitive to the minimum source halo mass if the reionization redshift is held fixed; and although the detailed numbers depend on uncertain source parameters (such as efficiency and duty cycle), we have shown that the overall conclusion remains the same: the full population of minihalos can substantially impede the overlap stage, while the bound minihalos themselves cannot.

Even after overlap, screening due to minihalos remains highly effective until the minihalos photo-evaporate; this screening keeps the ionizing intensity in the IGM low, and the Gunn-Peterson trough remains until even the bound minihalos are disrupted. This can produce a lag of $`200`$ Myr (equivalent to redshift interval $`z=7.2`$– 6) between overlap and the end of Gunn-Peterson absorption. This screening is so strong that it is highly effective even if molecular cooling occurs in most minihalos and only a few percent of them retain uncooled gas.

Intact minihalos may also directly absorb Ly$`\alpha `$ photons. For the full population at $`z=6`$, $`42\%`$ of a random line of sight is covered with $`\tau _{\mathrm{Ly}\alpha }>5`$. However, as noted above, a full population of minihalos screens ionizing sources and makes the IGM absorb Ly$`\alpha `$ photons much stronger than do the minihalos themselves. When only the bound minihalos remain, the ionizing intensity is still low in the IGM, but at this stage the direct absorption due to minihalos is already insignificant (e.g., only $`3.7\%`$ of the line of sight is covered at $`\tau _{\mathrm{Ly}\alpha }>5`$).

Current numerical simulations cannot resolve the population of minihalos during the reionization era. However, we have demonstrated in this paper that any conclusions regarding the timing of the overlap era of reionization and of the lifting of Gunn-Peterson absorption depend crucially on the behavior of the minihalos.

This work was supported in part by NASA grants NAG 5-7039, 5-7768, and NSF grants AST-9900877, AST-0071019 for AL. RB acknowledges the support of an Alon Fellowship at Tel Aviv University.


format: text

text:
# Isobaric multiplet yrast energies and isospin non-conserving forces

## Abstract

The isovector and isotensor energy differences between yrast states of isobaric multiplets in the lower half of the $`pf`$ region are quantitatively reproduced in a shell model context. The isospin non-conserving nuclear interactions are found to be at least as important as the Coulomb potential. Their isovector and isotensor channels are dominated by $`J=2`$ and $`J=0`$ pairing terms, respectively. The results are sensitive to the radii of the states, whose evolution along the yrast band can be accurately followed.

The electrostatic energy of a sphere of radius $`R`$ and charge $`Ze`$ is easily calculated to be $`E_C=3e^2Z^2/5R`$. It is under this guise that the Coulomb field enters the Bethe Weizäcker mass formula, and becomes a basic quantity in nuclear structure. Direct evidence of entirely Coulomb effects has long been available from displacement energies between mirror ground states (MDE), and in the last decade from differences in excitation energies of yrast bands in mirror nuclei (MED) cameron ; 49 ; 47-49 ; 51 ; 50fe 111As isospin non conserving nuclear forces play an important role we propose to replace the traditional acronym for Coulomb displacement energies—CDE—by MDE, where M stands for mirror. Similarly CED becomes MED..

The MDE energies range from few to tens of MeV. They are given mainly by $`E_C`$, but precise calculations were found to be unexpectedly inaccurate—the Nolen-Schiffer anomaly NS69 —and revealed the necessity to introduce charge symmetry breaking (CSB) nuclear potentials (mac01  and references therein). The anomaly is now under control to within shell effects bro00b ; duf02 , which we define as deviations from a Bethe Weiszäcker-type formula involving only number of particles ($`A`$), isospin ($`T`$), and its third component ($`T_z`$).

The MED are defined by ($`Z_>`$ and $`Z_<`$ are the largest and smallest $`Z`$ in the multiplet, and $`E_J`$ are the yrast excitation energies)

$$\mathrm{MED}_J=E_J(\overline{Z}+T)E_J(\overline{Z}T),\overline{Z}=\frac{Z_>+Z_<}{2}.$$
(1)
The observed MED are very small (of the order of 10-100 keV), and entirely due to shell effects. Recently, the experimental information on yrast bands has been extended to isospin triplets 50mn ; 46cr , thus determining new quantities, the TED given by

$$\mathrm{TED}_J=E_J(\overline{Z}+1)+E_J(\overline{Z}1)2E_J(\overline{Z}).$$
(2)
Both measurements are needed to achieve a clear understanding of the interplay between the Coulomb potential $`V_C`$, and $`V_B`$, the isospin breaking nuclear interaction.

To analyze them, we start by writing the isovector, $`\beta _𝒓^{(1)}`$, and isotensor, $`\beta _𝒓^{(2)}`$, contributions to $`V_B`$, as linear combinations of two body matrix elements in neutron-proton ($`\nu \pi `$) formalism ($`𝐫`$ $`r_1r_2r_3r_4`$, where $`r_i`$ is a subshell):

$$\beta _𝒓^{(1)}=V_{B𝒓}^{\pi \pi }V_{B𝒓}^{\nu \nu },\beta _𝒓^{(2)}=V_{B𝒓}^{\pi \pi }+V_{B𝒓}^{\nu \nu }2V_{B𝒓}^{\pi \nu }.$$
(3)
The isoscalar contribution $`\beta _𝒓^{(0)}=V_{B𝒓}^{\pi \pi }+V_{B𝒓}^{\nu \nu }+V_{B𝒓}^{\pi \nu }`$ is nil in $`V_B`$ , while for $`V_C`$ we have $`\beta _{C𝒓}^{(0)}=\beta _{C𝒓}^{(1)}=\beta _{C𝒓}^{(2)}=V_{C𝒓}^{\pi \pi }`$.

The MED are entirely of isovector origin and the first exact shell model calculations in the full $`pf`$ shell indicated that $`V_C^{\mathrm{ho}}`$, i.e., calculated in the harmonic oscillator (ho) basis fails to give a satisfactory description 47-49 . The way out proposed in this reference consisted in replacing the harmonic oscillator matrix elements $`V_{C𝐟_{\mathrm{𝟕}/\mathrm{𝟐}}}^{\mathrm{ho}}`$ by empirical ones derived from the $`A=42`$ spectrum which are very different. Therefore, it was hard to attribute the replacement to a renormalization of $`V_C`$, expected to be small. But it was equally hard to think in terms of CSB precisely because the effect was so large. Nevertheless the ansatz (or variants of it) worked quite well, and subsequent calculations incorporated it 51 ; 50fe  leading eventually to (almost) full quantitative agreement zuk01  for the MED in $`A=47`$, 49, 50 and 51. When the isotensor TED data came in, it became clear that both charge independence breaking mac01b  and CSB had to be invoked.

Ironically, this result is obvious from the, long known, $`A=42`$ spectra end90 . Assuming that the observed states are essentially $`f_{7/2}^2`$ configurations on top of the 40Ca core, these spectra define an interaction in the $`f_{7/2}`$ subshell. Therefore, by setting $`V_CV_{C𝐟_{\mathrm{𝟕}/\mathrm{𝟐}}}^{\mathrm{ho}}`$, the nuclear isovector and isotensor contributions can be extracted. They are shown in Table 1, where their centroids —$`_J(2J+1)V_{𝐟_{\mathrm{𝟕}/\mathrm{𝟐}}}^J/_J(2J+1)`$—have been subtracted for clarity. The assumption of $`f_{7/2}^2`$ dominance is not warranted, as—at least— the $`J=0`$ and 2 states are known to mix with core excitations. Therefore a safer procedure consists in replacing the lowest observed states by the $`f_{7/2}^2`$ centroids estimated from spectroscopic factors end90 . However, when this is done, no significant change obtains in Table 1, whose indications must therefore be taken very seriously.

A renormalized $`V_C`$ adapted to the $`f_{7/2}^2`$ space will remain a purely $`\pi \pi `$ force, and therefore the same for the isovector and isotensor channels. Furthermore, it is expected to be reasonably close to the bare $`V_C`$ in the first line of Table 1. Therefore, upon subtracting this bare $`V_C`$ from the observed data—second and third lines—we expect the same, reasonably small numbers. It is obvious that the corresponding numbers are neither equal nor small. The unavoidable conclusion is that the $`A=42`$ data indicate that the role of isospin non conserving nuclear forces is at least as important as that of the Coulomb potential in the observed MED and TED.

For the full description of these quantities in $`A=46`$-51 we rely on exact, isospin conserving shell model calculations ANTOINE  with single particle spectrum from 41Ca and the KB3G interaction. Very little changes are observed if the other standard interactions are used (KB3, FPD6, all defined in pov01 ). The energy differences are obtained in first order perturbation theory 222Diagonalizing instead of taking expectation values makes very little difference., as the sum of expectation values, in which we separate the monopole and multipole components of the Coulomb field $`V_C=V_{Cm}+V_{CM}`$ following Refs. duf02 ; zuk01 :

$`\mathrm{MED}_J`$ $`=`$ $`\mathrm{\Delta }_\mathrm{M}V_{Cm}_J+\mathrm{\Delta }_\mathrm{M}V_{CM}_J+\mathrm{\Delta }_\mathrm{M}V_B_J,`$ (4)
$`\mathrm{TED}_J`$ $`=`$ $`\mathrm{\Delta }_\mathrm{T}V_{CM}_J+\mathrm{\Delta }_\mathrm{T}V_B_J.`$ (5)
The monopole $`V_{Cm}`$ contains all terms quadratic in scalar products of Fermion operators $`a_i^+a_j`$. The non diagonal contributions ($`ij`$) lead to isospin mixing that demands second order perturbation theory. They will be considered here only through their influence on the radial wavefunctions, i.e., the Thomas Ehrman shift that depresses the single particle $`p_{3/2}`$ state in 41Sc by 200 keV below its analogue in 41Ca. The diagonal part ($`i=j`$) involves only proton number operators. It contains $`E_C`$ plus a single particle splitting induced by $`V_C`$ on the orbits of principal quantum number $`p`$ above harmonic oscillator (ho) closed shell $`Z_{cs}`$ 333This term can be derived from duf02 \[Eqs. 5,10-12\]. Hint: the overall coefficient in Eq. (6) comes from $`1.934\times 0.383\times (0.01)\times 1.5220.001`$.:

$$\epsilon _{Cl}=\frac{0.001Z_{cs}^{13/12}[2l(l+1)p(p+3/2)]}{A^{1/3}(p+3/2)}\mathrm{MeV}.$$
(6)
The effect of $`E_C`$ is proportional to the difference of (inverse) radii between a $`J`$-yrast and the ground state 50fe . The total radii depend on those of the individual orbits, and therefore—to good approximation—on the average neutron plus proton occupancies for each orbit, which we denote by $`m_k_J/2`$, with $`m_k=z_k+n_k`$ (number of neutrons plus number of protons in orbit $`k`$). We take averages relying on the near equality of proton radii in both members of a mirror pair duf02 . As it is reasonable to assume that orbital radii depend only on $`l`$, and the $`p_{1/2}`$ occupancy is always negligible, the whole radial effect will be taken to depend on the $`p_{3/2}`$ occupancy. Note that the single particle contribution from Eq. (6) is proportional to the difference of proton and neutron occupancies. It is important in $`A=41`$, but typically ten times smaller than the radial effect in $`A=`$47-51, so we neglect it and end up with $`\mathrm{\Delta }_MV_{Cm}_J=a_mm_{p_{3/2}}_J/2`$. The value of $`a_m`$ can be estimated by adding to the observed shift the single particle splitting (6) that depresses the $`l=3`$ orbits with respect to the $`l=1`$ ones (by 125 keV at $`A=40`$, $`Z_{cs}=20`$). Then, $`a_m.200+.125=0.325`$ MeV.

In the isotensor case the $`m_{p_{3/2}}`$ contributions cancel out.

The multipole contribution $`\mathrm{\Delta }V_{CM}_J`$ is calculated using oscillator Coulomb matrix elements in the $`pf`$ shell.

The only direct information on $`V_B`$ comes from Table 1. To make use of it we must explore the possibility of specifying an interaction acting in the full $`pf`$ shell, solely in terms of $`f_{7/2}`$ matrix elements. The idea turns out to be quite viable using the multiplicative prescription444This form proves sufficient for our present needs, but it is case-dependent and one can find examples that demand a term in $`am_{p_{3/2}}_J`$.

$$\mathrm{\Delta }V_{C\mathrm{𝐩𝐟}}^{ho}_J=b\mathrm{\Delta }V_{C𝐟_{\mathrm{𝟕}/\mathrm{𝟐}}}^{ho}_J,$$
(7)
as illustrated for the Coulomb potential in Fig. 1.

The same form efficiently relates the schematic pairing or quadrupole pairing forces in the $`pf`$ shell to the $`V_{𝐟_{\mathrm{𝟕}/\mathrm{𝟐}}}^{J=0\mathrm{or}J=2}`$ matrix elements. To minimize the number of parameters, for $`V_B`$ we retain only the leading term suggested by Table 1, and set $`\beta _{\mathrm{𝐩𝐟}}^{(1)}=\beta _1V_{𝐟_{\mathrm{𝟕}/\mathrm{𝟐}}}^{J=2}`$, $`\beta _{\mathrm{𝐩𝐟}}^{(2)}=\beta _2V_{𝐟_{\mathrm{𝟕}/\mathrm{𝟐}}}^{J=0}`$, where $`V_{𝐟_{\mathrm{𝟕}/\mathrm{𝟐}}}^J`$ is the matrix element with unit value. Collecting all the pieces we have

$`\mathrm{MED}_J`$ $`=`$ $`{\displaystyle \frac{1}{2}}a_mm_{p_{3/2}}_J+\mathrm{\Delta }_\mathrm{M}V_{C\mathrm{𝐩𝐟}}^{ho}+\beta _1V_{𝐟_{\mathrm{𝟕}/\mathrm{𝟐}}}^{J=2}_J,`$ (8)
$`\mathrm{TED}_J`$ $`=`$ $`\mathrm{\Delta }_\mathrm{T}V_{C\mathrm{𝐩𝐟}}^{ho}+\beta _2V_{𝐟_{\mathrm{𝟕}/\mathrm{𝟐}}}^{J=0}_J.`$ (9)
In Figs. 2 for the MED, $`V_{Cm},V_{CM}`$, and $`V_{BM}`$ stand respectively for the first, second and third terms in Eq. (8). The parameters are taken to be round numbers, $`a_m=300`$ keV and $`\beta _1=\beta _2=100`$ keV.

The reduction of $`V_B`$, for MED and TED, to a single matrix element is an oversimplification, but the results are so satisfactory that the need of extra terms is not felt. The only parameter-free alternative is to take matrix elements with the weights in Table 1. However, this choice is arbitrary because—from the discussion around Eq. (7) and Fig. 1—we expect a case by case (even matrix element by matrix element) renormalization. Nonetheless, though the agreement with experiment becomes less impressive, it remains acceptable. The conclusion is that the leading term in $`J=2`$ for MED is indeed dominant, and that in $`J=0`$ for TED very dominant.

The $`V_{Cm}`$, $`V_{CM}`$ and $`V_{BM}`$ contributions, shown separately in Fig. 2 for $`A=47`$ and 49 are quite far from the observed pattern, which is accurately reproduced only after these disparate terms are added. For $`A=50`$ and 51 we have replaced the $`V_{BM}`$ part by a variant of the full sum in which $`\beta _1`$ is halved. For $`A=50`$ the changes are insignificant, but there is a definite improvement in $`A=51`$ (remember again Eq. (7) and footnote before it).

It is especially worth noting in Fig. 2 that the strong signature effect in the $`A=49`$ band is erased in the MED by the out-of-phase $`V_{Cm}`$ and $`V_{CM}`$, while the signature staggering is enhanced in $`A=51`$.

The experimental TED patterns in Fig. 3 are quite nicely reproduced by the minimal $`\beta _{\mathrm{𝐩𝐟}}^{(2)}=\beta _2V_{𝐟_{\mathrm{𝟕}/\mathrm{𝟐}}}^{J=0}`$ choice. As mentioned, the inclusion of the $`J0`$ terms (third line of Table 1) makes little difference, and—interestingly enough—simply ignoring $`V_B`$ and doubling $`V_{CM}`$ (or the other way round) makes practically no difference. Which confirms the overwhelming dominance of $`J=0`$ pairing.

It can be hoped that a rigorous treatment calling upon state of the art CSB potentials mac01  will confirm the role of the $`J=2`$ pairing term for the isovector MED. The TED behaviour seems far simpler and our results are consistent with the findings in orm97  for $`\beta _𝐫^{(2)}`$ borne out in 46cr  for $`A=46`$. Therefore, here we may bet on—rather than hope for—confirmation by the charge independence breaking potentials mac01b .

The isovector channel raises a difficulty for $`A=46`$. In 46cr ; orm97  it was found that $`\beta _𝐫^{(1)}0`$ using the same functional form as for $`\beta _𝐫^{(2)}`$ with strong $`J=0`$ pairing, which does not square with our results. But in this case our results do not square with experiment either. The scheme that has been successful in $`A=47`$, 49, 50 and 51 fails in $`A=46`$: we are simply unable to do any better than in (46cr, , Fig. 3a)).

The trouble is no doubt due to the poor spectroscopy provided by full $`pf`$ shell diagonalizations for $`A46`$, at least when compared with the very high quality descriptions for the rest of the $`f_{7/2}`$ nuclei (i.e., $`A56`$). Fig. 4 illustrates the point: the calculated yrast energetics is wrong for the lowest states and, for the others, far less precise than the corresponding patterns in the heavier nuclei (see for example mar97  for $`A=47`$ and 49). This problem extends to transition rates and static moments. It was first noted and abundantly discussed in Ref. pov81  but its quantitative explanation remains a challenge. This unsatisfactory situation provides nonetheless a helpful clue: the TED may be unsensitive to details, but the MED demand accurate wave functions and could be taken as tests of their quality.

Within the $`A=46`$ proviso, our results make obvious something that may seem at first surprising: isospin non conserving potentials play a role that is at least as important as $`V_C`$ in explaining the MDE (and TDE, as found previously in 46cr ). In this respect, it is worth noting that direct evidence for charge symmetry breaking has been confined, so far, to the very light systems (basically $`A=2`$ and 3) mac01 . The mechanism plays an important part in resolving the Nolen Schiffer anomaly in the MDE, but the effects of $`V_C`$ remain much strongerduf02 . For the MED and TED, $`V_C`$ is at most as strong as $`V_B`$, for which we have shown that substantial quantitative information can be extracted from the data. To boot, the MED also provide a view of the evolution of yrast radii.

This work owes much to a stay of AZ at the UAM, made possible by a scholarship of the BBVA foundation. AP is supported by grant BFM2000-30 from MCyT (Spain).


text:
# 1 Introduction

## 1 Introduction

Given a set of positive integers $`A=\{a_1,\mathrm{},a_d\}`$ with $`\mathrm{gcd}(a_1,\mathrm{},a_d)=1`$, we call an integer $`n`$ *representable* if there exist nonnegative integers $`m_1,\mathrm{},m_d`$ such that

$$n=m_1a_1+\mathrm{}+m_da_d.$$

In this paper, we discuss the *linear diophantine problem of Frobenius*: namely, find the largest integer which is not representable. We call this largest integer the *Frobenius number* $`g(a_1,\mathrm{},a_d)`$.

One fact which makes this problem attractive is that it can be easily described, for example, in terms of coins of denominations $`a_1,\mathrm{},a_d`$; the Frobenius number is the largest amount of money which cannot be formed using these coins.

The following “folklore” theorem has long been known (probably at least since Sylvester ).

###### Theorem 1

$`g(a,b)=abab`$.

For $`d3`$, the quest for general formulas has so far been unsuccessful. For the case $`d=2`$, Sylvester  proved the following result.

###### Theorem 2 (Sylvester)

For $`A=\{a,b\}`$, exactly half of the integers between $`1`$ and $`(a1)(b1)`$ are representable.

Here we introduce and study a more general problem, a natural extension of the Frobenius problem.

###### Definition 1

We say that $`n`$ is *$`k`$-representable* if $`n`$ can be represented in the form

$$n=m_1a_1+\mathrm{}+m_da_d$$

in exactly $`k`$ ways.

In terms of coins, we can exchange the $`n`$ pennies in exactly $`k`$ different ways in terms of the given coin denominations. It is not hard to convince ourselves that—because the numbers in $`A`$ are relatively prime—for any $`k`$, eventually every integer can be represented in more than $`k`$ ways. Our extension of the Frobenius number is captured by the following definition:

###### Definition 2

$`g_k(a_1,\mathrm{},a_d)`$ is the smallest integer beyond which every integer is represented more than $`k`$ times.

This is a natural generalization of the concept of the Frobenius number, as $`g(a_1,\mathrm{},a_d)=g_0(a_1,\mathrm{},a_d)`$. As to be expected, the study of $`g_k`$ is extremely complicated for $`d3`$. There is an analogy here with $`k`$-representable integers and the classic problem of finding the number of representations of an integer as a sum of 4 squares, for example. However, the methods here are different. In this paper we concentrate on the case $`d=2`$ and present the following results. Let $`A=\{a,b\}`$.

###### Theorem 3

$`g_k(a,b)=(k+1)abab`$.

###### Theorem 4

Given $`k2`$, the smallest $`k`$-representable integer is $`ab(k1)`$.

###### Theorem 5

There are exactly $`ab1`$ integers which are uniquely representable. Given $`k2`$, there are exactly $`ab`$ $`k`$-representable integers.

Theorem 3 is a direct generalization of Theorem 1. Theorem 4 is meaningless for $`k=0`$ and trivial for $`k=1`$: the smallest representable integer is $`\mathrm{min}(a,b)`$. Theorem 5 extends Theorem 2 for $`k>0`$.

## 2 The restricted partition function

One approach to the Frobenius problem and its generalizations is through the study of the *restricted partition function*

$$p_A(n)=\mathrm{\#}\{(m_1,\mathrm{},m_d)^d:\text{ all }m_j0,m_1a_1+\mathrm{}+m_da_d=n\},$$

the number of partitions of $`n`$ using only the elements of $`A`$ as parts. In view of this function, $`g_k(a_1,\mathrm{},a_d)`$ is the smallest integer such that for every $`n>g_k(a_1,\mathrm{},a_d)`$ we have $`p_A(n)>k`$.

It is well known  that

$$p_A(n)=\frac{n^{d1}}{a_1\mathrm{}a_d(d1)!}+O\left(n^{d2}\right).$$

In particular,

$$p_{\{a,b\}}(n)=\frac{n}{ab}+c(n),$$

where $`c(n)=O(1)`$. In fact, \[11, p. 99\] gives a nice argument that $`c(n)`$ is periodic in $`n`$ with period $`ab`$, based on the generating function

$$\frac{1}{\left(1x^a\right)\left(1x^b\right)},$$

the coefficient of $`x^n`$ of which is equal to $`p_{\{a,b\}}(n)`$. This argument can be carried even further to give the following little-known formula.

###### Theorem 6 (Popoviciu)

Suppose $`a`$ and $`b`$ are relatively prime positive integers, and $`n`$ is a positive integer. Then

$$p_{\{a,b\}}(n)=\frac{n}{ab}\left\{\frac{b^1n}{a}\right\}\left\{\frac{a^1n}{b}\right\}+1.$$

Here $`\{x\}=xx`$ denotes the fractional part of $`x`$, and $`a^1a1(\text{mod }b)`$, and $`b^1b1(\text{mod }a)`$.

The earliest reference to this result that we are aware of is ; the formula has since been resurrected at least twice .

Instead of giving another proof of this theorem, we invite the reader to a scenic tour through the following modularized set of exercises. Consider the function

$$f(z)=\frac{1}{\left(1z^a\right)\left(1z^b\right)z^{n+1}}.$$

1. Compute the residues at all non-zero poles of $`f`$, and verify that Res$`(f(z),z=0)=p_{\{a,b\}}(n)`$.
2. Use the residue theorem to derive an identity for $`p_{\{a,b\}}(n)`$. (Integrate $`f`$ around a circle with center 0 and radius $`R`$, and show that this integral vanishes as $`R\mathrm{}`$.)
3. Verify that for $`b=1`$,

$`p_{\{a,1\}}(n)`$ $`=\mathrm{\#}\{(m_1,m_2):m_1,m_20,m_1a+m_2=n\}`$
$`=\mathrm{\#}\left([0,{\displaystyle \frac{n}{a}}]\right)={\displaystyle \frac{n}{a}}\left\{{\displaystyle \frac{n}{a}}\right\}+1.`$
4. Use this together with the identity found in 2. to obtain

$$\frac{1}{a}\underset{\lambda ^a=1\lambda }{}\frac{1}{(1\lambda )\lambda ^n}=\left\{\frac{n}{a}\right\}+\frac{1}{2}\frac{1}{2a}.$$
5. Verify that

$$\underset{\lambda ^a=1\lambda }{}\frac{1}{(1\lambda ^b)\lambda ^n}=\underset{\lambda ^a=1\lambda }{}\frac{1}{(1\lambda )\lambda ^{b^1t}}$$

and use this together with 4. above to simplify the identity found in 2.

Popoviciu’s beautiful and simple formula leads to very short proofs of the results stated in the introduction.

*First proof of Theorems 1 and 3.* We will show that $`p_{\{a,b\}}((k+1)abab)=k`$ and that $`p_{\{a,b\}}(n)>k`$ for every $`n>(k+1)abab`$. First, by the periodicity of $`\{x\}`$,

$`p_{\{a,b\}}((k+1)abab)=`$
$`={\displaystyle \frac{(k+1)abab}{ab}}\left\{{\displaystyle \frac{b^1((k+1)abab)}{a}}\right\}\left\{{\displaystyle \frac{a^1((k+1)abab)}{b}}\right\}+1`$
$`=k+2{\displaystyle \frac{1}{b}}{\displaystyle \frac{1}{a}}\left\{{\displaystyle \frac{b^1b}{a}}\right\}\left\{{\displaystyle \frac{a^1a}{b}}\right\}`$
$`=k+2{\displaystyle \frac{1}{b}}{\displaystyle \frac{1}{a}}\left\{{\displaystyle \frac{1}{a}}\right\}\left\{{\displaystyle \frac{1}{b}}\right\}`$
$`=k+2{\displaystyle \frac{1}{b}}{\displaystyle \frac{1}{a}}\left(1{\displaystyle \frac{1}{a}}\right)\left(1{\displaystyle \frac{1}{b}}\right)=k.`$
For any integer $`m`$, $`\left\{\frac{m}{a}\right\}1\frac{1}{a}`$. Hence for any positive integer $`n`$,

$$p_{\{a,b\}}((k+1)abab+n)\frac{(k+1)abab+n}{ab}\left(1\frac{1}{a}\right)\left(1\frac{1}{b}\right)+1=k+\frac{n}{ab}>k.$$

$`\mathrm{}`$

*Proof of Theorem 2.* We first claim that, if $`n[1,ab1]`$ is not a multiple of $`a`$ or $`b`$,

$$p_{\{a,b\}}(n)+p_{\{a,b\}}(abn)=1.$$
(1)
This identity follows directly from Theorem 6:

$`p_{\{a,b\}}(abn)={\displaystyle \frac{abn}{ab}}\left\{{\displaystyle \frac{b^1(abn)}{a}}\right\}\left\{{\displaystyle \frac{a^1(abn)}{b}}\right\}+1`$
$`=2{\displaystyle \frac{n}{ab}}\left\{{\displaystyle \frac{b^1n}{a}}\right\}\left\{{\displaystyle \frac{a^1n}{b}}\right\}`$
$`\stackrel{()}{=}{\displaystyle \frac{n}{ab}}+\left\{{\displaystyle \frac{b^1n}{a}}\right\}+\left\{{\displaystyle \frac{a^1n}{b}}\right\}`$
$`=1p_{\{a,b\}}(n).`$
Here, $`()`$ follows from the fact that $`\{x\}=1\{x\}`$ if $`x`$. This shows that, for $`n`$ between 1 and $`ab1`$ and not divisible by $`a`$ or $`b`$, exactly one of $`n`$ and $`abn`$ is not representable. There are

$$abab+1=(a1)(b1)=g(a,b)+1$$

integers between 1 and $`ab1`$ which are not divisible by $`a`$ or $`b`$. Finally, we note that $`p_{\{a,b\}}(n)>0`$ if $`n`$ is a multiple of $`a`$ or $`b`$, by the very definition of $`p_{\{a,b\}}(n)`$. Hence the number of non-representable integers is $`\frac{1}{2}(a1)(b1)`$. $`\mathrm{}`$

Note that we proved even more. By (1), every positive integer less than $`ab`$ has at most one representation. Hence, the representable integers in the above theorem are *uniquely* representable.

*First proof of Theorem 4.* Let $`n`$ be a nonnegative integer. Then

$`p_{\{a,b\}}(ab(k1)n)=`$
$`={\displaystyle \frac{ab(k1)n}{ab}}\left\{{\displaystyle \frac{b^1(ab(k1)n)}{a}}\right\}\left\{{\displaystyle \frac{a^1(ab(k1)n)}{b}}\right\}+1`$
$`=k{\displaystyle \frac{n}{ab}}\left\{{\displaystyle \frac{b^1n}{a}}\right\}\left\{{\displaystyle \frac{a^1n}{b}}\right\}.`$ (2)
If $`n=0`$, (2) equals $`k`$. If $`n`$ is positive, we use $`\{x\}0`$ to see that

$$p_{\{a,b\}}(ab(k1)n)k\frac{n}{ab}<k.$$

$`\mathrm{}`$

All nonrepresentable positive integers lie, by definition, in the interval $`[1,g(a,b)]`$. It is easy to see that the smallest interval containing all uniquely representable integers is $`[\text{min}(a,b),g_1(a,b)]`$. For $`k2`$, the corresponding interval always has length $`2abab+1`$, and the precise interval is given next.

###### Corollary 7

Given $`k2`$, the smallest interval containing all $`k`$-representable integers is

$$[g_{k2}(a,b)+a+b,g_k(a,b)].$$

*Proof.* By Theorems 3 and 4, the smallest integer in the interval is

$$ab(k1)=g_{k2}(a,b)+a+b.$$

The upper bound of the interval follows from the proof of Theorem 3. $`\mathrm{}`$

*First proof of Theorem 5.* First, in the interval $`[1,ab]`$, there are, by Theorems 2 and 4,

$$ab\frac{(a1)(b1)}{2}1$$

1-representable integers. Because of the periodic behavior of the partition function

$$p_{\{a,b\}}(n+ab)=p_{\{a,b\}}(n)+1,$$
(3)
which follows directly from Theorem 6, we see that there are

$$\frac{(a1)(b1)}{2}$$

1-representable integers above $`ab`$. For $`k2`$, the statement follows by similar reasoning. $`\mathrm{}`$

## 3 Elementary proofs

Although the proofs we have given so far are simple, they rely on Popoviciu’s formula (Theorem 6). In this section, we show how one can use an alternate recursive method. We will still use the philosophy behind the restricted partition function, but not its concrete form as given by Theorem 6. The basic idea is to view

$$p_A(n)=\mathrm{\#}\{(m_1,\mathrm{},m_d)^d:\text{ all }m_j0,m_1a_1+\mathrm{}+m_{d1}a_{d1}=nm_da_d\},$$

from which one obtains the recursion formula

$$p_{\{a_1,\mathrm{},a_d\}}(n)=\underset{m0}{}p_{\{a_1,\mathrm{},a_{d1}\}}(nma_n).$$
(4)
For the following, it will be useful to introduce the function

$$q_A(n)=\mathrm{\#}\{(m_1,\mathrm{},m_d)^d:\text{ all }m_j>0,m_1a_1+\mathrm{}+m_da_d=n\},$$

which counts those partitions of $`n`$ which uses only parts from $`A`$, where we additionally demand that each part gets used at least once. The functions $`p_A`$ and $`q_A`$ are intimitely related through

$$q_A(n)=p_A(na_1\mathrm{}a_d).$$
(5)
The recursion formula for $`q_A`$, corresponding to (4), is

$$q_{\{a_1,\mathrm{},a_d\}}(n)=\underset{m>0}{}q_{\{a_1,\mathrm{},a_{d1}\}}(nma_d).$$

*Second proof of Theorems 1 and 3.* For $`d=2`$, the last recursion simplifies to

$$q_{\{a,b\}}(n)=\underset{m>0}{}q_{\{a\}}(nmb).$$
(6)
Now

$$q_{\{a\}}(n)=\{\begin{array}{cc}1& \text{ if }a|n\text{ and }n>0\hfill \\ 0& \text{ otherwise. }\hfill \end{array}$$

Since $`\mathrm{gcd}(a,b)=1`$, the sum in (6) is larger than $`k`$ if we have $`(k+1)a`$ summands or more. There are $`\frac{n1}{b}`$ summands in (6). Hence for any $`n>(k+1)ab`$, $`q_{\{a,b\}}(n)>k`$. On the other hand, it is easy to see that $`q_{\{a,b\}}((k+1)ab)=k`$. Via (5), this translates into

$`p_{\{a,b\}}(n)>k\text{ if }n>(k+1)abab`$
$`p_{\{a,b\}}((k+1)abab)=k`$
$`\mathrm{}`$

*Second proof of Theorem 4.* We play a similar game as in the previous proof, now starting with (4), which gives for $`d=2`$

$$p_{\{a,b\}}(n)=\underset{m0}{}p_{\{a\}}(nmb).$$
(7)
Now

$$p_{\{a\}}(n)=\{\begin{array}{cc}1& \text{ if }a|n\text{ and }n0\hfill \\ 0& \text{ otherwise, }\hfill \end{array}$$

and since $`\mathrm{gcd}(a,b)=1`$, multiples of $`ab`$ give ‘peaks’ for the sum in (7) in the sense that

$$p_{\{a,b\}}(kab)>p_{\{a,b\}}(n)\text{ for all }n<kab.$$

The proof is finished by the easy observation that

$$p_{\{a,b\}}(kab)=k+1.$$

$`\mathrm{}`$

*Second proof of Theorem 5.* The first proof follows verbatim; the only place we have used the precise form of $`p_{\{a,b\}}`$ was for the periodicity property (3). But this identity follows easily from (7) and the fact that $`a`$ and $`b`$ are relatively prime. $`\mathrm{}`$

## 4 Final remarks

We note that for all $`d>2`$, generalized Dedekind sums  appear in the formulas for $`p_A(n)`$, which increases the complexity of the problem. The full details of these connections to Dedekind sums appear in .

We conclude with a few remarks regarding extensions of the above theorems to $`d>2`$. Although no ‘nice’ formula similar to the one appearing in Theorem 1 is known for $`d>2`$, there has been a huge effort devoted to giving bounds and algorithms for the Frobenius number . Secondly, we remark that Theorem 2 does not extend in general; however,  gives necessary and sufficient conditions on the $`a_j`$’s under which Theorem 2 does extend. The periodic behavior (3) of the partition function extends easily to higher dimensions . We leave the reader with the following “exercise”:

Unsolved problems: *Extend Theorems 3, 4, and 5 to $`d3`$.*

Acknowledgements. The authors are grateful to Mel Nathanson and Herb Wilf for helpful remarks and references on partition functions, and to Olivier Bordellès for pointing out the reference .

Department of Mathematical Sciences
Binghamton University
Binghamton, NY 13902-6000
matthias@math.binghamton.edu

Department of Mathematics
Temple University
Philadelphia, PA 19122
srobins@math.temple.edu