Does the structure of Pop III supernova ejecta affect the elemental abundance of extremely metal-poor stars?

We utilize the SPH/$N$-body hydrodynamics code gadget-2 (Springel, 2005), while considering non-equilibrium chemistry and cooling. The hydrodynamics of the gas component is solved with the standard SPH scheme with a fixed number of neighboring particles $N_{\rm ngb}=64\pm 8$. To alleviate spurious surface tension on contact discontinuities (Saitoh & Makino, 2013; Hopkins, 2013; Wadsley, Keller & Quinn, 2017), we adopt the shared timestep strategy, where the physical variables of all SPH particles are updated with a global timestep (Saitoh & Makino, 2009). The timestep is calculated as $\underset{i}{\min}\{\Delta t_{{\rm acc},i},\Delta t_{{\rm cou},i}\}$, where $\Delta t_{{\rm acc},i}$ and $\Delta t_{{\rm cou},i}$ are the acceleration and Courant timescales of a gas particle $i$, respectively.

We solve the chemical networks of 53 reactions for 15 species, e^-, H⁺, H, H^-, H${}_{2}^{+}$, H₂, D⁺, D, D^-, HD⁺, HD, He²⁺, He⁺, He, and HeH⁺. Our chemical model includes collisional ionization/recombination of H/He and the H^--process, H${}_{2}^{+}$-process, and three-body reactions for the formation of molecular hydrogen. We then calculate the cooling rates of inverse Compton, bremsstrahlung, ionization/recombination and collisional excitation of H/He, and ro-vibrational transition of H₂/HD molecules. We here ignore C and O fine-structure cooling, which is dominated by HD cooling for metallicities ${\rm[{C}/H]}\lesssim-3$ (Omukai et al., 2005; Jappsen et al., 2007).

We follow the dispersion of metals with Lagrangian tracer particles. The velocity ${\mbox{\boldmath$v$}}_{j}$ of a tracer particle $j$ are interpolated from the velocity ${\mbox{\boldmath$v$}}_{i}$ of neighboring SPH particles as

where $m_{i}$ and $\rho_{i}$ are the mass and density of a gas particle $i$, respectively, and $W(r_{ij},\tilde{h}_{j})$ is the cubic-spline kernel function of the distance $r_{ij}$ between particles $i$ and $j$ (Springel, 2005). The smoothing length $\tilde{h}_{j}$ is estimated so that a sphere with the radius $\tilde{h}_{j}$ contains four neighboring SPH particles. The number of neighboring particles for the velocity evaluation is smaller than the density evaluation to capture steep velocity gradients around SN shocks.

We create initial conditions for a cosmological zoom-in simulation with the music code (Hahn & Abel, 2011). We initialize the simulation at $z_{\rm ini}=99$ with cosmological parameters taken from Planck Collaboration et al. (2016). First, we run a dark matter simulation with $512^{3}$ particles in a $10h^{-1}$ Mpc (comoving) periodic box. We refine dark matter particles around MHs identified with a friends-of friends (FOF) algorithm. Then, we restart the simulation from $z_{\rm ini}$ including baryons. The minimum mass of dark matter and baryon particles are $m_{\rm p,gas}=4.45\ {\rm M_{\bigodot}}$ and $m_{\rm p,dm}=16.5\ {\rm M_{\bigodot}}$, respectively. A MH forms at redshift $z_{\rm form}=28.5$ with a virial radius $R_{\rm vir}=70.1$ pc and mass $M_{\rm vir}=3.33\times 10^{5}\ {\rm M_{\bigodot}}$. We cut out a spherical region with a radius 2.5 kpc centered on the center-of-mass of the MH, which contains 9,994,502 dark matter particles and 9,972,070 gas particles.

When the maximum density reaches $n_{\rm H,cen}=10^{3}\ {\rm cm^{-3}}$, we put a star particle with a mass $M_{\rm PopIII}=25\ {\rm M_{\bigodot}}$ at the center-of-mass of a region with densities $>n_{\rm H,cen}/3$. We solve radiation transport from the star with a scheme of Susa (2006). The emission rate of hydrogen ionization photons is set to $Q({\rm H})=7.58\times 10^{48}\ {\rm s}^{-1}$ (Schaerer, 2002). Fig. 1 shows the density, temperature, and radial velocity of the MH, and H ii abundance relative to hydrogen nuclei at the end of the lifetime of the star $t_{\rm life}=6.46$ Myr, corresponding to a redshift $z=27.3$. The virial mass and radius of the halo grow to $M_{\rm vir}=5.11\times 10^{5}\ {\rm M_{\bigodot}}$ and $R_{\rm vir}=87.6$ pc, respectively, through smooth mass accretion. The gas is partially ionized within $\sim 1$ pc from the star, where density, temperature, and H⁺ abundance are $n_{{\rm H}}\sim 300\ {\rm cm^{-3}}$, $T\sim 5000$ K, and $y({\rm H^{+}})\sim 0.1$, respectively. The radius of the H ii region is smaller than the halo virial radius because the expansion of ionizing front is halted by the gas accretion onto the MH (see also Kitayama et al., 2004).²²2The simulation is the same as the MH1-C25 run in C18 until the SN explosion occurs.

