© The Authors 2026

1. Introduction

Our Galaxy hosts a significant number of dwarf galaxies with diverse properties. These systems are commonly classified as dwarf spheroidal, dwarf elliptical, and dwarf irregular galaxies. The least luminous members of this population constitute the class of ultra-faint dwarf galaxies (UFDs). UFDs are defined as dwarf galaxies with an integrated V-band magnitude (M_V) > −7.7 (L ≲ 10⁵ L_⊙) (Simon 2019). The number of known UFDs has increased significantly since the Sloan Digital Sky Survey (SDSS) began identifying such objects (Willman et al. 2005). Over the subsequent two decades, the SDSS, the Dark Energy Survey (DES), and the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS) have discovered approximately eight times more dwarf satellites of the Milky Way (MW) than were known prior to the operation of these surveys. To date, the number of MW dwarf galaxies has exceeded 80 (Simon 2019). As of today, few extragalactic sources with properties similar to those of UFDs have been discovered. We report the discovery of PegasusW (Sand et al. 2022) in the Local Group, Leo M and Leo K (Battaglia et al. 2023) just outside of the MW halo, and six others orbiting NGC 253 (Sand et al. 2021) and NGC 300 (Mutlu-Pakdil et al. 2022), three for each galaxy. Note that these were detections just above the sky background limits and need confirmation (especially the megaparsec-distance source).

In the review by Simon (2019), it is implied that UFDs are comparable in luminosity (and stellar mass) to Galactic globular clusters but that their half-light radii are at least an order of magnitude larger. This definition is consistent with that adopted by Bullock & Boylan-Kolchin (2017). The main characteristics of UFDs, as summarized by Simon (2019), are as follows: (a) the stellar kinematics of UFDs indicate the presence of large dark matter (DM) fractions; (b) all dwarf galaxies, with the exception of the lowest-luminosity UFDs, are larger than any known globular cluster; (c) within individual UFDs, the abundances of Fe and α-elements exhibit substantial internal spreads, indicative of extended star formation and self-enrichment; (d) UFDs follow a luminosity–metallicity relation (Hidalgo 2017), whereas globular clusters do not; and (e) the abundance patterns of certain elements in UFDs resemble those observed in brighter dwarf galaxies and differ from the light-element abundance correlations characteristic of globular clusters.

In this work we focus primarily on point (a). Velocity dispersion measurements are available for a limited number of MW UFDs (26 of the 42 candidates reported in Simon 2019), as well as for Pegasus IV (reported in Cerny et al. 2023), yielding a total sample of 27 systems. We adopted the updated measurements of Chiti et al. (2022) for Grus I. The resulting mean velocity dispersion is ⟨σ⟩ = 4.66 km s⁻¹. In Simon (2019), the dynamical mass enclosed within the three-dimensional half-light radius (R_1/2), denoted M_1/2, was computed following Wolf et al. (2010) as

$$\begin{array}{c}\hfill M_{1/2}=930{\sigma }_{1\mathrm{D}}^2{R}_{1/2},\end{array}$$(1)

where σ_1D is the one-dimensional velocity dispersion in km s⁻¹ and R_1/2 is the projected two-dimensional half-light radius. The derived mean mass-to-light ratio for the UFD sample is ⟨M_1/2/L_V⟩ = 459.95 M_⊙/L_⊙, an extremely high value, approximately an order of magnitude higher than that of classical dwarf galaxies (M_V < −7.7), which have an average M/L ≈ 48.83 M_⊙/L_⊙. A similar trend is observed for M31 satellites, with ⟨M/L⟩ = 449.4 M_⊙/L_⊙ for UFDs and ⟨M/L⟩ = 55.34 M_⊙/L_⊙ for brighter dwarfs.

These results mirror earlier findings for Local Group dwarf spheroidal galaxies (dSphs; Amorisco & Evans 2011). These systems, typically located at Galactocentric distances exceeding 70 kpc and composed of old, metal-poor stellar populations, are therefore structurally similar to Galactic globular clusters. They exhibit no clear evidence of global rotation and are thus not rotationally supported, implying that mass estimates must rely on spectroscopic measurements of their velocity dispersions (Mateo et al. 1993).

Several studies of classical dSphs in the MW halo (e.g., Fornax, Sculptor, Ursa Minor, Draco, Leo I, Leo II, Sagittarius, Sextans, and LGS 3) have shown that the observed velocity dispersion, σ_obs, is significantly higher than the value expected if these systems were simple scaled-up versions of globular clusters, for which σ ∼ 1 − 3 km s⁻¹ (Mateo 1997). Subsequent work by Simon & Geha (2007) confirmed observed dispersions in the range 3.3 − 7.6 km s⁻¹, challenging earlier claims of an upper limit σ_obs ≤ 7 km s⁻¹ (Wyse et al. 2008) and motivating revisions of the inferred mass scale of these systems.

Several scenarios have been proposed to explain the unusually large velocity dispersions observed in dSphs and UFDs. Early studies (e.g., Aaronson et al. 1988) invoked the presence of substantial DM halos, consistent with predictions of the standard cosmological model. Alternative interpretations suggested that dSphs might be out of virial equilibrium due to ongoing tidal disruption. However, the importance of Galactic tides has been questioned based on the observed luminosity–metallicity relation (Kirby et al. 2008) and the lack of clear observational signatures of tidal features in many systems. For example, Geha et al. (2009) found no evidence of tidal tails or induced rotation in Segue I, rejecting the hypothesis that it is a disrupted globular cluster formerly associated with the Sagittarius stream (Belokurov et al. 2007).

