© The Authors 2026

1 Introduction

The 3D shape of the dark-matter (DM) halo serves as a fundamental test of the Λ cold dark matter (ΛCDM) model (Peter et al. 2013) and retains key information about a galaxy’s filamentary accretion history from the cosmic web and interactions with massive satellites (Arora et al. 2025). Functioning as an intermediary between the Galaxy and its satellites, the shape of the DM halo also plays a key role in deciphering the long-standing conundrum of a vertically aligned ‘plane of satellites’ with the Galactic disc, in the case of the Milky Way (Shao et al. 2016).

Determining the exact shape of the DM halo from observations, even for the Milky Way, is a challenge. One major method of probing the DM distribution is to model the motions of tracers, which reflect the underlying gravitational potential. Stellar streams, viewed as simple orbital pathways that trace the halo, played a significant role in constraining the shape of the Milky Way DM halo (Bovy et al. 2016; Bowden et al. 2015; Koposov et al. 2010; Palau & Miralda-Escudé 2023; Law & Majewski 2010). However, the most commonly used dynamically cold streams, for example Pal 5 and GD-1, trace only the inner halo at r ≲ 20 kpc. And a single stream is easily disturbed (Nibauer et al. 2024) and may lead to biassed results (Vasiliev et al. 2021).

Halo stars - with the full 6D phase-space information observed - span the entirety of space and, statistically, should exhibit strong power in constraining the underlying 3D DM distribution. However, these stars exist with complex orbital dynamics, necessitating accurate dynamical modelling. At present, the Milky Way halo is analysed using two primary types of dynamical models: the Jeans model and the distribution function (DF)-based dynamical model. These models yield divergent results regarding the shape of the DM halo, ranging from oblate with q ~ 0.4-0.7 (Loebman et al. 2014; Hattori et al. 2021 ; Li & Binney 2022), to nearly spherical with q ~ 0.9-1 (Wegg et al. 2019; Zhang et al. 2025), and to prolate with q ~ 1.3 (Posti & Helmi 2019). All these models are restricted by axisymmetric assumptions and incorporate ad hoc assumptions about the intrinsic dynamical structures. Despite extensive research over the years, the true shape of the Milky Way DM halo remains debated (Hunt & Vasiliev 2025).

Due to the large uncertainty surrounding the DM halo shape of the Milky Way from observations, its link to the spatial arrangement of satellite galaxies has not been seriously considered. The satellite galaxies in the Milky Way are not distributed isotropically but rather form a thin plane vertical to the Galactic disc (Lynden-Bell 1976; Kroupa et al. 2005; Pawlowski 2018; Sawala et al. 2023). The thickness of this plane of satellites has been extensively studied; it was previously thought to be rare but has now become more frequently observed through cosmological simulations (Sawala et al. 2023). The cause of its perpendicular alignment to the Galactic disc remains uncertain (Vasiliev et al. 2021). Typically, both the stellar disc and satellite system align parallel to the DM halo, sharing a similar net angular momentum direction, leading to a co-alignment of the stellar disc and satellite plane (Shao et al. 2016). However, in systems resembling the Milky Way, in which the satellite system lies perpendicular to the stellar disc, the DM halo might exhibit a twist, being co-aligned with the stellar disc in the inner regions but aligned vertically in the outer regions (Shao et al. 2021). Precisely measuring the DM shape as a function of radius could provide essential insights into the origin of the vertical satellite plane and elucidate its relationship with the large-scale matter distribution. Additionally, this suggests that the DM halo of the Milky Way may be triaxial, with a shape that changes from the inner to outer regions, necessitating dynamic models with more flexible DM halos.

Recently, a triaxial Schwarzwald model that numerically represents the DF of tracers by superposition of stellar orbits has been applied to the Milky Way halo, indicating a triaxial, nearly prolate halo with a large tilt angle (Dillamore & Sanders 2025). Nevertheless, there exists a strong degeneracy between the axis ratios of the DM halo and the tilt angle. This model was mainly constrained by the shape of a tilted stellar halo (Han et al. 2022a) associated with the Gaia-Sausage-Enceladus (GSE) (Helmi et al. 2018; Belokurov et al. 2018).

In a previous paper (Zhu et al. 2025), we presented and validated an empirical triaxial orbit-superposition model tailored for the Milky Way halo observed by the Large Sky Area Multi-Object Fibre Spectroscopic Telescope (LAMOST) and the Global Astrometric Interferometer for Astrophysics (Gaia). This approach accommodates a highly flexible triaxial DM halo, where the axis ratios are allowed to vary with radius, and it uses the full 6D phase-space information from the observations. Compared to traditional made-to-measure and Schwarzschild methods, this method does not theoretically sample the orbit libraries or determine the orbit weights by fitting the model to data; instead, the orbit library and orbit weights are entirely determined from the observations. The results are thus highly data-driven.

In this paper, we apply the empirical triaxial orbitsuperposition model to the Milky Way halo, utilising 6D phasespace data from LAMOST and Gaia, and robustly determine the 3D DM distribution of the Milky Way, extending to approximately 50 kpc. We then compare the findings with Milky Way analogues from the Auriga and TNG50 cosmological simulations and establish a coherent scenario linking the DM halo shape and satellite arrangement. The paper is organised in the following way: we describe the details of observational data in Section 2, introduce model construction in Section 3, and show the results in Section 4. We compare the 3D shape of the DM halo and the arrangement of the Milky Way’s satellite system with cosmological simulations in Section 5. We discuss in Section 6 and conclude in Section 7.

2 Observational data

2.1 Data preparation

We obtained the 3D positions and 3D velocities of halo stars by combining observations from LAMOST and Gaia. LAMOST has gathered spectroscopic data for approximately 10 million stars within the Milky Way. We used k giants from LAMOST DR9; the distance is accurately measured with a relative uncertainty of ~15% (Zhang et al. 2023). By cross-matching with Gaia DR3, we obtained a sample of 619 284 k-giants with 3D velocity measurements; the typical errors in the radial velocity, v_r, tangential velocity, v_φ, and v_θ are 10, 30, and 30 km/s, respectively. We excluded disc stars chemodynamically to avoid biassing the spatial distribution of the halo. We first removed all of the stars in the cold disc with λ_z > 0.8, then removed the group of stars with λ_z > 0.3 and [Fe/H] > -1.2 that are mostly within 20 kpc, not in dynamical equilibrium, and likely to be induced from the disc (see Fig. 1). The remaining stars were considered halo stars. We carefully excluded all substructures, including streams and overdensities, using the friend-of-friend method (Sun et al. 2025; Wang et al. 2022; Yang et al. 2019). We further removed stars with velocity errors of >100 km/s and performed 3σ clipping to remove the outliers. We preserved a refined halo sample comprising 14497 K-giants, with 11 415 in the northern hemisphere and 3082 in the southern hemisphere. The sample spans a vast area in the northern hemisphere, up to r ≲ 50 kpc.

Fig. 1 Density of stars in the parameter phase of circularity, λz, versus metallicity, [Fe/H]. We first removed the disc stars with a circularity of λz > 0.8, then also removed the group of stars with λz > 0.3 and [Fe/H] > -1.2. The rest were kept as halo stars in our sample.

2.2 Correction of the selection function

We corrected for the selection effects of the LAMOST survey following previous studies (Liu et al. 2017), in which the photometric data were assumed to be complete within their limiting magnitude. In the colour-magnitude plane (c, m) along each line of sight (l, b), the photometric density, ν_ph, is equal to the spectroscopic density, ν_sp, times the inverse of the selection function, S⁻¹, $${\nu }_{\mathrm{p}\mathrm{h}}(d|c,m,l,b)={\nu }_{\mathrm{s}\mathrm{p}}(d|c,m,l,b)S^{-1}(c,m,l,b).$$(1)

The spectroscopic density, ν_sp(d|c, m, l, b), was calculated using the kernel density estimation (KDE) method: $${\nu }_{\mathrm{s}\mathrm{p}}(d|c,m,l,b)=\frac{1}{\mathrm{\Omega }{d}^2}\sum _i^{{\nu }_{\mathrm{s}\mathrm{p}}(c,m,l,b)}{p}_i(d),$$(2)

where p_i(d) is the probability density function of distance for the i_th star, and Ωd² is the volume element between d and d + ∆d.

The selection function, S⁻¹(c, m, l, b), was evaluated from $$S^{-1}(c,m,l,b)=\frac{n_{\mathrm{s}\mathrm{p}}(c,m,l,b)}{n_{\mathrm{p}\mathrm{h}}(c,m,l,b)},$$(3)

where n_sp(c, m, l, b) and n_ph(c, m, l, b) are the numbers of stars with spectroscopic and photometric data, respectively, within each bin of (c, m, l, b).

The stellar density along a given line of sight is thus $${\nu }_{\mathrm{p}\mathrm{h}}(d|l,b)=\iint {\nu }_{\mathrm{s}\mathrm{p}}(d|c,m,l,b)S^{-1}(c,m,l,b)\mathrm{d}c\mathrm{d}m.$$(4)

The selection-corrected average density distribution of the stellar halo in the (R_gc, z_gc) plane (ν_ph(R_gc, z_gc)) could then be constructed (Xu et al. 2018; Yang et al. 2022).

We binned the stars observed spectroscopically in the (R_gc, z_gc) plane and obtained the number of stars in each bin n_sp(R_gc, z_gc). For each of the stars in the bin, i, we calculated their photometric weight: $${\omega }_i=\frac{\overline{{\nu }_{\mathrm{p}\mathrm{h}}}(R_i,z_i)\times V}{n_{\mathrm{s}\mathrm{p}}(R_i,z_i)},$$(5)

