© The Authors 2025

1. Introduction

Dwarf galaxies have been extensively observed in the Local Group, where more than 70 objects have been detected as satellites orbiting around the Milky Way (MW) and Andromeda (McConnachie 2012). Among these satellites, the ultra-faint dwarf galaxies (UFDs; M ≲ 10⁵ M_⊙, Simon 2019) are dominated by ancient, > 10 Gyr (Brown et al. 2012); Gallart:2021, and very metal-poor stars (Kirby et al. 2008); Spite:2018 and, similarly to the classical dwarf spheroidal galaxies, have little to no gas mass (e.g., Westmeier et al. 2015). The observed properties of UFDs, which include the metallicity distribution function and the fraction of carbon-enhanced metal-poor stars (CEMP-no, e.g., Lai et al. 2011; Yoon et al. 2019) strongly support the idea that UFDs are living relics of the first star-forming mini-halos (e.g., Ricotti & Gnedin 2005; Salvadori & Ferrara 2009; Jeon et al. 2017; Barnes & Hut 1986), which formed before the end of reionization and are hosted by ∼10⁶ − 10⁸ M_⊙ dark matter halos. Their early formation and low mass make them the best candidates to host the first stars’ descendants (e.g., Salvadori et al. 2015; Magg et al. 2018; Rossi et al. 2024) and perfect laboratories in which to study the effects of feedback processes on galaxy formation (e.g., Jeon et al. 2015; Agertz et al. 2020; Gutcke et al. 2022). Moreover, UFDs are the most dark-matter-dominated stellar systems known, with mass-to-luminosity ratios of M/L > 100 M_⊙/L_⊙, and studying them is expected to play a key role in helping us understand the nature of dark matter (e.g., Bullock & Boylan-Kolchin 2017; Strigari 2018; Battaglia & Nipoti 2022).

Identifying UFD stars becomes increasingly difficult as the distance from the galactic center increases due to their faint nature and the foreground contamination of MW stars. Nonetheless, recent observations of UFDs have started to provide evidence for stars at large distances from the center. This is the case for Tucana II, where stars have been identified up to a half-light radius of ∼9 R_1/2 (Chiti et al. 2021), Hercules I, with stars detected up to ∼9.5 R_1/2 (Garling et al. 2018; Longeard et al. 2023; Ou et al. 2024), Tucana V, up to ∼10 R_1/2 (Hansen et al. 2024), Ursa Major I up to ∼4 R_1/2, Boötes I up to ∼4 R_1/2, and Coma Berenice up to ∼2 R_1/2 (Waller et al. 2023). In addition, multiple studies identify new and distant stars exploiting Gaia data. Tau et al. (2024) explore the outskirts of 30 UFDs using already known and new RR Lyrae stars and show that at least ∼30% present extended stellar populations. By considering ∼60 MW dwarf galaxies’ satellites, Jensen et al. (2024) systematically identify candidate member stars based on spatial, color-magnitude, and proper motion information; they report an extended secondary outer profile for nine dwarf MW satellites, seven of which are UFDs: Boötes I, Boötes III, Draco II, Grus II, Segue I, Tucana II, and Tucana III.

The formation channels for such extended features can be i) tidal stripping, ii) in situ formation via supernova (SN) feedback, and iii) merging events. In proximity to the MW, its gravitational pull can deform the UFD’s morphology, forming an extended halo, as is shown by Fattahi et al. (2018). Signatures of tidal interaction are signaled by the presence of velocity gradients (Li et al. 2018) and deviations from the expected exponential profile that typically describes the stellar density distribution (Peñarrubia et al. 2008). Early SN feedback can induce an inside-out stellar formation, with the successive generation of stars forming at larger radii. This scenario is plausible, since UFDs may have a bursty star formation history (Wheeler et al. 2019) – that is, high rates of star formation lasting a short period of time – as is expected for galaxies that formed most of their stars prior to reionization (Gelli et al. 2020); Gelli:2023; pallottini:2023; sun:2023. However, compared to more massive galaxies, the star formation rate (SFR) of UFDs is expected to be low – < 10⁻⁴ M_⊙ yr⁻¹ (Salvadori et al. 2014) – and thus the stellar feedback scenario requires energetic SNe to compensate for the low SN rate. Signs of SNe can be retrieved in the chemical composition of stars, and thus high-resolution spectroscopy is the primary way to determine if the UFD stellar halos have an in situ origin (e.g., Chiti et al. 2023). Merger events between UFDs are expected (Deason et al. 2014); Revaz:2023 and they can form extended stellar components beyond 5 − 20 R_1/2 (Deason et al. 2022); Ricotti:2022. Determining if a stellar halo formed via merging events requires a combination of spatial and chemical information to rule out the other scenarios. Currently, only in Tucana II does the stellar halo seem consistent at a high confidence level with the merger scenario (Chiti et al. 2021); chiti:2023.

Tucana II is a UFD with a stellar mass of $M_{\star }^{\mathrm{TucII}}=3_{-2}^{+7}\times {10}^3{M}_{\odot }$ and a half-light radius of $R_{1/2}^{\mathrm{TucII}}=120\pm 30$ pc (Bechtol et al. 2015), and it has one of the lowest average metallicities in the Local Group, ⟨[Fe/H]⟩ = −2.7 (Ji et al. 2016; Chiti et al. 2018). When observing the outskirts of Tucana II, Chiti et al. (2021) identified five and two giant stars located at distances of > 2 R_1/2 and > 5 R_1/2, respectively. Chiti et al. (2021, 2023) rule out the tidal interaction scenario for Tucana II because of: i) the orthogonality of the stellar extension with respect to the predicted direction of the tidal debris; and ii) the absence of a velocity gradient in the extended stellar component. The position of the pericenter (around ∼40 kpc; Battaglia & Nipoti 2022; Pace et al. 2022) is distant enough to corroborate a scenario of low impact from tidal interactions, even though recent studies of cosmological zoom-in simulations suggest that satellites experience tidal disruption at comparable pericenters (e.g. Shipp et al. 2023, 2024). The SN feedback scenario for the formation of the stellar halo in Tucana II is disfavored because of the low SFR and for two additional reasons. First, Chiti et al. (2023) report that 75% of the observed metal-poor stars ([Fe/H] < −2.9) are carbon-enhanced (i.e., the so-called CEMP-no, [C/Fe] > 0.7), which suggests enrichment by low-energy primordial “faint SNe” (e.g., Salvadori et al. 2015; Jeon et al. 2017; Rossi et al. 2023). Second, the chemical abundances of stars in Tucana II at various [Fe/H] are consistent with the ones observed in other UFDs. Therefore, if SN feedback is the main formation channel of stellar halos, most UFDs should present an extended feature, which instead is not observed. Ruling out the other possibilities, the remaining formation scenario for Tucana II stellar halo is the merging of two UFDs.

A result supporting such a scenario comes from the cosmological hydrodynamical simulations by Tarumi et al. (2021). They report that a galaxy closely resembling Tucana II in terms of the fraction of distant stars, metallicity, and metallicity gradient can be formed by a major merger; that is, a merger for which the mass ratio between the two galaxies is within 1:1 and 10:1. There have been other studies analyzing the impact of dwarf galaxy mergers on the formation of stellar halos. Using N-body simulations, Deason et al. (2022) find that intermediate merger ratios (∼5:1) maximize the growth of extended stellar halos regardless of the stellar mass–halo mass relation. By preparing cosmological N-body simulations, Ricotti et al. (2022) show that the shape and mass of the stellar halos can be described as the sum of nearly exponential profiles, with scale radii determined by the mass ratio of the mergers building it. Goater et al. (2024) analyzed tidally isolated UFDs in zoom-in hydrodynamical simulations, taken from the EDGE cosmological simulation suite (Agertz et al. 2020), finding that many galaxies exhibit anisotropic extended stellar halos mimicking tidal tails, but such halos are actually formed by late-time mergers. These studies primarily focus on the influence of the merger mass ratio on the stellar halo formation, with limited exploration of other potential merger or progenitor characteristics, leaving open questions such as whether the merger trajectory has any impact on the stellar halo formation, whether the dark matter total mass influences the size of the post-merger galaxy or the stellar halo formation, and what properties of the merger we can infer using present-day observations.