Another possibility is that the elevated velocity dispersions arise from binary-star orbital motion. Although the tidal disruption scenario has largely been disfavored, the role of binaries remains an active topic of investigation. Mateo (1997) argued that unresolved binaries are unlikely to fully account for the inflated dispersions observed in classical dSphs, which are therefore still considered DM-dominated. However, their impact may be non-negligible in UFDs, the low-luminosity counterparts of classical dSphs (McConnachie & Côté 2010; Spencer et al. 2018). Despite recent expansions of the UFD sample (Massari & Helmi 2018), small-number statistics and the lack of extensive multi-epoch spectroscopic observations continue to hinder robust constraints on binary fractions and orbital parameter distributions (McConnachie & Côté 2010; Spencer et al. 2018). In this context, it is worth mentioning the case of Segue II, for which the inflation of the velocity dispersion has been extensively debated (Belokurov et al. 2009; Kirby et al. 2013). Unfortunately, only in a limited number of cases do the available spectroscopic data allow meaningful constraints to be placed on the binary fraction (e.g., the UFD galaxy Reticulum II; Minor et al. 2019). Other recent studies also suggest the importance of wide binaries in UFDs (El-Badry & Rix 2019), and a resolved stellar population will allow us to place constraints on the binary fractions and eventually the DM abundance (Mutlu-Pakdil et al. 2023). For this reason, modeling approaches based on Monte Carlo simulations and Bayesian analyses have been widely adopted. To date, most models have aimed to reproduce the observed velocity dispersion of classical dSphs by varying both the binary fraction and the binary orbital parameters, and subsequently comparing the results with spectroscopic data to estimate the contribution of binaries to the observed velocity dispersion, σ_obs, in UFDs (Spencer et al. 2017; Massari & Helmi 2018). Nevertheless, assumptions regarding the distributions of orbital parameters – particularly orbital periods and semimajor axes – represent a significant limitation in this framework. This motivates the development of theoretical models capable of making inferences about the binary populations in these systems that are as reliable and general as possible. Assuming that both dSphs and UFDs are DM-dominated systems, which is currently the most widely supported scenario, we present in this paper an N-body study aimed at exploring the effects of different choices of binary orbital parameters using different binary fractions on the observed velocity dispersion of such galaxies.

The ultimate purpose of this work is to investigate the impact of binaries on the determination of the dynamical mass in the faintest MW satellites, with particular reference to Rastello et al. (2020) regarding the methodology used to compute velocity dispersions. Finally, recent studies (Rostami-Shirazi et al. 2025) suggest that baryonic explanations may be possible in some cases, although additional components (e.g., a central black hole) may be required to match the inferred dynamical mass in specific systems.

This paper is organized as follows. In Sect. 2 we describe the methodology and the simulation set-up. In Sect. 3 we present our results, and in Sect. 4 we summarize our main conclusions.

2. Methodology and initial conditions

Of the dwarf galaxies known in the local Universe, UFDs are the faintest. They are satellites of the MW that populate the lowest-luminosity tail of the galaxy luminosity function (L ≲ 10⁵ L_⊙) and are commonly considered to be the oldest, most metal-poor, least chemically evolved, and most DM-dominated stellar systems known.

The distinction between dSphs and UFDs is somewhat arbitrary. Here, we adopted the classification proposed by Bullock & Boylan-Kolchin (2017) and subsequently used by Simon (2019), defining UFDs as systems with integrated magnitudes M_V ≥ −7.7 (L ≲ 10⁵ L_⊙).

In this paper we present a set of direct N-body simulations of a prototype UFD composed of both single stars and binaries. The dynamical evolution of each model was followed with high precision using the NBODY6++GPU code, from time zero (after gas dissipation) up to one Hubble time.

The main goals of this paper are: (i) to investigate the long-term evolution of a UFD-like system in the presence of low, intermediate, and high binary fractions, thereby testing the assumption that such systems are effectively collisionless; and (ii) to quantify how binaries affect both the internal dynamics of the galaxy and the observational determination of its virial mass. This latter aspect extends previous studies (e.g., El-Badry et al. 2018; Rastello et al. 2020; Pianta et al. 2022), in which the impact of binaries on mass estimates was explored using statistical or semi-analytical approaches.

It is well established that the binary fraction increases toward lower metallicities (Moe et al. 2019). Several binaries have been identified in UFDs (e.g., Spencer et al. 2018), and observational evidence suggests that the binary fraction may range from a few percent up to approximately 70%. Combined with the growing evidence that binaries can significantly inflate the measured velocity dispersion in faint galaxies, this motivates a detailed investigation of their dynamical role, especially given the claimed extreme DM content of these systems (Irwin & Hatzidimitriou 1995).

Binary-induced velocity dispersion inflation has been observed, for example, in Very Large Telescope studies of Boötes I, Leo IV, and Leo V (Jenkins et al. 2021). In Boötes I and Leo IV, the one-dimensional velocity dispersion decreases by approximately 0.2 − 0.3 km s⁻¹ once the spurious binary contribution is removed. However, we believe that due to the small-number statistics, observational data alone remain insufficient to robustly quantify the binary contribution to σ, which enters quadratically into virial mass estimates.

For this reason, numerical approaches have been adopted to assess the influence of binaries on velocity dispersion measurements. Rastello et al. (2020) conducted a full, high-precision N-body study of low-mass stellar systems comparable to open clusters, while Pianta et al. (2022) employed a statistical framework to study systems with sizes comparable to dSphs and UFDs. These authors found that, particularly in the faintest galaxies, binary-induced inflation of σ can be substantial, potentially leading to virial mass overestimates of a factor of two or more for binary fractions of ∼40% when multi-epoch observations are unavailable. Actually, multi-epoch are extremely rare and limited for UFDs (e.g., Kirby et al. 2015; Spencer et al. 2017, 2018; Jenkins et al. 2021; Buttry et al. 2022; Simon et al. 2011, 2015).

