WEAKLY INTERACTING MASSIVE PARTICLE DARK MATTER AND FIRST STARS: SUPPRESSION OF FRAGMENTATION IN PRIMORDIAL STAR FORMATION

We perform the calculations for this paper using the smoothed particle hydrodynamics (SPH) code gadget 2 (Springel 2005). We have substantially modified this code to include a fully time-dependent chemical network, details of which can be found in the Appendix of Clark et al. (2011a). Our treatment includes: H₂ line cooling (Glover & Abel 2008), which is treated in the optically thick regime using the Sobolev approximation (Yoshida et al. 2006; Clark et al. 2011a); collision-induced emission from very high density H₂ (Ripamonti & Abel 2004); cooling due to the collisional ionization and recombination of hydrogen and helium (Glover & Mac Low 2007); and heating and cooling arising due to changes in the molecular fraction, shock heating, compressional heating, and cooling due to rarefactions. Turk et al. (2011) showed that the choice of the H₂ three-body formation rate coefficient influences the resulting dynamics of the gas within the halo. In this work we use the three-body H₂ formation rate of Glover (2008) which is intermediate within the range of the published rates.

To treat gas which collapses to scales smaller than the resolution limit of our simulations, we use a sink particle approach. We create sinks using the standard Bate et al. (1995) algorithm. Our particular implementation of this is the same as that used in Jappsen et al. (2005). The minimum number density for sink particle creation is 10¹⁶ cm⁻³, but it is important to note that candidate regions must also satisfy a series of tests to ensure that they are unambiguously collapsing. The ratio of thermal energy and rotational energy to the magnitude of the gravitational potential energy (α and β, respectively) for the particles that will be turned into the sink must satisfy the conditions

$$\begin{equation} \alpha \le \frac{1}{2}, \end{equation} \tag{ 1 }$$

$$\begin{equation} \alpha + \beta \le 1. \end{equation} \tag{ 2 }$$

In addition, the total energy of the particles must be negative, and the divergence of their accelerations must be zero in order to prevent structures that are likely to bounce or be tidally disrupted being turned into sinks. Consequently sink particle formation often occurs above this density threshold.

We set the outer accretion radius of the sink particles, r_out, to 6 AU and the inner accretion radius, r_in, to 4 AU. SPH particles that come within a distance r_out > r > r_in of a sink are accreted only if they are gravitationally bound to that sink. SPH particles that come closer than r_in are always accreted. Once sinks have formed, we account for the accretion luminosity produced by the ongoing accretion of gas onto the sinks using the scheme described in Smith et al. (2011), which assumes that all of the protostars that form are normal Population III protostars, and not dark stars. We would expect the heating produced by accretion onto dark stars to be smaller, owing to their larger radii, and so this procedure gives us an upper limit on the effect of the accretion luminosity.

Most current simulations have insufficient resolution to follow the DM contraction within a primordial halo with the same resolution as that used for the baryons. For example, Abel et al. (2002) had a DM mass resolution of 1.1 M_☉ which meant they could not determine the DM profile at radii smaller than around 0.05 pc. Instead, the method of adiabatic contraction developed by Blumenthal et al. (1986) has been used by various authors (e.g., Spolyar et al. 2008; Ripamonti et al. 2010) to calculate how the DM distribution will respond to changes in the gravitational potential of the gas. Spolyar et al. (2008) applied this method to the collapsing halo simulations of Abel et al. (2002) and Gao et al. (2007) and found that the DM density at the outer edge of the baryonic core after adiabatic contraction was well fit by the expression

$$\begin{equation} \rho _{xc} \approx 5 n_{\rm p}^{0.81} \: {\rm GeV \: cm^{-3}}, \end{equation} \tag{ 3 }$$

where ρ_xc is the DM density at the edge of the core and n_p is the number density of hydrogen nuclei in the core.⁵

Given the considerable technical challenges involved if one attempts to follow the evolution of the DM density profile directly, we do not attempt to do so. Instead, we parameterize its effects with the help of the adiabatic contraction results of Spolyar et al. (2008). The DM density profile is calculated within the code as follows. At each time step, the point of maximum gravitational potential energy in the halo is found, which we assume to be the point at which the DM density has its maximum value. Next, the maximum baryon density of the halo is found and the DM density at the outer edge of the core calculated from it using Equation (3). Its radial dependence is then described by a two-part power-law fit, using power-law slopes drawn from the calculations of Ripamonti et al. (2010). In the outer halo, at distances r > r_c from the center of the density profile, we have

