A Line Profile with Rayleigh Scattering

# The Formation of the First Stars II. Radiative Feedback Processes and Implications for the Initial Mass Function DRAFT: March 6, 2018

## Abstract

We consider the radiative feedback processes that operate during the formation of the first stars, including the photodissociation of H, Lyman- radiation pressure, formation and expansion of an H ii region, and disk photoevaporation. These processes may inhibit continued accretion once the stellar mass has reached a critical value, and we evaluate this mass separately for each process. Photodissociation of H in the local dark matter minihalo occurs relatively early in the growth of the protostar, but we argue this does not affect subsequent accretion since by this time the depth of the potential is large enough for accretion to be mediated by atomic cooling. However, neighboring starless minihalos can be affected. Ionization creates an H ii region in the infalling envelope above and below the accretion disk. Lyman- radiation pressure acting at the boundary of the H ii region is effective at reversing infall from narrow polar directions when the star reaches , but cannot prevent infall from other directions. Expansion of the H ii region beyond the gravitational escape radius for ionized gas occurs at masses , depending on the accretion rate and angular momentum of the inflow. However, again, accretion from the equatorial regions can continue since the neutral accretion disk has a finite thickness and shields a substantial fraction of the accretion envelope from direct ionizing flux. At higher stellar masses, in the fiducial case, the combination of declining accretion rates and increasing photoevaporation-driven mass loss from the disk act to effectively halt the increase in the protostellar mass. We identify this process as the mechanism that terminates the growth of Population III stars (i.e., stars with primordial composition) that have not been affected by prior star formation (Population III.1 stars). We discuss the implications of our results for the initial mass function of these stars. In the Appendix we develop approximate solutions to a number of problems relevant to the formation of the first stars: the effect of Rayleigh scattering on line profiles in media of very large optical depth; the intensity of Lyman- radiation in very opaque media; an approximate determination of the radiative acceleration in terms of the gradient of a modified radiation pressure; the determination of the flux of radiation in a shell with an arbitrary distribution of opacity; and the vertical structure of an accretion disk supported by gas pressure with constant opacity.

stars: formation — early universe — cosmology: theory

## 1 Introduction

There has been substantial recent progress in our theoretical understanding of how the first stars formed (Bromm & Larson 2004). In marked contrast to the case for contemporary star formation, the initial conditions for the formation of the first stars are believed to be relatively well understood: they are determined by the growth of small-scale gravitational instabilities from cosmological fluctuations in a cold dark matter universe. The first stars are expected to form at redshifts in dark matter halos of mass (Tegmark et al. 1997). In the absence of any elements heavier than helium (other than trace amounts of lithium) the chemistry and thermodynamics of the gas are very simple (Abel, Bryan, & Norman 2002, hereafter ABN; Bromm, Coppi, & Larson 2002, hereafter BCL02). There are no dust grains to couple the gas to radiation emitted by the protostar. There are no previous generations of stars to roil the gas out of which the stars form, nor is there any radiation other than the cosmic background radiation. Existing calculations have assumed that there were no significant primordial magnetic fields, thereby eliminating a major complication that occurs in contemporary star formation. However, even in the absence of significant primordial fields, it is possible that magnetic fields could have been generated in the accretion disk surrounding a primordial protostar (Tan & Blackman 2004), although even in this case the magnetic fields become dynamically important only after the star formation process is well underway, and they do not affect the initial conditions. Given this relative simplicity, there is some confidence in the results of numerical simulations that have followed the collapse of cosmological scale perturbations down to almost stellar dimensions (ABN; Bromm, Coppi, & Larson 1999; Yoshida et al. 2006; O’Shea & Norman 2007). This confidence is strengthened by the fact that it appears to be the microphysics of cooling that determines the types of baryonic structures that are formed, and not, for example, the details of the initial power spectrum of fluctuations in dark matter density. The results of these simulations suggest that the initial gas cores out of which stars form are quite massive, .

The observational imprint of the first stars depends on their mass: These stars were likely to be of critical importance in reionizing the universe, in producing the first metals, and in creating the first stellar-mass black holes. The number of ionizing photons emitted per baryon depends on the stellar mass for (Bromm, Kudritzki, & Loeb 2001). The hardness of the radiation field is also sensitive to the mass (e.g. Tumlinson & Shull 2000; Schaerer 2002), so that He reionization can be affected. The effectiveness of the first stars in enriching the intergalactic medium (IGM) with metals and in producing the first stellar-mass black holes also depends sensitively on the mass of the star. A potential simplification in assessing these effects is that massive primordial stars are thought to have much smaller mass-loss rates than contemporary massive stars (Kudritzki 2002), so that the mass at core collapse should be quite similar to the initial mass. However, Meynet, Ekström, & Maeder (2006) have argued that if rotation is allowed for, then mass loss can be significant. Heger & Woosley (2002) showed that stars exploding as supernovae above about and between 40 and should collapse directly to black holes, and they argued that such stars would providing relatively little metal enrichment. However, Ohkubo et al. (2007) followed the collapse and explosion of and stars in two dimensions and concluded that substantial amount of metals could be ejected; they proposed that such supernovae could produce intermediate-mass black holes. Heger & Woosley (2002) also showed that for , the pair instability leads to explosive O and Si burning that completely disrupts the star, leaving no remnant and ejecting large quantities of heavy elements. Such supernovae produce a dramatic odd-even effect in the nuclei produced. Stars below are expected to form neutron stars, with more normal enrichment rates. In principle, metallicity determinations from high-redshift absorption line systems (Schaye et al. 2003; Norman, O’Shea, & Paschos 2004) or from very metal poor local stars (Beers & Christlieb 2005) can constrain the initial mass function (IMF) of the early generations of stars. Indeed, based on observations such as these, Daigne et al. (2004) argue that the stars responsible for reionizing the universe mostly likely had masses , and Tumlinson (2006) concludes that stars above could have contributed at most 10% of the iron observed in extremely metal poor stars (those with [Fe/H] —Beers & Christlieb 2005).