The full dynamical evolution presented in this paper should shed light on the debated question of the selective effect of the binary content of a UFD on the overall dynamics. Although UFDs are commonly treated as collisionless systems, the presence of a significant binary population complicates this picture. Interactions among different stellar components can alter the binary population, preferentially disrupting wide binaries while leaving hard binaries intact. This evolution has direct consequences for mass estimates derived from velocity dispersion measurements. For the evolution of wide binaries in dwarf galaxies, we refer to Peñarrubia et al. (2016) and Livernois et al. (2023).

Our data analysis follows the methodology described in Pianta et al. (2022) and includes several definitions of the velocity dispersion: (1) a global dispersion, σ_g, computed by treating all stars as if they were single; (2) a dispersion σ_g, cm, in which binaries contribute only through their center-of-mass velocities; and (3) a luminosity-weighted dispersion, σ_L, in which single stars are treated normally and binaries are weighted according to

$$\begin{array}{c}\hfill v=\frac{L_A{v}_A+L_B{v}_B}{L_A+L_B},\end{array}$$(2)

where L_A, v_A, and L_B, v_B are the luminosities and speeds of the binary components (see Rastello et al. 2020).

The uncertainty in our work arises from both systematic and statistical sources. Systematic uncertainties include numerical integration errors as well as uncertainties in the photometric output derived from the mass–luminosity relations implemented in the code. For a detailed discussion of these effects, we refer to Kamlah et al. (2022a).

2.1. UFD initial conditions

Table 1 summarizes the settings of our UFD models, which are largely inspired by Pianta et al. (2022). We sampled a Plummer model (Plummer 1911) with 82 000 stars, adopting different binary fractions f_b = N_b/N, where N_b is the number of binary systems and N is the total number of stars. We explored values of f_b from 0 to 0.5, in steps of 0.1.

The system was initially contained within a virial radius of r_vir = 50 pc. Stellar masses were drawn from a Kroupa (2001) initial mass function over the range 0.08 − 150 M_⊙, yielding a total initial stellar mass of 6.24 × 10⁴ M_⊙. The model was initialized in virial equilibrium and evolved dynamically for a time comparable to one Hubble time (13.7 Gyr).

Ultra-faint dwarf galaxies are generally regarded as collisionless stellar systems due to their low densities and their characteristic radii, which can be either small or relatively large. This, despite the low velocity dispersion, implies very long two–body relaxation times. The local relaxation time can be estimated as (Binney & Tremaine 2008, Eq. 7.106)

$$\begin{array}{c}\hfill t_{\mathrm{rel}}=\frac{0.34{\sigma }^3}{G^2{m}_{\ast }\rho ln\mathrm{\Lambda }},\end{array}$$(3)

where σ and ρ are the local one-dimensional velocity dispersion and mass density, m_* is the mean stellar mass, and lnΛ is the Coulomb logarithm. The Coulomb logarithm represents the cumulative effect of numerous weak gravitational encounters between stars. Formally, it arises from integrating over the range of impact parameters, b, that contribute to the scattering between two bodies, and is defined as lnΛ = ln(b_max/b_min). The lower cutoff, b_min, is typically taken to be on the order of the impact parameter corresponding to a 90° deflection, b₉₀ = 2Gm_*/σ², while the upper cutoff, b_max, is set by the characteristic size of the system, usually the local scale length of the density distribution or the half-mass radius (e.g., Spitzer 1987; Binney & Tremaine 2008). Because the exact values of these cutoffs are somewhat ambiguous and system-dependent, lnΛ is often approximated as a constant in the range 5−15, depending on the stellar system under consideration. For globular clusters, lnΛ ≃ 10 is commonly adopted, as it yields accurate order-of-magnitude estimates of the relaxation time. Thus, following standard practice, we adopted lnΛ = 10.

Under this assumption, we estimate a half-mass relaxation time of t_rel, hm ≈ 43 000 Myr and a core relaxation time of t_rel, c ≈ 6000 Myr. Here, the core radius is defined as the three-dimensional radius at which the stellar density drops to half its central value. Although the relaxation time exceeds a Hubble time in most regions, the innermost regions experience relaxation on timescales shorter than the simulation duration. Consequently, we expect mild but non-negligible collisional effects to develop over the course of the evolution.

Binary systems segregate on the dynamical friction timescale (e.g., Binney & Tremaine 2008),

$$\begin{array}{c}\hfill t_{\mathrm{df}}=0.66\frac{m_{\ast }}{m_{\mathrm{bin}}}{t}_{\mathrm{rel}}\approx 0.33t_{\mathrm{rel}},\end{array}$$(4)

assuming a typical binary mass m_bin ≈ 2m_*. Figure 1 shows the relaxation time, t_rel, and the dynamical friction timescale, t_df (discussed below), as functions of radius for our model. The relaxation time greatly exceeds a Hubble time in the outer regions of the system.

Fig. 1. Initial relaxation (blue) and dynamical friction (orange) times of our UFD Plummer model as functions of radius. The vertical dashed red and black lines denote the initial core and half-mass radius of the adopted Plummer model. For reference, the horizontal gray line marks the Hubble time.