In this work, we explore the parameter space of UFD mergers, considering five key properties – the specific angular momentum, specific kinetic energy, merger mass ratio, total UFD mass, and initial radius – to model the trajectory and mass ratio of the merger and the mass and radius of the progenitors. Our goal is to investigate how these parameters influence the formation of the stellar halo in UFDs. We developed a suite of N-body simulations of isolated mergers between two possible Tucana II progenitors, to perform a wide (> 1 dex) parameter space exploration. Ideally, we would need: i) a large number of simulations, to perform a high-frequency sampling of the parameter space, and ii) single-stellar mass resolution, to have good precision in the identification of merger debris. To achieve both features, we developed an emulator that can efficiently generate the properties of the post-merger galaxies, once trained with simulation data. To validate our findings, we conducted a comparative analysis of simulated stellar surface density profiles and the one observed in Tucana II.

This article is structured as follows. First, we describe the simulation setup, the parameter space explored, and how we constructed the emulator (Sect. 2). Then, the results of the parameter space exploration and the comparison with observations are presented (Sect. 3). Finally, we summarize and discuss our findings in (Sect. 4).

2. Methods

Here, we present the suite of simulations for isolated merging UFD galaxies and detail our strategy for exploring the parameter space of the possible merger configurations.

We start by reviewing the base assumptions and the schemes used to generate the initial conditions (ICs) for the progenitors (Sect. 2.1). Then, we discuss the intervals adopted for the ICs (Sect. 2.2) and show how we can explore them with a Latin hypercube sampling of the merger parameter space (Sect. 2.3). Finally, we describe the numerical code used to follow the time evolution of the merger (Sect. 2.4) and we explain how we trained an emulator by using the results from our suite (Sect. 2.5).

2.1. Initial conditions

Our merging simulations were set up at redshift z = 6, after the reionization of the Universe, when most of the gas of UFDs – that is, minihalos – is expected to be removed or heated by radiation (Barkana & Loeb 1999; Sobacchi & Mesinger 2013; Gutcke et al. 2022). For this reason, we neglected the gas component and focused on dry mergers; that is, mergers of galaxies with only dark matter and stars.

Throughout the simulation suite, we varied the main properties of the dark matter and stellar component of the two merging galaxies together with the parameters defining the kinematic of the merger. Specifically, we explored: i) the dark matter mass, M_DM; ii) the stellar mass, M_⋆; iii) the scale radius of the progenitors; and iv) the trajectory of the merger. A detailed discussion of the ranges adopted for these quantities is given in Sect. 2.2 (see Table 1 for a summary), while, in the following, we describe the general setup and parametrization of the ICs.

2.1.1. Dark matter profile

We assume the dark matter density profile, ρ_DM, to be a spherically symmetric Navarro-Frank-White (Navarro et al. 1996):

$$\begin{array}{cc}\hfill {\rho }_{\mathrm{DM}}(r)& =\frac{3M_{\mathrm{DM}}}{4\pi R_{\mathrm{vir}}^3}{(ln(1+c)-\frac{c}{1+c})}^{-1}\times \hfill \\ & \times {\frac{r}{R_{\mathrm{vir}}}}^{-1}{(\frac{1}{c}+\frac{r}{R_{\mathrm{vir}}})}^{-2},\hfill \end{array}$$(1)

where r is the distance from the galaxy center, R_vir is the virial radius, and c is the concentration parameter. Considering the virial radius, defined as the radius where the density is Δ_c = 200 times the critical density of the Universe, we can express it as:

$$\begin{array}{c}\hfill R_{\mathrm{vir}}={\left[\frac{(M_{\mathrm{DM}}+M_{\star })G}{100H^2}\right]}^{1/3},\end{array}$$(2)

with G being the gravitational constant and H = H(z) the Hubble parameter¹; that is,

$$\begin{array}{c}\hfill H=H_0{[{\mathrm{\Omega }}_m{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{\Lambda }}]}^{1/2}.\end{array}$$(3)

The concentration is c = 3.5, assuming the results found by Diemer & Kravtsov (2015) for galaxies hosted in M_DM = 10⁸ M_⊙ halos at z = 6.

2.1.2. Stellar profile

We adopted an exponential profile, which is commonly utilized to fit stellar density profiles of UFDs (e.g., Belokurov et al. 2006; Martin et al. 2008; Muñoz et al. 2018). The profile is defined by

$$\begin{array}{c}\hfill {\rho }_{\star }(r)={\rho }_0{e}^{-r/(0.5R_{1/2})},\end{array}$$(4)

where ρ₀ is the central density, R_1/2 is the projected radius that contains half of the mass of the galaxy, in other words the half-mass radius, and the 0.5 factor arises integrating the exponential profile, factor 0.351, and projecting the three-dimensional scale radius to the projected two-dimensional R_1/2, factor 1.33 (Wolf et al. 2010). In the suite, we varied R_1/2 to account for different initial stellar densities, as is discussed in Sect. 2.2.

2.1.3. Discretization of the progenitors

We used the public code DICE (Perret 2016) to perform a Lagrangian sampling of the particles composing each progenitor in order to obtain the selected density profiles for dark matter (Eq. (1)) and stars (Eq. (4)). DICE uses an N-try Markov chain Monte Carlo method to extract the particles’ position, x_s, from the profiles and integrates the Jeans equations to assign the velocity, v_s, for each particle, s; this is done separately for the progenitors.

To match the observed stellar mass of Tucana II, we fixed the sum of the stellar mass of the two progenitors to M_⋆,tot = 6000 M_⊙. To compare our results with star-by-star observations, we opted for a single-stellar mass resolution; in all the simulations, the mass resolution of the stellar component is m_⋆ = 1 M_⊙. The number of dark matter particles was adjusted to have a mass resolution of m_DM = 100 M_⊙. The dark matter to stellar mass resolution ratio adopted is consistent with the one used in Ricotti et al. (2022) and satisfies the trade-off between computational cost and spurious numerical heating (e.g., Ludlow et al. 2023; Di Cintio 2024). The numerical heating mostly affects the size of the post-merger galaxy, while it is negligible for the stellar halo formation, as is shown in Appendix A.

As an example of the procedure, in Fig. 1 we report the stellar density profile (left panel) and the stellar surface density map on the XY plane (right panel) generated by DICE for a progenitor with M_⋆ = 5.48 × 10³ M_⊙ and R_1/2 = 118 pc. In the left panel, the solid red line shows the exponential profile fit to the progenitor density profile. Since DICE computes galaxies at Jeans’ equilibrium, the initial progenitors are not at virial equilibrium and must evolve in isolation to virialize. We address this aspect in the next paragraph (Sect. 2.1.4).

Fig. 1. Example of the virialized stellar distribution of a progenitor with M⋆ = 5.48 × 103 M⊙ and R1/2 = 118 pc generated using DICE (see Sect. 2.1.3). Left panel: Stellar density (ρ⋆) extracted by spherically averaging the mass in a simulation shell (points) and the fit exponential profile (solid line, Eq. (4)) as a function of the radius (r). Right panel: Stellar surface density map, Σ⋆ integrated along the line of sight; the spatial scale is indicated as an inset. In both panels, the dashed lines mark the location of spheres with radius R1/2, 3 R1/2, and 5 R1/2.

2.1.4. Merger setup

Using DICE, we placed the two progenitors in the simulation box, assigned them a relative velocity, set the system’s center of mass at rest in the center, and placed the progenitors at a distance of D = 8 kpc along the X axis, which ensured that they virialized while the gravitational interaction was still negligible. Hereafter, the subscript 1 refers to the more massive progenitor and 2 to the other one.

We set the progenitors on the same Z axis, assuming that the velocity differs from zero only along the X axis. The initial position of the progenitors can be written as

$$\begin{array}{cc}\hfill X_{1,2}& =\mp D\frac{M_{2,1}}{M_1+M_2},\hfill \end{array}$$(5a)

$$\begin{array}{cc}\hfill Y_{1,2}& =\mp \frac{l}{|V_1|+|V_2|}\frac{M_1+M_2}{M_{1,2}},\hfill \end{array}$$(5b)

with their velocity along the X axis being

$$\begin{array}{c}\hfill V_{1,2}=\pm \sqrt{2k\frac{M_{2,1}}{M_{1,2}}}.\end{array}$$(5c)

In Eq. (5), M_i is the total mass of the i-th progenitor (M_i = M_i,DM + M_i,⋆), l is the specific angular momentum (with regard to the center of mass),