$$\begin{equation} \rho _{x} = \rho _0 \left(\frac{r}{1 \,\mathrm{pc}} \right)^{-1.8} \: {\rm GeV \: cm^{-3}}, \end{equation} \tag{ 4 }$$

where r is the radial distance from the center and ρ₀ = 5 × 10⁴ GeV cm⁻³ is a normalizing constant that is equal to the DM density at r = 1 pc. The size of the central DM core, r_c, is determined by finding the radius at which the DM densities given by Equations (3) and (4) are equal, which yields

$$\begin{equation} r_{c} = 16.7 \left(\frac{n_{\rm p}}{10^{14} \: {\rm cm^{-3}}} \right)^{-0.81/1.8} \: {\rm AU}. \end{equation} \tag{ 5 }$$

Finally, the DM density profile within the core, at radii r < r_c, is given by

$$\begin{equation} \rho _{x} = \rho _{xc} \left(\frac{r}{r_c} \right)^{-0.5}. \end{equation} \tag{ 6 }$$

It is unclear how efficient adiabatic contraction will be within the central regions of the DM halo, given that gravitational collapse of the gas in these regions occurs rapidly. Our adoption of a peaked DM profile within the core may therefore be an overestimate; see, for instance, the discussion of the validity of the adiabatic contraction approximation given in Ripamonti et al. (2010).

Given the DM density profile, the contribution of DMA to the halo energy budget is computed as follows. We follow Spolyar et al. (2008) and adopt the standard DMA cross section of 〈σv〉 = 3 × 10⁻²⁶ cm³ s⁻¹. For the DM particle mass, m_x, we note that astrophysical constraints require that m_x ⩾ 10 GeV (see, e.g., Ackermann et al. 2011; Galli et al. 2011). Additionally the desire to avoid fine-tuning while still explaining the gauge hierarchy problem argues for an upper limit m_x ∼ a few TeV. We therefore adopt a DM particle mass of m_x = 100 GeV for our fiducial case, close to the middle of this range of values, but we also explore the effect of varying m_x. Following Valdés & Ferrara (2008), we assume that one third of the energy from DMA is immediately carried away by neutrinos and that the remaining energy is split between direct heating of the gas, ionization and dissociation of the constituents of the gas, and the excitation of H, He, and H₂. The fraction that is deposited as heat is given by

$$\begin{equation} f_h=1-0.875\big(1-x_{e}^{0.4052}\big), \end{equation} \tag{ 7 }$$

and the fraction that goes into ionization is given by

$$\begin{equation} f_i=0.384\big(1-x_{e}^{0.542}\big)^{1.1952}, \end{equation} \tag{ 8 }$$

where x_e ≡ n_e/n is the fractional abundance of electrons. The remaining energy is radiated away in electronic transitions from H, He, and H₂ and plays no direct role in the evolution of the gas. The total power per unit volume injected by DMA is

$$\begin{equation} W_{\rm ann}=2 m_x \dot{N}_{\rm ann}, \end{equation} \tag{ 9 }$$

where $\dot{N}_{\rm ann}$ is the number of DMAs per unit volume per unit time. This is given by

$$\begin{equation} \dot{N}_{\rm ann}=\frac{1}{2} n_x^2 \langle \sigma v\rangle , \end{equation} \tag{ 10 }$$

where n_x = ρ_x/m_x is the number density of DM particles. The resulting energy injection rates per unit volume for heating (Q_h) and ionization (Q_i) are

$$\begin{equation} Q_h=2 f_a \dot{N}_{\rm ann} (1-e^{-\tau _x})f_h m_x, \end{equation} \tag{ 11 }$$

$$\begin{equation} Q_i=2 f_a \dot{N}_{\rm ann} (1-e^{-\tau _x})f_i m_x, \end{equation} \tag{ 12 }$$

where f_a = 2/3 is the fraction of the total energy that affects the gas (one third of the energy escapes in the form of neutrinos).