For the bulk of the system (r ≲ r_h), this confirms the validity of the collisionless approximation. However, in the innermost regions (r ≲ r_c), the local relaxation time is shorter than 13 700 Myr, implying that the core experiences approximately two relaxation times by the end of the simulation. As a result, we expect small but non-negligible collisional effects to develop over the course of the simulation, such as mild core contraction or mass segregation in the central regions of the UFD. In particular, binary systems are expected to segregate on the dynamical friction timescale.

2.2. Binary initial conditions

Binary pairing is assigned randomly following Weidner et al. (2009) and Wang et al. (2015). The logarithm of the semimajor axis is drawn from a uniform logarithmic distribution over the range 0.1−50 AU, ensuring that the primordial population is dominated by hard binaries, following the conclusion on the contribution of hard binaries to velocity dispersion of Pianta et al. (2022). Eccentricities are drawn from a thermal distribution, f(e) = 2e. Since the 0% and 50% models are primarily intended as proof-of-concept cases, we adopted a 30% binary fraction for most of our plots. This choice is motivated by several considerations: (i) the system is effectively collisionless, with negligible direct stellar collisions (a point that will be demonstrated by our results below); (ii) binary stars play a critical role in shaping the measured velocity dispersion; and (iii) binaries serve as an important proxy for stellar evolution, through processes such as Roche lobe overflow and other binary evolutionary channels. The models are labeled following the nomenclature reported in Table 1.

2.3. NBODY6++GPU

High-precision simulations of large stellar systems are performed using NBODY6++GPU, a direct N-body integrator optimized for GPU acceleration (Nitadori & Aarseth 2012; Wang et al. 2015) and parallel computing (Spurzem & Kamlah 2023). The code is based on the original N-body family developed by Sverre Aarseth, able to integrate systems with up to more than one million particles on modern supercomputing architectures (see, Wang et al. 2016; Arca Sedda et al. 2023, 2024a,b).

The NBODY6++GPU code uses Kustaanheimo-Stiefel regularization (Kustaanheimo & Stiefel 1965), to efficiently handle close binaries. Single and binary stellar evolution is also treated with the “Binary Stellar Evolution” and “Single Stellar Evolution” packages of Hurley et al. (2000, 2001, 2002), which have been refined and adapted to the current discoveries in stellar evolution (for more information, see Kamlah et al. 2022a; Arca Sedda et al. 2024a, 2023, 2024b). Although NBODY6++GPU is most commonly applied to star clusters, it is well suited for modeling low-density systems such as UFDs. In the present work, DM is not included, but its incorporation through analytic density profiles will be explored in future studies.

3. Results

3.1. Dynamical evolution of the UFD

The evolution of the Lagrangian radii of the UFD shown in Fig. 2 illustrates the behavior of both the single-star and binary components. With the exception of model U0 – which contains no binaries – all models exhibit a similar evolutionary trend in their binary populations. In particular, in the innermost regions, models with lower binary fractions show slightly earlier mass segregation, since binaries tend to delay core collapse (see Appendix A, where we display and compare the Lagrangian radii of models U0, U20, and U40).

Fig. 2. Lagrangian radii evolution of the UFD model U30 over one Hubble time, using the time-updated value of the total mass. Top: For both single and binary stars. Middle: For single stars. Bottom: For binary stars.

Figure 2 highlights the nearly collisionless nature of UFDs. Although the system loses mass from its most massive stars – which evolve into black holes or neutron stars within the first 100 Myr – the Lagrangian radii are only mildly affected, and only in the innermost regions. This behavior is a consequence of the low stellar density, even at small radii, which results in rare stellar encounters. The outer regions also expand slowly, since the UFD is treated as an isolated system.

Although the cluster is considered isolated, we adopted a mass-loss criterion defined by r_esc = 20 r_scale, where r_scale is the scale radius of the system, which we set equal to the half-mass radius. Thus, if the UFD expands to 100 pc, the escape radius becomes 2000 pc. All stars that move beyond r_esc are considered lost from the system.

The UFD does not show any significant structural variation until approximately 3000 Myr, when bumps appear in the inner Lagrangian radii and the outer shells expand. This behavior is a consequence of the onset of mass segregation, which increases the encounter rate in the central regions and triggers a mild core contraction. As stellar evolution proceeds, white dwarfs form, and the associated mass loss halts further core contraction once the majority of red giants have evolved into white dwarfs.

We therefore conclude that, prior to the onset of mass segregation, binary disruptions and stellar ejections are limited to processes related to binary evolution (e.g., Roche lobe overflow). This outcome is also a consequence of our choice of initially hard binaries, which is well motivated by the results of Pianta et al. (2022) and appropriate for the scope of this study.

The middle panel of Fig. 2 shows the evolution of the Lagrangian radii for the single-star population of model U30. Since single stars constitute 70% of the total stellar population in this model, they play a dominant role in determining the long-term evolution of the UFD, which remains largely collisionless. After 3000 Myr, the positions of single stars are more strongly perturbed in the inner regions as a result of mass segregation, which drives lower-mass objects toward larger radii.

The bottom panel of Fig. 2 shows a similar pattern for the Lagrangian radii of the binary population. Both single stars and binaries are driven toward the outer regions, but the effect is more pronounced for binaries. In particular, the 99% Lagrangian radius of the binary component increases rapidly, whereas the corresponding single-star shell expands more steadily. This behavior is induced by mass segregation and core contraction, which preferentially drive low-mass single stars and low-mass binaries outward. As a result, a prominent expansion peak appears around the time when the system ceases segregating, at approximately 5 Gyr for model U30.

By 3 Gyr, the UFD has already expanded beyond 1000 pc, a scale at which tidal stripping by the MW would become significant in a more realistic, non-isolated scenario. However, since the UFD is modeled as isolated in our simulations, it does not experience evaporation in the denser regions until several relaxation times have elapsed, well beyond a Hubble time.