where V = 2πR_idRdz is the volume of the bin. We considered all the stars in a 1.0 × 1.0 kpc bin to have the same weight.

We present our sample of k-giants in the northern hemisphere in Fig. 2. The observations close to the disc plane and at r ≳ 30 kpc are less complete; thus, the stars have higher weights. With weights correcting for their bias from observational selection, the sample in the northern hemisphere can be considered representative of the halo stars in 6D position-velocity distributions.

Fig. 2 K giants in our final smooth halo sample in the northern hemisphere, coloured by the weights of particles obtained by correction of the selection function. The sample can be taken as complete in the Rgc-zgc space after correction of the selection function.

2.3 Correction of the LMC bias

The influence of LMC on the stellar halo at r < 50 kpc is small. We binned the data in the R_gc-z_gc plane with a bin size of 10 × 10 kpc and evaluated their mean velocities, v_x, v_y, and v_z, within each bin. We observe a minor net velocity in v_z and v_x, around ~10 - 20 km/s, whereas a system in perfect dynamical equilibrium should have these velocities at 0. These offsets are presumably caused by the LMC’s influence (Erkal et al. 2020). We adjusted for this by deducting the mean v_z and v_x velocities from each star’s velocity, ensuring zero mean velocities in every bin. Although the LMC’s effect is limited in the areas under examination, correcting for it does not significantly alter our final results in constraining the 3D distribution of the DM halo. However, the agreement between the velocity distributions of the model and the data improves slightly when this correction is applied, as was expected.

3 Model construction

Our sample serves as a good approximation of being stationary, i.e. the stars’ DF f (x, u) will statistically not change when orbiting them in the correct gravitational potential. Driven by data insights, we presented a novel dynamical model (Zhu et al. 2025) aimed at constraining the gravitational potential by directly comparing the observed data distribution, f_data(x, u), with the distribution derived from a model superposed by the orbits of these stars integrated within the potential f_model(x, u|potential). This approach, named the empirical triaxial orbit superposition method - based on minimal assumptions of dynamical models that the system is stationary - overcomes key limitations of previous models. The method has been validated and proves effective at uncovering the 3D DM distribution in galaxies similar to the Milky Way. We have applied it to mock data created from Auriga simulations, selecting three galaxies with different 3D shapes of the underlying DM halo. In all three cases, our model accurately recovers the 3D DM distribution, including the axis ratios and the radial density profile. For one case with a twisted DM halo from the inner to outer regions, our model succeeds in capturing the radial variation in the DM halo shape.

3.1 Empirical triaxial orbit superposition method

Here we briefly introduce the main steps to create the empirical orbit superposition model: (1) we built a model for the gravitational potential with a few free parameters; (2) given one set of parameters, we used the 6D phase space information of stars in the northern hemisphere as initial conditions, and calculated their stellar orbits in the gravitational potential; (3) we superposed the stellar orbits with their weights given by the selection function correction and constructed stellar density distributions and velocity distributions from the orbit superposition model; (4) we compared the data and model to calculate the likelihood or χ² of each model; (5) we explored the parameter grid of the gravitational potential and found the best-fitting models with the maximum likelihood or least χ².

We constructed a model of the gravitational potential by including a barred bulge, a thin disc, a thick disc, and a triaxial DM halo with the following density profiles: $${\rho }_{\mathrm{b}\mathrm{u}\mathrm{l}\mathrm{g}\mathrm{e}}=\frac{105M}{32\pi pq{a}^3}{[1-{\left(\frac{\tilde{r}}{a}\right)}^2]}^2;$$(6) $${\rho }_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{c}}=\frac{{\mathrm{\Sigma }}_0}{2h}\times \exp (-{\left[\frac{R}{R_d}\right]}^{1/n}-\frac{R_{\mathrm{c}\mathrm{u}\mathrm{t}}}{R})\times \exp (-\frac{|z|}{h}).$$(7)

We fixed the density distribution of the bulge and discs obtained from previous studies (Bland-Hawthorn & Gerhard 2016): with M_*,bulge = 1.6 × 10¹⁰M_⊙, a scale radius of a = 3.5 kpc, axis ratios of p = 0.44, q = 0.31 for the bulge, M_*,thin = 3.16 × 10¹⁰ M_⊙, a scale radius of R_d = 2.6 kpc, a scale height of h = 0.3 kpc, an inner cut-off radius of R_cut = 7 kpc (Lian et al. 2024), a Sersic index of n = 1 for the thin disc and M_*,thick = 6 × 10⁹ M_⊙, a scale radius of R_d = 2.0 kpc, a scale height of h = 0.9 kpc, and a Sersic index of n = 1 for the thick disc¹.

We adopted a flexible triaxial generalised NFW model for the DM halo, allowing free orientations of the DM halo and a variable 3D shape as a function of radius, generally following (Vasiliev et al. 2021): $${\rho }_{\mathrm{h}\mathrm{a}\mathrm{l}\mathrm{o}}={\rho }_0{\left(\frac{\tilde{r}}{r_s}\right)}^{-\gamma }{[1+{\left(\frac{\tilde{r}}{r_s}\right)}^{\alpha }]}^{\frac{\gamma -\beta }{\alpha }}\times \exp -{\left(\frac{\tilde{r}}{r_{\mathrm{c}\mathrm{u}\mathrm{t}}}\right)}^{\xi },$$(8)

in which $\tilde{r}=\sqrt{X^2+(Y/p_{\mathrm{D}\mathrm{M}}{)}^2+(Z/q_{\mathrm{D}\mathrm{M}}{)}^2}$, p_DM, and q_DM are the ratios of the axes. The co-ordinates X ≡ {X, Y, Z} are related to the usual Galactocentric Cartesian co-ordinates x_gc {x_gc,y_gc, z_gc} by a rotation matrix, X = Rx_gc, parameterised by three Euler angles, α_q, β_q, and γ_q. The Z-axis is tilted by the angle β_q relative to the z_gc-axis. Meanwhile, the angles α_q (γ_q) define the rotation between the x_gc(X)-axis and the intersection line of the x_gcy_gc and XY planes.

The tilt angle β_q highly degenerates with the vertical axis ratio, q_DM. We fixed β_q = 0° and allowed q_DM > 1 in our model. The Z-axis is thus always defined as perpendicular to the stellar disc. Because the direct observations are not complete and we have averaged the information in the azimuthal direction, our model cannot constrain the azimuthal orientation of the DM halo. We fixed α_q = 0° in our fiducial model. In this model with fixed β_q and α_q, the axis ratio q_DM and p_DM will reflect the DM halo orientation.

For the radial profile, we fixed α = 1, β = 3, and chose the outer cut-off radius r_cut = 500 kpc and the cut-off strength ξ = 5. We were left with five free parameters in the DM halo: ρ₀, r_s, γ, ρ_DM, and q_DM.

When needed, we allowed p_DM and q_DM to vary as a function of radius, with $$p_{\mathrm{D}\mathrm{M}}(r)=(p_{\mathrm{i}\mathrm{n}}+p_{\mathrm{o}\mathrm{u}\mathrm{t}}{\left(\frac{r-10}{r_q}\right)}^2)/(1+{\left(\frac{r-10}{r_q}\right)}^2)$$(9) $$p_{\mathrm{D}\mathrm{M}}(r<10)=p_{\mathrm{i}\mathrm{n}}$$(10) $$q_{\mathrm{D}\mathrm{M}}(r)=(q_{\mathrm{i}\mathrm{n}}+q_{\mathrm{o}\mathrm{u}\mathrm{t}}{\left(\frac{r-10}{r_q}\right)}^2)/(1+{\left(\frac{r-10}{r_q}\right)}^2)$$(11) $$q_{\mathrm{D}\mathrm{M}}(r<10)=q_{\mathrm{i}\mathrm{n}}$$(12)

where $r=\sqrt{X^2+Y^2+Z^2}$, p_in, p_out, q_in, q_out, and the scale radius, r_q, were allowed to be free parameters.

Given each set of parameters in the gravitational potential, we integrated the stellar orbits using the publicly released AGAMA package (Vasiliev 2019)² for the 11415 stars in the northern hemisphere. We integrated ten orbital periods for each star and withdrew 1000 particles from each orbit with equal time steps, and each particle was given the weight of the star that initialises the orbit. We superposed particles drawn from the orbits together, obtaining an orbit superposed model that numerically represents the DF in 6D phase space of the Milky Way stellar halo.

3.2 Evaluation of the goodness of the model

We detail particular components of the DF: the spatial density distribution and velocity distributions at different positions, and we evaluate these for the model against the data. We calculated the stellar density distribution in the R_gc-z_gc plane, dividing it into 25 × 25 bins within 50 kpc. The density in each bin was calculated by N/(2πRdRdz) in units of N/kpc³. We calculated the density directly for both the observed data, p_data, and the model, p_model. The uncertainty of the data, dρ_data, was assessed by bootstrapping the distances of stars along with their uncertainties. We calculated ${\chi }_{\mathrm{d}\mathrm{e}\mathrm{n}}^2$ by directly comparing the data and the model: $${\chi }_{\mathrm{d}\mathrm{e}\mathrm{n}}^2=\sum ({\rho }_{\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{a}}-{\rho }_{\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l}}{)}^2/d{\rho }_{\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{a}}^2.$$(13)

To analyse velocity distributions, we categorised model particles into N_bin = 7 × 5 bins in the direction of r and θ as is shown in Fig. 3. We constructed the velocity distributions (v_r, v_φ, v_θ) within each bin, j, and represent the distribution of each velocity component by a histogram using intervals of M, labelled ($v_{m,j}^k$, $h_{m,j}^k$), where k corresponds to the velocities v_r, vφ, and v_θ.

