Metal Mixing in the R-Process Enhanced Ultra-Faint Dwarf Galaxy Reticulum II

The ultra-faint dwarf galaxy Reticulum~II was enriched by a single rare and prolific r-process event. The r-process content of Reticulum~II thus provides a unique opportunity to study metal mixing in a relic first galaxy. Using multi-object high-resolution spectroscopy with VLT/GIRAFFE and Magellan/M2FS, we identify 32 clear spectroscopic member stars and measure abundances of Mg, Ca, Fe, and Ba where possible. We find $72^{+10}_{-12}$% of the stars are r-process-enhanced, with a mean $\left\langle\mbox{[Ba/H]}\right\rangle=-1.68~\pm~0.07$ and unresolved intrinsic dispersion $\sigma_{\rm [Ba/H]}<0.20$. The homogeneous r-process abundances imply that Ret~II's metals are well-mixed by the time the r-enhanced stars form, which simulations have shown requires at least 100 Myr of metal mixing in between bursts of star formation to homogenize. This is the first direct evidence of bursty star formation in an ultra-faint dwarf galaxy. The homogeneous dilution prefers a prompt and high-yield r-process site, such as collapsar disk winds or prompt neutron star mergers. We also find evidence from [Ba/H] and [Mg/Ca] that the r-enhanced stars in Ret~II formed in the absence of substantial pristine gas accretion, perhaps indicating that ${\approx}70$% of Ret~II stars formed after reionization.