$$\begin{array}{c}\hfill l=\frac{{\displaystyle {\sum }_s}{m}_s{\mathit{v}}_s\times {\mathit{x}}_s}{{\displaystyle {\sum }_s}{m}_s},\end{array}$$(6a)

where the sum over s is extended to all particles, and k is the specific kinetic energy,

$$\begin{array}{c}\hfill k=\frac{{\displaystyle {\sum }_s}{m}_s{\mathit{v}}_s^2}{{\displaystyle {\sum }_s}{m}_s}.\end{array}$$(6b)

We varied the values of l and k to probe different merging trajectories. Figure 2 shows an example of the merger’s ICs, in which the two progenitors have total masses of M₁ = 4.8 × 10⁷ M_⊙ and M₂ = 4.6 × 10⁶ M_⊙, the specific merger kinetic energy is k = 0.72 pc² Myr⁻², and the angular momentum is l = 1.1 × 10³ pc² Myr⁻¹.

Fig. 2. Example of merger ICs prepared with DICE (see Sect. 2.1.4). The density map in the X − Y plane shows the total (dark matter+stars) density (ρ) that has been mass-averaged along the Z axis. The two progenitors were placed in the simulation box at a distance of D = 8 kpc along the X axis and with an impact parameter of |Y1|+|Y2| (Eq. (5b)). The arrows mark the direction of the initial velocity of each progenitor (Eq. (5c)).

2.2. Parameterization of the merger

We focus our exploration on the five parameters: the merger mass ratio, M₁/M₂, the dark matter content, M_DM/M_⋆, the stellar radius of the progenitors, R_1/2, the specific angular momentum, l, and the specific kinetic energy, k. In the following, we first describe the interval explored for each of these properties and then present how we explore the parameter space via a Latin hypercube sampling. Table 1 summarizes the properties and relative intervals explored.

2.2.1. Merger mass ratio

Defining the merger mass ratio between the two progenitors as M₁/M₂, we note that previous works (e.g., Deason et al. 2022) identify a peak in the fraction of stars in the stellar halo when M₁/M₂ ≃ 5. Therefore, to fully capture its dependency in UFD mergers, we explored the interval M₁/M₂ ∈ [1, 30].

It is convenient to use the merger mass ratio to define the stellar mass of each progenitor:

$$\begin{array}{cc}\hfill M_{1,\star }& =M_{\mathrm{tot},\star }\frac{M_1/M_2}{1+M_1/M_2}\hfill \end{array}$$(7a)

$$\begin{array}{cc}\hfill M_{2,\star }& =M_{\mathrm{tot},\star }\frac{1}{1+M_1/M_2},\hfill \end{array}$$(7b)

where M_tot, ⋆ = M_1, ⋆ + M_2, ⋆ is the total stellar mass in the simulation, which was set to M_tot, ⋆ = 6 × 10³ M_⊙ to ensure that, even in the case of stellar mass loss during the merger, the post-merging galaxy stellar mass is consistent with Tucana II estimations, $M_{\star }^{\mathrm{TucII}}=3_{-2}^{+7}\times {10}^3{M}_{\odot }$ (Bechtol et al. 2015).

2.2.2. Dark matter to stellar mass ratio

The dark matter content of UFDs is poorly constrained. Velocity dispersion observations can determine the dark matter mass within the luminous part of the galaxy; however, extrapolating the total halo masses from such a measurement is an uncertain endeavor (e.g., Errani et al. 2018). For this reason, we explored an interval of values of M_DM, investigating whether the dark matter abundance impacts the stellar halo formation. To this end, we varied the stellar to dark matter mass ratio, M_DM/M_⋆, which determines the dark matter mass of each progenitor once the merger mass ratio sets the stellar mass.

The interval explored is M_DM/M_⋆ ∈ [8, 100]×10³, and, in Fig. 3, we report it in an M_⋆ - M_DM plane with a shaded region. The selected range includes the latest dark matter mass estimate within 1 kpc for Tucana II, $M_{\mathrm{DM}}^{\mathrm{TucII}}\sim 2.5\times {10}^7{M}_{\odot }$ (Chiti et al. 2021). Further, the explored range roughly brackets the predictions by abundance-matching methods performed by Behroozi et al. (2013) for galaxies between z = 6 and z = 0 (green lines in Fig. 3). Recently, a substantial effort has been carried out to constraint the stellar mass for smaller system both from local galaxy observations (e.g., Jethwa et al. 2018; Nadler et al. 2020) and simulations (e.g., Munshi et al. 2021; Kim et al. 2024). The explored range in this work encompasses most of the scatter measured at z = 0 in simulated dwarf galaxies by Kim et al. (2024, green-shaded area in Fig. 3).

Fig. 3. Parameter sampling in the dark matter (MDM) vs. stellar mass (M⋆) plane. Each pink (purple) point represents the least (more) massive progenitor in a merging simulation. The shaded blue region indicates the dark matter to stellar mass (MDM/M⋆) interval explored by the suite. The gold star shows the estimate of Tucana II from Chiti et al. (2021) observations, while the green lines are the extrapolation from the abundance matching model from Behroozi et al. (2013) at redshift z = 0 and 6. The shaded green region represents the extrapolation of the scatter in the stellar-mass to halo-mass relation measured by Kim et al. (2024) for dwarf galaxies.

2.2.3. Half-mass radius

Observations show that UFDs with a stellar mass similar to Tucana II have a half-light radius – the radius that contains half of the galaxy’s luminosity – ranging from tens to hundreds of parsecs (Simon 2019). The half-light radius is equal to the half-mass radius, assuming a constant luminosity-to-mass ratio, and thus we decided to explore an interval of values, R_1/2 ∈ [15, 150] pc.

We assigned the value of R_1/2 to the more massive progenitor, and we defined the mean stellar density as

$$\begin{array}{c}\hfill {\overline{\rho }}_{\star }=\frac{3}{8\pi }\frac{M_{1,\star }}{R_{1/2}^3}.\end{array}$$(8)

To set the scale radius of the other progenitor, we imposed that the two galaxies have the same ${\overline{\rho }}_{\star }$; that is, the scale radius of the less massive was determined by R_1/2(M₂/M₁)^1/3. This approach prevented us from exploring the dependency of the stellar halo formation on the density difference between the two progenitors. However, by avoiding double sampling the R_1/2 parameter, we reduced the number of dimensions explored, thereby improving the robustness of our results.

In Fig. 4, we report the sampling of the M_⋆-R_1/2. We considered the stellar mass and half-light radii of all the UFDs in the Local Universe Database (Pace 2024) and performed a fit with a polynomial function. The fit ±1 − σ range is reported as the green-shaded area in Fig. 4. The overlap between the sampled values and the observations in the M_⋆-R_1/2 plane confirms that the adopted scaling and range of the radii are adequate to explore the physical parameter space.

Fig. 4. Parameter sampling in the stellar mass (M⋆) vs. scale radius (R1/2) plane. Each pink (purple) point represents the least (more) massive progenitor in a merging simulation and the shaded blue area indicates the scale radius interval explored in the suite. The shaded green region represents the polynomial fit ±1σ of all the available UFD data in the Local Volume Database (Pace 2024).

2.2.4. Specific angular momentum and kinetic energy

Since we want to study the merging event between galaxies, we limited our explorations to values of l and k that are highly likely to produce the merging of the two progenitors. Specifically, the intervals explored are l ∈ [100, 3000] pc² Myr⁻¹ and k ∈ [0.1, 3] pc² Myr⁻². We adopted these intervals after testing them with (i) the solution of the two-body problem, assuming both galaxies as point masses, and (ii) a low-resolution version of the simulation suite, as is further explained in the following.

2.3. Parameter space sampling

Sampling the parameter space is a delicate endeavor. The more samples we have, the more finely we explore the parameter space, with increasing computational cost. We used a five-dimensional Latin hypercube sampling (McKay et al. 1979) to maximize the exploration, while keeping the number of samples (simulations) low. The implemented sampling consists of 47 nodes (sampled points) with each node being a unique set of the five parameters defining the merger; that is,

$$\begin{array}{c}\hfill \mathit{\kappa }=(l,k,M_1/M_2,M_{\mathrm{DM}}/M_{\star },R_{1/2},).\end{array}$$(9)