In bin j, we have N_j stars, and for a star, i, with a velocity and measurement uncertainty ($v_i^k$, ${\sigma }_i^k$) within that bin, we determined its likelihood compared to the model for each of the three velocity components as follows: $$P_{i,j}^k=\frac{1}{\sqrt{2\pi }\times {\sigma }_i^k}\times \frac{\sum _{m=1}^M{h}_{m,j}^k\times \exp [-(v_i^k-v_{m,j}^k{)}^2/(2({\sigma }_i^k{)}^2)]}{\sum _{m=1}^M(h_{m,j}^k)}.$$(14)

The collective likelihood of all Nj stars in this bin is given by $L_j^k=\sum _{i=1}^{N_{\mathrm{j}}}\log (P_{i,j}^k)$, and the likelihood of combining all bins is represented by $L^k=\sum _{j=6}^{N_{\mathrm{b}\mathrm{i}\mathrm{n}}}{L}_j^k$.

Our data for integrating stellar orbits primarily cover the range 4 < r < 50 kpc, meaning that we lack data for orbits starting at r < 4 kpc. These orbits could significantly influence the DF of stars within r < 8 kpc. Consequently, the estimated DF from the model remains incomplete for r < 8 kpc, whereas it is more than 90% complete for the range 8 < r < 50 kpc. Thus, we excluded the five bins in r = [4, 8] kpc from our likelihood calculations. Although our observations for the southern hemisphere are incomplete, the 3082 stars available should conform to the same DF. We incorporated these stars into the northern spatial bins and included them in our likelihood calculations.

The likelihood across the three velocity components was merged to derive the constraints from the velocity distributions: $${\chi }_{\mathrm{v}\mathrm{d}\mathrm{i}\mathrm{s}}^2=-2(L^{v_r}+L^{v_{\varphi }}+L^{v_{\theta }}).$$(15)

Finally, we combined the χ² values from both the density distribution and the velocity distributions to obtain $${\chi }_{\mathrm{t}\mathrm{o}\mathrm{t}}^2=-2(L^{v_r}+L^{v_{\varphi }}+L^{v_{\theta }})+{\chi }_{\mathrm{d}\mathrm{e}\mathrm{n}}^2.$$(16)

Fig. 3 Spatial binning scheme we used for comparing the velocity distributions of the model and data. We divided the model into 7 × 5 spatial bins along r × θ, and calculated the likelihood of stars located within each bin by comparing to the velocity distributions of the model in the corresponding bin.

Fig. 4 Model constraints on the parameters of the underlying DM distribution, ρ0, rs, γ, qDM , and pDM. The three parameters determining the radial distribution, ρ0, rs, and γ, are degenerate due to the lack of data points at r < 8 kpc and r > 50 kpc. However the 3D DM distribution within our data coverage (r < 50 kpc) are well constrained by our model, with a significant role played by the density distribution, ${\chi }_{\mathrm{d}\mathrm{e}\mathrm{n}}^2$.

Fig. 5 Enclosed mass profiles and rotation curve that we obtained for the Milky Way. The solid line is the average of models within the 1σ confidence level, while the dashed curves are the upper and lower boundary.

4 Model results

4.1 The 3D dark matter distribution of Milky Way

We investigated the parameter space of ρ₀, r_s, γ, p_DM, and q_DM within the gravitational potential. A parameter grid was established with intervals of 0.1, 10 kpc, 0.02, 0.04, and 0.02 for the five parameters, respectively. An iterative approach was used to find the optimal models. The process began with an initial model, and after developing initial models, iterative refinement followed. During each iteration, we identified the optimal models using the criterion χ² - min(χ²) < χ²_s, where X = 100. We then generated new models by walking two steps in every direction of the parameter grid of each optimal model. This approach guided the search towards the lower χ² within the parameter grid and halted when the model with the minimum χ² was identified. The relatively high value of χ²_s was selected to ensure that all models within the confidence level of 1σ were considered before the end of the iteration. Ultimately, we determined the best-fitting models by achieving the minimum ${\chi }_{\mathrm{t}\mathrm{o}\mathrm{t}}^2$.

The 1σ confidence level was determined by bootstrapping. We took the gravitational potential of the best-fitting model and perturbed the position and velocity of observed stars with their uncertainty for 100 trials. In each trial, we recomputed the model and evaluated the resulting ${\chi }_{\mathrm{d}\mathrm{e}\mathrm{n}}^2$, ${\chi }_{\mathrm{v}\mathrm{d}\mathrm{i}\mathrm{s}}^2$, and ${\chi }_{\mathrm{t}\mathrm{o}\mathrm{t}}^2$. The standard deviation of these 100 values of χ² was taken as the 1σ confidence level of the model $\mathrm{\Delta }{\chi }_{1\sigma }^2$.

The parameter grid we explored for the Milky Way is illustrated in Fig. 4 with about 15 000 models calculated. Both the density distribution, ${\chi }_{\mathrm{d}\mathrm{e}\mathrm{n}}^2$, and the velocity distributions, ${\chi }_{\mathrm{k}\mathrm{i}\mathrm{n}}^2$, yield consistent results of the gravitational potential parameters, with ${\chi }_{\mathrm{d}\mathrm{e}\mathrm{n}}^2$ playing a key role in determining the DM shape parameters, p_DM and q_DM. We established the DM halo shape parameters as q_DM = 0.92 ± 0.08 and p_DM = 0.8 ± 0.2, a nearly oblate DM halo, but with the long-intermediate axis plane vertically aligned with the stellar disc (q_DM > p_DM) that is strongly preferred by our model.

The three parameters determining the radial distribution, ρ₀, r_s, and γ, are still degenerate. However, the enclosed DM mass profile within the data coverage is reasonably constrained. We determined the DM mass to be M_DM(< 50kpc) = (5.4 ± 1.6) × 10¹¹ M_⊙. The circular velocity at the solar location (r = 8.2 kpc (Bland-Hawthorn & Gerhard 2016)) was found to be $V_c={220}_{-23}^{+29}$ km/s, which agrees well with previous findings (McMillan 2017). Our derived enclosed mass profile and rotation curve for the Milky Way are presented in Fig. 5.

4.2 The best-fitting models

To illustrate how well the models match the data, we show the stellar density distribution in the R_gc-z_gc plane in Fig. 6. The stellar density distributions of the data and model were constructed in the same manner. We further extracted the flattening of the halo stars, q_star, along the iso-density contours of log(ρ_star) = [−3, −2, −1,0,1,2,3] and compared that from the data and the model. As is shown in Fig. 6, the model with an underlying DM halo of q_DM = 0.92, p_DM = 0.80 constructed a stellar halo density that matches the data well. Variations in the DM halo shape yield significant deviations in the stellar density distribution of the orbit superposition model from the actual data, especially the flattening q_star at different radii. In other words, the density distribution of the halo stars, especially the flattening q_star, places strong constraints on the shape of the underlying DM halo; our results for the 3D shape of DM halo are highly data-driven.

To analyse velocity distributions, we categorised model particles into bins of N_bin = 7 × 5 in r and θ. The divisions are r = [4,8,10,12,15,20,30,50] kpc and θ = [0,15,30,45,60,90°]. We then constructed the velocity distributions (v_r, v_φ, v_θ) within each bin for the data and model, as is illustrated in Fig. 7. Likewise, the model with DM mass of log(ρ₀[M_⊙/kpc³]) = 6.2, r_s = 70 kpc, γ = 1.0 (red) matches the data well, while an altered density profile results in notable differences in the model’s velocity distributions. The velocity distribution of stars, especially v_r, strongly constrains the radial distribution of the underlying DM distributions. We also note that our best-fitting model, which constructed both stellar density and the 3D velocity distributions, agrees with the data generally well, supporting our hypothesis of a ‘stationary system’.

Fig. 6 Comparison of the observational data and the empirical orbit superposed model: the stellar density distribution. The top panels from left to right are the 2D density distribution ρstar in Rgc versus zgc constructed by data, our best-fitting model (qDM = 0.92, pDM = 0.80), a model with a flatter DM halo, and a model with a rounder DM halo. The bottom panels show the relative error of the data, and the corresponding model residuals. The inset panel in the right shows the flattening of the stellar halo, qstar ≡ zgc/Rgc, as a function of the radius, Rgc, extracted from the iso-density contours of the 2D density distribution. The star symbols with the error bar are data. The solid red, light blue, and black curves indicate the three models shown on the left. The dashed red curve indicates a model allowing qDM and pDM to vary as a function of radius. The best-fitting model constructed stellar density distribution that matches the observational data well, while models with flatter or rounder DM halos constructed stellar halos that are either too flat or too round. The density distribution of the halo stars, especially the flattening qstar, has strong constraints on the shape of the underlying DM halo.

4.3 Models with variable p_DM(r), q_DM(r)

The best-fitting models with constant p_DM, q_DM generally succeed in replicating the full DF of the data. There is a small deviation between the 2D stellar density distribution constructed from the data and the model. This deviation is not statistically significant, but it becomes noticeable when we evaluate the flattening of the stellar halo (q_star) along the iso-density contours, as is shown in the right panel of Fig. 6.