As discussed above, the key feature revealed by these plots is the presence of mass segregation, which is quantified by the average stellar mass within the Lagrangian shells and shown in the top panel of Fig. 3. Mass segregation remains negligible up to approximately 1 Gyr, after which a pronounced segregation of higher-mass stars develops in the central regions. Nevertheless, its impact on the overall dynamical evolution of the inner regions is limited, owing to the large spatial extent of the UFD at the time segregation becomes effective, which results in long segregation timescales.

Fig. 3. Average mass within the different Lagrangian radii of the U30 model over one Hubble time.

The formation of new binaries – amounting to only ∼0.004% of the total stellar population – does not significantly affect the evolution of the Lagrangian radii of either the single-star or binary components. As we show in the next section, these newly formed systems are wide binaries that are rapidly disrupted (within ≲3 Myr), rendering their contribution to the global dynamical and luminosity evolution of the UFD entirely negligible.

3.2. Mass loss and binaries lost in the UFD

The total mass ejected is extremely small compared to that of denser systems, as shown in the top panel of Fig. 4. The total mass lost amounts to 9% of the initial mass. Notably, this loss is driven primarily by stellar evolution (∼6%, corresponding to 66% of the total mass loss) and, to a lesser extent, by ejections (∼3%, corresponding to 33% of the total mass loss). Approximately 6% of the mass loss originates from single stars, while the remaining 3% comes from binary systems, with percentages computed relative to the initial total mass of the UFD.

Fig. 4. Top: Fraction of mass present in the UFD, compared to the total mass of each component over the simulation time. Bottom: Evolution of the members of dynamical binaries the survive in the UFD (blue, right-side ordinate) and the fraction of primordial binaries that survive in the UFD (orange, left-side ordinate).

The first episode of mass loss in the UFD occurs during the early evolutionary phase, when the most massive stars rapidly evolve into neutron stars or black holes within the first ∼25 − 100 Myr. In star cluster systems, neutron stars are typically ejected (e.g., Flammini Dotti et al. 2025); however, the large spatial extent of the UFD results in a negligible amount of mass loss associated with neutron-star ejections, which instead occur much later in the evolution (see the next section). This behavior is also a consequence of treating the system as isolated in our models.

As noted above, less massive objects tend to migrate toward the outer regions due to two-body relaxation after approximately 3000 Myr. In the bottom panel of Fig. 4, we show the fraction of disrupted and ejected primordial binaries, together with the absolute number of dynamically formed binaries. The majority of dynamically formed binaries are extremely wide systems (with semimajor axes a > 10³ AU) and are therefore disrupted on short timescales compared to hard primordial binaries. As is well known, most dynamically formed binaries become bound at large separations from their host stars. In our simulations, these dynamically formed binaries are consequently very short-lived.

As a result, the dynamical binary population contributes negligibly to both the luminosity and the velocity dispersion of the system, accounting for only ∼0.004% of the total number of binaries. In contrast, primordial binaries are lost at the ∼10% level before 3000 Myr. Most of these losses occur within the first 100 Myr due to stellar evolution, particularly the formation of compact objects, which can unbind the binary system. The remaining ∼6% of primordial binaries are lost between 3000 Myr and a Hubble time. As the UFD continues to expand, stellar crossing times increase significantly, making further ejections progressively less likely.

By a Hubble time, the UFD has lost approximately 16% of its original primordial binary population. Of these lost binaries, ∼3% are removed through ejections, while ∼13% are lost through binary disruption.

Finally, we examined the evolution of the semimajor axis of binaries in model U30, both including dynamically formed binaries (shown in blue) and considering only primordial binaries. The average semimajor axis in both cases remains similar during the first ∼100 Myr, after which it gradually increases. This evolution is driven by stellar evolution, particularly the formation of white dwarfs and giant stars, which can significantly alter binary orbital parameters. The peak observed at the end of the curve is caused by the inclusion of dynamically formed binaries; when these systems are excluded, the evolution of the semimajor axis is considerably smoother (see Fig. 5).

Fig. 5. Evolution of the average semimajor axis for the binaries in model U30. The orange curve corresponds to the sample of primordial binaries; the blue curve also takes the dynamical binaries into account.

3.3. Impact of stellar evolution

The role of stellar evolution is extremely important in a nearly collisionless system. For this reason, a detailed examination of the stellar population content of the simulated system is particularly informative. Table 2 reports the fraction (in percent) of the different stellar types populating the U30 UFD model over its evolution, from t = 0 to one Hubble time. The fraction of helium stars is sufficiently small that their overall contribution is negligible; consequently, they are not included in Table 2.

We considered the following specific evolutionary timestamps:

1. The 25 Myr timestamp is relevant because many massive stars begin rapid evolution and become neutron stars or black holes;

2. The 100 Myr timestamp corresponds to the formation of the first white dwarfs, primarily in binary systems;

3. The 1000 Myr timestamp approximately marks the onset of mass segregation, when the average stellar mass in the inner regions begins to increase (see Fig. 3);

4. The 5000 Myr timestamp corresponds to a phase in which the number of white dwarfs has nearly doubled relative to the previous timestamp;

5. The 13 700 Myr timestamp marks the end of the simulation, corresponding to one Hubble time.

At a Hubble time, white dwarfs constitute approximately 11% of the initial total number of stars and about 13% of the stellar population remaining at the end of the simulation. In Fig. 6 we show the number of stars of each stellar type as a function of time, providing an overview of the stellar population evolution in the UFD. One of the most prominent features is the peak in the formation of giant stars and the rapid increase in the number of white dwarfs, which become the dominant compact object population after ∼100 Myr. Both of these characteristics were first identified in dwarf galaxies by Angeletti et al. (1980).