To compute the node values with the LHS prescription, we divided the range of each individual parameter into a set of equally probable intervals and the algorithm selected one sample from each interval, ensuring that the samples were evenly spaced out across the entire range. We adopted a logarithmic scaling for the five components in the node κ.

We verified that the progenitors were on merging trajectories before evolving the ICs generated by the 47 nodes at full resolution. To do so, we generated the ICs at lower dark matter resolution, $m_{\mathrm{DM}}^{\mathit{low}}={10}^3{M}_{\odot }$, and evolved them. The coarser resolution makes these simulations run at around a tenth of the computational cost, enabling us to evolve each one for 5 Gyr and check whether or not the progenitors merge. In addition, having the same ICs simulated using two different resolutions allows us to evaluate the resolution effects, as is detailed in Appendix A.

2.4. Time evolution

We used the public version of the moving-mesh hydrodynamic code AREPO (Springel 2010; Weinberger et al. 2020) to evolve the ICs. AREPO performs the time evolution with a second-order accurate leapfrog scheme, using a tree-particle-mesh method (Bagla 2002) to compute the gravitational interaction between particles. Having just two galaxies in the simulation box, we adopted only the oct-tree grouping scheme (Barnes & Hut 1986) for the force evaluation, using nonperiodic boundary conditions. The gravitational softening length of dark matter and stellar particles is ϵ = 1 pc, and the simulation box has side L = 100 kpc.

Each simulation evolved until the two progenitors merged into a single, virialized galaxy. The merging instant was identified using the on-the-fly friends-of-friends group finder (FoF) available within AREPO and originally described in Springel et al. (2001). The algorithm groups together particles within a distance of 0.2 times the mean interparticle separation. We defined the merging instant as the point at which the FoF algorithm finds a single galaxy. After the merging instant, each simulation evolves for an additional ∼1 Gyr, to ensure that the two progenitors fully merge and that the virial ratio is 2T/U ≃ 1 for several hundred million years, with T being the kinetic energy and U the potential energy.

2.5. Merger emulator

Having a high computational cost is a common limitation in numerical simulations, especially when a parameter space exploration is involved; a possible solution is the implementation of an emulator, to either reduce the number of simulations needed (e.g., Bird et al. 2019; Rogers et al. 2019; Brown et al. 2024) or to speed up the computation (e.g., Spurio Mancini et al. 2022; Branca & Pallottini 2024; Robinson et al. 2024; Bartlett et al. 2024). In this work, we have adopted a neural network (NN) as an emulator to enhance the parameter space exploration. The NN approximates a function that, taking the merger configuration and progenitor ICs as input, returns the properties of the post-merger galaxy without requiring new simulations. Indeed, thanks to the universal approximation theorem (UAT, Cybenko 1989), a trained NN can approximate with arbitrary precision the nonlinear function that maps the input parameters, κ, into outputs characterizing the galaxy and the stellar halo. For the present work, we set two outputs,

$$\begin{array}{c}\hfill (R_{\star },f_5)=\mathrm{NN}(\mathit{\kappa }),\end{array}$$(10)

where with R_⋆ we refer to the half-mass radius of the post-merger galaxy to distinguish it from the massive progenitor’s half-mass radius, R_1/2, which is one of the inputs of the NN. f₅ is the stellar halo mass fraction; that is, the fraction of stars at r > 5R_⋆. Such a choice of outputs enables us to directly compare our results with Tucana II observations (see Sects. 3.2 and 3.3).

We adopted a regularized loss function to prevent overfitting and trained the emulator on the simulation outputs. The number of training points is small compared to the dimensionality of the problem and we overcame the sparse sampling issue by adopting a Latin hypercube sampling (see Sect. 2.3), which ensures uniform coverage of the parameter space. Furthermore, deep NNs are relatively unaffected by the growing difficulty of analyzing data as the number of dimensions increases; that is, the “curse of dimensionality”, when approximating compositional functions (Poggio & Liao 2018). We refer to Appendix B for a complete description of the NN.

3. Results

First, we present the evolution of a single simulation and its analysis to overview the merging process (Sect. 3.1). Then, we repeat the same analysis on all the simulations, use the results to train the emulator, and explore the parameter space to determine the dependency of the stellar halo on the ICs (Sect. 3.2). Finally, we directly compare our findings with the observations of Tucana II (Sect. 3.3).

3.1. Merging of ultra-faint dwarfs: An overview

To showcase the merging process, we selected simulation SIM_47, which produces a post-merger galaxy with a stellar surface density consistent with the one observed in Tucana II (see Sect. 3.3). Simulation SIM_47 has a merger mass ratio of M₁/M₂ ≃ 10, scale radius of R_1/2 ≃ 120 pc, mass-to-light ratio of M_DM/M_⋆ ≃ 10⁴, specific angular momentum of l ≃ 1.1 × 10³ pc² Myr⁻¹, and specific kinetic energy of k ≃ 0.72 pc² Myr⁻². We note that this is the set of parameters already used in Figs. 1 and 2, which show the ICs and merger setup, respectively.

To visualize the spatial evolution in the simulated volume, we report the stellar surface density for the various phases of the merger in Fig. 5. The top left panel is at t = 0.2 Gyr when the two progenitors are virialized. The top right panel shows the gravitational interaction between the two progenitors during the third and last passage before merger completion. The bottom panels are at the end of the simulation; that is, at t = 4.01 Gyr, when the post-merger galaxy is virialized, and when the virial ratio is around 1. The effect of the merger on the morphology of the galaxy is evident when comparing the post-merger galaxy with one of its progenitors (the right panel in Fig. 1). To help the comparison, we mark 1, 3, and 5 R_⋆ in both panels of Fig. 1 and Fig. 5. The progenitor is mostly contained within 5 R_1/2, while the post-merger galaxy is wider, extending way beyond 5 R_⋆. The extended stellar distribution at large radii – the stellar halo – forms as a result of the merger event and is bar-shaped. Therefore, contrary to the progenitors, the post-merger galaxy is asymmetric and we account for the asymmetry using an ensemble average² to estimate R_⋆, which we recall is the half-mass radius of the 2D projection of the galaxy. We address the differences between the 3D and 2D half-mass radius in Appendix C.

Fig. 5. Example of merger evolution. The different panels report the stellar surface density (Σ⋆) maps at various evolutionary stages for progenitors with a M1/M2 ≃ 10 (profile in Fig. 1) and the merger setup in Fig. 2. The bottom right panel is a zoom-in on the virialized post-merger galaxy, and the three dashed circles have radii of one, three, and five half-mass radii. After the merging, the resulting fraction of stars in the outskirts is f5 ∼ 5.5%.

To analyze the difference between pre- and post-merger galaxies’ stellar distribution, in Fig. 6 we compare the stellar density of the two progenitors (shades of purple) and the post-merger galaxy (orange). The vertical lines are the half-mass radius for each galaxy, while the solid lines are the fit profiles: exponential (Eq. (4)) for both progenitors and the sum of two exponential profiles for the post-merger galaxy. The double exponential profile models the formation of a dense core and a sparse stellar halo and better fits the post-merger stellar distribution. The difference between the post-merger galaxy and its massive progenitor stellar distribution lies in the outskirts, where the post-merger galaxy has a higher mass fraction due to the formation of a stellar halo. The higher mass fraction affects the half-mass radius, going from R_1/2 = (114 ± 1) pc of the massive progenitor to $R_{\star }={130}_{-1}^{+3}$ pc of the post-merger galaxy. The error distribution in R_⋆ is due to the asymmetry of the stellar halo formed via a merger.

Fig. 6. Example of pre- vs. post-merging stellar distributions. We show the profile for the simulation selected for Fig. 5 by reporting ρ⋆ for the two progenitors with shades of purple and the post-merger galaxy in orange. The solid lines are the fit profile, which is exponential (Eq. (4)) and a sum of exponential for the progenitors and the post-merger galaxy, respectively. The vertical lines mark the half-mass radius for each galaxy.

Tidal interactions during the merger strip stars from the less massive progenitor and form the stellar halo, as is shown in Fig. 7, where we report the fraction of stars of the post-merger galaxy coming from each progenitor at different radii. We adopted the same progenitor colors of Fig. 6, and the dashed lines are at 1, 3, and 5 R_⋆. While the center of the post-merger galaxy is dominated by stars of the more massive progenitor, beyond 3 R_⋆ the majority of the stars come from the smaller progenitor. This distribution shows that stars of the low-mass progenitor compose the stellar halo, and this is key to understanding the trends we observe exploring the parameter space (Sect. 3.2).