We further created a set of models for the Milky Way with p_DM(r), q_DM(r) following Equations (9)-(12). We fixed r_q = 20 kpc and allowed p_in, p_out, q_in, and q_out to be free parameters. The parameters determining the radial DM profile, ρ₀, r_s, and γ, do not degenerate significantly with the 3D shape. To optimise the efficiency of the model, we fixed them at log(ρ₀[M_⊙/kpc³]) = 6.3, r_s =40 kpc, and γ = 1.4. We created about 2000 models to explore the parameter space of p_in, p_out, q_in, and q_out. We obtained models in which q_in < q_out and p_in > p_out are preferred, although statistically not significant. The flattening of the stellar halo (q_star) of a best-fitting model is shown in Fig. 6, which indeed matches the data slightly better. We extracted p_DM and q_DM at r = 15, 30, 50 kpc from each model and compare these with cosmological simulations in the next section.

5 Comparison with cosmological simulations

5.1 Shape of the dark matter halo

We took 198 Milky Way-like or M31-like galaxies (Pillepich et al. 2024) from the TNG50 simulations (Nelson et al. 2019b,a; Pillepich et al. 2019) and 30 galaxies from Auriga (Grand et al. 2017), and we performed further selections to enhance their comparability with the Milky Way. We calculated the flattening of the stellar component, q_star = z_gc/x_gc, of each simulation and defined the disc fraction by considering stars with circularity λ_z > 0.5 as part of the disc. Note here that we did not define a dynamically cold disc with λ_z > 0.7 (Genel et al. 2015) or λ_z > 0.8 (Zhu et al. 2018). Instead, we included some stars on dynamically warm orbits to make the disc comparable to the Milky Way disc that was defined by fitting exponential profiles to the stellar number density distribution (Grand et al. 2017). The Milky Way has f_disc ~ 0.7 and global q_star ~ 0.15 (Bland-Hawthorn & Gerhard 2016). We selected galaxies using a slightly loose threshold with f_disc > 0.6 and q_star < 0.3 to retain a relatively large sample. This selection removed galaxies that have experienced relatively recent major mergers, which create a more prominent bulge or stellar halo than that of the Milky Way. We further removed a few galaxies with ongoing major mergers that have not stabilised their stellar components. This results in 94 TNG50 galaxies and 11 Auriga galaxies.

We calculated the intrinsic shape of the DM halo in the co-ordinate X ≡ {X, Y, Z}, defined exactly the same way as we defined it for the Milky Way in our model. We first defined the Z-axis as being vertical to the stellar disc. Then, we projected the DM halo along the Z-axis and defined the long and short axes in the disc plane as X and Y, respectively. We measured the intrinsic shape of the DM halo in the fixed co-ordinate system and measured the axis ratios p_DM = Y/X and q_DM = Z/X at r = 15,30, 50 kpc, respectively. In this definition, p_DM is restricted to be less than unity, while q_DM is allowed to be either lower or higher.

We took the average of p_DM and q_DM measured at r = 15,30,50 kpc for all 94 TNG50 and 11 Auriga galaxies and compared them with the Milky Way DM axis ratio obtained by our default model with constant p_DM and q_DM. As is shown in Figure 8, this particular DM halo shape of the Milky Way is not typically anticipated by the ΛCDM model for galaxy formation and is a rarity in both TNG50 and Auriga. Nevertheless, a total of 13 out of the 105 TNG50/Auriga galaxies are found within the 3σ confidence interval of our model results for the Milky Way. The three within 1σ, TNG50 501208, TNG50 546474, and Auriga 12, have intrinsically oblate DM halos, but with the long-intermediate axis plane vertically aligned with the stellar disc.

We further show the DM axis ratios of TNG50/Auriga galaxies measured at different radii in Fig. 9 and compare them with the Milky Way DM axis ratio obtained by the model with varying p_DM(r) and q_DM(r). The above three galaxies, with average axis ratios within 1σ, actually have twisted DM halos that are more oblate in the inner regions and become more vertically aligned in the outer regions. TNG50 483594, outside of 1σ but within 3σ, exhibits a similar trend but is intrinsically rounder. Although not statistically significant, the Milky Way DM halo is likely to show a similar trend, as is suggested by our model with varying p_DM(r) and q_DM(r). The radial variation in the DM halo shape that we found, although tentative, is consistent with recent results from the modelling of the Sagittarius stream (Vasiliev et al. 2021).

Fig. 7 Comparison of the observational data and the empirical orbit superposed model: the 3D velocity distributions vr, vφ, and vθ in different spatial bins. We divided the system into 7 × 5 bins in r versus θ, but eliminating the first row with r < 8 kpc and only showing the odd number of columns here. In each sub-panel, the solid and dashed grey curves are the mean and 1σ uncertainty of the distribution constructed from observational data. The red, dark blue, and light blue curves were constructed by models with different radial distributions of the underlying DM mass. The model with a DM mass of log(ρ0[M⊙/kpc3]) = 6.2, rs = 70 kpc, and γ = 1.0 (red) matches the data well. The velocity distribution of stars, especially vr, strongly constrains the radial distribution of the underlying DM distributions.

5.2 Arrangement of the satellite system

The so-called ‘satellite plane’ of the Milky Way is especially notable in the orbits of the most massive satellites (Sawala et al. 2023). We first illustrate the orbital configuration of the massive satellite galaxies in the Milky Way, as well as those in the 13 TNG50/Auriga galaxies whose DM shapes fall within the 3σ confidence level of the Milky Way.

For the Milky Way, we selected the ten brightest satellites and used their observed 3D positions and velocities as initial conditions, as per Pace et al. (2022), to compute their orbital paths within our derived optimal gravitational potential. For each of the TNG50/Auriga galaxies, we took the ten most massive satellites and used their actual trajectories back approximately one orbital period from z = 0 in most cases. For the satellites that have been orbited for many periods and are already very close to the galaxy centre at z = 0, we took an orbital period from 3 Gyr ago for better visualisation.

In Fig. 10, we show these four galaxies that have a nearly oblate outer DM halo vertical to the stellar disc: TNG50 501208, TNG50 546474, Auriga 12, and TNG50 483594. They exhibit vertical satellite planes similar to the Milky Way; the first three are the only galaxies with DM shapes that fall within the 1σ confidence level of the Milky Way. There are no clear satellite planes in the remaining nine galaxies within 3σ (see Appendix A).

In general, a correlation is evident between the DM halo shape and the arrangement of the satellite systems in TNG50 galaxies (see Appendix A). Most TNG50 galaxies have oblate DM halos with long-intermediate axes aligned with the stellar disc and paired with coplanar satellite systems. While galaxies with almost spherical DM halos display satellites distributed in a nearly isotropic fashion, those with DM halos aligned vertically to the stellar disc tend to have satellite systems organised vertically as well. This correlation is expected as a result of the DM halo assembling over time through the accretion of satellite galaxies with similar net angular momentum, consistent with previous findings from other cosmological simulations (Shao et al. 2016).

Fig. 8 DM halo shape, qDM = Z/X versus pDM = Y/X, of the Milky Way in comparison with TNG50 and Auriga galaxies. The dashed red contours indicate the 1σ and 3σ confidence level of the Milky Way DM axis ratios obtained by our models. The red star marks the best-fitting model with qDM = 0.92 and pDM = 0.8. Each black dot is a Milky Way-like galaxy from TNG50 simulations. Each black multiple is a galaxy from Auriga simulations. The plus at the bottom left indicates the typical uncertainty of DM halo shapes measured for TNG50/Auriga galaxies. The four simulations highlighted by red/magenta circles/multiple are those associated with a vertical satellite plane. As we show later, those highlighted by black circles are all simulations lying within a 3σ range of the Milky Way but without significant vertical satellite planes.

5.3 A vertically orientated dark matter halo in alignment with a satellite plane results from a flip of the stellar disc

We further traced the evolutionary history of the four Auriga/TNG50 galaxies, TNG50 501208, TNG50 546474, TNG50 483594, and Auriga 12, and found that the distinct vertical alignment of the DM halo and satellite plane is due to the tilt of the stellar disc with a large angle, which we call a ‘flip of the disc’. To check the evolutionary history, we paid special attention to the merger history, the evolution of the stellar disc, and the shape of the DM halo. We calculated the spin direction of all stars at r < 20 kpc at each snapshot in the history and defined θ_disc as the angle between the spin vector and its projection line in the XY plane. In this definition, we have θ_disc,z=0 = 90°. We then calculated the DM axis ratios at 50 kpc (p_DM = Y/X, q_DM = Z/X) in different snapshots. XYZ is the coordinate defined at z = 0, with the Z -axis perpendicular to the disc, as we defined previously.

As is illustrated in Fig. 11 and quantified in Fig. 12, the last major merger of TNG50 501208 ended t ~ 8.5 Gyr ago. A stellar disc was formed gradually after the major merger. At t ~ 6-7 Gyr ago, the stellar disc was almost aligned in the long-intermediate axis plane of the DM halo; if not perturbed by other mergers, the system would stabilise in this configuration. The galaxy was followed by a series of minor mergers, the most massive one having M_max = 7.1 × 10¹⁰ M_⊙ coming in on an orbit that is retrograde to the disc rotation. At t = 6 Gyr ago, when the massive minor merger reached its pericenter for the first time, the disc started to tilt, initiating continuous tilting in the subsequent 6 Gyr. The disc has since flipped about 80^° until z = 0, with an average rate of ~13° per Gyr, while it was not significantly perturbed. Although the satellite perturbed the DM halo’s shape at pericenter, the DM halo quickly re-stabilised. The mass accreted from this satellite and other minor mergers continuously built the DM halo along the orbit, consistent with the pre-existing DM halo shape; the axis ratios of the DM halo measured at r = 50 kpc remained largely unchanged. The disc flip resulted in the DM halo and stellar disc’s vertical alignment at z = 0.