In discussing the formation of the first stars, the terms “first stars” and “Population III stars” are often used interchangeably, but this can lead to confusion. To be precise, we shall follow the conventions suggested by one of us at the First Stars III conference (O’Shea et al. 2008) and define Population III stars as those stars with a metallicity sufficiently low that it has no effect on either the formation or the evolution of the stars. The value of the critical metallicity for star formation—i.e., the value below which the metals do not influence star formation—is uncertain, with estimates ranging from if the cooling is dominated by small dust grains that contain a significant fraction of the metals (Omukai et al. 2005) to if the cooling is dominated by the fine structure lines of C and O and there is negligible H (Bromm & Loeb 2003); if H cooling is included, Jappsen et al. (2007) argue that there is no critical metallicity for gas-phase metals. It is possible that values of the metallicity even less than could significantly affect the evolution of primordial stars (Meynet, private communication). Among Population III stars, we distinguish between the first and second generations, termed Population III.1 and III.2, respectively: The initial conditions for the formation of Population III.1 stars are determined solely by cosmological fluctuations whereas those for Population III.2 stars are significantly affected by other stars. It is likely that Population III.1 stars have a primordial composition, since it is hard to see how the gas out of which they form could be contaminated by even small trace amounts of metals without having been affected by radiation from the star that produced the metals. Stars affect the initial conditions for the formation of Population III.2 stars primarily by radiation, both ionizing radiation and Lyman-Werner band radiation that destroys H molecules. The latter reduces the cooling efficiency of the gas, allowing compression to heat up to the point that it begins to become collisionally ionized; shocks associated with H ii regions and supernova remnants can also ionize the gas. Once the gas has been partially ionized, HD cooling can become important, reducing the characteristic star-formation mass (Uehara & Inutsuka 2000). Greif & Bromm (2006) use the term “Population II.5” to describe stars that form from gas in which HD cooling is important, whereas in our terminology these would be Population III.2 stars, but we believe that it is better to describe all essentially metal-free stars as “Population III.” It should be noted that our definition of Population III.2 stars includes all Population III stars that were significantly affected by previous generations of star formation, even if that did not result in significant HD production. Those Population III.2 stars that form out of gas that is enriched in HD will typically be less massive than Population III.1 stars—by about a factor 10 according to Greif & Bromm (2006). They infer that Population III.1 stars are relatively rare, constituting about 10% by mass of all Population III stars.

In this paper we wish to estimate the characteristic mass of the first generation of stars (Population III.1). Even if they are relatively rare, they are critical in setting the initial conditions for the star formation that followed, and therefore in determining the reionization of the universe and the production of the first metals and the first stellar mass black holes. For contemporary star formation, it is believed that the IMF is set by a combination of gravitational fragmentation in a turbulent medium (Elmegreen 1997; Padoan & Nordlund 2002) and feedback effects. The characteristic stellar mass is of order the Bonnor-Ebert mass, . However, not all of the initial core mass is incorporated into the final star, since contemporary protostars have powerful outflows that eject some of the core mass (Nakano, Hasegawa, & Norman 1995; Matzner & McKee 2000). There are a number of feedback effects that occur in contemporary massive star formation (Larson & Starrfield 1971), particularly radiation pressure on dust and photoionization associated with the growth of an H ii region. It remains unclear whether the upper limit on the contemporary IMF is set by feedback or by instabilities that afflict very massive stars.

The clumps out of which the first stars form have total masses, including dark matter, of order ; these objects are typically referred to as minihalos. Cooling by trace amounts of H generally leads to the formation of a gravitationally unstable core of baryonic mass (ABN). In contrast to contemporary star-forming regions, the turbulence in this gas is subsonic, due to the higher temperatures and the lack of internal and external sources of turbulence. As a result, gravitational fragmentation is much less effective—indeed, numerical simulations show no evidence for it (e.g., ABN, Yoshida et al. 2006, but see Clark, Glover, & Klessen 2008), and analytic calculations (Ripamonti & Abel 2004), including those that consider disk fragmentation (Tan & Blackman 2004), show no evidence for fragmentation either. It therefore appears that the mass of the first stars is likely to be set by feedback effects.

Feedback effects can be classified as either radiative or kinetic. Kinetic feedback includes protostellar outflows and main-sequence winds. In contemporary star formation, outflows are believed to be hydromagnetically driven; in this paper, we assume that the magnetic fields associated with the protostar are too weak or too tangled to drive a strong outflow (see Tan & Blackman 2004 for a more extensive discussion of the effect of these outflows). Due to the absence of metals, main-sequence winds are very weak in the absence of rotation (Kudritzki 2002), although they could become important in the later stages of evolution if the star is rapidly rotating. Since we are primarily interested in the early stages of evolution of the star, we neglect kinetic feedback (i.e., outflows and winds).

Radiative feedback includes radiation pressure, photoionization heating, and photodissociation. Radiation pressure can be due to continuum radiation or to resonance line scattering; the continuum radiation pressure can be due to electron scattering or to photoionization. Photoionization also leads to heating, which unbinds gas beyond the gravitational radius , where is the isothermal sound speed of the gas. If the gas is initially in a disk, this process is termed photoevaporation. Finally, photodissociation destroys H  the dominant coolant in neutral primordial gas.

Most previous work has focused on the effects of continuum radiation pressure and photoionization heating in limiting the mass of primordial stars. Omukai & Palla (2001; 2003) focused on electron scattering, which leads to the Eddington limit on an accreting mass. For the case of spherical accretion at a rate of , radiation pressure first becomes important at around , leading to a dramatic swelling of the stellar surface. This, however, is a transient effect because an important part of the luminosity is due to accretion, which is reduced by the increase in the stellar radius. Only at masses around does the internal luminosity become very close to the Eddington value, leading to runaway expansion of the star and, presumably, the end of accretion. Omukai & Palla (2003) considered a range of accretion rates. They found that if the accretion rate is smaller than , then the total luminosity remains sub-Eddington and the star continues to grow along the main sequence to arbitrarily large masses. On the other hand, if the accretion rate is somewhat larger than this critical rate, the Eddington limit becomes important at around . Accretion at a rate based on the settling motions in the core of ABN is slow enough that the Eddington limit does not affect the final mass. However, these models ignored the influence of other protostellar feedback mechanisms on the infalling envelope. These models also assumed spherical symmetry, which leads to much larger photospheric radii and thus a softer radiation field than in the more realistic case of disk accretion (see §7 in Paper I).