Fig. 7. Origin of the stars in the post-merger galaxy. For the simulation selected for Fig. 6, we show the radial distribution of the fraction of stars that originally belonged to the high-mass progenitor (purple) and the low-mass progenitor (pink). The dashed vertical lines are at 1, 3, and 5 R⋆.

To quantify the presence of a stellar halo, we followed Tarumi et al. (2021) and calculated f₅, which we recall is the fraction of stars beyond 5 R_1/2. Since f₅ depends on the half-mass radius, we performed the same ensemble average considering all possible projection planes. For each projection plane, we randomly sampled the stellar surface density distribution using the same number of stars of Chiti et al. (2021, i.e. 19) and considered the mean and standard deviation as f₅ and its error. We refer the reader to Appendix C for a discussion on the differences between 3D and projected properties. We estimate f₅ = 0.16 ± 0.4% and f₅ = 0.19 ± 0.03% for the two progenitors. These values are consistent with the ones expected in isolated galaxies with an exponential profile, where f₅ tends to zero. Figure 5 shows the formation of a stellar halo after the merger, and as a result of the stellar halo formation, the fraction of outskirts stars in the post-merger galaxy increases to f₅ = 5.4 ± 1.8%. For comparison, the observed value in Tucana II is f₅^TucII = 2/19 = 10.5%.

3.2. Exploration of the parameter space with the emulator

Using the radius, R_⋆, and the fraction of stars in the outskirts, f₅, of each post-merger galaxy, we built the dataset to train the emulator. We recall that (Sect. 2.5) the input values of the emulator are the explored initial parameters, l, k, M₁/M₂, M_DM/M_⋆, and R_1/2, and the output values are the median of the two post-merger quantities, f₅ and R_⋆ (see Eq. (10); for details of the NN architecture and its training, we refer to Appendix B).

Figure 8 shows the distribution of the emulated outputs (f₅ top row, R_⋆ bottom row) of 10⁶ different ICs log-spaced distributed. The orange dots represent the results of the training simulations, and their overlap with the emulated results confirms the emulator’s robustness. The last column of Fig. 8 shows the distribution of output values in the parameter space explored. The distribution of f₅ peaks around 5%, the value found by Tarumi et al. (2021) for their dry merger. The R_⋆ distribution increases at low values. This trend is primarily due to the distribution of the training data, which is affected by numerical heating at low R_1/2, as is explained in detail in Appendix A. Consequently, the emulator predicts a higher frequency of galaxies with low R_⋆ due to the higher number of simulations with R_⋆ ≃ 40 pc.

Fig. 8. Emulation of the post-merging f5 (upper panels) and stellar radius (R⋆, lower panels) as a function of merger setup parameters, i.e., specific angular momentum (l, Eq. (6a)), specific kinetic energy (k, Eq. (6b)), mass ratio (M1/M2), dark to stellar mass ratio (MDM/M⋆), and pre-merger stellar size (R1/2, Eq. (4)). Each of the first five columns gives the distribution of outputs of the merging process as a function of a single parameter. The distribution is obtained by emulating 106 IC log-spaced distributed (see Table 1), and the gold line represents its median. The orange points are the individual simulations used to train the emulator. The red line (shaded region) gives the average (dispersion) of the observed values for Tucana II (Bechtol et al. 2015; Chiti et al. 2021). The last column shows the distribution of the emulated output normalized to the total number of inferences.

The red-shaded areas are the observed values for Tucana II with the relative error. There are several estimates of the angular size of Tucana II (Bechtol et al. 2015; Koposov et al. 2015; Moskowitz & Walker 2020) and distances (Bechtol et al. 2015; Koposov et al. 2015; Vivas et al. 2020) each resulting in different half-light radii. To ensure consistency with Chiti et al. (2021), we adopted the half-light radius $R_{1/2}^{\mathrm{TucII}}=(120\pm 30)$ pc (Bechtol et al. 2015). This choice allows for a direct comparison of our simulated stellar distribution with their observational data. We computed the fraction of distant stars, f₅^TucII = 10 ± 5%, from the distribution of individual stars observed in Chiti et al. (2021). The error assigned to f₅ considers that one of the two stars in the outskirts can lie within 5 R_1/2, changing the fraction to 5%.

The emulated results, whose median is given by the gold line, show little to no dependency on the specific angular momentum, l, kinetic energy, k, and dark matter content, M_DM/M_⋆. Clear trends are instead present for the merger mass ratio, M₁/M₂, and the progenitor half-mass radius, R_1/2; specifically, R_⋆ correlates almost linearly with R_1/2, while f₅ increases with M₁/M₂, peaks around M₁/M₂ ≃ 6, and then decreases. The correlation between the sizes of the progenitor and the post-merger galaxy is expected, since the larger the progenitor is, the larger the post-merger galaxy will be. Meanwhile, the correlation between f₅ and the merger mass ratio emerges because, as is discussed in Sect. 3.1, the stellar halo stars primarily originate from the less massive progenitor (the “satellite”), leading to two distinct regimes. When the two progenitors have comparable masses, they mix, forming a post-merger galaxy with a minimal stellar halo and a small f₅. In this regime, increasing the merger mass ratio weakens the depth of the potential well of the satellite, making it easier for its stars to be stripped during the merger and form a more prominent stellar halo. As a result, f₅ initially increases with the M₁/M₂ ratio. However, increasing M₁/M₂ reduces the number of stars in the “satellite”, leading to a decrease in the number of stars available to form the stellar halo. For this reason, beyond a certain threshold (approximately M₁/M₂ ≃ 6), f₅ begins to decline.

To quantify the dependencies of f₅ and R_⋆ from the ICs (κ), we performed a principal component analysis (PCA). A PCA studies the linear transformation (component) of the parameter space of the data. Such components are ranked so that the variance decreases. Thus, the first (principal) components can explain most of the variation in the original dataset. We performed the PCA separately for the two outputs; that is, one for (f₅, κ) and one for (R_⋆, κ). For each PCA, we considered 10⁶ points in six dimensions, normalizing each dimension such that its mean and standard deviation are 0 and 1, respectively. The first PCA component for f₅ (R_⋆) is dominated by f₅ (R_⋆) and M₁/M₂ (R_1/2), confirming the trends observed in Fig. 8.

To determine the secondary dependencies, we performed two additional PCAs considering a reduced number of points. Specifically, when analyzing f₅, we considered the points with M₁/M₂ ∈ [5, 7], where f₅ peaks, while for R_⋆, we considered the points with the progenitor radius R_1/2 ∈ [90, 110] pc; that is, the progenitor radius producing a post-merger galaxy radius consistent with Tucana II observations. The first reduced PCA component for f₅ (R_⋆) is dominated by f₅ (R_⋆) and R_1/2 (M₁/M₂), revealing that the principal dependency of one output is the secondary dependency of the other.

Once we removed the dependency on M₁/M₂ and R_1/2, selecting points simultaneously in the intervals M₁/M₂ ∈ [5, 7] and R_1/2 ∈ [125, 145] pc, we noticed that f₅ and R_⋆ marginally depend on the other three parameters, l, k, and M_DM/M_⋆. The relative importance depends on the intervals considered, but generally M_DM/M_⋆ is the least important for both f₅ and R_⋆.

3.3. Comparison with Tucana II observations

We can now compare the emulated f₅ and R_⋆ with the observations to infer the properties of Tucana II progenitors. As is described in Sect. 3.2, we adopted $R_{1/2}^{\mathrm{TucII}}=(120\pm 30)$ pc and f₅^TucII = 10 ± 5% as observed values of half-mass radius and fraction of stellar halo stars in Tucana II, respectively.

