Origin of highly r-process-enhanced stars in a cosmological zoom-in simulation of a Milky Way-like galaxy

ABSTRACT

The r-process-enhanced (RPE) stars provide fossil records of the assembly history of the Milky Way (MW) and the nucleosynthesis of the heaviest elements. Observations by the R-Process Alliance (RPA) and others have confirmed that many RPE stars are associated with chemo-dynamically tagged groups, which likely came from accreted dwarf galaxies of the MW. However, we do not know how RPE stars are formed. Here, we present the result of a cosmological zoom-in simulation of an MW-like galaxy with r-process enrichment, performed with the highest resolution in both time and mass. Thanks to this advancement, unlike previous simulations, we find that most highly RPE (r-II; [Eu/Fe] > +0.7) stars are formed in low-mass dwarf galaxies that have been enriched in r-process elements for [Fe/H] |$\lt -2.5$|⁠, while those with higher metallicity are formed in situ, in locally enhanced gas clumps that were not necessarily members of dwarf galaxies. This result suggests that low-mass accreted dwarf galaxies are the main formation site of r-II stars with [Fe/H] |$\, \lt -2.5$|⁠. We also find that most low-metallicity r-II stars exhibit halo-like kinematics. Some r-II stars formed in the same halo show low dispersions in [Fe/H] and somewhat larger dispersions of [Eu/Fe], similar to the observations. The fraction of simulated r-II stars is commensurate with observations from the RPA, and the distribution of the predicted [Eu/Fe] for halo r-II stars matches that observed. These results demonstrate that RPE stars can be valuable probes of the accretion of dwarf galaxies in the early stages of their formation.

methods: numerical, stars: abundances, Galaxy: abundances, Galaxy: formation, Galaxy: halo, Galaxy: kinematics and dynamics

1 INTRODUCTION

The chemical abundances of stars potentially provide clues to understanding galaxy formation and nucleosynthesis. In particular, elements synthesized by the rapid neutron-capture process (r-process elements), such as Eu, Th, and U, are crucial elements for resolving early Galactic chemical evolution (e.g. Beers & Christlieb 2005; Sneden, Cowan & Gallino 2008; Frebel & Norris 2015; Thielemann et al. 2017; Frebel 2018; Cowan et al. 2021), using low-metallicity stellar probes in the Milky Way (MW; e.g. Spite & Spite 1978; Truran 1981; McWilliam 1997).

Recognition of the highly r-process-enhanced (RPE) star CS 22892-052 (Sneden et al. 1994; [Eu/Fe]¹ = +1.6 at [Fe/H] = −3.1), discovered in the HK survey (Beers, Preston & Shectman 1985, 1992), made it clear that the abundances of r-process elements in the early epochs of star formation in the MW were not homogeneous. Later observations of this star have been reported by many authors, including Sneden et al. (1996, 2003) and Barklem et al. (2005).

The general definition of RPE stars are those exhibiting enhancement of heavy r-process elements thought to be primarily produced by the r-process such as europium (Eu; Z = 63), [Eu/Fe] |$\gt \, +0.3$|⁠. Subsequent to the recognition of CS 22892-052, r-process-element abundance measurements for additional RPE stars were reported: e.g. CS 31082-001 (Cayrel et al. 2001; Hill et al. 2002; Plez et al. 2004), CS 22183-031 (Honda et al. 2004), and HE 1523-0901 (Frebel et al. 2007), among many others. Furthermore, a star-to-star scatter of [Eu/Fe] (and [Sr/Fe]) ratios (more than two orders of magnitude) has been confirmed for stars with |$\mathrm{[Fe/H]}\, \lt -$|2.5 (e.g. McWilliam et al. 1995; Ryan, Norris & Beers 1996; François et al. 2007; Ishigaki, Aoki & Chiba 2013; Roederer et al. 2014). Most recently, ground- and space-based observations of the bright metal-poor (V = 9, [Fe/H] = −1.46) highly RPE ([Eu/Fe] = +1.32) horizontal-branch star HD 225925 has produced the largest inventory of measured elemental abundances for any star other than the Sun (Roederer et al. 2022), and provided a new ‘template star’ for the observed distribution of r-process elements.

1.1 The classification of RPE stars and their observed environments

RPE stars are conventionally sub-classified by their [Eu/Fe] ratios, originally defined by Beers & Christlieb (2005) in order to separate moderately and highly enhanced stars (r-I: + 0.3 |$\lt $| [Eu/Fe] |$\le $| + 1.0, [Ba/Eu] |$\lt$| 0; r-II [Eu/Fe] |$\gt$| + 1.0, [Ba/Eu] |$\lt$| 0). The condition [Ba/Eu] |$\lt$| 0 is used in order to select stars dominantly enriched by the r-process, rather than the s-process.² Holmbeck et al. (2020) recently statistically reclassified RPE stars, based on the large sample of such stars reported by the RPA. They defined the r-I and r-II boundaries, using a k-medoids partitioning technique (Kaufman & Rousseeuw 1990) and a mixture-model analysis. Their revised classification of r-I and r-II stars is

$$\begin{eqnarray*} r\text{-I}&:& +0.3\, \lt \, \mathrm{[Eu/Fe]}\, \le \, +0.7,\ \mathrm{[Ba/Eu]}\, \lt \, 0,\\ r\text{-II}&:& \mathrm{[Eu/Fe]}\, \gt \, +0.7,\ \mathrm{[Ba/Eu]}\, \lt \, 0. \end{eqnarray*} $$

This classification was defined to minimize the distance between the members of three groups (the third group being those stars with [Eu/Fe] ≤ +0.3) to determine the group centres in the [Eu/Fe] distribution. They also found that stars classified into these three groups more appropriately represented the sample [Eu/Fe] distribution than those classified into two groups, based on the Akaike information criterion (Akaike 1973). Recently, Cain et al. (2020) defined r-III stars for the handful of stars with [Eu/Fe]|$\, \gt \,$| + 2.0 and [Ba/Eu]|$\, \lt -0.5$|⁠.

Since the r-process-element abundances of r-II stars are well above the average values in the MW, they are ideal probes to investigate their formation environment. As we point out below, the formation of r-I stars in our simulation is heavily affected by the assumption of the appropriate yields, which is highly uncertain. Thus, in this paper, we concentrate on the r-II stars (for simplicity, we combine these stars with the few r-III stars known and simply refer to their union as ‘r-II’ stars). The r-I stars will be considered in detail in a forthcoming paper.

RPE stars are also found in dwarf galaxies. The ultrafaint dwarf (UFD) galaxy Reticulum II is a notable example (Ji et al. 2016a, b, 2022; Roederer et al. 2016). Seven out of nine observed stars are classified as r-II stars in this galaxy. Ji et al. (2016a) argued that a nucleosynthetic event with a europium yield of |$10^{-4.5}\, \mathrm{M}_{\odot }$| could enhance the r-process-element abundances of gas with a mass |${\sim}10^{6}\, \mathrm{M}_{\odot }$| to the level required for r-II stars. The Tucana III UFD is another example. An r-I star, DES J235532.66-593114.9, with [Eu/Fe] = + 0.60, is confirmed in this galaxy, along with at least two other less-enhanced r-I stars (Hansen et al. 2017; Marshall et al. 2019). Several RPE stars are also found in more massive dwarf galaxies: Fornax (Letarte et al. 2010; Reichert, Hansen & Arcones 2021), Ursa Minor (Shetrone, Côté & Sargent 2001; Sadakane et al. 2004), Carina (Shetrone et al. 2003; Norris et al. 2017), Draco (Cohen & Huang 2009), Sculptor (Shetrone et al. 2003; Geisler et al. 2005) dwarf spheroidal galaxies (dSphs), and the Large/Small Magellanic Clouds (LMC/SMC; Reggiani et al. 2021). These findings clearly suggest a link between the formation of the MW and dwarf galaxies of a variety of masses.

1.2 The nature of the stellar halo of the MW

Studies of the nature of the MW halo have a long, rich history (see e.g. Beers & Christlieb 2005 for a brief overview). Two decades ago, Chiba & Beers (2000) analysed the motion of ∼1200 stars with |$\mathrm{[Fe/H]}\, \lt -0.6$|⁠, based on the kinematically unbiased sample of over 2000 metal-poor stars from Beers et al. (2000), and found that there is no correlation between stellar metallicities and their orbital eccentricities. This finding clearly indicated, for the first time, that the correlation seen in Eggen, Lynden-Bell & Sandage (1962) that led to a monolithic collapse scenario for MW formation was due to the kinematic selection bias employed to assemble their sample. Chiba & Beers concluded that the halo was complex and not a single structure with a simple formation history, initiating the discussion of the possible existence of a flattened ‘inner halo’ and more spherical ‘outer halo’ with differing formation mechanisms.

The series of papers by Carollo et al. (2007, 2010) and Beers et al. (2012) amplified these conclusions based on the much larger sample of metal-poor stars obtained during the Sloan Digital Sky Survey (York et al. 2000), in particular the stellar-specific sub-survey Sloan Extension for Galactic Understanding and Exploration (SEGUE; Yanny et al. 2009), and clearly demonstrated ‘the halo’ was, in fact, a superposition of at least two primary components – an inner and an outer halo – that exhibit different spatial density profiles, stellar orbits, and stellar metallicities. These properties indicate that the individual halo components likely formed in fundamentally different ways, through successive dissipational and dissipationless mergers and the tidal disruption of proto-Galactic clumps, refining the original suggestions from Chiba & Beers (2000).

In the context of this scenario, stars in the MW can be broadly classified as members of in situ (primarily inner-halo) and accreted (primarily outer-halo) components. In situ stars are thought to have formed in what is sometimes referred to as the ‘main halo’ of the MW (dominated by the inner-halo population), while accreted stars are expected to have formed in the halo outside of the MW’s main halo and accreted by the MW at a later time, and are expected to predominantly populate the outer halo.

1.3 Theoretical expectations and halo substructure

According to the hierarchical structure formation scenario in the Lambda cold dark matter (ΛCDM) paradigm, galaxies could be formed by the clustering of smaller systems (e.g. White & Rees 1978; Blumenthal et al. 1984). Therefore, it is crucial to find signatures of the accretion events of dwarf galaxies to the MW. Relaxation of the initial dynamical phase-space distribution is thought to exceed the age of the Universe in order to fully mix the stellar orbital energies and angular momenta. Thus, the observed phase-space distribution of stars has been used to identify substructures in the MW (and references therein Helmi 2020). From a sufficiently large sample of stars with available dynamical information, the assembly of the MW can be reconstructed. If stellar metallicities are also available, one can use this to make a rough inference of the ages of the stars involved. Results based on precision age-determination techniques (in particular, for subgiant stars) have been recently reported for halo and disc-system stars by Xiang & Rix (2022).

The astrometric data provided by the Gaia satellite mission have already revolutionized our understanding of the building blocks of the MW (Gaia Collaboration 2016a, b, 2018, 2021). Several substructures have been identified in the dynamical phase space, which could be remnants of disrupted dwarf galaxies accreted to the MW (e.g. Belokurov et al. 2018; Haywood et al. 2018; Helmi et al. 2018; Koppelman, Helmi & Veljanoski 2018; Myeong et al. 2018; Koppelman et al. 2019; Naidu et al. 2020; Yuan et al. 2020). Chemical-abundance analyses of member stars in each substructure show that some of these have chemical abundances similar to presently observed MW satellites (e.g. Matsuno, Aoki & Suda 2019; Aguado et al. 2021; Limberg et al. 2021b; Matsuno et al. 2021, 2022; Carrillo et al. 2022; Horta et al. 2022). Numerous stellar debris streams have also been reported (e.g. Mateu 2022); their full chemical compositions are only just now beginning to be studied. Notably, Gull et al. (2021) have reported the detection of stars that exhibit clear r-process patterns in several streams, including a number of RPE stars.

