© The Authors 2025

1 Introduction

The three-dimensional (3D) distribution of dark matter (DM) within the Milky Way (MW) is crucial to understanding its formation history. Despite extensive studies, there remain significant uncertainties regarding its total mass, density profile, and, in particular, its debated 3D shape.

Streams, viewed as simple orbital pathways that trace the halo, played a significant role in constraining the shape of the MW DM halo. Several approaches have aimed to constrain the flattening of the DM using the two cold streams, Pal 5 and GD-1, which generally suggest a consistent result of q = 0.8–0.95, assuming an axisymmetric halo (Bovy et al. 2016; Bowden et al. 2015; Koposov et al. 2010; Küpper et al. 2015; Malhan & Ibata 2019). These streams trace the inner halo at r ≲ 20 kpc. However, some other streams in similar regions indicate a preference for a more prolate halo. Recent research involving the Heimi stream has identified a mildly triaxial DM halo, with parameters p = 1.013 and q = 1.204 within the inner 20 kpc (Woudenberg & Helmi 2024). Another stream, M68 (Palau & Miralda-Escudé 2019), might also favour a prolate halo. The combination of data from NGC 3201, M68, and Palomar 5 fit well with a halo axis ratio of q = 1.06. (Palau & Miralda-Escudé 2023).

In the outer halo, the Sagittarius stream, being the most notable structure, is analysed using a variety of methods, which results in different outcomes, with q ranges from 0.44 to 1.3 (Law & Majewski 2010; Johnston et al. 2005; Helmi 2004; Ibata et al. 2001). These inconsistent results could be mainly due to different parametrisations of the potential and treatment of the influence of the LMC on the outer portion of the stream (Vera-Ciro & Helmi 2013). Recent research shows that when modelling the Sagittarius stream in a dynamic potential, acknowledging our Galaxy’s motion towards the massive in-falling LMC, the shape of the DM halo is not strongly constrained by the Sagittarius stream (Vasiliev et al. 2021).

Conversely, halo stars - with full 6D phase-space information observed (3D position and 3D velocities) – span the entire space and, statistically, should exhibit stronger power in constraining the underlying DM distribution. However, these stars exist with complex orbital dynamics, necessitating accurate dynamical modelling. The current dynamical models applied to the MW halo usually impose ad hoc assumptions that might bias the results.

A commonly used dynamical model for the MW halo is the Jeans model, which often assumes either a spherical or axisymmetric configuration. Spherical Jeans models, which utilise giant stars in the halo extending to large radii of r ~ 100 kpc, provide strong constraints on the MW halo mass (e.g. Xue et al. 2008; Bird et al. 2022), but do not restrict its shape by definition. The axisymmetric assumptions in the Jeans model generally yield mass estimates for the MW that are consistent with those of the spherical model (see review by Wang et al. 2020). However, the results diverge with respect to the DM halo’s shape. Using RR Lyrae stars, research by Loebman et al. (2014) identified an oblate halo with q = 0.4 ± 0.1 at r < 20 kpc, but Wegg et al. (2019) found a nearly spherical halo with q = 1.00 ± 0.09 at 5 < r < 20 kpc. The alignment of the velocity ellipsoid in the traditional axisymmetric Jeans model differs from observations of MW halo stars, potentially leading to bias. Recently, an axisymmetric Jeans model that allows for a generic shape and radial orientation of the velocity ellipsoid was developed (Nitschai et al. 2020) and applied to MW BHB and K-giant stars, resulting in a near spherical DM halo with q = 0.9−0.95 (Zhang et al. 2025). Nevertheless, Jeans models are still limited in axisymmetric assumption and do not use the full 6D phase-space information, as they often only fit the first and second velocity moments and assume a specific tracer density.

To effectively use the 6D phase-space data, distribution functions (DFs) based dynamical models have been developed and implemented for the MW. A DF-based dynamical model applied to globular clusters (GCs) by Posti & Helmi (2019) identifies a prolate halo with q = 1.30 ± 0.25. In contrast, using DF-based axisymmetric models on halo stars reveals an oblate halo with q > 0.963 within r < 30 kpc, as determined by RR Lyrae stars (Hattori et al. 2021; Li & Binney 2022). These models presume an analytically assumed DF. A new approach named the empirical DF was recently proposed; it uses the observed 6D data of the tracers to construct the DF through appropriate smoothing of the observational data based on theoretical insights, making the DF data driven, with fewer inherent assumptions. Nevertheless, this model is restricted to spherical symmetry (Han et al. 2016; Li et al. 2025). Recently, a novel method was developed to determine the gravitational potential of the MW using 6D data, using deep-learning tools (Green et al. 2023), making minimal physical assumptions that the system is stationary. This method learns the DF from the data and does not restrict the gravitational potential; however, it has not been applied to the MW halo yet, as it may require × 10⁵ data points to construct the model.

Currently, there is no triaxial dynamical model that properly represents the MW halo. However, triaxial, particle-based made-to-measure (Portail et al. 2017; Long et al. 2013) and orbit-based Schwarzschild models (Khoperskov et al. 2025a) are widely used for nearby galaxies and the MW bulge. In such models, the DF of the tracer is numerically constructed by superposition of stellar orbits. For nearby galaxies, the observational information is integrated along the line of sight, without resolution of single stars; the true stellar orbit distribution is unknown. In the traditional Schwarzschild orbit-superposition method, a comprehensive orbit library is created based on theoretical assumptions, and orbit weights are determined by fitting to the observational data, i.e. line-of-sight velocity distribution across the 2D projected plane (e.g., van den Bosch et al. 2008; Zhu et al. 2018; Thater et al. 2022; Tahmasebzadeh et al. 2024; Vasiliev & Athanassoula 2015). The limited observational data restrict the precision of the DF composed of stellar orbits and the ability to determine the DM shape of nearby galaxies. In the MW, we have measurements of the 6D phase space for individual stars. However, the bulge areas are significantly affected by dust extinction, leading to incomplete observations. Analytical models continue to be used in representing the stellar-density distribution and as model input, and libraries of orbits based on theoretical assumptions, along with the determination of orbit weights through data fitting, remain necessary in modeling the MW bulge (Khoperskov et al. 2025a; Portail et al. 2017; Long et al. 2013; Khoperskov et al. 2025c) and disc (Khoperskov et al. 2025b; Bovy et al. 2018).

In the MW halo, using data from LAMOST and Gaia, a large sample of K-giant stars has been constructed, with distances measured with an accuracy level within ≲15%, line-of-sight velocities with a typical error of ~10 km/s, and accurate proper motions with a typical error of ~30 km/s for tangential velocities. These stars cover a substantial portion of the halo, especially in the northern hemisphere, and the selection function is accurately calibrated for volumes up to 50 kpc (Liu et al. 2017; Yang et al. 2022). The DF function of the stellar halo is thus observationally known.