INTRODUCTION
Ultra-faint dwarf galaxies (UFDs) are Milky Way satellite galaxies with luminosities M V > −7.7 (stellar masses 10 5 M , Simon 2019). UFDs appear to form all their stars in the first 1-2 billion years, before their arXiv:2207.03499v2 [astro-ph.GA] 23 Feb 2023 star formation is cut off by reionization (Benson et al. 2002;Brown et al. 2014). UFDs probe the extreme lowmass end of galaxy formation, where star formation is inefficient and massive stars form stochastically, resulting in intermittent feedback and incomplete sampling of nucleosynthetic sources (Koch et al. 2008(Koch et al. , 2013Simon et al. 2010;Frebel et al. 2010;Frebel & Norris 2015;Ji et al. 2016aJi et al. , 2019c. UFDs are also relics of early galaxy formation, providing a unique window into the first stars and galaxies in a pre-reionization universe, as well as a clean probe of the first metal-free Population III stars (e.g., Bovill & Ricotti 2009;Salvadori & Ferrara 2009;Frebel & Bromm 2012;Ji et al. 2015). Since many halo stars with [Fe/H] < −2.5 likely form in UFD-like environments, even if they later grow into or accrete into larger systems (Brauer et al. 2019), it is crucial to understand the star formation conditions for UFDs to interpret the most metal-poor stars. To understand these early properties, the red giant branch stars in UFDs have been the subject of intense spectroscopic study. The last 15 years have resulted in high-resolution spectra of 100 stars across ∼20 UFDs with detailed elemental abundances (see Frebel & Norris 2015;Simon 2019;Ji et al. 2019cJi et al. , 2020a for a description of the basic characteristics and chemical evolution trends).
Reticulum II (Ret II) is a UFD discovered in the Dark Energy Survey (Bechtol et al. 2015;Koposov et al. 2015a), located only 32 kpc away. Initial followup spectroscopy showed that its velocity dispersion, mean metallicity, and metallicity dispersion were consistent with typical UFDs (Simon et al. 2015;Koposov et al. 2015b;Walker et al. 2015, henceforth S15; K15; W15). Subsequent high-resolution spectroscopy surprisingly showed that most Ret II stars displayed some of the highest r-process enhancements known (Ji et al. 2016a,b;Roederer et al. 2016b, henceforth J16;R16). By comparing to other UFDs (which display unusually low neutron-capture element abundances, Frebel & Norris 2015;Ji et al. 2019c), the clear conclusion is that Ret II experienced enrichment from a single rare and prolific r-process event. The source of the r-process elements is still debated, as it could be consistent with rprocess nucleosynthesis in a prompt neutron star merger or rare core-collapse supernova (e.g., Ji et al. 2016a;Beniamini et al. 2016;Safarzadeh & Scannapieco 2017;Safarzadeh et al. 2019a;Ojima et al. 2018;Siegel et al. 2019;Tarumi et al. 2020;Molero et al. 2021;Jeon et al. 2021;Cowan et al. 2021).
The single r-process event in Ret II provides a unique opportunity to probe metal mixing in a UFD. Because all the r-process elements (including barium and europium) were deposited in a single enrichment event, the distribution of r/H ratios in Ret II stars depends only on the overall amount of enriched gas and the and homogeneity of metal mixing into the gas. This contrasts with elements synthesized by more common sources like supernovae or asymptotic giant branch (AGB) stars, since the frequency of element production interacts with metal mixing to produce the distribution of stellar abundances (e.g., Krumholz & Ting 2018;Emerick et al. 2019Emerick et al. , 2020. A direct constraint on metal mixing by measuring the [r/H] distribution could play a major role in interpreting UFD abundances and formation histories (e.g., Frebel & Bromm 2012;Ji et al. 2015;Webster et al. 2016;Tarumi et al. 2020).
We thus present a detailed spectroscopic study of Ret II chemical abundances obtained with multiobject spectroscopy using VLT/FLAMES and Magellan/M2FS. We find 32 clear member stars and 8 more candidates, the most spectroscopically confirmed members to date in Ret II. About half the stars have Ba and Fe constraints, while a third have Mg and/or Ca measurements as well. Our primary focus is measuring the distribution of [Ba/H] in these stars, which is a tracer of r-process enrichment in Ret II due to the high r-process enhancement and negligible s-process contribution in Ret II (Ji et al. 2016a). Since this is the largest spectroscopic sample of members yet, we also explore more general kinematics, binarity, chemical evolution, and spatial gradients. Section 2 presents the spectroscopic observations and data reduction. Section 3 describes the velocity and chemical abundance analysis methods, as well as membership determination including auxiliary information from the Dark Energy Survey (DES, DES Collaboration et al. 2018) and Gaia EDR3 (Gaia Collaboration et al. 2016Lindegren et al. 2021). Section 4 gives our results for the Ret II radial velocity distribution, chemical abundance trends, Fe and Ba distributions, and radial gradients. Section 5 discusses the implications of our measurements on metal mixing in dwarf galaxies, the origin of the r-process elements, and chemical evolution in Ret II. We summarize and conclude in Section 6. Multi-epoch velocities are provided in Appendix A. A major systematic for our Ba results is microturbulence, which is discussed extensively in Appendices B and C.

OBSERVATIONS AND DATA REDUCTION
We observed Reticulum II with VLT/FLAMES in October 2017 (Pasquini et al. 2002), and with Magellan/M2FS in September 2016 at medium-resolution and November 2017 at high-resolution (Mateo et al. 2012). Table 1 contains details about which stars were observed at which settings. Note that all signal-to-noise ratios (SNR) quoted in this paper refer to the SNR per pixel.

VLT/FLAMES, GIRAFFE
The FLAMES/GIRAFFE setup on the VLT UT2 provides high-resolution spectra of ∼100 stars over a field of view of diameter 0.4 degrees. Observations were taken in visitor mode on 26-27 October 2017 with excellent weather. We used the HR14A setting, covering one order from 6300−6500Å with R ∼ 18000. Targets were selected based on our own photometry of public DES Y1 images following Koposov et al. (2015a). We chose targets near the fiducial CMD within the single field we targeted. The total FLAMES exposure time was 11.8h, with most exposures being 3000s but a few exposures of 2400s and 3600s at the end of the night.
Data were reduced with the standard ESO pipeline, which provides flat-corrected and wavelength calibrated 1D flux and error spectra. The 1D spectra are extracted to a common rebinned dispersion without cosmic ray rejection or sky subtraction. For each object, we removed cosmic rays in 1D by normalizing individual exposures by their median flux, then masking pixels with >5σ deviations from the combined median spectrum. Care was taken not to mask pixels associated with variable sky lines. Sky subtraction was performed in 1D mostly following Battaglia et al. (2008). For each exposure, we constructed a master sky spectrum from ∼15 sky fibers using an inverse-variance weighted mean. The master sky flux was split into two components, a sky emission line and a continuum component. The line component was used to identify wavelength bins associated with emission lines. Then for each object spectrum, we also split the flux into emission lines and continuum, rescaled the master sky line flux to match the object emission line flux by minimizing the L1 norm (total absolute deviation at wavelengths associated with sky emission lines), applied the same scaling factor to the sky continuum, and subtracted the rescaled master sky from the object spectrum. Visual inspection of the sky-subtracted spectra suggests this procedure was generally effective, with no correction to the line spread function or wavelength recalibration needed. Still, there are sometimes sky subtraction residuals from spatially variable sky lines, which does impact our Ba line of interest (see Section 2.4). Final coadded spectra were obtained using an inversevariance weighted average of individual exposures.

Magellan/M2FS HiRes
We obtained high-resolution spectra of Ret II stars with M2FS on 16-17 November 2017. We used the HiRes mode with 180µm slits, providing R ∼ 18000.
The detectors were binned 2x2 with 4 amplifier slow readout. Two different blocking filters were used to observe the targets (one on each M2FS channel), based on a visual examination of the VLT spectra. For fainter targets with unclear Ba detections or upper limits, we used the BulgeGC1 filter, which includes 24 fibers covering 6 orders from 6100 − 6700Å, including the Ba line at 6496.7Å. The Ba line at 6141Å is on the blue end of the filter cutoff and cannot be used. For brighter targets that already had clear Ba detections in the VLT data or upper limits, we instead used the MgWide filter, which includes 28 targets covering 4 orders from 5150 − 5400Å. 6 sky fibers were allocated for each arm. The total exposure time was 14h.
The data were reduced with a custom pipeline 1 . Each of 4 amplifier images was bias subtracted using the overscan and stitched into one image, then had dark current subtracted. Every science frame was associated with a single arc and flat obtained closest in time to the science frame. The object trace was fit to each flat using a 5th order Legendre polynomial. Scattered light was subtracted from every flat and science frame by fitting the inter-object regions with a 2D Legendre polynomial of degree 5 in either direction. Twilight flats were used to determine throughput corrections for each fiber.
The wavelength calibration was motivated by Kelson (2003) and adapted for fiber spectroscopy. An initial feature identification was done once by hand, extracting all orders of each fiber and using the IRAF identify command to manually identify positions of 50 − 70 arc lines in each order of each fiber in the X (wavelength) direction on the CCD (Tody 1986(Tody , 1993. These identifications were then turned back into 2D coordinates using the trace functions. The actual wavelength calibration was performed in 2D, finding sources in each arc frame using Source Extractor (Bertin & Arnouts 1996), and matching the detected sources to identified lines using a KD tree 2 . We then fit a 5th order Legendre polynomial for the wavelength solution, iteratively rejecting outliers, and using lines from all object fibers to fit the overall distortion. A single X and Y pixel offset is allowed for each fiber (but not orders within a fiber) to account for any movement of the fibers in the pseudo-slit. In total, 40 − 50 lines were identified and used in each order for the MgWide filter, and 10 − 30 lines were identified and used for the BulgeGC1 filter (fewer due to the satu-1 https://github.com/alexji/m2fs reduction 2 2D wavelength calibration was necessary for the BulgeGC1 filter because the ThArNe arcs taken for this setting were extremely saturated, introducing many spurious features in 1D extracted arcs. rated arcs). The final wavelength solution has a typical RMS <0.01Å in both arms. Data were then extracted using flat-relative optimal extraction (Zechmeister et al. 2014), which we found performed better than fitting a functional form to the object profile.
To perform sky subtraction, we linearly rebinned the extracted spectra onto a uniform wavelength grid, then followed essentially the same sky subtraction procedure as the VLT data. The main differences were that the M2FS data has multiple orders, so sky subtraction was done independently for each order; and since the MgWide filter has few sky lines, we did not rescale the master sky spectrum to match line strengths, instead just directly subtracting the throughput-corrected master sky. There were clear differences in the line spread function for different fibers, resulting in residuals around sky lines. We thus rejected data around sky lines. Different exposures were coadded order-by-order with an inversevariance weighted average. Coadded orders were then continuum-normalized separately in smhr 3 before being stitched into a single spectrum.

Magellan/M2FS MedRes
We conducted two sets of medium resolution observations of Ret II stars using M2FS on 2016 September 6 and 10, totaling 6.72 hr of integration time. We used the MedRes grating on the 'R' spectrograph, 95 µm slits, 2x2 binning, and the MedRes Ba(23) filter, which transmits one order with 4450 ≤ λ ≤ 4615Å. This setup yields R ∼ 9, 000, as measured from individual Ar or Th emission lines in the comparison lamp spectra. We performed data reduction, extraction, wavelength calibration, sky subtraction, co-addition, and continuum normalization following the procedures described in Roederer et al. (2016a), modified for use with MedRes spectra.
We extracted the two sets of Ret II MedRes spectra (one set from each night) separately. We measured radial velocities of each star in each observation by cross-correlating (using the IRAF fxcor task) its spectrum against a synthetic metal-poor template spectrum smoothed to the same spectral resolution. Repeat observations of probable Ret II members show a standard deviation of 1.7 km s −1 , which we regard as the uncertainty of an individual measurement. We also observed two comparison stars, CD −24 • 1782, and CS 31082-001 using the same M2FS MedRes setup. Their radial velocities, after applying Heliocentric corrections computed using the IRAF rvcorrect task, agree with pub-3 https://github.com/andycasey/smhr, originally described in Casey (2014) and expanded in Ji et al. (2020b) lished values (Roederer et al. 2014;Hansen et al. 2015) to better than 1.7 km s −1 . which we regard as the systematic uncertainty of our measurements. Combining the individual and systematic uncertainty, we adopt a total velocity uncertainty of 2.4 km s −1 for the MedRes velocities. However, when later comparing repeat velocity measurements of Ret II stars, we found a systematic offset of ≈10 km s −1 . We thus decided not to use the M2FS MedRes velocities in this work, though we report their values in Appendix A.

Comment on Sky Subtraction and Ba Lines
Our primary Ba line at 6496.7Å is very close to a strong sky line at 6498.7Å, and our Ba abundances are potentially susceptible to sky subtraction residuals. For the VLT data with SNR > 25, we were able to obtain good simultaneous fits to the Ba lines and the sky line residuals. We decided stars with SNR < 25 were too strongly impacted by sky subtraction to have reliable Ba measurements, especially since a small error can make a big asymmetric difference in the abundances for saturated lines (Section 3.4.2). For the M2FS data, the signal-to-noise was lower and the heliocentric correction moved the Ba line right into the sky line, so none of the Ba 6496.7Å line measurements from M2FS were reliable. There was also a 6141Å Ba line located on the filter cutoff, but after investigation we decided it was adversely affected by scattered light subtraction systematics and thus not sufficiently reliable for abundance measurements.

ANALYSIS
We assigned each star a numerical ID. IDs 1 − 26 were red giant branch stars identified as members or candidate members by Simon et al. (2015) (a superset of Walker et al. 2015;Koposov et al. 2015b), sorted by magnitude. IDs 97 − 200 were assigned arbitrarily to other observed stars.

Photometry and Astrometry
For all observed stars, we queried the DES DR1 catalog for griz magnitudes using NOAO datalab, which were dereddened using the DES DR1 reddening coefficients (DES Collaboration et al. 2018). The exception is star ID 1, which is saturated in the DES data release but had photometry determined from an individual DECam exposure in Simon et al. (2015). The subscript 0 (e.g., g 0 ) indicates dereddened magnitudes. RA and Dec were taken from the DES catalog (differing by ∼0. 1 compared to Gaia). We also computed the elliptical distance r e for each star, assuming Ret II is centered at 03:35:47.83 −54:02:47.8 with a position angle of 68 • , ellipticity of 0.6, and half-light radius 6.3 arcmin (Mutlu-Pakdil et al. 2018). We adopted a distance modulus of 17.5 whenever needed (Mutlu-Pakdil et al. 2018). Proper motions µ α cos δ and µ δ were obtained by cross-matching with Gaia EDR3 (Lindegren et al. 2021). Our spectroscopic target selection did not use proper motions, as it was performed before Gaia DR2 was released.

Radial Velocities
Radial velocities were measured with weighted crosscorrelation against a high-resolution template spectrum of HD 122563 obtained with the MIKE spectrograph and shifted to rest-frame (Bernstein et al. 2003). Each spectrum was first normalized with a 3rd order polynomial. Then, for a range of velocities, we shifted the template spectrum by that velocity and calculated the χ 2 of pixels within a specific wavelength interval. For the VLT data and the M2FS BulgeGC1 data, we used Hα to measure velocities, cross-correlating from 6550 − 6575Å. For the M2FS MgWide data, we used the Mg triplet, cross-correlating from 5150−5200Å. Velocities were then found as the minimum of the χ 2 contour, and 1 σ statistical velocity errors were determined by ∆χ 2 = 1 (e.g., Martini et al. 2006). For the M2FS data, we used just one of the orders for velocity measurement: order 54 for BulgeGC1 and order 69 for MgWide. We did not attempt to determine very detailed radial velocities for the MedRes M2FS data, given its lower resolution.
We determined systematic velocity uncertainties from repeated velocity measurements. For both the VLT and M2FS run, every star was observed with 7-15 individual exposures over two nights. We thus measured the velocity and statistical uncertainty for each individual exposure. Then, following Li et al. (2019), we took all pairs of repeated velocity measurements with σ v < 30 km s −1 and fit the following Gaussian-plus-outlier model: where v i , v j , σ v,i , and σ v,j are the individual velocities and uncertainties; f is the fraction of pairs that are good; σ outlier is a large value characterizing a broad background outlier model; and F (σ v ) = σ 2 v,floor + (k × σ v ) 2 is a rescaling of the velocity uncertainty with a scale factor k and a systematic floor σ v,floor . The VLT data had f = 0.98, σ v,floor = 0.69 km s −1 , and k = 1.02. The M2FS MgWide data had f = 0.66, σ v,floor = 0.74 km s −1 , and k = 1.64. The M2FS BulgeGC1 data had f = 0.93, σ v,floor = 3.20 km s −1 , and k = 0.11. In this case, having k < 1 suggests that the individual velocities are dominated by systematic errors, which is primarily due to the low SNR of individual exposures causing poor continuum fits and template matches. We thus conservatively take k = 1 for the BulgeGC1 data. These values of k and σ v,floor are applied to generate our final velocity uncertainties. Figure 1 shows the individual velocity measurements for all observed stars. The left panel shows the 26 likely members from previous observations, while the right panel shows other observed stars. The red bars on the lower axis indicate stars that are clear members, while orange bars indicate candidate members (see Section 3.3 for the definitions). For comparison, in the left panel we also plot velocities derived by S15, K15, J16, and R16. Since many literature velocities were measured using the same spectra, duplicate velocity measurements from W15 and S15 were removed from this figure.
For the 40 members and candidate members (determined in Section 3.3), we determined velocities by averaging our VLT and M2FS velocities with the literature velocities from S15, K15, J16, and R16, for a total of up to 6 independent velocity measurements. We combine all available literature velocities with an inverse-variance weighted mean. Details are given in Appendix A, but this includes estimating a systematic velocity offset for each data sample using the S15 velocities as a reference, and identifying binary star candidates using a chisquared test. Five candidate binary stars are identified in Table 1 and as grey bars on the top axis of Fig 1. These binary star velocities were also found with a weighted mean, but the uncertainty indicates the range of velocities. Binary candidates are excluded from the velocity and velocity dispersion calculations. Our focus in this paper is on chemical abundances, so this assessment of binarity and velocity systematics is incomplete, and a more comprehensive future study is warranted.

Membership
We used four criteria to establish membership within our selected targets: radial velocities in the range 53 km s −1 < v r,hel < 74 km s −1 ; proper motion within 2 units of Mahalanobis distance from the Ret II mean proper motion of (µ * α , µ δ ) = (2.39, −1.36) (McConnachie & Venn 2020) 4 ; position within 2 elliptical half light radii of the Ret II center 5 ; and g 0 − r 0 color between 13 Gyr, α-enhanced Dartmouth isochrones (Dotter et al. 2008) with [Fe/H] = −2.5 and [Fe/H] = −1.5, allowing a 0.03 mag buffer on each side of the isochrone given an expected reddening uncertainty. Since we used hard cuts, edge cases near these boundaries were inspected individually.
The results are shown in Figure 2. 70 of the 129 total spectroscopic targets were rejected as being outside of our radial velocity limits (small blue points). Of the remaining 59 stars, 18 stars had consistent radial velocities but were rejected based on clearly discrepant spatial location, CMD position, or proper motion (cyan crosses). One additional star (ID 191) was also rejected as having [Fe/H] ∼ −0.5 from our later analysis, too high to be part of Ret II. This leaves 40 stars, of which 32 stars were confident members, clearly matching at least three of the four criteria (large red points). In the rest of this paper, we refer to these stars as "clear members." The remaining 8 stars were highlighted as uncertain members (small orange squares), which we refer to as "candidate members" in the rest of this paper. Stars 188 and 143 are CMD members but at somewhat large distance and proper motion away from the galaxy core to be considered very confident members. Stars 14, 97, and 142 are redward of the isochrones, requiring high metallicities [Fe/H] −1.5 (or perhaps large carbon bands) to be considered part of Ret II. Our subsequent analysis finds none of these stars has such a high metallicity. Star 16 is blueward of the isochrone and has an unusually shallow Hα line shape suggesting it is a hotter star. Finally, stars 151 and 154 have consistent kinematics and would be CMD non-members, but they are located where evolved blue stragglers or an unusually young stel-lar population might be found, as indicated by a 10Gyr isochrone in the top-left panel of Figure 2. In other dwarf spheroidal galaxies, blue stragglers tend to appear even younger (2.5 ± 0.5 Gyr, Santana et al. 2013) so it is possible these stars could indicate a younger population of stars in Ret II. However, the spectra of these potential member stars are too low SNR to reliable measure any abundances. Stars 100 and 102 are known BHB member stars (Simon et al. 2015). Other than rejecting the clear non-member star ID 191, we have avoided using metallicity information in the membership selection so as to remain unbiased in final abundance distributions.
Note that several Ret II members have slightly lower µ α and µ δ than most of the galaxy (at (µ α , µ δ ) = (+2.0, −2.0) compared to (+2.4, −1.3)). These 5 stars were also offset in the same direction in Gaia DR2. However, they have larger proper motion uncertainties individually consistent with the bulk of the galaxy, and they do not stand out in radial velocity. Future Gaia data releases can test whether these represent a true feature.

Chemical Abundances
We derived most abundances using a standard analysis using 1D Castelli & Kurucz (2004) model atmospheres and Local Thermodynamic Equilibrium (LTE) radiative transfer with MOOG (Sneden 1973) including scattering  and Barklem et al. (2000) damping 6 . Since Ba is the element of most interest, we analyzed Ba with both LTE and non-LTE, see Section 3.4.5. Our stellar parameters are given in Table 2, atomic data in Table 3, and abundance results in Table 4.

Stellar Parameters
Since our spectra cover a very small wavelength range with few lines, we determined effective temperatures for member stars from DES photometry and Dartmouth isochrones (Dotter et al. 2008). We adopted a 13 Gyr, [Fe/H] = −2.5, [α/Fe]= +0.4 isochrone as our fiducial isochrone. The isochrone was used to fit T eff as a function of g 0 − r 0 , and then we applied this to every member star. The brightest known member star (ID=1) is near the DES saturation limit with obviously incorrect r and i magnitudes in DES DR1. We adopted g 0 −r 0 = 0.80 from Simon et al. (2015) to determine this star's stellar parameters. Statistical uncertainties were calculated assuming a fixed error in g 0 − r 0 = 0.02 mag (the typical reddening uncertainty) resulting in a typical temperature error of 50-100 K. Systematic uncertainties were estimated by taking the largest difference be- The total temperature uncertainty was the quadrature sum of these two uncertainties.
The surface gravity log g was determined photometrically using the equation log g = 4.44 + log M + 4 log T eff /5780K+0.4(g 0 −µ+BC(g)−4.75) (Venn et al. 2017) where M = 0.75 ± 0.1M is the typical mass of an old red giant branch star, g 0 is the dereddened DES g magnitude, µ = 17.5 is the distance modulus to Ret II, and BC(g) is the Casagrande & VandenBerg (2014) bolometric correction. Besides the M and temperature uncertainties, we propagated a conservative 0.2 mag total uncertainty for the distance modulus and dereddening, as well as 0.03 mag uncertainty in the bolometric correction. The total log g uncertainty is about 0.2 dex for all stars. Using r instead of g led to the same log g within 0.05 dex, and using different metallicities for the bolometric correction led to differences within 0.02 dex.
The T eff and log g in Table 2 agree with two previous works that studied a total of 9 stars in Ret II using high-resolution spectroscopy (J16; R16). The stellar parameters in these two studies were derived using standard 1D-LTE methods, i.e., by balancing abundances of Fe lines with respect to excitation potential, ionization, and line strength, including a temperature recalibration to a photometric scale (Frebel et al. 2013). All nine stars agree within the 1σ stellar parameter uncertainties with no mean offset.
Our primary line of interest, the 6496Å Ba line, is saturated in essentially all of our stars with a detected Ba line. Microturbulence (ν t ) thus plays a major role in determining the final abundance and uncertainty, espe-cially for the two coolest stars in our sample (IDs 1 and 97). As we do not have enough lines to determine ν t self-consistently in our stars, we instead adopt a microturbulence relation based on log g from the metal-poor giants in Roederer et al. (2014): where the typical scatter around the relation is 0.13 km s −1 . We note that this relation is quite different from the direct ν t measurements in J16 and R16, but it eliminates a trend in Ba abundance with effective temperature that would otherwise be present. We have decided to adopt this relation for the rest of the main paper, and a complete investigation of microturbulence choices and the impact on our results is discussed in Appendices B and C. As we do not have Fe or α-element constraints for most stars, we assume everywhere a model metallicity of [M/H] = −2.5 and [α/Fe] = +0.4. We verified in stars with [Fe/H] and [Mg/Fe] measurements that changing these parameters produces negligible differences to the resulting abundances. The final stellar parameters and uncertainties for all members or candidate members are given in Table 2. We do not determine stellar parameters for the BHB stars (IDs 100 and 102).

VLT Data
In the VLT data, most Ret II member stars only have 1-5 significant absorption features: H-α, the Ba ii line at 6496.897Å, an Fe i line at 6494.98Å, a Ca i line at 6439.08Å, and sometimes a Ca i line at 6493.78Å. We use this data to derive Ba, Fe, and Ca abundances where possible. Atomic data used are given in Table 3.
Stellar parameter uncertainties were found by calculating new curves of growth for 1σ differences in stellar parameters, redoing the above calculations, and taking the difference. These are added in quadrature to the statistical uncertainties.
As part of the above procedure we fit the Ca line at 6493.8Å. However, there is a stronger Ca line at 6439Å that is more often detected and also less susceptible to non-LTE effects (Mashonkina et al. 2017). The Ca abundance is thus derived using this 6439Å line with smhr. This includes the normalization, equivalent width measurement, and stellar parameter uncertainty propaga-tion . Formal 4 σ upper limits for Ca were also calculated in smhr by synthesizing a Ca line such that ∆χ 2 = 16.

M2FS HiRes Abundances
We used smhr to normalize and stitch coadded orders, fit equivalent widths, measure element abundances from MOOG, and propagate stellar parameter uncertainties (see description in Casey 2014; Ji et al. 2020b). The primary use of this data is to measure the Mg abundances, which we derive from equivalent widths of the Mg b lines at 5172 and 5183Å. Useful abundances were ultimately only derived for the MgWide arm, since the BulgeGC1 arm only had fainter and warmer stars, and the 6496Å Ba line was completely blended with a sky line.

M2FS MedRes Abundances
The primary line of interest in these data is the Ba 4554Å line. We measure the abundance from this line with a procedure identical to that used for the VLT data, except that when fitting the equivalent width we manually define valid continuum wavelength ranges instead of modeling the many absorption lines. The results are provided in Table 4, and they are consistent with the Ba abundances derived from the 6496Å line in the VLT data but with larger uncertainties. We thus do not use these abundances further. However, in the future this medium-resolution mode may be an efficient way to search for strong Ba lines in other UFDs.

Barium NLTE Abundances
Ba can be substantially affected by non-LTE effects (NLTE, e.g., Bergemann & Nordlander 2014;Mashonkina et al. 2014;Gallagher et al. 2020) especially when the lines are near saturation. For metal-poor giants, NLTE makes the Ba 6497 line stronger, resulting in lower inferred Ba abundances when including NLTE. To determine this quantitatively, we computed NLTE abundances for the Ba 6497 line using an updated version of the MULTI radiative transfer code (Carlsson 1986;Bergemann et al. 2019;Gallagher et al. 2020) and MARCS model atmospheres (Gustafsson et al. 2008). Like the MOOG/ATLAS analysis, we adopted a metallicity of [M/H] = −2.5 and [α/Fe] = +0.4 for all model atmospheres, and used the model atom presented in Gallagher et al. (2020) including Barklem et al. (2000) damping. After pre-computing a grid of NLTE curves of growth at many T eff and log g values, we use Delaunay triangulation and linear interpolation to find NLTE abundances as a function of stellar parameters and equivalent widths (using scipy.interpolate.LinearNDInterpolator).  . VLT spectra with best-fit models of the Ba line. From low to high wavelength the primary absorption features are due to Ca, Fe, and Ba. The Ba abundance measurement or limit is indicated. Stars are sorted by SNR, membership, and effective temperature. Top section: stars with SNR > 25. Bottom section: stars with SNR < 25. The sky residuals are clearly substantial for the lower SNR stars. The black line shows the data, while the grey shaded region indicates the 1σ spectrum noise. Best-fit models are shown in red and orange for clear and candidate members, respectively. The solid colored lines are the model including the sky line, while the dashed colored lines are the model fit in the absence of the fitted sky line residual. In some cases (e.g. star 020, 134, 157) the Ba line is completely degenerate with the sky emission line, so no useful constraint is obtained. The shaded colored regions indicate the 16-84th percentile range of model parameters. Note that we have not explicitly plotted the upper limit model.
In the rest of this paper, we adopt the NLTE [Ba/H] abundances as our fiducial abundances. However in Appendix C we show a comparison between the MOOG LTE/ATLAS and MULTI NLTE/MARCS abundances. Overall, the NLTE effects are approximately −0.3 dex, resulting in lower abundances compared to LTE modeling.

Abundance Summary
In summary, the main abundance results are Ba, Fe, Ca, and Mg measurements or upper limits. The VLT data are used to measure Ba, Fe, and Ca. The M2FS high-resolution data are used to measure Mg and to verify Fe and Ca. The M2FS medium-resolution data are used to verify Ba (Table 4). For the brightest stars, we use the M2FS HiRes data to derive abundances for Cr, Ti, and Nd, which we report in Appendix D. We do not discuss these elements more, as they are consistent with the discussion in J16 and R16. We adopt NLTE abundances for Ba and LTE abundances for other elements. Figure 4 shows the radial velocities of all stars in our sample (solid blue histogram). There is a clear peak at ≈65 km s −1 associated with Ret II. Stars considered clear Ret II members are shown as the solid red histogram, while all Ret II members including candidates are the open orange histogram (see Section 3.3 for details). For comparison, we show the radial velocity distribution of a smooth background halo from the Besançon model (Robin et al. 2003). We restrict the background to a CMD region surrounding our targets defined by four points (g − r, g) = (0.2, 21.2), (0.7, 21.2), (0.45, 16.0), (1.0, 16.0). We query a large area for statistics and rescale the distribution to the FLAMES field of view.

Radial Velocity Distribution
The residual background after removing the clear and candidate members is shown as a purple histogram. There continues to be an excess of stars near the velocity of Ret II relative to the Besançon model. Detailed investigation of the proper motions and spatial position shows that these residual stars cannot be members of Ret II. There may potentially be additional structure in this region of the sky that is not part of Ret II, though we do not see any clear spatial or proper motion trends.
We determine the mean velocityv and velocity dispersion σ v of Ret II with a Gaussian scatter model. Each star is assumed to have a true velocity distributed according to a Gaussian with mean and standard deviationv and σ v , and the actual observed velocity has a Gaussian noise added to it with the individual velocity uncertainty in Table 1. The prior on the mean velocity is uniform with no bounds, and the prior on the scatter is uniform in log space from σ v ∈ [0.1, 10] km s −1 . The posterior is sampled using Stan (Carpenter et al. 2017, following Casey et al. 2021). All five likely binary stars are removed. Using only the 27 clear likely single members of Ret II, we obtainv = 63.9 ± 0.5 km s −1 and σ v = 2.97 +0.43 −0.35 km s −1 . Adding the 8 candidate members givesv = 64.0 ± 0.5 km s −1 and σ v = 2.96 +0.44 −0.36 km s −1 , which is the same within the uncertainties. Our velocity dispersion is consistent both with previous measurements of ≈3.3 km s −1 (S15; W15; K15) and with the lower value of 2.8 +0.7 −1.2 inferred by Minor et al. (2019) who statistically accounted for binaries. Our velocity dispersion is about two times more precise given more stars and additional independent velocity measurements, but we caution that our velocity dispersion uncertainties may be overly optimistic given the simple treatment of possible systematics (see Appendix A).    Note-Rows where the third column has a prior are measured from this work.
is an extremely carbon-enhanced star (e.g., Koposov et al. 2018 Figure 6. The number of stars is low and the error bars are large, but the left two panels show that both [Mg/Fe] and [Ca/Fe] broadly decline with increasing [Fe/H]. Given the low statistics, there is no clear "αknee" as seen in the Milky Way and some dwarf galaxies (Venn et al. 2004;Tolstoy et al. 2009;Hill et al. 2019), though we note that most dwarf galaxies do not have such a clear knee (Kirby et al. 2020;Theler et al. 2020 Li et al. 2018), and sampled using Stan in a manner similar to the velocity dispersion. We adopt a log-uniform prior for the intrinsic scatter of 10 −2 − 10 0 . The results are a mean [Fe/H] = −2.64 with dispersion 0.32 dex, matching previous results (Simon et al. 2015;Walker et al. 2015;Koposov et al. 2015b). Adding Fe measurements from the two candidate members or including upper limits does not substantially affect these results.

Intrinsic Scatter in [Ba/H]
We detected [Ba/H] in 16 of our 32 member stars (21 out of 40 stars including candidate members). However, only [Ba/H] measurements from stars with SNR > 25 were considered reliable due to sky subtraction residuals (see Figure 3 and Appendices B and C), restricting the [Ba/H] measurements to 11 out of 16 member stars (13 out of 19 candidate members).  from 10 −2 − 10 0 , a flat prior from 0 to 1 for the mixing fraction, and no constraint on the means (except that one has to be larger than the other). For the abundance likelihoods, for stars with [Ba/H] detections we used the usual Gaussian likelihood using the [Ba/H] measurement and uncertainty as the mean and standard deviation for that star. For the stars with [Ba/H] 99% upper limits, we adopted a step-function likelihood where 99% of the probability is uniform between [Ba/H]= −5 to the measured upper limit, and 1% of the probability is uniform from the measured upper limit to [Ba/H]= +1. This model was implemented in Stan and sampled using the NUTS sampler with 4 chains and 10 6 steps. Model parameters were initialized near likely values from initial visual examination (i.e., mixing fraction based on the number of limits vs measurements, intrinsic scatters of 0.1 dex, component means of −4.0 and −1.5). Our The arrows indicate the 99% upper limits for stars with undetected Ba lines, which are included in our modeling procedure. The best-fit two-component models are shown as shaded regions (red for clear members, orange for candidate members).
Overall, the detected Ba abundances for the clear member stars are clearly unimodal, and the per-star [Ba/H] uncertainty can explain the observed scatter in [Ba/H] abundances. We obtain a 95% upper limit σ [Ba/H] < 0.20 dex. If adding the candidate members, we obtain a weaker upper limit σ [Ba/H] < 0.31 dex, which is primarily driven by star 14 that has a weak Ba line. As discussed previously, if this star is a member it is not likely that its stellar parameters can be determined with photometry, so the [Ba/H] inference for this star is likely inaccurate.
As a note of caution, we believe the choices made for our fiducial measurement are the most appropriate for quantifying the [Ba/H] intrinsic scatter in Ret II, as they minimize the impact from known systematic issues. But for completeness, in Appendix C we show the effect of all choices (i.e., SNR cut, membership, ν t − log g relation) on the intrinsic scatter measurement. Changing the log g −ν t relation changes the limit to σ [Ba/H] < 0.28 dex for clear members and σ [Ba/H] < 0.37 dex for candidate members. This effect is driven almost entirely by the two cool stars (ID 1 and 97), which have low statistical uncertainty due to their brightness but have a strong dependence on microturbulence. with radius for confirmed and candidate members. Overall, there are no strong abundance gradients for Fe and Ba. This lack of trend within this radius is expected, as any initial abundance gradients produced when Ret II was forming its stars at z 6 will likely be dynamically mixed away by z = 0: for 10 4 stars within the half-light radius of 58 pc and typical velocity of 3 km s −1 , the two-body relaxation time is about 2.5 Gyr (Binney & Tremaine 2008). Note a metallicity gradient could be present at larger radii, as illustrated by the clear gradient in Tucana II (Chiti et al. 2021), but our observations only reach two half-light radii.

Abundance Trends with Radius
For Mg and Ca, there may be a small abundance gradient in [Mg/Fe] and [Ca/Fe], perhaps decreasing slightly at larger distances. For Mg, this could easily be due to small number statistics rather than an abun-  The symbols are as in Figure 5, with large red points indicating high-SNR stars with confident membership. There are no significant abundance gradients. dance gradient; and for Ca, the scatter around a mean trend would be comparable to or larger than the size of the gradient itself. Thus, we do not consider there to be evidence for any abundance gradients out to two half-light radii.

DISCUSSION
We have carefully determined the [Ba/H] distribution in Ret II. Figures 5 and 7 and Table 6 show that ∼30% of the stars in Ret II are relatively metal-poor with no detected r-process elements, while ∼70% of the stars in Ret II have identical r-process abundances of [Ba/H] = −1.7, precise to 0.2 dex. We now consider the implications of this measurement for galaxy formation and metal mixing (Section 5.1), the r-process site (Section 5.2), and chemical evolution in Ret II (Section 5.3). In most of the discussion we assume that the high [Ba/H] plateau in Ret II originates predominantly from a single r-process event. However, we consider (and dismiss) the possibility of Ba contamination from multiple r-process sources or non-r-process sources in Section 5.4.
Previous studies indicated that the r-process material in Ret II likely originates from a single r-process event, rather than multiple r-process events occurring within Ret II (Ji et al. 2016a;Roederer et al. 2016b). To update those conclusions with more recent results, out of 20 UFDs with high-resolution chemical abundances measured, Ret II has r-process abundances 2-3 orders of magnitude higher than the neutron-capture element abundances in other UFDs. It is thus extremely unlikely that Ret II would be enriched by multiple prolific r-process events, while other UFDs have no similarly prolific events. To estimate this probability, we ran Monte Carlo simulations assuming that the number of r-process events in each of 20 galaxies is distributed as a Poisson random variable. We assume each UFD is equally likely to host an r-process event and run simulations where the intrinsic rate of r-process events is between 0.00 to 0.50 events per galaxy. We then conservatively take the highest probability. The probability that any one galaxy gets ≥ 2 r-process events while every other galaxy gets 0 r-process events is < 1.5%. Two other ultra-faint satellites also appear to display a lower level of r-process enhancement (Tucana III and Grus II; Hansen et al. 2017Hansen et al. , 2020Marshall et al. 2018). If we assume that both of these systems are dwarf galaxies and hosted an r-process event, the probability that three UFDs have at least one r-process event and one has at least two r-process events is < 6%. However, the classification of Tucana III and Grus II as dwarfs is currently uncertain ; if neither is a dwarf galaxy then the probability of multiple r-process events in Ret II is again 1.5%. Thus, it is unlikely that Ret II was enriched by more than one r-process event.

Well-mixed [Ba/H] Implies Bursty Star Formation in Ret II
We first interpret the distribution of detected [Ba/H], which has a low intrinsic dispersion σ Ba < 0.20. Because all this r-process material is deposited at a single time in Ret II, the observed variation in [Ba/H] is entirely due to variations in metal mixing. The unresolved Ba dispersion thus implies that the r-process material in Ret II must have been very well-mixed by the time the high-Ba stars in Ret II have formed. To interpret this homogeneous Ba, consider a parcel of gas with a fresh source of metals. There are two crucial timescales: τ mix , the time for the metals to completely mix within that gas, and τ sf , the time to turn that parcel of gas into stars. Stars forming from this gas will be chemically homogeneous if mixing is faster than star formation, i.e. τ mix < τ sf . The homogeneous Ba abundances in Ret II thus provide a lower limit on the time between when the r-process event occurs and when the next generation of stars forms. In other words, the mixing timescale provides a constraint on the burstiness of star formation. We argue in this section that both analytic estimates and simulations have shown that τ sf > τ mix > 100 Myr, and thus there must be at least a 100 Myr gap in the star formation history of Ret II. This is the first observational evidence of bursty star formation in an ultra-faint dwarf galaxy.
We first describe some basic physical ingredients determining τ mix . Metal mixing in dwarf galaxies can be approximated as proceeding in two phases: the initial explosion remnant and subsequent turbulent mixing (e.g., Karlsson et al. 2008;Emerick et al. 2020). For an individual explosion, the initial remnant is dominated by a momentum-driven snowplow, lasts only ∼10 5 yr, and sweeps up ∼10 5 M of gas (for a 10 51 erg explosion; e.g., Cioffi et al. 1988;Ryan et al. 1996;Greif et al. 2007;Ji et al. 2015;Macias & Ramirez-Ruiz 2019;Magg et al. 2020). This mass is insignificant compared to a dwarf galaxy's ISM. The dominant process determining τ mix is thus turbulent mixing, where large-scale energy injection cascades to small-scale velocity eddies that enable microscopic diffusion (e.g., Pan et al. 2013). Turbulent mixing is typically modeled as a diffusion process (Klessen & Lin 2003;Karlsson et al. 2008;Greif et al. 2010;Ji et al. 2015;Krumholz & Ting 2018;Beniamini & Hotokezaka 2020;Tarumi et al. 2020), R 2 t ∝ D t τ , where R t is the turbulent diffusion distance, D t is the turbulent diffusion coefficient, and τ is the time since material is deposited. Drawing an analogy to mixing length theory, the diffusion coefficient can be estimated by a typical length and velocity scale driving turbulence, D t ∼ R turb v rms . Complete mixing occurs when R t reaches a length scale associated with the full size of the galaxy, R t = R gal . Thus, τ mix ∼ R 2 gal /(R turb v rms ) (Pan et al. 2013). In early dwarf galaxies, turbulence is primarily driven by gravitational gas accretion or mergers (Wise & Abel 2007;Greif et al. 2008;Klessen & Hennebelle 2010;Safranek-Shrader et al. 2012;Ritter et al. 2015), which can be used to estimate a turbulent diffusion coefficient (Karlsson et al. 2008;Ji et al. 2015). Putting in typical values for a dwarf galaxy (R gal ∼ 500 pc the extent of the star forming gas, R turb ∼ 100 pc the size of a typical explosion remnant, and v rms ∼ 10 km s −1 the turbulent velocity of ambient gas, from Tarumi et al. 2020) gives τ mix ∼ 250 Myr.
While intuitive, the mixing length formalism fails to describe the full physics of the complex, anisotropic, and multi-phase metal-mixing process in dwarf galaxies. For example, it is well known that mixing depends on the temperature of the ISM phase that the metals reside in, such that hot ISM phases mix much more efficiently than cold ISM phases (e.g., Kobulnicky & Skillman 1997;de Avillez & Mac Low 2002;Emerick et al. 2018Emerick et al. , 2019Emerick et al. , 2020; and the anisotropic topology of cold gas clumps and filaments in early dwarf galaxies affects where metals get deposited and new stars form (e.g., Webster et al. 2016;Chiaki et al. 2018;Magg et al. 2022). It is thus crucial to study metal mixing with hydrodynamic galaxy formation simulations. A few recent simulations have explicitly studied metal mixing in dwarf galaxies (e.g., Webster et al. 2014Webster et al. , 2015Hirai et al. 2015Hirai et al. , 2017Revaz et al. 2016;Escala et al. 2018;Emerick et al. 2019Emerick et al. , 2020Tarumi et al. 2020;Jeon et al. 2021). In the vast majority of simulations, abundance scatter from individual metal sources is typically very large, ranging from 0.4-2.0 dex (e.g., Safarzadeh & Scannapieco 2017;Emerick et al. 2020;Applebaum et al. 2021) 9 . Note that many simulations compare their simulation abundance scatter directly to the observed abundance data, without deconvolving the abundance uncertainties. Tarumi et al. (2020) performed a direct simulation of r-process enrichment in Ret II, simulating a UFD and injecting r-process elements from a single NSM to see the resulting r-process abundance spread. In order to homogenize the gas to a level consistent with that observed in Ret II, they found that the gas had to mix for a time period of a few hundred Myr before forming stars. They measure an effective diffusion coefficient in their simulation of D t ≈ 10 −3 kpc 2 Myr −1 , resulting in a complete mixing timescale of ≈250 Myr. Thus, we expect 100 Myr < τ mix < τ sf , and the timescale between bursts of star formation in Ret II should be over 100 Myr, a significant fraction of a Hubble time at z > 6.
There are a few important caveats to the Tarumi et al. (2020) simulation. First, the mixing time can be affected by stochastic events. When Tarumi et al. (2020) exploded the NSM at the outskirts of the galaxy rather than the center (to mimic a velocity kick, also see Safarzadeh & Scannapieco 2017;Safarzadeh et al. 2019b;Bonetti et al. 2019), the lower diffusion coefficient resulted in less efficient mixing of the r-process elements. Emerick et al. (2019Emerick et al. ( , 2020 also emphasize that the exact timing and location of r-process production relative to other stellar feedback sources that drive turbulence can both increase and decrease the mixing time. Second, the abundances in the Tarumi Tarumi et al. (2020) indeed have a gas-rich merger that helps homogenize Ba and Fe but causes a dispersion in [X/H] at fixed time (Y. Tarumi, private communication). Finally, most galaxy formation simulations stay at relatively low resolutions, e.g. Tarumi et al. (2020) adopt the ISM model from the Auriga Project that uses an effective equation of state model below 0.1 cm −3 (Grand et al. 2017). This may resolve mixing in the large-scale ISM, but it may not resolve small-scale inhomogeneous mixing (e.g., Pan et al. 2013;Chiaki et al. 2018;Magg et al. 2022). Tarumi et al. (2020) is the only simulation directly comparable to Ret II, but the overall timescale of > 100 Myr to mix matches results from other simulations (Hirai et al. 2017;Emerick et al. 2019) and estimates based on the mixing length scaling relation (Karlsson 2005;Pan et al. 2013;Ji et al. 2015). Thus, the mixing timescale of > 100 Myr seems robust.
In summary, the mixing time in UFDs is likely larger than 100 million years. The homogeneous r-process abundances in Ret II thus indicate that at least two early bursts of star formation occurred in Ret II, separated by at least a few hundred million years. The first burst produced r-process elements that enriched stars born in the second burst (which should be relatively extended due to the presence of Type Ia enrichment, see Section 5.3). Given that we find ≈70% of Ret II stars are r-process enhanced, a concrete prediction is that star formation histories of Ret II with a precision of 100Myr should show 30% of stars forming first, a gap of > 100 Myr, and then the other 70% of star formation. Resolving bursty star formation on 100 Myr timescales is currently out of reach of direct star formation history measurements (e.g., Brown et al. 2014;Weisz et al. 2014a), but it qualitatively matches predictions from several UFD simulations (e.g., Wheeler et al. 2019;Jeon et al. 2017). These simulations achieve burstiness due to strong clustered supernova feedback, which purges the galaxy of starforming gas and requires a long recovery time for gas to cool and form the next generation of stars (e.g., Jeon et al. 2014). The r-process elements can thus homogenize during this extended recovery time. The mass ratio M Eu /M r is 10 −3.0 , assuming M r is elements with A ≥ 80 from the solar r-process pattern (Sneden et al. 2008;Côté et al. 2018;Ji et al. 2019b). Thus, we find that M r /M H ≈ 10 −7.2±0.1 , where M H is an effective dilution mass of hydrogen.

A Prompt, High-Yield R-Process Site
Inferring the yield M r of the r-process site depends on what is assumed for M H . Ji et al. (2016a) previously argued that expected dilution masses range from 10 5 − 10 7 M . The lower end is set by the initial explosion remnant (also see simulations by Magg et al. 2022), and the upper end is set by the total available gas in a 10 8 M dark matter halo that is likely to host Ret II at high redshift. This would correspond to M r ∼ 10 −2.2 − 10 −0.2 M . However, simulations show that the effective dilution masses are near the minimum of 10 5 M only when metals are inhomogeneously mixed (e.g., Chiaki et al. 2018;Jeon et al. 2021;Magg et al. 2022). The homogeneity of r-process elements in Ret II thus suggests that the metals have diluted into a larger hydrogen mass, which should be in the range 10 6 − 10 7 M . The higher dilution mass would imply a higher r-process yield of M r ∼ 10 −1.2 − 10 −0.2 M .
An alternate way to estimate the r-process yield is to count up the amount of r-process elements locked into stars, then apply a correction factor for how much r-process material would not be captured in stars. A large fraction of r-process can be lost to the intergalactic medium due to the low gravitational potential of early dwarf galaxies (e.g., Beniamini et al. 2017;Brauer et al. 2021). Both empirically and theoretically, only 10 −2 of metals are retained in dwarf galaxies (Dekel & Woo 2003;Robertson et al. 2005;Kirby et al. 2011;McQuinn et al. 2015). We can estimate the total mass of r-process in Ret II using its present-day stellar mass of ≈3300M (Mutlu-Pakdil et al. 2018). Assuming a hydrogen mass fraction of 0.75, the total r-process mass contained in Ret II today is 10 −3.8 M . Thus the expected yield of the r-process site should be M r 10 −1.8 M , consistent with our previous estimate M r ∼ 10 −1.2 − 10 −0.2 M . Note that the dilution masses described in Ji et al. (2016a) implicitly include this metal loss to the IGM, as the higher effective dilution masses can be thought of as corresponding to lower metal retention (also see Fig  11 of Magg et al. 2022).
The homogeneous [Ba/H] also suggests that the rprocess site has to be fairly prompt in Ret II. Ji et al. (2016c) originally argued that the recovery times in UFDs (i.e., the time for gas to recollapse to the center of a halo after feedback) were longer than 10 − 100 Myr (Jeon et al. 2014;Ji et al. 2015;Bland-Hawthorn et al. 2015). This allowed a significant fraction of ordinary neutron star mergers with 10 − 100 Myr delay times to merge and enrich the gas before star formation. However, if the r-process material must then subsequently mix for an additional 100 Myr to homogenize, this puts strong pressure on the r-process site to be very prompt in order to mix fully before turning into stars. Additionally, Simon et al. (2022) recently obtained a precise star formation history for Ret II using Hubble colormagnitude diagrams. Combined with the result from this paper that 70% of Ret II stars are r-process enhanced, they concluded that the r-process time delay in Ret II must be shorter than 500 Myr.
Together, the higher r-process yield and more prompt r-process event implied by homogeneous r-process mixing slightly favor rare core-collapse supernovae over neutron star mergers as the source of r-process elements in Ret II. Our higher expected r-process yield of 10 −1.2 −10 −0.2 M is a better match to the 0.08−0.3M of r-process produced in collapsar disk winds (Fryer et al. 2006;Surman et al. 2006;Siegel et al. 2019;Miller et al. 2020), but also consistent with magnetorotation-ally driven jets (10 −2.5 − 10 −1.5 M of r-process, Mösta et al. 2018) or neutron star mergers (10 −3 − 10 −1 M , Wu et al. 2016;Radice et al. 2018; note GW170817 had r-process mass ≈10 −1.5±0.3 M , Drout et al. 2017;Kilpatrick et al. 2017;Tanaka et al. 2017;Tanvir et al. 2017;Chornock et al. 2017), or common envelope jet supernovae (e.g., Grichener & Soker 2019;Grichener et al. 2022). The fact that there is a few hundred Myr delay after r-process enrichment also favors core-collapse supernovae. However, very prompt and high-yield neutron star mergers are still on the table (Beniamini et al. 2016;Beniamini & Piran 2019;Safarzadeh et al. 2019a). Figure 5 shows that the [Ba/H] abundance of the rprocess rich stars stays very flat over an extended range of [Fe/H]. This can be seen quantitatively by comparing the metallicity (Fe) dispersion of 0.32 +0.10 −0.07 dex to the r-process (Ba) dispersion of < 0.20 dex. The simplest interpretation of the larger Fe dispersion is that the rprocess stars formed over some extended period of time where Ret II was able to self-enrich with iron from supernovae, as expected for a dwarf galaxy (Willman & Strader 2012). The flat [Ba/H] abundance would then clearly indicate that there is no pristine gas accretion nor any significant r-process production during the last 70% of Ret II's stellar mass growth. If there were significant pristine gas accretion during this time, it would reduce [Ba/H] at high [Fe/H] 10 .

No Gas Accretion During Most Ret II Star Formation
This gas cutoff scenario also could explain the [Mg/Ca] trend seen in Figure 6 through the integrated galactic initial mass function (IGIMF) (Weidner et al. 2013). In this model, a gas-poor galaxy is unable to create the densest and largest molecular clouds, introducing an effective upper mass limit to stars formed. Since Mg is predominantly produced in the most massive core-collapse supernovae and Ca is produced in all supernovae, a restricted gas supply will result in lower [Mg/Ca] abundances relative to a fully sampled IMF (McWilliam et al. 2013;Ji et al. 2020a;Lacchin et al. 2020). Thus, a lack of gas accretion could explain both the flat [Ba/H] and the declining [Mg/Ca] observed in Ret II. This observation should simplify chemical evolution models aiming to reproduce the r-process abundance trends of Ret II (e.g., Komiya & Shigeyama 2016;Ojima et al. 2018;Molero et al. 2021;Cavallo et al. 2021).
A lack of gas accretion may indicate something about the broader formation environment of Ret II. In particular, it is expected that UFDs like Ret II are ultimately quenched by reionization (e.g., Bullock et al. 2000;Benson et al. 2002;Brown et al. 2014;Rodriguez Wimberly et al. 2019), but it is not yet clear whether reionization immediately removes cold gas from halos or just restricts gas inflow (e.g., Okamoto et al. 2008;Weisz et al. 2014b;Jeon et al. 2017;Bose et al. 2018;Wheeler et al. 2019). Since >70% of Ret II stars form in the absence of significant gas accretion, it may be that it formed all these stars after reionization.
In this vein, it is interesting to note that Ret II is a satellite of the Large Magellanic Cloud (Patel et al. 2020;Erkal & Belokurov 2020;Battaglia et al. 2022). Sacchi et al. (2021) tentatively find that the star formation histories of LMC UFD satellites (including Ret II) take longer to complete the last 10% of their star formation history compared to Milky Way UFD satellites. This could suggest that LMC UFD satellites like Ret II were relatively isolated when they formed compared to Milky Way UFD satellites, and thus experienced delayed reionization. If so, it would support the concept of patchy reionization at the smallest galactic scales (e.g., Lunnan et al. 2012;Aubert et al. 2018). More recently, Simon et al. (2022) have determined the star formation history of Ret II with additional HST photometry. They find that Ret II likely took ∼2.5 Gyr to form all its stars, largely corroborating the arguments from chemical evolution here.
For completeness, we note that the above discussion has implicitly assumed that τ mix and τ sf are the same for both Fe and Ba. One can imagine scenarios where the timing of stellar feedback causes an early source of Ba to be mixed more than a later source of Fe (e.g., Ritter et al. 2015;Schönrich & Weinberg 2019). In this case, stars at all metallicities would form simultaneously. This scenario seems unlikely for Ret II given the coherent evolution of Mg and Ca with [Fe/H], but it provides motivation to obtain more precise star formation histories in Ret II.

Contamination of Ba by Other Sources
We have interpreted our Ba measurements in Ret II as tracing pure r-process, based on the pure r-process patterns found in Ji et al. (2016a) and Roederer et al. (2016b). However in principle there could be three possible contaminating sources of Ba that are empirically found in UFDs: (1) a low-yield source of Ba observed in most UFDs, of unknown origin (Frebel & Norris 2015;Roederer 2017;Ji et al. 2019c), but possibly attributed to r-process in neutrino driven winds (J16, Simon 2019) or s-process in rotating massive stars (Frischknecht et al. 2016;Limongi & Chieffi 2018;Tarumi et al. 2021); (2) late-time AGB enrichment in the ISM Ji et al. 2020a); or (3) mass transfer of s-process Ba from a binary companion (Frebel et al. 2014).
None of these possible contaminants will impact our conclusions. The impact of the first two sources is much less than the r-process content of Ret II, and it can be estimated by considering Ba abundances in other UFDs. The low-yield Ba source produces typical [Ba/H] ∼ −4 (Ji et al. 2019c). Ba in the ISM from AGB stars is not often seen in UFDs given their short star formation durations, but where it is seen it reaches [Ba/H] LTE ∼ −2.5 Ji et al. 2020a). In both cases, the amount of contamination is at most 1/10 of the Ba in Ret II, too low to make a significant perturbation. For the third source, AGB mass transfer tends to produce [Ba/Fe] ∼ +2, much more Ba than is observed in these stars (e.g., Frebel et al. 2014;Hansen et al. 2016). Thus, we find it extremely unlikely that Ba is tracing anything other than the single r-process event in Ret II.
However, we note that one candidate member, Star 97, is quite red both in DES and Gaia photometry. It is possible this is due to large amounts of carbon on the surface of the star, in which case it may have experienced mass transfer possibly including Ba. If so, this further justifies excluding star 97 from our main results.

CONCLUSION
We have obtained multi-object spectroscopy of red giant branch members in the ultra-faint dwarf galaxy Reticulum II using VLT/GIRAFFE and Magellan/M2FS. Our redetermination of the velocity and metallicity dispersion is consistent with past results, and we detect no significant spatial gradients in the element abundances.
Ret II is of special interest due to its enrichment by a single r-process event. Our primary new result is a quantitative measurement of the [Ba/H] distribution (Figure 7, Table 6), which is a unique probe of gas dynamics and metal mixing within a faint, currently gas-free dwarf galaxy. Approximately 30% of Ret II stars have no detected r-process material, while the other 70% are enriched to a high enhancement. We place an upper limit of σ [Ba/H] < 0.20 dex on the intrinsic [Ba/H] dispersion of the high-Ba stars, which implies that the initial r-process enrichment needs to turbulently mix and homogenize for at least 100 Myr before stars form. This is the first direct evidence of bursty star formation in a UFD. The long mixing time also favors an r-process site that is very prompt and produces a high r-process yield ( 10 −1.5 M ). We thus slightly favor rare corecollapse supernovae as the source of r-process elements in this galaxy due to their higher r-process yield, though prompt high-yield neutron star mergers are allowed as well.
Examining the chemical evolution in Ret II, we find an overall declining [α/Fe] vs [Fe/H] pattern as expected in dwarf galaxies. Since [Ba/H] is flat over an extended [Fe/H] range, this suggests that Ret II did not accrete significant gas during the last 70% of its star formation. This is consistent with the observed declining [Mg/Ca] ratio if Ret II was too gas-poor to form the most massive core-collapse supernovae. The chemical evolution of Ret II thus suggests that it may have formed in an underdense environment, consistent with its status as a satellite of the Large Magellanic Cloud.
These constraints on UFD formation and the r-process site demonstrate the power of dwarf galaxy archaeology. By finding stars in a common formation environment, it becomes possible to ask questions that could not be answered if these same stars were found individually scattered through the Milky Way. Reticulum II is unusually nearby and thus currently accessible for this type of study, but as the next generation of extremely large telescopes comes online, it will become possible to extend similar techniques to study the chemistry of ultra-faint dwarf galaxies throughout the Milky Way (Ji et al. 2019a This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France (Wenger et al. 2000). This research has made use of NASA's Astrophysics Data System Bibliographic Services. This work has made use of the VALD database, operated at Uppsala University, the Institute of Astronomy RAS in Moscow, and the University of Vienna.
This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www. cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www. cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research uses services or data provided by the Astro Data Lab at NSF's National As a relatively nearby UFD of great scientific interest, Ret II has obtained many different epochs of radial velocities. We have collected all currently available literature velocities in Table 7. The literature velocities are mostly derived from coadding spectra taken across 1-4 adjacent nights, so the MJD reported here is only accurate to 2 days of precision. These velocities are not homogeneous and may suffer from systematic effects.
In a first attempt to calibrate the systematic effects, we adopt the S15 velocities as a reference velocity scale, as they have the most stars with common velocities compared to other literature sources. For matched stars in each sample, we calculate a weighted mean velocity offset. After removing this offset, stars with velocity variations inconsistent with a chi-squared test with p < 0.01 are identified as likely binaries (e.g., Chiti et al. 2022). We iterate this process until convergence, resulting in five binary stars: Star 2, 13, 18, 19, and 21. The final mean velocity offsets relative to S15 are 1.01 km s −1 for K15, 0.63 km s −1 for J16, 1.02 km s −1 for Roederer et al. (2016b), 2.79 km s −1 for our VLT spectra, and 2.37 km s −1 for our HighRes M2FS spectra. Our MedRes M2FS data have a very large offset of −8.9 km s −1 , and we thus decided to exclude it from any velocity studies.

B. EFFECT OF MICROTURBULENCE ON BARIUM ABUNDANCES
Microturbulence (ν t ) is a parameter introduced to 1D stellar atmospheres to account for unmodeled 3D atmospheric effects. It affects lines at the saturated part of the curve of growth, where a higher microturbulence effectively desaturates strong lines by adding some extra doppler broadening (e.g., Gray 2008). When lots of Fe I lines are available, microturbulence is usually found by balancing Fe I abundance as a function of line strength. Empirical measurements of microturbulence show that red giants with lower surface gravities (and temperatures) tend to have higher microturbulence values (e.g., Barklem et al. 2005).
In our VLT spectra, the Ba 6496Å line is at or near saturation when detected. For our coolest giants, increasing ν t by 0.2 km s −1 reduces [Ba/H] by 0.15 dex. Since we aim to resolve abundance scatter on the order of 0.20 dex, this systematic effect is often the dominant uncertainty, especially for the brightest giants where the equivalent width is well-measured. There are not enough Fe lines to self-consistently measure ν t in our stars, so we must use existing correlations between log g and ν t . Here, we investigate five different datasets that measured ν t using high-resolution spectroscopy of metal-poor red giant stars, examining systematic differences in ν t −log g relations, as well as the scatter around those relations. The effect of these choices on our barium abundance results is investigated in Appendix C. Figure 9 shows the result of our investigation. The left column in Figure 9 plots log g vs ν t for data from Barklem et al. (2005) (R14g). When available, we use the stellar parameters as tabulated in JINAbase (Abohalima & Frebel 2018). Horizontal branch stars have higher microturbulence than RGB stars, so they are removed with a cut log g > 0.00286 T eff − 12.7. To each dataset, we fit a linear and quadratic polynomial for ν t as a function of log g, using a robust fitter based on the routine robust poly fit in the AstroIDL library. The scatter around the fit is measured with the biweight scale of the residuals (a robust standard deviation). The coefficients and scatter around each relation are given in Table 8.
Note that the right column of Figure 9 shows that our stars with SNR < 25 have a much larger scatter than stars above that threshold. Visual examination of the spectra (Figure 3) suggests that these stars are adversely affected by inaccurate sky subtraction that is blended with the Ba line. We thus have decided to exclude the low-SNR stars from most analyses. Also note that we used MOOG LTE/ATLAS [Ba/H] abundances for this figure, though the conclusions are robust if using NLTE instead.

B.1. Mean Relations
There are clear differences in the average ν t − log g relation across different literature samples: at low log g, the first three rows (B05, C13, J15) have systematically higher ν t than the last three rows (M08, R14, R14g). The origin of this difference is not clear. One possibility is the spectra for M08 and R14 have SNR ∼ 100, substantially higher than B05 and J15 which typically have SNR ∼ 30. The low SNR of weak iron lines could bias microturbulence too high (as 0.0175 log g 2 + −0.3242 log g + 2.009 0.08 −0.2545 log g + 1.944 0.08 R14 0.0471 log g 2 + −0.3474 log g + 1.969 0.18 −0.1201 log g + 1.764 0.18 R14 giants 0.0386 log g 2 + −0.3313 log g + 1.960 0.13 −0.2247 log g + 1.897 0.10 our final inferred Ba scatter. Figure 10 shows this specifically for the high SNR and clear member stars for the R14g and B05 relations. The R14g ν t − log g relation clearly has less trend with T eff than B05. Note that the trends are primarily driven by the two coolest and brightest stars (T eff ∼ 4500 K). Because these two stars have low statistical uncertainty on their [Ba/H] abundances, the systematic effect of microturbulence can make a large difference in the inferred [Ba/H] scatter.
We have decided to adopt the quadratic ν t − log g relation from the giant stars in Roederer et al. (2014) (R14g) as our fiducial results. R14g have the highest SNR spectra of metal-poor giants with the largest wavelength coverage out of all these data samples. The results of this paper would not change if we used the very similar M08 or R14 relations instead. However, this ν t − log g relation is different from most previous studies of dwarf galaxy stellar abundances (Ji et al. 2016b;Ji & Frebel 2018;Roederer et al. 2016b). Thus in Appendix C, we give all results using the ν t − log g relation from B05 as well, which matches those previous abundance studies. We note that adding NLTE effects for Ba exacerbates the trend for [Ba/H] vs T eff when using the B05 ν t -log g relation, because the NLTE corrections are larger (more negative) when microturbulence is higher. It is also worth noting that past studies of Ret II had clear [Ba/H] trends with temperature that can be explained by microturbulence.

B.2. Microturbulence Scatter
Typically, a systematic uncertainty of ≈0.2 km s −1 is adopted for microturbulence, which accounts for the systematic mean differences described above. However, because we are interested in the abundance scatter within Ret II, another crucial value is the intrinsic scatter in microturbulence around a "true" ν t − log g relation, i.e., changes in the atmospheric structure that are unmodeled by log g alone. This error is likely smaller than the observational scatter, because the microturbulence measurements themselves are noisy.
Examining our five datasets, two have a scatter of ∼0.2 km s −1 (C13, J15), while three have a scatter of ∼0.1 km s −1 (B05, M08, R14). Using a smaller intrinsic ν t scatter would increase the significance of our scatter detections, while a larger intrinsic ν t scatter reduces the significance. Since we have adopted the mean ν t − log g relation from R14g, we also decide to adopt the intrinsic scatter of 0.13 km s −1 from that data sample. In our systematic investigations using the B05 sample, we adopt the corresponding intrinsic scatter of 0.12 km s −1 .

C. SYSTEMATIC EFFECTS ON BARIUM SCATTER
Here we explore the effect of different data subsets and microturbulence relations on the main result of this paper, the [Ba/H] mean and scatter. For the data samples, we consider permutations of membership (clear members only vs including candidate members) and MULTI NLTE/MARCS vs MOOG LTE/ATLAS. For the microturbulence relations, we use the fiducial R14 giants (R14g) relation, as well as the B05 relation that has a higher microturbulence for the coolest/lowest gravity giants. For each of these permutations, we fit the two-component Ba scatter model described in Section 4.4. Note that while our fiducial model is run with a very large number of steps, for the other models we only sampled to reach 100 effective samples, and thus the uncertainties and limits on the parameters will be less accurate. Table 9 gives the results of the model fits. The first row is our fiducial value, while the other rows show various data permutations. µ 1 and σ 1 are the most important values, indicating the mean and intrinsic spread on the detected [Ba/H] abundances. µ 2 and σ 2 are the mean and scatter of the undetected [Ba/H] component, which is not wellconstrained given that no low [Ba/H] abundances were detected. p 2 is the fraction of stars in the undetected [Ba/H] component, i.e., 1 − p 2 is the fraction of r-enhanced stars. The uncertainties are 1σ, and the limit on σ 1 is a 95% limit.  We point out three main conclusions of Table 9. First, the MULTI NLTE/MARCS mean [Ba/H] abundances (µ 1 ) are typically lower than the MOOG LTE/ATLAS abundances by about 0.3 dex. Figure 11 shows that the typical [Ba/H] correction going from MOOG to MULTI is −0.27 ± 0.04 dex in a way that is fairly close to a constant offset. This is a result both of the different model atmospheres as well as the effect of NLTE. Second, when considering just the clear member stars, none of the models detect a significant [Ba/H] dispersion σ 1 . The constraint is stronger when using NLTE, but weaker when using the B05 microturbulence relation instead of the R14g microturbulence relation. This is driven primarily by the coolest Ret II star (ID 1), which is most affected by the different microturbulence relations (Figure 9). Third, when including the candidate member stars, all the upper limits on σ 1 get looser, and actually in one case (B05 MOOG with candidates) the intrinsic dispersion is resolved at 2σ. This is predominantly because of the outlier star 14, which has a weak but detected Ba line. If this star is actually part of Ret II, then our two-component model for [Ba/H] is likely insufficient to describe the data because star 14 is well outside of the main peak of [Ba/H] detections, but well above the more stringent [Ba/H] upper limits.
After examining Table 9 and Figure 9, we decided that using the R14g MULTI/NLTE results with only clear members is the most reliable measurement. It is clear that using NLTE and definite members will result in a better measurement, and eliminating the trend with stellar parameters discussed in Appendix B justifies using R14g instead of B05. However, for completeness, we show several permutations of best-fit [Ba/H] distributions in Figure 12. The top-left panel shows the R14g and NLTE abundances used in the main paper, but plotting all low-SNR detections as small data points and low-SNR upper limits as grey arrows. As also seen in the right column of Figure 9, the low-SNR data are skewed towards higher [Ba/H] abundances primarily due to bad sky subtraction (Figure 3). The top-right panel shows three alternate fits to different permutations of data used (members and candidates; high-and low-SNR data). The bottom left panel shows the effect of changing the radiative transfer, and the bottom-right panel shows the effect of changing the microturbulence relation. These differences make a relatively small change to the Ba dispersion (which is not resolved) but a fairly large change to the mean abundance.