Statistical analyses based on the combination of spectroscopic and Gaia astrometric data for a variety of stellar samples have been used to identify dynamically tagged groups (DTGs) in the phase-space distribution (Yuan et al. 2020; Limberg et al. 2021a; Shank et al. 2022b, c). In the largest such sample studied to date, Shank et al. (2022b) applied an unsupervised learning algorithm Hierarchical Density-Based Spatial Clustering of Applications with Noise (HDBSCAN; Campello, Moulavi & Sander 2013) for orbital energy and cylindrical actions to some 8000 stars selected from the RAdial Velocity Experiment (RAVE) DR6 (Steinmetz et al. 2006, 2020a, b); they identified 179 DTGs. Most of these DTGs are associated with known substructures such as Gaia-Sausage-Enceladus (GSE; Belokurov et al. 2018; Helmi et al. 2018), the Metal-Weak Thick Disc (Morrison, Flynn & Freeman 1990; Beers et al. 2014; Carollo et al. 2019), the Splashed Disc (Di Matteo et al. 2019; Belokurov et al. 2020), Thamnos (Koppelman et al. 2019), the Helmi Stream (Helmi et al. 1999; Chiba & Beers 2000), and LMS-1/Wukong (Naidu et al. 2020; Yuan et al. 2020). They also found 22 DTGs that are associated with RPE stars.

1.4 RPE stars as probes of halo substructure

In a seminal investigation, Roederer, Hattori & Valluri (2018) showed that RPE stars tend to be associated with groups of stars with similar orbits, based on an analysis of the dynamics of 35 r-II stars using Gaia DR2 (Gaia Collaboration 2018) and previously published radial velocities and elemental abundances (see also Hattori, Okuno & Roederer 2022). They identified eight groups of 2–4 stars each, with an indication that the stars within each group had more similar abundances of [Fe/H] and [Eu/Fe] than would have been expected from random draws. Their orbital analysis also showed that their sample stars all exhibited halo-like kinematics. These results suggested that RPE stars may have arisen from disrupted dwarf galaxies similar to UFDs or low-luminosity dSphs.

The ongoing R-Process Alliance (RPA) effort has greatly increased the number of confirmed RPE stars (Hansen et al. 2018; Sakari et al. 2018; Ezzeddine et al. 2020; Holmbeck et al. 2020). Holmbeck et al. (2020) reported a total of 72 r-II stars and 232 r-I stars based on this survey alone. Gudin et al. (2021) analysed the dynamics of 446 RPE stars in the RPA and additional literature samples, again using HDBSCAN, to identify 30 dynamical groups of RPE stars with between 3 and 12 members each. Since they performed the clustering analysis exclusively for chemically peculiar stars (RPE stars), they call these groups chemo-dynamically tagged groups (CDTGs). They found that each CDTG has smaller than expected dispersions in metallicity ([Fe/H]), carbonicity ([C/Fe]), and neutron-capture-element abundance ratios ([Sr/Fe], [Ba/Fe], and [Eu/Fe]). These results strongly suggested that RPE stars within each CDTG experienced a common chemical evolution, supporting the scenario that RPE stars in CDTGs were from disrupted dwarf galaxies or globular clusters. Most recently, Shank et al. (2022a) have compiled a sample of ∼1800 RPE stars based on the RPA, GALAH DR3 (Buder et al. 2021), and other literature sources and described the abundance dispersion results for 38 CDTGs with between 5 and 20 members each.

1.5 The origin of RPE stars and their birth environments

The enrichment of r-process elements in the MW and dwarf galaxies with chemical-evolution models has been studied for over 30 yr (e.g. Mathews & Cowan 1990; Mathews, Bazan & Cowan 1992; Ishimaru & Wanajo 1999; Tsujimoto, Shigeyama & Yoshii 2000; Argast et al. 2004; Ishimaru et al. 2004; Cescutti & Chiappini 2014; Komiya et al. 2014; Matteucci et al. 2014; Tsujimoto & Shigeyama 2014; Cescutti et al. 2015; Hirai et al. 2015, 2017; Ishimaru, Wanajo & Prantzos 2015; Wehmeyer, Pignatari & Thielemann 2015; Côté et al. 2018, 2019; Ojima et al. 2018; Safarzadeh, Sarmento & Scannapieco 2019; Schönrich & Weinberg 2019; Wanajo, Hirai & Prantzos 2021). Most of these studies focused on constraining the astrophysical sites of the r-process elements by explaining the trend and dispersions of [Eu/Fe] ratios as a function of [Fe/H]. Argast et al. (2004) preferred core-collapse supernovae (CCSNe) as the astrophysical site rather than binary neutron star mergers (NSMs) in terms of the rates and delay times. On the other hand, Hirai et al. (2015) have shown that NSMs can be a major contributor of r-process elements in dwarf galaxies because of their slow chemical evolution. Evidence for the (slightly) delayed production of the r-process elements, which would come from NSMs, has been seen in dwarf galaxies and accreted components (Duggan et al. 2018; Skúladóttir & Salvadori 2020; Matsuno et al. 2021; Naidu et al. 2022).

Nucleosynthetically, NSMs are one of the most promising candidates for the astrophysical sites of the r-process (e.g. Lattimer & Schramm 1974; Symbalisty & Schramm 1982; Eichler et al. 1989; Meyer 1989; Freiburghaus, Rosswog & Thielemann 1999; Goriely, Bauswein & Janka 2011; Korobkin et al. 2012; Bauswein, Goriely & Janka 2013; Wanajo et al. 2014; Goriely et al. 2015; Sekiguchi et al. 2015, 2016; Radice et al. 2016; Fujibayashi et al. 2022). For other sites, CCSNe, driven by magneto-rotational instability (Winteler et al. 2012; Nishimura, Takiwaki & Thielemann 2015; Nishimura et al. 2017) and collapsars (Siegel, Barnes & Metzger 2019) are recently proposed, while CCSNe driven by neutrino winds have been found to have difficulty synthesizing elements heavier than A ∼ 110 (Wanajo, Janka & Müller 2011; Wanajo 2013; Wanajo et al. 2018). Although contributions from each site to the r-process enrichment are still debated, many nucleosynthesis and chemical-evolution studies support that NSMs play an important role as the astrophysical site of the r-process. Observations of the kilonova associated with the gravitational wave event GW170817 have shown a definitive astrophysical source of heavy elements created by the r-process in NSMs (e.g. Abbott et al. 2017a,b; Arcavi et al. 2017; Drout et al. 2017; Nicholl et al. 2017; Shappee et al. 2017; Tanaka et al. 2017; Watson et al. 2019; Domoto et al. 2021).

Prolific amounts of r-process elements from NSMs could form RPE stars in UFDs (Safarzadeh & Scannapieco 2017; Tarumi, Yoshida & Inoue 2020; Jeon, Besla & Bromm 2021). Chemical-evolution models following the hierarchical structure formation scenario or cosmological simulations implied that RPE stars in the MW are formed in UFD-size galaxies (Ojima et al. 2018; Brauer et al. 2019; Wanajo et al. 2021). However, the mechanism and environments of RPE star formation in the MW are not yet well-understood due to the lack of detailed simulations. In order to understand the formation of RPE stars in the MW, it is necessary to perform high-resolution cosmological zoom-in simulations that can resolve dwarf galaxies with stellar masses down to |${\lesssim}10^5\, \mathrm{M}_{\odot }$| with models of NSMs, metal diffusion, and sufficient numbers of snapshots to follow RPE star formation. Although sophisticated galaxy-formation simulations have been published recently (e.g. Guedes et al. 2011; Sawala et al. 2016; Wetzel et al. 2016; Grand et al. 2017; Buck et al. 2020; Font et al. 2020; Agertz et al. 2021; Applebaum et al. 2021), it is still challenging to perform simulations considering all of these requirements due to the high computational and storage costs.

This study has performed the highest stellar mass and time resolution simulation of an MW-like galaxy with enrichment of r-process elements, which is able to distinguish between r-process enrichment in the UFD-like low-mass building blocks and in situ (local) enrichment. Thanks to these advancements, we can study the kinematics of the simulated stars in order to compare with observed CDTGs.

Our paper is outlined as follows. In Section 2, we summarize and compare previous cosmological simulations of MW-like galaxies with r-process enrichment with our study. In Section 3, we describe our code, adopted models, and initial conditions. Section 4 describes the chemical abundances and kinematics of our simulated stars. Section 5 considers the formation and chemo-dynamical properties of r-II stars. Section 6 summarizes our conclusions and provides some perspectives on future efforts.

2 COSMOLOGICAL SIMULATIONS OF MILKY WAY-LIKE GALAXIES WITH r-PROCESS ENRICHMENTS

Several previous studies have performed cosmological simulations of MW-like galaxies with r-process enrichment (Shen et al. 2015; van de Voort et al. 2015, 2020, 2022; Naiman et al. 2018; Haynes & Kobayashi 2019). Shen et al. (2015) investigated the enrichment of r-process elements in the Eris simulation (Guedes et al. 2011) using a post-processing model for NSMs, and showed that NSMs can explain the level of Eu abundance of metal-poor stars in an MW-like galaxy. van de Voort et al. (2015) computed the enrichment of r-process elements in the FIRE simulation (Hopkins et al. 2014). Their results suggested that large-scale hydrodynamic mixing processes play an essential role in explaining the observed scatter of [r-process/Fe] abundance ratios. Naiman et al. (2018) discussed the distribution of Eu in MW-like galaxies formed in the IllustrisTNG simulations (Springel et al. 2018), and found that the scatter of [Eu/Fe] is sensitive to gas properties in the redshift range z = 2–4.

These studies mainly focused on understanding how to reproduce the observed scatter of [Eu/Fe] in low-metallicity stars. However, it has proven difficult to detail the formation mechanism of stars with r-process elements with these simulations. As described in Section 1, it is important to resolve UFD-size galaxies in cosmological simulations of MW-like galaxies. In addition, as shown in Jeon et al. (2021), it is necessary to resolve r-process enrichment with a time-scale of ∼10 Myr to clarify the formation history of RPE stars. To date, no simulations have had sufficient mass and time resolution to discuss RPE star formation in the context of the formation and assembly history of the MW.

The main aims of this study are to (1) understand the mechanism and environments of RPE star formation and (2) provide a connection between the chemo-dynamical properties of RPE stars and the assembly history of the MW. Here we compute the enrichment of r-process elements in the highest stellar mass and time resolution cosmological zoom-in simulation of an MW-like galaxy, which can resolve dwarf galaxies down to stellar masses of |${\sim}10^5\, \mathrm{M}_{\odot }$|⁠, with a snapshot interval of ∼10 Myr. This simulation, with in total 42 TB outputs, makes it possible to extract information about RPE star formation and chemo-dynamical properties. As noted above, we focus on r-II stars ([Eu/Fe]|$\, \gt \, +0.7$| and [Ba/Eu]|$\, \lt \, 0$|⁠) in the solar neighbourhood, defined in Section 3.3, and regard europium (Eu) as representative of the r-process elements in order to compare with observations.