We estimate the optical depth of the gas to the annihilation products, τ_x, by assuming that the baryonic density profile is flat within the DM core and declines as a power-law ρ∝r⁻² at larger radii. For a profile of this form, the column density of gas, measured radially outward from a point r, is given at r ⩽ r_c by

$$\begin{equation} \Sigma (r) = \rho (r) \left[(r_{c} - r) + \left(r_{c} - \frac{r_{c}^2}{r_{h}}\right) \right] \end{equation} \tag{ 13 }$$

and at r > r_c by

$$\begin{equation} \Sigma (r) = \rho (r) \left(r - \frac{r^2}{r_{h}}\right), \end{equation} \tag{ 14 }$$

where r_h is the virial radius of the halo, which we take to mark the edge of the gas distribution. If we follow Ripamonti et al. (2010) and adopt a constant gas opacity of κ = 0.01 cm² g⁻¹, and also restrict our attention to scales small compared to r_h, then we can write the optical depth as

$$\begin{equation} \tau _x \equiv \kappa \Sigma (r) = \left\lbrace \begin{array}{@{}lr}\kappa \rho (r) (2r_{c} - r) & r \le r_{c} \\ \ms \kappa \rho (r) r & r > r_{c} \end{array} \right.. \end{equation} \tag{ 15 }$$

In practice, this expression for τ_x is an underestimate, for a couple of reasons. First, it assumes that the optical depth in any direction from point r is the same as the value that we measure along a ray passing radially outward, when in reality the optical depths in other directions will be somewhat larger. Second, the baryonic density profiles that we recover in our simulations do not completely match up with the profile we assumed to derive τ_x. Our assumed profiles and the real profiles both have power-law slopes ρ∝r⁻² in their outer regions, but the real profiles do not have the completely flat core that we assume within r_c. Our estimate of τ_x is therefore too small at distances r ≪ r_c, meaning that our derived heating and ionization rates may also be too small. In practice, we do not expect this to be a major source of error, as in the regime where DMA becomes important, gas with r ≪ r_c will typically have τ_x ≫ 1, meaning that the heating and ionization rates are insensitive to the precise value of τ_x.

To convert from Q_i, the total energy per unit time per unit volume that goes into ionizations, to the ionization rates of the individual chemical species present in the gas, we follow the procedure outlined in Ripamonti et al. (2007). They identified seven main reactions that can be caused by DMA: the ionization of H, He, He⁺, and D, and the dissociation of H₂, HD, and H⁺₂. Between them, these seven reactions account for the bulk of the energy lost in the form of ionizations or dissociations, and we follow Ripamonti et al. (2007) and split the energy available for ionization between these seven species according to their relative abundances.

In Figure 1, we show the density profiles of gas and DM for simulations H1, H1-lm, and H2, plus the two reference models. The simulations are compared when the hydrogen nuclei number density at the center of the halo first reaches 5 × 10¹⁴ cm⁻³, which corresponds to a DM core radius of r_c ∼ 8 AU and a time of roughly six years before the formation of the first protostar. It is immediately apparent that the gas in each halo is able to collapse to high densities, regardless of the strength of the DMA feedback. This is true even in our maximal feedback model, H1-lm, where the DM particle mass was only 10 GeV. For reference, in Spolyar et al. (2008) it was predicted that collapse would stop at densities of 10¹³ cm⁻³ for a DM mass of 100 GeV, and at densities of 10⁹ cm⁻³ for a DM mass of 10 GeV. We find no evidence for this in our simulations. On the other hand, our results are in good agreement with the one-dimensional results of Ripamonti et al. (2010), who found that the gas could collapse to densities of 10¹³–10¹⁴ cm⁻³ in all of their models. We also found collapse up to high densities in run H1-hm, performed with a DM particle mass of 1000 GeV, but as the effects of DMA in this case are much smaller than in our fiducial case, we do not discuss this run further in this section.

Zoom In Zoom Out Reset image size

The main effect that DMA has on the density profile of the gas is the clear density enhancement at radii of 100–1000 AU. We also find that in the runs with DMA, the time taken for the gas to collapse to protostellar densities is shorter than in the reference runs (see Table 2). Both of these effects can be understood by an analysis of the temperature structure of the gas. In Figure 2, we show how the temperature of the gas varies as a function of the gas number density. In the outer regions of the halo, the gas is cooler than in the reference case. At higher densities, however, there is a sharp rise in the gas temperature, taking it above the value in the reference case. Comparison of Figures 1 and 2 shows that the density at which this temperature increase occurs is the same as the density at which we first see the bump in the density profile.