In this paper, we introduce an empirical, triaxial orbit-superposition model for the MW halo. Compared to traditional made-to-measure and Schwarzschild methods, we did not theoretically sample the orbit libraries or determine the orbit weights by fitting the model to data. Instead, we integrated the orbits of stars with 6D phase-space information observed and took the weights of stars from the correction of the selection function; i.e. the orbit library and orbit weights are entirely determined from the observations. We took the minimum physical assumption that the stellar halo is stationary, meaning the DF of stars integrated in a correct potential matches that of stars directly observed instantaneously; our model results will thus be highly data driven. However, this method requires high-quality observations and is confined to the Milky Way.

We applied the method to mock datasets that mimic the LAMOST and Gaia observations. Constrained by halo stars in 6D phase space, we show that the 3D shape and the radial density distribution of the Milky Way DM halo can be determined simultaneously. The paper is organised in the following way: we describe the details of mock data created from Auriga simulations in Section 2, we introduce the method and its application to the mock data in Section 3, we show the results in Section 4, and we discuss them in Section 5. We conclude in Section 6.

2 Simulations and mock data

We first created mock data mimicking the LAMOST and Gaia observations from Auriga simulations. Then, we processed and refined the data to prepare them for dynamical modelling, following the same procedures applied to real observations.

2.1 Auriga simulations

The simulations used for our study are taken from the Auriga project (Grand et al. 2017, 2019), which is a suite of 40 cosmo-logical magneto-hydrodynamical simulations for the formation of the Milky Way-mass halos. These simulations were performed with the AREPO moving-mesh code (Springel 2010) and follow many important galaxy formation processes such as star formation, a model for the ionising UV background radiation, a model for the multi-phase interstellar medium, mass loss and metal enrichment from stellar evolutionary processes, energetic super-novae and AGN feedback, and magnetic fields (Pakmor et al. 2017). We refer the reader to Grand et al. (2017) for more details. In this study, we selected three galaxies from the Auriga simulation suite at a mass resolution of ~5 × 10⁴ M_⊙ for baryons. The comoving gravitational softening length for the star particles and high-resolution DM particles is set to 500 h⁻¹ cpc. The physical gravitational softening length grows with the scale factor until a maximum physical softening length of 369 pc is reached. This corresponds to z = 1, after which the softening remains constant.

2.2 Coordinates and measurement of true DM shape

Initially, we found the principal axis of the stellar component. The galactocentric coordinate x_gc ≡ {x_gc, y_gc, z_gc} is established such that the z_gc-axis is perpendicular to the disc, while the x_gc and y_gc axes are set somewhat arbitrarily within the plane of the disc, similarly to the Milky Way’s galactocentric coordinate.

The intrinsic shapes of the DM halos in the Auriga simulation were determined by Prada et al. (2019), using the axis ratio measured along the principal axis of the DM halo. In this analysis, we used Auriga 23, 5, and 12 to test our models. These three galaxies exhibit considerable variability in the intrinsic shapes of their DM halos, allowing us to assess our model’s adaptability to various DM-halo configurations.

The DM halo does not typically align precisely with the stars. The principal axes of the DM halo, denoted X ≡ {X, Y, Z}, may be connected to the galactocentric coordinate through a rotation matrix, characterised by three Euler angles, α_q, β_q, and γ_q: $$\begin{array}{l}X=\text{R}{x}_{\text{gc}},\hfill \\ \text{R}\equiv \left(\begin{array}{ccc}{c}_{\alpha }{c}_{\gamma }-s_{\alpha }{c}_{\beta }{s}_{\gamma }& s_{\alpha }{c}_{\gamma }+c_{\alpha }{c}_{\beta }{s}_{\gamma }& s_{\beta }{s}_{\gamma }\\ -c_{\alpha }{s}_{\gamma }-s_{\alpha }{c}_{\beta }{c}_{\gamma }& -s_{\alpha }{s}_{\gamma }+c_{\alpha }{c}_{\beta }{c}_{\gamma }& s_{\beta }{c}_{\gamma }\\ s_{\alpha }{s}_{\beta }& -c_{\alpha }{s}_{\beta }& c_{\beta }\end{array}\right)\hfill \end{array},$$(1)

where c_○, s_○ stands for cos ○, sin ○. The z-axis is tilted by the angle β_q relatively to the z_gc-axis. Meanwhile, the angles α_q (γ_q) define the rotation between the x_gc (X) axes and the intersection line of the x_gcy_gc and XY planes.

In practice, when DM distributions are derived from observational data, the orientation angles and intrinsic shapes of DM are highly degenerate. To simplify this, we set Z = z_gc by assuming β_q = 0^°. The orientation of the DM halo in the x_gcy_gc plane is left unrestricted. By projecting DM particles along the z-axis, we determine the principal axis in the XY plane using particles within r < 100 kpc. The long axis in this plane is labelled as X and the short axis as Y. We define the axis ratio as p_DM = Y/X, and q_DM = Z/X. In this definition, p_DM is restricted to be <1, while q_DM is allowed to be >1.

The intrinsic shape of Auriga halos 23, 5, and 12 measured in such a coordinate XYZ is shown in Fig. 1. The intrinsic shapes measured in our coordinate XYZ are slightly different but close to the intrinsic shapes measured along the principal axes of the DM halo shown in Prada et al. (2019).

Fig. 1 Shapes of DM halos in Auriga 23, 5, and 12, measured in the coordinate system where the z-axis is perpendicular to the stellar disc, and X and Y are the long and short axes of the DM halo in the disc plane. The DM halos of Auriga 23 and 5 are oblate, aligned with the stellar disc (pDM > qDM), and vary little with radius. The DM shape of Auriga 12 varies from being oblate and aligned with the stellar disc (pDM ~ 1 and pDM > qDM) in the inner regions to being vertically aligned with the disc (pDM < qDM and qDM ~ 1) at r ≳ 20 kpc.

Fig. 2 Left: MW observations from LAMOST + Gaia. Right: mock data created for Auriga 23. Each dot represents one star/particle coloured by weight obtained from selection-function correction. There are about 15 000 k giants in the final sample of the smooth halo of the Milky Way. The mock data of Auriga 23 include about 20 000 particles, which generally match the spatial distribution of the MW sample.

2.3 Mock data

By cross-matching LAMOST with Gaia DR3, we achieve a sample of 619 284k giants; with 3D velocity measured, the typical errors of radial velocity υ_r and tangential velocities υ_ϕ and υ_θ are 10, 30, and 30 km/s, respectively. The typical distance error is about 15%. We excluded disc stars chemodynamically to avoid biasing the spatial distribution of the halo. We carefully excluded all substructures, including streams, overdensities, and GES using the friend-of-friend method (Wang et al. 2022; Yang et al. 2019), and we further performed 3σ clipping to remove outliers. In the end, we had 15 000 k giants in the smooth halo with the spatial distribution shown in Fig. 2.