3 METHOD AND MODELS

3.1 Code

Here, we briefly describe the code employed in this study. Details of the code are described in Saitoh et al. (2008), Hirai et al. (2018), and Hirai, Wanajo & Saitoh (2019). In this study, we adopt the N-body/smoothed particle hydrodynamics (SPH) code asura (Saitoh et al. 2008, 2009). Gravity is computed with a tree method (Barnes & Hut 1986) parallelized following Makino (2004). Hydrodynamics are computed using the density-independent SPH (DISPH) method, which can correctly treat instabilities in contact discontinuity (Saitoh & Makino 2013, 2016). A Fully Asynchronous Split Time-Integrator (FAST) for a Self-Gravitating Fluid scheme is adopted to reduce the computational cost of the self-gravitating fluid systems in SPH (Saitoh & Makino 2010). A time-step limiter, which enforces the time-step difference among local particles small enough to solve the propagation of shocks, is also implemented (Saitoh & Makino 2009). We adopt the metallicity dependent cooling and heating functions from 10 to 10⁹ K generated by cloudy version 13.5 (Ferland et al. 1998, 2013, 2017). Effects of self-shielding (Rahmati et al. 2013) and the ultraviolet background radiation field (Haardt & Madau 2012) are also implemented in our simulations.

Star particles are stochastically created following the Schmidt law (Schmidt 1959). Here, we implement the models for star formation described in Okamoto et al. (2003) and Saitoh et al. (2008). In this model, gas particles need to satisfy conditions for star formation as follows: (1) the number density is higher than 100 cm⁻³, (2) the temperature is lower than 1000 K, (3) the divergence of the velocity is less than zero, and (4) the particle is not heated by supernovae. The threshold density for star formation is similar to the mean density of giant molecular clouds (e.g. Solomon et al. 1987; Heyer et al. 2009).

When a gas particle satisfies all of these conditions, we compare a randomly generated number (⁠|$\mathcal {R}$|⁠) between 0 and 1 with the probability (p_*) in a given time-step (Δt):

$$\begin{eqnarray*} p_{*} = \frac{m_{\mathrm{gas}}}{m_{*}} \left\lbrace 1-\exp \left(-c_{*} \frac{\Delta t}{t_{\mathrm{dyn}}}\right)\right\rbrace , \end{eqnarray*} $$

(1)

where m_gas, m_*, t_dyn, and c_* are the mass of one gas particle, the mass of one star particle, the dynamical time of the star-forming region, and a dimensionless star-formation efficiency. Here, we set c_* = 0.5, motivated by observations of star formation efficiency in star clusters (e.g. Lada & Lada 2003) and previous simulations (e.g. Tasker & Bryan 2006, 2008; Saitoh et al. 2008). If |$\mathcal {R}\, \lt \, p_{*}$|⁠, one-third of the mass of gas particles are converted into star particles, following Okamoto et al. (2003, 2005) and Saitoh et al. (2008), in order to prevent an abrupt decrease of the number of cold gas particles due to star formation.

Star particles are approximated as simple stellar populations (SSPs). We adopt the initial mass function (IMF) of Chabrier (2003) from 0.1 to 100 M_⊙. The metallicity-dependent stellar lifetime table is taken from Portinari, Chiosi & Bressan (1998). Since mega metal-poor (⁠|$\mathrm{[Fe/H]}\, \lt \, -6$|⁠) stars are mainly formed by molecular hydrogen cooling, rather than the fine-structure line cooling by heavy elements (Omukai et al. 2005), Population III star formation simulations predict top-heavy IMFs (e.g. Hirano et al. 2015; Stacy, Bromm & Lee 2016). Therefore, we adopt the top-heavy IMF from 0.6 to 300 M_⊙ estimated from simulations of Population III star formation (Susa, Hasegawa & Tominaga 2014) for stars with |$Z\, \lt \, 10^{-5}\, Z_{\odot }$|⁠. A larger number of low-mass stars are predicted to be formed towards higher metallicity for |$Z\, \gtrsim \, 10^{-5}\, Z_{\odot }$| (Chon, Omukai & Schneider 2021). The number of Lyman α photons is evaluated with pégase (Fioc & Rocca-Volmerange 1997) in order to heat gas around massive stars with ages less than 10 Myr to 10⁴ K as an H ii region. The Chemical Evolution Library (celib) is adopted for the implementation of different nucleosynthetic sources (Saitoh 2016, 2017; Hirai, Fujii & Saitoh 2021). Stars from 13 to 40 M_⊙ for |$Z\, \ge \, 10^{-5}\, Z_{\odot }$| (from 13 to 300 M_⊙ for |$Z\, \lt \, 10^{-5}\, Z_{\odot }$|⁠) are assumed to explode as CCSNe. We adopt the yields from Nomoto, Kobayashi & Tominaga (2013). We also adopt a model of type Ia supernovae (SNe Ia), assuming a power-law delay-time distribution with an index of −1 and a minimum delay time of 40 Myr (e.g. Totani et al. 2008; Maoz & Mannucci 2012; Maoz, Mannucci & Brandt 2012). We adopt the yields of SNe Ia from the N100 model of Seitenzahl et al. (2013). The minimum delay time for SNe Ia is taken from the lifetime of the likely most massive progenitors (8 M_⊙) of white dwarfs. We note that the short end of the delay-time distribution of SNe Ia is highly uncertain (Strolger et al. 2020). Several chemical-evolution studies prefer a minimum delay time of SNe Ia longer than 500 Myr (Homma et al. 2015; Wanajo et al. 2021).

We assume that Eu is synthesized exclusively by NSMs. We adopt a model of NSMs with a delay-time distribution with power-law index −1, with a minimum delay time of 10 Myr, following population-synthesis calculations (e.g. Dominik et al. 2012). Recent studies have shown that the power-law index of −1 in NSMs cannot reproduce the decreasing trend of [Eu/Fe] towards higher metallicity seen in the MW disc (Hotokezaka, Beniamini & Piran 2018; Côté et al. 2019). This problem can be resolved if the delay times of NSMs are sufficiently shorter than those of SNe Ia (Wanajo et al. 2021). In this study, we assume that 0.2 per cent of stars with 8–20 M_⊙ lead to NSMs. This rate corresponds to ∼10 events Myr⁻¹ in MW-equivalent galaxies. We also assume that each NSM produces 1.0 × 10⁻⁴ M_⊙ of Eu to reproduce the observed peak of [Eu/Fe] distribution ([Eu/Fe] ∼+0.4; Holmbeck et al. 2020). This value corresponds to the total ejecta mass of 0.05 M_⊙ of one NSM, which is consistent with the values estimated from observations of GW170817 (e.g. Abbott et al. 2017a; Arcavi et al. 2017; Nicholl et al. 2017; Tanaka et al. 2017).

Supernova-feedback models are implemented following Hopkins et al. (2018a). Thus, we transfer the terminal momentum of supernova remnants to surrounding gas particles when a supernova explodes. We assume that CCSNe and SNe Ia explode with 10⁵¹ erg. We also assume that 5 per cent of stars with 20–40 M_⊙ explode as broad-line type Ic supernovae (hypernovae) with energies of 10⁵² erg. This model does not include the energy feedback from stellar winds. The energy of NSMs is not taken into account, as in Shen et al. (2015) and van de Voort et al. (2015), because the rate is significantly lower than CCSNe, and the estimation of energy from NSMs is highly uncertain (e.g. Hotokezaka et al. 2013). Safarzadeh & Scannapieco (2017) confirmed that the energy of NSMs does not largely affect the abundances of Eu in UFDs.

In this study, we implement the turbulence-induced metal-diffusion model (Shen, Wadsley & Stinson 2010; Hirai & Saitoh 2017; Saitoh 2017). The ith elements (Z_i) injected into the interstellar medium (ISM) are diffused to the surrounding gas particles following the equation below:

$$\begin{eqnarray*} \frac{{\rm {d}}Z_i}{{\rm {d}}t} = \nabla (D\nabla {Z_i}),~~ D = C_{\rm {d}}|S_{ij}|h^2, \end{eqnarray*} $$

(2)

where h is the smoothing length of the SPH particle, S_ij is the trace-free shear tensor, and C_d is a scale factor for metal diffusion. In this study, we adopt C_d = 0.01. This value is calibrated using the observed r-process abundances in dwarf galaxies (Hirai & Saitoh 2017). Their study and Hirai et al. (2018) have shown that models with C_d = 0.01–0.1 can produce the scatter of [Ba/Fe] and [Zn/Fe] consistent with observations.

3.2 Initial conditions

In this study, we performed a cosmological zoom-in simulation of an MW-like galaxy. A low-resolution pre-flight cosmological simulation of structure formation was performed with the N-body code gadget-2 (Springel 2005). The initial condition was generated by music (Hahn & Abel 2011). The box size is (36 h⁻¹ Mpc)³. Density fluctuations on this scale are in the linear regime at z = 0 (Okamoto et al. 2003). We set the cosmological parameters as follows: |$\Omega _{\rm {m}} = 0.308$|⁠, |$\Omega _{\Lambda } = 0.692$|⁠, |$\Omega _{\rm {b}} = 0.0484$|⁠, |$H_{0} = 67.8$| km s⁻¹ Mpc⁻¹, |$\sigma _{8} = 0.815$|⁠, and |$n_{\rm {s}} = 0.968$| (Planck Collaboration XIII 2016).

We have selected the haloes from this simulation that can form MW-like galaxies. The selection criteria is as follows: (1) the total mass is from 5 × 10|$^{11}\, \mathrm{M}_{\odot }$| to 2 × 10|$^{12}\, \mathrm{M}_{\odot }$|⁠, (2) over half of the halo’s total mass is developed by z = 2, (3) spin parameters of the haloes are from 0.02 to 0.07, and (4) there are no haloes larger than 10|$^{13}\, \mathrm{M}_{\odot }$| within 3 h⁻¹ Mpc from the simulated halo (e.g. Ishiyama et al. 2013; Griffen et al. 2016). We use the AMIGA halo finder (ahf) to find haloes and construct merger trees (Gill, Knebe & Gibson 2004; Knollmann & Knebe 2009). We select high-resolution particles following the convex-hull structure to minimize the computational cost. High-resolution particles are collected within three times the virial radius of the main halo at z = 0. The total number of high-resolution dark matter and gas particles is 8.0 × 10⁷. High-resolution zoom-in hydrodynamics simulations were performed with asura (see Section 3.1 for details). We have confirmed that the contamination from low-resolution particles is less than 1 per cent within the virial radius of the main halo (Oñorbe et al. 2014). Table 1 lists the parameters of the model halo adopted in this study.

Table 1.

Parameters of the selected halo. From left to right, columns show the name of the model, the total number of particles (N), the virial mass of the main halo (M_vir), the virial radius (R_vir) at z = 0, the mass of one dark matter particle (m_DM), the initial mass of one gas particle (m_gas), the initial mass of one star particle (m_*) formed from gas particles with a mass of 1.3 × 10|$^{4}\, \mathrm{M}_{\odot }$|⁠, and the gravitational softening length for dark matter (ϵ_DM) and gas particles (ϵ_g).