Figure 9 shows the distribution of the output values consistent with Tucana II observations, and based on the PCA results, we limited the study to two initial parameters that define the output values; that is, M₁/M₂ and R_1/2. The top row shows the distribution of these initial parameters, with the dashed lines indicating the median and 16-th and 84-th percentiles, while the four orange points represent four simulations in the suite consistent with Tucana II observations. Consistency with observations is achieved for ICs with a progenitor half-mass radius of $R_{1/2}={97}_{-18}^{+25}$ pc and a merger mass ratio of $M_1/M_2=8_{-3}^{+4}$. Therefore, we determine that an intermediate merger (5 < M₁/M₂ < 12) forms a stellar halo similar to the one observed for Tucana II. Such a finding is consistent with Deason et al. (2014), who report that 40% of z = 0 UFDs in their simulations underwent at least one major merger in the last 12 Gyr. Among the four post-merger simulations consistent with observations, SIM_41 and SIM_47 are consistent within 1.1σ with³ the inferred values of M₁/M₂ and R_1/2, SIM_17 has a slightly larger progenitor radius, > 1.1σ, and SIM_39 has both M₁/M₂ and R_1/2 values > 1.1σ. We note that the errors on the inferred properties highly depend on the value and error adopted for f₅, and better constraining it helps to decrease the uncertainties on M₁/M₂ and R_1/2. For example, lowering f₅ error by a factor of two – ∼ ± 2.5% – results in a decrease of a factor of three on the error of the ICs consistent with Tucana II.

Fig. 9. Distributions of the input couple R1/2 and M1/M2 that give f5 and R⋆ simultaneously consistent with the observations from Tucana II. Similarly to Fig. 8, the red line gives the average value inferred from observations, but the ranges for f5 and R⋆ were cut to include only the observational uncertainty, i.e., f5TucII = 10 ± 5% (Chiti et al. 2021) and $R_{1/2}^{\mathrm{TucII}}=(120\pm 30)$ pc (Bechtol et al. 2015), respectively. The orange points are the post-merger galaxies consistent with observations from Tucana II, and are numbered with their simulation ID.

One method to reduce the error on f₅ is to increase the number of identified stars; for example, using Gaia data. In recent work, Jensen et al. (2024) exploit Gaia data, employing a selection algorithm based on spatial, color-magnitude, and proper motion information targeting low brightness features of MW dwarf satellites. They are able to identify the second low-density outer profile of nine satellites, including Tucana II, and for each of them, to derive the stellar surface density profile. They model it as the sum of two exponential profiles:

$$\begin{array}{c}\hfill \mathrm{\Sigma }(r)={\mathrm{\Sigma }}_0(e^{-r/r_e}+B{e}^{-r/r_s}),\end{array}$$(11)

where Σ₀ is the central surface density, r_e is the exponential scale radius, B is the normalization of the outer component (0 ≤ B ≤ 1), and r_s the scale radius of the outer component (1.68r_s > R_1/2). They derive r_e from the literature value of the half-light radius such that R_1/2 = 1.68r_e. Unfortunately, we cannot use their star candidates list to evaluate f₅, since there are no spectroscopic follow-ups to verify the membership of each star. However, we can test our results by comparing the surface density of our simulated post-merger galaxies with the one derived by them.

Figure 10 shows the derived stellar surface density profile of Tucana II (green) and our simulated post-merger galaxies (orange), and each of the four columns refers to a simulation that is simultaneously consistent with the observed values of f₅ and R_1/2. As we can see, most of these four simulations reasonably match the observed stellar halo well beyond 1 kpc and only SIM_39 is not consistent with the observed data. We quantified the accordance between simulations and observations using the standard χ² procedure:

Fig. 10. Comparison of the observed stellar number density distribution of Tucana II and the simulations with the f5 and R1/2 parameters closest to the observed values (f5TucII = 10 ± 5% from Chiti et al. 2021 and $R_{1/2}^{\mathrm{TucII}}/\mathit{pc}=120\pm 30$ from Bechtol et al. 2015). Each panel shows a different simulation, as is indicated in the inset, ordered by decreasing χ2 value. The green line (shaded area) gives the mode (error) of the profile for Tucana II, observed by Jensen et al. (2024). The orange line (points with bars) gives the best fit (average and standard deviation of the data points) of the profile of a simulated post-merger galaxy, obtained by adopting a two-exponential profile (Eq. (11)). The solid (dashed) gray line shows the R⋆ (5 R⋆) post-merger radius.

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _{i=1}^{10}}\frac{{({\mathrm{\Sigma }}_{i,\mathrm{obs}}-{\mathrm{\Sigma }}_{i,\mathrm{sim}})}^2}{{\sigma }_{i,\mathrm{obs}}^2+{\sigma }_{i,\mathrm{sim}}^2},\end{array}$$(12)

where Σ_i, obs (σ_i, obs) and Σ_i, sim (σ_i, sim) are the stellar surface density profiles (standard errors) derived from observations and simulations, respectively. We used ten radius bins to be consistent with Jensen et al. (2024), and for each bin we computed σ_i, obs propagating R_1/2, B, and r_s errors, and we computed σ_i, sim as the Poisson error. The three simulations within 1.1σ from the inferred M₁/M₂ – SIM_17, SIM_41, and SIM_47 – have a good agreement; that is, χ² = 0.5 − 2. It should be noted that χ² of SIM_17 is larger, since the accordance with R_1/2 is worse. All three simulations show the two-component profiles with differences in the central density Σ₀ and the transition to the second exponential profile. Their ICs span almost all the angular momentum values, confirming that the formation of a stellar halo does not depend on angular momentum in the interval explored in this work. We recall that we considered l values that result in a merger between the two progenitors. For SIM_39, the surface density significantly diverges from the observed one (χ² = 12.0). Although the fraction of distant stars is consistent with observations, $f_5=6.4_{-4}^{+2}$, the simulation features an inner core that is 2× denser and 0.5× smaller compared to observations. This is due to the offset from M₁/M₂ and R_1/2, which are at > 1.1 σ from the emulator predictions.

4. Summary and discussion

In this work, we investigate merger events as the formation scenario of stellar halos around UFDs. We study the impact of different merger and progenitor properties and focus on reproducing the stellar halo observed in Tucana II by Chiti et al. (2021, 2023), Jensen et al. (2024).

To this end, we developed a suite of 47 N-body simulations of isolated dry mergers between two idealized progenitors of Tucana II with a stellar resolution of a single solar mass. We used DICE (Perret 2016) to generate the ICs, assuming an exponential profile for the stellar component of both progenitors, and we used AREPO (Weinberger et al. 2020) to evolve the simulation until the post-merger galaxy virializes. The suite explores different merger trajectories and progenitors properties, assigning to each simulation a unique set of five parameters: i) the specific angular momentum, l, ii) the specific kinetic energy, k, iii) the merger mass ratio, M₁/M₂, iv) the dark matter content, M_DM/M_⋆, and v) the half-mass radius of the more massive progenitor, R_1/2. For each post-merger galaxy, we quantified the formation of a stellar halo using two properties: the half-mass radius, R_⋆, and the fraction of stars, f₅, at high radii (> 5R_⋆).

We developed a NN-based emulator to avoid the extreme computational cost expected from a high-resolution sampling of the selected five-dimensional parameter space. Via the emulator, we quickly obtained the properties of one million post-merger galaxies – R_⋆ and f₅ – and studied the dependency of the stellar halo formation on the ICs of the merger and progenitors. By comparing the emulated properties to Tucana II observations, we inferred consistent merger and progenitor characteristics. Additionally, we compared the stellar surface density of simulated post-merger galaxies and Tucana II to further validate the merger scenario.

We can summarize our most important results as follows:

• The merger mass ratio, M₁/M₂ (progenitors’ size, R_1/2), drives the formation of a stellar halo, f₅ (the size of the galaxy, R_⋆), in the aftermath of a merger, with the progenitors’ size, R_1/2 (merger mass ratio, M₁/M₂), being the secondary dependency. In the interval explored, the specific angular momentum, l, specific kinetic energy, k, and dark matter abundance, M_DM/M_⋆, have a negligible impact on both the stellar halo formation and the size of the post-merger galaxy.

• Using the observed fraction of distant stars, f₅^TucII = 10.5% (Chiti et al. 2021), and half-light radius, R_⋆^TucII = 120 pc (Bechtol et al. 2015), of Tucana II, we predict that the merger mass ratio and massive progenitor size consistent with observations are $M_1/M_2=8_{-3}^{+4}$ and $R_{1/2}={97}_{-18}^{+25}$ pc. The merger that can form Tucana II’s stellar halo is intermediate, while the size of present-day Tucana II is similar to the radius of its massive progenitor.

