Impact of Gaia-Sausage-Enceladus Debris on Milky Way Escape-velocity Estimates in TNG50

Introduction

In the $\Lambda$CDM paradigm, galaxies grow through hierarchical accretion and mergers of lower-mass systems ((Press & Schechter 1974; White & Rees 1978). The Milky Way (MW) is no exception to this scenario, as its stellar halo preserves tidal debris from multiple accretion events, including disrupted dwarf galaxies (e.g., [8, 5] (Helmi et al. 1999)). Therefore, the kinematics of these remnants provides a powerful probe of MW dynamics (for recent reviews, see [12] (Helmi 2020; Hunt & Vasiliev 2025)). This context motivates a reexamination of methods used to infer MW properties, specifically the escape velocity.

Corresponding author: Haining Li

[email protected]

The escape velocity, defined as the maximum velocity for a star bound to the gravitational field at a given location, is a measure of the potential depth of the MW. The seminal work of Leonard & Tremaine (1990), hereafter LT90, proposed a standard method to measure the escape velocity by modeling the velocity distribution of high-velocity stars as a truncated power law. Combined with a mass model of the MW, this formalism has been extensively applied to estimate the escape velocity and thereby constrain the dark matter (DM) mass and the overall mass distribution (e.g., [13] (Smith et al. 2007; Piffl et al. 2014; Williams et al. 2017; Monari et al. 2018; Koppelman & Helmi 2021; Necib & Lin 2022a; Prudil et al. 2022; Roche et al. 2024; Wu et al. 2025)). However, the theoretical assumptions underlying the LT90 method have been increasingly challenged. One of the most critical issues is that the LT90 method relies on a smooth, well-mixed stellar distribution, which may be violated by the complex, unrelaxed Galactic stellar halo ([23]). Thanks to the massive astrometric data from the Gaia mission ([22] (Lindegren et al. 2016)), combined with other large spectroscopic surveys such as RAVE ((Steinmetz et al. 2006)), the Sloan Digital Sky Survey (SDSS, York et al. 2000), LAMOST ([10] (Zhao et al. 2012; Yan et al. 2022)), APOGEE (Majewski et al. 2017), GALAH ([11]), numerous dynamical substructures have been identified (e.g., [4] (Haywood et al. 2018; Helmi et al. 2018; Koppelman et al. 2019; Matsuno et al. 2019; Myeong et al. 2019; Naidu et al. 2020; Yuan et al. 2020; Malhan 2022; Malhan & Rix 2024)). These substructures strongly overlap in the local phase space (e.g., [14] (Naidu et al. 2020; Horta et al. 2023; Zhang et al. 2024; Xie et al. 2025)), introducing kinematic complexity. Moreover, systems with substructures exhibit high-velocity tails that differ from those of a smooth, fully relaxed halo (e.g., [23] (Koppelman & Helmi 2021; Necib & Lin 2022b; Roche et al. 2024)).

Among the various substructures, Gaia-Sausage-Enceladus (GSE; [4] (Haywood et al. 2018; Helmi 2020)) stands out as a dominant component of ex situ stars in the MW (e.g., Nelson et al. 2019a; Naidu et al. 2020, 2021). Such debris, with highly radial orbits and high anisotropy, can populate the high-velocity tail ([13] (Necib et al. 2019)). Since escape-velocity estimates are sensitive to the shape of the high-speed tail, it is essential to validate the LT90 method in systems with GSE structures.

Motivated by recent findings of GSE-like debris in TNG50 MW analogs ([20]), we investigate how GSE-like structures bias escape-velocity estimates using MW analogs in Pillepich et al. (2024). In Section 2, to construct our mock solar neighborhood samples, we isolate GSE-like stars using purely kinematic criteria (e.g., [9] (Naidu et al. 2020; Lane et al. 2022)) and incorporate LAMOST magnitude and footprint selection ([10] (Zhao et al. 2012; Yan et al. 2022)) to represent realistic survey data. We introduce an adaptive selection of the high-velocity tail and the modeling process of $v_{\mathrm{esc}}$ in Section 3. By comparing the modeled $v_{\mathrm{esc}}$ with the true $v_{\mathrm{esc}}$ of the simulation, we rigorously evaluate the extent to which removing these GSE-like structures can recover the true $v_{\mathrm{esc}}$ and total mass, as described in Section 4 and Section 5. Finally, we summarize our findings in Section 6.

Data

TNG50 Simulations and MW Analogs

For our study, we used the high-resolution TNG50 simulation (Pillepich et al. 2019; Nelson et al. 2019b) of the IllustrisTNG project (Marinacci et al. 2018; Naiman et al. 2018; Pillepich et al. 2018a; Springel et al. 2018; Nelson et al. 2019a). The simulation is based on the Planck Collaboration et al. (2016) cosmology, with a Hubble constant of $H_0 = 100h\ \mathrm{km\ s^{-1}\ Mpc^{-1}} = 67.74\ \mathrm{km\ s^{-1}\ Mpc^{-1}}$, and evolved using the moving mesh code AREPO (Springel 2010). The initial conditions are generated in the early Universe (redshift $\approx127$) and evolved to the present day (redshift $= 0$), following a comprehensive model of galaxy formation (Weinberger et al. 2017; Pillepich et al. 2018b). The TNG50 simulation operates in a periodic cubic volume with a side length of 51.7 cMpc, containing $2160^3$ DM particles and an equal initial number of gas cells. The mass resolutions are $4.5 \times10^5\ M_\odot$ and $8.5 \times10^4\ M_\odot$ for DM and baryonic particles (Pillepich et al. 2024), respectively, which are crucial for properly resolving low-mass merger remnants with stellar masses of order $10^7$–$10^9\ M_\odot$ and provide new opportunities to constrain the formation and assembly histories of MW-mass galaxies.

We made use of the catalogs proposed by Pillepich et al. (2024), who identify 198 MW and M31 analogs at redshift $= 0$ according to stellar mass ($M_\star$), disk morphology, and environment. Since we focus only on MW analogs, an additional stellar mass cut ($M_\star< 10^{10.9}\ M_\odot$) was applied, resulting in 159 MW analogs. For these analogs, we extracted essential data for all particles and cells, such as coordinates, velocities, masses, and stellar photometry (Nelson et al. 2019a). The coordinates and velocities were converted to the coordinate system of the main galaxy, with more details in Pillepich et al. (2024). We excluded particles belonging to satellite galaxies, which have been studied in [18] and [17], and focused only on the remaining stars in the main galaxies for the subsequent dynamical and stellar property calculations.

To separate stellar populations based on their origins, we utilized SUBLINK merger trees (Rodriguez-Gomez et al. 2015) to trace particles from redshift $= 0$ back to their formation time. Following the definitions described in Rodriguez-Gomez et al. (2016), stars formed within the main progenitor branch are classified as in situ, which can be effectively distinguished from ex situ stars that formed in satellite galaxies and were subsequently accreted into the host.

Properties of MW Analogs with GSE-like Features

The methodology for deriving dynamics from TNG50 simulations is more straightforward than in observations, as the mass distribution of each type of particle is directly accessible. In this study, we adopted two right-handed coordinate systems: a cylindrical system $(R, \phi, z)$ and a spherical system $(r, \theta, \phi)$, both centered on the main galaxy. The gravitational potential $\Phi(r)$ of the MW analogs was constructed with the AGAMA software package Vasiliev (2019) under the assumption of approximate spherical symmetry, using the positions and masses of all DM, gas and stellar particles. A detailed discussion on the impact of adopting different symmetry assumptions is provided in Appendix § 8. We then computed the orbital parameters, including orbital energy $(E)$, azimuthal action $(J_\phi= -Rv_\phi)$, circularity $(\eta= J_\phi/|J_{\rm circ}(E)|)$, orbital apocenter $(r_{\rm apo})$ and pericenter $(r_{\rm peri})$, and eccentricity $(e = \frac{r_{\rm apo}-r_{\rm peri}}{r_{\rm apo}+r_{\rm peri}})$. Here, $J_{\rm circ}(E)$ is the angular momentum of a circular orbit with the same energy $E$. We define negative $v_\phi$ as corresponding to disk rotation; thus, positive $J_\phi$ corresponds to prograde orbits.

While a variety of chemical and dynamical criteria have been employed in the literature to identify GSE stars (e.g., Feuillet et al. [19]; Naidu et al. 2020; Lane et al. 2022; Horta et al. 2023; Lane et al. 2023; Carrillo et al. [9]; Zhang et al. 2024), selections relying on energy or angular momentum are inherently sensitive to the varying potential depths of host galaxies. In this study, we select GSE-like stars by combining intrinsic ex situ classification with present-day kinematics, i.e., ex situ stars with high eccentricity $(e > 0.7)$ and low azimuthal velocity $(|v_\phi| < 100\ {\rm km\ s}^{-1})$, in order to select a more complete GSE-like sample (e.g., Lane et al. 2022; Carrillo et al. [9]). Furthermore, to ensure prominent GSE-like structures, we selected analogs in which the mass fraction of GSE-like stars to ex situ stars is $\gtrsim50\%$ (e.g., Naidu et al. 2020, 2021). The ex situ classification is an intrinsic simulation label and is not directly observable. In MW observations, accreted populations can instead be identified using chemical information, such as the [Mg/Mn] versus [Al/Fe] plane (e.g., Horta et al. 2023), combined with kinematics (e.g., Matsuno et al. 2021; Zhang et al. 2024; Xie et al. 2025). However, applying the same criteria to TNG50 is nontrivial because the predicted abundance distributions depend on the chemical-evolution model, the adopted selection boundaries, and observational uncertainties. We therefore use intrinsic ex situ labels to define a clean, idealized baseline for isolating the dynamical effect of accreted debris.

Furthermore, we impose no explicit constraint on the infall times of the GSE-like progenitors. Accretion times cannot be measured directly from observations, and systems selected by their present-day GSE-like kinematics span lookback times of $\sim3$–11 Gyr (see also Folsom et al. [20]). This range extends to more recent epochs than the commonly inferred 8–11 Gyr lookback-time interval for the MW’s GSE accretion event (e.g., Helmi et al. 2018; Naidu et al. 2021; Donlon et al. [15]). Nevertheless, our subsequent goodness-of-fit criterion naturally disfavors systems containing strongly unrelaxed debris from very recent mergers.

We note that our selected GSE-like stars may originate from multiple merger events, consistent with recent studies suggesting that GSE stars may comprise debris from more than one progenitor (e.g., Orkney et al. 2023; Folsom et al. [20]). Regardless of their accretion histories, these stars occupy phase-space regions similar to observed GSE populations, making them suitable proxies for GSE debris when evaluating their effects on the escape velocity. Our final sample contains 28 MW analogs with prominent GSE-like debris, 19 of which overlap with the sample of Folsom et al. [20]. Table 1 summarizes the selection criteria and the number of the resulting galaxies. In practice, not all MW analogs have continuous high-velocity tails that allow statistically reliable $v_{\rm esc}$ fits. We therefore restrict our analysis to the valid subset of 16 galaxies, as detailed in § 4.2 and Appendices § 9 and § 10.

Selection Criteria	Number of Galaxies
MW/M31 analogs from Pillepich et al. (2024)	198
$M_{\star} < 10^{10.9}\,M_{\odot}$	159
Mass fraction of the GSE-like stars to ex situ stars $\gtrsim50\%$	28
Reliable escape-velocity modeling	16

As noted in previous work, the TNG50 MW analogs typically exhibit more extended disks (scale lengths $R_{\rm d} = 4.4^{+3.2}_{-1.8}\ {\rm kpc}$; e.g., Sotillo-Ramos et al. 2022; Pillepich et al. 2024) than the MW, which has $R_{\rm d,MW} = 1.7$–2.9 kpc (e.g., Hammer et al. [24]; Jurić et al. 2008; Bovy & Rix [7]; Bland-Hawthorn & Gerhard [6]). Local kinematic and mass distributions are poorly reproduced at a fixed physical radius of 8.2 kpc (Bland-Hawthorn & Gerhard [6]). Consequently, we compared the circular velocity curves calculated by $v_{\rm circ}(r) = \sqrt{GM(<r)/r}$ with observations to determine the physical solar neighborhood.

The left panel of Figure 1 presents the circular velocity curves for 28 MW analogs with prominent GSE-like structures in our study. In most galaxies, the circular velocity curves become flat at the solar radius, $r_{\odot}/r_{200} \approx0.043$ in Eilers et al. [16].¹ Therefore, we adopted the regime corresponding to $|r - 8.2\ \mathrm{kpc}| < 1.0\ \mathrm{kpc}$ in the MW, i.e., $0.038 < r/r_{200} < 0.048$, as the physically defined solar neighborhood.² We did not impose an additional cut on $|z|$. As a robustness check, we repeated the analysis for the representative well-mixed system (Subfind ID: 505100) using the cuts of $|z| < 2\ \mathrm{kpc}$ and $|z| < 5\ \mathrm{kpc}$. The inferred escape velocities are consistent with the result obtained without a vertical cut, although the vertical cuts lead to broader constraints due to the reduced sample sizes.

Vector (PDF)

The local escape velocity at the solar volume is widely defined as

where $r_{\mathrm{cut}}$ is an arbitrary limiting radius for realistic, nonisolated galaxies like the MW. Unlike an isolated system, the gravitational potential in cosmological simulations does not vanish at infinity, implying that the escape velocity should be defined relative to a finite reference radius. Previous studies have generally taken $r_{\mathrm{cut}}$ as $2r_{200}$ (e.g., [13] (Necib & Lin 2022a; Roche et al. 2024; Wu et al. 2025)) or $3r_{340}$ (e.g., Piffl et al. 2014; Monari et al. 2018; Koppelman & Helmi 2021).

To assess the sensitivity of the local escape velocity to the reference radius, we calculated $v_{\mathrm{esc}}(r_{\odot})$ over a range of $r_{\mathrm{cut}}$ values and compared with the most recent measurements of $v_{\mathrm{esc}}(r_{\odot}) = 523.74^{+12.83}_{-13.47}\ \mathrm{km\ s^{-1}}$ (Wu et al. 2025), as illustrated in the right panel of Figure 1. The comparison indicates that the median simulated escape velocity is consistent with the observed value at $r_{\mathrm{cut}} \approx2.16 \times r_{200}$. This radius is consistent with the region where the simulated $v_{\mathrm{esc}}(r_{\odot})$ profile begins to flatten. Beyond this radius, the local escape velocity changes little with increasing $r_{\mathrm{cut}}$ ($\lesssim10\ \mathrm{km\ s^{-1}}$ within the range of $2$–$3r_{200}$), indicating that matter at larger radii contributes only weakly to the local gravitational potential. This convergence provides independent support for our adopted $r_{\mathrm{cut}}$, allowing it to reasonably serve as a physical boundary of the simulated galaxies. Different choices of $r_{\mathrm{cut}}$ would shift the absolute recovery fractions but would not affect the relative comparison between the “with GSE” and “no GSE” scenarios.

Samples for Local Escape-velocity Estimates

We then constructed samples for local escape velocity estimates. Since the LT90 method focuses on halo stars, we identified the stellar halo of the MW analogs using the Toomre diagram (see, e.g., Helmi et al. (2018); Koppelman et al. (2019); Dodd et al. [14]), defined as:

where $v_{\mathrm{LSR}}$ is the velocity of the local standard of rest (LSR). Since our sample of MW analogs spans a very narrow circular velocity range ($v_{\mathrm{circ}} \sim200$–$250\ \mathrm{km\ s^{-1}}$, see Table A1), adopting a fixed velocity threshold yields results nearly identical to a scaled cut while maintaining direct consistency with MW observational works. Figure 2 displays the distribution of local stars in kinematic space for a typical MW analog (Subfind ID: 505100), with $v_{\mathrm{circ}}(r_{\odot}) = 215.05\ \mathrm{km\ s^{-1}}$.

Index	Subfind ID	$M_\star$	$r_{200}$	$M_{200}$	$v_{\mathrm{circ}}(r_\odot)$	$v_{\mathrm{esc}}^{\mathrm{sph}}(r_\odot)$	$v_{\mathrm{esc}}^{\mathrm{axi}}(r_\odot)$	$f_{\mathrm{local,GSE}}$	Status/Label
		$[10^{10}\ M_\odot]$	$[\mathrm{kpc}]$	$[10^{12}\ M_\odot]$	$[\mathrm{km\ s}^{-1}]$	$[\mathrm{km\ s}^{-1}]$	$[\mathrm{km\ s}^{-1}]$
1	446665	3.24	253.11	1.74	212.96	537.56	540.16	0.61	Valid
2	454172	3.63	248.18	1.62	215.45	512.41	514.88	0.49	Valid
3	465255	6.17	285.78	2.51	258.43	625.70	628.79	0.72	Valid
4	492876	6.92	263.10	1.95	243.53	581.70	586.41	0.66	Valid
5	496186	4.90	262.50	1.95	249.81	591.46	594.05	0.66	Valid
6	498522	7.59	259.81	1.86	249.87	589.79	595.49	0.42	Valid
7	504559	3.98	257.82	1.82	239.57	567.93	570.55	0.62	Valid
8	505100	3.98	249.72	1.66	215.05	501.28	503.16	0.63	Valid
9	515296	5.01	245.63	1.58	241.35	558.94	561.17	0.35	Valid
10	528322	3.16	207.45	0.95	204.42	469.67	470.54	0.64	Valid
11	530330	4.68	234.77	1.38	222.20	524.12	528.24	0.70	Valid
12	538905	4.27	218.43	1.12	205.39	490.68	496.33	0.52	Valid
13	540082	4.07	207.68	0.95	217.06	484.24	486.77	0.64	Valid
14	540920	3.39	222.47	1.17	211.07	496.24	503.99	0.80	Valid
15	541218	3.31	223.16	1.17	213.44	493.88	496.57	0.80	Valid
16	552414	4.57	203.72	0.91	195.54	485.98	490.28	0.62	Valid
17	422754	6.46	305.66	3.02	257.03	642.79	645.41	0.63	A
18	440407	6.61	258.00	1.82	237.28	562.59	564.95	0.64	A
19	511920	6.31	245.69	1.58	253.67	581.63	584.24	0.74	A
20	512425	3.55	240.28	1.48	209.41	524.01	527.63	0.56	A
21	534628	4.68	223.36	1.20	223.92	516.35	521.10	0.55	B
22	537941	4.37	207.18	0.95	195.73	485.63	492.34	0.61	C
23	543114	4.68	216.65	1.10	247.63	532.83	536.26	0.38	B
24	543729	4.07	208.18	0.95	184.64	478.24	485.16	0.61	A
25	550149	3.55	204.46	0.91	200.73	466.43	470.35	0.64	A
26	552581	3.47	202.85	0.89	217.89	486.09	488.22	0.67	B
27	565251	2.69	188.43	0.71	178.28	438.87	442.52	0.29	A
28	571633	3.47	182.58	0.65	179.88	438.47	445.29	0.34	C

The left panel shows that the halo and disk populations are distinctly separated by the Toomre diagram in velocity space. The middle and right panels present the selected GSE-like stars in the $v_{\phi}$–$v_{r}$ and $E$–$J_{\phi}$ planes, which confirm that the GSE-like particles cluster around $J_{\phi} \sim0$ while spanning a wide range of orbital energies (see also Naidu et al. (2020); Dodd et al. [14]; Horta et al. (2023); Zhang et al. (2024); Xie et al. (2025)). In addition to the dominant halo and GSE-like structures, we can also identify a less prominent retrograde disk component with $\eta< -0.7$. We note that this component has a negligible impact on the validity of our subsequent escape velocity measurements for this galaxy.

In order to mimic observational effects, we also took LAMOST as an example and generated subsets with its selection strategy. We converted the simulated intrinsic photometry to apparent magnitudes in the SDSS $r$ band (Nelson et al. 2019a).³ Applying the standard selection criteria of the LAMOST survey, decl. $> -10^\circ$ and $9.0 < r_{\mathrm{mag}} < 17.5$ (Cui et al. [10]; Zhao et al. (2012); Yan et al. (2022)), we defined our “LAMOST sample”, which allows for a direct comparison with the complete “All sample”. Here, the apparent magnitude of each star particle is used as a practical proxy for the detectability of the stellar population it represents, since our intention is to assess the statistical impact caused by the LAMOST-like observational selection, rather than constructing a fully realistic synthetic survey via full sampling.

Methodology

In this section, we present our framework for modeling $v_{\mathrm{esc}}$ and investigating the impact of GSE-like debris. First, we introduce a statistical method to evaluate when this model is valid. Then, we isolate the impact of GSE by comparing fits from stellar subsets with and without GSE-like stars. We also consider observational selection effects to provide practical insights into observations.

Escape-velocity Modeling

As introduced in § 1, the standard LT90 approach assumes velocity distributions following

for $v_{\min} \leq v < v_{\mathrm{esc}}$, where $v_{\mathrm{esc}}$ is the escape velocity, $k$ is the exponent, and $v_{\min}$ is the low-velocity threshold. Since the fraction of unbound stars in the TNG50 simulated galaxies is negligible, we did not include the outlier component adopted in recent observational works (e.g., Williams et al. 2017; Necib & Lin 2022a; Roche et al. 2024; Wu et al. 2025). Therefore, the total likelihood function is given as:

We utilized a Bayesian framework to infer the free parameters $v_{\mathrm{esc}}$ and $k$. Our baseline setup assumed a uniform prior on $k$, $\mathcal{U}(0.1 < k < 10)$, and a widely used prior $\propto1/v_{\mathrm{esc}}$ over $v_{\min} < v_{\mathrm{esc}} < 1000\ \mathrm{km\ s^{-1}}$, similar to previous works (Koppelman & Helmi 2021; Necib & Lin 2022a; Prudil et al. 2022). To estimate the posterior distributions of $v_{\mathrm{esc}}$ and $k$, the Markov Chain Monte Carlo (MCMC) sampling was performed in EMCEE (Foreman-Mackey et al. [21]), using 50 walkers, 1000 steps as a burn in, and additional 2000 steps for further sampling.

However, one of the most critical steps in this Bayesian inference is to determine an appropriate $v_{\min}$. While the large volume of the simulation data allows for well-converged posterior distributions (Piffl et al. 2014), the validity of the fitted parameters depends on whether the adopted $v_{\min}$ correctly identifies the regime where the truncated power law holds. Unlike observational works that often choose a fixed $v_{\min}$ (e.g., Monari et al. 2018; Deason et al. [13]; Wu et al. 2025) for the solar neighborhood, such an assumption is not suitable for our study. The velocity distribution in simulated galaxies varies from galaxy to galaxy and may differ from that of the real MW. The reliability of the truncated power law for TNG50 MW analogs therefore warrants further verification.

To address this, we implemented an adaptive $v_{\min}$-selection strategy. We scanned candidate $v_{\min}$ values over a 100 km s$^{-1}$ interval with a step size of 2 km s$^{-1}$, setting the upper limit between 320 and 380 km s$^{-1}$ to retain a statistically large number of high-velocity GSE-like stars. For each trial $v_{\min}$, we modeled the velocity distribution of local halo stars with $v > v_{\min}$ to derive the parameters $v_{\mathrm{esc}}$ and $k$. We then evaluated the goodness of fit by calculating the Kolmogorov–Smirnov (KS) $p$-value, comparing the modeled distribution with that of the simulated data. As $v_{\min}$ increases and the velocity distribution enters the power-law decline regime,

the $p$-value is expected to undergo a sharp initial increase. Our fiducial $v_{\min}$ was then chosen as the velocity at which this initial rise ends, which can be identified as the first prominent peak with the $p$-value $> 0.05$, as detailed in § 4.1. This strategy retains the largest possible sample while ensuring that the high-velocity tail has already entered a stable power-law decline, balancing a valid fit with sufficient statistics. Since observed high-velocity tails consistently exhibit a clear power-law decline (e.g., Monari et al. 2018; Koppelman & Helmi 2021; Prudil et al. 2022; Necib & Lin 2022b; Roche et al. 2024; Wu et al. 2025), systems that produce valid fits are considered closer to the MW observations.

During the search process, the fitted parameters $v_{\mathrm{esc}}$ and $k$ typically show a plateau near the fiducial $v_{\min}$, with variations of $\Delta v_{\mathrm{esc}} < 5$ km s$^{-1}$ and/or $\Delta k < 0.1$ in most well-modeled galaxies. This behavior implies that our method can indeed identify a stable power-law decline of high-velocity tails. For systems with sparse high-velocity tails, however, posterior distributions can become strongly degenerate at high $v_{\min}$. In these cases, we obtained the unified inference using the posterior of $k$ at the lower $v_{\min}$ as the prior at the higher $v_{\min}$.

Comparative Analysis

To systematically evaluate the impact of GSE-like structures and observational limitations, we considered four distinct scenarios: retaining (“with GSE”) versus removing (“no GSE”) the GSE-like stars, applied across both the complete “All sample” and the observationally motivated “LAMOST sample”. Applying our pipeline to these scenarios, we obtained the posterior distribution of $v_{\mathrm{esc}}$ to directly compare with the true local $v_{\mathrm{esc}}$ values calculated in § 2.2.

Furthermore, we propagated these $v_{\mathrm{esc}}$ estimates into a mass modeling to assess the impact of GSE on mass inference. We modeled the total potential by combining a fixed baryonic component with a Navarro-Frenk-White (NFW; Navarro et al. 1996) DM halo. The baryonic potential of each simulated galaxy was constructed directly via multipole expansion in the AGAMA software package (Vasiliev 2019), using the positions and masses of all gas and stellar particles. MCMC sampling was then performed using EMCEE (Foreman-Mackey et al. [21]) to constrain the DM virial mass ($M_{\mathrm{DM,200}}$) and the concentration ($c_{\mathrm{DM,200}}$), using the full posterior distribution of $v_{\mathrm{esc}}$. During the MCMC sampling, we fixed the contribution of the baryonic component and varied only the parameters of the NFW DM halo. The total virial mass was then calculated by summing the baryonic mass and the DM virial mass for each MW analog. Given the significant deviation of the inner DM halo from an NFW profile in the TNG cosmological simulation (e.g., Lovell et al. 2018), we added the circular velocity at $0.1r_{200}$ instead of the local circular velocity as an additional constraint in the likelihood evaluation.

Results

Validation on a Single Galaxy

To validate our adaptive $v_{\min}$-selection strategy, we present the results for a typical MW analog with a relatively smooth velocity distribution (SUBFIND ID: 505100) in Figure 3. The left panels show the results of our adaptive $v_{\min}$ selection. For both the “All sample” and “LAMOST sample”, the $p$-value exhibits an initial rise followed by a drop as $v_{\min}$ increases. Since this trend remains highly consistent regardless of whether the sample contains GSE-like stars, we adopted the same fiducial $v_{\min}$ for the “with GSE” and “no GSE” scenarios.

Vector (PDF)

As described in § 3.1, the initial rise indicates that the velocity distribution has entered the region governed by a power-law decline. The subsequent drop indicates that when $v_{\min}$ is too high, the statistical reliability is compromised due to the significant decrease in the sample size. The extreme high-velocity tail is particularly sensitive to discrete, unrelaxed substructures. Even in this well-mixed system, the distribution fluctuates more severely at extremely high velocities, which reduces the goodness of fit to a smooth power-law model.

In the overall velocity distributions (middle panels of Figure 3), our fiducial $v_{\min}$ corresponds to the onset of a stable power-law decline, i.e., beyond this threshold, the exponent $k$ remains relatively constant. The final modeling results shown in the right panels reproduce the decreasing trend of the high-velocity tails well, indicating that our $v_{\min}$-selection strategy identifies an appropriate $v_{\min}$ consistent with the definition in LT90.

Results of Multiple Milky Way Analogs

We then expanded our adaptive $v_{\min}$-selection pipeline to the full suite of 28 MW analogs (see Appendices § 9 and § 10 for details). Among these systems, 16 galaxies with relatively continuous and smooth high-velocity tails yield reliable escape velocity fits. By contrast, our pipeline naturally rejects the remaining 12 systems with high-velocity tails that have obvious clumpy features and thereby deviate from the single power-law model, as shown in Appendix § 9. For such systems with clumpy velocity distributions, an alternative approach is to introduce multiple power-law components to account for distinct substructures (e.g., Necib & Lin 2022b). However, since our primary objective is to evaluate the reliability of the standard LT90 method, which assumes a smooth, well-mixed stellar halo, we do not introduce additional power-law components, allowing our pipeline to naturally filter out unrelaxed systems. Taken together, these results indicate that our pipeline provides a robust diagnostic for identifying systems with dynamically smooth stellar halos.

In Appendix § 10, even among the 16 galaxies with reliable escape-velocity modeling, the accuracy of the LT90 method depends sensitively on the smoothness of the high-velocity tails (see also Grand et al. [23]). In some galaxies (e.g., SUBFIND IDs: 498522, 538905, and 540082), the true $v_{\mathrm{esc}}$ lies in the $1\sigma$ confidence intervals of the fitted $v_{\mathrm{esc}}$ for both the “with GSE” and “no GSE” scenarios. This behavior suggests that the velocity distribution of GSE-like components is consistent with that of the underlying smooth halos in these systems, i.e., GSE-like structures may be well-mixed in the host halos.

However, in other systems whose velocity distributions exhibit multiple small bumps introduced by substructures (e.g., SUBFIND IDs: 492876, 504559, and 530330), GSE-like structures substantially reshape the high-velocity tails, producing more extended tails and hence lower $k$ values (see also [13] (Necib & Lin 2022a; Roche et al. 2024)). Removing the GSE-like stars leads to the escape-velocity estimates which are more consistent with the true values. A detailed discussion of the $v_{\mathrm{esc}}$-modeling bias caused by GSE-like structures is provided in § 5.

Comparing the modeling results for the “All sample” and the “LAMOST sample”, we note that observational selection effects may randomly change the smoothness of the high-velocity tails (e.g., Subfind IDs: 454172 and 530330). However, the systematic bias caused by GSE-like structures and the improvement from removing GSE-like stars still persist.

Discussion

To quantify the systematic bias induced by GSE-like structures, we construct joint posterior distributions for 16 final galaxies and compare the fitting results between the “with GSE” and the “no GSE” scenarios. Figure 4 shows the comparison of the recovery fraction, given by $v_{\mathrm{esc,fitted}}/v_{\mathrm{esc,true}}$, with the $1\sigma$ confidence intervals defined by the 16th and 84th percentiles. Among these galaxies, the proportion of GSE-like stars within high-velocity tails spans a broad range, from $\sim10\%$ to 50%. We note that the recovery fraction is insensitive to the proportion of GSE-like stars, confirming that the systematic bias is driven by differences in the intrinsic dynamical properties rather than by the relative fraction of the unrelaxed substructures.

As shown in the left panels of Figure 4, the marginalized distribution for the “with GSE” scenario of the “All sample” displays a broad, noticeably left-skewed profile, yielding a systematic underestimate with a recovery fraction of $97.5^{+2.1\%}_{-2.9\%}$ and a deviation of $1.2\sigma$ from the true values. This underestimate can be physically attributed to the dynamical properties of the substructures. Stars from a single merger event, such as the GSE-like merger, possess distinct binding energies and phase-space distributions, failing to populate high-velocity tails efficiently and continuously (e.g., [23] (Necib et al. 2019; Necib & Lin 2022b)). Such strongly localized debris that produces localized bumps or discontinuities is therefore generally rejected by our goodness-of-fit criterion. However, early merger debris may become phase mixed sufficiently to produce a smooth velocity distribution while still retaining dynamical properties distinct from those of the underlying stellar halo. In this case, the debris can contribute continuously over a broad velocity range without causing the single power-law model to be rejected. The fitting algorithm accommodates these components by adopting a flatter decline (lower $k$), which inherently biases the best-fit $v_{\mathrm{esc}}$ towards lower values. This explains the observed skewed distribution in the “with GSE” scenario even after systems with the most prominent bumpy substructures have been excluded.

Similar underestimates have also been reported in the Auriga simulations (see, e.g., Grand et al. [23]; Koppelman & Helmi (2021)). Our inferred offset is smaller ($\sim2\%$) compared to $\sim7\%$ in Grand et al. [23] or $\sim10\%$ in Koppelman & Helmi (2021). In addition to differences in simulation resolution and in the exact truncation location of high-velocity tails, this reduction likely reflects the influence of our pipeline, which selects only systems with robust $v_{\mathrm{esc}}$ fits.

Crucially, as discussed in § 4.2, removing GSE-like structures effectively mitigates this bias. For the “no GSE” scenario of the “All sample”, the right marginalized distribution in the left panel of Figure 4 shows a more symmetric shape. The recovery fraction of the true $v_{\mathrm{esc}}$ increases to $99.4^{+2.1}_{-2.5}\%$, with a reduced deviation of $0.34\sigma$. This confirms that GSE-like structures can bias inferred escape velocities under the well-mixed assumption of the LT90 method and serve as an important and likely leading source of systematic bias among the systems considered here. Indeed, our findings are consistent with the established physical picture that dominant GSE structures with highly radial orbits (e.g., Belokurov et al. [4]; Haywood et al. (2018); Helmi et al. (2018); Naidu et al. (2020, 2021)) can populate high-velocity tails (Deason et al. [13]; Necib et al. (2019)).

The corresponding recovery fraction for the “LAMOST sample” is presented in the right panel of Figure 4. Observational selection effects manifest as broadening and flattening the overall posterior distribution, thereby erasing some intrinsic statistical features. However, despite this flattening effect, the improvement in the removal of GSE-like structures remains significant. The recovery fraction increases from $96.1^{+3.7}_{-3.6}\%$ to $98.7^{+3.7}_{-3.9}\%$, with a deviation decreasing from $1.06\sigma$ to $0.41\sigma$. Therefore, our pipeline, which combines the adaptive $v_{\min}$-selection strategy with the removal of GSE-like stars, remains effective under this simplified observational selection. It effectively mitigates intrinsic physical biases and remains insensitive to observational selection effects, providing a reliable framework for accurately estimating $v_{\mathrm{esc}}$.

Furthermore, we demonstrate the posterior distributions of the mass recovery fraction, $M_{200,\mathrm{fitted}}/M_{200,\mathrm{true}}$, in Figure 5. To better demonstrate the posterior distribution of the mass, we also modeled it using an asymmetric Gaussian model and derived the best-fit median and $1\sigma$ credible intervals. Comparing the results from the “with GSE” and the “no GSE” scenarios, the mass estimates show substantial overlap within the $1\sigma$ confidence intervals, both consistent with the values of the true total mass. Nevertheless, the results of the “no GSE” scenarios can better recover the true total mass. The $\sim2\%$ underestimation of $v_{\mathrm{esc}}$ results in an $\sim6\%$ underestimation

Vector (PDF)

of the total mass, indicating the necessity to consider the impact of GSE more carefully when estimating the total mass of the MW by the escape-velocity method.

Conclusions

In this work, we used TNG50 MW analogs to quantify how GSE-like debris biases local escape-velocity measurements based on the LT90 method. From 159 MW analogs in Pillepich et al. (2024), we selected 28 systems with prominent GSE-like structures, and 16 of them yielded statistically valid fits with continuous high-velocity tails. Our main conclusions are as follows.

1. We develop an adaptive $v_{\min}$-selection strategy that combines statistical reliability with sample size. This approach robustly identifies the velocity range in which the truncated power-law model is valid, and provides a reliable diagnostic for dynamically smooth stellar halos. This step is essential for identifying the subset of halos for which LT90-based inference is meaningful.

2. The accuracy of the inferred escape velocity is intrinsically related to the smoothness of the stellar halo’s velocity distribution. Complex assembly and evolution histories of MW analogs result in unrelaxed substructures and introduce clumpy features in the stellar velocity distribution, thereby invalidating the single power-law model.

3. Although multiple substructures may have an impact on high-velocity tails, GSE, as the most kinematically extreme and dominant substructure, is an important source of systematic underestimation in escape-velocity measurements for MW-like systems. Specifically, removing GSE-like stars alone significantly increases the recovery fraction of the true escape velocity ($v_{\mathrm{esc,fitted}}/v_{\mathrm{esc,true}}$) from $97.5^{+2.1\%}_{-2.9\%}$ to $99.4^{+2.1\%}_{-2.5\%}$ for the complete “All sample”.

4. The same trend persists after applying observational selection effects, such as the LAMOST-like magnitude and footprint selection. Despite broader posteriors from smaller sample sizes, removing GSE-like stars still improves the recovery fraction from $96.1^{+3.7\%}_{-3.6\%}$ to $98.7^{+3.7\%}_{-3.9\%}$ for the “LAMOST sample”. Since the results of the “All sample” have higher recovery fractions and nar-

rower $1\sigma$ confidence intervals, it is necessary to obtain more accurate escape velocity measurements using all-sky coverage data.

5. We also investigated the impact of GSE on the mass estimates. In this study, an $\sim2\%$ underestimate of the escape velocity leads to an $\sim6\%$ underestimate of the total mass for MW-like systems. Therefore, accounting for GSE substructures is essential for an accurate recovery of the Galactic potential and total mass.

Acknowledgements

We thank the anonymous referee for very efficient and helpful suggestions to improve the draft. This work is supported by the National Natural Science Foundation of People’s Republic of China (NSFC) grant No. 12588202, 12633005, and 12273053, the National Key $R\&D$ Program of People’s Republic of China Nos. 2024YFA1611902, 2024YFA1611903 and 2023YFE0107800, and the CAS Project for Young Scientists in Basic Research (Nos. YSBR-092 and YSBR-062), the Strategic Priority Research Program of Chinese Academy of Sciences grant No. XDB1160100, the science research grants from the People’s Republic of China Manned Space Project with Nos. CMS-CSST-2021-B03 and CMS-CSST-2025-A11, and National Astronomical Observatories Chinese Academy of Sciences Nos. E5ZB0301, E4ZB0301, E6510301 and E6510302. The TNG50 simulation was run with compute time granted by the Gauss Center for Supercomputing (GCS) under Large-Scale Projects GCS-DWAR on the GCS share of the supercomputer Hazel Hen at the High Performance Computing Center Stuttgart (HLRS).

Software: AGAMA ((Vasiliev 2019)), astropy ([1, 2, 3]), emcee ([21]), matplotlib ((Hunter 2007)), numpy ((Harris et al. 2020)), scipy ((Virtanen et al. 2020)), TNG access scripts ((Nelson 2021)).

Impact of Different Potential Approximations

Table A1 summarizes the fundamental properties of the 28 MW analogs with prominent GSE-like structures. To check on the impact of different potential approximations, we construct the gravitational potential of the MW analogs under the axisymmetric approximation. The procedure is implemented using multipole expansion for the positions and masses of all DM, gas and stellar particles, with the AGAMA software package ((Vasiliev 2019)). Then we calculate the local escape velocity using the modification of Equation (1), given by:

where $\Phi(r_{\mathrm{cut}})$ represents the average potential of the spherical shell at $r_{\mathrm{cut}}$. The local position is set to $(R,z)=(r_{\odot},0)$, and the results are also listed in Table A1.

We note that the typical error introduced by our spherical potential approximation is $<1\%$, and the spherical approximation tends to slightly underestimate the local escape velocity. Furthermore, when evaluating the recovery of the true escape velocity under the axisymmetric approximation, the recovery fraction of the escape velocity in the all-sky scenario increases from $96.8\%$ to $98.7\%$ after removing the GSE-like structures. Similarly, when considering LAMOST-like selection effects, the recovery fraction improves from $95.4\%$ to $98.0\%$. Since the improvement of $\sim2\%-3\%$ is consistent with that under the spherical approximation, we can conclude that different potential approximations have a negligible impact on our final results.

MW analogs with Invalid Fits

Table A1 also summarizes the fitting status of the 28 MW analogs. Among these, 16 galaxies exhibit continuous and relatively smooth high-velocity tails, yielding reliable $v_{\mathrm{esc}}$ fits under the LT90 formalism, i.e., the $p$-value $> 0.05$ at the fiducial $v_{\min}$. The remaining 12 systems fail to obtain valid fitting results. Figure B1 presents representative examples of these invalid cases, which can be attributed primarily to three distinct kinematic origins.

Vector (PDF)

First, the dominant reason is unrelaxed kinematic substructures, as shown in the left panel of Figure B1.

Note— Columns from left to right are: (1) index, (2) Subfind ID from the TNG50 simulation at redshift = 0, (3) stellar mass, (4) virial radius, (5) virial mass, (6) local circular velocity, (7) and (8) local escape velocity for spherical and axisymmetric approximations, respectively, (9) mass fraction of GSE-like stars to ex situ halo stars within the solar volume, and (10) fit status or the main reason for invalid fits using the LT90 method. Systems with reliable fits are marked as “Valid”. For systems with invalid fits, the labels are defined as follows: (A) clumpy high-velocity tails caused by recent massive mergers; (B) prominent contamination from retrograde disk-like structures; and (C) highly disk-dominated systems (disk fraction > 70%) with insufficient high-velocity halo stars.

Some galaxies have experienced relatively recent massive mergers ($\lesssim4$ Gyr);⁴ thus, the accreted components are not fully mixed in the stellar halo. Unrelaxed accreted components may produce clumpy, discontinuous high-velocity tails that violate the smooth-halo assumption of the LT90 method. Notably, several systems with recent massive mergers, e.g., Subfind ID 530330, still exhibit relatively smooth halo velocity distributions that are sufficient to yield valid fits.

Second, the middle panel of Figure B1 exhibits the velocity distribution of a MW analog with a prominent retrograde disk-like structure, which is kinematically defined as circularity $< -0.7$. Since disk and halo populations possess fundamentally different velocity distributions, prominent retrograde disk-like structures significantly distort the morphology of the high-velocity tail. Consequently, future precise measurements of the MW’s escape velocity may require careful consideration of these retrograde structures.

Finally, a few MW analogs are highly disk-dominated, in which over 70% of the stars belong to the disk population. These systems lack a prominent stellar halo; thus, we are unable to extract a statistically significant sample of local high-velocity halo stars. Furthermore, as shown in the right panel of Figure B1, the fastest halo stars do not reach the $\sim0.9 \times v_{\mathrm{esc}}(r_{\odot})$ threshold, which also limits the escape velocity modeling.

Escape Velocity Fitting Results

Figure C1 and Figure C2 show $v_{\mathrm{esc}}$ fitting results for the 16 final MW analogs with reliable fits, corresponding to the “All sample” and “LAMOST sample”, respectively. As we can see, these galaxies exhibit relatively smooth and continuous high-velocity tails and can be well described by a single power law.

Vector (PDF)

Vector (PDF)