Fig. 6. Number of stars in different evolutionary stages over time (see Table 2).

Black holes are largely retained by the galaxy, whereas neutron stars escape in significant numbers due to supernova natal kicks (see stellar evolution recipe C in Kamlah et al. 2022a). However, neutron-star escape occurs predominantly at later times, which can be explained by the crossing timescale of the system, t_cr ≈ 62 Myr. The fastest neutron stars escape first, followed by progressively slower ones. White dwarfs, by contrast, are almost entirely retained by the UFD.

In Fig. B.1 we show the Hertzsprung–Russell diagram of the U30 model in the L_bol − T_eff plane at several key evolutionary times, including 2 Gyr to capture the system during the onset of mass segregation. The HR diagrams are generated using the stellar and binary evolution prescriptions of Hurley et al. (2005) and their updated implementations (Kamlah et al. 2022b; Spurzem & Kamlah 2023; Vergara et al. 2025).

Figure 7 shows the luminosity evolution of the different stellar populations in the UFD. Black holes are excluded, as their luminosity contribution is negligible. The most massive stars evolve into neutron stars or black holes within the first 20−25 Myr, forming compact remnants with low luminosities. As shown also in Figure 8, main-sequence stars dominate the luminosity budget during the first 100 Myr, while red giant branch stars become comparable contributors at ∼1000 Myr and subsequently dominate the luminosity of the UFD. This occurs as stars near the upper end of the main sequence ignite hydrogen shell burning and evolve into red giants, leading to a substantial increase in the total luminosity of the system.

Fig. 7. Total bolometric luminosity of stars in different evolutionary stages over time (model U30; see Table 2).

Fig. 8. Time evolution of the total bolometric luminosity (blue) of model U30 and of the luminosity of single stars (orange) and binaries (green).

Blue stragglers are extremely rare in our simulations, with only a few forming, consistent with previous observational findings (Momany 2015). The luminosity contribution of neutron stars remains entirely negligible, even at their peak luminosities of ∼10¹ L_⊙. By the end of the simulation, white dwarfs constitute the third most significant contributors to the total luminosity. Since not all main-sequence stars have evolved into white dwarfs by a Hubble time, main-sequence stars continue to contribute more strongly than white dwarfs even at late times.

In summary, the luminosity of the UFD is ultimately dominated by red giant stars. As the system begins to segregate, an increasing number of red giants form, enhancing their contribution to the total luminosity. At this stage, red giant stars become approximately five times more luminous than main-sequence stars, significantly shaping the integrated light of the UFD.

3.4. Velocity dispersions of UFD models

As anticipated in the Introduction, one of the primary motivations of this study is to assess the role of binary stars in estimating the virial mass of low-mass stellar systems, such as dSphs and UFDs. In this way, observational measurements of this fundamental quantity can be more accurately interpreted. Although the effect of binaries in inflating the velocity dispersion of stellar systems has already been explored analytically (see, e.g., Pianta et al. 2022), in this work we present, for the first time, a fully self-consistent dynamical N-body model of a UFD-sized galaxy.

In particular, we simulated the evolution of a UFD-like system from the initial burst of star formation up to a Hubble time, adopting a 1:1 star–particle representation. Compared to previous studies, this approach allows us to follow the complete dynamical evolution of the system while treating binary stars as living components of the stellar population, rather than as a static or statistical correction.

To evaluate the velocity dispersion – which represents both the true kinetic support against gravity and the spurious contribution from binary orbital motions – we considered different prescriptions for computing velocity dispersions in our simulated UFDs. We refer the reader to Rastello et al. (2020) and Pianta et al. (2022) for a detailed discussion of these approaches, which are designed to reproduce, with varying degrees of realism, quantities comparable to those measured in spectroscopic observations. In this study, we adopted both a direct evaluation of the velocity dispersion and a luminosity-weighted one.

To disentangle the contribution from binary orbital motion, we compared: (i) the true virial velocity dispersion, which provides the kinetic support against gravitational contraction and, for binaries, accounts only for the center-of-mass motion (thereby excluding internal orbital motion); and (ii) a blind determination, in which all stars are treated as single objects, implicitly including the contribution from internal binary motions. The latter clearly represents an overestimate of the true kinetic support against gravity. For the former, we also included the binary contribution from the weighted luminosity distribution of binaries. We can describe the classical equation for velocity dispersion as

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{g}}=\sqrt{\frac{\sum _{i=1}^N{(v_i-⟨v⟩)}^2}{N}},\end{array}$$(5)

where i is the single star component, while ⟨v⟩ is the average velocity of all the stars. We treated each star as a single star, marking the blind approach. In order to analyze the binary center of mass and luminosity weighted approach, we changed the factor v_i for the stars in binaries. For the center of mass binaries, we wrote this as

$$\begin{array}{c}\hfill v_{\mathrm{cm},\mathrm{j}}=\frac{m_{A,j}{v}_{A,j}+m_{B,j}{v}_{B,j}}{M_j},\end{array}$$(6)

where A and B are the binary components indexes, while M_j = m_A, j + m_B, j is the sum of the two masses. A similar approach is done for the luminosities:

$$\begin{array}{c}\hfill v_{\mathrm{L},\mathrm{j}}=\frac{L_{A,j}{v}_{A,j}+L_{B,j}{v}_{B,j}}{L_j},\end{array}$$(7)