Note that the radius of the ionized region in this simulation is smaller than the Strömgren radius $R_{\rm St}=36\left(Q({\rm H})/8\times 10^{48}\ {\rm s}^{-1}\right)\left(n_{\rm e}/1\ {\rm cm^{-3}}\right)^{-2/3}\ {\rm pc}$ because of the resolution limit as we discuss in Section 4.2.1. Also, we here consider the formation of a single Pop III star in the MH although the fragmentation of accretion disks can lead to multiple Pop III star formation (e.g., Turk et al., 2009; Clark et al., 2011; Greif et al., 2012). In addition, we prevent star formation in the other MHs by switching off gas cooling in dense regions with $>10^{3}\ {\rm cm^{-3}}$ at distances $>100$ pc from the central MH. Lastly, we do not consider a relative velocity offset between dark matter and baryons (“streaming velocity”). This would play an important role in the structure formation and star formation (Chiou et al., 2018, 2019; Druschke et al., 2019). The validity of these assumptions is discussed in Section 4.2.3.

From the snapshot at $t_{\rm life}$ (Fig. 1), we run four simulations with explosion energies $E_{{\rm SN}}=0.7$, $1$, $5$, and $10$ B ($1\ {\rm B}=10^{51}\ \rm erg$), hereafter called as E0.7, E1, E5, and E10, respectively. We replace the star particle with a Pop III remnant particle with a mass $M_{\rm rem}$ (Table 1) and uniformly inject the explosion energy $E_{{\rm SN}}$ to central 200 SPH particles in a sphere with a radius $R_{\rm eje,sim}=2.77$ pc. We deposit all of the explosion energy as thermal energy.

In this region, we insert $10^{6}$ metal particles, and the elemental mass fractions of each particle are remapped from one-dimensional stellar evolution/nucleosynthesis models of Tominaga et al. (2014). In their models, heavier elements than He are contained in only a small part of the ejecta. The metal-rich region with $X({\rm{H}})+X({\rm{He}})<0.95$ (hereafter called “CO core”) lies only within the radius $R_{\rm met}$, at least $\sim 4$% of entire ejecta radius $R_{\rm eje}$ (Table 1), where $X({\rm{M}})$ is the mass fraction of an element M. If we include the all part of ejecta, the velocity of metal particles in the CO core would be interpolated from only $200\times 0.04^{3}=0.01$ SPH particles (Table 1). Therefore, we consider only the CO core in our simulations. Fig. 2 shows the mass fractions of major elements in the CO core calculated by Tominaga et al. (2014).

We note that 2.34%, 13.1%, 13.1%, and 14.8% of C is produced in the outer region containing mainly primordial elements (called “hydrogen envelope”) for E0.7, E1, E5, and E10, respectively. This means 0.01–0.07 dex loss of carbon mass. We also note that only $1.18\times 10^{-5}$, $2.09\times 10^{-4}$, $3.65\times 10^{-4}$, and $8.66\times 10^{-4}$ of N is produced in the CO core. Thus, even without the hydrogen envelope, the radial velocity of the CO core is not overestimated because of the deceleration from the swept up materials.

Table 1 shows the CO core mass $M_{\rm met}$ and C, O, and Fe mass ($M_{\rm C}$, $M_{\rm O}$, and $M_{\rm Fe}$, respectively). The dominant element is O for all $E_{{\rm SN}}$. The mass of lighter elements decreases while the mass of heavier elements increases with increasing $E_{{\rm SN}}$. The elemental abundance ratio of the uniform ejecta, which is often used when fitting EMP stars, are calculated from the total mass of each element produced in the ejecta. We hereafter attach the index ‘0’ to depict the average elemental abundance ratios. In our normal CCSN model, the average carbon-to-iron abundance ratio ${\rm[{C}/Fe]}_{0}$ increases from $-0.65$ for E10 to $0.48$ for E0.7, below the definition of CEMP stars (${\rm[{C}/Fe]}=0.7$). Note that 80% of EMP stars are C-enhanced stars with metallicities ${\rm[{Fe}/H]}<-4.0$ (Yoon et al., 2018; Norris & Yong, 2019). We discuss the metal enrichment scenario from a different SN model which can explain C-enhanced star formation in Section 4.2.4.