Omukai & Inutsuka (2002) considered the the combined effects of photoionization heating and continuum radiation pressure due to photoionization in the infalling envelope. They show that there is a critical stellar mass at which the hydrogen ionizing luminosity is sufficient to create an H ii region, which rapidly spreads to large distances where its thermal pressure becomes dynamically important in slowing infall. However, the ionizing radiation force decelerates the inflowing gas, raising the gas density and therefore reducing the radius of the H ii region. For spherical inflow, this mechanism is so effective that the radius of the H ii region remains well below the gravitational radius , stopping any mass loss. They concluded that with this effect, there was no limit on protostellar masses below . Without this radiation force, they predicted a mass limit of order . As we shall see below, these conclusions are sensitive to the assumption of spherical accretion.

In Paper I (Tan & McKee 2004) we modeled the growth of a primordial protostar from very small masses to large. We included the effects of rotation of the infalling gas, which led to the formation of an accretion disk around the protostar. The goal of this paper is to determine when the energy output from the protostar is sufficient to halt accretion and set the final stellar mass. This is an extremely complicated problem, the full solution of which requires three dimensional hydrodynamical simulations that include the generation, propagation, and dynamical influence of radiation. Furthermore these simulations must resolve a large range of scales: the protostar is of order , while the size of the quasi-hydrostatic core that encloses is of order 1 pc, several million times greater. The demands on the time scale are even greater: the simulation may have to follow the evolution of the star over its lifetime Myr (Schaerer 2002) while at the same time following the dynamics of an accretion disk with a characteristic time scale as short as s. The numerical simulations of ABN were able to resolve an even greater range of radii, but it will be some years before it is possible to meet the required dynamical range in time scales. As a result, we shall develop simple analytic models for the feedback interaction that we hope will provide a useful first step.

We begin our discussion with a review the results of Paper I in §2. Feedback effects are then considered in the approximate order in which they become manifest. In §3, we briefly discuss the effects of photodissociation. Lyman- radiation pressure feedback, is considered in §4 and in several Appendixes. The feedback from ionizing photons that can create an H ii region is considered in §5. After discussing shadowing by accretion disks in §6 and an Appendix, the feedback due to disk photoevaporation is considered in §7. Figure 1 gives an overview of these feedback processes occurring near the protostar. Finally, our conclusions are summarized in §8.

## 2 Review of Paper I: Properties and Evolution of Primordial Protostars

The radiative output from a protostar depends on the temperature and luminosity of its emitting components, which are the star itself (stellar photosphere), the boundary layer of the accretion disk with the star, and the larger scale accretion disk. The luminosity of the star depends mostly on its mass. The size of the star then determines its temperature. The size of the star and the accretion rate determine the radiative properties of the both the boundary layer and the accretion disk, since emission from the latter is dominated by the inner regions.

The size of the star depends on the accretion rate during its formation history. At lower masses there is a balance in the size set by the need to radiate the luminosity, which is mostly due to accretion, with a photospheric temperature that is likely to be close to  K because the opacity due to rapidly declines in this temperature regime. Under the assumption of spherical accretion, Stahler, Palla, & Salpeter (1986) found the protostellar radius to be

 r∗≃90.8m0.27∗,2˙m0.41∗,−3   R⊙     (m∗,2≲0.1), (1)

where and . For the accretion rates typical of primordial star formation we see that the size is very large. For more massive protostars there is a transition once the star is about as old as its local Kelvin-Helmholtz time, then contraction proceeds towards the main sequence size, where accretion can continue. In Paper I, we found that for typical conditions, the protostar reached its main sequence radius at about . According to Schaerer (2002), this radius is

 r∗≃4.3m0.55∗,2 R⊙     (main sequence) (2)

to within 6% for .

Thus almost all the radiative stellar properties depend on just two parameters: the stellar mass and the accretion rate. Note that in principle these properties also depend on the angular momentum of the infalling gas, since if there was no rotation, then spherical accretion implies very high gas densities near the protostar, affecting the location of its photosphere. However, for any realistic amount of angular momentum, a disk forms whose size is much larger than , and the star’s properties no longer depend on the rotation of the core from which the star forms.

The accretion rate of Population III protostars depends on the density structure of the gas core at point when the star starts to form. This density structure is set by the balance of thermal pressure and self-gravity, which in turn depends primarily on the cooling properties of molecular hydrogen. This cooling creates almost isothermal cores at  K with an outer bounding density of about , which is the critical density of cooling transitions (for molecules interacting with atomic H). In fact the temperatures increase to several hundred K in the inner parts of the core because of the reduced cooling efficiency above the critical density. These basic features have been confirmed by numerical studies (ABN; BCL02; O’Shea & Norman 2007). The trigger for dynamical collapse is thought to be the rapid formation of by three-body collisions at high densities  cm, since this then dramatically increases the cooling rate in this region.