where instead of the mass, we use the luminosity. Thus, we obtain

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{g},\mathrm{cm}}=\sqrt{\frac{1}{N_s+N_b}({\displaystyle \sum _{i=1}^{N_s}}{v}_i^2+{\displaystyle \sum _{j=1}^{N_b}}{v}_{\mathrm{cm},j}^2)},\end{array}$$(8)

where N_b is the number of binaries. Similarly for the luminosity-weighted case:

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{L}}=\sqrt{\frac{1}{N_s+N_b}({\displaystyle \sum _{i=1}^{N_s}}{v}_i^2+{\displaystyle \sum _{j=1}^{N_b}}{v}_{\mathrm{L},j}^2)}.\end{array}$$(9)

The early evolution of both σ_g and σ_g, cm in the top panel of Fig. 9 is dominated by stellar evolution. The most massive stars evolve into neutron stars and black holes within the first 20−25 Myr. The natal kicks imparted to neutron stars significantly increase their velocities, producing a sharp rise in the velocity dispersion. This effect is visible in both σ_g and σ_g, cm, since both measures are sensitive to changes in stellar velocities.

Fig. 9. Model U30. Top: σg (blue), σg, cm (orange), and σL (green) evolution over time. Bottom: σL for models with initial binary fractions in the 0 − 50% range.

At later times, around 3000 Myr, the velocity dispersion increases again as a consequence of core contraction associated with the onset of mass segregation. The increased encounter rate in the central regions raises the velocities of stellar members, leading to a modest enhancement of the dispersion.

In contrast, the luminosity-weighted velocity dispersion, σ_L, is much less sensitive to mass loss and remains nearly constant throughout the simulation, aside from a mild perturbation around 3000 Myr. This behavior reflects the dominant contribution of red giant branch stars to the total luminosity during this phase. Once the number of red giants declines, σ_L returns to its previous value. Overall, σ_L is effectively time-independent, with variations limited to at most ∼1% of its initial value. In model U0, which contains no primordial binaries and only a negligible number of dynamically formed binaries, variations in σ_L are correspondingly minimal.

In the bottom panel of Fig. 9, the velocity dispersion remains approximately constant in time for all models, consistent with the largely collisionless nature of the system and the relatively stable binary fraction over the simulation.

The simulation of model U50 was not completed to a full Hubble time due to computational constraints on the LEONARDO booster system. However, this does not affect our conclusions, as the velocity dispersion in this model clearly reaches a stable plateau, consistent with the behavior observed in the fully evolved models.

As expected, the velocity dispersion increases significantly with increasing binary fraction, reaching values up to ∼10.58 times larger than those measured in the absence of binaries. The data shown in Fig. 9 confirm the expected approximately linear dependence of σ² on the binary fraction f_b. The best-fitting linear relation, which remains nearly unchanged over time, is (for f_b ≥ 0.1)

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{L}}^2=201.39f_b-12.69{\mathrm{km}}^2{\mathrm{s}}^{-2},\end{array}$$(10)

with an rms scatter of 14.80. A better fit, still for f_B ≥ 0.1, is obtained with a quadratic formula

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{L}}^2=768.0f_b^2-218.321f_b+26.0{\mathrm{km}}^2{\mathrm{s}}^{-2},\end{array}$$(11)

giving an rms scatter of 4.345. Since the virial mass estimate scales with the square of the velocity dispersion, this implies an overestimate of the virial mass by a factor of ∼27. For completeness, we also present a zoomed-in view of σ_L for each model in Fig. 10. The variations are small in all cases, consistent with the collisionless character of the UFD. In model U0, the absence of binaries leads to only a minor change in σ_L during the onset of mass segregation. This behavior further reflects the modest overall mass loss experienced by the system, confirming that models with lower numerical noise provide an accurate representation of UFD dynamics. Finally, a visual comparison with Fig. 2 of Simon (2019) of ours models is present in Appendix C. Our models covers the observed range from Simon (2019) in the velocity dispersion quite well. Note that our simulated UFDs are on the high luminosity edge of what Simon (2019) classified as UFDs (dwarfs with M_V > −7).

Fig. 10. Velocity dispersions (σL) for all models, from U0 (top left) to U50 (bottom right).

4. Conclusions

We computed six full N-body models representing UFDs populated by both single and binary stars, spanning a range of binary fractions, as described in Sect. 2.2. The main characteristics of the models are summarized in Table 1. Each model was evolved for a Hubble time following the initial burst of star formation, which allowed us to follow the complete dynamical evolution of the systems.

Ultra-faint dwarf galaxies are generally expected to be quasi-collisionless systems. We consistently find a broadly similar global dynamical evolution across all models, with differences arising primarily from the varying binary fractions and from stellar evolution. In particular, the luminosity-weighted velocity dispersion – which is the quantity most relevant for observational mass estimates of dwarf galaxies – is significantly affected by the fraction of binaries present in the system.

In the following, we summarize the main conclusions of this study:

1. The long-term evolution of a UFD can be divided into two main phases. During the first phase, lasting until approximately 3000 Myr, the system remains nearly stationary. In this early stage, some neutron stars and black holes are ejected following supernova explosions, but these events have only a minor impact on the global dynamical evolution. Afterward, the system shows a weak collisional behavior, reaching partial relaxation at around 5000 Myr. Mass segregation develops primarily in the most compact central regions (within ∼10 pc).

2. The role of binaries in the long-term evolution of the system is mainly to moderate the contraction of the inner regions. Larger binary fractions tend to delay the onset of mass segregation, leading to a slightly more extended core structure.