Zoom In Zoom Out Reset image size

In Figure 3, we show how the fractional abundances⁶ of H₂ and e⁻ vary as a function of density within our different models. At low densities, DMA-induced ionization enhances the number of electrons by two or more orders of magnitude. This increases the rate of H₂ formation relative to the case without DMA by increasing the effectiveness of the following reaction chain:

$$\begin{eqnarray} {\rm H}+e^- & \rightarrow & {\rm H}^- + \gamma , \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} {\rm H}^- + {\rm H} & \rightarrow & {\rm H}_2 +e^{-}. \end{eqnarray} \tag{ 17 }$$

As H₂ is the main coolant of the gas at these densities, the enhanced H₂ abundance leads to more efficient cooling. We therefore find, counterintuitively, that this lower density gas is cooler in the case with DMA than in the case without DMA. This is in agreement with the results of the one-dimensional simulations of Ripamonti et al. (2010) and provides an explanation for the shorter collapse times of the halos in which DMA is occurring, as for the majority of its lifetime, the gas has less thermal support.

Zoom In Zoom Out Reset image size

Once the gas reaches a density of around 10¹² cm⁻³, however, the H₂ fraction decreases sharply in the runs with DMA. By comparing this figure with the temperature distribution shown in Figure 2, we see that this sharp decrease occurs once the gas temperature reaches T ∼ 2000 K. It occurs because at this temperature, collisional dissociation of H₂ by the reactions

$$\begin{eqnarray} {\rm H_{2}} + {\rm H} & \rightarrow & {\rm H} + {\rm H} + {\rm H}, \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray} {\rm H_{2}} + {\rm H_{2}} & \rightarrow & {\rm H} + {\rm H} + {\rm H_{2}}, \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} {\rm H_{2}} + {\rm He} & \rightarrow & {\rm H} + {\rm H} + {\rm He} \end{eqnarray} \tag{ 20 }$$

and by the charge transfer reaction

$$\begin{equation} {\rm H_{2}} + {\rm H^{+}} \rightarrow {\rm H_{2}^{+}} + {\rm H} \end{equation} \tag{ 21 }$$

becomes effective. The increased temperature leads to an increased thermal pressure, which slows the collapse of the gas at these densities and leads to a "pile-up" of material visible as a pronounced bump in the gas density profile. It is also clear from Figure 3 that not all of the H₂ is destroyed at this density. The H₂ that survives in the higher density gas enables it to dissipate much of the energy introduced into the gas by the DMA, and hence allows it to maintain its temperature at close to 2000 K. This allows gravitational collapse to continue past the point at which DMA heating first outstrips H₂ cooling. We examine the physics of this in more detail in the next section.

Finally, we note that one must take care in interpreting the behavior of the H₂ abundances shown in Figure 3. It is all too easy to think of the evolution of $x_{\rm H_{2}}$ with density as a time history of the H₂ abundance. In this picture, the H₂ abundance first increases as the gas collapses, then sharply decreases once the number density reaches 10¹² cm⁻³, before increasing once more at higher gas densities. However, this is incorrect, as Figure 4 demonstrates. The rise in the H₂ abundance in the densest gas that we see if we look at a single snapshot is not due to the reformation of H₂ in this gas. Instead, it occurs because the amount of H₂ that the gas loses once it reaches a density of around 10¹² cm⁻³ increases as the collapse proceeds. The first gas to collapse—the material that is now at densities above 10¹⁴ cm⁻³—loses only a relatively small fraction of its H₂, while the material falling in later loses much more of its molecular content. In part, this behavior is due to an increase in the DM number density associated with the n = 10¹² cm⁻³ gas, caused by the increasing concentration of the halo. However, an additional contribution comes from the fact that gas falling in at later times will generally have a higher infall velocity, and hence will shock more strongly once it reaches this density.

Zoom In Zoom Out Reset image size

4.2.1. Fiducial Case

