{'fraction_non_alphanumeric': 0.03324349960101778, 'fraction_numerical': 0.018323070205814602, 'mean_word_length': 4.710110041265475, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 7, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Non-Parametric Spherical Jeans Mass Estimation with B-splines

Spherical Jeans modeling is widely used to estimate mass profiles of systems from star clusters to galactic stellar haloes to clusters of galaxies. It derives the cumulative mass profile, M(<r), from kinematics of tracers of the potential under the assumptions of spherical symmetry and dynamical equilibrium. We consider the application of Jeans modeling to mapping the dark matter distribution in the outer reaches of the Milky Way using field halo stars. We present a novel non-parametric routine for solving the spherical Jeans equation by fitting B-splines to the velocity and density profiles of halo stars. While most implementations assume parametric forms for these profiles, B-splines provide non-parametric fitting curves with analytical derivatives. Our routine recovers the mass profiles of equilibrium systems with flattened haloes or a stellar disc and bulge excellently (<~ 10% error at most radii). Tests with non-equilibrium, Milky Way-like galaxies from the Latte suite of FIRE-2 simulations perform quite well (<~ 15% error for r<~ 100 kpc). We also create observationally motivated datasets for the Latte suite by imposing selection functions and errors on phase space coordinates characteristic of Gaia and the DESI Milky Way Survey. The resulting imprecise and incomplete data require us to introduce an MCMC-based subroutine to obtain deconvolved density and velocity dispersion profiles from the tracer population. With these observational effects taken into account, the accuracy of the Jeans mass estimate remains at the level 20% or better.


INTRODUCTION
The European Space Agency's Gaia satellite (Perryman et al. 2001;Lindegren et al. 2016;Gaia Collaboration et al. 2016, launched in 2013, has released 3D positions and proper-motions for more than a billion Milky Way stars. In conjunction with lineof-sight velocities and chemical abundances for stars obtained with ground-based spectroscopic surveys (e.g. RAVE, Steinmetz et al. 2006;LAMOST, Hu & Jiang 2005; APOGEE, Allende Prieto et al. 2008;GALAH, De Silva et al. 2015; Gaia-ESO, Gilmore et al. 2012, DESI Allende Prieto et al. 2020,various components of the Milky Way's stellar halo are being characterized: field halo stars, stars in individual tidal streams, globular clusters, and satellites. The dynamics of these halo objects are important probes of the Milky Way's dark matter distribution and, in principle, allow us to determine several properties of the Milky Way's dark matter halo, most fundamentally its mass density profile. The mass of the Milky Way is a fundamental quantity of interest for comparisons with cos-of sophisticated techniques that incorporate the observational uncertainties via forward modeling. Starting with assumed parametric forms for the gravitational potential and distribution functions of tracers (star or globular clusters), these techniques impose observational errors and selection functions on the models to derive, often using Bayesian inference, the best fit estimates and their confidence intervals for a variety of potential parameters, including the mass of the halo. The most recent example of such distribution function modeling applied to Gaia (E)DR3 data (Deason et al. 2021) is based on a sample of 665 halo stars (98% with full 6D phase space data) in the distance range 50-100 kpc (including K giant stars, Blue Horizontal Branch (BHB) stars, RR Lyrae stars, and Blue Stragglers) and attempts to correct for the dynamical effect of the LMC.
We are on the verge of an era where the numbers of stars with full 6D phase space data is going to see a dramatic increase from hundreds to hundreds of thousands or even millions. The Dark Energy Spectroscopic Instrument (DESI) (DESI Collaboration et al. 2016a,b;Allende Prieto et al. 2020) is a 5000 fiber spectrograph on the Mayall 4-meter telescope at Kitt Peak National Observatory, Arizona, and is expected to obtain spectra of 5-8 million Gaia stars (down to = 19.5) in the next five years. In addition, the 4MOST (de Jong et al. 2012) and WEAVE (Dalton et al. 2014) surveys are similar multi-fiber spectrographic surveys which will together obtain spectra for millions of stars. These spectroscopic surveys will deliver line-of-sight velocities and spectro-photometric distances to stars out to 100 kpc, enabling -for the first time -the assembly of samples of stars with 6D phase space coordinates that are orders of magnitude larger than the best samples currently available.
In this work we implement a new version of one of the oldest and conceptually simplest dynamical modeling tools -the spherical form of the Jeans equations (Jeans 1915;Binney 1980;Binney & Tremaine 2008) -under the assumption that very soon we will have samples of 10 4 − 10 5 halo stars out to ∼ 100 kpc with full 6D phase space coordinates. Our goal in this paper is to assess the effects of various factors on the derived cumulative mass profiles obtained with the spherical Jeans equation: (a) the underlying assumptions of the modeling approach: spherical symmetry and dynamical equilibrium; (b) realistic observational errors on proper motions, line-of-sight velocities and distances to halo stars and (c) the effects of limited survey volumes inherent to all surveys.
Most previous implementations of the spherical Jeans equation, applied to the Milky Way halo either bin the data (Kafle et al. 2018;Wang et al. 2018) and/or use analytic functions (e.g. power laws) to describe the density and velocity profiles and velocity anisotropy of tracer stars (Gnedin et al. 2010). Binning the data necessitates the computation of numerical derivatives which can be noisy (Kafle et al. 2018). For this reason, the use of analytic functions has been favored because they are relatively easy to implement and, in some cases, align well with theoretical predictions. For example, Gnedin et al. (2010) assign a power law relation to the tracer density profile of BHB stars beyond 25 kpc, which is broadly consistent with predictions for the "broken power-law" density profile of the stellar halo from cosmological simulations (Bullock & Johnston 2005;Johnston et al. 2008;Cooper et al. 2011). Other studies use parametric fits to the velocity distribution and velocity anisotropy profile. These approaches have their limitations since parametric fitting curves may not adequately represent the true velocity or density profiles. Many studies construct parametric fits from binned data, which is sensitive to how bins are chosen, although discretekinematic Jeans models that avoid binning have also been developed (Watkins et al. 2013). Since the velocity anisotropy profile, in particular, requires accurate determination of proper motions (which have only become available in large numbers in the past 3 years, thanks to Gaia), it has historically been common to assume either a constant value or a few functional forms (e.g. Battaglia et al. 2005), or to use functional forms motivated by cosmological simulations (Xue et al. 2008;Gnedin et al. 2010). Only in recent years has it become possible to use the Gaia proper motions, sometimes even without line-of-sight velocities, to estimate the velocity dispersion tensor within ∼ 10 kpc from the Sun and to constrain the Galactic potential from the axisymmetric Jeans equations (Wegg et al. 2019;Nitschai et al. 2020).
In the context of external dwarf spheroidal galaxies, sophisticated Bayesian methods have been developed to derive parametric fits to the tracer density, line-of-sight velocity distributions and velocity anisotropy profiles (e.g. G S , by Read & Steger 2017). Diakogiannis et al. (2017) have developed a hybrid approach and use B-splines to construct the radial velocity dispersion profile 2 , while using a parametric mass model. They apply their method to line-of-sight velocities in the Fornax dwarf spheroidal galaxy. There have been several recent studies to assess how well the spherical Jeans equation performs on dwarf spheroidal galaxies when the various underlying assumptions are broken (e.g. spherical symmetry, dynamical equilibrium). For example, Evslin & Del Popolo (2017) use tracer particles drawn from N-body simulations to evaluate the fourth-order spherical Jeans equation when applied to dwarf spheroidal galaxies. Genina et al. (2020) assess how G S fares on mock dwarf spheroidal galaxies drawn from cosmological ΛCDM and SIDM simulations. El-Badry et al. (2017) show that, using mock data sets from cosmological hydrodynamical simulations, dwarf galaxies with recent episodes of star formation and feedback are not in dynamical equilibrium and Jeans modeling results in an overestimate of the mass of the dwarf galaxy. We do not discuss these methods for modeling external dwarf galaxies any further since the nature of the data, both in sample size and phase space dimensionality (typically only sky-positions and line-of-sight velocities for < 5000 stars) are completely different from the situation considered in this work.
We present a novel B-spline based routine for performing spherical Jeans modeling for the Milky Way. Prior to applying it to real data, it is important to assess the modeling biases introduced by any spherical Jeans modeling code. This is pertinent because real potentials of Milky Way-like disc galaxies are not spherically symmetric nor are they in dynamical equilibrium. In addition, it is helpful to have a quantitative understanding the effect that breaking these assumptions has on the resulting mass profile. For this purpose we construct a diverse suite of mock datasets to validate our routine and quantify the errors in the derived cumulative mass estimate when the standard assumptions are broken. Although there have been several recent works focused on understanding the errors introduced by breaking various assumptions of the spherical Jeans equation by testing the method on simulations (e.g. El-Badry et al. 2017;Kafle et al. 2018;Wang et al. 2018), such tests are required to validate each code. Consequently, the mock datasets in our suite include smooth, self-consistent equilibrium distribution functions generated using A (Vasiliev 2019) -an efficient all-purpose galactic modeling package. We also generate mock stellar haloes from three galaxies of the Latte suite of FIRE-2 cosmological hydrodynamic zoom-in simulations. (Wetzel et al. 2016;Sanderson et al. 2020).
For the Latte galaxies, we also consider the effect of imposing survey selection functions and mock observational errors. For example, we exclude star particles outside the fiducial footprint of the Dark Energy Spectroscopic Instrument (DESI) (DESI Collabora-tion et al. 2016a,b;Allende Prieto et al. 2020). Simultaneously, we impose mock observational errors on proper-motions, line of sight velocities and distances to tracer stars. Since selection functions and measurement errors lead to poor Jeans mass estimates if not accounted for, we attempt to alleviate this by developing a Bayesian MCMC subroutine to deconvolve observational effects from the underlying data. We demonstrate the efficacy of this deconvolution process with our Latte-based mocks. Once DESI data are available, we will apply our Jeans routine to Gaia proper motions and distances and DESI line-of-sight velocities to construct a mass profile of the Milky Way out to ∼ 80 − 100 kpc.
The layout of this paper is as follows. Section 2 gives a brief overview of spherical Jeans modeling and our B-spline implementation in the case of perfect and complete data (Sec. 2.1) and imprecise, incomplete data (Sec. 2.2). Section 3 describes the various mock datasets we test the implementation on. Section 4 describes the results of these tests. We summarize and conclude in Section 5, and Appendix A gives an introduction to B-splines.

ANALYSIS METHODS
The Jeans equations are obtained from the collisionless Boltzmann equation by computing second order moments of the distribution function (Binney & Tremaine 2008). Under the assumption of spherical symmetry, they reduce to a single equation which can be used to derive the enclosed mass at any radius (< ), from which the full gravitational potential Φ can be obtained using Φ( ) = (< ) . The form of the spherical Jeans equation we use here is given by The velocity anisotropy parameter is defined as The terms 2 , 2 and 2 are the means of the squares of the radial, azimuthal, and polar velocities respectively; is the number density of tracers (stars, globular clusters, etc.); is the spherical Galactocentric radius; is the gravitational constant. < 0 indicates a tangentially biased velocity distribution, = 0 indicates velocity isotropy, and > 0 indicates a radially biased velocity distribution.
Note that in the above we use the full second moments 2 , 2 and 2 rather than velocity dispersions 2 , 2 and 2 , since the former are more suitable for B-spline evaluation and these quantities are more relevant from the dynamical standpoint. However, in the cases with imprecise and incomplete data (Sec. 2.2), we assume zero mean tangential streaming velocity ( = = 0) to reduce the number of free parameters, in which case 2 ... = 2 ... . The spherical Jeans equation can be applied in two complementary regimes: (1) solving a "dynamical inverse problem" where one measures the dispersion of sky-plane and line-of-sight velocities directly from observations (with errors) and then "inverts" it to derive the radial velocity dispersion and anisotropy profiles (or assume ), and then infer the potential or cumulative mass distribution (Eq. 1); (2) "forward modeling," on the other hand, assumes some parametric form for the gravitational potential and tracer distribution function and then computes the expected velocity dispersion and anisotropy profiles and convolves them with expected observational uncertainties and selection functions and compares with observations.
Although dynamical inverse modeling was the original formulation of the method, inverse problems in general have difficulty dealing properly with errors and noise from small samples. Furthermore, it has been shown that in the context of the Milky Way halo, inferring the anisotropy profile from the line-of-sight velocity distribution results in a biased profile at Galactocentric distances beyond ∼20 kpc (Hattori et al. 2017). Most current applications therefore follow the second route because it is easier to account for observational errors. In this work, we assess how well the "inverse modeling" route behaves when full 6D phase space coordinate information is available (i.e. does not need to be inferred from los alone), for samples that are orders of magnitude (>10 4 ) larger than have previously been available. This method gives the potential (or rather the cumulative mass) more directly, but also suffers from errors more directly.
We split the analysis workflow into two stages. First, we reconstruct the radial profiles of tracer density ( ) and components of second moment of velocity 2 , , (or the velocity dispersion 2 , , , depending on the method) from the observed distribution of tracers. Then we substitute these functions into the Jeans equation (Eq. 1) to infer the mass profile. The first stage could be used in a more general context, independently from the second stage, and is implemented in two variants, described in the following sections: neglecting the observational errors and incomplete spatial coverage, or accounting for them. In both variants, these functions are represented as B-splines -a category of piecewise polynomials defined by an array of grid points (knots) and amplitudes of basis functions.

The case of precise and complete measurements
The first approach is fairly simplistic and suitable only for situations in which the observational errors can be neglected and the spatial coverage is uniform. In this case, the logarithms of the density and the second moments of velocity 2 as functions of ln are determined using out-of-the-box penalized spline fitting and log-density estimation routines implemented in A . The mathematical formalism is outlined in Appendix A. Specifically, we perform penalized spline regression to fit B-splines to 2 , 2 , 2 and perform a penalized spline density estimate to calculate ln( ). An equally valid alternative is to construct an estimate of the kinetic energy density in each dimension ln( 2 ... /2) and obtain the dispersions by dividing the kinetic energy density by the mass density.
There is no universally optimal procedure for determining the parameters of the knots: their number, spacing, minimum and maximum bounds. Diakogiannis et al. (2017) use a sophisticated but more computationally intensive evolutionary modeling algorithm to optimize and adapt their choice of knots when applying the JEAnS code to dwarf spheroidal galaxies. In contrast, we empirically determine our knot configurations, placing significant emphasis on keeping their parameters (minimum and maximum bound, and count) consistent over the variety of datasets to guard against overfitting.
In these fits, we use logarithmically spaced knots in radius (equally spaced in ln ) because we find that this best captures the radial distribution of tracer particles. We adopt 6 knots for datasets described in Section 3.1 and 5 knots for the datasets described in Section 3.2 and 3.3. The minimum and maximum radial bound of the knots is similarly dataset dependent, but we have arranged them such that the knots are identical across all datasets presented in a given figure.
B-splines are an attractive choice for multiple reasons: chiefly, they provide analytical derivatives and non-parametric fits. The radial derivative terms in Eq. 1 on ln( ) and ln( 2 ) are computed trivially by utilizing the B-splines' analytical derivatives, which are continuous across grid knots for B-splines of degree 2 and higher (we use degree 3 B-splines, which are equivalent to the familiar cubic splines). Further, the B-spline's non-parametric nature means that they do not restrict the velocity or density profiles to a specific shape as power-laws -which are often used for Jeans modelingwould.

The case of imprecise and incomplete data: deconvolving observational effects from underlying distributions
The second approach is designed to deal with the more complicated case when the observational errors are significant and the spatial coverage is non-uniform. We do this by forward modelling of the intrinsic density and velocity dispersion profiles, accounting for observational effects before comparing them to the actual data. The distribution function of tracers is written in a factorized form, where ( ) is the density and N | , 2 is the Gaussian velocity distribution for each of the three orthogonal velocity components in spherical coordinates (one radial and two tangential, the latter sharing the same velocity dispersion profile). As in the previous section, the logarithm of the density and the two components of the velocity dispersion , are represented by B-splines in ln with fixed knots and free amplitudes varied during the fit. For simplicity, we assume zero mean streaming velocity ( = = 0), but in principle one could model it as another spline function of ln with a few more free parameters. To account for the limited spatial coverage of the observed dataset, we normalize ( ) to have unit integral over the survey footprint and over the range of heliocentric distances corresponding to the allowed range of apparent magnitudes (including the uncertainty in the distance modulus, as described in Section 4.3.1 of Hattori et al. 2021). This normalization needs to be computed with a relative error 1/ stars in order not to affect the likelihood, and we optimize the procedure for the specific case of the DESI survey footprint, which covers two spatially disjoint regions on the sky. The 3D integration across the sky coordinates and distance modulus is performed separately in each region with a suitable variable substitution transforming it into a unit cube.
The likelihood of the model against the observed dataset is given by convolving the distribution function (Eq. 3) with observational uncertainties of each datapoint. Hattori et al. (2021) performed this convolution with the Monte Carlo method, splitting each star from the catalogue into samples ∼ 100 points drawn from the joint error distribution of all relevant coordinates (distance, proper motion and line-of-sight velocity) and converting them to Galactocentric position and velocity; these samples are kept fixed in the course of MCMC analysis to minimize the impact of the Poisson noise on the likelihood. In their case, the model parameters specify both the potential and the distribution function of stars. To evaluate the likelihood of each model, the position and velocity of all samples are converted to actions computed in the model potential, and then the action-space distribution function of the model is averaged over all samples belonging to each star. In our case, the general approach is similar, but with two important differences. First, the likelihood of a model is determined by the phase-space (rather than action-space) distribution function (Eq. 3), and the corresponding potential of each model is computed from the Jeans equation at a later step. Second, we benefit from the fact that the error convolution can be performed analytically for all three velocity dimensions thanks to the Gaussian form of the distribution function, and use the Monte Carlo integration only to account for the distance uncertainty. Specifically, we sample samples = 20 points from each star's distance error distribution and convert the 3D coordinates of these samples into the Galactocentric radii , where = 1 . . . stars enumerates stars and = 1 . . . samples iterates over samples for each star. These radii are kept fixed throughout the entire MCMC run. We also precompute auxiliary matrices for converting the two Galactocentric velocity dispersions , ( ) into the dispersions of three heliocentric kinematic quantities -two proper motion components ; , ; and line-of-sight velocity los; . Since the mean velocity in the model is assumed to be zero, the mean values of these heliocentric quantities are determined entirely by the Solar velocity and, in the case of proper motion, the distance to each sample, hence also precomputed in advance. Then for each choice of model parameters, we construct the splines for , ( ), evaluate them at each sample's radius, convert into the Heliocentric coordinates, and add in quadrature the measurement uncertainties of each star los; , ; , ; to obtain error-broadened dispersions of all three kinematic measurements (different for each sample). The likelihoods of the density, line-of-sight velocity, and both proper motion components are multiplied, averaged over all samples for each star, and the logarithms of these values are summed up to obtain the likelihood of the entire dataset: . (4) We first identify the best-fit model by maximizing the likelihood function, using the Nelder-Mead deterministic optimization algorithm. After this, we explore the parameter space with the MCMC search, using the emcee package (Foreman-Mackey et al. 2013). This gives the joint distribution of model parameters, and for each model in the chain, we compute the cumulative mass profile from the Jeans equation (Eq. 1), as before, with the only difference that the second velocity moments are replaced by velocity dispersions.
We stress that even if we allowed for a nonzero mean velocity and used the full second velocity moment 2 = 2 + 2 in the equations, this approach would not be equivalent to the one described in the previous section in the limit of complete sky coverage and negligible errors. The difference is subtle but might be important in some cases (e.g., a strongly anisotropic system): whereas in the first approach, we compute the kinetic energy directly, without assuming anything about the velocity distribution function, in the second approach we explicitly model it as a Gaussian function. In any case, we need to assume a certain functional form for the velocity distribution in order to perform the error convolution, and a Gaussian is by far the simplest and still fairly realistic case.

MOCK DATA
We now describe the various mock datasets used to validate our B-spline-based routine.
In Section 3.1 we describe the equilibrium mock datasets that were generated using A . These come in two categories: "Haloalone" models (hereafter HA-models) and "Halo-Disc-Bulge" models (hereafter HDB-models). The purpose of these tests is to verify the accuracy of the routine and to assess the effects of realistic deviations from sphericity and, in the case of the HDB-models, the effect of contamination of halo tracers by disc particles.
In Section 3.2 we describe more realistic mock datasets generated from star particles drawn from three galaxies (m12f, m12i, and m12m) of the Latte suite of FIRE-2 cosmological hydrodynamic zoom-in simulations (Wetzel et al. 2016). These models introduce significantly more complexity and realism, including the effects of halo substructure like streams and clustered halo stars, and disequilibrium arising from tidal interactions with nearby satellites. They also present the challenge of accurately separating disc and halo stars.
In Section 3.3 we describe how we impose observational selection functions and errors onto the mocks described in 3.2. The selection functions we consider are for the fiducial DESI footprint and survey magnitude limits of the Gaia and DESI surveys. The mock errors are imposed on proper motions, line-of-sight velocities, and distance moduli. These mocks also introduce a dependence on the choice of solar position, which we vary between three locations in the disc.

Self-consistent equilibrium mock datasets
The A dynamics package is used to generate self-consistent equilibrium distribution functions with all particles having the same mass. In practice, the halo tracers provided to the Jeans routine from these datasets are dark matter particles, but we assume in these simplified models that halo stars would be similarly distributed.
The HA-models contain only a spheroidal halo with axis ratio varying from = 1.0 (spherical) to = 0.6. A uses the double-power-law distribution function in action space from Posti et al. (2015) in which we tailor the mixing coefficients to produce the desired oblate axisymmetric shape. The mixing coefficients in the linear combination of actions in the distribution function are responsible for both the flattening and the velocity anisotropy, and it may not be straightforward to construct a model that follows a particular spatial density profile precisely. We determine the parameters of the distribution function (power-law indices, mixing coefficients, etc.) that approximately reproduce the Hernquist (1990) density profile. We then tune the mixing coefficients to create a desired amount of flattening, which we quantify by the axis ratio of the inertia tensor, as described in Zemp et al. (2011, their method E1). The models are fully specified by the distribution function, and the corresponding potential-density pair is determined via an iterative procedure (Binney 2014). Once generated, the halo mass and scale radius are scaled to halo = 1.3 × 10 12 and halo = 25 kpc, roughly matching the Galactic virial mass and halo density break radius quoted in Bland-Hawthorn & Gerhard (2016).
The HDB-models have a spherical halo, a disc, and a spherical bulge. The HDB's initial halo density function adds an exponential cutoff to the HA density function but scales to the same mass halo = 1.3 × 10 12 with the same scale radius halo = 25 kpc. The double-exponential disc is a combination of the thin and thick disc parameters quoted in Bland-Hawthorn & Gerhard (2016). Summing the thin and thick disc's stellar mass yields disc = 4.1 × 10 10 and a mass weighted average of the radial scalelengths and scaleheights for the thin and thick discs yields disc = 2.5 kpc and ℎ disc ≈ 375 pc. The spherical bulge follows the Spheroid density profile with parameters taken from McMillan (2017): total mass bulge = 8.9 × 10 9 , scale radius bulge = 0.075 kpc, outer cutoff radius cut = 2.1 kpc, and powerlaw indices ( , , ) = (1, 1.8, 0).
While we create variations of the HA-models with different axis ratios, the halo in the HDB-model remains spherical. However, we make additional HDB-models with other alterations. In one variant, we inject a random sample containing 1/4 of the system's disc particles into the halo tracer population to test the routine's sensitivity to contamination of the halo tracer population with disc particles. In another mock, we assess the effect of a radially varying velocity distribution for halo particles by creating a halo particle distribution with the Cuddeford-Osipkov-Merritt anisotropy profile (Osipkov 1979;Merritt 1985a,b;Cuddeford 1991), using ( = 0) = 0.2 and anisotropy radius = 100 kpc. The HA-models and HDB-models are generated with 300 000 particles in their haloes. The HDB-model with disc contamination has an additional 40 000 particles used as input, which come from the 160 000 particles comprising its disc.

Mock datasets from cosmological hydrodynamic simulations
The Latte suite of FIRE-2 cosmological simulations are performed using mesh-free hydrodynamics from the GIZMO code (Hopkins 2015) and the Feedback in Realistic Environments model 2 (FIRE-2) (Hopkins et al. 2014;Wetzel et al. 2016). The Latte galaxies are a collection of simulated Milky Way-mass galaxies at uniquely high resolution that are well suited for generating mock Milky Way-like datasets. We use the redshift = 0 public snapshots of the m12f, m12i, and m12m galaxies released with the Ananke synthetic Gaia catalog (Sanderson et al. 2020). We create datasets for each galaxy by selecting only the star particles with metallicities [M/H] < −1.5 and ages age > 8 Gyr. Note that we do not impose a geometric cut to exclude disc stars since this complicates the estimation of the density profile, we assume that our metallicity and age selection largely suffices to exclude disc stars. This is a proxy for selecting metal-poor accreted halo stars since the public FIRE snapshots do not contain information about whether stars were accreted from satellites or formed in situ in the host galaxy. Since information about whether a star was accreted is also not available a priori in real data, this provides an additional level of realism to our mock datasets. We calculate the [M/H] ratio and the hydrogen mass fraction ( ) from the mass fractions of He ( ) and metals ( ) provided in the public Latte data as follows: These mocks and the mocks in Section 3.3 operate under the simplification considering each Latte star particle to be a single tracer star. We consider our tracers to be RR-Lyrae stars (RRLs) because they serve as standard candles. Each metal poor halo star particle (∼7070 ) can be considered to represent the entire stellar population of approximately a single star cluster. It is generally considered reasonable to assume that on average a star cluster of this mass would contain one RRL. 1 After selecting tracers by metallicity and age, there are approximately 200 000 remaining particles in each dataset (m12f : 187 457, m12i: 229 347, m12m: 217 286).
We note in passing that we chose not to use mock stars from the Ananke mock Gaia catalogs (Sanderson et al. 2020) of the same Latte galaxies for two reasons. First, although Ananke produces stellar populations and phase space distributions of stars with uncertainties, the currently public catalogs are only for Gaia DR2 uncertainties, while we wish to generate mock catalogs for Gaia (E)DR3 and DR5. Second, the process employed to create the Ananke catalogs "explodes" each star particle from the simulation into a stellar population. While the stellar population particles preserve the local phase space density around the original star particle, this population is isotropic in velocity space, i.e. it does not preserving the local velocity anisotropy in the vicinity of the original star particle. Since the cumulative mass estimate (Eq. 1 & Eq. 2) depends explicitly on the velocity anisotropy of the stellar tracers, the mock stars from Ananke, with their isotropic velocity distribution would bias our results.
The three Latte galaxies, while all Milky Way mass, have varying amounts of halo substructure in the form of tidal streams and shells and a fairly broad range of satellite properties. Their haloes are also in different states of dynamical equilibrium, have a range of realistic morphologies, have experienced a varied assembly history (Garrison-Kimmel et al. 2018), and have realistic satellite populations (Wetzel et al. 2016). In addition, these galaxies have halo stars with moderate amounts of tangential streaming motions (azimuthal and polar ). Figure 1 shows heatmaps of star particles in these three galaxies from the present-day (redshift zero) snapshot, color coded by | | (upper row) and (lower row). All three galaxies, show significant streaming motions and substructure but m12m has the most dramatic streaming motion, so we expect our Jeans estimation to perform the worst on the mock datasets derived from it.
In addition to substructure in the form of streams, shells and over-densities, there is recent observational evidence that the inner part of the Milky Way, including the Sun, is moving with respect to the outer stellar halo as a result of on-going tidal interaction with the Large Magellanic Cloud (LMC) (Petersen & Peñarrubia 2021;Erkal et al. 2021) and that the stellar halo density distribution shows a distinct polarization pattern also likely caused by the LMC (Conroy et al. 2021). Erkal et al. (2020) and Deason et al. (2021) demonstrated that these perturbations lead to an overestimate of the Milky Way mass profile with traditional dynamical modelling methods relying on the equilibrium assumption, and Correa Magnus & Vasiliev (2022) presented a method for compensating these perturbations in mass modelling. Although the main effect of the LMC is kinematical, by coincidence, the magnitude of Milky Way mass overestimate roughly matches the mass of the LMC itself, if we were to add it to the virial mass of our Galaxy, as noted in Li et al. (2020). In this paper we do not consider the effects of significant perturbations, such as that arising from tidal interaction with the LMC.

Adding observationally motivated selection functions and errors
One must choose a solar position or a Local Standard of Rest (LSR) before imposing observationally motivated selection functions (SFs) or errors on mock halo stars from simulated galaxies. As part of the construction of the Ananke synthetic Gaia catalog (Sanderson et al. 2020), the Latte collaboration provided three LSRs with precomputed solar velocities for each of the three galaxies m12f, m12i, and m12m. The first LSR position, LSR0, places the sun at ( , , ) = (0, 8.2, 0) kpc and the other two, LSR1 and LSR2, are located at the same Galactocentric radius, in the galactic midplane, and rotated 120 • in either direction from LSR0. We transform the coordinates of the star particles to the equatorial and Galactic coordinate systems with A when applying the following selection functions and errors.
We start by imposing a sky-selection function similar to the fiducial DESI footprint with a declination cut −35 • ≤ ≤ 90 • and a galactic latitude cut 30 • ≤ | | ≤ 90 • . We also impose a magnitude selection function in the Gaia -band (16.0 < < 20.7), the lower bound determined by the DESI Milky Way Survey (MWS) (Allende Prieto et al. 2020) and the upper bound by Gaia's expected performance 2 . The tracer star particles from the Latte galaxies represent metal-poor RR Lyrae stars (RRL), for which the distance modulus D can be computed from the apparent magnitude , assuming the absolute -band magnitude = 0.58 with an intrinsic scatter 0.24 mag, corresponding to about a 10% distance error (Iorio & Belokurov 2019). We compute the distance modulus for our mock stars as follows: D = 5 * log( /kpc) + 10 (7) = + D The mock errors are randomly generated from independent Gaussian distributions on each quantity for each star particle with standard deviations determined by characteristic or expected errors for Gaia and DESI observations. We use a very conservative line-of-sight velocity error for the DESI Milky Way Survey los = 10 km s −1 (Allende Prieto et al. 2020), which is reasonable for RRLyrae which are pulsating stars, although the los errors are expected to be 2-5 km s −1 for other tracer populations such as luminous red giants. The Gaia proper motion errors ( * , ) come from PyGaia 3 . PyGaia allows for calculations of errors from Gaia (E)DR3, and expected errors from DR4 (containing the first 5 1 / 2 years of observations) and DR5 (containing up to 10 years of observations). The uncertainty in distance modulus D is the same for all particles and is caused by the intrinsic scatter in the RRL period-luminosity relation. We add distance modulus errors first, because the proper motion errors are magnitude dependent. Lastly, we assume the right ascension ( * ) and the declination ( ) have no error. Figure 2 shows the 3D velocities squared of tracers in m12f after imposing the listed observational selection functions and errors at LSR0 over three Gaia data releases. Because the Gaia proper motion errors increase with increasing , tracers at large Galactocentric radii tend to have their velocities dramatically boosted, often above the expected escape velocity curve. These extremely high velocity stars will cause enclosed mass estimate from Jeans modeling to blow up at large radii if not corrected for by our more advanced modeling described in Section 2.2. This effect makes the datasets described in this section practically unusable without deconvolution especially when generated with proper motion errors corresponding to earlier Gaia releases: (E)DR3 & DR4. The forward modeling method described in Sec. 2.2 also allows us to correct for the volume removed by the DESI footprint when calculating density.
The mock datasets described here are LSR-dependent, so imposing these SFs and errors produces three datasets for each Latte galaxy -one for each of the three LSRs -and have roughly 20 000-30 000 remaining tracers each. Further, for each Gaia release -(E)DR3, DR4, DR5 -we make a separate mock dataset, totaling to 27 possible mocks. We only show results for (E)DR3 and DR5 as the DR4 results lie between the other two as expected.

RESULTS
Each of the mock datasets described in Section 3 is provided as input to our routine, the B-spline functions for ln ( ), 2 , , ( ) are generated, and the resulting (< ) is computed using Eq. 1. This estimated mass profile is compared with the system's true mass profile over the relevant radial range. The true mass profile is computed for the entire radial extent of the system by cumulatively Mass profiles of three "Haloalone" (HA) models with flattening quantified by the axis ratio = 1.0 (spherical), 0.8, and 0.6 (oblate axisymmetric). Increasing flattening causes minor fluctuations of no more than 4% from the true mass at most radii. The larger error at small radii is most likely due to small numbers of tracers. These suggest that sphericity is not absolutely required for very good estimations with spherical Jeans modeling.
summing the masses of all particles (stars, gas, and dark matter where relevant) within a given Galactocentric radius.
The estimated mass profile from each mock is judged for accuracy by computing the percent error relative to the corresponding true mass profile. Figures 3-5 and 7 show the estimated mass profiles (solid lines), their percentage errors (lower panels), and the true mass profiles (dashed lines) as a function of Galactocentric radius for a collection of mock datasets. The mocks described in each subsection of Section 3 correspond to a matching subsection in Section 4.

Initial tests with self-consistent equilibrium mocks
In this section we describe the results from tests with the selfconsistent equilibrium mock datasets described in Section 3.1. Figure 3 shows the results for the three HA-models, whose B-spline fits are each run with knots = 6 knots logarithmically spaced between min = 1 kpc and max = 70 kpc. Here, we also remove particles beyond the bounds of the knots when performing the spline fits. The blip in error within < 5 kpc is merely an artifact of small number statistics at small radii. Aside from this, the progression of flattening from = 1.0 to = 0.6 only has a very small effect on the accuracy of the estimation. The mass profile from the = 0.6 mock oscillates around the true profile and is noticeably less stable than that of the = 1.0 and = 0.8 mocks. For all three mocks however, the error is less than 4% beyond = 5 kpc. The very small error in the estimations on the = 0.8 and = 0.6 mocks shows that breaking the sphericity assumption does not necessarily mean the Jeans estimate will be poor. We also run tests on HA models with = 0.9 and = 0.7, which perform in-line with the results shown here.
In MW-like galaxies, the sphericity assumption is also broken by the presence of a disc. Figure 4 shows the results of tests with the original HDB-model and its variants. For these models, the Bspline fits are all performed with knots = 6 knots logarithmically spaced between min = 5 kpc and max = 80 kpc. The log-scaled   Fig. 3, dashed curves show true distributions and solid curves show our Jeans estimates for mass profiles of the "Halo-Disc-Bulge" (HDB) models including variants with contamination from disc particles and a radially varying Cuddeford-Osipkov-Merrit velocity anisotropy profile. All variants perform very well: error ≤10% at most radii. The inclusion of a disc breaks the assumption of sphericity but has little effect on the accuracy, even when the tracer population is contaminated with disc particles (yellow curve). A radially varying anisotropy profile also imparts very little error beyond = 30 kpc where the disc ends. spline functions in log-radius are linearly extrapolated beyond the grid endpoints, corresponding to a power-law extrapolated density and kinetic energy profiles. The error in the estimation on the HDBmodels is within 10% beyond = 10 kpc. Notably, the significant deviations between the estimation on the original HDB model and its variants are confined to roughly < 30 kpc. Beyond that, the three mass profiles are quite similar. The differences at small radii between the disc contamination variant and the others are expected as the injection of rapidly orbiting and spatially flattened disc particles are significant changes to that region of the system.

Early tests with metal-poor particles from cosmological simulations
In this section we describe tests of our Jeans modeling routine with mock halo stars drawn from galaxies belonging to the Latte suite of FIRE-2 cosmological simulations: m12f, m12i, and m12m. The datasets we generate using these three galaxies select a population of metal-poor, old star particles with [M/H] < −1.5 and age > 8 Gyr meant to represent a stellar halo RR Lyrae population. The B-spline fits on these datasets are each run with knots = 5 knots logarithmically spaced between min = 5 kpc and max = 80 kpc. We reduce the number of knots from the tests with equilibrium mocks (Section 4.1) to increase the smoothness in the B-spline fits. Figure 5 shows the results of the tests for these three galaxies when all metal poor stars are selected. The accuracy with which our Jeans modeling routine recovers the input mass distribution is, of course, worse for these mocks than for the idealized equilibrium mocks. The resulting mass profiles all have very similar shapes with inflection points in the same radial regions. They are all within 15% error nearly everywhere and often much less. Like the HDB-models, they uniformly exhibit underestimates at roughly > 40 kpc in either set. Here, we only show estimates out to = 100 kpc for ease of comparison with the observationally motivated results in the next section (the DESI and Gaia magnitude limits will yield too few stars and large velocity errors beyond this distance). Nonethe- -10% 0 +10% % Error Figure 5. Mass profiles of mock datasets generated using three galaxies from the Latte FIRE-2 cosmological simulations: m12f (purple), m12i (orange), and m12m (teal). For each, we emulate a stellar halo RR Lyrae population by selecting star particles with [M/H] < −1.5 and age > 8 Gyr. All three resulting mass profiles are within 15% of the truth beyond = 5 kpc. These tests with advanced cosmological simulations illustrate our B-spline routine's ability to smooth out disequilibrium like tidal streams, halo substructure, and tangential streaming motion to recover the underlying mass distribution quite well. less, we maintain a very similar level of accuracy for our cumulative mass estimate ( 15% error) when we use error-free Latte data with ( min , max , knots ) = (5 kpc, 100 kpc, 5) evaluated out to = 200 kpc. We also obtain a similar level of accuracy ( 20% error) when using variants of these mocks with different metallicity thresholds ([M/H] < -1.0, -2.0) but the same age threshold.
In the following section, it is important to keep in mind that the three Latte galaxies used here have triaxial dark matter haloes, show significant halo substructure (shells and tidal streams), show varying amounts of tangential streaming (see of Fig. 1), and have satellites within 100 kpc. Despite these various sources of disequilibrium and spherical asymmetry, Figure 5 shows that, with error-free and complete data, the cumulative mass profiles are still well recovered. Further, we observe that the introduction of these sources of disequilibrium are more detrimental to the accuracy than the introduction of spherical asymmetries on our equilibrium models shown in Sec. 4.1. The distribution function fitting performed on mock stellar halo catalogues in Wang et al. (2015) made the same assumptions that we do here, spherical symmetry and dynamical equilibrium, and similarly found that errors from disequilibrium are dominant over those from asymmetry.

Tests with observational errors and selection functions on cosmological simulations
All of our previous tests have assumed that the data available to us are error-free, full 6D phase space coordinates of halo tracers throughout the volume of interest. This simplification enabled us to assess the systematic error introduced by the code itself, but did it not allows us to assess how realistic observational selection functions and errors affect the Jeans mass estimates. In reality, practically all observations come with errors and spatial selection functions.
In this section, we describe tests on the observationally motivated mock datasets described in Section 3.3 and the more advanced analysis we apply to them described in Section 2.2. These mocks better represent the data we expect to obtain from Gaia and DESI by imposing their observational selection functions and uncertainties in phase space coordinates of tracers. We now model the density and velocity distributions in the halo by deconvolving the observational effects of incomplete and imprecise data from the true underlying data. Figure 6 shows example output of our deconvolution routine on the three Latte galaxies in question with mock Gaia (E)DR3 and DESI errors imposed at LSR0. The MCMC fits (solid colored lines) recover the true density and velocity dispersion profiles (solid black lines) despite the input data (dashed colored lines) deviating from them significantly. In the right column the proper motion errors on the data clearly blow-up at large radii for the observed tangential velocity dispersion profile (dashed blue line), but the effect is entirely corrected for by deconvolution in all three galaxies (solid blue line). Properly accounting for the selection function also largely mitigates the effect of the incomplete spatial coverage of the survey on the tracer density profile, as seen most clearly as a sharp drop-off in density (dashed red curve) at < 10 kpc in the left column. The accuracy of the corrected profile is not uniform, however. While m12f has its density well recovered by the deconvolution routine, m12i and m12m still have large errors within = 10 kpc. We attribute this to their elevated amounts of substructure and streaming motion as seen in Figure 1. We also note that the true density profile and error-imposed density (red dashed curve) of m12m line up nearly exactly beyond = 80 kpc. This is undoubtedly the result of an enormous stream, visible in Fig.1, dominating the halo density at that radius. We discuss m12m in light of this stream and its relevance to our future MW application later in this section.
We can now plug the deconvolved density and velocity profiles into the spherical Jeans equation (Eq. 1) to calculate the enclosed mass profile (< ). Figure 7 shows results of the deconvolution run on Latte galaxies m12f, m12i, and m12m with mock Gaia DR5 and (E)DR3 errors. As described in Section 3.3, we consider three LSRs for each galaxy whose respective mass profiles are distinguished as separate solid colored lines. The colored uncertainty bands above and below the estimated mass profiles are 16-84 percentile intervals of models in the MCMC series.
From Figure 7, there are a couple key takeaways. In the DR5 panel (left), the small radii region ( < 20 kpc) has large errors and dramatic swings between over-and underestimating the enclosed mass across all three galaxies. This is not particularly surprising as much of this region had all the star particles removed by the observational selection functions, so we expect low confidence in this region. Beyond about = 20 kpc, the mass profiles for m12f and m12i are stable and yield impressive accuracy. It is only at about = 85 kpc when some of the m12f mass profiles reach more than 20% error while all the m12i mass profiles stay within 20% out to = 100 kpc. The results are quite different, however, for m12m. After the dramatic swings within 20 kpc, the m12m mass profiles continue exhibiting oscillatory motion around their true mass profile. Such poor performance is expected considering the tremendous streaming motion in m12m (which we ignored in our deconvolution algorithm) in comparison to m12f and m12i (Fig. 1). We note that even with perfect data over the entire volume m12m showed the largest deviation from the true mass distribution (see Fig. 5), but since the formulation of the Jeans mass estimate with perfect data shown in Fig. 5 did not assume , , = 0, the results were better than the deconvolved estimates which assume no streaming motion. In future applications, we may need to add more flexibility in the models by allowing a nonzero streaming motion around a particular axis (not necessarily perpendicular to the disc plane). It must be noted, however, that m12m should be considered a worst-case-  -20% 0 +20% m12m Figure 7. Results of Jeans estimation on Latte galaxies m12f (purple), m12i (orange), and m12m (teal) after deconvolving observational effects from underlying data. Mocks with Gaia DR5 proper motion errors are shown in the left column and (E)DR3 proper motion errors in the right column. The bands above and below each Jeans mass estimate represent the 16-84 percentile range of the models in the MCMC series. The three lower panels show the errors on the Jeans mass estimates for the three galaxies. The accuracy at = 100 kpc is usually within 20% for either Gaia release, but the uncertainty bands and the spread between LSRs are notably larger for the (E)DR3 mocks than for the DR5 mocks scenario. When we progress to real data we will not encounter many of the reasons that m12m produced a poor fit. There is not nearly as much streaming motion in the Milky Way as there is in m12m, and the well defined streams that are present, such as Sagittarius, are often easily removed from observational datasets. For these reasons, we expect results using MW data to be more akin to that from m12f and m12i rather than m12m.
We emphasize here that while non-monotonic cumulative mass profiles are completely unphysical they arise from the fact that the enclosed mass estimate at a given radius depends only on the tracer density profile and velocity distribution at that radius. We did not impose any relation between density and velocity dispersion profiles in the model from the outset to preclude non-monotonic mass profiles. It is straightforward to add this constraint in the fitting procedure (disqualifying non-monotonic models from the MCMC chain), but it is unlikely to mitigate the underlying physical reasons for the poor fit.
In the (E)DR3 panel (right), the results are remarkably similar to the DR5 panel (left) despite the much greater effect of larger proper motion errors on pre-deconvolution data (Fig. 2). The m12f results are nearly identical across (E)DR3 and DR5 runs. The m12i results are also nearly identical until about = 50 kpc, beyond where the spread between different LSRs and the uncertainty bands are much greater for the (E)DR3 than the DR5 runs. The (E)DR3 m12m results have the same issues as those of DR5 but, similarly to m12i, with increased spread between LSRs and increased uncertainty band size. The similarity in accuracy is a testament to the efficacy of the deconvolution routine when dealing with large or small proper motion errors.

SUMMARY & CONCLUSION
We have presented a new non-parametric spherical Jeans modeling code that uses B-splines to fit the enclosed mass distribution within the Milky Way halo. We present tests of this code with a variety of mock datasets to assess the effects of breaking various assumptions made in deriving the spherical Jeans equations. Our main results are listed below: • With our simplest "Halo-alone" (HA) models we show that going from a spherical mass distribution to a flattened mass distribution with axis ratio = 0.6 only slightly increases the error on the estimate of enclosed mass but the error is within 4% beyond = 5 kpc (see Fig. 3).
• Our "Halo-Disc-bulge" (HDB) models are constructed to resemble the Milky Way halo, disc, and bulge but consider only a spherical halo. We also include the additional test of considering both isotropic and anisotropic velocity distributions for halo particles, and the test of having "contamination" from disc stars in the halo. Our Jeans modeling code performs very well with these HDB models with maximum errors on the enclosed mass almost always within 10% and less than 5% at most radii (see Fig. 4).
• The mock datasets with the next level of complexity are constructed from the Latte FIRE-2 zoom-in cosmological hydrodynamical simulations of three galaxies, m12f, m12i, and m12m. These galaxies were considered Milky Way analogs and contain realistic amounts of substructure and display varying amounts of disequilibrium (see Fig. 1). With error-free data over the whole sky, the code still gives better than 15% error out to 100 kpc in all three Latte galaxies (see Fig. 5).
• The most realistic datasets build upon the previous by considering the effects of observational selection functions and errors from Gaia and the DESI Milky Way Survey (Allende Prieto et al. 2020). The Gaia component includes a -band magnitude upperlimit, proper motion errors from the end of (E)DR3 and DR5, and errors on spectro-photometric distances. The DESI component includes its observational footprint (declination and galactic latitude bounds) and errors on line-of-sight velocity. These errors and selection functions are imposed on the same three Latte (FIRE-2) galaxies and consider three possible solar positions. Paired with the addition of observational effects is the inclusion of a more advanced method for modeling the velocity and density profiles. This method is designed to deconvolve the effects resulting from the observational errors and selection functions causing the data to be imprecise and incomplete. The two primary observational effects are on the density at small Galactocentric radii, due to the removal of volume by the DESI footprint, and on the tangential velocity dispersion at large radii, due to Gaia proper motion errors scaling up with heliocentric distance. Figure 6 shows the 2 and profiles of the underlying incomplete and imprecise dataset as well as their corrected counterparts after applying deconvolution. Jeans estimates performed with this setup yield very encouraging results (see Fig. 7). When run with DR5 proper motion errors, the error on (< 100 kpc) is around 20% and often better. The (E)DR3 runs give very similar results for m12f but have greater spread between LSRs in m12i, and m12m. The resulting mass profile for m12m has dramatic oscillations above and below true mass profile and violates monotonicity for cumulative mass profiles. We expect this breakdown of our Jeans routine is due to the very large streaming motions present in m12m which is not accounted for in our deconvolution routine (see Fig. 1).
The spherical Jeans equation method is a fairly simple, but popular method for recovering the mass distribution of (nearly) spherical potentials. A fundamental feature of this method is that it depends on the radial derivatives of the tracer mass distribution as well as the radial velocity dispersion profile. Since derivatives of observational data profiles are extremely noisy (Kafle et al. 2018), most authors model the tracer density and velocity dispersion profiles with analytic functions (often power-laws) to enable convenient computation of the derivatives. In our implementation of this method, we have used B-splines to estimate the density profile and velocity dispersion profiles of the tracer particles, enabling us to avoid the need to assume a functional form but still allowing for the rapid computation of the enclosed mass profile 4 .
Other implementations of the spherical Jeans equation formalism applied to full 6D phase space data (e.g. Wang et al. 2018) have also been tested on several hundred isolated and binary haloes from the dark-matter only Millenium II simulations (Boylan-Kolchin et al. 2009) and 12 intermediate resolution cosmological hydrodynamical 4 The routine's run time without deconvolution is a few seconds. zoom-in simulations of Local Group analogues (the APOSTLE simulations, Fattahi et al. 2016;Sawala et al. 2016). Wang et al. (2018) find that their Jeans estimates of 200 and the halo concentration parameter 200 yielded errors of around 25% when dark matter particles are used as tracers but the errors can go up to 200-300% when using star particles as tracers (even with error-free phasespace data). The APOSTLE datasets used by these authors are most directly comparable to our error-free Latte datasets (Sec. 3.2), on which we achieve errors on (< 100 kpc) of ∼ 5−15% (Fig. 5). We believe this discrepancy is likely due to three factors. First, although not explicitly stated in their paper we infer that Wang et al. (2018) used binned data for the density ( ), velocity dispersion ( ) and anisotropy ( ) profiles which were likely quite noisy. Our early experiments with binned data showed that numerical derivatives of such data are extremely noisy and do not have the desirable qualities B-splines provide our implementation: smooth non-parametric fits and analytical derivatives. Second, Wang et al. (2018) cut out the inner 20 kpc to exclude the in situ stellar component but include tracers out to the virial radius. We instead selected old metal poor stars with [M/H] < −1.5 and age > 8 Gyr in an attempt to capture a more virialized halo population and placed our maximum knot at = 80 kpc, within the ∼ 100 kpc limit where we expect a significant sample of Gaia+DESI data. Our experiments selecting more metal rich stars do show a small increase in the estimated errors on (< 100 kpc) but only up to 20% (not 300%). Increasing the radius of our maximum knot to 100 kpc and evaluating our B-splines out to 200 kpc with extrapolation maintains our great accuracy at 15% error on (< 200 kpc). Finally, while (Wang et al. 2018) fit their derived potentials to NFW potentials and report the errors on the best-fit 200 , 200 we simply report the errors on the actual cumulative mass profiles of our simulated haloes, which we do not compare with NFW profiles. While other factors like the specific implementation of baryonic microphysics in the different simulations might cause some differences, the three galaxies from the Latte suite of FIRE-2 simulations we have chosen in this study have global properties quite similar to the Milky Way. Our tests in this article suggest that, for real Gaia+DESI data, the prospects for obtaining errors of no more than 25% are quite favorable and would be a significant improvement over the current factor of 2 uncertainty in halo mass.
Other methods such as distribution function fitting (Watkins et al. 2010;Eadie et al. 2017Eadie et al. , 2018Eadie & Jurić 2019;Hattori et al. 2021;Deason et al. 2021;Shen et al. 2021) use Bayesian MCMC analysis to properly account for the errors on the mass estimates arising from the error on the data. These methods, while much more sophisticated, are computationally and technically challenging, and for simplicity, most authors assumed very simple power-law profiles for the tracer density and the potential. Our goal in developing this non-parametric B-spline Jeans modeling code and testing it with mock data is to provide the community with a fast and easy-to-use code with well quantified systematic errors. We dub the code (Non-parametrIc jeans Modeling with B-spLinEs) and have made it publicly available on GitHub (https://github.com/nabeelre/ nimble). Included in the repository are the routines for performing the analysis described in Sec. 2.1 and Sec. 2.2 as well as code to generate all the mock datasets described in Sec. 3. Generally, it can be applied to a sample of observed 6D coordinates for tracer stars to construct B-spline estimates of the radial density profile, velocity dispersion profiles and the velocity anisotropy profile, which are then used to evaluate Eq. 1.
In future, we aim to apply this code to Gaia DR4 and Gaia DR5 data with line-of-sight velocities from DESI, as well as future data expected from WEAVE (Dalton et al. 2014) and 4MOST (de Jong et al. 2012).

DATA AVAILABILITY
The data in this article come from two primary sources.
The data described in Sec. 3.1 are generated using the A package (Vasiliev 2019) using scripts provided in the NIMBLE repository, which can be accessed at https://github.com/ nabeelre/nimble.
The data used to make the mocks described in Sec. 3.2 and 3.3 are of the Latte suite of FIRE-2 cosmological zoom-in simulations (Wetzel et al. 2016)

APPENDIX A: FITTING DENSITY AND VELOCITY PROFILES WITH B-SPLINES
In this section we provide a brief overview of the B-spline formalism in the context of fitting and density estimation. For a more complete mathematical background see Appendix A2 in Vasiliev (2018).
B-splines are a set of piecewise-polynomial basis functions ( ) defined by knots of a grid , = 1.. . B-splines of polynomial degree are nonzero on at most + 1 consecutive grid segments, their − 1'th derivatives are continuous at grid knots, and a grid with knots ( − 1 segments) generates + − 1 basis functions. A B-spline representation of a curve is a linear combination of basis functions with some coefficients: ( ) = ( ). The most familiar case is a histogram, which is none other than a 0th degree B-splines basis; however, since even the curve itself is discontinuous, it is not suitable for evaluating derivatives. In this study we use 3rd degree B-splines, which are mathematically equivalent to a clamped cubic spline curve (clamped meaning that its 2nd derivatives at endpoints are not necessarily zero).
In the present context, we use B-splines to estimate the density from discrete samples. Namely, given a set of points with weights , = 1 . . . data , we seek to construct a density function ( ) such that ∫ ( ) ≈ ( ; , ), where the indicator function ( ; , ) = 1 when ≤ ≤ and zero otherwise. The density is normalized so that the integral over the entire interval containing all points equals . It turns out that B-splines are very useful for representing the logarithm of , ensuring the positivity constraint on density. After choosing the interval and grid knots, the coefficients of the B-spline representation ln ( ) = ( ) can be found by maximizing the penalized log-likelihood of the input sample: { } = arg max ln ( ) − ∫ 2 ln ( )/ 2 2 , where the second term biases the solution towards smooth curves, and the roughness penalty is automatically determined by crossvalidation. With the optimal choice of , the number of grid knots has little effect on the resulting curve (as long as it is sufficient to resolve key features in the distribution). Note that since the density is normalized to unity, the number of free parameters in the fit is one fewer than the number of spline knots.
For the Jeans equation, we need the density and the kinetic energy, which are estimated as follows. The distribution of points in spherical radius (the number of points in the radial interval .. + ) is (< )/ = 4 2 ( ), and it is convenient to use ≡ ln instead of radius. We construct a B-spline estimate for the logarithm of the density of points in , using the masses of particles in the input snapshot as the weights associated with points: 1 ( ) ≡ ln / ln = ln 4 3 ( ) ; hence the mass density is given by ( ) = (4 ) −1 exp 1 ( ) − 3 . Likewise, when using the masses multiplied by squared radial and tangential velocity components 2 , 2 + 2 of each particle as the weights, we obtain the estimates of the corresponding terms in the kinetic energy tensor: 2 ( ) ≡ ln 4 3 ( ) 2 , etc. The logarithmic derivative of the radial kinetic energy ln 2 / ln appears in the Jeans equation directly, and the other terms contains the ratios of exp 1 , exp 2 , etc. This paper has been typeset from a T E X/L A T E X file prepared by the author.