ABN carried their calculations almost to the point of protostar formation, and at this time gas was flowing inward subsonically almost everywhere (except for , where the inflow was slightly supersonic). Shu’s (1977) expansion wave solutions for protostellar accretion are based on the assumption that the inflow velocity at this time is zero. Hunter (1977) generalized these solutions and showed that there is a discrete set of self-similar solutions that begin at rest at and have a constant infall velocity at the time of protostar formation (). One of these solutions, the Larson-Penston solution (Larson 1969; Penston 1969), has supersonic inflow (Mach number = 3.3 at ; this solution is clearly inconsistent with the numerical results. In fact, the accretion flow appears to be a settling solution regulated by H cooling. Only one of Hunter’s solutions corresponds to mildly subsonic inflow (Mach number =0.295 at ), comparable to that found by ABN, and this is the solution adopted in Paper I. This solution has a density that is 1.189 times greater than a singular isothermal sphere at , and the accretion rate is 2.6 times greater.

Hunter’s (1977) solution applies to an isothermal gas. Omukai & Nishi (1998) and Ripamonti et al. (2002) have numerically calculated accretion rates for primordial protostars, and showed that the accreting gas is isentropic with an adiabatic index due to H cooling; i.e., each mass element satisfies the relation with the “entropy parameter”  const. In hydrostatic equilibrium, such a gas settles into a polytropic configuration, which in general has . For an isentropic gas, we have and . In Paper I, we presented an analytic model for the protostellar accretion rate for isentropic gas. We allowed for the existence of an accretion disk around the protostar with a significant fraction of the stellar mass,

 m∗d=m∗+md≡(1+fd)m∗, (3)

with a fiducial value for the disk mass fraction . Following McKee & Tan (2002, 2003), we wrote the accretion rate as

 ˙m∗d=ϕ∗m∗dtff, (4)

where is a numerical constant of order unity and is the free-fall time of the gas (measured at ). For gas that is in hydrostatic equilibrium at , McKee & Tan (2002) showed that to within about 1% for ; we have since verified that this is valid for . To our knowledge, Hunter’s self-similar solutions starting at have not been generalized to the non-isothermal case. 1 In Paper I, we therefore assumed that the accretion rate for the case is 2.6 times greater than that for hydrostatic initial conditions, just as in the isothermal case.

Feedback from the star, whether due to winds, photoionization, or radiation pressure, can reduce the accretion rate onto the star. We define a hypothetical star-disk mass, , and accretion rate, , in the absence of feedback. In this case, the star-disk mass equals the mass of the core out of which it was formed, . The instantaneous and mean star formation efficiencies are

 ϵ∗d ≡ ˙m∗d˙m∗d,0, (5) ¯ϵ∗d ≡ m∗dm∗d,0=m∗dM. (6)

In our previous work, we assumed that the star formation efficiency was constant, so that . in the present work, we shall find that significant feedback does not set in until the star is fairly massive, so that we must distinguish the instantaneous and mean values. In Paper I, we found that the accretion rate onto the star-disk system is

 ˙m∗d=0.026[ϵ∗dK′15/7(M/M⊙)3/7]     M⊙ yr−1, (7)

where

 K′≡P/ργ1.88×1012 cgs=(T′eff300 K)(104 cm−3nH)0.1 (8)

is a measure of the entropy of the accreting gas. Here is an effective temperature that includes the effect of turbulent motions; we have added a prime to the defined in Paper I to distinguish it from the effective temperature of a radiating atmosphere. Expressing the accretion rate in terms of the stellar mass, which equations (3) and (6) imply is , we find

 ˙m∗d=0.026⎡⎣ϵ∗d¯ϵ3/7∗dK′15/7(1+fd)3/7(m∗/M⊙)3/7⎤⎦     M⊙ yr−1, (9)

With , this result is in good agreement with the results of Omukai & Nishi (1998) and Ripamonti et al. (2002); since their 1D calculations did not allow for disks, this comparison is made for . Note that this agreement validates our use of the Hunter mildly subsonic solution to infer the accretion rate. If (i.e., no feedback) and then the accretion rate onto the star disk is

 ˙m∗d,−3→3.20⎛⎜⎝K′15/7m3/7∗,2⎞⎟⎠, (10)

where henceforth it will be understood that numerical evaluations denoted by “” have and . In this case the accretion rate onto the star (which may be primarily through the disk) is 3/4 of this [since ].

Our estimate of the accretion rate is somewhat above that estimated by ABN, but this is to be expected since their calculation stopped prior to the formation of the protostar. Indeed, at the time at which the protostar first forms (), the accretion rate at any finite radius [i.e. ] in a self-similar, isothermal collapse is smaller than the value it has at times , where is the isothermal sound speed. Equivalently, at a given time, the accretion rate at radii is less than that at small radii, . In Shu’s (1977) solution for the collapse of a singular isothermal sphere, the accretion rate at a given time is zero outside the expansion wave at ; inside the expansion wave, the accretion rate smoothly increases to as . For the Larson-Penston solution, the accretion rate at a given time increases from at large radii () to at small radii (). For Hunter’s mildly subsonic solution, which we have suggested is closest to the numerical simulations, the accretion rate increases from at large radii to at small radii (Hunter 1977), an increase of a factor 3.7. This demonstrates that caution should be exercised in inferring accretion rates at late times from those measured at early times, which is commonly done in simulations (e.g., ABN, Yoshida et al. 2006, O’Shea and Norman 2007).

The age of the star when it reaches a mass is (Paper I)

 tyr=27(1+fd¯ϵ∗d)10/7K′−15/7(m∗M⊙)10/7→2.92×104K′−15/7m10/7∗,2, (11)

where  yr) and where it is the mean star formation efficiency that enters. The resulting stellar mass is

 m∗→0.075K′1.5t0.7yr  M⊙. (12)

Bromm & Loeb (2004) have carried out a 3D simulation of the accretion onto the protostar for the first yr, and for our result is within a factor of theirs for this time period. (However, it should be noted that an extrapolation of their result to times beyond yr gives a mass less than half our estimate of the mass of the star plus disk. It remains to be determined whether such an extrapolation is valid.)

With this estimate of the protostellar mass, it is possible to calculate the maximum possible mass a primordial star could have. Schaerer (2002) has calculated the main sequence lifetimes of primordial stars with no mass-loss for , and his results are accurately described by the expression  Myr for . If we assume that the accretion is not limited by any feedback processes (), that Schaerer’s results can be extrapolated to higher masses, and that accretion during the relatively short post main-sequence phase is negligible, then we find

 m∗,max=∫tms0˙m∗dt≃1900(11+fd)0.86K′1.28   M⊙→1500K′1.28   M⊙. (13)

The maximum possible stellar mass is thus controlled by the value of the entropy parameter of the core.

In Paper I we also considered the effect of rotation. Rotation of the infalling gas has a dramatic effect on the evolution of the protostar, since it leads to much smaller photospheric radii and correspondingly higher temperatures and accretion luminosities. We parameterized the rotation in terms of

 fKep≡vrot(rsp)vKep(rsp), (14)

the ratio of the rotational velocity to the Keplerian velocity measured at the sonic point at . ABN found independent of radius, so we adopt this as a fiducial value. We then showed that the accreting gas formed a disk with an outer radius

 rd = 1.92×1016(fKep0.5)2(m∗d,2¯ϵ∗d)9/7K′−10/7     cm, (15) → 2.78×1016(fKep0.5)2m9/7∗,2K′10/7     cm. (16)

## 3 Photodissociation Feedback

As the protostar grows in mass, it begins to emit copious amounts of non-ionizing ultraviolet radiation (far-ultraviolet, or FUV, radiation), as shown in Figure 2. This radiation can photodissociate the H that is critical for cooling the accreting gas (Omukai & Nishi 1999; Glover & Brand 2001); its dynamical effects are considered in the next section.