Model	N	Mvir	Rvir	mDM	mgas	m*	ϵDM	ϵg
		(M⊙)	(kpc)	(M⊙)	(M⊙)	(M⊙)	(pc)	(pc)
220lv12	9.2 × 107	1.2 × 1012	1.5 × 102	7.2 × 104	1.3 × 104	4.5 × 103	85	82

3.3 Definition of solar neighbourhood stars in the simulation

We define the position of the ‘sun’ in the simulation to select star particles in a region comparable to the observations. Following the Sun’s location and velocity in the Galaxy (Bland-Hawthorn & Gerhard 2016), we set the ‘sun’ in the simulation to the centre of the mass of stars located in (1) 5–11 kpc from the Galactic Centre, (2) −50 to 50 pc from the disc plane, and (3) 150–350 km s⁻¹ of the azimuthal velocity. Following these criteria, the position of the ‘sun’ in model 220lv12 is (⁠|$x,\, y,\, z$|⁠) = (4.4 kpc, −5.0 kpc, 0.0 kpc). This study mainly analyses stars within 5 kpc from this position to compare with observations. We also ignore stars within 3 kpc from the Galactic Centre to remove bulge stars, which are not typically targeted in observational surveys. Here, we define star particles satisfying these conditions as solar neighbourhood stars. Although we do not explicitly categorize stars in this study, this region contains stars belonging to the thin, thick discs, and stellar halo.

We note that all of the above assumptions for our simulations were made in advance of any comparison with the observed stellar kinematic and abundance information described below. That is, no ‘tweaking’ of the assumptions has been made, although this might be explored in the future.

4 RESULTS

4.1 Formation of a Milky Way-like galaxy

The simulated galaxy has similar properties to the MW. Fig. 1 shows face-on (Fig. 1a) and edge-on (Fig. 1b) views of the stars and gas in 220lv12, indicating that a disc galaxy is formed in this simulation. The stellar and halo masses within the virial radius at |$z = 0$| are |$8.44\, \times \, 10^{9}\, \mathrm{M}_{\odot }$| and |$1.20\, \times \, 10^{12}\, \mathrm{M}_{\odot }$|⁠, respectively. The stellar mass–halo mass ratio at the time of the peak stellar mass is 0.007. This value is similar to the stellar mass–halo mass ratios of the MW (Bland-Hawthorn & Gerhard 2016) and other cosmological simulations of MW-like galaxies (e.g. Guedes et al. 2011; Grand et al. 2017; Hopkins et al. 2018b; Font et al. 2020; Agertz et al. 2021; Applebaum et al. 2021). Although the final stellar mass is smaller than the MW (≈5 |$\pm \, 1\, \times$| 10|$^{10}\, \mathrm{M}_{\odot }$|⁠), we minimize the effects of the difference in the stellar mass by adjusting the position of the ‘sun’ to avoid selecting outer-disc stars, which are not presently searched for RPE stars. In Section 5.4, we compare our results to previous simulations of MW-like galaxies. We compare the results from this simulation to the observations of pertinent elemental-abundance ratios and the stellar kinematics of RPE stars.

Figure 1.

The (a) face-on and (b) edge-on view of the stellar and gas distribution in the simulated MW-like galaxy at z = 0. The colour scale represents the log-scale mass fraction of stars and gas in each grid of 0.47 × 0.47 kpc².

The metallicity distribution function (MDF) of star particles within 5 kpc from the ‘sun’ and outside of 3 kpc from the Galactic Centre in our simulation shown in Fig. 2 is similar to that of the MW. Fig. 2(a) compares simulated and observed MDF taken from GALAH DR3 for |$\mathrm{[Fe/H]}\, \gt \, -2$| (Buder et al. 2021). Since GALAH data do not obtain accurate metallicity estimates for stars with |$\mathrm{[Fe/H]}\, \lt -$|2, we compare the MDFs for stars with |$\mathrm{[Fe/H]}\, \gt -2$|⁠. The simulated MDF of solar neighbourhood stars for |$\mathrm{[Fe/H]}\, \gt \, -2$| exhibits a median [Fe/H] = −0.17, with an interquartile range of 0.44 dex. By comparison, the MDF based on GALAH DR3 exhibits a median [Fe/H] = −0.16 and interquartile range of 0.35 dex (Buder et al. 2021).

Figure 2.

MDFs of (a) the MW’s solar neighbourhood stars and (b) halo stars. Blue and orange histograms represent simulated and observed stars, respectively. The simulated data are the MDF of star particles with distances within 5 kpc from the ‘sun’ and outside of 3 kpc from the Galactic Centre in the simulation. Observed data in panel (a) are taken from GALAH DR3 (Buder et al. 2021) with |$\tt {flag\_sp}=0$|⁠, |$\tt {snr\_c3\_iraf}\gt 50$|⁠, |$\tt {flag\_fe\_h}=0$|⁠. Observed data in panel (b) are taken from the Galactic halo stars of Youakim et al. (2020). See text for details.

The lower metallicity end of simulated MDFs is also similar to observations. Fig. 2(b) compares simulated and observed halo star MDFs (Youakim et al. 2020) for stars with |$\mathrm{[Fe/H]}\, \lt \, -2$|⁠. Although a two-sample Kolmogorov–Smirnov (KS) test for the observed and simulated MDFs rejects the null hypothesis that the two distributions are drawn from the same parent populations at a significance level of 0.05 (p = 0.04), both data exhibit decreasing trends towards lower metallicity. Peaks in the simulated MDFs at [Fe/H] = −2.0 and −2.6 are possibly due to the contribution of the relatively massive accreted components (e.g. haloes 3, 4, and 5 in Table 2).

Table 2.

List of haloes where r-II stars are formed. From left to right, the columns show the halo number, the time when the halo formed (t_form), the time when the first r-II star formed in the halo (t_rII), the time when the halo is accreted into the main halo (t_acc), and the total mass (M_tot), gas mass (M_gas), and stellar mass (M_*) at the time of the first r-II star formation in a given halo. The time listed in this table is the one from the beginning of the simulation. Halo 0 is the main halo (in situ component) of this simulation. The halo number is assigned based on the order of the snapshots in which the first r-II star is formed.

Halo No.	tform	trII	tacc	Mtot	Mgas	M*
	(Gyr)	(Gyr)	(Gyr)	(⁠|$10^9\, \mathrm{M}_{\odot }$|⁠)	(⁠|$10^8\, \mathrm{M}_{\odot }$|⁠)	(⁠|$10^6\, \mathrm{M}_{\odot }$|⁠)
0	0.24	0.43	–	5.02	3.60	7.34
1	0.17	0.40	0.73	0.36	0.47	0.05
2	0.18	0.39	0.49	0.38	0.36	0.12
3	0.18	0.44	0.76	6.49	4.66	4.28
4	0.15	0.45	1.27	2.91	1.88	2.26
5	0.13	0.53	1.12	3.65	3.31	0.96
6	0.36	0.53	0.57	1.29	1.14	0.37
7	0.15	0.55	0.68	0.30	0.21	0.24
8	0.13	0.58	2.46	2.55	2.87	4.03
9	0.18	0.58	8.44	2.85	1.68	1.69
10	0.13	0.59	2.48	2.53	2.85	1.64
11	0.39	0.59	0.94	1.04	2.11	0.45
12	0.44	0.71	1.69	3.57	0.68	0.50
13	0.39	0.93	3.43	3.08	1.92	6.12
14	0.47	1.07	4.57	4.24	6.49	1.75
15	0.33	1.23	1.64	1.37	0.73	0.50
16	0.38	1.67	9.78	33.99	20.27	14.38

Our simulation also exhibits general trends of [Mg/Fe] ratios that are similar to the observations. Fig. 3 shows the [Mg/Fe] ratios, as a function of [Fe/H], for solar neighbourhood stars simulated in 220lv12. For [Fe/H] < −2, the [Mg/Fe] ratios are constant (∼+0.5) with a small scatter (⁠|$\sigma = 0.19$| dex). Most of these stars are associated with the stellar halo. As is expected, the iron production from SNe Ia decreases the [Mg/Fe] ratios for [Fe/H] > −2. Most stars with |$\mathrm{[Mg/Fe]}\, \lt \,$|0 and |$-2\, \lt \, \mathrm{[Fe/H]}\, \lt -1$| come from the accreted component, halo 5, discussed in Section 5. Note that our assumption of the minimum time-scale for SNe Ia (40 Myr) would shift the metallicity of the ‘knee’ in [Mg/Fe]. This assumption is likely responsible for the constant [Eu/Fe] ratios at high metallicity due to small differences in the delay time distribution of SNe Ia and NSMs. This issue is discussed further in Section 5.5.

Figure 3.

[Mg/Fe] versus [Fe/H] at z = 0 in the simulation. Stars in the simulation are within 5 kpc distance from the ‘sun’ and outside of 3 kpc from the Galactic Centre. The colour scale represents the log-scale mass fraction of stars in each grid of 0.02 × 0.01 dex², from 10⁻⁶ (purple) to 10⁻⁴ (yellow). The orange represent observed abundances in the MW taken from the SAGA data base (Suda et al. 2008, 2011, 2017; Yamada et al. 2013) for |$\mathrm{[Fe/H]}\, \lt -1.5$| and GALAH DR3 (Buder et al. 2021) with |$\tt {flag\_sp}=0$|⁠, |$\tt {snr\_c3\_iraf}\gt 50$|⁠, log g < 1.9, |$\tt {flag\_fe\_h}=0$|⁠, and |$\tt {flag\_Mg\_fe}=0$|⁠. Note that RPE stars from Shank et al. (2022a) are not included in this figure.

4.2 The enrichment of Eu

Fig. 4 shows [Eu/Fe], as a function of [Fe/H], for solar neighbourhood star particles. As shown in this figure, there exists a star-to-star scatter in the [Eu/Fe] ratios of over 2 dex, similar to the observations for [Fe/H] ≲ −2. Here, we define stars with [Eu/Fe] > +0.7 as r-II stars. The fraction of r-II stars is 8.7 per cent in this simulation, which is similar to the fraction reported in the RPA (8.3 per cent; Holmbeck et al. 2020) – note that both of these values are computed for stars with |$\mathrm{[Fe/H]}\, \lt -2.5$| and [Eu/Fe]|$\, \gt -0.8$|⁠. SNe Ia dominantly start to contribute for |$\mathrm{[Fe/H]}\, \gt \, -2$| in this simulation. Owing to the delay times adopted for NSMs, the RPE stars in the simulation are formed with [Fe/H] ≳ −3.5.

Figure 4.

[Eu/Fe] as functions of [Fe/H] at z = 0. Stars in the simulation are taken within 5 kpc from the ‘sun’ and outside of 3 kpc from the Galactic Centre. The colour scale represents the log-scale mass fraction of stars in each grid of 0.02 × 0.01 dex² from 10⁻⁶ (purple) to 10⁻² (yellow). Orange plots represent observed abundances for stars in the MW taken from RPA (Holmbeck et al. 2020), GALAH DR3 (Buder et al. 2021) with |$\tt {flag\_sp}=0$|⁠, |$\tt {snr\_c3\_iraf}\gt 50$|⁠, log g < 1.9, |$\tt {flag\_fe\_h}=0$|⁠, |$\tt {flag\_Eu\_fe}=0$| and Shank et al. (2022a) for stars with [Fe/H] < −0.8 and [Eu/Fe] >+0.3. Shank et al. (2022a) are composed of the literature data taken from JINAbase (Abohalima & Frebel 2018) and GALAH DR3. The upper left region indicated by the red dashed line is plotted in Fig. 5.