We terminate the simulations at the time $t_{\rm recol}$ when the maximum density of an enriched gas cloud reaches $n_{\rm H,cen}=10^{3}\ {\rm cm^{-3}}$. We compare the elemental abundance within half of a Jeans length

with the abundances from uniformly mixed ejecta.

Figs. 3 and 4 show the gas density at 0.1 and 1 Myr after the SN explosion for E0.7 and E10. The blue and red dots depict the distribution of metal particles with $X({\rm{C}})>X({\rm{Fe}})$ and $X({\rm{C}})<X({\rm{Fe}})$, respectively. In all progenitor models, a part of ejecta falls back mainly along the cosmological filaments after the SN shells lose its thermal energy through radiative cooling in the dense contact surfaces with the filaments (Fig. 1). The gas falling back into the MH begins to collapse again, i.e., internal enrichment (IE) occurs at $t_{\rm recol}=5.69$–37.2 Myr for different explosion energies (Table 2). The darker colors in Figs. 3 and 4 depict the particles which eventually fall back into the MH. These particles are confined in the contact surfaces.

Fig. 5 shows the radial distribution of density, $A({\rm{C}})$, ${\rm[{Fe}/H]}$, and ${\rm[{C}/Fe]}$ of the gas particles as a function of distance from the density maximum of the enriched clouds. We also show half of the Jeans length $R_{\rm J}$ (black dashed line) and the resolution limit defined as twice the smoothing length (red dotted line). The carbon and iron abundances show large deviations at distances $>10$ pc because of the non-uniform metal mixing in the voids. Toward the higher density region at distances $<3$ pc, the metal abundances converge. In Fig. 6 we compare the timescale for metal mixing due to turbulence, $t_{\rm turb}=R/\sigma_{v}$, with the free-fall timescale, $t_{\rm ff}=(3\pi/32G\bar{\rho})^{1/2}$, where $\sigma_{v}$ and $\bar{\rho}$ are the velocity dispersion and mean density within a radius $R$, respectively. Within $R\lesssim 3$ pc the mixing timescale becomes longer than the collapsing timescale, which indicates that further metal mixing does not occur as pointed out by earlier studies (Smith et al., 2015; Ritter et al., 2015; Chiaki & Wise, 2019).

We measure the average abundance of each element within $R_{\rm J}\sim 3$ pc, regarding it as the characteristic abundance of the clouds. The red and black curves in Fig. 7 show the elemental abundance ratios ${\rm[{M}/Fe]}$ for the stratified and uniform ejecta, respectively. The difference $\Delta{\rm[{M}/Fe]}$ of these abundances is within $\pm 0.4$ dex for all elements for all progenitor models (bottom panels of Fig. 7). We can conclude that the stratified structures of the ejecta hardly affect on the elemental abundance ratio of the internally enriched clouds, which we will further analyze in Section 3.2.1 (i).

For E10, metals reach the neighboring halo, and external enrichment (EE) occurs (Fig. 8) as well as IE. The blue curve of Fig. 9 shows that only the outer layers of ejecta with mass coordinates $M_{r}>2.5\ {\rm M_{\bigodot}}$ reach the neighboring halo. Since these layers mainly contain lighter elements (bottom panel of Fig. 7), the relative abundances of lighter elements ($<{\rm Ca}$) are enhanced more than $\Delta{\rm[{M}/Fe]}>1$. This indicates that stars with the C-enhanced abundance (${\rm[{C}/Fe]}=1.49$) can form in externally enriched clouds from normal CCSNe.

If low-mass stars ($M_{*}<0.8\ {\rm M_{\bigodot}}$) form in the enriched clouds, they can survive for the Hubble time (Smith et al., 2015; Chiaki & Wise, 2019). Then, if some of the stars are accreted onto our Milky Way halo or local dwarf galaxies through cosmic structure formation, they can be observed as EMP stars (Frebel & Norris, 2015, and references therein). The colored symbols in Fig. 10 show the distribution of the stars on the $A({\rm{C}})$-${\rm[{Fe}/H]}$ plane. The closed and open symbols depict the abundances of stars forming in the internally and externally enriched clouds (hereafter called IE-stars and EE-stars), respectively. The small and large symbols depict the abundances for the uniform and stratified ejecta, respectively. Interestingly, the Fe abundances are shifted while C abundances hardly change. The mass of the outer layers is larger than that of the inner layers, where the fall back fraction fluctuates, and the mass-weighted average fraction $\langle f_{\rm enr}\rangle$ is closer to the fraction $f_{\rm enr}(M_{\rm r})$ in the outer layers (Fig. 9).