Once the molecular coolants in the accreting gas are destroyed, the adiabatic index rises from to . If the gas were able to continue contraction, it would heat up until the temperature became high enough ( K) to excite the Ly- transition. For K, the adiabatic index would then drop back to about 1.

Can the protostar continue to accrete when ? If one considers the related problem of the gravitational stability of polytropic gas spheres, one might be led to conclude that accretion would stop: polytropic stars are stable against gravitational collapse for (Chandrasekhar 1939), and the same is true for polytropic gas clouds even if (McKee & Holliman 1999). However, there is a crucial distinction between collapse onto a protostar and the contraction of a gas cloud prior to protostar formation, and that is the presence of the central protostar, which is effectively a mass singularity. The stability analyses cited above assumed that there was no mass singularity at the origin. When one is present, the problem is analogous to that of Bondi accretion, the accretion of non-self gravitating gas onto a star; this can occur for (e.g., Shapiro & Teukolsky 1983). The problem of protostellar accretion, which includes the self-gravity of the gas, has been considered by Fatuzzo et al. (2004) for a wide range of conditions. They showed that gravity dominates over pressure close to the protostar, so that accretion can occur, provided . It is straightforward to see why: In supersonic inflow, the density scales as , so that the temperature rises more slowly than the kinetic energy per unit mass provided . They demonstrated that the accretion rate for a singular, initially isothermal sphere with is only slightly smaller than for the case in which .

The argument of Fatuzzo et al. applies to the inner, supersonic region of infall. What about the outer, subsonic region? There the density varies as a higher power of radius (e.g., for = const., ), so that pressure can overcome gravity at a lower value of . To see this more quantitively, consider the case of a singular isothermal sphere with . Assume that the initial density of the sphere is times greater than the equilibrium value. If the gas is flowing inward at a velocity far from the protostar (i.e., at large values of the similarity variable , which is just in the isothermal case), then varies as

 v=−v∞−2(Λ−1)x−(4−3γ)v∞x2+⋯ . (17)

(Fatuzzo et al. 2004; we have corrected a typo in the last factor). Thus, for , pressure forces will tend to decelerate the flow; however, this can be overcome by a suitable overdensity . We have confirmed this by numerical integration of the equations given by Fatuzzo et al.: For and , accretion is possible for ; for , accretion is possible provided exceeds some threshold. For primordial star formation, we estimate and ; accretion can occur in this case for for , respectively. These overdensities are quite modest (for example, Hunter’s 1977 subsonic infall solution has and ), so we anticipate that photodissociation should not prevent protostellar accretion. The value of the overdensity is likely to vary from one protostar to another, however, so it is possible that in some cases it would be too small to permit accretion. In such cases the infalling gas would decelerate; once it is stationary, however, it could resume accretion when it is overtaken by an expansion wave, as shown by Fatuzzo et al. Our numerical calculations show that the increase in from 1.1 to has only a minor effect on the accretion rate, diminishing it by less than 20%. We conclude that photodissociation of molecular coolants by the protostar does not have a significant effect on its accretion rate.

On the other hand, collapsing cores that do not contain a protostar, but that are close enough to a protostar that their molecular coolants are destroyed, will cease collapsing if their central temperature is low enough ( K) that exceeds . Thus, FUV emission from the first stars is potentially quite effective at suppressing star formation in their vicinity. We can estimate the distance over which star formation is suppressed from the work of Glover & Brand (2001). As in Paper I, we assume that the core is in approximate hydrostatic equilibrium and is characterized by an entropy parameter . We find that the time to dissociate H is less than the free-fall time if the core is within a distance

 D=24[(SLW1049 s−1)(10−3x2)(fabsfdiss0.01)]1/21¯n21/404K′1/4     pc, (18)

of the protostar, where is the photon luminosity in the Lyman-Werner bands, is the fractional abundance of H, is the fraction of the Lyman-Werner flux absorbed by the H, is the fraction of absorptions that result in dissociation, and is the mean density of H nuclei. Thus, even a star, which has  s, can suppress star formation in an existing core only if the core is relatively nearby. A more detailed analysis by Susa (2007) comes to the same conclusion. Ahn & Shapiro (2007) model both dissociation and ionizing feedback and also find a relatively weak suppression of PopIII.2 star formation by neighboring PopIII.1 stars. Whalen et al. (2008) have presented multi-dimensional numerical simulations of these processes.

The second feedback effect of FUV radiation is radiation pressure acting on the Lyman absorption series in the infalling neutral gas. This effect has been considered previously by Oh & Haiman (2002), who studied feedback effects in halos with virial temperatures above K, which are more massive than those we consider. They concluded that Lyman- radiation pressure could be important, but did not find any constraint on the mass of the star that could form. Our work improves upon theirs in several respects: we include stellar continuum photons injected away from line center as well as Lyman- photons emitted in the H ii region; we allow for Rayleigh scattering; we include the limitations on the radiation pressure set by two-photon emission and by the blackbody constraint; and we allow for the effects of rotation in the infalling gas.

Since conditions are very opaque, the Lyman- radiation can be considered to be isotropic. The Lyman- radiation pressure is then

 Pα=13uα=4πJα3c, (19)

where and are the energy density and mean intensity of the Lyman- radiation. The estimation of is complicated by the fact that Lyman- photons can diffuse in frequency as well as in space, and that at the optical depths we are considering, the transfer is dominated by the damping wings of the line profile (Adams 1972). This problem is far too difficult to treat analytically for the geometry and dynamics that we are using to model the protostellar accretion. We therefore make the following substantial approximations when evaluating at the outer boundary of the H ii region, (which may be at the surface of the protostar), and at a particular polar angle. (1) The axially symmetric geometry can be replaced by an equivalent slab geometry. The effects of spherical divergence are incorporated by normalizing the mean intensity to the flux at . The slab column is set equal to 20% of that in the infalling gas based on the discussion in Appendix C. (2) The anisotropy in the optical depth can be accounted for by taking the harmonic mean of the opacity,

 1¯τeff=1A∫dAτ(r) (20)

(see Appendix D). In practice, the escape of photons is primarily controlled by the minimum optical depth, which is in the polar direction, so in our numerical models we evaluate with a column density that is 20% of the column in the vertical direction from the point of interest. (3) Finally, we assume that the effect of the velocity field can be approximated by a Doppler line profile of suitable width (see below).

