COMPACT STELLAR BINARY ASSEMBLY IN THE FIRST NUCLEAR STAR CLUSTERS AND r-PROCESS SYNTHESIS IN THE EARLY UNIVERSE

The basic assumption of our modeling is that compact star clusters can be formed at the center of high-z dark matter (DM) halos with a virial temperature ${{T}_{{\rm vir}}}\sim (1-2)\times {{10}^{4}}\;{\rm K}$ and non-zero but very low metallicity ${{{\rm log} }_{10}}(Z/{{Z}_{}})\sim -4$. In fact, under these conditions disk fragmentation can be suppressed and gas efficiently funneled by the Toomre instability in the innermost few parsecs where it forms a star cluster (see Devecchi & Volonteri 2009 for details).

To estimate the formation rate of early forming compact and high-density star clusters, we derive the formation rate of DM halos meeting the criteria above from the simulations of Trenti & Shull (2010), which focus on early chemical enrichment of three Milky-Way-like halos with DM mass ${{M}_{{\rm DM}}}(z=0)\approx {{10}^{12}}\;{{M}_{\odot }}$. The simulations have been carried out using ${{N}_{{\rm DM}}}={{1024}^{3}}$ DM particles in a comoving volume of ${{10}^{3}}\;{\rm Mp}{{{\rm c}}^{3}}$, with DM mass resolution of $3.4\times {{10}^{4}}\;{{M}_{\odot }}$. The DM only simulation is post-processed with the star formation and chemical enrichment model described in Trenti & Shull (2010) and Trenti et al. (2009), which includes treatment of metal outflows and radiative feedback in the Lyman–Werner bands to regulate Population III formation (Trenti & Stiavelli 2009).

Because, as shown in Section 2.2, the encounter rate in the core of the first star clusters declines significantly at redshifts $z\lt 10$ (mainly due to a decrease in the central density), we focus only on systems formed before then, which are shown in Figure 1. By z = 10, a total of $\langle {{N}_{c}}\rangle =515$ clusters are expected to be formed in a region that will collapse by z = 0 in a Milky-Way-like halo. This number is higher by a factor of about two when compared to the formation rate of similar objects in a random region of the universe, owing to the enhanced DM halo formation rate induced by the over dense environment (see, e.g., Trenti & Shull 2010).

Zoom In Zoom Out Reset image size

Overall our simplified approach provides a robust estimate of the formation rate of halos capable of hosting dense star clusters. For example, feedback at $z\gt 10$ is unlikely to affect gas in halos with ${{T}_{{\rm vir}}}\gt {{10}^{4}}$ K. This leaves the leading source of uncertainty in the details of star formation in such halos, which are very challenging to model from first principles (e.g., Wise et al. 2012). Thus, as a first characterization, we describe the internal properties of first stellar clusters at z ≈ 10 following Devecchi & Volonteri (2009). From Figure 4 in Devecchi & Volonteri (2009) we derive the characteristic stellar mass ${{M}_{{\rm tot}}}$ and half-mass radius ${{r}_{{\rm hm}}}$ for the central star clusters as a function of their formation redshift. While the core density during initial stages of cluster evolution is fairly uncertain and dependent upon the initial conditions, as a basic estimate, based on direct N-body modeling experiments (Trenti et al. 2010), we assume a typical core to half mass radius ratio ${{r}_{c}}/{{r}_{{\rm hm}}}=0.1$ and a core to characteristic stellar mass ratio ${{M}_{c}}/{{M}_{{\rm tot}}}=0.04$, corresponding to a pre-core collapse cluster with a W₀ = 8.5 King (1966) profile. With the model (in virial equilibrium), we derive the central velocity dispersion from the cluster mass and radius:

$$\begin{eqnarray}&&{{\sigma }_{c}}=18.81\;{\rm km}\;{{{\rm s}}^{-1}}\left( \frac{{{M}_{{\rm tot}}}}{{{10}^{5}}\;{{M}_{\odot }}} \right)\left( \frac{1\;{\rm pc}}{{{r}_{{\rm hm}}}} \right).\end{eqnarray} \tag{ 1 }$$

Having an estimate for the mass and velocity dispersion of a cluster (Figure 2), the number of NSs can be calculated by assuming the star clusters have a Kroupa (2001) IMF $\xi (m)={{\xi }_{0}}{{(m/0.5\;{{M}_{\odot }})}^{\alpha (m)}}$, with $m\in [{{m}_{{\rm min} }}:{{m}_{{\rm max} }}]$ and either a cut-off at ${{m}_{{\rm min} }}=0.2\;{{M}_{\odot }}$ or at ${{m}_{{\rm min} }}=1\;{{M}_{\odot }}$ (top-heavy IMF). ${{m}_{{\rm max} }}=100\;{{M}_{\odot }}$, $\alpha (m)=-2.3$ for $m\gt 0.5\;{{M}_{\odot }}$ and $\alpha (m)=-1.3$ for $0.2\;{{M}_{\odot }}\leqslant m\leqslant 0.5\;{{M}_{\odot }}$. By assuming, for example, that NSs are produced in the mass range $[8:25]\;{{M}_{\odot }}$ (Hurley et al. 2000), the fraction of NSs is ${{f}_{{\rm ns}}}=9\times {{10}^{-3}}$ for the standard IMF and ${{f}_{{\rm ns}}}=5\times {{10}^{-2}}$ for the top heavy IMF.