We took three galaxies from the Auriga simulations, Auriga 23, 5, and 12. For each one, we took the following steps to create mock data to mimic LAMOST + Gaia observations for the MW halo.

First, we treated each particle as a single star, transferring their positions and velocities into the ICRS coordinate. We adopted a relative error of 15% for distance, an error of 30 km/s for the tangential velocity in the RA and DEC directions, and an error of 20 km/s for the LOS velocity. The positions and velocities of stars were then perturbed by random numbers, which are normally distributed with dispersions equal to the observational errors.

Second, we performed spatial selection of the stars. We divided the particles into an original bulge with r_gc < 4 kpc, an original disc with z_gc < 4 kpc, r_gc > 4 kpc and Z/Z_⊙ > 0.5, and an original halo with the rest. We first randomly selected 1/100 particles from the disc, 1/3 particles from the halo (1/4 for Auriga 12), and no particles from the bulge, and we then performed further spatial selection by only including stars in the northern hemisphere and within a spatial volume similar to the MW LAMOST + Gaia observations, as shown in Fig. 2. This results in ~30 000 stars with the full 6D phase-space information measured. The sample size is much smaller than the original LAMOST + Gaia sample because we took a sparse sampling of the disc stars.

2.4 Data preparation and correction of selection function

With mock data created in an observational manner, we converted the stars into the galactocentric coordinate and obtain the uncertainty of position and velocity in the new coordinate via a bootstrapping process. We cleaned the data and corrected the select function to make it ready for the model construction, in a way similar to how we dealt with the real observations.

First, we used a position-velocity clustering estimator method to identify substructures; we selected the smooth halo by removing all stars identified in substructures with a link length of 0.84 in the 6D distance (Yang et al. 2019). The stellar halos of Auriga 5 and 23 are rather smooth, with no significant structures; we only performed this step on the mock data of Auriga 12. There are fewer substructures in the halo of Auriga simulations than in the real Milky Way galaxy due to the limited mass resolution in the simulation.

Second, we removed the stars from the disc. We integrated the stellar orbits into an initially estimated potential and calculated the circularity, λ_z, of each orbit. λ_z is defined as the angular momentum around the z-axis normalised by the maximum allowed by circular orbits with the same binding energy: ${\lambda }_z=\overline{L_z}/\left(r{V}_c\right)$, where $\overline{L_z}=\overline{x{\upsilon }_y-y{\upsilon }_x},r=\sqrt{\overline{x^2+y^2+z^2}}$, and $V_c=\sqrt{\overline{{\upsilon }_x^2+{\upsilon }_y^2+{\upsilon }_z^2}}$, taking the average of the points along each orbit. A star in a near circular orbit has λ_z ~ 1, and on a radial orbit it has λ_z ~ 0. We removed the disc stars in a cut that combines circularity and metallicity; all the stars with λ_z > 0.8 or λ_z > 0.5 and Z/Z_⊙ > 0.5 were removed. As shown in Fig. 3, the remaining halo stars have an almost symmetric distribution in λ_z centred at λ_z = 0.

Third, we cleaned the data by 3σ clipping in the velocity distributions. Although we removed the substructures identified by the clustering estimator method, there are still some particles with high tangential velocities, which are unlikely to be part of a smooth halo in dynamical equilibrium. We also performed a 3σ clipping in υ_ϕ versus υ_θ to remove the outliers.

Lastly, we corrected for the selection function. For this proof of concept, we used a simple approach by correcting the selection function in the R_gc – z_gc plane. We divided the R_gc – z_gc plane into 50 × 50 bins in R_gc = [0, 50] kpc and z_gc = [0, 50] kpc. We calculated the density of stars in each volume of 2πR_gcdR_gcdz_gc constructed from our sample (ρ_select) and that of the original halo before performing spatial selection (ρ_original). Any star, i, in the bin, j, will have a weight of w_i = ρ_original,j/ρ_select,j. Such a selection correction is less accurate at z_gc < 4 kpc, where disc stars are mixed in the original halo.

In the end, we were left with about 20 000 stars in the smooth halo; these stars are all located in the northern hemisphere and have the selection function corrected. The final sample of stars with their weights is shown in Fig. 2.

The real selection function of LAMOST is complicated and is corrected by comparing with photometric data, which are considered complete with the magnitude limit (Liu et al. 2017; Yang et al. 2022), resulting in the correction of the selection function in the 3D space, more accurate than what we did here. It is beyond the scope of this paper to test the method of selection function correction. We just wanted the spatial sampling of our mock data to match some major properties of the real observations: (1) the data are not complete in different azimuthal angles, and we have fewer data in the direction of the Galactic centre; (2) the observations at larger distances or with small polar angles are less complete; thus, we obtained larger weights for the stars at these positions. As shown in Fig. 2, the particle weight distribution in our mock data set generally matches the trend in the LAMOST + Gaia MW data.

Fig. 3 Mock data for Auriga 23. The blue dots are disc stars that we excluded, i.e. stars with λz > 0.8 or stars with λz > 0.5 and Z/Z⊙ > 0.5. The black dots are taken as halo stars and kept in our sample. We divided the data into 7 × 6 bins in rgc versus θ with the intervals of rgc = [5, 8, 12, 16, 22, 30, 40, 50] kpc, θ = [0, 10, 20, 30, 40, 55, 90°]; here, we only show three columns to illustrate the data.

3 Model construction

We have a sample of stars in the halo, which can be considered complete in the northern hemisphere at r_gc ≲ 50 kpc after correction of the selection function. That is, we have the DF of halo stars in 6D phase space f (x, υ) numerically known from the observations. We assumed that the smooth halo is stationary and that the DF of stars integrated in the correct potential matches that directly observed instantaneously. Under this assumption, the DF constructed from the northern hemisphere should be representative of the whole system.

We constructed a dynamical model called the empirical orbit-superposition model according to the following steps: (1) we constructed an adaptable gravitational potential for the MW with a few free parameters; (2) in a gravitational potential with given parameters, we calculated the stellar orbits using stars’ positions and velocities measured as initial conditions and constructed a model by superposing the stellar orbits using the weights of stars given by the selection function correction; (3) we extracted stellar-density distribution and velocity distributions from the orbit-superposition model and evaluated the goodness of the model by calculating the likelihood (χ²) of each model against the data; (4) we explored the parameter grid of the gravitational potential to find the best-fitting models with the maximum likelihood (or least χ²). This method, in principle, is similar to the traditional Schwarzschild orbit-superposition method, but we have the true DF (orbit library and orbit weights) known from the observations. The only assumption we make in this model is that the system is stationary, which is the minimum assumption for all dynamical models. We confirmed the feasibility of this assumption for the smooth halo with a good match between the model and the data, as shown in Section 3.4.

3.1 Gravitational potential