• Four simulated post-merger galaxies are consistent with f₅^TucII and R_⋆^TucII. Three simulations have M₁/M₂ (R_1/2) within 1.1σ (1.2σ) from the predicted values: their stellar surface density reproduce (χ² ≈ 0.5 − 2) the observed one (Jensen et al. 2024) well. The simulation with both M₁/M₂ and R_1/2 at > 1.2σ from the predictions deviates the most from observations (χ² = 12). This underlines the importance of M₁/M₂ and R_1/2 in determining the post-merger galaxy morphology.

Our results confirm previous findings about the importance of the merger mass ratio on the formation of a stellar halo in UFDs (Deason et al. 2022; Ricotti et al. 2022), show the impact of the progenitor size on the radius of the post-merger galaxy, and highlight the possibility of inferring progenitor and merger properties using present-day observations. The agreement between the results from our simulations and emulations with observations in terms of f₅, R_⋆, and stellar surface density further corroborates the merging scenario for the formation of Tucana II’s stellar halo, suggested by observations (Chiti et al. 2021, 2023) and simulations of two specific UFD progenitors (Tarumi et al. 2021).

The NN-based emulator proves to be a fundamental tool to perform the parameter space exploration, especially when computational resources are limited, and it can be easily adapted to infer any scalar property from simulations; for example, the post-merger galaxy’s total mass. To emulate functional properties, such as the parametrical dependence of density profiles or surface brightness, and directly compare them with observations (e.g. Jensen et al. 2024; Conroy et al. 2024), it is required to change the architecture to a neural operator-based model (Lu et al. 2019). Such a possibility highlights the versatility and power of emulation techniques in astrophysical research.

We expect ongoing and future spectroscopic surveys targeting UFDs to have a high impact on the search for stellar halos. For example, 4DWARFS (Skúladóttir et al. 2023) will provide the chemical abundances for more than 2000 stars in ∼40 UFDs (Skúladóttir, private communication). The higher statistic will produce precise measurements of f₅ in Tucana II and possibly other UDFs, revealing the frequency of stellar halos around UFDs. Using our method on this upcoming data will allow us to uncover the past assembly history of the faintest companions of our MW.

Acknowledgments

This project received funding from the ERC Starting Grant NEFERTITI H2020/804240 (PI: Salvadori). We acknowledge the CINECA award under the ISCRA initiative for the availability of high-performance computing resources and support from the Class C project UFD-SHF, “UFD mergers: Stellar Halo Formation”, HP10CM106H (PI: Querci). We acknowledge the usage of AREPO (Weinberger et al. 2020), ASTROPY (The Astropy Collaboration 2022), DEEPXDE (Lu et al. 2021), DICE (Perret 2016), MATPLOTLIB (Hunter 2007), NUMPY (Harris et al. 2020), PYNBODY (Pontzen et al. 2013), PYTHON3 (Van Rossum & Drake 2009), SCIPY (Virtanen et al. 2020), and SCIKIT LEARN (Pedregosa et al. 2011).

References

Appendix A: Dynamical heating from dark matter resolution

Two species simulations, i.e. dark matter and stars, are affected by numerical heating that depends on the resolution mass ratio μ = m_DM/m_⋆ and on the number of particles (Ludlow et al. 2023). High μ or low number of particles induce numerical heating, resulting in spurious heating of the stellar component and a larger galaxy. To test the impact of numerical heating on our results, we compare the properties of the post-merger galaxies of three suites, varying the dark matter and stellar resolution. The dark matter particle of the low-resolution (high-resolution) suite has a mass of $m_{\mathrm{DM}}^{\mathit{low}}={10}^3{M}_{\odot }$ ($m_{\mathrm{DM}}^{\mathit{high}}={10}^2{M}_{\odot }$) and the resolution mass ratio is μ^low = 10³ M_⊙ (μ^high = 10²). The mass of a single dark matter particle in the mid-resolution suite is $m_{\mathrm{DM}}^{\mathit{mid}}=m_{\mathrm{DM}}^{\mathit{low}}$ while μ^mid = μ^high. The adopted μ^low is similar to the value used by Agertz et al. (2020), while μ^high is consistent with the value in Ricotti et al. (2022).

In Fig. A.1, we show the distribution of the post-merger galaxies properties, f₅ top row and R_⋆ bottom row, in the parameter space explored. The high-resolution simulations are in orange, the mid-resolution in blue, and the low-resolution in green. Comparing the distributions of the fraction of distant stars, we note that they overlap, with the higher resolution suite peaking at a higher value (∼13%) compared to the mid-resolution (∼8.6%) and the low-resolution (∼8.2%). On the contrary, the half-mass distributions do not overlap, with the mid-resolution R_⋆ systematically between the low- and high-resolution one.

Considering the typical error estimating f₅, i.e., of the order of 20%, the f₅ distributions at different resolutions are consistent. Such consistency confirms that the simulations capture the stellar halo formation by a dry merger, even with a resolution as low as $m_{\mathrm{DM}}^{\mathit{low}}$. On the contrary, the half-mass radius distributions do not overlap, even considering the error associated with the measure, showing the impact of spurious numerical heating and forcing a different method to test the convergence of the simulations. The size of the post-merger galaxy (orange points) is larger than the progenitor radius (the dashed line), as expected in N-body dry-merger simulations (Oogi et al. 2015) and evident from the bottom right panel of Fig. A.1. For R_1/2 > 50 pc, the two radii are proportional, while below 50 pc R_⋆ saturates at a minimum value, i.e., the radius resolution limit, suggesting that the numerical heating can be neglected in simulations with R_1/2 > 50 pc. To further verify that numerical heating is negligible, we limit R_1/2 to the proportional regime and determine the impact of each parameter on R_⋆. We find the same hierarchy as when using the whole dataset, with R_1/2 being the dominant parameter and M₁/M₂ playing a secondary role, showing that numerical heating does not affect our results. Moreover, the values of R_⋆ derived at R_1/2 ≃ 100 pc, i.e., the expected value of the massive progenitor radius consistent with observations, for different resolutions are consistent when considering the observational errors are of the order of ∼30 pc.

For the reasons stated above, our results on the stellar halo formation and size of the post-merger galaxies are well captured by the dark matter resolution $m_{\mathrm{DM}}^{\mathit{high}}={10}^2{M}_{\odot }$, and there is no need for increasing the resolution.

Fig. A.1. Comparison of the distribution of the post-merger galaxies properties, i.e., half-mass radius R⋆ (bottom row) and the fraction of stars f5 beyond 5R⋆ (top row), changing the dark matter and stellar resolution. Each column represents an IC property explored in this work, from left to right are: specific angular momentum l, specific kinetic energy k, merger mass ratio M1/M2, dark matter abundance MDM/M⋆, and progenitor half-mass radius R1/2. The dashed line in the bottom right corner is R⋆ = R1/2. The low-resolution simulations (mDM = 103 M⊙, mDM/m⋆ = 103 M⊙) are in green, the mid-resolution simulations are in blue (mDM = 103 M⊙, mDM/m⋆ = 102 M⊙), while the high-resolution simulations ($m_{\mathrm{DM}}^{\mathit{high}}={10}^2{M}_{\odot }$, mDM/m⋆ = 102 M⊙) are in orange.

Appendix B: Using neural networks as an emulator

In this work, we adopt and train a NN as an emulator to infer an approximate value of f₅ and R_⋆ of post-merger galaxies based on the ICs of the merger. The NN has an input layer with 5 neurons, one for each parameter explored, i.e., l, k, M₁/M₂, M_DM/M_⋆, and R_1/2, 5 hidden layers with 128 dense neurons, and an output layer with a neuron for f₅ and R_⋆, each.

The number of simulations available for training and testing the NN is limited, and each simulation carries information in a specific interval of ICs. To avoid losing information, we train the NN on all the available data, and we use a l2 regularization to control overfitting (Bühlmann & van de Geer 2011). We tested our approach by training nine identical emulators on 90%, 80%, and 70% of the dataset and testing them on the remaining fraction.

Figure B.1 shows the median of the inferred values of all emulators, colored by the training data size. All medians do not deviate more than 0.1 dex between each other, 0.05 dex in some cases, confirming that using all the available data to train the NN does not affect the results. Moreover, the primary and secondary dependencies on the properties of the IC and the inferred value of M₁/M₂ and R_1/2 are consistent between all emulators, corroborating that our results do not critically depend on the training set of the NN.

Fig. B.1. Same as Fig. A.1 considering only the high resolution simulation data in orange. Each line represents the median distribution of the reference emulator, black, or the test emulator, colored in shades of blue based on the size of the training set.