Zoom In Zoom Out Reset image size

Next, we estimate the fraction of these NSs that are retained in the cluster core. The core retention fraction is affected by NS natal kicks and binary and mass segregation. We assume that the progenitors of NSs are generally in binaries since massive stars are expected to form with high binarity (Krumholz et al. 2009), and use the model of Pfahl et al. (2002, Figure 12) to estimate the retention of NSs formed in clusters with different escape velocities. We consider a Maxwellian kick distribution, $p(v)=\sqrt{2/\pi }{{v}^{2}}/\sigma _{{\rm kick}}^{3}{{e}^{-{{v}^{2}}/(2\sigma _{{\rm kick}}^{2})}}$, with ${{\sigma }_{{\rm kick}}}=100\;{\rm km}\;{{{\rm s}}^{-1}}$ and ${{\sigma }_{{\rm kick}}}=200\;{\rm km}\;{{{\rm s}}^{-1}}$, which are the fastest kick distribution studied by Pfahl et al. (2002). Since that study concludes that even with their optimistic assumptions they might underestimate the actual retention fraction (Section 7 in Pfahl et al. 2002), our approach is fairly conservative. We assume a cluster escape speed of ${{v}_{{\rm esc}}}(0)=4.12{{\sigma }_{c}}$ as derived from the King density profile. In what follows we will assume that 50% of the retained NSs are unbound from their parent binary while 50% remain bound, as roughly found by Pfahl et al. (2002) for the range of cluster's escape velocities derived here. Mass segregation can enhance the number of NSs and NS progenitors in the core relative to the cluster mean. From N-body simulations with stellar evolution (M. MacLeod et al. 2015, in preparation), we see that mass segregation of the NS progenitors is efficient (core density increased by ≈4), but the core density of NSs is increased only by ≈50%, owing to the redistribution following natal kick. We thus apply this latter correction factor (that is 1.5×) to map the global to core retention fraction (${{\zeta }_{{\rm ret}}}$), which is then shown in panel (d) of Figure 2. ${{\zeta }_{{\rm ret}}}\sim 0.2$ for fast kicks (${{\sigma }_{{\rm kick}}}=200\;{\rm km}\;{{{\rm s}}^{-1}}$) at all redshifts, while for the distribution with slower kicks ${{\zeta }_{{\rm ret}}}\sim 0.8$ at z = 14, decreasing to ${{\zeta }_{{\rm ret}}}\sim 0.55$ by $z\sim 10$.

To calculate the encounter rate within the cluster we assume NSs are distributed homogeneously within the core, with fractional number ${{f}_{{\rm ns}}}$, total number density ${{n}_{c}}$, and binary fraction ${{b}_{{\rm ns}}}$. We further assume that the stars follow a Maxwellian velocity distribution function with dispersion ${{\sigma }_{c}}$. We calculate a fiducial encounter rate ${{\xi }_{{\rm ns}}}$ for individual clusters provided a distance of closest approach ${{R}_{{\rm min} }}$ (e.g., Lee et al. 2010)

$$\begin{eqnarray}\begin{array}{rcl} {{\xi }_{{\rm ns}}} & = & 3\times {{10}^{-3}}\ {\rm Gy}{{{\rm r}}^{-1}}f_{{\rm ns}}^{2}{{\left( \frac{{{n}_{c}}}{{{10}^{6}}\ {\rm p}{{{\rm c}}^{-3}}} \right)}^{2}}{{\left( \frac{{{r}_{c}}}{0.1\ {\rm pc}} \right)}^{3}} \\ {} & {} & \times {{\left( \frac{{{\sigma }_{c}}}{20\ {\rm km}\;{{{\rm s}}^{-1}}} \right)}^{-1}}\left( \frac{{{M}_{{\rm ns}}}}{1.4\ {{M}_{\odot }}} \right)\left( \frac{{{R}_{{\rm min} }}}{10\ {\rm km}} \right). \\ \end{array}\end{eqnarray} \tag{ 2 }$$

To estimate the rate of single–single NS encounters, we multiply by the single NS fraction squared, $1/2{{(1-{{b}_{{\rm ns}}})}^{2}}$, where the factor of $1/2$ avoids double counting pairwise encounters. To estimate binary–single encounters, we must include a factor ${{b}_{{\rm ns}}}(1-{{b}_{{\rm ns}}})$. To form a binary by transferring orbital energy to internal stellar oscillations, we adopt the cross section for tidal capture formalism described in Lee & Ostriker (1986) and Kim & Lee (1999), which in the case of NSs gives