We constructed a model of the gravitational potential by including a bulge, a disc, and a DM halo.

We adopted a Sérsic bulge, $${\Sigma }_{\text{bulge}}={\Sigma }_0\exp \left(-b_n{\left(R_{\text{gc}}/a_{\text{bulge}}\right)}^{1/n}\right),$$(2)

and a MiyamotoNagai disc, $${\Phi }_{\text{disk}}=-\frac{M_{\text{disk}}}{\sqrt{R_{\text{gc}}^2+{\left(a_{\text{disk}}+\sqrt{z_{\text{gc}}^2+b_{\text{disk}}^2}\right)}^2}},$$(3)

in which $R_{\text{gc}}=\sqrt{x_{\text{gc}}^2+y_{\text{gc}}^2}$. Our data in the halo do not have strong constraints on the bulge and disc. We fixed the density distribution of the bulge and discs via the true values of the simulation (Grand et al. 2017). We have M_*,bulge = 3.2 × 10¹⁰ M_⊙, a scale radius of a_bulge = 1.71 kpc, and a Sérsic index of n = 1.44, M_*,disk = 5.6 × 10¹⁰ M_⊙. Moreover, we obtain a scale radius of a_disk = 4.9 kpc and a scale height of b_disk = 0.25 kpc for Auriga 23. The disc and bulge masses still have some uncertainty for the MW. We tried models with a 20% higher or lower disc mass for Auriga 23 and found no difference in our results. The true values parametrising the bulge and disc of Auriga 5 and Auriga 12 in Grand et al. (2017) are used in our models.

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

in which $\tilde{r}={\left(p_{\text{DM}}{q}_{\text{DM}}\right)}^{1/3}\sqrt{X^2+{\left(Y/p_{\text{DM}}\right)}^2+{\left(Z/q_{\text{DM}}\right)}^2}$, p_DM, and q_DM are the axis ratios of the DM halo.

The coordinate X ≡ {X, Y, Z} is related to the galactocentric Cartesian coordinate x_gc ≡ {x_gc, y_gc, z_gc} by a rotation matrix parametrised by three Euler angles, α_q, β_q, and γ_q, as shown in Section 2.3. In principle, all three angles can be left free. Taking into account the degeneracy between the DM orientation and the 3D shape, we fixed β_q = 0. The combination of α_q and γ_q determines the DM halo orientation in the XY plane; we fixed γ_q = 0 and allowed α_q to be free in order to test the ability of constraining it. For default models, we also fix α_q = 0.

When it is needed, we allowed p_DM and q_DM to vary as a function of radius with $$p_{\text{DM}}\left(r\right)=\left(p_{\text{in,DM}}+p_{\text{out,DM}}{\left(\frac{r-10}{r_q}\right)}^2\right)/\left(1+{\left(\frac{r-10}{r_q}\right)}^2\right),$$(5) $$p_{\text{DM}}\left(r<10\right)=p_{\text{in,DM}},$$(6) $$q_{\text{DM}}\left(r\right)=\left(q_{\text{in,DM}}+q_{\text{out,DM}}{\left(\frac{r-10}{r_q}\right)}^2\right)/\left(1+{\left(\frac{r-10}{r_q}\right)}^2\right),$$(7) $$q_{\text{DM}}\left(r<10\right)=q_{\text{in,DM}},$$(8)

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

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 started with models for all three mock galaxies with fixed DM orientation (α_q = 0, β_q = 0 and γ_q = 0) and constant p_DM and q_DM. The five parameters in the DM halo, that is, the scale density ρ₀ (or the mass M_halo), the scale radius r_s, the inner density slope γ, p_DM and q_DM, were left as free parameters. For Auriga 12, we further constructed a model with different orientation angles, α_q, and a model that allows p(r) and q(r) to vary as a function of radius.

3.2 Orbit integration

Orbit integration was performed using the publicly released AGAMA package (Vasiliev 2019)¹. We have about 20 000 stars in the cleaned sample of the smooth halo, with their positions and velocities measured from observations, and the weight of each star obtained from the correction of selection. We started with the position and velocity of the stars from observations and integrated ten orbital periods for each star. We withdrew 1000 particles from each orbit with equal time steps, and each particle was given the weight as the star that initialises the orbit. We superposed particles drawn from the orbits together, thus obtaining an orbit-superposed model, which numerically represents the DF of the stellar system.

3.3 Comparing model and data

The DF of the orbit-superposed model, f_model(x, υ|potential), is based on data but incorporates the gravitational potential. If the system is in dynamical equilibrium, the DF of the model with the accurate gravitational potential should be statistically the same as that of the stars that initialised it ƒ_data(x, υ). We detail particular components of the DF: the spatial density distribution and velocity distributions at different positions. We evaluated the goodness of the model by comparing these with the data.

We calculated the stellar-density distribution in the R_gc–z_gc plane, by dividing it into 25 × 25 bins within 50 kpc. The density in each bin is calculated by N/(2πRdRdz) in units of N/kpc³. We calculated the density directly for both the observed data ρ_data and the model ρ_model. The uncertainty of the data dρ_data is assessed by bootstrapping the position of stars with their uncertainties. The density distributions of the data and model are shown in Fig. 4; we calculated ${\chi }_{\text{den}}^2$ by directly comparing the data and model: $${\chi }_{\text{den}}^2={\displaystyle \sum {\left({\rho }_{\text{data}}-{\rho }_{\text{model}}\right)}^2/d{\rho }_{\text{data}}^2.}$$(9)

To analyse velocity distributions, we categorised model particles into bins N_bin = 7 × 6 in r and θ. The divisions are r = [5, 8, 12, 16, 22, 30, 40, 50] kpc and θ = [0, 10, 20, 30, 40, 55, 90°]. We then constructed the velocity distributions (υ_r, υ_ϕ, υ_θ) within each bin, as illustrated in Fig. 3. For each velocity component within bin j, we represent its distribution by a histogram using intervals M, labelled as (${\upsilon }_{m,j}^k,h_{m,j}^k$), where k corresponds to the velocities υ_r, υ_ϕ, and υ_θ.

In bin j, we have N_j stars, and for a star i with a velocity and measurement uncertainty (${\upsilon }_i^k,{\sigma }_i^k$) within that bin, we determine its likelihood in comparison 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{{\Sigma }_{m=1}^M{h}_{m,j}^k\times \exp \left[-{\left({\upsilon }_i^k-{\upsilon }_{m,j}^k\right)}^2/\left(2{\left({\sigma }_i^k\right)}^2\right)\right]}{{\Sigma }_{m=1}^M\left(h_{m,j}^k\right)}.$$(10)