Figure 5 shows the main heating and cooling processes acting in run H1 shortly before the formation of the first protostar, at the time when the central number density first reached 5 × 10¹⁴ cm⁻³. At low densities, compressional heating and H₂ formation heating are the two most important heat sources, although even at densities as low as 10⁶ cm⁻³, DMA heating is beginning to contribute to the total heating rate at the level of a few percent or more. The cooling of the gas at these densities is dominated by H₂ line cooling. Once the gas density reaches n ∼ 10⁸ cm⁻³, compressional heating starts to become unimportant, and DMA heating catches up with H₂ formation heating, with both subsequently playing important roles in the thermal balance of the gas. The heating produced by these two processes raises the gas temperature, but this enables H₂ line cooling to become more effective, and this is able to offset much of the additional heat input from the DMA at densities n < 10¹² cm⁻³. Above this density, H₂ line cooling becomes relatively ineffective, both because the H₂ is starting to dissociate and also because the H₂ emission lines become optically thick at high densities. The heating rate due to DMA therefore outstrips the H₂ line cooling rate, in agreement with the prediction of Spolyar et al. (2008). However, we see from Figure 5 that this is not the whole story. As the temperature rises, H₂ collisional dissociation (reactions 18–20) and the destruction of H₂ by charge transfer with H⁺ (reaction 21) become increasingly effective. These reactions are endothermic, and hence remove thermal energy from the gas. If H₂ subsequently reforms, then the associated H₂ formation heating will return some of this energy to the gas, but if more H₂ is destroyed than can reform, then the overall effect is to dissipate energy. This effect should not be regarded as "cooling" in the same sense as H₂ line cooling or collisionally induced emission cooling, as it will not bring about a decrease in the gas temperature. Instead, it acts more like a thermostat, preventing the temperature from increasing significantly until all of the H₂ has been consumed.

Zoom In Zoom Out Reset image size

The importance of this effect can be appreciated if we compare the time required to destroy all of the H₂ at the center of the halo with the free-fall collapse time of the gas. By considering the energy budget of the gas in this way we can obtain a rough estimate of whether collapse is likely to be halted at densities close to our resolution limit. From Figure 5, we see that at a density n = 5 × 10¹⁴ cm⁻³, the DMA heating rate per unit volume at the center of the halo is roughly 5 × 10⁻⁶ erg s⁻¹ cm⁻³. If we take the H₂ fractional abundance in this gas to be $x_{\rm H_{2}} = 0.4$, which is a conservative estimate, then the total amount of energy per unit volume that can be dissipated by dissociating H₂ is

$$\begin{eqnarray} E_{\rm H_{2}} & = & 4.48\, {\rm eV} \times x_{\rm H_{2}} n, \hspace{92pt} \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} & \simeq & 1400 \left(\frac{n}{5 \times 10^{14} \: {\rm cm^{-3}}}\right) \: {\rm erg} \: {\rm cm^{-3}}, \end{eqnarray} \tag{ 23 }$$

where 4.48 eV is the binding energy of a single H₂ molecule. The time required to destroy all of the H₂ is therefore t_dis = 1400/5 × 10⁻⁶ ∼ 3 × 10⁸ s, or around 10 years. In comparison, the free-fall time of the gas at this density is around two years. The energy input from the DM therefore cannot destroy the H₂ rapidly enough to prevent the gas from collapsing. At even higher densities, the DMA heating rate will be larger. However, the heating rate per unit volume within the core of the DM density profile scales with the central gas density as n^1.62, and hence t_dis∝n^−0.62. In comparison, the free-fall timescale scales as t_ff∝n^−0.5. Therefore, a large increase in density is required in order to significantly alter the ratio of the H₂ dissociation timescale to the free-fall timescale.

In reality collapse of the dense gas is likely to occur on a timescale that could be a factor of a few longer than the free fall timescale. However, even in this case the collapse timescale is shorter than the H₂ dissociation timescale, and will remain so until the density increases significantly. It therefore seems unlikely that collapse will halt at a density just above our sink creation threshold, although we cannot say at exactly what point the collapse will halt without performing much higher resolution simulations.