$$\begin{eqnarray}&&R_{{\rm min} }^{{\rm tidal}}=13.37{{R}_{{\rm ns}}},\end{eqnarray} \tag{ 3 }$$

(Lee et al. 2010). To compute the cross section for binary formation by GW radiation, we follow Lee (1993). That is:

$$\begin{eqnarray}&&R_{{\rm min} }^{{\rm gw}}=1458{{\left( \frac{{{\sigma }_{c}}}{20\;{\rm km}\;{{{\rm s}}^{-1}}} \right)}^{-4/7}}{\rm km}.\end{eqnarray} \tag{ 4 }$$

This is significantly larger than $R_{{\rm min} }^{{\rm tidal}}=160({{R}_{{\rm ns}}}/12\;{\rm km})\ {\rm km}$. Although, as we show in Section 3, tidal capture encounters eject significantly more r-process material.

Finally, we consider encounters between one binary containing a NS with a single NS leading to an exchange, and subsequent formation of a NS–NS pair. Using exchange cross sections from Heggie et al. (1996), the cross section for a $1.4\;{{M}_{\odot }}$ NS to exchange into a binary given an equal mass companion is ${{\sigma }_{{\rm ex}}}=32.2\;{\rm A}{{{\rm U}}^{2}}(a/{\rm AU}){{(v/20\;{\rm km}\;{{{\rm s}}^{-1}})}^{-2}}$. If the eccentricity is sufficiently high, the pair will merge on a short timescale (Samsing et al. 2014). To calculate the cross section for the formation of NS binaries with lifetime less than 1 Gyr, we use the lifetime estimate ${{t}_{{\rm life}}}\approx 2.9\times {{10}^{17}}{{(1-{{e}^{2}})}^{7/2}}\;{\rm yr}$ and assume thermal distribution of eccentricity, $p(e)=2e$, following the exchange (Samsing et al. 2014). Thus, for any binary semimajor axis, a, there is a critical eccentricity ${{e}_{{\rm crit}}}$ above which newly formed binaries will merge in less than 1 Gyr. The fraction of binaries with $e\gt {{e}_{{\rm crit}}}$ is thus $1-e_{{\rm crit}}^{2}$ and scales as ${{a}^{-8/7}}$. The resultant cross section scales weakly with binary semimajor axis, ${{a}^{-1/7}}$ (e.g., Samsing et al. 2014, Figure 15). For a representative initial binary population with a binary semimajor axis distribution with $p(a)\propto {{a}^{-1}}$ and $a\in [10\;{{R}_{\odot }}:1\;{\rm AU}]$, the integrated cross-section can be expressed as

$$\begin{eqnarray}&&R_{{\rm min} }^{3-{\rm body}}=1.6\times {{10}^{6}}\;{\rm km},\end{eqnarray} \tag{ 5 }$$

in Equation (2).¹⁰

Equipped with $R_{{\rm min} }^{{\rm tidal}}$, $R_{{\rm min} }^{{\rm gw}}$, and $R_{{\rm min} }^{3-{\rm body}}$, we can calculate the binary assembly rate ${{\xi }_{{\rm ns}}}$ for individual clusters (Equation (2)) as a function of their formation redshift, which, as shown in Figure 2, sets their structural properties. The rates for tidal, GW captured, and dynamically formed binaries within individual clusters are given in panel (a) of Figure 3, where we assume ${{b}_{{\rm ns}}}=0.5$. As expected, tidal encounters are very rare and the rate is dominated by binary–single interactions, which have been proposed as viable merger channel by several studies (Sigurdsson & Phinney 1993; Grindlay et al. 2006; Samsing et al. 2014). In this work we ignore binary–binary encounters and their exchange products, which will represent, at most, a contribution of order unity to the binary–single encounter rate.

Zoom In Zoom Out Reset image size

An estimate of the integrated encounter rate ${{\eta }_{{\rm ns}}}$ over all the Milky Way central star clusters present at high z can then be made under the assumption that each cluster contributes to the encounter rate from its redshift of formation until its dissolution. Conservatively, we assume that dissolution happens 0.5 Gyr after cluster formation, although some clusters may well survive for longer and therefore provide additional contribution to the encounter rate.¹¹ Our estimate of ${{\eta }_{{\rm ns}}}$ for dynamically formed binaries as a function of z is shown in panel (b) of Figure 3. From this panel it is possible to compute the total number of mergers (which will follow within 1 Gyr after binary capture) in the Milky Way by integrating the rate over time. In panel (c) of Figure 3 we also plot for clarity the fraction of high-z star clusters that experienced a merger as a function of z derived by dividing the number of expected mergers by the number of clusters present. To calculate the total amount of r-process material synthesized in these clusters we need to combine the rate calculations presented here with an estimate of the mass production rate per event. It is to this issue that we now turn our attention.