The collective likelihood of all N_j stars in this bin is given by $L_j^k={\displaystyle {\sum }_{i=1}^{N_j}\log \left(P_{i,j}^k\right)}$, and the overall likelihood across all the bins is represented by $L^k={\displaystyle {\sum }_{j=7}^{N_{\text{bin}}}{L}_j^k}$. Due to a deficit of stars at r < 4 kpc, which contributes to the DF of the inner bins, the DF for the innermost bins of our model remains incomplete. Therefore, we omit the first six bins within the r = [5, 8] range in our likelihood calculations. The likelihood across the three velocity components is merged to derive the constraints from the velocity distributions: $${\chi }_{\text{vdis}}^2=-2\left(L^{{\upsilon }_r}+L^{{\upsilon }_{\varphi }}+L^{{\upsilon }_{\theta }}\right).$$(11)

Finally, we combine the χ² values from both the density distribution and velocity distributions to obtain $${\chi }_{\text{tot}}^2=-2\left(L^{{\upsilon }_r}+L^{{\upsilon }_{\varphi }}+L^{{\upsilon }_{\theta }}\right)+{\chi }_{\text{den}}^2.$$(12)

Fig. 4 Comparison of stellar-density distribution from data and model for Auriga 23. The top panels from left to right are composed of observational data and several models with different DM axis ratios of qDM and pDM. The colours represent the number density as indicated by the colour bar. and the dashed curves are iso-density contours. The second row shows the uncertainty of the data, and model residuals. The data are well matched by the model with DM qDM = 0.64 and pDM = 1, which are well consistent with the ground truth. The stellar-density distribution constructed by models are either too round or too flat if the shape of underlying DM changes. The stellar-density distribution strongly constrains the 3D shape of the underlying DM halo.

3.4 Best-fitting models

We investigated the parameter space of ρ₀, r_s, γ, p_DM, and q_DM within the gravitational potential. A parameter grid is established with intervals of 0.1, 5 kpc, 0.02, 0.04, and 0.04 for ρ₀, r_s, γ, p_DM, and q_DM, respectively. An iterative approach is used to find the optimal models. The process begins with an initial model, and, after developing initial models, iterative refinement follows. During each iteration, we identify the optimal models using the criterion ${\chi }^2-\min \left({\chi }^2\right)<{\chi }_s^2$, where ${\chi }_s^2=100$. We then generate new models by walking two steps in every direction of the parameter grid of each optimal model. This approach guides the search towards the lower χ² within the parameter grid and halts when the model with the minimum χ² is identified. The relatively high value of ${\chi }_s^2$ is selected to ensure that all models within a 1σ confidence level are considered before the end of the iteration. Ultimately, we determine the best-fitting models by achieving the minimum ${\chi }_{\text{tot}}^2$.

We illustrate how the model works using Auriga 23 as an example. In Fig. 4, it is shown that the best-fitting model accurately reproduces the stellar density across the R_gc – z_gc plane. Furthermore, it imposes strict constraints on the 3D shape of the underlying DM halo. The shape of the DM halo directly affects the direction of stellar acceleration and thus the shape of the stellar orbits. Variations in the DM halo shape can cause the constructed stellar halo to appear overly flat or excessively round. As shown in Fig. 5, the best-fitting model also reproduces the velocity distributions well in all three components: υ_r, υ_ϕ, and υ_θ. These velocity distributions provide strong constraints on the radial density distribution of the DM halo which directly affects the strength of the acceleration, especially in the radial direction. The different radial profiles of the DM halo lead to notable discrepancies in the velocities between the observational data and the model, particularly in υ_r. Of course, both the stellar-density and velocity distributions are determined by the overall 3D DM distribution, which will be better constrained by combining both.

Our best-fitting model generally reproduces both the density and velocity distributions at various locations in the R_gc – z_gc plane, thus validating our assumption of an overall stationary system with the agreement between the model and the data.

However, these stellar systems are never in perfect equilibrium. Although we have tried our best to clean the substructures, there will inevitably be some residual substructures or non-equilibrium features left in the data. For example, in Fig. 5, there are some peaks/bumps in the velocity distributions that are not matched by the model; some of them can indicate residual sub-structures. The υ_ϕ distributions in the bin at 8 < r < 12 kpc close to the disc plane are not perfectly matched by the model either, which can indicate some non-equilibrium features.

In Fig. 6, we illustrate the parameter space explored for Auriga 23. The points denote models within a confidence level of 3σ, coloured according to the normalised $\Delta {\chi }_{\text{tot}}^2=\left({\chi }_{\text{tot}}^2-\min \left({\chi }_{\text{tot}}^2\right)\right)/\Delta {\chi }_{1\sigma }^2$, where $\Delta {\chi }_{1\sigma }^2$ is the confidence level 1σ determined by bootstrapping². As shown in Fig. 6, the combination of density distribution and velocity distribution strongly constrains the 3D shape of DM p_DM and q_DM. Due to the still limited data coverage, there are still strong degeneracies between the three parameters (ρ₀, r_s, and γ) determining the radial profile.

Nevertheless, the radial density profile within the data coverage is well constrained, as we detail in the next section.

Fig. 5 Comparison of velocity distributions from data and model for Auriga 23. The velocity distributions, υr, υϕ, and υθ are calculated in 7 × 6 bins in r versus θ, but we only show three columns here as labelled. In each panel, the grey solid and dashed curves are the velocity distribution and uncertainty of observational data. The red solid curves are from the best-fitting model. Blue and magenta dashed curves show models with different radial distribution of the underlying DM mass. The velocity distributions, especially υr, strongly constrain the underlying DM radial density distribution. The best-fitting model (red solid) reproduces the velocity distributions in all three components well: υr, υϕ, and υθ.

Fig. 6 Constraints on parameters of underlying DM distribution for Auriga 23. Each dot represents one model coloured by the relative $\Delta {\chi }_{\text{tot}}^2$ compared to the least χ2 model. The parameters are sampled on grid with fixed step sizes. We slightly dithered the positions of the points to make them visible. The solid red, black, and dashed purple contours are the 1σ confidence levels determined by $\Delta {\chi }_{\text{tot}}^2,\Delta {\chi }_{\text{vdis}}^2$ and $\Delta {\chi }_{\text{den}}^2$, respectively. The black horizontal and vertical dashed lines indicate the ground truth. The three parameters determining the radial distribution, p0, rs, γ, are degenerated. The shapes of the DM halo, qDM and pDM, are well constrained and consistent with the ground truth.

4 Results

We applied the model to mock data created from Auriga 23, Auriga 5, and Auriga 12. In this section, we show our recovery of the DM radial density profile and 3D shape for these three galaxies.

4.1 Recovery of DM density profile

We demonstrate how effectively our models recover the DM radial-density profile and the enclosed DM mass in Fig. 7. Despite significant degeneracy between p₀, r_s, and γ, our models accurately capture the radial-density profile within the observational range of 4–50 kpc. The uncertainty is larger at r < 8 kpc due to the absence of data.