Figure 6 shows the main heating and cooling processes acting in the reference case, run H1-ref. In the absence of DMA heating, the main heating term is compressional heating over most of the range of densities examined, with H₂ formation heating becoming important at densities between 10¹⁰ cm⁻³ and 10¹² cm⁻³, and at n > 10¹⁴ cm⁻³. The main source of cooling at n < 10¹⁴ cm⁻³ is H₂ line cooling, while at higher densities, H₂ collisional dissociation plays an important role in regulating the temperature of the gas. This is in reasonable agreement with the results of other models; for instance, Yoshida et al. (2006) find that H₂ dissociation becomes significant at densities of ∼10¹⁵ cm⁻³, at which point the gas temperature is roughly 2000 K. The main effect of the DMA heating seems to be simply to bring the gas to this state at an earlier point in its evolution.

Zoom In Zoom Out Reset image size

4.2.2. Maximal Case

The previous analysis is for the fiducial case where the DM particle mass was 100 GeV. However, even in our maximal case, where the DM particle mass was 10 GeV, we find broadly similar behavior. In Figure 7, we show the rates of the main heating and cooling processes in run H1-lm shortly before the formation of the first protostar. In this case, we see that DMA heating becomes important at an earlier time and that H₂ dissociation becomes the dominant energy dissipation mechanism once the gas density reaches a lower density, n ∼ 10¹⁰ cm⁻³, than in the fiducial model. This is consistent with the results plotted in Figure 2, which show that the gas reaches T ∼ 2000 K slightly earlier in its evolution. However, the subsequent behavior of the gas is very similar in both models. H₂ dissociation is such an effective thermostat that an increase in the heating rate by a factor of 10 produces only a small increase in the gas temperature. Repeating our previous analysis of the H₂ dissociation timescale, we find that in this case t_dis ∼ t_ff, meaning that the gas will probably lose most of its H₂ before it collapses to protostellar densities. It is therefore possible that in this extreme case, a "dark star" will form, but if so, then its size will be smaller than our minimum spatial resolution of a few AU. For comparison, Spolyar et al. (2008) predict that even in the 100 GeV case, a dark star of size ∼20 AU should form, and expect that reducing the DM particle mass will lead to an even larger dark star. We find no evidence for such large DMA-supported structures in our simulations.

Zoom In Zoom Out Reset image size

In their original study, Spolyar et al. (2008) made the simple assumption that the gravitational collapse of the gas would come to a halt as soon as the heating rate produced by DMA exceeded the H₂ cooling rate. As a result, they predicted the formation of large (∼20 AU or more) protostellar objects supported by DMA heating, the so-called dark stars. However, their assumption is incorrect. Ripamonti et al. (2010) first showed, using one-dimensional spherically symmetric models, that the collapse does not halt once the DM heating rate exceeds the H₂ cooling rate, and in our present study we confirm their results. The key factor that allows collapse to continue is the collisional dissociation of H₂. This acts as a thermostat, preventing the gas temperature from increasing significantly above 2000 K until all of the H₂ is destroyed. A simple estimate of the timescale on which the H₂ is destroyed shows that even in the most extreme case, it is comparable to the free-fall time of the gas, meaning that the gas can reach much higher densities (and hence smaller length scales) than anticipated in the Spolyar et al. (2008) study. As a result, we find no evidence for the formation of dark stars within our simulations, although we note that we cannot rule out the formation of such objects if they form above densities of 10¹⁶ cm⁻³, have initial masses of less than 0.01 M_☉ and radii of less than an AU, as this is below our resolution limit.

Perhaps the most striking finding of our work is the large reduction of the level of fragmentation within the protostellar accretion disk. Of our four simulations with DM, secondary fragmentation occurred in only two cases and then only at large radii. Much of the recent focus of studies of Population III star formation has been on the fragmentation of protostellar disks (see, e.g., Clark et al. 2011b; Greif et al. 2011; Smith et al. 2011). Greif et al. (2012) have shown that fragmentation occurs on scales of 10 AU or smaller and that mergers between protostars should be common. The possibility of mergers and interactions between protostars was also highlighted by Smith et al. (2012), who show that the high accretion rates experienced by the young protostars lead to large, extended, "fluffy" objects that have a higher probability of interacting than more conventional protostars.