The satellite galaxy that spiralled inwards brings a substantial net angular momentum. It is important to note that a tilt rate of ~13° per Gyr corresponds to about 2 km/s at a radius of r = 10 kpc, indicating a very minor rigid motion of the disc. This scenario requires torques that transfer less than ≲1% of the angular momentum of the satellite galaxy, which has a total mass comparable to that of the stellar disc. In our analysis of 13 TNG/Auriga galaxies within a 3σ threshold, eight galaxies exhibit significant disc tilts, with tilt angles ranging from 30° to 100°. Seven of them are likely caused by minor mergers with M_max > 10¹⁰ M_⊙, mostly spiralled in (see Appendix A). Disc flips are fairly common in cosmological simulations (Dillamore et al. 2022; Earp et al. 2017; Bett & Frenk 2012). Direct N-body simulations show that the torque on the stellar disc from a satellite galaxy is a promising mechanism that induces the disc tilt (Dodge et al. 2023). Previous analysis of Auriga galaxies shows that disc tilt is correlated with the time duration of close satellite interactions (Gómez et al. 2017a), consistent with our results. It is also found that even low-mass satellites that penetrate the outer regions of a galaxy can significantly perturb and tilt a host galactic disc (Gómez et al. 2017b), due not only to direct tidal perturbation but also to the generation of asymmetric features in the DM halo that can be efficiently transmitted to its inner regions.

In a similar pattern, the evolution of TNG50 546474 is presented in the right panel of Fig. 12. This galaxy underwent a major merger that ended around t ~ 7.5 Gyr ago, immediately followed by the most massive minor merger with M_max = 1.4 × 10¹⁰ M_⊙. At t ~ 7 Gyr ago, after the reformation of the disc post-major merger and stabilisation of the DM halo, the disc and DM halo were nearly coplanar and aligned. Subsequently, the DM halo’s shape saw minimal change, although it gained about 30% of its mass in the upcoming 7 Gyr through minor mergers on orbits consistent with the pre-existing DM halo shape. At the same time, the disc has tilted approximately 60° with an average rate of ~8° per Gyr.

Similar figures for Auriga 12 and TNG50 483594 are shown in Figure 13. Auriga 12 has had more frequent recent mergers. Its stellar disc was coplanar aligned with the DM halo at t ~ 8-10 Gyr ago; the torques from the combination of a few subsequent minor mergers have tilted the disc by ~90°. The DM halo was influenced but has maintained a similar shape since about 10 Gyr ago.

The DM shape of TNG50 483594 is outside the 1σ but within the 3σ confidence level of the Milky Way. This galaxy is special in that it was quickly assembled with many head-on mergers at t ≳ 8 Gyr ago, but without a GSE-like major merger, which resulted in a rather spherical DM halo. Due to frequent mergers in the early Universe, it has quickly assembled stellar mass and halted star formation. With a near-spherical DM halo, its stellar disc was arbitrarily orientated at t = 6-10 Gyr ago. The subsequent minor mergers exhibit a satellite plane whose contribution has not significantly changed the DM halo shape until t ~ 6 Gyr ago, along with a massive minor merger with M_max = 4.5 × 10¹⁰ M_⊙. The minor mergers built the halo along their orbits and tilted the disc by about 30° in another direction, leading to the vertical alignment of the DM halo and the disc at z = 0. In this case, the tilt of the disc is not the only reason; it strengthens the misalignment and contributes to the vertical alignment of the stellar disc and the satellite distribution observed at z = 0. Even without the disc tilt, there is still a misalignment of the stellar disc and the satellite plane due to variations in the filament network along with mass accreted into the galaxy, as is shown in Shao et al. (2021).

Despite the tilting of the disc, the DM halo of the inner regions (r ≲ 20 kpc) was tilting with the stellar disc. As a result, the DM halos of these galaxies tend to be twisted, being oblate in the inner regions, and more vertically aligned with the disc in the outer regions, as has already been shown in Figure 9, consistent with the trend we obtained for the Milky Way.

The three 1σ galaxies, TNG50 501208, TNG50 546474, and Auriga 12, would transform into a typical galaxy with an oblate DM halo by hypothetically tilting their stellar disc (and inner DM halo at r ≲ 20 kpc) back to the long-intermediate axis plane of the DM halo in the outer regions. Note that TNG50 403894 also has a twisted DM halo; however, its transition radius, similar to other cases of twisted halos caused by variations in filament accretion (Shao et al. 2021), is larger than the three cases caused purely by the tilt of the disc (see Fig. 9).

Fig. 9 DM axis ratio measured at three different radii, r = 15,30,50 kpc, for TNG50 and Auriga galaxies, in comparison with that of the Milky Way obtained by our model allowing pDM(r) and qDM(r) vary as a function of radius. This model has a large uncertainty and we only weakly constrain the DM axis ratios at different radii, which shows a trend to be more oblate in the inner regions and to be more vertically aligned in the outer regions, consistent with the trend shown in the Milky Way analogues, TNG50 501208, TNG50 546474, Auriga 12, and TNG50 483294. For these four galaxies, the arrows point to their intrinsic axis ratios measured at the principle axes of the DM halo. These four galaxies have twisted DM halos and are vertically aligned at r > 20 kpc, as we obtained for the Milky Way, bur their intrinsic DM shapes are common compared to the general TNG50/Auriga galaxies if we hypothetically tilt their disc and inner DM halo to the long-intermediate axes plane of the outer DM halo.

Fig. 10 Orbit configuration of satellite galaxies for the Milky Way and the four analogues highlighted in Figs. 8 and 9. We show orbits of the ten brightest or most massive satellites in YZ planes. The blue disc at the top left indicates the orientation of the stellar disc. The analogues exhibit striking similarities with the Milky Way that the DM halo coincides with a plane of satellites and vertically aligns with the Galactic disc.

Fig. 11 Spiral-in minor merger flips the stellar disc illustrated by TNG50 501208. The left column show orbits of the ten most massive mergers coloured their maximum history mass. The dashed red curve with Mmax = 5.3 × 1010 M⊙ is a major merger finished at about 8.5 Gyr ago, and the solid red curve with Mmax = 7.1 × 1010 M⊙ is the most massive minor merger followed. The rest columns show the DM distribution (the grey shadow), the projection of the stellar disc (the blue ellipse with size enlarged five times for visualisation), and the orbit progress of the most massive minor merger (solid red curve) at different snapshots in the history. The axis ratios of the DM halo measured at r = 50 kpc are labeled. A thin stellar disc was formed right after the major merger finished at about 8.5 Gyr ago. At about t = 6-7 Gyr ago the disc was aligned almost co-planar with the long-intermediate axis plane of the DM halo and orbital plane of the satellites. The disc was flipped in the last seven gigayears, associated with the spiral-in minor merger (solid red curve), while the general shape of DM halo (illustrated by qDM and pDM measured at 50 kpc in a fixed co-ordinate) and satellite arrangement remain unchanged, resulting in the stellar disc vertical to the DM halo and satellite orbit plane at present (also see Movie S1).

6 Discussion

6.1 Significance of the results

In our default model with a constant axis ratio of q_DM = Z/X and p_DM = Y/X, the axis ratio q_DM is constrained to be 0.92 ± 0.08, consistent with most previous studies (e.g., Wegg et al. 2019; Huang et al. 2024; Zhang et al. 2025). However, we have also shown that the Milky Way DM halo is triaxial, with q_DM > p_DM strongly preferred. Although there is still some degeneracy between the two parameters, with a relatively large uncertainty on the axis ratio p_DM from our dynamical model, the allowed region in q_DM versus p_DM is already narrow enough for us to distinguish the galaxy formation scenarios. Except for the four Milky Way analogues aforementioned, TNG50 501208, TNG50 546474, TNG50 483594, and Auriga 12, all the other TNG50 galaxies fall within the 3σ confidence level of our model, displaying nearly isotropic or coplanar satellite distributions, as was expected from their intrinsic near-spherical or normally oblate DM halos (see Appendix B).

The presence of a vertical satellite plane thus helps to pin down the shape of the DM halo: a nearly oblate but vertically orientated DM halo is the only type allowed within the 3σ threshold of the Milky Way that is associated with a vertical satellite plane. In three of the four galaxies, the vertical alignment of the DM halo and satellite plane is purely due to the identical physical mechanism: a flip of the stellar disc, likely torqued by minor mergers that spiralled inwards. They are also the only three galaxies with DM halo shapes lying within 1σ range of the Milky Way. In the remaining case, TNG50 483594, its vertical alignment was initialised by variations in the filament network along with mass accreted into the galaxy, but was strengthened by disc tilt. Its DM halo is intrinsically rounder, out of 1σ range of the Milky Way. What is more, it has a very different assembly history from the Milky Way. We therefore strongly propose that the Milky Way arises from the same scenario as the three 1σ galaxies; otherwise, a system combining such a DM halo shape and vertical satellite plane would be highly impossible in either the TNG50 or Auriga simulations.

The merger history of the Milky Way closely resembles that of two of the aforementioned 1σ galaxies, TNG50 501208 and TNG50 546474. It is believed that the stellar disc formed after the last major merger (Belokurov et al. 2020) - the Gaia-Enceladus-Sausage (GSE) merger - occurred 8-10 Gyr ago (Helmi et al. 2018; Belokurov et al. 2018), comparable to the major mergers that happened in these two galaxies 9-10 Gyr ago (highlighted in Fig 12). The prominent Sagittarius stream (Majewski et al. 2003; Belokurov et al. 2006) imprints several spiral-in orbital passages in the halo over the past 6-8 Gyr (Dierickx & Loeb 2017). Recent studies suggest that the progenitor mass of Sagittarius is likely to be 5-6 × 10¹⁰ M_⊙ or even larger (Gibbons et al. 2017; Read & Erkal 2019), comparable to the minor mergers that induced the disc flip in TNG50 501208 and TNG50 546474.