The enclosed DM mass profile is also well recovered, with a relative uncertainty of about 10%. It should be noted that the DM mass we obtained in the outer regions is 5%–10% higher than the true values but aligns well with the enclosed mass that combines DM and gas. The circumgalactic medium (CGM) in the halo of the simulation mirrors the DM’s density distribution. The mass of the CGM is compensated by an overestimation of the DM mass without a separate gas contribution in our model of the gravitational potential.

4.2 Recovery of the 3D shape of DM

4.2.1 Models with constant q_DM and p_DM

The three selected galaxies exhibit varied DM halo structures, as illustrated in Fig. 1. The DM halo of Auriga 23 is predominantly oblate, showing minimal change from the inner to the outer regions. The DM halo of Auriga 5 is more triaxial, displaying slight variation. On the other hand, the DM halo of Auriga 12 undergoes notable variations, transitioning from oblate within the inner 20 kpc to a vertical orientation at r > 20 kpc.

Initially, we developed models with constant q_DM and p_DM for all three galaxies. As depicted in Fig. 8, our model successfully recovered both p_DM and q_DM for Auriga 23 and Auriga 5, with relatively small uncertainties. For Auriga 12, the model captured the outer configuration of the DM halo, and there is a significant degeneracy between p_DM and q_DM when the halo is vertically aligned; this is characterised by p_DM < q_DM.

The stellar-density distribution in the R_gc – z_gc plane, especially the flattening of the stellar halo, is crucial to constraining the DM shape. We measure the stellar-halo flattening, q_star, across various radii using isodensity contours and compared these from the data and several models in Fig. 9. The stellar halo of Auriga 23 is notably flat, becoming increasingly so from the centre to the outskirts, and our best-fitting model successfully mimicked q_star at all radii. The stellar halo of Auriga 5 appears rounder, with a tendency to become rounder moving outward; our model approximates the radial variation q_star(r), albeit with moderate inner-region discrepancies. The stellar halo of Auriga 12 is the most spherical, notably rounding from inner to outer areas, with q_star(r) challenging to replicate by a model with a constant DM shape. The failure to reproduce the tendency of q_star(r) suggests a variation of DM shapes from the inner to outer regions.

Fig. 7 Recovery of DM radial-density distribution. The radial-density profile and enclosed mass profile we obtained for Auriga 23, Auriga 5, and Auriga 12, which are well consistent with the ground truth within 1σ uncertainty, are shown.

Fig. 8 Recovery of DM intrinsic shapes by models with constant qDM and pDM. The panels from left to right are Auriga 23, Auriga 5, and Auriga 12. In each panel, the solid red, black dashed, and purple dashed curves indicate the model 1σ confidence level determined by ${\chi }_{\text{tot}}^2,{\chi }_{\text{vdis}}^2$, and ${\chi }_{\text{den}}^2$, respectively. The dots are the ground truth measured from the DM halo of the simulation at different radii. The colours represent the radii where pDM and qDM are measured as indicated by the colour bar. The DM halos of Auriga 23 and Auriga 5 are oblate, with little variation from the inner to outer regions, which are well recovered by our model. The DM halo of Auriga 12 varies from oblate in the inner region to vertically aligned in the outer regions, the current models capture the shape of the DM at the outer regions.

Fig. 9 Reproduction of stellar halo flattening qstar as a function of radius Rgc for models with different underlying DM halos. The three panels from left to right are Auriga 23, 5, and 12. The DM halo shapes of Auriga 23 show little variation within our data coverage. The qstar(r) is well reproduced by a model with constant pDM and qDM close to the ground truth. Auriga 5 shows mild variation, and Auriga 12 shows relatively strong variation of the DM halo shape, from r ~ 10 kpc to r ~ 50 kpc. Especially for Auriga 12, the stellar-halo flattening, qstar(r), at all regions can hardly be reproduced by a model with constant pDM and qDM; in contrast, it can be reproduced well by a model with pDM(r) and qDM(r), which vary with radius, close to its ground truth. The stellar-halo flattening qstar(r) has strong constraints on the underlying DM shape.

4.2.2 Models allowing DM shape to vary from inner to outer regions

We create models with variables p_DM(r) and q_DM(r) as specified in Eqs. (5)–(8) for Auriga 12. We fixed the transition radius r_q = 20 kpc, while the remaining four parameters, q_in,DM, q_out,DM, p_in,DM, and p_out,DM were treated as free parameters. The parameters determining the radial DM profile do not degenerate significantly with the 3D shape. To optimise the efficiency of the model, we fixed them at log(ρ₀[M_⊙/kpc³]) = 6.2, r_s = 40 kpc, and γ = 1.6.

The parameter grid we explored for this model is shown in Fig.10. Despite significant uncertainties, a twisted DM halo that is oblate (p_in,DM > q_in,DM and p_in,DM ~ 1) in the inner regions and vertically aligned in the outer regions (p_out,DM < q_out,DM and q_out,DM ~ 1) is strongly favoured.

We derived p_DM(r) and q_DM(r) from our models within the 1σ confidence level and compared those with the true profiles derived from the simulation. Our model generally recovers the true tendency of the DM shape variation with radius, as shown in the right panels of Fig. 10.

Such models allowing variables p_DM(r) and q_DM(r) are indeed able to reproduce the stellar-halo-density distribution. We took one of the best-fitting models with p_in = 1.04, p_out = 0.84, q_in = 0.6, and q_out = 1.32, which reproduces the surface density distribution well from the inner to outer regions, with the residuals significantly reduced compared to the model with constant p_DM and q_DM as shown in Fig. 11. The tendency of q_star(r) derived from this model matches the data well (the dashed black curve in Fig. 9).

Fig. 10 Left: constraints on qin,DM, pin,DM, qout,DM, and pout,DM for model allowing DM axis ratios to vary as a function of radius for Auriga 12. Right: comparison of pDM(r) and qDM(r) profiles obtained by the model (red) and the ground truth (black). The solid and dashed red curves represent the mean and 1σ uncertainty of the results obtained by the model.

Fig. 11 Comparison of stellar-density distribution constructed from data and different models for Auriga 12. The columns from left to right are the data, best-fitting model with constant pDM and qDM, and best-fitting model with pDM(r) and qDM(r) varying as a function of radius, r. The model with constant pDM and qDM cannot match the density distribution in all radial regions; there are significant residuals in the inner 20 kpc. The consistency between data and model improves with the model allowing variable pDM(r) and qDM(r).

5 Discussion

We fixed the DM orientation with the tilt angle β_q = 0. In reality, the DM halo might not be perfectly aligned with the stellar disc. The three galaxies we tested have DM tilt angles of β_q ~ 10–15°.

The axis ratio, q, recovered from our model is close to the true axis ratio measured along the fixed z_gc axis, not the intrinsic q measured along the principal axes of the DM halo, although the latter two are close to each other when β_q is small. For galaxies with a larger tilt angle, β_q, we expect the model to still recover q measured along the fixed z_gc in a similar way; in this case, it may be significantly different from the intrinsic q measured along the DM principal axes.