The simulated distribution of [Eu/Fe] for metal-poor r-II stars provides a good match to that of the observations for r-II stars, as shown in Fig. 5 for r-II stars with |$\mathrm{[Fe/H]}\, \lt -0.8$|⁠. The metallicity cut (⁠|$\mathrm{[Fe/H]}\, \lt -0.8$|⁠) is based on Shank et al. (2022a) to be the same as the data taken by the RPA (Holmbeck et al. 2020). As can be appreciated from inspection, the simulated [Eu/Fe] distribution of metal-poor r-II stars reproduces the decreasing number fraction of these stars towards higher [Eu/Fe] ratios. The skewness of the [Eu/Fe] distribution for the simulation and the observations are 1.9 and 2.2, respectively, indicating that the two distributions are quite similar. From the application of a two-sample KS test for the observed and simulated [Eu/Fe] distributions of r-II stars, we are unable to reject the null hypothesis that the two distributions are drawn from the same parent populations (p = 0.32).

Figure 5.

The [Eu/Fe] distribution of the observed (orange) and simulated (blue) r-II stars with |$\mathrm{[Fe/H]}\, \lt -0.8$|⁠. The simulated data are the solar neighbourhood stars defined in Section 3.3. Observed data are taken from Shank et al. (2022a) and GALAH DR3 (Buder et al. 2021).

4.3 Kinematics of r-process-enhanced stars

Few r-II stars with [Fe/H] < −1 exhibit disc-like orbits. Fig. 6 is the Toomre diagram for stars around the ‘sun’ in our simulation. In the figure, we plot the solar neighbourhood disc and r-II stars with [Fe/H] < −1 within 0.2–5 kpc from the ‘sun’ in the simulation, respectively, to compare with the observations (Shank et al. 2022a). In the simulation, 92 per cent of stars lie beyond a 100 km s⁻¹ velocity difference from the ‘sun’. This result indicates that most RPE stars are kinematically associated with the stellar halo system. The fraction of simulated r-II stars with prograde and retrograde orbits are 55 per cent and 45 per cent, respectively, similar to the observations. Shank et al. (2022a) find that the fractions of prograde and retrograde r-II stars are 65 per cent and 35 per cent, respectively. They also find that only 8 per cent of r-II stars exhibit disc-like orbits.

Figure 6.

Toomre diagram for the observations (orange filled circles), simulated stars with a distance of 0.2 kpc from the ‘sun’ (green squares), and simulated r-II stars with [Fe/H] < −1 (blue filled circles). The observed data are taken from Shank et al. (2022a). The velocities are plotted with a cylindrical coordinate system (V_R, V_ϕ, V_z), where V_ϕ is the azimuthal velocity and V_⊥ is defined as (V_R + V_z)^1/2.

For [Fe/H] > −1, all of the simulated r-II stars are formed in situ. In this metallicity range, 62 per cent of the simulated r-II stars lie within a 100 km s⁻¹ velocity difference from the ‘sun’, i.e. there are a larger fraction of stars with disc-like kinematics compared to stars with [Fe/H] < −1. Moreover, the fraction of stars with prograde orbits is larger than for the lower metallicity stars. For [Fe/H] > −1, 84 per cent of the stars exhibit prograde orbits. If this were the case, a greater fraction of r-II stars with prograde orbits would be observed at high metallicity. As we discuss in Section 5.5, the fraction of r-II stars could be overestimated in this simulation due to the assumed level of metal mixing and the delay-time distribution of NSMs. Previous observations suggest that relatively few r-II stars are found in this metallicity range. Therefore, we do not include these stars in Fig. 6.

5 DISCUSSION

Here, we discuss the formation of r-II stars and their chemo-dynamical properties.

Table 2 shows the list of the haloes that formed r-II stars and were accreted into the main halo (halo 0) in our simulation. These are ordered by the snapshot in time when the first r-II star is detected in each.

5.1 The formation of r-II stars

In our simulation, most of the r-II stars are formed in the early epochs of galaxy evolution. Fig. 7 shows [Fe/H] as a function of the time from the beginning of the simulation. The discontinuous distribution of stars seen around 2–6 Gyr is due to the suppression of star formation by supernova feedback. Star formation in galaxies is also affected by a cold accretion flow (e.g. Kereš et al. 2005). As seen in this figure, over 90 per cent of the r-II stars are formed within 4 Gyr from the beginning of the simulation (⁠|$z\, \gt \, 1.6$|⁠). All of the stars with [Fe/H] < −1.5 are formed at |$z\, \gt \, 1.6$|⁠. This epoch corresponds to the time when the stellar halo is formed, suggesting that most of these stars come from dwarf galaxies that are later accreted to the main halo. Since the Galactic disc is formed in the later stage of the galaxy evolution, this result explains the halo-like kinematics seen in Fig. 6.

Figure 7.

[Fe/H], as a function of the formation time, of the solar neighbourhood simulated stars. The colour scale represents the log-scale mass fraction of stars in each grid of 0.05 Gyr × 0.02 dex, from 10⁻⁶ (purple) to 10⁻³ (yellow). The red dots denote the r-II stars formed in the simulation.

The r-II stars with |$\mathrm{[Fe/H]}\, \ge -1.5$| are formed throughout the entire history of galaxy evolution. Most of these stars are formed in the main progenitor halo or in the disc, where they were locally enhanced in r-process elements. Since NSMs eject prolific amounts of r-process elements, some RPE stars are still formed in the later phases of galaxy evolution.

The low-metallicity r-II stars tend to be formed in accreted components, as shown in Fig. 8(a). The number fraction of stars formed in accreted components, as a function of [Fe/H], demonstrates that there is an increasing trend of accreted fractions towards lower metallicity. Over 80 per cent of the r-II stars with [Fe/H] < −2 are formed in accreted components, while all r-II stars with [Fe/H] > −1 are formed in situ.

Figure 8.

The number fraction of stars from accreted components, as a function of (a) [Fe/H] (blue solid curve: stars with [Eu/Fe] > +0.7; orange dashed curve: stars with [Eu/Fe] ≤ +0.7) and (b) [Eu/Fe] (blue solid curve: all r-II stars, magenta dashed curve: r-II stars with [Fe/H] < −1).

The accreted fraction of r-II stars is larger than that of non-RPE stars in the metallicity range |$-2.5\, \lt \, \mathrm{[Fe/H]}\, \lt -1.5$|⁠. At [Fe/H] = −2.2, 96 per cent (71 per cent) of r-II (non-RPE) stars are formed in the accreted components. This result is consistent with van de Voort et al. (2020), who found that 91 per cent of r-II stars with [Fe/H] < −2 were formed in accreted components, using data from the Auriga simulation (Grand et al. 2017). This fraction was larger than for all very metal-poor stars (78 per cent).

Fig. 8(b) shows the accreted fractions as a function of [Eu/Fe]. From inspection, there exists an increasing trend of the accreted fraction towards higher [Eu/Fe] ratios for all r-II stars. This trend is mainly due to the contribution of in situ r-II stars with [Fe/H] > −1. These stars decrease the accreted fraction of r-II stars with [Fe/H] < −1.5. Excluding these stars increases the accreted fraction (the magenta dashed curve in Fig. 8b).

The r-II stars are formed in gas clumps enhanced with r-process elements. Fig. 9(a) shows the time evolution of the average gas-phase [Eu/Fe]. The time variation of [Eu/Fe] depends on the gas mass of the building blocks. In the case of small building blocks (e.g. haloes 1 and 2), the mean [Eu/Fe] largely varies owing to the frequency of NSM events. The maximum mean [Eu/Fe] in halo 1 is +2.11 at the time from the beginning of the simulation (⁠|$t = 0.29\,$|Gyr). At the time of r-II star formation (t = 0.40 Gyr), it is still [Eu/Fe] = +0.73. At t = 0.29 Gyr, the gas mass of this galaxy is 4.7|$\, \times \, 10^6\, \mathrm{M}_{\odot }$| (Fig. 10a). In such a small system, one NSM can enhance most of the gas contained in the halo; stars formed in this phase may have large enhancements of r-process elements as a result.

Figure 9.

Time evolution of the mean gas-phase [Eu/Fe] ratios. Panel (a) represents halo 1 (orange dashed curve) and halo 2 (green dotted curve), while panel (b) shows halo 0 (main progenitor halo, blue solid curve) and halo 3 (red dash–dotted curve). Arrows indicate the formation time of r-II stars in each component (haloes 0: blue, 1: orange, 2: green, 3: red).

Figure 10.

Similar to Fig. 9 but for gas mass, as a function of time, in each halo.

This mechanism is similar to r-II star formation in UFDs. Cosmological zoom-in simulations have shown that the small gas mass of UFDs is sufficient to enhance the r-process abundance after an NSM to form r-II stars (Safarzadeh & Scannapieco 2017; Tarumi et al. 2020), due to the lack of significant dilution by the gas reservoir. Jeon et al. (2021) have shown that the most important parameter to form r-II stars is the time-scale of r-II star formation after the Eu enrichment. Their results suggested that r-II stars need to be formed 10–100 Myr after an NSM event to keep the gas-phase Eu abundance ratio over [Eu/Fe] > +1. The condition to form r-II stars in our simulation is similar. As shown in the orange dashed curve in Fig. 9(a), the average [Eu/Fe] drops within |${\approx}$|100 Myr because of SN feedback, metal diffusion, and newly accreted gas.

The r-II stars found in relatively large building blocks and the main progenitor halo are formed in gas clumps that are locally enhanced in r-process elements. The red dash–dotted curve in Fig. 9(b) shows the mean gas-phase [Eu/Fe] ratio in halo 3. The gas mass of this galaxy reaches |$4.7\, \times \, 10^{8}\, \mathrm{M}_{\odot }$| at the time of the formation of r-II stars (Fig. 10b). At the time of the r-II star formation, the total mass of this halo is larger than the main progenitor halo (halo 0, Table 2). Most of the mass of this halo is stripped when it accretes to halo 0. In this halo, the average gas-phase [Eu/Fe] is not highly enhanced when r-II stars are formed. This result suggests that an NSM cannot solely enhance r-process elements in haloes with a gas mass larger than |${\sim }10^{8}\, \mathrm{M}_{\odot }$|⁠.

Fig. 11 shows the gas-phase [Eu/Fe], as a function of [Fe/H], within 200 pc from a progenitor gas particle of in situ r-II stars. We point to the most r-process enhanced in situ star ([Fe/H] = −1.75, [Eu/Fe] = +1.76) from our sample for this analysis. Before star formation initiates, all gas particles have [Eu/Fe] < 0 (t = 623 Myr). After that, some gas particles show [Eu/Fe] > +1 in this region owing to a nearby NSM (t = 637 Myr); the [Fe/H] ratios are increased due to a nearby supernova. The number of gas particles is increased compared to the previous snapshot because the gas particles are in a converging flow that enters this region. After star formation initiates, several gas particles are consumed to form stars, and the [Eu/Fe] ratios are decreased due to metal diffusion.

Figure 11.