Direct N-body simulations (Dodge et al. 2023) show that the GSE could effectively tilt the Galactic disc, if not destroy it. Although GSE-like major mergers usually result in the coalignment of the disc and DM halo (Dillamore et al. 2022), the ultimate orientation of the disc depends on the orbital properties of the satellite. We find that satellites coming in on prograde orbits usually result in the co-alignment of the disc and orbital plane, while satellites coming in on retrograde orbits could result in a disc vertical to the satellite orbital plane (Cai & Zhu et al., in preparation). The latter is consistent with the case in the MW that GSE remnants are with retrograde rotations (Helmi et al. 2018). We suggest that Sagittarius could also be a promising candidate that induced the disc tilting to vertical alignment, similar to the tilt that began about 6 Gyr ago in TNG50 501208. In real cases, the torques of multiple mergers could have worked together, along with the accretion of unbound DM (Genel et al. 2010).

A diagram is included in Fig. 14 to illustrate the structural evolution we propose for the Milky Way. The current disc-halo orientation of the Milky Way may not be fully stable; the disc could still be tilting. Although there are currently no direct observational constraints, the model of the stream Pal 5’s tidal tails favours a tilting rate of 15° per Gyr for the Milky Way disc (Nibauer et al. 2024), supporting our scenario.

Fig. 12 Evolution history of TNG50 501208 (left) and TNG50 546474 (right). We describe the four panels with TNG50 501208. Top : assembly history of the galaxy. The solid thick curves represent two of the most massive mergers the galaxy experienced. The dashed curves represent small mergers. TNG50 501208 has assembled 80% of its DM mass and 72% of its stellar mass at t = 6 Gyr ago. Second row : distance of the satellite to the galacetic center as a function of time. The solid black curve is a major merger ended at t ~ 8.5 Gyr ago. The red curve is the following most massive minor merger with Mmax = 7.1 × 1010 M⊙. The vertical dashed red line indicates the time at which it reached the pericenter for the first time. Third row : variation in the disc spin direction as a function of time. The disc was co-aligned with the long-intermediate axis plane of the DM halo at t ~ 6 Gyr ago, and it has tilted 80° in the past 6 Gyr. Bottom: evolution of the DM axis ratios, with pDM = Y/X and qDM = Z/X. The DM shape remain largely unchanged in the past 8 Gyr. The four panels in the right show a similar evolution history of TNG50 546474.

Fig. 13 Similar to Fig. 12, but for Auriga 12 (left) and TNG50 483594 (right).

Fig. 14 Diagram illustrating the structural evolution that we have proposed for the Milky Way. The grey ellipse shows the DM halo, with the darker region representing the inner halo associated with the stellar disc. The black dots denote satellite galaxies, while the red curve represents orbit of a massive satellite retrograde to the rotation of stellar disc. Several gigayears ago, the Milky Way featured an oblate DM halo, with the stellar disc, DM halo, and satellite galaxies co-aligned. Around 6-8 Gyr ago, the disc and inner halo began tilting, and this tilt continued over the subsequent gigayears, eventually leading to a twisted DM halo, with the outer DM halo and satellite plane becoming perpendicular to the stellar disc.

6.2 The effects of a tilting disc on dynamical models

If the Galactic disc is still tilting at z = 0, a tilting disc might affect our dynamical models in two aspects: first, the overall potential is evolving, and thus the stellar halo, especially the inner halo, is not necessarily in a stationary state; second, the inner halo may be tilting with the disc and cause a velocity shift and finally, a twist compared to the outer halo, as suggested by the simulations. We evaluate the significance of these effects, considering a disc tilt rate of 15° per Gyr.

The cross time of stars in the inner stellar halo is 0.1-0.2 Gyr. During one period of cross time, the disc could tilt about two degrees, which has negligible effects on the potential. In other words, the potential evolves such that the motion of stars will be affected, but the evolution is so slow that the halo stars can update their motions efficiently to reach dynamical equilibrium in a timely manner. Our assumption of dynamical equilibrium should still hold with any momentum. Moreover, the dynamical equilibrium of the system is indeed supported by the good match between the data and our best-fitting model.

If there is a velocity shift between the inner and outer halo caused by the tilt of the inner halo with the disc, it will lead to some non-equilibrium features in the outer halo. A disc tilt rate of 15° per Gyr equals 0.05 mas/yr, and corresponds to 8 km/s at r = 30 kpc and 13 km/s at r = 50 kpc, which are less than the net velocities we subtracted for the correction of LMC effects. The effects of tilt to the outer halo, if they exist, will be mixed with the LMC effects we subtracted. As we tested, imposing the LMC correction or not does not make a noticeable difference in our final results. We thus do not expect a noticeable difference in our final results if there is a tilt of the inner halo with the disc. In summary, the disc tilt, if it is still ongoing at z = 0, is actually very slow and has minor effects on the kinematics of the stellar halo, which should not affect the results of the dynamical models.

6.3 The tilt angle of the dark matter halo

We fixed the tilt angle of the DM halo, β_q, to be either 0° or 90° compared to the disc, by fixing the Z-axis as perpendicular to the stellar disc. Our current approach allows for a fair comparison with the cosmological simulations by defining exactly the same co-ordinates for measuring the DM shape of the simulated galaxies. In reality, the tilt angle of the DM halo is not likely to be exactly 0° or 90°; a disc vertical to the intermediate axis of the DM halo was suggested to be unstable, while a disc can persist off one of the principal planes of the potential (Debattista et al. 2013). However, the disc tilt angle is strongly degenerated with the axis ratio q_DM for the measurement. A tilt angle of β_q ~ 15° of the DM halo compared to the stellar disc was preferred in creating the disc warp (Han et al. 2023) and the presence of a tilted stellar halo dominated by GSE (Han et al. 2022b). The tilted stellar halo was used to constrain an orbitsuperposition model (Dillamore & Sanders 2025), in which the DM tilt angle is allowed to be free. They obtain an average tilt angle of β_q ~ 45° for the DM halo, but there is a strong degeneracy between the two parameters such that β_q ~ 30-80° and q_DM ~ 0.7-0.95 are allowed within their 1σ confidence level. Note that all the above results are mainly constrained by GSE at r ≲ 30 kpc. A tilted DM halo in the inner regions is generally consistent with our results; we find tentative evidence for a radially varying shape of the DM halo, which being less tilted in the inner regions and becomes vertically aligned in the outer regions, consistent with Vasiliev et al. (2021). Data with larger spatial coverage are needed to further constrain the twist of the DM halo. Interestingly, a recent study has constructed the 3D distribution of stellar halo out to about 50 kpc using Desi observations, which found the Milky Way stellar halo is twisted and becomes vertical to the stellar disc in the outer regions (Li et al. 2025), which turns out to be consistent with the shape of the DM halo we found.

7 Summary

We determined the 3D DM distribution of the Milky Way by utilising the empirical triaxial orbit-superposition model applied to the 6D phase-space data of halo stars provided by LAMOST and Gaia. This model accommodates a highly flexible triaxial DM halo, where the axis ratios may change with varying radii. Both the radial density profile and the 3D shape of the DM halo are constrained simultaneously. Furthermore, we benchmark our findings against 105 galaxies resembling the Milky Way from cosmological simulations. Our major findings are as follows.

• We discover a triaxial, nearly oblate DM halo with qDM = Z/X = 0.92 ± 0.08, p = Y/X = 0.8 ± 0.2 averaging r ≲ 50 kpc, where the Z-axis is defined perpendicular to the stellar disc. The axis ratio qDM > pDM is strongly preferred; the long-intermediate axis plane of the DM halo is unexpectedly vertical to the Galactic disc;

• We determine the DM mass to be M_DM(< 50kpc) = (5.4 ± 1.6) × 10¹¹ M_⊙, and the circular velocity at the solar location to be $V_c={220}_{-23}^{+29}$ km/s, which agrees well with previous research;

• This particular DM halo shape of the Milky Way is not typically anticipated by the ΛCDM model for galaxy formation and is a rarity in both TNG50 and Auriga. Among the 105 galaxies in TNG50/Auriga, only 13 are within the 3σ confidence level of our Milky Way findings, and just three fall within the 1σ range. Notably, all three 1σ galaxies display a plane of satellites that aligns with the long-intermediate plane of the DM halo and is perpendicular to the stellar disc, mirroring the Milky Way;

• This striking configuration of the DM halo and satellite system strongly suggests that the Galactic disc has flipped, torqued by minor mergers, from its original alignment with the DM halo and satellite plane, as is evidenced by Milky Way analogues we found in Auriga and TNG50;

• We find tentative evidence that the Milky Way DM halo is twisted, consistent with alignment with the disc in the inner r ≲ 20 kpc, and becomes vertically orientated in the outer regions, consistent with the prediction of disc flip scenario.

With the 3D shape of the DM halo uncovered, we establish a coherent scenario within the ΛCDM framework for the origin of the vertically orientated DM halo and the satellite plane of the Milky Way as a result of the disc flip. This addresses the vertical alignment of the satellite plane problem in a nontrivial but straightforward way. At this point, our Milky Way is special but not an outlier in the ΛCDM framework. Additionally, the DM shape and disc orientation have broad implications; they potentially impact the development of other Milky Way formations, including disc warps (e.g., Yaaqib et al. 2025) and stellar streams (e.g., Nibauer et al. 2024). Our work may trigger a significant shift in understanding the formation and evolution of the Milky Way.

Data availability

The movie associated to Fig. 11 is available at https://www.aanda.org