The orientation of the DM halo in the x_gcy_gc plane remains largely undetermined. Ideally, if the DM orientation α_q can be constrained, we converge two solutions: one with p < 1 and the correct α_q and one with p′ = 1/p when α_q is misaligned with 90°. In our exploration of the Auriga 12 model that allows a free orientation (α_q) of the DM halo in the x_gcy_gc plane, we failed to find significant constraints on α_q. In particular, the stellar-density distributions created by various α_q values of the DM halo appeared similar, resulting in no constraints on α_q from ${\chi }_{\text{den}}^2$, and the same is true for the radial velocity distribution υ_r. The velocities υ_ϕ and υ_θ are somewhat influenced by the orientation of the DM α_q, but this effect occurs in a somewhat unpredictable manner.

This arises because we averaged out information along the azimuthal angle when correcting the selection function and comparing the data with the model. We had to do so because observational data at various azimuthal angles remain incomplete, particularly with a lack of observations directed towards the Galactic centre. Consequently, this process results in a significant level of uncertainty for p_DM = Y/X, especially when p_DM < q_DM, as demonstrated by Auriga 12. Despite this, the relative positioning of p_DM and q_DM is effectively constrained. In our models for Auriga 12, p_DM < q_DM is strongly preferred in the outer regions, in contrast to Auriga 23 and 5, where p_DM > q_DM is strongly preferred.

The stellar halo of the three galaxies we tested is generally in dynamical equilibrium, but inevitably with some residual substructures and non-equilibrium features; they did not introduce obvious systematic bias according to our test results. The LMC, which fell in recently, may have perturbed the Milky Way halo and will affect the results of dynamical modelling. However, the LMC mainly induces a bulk motion in the outskirts of the halo, while it only mildly affects the inner halo, causing a velocity shift in V_z with ~10 km/s at r < 50 kpc (Erkal et al. 2020). Among the Auriga simulations, Auriga 25 had a recent minor merger with mass and orbits similar to the LMC, but its perturbation of the disc and inner halo is significantly stronger than that of the LMC of the Milky Way. We applied the method to Auriga 25 and showed that even with some obvious non-equilibrium in the data, the DM distribution can still be roughly recovered but with larger statistical uncertainty (see the appendix).

6 Summary

We present an empirical, triaxial orbit-superposition model tailored for the MW halo. In this method, we take the minimum assumption of dynamical models whereby the stellar halo is stationary, meaning that the DF in the 6D phase space, f (x, υ), does not change when the stars orbit in the correct gravitational potential. The feasibility of this assumption can be gauged by the final goodness of the match of the data and the model.

The gravitational potential of the model is a combination of disc, bulge, and DM halo with adaptable parameters. Our approach employs a highly flexible triaxial NFW halo, which allows variations in the DM axis ratios (p_DM = Y/X, q_DM = Z/X) with radius. Within each specified gravitational potential, we take the 6D phase-space information of stars observed as initial conditions, integrate their orbits in the gravitational potential, and form a model by superposing the orbits together. The orbits are weighted by the weights assigned to the stars when correcting the observation selection function. We evaluated the goodness of the model by comparing the DF of the model with that directly constructed from the data, which is expressed in terms of density distribution and velocity distributions (υ_r, υ_ϕ, υ_θ) across the R_gc – z_gc plane.

This methodology was validated using mock observational data of the MW halo from LAMOST + Gaia with three Auriga galaxies, halo 23, 5, 12, with different intrinsic shapes of the DM halo. The mock data were created with similar spatial selection and observational errors as in LAMOST + Gaia. We cleaned the data by subtracting the disc, substructures, and so on, in exactly the same way as we dealt with the real observational data.

We created models with constant DM axis ratios from all three mock galaxies and extra models allowing p_DM(r), q_DM(r) to vary with radius for Auriga 12. We find the following:

• The density distribution and velocity distributions of the stellar halo provide consistent constraints on the DM distribution. Specifically, the stellar-density distribution has stronger constraints on the DM shape, while velocity distributions better restrict the DM radial-density profile;

• The DM radial-density profile and the enclosed DM mass of the three galaxies within the data coverage are well recovered, with a relative uncertainty of about 10%. We still have relatively large uncertainty in each of the three parameters that determine the radial profile (ρ₀, r_s, and γ) with limited data coverage;

• Auriga 23 and Auriga 5 have oblate-to-triaxial DM halos that vary only mildly with radius. The DM axis ratios (p_DM, q_DM) are well recovered by the models with constant p_DM and q_DM, and the best-fitting models reproduce the stellar-halo-density distribution reasonably well;

• The DM halo of Auriga 12 is twisted, oblate (p_DM > q_DM) in the inner 20 kpc, and vertically aligned (p_DM < q_DM) in the outer regions. Our model with constant p_DM and q_DM recovers those in the outer regions, and it does not reproduce the stellar-halo-density distribution in the inner regions well. We recover p_DM(r), q_DM(r) of Auriga 12 as a function of radius generally well with a model that allows them to vary, which also better reproduces the stellar-halo-density distribution from inner to outer regions.

Acknowledgements

We thank David Hogg, Jianhui Lian, Zhaozhou Li for useful discussions. LZ acknowledges the support from the CAS Project for Young Scientists in Basic Research under grant No. YSBR-062 and the National Key R&D Program of China No. 2022YFF0503403. X.X.X. acknowledges the support from the National Key Research and Development Program of China No. 2024YFA1611902, CAS Project for Young Scientists in Basic Research grant No. YSBR-062 and YSBR-092, and the science research grants from the China Manned Space Project with No. CMS-CSST-2025-A11.

Appendix A The recovery of DM distributions for Auriga 25

Auriga 25 has a recent minor merger which significantly perturbed its stellar disk and inner halo. The infall orbit and mass of this merger are similar to those of the LMC in the Milky Way, but its perturbation to the disk and inner halo are significantly stronger than the LMC did to the Milky Way.

We orient Auriga 25 in the way that the LMC-like satellite is located in the southern hemisphere, and we take the particles in the north to create the mock data; the stars in the north look closer to dynamical equilibrium with less asymmetrical structures. The mock data are created in a similar way as for the other three galaxies. As shown in Figure A.1, there are significant substructures in the inner halo. We removed the big tail at the end of the disk, but not the other substructures.

We applied the dynamical model to this galaxy in the same way as we did to the three galaxies we showed in the paper. In Figure A.2, we show the recovery of the DM radial density distribution and the 3D shape of the DM halo. In this case, we have larger statistical errors in the results. However, the true DM distribution is still generally recovered by our model within the 1σ confidence level.