Evolution of gas-phase [Eu/Fe], as a function of [Fe/H], within 200 pc from the formation site of an in situ r-II star. The blue circle, orange triangles, and green squares represent gas particles in this region at t = 623, 637, and 651 Myr, respectively.

One NSM can cause local enhancements of r-process elements in gas clumps. In this simulation, an NSM typically distributes r-process elements to gas clumps with gas mass |${\sim}10^6\, \mathrm{M}_{\odot }$|⁠. Inside these clumps, there are variations of [Eu/Fe] ratios (Fig. 11). Gas particles closer to the NSM are more enhanced in r-process elements than more distant particles. This coincidence is highly rare. Therefore, this mechanism of forming r-II stars is only possible in relatively massive building blocks that can host sufficient numbers of NSMs. This result implies that if stars are formed in a highly RPE region near an NSM, high metallicity r-II stars could be formed.

The formation of r-II stars in relatively massive building blocks found in this study could be used to understand the presence of r-II stars found in massive satellite galaxies (e.g. LMC, SMC, Fornax). Reggiani et al. (2021) found that six out of nine (two out of four) metal-poor LMC (SMC) giants are classified as r-II stars. In the Fornax dSph, three r-II stars with +1.25 ≤ [Eu/Fe] ≤ +1.45 have been confirmed (Reichert et al. 2021). Although these stars have relatively high metallicity (⁠|$-1.3\, \le$| [Fe/H] ≤−0.8), their chemistry clearly shows the typical r-process abundance pattern. These r-II stars could be formed in locally RPE gas clumps.

Lower mass building blocks tend to have larger fractions of r-II stars. Fig. 12 shows the number fraction of r-II stars in building-block galaxies. As seen in this figure, there are clear increasing trends of this fraction towards the lower stellar mass. Since lower mass galaxies tend to have lower gas mass, an NSM enhances the abundances of r-process elements significantly due to the lack of dilution (Fig. 9). Subsequent star formation from r-process-depleted gas decreases the fraction of r-II stars in each halo. The scatter is mainly caused by a newly contributing NSM in each halo. The r-II stars formed around the ejecta of the NSM enhance the fraction of r-II stars and stellar mass of the building block galaxy.

Figure 12.

The number fraction of r-II stars, as a function of the stellar mass of building-block galaxies. The same colour and symbol indicate the building-block galaxy in different snapshots from 0.13 to 1.06 Gyr from the beginning of the simulation; the snapshots are separated by 13.8 Myr. The circles represent haloes from 1 to 10, and the inverted triangles represent haloes from 11 to 12. See legend for details.

The above results could explain the observed trend of r-process enhancement in dwarf galaxies. In the Reticulum II UFD, seven out of nine stars observed in this galaxy are enhanced in r-process elements (Ji et al. 2016a). In contrast, classical dSphs tend to have smaller fractions of RPE stars than in Reticulum II (e.g. Frebel & Norris 2015). There are 41 out of 165 r-II stars in Fornax (Letarte et al. 2010), Ursa Minor (Shetrone et al. 2001; Sadakane et al. 2004), Carina (Shetrone et al. 2003; Norris et al. 2017), Draco (Cohen & Huang 2009) and Sculptor (Shetrone et al. 2003; Geisler et al. 2005) dSphs (counted from the SAGA data base; Suda et al. 2008, 2011, 2017; Yamada et al. 2013). These results imply that r-II stars in classical dSphs are formed by local inhomogeneities, while those in UFDs are formed in a halo entirely enhanced in r-process elements. Another possibility is that r-II stars formed in UFD-mass galaxies that are accreted by other dwarf galaxies. This could happen in relatively large dwarf galaxies, such as LMC and Fornax (e.g. Reichert et al. 2021). Detailed studies of the spatial distribution and kinematics of r-II stars in surviving dwarf galaxies could help resolve this issue.

5.2 Chemo-dynamical properties of r-II stars

In order to extract information on the birth environments of r-II stars from the observations, it is necessary to clarify the chemo-dynamical properties of r-II stars. Integral of motions such as energy and angular momentum are conserved on long time-scales and have been used to identify accreted components (e.g. Helmi 2020). Fig. 13 shows a Lindblad Diagram (the total specific energy, E, as a function of the orbital angular momentum, L_z) for r-II stars. The total specific energy is the sum of the kinetic and potential energy per stellar mass. The kinetic energy is computed based on the velocity of the particle, and the potential energy is calculated using the distance and mass of all particles involved in the force calculation. The orbital angular momentum is defined to be L_z = RV_ϕ, where |$R=\sqrt{x^2+y^2}$|⁠. Stars with low E (⁠|$E\, \lt -4.5\, \times \, 10^5\,$|km² s⁻²) and high L_z (⁠|$L_{\rm {z}}\, \gt \, 0.7\times \, 10^3\,$|kpc s⁻²) are newly formed stars with ages ≲ 50 Myr around the solar neighbourhood. Since these stars are not observed in the MW and are not generally r-process enhanced, we exclude these stars from the discussion.

Figure 13.

Distribution of star particles in the Lindblad diagram; total energy (E) versus the vertical angular momentum (L_z) for (a) in situ r-II stars, (b) stars formed in haloes 1–5, (c) 6–10, and (d) 11–16. The grey-scale represents the mass fraction of all stars within 5 kpc from the ‘sun’ and outside of 3 kpc from the Galactic Centre. The coloured symbols show r-II stars from the same haloes. The halo numbers are indicated in the legend.

As shown in Fig. 13(a), most in situ r-II stars exhibit disc-like kinematics. The solar neighbourhood disc stars reside in the parabola of the E versus L_z diagram with |$L_{\rm {z}}\, \gt \, 0$|⁠. These stars are formed in the Galactic disc’s locally RPE regions. Shank et al. (2022a) find 30 r-II stars with disc-like kinematics in the RPA sample; such stars are candidates for stars formed by local inhomogeneities of r-process enhancement in the Galactic disc.

On the other hand, a few in situ stars have retrograde orbits, high energy, or both. Most of these stars are formed in the primordial main progenitor halo before the disc formation. The r-II stars with the lowest binding energy (⁠|$E\, \lt -4.2\, \times \, 10^5$| km² s⁻²) are also old (ages with ≳11 Gyr). These stars are associated with the main progenitor halo. The kinematics of these r-II stars are similar to the in situ halo stars, as defined in Carollo & Chiba (2021). These authors pointed out that in situ halo stars are located in the lowest binding energy or along a parabola in the E versus L_z diagram. These results suggest that r-II stars with the lowest binding energy or located along the parabola of the Lindblad diagram could be formed in the main progenitor halo.

Figs 13(b)–(d) show E versus L_z for r-II stars formed in accreted components. As expected from Helmi & de Zeeuw (2000), stars formed in the same halo tend to reside in a similar region in their integrals of motion space. Roederer et al. (2018) found eight groups for r-II stars in the phase space by applying the clustering methods. Shank et al. (2022a) identify 38 CDTGs among ∼1800 RPE stars. The results of these studies indicate that r-II stars clustered in dynamical phase space could come from accreted components. According to Fig. 13, some groups of r-II stars from the same halo are located in similar regions in the E versus L_z diagram (e.g. haloes 8, 10, and 11). For example, the standard deviation around the cluster central value (⁠|$\sigma _{\rm {d}} = \sqrt{\sigma _{Lz}^2+\sigma _{E}^2}$|⁠, where σ_Lz and σ_E are standard deviations of L_z and E, respectively) of halo 11 is |$\sigma _{\rm {d}}\, =0.23$| while that of halo 0 is |$\sigma _{\rm {d}} = 0.57$|⁠. These results support the scenario that r-II stars clustered in their dynamical phase space come from accreted components.

The r-II stars from more massive haloes tend to be broadly distributed in the dynamical phase space. For example, haloes 4 and 5 have the peak stellar mass of 9.2|$\, \times \, 10^7\, \mathrm{M}_{\odot }$| and 9.3|$\, \times \, 10^7\, \mathrm{M}_{\odot }$|⁠, respectively. Stars in these haloes exist in |$-1.0\, \times \, 10^3$| kpc km s|$^{-1}\, \lesssim L_{\rm {z}}\, \lesssim 2.0\, \times \, 10^3$| kpc km s⁻¹ and |$-4.6\, \times \, 10^5$| km² s|$^{-2}\, \lesssim E\lesssim -4.0\, \times \, 10^5$| km² s⁻² (Fig. 13b). In contrast, halo 11 has a peak stellar mass of 2.4|$\, \times \, 10^6\, \mathrm{M}_{\odot }$|⁠. This halo is highly clustered in the phase space (Fig. 13d), implying that r-II stars that are tightly clustered could originate from low-mass building blocks.

The r-II stars formed in haloes that were accreted earlier tend to have low energy. For instance, haloes 6 and 11, which reside in the lowest energy in the phase space, are examples of haloes that were accreted early in the simulation. From Table 2, these haloes are accreted into the main halo at t_acc = 0.57 Gyr (halo 6) and 0.94 Gyr (halo 11). In contrast, stars from haloes 9 and 16 have high energy and much later accretion times (t_acc = 8.44 and 9.78 Gyr, respectively). Fig. 14 shows [Mg/Fe], as a function of [Fe/H], grouped by r-II stars from the same haloes. As shown in this figure, early accreted haloes (haloes 6 and 11) are characterized by low metallicity (⁠|$\mathrm{[Fe/H]}\, \lt -2$|⁠) and high [Mg/Fe] ratios (⁠|$\mathrm{[Mg/Fe]}\, \gt \, 0$|⁠). In contrast, later accreted haloes (haloes 9 and 16) exhibit higher metallicity (⁠|$\mathrm{[Fe/H]}\, \gt -2$|⁠) and lower [Mg/Fe] ratios (⁠|$\mathrm{[Mg/Fe]}\, \sim 0$|⁠) compared to earlier accreted haloes, indicating that the latter haloes experienced contributions of iron from SNe Ia.

Figure 14.

[Mg/Fe], as functions of [Fe/H], for (a) in situ r-II stars, (b) stars formed in haloes 1–5, (c) 6–10, and (d) 11–16. Colour plots show r-II stars from the same haloes. The halo numbers are indicated in the legend.

Some groups of r-II stars exhibit a decreasing trend of [Mg/Fe] ratios, typically seen in dwarf galaxies and accreted components. For in situ stars, there is an overall decreasing trend of [Mg/Fe] for |$\mathrm{[Fe/H]}\, \gt -2$| (Fig. 14a). Note that the trend is unclear because we adopt the short minimum delay time for SNe Ia (t_min = 40 Myr; Section 3.1). For accreted stars, some haloes (e.g. haloes 4 and 5) clearly show decreasing trends of [Mg/Fe] at |$\mathrm{[Fe/H]}\, \sim -2$|⁠. On the other hand, relatively low-mass haloes (e.g. haloes 1, 2, and 11) do not exhibit decreasing trends. These haloes halted star formation before SNe Ia started contributing to the Fe enrichment. SNe Ia begin to contribute to the Fe enrichment after 40 Myr from the first star formation in each halo in this simulation.