For comparison, we also plot the C and Fe abundances of observed C-normal and C-enhanced stars with blue and red dots, respectively (taken from the SAGA database; Suda et al., 2008). The stars with below and above ${\rm[{C}/Fe]}=0.7$ (dashed line) are defined as C-normal stars and CEMP stars, respectively (Aoki et al., 2007). C-normal stars are distributed around the line ${\rm[{C}/Fe]}=0$ (dot-dashed line). CEMP stars with ${\rm[{Fe}/H]}<-3$ are classified into two sub-groups according to their distributions (Yoon et al., 2016). Group III stars are distributed around the horizontal line with a constant $A({\rm{C}})$ while Group II stars are distributed along the line with a constant ${\rm[{C}/Fe]}$ and just above the region where C-normal stars are distributed. The $A({\rm{C}})$ show the large star-to-star scatter of $\sigma(A({\rm{C}}))\sim 1$ dex with a fixed metallicity for all groups. This can be considered to reflect the variation of mass and explosion energy of few progenitor stars in the early stage of galactic chemical evolution (Ryan et al., 1996; Cayrel et al., 2004).

From our simulations, the IE-stars are distributed around the line ${\rm[{C}/Fe]}=0$ (dot-dashed line) with $\sim 1$ dex deviation. This is consistent with the distribution of observed C-normal and Group II stars. For $E_{{\rm SN}}\geq 5$ B (green circle and blue square), the IE-stars have low ${\rm[{C}/Fe]}$ values (bottom right in Fig. 10) because ${\rm[{C}/Fe]}$ for a uniform ejecta ($-0.18$ and $-0.65$) is smaller than the solar abundance ratio, and is unchanged or depleted by the effect of the stratified ejecta. For $E_{{\rm SN}}\leq 1$ B (red diamond and orange triangle), ${\rm[{C}/Fe]}$ is larger (0.48 and 0.31) for the uniform ejecta and further enhanced for the stratified ejecta (see Section 3.2.1 ii). The stars are distributed around the critical line for CEMP stars ${\rm[{C}/Fe]}=0.7$ (dashed line). For $E_{{\rm SN}}=0.7$ B, the stars would be classified as CEMP stars (red diamond). In this work, we fix the progenitor mass to $25\ {\rm M_{\bigodot}}$, but the range of explosion energies (0.7–10 B) can explain the distribution of observed C-normal and CEMP Group II stars with a scatter of $\sim 1$ dex.

For EE-stars, both $A({\rm{C}})=2.48$ and ${\rm[{Fe}/H]}=-7.44$ are smaller than the ones of IE-stars by two orders of magnitude because metals are mixed with a large mass of pristine gas during a merger of three host halo progenitors which induces EE (Section 3.2.2). As a result, $A({\rm{C}})$ and ${\rm[{Fe}/H]}$ of the EE-stars exist in the low metallicity region where stars have not been observed (grey shaded region in Fig. 10). We can consider several reasons why these EE-stars with extremely small metallicities have not been observed. First, because the metallicity for IE is larger than that for EE, the C-enhanced abundance pattern resulting from EE may be washed out by the IE of the neighboring cloud itself. Halos only externally enriched may be so rare that EE-stars have not so far been observed among the current samples of EMP stars (Hartwig et al., 2018). Second, clouds with insufficient metal content ($A({\rm{C}})\lesssim 6$) should contain a negligible amount of dust grains, which is the main coolant that can induce cloud fragmentation. Chiaki et al. (2017) estimate the elemental abundances below which dust cooling works inefficiently as $10^{{\rm[{C}/H]}-2.30}+10^{{\rm[{Fe}/H]}}=10^{-5.07}$ (purple curve in Fig. 10). The EE-stars are plotted below this curve, where clouds can collapse stably against fragmentation, and a single massive star ($>10\ {\rm M_{\bigodot}}$) is likely to form. Such massive stars can not be observed because they will not survive for 13.6 Gyr between the recollapsing redshift $z_{\rm recol}=22.5$ and the present day.

In this work, we perform numerical simulations of the transition from Pop III stars to Pop II stars, considering the non-uniform structure of SN ejecta with realistic nucleosynthesis models. Because of computational limits, there are several caveats in the simulations.