The propagation of resonance line photons in very opaque media has been treated by a number of authors (Adams 1972; Harrington 1973; Hummer & Kunasz 1980; Neufeld 1990). Let the mean optical depth in the line be

 ¯τ=1ΔνD∫τνdΔν. (21)

In terms of the normalized frequency , we have , where is the line profile. In the Doppler core, the line profile is ; in the damping wings it is , where is the ratio of the natural line width to the Doppler width. In applications, we generally have , and in that case the optical depth at line center, , is related to mean optical depth by . The opacity is , and the mean free path is .

As shown by Adams (1972), resonance photons escape in a single longest excursion from line center. After scatterings, the escaping photon has a frequency shift and it has traveled a distance , where , etc. In order for the photon to escape, this distance must be the size of the region, . This implies and , which in turn gives the characteristic frequency of the escaping photons as . The total path length traversed by the escaping photons is about . As a result, we expect the mean intensity to exceed the incident intensity by a factor of about .

The velocity field has contributions from thermal motions, turbulent motions, and the overall flow. Thermal and microturbulent velocities are naturally included in . In the cases we shall consider, the overall flow is highly opaque, so it generally does not contribute to the random walk of the photons. If the velocity width of the flow (including macroturbulence) is small compared to the line width of the escaping photons, (from Neufeld 1990), then the flow has negligible effect on the escape of the photons. On the other hand, if , then the effective column density of the gas is reduced. For example, in the simple case in which the velocity varies linearly with the column density, photons will interact with only a fraction of the gas. If is the line width in the absence of any flow velocity, then one can show that the effective column density is reduced by a factor of about from the total value. In our numerical models we always set , the value appropriate for gas, the assumed temperature of the infalling neutral gas near the protostar. We set equal to half the difference in radial velocities of the inner and outer edges of the slab. If , which is not usually the case, we reduce the effective column by the factor .

Appendix B describes the general enhancement in the intensity of photons that are trapped by the Lyman- damping wings, and, if the columns are very large, by the opacity due to Rayleigh scattering. Thus, photons from the protostellar continuum, outside the frequency interval defined by , can contribute to the radiation pressure. The enhancement in intensity leads to an increase in the radiation pressure so that the momentum transferred from the radiation to the gas is in each optical depth. As shown in Appendix B, the isotropic component of the radiation pressure is

where and is the intensity of a blackbody with a temperature equal to that at the protostellar surface. When this limit is evaluated for the reprocessed Lyman- photons, the intensity is limited to that of a blackbody at the temperature of the ionized gas in the H ii region (see Appendix B.3). This expression is valid provided , corresponding to  cm for Lyman .

What is the condition for the radiation pressure to halt the accretion? We assume that the accreting gas is inside the sonic point, so that the gas pressure is negligible. For steady, radial flow, the equation of motion of the gas is

since we have shown in Appendix C that the radiative force can be represented by the gradient of the isotropic component of the radiation pressure. We assume that the radiation pressure builds up rapidly over a distance small compared to the radius ; this is justified below. Then constancy of the mass flux implies  const. If the radiative force is to stop the flow in a small distance, then the gravitational force must be negligible in comparison. We then have

so that  const in the deceleration region. When gas enters this region, the radiation pressure is small and , the free-fall velocity; as the gas decelerates, the radiation pressure rises and drops. The inflow will be halted if the radiation pressure at the inner edge of the deceleration region is , where is the density in the freely falling gas.2 We have verified this simple argument by solving for the flow in the case that the flux varies as , with ; such a faster than spherical fall-off in the flux is expected when the density distribution is not spherically symmetric, so that flux will escape into regions of lower opacity (as in the case of the “flashlight effect”—Yorke and Bodenheimer 1999). We find that the radiation pressure required to halt the infall is within about 10% of for In order to reverse the inflow and eject the matter, the radiation pressure must be twice this, . Of course, a steady radial flow cannot reverse direction. To see that the factor 2 is required, one can imagine that the flow is inward over half the sky and outward over the other half; to maintain the same accretion rate, the density would have to be twice as large. If the gas is initially in a disk, there is no infall to start with and the pressure required for ejection is .

We evaluate this criterion along the polar axis at the edge of the H ii region, which is where the ram pressure of infalling gas is a minimum and where breakout should occur first (Figure 3). At low values of the breakout does not occur until the star has reached several hundred solar masses as the photosphere is very large and cool, producing little FUV flux. At reasonable values of , polar breakout can occur relatively early, at . This is the point in the protostellar evolution when the star is starting its rapid contraction to the main sequence, and the surface temperature and luminosity are thus rising. For these values of the densities and ram pressures become significantly greater as the sight line moves away from the pole. By the time that the radiative flux from the star is large enough to reverse the flow in these directions, a polar cavity would have been blown out, thus reducing the enhancement in the radiation pressure due to trapping of photons. Thus although Lyman- radiation pressure can act to reduce the efficiency of accretion, we expect it to be unable to stop it. Even the reduction of the accretion efficiency is likely to be relatively modest, since even a small polar cavity can dramatically reduce the radiation pressure in the H ii region. In the following sections we consider other feedback mechanisms that are more effective at limiting accretion, although they do so at higher masses. When following the stellar evolution to these masses we shall assume the reduction in accretion efficiency due to Lyman- feedback is negligible.

## 5 Ionizing Feedback and Breakout of the H ii Region

### 5.1 Photoionization Heating

Extreme ultraviolet ( eV) radiation from the protostar can ionize infalling neutral gas, creating an H ii region. The temperature of the ionized gas is about  K, based on the models of Giroux & Shapiro (1996) and Shapiro, Iliev and Raga (2004) with stellar spectra. The thermal pressure of the ionized region is typically much greater than that in neutral gas of the same density because of the elevated temperatures and sound speeds: . The pressure gradients created at this ionized-neutral boundary can become steep enough to cause the H ii region to expand and to dramatically reduce the accretion of gas to the star. In this section we attempt to calculate the point in the protostellar evolution at which this transition occurs. This problem has been considered previously by Omukai & Inutsuka (2002). The new feature in our treatment is that we allow for the rotation of the infalling gas, which can significantly reduce the density near the protostar. As we shall see, this completely changes the evolution of the H ii region.