The observed r-II stars generally have similar [Mg/Fe] ratios, as found in this study. Xing et al. (2019) describe the star LAMOST J112456.61+453531.3 with [Fe/H] = −1.27, [Eu/Fe] = +1.1, and [Mg/Fe] = −0.31 (see also Sakari et al. 2019). They argued that this star could come from a disrupted dwarf galaxy. In our simulation, there is a star with [Fe/H] = −1.49, [Eu/Fe] = +1.27, and [Mg/Fe] = −0.35 formed in halo 5. This star is formed at t = 6.64 Gyr when the stellar mass of halo 5 is 4.10|$\, \times \, 10^7\, \mathrm{M}_{\odot }$|⁠. We also find a similar in situ halo star. This star has [Fe/H] = −1.26, [Eu/Fe] = +1.75, and [Mg/Fe] = −0.46, formed at t = 6.68 Gyr when the stellar mass of the halo is 9.20|$\, \times \, 10^7\, \mathrm{M}_{\odot }$|⁠. These results suggest that r-II stars with high metallicity and low [Mg/Fe] ratios can arise from the RPE gas clumps formed in relatively massive building blocks.

Many of the r-II stars with |$\mathrm{[Fe/H]}\, \lt -2$| exhibit |$\mathrm{[Mg/Fe]}\, \gt 0$|⁠. According to Fig. 14, most r-II stars have [Mg/Fe] > 0. For example, Shank et al. (2022a) find that most r-II stars in the RPA sample have |$\mathrm{[Mg/Fe]}\, \gt 0$|⁠. These results suggest that r-II stars with |$\mathrm{[Fe/H]}\, \lt -2$| are mainly formed in haloes without significant contributions from SNe Ia, i.e. those in an early stage of their chemical evolution.

The r-II stars in each accreted component exhibit low dispersions in [Fe/H] and larger variations of their [Eu/Fe] ratios. For example, Fig. 15 depicts [Eu/Fe], as a function of [Fe/H], for in situ r-II stars and those in each halo. This figure shows that r-II stars in the same accreted component have low scatter in [Fe/H], while in situ r-II stars are seen over a wide metallicity range. All of the accreted components with more than three stars have standard deviations of [Fe/H] less than 0.5 dex. This result is consistent with observed r-II stars in CDTGs (Roederer et al. 2018; Gudin et al. 2021). Shank et al. (2022a) show that the standard deviations of [Fe/H] in CDTGs with more than five stars are less than 0.33 dex, and those in [Eu/Fe] are less than 0.16 dex. Note that this estimate includes r-I stars. Scatters of [Eu/Fe] in accreted components, including r-I stars, will be examined in our forthcoming paper. The short time-scale for r-II star formation results in the low scatter in [Fe/H]. The longest time separation of the first and the last r-II star formation is 620 Myr in halo 12. Most haloes form r-II stars within 100 Myr of one another.

Figure 15.

Same as Fig. 14 but for [Eu/Fe] as a function of [Fe/H].

The dispersion of [Eu/Fe] ratios seen in accreted components could explain the r-II star found in the Indus stellar stream. Hansen et al. (2021) reported that a star Indus_13 shows [Eu/Fe] = +1.81, while other Indus member stars have [Eu/Fe] ≈ +0.5. Halo 1 in our simulation shows a similar trend. In this halo, the oldest and most RPE star has [Eu/Fe] = +1.82, but other stars have +0.71|$\, \le \,$|[Eu/Fe]|$\, \le \,$| + 0.86. The lack of stars between these values is because we only select solar neighbourhood stars. The scatter seen in accreted components can be caused by these haloes’ inhomogeneous [Eu/Fe] distribution, as even if their average [Eu/Fe] is enhanced, there are large scatters of the gas-phase [Eu/Fe] ratios from [Eu/Fe]|$\, \lt -1$| to [Eu/Fe]|$\, \gt \, +2$|⁠.

5.3 Requisite mass resolution

Here, we discuss the resolution required to resolve the scatter of [Eu/Fe] abundance ratios caused by the building blocks and in the in situ component. In order to accurately compute the chemo-dynamical properties of the building blocks, the dark matter, gas, and stellar components need to be resolved. The typical mass and the virial radius of dark matter haloes that form UFD-mass galaxies is |${\sim}10^9\, \mathrm{M}_{\odot }$| and ∼10 kpc, respectively (e.g. Wheeler et al. 2019). In such haloes, most stars are formed within ∼1 kpc. This means that it is necessary to resolve the inner structures of haloes down to at least 10 per cent of the virial radius. According to Power et al. (2003), ∼10³ particles need to be contained within 10 per cent of the virial radius to estimate the density profile of haloes robustly (see their fig. 14). In this case, ∼10⁴ dark matter particles should be contained within the virial radius. Therefore, simulations with |$m_{\rm {DM}}\lesssim 10^5\, \mathrm{M}_{\odot }$| in a zoomed-in region can resolve haloes that are able to form UFDs.

Gas clumps of |${\sim}10^6\, \mathrm{M}_{\odot }$| need to be resolved to discuss the formation of RPE stars. Ji et al. (2016a) estimated that an NSM could explain the value of [Eu/H] seen in Reticulum II if the ejecta from the NSM is diluted by a gas cloud with |${\sim}10^6\, \mathrm{M}_{\odot }$|⁠. SPH simulations evaluate physical quantities using ∼100 nearest neighbour particles, thus the smallest gas-mass scale for which one can evaluate physical quantities is ∼100 m_gas. Therefore, simulations with |$m_{\rm {gas}}\, \lesssim \, 10^4\, \mathrm{M}_{\odot }$| can resolve the formation of r-II stars in both the accreted and in situ components.

For the stellar components, recent very high-resolution simulations of MW-like galaxies have shown that the size and velocity dispersions of satellite galaxies with more than 10 star particles successfully reproduce the size and the velocity dispersion of the observed UFDs (Applebaum et al. 2021; Grand et al. 2021). This result indicates that we can obtain reliable stellar properties for galaxies with |${\sim}10^5\, \mathrm{M}_{\odot }$| in simulations with star particle masses of |$m_{*}\lesssim \, 10^4\, \mathrm{M}_{\odot }$|⁠. As shown in Table 1, all of these requirements are satisfied by our simulation. We note that these are minimum requirements, and further finer resolutions improve analyses. Thanks to this high resolution, our simulation can resolve the scatter of [Eu/Fe] in each RPE star-forming region and compare the kinematics of r-II stars with observed CDTGs.

5.4 Comparison to previous studies

As discussed in Section 2, some previous studies have considered the enrichment of r-process elements in their simulations of MW-like galaxies (Shen et al. 2015; van de Voort et al. 2015, 2020, 2022; Naiman et al. 2018; Haynes & Kobayashi 2019). Compared with these, our simulation has the highest stellar-mass resolution. Table 3 compares the mass-resolution and metal-mixing models of previous simulations. The mass resolution affects the star-to-star scatter of [Eu/Fe] at low metallicity. All simulations with |$m_{\rm {gas}}\, \lesssim \, 10^4\, \mathrm{M}_{\odot }$| produce such scatter. van de Voort et al. (2015) have shown that their low-resolution model underestimates the observed scatter. Haynes & Kobayashi (2019) could not reproduce the [Eu/Fe] scatter by their NSM models. Since the number of stars with |$\mathrm{[Fe/H]}\, \lt \, -$|2.5 is much lower than for more metal-rich stars, and they tend to form in accreted components, it is difficult to correctly follow the formation of RPE stars in simulations with mass resolution lower than the requisite resolution (⁠|$m_{\rm {DM}}\, \lesssim \, 10^5\, \mathrm{M}_{\odot }$|⁠, |$m_{\rm {gas}}\, \lesssim \, 10^4\, \mathrm{M}_{\odot }$|⁠, |$m_{*}\, \lesssim \, 10^4\, \mathrm{M}_{\odot }$|⁠), as discussed in Section 5.3.

Table 3.

Mass resolution of this work and previous studies. From left to right, the columns show the name of the simulation, the mass of one dark matter particle (m_DM), the initial mass of one gas particle or cell (m_gas), the initial mass of one star particle (m_*), the final simulation redshift (z_final), metal-mixing model, and references. Acronyms used for the metal-mixing models are HS: turbulence-induced metal mixing model (Shen et al. 2010; Hirai & Saitoh 2017; Saitoh 2017), SF: metal mixing in star-forming regions (Shen et al. 2015; Hirai et al. 2015, 2017), and MM: metal mixing computed together with mass transfer among moving mesh (Springel 2010).

Simulations	mDM	mgas	m*	zfinal	Metal mixing	References
	(M⊙)	(M⊙)	(M⊙)
This work	7.2 × 104	1.3 × 104	4.5 × 103	0	HS	–
Eris	9.8 × 104	2.0 × 104	6.0 × 103	0	SF	Shen et al. (2015)
FIRE (Low)	2.3 × 106	4.5 × 105	4.5 × 105	0	–	van de Voort et al. (2015)
FIRE (Fiducial)	2.8 × 105	5.7 × 104	5.7 × 104	0	–	van de Voort et al. (2015)
FIRE (High)	3.5 × 104	7.1 × 103	7.1 × 103	2.4	–	van de Voort et al. (2015)
IllustrisTNG100-1	7.5 × 106	1.4 × 106	1.4 × 106	0	MM	Naiman et al. (2018)
IllustrisTNG100-2	6.0 × 107	1.1 × 107	1.1 × 107	0	MM	Naiman et al. (2018)
Aquila (G3-CK)	2.2 × 107	3.5 × 106	1.8 × 106	0	–	Haynes & Kobayashi (2019)
Auriga	4.3 × 104	8.1 × 103	8.1 × 103	0	MM	van de Voort et al. (2020)
Auriga	3.6 × 104	6.7 × 103	6.7 × 103	0	MM	van de Voort et al. (2022)

Metal mixing also plays a vital role in correctly computing the enrichment of r-process elements. Shen et al. (2015) improved the agreement between the predicted and observed scatter of [Eu/Fe] abundance ratios when they implemented a model for metal mixing in the star-forming region. Hirai & Saitoh (2017) found that models without metal mixing overproduce RPE stars in dSphs. Thus, they calibrated the diffusion coefficient for metal mixing to suppress RPE star formation in dSphs. This study uses their value to compute the fraction of r-II stars in low-metallicity environments.

Both the simulations of this study and van de Voort et al. (2020) predict that very metal-poor stars are formed in the early phase of galaxy formation. van de Voort et al. (2020) reported that the median formation redshifts were z = 4.7 (12.6 Gyr ago, |$\mathrm{[Fe/H]}\, \lt \, -$|2) and z = 6.2 (12.9 Gyr ago, |$\mathrm{[Fe/H]}\, \lt \, -$|3). They also mentioned that the ages of outliers (in [r-process/Fe]) are older than the median in the same metallicity bin. Our simulation also finds that very and extremely metal-poor stars are formed at high-z: z = 7.530 (13.18 Gyr ago, |$\mathrm{[Fe/H]}\, \lt \, -$|2) and z = 7.844 (13.22 Gyr ago, |$\mathrm{[Fe/H]}\, \lt \, -$|3). We also confirm that the median age of r-II stars is older than all-stars: z = 7.836 (13.22 Gyr ago, |$\mathrm{[Fe/H]}\, \lt \, -$|2) and z = 10.44 (13.41 Gyr ago, |$\mathrm{[Fe/H]}\, \lt \, -$|3). The older ages in our study compared to van de Voort et al. (2020) could be due to differences in the efficiency of metal mixing and initial conditions. Nevertheless, these results suggest that the median ages of very metal-poor stars are over 10 Gyr, and those of r-II stars are even older.