Acknowledgements

The authors thank David Hogg, Jianhui Lian, Zhaozhou Li, Wenting Wang for useful discussions. The authors acknowledge the support from the National Key R&D Program of China No. 2025YFF0511002(LZ), 2022YFF0503403 (LZ), 2022YFA1602903 (XK), and No. 2024YFA1611902 (X.X.X), the CAS Project for Young Scientists in Basic Research under grant No. YSBR-062 (LZ, X.X.X, LZhang), and the science research grants from the China Manned Space Project with NO. CMS-CSST-2025-A11 (X.X.X., LZ, CY, LZhang), the Strategic Priority Research Program of Chinese Academy of Sciences grant No. XDB1160102 (X.X.X), and National Natural Science Foundation of China (NSFC) No. 12588202 (X.X.X.).

References

Appendix A Quantifying the spatial arrangement of satellites

We further quantify the distribution of the overall population of satellite galaxies for the 94 selected Milky Way-like galaxies from TNG50 and the 11 galaxies from Auriga. We include all satellites at z = 0 that are more massive than 10⁸ M_⊙ and located within 200 kpc of the halo, and we characterise the arrangement of the satellite system at z = 0 by defining two parameters: the fraction of satellites in vertically aligned orbits, denoted as f_Sat,vertical, and the anisotropy of their distribution in the azimuthal angle, represented by ξ_Sat,φ.

Fig. A.1 The satellite distribution in θL and φr for the Milky Way (Top) and TNG50 501208 (Botton) illustrating the definition of fSat,vertical and ξSat,φ. The parameter fSat,vertical is defined as the fraction of satellites with θL < 45°. We fit the distribution of φr with a sine function n = n1 sin(2φr + φ0) + n0. The parameter ξsat,φ = n1/n0 is defined to describe the anisotropy of satellite galaxies distributed in azimuthal angle φr.

First, we compute their angular momentum; the orientation of the angular momentum θ_L is determined by the angle between the spin direction and its projection line on the disc plane. A satellite with θ_L ≈ 0 is in an orbit vertically aligned with the stellar disc. The f_Sat,vertical is defined as the fraction of satellites with θ_L < 45°. We evaluated the azimuthal anisotropy of the satellite distribution by calculating their azimuthal angle φ_r in the disc plane. The distribution of all satellites’ φ_r is plotted, and we fit the histogram with a sine function n = n₁ sin(2φ_r + φ₀) + n₀. An anisotropy parameter, ξ_Sat,φ = n₁/n₀, is defined. The definition of these two parameters is illustrated in Fig. A.1. Under this definition, a perfectly isotropic satellite system has f_Sat,vertical = 0.707 and ξ_Sat,φ = 0.

For the Milky Way, we took 52 satellite galaxies with measured 3D position and velocity (Pace et al. 2022), and added LMC to the data, thus with 53 in total. As shown in Figure A.1, the satellite galaxies of the Milky Way are vertically orbiting with f_Sat,vertical = 0.87 and tend to align in a plane with high azimuthal anisotropy ξ_Sat,φ = 0.55.

The results of these 94 TNG50 and 11 Auriga galaxies are presented in Fig. A.2. As shown in the figure, the distribution of satellite galaxies around most galaxies is not entirely isotropic. Note that these parameters are affected by the number of satellites and there are still certain scatters in describing the satellite arrangement. However, we notice that there is a general correlation between the shape of the DM halo and the arrangement of the satellite system; the correlation is stronger for the DM shape measured at a larger radius. In galaxies with oblate DM halos (q_DM < p_DM), satellite galaxies tend to align coplanarly with the stellar disc, characterised by f_Sat,vertical < 0.707 and ξ_Sat,φ ≲ 0.3. This configuration is typical for systems with non-zero total angular momentum, where the orientation of the stellar disc, DM halo, and satellite system are horizontally synchronised, consistent with previous results (Shao et al. 2016). There are a number of galaxies with satellites distributed highly coplanar with f_Sat,vertical < 0.5; a plane of satellites like the Milky Way (not considering its orientation) is thus not so rare, as also revealed by previous studies (Sawala et al. 2023). However, galaxies with triaxial or near-spherical DM halos (q_DM ≃ p_DM) show a more isotropic distribution of satellites (f_Sat,vertical ~ 0.707). Meanwhile, in the rare cases of galaxies with DM halo vertically aligned with the stellar disc (q_DM > p_DM), their satellite systems appear vertically orientated f_Sat,vertical > 0.707, and at the same time show considerable anisotropy in the azimuthal direction (large value of ξ_Sat,φ).

For Milky Way and the 13 TNG50/Auriga galaxies within the 3σ confidence level, this analysis of the overall satellite populations with f_Sat,vertical and ξ_Sat,φ is generally consistent with the orbital configurations of the ten brightest/most massive satellites (see Fig. A.3). The Milky Way, TNG50 501208, TNG50 546474, and TNG50 483594 show clear planes in the orbits of their ten most massive satellites, while their overall satellite populations are also vertically aligned and highly anisotropic. Auriga 12 also exhibits a plane formed by the ten most massive satellites; its overall satellite populations are also vertically aligned, whereas the azimuthal anisotropy is not so significant. Among these, TNG50 501208 and TNG50 546474 are the closest to the Milky Way when considering either the arrangement of the ten most massive satellites or the overall satellite populations. For the nine remaining galaxies located in 3σ, there are no clear vertical satellite planes formed by the ten most massive satellites, and their overall satellite populations are generally distributed rather isotropically or coplanar aligned, as revealed by f_Sat,vertical and ξ_Sat,φ.

Appendix B The other TNG50 and Auriga galaxies located within 3σ confidence level

Despite the four galaxies with DM halos and satellite planes vertically aligned due to disc tilt, there are another nine galaxies, TNG50 392277, 402555, 530852, 539333, 550149, 550475, 566365, 571454, and 580406 with DM shape outside of the 1σ but within the 3σ confidence level of the Milky Way. These galaxies have intrinsically near-spherical DM halos. Four of them, TNG50 530852,550149,571454, 580406 also experienced significant disc tilt with tilt angles from 30° to 100°. The disc tilts of the former three galaxies are associated with and likely induced by a massive minor merger. Such events occurred quite often and were not necessarily associated with systems exhibiting satellite planes. The disc of TNG50 580406 was tilted about 30° and was not associated with any significant mergers.

Fig. A.2 Arrangement of satellite system in the Milky Way and the Milky Way-like simulated galaxies, and relation with the 3D shape of DM halo. Left: the configuration of satellite system of each galaxy is characterised by two parameters: the fraction of satellites moving on vertical orbits fSats,verticai, and the azimuthal anisotropy of the satellite distribution ξSat,φ. Each dot is one galaxy coloured by their DM halo shape qDM - pDM (< 50 kpc). The red star represents the Milky Way. The 13 TNG50/Auriga galaxies having DM halo shapes within 3σ confidence level of the Milky Way are highlighted. Right: correlation of the DM halo shape qDM - pDM versus satellite distribution fSats,vertical of the 94 TNG50 galaxies, R in the legend is the Pearson correlation coefficiant. The correlation becomes stronger when increasing the number of satellites by taking a lower mass threshold and measuring DM shape at a larger radius.

However, all nine of these galaxies have satellites near isotropically, as shown in the right panels of Figure A.3, distinguished from the vertical satellite plane as observed in the Milky Way. This is consistent with the correlation between the DM shape and satellite arrangement we showed in the right panel of Figure A.2; satellite galaxies in near-spherical halos tend to be isotropically distributed.

Such galaxies with a near-spherical DM halo, although lying within 3σ confidence of our model for the Milky Way, are not likely the case for the Milky Way considering their satellite arrangements.

Fig. A.3 Orbital configuration of satellites in the Milky Way and TNG50/Auriga galaxies within the 3σ confidence level. For each galaxy, we show the orbits of the ten brightest/most massive satellites. For the Milky Way, the pluses represent the orbit of the Sagittarius stream. The colours represent the relative mass of the satellites in each galaxy, red is the most massive one; the thickness of the curves are scaled with the maximum history mass of the satellites. The blue ellipses in the top panels indicate the orientation of stellar disc. The first four simulated galaxies in the topleft panels are the four galaxies we discussed, they show similar vertical satellite plane as the Milky Way. While the rest nine simulated galaxies show rather isotropic distribution of the satellite galaxies. Eight of the 13 galaxies (highlighted by red boxes), experienced disc tilt from 30° to 100°, most of them are likely induced by a massive minor merger (red curves), mostly spiraled inward on tangential orbits, except 580406.

All Figures

	Fig. 1 Density of stars in the parameter phase of circularity, λz, versus metallicity, [Fe/H]. We first removed the disc stars with a circularity of λz > 0.8, then also removed the group of stars with λz > 0.3 and [Fe/H] > -1.2. The rest were kept as halo stars in our sample.
In the text

	Fig. 2 K giants in our final smooth halo sample in the northern hemisphere, coloured by the weights of particles obtained by correction of the selection function. The sample can be taken as complete in the Rgc-zgc space after correction of the selection function.
In the text

	Fig. 3 Spatial binning scheme we used for comparing the velocity distributions of the model and data. We divided the model into 7 × 5 spatial bins along r × θ, and calculated the likelihood of stars located within each bin by comparing to the velocity distributions of the model in the corresponding bin.
In the text

Fig. 4 Model constraints on the parameters of the underlying DM distribution, ρ0, rs, γ, qDM , and pDM. The three parameters determining the radial distribution, ρ0, rs, and γ, are degenerate due to the lack of data points at r < 8 kpc and r > 50 kpc. However the 3D DM distribution within our data coverage (r < 50 kpc) are well constrained by our model, with a significant role played by the density distribution, ${\chi }_{\mathrm{d}\mathrm{e}\mathrm{n}}^2$.