As in Paper I, we approximate the density distribution of the infalling gas by the Ulrich (1976) solution. The gas is assumed to be spherically symmetric far from the protostar, and each mass element conserves its angular momentum as it falls ballistically toward the star. Terebey, Shu, & Cassen (1984) showed that this solution matches on to an expansion-wave solution for the gravitational collapse of a singular isothermal sphere. The resulting density distribution is

 ρ=˙m∗dψ(μ,r)4πr3/2(2Gm)1/2, (25)

where ,

 ψ(μ,r)=⎛⎜ ⎜ ⎜⎝21+μμ0⎞⎟ ⎟ ⎟⎠1/21μμ0+2μ20(rdr), (26)

and is the value of far from the protostar. The two angles are related by

 rdr=μ0−μμ0(1−μ20), (27)

which shows that : the gas converges toward the disk plane. Ulrich assumed that the disk had negligible mass, so that in equation (25). In our case, varies smoothly from to as shrinks from being much greater than to being much less than . This variation in the mass acting on the infalling gas leads to small, unknown deviations from the Ulrich solution. In view of the necessarily approximate nature of the solution and the fact that and differ by only a factor in the fiducial case, we shall set in applying equation (25).

The density factor depends on both the current direction cosine, , as well as the initial one, , with the two being related by the cubic equation (27). In our analysis, it is convenient to have an approximation for in which the dependence on is eliminated:

 ψ≃[21+Max(μ2/3,1−ζ)]1/210.5(ζ−1+3|ζ−1|)+3μ2/3Min(1,ζ), (28)

where . This is exact for all at , where , and at , in the plane of the disk. It is also exact at for all . For it is accurate to better than 20%; for , it is accurate to within a factor 2. To simplify our results, we approximate this further for and take

 ψ∼(21+μ)1/2r2rd,     (r≲0.5rd). (29)

This approximation is quite accurate at (better than 20% for ), but it deteriorates away from there, underestimating the density by a factor in the plane for (the accuracy improves as decreases). Nonetheless, it suffices to give an analytic estimate for the behavior of the H ii region.

#### Early Evolution of the H ii Region

H ii regions are bounded by ionization fronts. Ionization fronts that move faster than about with respect to the neutral gas, where is the isothermal sound speed of the ionized gas, are termed “R-type;” such fronts have little dynamical effect (Spitzer 1978). However, if the velocity of the front slows below , a shock forms in front of the ionization front and the velocity of the front into the shocked gas falls to , where is the isothermal sound speed of the shocked neutral gas; such ionization fronts are termed “D-type.” When the H ii region first forms, it is embedded in gas falling inward with a velocity . As a result, the ionization front is initially R-type, and the radius of the H ii region, , is determined by ionization balance.

Since the density of the infalling gas depends on the angle relative to the axis of rotation, depends on angle also. We determine this radius in the sector approximation, in which ionizations balance recombinations in each element of solid angle:

 dSdΩ=S4π=∫rHIIr∗r2α(2)nenpdr, (30)

where is the rate of emission of ionizing photons and  cm s is the recombination rate to the excited states of ionized hydrogen. In writing equation (30), we have made three approximations. First, we have assumed that the rate of emission of ionizing photons is much greater than the rate of accretion of hydrogen atoms so that ionizations and recombinations are very nearly in balance (note that advection of neutral H into the H ii region is allowed for in the numerical models). For a mass accretion rate of yr, the hydrogen accretion rate is about  s. The mass at which the ionizing photon luminosity exceeds this value depends on ; for the fiducial case of , this occurs at about . Second, we have assumed that in the outer parts of the H ii region, where the helium is singly ionized for stellar temperatures  K, each recombination of He results in one H ionization, whereas it actually results in about 2/3 of an ionization at the relevant densities ( cm; Osterbrock 1989). In fact, the abundance of He is sufficiently small () that we shall henceforth neglect it in our analytic estimates (however, we do not neglect its contribution to the mass density, nor is it neglected in the numerical calculations). Finally, we have ignored photoionization from the level of H, so that our calculation somewhat underestimates the size of the ionized region, although this is not very important at the densities resulting from realistic values of .

With the density distribution given by equation (25), equation (30) becomes

 S=α(2)˙m2∗dI8πμ2HGm∗d≡ScrI, (31)

where  g is the mass per hydrogen and

 I≡∫rHIIr∗ψ2(μ,r)r dr. (32)

We have set in equation (25) in accord with the discussion following equation (27). Equation (31) reduces to that of Omukai & Inutsuka (2002) for (and if is replaced by , by and by ). As shown by Omukai & Inutsuka, an H ii region in an density profile is confined to the vicinity of the star for and expands to exponentially large distances as increases above . Numerically, we have

 Scr=3.07×1050(2.5T4)0.8˙m2∗d,−3m∗d,2    ph s−1→2.36×1051(2.5T4)0.8K′30/7m13/7∗,2    ph s−1. (33)

By comparison, the ionizing luminosity of a Pop III star is

 S≃7.9×1049ϕSm1.5∗,2    ph s−1, (34)

which for is a fit to Schaerer’s (2002) results for the ionizing luminosity of main sequence primordial stars; the fit is accurate to within about 5% for . As shown in Paper I, the ionizing luminosity can be less than the main sequence value () when the star is still contracting toward the main sequence, and it can be greater when it is accreting while on the main sequence; for the case illustrated in Paper I, . If the accretion rate is not reduced by feedback effects, would not exceed until for . However, as we shall see, rotation makes the factor small and allows the H ii region to expand until it is almost as large as the disk even when the mass is less than this.

At early times we have so that we can use approximation (29) for the density. As a result, we find

 SScr=I≃14(1+μ)(rHIIrd)2. (35)

With equations (33) and (34), we then obtain

 rHIIrd = 1.01(1+μ)1/2(1+fd)1/2ϕ1/2S(T42.5)0.4m1.25∗,2˙m∗d,−3, (36) → 0.37(1+μ)1/2ϕ1/2S(T42.5)0.4m47/28∗,2K′15/7; (37)

recall that “” indicates that we have taken . We see that for we have , so that our approximation for the density, equation (29), is reasonably accurate at early times. The radius of the H ii region is then

 rHII=5.40×1015(1+μ)1/2ϕ1/2S(T42.5)0.4(fKep0.5)2(1+fd)31/14m3∗,2K′25/7   cm      (m∗,2≲1), (38)