While this study and van de Voort et al. (2020) predict a higher fraction of RPE stars formed in accreted components (Fig. 8a) than that of non-RPE stars, Naiman et al. (2018) could not find any correlation with the Eu abundance and assembly history. This discrepancy could come from the difference in the resolution. The average mass gas cell in Naiman et al. (2018) is |$1.4\times 10^6\, \mathrm{M}_{\odot }$|⁠. With this resolution, it is not possible to resolve UFD-like building blocks, which are key to forming r-II stars (Section 5.1). Most of the r-II stars in their simulation would be formed in local inhomogeneities within the disc or haloes.

5.5 Caveats

For [Fe/H] ≳ −1, the simulated [Eu/Fe] abundance ratios in Fig. 4 are constant towards higher metallicity, contrary to the observations, which exhibit a decreasing trend of [Eu/Fe]. This result is because we adopt a power-law index of −1 for the delay-time distributions of NSMs, which is the same power-law index adopted for SNe Ia. The effects of increasing [Eu/Fe] ratios by NSMs and the expected decrease due to SNe Ia offset each other under this assumption (Hotokezaka et al. 2018). As a result, the fraction of r-II stars at high metallicity is expected to be overestimated. This problem could be resolved by considering r-process production from short time-scale events (Côté et al. 2019), a longer minimum delay time for SNe Ia (Wanajo et al. 2021), natal kicks, and inside-out disc evolution (Banerjee, Wu & Yuan 2020), possibly in combination. These effects start to dominate for |$\mathrm{[Fe/H]}\, \gtrsim -1$|⁠. Therefore, the lack of modelling of these effects does not largely affect the [Eu/Fe] distribution with |$\mathrm{[Fe/H]}\, \lesssim -1$|⁠.

This simulation tends to produce more r-II stars in this metallicity range compare to observations. There are two reasons for this tendency: (1) Since this simulation assumes the similar delay-time distributions for SNe Ia and NSMs, the average [Eu/Fe] ratios are higher than the observations for |$\mathrm{[Fe/H]}\, \gt \, -1$|⁠, as discussed above. Due to the high average [Eu/Fe] abundance ratios, the [Eu/Fe] in RPE gas are not significantly diluted by the surrounding gas; and (2) the value of the diffusion coefficient assumed in this simulation could be insufficient. In this study, we adopt a lower value of the scaling factor for the metal diffusion (⁠|$C_{\rm {d}}\, =0.01$|⁠). However, Hirai & Saitoh (2017) could not constrain the value of C_d larger than 0.01. Adopting a higher value of C_d would suppress the formation of r-II stars.

The existence and persistence of the local inhomogeneity caused by an NSM (Fig. 11) would depend on the resolution and the efficiency of the dilution of metals. As discussed in Section 5.1, the mass of the gas clump is typically |${\sim}10^6\, \mathrm{M}_{\odot }$|⁠, which is the smallest mass scale to evaluate the physical quantities in our SPH simulation. The time-scale of the dilution is also affected by the metal-diffusion coefficient (Hirai & Saitoh 2017). This should be tested in future higher-resolution simulations.

The dispersions seen in Fig. 15 could also be due to insufficient metal mixing in our simulation. Since metal mixing is caused by turbulence in the ISM, which is not resolved in the simulation, the efficiency of metal mixing on a galactic scale is unknown. Hirai & Saitoh (2017) have shown that the scatter of [Eu/Fe] can be artificially caused in SPH simulations with low efficiencies of metal mixing. Although we adopt the recommended value from their study, which is consistent with other estimates on different scales (e.g. Horiuti 1987; Meneveau & Katz 2000), there is still the possibility that our simulation does not have sufficient efficiency to thoroughly mix metals. Determining the dispersions of [Eu/Fe] from high-resolution spectroscopic observations in accreted components and dwarf galaxies could help empirically resolve this issue.

In this study, we have concluded that r-II stars can be formed in low-mass building blocks enhanced in r-process elements or from locally enhanced gas clumps in relatively massive building-block galaxies and the main halo. The birthplaces of r-II stars depend on metallicity. For |$\mathrm{[Fe/H]}\, \lt -2$|⁠, most of the r-II stars are formed in accreted components (Fig. 8a). Especially for |$\mathrm{[Fe/H]}\, \lt -2.5$|⁠, r-II stars are formed in low-mass dwarf galaxies with stellar mass |${\sim}10^5\, \mathrm{M}_{\odot }$| when their gas was enhanced in r-process elements (haloes 1, 2, 6, 9, 11, 13, and 14). NSMs are allowed to occur at low metallicity in such low-mass systems because of their slow chemical evolution. For higher metallicity, r-II stars are formed in local inhomogeneities of r-process abundances in massive building blocks or the main halo. Although the local enhancement of r-process abundances is erased within a relatively short time-scale (∼10 Myr), r-II stars can be formed if the system is sufficiently massive to host NSMs that lead to star formation in the RPE gas clumps (Ji et al. 2016a; Roederer et al. 2016; Ojima et al. 2018; Brauer et al. 2019; Wanajo et al. 2021). Simultaneously, it is possible to explain the existence of high-metallicity RPE stars by star formation in local inhomogeneities of the gas-phase r-process abundances.

6 CONCLUSIONS

In this paper, we perform a cosmological zoom-in simulation of an MW-like galaxy to investigate the birth environments and chemo-dynamical properties of r-II stars, one of the first with sufficient mass (⁠|$m_{\rm {gas}} = 1.3\, \times \, 10^4\, \mathrm{M}_{\odot }$|⁠) and time resolution (∼10 Myr) to accomplish this purpose. This simulation produces star-to-star dispersions of [Eu/Fe] ratios at |$\mathrm{[Fe/H]}\, \lt -2$| (Fig. 4), similar to the observations, thanks to appropriate modelling of NSMs and metal diffusion. We confirm that most r-II stars with |$\mathrm{[Fe/H]}\, \lt -1$| exhibit halo-like orbits reported in previous observations (Roederer et al. 2018; Gudin et al. 2021; Shank et al. 2022a).

We find that over 90 per cent of r-II stars are formed in the early epochs of the galaxy formation (⁠|$t\, \lesssim 4\,$|Gyr), suggesting that the majority of r-II stars are formed in small building blocks and accreted into the main halo later (see Fig. 7). We also find an increasing trend in the fraction of accreted components towards low metallicity (see Fig. 8) and that the accreted fraction of r-II stars is higher than that of non-RPE stars. These results indicate that the majority of r-II stars were formed in disrupted dwarf galaxies over 10 Gyr ago.

The r-II stars are formed in gas clumps enhanced in r-process elements. For |$\mathrm{[Fe/H]}\, \lt -2.5$|⁠, r-II stars are formed in dwarf galaxies with |$M_*\lesssim 10^5\, \mathrm{M}_{\odot }$|⁠. Because of the small gas mass (⁠|$M_{\rm {gas}}\, \lesssim 10^7\, \mathrm{M}_{\odot }$|⁠) of these galaxies, the average gas-phase [Eu/Fe] is enhanced once an NSM occurs. In such a system, most stars are formed as RPE stars. We find an increasing trend in the fraction of r-II stars towards haloes with lower stellar mass. For |$\mathrm{[Fe/H]}\, \gt -2.5$|⁠, r-II stars are mainly formed in locally RPE gas clumps in relatively massive accreted components or the main halo. This formation channel could account for the r-II stars found in classical dSphs (e.g. Reichert et al. 2021) and the MW star with |$\mathrm{[Fe/H]}\, =-1.27$| (Xing et al. 2019).

The chemo-dynamical properties of r-II stars also provide valuable information on their origin. The r-II stars from accreted components tend to be clustered in their dynamical phase space, while the majority of r-II stars formed in situ have disc-like kinematics. We find that some dynamical groups of r-II stars exhibit decreasing [Mg/Fe] towards higher [Fe/H], similar to present dwarf galaxies. We also find that r-II stars in the same halo exhibit low dispersions of [Fe/H] (since r-II stars are formed on relatively short time-scales, ∼100 Myr) and somewhat larger dispersions of [Eu/Fe] in each halo, as also seen from the observations. Furthermore, the fraction of halo r-II stars found in our simulation is commensurate with observations from the RPA, and crucially, the distribution of the predicted [Eu/Fe] for halo r-II stars well matches that observed. The dispersion of [Eu/Fe] in each halo could also explain the r-II star in the Indus stellar stream (Hansen et al. 2021). These results clearly indicate that observations of r-process elements, together with α-elements and kinematics, help distinguish metal-poor stars’ formation environments.

Efforts to identify and analyse RPE stars and perform refined high-resolution galaxy formation simulations will greatly improve our understanding of the MW’s formation in the earliest epochs. As shown in this study, most r-II stars are formed in the accreted components in the early Universe. Detailed comparisons with the chemo-dynamical properties of RPE stars and galaxy-formation simulations can constrain the evolutionary histories of the building blocks of the MW and the astrophysical sites responsible for the production of the heaviest elements. Observationally, the RPA and Gaia are greatly increasing the number and accuracy of the derived chemo-dynamical properties of RPE stars. Eventually, individual stellar abundances and kinematics will be used to compare with future galaxy-formation simulations resolving individual stars (e.g. Hu et al. 2017; Hu 2019; Lahén et al. 2020; Fujii et al. 2021a,b; Gutcke et al. 2021, 2022; Hirai et al. 2021). These efforts will open up a new window to reveal galaxy formation and nucleosynthesis in the early Universe.

ACKNOWLEDGEMENTS

We appreciate the anonymous referee for providing valuable suggestions to improve the manuscript. We also thank Ian U. Roederer for his careful reading of the draft. We are grateful for fruitful discussions with Alexander P. Ji and Ting S. Li. This work was supported in part by JSPS KAKENHI grant numbers JP21J00153,JP20K14532,JP21H04499,JP21K03614,JP22H01259,JP21H05448,JP18H05437,JP21H00055,JP21H04500,JP21K03633,JP20H05861,JP21H04496,JP19H01931,JP20H01625, and JP21K11930, MEXT as ‘Program for Promoting Researches on the Supercomputer Fugaku’ (Toward a unified view of the Universe: from large-scale structures to planets, Grant No. JPMXP1020200109), JICFuS, grants PHY-1430152; Physics Frontier Center/JINA Center for the Evolution of the Elements (JINA-CEE), and OISE-1927130: The International Research Network for Nuclear Astrophysics (IReNA), awarded by the US National Science Foundation. Numerical computations and analysis were carried out on Cray XC50 and computers at the Center for Computational Astrophysics, National Astronomical Observatory of Japan, and the Yukawa Institute Computer Facility. This research has made use of NASA’s Astrophysics Data System and astropy,³ a community-developed core python package for Astronomy (Astropy Collaboration 2013, 2018).

DATA AVAILABILITY

Data will be available with a reasonable request to the corresponding author. celib is available at https://bitbucket.org/tsaitoh/celib.

Footnotes

1

[A/B] = log₁₀(N_A/N_B) − log₁₀(N_A/N_B)_⊙, where N_A and N_B are the number densities of elements A and B, respectively.

2

According to Burris et al. (2000), the pure s- and r-process ratios of [Ba/Eu] are [Ba/Eu]|$_{s} = +1.5$| and [Ba/Eu]|$_{r}\, = -0.82$|⁠, respectively.

3

http://www.astropy.org

Author notes

JSPS Research Fellow.