Such scenarios are far less likely when the effects of DMA are taken into account. In this picture either a single massive star forms at the center of the halo, or a secondary forms in a wide binary with a separation of around 1000 AU. Fragmentation on such scales was also seen by Stacy et al. (2010) in lower resolution studies which focused on scales larger than the inner disk around the protostar. In our model, when secondary fragmentation begins the mass of the primary protostar is already greater than 10 solar masses. McKee & Tan (2008) find that for an accretion rate of 10⁻³ M_☉ yr⁻¹ an ionized H ii region will form around the protostar after it obtains a mass of around 20 M_☉, depending upon its rotation. Such an ionized region is likely to form even earlier in our model due to the additional ionization of hydrogen from DMA. Thus, it is possible that at the point when secondary stars form the original protostar will already be surrounded by an H ii region and any further growth of the primary will be determined by the balance of radiative transfer effects. Hosokawa et al. (2012) find that UV radiation from primordial protostars would photo-evaporate their accretion disks once they achieved a mass of around 43 M_☉. We note however that similar claims for present-day high-mass star formation (Yorke & Sonnhalter 2002) have not borne out in detailed three-dimensional calculations (Krumholz et al. 2009; Peters et al. 2010; Kuiper et al. 2012).

The biggest assumption in our model is the adoption of an analytic DM halo rather than a live halo. No calculation has yet been able to fully follow the contraction and fragmentation of baryonic matter down to AU scales with a comparable resolution in DM. Simulations by Abel et al. (2002) followed the collapse of a minihalo with a mass resolution of 1.1 M_☉ for the DM component and confirmed that the DM had a peaked profile to radii as small as ∼1, 500 AU. Previous work by Ripamonti et al. (2010) and Spolyar et al. (2008) adopts the method proposed by Blumenthal et al. (1986) to model the contraction of a DM halo due to baryonic collapse, and upon these results we base our DM profile. In the Appendix of Ripamonti et al. (2010), the applicability of this method is discussed and found to be in better agreement with the current numerical findings than alternative models Steigman et al. (e.g., 1978). We are therefore confident that our model for increasing the DM density with increasing baryon density is reasonable. If anything such a profile may give a slight overestimate of DM density at the center of the halo. This would further reduce the effective of DMA heating and make it even less likely that any "dark stars" are formed.

A bigger potential uncertainty in our adopted DM profile, however, is the assumption of spherical symmetry. Before the formation of the first sink particle this is a reasonable simplification, as the baryon profile itself is smooth and centrally concentrated. However, after the first sink particle is formed, this rapidly changes within the central regions of the halo. A disk quickly develops, and strong spiral arm features develop within the disk. It is possible that interactions between these high asymmetric features and inhomogeneities in the DM density distribution could lead to a transfer of energy to the DM component and a flattening of the central DM distribution, in a similar manner to that thought to occur from a rotating bar in disk galaxies (Weinberg & Katz 2002). To fully resolve these issues, calculations with a live three-dimensional halo component would be needed.

Another important issue is whether the central protostar remains within the center of the DM distribution, as this affects the amount of DM that can be captured by the protostar by particle scattering, which strongly influences the lifetime of any dark star phase (Iocco et al. 2008; Taoso et al. 2008a; Yoon et al. 2008). In our calculations, we find that the central protostar moves by a significant amount after its formation, as a result of the fact that accretion onto the protostar does not occur in a perfectly symmetrical fashion, and also because it is being acted on by torques from the asymmetrical gas distribution in the protostellar accretion disk. We have assumed that the DM responds in a similar fashion and that the peak in the DM distribution follows the central protostar. However, if an offset were to develop between the DM cusp and the baryon peak, dynamical heating of the DM would occur and the resulting reduction in the DM density would reduce the amount of annihilation.

This issue has been addressed recently by Stacy et al. (2012), who performed a calculation in which they resimulated the inner regions of a halo with a live DM halo of similar resolution to that of the gas after the formation of an initial sink particle. The effects of DMA were not, however, included. They formed additional sink particles at around 1000 AU and found that the motion of the star–disk system became displaced from the DM density peak and that after 5000 years there was insufficient DM to influence the formation and evolution of the protostars. However, Stacy et al. (2012) do not resolve protostar formation in the inner disk around the central object and do not include DMA and so cannot comment on fragmentation in this regime. In our runs without DMA, secondary fragmentation occurred on timescales of only a few hundred years. It is therefore likely that there will be sufficient DM during the period in which the inner disk is prone to fragmentation to maintain the findings of this paper.