Appendix C: Projecting 3D properties

Since the merger orbital plane is fixed in the suite, projecting both R_1/2 and f₅ on the XY plane can have systematic effects on the measures. To minimize such effects and to account for the asymmetry of the post-merger galaxies, we project the galaxy into 8000 different planes uniformly distributed in the 3D space and measure R_1/2 and f₅ for each projection. We take the median and 16%-th and 84%-th percentiles of the resulting distributions as the value and errors, respectively. The results of the ensemble average are shown in Fig. C.1, where each solid line represents a simulation and spans the interval between the 16%-th and the 84%-th percentile, and the dashed line is the plane’s bisector. As reported by Wolf et al. (2010), the ratio between the 2D and 3D half-mass radii can be approximated to 3/4 for many profiles, which is also true in our post-merger simulations. The 3D values of f₅ are within the error bars of the projected f₅, suggesting that the parameter is negligibly affected by the projection procedure.

Fig. C.1. Distribution of the half-mass radius (left) and f5 (right) of the post-merger galaxy on the 3D-2D plane. Each solid line represents a simulation in the suite, spanning from the 16%-th to the 84%-th percentiles. The dashed lines in the two panels are the plane bisectors. The 3D half-mass value is multiplied by 3/4 to account for the projection, as reported by Wolf et al. (2010).

All Tables

All Figures

Fig. 1. Example of the virialized stellar distribution of a progenitor with M⋆ = 5.48 × 103 M⊙ and R1/2 = 118 pc generated using DICE (see Sect. 2.1.3). Left panel: Stellar density (ρ⋆) extracted by spherically averaging the mass in a simulation shell (points) and the fit exponential profile (solid line, Eq. (4)) as a function of the radius (r). Right panel: Stellar surface density map, Σ⋆ integrated along the line of sight; the spatial scale is indicated as an inset. In both panels, the dashed lines mark the location of spheres with radius R1/2, 3 R1/2, and 5 R1/2.

In the text

Fig. 2. Example of merger ICs prepared with DICE (see Sect. 2.1.4). The density map in the X − Y plane shows the total (dark matter+stars) density (ρ) that has been mass-averaged along the Z axis. The two progenitors were placed in the simulation box at a distance of D = 8 kpc along the X axis and with an impact parameter of |Y1|+|Y2| (Eq. (5b)). The arrows mark the direction of the initial velocity of each progenitor (Eq. (5c)).

In the text

Fig. 3. Parameter sampling in the dark matter (MDM) vs. stellar mass (M⋆) plane. Each pink (purple) point represents the least (more) massive progenitor in a merging simulation. The shaded blue region indicates the dark matter to stellar mass (MDM/M⋆) interval explored by the suite. The gold star shows the estimate of Tucana II from Chiti et al. (2021) observations, while the green lines are the extrapolation from the abundance matching model from Behroozi et al. (2013) at redshift z = 0 and 6. The shaded green region represents the extrapolation of the scatter in the stellar-mass to halo-mass relation measured by Kim et al. (2024) for dwarf galaxies.

In the text

	Fig. 4. Parameter sampling in the stellar mass (M⋆) vs. scale radius (R1/2) plane. Each pink (purple) point represents the least (more) massive progenitor in a merging simulation and the shaded blue area indicates the scale radius interval explored in the suite. The shaded green region represents the polynomial fit ±1σ of all the available UFD data in the Local Volume Database (Pace 2024).
In the text

Fig. 5. Example of merger evolution. The different panels report the stellar surface density (Σ⋆) maps at various evolutionary stages for progenitors with a M1/M2 ≃ 10 (profile in Fig. 1) and the merger setup in Fig. 2. The bottom right panel is a zoom-in on the virialized post-merger galaxy, and the three dashed circles have radii of one, three, and five half-mass radii. After the merging, the resulting fraction of stars in the outskirts is f5 ∼ 5.5%.

In the text

Fig. 6. Example of pre- vs. post-merging stellar distributions. We show the profile for the simulation selected for Fig. 5 by reporting ρ⋆ for the two progenitors with shades of purple and the post-merger galaxy in orange. The solid lines are the fit profile, which is exponential (Eq. (4)) and a sum of exponential for the progenitors and the post-merger galaxy, respectively. The vertical lines mark the half-mass radius for each galaxy.

In the text

	Fig. 7. Origin of the stars in the post-merger galaxy. For the simulation selected for Fig. 6, we show the radial distribution of the fraction of stars that originally belonged to the high-mass progenitor (purple) and the low-mass progenitor (pink). The dashed vertical lines are at 1, 3, and 5 R⋆.
In the text

Fig. 8. Emulation of the post-merging f5 (upper panels) and stellar radius (R⋆, lower panels) as a function of merger setup parameters, i.e., specific angular momentum (l, Eq. (6a)), specific kinetic energy (k, Eq. (6b)), mass ratio (M1/M2), dark to stellar mass ratio (MDM/M⋆), and pre-merger stellar size (R1/2, Eq. (4)). Each of the first five columns gives the distribution of outputs of the merging process as a function of a single parameter. The distribution is obtained by emulating 106 IC log-spaced distributed (see Table 1), and the gold line represents its median. The orange points are the individual simulations used to train the emulator. The red line (shaded region) gives the average (dispersion) of the observed values for Tucana II (Bechtol et al. 2015; Chiti et al. 2021). The last column shows the distribution of the emulated output normalized to the total number of inferences.

In the text

Fig. 9. Distributions of the input couple R1/2 and M1/M2 that give f5 and R⋆ simultaneously consistent with the observations from Tucana II. Similarly to Fig. 8, the red line gives the average value inferred from observations, but the ranges for f5 and R⋆ were cut to include only the observational uncertainty, i.e., f5TucII = 10 ± 5% (Chiti et al. 2021) and $R_{1/2}^{\mathrm{TucII}}=(120\pm 30)$ pc (Bechtol et al. 2015), respectively. The orange points are the post-merger galaxies consistent with observations from Tucana II, and are numbered with their simulation ID.

In the text

Fig. 10. Comparison of the observed stellar number density distribution of Tucana II and the simulations with the f5 and R1/2 parameters closest to the observed values (f5TucII = 10 ± 5% from Chiti et al. 2021 and $R_{1/2}^{\mathrm{TucII}}/\mathit{pc}=120\pm 30$ from Bechtol et al. 2015). Each panel shows a different simulation, as is indicated in the inset, ordered by decreasing χ2 value. The green line (shaded area) gives the mode (error) of the profile for Tucana II, observed by Jensen et al. (2024). The orange line (points with bars) gives the best fit (average and standard deviation of the data points) of the profile of a simulated post-merger galaxy, obtained by adopting a two-exponential profile (Eq. (11)). The solid (dashed) gray line shows the R⋆ (5 R⋆) post-merger radius.

In the text

Fig. A.1. Comparison of the distribution of the post-merger galaxies properties, i.e., half-mass radius R⋆ (bottom row) and the fraction of stars f5 beyond 5R⋆ (top row), changing the dark matter and stellar resolution. Each column represents an IC property explored in this work, from left to right are: specific angular momentum l, specific kinetic energy k, merger mass ratio M1/M2, dark matter abundance MDM/M⋆, and progenitor half-mass radius R1/2. The dashed line in the bottom right corner is R⋆ = R1/2. The low-resolution simulations (mDM = 103 M⊙, mDM/m⋆ = 103 M⊙) are in green, the mid-resolution simulations are in blue (mDM = 103 M⊙, mDM/m⋆ = 102 M⊙), while the high-resolution simulations ($m_{\mathrm{DM}}^{\mathit{high}}={10}^2{M}_{\odot }$, mDM/m⋆ = 102 M⊙) are in orange.

In the text

	Fig. B.1. Same as Fig. A.1 considering only the high resolution simulation data in orange. Each line represents the median distribution of the reference emulator, black, or the test emulator, colored in shades of blue based on the size of the training set.
In the text

	Fig. C.1. Distribution of the half-mass radius (left) and f5 (right) of the post-merger galaxy on the 3D-2D plane. Each solid line represents a simulation in the suite, spanning from the 16%-th to the 84%-th percentiles. The dashed lines in the two panels are the plane bisectors. The 3D half-mass value is multiplied by 3/4 to account for the projection, as reported by Wolf et al. (2010).
In the text