where we have set the star formation efficiencies and equal to unity and where we have made the approximation in the exponent of . As the stellar mass increases above , the approximation for the density, equation (29), increasingly overestimates the density except near the equator; as a result the radius of the H ii region, , becomes larger than the value given in equation (38) except near the plane of the disk, where the high density traps the H ii region. As remarked above, for , which occurs for if the accretion continues unabated by feedback, increases exponentially with (Omukai & Inutsuka 2002).

#### Later Evolution of the H ii Region: Suppression of Accretion

According to equation (38), the radius of the H ii region expands on the protostellar evolution timescale yr, which is far longer than the dynamical time  yr. As a result, the velocity of the ionization front relative to the infalling gas is very nearly equal to the free-fall velocity. The first phase of evolution of the H ii region ends when it expands to the point that the radius becomes comparable to the gravitational radius,

 rg≡GϕEddm∗dc2i=3.92×1015ϕEdd(2.5T4)m∗d,2   cm, (39)

where we have taken the gravitating mass to be . Here we have allowed for the decrease in the radiation pressure due to electron scattering through the factor

 ϕEdd≡1−LLEdd, (40)

where is the Eddington limit. In Paper I we found that for , which corresponds to . Equations (38) and (39) relate the protostellar mass to ,

 m∗,2=0.85⎡⎣ϕ1/2EddK′25/14(1+μ)1/4(1+fd)17/28ϕ1/4S⎤⎦(2.5T4)0.7(0.5fKep)(rHIIrg)1/2. (41)

Keep in mind that this relation is valid only for , so that cannot be small. This condition is satisfied insofar as the simulations of ABN are representative of the angular momentum of the accreting gas, since they give .

When the H ii region expands to the point that , a shock front forms and the ionization front becomes D-type; this occurs at . The accretion rate through the H ii region will begin to decrease at this point. Since the shocked neutral gas is denser than the ionized gas in the H ii region, the accelerating expansion of the shocked shell is subject to the Rayleigh-Taylor instability, and as a result it is difficult to estimate by how much the accretion is reduced. While the shell is moving slowly, it can fall onto the disk and accrete that way. However, once the shocked shell is moving faster than the local free-fall velocity it seems unlikely than any significant further accretion can occur. To obtain an approximate upper limit on the accretion through the H ii region, we assume that, from a given direction, the accretion is unimpeded until the H ii region has expanded to a radius equal to . Because of the declining density distribution in the infall envelope, the H ii region is expanding relatively rapidly at this stage and so soon ionizes a large region beyond , which we expect leads to a substantial reduction in accretion rate from the affected directions. This approximation needs to be checked with multi-dimensional radiation-hydrodynamical simulations. It is important to bear in mind that this suppression of the accretion occurs only in the ionized gas. Gas in the shadow of the accretion disk around the star can continue to accrete, as discussed below.

In figure 4 we show the geometry of the H ii region near the point of polar breakout of the ionized gas beyond . In this calculation the protostellar evolution has been followed as described in Paper I, but now including the effects of a reduction in accretion rate once the H ii region breaks out beyond . This has only just started to occur at the point of the evolution shown in the figure. We have assumed there is negligible reduction in the accretion rate because of the Lyman- radiation feedback since we expect its effects to be limited to relatively small angles around the polar axis. The extent of the H ii region is calculated using the sector approximation using the density distribution model of Ulrich (1976) as described above. We include the effect of electron scattering, but not force due to photoionization, which is discussed below. The effect of radiation pressure due to photoionization is strongest for purely radial infall, so its neglect is not critical for the models presented here. We do allow for advection of neutrals into the H ii region, though they are not very important by the time the protostar is . A temperature of was adopted for the ionized gas, which affects the value of . Note that in figure 4 we have assumed an infinitely thin accretion disk. The polar and equatorial breakout conditions are shown as a function of in figure 5. Once the protostar has reached the masses indicated by the “Equatorial” line in this figure, we do not expect accretion to be able to proceed from directions that have a direct line of sight to the protostar, i.e. those directions that are not shielded by the accretion disk. Thus in most cases we do not expect H ii region breakout to set the final mass of the star, but rather to cause a decrease in accretion efficiency that starts to become important at about in the fiducial case. The actual reduction in accretion efficiency depends on the thickness of the accretion disk, to be discussed below (§6).

We can compare the analytic prediction for H ii region breakout (eq. 41) with our numerical calculation, which for the fiducial model (, , , ) finds breakout in the polar direction at a mass of 45.3 . At this point the total H-ionizing luminosity is so that and the bolometric luminosity is so that . With these values, the analytic estimate for the mass at which polar () H ii region breakout () occurs (eq. 41) is 44.7 , in excellent agreement with the numerical value. In Fig. 5 we see that the mass at which the H ii region breaks out does not scale as a simple power of . This is because and vary with stellar mass, especially for .

### 5.2 Radiation Pressure due to Photoionization

Continuum radiation is dynamically coupled to the gas in the H ii region, both through Thomson scattering and through photoionization. Since the H ii region is optically thin to Thomson scattering, it effectively reduces the force of gravity by a factor , which as discussed above is for . At distances large enough that the mass acting on the gas is , is reduced by a factor . Keep in mind that the decrease in the effective gravity due to electron scattering reduces the accretion rate of ionized gas only; it does not affect the accretion of neutral gas outside the H ii region onto the disk.

Every photoionization results in a transfer of momentum to the gas, where is the mean energy of the photons that ionize the gas. The importance of radiation pressure associated with photoionization has long been appreciated in studies of active galactic nuclei and X-ray sources (e.g., Tarter & McKee 1973); Haehnelt (1995) pointed out its importance in the formation of the first galaxies, and Omukai & Inutsuka (2002) discussed its role in the formation of the first stars. They showed that, under the assumption of perfect spherical symmetry, radiation pressure would have the counter-intuitive effect of reducing the size of the H ii region, thereby eliminating any feedback effect on the growth of the protostar. Since photoionizations are balanced by recombinations, the radiative force is given by . Generalizing their treatment to include electron scattering, we find that this radiative force balances the effective gravity at a critical density given by

 α(2)n2cr(hνic)=ϕEddρcrGm∗dr2, (42)

so that

 ncr=2.15×106ϕEdd(T42.5)0.8m∗d,2r2