In the text

	Fig. 5 Enclosed mass profiles and rotation curve that we obtained for the Milky Way. The solid line is the average of models within the 1σ confidence level, while the dashed curves are the upper and lower boundary.
In the text

Fig. 6 Comparison of the observational data and the empirical orbit superposed model: the stellar density distribution. The top panels from left to right are the 2D density distribution ρstar in Rgc versus zgc constructed by data, our best-fitting model (qDM = 0.92, pDM = 0.80), a model with a flatter DM halo, and a model with a rounder DM halo. The bottom panels show the relative error of the data, and the corresponding model residuals. The inset panel in the right shows the flattening of the stellar halo, qstar ≡ zgc/Rgc, as a function of the radius, Rgc, extracted from the iso-density contours of the 2D density distribution. The star symbols with the error bar are data. The solid red, light blue, and black curves indicate the three models shown on the left. The dashed red curve indicates a model allowing qDM and pDM to vary as a function of radius. The best-fitting model constructed stellar density distribution that matches the observational data well, while models with flatter or rounder DM halos constructed stellar halos that are either too flat or too round. The density distribution of the halo stars, especially the flattening qstar, has strong constraints on the shape of the underlying DM halo.

In the text

Fig. 7 Comparison of the observational data and the empirical orbit superposed model: the 3D velocity distributions vr, vφ, and vθ in different spatial bins. We divided the system into 7 × 5 bins in r versus θ, but eliminating the first row with r < 8 kpc and only showing the odd number of columns here. In each sub-panel, the solid and dashed grey curves are the mean and 1σ uncertainty of the distribution constructed from observational data. The red, dark blue, and light blue curves were constructed by models with different radial distributions of the underlying DM mass. The model with a DM mass of log(ρ0[M⊙/kpc3]) = 6.2, rs = 70 kpc, and γ = 1.0 (red) matches the data well. The velocity distribution of stars, especially vr, strongly constrains the radial distribution of the underlying DM distributions.

In the text

Fig. 8 DM halo shape, qDM = Z/X versus pDM = Y/X, of the Milky Way in comparison with TNG50 and Auriga galaxies. The dashed red contours indicate the 1σ and 3σ confidence level of the Milky Way DM axis ratios obtained by our models. The red star marks the best-fitting model with qDM = 0.92 and pDM = 0.8. Each black dot is a Milky Way-like galaxy from TNG50 simulations. Each black multiple is a galaxy from Auriga simulations. The plus at the bottom left indicates the typical uncertainty of DM halo shapes measured for TNG50/Auriga galaxies. The four simulations highlighted by red/magenta circles/multiple are those associated with a vertical satellite plane. As we show later, those highlighted by black circles are all simulations lying within a 3σ range of the Milky Way but without significant vertical satellite planes.

In the text

Fig. 9 DM axis ratio measured at three different radii, r = 15,30,50 kpc, for TNG50 and Auriga galaxies, in comparison with that of the Milky Way obtained by our model allowing pDM(r) and qDM(r) vary as a function of radius. This model has a large uncertainty and we only weakly constrain the DM axis ratios at different radii, which shows a trend to be more oblate in the inner regions and to be more vertically aligned in the outer regions, consistent with the trend shown in the Milky Way analogues, TNG50 501208, TNG50 546474, Auriga 12, and TNG50 483294. For these four galaxies, the arrows point to their intrinsic axis ratios measured at the principle axes of the DM halo. These four galaxies have twisted DM halos and are vertically aligned at r > 20 kpc, as we obtained for the Milky Way, bur their intrinsic DM shapes are common compared to the general TNG50/Auriga galaxies if we hypothetically tilt their disc and inner DM halo to the long-intermediate axes plane of the outer DM halo.

In the text

Fig. 10 Orbit configuration of satellite galaxies for the Milky Way and the four analogues highlighted in Figs. 8 and 9. We show orbits of the ten brightest or most massive satellites in YZ planes. The blue disc at the top left indicates the orientation of the stellar disc. The analogues exhibit striking similarities with the Milky Way that the DM halo coincides with a plane of satellites and vertically aligns with the Galactic disc.

In the text

Fig. 11 Spiral-in minor merger flips the stellar disc illustrated by TNG50 501208. The left column show orbits of the ten most massive mergers coloured their maximum history mass. The dashed red curve with Mmax = 5.3 × 1010 M⊙ is a major merger finished at about 8.5 Gyr ago, and the solid red curve with Mmax = 7.1 × 1010 M⊙ is the most massive minor merger followed. The rest columns show the DM distribution (the grey shadow), the projection of the stellar disc (the blue ellipse with size enlarged five times for visualisation), and the orbit progress of the most massive minor merger (solid red curve) at different snapshots in the history. The axis ratios of the DM halo measured at r = 50 kpc are labeled. A thin stellar disc was formed right after the major merger finished at about 8.5 Gyr ago. At about t = 6-7 Gyr ago the disc was aligned almost co-planar with the long-intermediate axis plane of the DM halo and orbital plane of the satellites. The disc was flipped in the last seven gigayears, associated with the spiral-in minor merger (solid red curve), while the general shape of DM halo (illustrated by qDM and pDM measured at 50 kpc in a fixed co-ordinate) and satellite arrangement remain unchanged, resulting in the stellar disc vertical to the DM halo and satellite orbit plane at present (also see Movie S1).

In the text

Fig. 12 Evolution history of TNG50 501208 (left) and TNG50 546474 (right). We describe the four panels with TNG50 501208. Top : assembly history of the galaxy. The solid thick curves represent two of the most massive mergers the galaxy experienced. The dashed curves represent small mergers. TNG50 501208 has assembled 80% of its DM mass and 72% of its stellar mass at t = 6 Gyr ago. Second row : distance of the satellite to the galacetic center as a function of time. The solid black curve is a major merger ended at t ~ 8.5 Gyr ago. The red curve is the following most massive minor merger with Mmax = 7.1 × 1010 M⊙. The vertical dashed red line indicates the time at which it reached the pericenter for the first time. Third row : variation in the disc spin direction as a function of time. The disc was co-aligned with the long-intermediate axis plane of the DM halo at t ~ 6 Gyr ago, and it has tilted 80° in the past 6 Gyr. Bottom: evolution of the DM axis ratios, with pDM = Y/X and qDM = Z/X. The DM shape remain largely unchanged in the past 8 Gyr. The four panels in the right show a similar evolution history of TNG50 546474.

In the text

	Fig. 13 Similar to Fig. 12, but for Auriga 12 (left) and TNG50 483594 (right).
In the text

Fig. 14 Diagram illustrating the structural evolution that we have proposed for the Milky Way. The grey ellipse shows the DM halo, with the darker region representing the inner halo associated with the stellar disc. The black dots denote satellite galaxies, while the red curve represents orbit of a massive satellite retrograde to the rotation of stellar disc. Several gigayears ago, the Milky Way featured an oblate DM halo, with the stellar disc, DM halo, and satellite galaxies co-aligned. Around 6-8 Gyr ago, the disc and inner halo began tilting, and this tilt continued over the subsequent gigayears, eventually leading to a twisted DM halo, with the outer DM halo and satellite plane becoming perpendicular to the stellar disc.

In the text

Fig. A.1 The satellite distribution in θL and φr for the Milky Way (Top) and TNG50 501208 (Botton) illustrating the definition of fSat,vertical and ξSat,φ. The parameter fSat,vertical is defined as the fraction of satellites with θL < 45°. We fit the distribution of φr with a sine function n = n1 sin(2φr + φ0) + n0. The parameter ξsat,φ = n1/n0 is defined to describe the anisotropy of satellite galaxies distributed in azimuthal angle φr.

In the text

Fig. A.2 Arrangement of satellite system in the Milky Way and the Milky Way-like simulated galaxies, and relation with the 3D shape of DM halo. Left: the configuration of satellite system of each galaxy is characterised by two parameters: the fraction of satellites moving on vertical orbits fSats,verticai, and the azimuthal anisotropy of the satellite distribution ξSat,φ. Each dot is one galaxy coloured by their DM halo shape qDM - pDM (< 50 kpc). The red star represents the Milky Way. The 13 TNG50/Auriga galaxies having DM halo shapes within 3σ confidence level of the Milky Way are highlighted. Right: correlation of the DM halo shape qDM - pDM versus satellite distribution fSats,vertical of the 94 TNG50 galaxies, R in the legend is the Pearson correlation coefficiant. The correlation becomes stronger when increasing the number of satellites by taking a lower mass threshold and measuring DM shape at a larger radius.

In the text

Fig. A.3 Orbital configuration of satellites in the Milky Way and TNG50/Auriga galaxies within the 3σ confidence level. For each galaxy, we show the orbits of the ten brightest/most massive satellites. For the Milky Way, the pluses represent the orbit of the Sagittarius stream. The colours represent the relative mass of the satellites in each galaxy, red is the most massive one; the thickness of the curves are scaled with the maximum history mass of the satellites. The blue ellipses in the top panels indicate the orientation of stellar disc. The first four simulated galaxies in the topleft panels are the four galaxies we discussed, they show similar vertical satellite plane as the Milky Way. While the rest nine simulated galaxies show rather isotropic distribution of the satellite galaxies. Eight of the 13 galaxies (highlighted by red boxes), experienced disc tilt from 30° to 100°, most of them are likely induced by a massive minor merger (red curves), mostly spiraled inward on tangential orbits, except 580406.

In the text