Fig. A.1 The surface brightness of Auriga 25 projected with different directions. Its disk was significantly destroyed by a recent merger. Both its stellar disk and its halo are not in fully dynamical equilibrium.

Fig. A.2 The recovery of DM radial density distribution and DM 3D shape of Auriga 25, similar to that shown in Figure 7 and Figure 8 for other galaxies.

References

All Figures

Fig. 1 Shapes of DM halos in Auriga 23, 5, and 12, measured in the coordinate system where the z-axis is perpendicular to the stellar disc, and X and Y are the long and short axes of the DM halo in the disc plane. The DM halos of Auriga 23 and 5 are oblate, aligned with the stellar disc (pDM > qDM), and vary little with radius. The DM shape of Auriga 12 varies from being oblate and aligned with the stellar disc (pDM ~ 1 and pDM > qDM) in the inner regions to being vertically aligned with the disc (pDM < qDM and qDM ~ 1) at r ≳ 20 kpc.

In the text

Fig. 2 Left: MW observations from LAMOST + Gaia. Right: mock data created for Auriga 23. Each dot represents one star/particle coloured by weight obtained from selection-function correction. There are about 15 000 k giants in the final sample of the smooth halo of the Milky Way. The mock data of Auriga 23 include about 20 000 particles, which generally match the spatial distribution of the MW sample.

In the text

Fig. 3 Mock data for Auriga 23. The blue dots are disc stars that we excluded, i.e. stars with λz > 0.8 or stars with λz > 0.5 and Z/Z⊙ > 0.5. The black dots are taken as halo stars and kept in our sample. We divided the data into 7 × 6 bins in rgc versus θ with the intervals of rgc = [5, 8, 12, 16, 22, 30, 40, 50] kpc, θ = [0, 10, 20, 30, 40, 55, 90°]; here, we only show three columns to illustrate the data.

In the text

Fig. 4 Comparison of stellar-density distribution from data and model for Auriga 23. The top panels from left to right are composed of observational data and several models with different DM axis ratios of qDM and pDM. The colours represent the number density as indicated by the colour bar. and the dashed curves are iso-density contours. The second row shows the uncertainty of the data, and model residuals. The data are well matched by the model with DM qDM = 0.64 and pDM = 1, which are well consistent with the ground truth. The stellar-density distribution constructed by models are either too round or too flat if the shape of underlying DM changes. The stellar-density distribution strongly constrains the 3D shape of the underlying DM halo.

In the text

Fig. 5 Comparison of velocity distributions from data and model for Auriga 23. The velocity distributions, υr, υϕ, and υθ are calculated in 7 × 6 bins in r versus θ, but we only show three columns here as labelled. In each panel, the grey solid and dashed curves are the velocity distribution and uncertainty of observational data. The red solid curves are from the best-fitting model. Blue and magenta dashed curves show models with different radial distribution of the underlying DM mass. The velocity distributions, especially υr, strongly constrain the underlying DM radial density distribution. The best-fitting model (red solid) reproduces the velocity distributions in all three components well: υr, υϕ, and υθ.

In the text

Fig. 6 Constraints on parameters of underlying DM distribution for Auriga 23. Each dot represents one model coloured by the relative $\Delta {\chi }_{\text{tot}}^2$ compared to the least χ2 model. The parameters are sampled on grid with fixed step sizes. We slightly dithered the positions of the points to make them visible. The solid red, black, and dashed purple contours are the 1σ confidence levels determined by $\Delta {\chi }_{\text{tot}}^2,\Delta {\chi }_{\text{vdis}}^2$ and $\Delta {\chi }_{\text{den}}^2$, respectively. The black horizontal and vertical dashed lines indicate the ground truth. The three parameters determining the radial distribution, p0, rs, γ, are degenerated. The shapes of the DM halo, qDM and pDM, are well constrained and consistent with the ground truth.

In the text

	Fig. 7 Recovery of DM radial-density distribution. The radial-density profile and enclosed mass profile we obtained for Auriga 23, Auriga 5, and Auriga 12, which are well consistent with the ground truth within 1σ uncertainty, are shown.
In the text

Fig. 8 Recovery of DM intrinsic shapes by models with constant qDM and pDM. The panels from left to right are Auriga 23, Auriga 5, and Auriga 12. In each panel, the solid red, black dashed, and purple dashed curves indicate the model 1σ confidence level determined by ${\chi }_{\text{tot}}^2,{\chi }_{\text{vdis}}^2$, and ${\chi }_{\text{den}}^2$, respectively. The dots are the ground truth measured from the DM halo of the simulation at different radii. The colours represent the radii where pDM and qDM are measured as indicated by the colour bar. The DM halos of Auriga 23 and Auriga 5 are oblate, with little variation from the inner to outer regions, which are well recovered by our model. The DM halo of Auriga 12 varies from oblate in the inner region to vertically aligned in the outer regions, the current models capture the shape of the DM at the outer regions.

In the text

Fig. 9 Reproduction of stellar halo flattening qstar as a function of radius Rgc for models with different underlying DM halos. The three panels from left to right are Auriga 23, 5, and 12. The DM halo shapes of Auriga 23 show little variation within our data coverage. The qstar(r) is well reproduced by a model with constant pDM and qDM close to the ground truth. Auriga 5 shows mild variation, and Auriga 12 shows relatively strong variation of the DM halo shape, from r ~ 10 kpc to r ~ 50 kpc. Especially for Auriga 12, the stellar-halo flattening, qstar(r), at all regions can hardly be reproduced by a model with constant pDM and qDM; in contrast, it can be reproduced well by a model with pDM(r) and qDM(r), which vary with radius, close to its ground truth. The stellar-halo flattening qstar(r) has strong constraints on the underlying DM shape.

In the text

	Fig. 10 Left: constraints on qin,DM, pin,DM, qout,DM, and pout,DM for model allowing DM axis ratios to vary as a function of radius for Auriga 12. Right: comparison of pDM(r) and qDM(r) profiles obtained by the model (red) and the ground truth (black). The solid and dashed red curves represent the mean and 1σ uncertainty of the results obtained by the model.
In the text

Fig. 11 Comparison of stellar-density distribution constructed from data and different models for Auriga 12. The columns from left to right are the data, best-fitting model with constant pDM and qDM, and best-fitting model with pDM(r) and qDM(r) varying as a function of radius, r. The model with constant pDM and qDM cannot match the density distribution in all radial regions; there are significant residuals in the inner 20 kpc. The consistency between data and model improves with the model allowing variable pDM(r) and qDM(r).

In the text

	Fig. A.1 The surface brightness of Auriga 25 projected with different directions. Its disk was significantly destroyed by a recent merger. Both its stellar disk and its halo are not in fully dynamical equilibrium.
In the text

	Fig. A.2 The recovery of DM radial density distribution and DM 3D shape of Auriga 25, similar to that shown in Figure 7 and Figure 8 for other galaxies.
In the text