3. The final state of the simulated systems is characterized by a luminosity budget dominated by red giant stars, in agreement with observational evidence for UFDs (Simon 2019), although they make up, in number, less than 1% of the total.

4. With our detailed numerical simulations, we confirm earlier statistical indications from Pianta et al. (2022) that the presence of unresolved binaries significantly inflates the estimates of the system velocity dispersion, thereby leading to a substantial overestimate of the virial mass. These results reduce the need to invoke large amounts of DM in very small, low-luminosity systems such as UFDs.

In conclusion, our study demonstrates that UFDs are not perfectly collisionless systems and that their dynamical evolution is nontrivial, being partly governed by their binary content, which remains significant even after a Hubble time. Moreover, we validate previous statistical estimates of the impact of internal binary motions on observational determinations of galaxy velocity dispersions. In particular, sufficiently hard binaries survive throughout the evolution of the galaxy, and their rapid internal motions effectively inflate the Doppler-derived velocity dispersion. The next step of this research would be to include a possible halo of DM, to test its role in the overall dynamics of these small galaxies.

Acknowledgments

This work was possible thanks to the use of the CINECA Leonardo booster, in the frame of the project “The Case of Dark Matter Overabundance in Small Galaxies”, which was awarded 80 000 core hours. This material is based upon work supported by Tamkeen under the NYU Abu Dhabi Research Institute grant CASS. R.S. acknowledges support from NAOC International Cooperation Office in 2023, 2024, and 2025, Chinese Academy of Sciences President’s International Fellowship Initiative for Visiting Scientists (PIFI, 2026PVA0089), and the National Natural Science Foundation of China (NSFC) under grant No. 12473017.

References

Appendix A: Lagrangian radii of models U0, U20, and U40

The figure represents the time evolution of the Lagrangian radii containing from 0.1% to 99% of the actual total mass for models U0, U20, and U40.

Fig. A.1. Lagrangian radii for all cluster components (top row), single components (middle row), and binary components (bottom row) for models U0 (left column), U20 (middle), and U40 (right). In the bottom-left panel, we show only the 100% Lagrangian radius, due to the presence of extremely few dynamical binaries. As expected, there are no major differences in the evolution of the components in the mass segregation and expansion of the UFD.

Appendix B: HR diagrams of model U30

Figure B.1 presents HR diagrams of model U30 at different times, as labeled.

Fig. B.1. HR diagrams for 0, 25, 100, 500, 1000, 2000, 5000, 10000, and 13700 Myr (as labeled) for model U30.

Appendix C: σ_L comparison between observations and our models

In Fig. C.1 we compare the observations of UFD present in Supplement Table 1 of Simon (2019), where M_V > −7.7 as a condition for identification as UFD is used. We include our six models, after properly deriving the integrated V magnitude from our bolometric magnitudes. Our data are taken at the end of the simulation (Hubble time), except for model U50 (stopped at around 3 000 Myr). Our UFD models are in the high luminosity tail of the plot.

Fig. C.1. Velocity dispersion of UFDs as defined by Simon (2019, red dots) and our models (blue stars; U0 is the lowest and U50 the highest). Model U50 stopped at 3000 Myr.

All Tables

All Figures

	Fig. 1. Initial relaxation (blue) and dynamical friction (orange) times of our UFD Plummer model as functions of radius. The vertical dashed red and black lines denote the initial core and half-mass radius of the adopted Plummer model. For reference, the horizontal gray line marks the Hubble time.
In the text

	Fig. 2. Lagrangian radii evolution of the UFD model U30 over one Hubble time, using the time-updated value of the total mass. Top: For both single and binary stars. Middle: For single stars. Bottom: For binary stars.
In the text

	Fig. 3. Average mass within the different Lagrangian radii of the U30 model over one Hubble time.
In the text

	Fig. 4. Top: Fraction of mass present in the UFD, compared to the total mass of each component over the simulation time. Bottom: Evolution of the members of dynamical binaries the survive in the UFD (blue, right-side ordinate) and the fraction of primordial binaries that survive in the UFD (orange, left-side ordinate).
In the text

	Fig. 5. Evolution of the average semimajor axis for the binaries in model U30. The orange curve corresponds to the sample of primordial binaries; the blue curve also takes the dynamical binaries into account.
In the text

	Fig. 6. Number of stars in different evolutionary stages over time (see Table 2).
In the text

	Fig. 7. Total bolometric luminosity of stars in different evolutionary stages over time (model U30; see Table 2).
In the text

	Fig. 8. Time evolution of the total bolometric luminosity (blue) of model U30 and of the luminosity of single stars (orange) and binaries (green).
In the text

	Fig. 9. Model U30. Top: σg (blue), σg, cm (orange), and σL (green) evolution over time. Bottom: σL for models with initial binary fractions in the 0 − 50% range.
In the text

	Fig. 10. Velocity dispersions (σL) for all models, from U0 (top left) to U50 (bottom right).
In the text

Fig. A.1. Lagrangian radii for all cluster components (top row), single components (middle row), and binary components (bottom row) for models U0 (left column), U20 (middle), and U40 (right). In the bottom-left panel, we show only the 100% Lagrangian radius, due to the presence of extremely few dynamical binaries. As expected, there are no major differences in the evolution of the components in the mass segregation and expansion of the UFD.

In the text

	Fig. B.1. HR diagrams for 0, 25, 100, 500, 1000, 2000, 5000, 10000, and 13700 Myr (as labeled) for model U30.
In the text

	Fig. C.1. Velocity dispersion of UFDs as defined by Simon (2019, red dots) and our models (blue stars; U0 is the lowest and U50 the highest). Model U50 stopped at 3000 Myr.
In the text