Charting Galactic Accelerations with Stellar Streams and Machine Learning

We present a data-driven method for reconstructing the galactic acceleration field from phase-space measurements of stellar streams. Our approach is based on a flexible and differentiable fit to the stream in phase-space, enabling a direct estimate of the acceleration vector along the stream. Reconstruction of the local acceleration field can be applied independently to each of several streams, allowing us to sample the acceleration field due to the underlying galactic potential across a range of scales. Our approach is methodologically different from previous works, since a model for the gravitational potential does not need to be adopted beforehand. Instead, our flexible neural-network-based model treats the stream as a collection of orbits with a locally similar mixture of energies, rather than assuming that the stream delineates a single stellar orbit. Accordingly, our approach allows for distinct regions of the stream to have different mean energies, as is the case for real stellar streams. Once the acceleration vector is sampled along the stream, standard analytic models for the galactic potential can then be rapidly constrained. We find our method recovers the correct parameters for a ground-truth triaxial logarithmic halo potential when applied to simulated stellar streams. Alternatively, we demonstrate that a flexible potential can be constrained with a neural network, though standard multipole expansions can also be constrained. Our approach is applicable to simple and complicated gravitational potentials alike, and enables potential reconstruction from a fully data-driven standpoint using measurements of slowly phase-mixing tidal debris.


INTRODUCTION
Stellar streams are the remnants of tidal disruption, a phenomenon which occurs when an ensemble of stars becomes tidally stripped in the underlying galactic field. Stream progenitors range from satellite dwarf galaxies (e.g. Majewski et al. 2003;Malhan et al. 2021;Panithanpaisal et al. 2021) to globular clusters (e.g. Odenkirchen et al. 2001a;Yuan et al. 2020;Alabi et al. 2020). During disruption, stars lost from the progenitor follow galactocentric orbits similar to that of the progenitor itself. This results in long stream-like structures extending many degrees on the sky (see, e.g., Shipp et al. 2018), sometimes wrapping around the parent galaxy multiple times (Koposov et al. 2012;Belokurov et al. 2014;Ramos et al. 2020).

jnibauer@princeton.edu
Kinematically cold, metal-poor streams are particularly useful for studies of the galactic gravitational potential, as they provide a snapshot of an extended stellar orbit in phase-space. There are several independent methods for constraining the galactic potential: fitting orbits directly to stream measurements (e.g. Johnston et al. 1999;Fardal et al. 2006;Koposov et al. 2010;Varghese et al. 2011), particle ejection methods (Küpper et al. 2012;Fardal et al. 2015;Bonaca et al. 2014;Gibbons et al. 2014), action-angle tracks (Sanders 2014;Bovy 2014), clustering in integrals-of-motion space (Sanders & Binney 2013a;Sanderson et al. 2015;Reino et al. 2022), or, finally, direct N -body simulations (e.g. Dehnen et al. 2004;Law & Majewski 2010). While Nbody simulations provide the most accurate kinematic depiction of stellar streams, they are too costly to search the parameter space efficiently for even a semi-realistic galactic potential model. Alternatively, orbit fitting of streams is attractive in its relative computational effi-ciency, but while the track of a stream system is similar to the orbit of the progenitor, stellar streams do not generally delineate a perfect orbit (Sanders & Binney 2013b). Particle-release methods attempt to model the formation of stream systems by releasing particles in a trial orbit according to the (analytic) tidal radius of the progenitor system (Küpper et al. 2012;Gibbons et al. 2014;Fardal et al. 2015;Bowden et al. 2015). After a given integration time, the final generated stream can be compared to the observed system. This method is far less costly than direct N -body simulations, though it still scales poorly to complicated potentials that require many model parameters. Action-angle coordinates provide a natural basis for stream formation, and can be used to obtain the model parameters of a given potential that maximize clustering in action space (Sanders & Binney 2013a;Sanderson et al. 2015;Reino et al. 2022). Related methods include modeling the action-angle distribution of streams directly. For instance, Sanders (2014) constructs a generative model for stream formation in action-angle space. Conversion from actionangle space to phase (i.e., position-velocity) space relies upon the choice of the potential. Similarly, Bovy (2014) models the initial action-angle distribution of tidal debris, which can easily be integrated forward. In this framework, the orbit of the progenitor is estimated and the mean difference in orbital frequencies along the stream relative the progenitor is characterized. The mean stream track in action-angle space is converted to phase-space using an approximate linear interpolation method. Both of these methods require a prescription for stream formation in action-angle space, and assume a uniform distribution in stripping times for stars that become unbound from the progenitor. In general, analytic transformations from phase-space coordinates (x, v) to action-angle coordinates (θ, J ) are known for only a limited number of potentials for which the Hamilton-Jacobi equation is separable. Action-angle methods therefore require relatively simple models for the underlying potential, the most flexible of which is the two-component Stäckel model (Stäckel 1891), which assumes that all orbits are defined by the same foci. This constraint has been shown to be incorrect for real galaxies (Kuijken & Gilmore 1989;Binney 2012). The Stäckel-Fudge method (Binney 2012;Sanders & Binney 2015) enables greater flexibility, by treating the true potential as locally similar to a Stäckel potential for which action-angle coordinates can be derived. This approach is typically implemented in the previously-mentioned works that model the action-angle distribution of tidal streams. While useful, the Stäckel-Fudge approach is approximate and will suffer if the local potential is poorly described by a Stäckel model.
Regardless of the method and trade-offs therein, most studies using stellar streams to constrain the galactic potential rely on a narrow class of analytic models or truncated basis function expansions to represent the potential. However, the true morphology of the dark matter halo in the Galaxy could deviate from simplified models (see, e.g., Prada et al. 2019;Garavito-Camargo et al. 2021;Shao et al. 2021). Adopting an incorrect model for the galactic potential can lead to biased parameter constraints, regardless of the accuracy of the adopted stream-modeling approach (Bonaca et al. 2014). A method which does not place functional priors on the underlying galactic potential could therefore provide an important advancement in mapping the detailed distribution of dark matter in e.g. the Milky Way stellar halo. Such constraints could then be compared to cosmological simulations, for which functional priors on the dark matter distribution are not adopted (e.g. Diemand et al. 2008).
Recent work (Bonaca & Hogg 2018) has demonstrated that more flexible potential models (e.g. basis function expansions) fit to individual stellar streams probe the local properties of the potential in the neighborhood of the given stream, but are not necessarily representative of the global potential. Consequently, a population of streams constrains different aspects of the global potential. Bonaca & Hogg (2018) therefore highlight the need for a flexible method that can constrain the global potential from multiple streams simultaneously, while accurately reproducing the local features of the potential probed by individual systems. Motivated by a Fisherinformation analysis of simulated stellar streams, these authors suggest to analyze each stream independently using a flexible model, thereby maximizing the learned information content of the given stream. As a postprocessing step, a global potential can be constrained by interpolating the individualized stream fits across multiple systems.
In this work we consider a fully flexible and differentiable approach to directly estimate the galactic acceleration field in the neighborhood of a given stream. Our method does not require a model for the underlying galactic potential, nor does it require streams to delineate isoenergy curves in phase-space (i.e. stellar orbits). The output of our analysis is similar to that of Naik et al. (2022), since we constrain galactic accelerations rather than a latent potential model. We assume that the potential is mostly static, or changing slowly. Once the acceleration field in the neighborhood of a number of streams is measured, we demonstrate that a global,  Stars are color-coded by the phase angle γ, which increases monotonically along the stream from the trailing to the leading tail. The positions and velocities of stars along the stream are then fed to a neural network (C), which parametrizes a mean track, x θ (γ), through the position of stream stars along with a track speed. Both the stream track and track speed are parametrized in terms of the γ-parameter. These fits can be differentiated and combined through Eq. 6 to estimate the cartesian components of the acceleration vector along the stream directly. An intermediate model for the galactic potential never needs to be specified. During training of the neural network, we encourage the condition that the tangent vector to the stream track falls roughly parallel to the local trajectory of stream stars. We also encourage a smooth acceleration field, by penalizing neural-network parameters that give rise to large local changes in accelerations along the stream track. An example output is illustrated in the bottom right panel (D), where each cartesian component of the stellar stream is "unwrapped" along the γ-axis, and color-coded by the inferred accelerations.
flexible gravitational potential can be constrained. The resulting potential is consistent with the local properties of the acceleration field where each stream resides, and is fully data-driven. Alternatively, interpretable analytic models for the potential can be fit to the estimated accelerations along each stream to recover physical parameters. The paper is organized as follows. In §2 we introduce our method for estimating accelerations from stellar streams. We provide a probabilistic framework to connect our model to the data in §2.2, and discuss regularization of the inferred acceleration field in §2.3. In §3, we apply our method to simulated stellar stream data to demonstrate an accurate reconstruction of the accelera-tion field in a known ground truth potential. We compare our approach to fitting stellar streams with orbits in §4, and demonstrate accurate potential reconstruction in §5. We provide a discussion of our method in §6, along with limitations of this work and future directions.

METHOD
We now introduce our approach for estimating the galactic acceleration field from 6D phase-space measurements of stellar streams. Several streams in the Milky Way have 6D phase-space measurements available (e.g., Koposov et al. 2010;Antoja et al. 2020), though the majority have only partial coverage. As a first analysis, we utilize the full 6D phase-space distribution of stellar streams. We will consider the case of missing phase-space observations in a future work. We provide a brief summary of our approach below, and a more technical description in §2.1. Stellar streams are populated by stars on proximate but distinct orbits. Diversity among the orbits depends upon the potential of the progenitor, its history of mass loss, the mean potential of the galactic host, and perhaps on substructure within the host's halo. Absent the last of these, one expects the local mean of the orbits to vary slowly along each arm of the stream. Motivated by these considerations, in this work we treat stellar streams as a mixture of orbits, with neighboring segments of the stream populated by stars with a similar mixture. However, distinct segments of the stream may have different mixtures of orbits. Accordingly, we do not assume that streams delineate a single orbit in the correct galactic potential. Instead, we parametrize a path through a given stream in phase-space which characterizes the average local trajectory of the stellar orbital mixture. This path is flexible and differentiable, and does not necessarily coincide with any single orbit in the galactic potential. By applying the chain rule to eliminate explicit time-dependence, an acceleration vector can be estimated at each evaluation point along the stream from the parametrized path. An intermediate model for the galactic potential never needs to be specified.
Our analysis workflow is illustrated conceptually in Fig. 1. In the bottom left panel (A), we plot a projection of a mock stream in angular on-sky coordinates (φ 1 , φ 2 ). Arrows indicate velocity vectors. The track of the stream is determined by linking nearest neighbors in 6D phase-space. This induces an ordering along the stream, which we encode with a scalar parameter γ ∈ [−1, 1] (B). This enables us to "unwrap" a stream along the γ-axis, similar to unwrapping a particle's trajectory, (x(t), y(t)), in time. In practice, the parameter can be determined using an on-sky position angle and velocities to distinguish the leading arm from the trailing arm. We delay a more general discussion of the γ-parameter and its determination to §2.4.
Once the data are ordered, a neural network is used to parametrize a flexible and differentiable track along the observed stream ( Fig. 1; panel C). This neural network takes γ as an input, and outputs the three cartesian coordinates (x, y, z) and speed v of the track. We denote the position-space stream track with x θ (γ), where θ are the parameters of the neural network. The mean speed as a function of γ is v ϕ (γ), where ϕ are the parameters of an independent neural network that characterizes the track speed. In Fig. 1 we consolidate these neural networks down to one system, though the usage of two independent neural networks for the stream track and track speed is only an implementation detail.
We assume that locally, stream stars are characterized by a similar mixture of orbits at a slightly different orbital phase. In order to maximize the validity of this assumption, neural-network training is carried out so that the stream track is encouraged to point along the local mean motion of stars along the stream. This condition helps the neural network determine a path through the stream for which our assumptions are most completely satisfied. In addition, we adopt a prior on the neural-network parameters θ that encourages smooth estimates of the acceleration field as a function of the scalar parameter γ. This prior is implemented by penalizing neural-network parameters that give rise to large derivatives of estimated accelerations with respect to γ. Parameters of the stream-track neural network are adjusted so as to satisfy the tangent and smoothness conditions while also remaining compatible with the observed phase-space data.
An example output of our analysis is shown in the bottom right panel of Fig. 1 (D), where stream stars are "unwrapped" along the γ-axis and color-coded by their inferred accelerations in each cartesian component. In summary, our analysis utilizes a neural network to model the track of a stream in 6D phase-space. By applying the chain rule, derivatives of position and velocity along this track provide direct estimates of galactic accelerations. The required inputs are the positions and velocities of stream-members.

Technical Description
Any given stream consists of a collection of stars on somewhat different orbits. In the space of Integrals of Motion (IoM), stream tracers are clustered with locally similar IoM at slightly different phase angles (see, e.g., Sanders & Binney 2013a;Fardal et al. 2015;Reino et al. 2021). However, at a global level stream tracers span a range of energies. We demonstrate this property in Fig. 2, where we plot an N -body realization of a stellar stream generated in an axisymmetric potential after ∼ 3.5 Gyr of evolution (left panel). For an axisymmetric potential, the energy (E) and z-component of angular momentum (L z ) are IoM. We plot the N -body stream in the space of binned mean energy and binned mean L z in the right panel of Fig. 2, where these quantities are relative to the progenitor cluster. Points are colorcoded by the average angular coordinate φ 1 in the left panel, which increases monotonically along the stream. The progenitor cluster is shaded in gray. In the space of IoM, we see a distinctive "bow-tie" feature, with separate wings corresponding to the leading and trailing In the left panel, we plot a projection of the N -body stream in angular coordinates. In the right panel, the stream is plotted in the space of energy and the z-component of angular momentum, both of which are integrals of motion in a stationary axisymmetric potential. These quantities are plotted relative to the progenitor cluster, shaded in gray. Outside of the progenitor, tracers are color-coded by their position angle, φ1, in the left panel. This angle increases monotonically along the stream, demonstrating that stream stars have locally similar integrals of motion at slightly different phase angles. This figure is provided only for conceptual purposes; our method does not model the energy distribution of stellar streams directly, nor does it work in the space of integrals of motion.
arms of the tidal stream. The "bow-tie" is a result of the progenitor oscillating between pericenter and apocenter while losing stars to the galactic tidal field. Once stars become unbound from the progentior, they are "lockedin" to the potential of the galaxy and no longer oscillate with the progenitor. Further discussion is provided in Gibbons et al. (2014) and Yoon et al. (2011). From this figure, it is clear that tracers are-on average-sorted by their energy and angular momentum along the stream. That is, stream stars farthest from the progenitor typically have the largest offsets in E and L z from the progenitor, whereas unbound stars near the progenitor have smaller offsets in these quantities (see, e.g., Johnston et al. 2001;Gibbons et al. 2014).
For a given band in φ 1 , Fig. 2 illustrates that the stream is populated by stars each with slightly different orbits. However, the local mixture of orbits changes only gradually along φ 1 from one segment of the stream to the next. Consequently, while stellar streams do not trace a single orbit in the galactic potential, they typically consist of an ensemble of stars characterized by a locally similar mixture of orbits. At a global level, however, any two segments of the stream separated by a substantial phase angle can have larger differences in orbits.
Our analysis exploits this property of stellar streams: that is, we assume adjacent small segments of the stream are populated by an ensemble of stars with a similar distribution in the space of IoM. The properties of this distribution are allowed to evolve along the stream, albeit slowly. If two local segments of the stream have a similar distribution in the space of IoM, then the mean of the IoM distribution will not vary significantly with small changes in phase angle along the stream. Under this view, the mean path of a stream is characteristic of many small orbital segments which have locally similar IoM.
The local clustering of IoM in Fig. 2 (bottom right) as a function of position angle φ 1 motivates us to fit the mean dynamical properties of the stream as a function of some underlying phase. We refer to this fit as the stream track, which characterizes the local position and motion of stream tracers as a function of e.g. position angle φ 1 . The stream track forms a curve in 6-dimensional phasespace. We parametrize this curve in terms of the scalar parameter γ: where x and v are position and velocity vectors, respectively. In practice, we estimate x(γ) and v(γ) using a neural network with parameters θ. We discuss estimating these quantities from phase-space data of stellar streams in §2.2. However, in this section we suppress dependence on θ for simplicity. Typically, γ is taken to be a position angle on the sky, increasing from the trailing to leading arm of the stream (e.g., φ 1 in Koposov et al. 2010 andFig. 2, or Λ in Vasiliev et al. 2021).
Our central assumption is that stellar streams have locally similar mixtures of IoM at slightly different phase angles. That is, the mean instantaneous motion of a small stream segment will fall roughly parallel to the local direction of the stream track. Mathematically, this assumption implies the existence of a relationship between the parametrized stream track in position space, x(γ), and the parametrized track velocity, v(γ). In particular, we define the unit vector tangent to the stream track at scalar parameter γ as where x (γ) ≡ dx(γ)/dγ. If the velocity of an infinitesimal segment of the stream track follows the local morphology of the stream, this means where v(γ) is the local mean speed of stream tracers at scalar parameter γ. In practice, this relationship (Eq. 3) between x(γ) and v(γ) is enforced through our modeling, such that the best fit stream track is the one for which our assumptions are most completely satisfied given the observed data (see §2.2 for details). We emphasize that while Eq. 3 is valid for a stellar orbit in the galactic potential, it does not necessarily imply that the parametrized track (x(γ), v(γ)) is that of any particular orbit in the underlying potential. In fact, provided that the mean IoM of neighboring small patches of a stream are locally similar, the parametrized track as a whole is not required to delineate a path of constant energy in any potential whatsoever. This is because while stream tracers typically fall along similar orbits locally, distinct segments of a stream separated by substantial phase angles may have significantly different energies (as well as other IoMs) and therefore different orbits. We demonstrate this major advantage of our analysis in §4.
We next define the scalar path length s, which represents the path of the stream track with dimensions of length. Explicitly, this quantity can be expressed as a function of the scalar parameter γ as follows: where we have used x = (x, y, z) andγ is the integration variable. Following the discussion above, we identify x(γ) locally with an average particle trajectory x(t), the velocity v(γ) with dx/dt, and ds/dt with v(t) . Applying the chain rule, the local correspondence between γ and time is then With this correspondence, the local particle acceleration is obtained as In other words, whereγ is an evaluation point, and the term in parentheses is from Eq. 4. Consequently, given a differentiable parametrization of the stream track in phase-space, the acceleration vector along the track can be estimated using Eq. 6 without having to formulate a model for the potential. This enables a substantial increase in flexibility over previous methods, which require an analytic form for the potential to be adopted. In this way, our approach is agnostic to the functional form of the potential. Furthermore, Eq. 6 is inherently local and does not assume that distinct segments of a stream separated by a substantial phase fall along the same orbit. We demonstrate this quality in subsequent sections. With Eq. 6, stellar streams can be used to sample the galactic acceleration field directly. Samples of the galactic acceleration field can be utilized to infer the potential, Φ, through the relation

Fitting Stellar Streams
From a set of 6D phase-space observations (that is, positions and velocities of stars belonging to the stream) a number of analytical and data-driven methods can be employed to estimate a differentiable parametrized curve (x(γ), v(γ)) in phase-space. These range from polynomial fits, splines, to fully flexible neural network representations. In the subsequent sections we consider the latter of these methods, since a neural network representation allows us to flexibly describe the data while encouraging smooth estimates of the underlying acceleration field. In this section, we also assume that a suitable γ for each tracer has been assigned, though we discuss a more general determination of the γ parameter in §2.4. The set of measurements required for the present analysis is then D = {(x n , v n , γ n )} N n=1 , where x n and v n denote the position and velocity of the n th streammember, and γ n encodes where the stream-member belongs along the unwrapped stream. There are N total stream members.
Given the separability of the tangent vector T from the stream track speed v in Eq. 3, we fit the stream track x θ (γ) and speed v ϕ (γ) separately, where θ and ϕ are the parameters of the two independent neural network representations of these quantities. We use a standard class of neural networks called multilayer perceptron (MLP) with a fully connected feed-forward architecture. The neural network x θ maps a scalar input, γ, to a 3-dimensional position output (x, y, z). The network v ϕ (γ)-which we refer to as the speed neural network, or track speed-maps the same scalar input γ to a scalar output, representing the local speed of stream tracers. We discuss the specific architecture of these neural networks in Appendix B.
Fitting a MLP involves tuning the parameters of the neural network until an objective function is optimized, typically analogous to χ 2 fitting. For the case of the speed neural network, v ϕ (γ), we adopt a mean-squarederror loss function of the form where w n is a weight for each measurement that could be informed by estimated signal-to-noise of a given data point. The term ϕ 2 2 denotes the L 2 norm of the neural network parameters ϕ, and h is a hyper-parameter of the model. Practically, the L 2 norm works to enforce smoothness in the space of v ϕ , and reduce overfitting by penalizing particularly large parameter values. In the machine-learning literature, this regularization technique is commonly referred to as weight decay. We adopt h ∼ 10 −6 , where h is the weight decay parameter. The optimal parametersφ are determined using standard back-propagation routines with the Adam optimizer (Kingma & Ba 2014), until the loss function in Eq. 8 is minimized. All neural-network training in this analysis is performed in mini-batches, such that model parameters are only updated after iterating through several observations (see Masters & Luschi 2018 for a review on mini-batch training). This is standard practice when training neural networks, and reduces computational load while increasing the efficiency of training.
We next turn to our determination of the stream track, x θ (γ). Eq. 3 implies that the stream track is parallel to the mean velocity of the local stream stars. This follows from our assumption that a given segment of the stream has a similar mixture of IoM compared to an adjacent segment. We build this property into our modeling, in order to ensure that the best-fit stream track is the one that most completely satisfies our assumptions. We first define the quantity which represents the unit-vector tangent to the stream track at scalar parameter γ. We subscript this quantity with θ, since it depends on the parameters of the stream track neural network. The observed unit velocity vector for the n th stream star is T n ≡ v n / v n . Generally, the objective function can be thought of as a likelihood for the data D given a set of model parameters (e.g., θ). Indeed, Eq. 8 can be equivalently expressed as the combined negative log-likelihood of Gaussian distributed random variables, though presenting the objective function as a weighted mean-squared-error loss as we have done is more typical in the machine-learning literature. For the case of the stream track neural network x θ , we obtain an objective function using the more general Bayesian framework of likelihoods and priors, since we later incorporate a regularization prior in our analysis. We connect the data D to the model parameters θ through a multivariate Gaussian distribution. The likelihood for the n th stream star is where N (x|µ, C) is a multivariate Gaussian distribution with mean vector µ and covariance matrix C, evaluated at x. The symbols I 3×3 denote 3 × 3 identity matrices. We treat τ 1 and τ 2 as constant hyper-parameters of the model, which act as weights on the real space position and trajectory of the parametrized track, respectively. Similar to the weights w n in Eq. 8, the likelihood in Eq. 10 provides a pathway to incorporate errors on the phase-space measurements used in our fits through the covariance term, analogous to the likelihood in Bonaca et al. (2014);Hogg (2018). We postpone the inclusion of errors to a future work.
Assuming that the measurement of position and velocity for each star is statistically independent, we may write the likelihood for the data as a product over all stars: The logarithm of this quantity can be taken as an objective function to be maximized over the model parameters θ. We derive the explicit objective function used during model training in §2.3, where we include the addition of regularization prior which encourages smooth estimates of the acceleration field.
While fitting for both the position and trajectory of the stream track simultaneously in Eq. 11 might appear redundant, it has the beneficial feature of self-regulation. That is, we fit for a parametrized curve in phase-space with the property that the tangent vector at each point along the curve is parallel to the local instantaneous motion of stream tracers in 3D space. Therefore, the best fit stream track is not required to pass through the the centroid of the stream (i.e., the central region perpendicular to the elongation of the stream) as one would find when fitting the 3D positions of stream stars alone. Instead, the best fit stream track is the one which will satisfy our assumptions most completely. We illustrate this aspect of our analysis in §3.

Regularization of Inferred Accelerations
Poisson's equation, ∇ 2 x Φ = 4πGρ, connects the Laplacian of the gravitational potential to the underlying matter density, ρ. From this relation, in order for the acceleration field to be physical it must have the property ∇ x · a ≤ 0. Furthermore, because the force of gravity is conservative we also expect ∇ x × a = 0; this ensures the existence of a scalar potential Φ such that a = −∇ x Φ. We have experimented with incorporating similar constraints as physically motivated priors on the neural network parameters that characterize the stream track and track speed. However, the method presented in this paper works by estimating changes in position and velocity along the 1-dimensional track of a stream, which we relate to the underlying galactic acceleration field. Our estimate of the acceleration field, while 3dimensional, is also confined to this track. As a result, we infer what is effectively a "slice" through the full acceleration field, parametrized by x θ (γ). Therefore, from our method alone and without further assumptions on the functional form of the potential, we only have access to derivatives of the acceleration field along the stream track. We cannot in general calculate ∇ x ·a at each evaluation point along the stream track. With regards to enforcing a physical mass density ρ ≥ 0 along the stream track, in Appendix C we show that this condition can always be satisfied independent of the shape, location, or velocity of the track. This is again a consequence of our analysis providing a slice through the acceleration field, rather than a full spatial map. However, estimates of the acceleration field along one or more streams can be combined to constrain a flexible potential. From this potential, one can compute ∇ x · a and other expressions involving full-spatial derivatives. We demonstrate this aspect of our analysis in §5.
As an alternative and more feasible regularization, we adopt a prior on the neural network parameters that penalizes noisy estimates of the acceleration field. That is, we down-weight neural network parameters for which da θ /dγ is large. The acceleration vector inferred from Eq. 6 depends on two independent neural networks: namely, the stream track x θ (γ), and the track speed v ϕ (γ). While it is possible to train both of these neural networks jointly so that da θ /dγ is typically small, we have found this to be unnecessary and computationally expensive. Instead, we fit the speed neural network v ϕ (γ) separately, and fix its model parameters to the optimal valuesφ. We then train x θ (γ), penalizing large derivatives of the inferred accelerations along the stream with ϕ set to the optimal values,φ. We find that the stream track will sometimes induce noise artifacts in the inferred accelerations, presumably because it is a more complicated neural network which parametrizes a curve in 3-dimensional space. We find that such artifacts are uncommon for the simpler speed neural network, v ϕ , provided that we regularize the model parameters of this network with a weight decay term (Eq. 8).
In the context of Bayesian inference, we promote smoothness of the acceleration field by adopting an improper exponential prior over the stream track parameter θ, given the optimal track speed parametersφ as follows: where τ 3 > 0 is a hyper-parameter of the model, and a θ|φ is the acceleration vector along the stream inferred with θ free to vary and ϕ fixed to the optimal values, ϕ. Applying Bayes' theorem, the posterior over the parameters θ is proportional to the product between the likelihood and the prior, where the likelihood is from Eq. 11. We utilize Eq. 13 as an objective function over θ. In particular, up to a constant, we minimize the negative logarithm of the posterior density over θ, where we have defined λ 1 ≡ τ −1 1 , λ 2 ≡ τ −1 2 , and λ 3 ≡ τ −1 3 . Minimization of Eq. 14 can be performed efficiently using standard back-propagation routines. In this work, we use the Adam optimizer implemented in Py-Torch (Kingma & Ba 2014;. In practice, we set λ 1 = 1 and initialize training with λ 2 = 1 and λ 3 = 0. After ∼ 50 initial training epochs, we set λ 2 to a larger value, such that the order of magnitude of the vector norm λ 2 T n − T θ (γ n ) 2 is at least the same as that of x n − x θ (γ n ) 2 . This ensures that the first two terms in the loss function of Eq. 14 receive similar weight. After an additional ∼ 50 training epochs, we then choose λ 3 such that λ 3 da θ|φ /dγ 2 has roughly the same order of magnitude compared to the first penalty term on average. Our analysis relies heavily on the tangent vector to the stream track, T θ . This can be seen from our expression for the acceleration vector along the stream track, Eq. 6, and our expression for the track velocity, Eq. 3. While Eq. 6 depends on dx/dγ explicitly, there is an additional implicit dependence on this quantity through dv/dγ and Eq. 3. Meanwhile, the position of the track in Eq. 6 is less relevant, since only its derivatives are involved when estimating accelerations. However, the position of the track is important when we need to assign a location to our estimates of the acceleration field. Because of this, we ultimately give a larger weight to the λ 2 , tangent vector term in Eq. 14 due to its relative importance in estimating accelerations. However, the position of the stream track cannot be neglected, and must still remain compatible with the data. To strike a balance, we typically choose a final value of λ 2 that is a few times the order of magnitude of the mean λ 1 term in Eq. 14. This ensures that we find a path which characterizes the mean trajectory of the stream, while also remaining compatible with the measured positions of stars along the stream. A more optimal algorithm might maximize λ 2 in Eq. 14 while simultaneously minimizing θ through the other terms in the expression. For the present analysis, however, providing a slight advantage to the λ 2 term in Eq. 14 is found to provide robust estimates of the acceleration field while remaining compatible with the data. A simple diagnostic can be performed by visually inspecting the stream track x θ (γ) compared to the position of stars along the stream. In addition, the unit vector trajectories of the stars can also be compared to the track tangent vector, T θ (γ). A successful fit will characterize the stars in both of these sub-spaces, as illustrated in Fig. 3. In this way hyperparameter tuning is not arbitrary, since the data can be used to inform proper hyper-parameter values. Indeed, a visual inspection of the stream fits is used in this analysis when tuning hyper-parameter values.
We find that typically, λ 2 will be set to relatively large values, ∼ 10 2 − 10 3 , partly due to the difference in units between the λ 1 and λ 2 terms in Eq. 14 and due to the relative importance assigned to the λ 2 term. For the third penalty term in Eq. 14, we find λ 3 ∼ 10 is typically sufficient to ensure a smooth acceleration field devoid of high frequency noise artifacts from the stream track fit. We also adopt a weight decay value of 10 −3 for the stream track neural network. This limits potential over-fitting. Hyper-parameter tuning can be automated using hyper parameter optimization libraries such as RayTune implemented in PyTorch (Liaw et al. 2018) to yield the best fit stream track. However, we find our analysis is robust to hyper-parameter values, provided that the individual terms in Eq. 14 receive roughly similar non-zero weight overall, with λ 2 set to somewhat larger values relative to the other two terms. Neural networks are trained for roughly 3000 epochs in mini-batches consisting of ∼ 50 − 100 tracers.

Determining γ
Our modeling introduced in §2.1 depends on a scalar γ, which acts as a phase parameter that increases monotonically along the stream. For the majority of streams discovered using e.g. Gaia so far, an on-sky position angle can be used to inform γ since these systems do not subtend large angles, nor do they wrap back onto themselves ). However, for more complicated streams like Sagittarius, an on-sky position angle is not sufficient since the stream wraps around the Milky Way several times (e.g. Koposov et al. 2012;Belokurov et al. 2014;Antoja et al. 2020;Ramos et al. 2020). In the context of simulated data, Gibbons et al. (2014) devised a method to "unwrap" a given stream using the known time evolution of each tracer. This method is not applicable to real data, for which we only have access to a kinematic snapshot of any given stream. To devise an alternative approach, we first note that there is a "gauge freedom" in the choice of γ, in the sense that no physically significant quantity is affected by the replacement γ → f (γ) provided the function f (·) is smooth and monotonic. This can be readily seen if we consider γ to be an on-sky position angle, e.g., φ 1 in Fig. 2. Each tracer has its own unique angular coordinate φ 1 , which can be mapped to its 6d phase-space position. An equally valid mapping can be achieved in a frame which is scaled and rotated relative to Fig. 2, with new angular coordinatesφ 1 ,φ 2 . Consequently, γ simply assigns an ordering to stream tracer particles, with the property that a continuous change in the parameter corresponds to a continuous change in position along the stream track.  . We illustrate this parametrization with the colored curve, color-coded by the γ value along the stream. Bottom Row: For the same mock stream, we plot the unit-vector trajectories of tracers in the x, y, and z direction (black scatter points). The red curve is our parametrization of these points, taken as the derivative of the colorful curve in the top panel with respect to γ, and normalized. The parametrized curves in this figure are generated from a single neural network. Streams are unwrapped using the method described in Appendix A.
Because of this gauge invariance, once stream tracers have been ordered, we have the freedom to assign any order-preserving scalar γ to the stream particles. For complicated stream morphologies with one or several loops, we use a nearest-neighbors algorithm and an auto-encoder neural network to assign an ordering to stream stars. We provide a technical discussion of this approach in Appendix A, and we briefly summarize the method below.
Given an unordered set of 6D phase-space measurements of stream tracers, we first construct a graph in phase-space which connects nearby particles. We impose the condition that the ray connecting neighboring segments in position space is roughly aligned with the local trajectory of stream particles. We refer to this method as nearest neighbors with momentum, since the algorithm moves along the stream following the local motion of tracers. This algorithm connects adjacent particles in phase-space, allowing us to assign an ordering to particles. Due to the momentum condition, the algorithm inevitably passes over some stream particles without incorporating them into the nearest neighbors graph. To fill in these gaps, we require a flexible interpolation function which maps a tracer particle in 6d phase-space to an ordered scalar parameter γ. We adopt an auto-encoder neural network to assign a suitable ordering to tracers that were left out of the graph. This neural network maps a phase-space coordinate (x i , v i ) to a scalar γ i , with the arbitrary condition that the output of the auto-encoder, γ i , is in the interval [−1, 1]. This method is found to successfully assign a scalar parameter γ to a variety of stream morphologies. While this approach is not strictly necessary for simple stream morphologies, it provides an automated means to "unwrap" a stellar stream and measure derivatives along its track. We emphasize that the unwrapping of a stream occurs as a pre-analysis step: in order to measure derivatives along the track of a stream, we must first assign an ordering to stars along the stream. For the remainder of this work, all streams are unwrapped using the automated method and the full 6D phase-space measurements for the given stream.
We compared our nearest-neighbor-with-momentum technique to classical dimensionality-reduction methods such as kernel PCA (Schölkopf et al. 1999) with a radialbasis-function kernel, and found that traditional techniques could not successfully unwrap stellar streams.
In this section we generate multiple stream systems in a ground-truth gravitational potential, and fit these streams directly using the methods introduced in §2 to sample the underlying galactic acceleration field.

Generating Mock Stellar Streams
We focus primarily on dynamically cold streams with globular-cluster progenitors. The disruption of globular clusters has been studied extensively in direct N -body simulations (see, e.g., Heggie & Hut 2003;Küpper et al. 2010), which have shown that clusters tend to experience mass loss driven by the tidal forces of the host galaxy and two-body relaxation. Stars which escape the progenitor's tidal radius have slightly different energies from that of the progenitor, and therefore somewhat different orbits in general (Eyre & Binney 2011;Sanders & Binney 2013b).
We generate mock streams in a given potential using the "particle-spray" method from Fardal et al. (2015) implemented in the Gala package (Price-Whelan 2017a). This method is not a direct N −body simulation of globular cluster disruption, though it provides a prescription to reproduce streams generated in such simulations efficiently. In brief, particles are released from the progenitor near its Lagrange points for each step of an orbit integration. Released particles are then integrated forward in the background potential of the galaxy, and that of the (evaporating) progenitor system. This prescription has been shown to reproduce several detailed features of streams generated in realistic N −body simulations (e.g., stream fanning and morphology; Fardal et al. 2015;Pearson et al. 2015), and has been utilized extensively in other studies of stellar streams (e.g., Valluri et al. 2021;El-Falou & Webb 2022;Qian et al. 2022).
For the galactic potential, we adopt a triaxial logarithmic potential of the form We choose this potential not because it is necessarily representative of the actual Milky Way halo. Rather, if no two of the q i are equal, the potential will not be well matched to commonly-used axisymmetric models. We adopt q 1 = 1, q 2 = 1.3, q 3 = 0.9, R h = 2 (the lengths R h , x, y, z being measured in kpc) and v c = 150 km/s. We model the progenitor system as a Plummer-sphere potential with a mass of 2.5 × 10 4 M and a Plummer radius of 4 pc. We initialize the progenitor system with a randomly sampled position in a 3-dimensional box with a side length of 30 kpc. Initial velocities are sampled such that the initial speed is less than or equal to the local circular velocity. We generate 6 streams using this method, giving rise to a diverse range of tidal debris. The progenitor's orbit is integrated forward for 6 Gyr in timesteps of dt = 1 Myr, with two particles being released near the progenitor's Lagrange points at each time-step. An example of one generated stream is illustrated in the top panel of Fig. 6, in a rotated angular coordinate frame. Stream particles are colored by their energy relative to that of the progenitor.
Our method assumes that the local elongation of a given stream is tangent to the local motion of its tracers. For systems where the progenitor globular cluster is still intact, particles released from opposing lagrange points lead to sharp morphological features for which our assumptions are challenged (for example, the progenitor cluster has not fully disrupted in Fig. 2; φ 1 ≈ −35.5 deg). We therefore remove the progenitor from our analysis on simulated data to avoid fitting this feature. For real stellar streams, only a few systems have progenitors that can still be observed (see, e.g., Pal-5, Odenkirchen et al. 2001b;Erkal et al. 2017;Pal-13, Shipp et al. 2020;Globular cluster NGC 5466, Belokurov et al. 2006;Grillmair &Johnson 2006, andIbata et al. 2021). The progenitors of the vast majority of streams discovered so far appear to have dissolved Shipp et al. 2019;Riley & Strigari 2020;Bonaca et al. 2021). For systems with a prominent progenitor, we anticipate removing this segment of the stream to avoid bias. However, for the majority of observed streams we do not expect the progenitor to have substantial remnants. After removing progenitors, we also remove ∼ 1.5% of tracers from the extended tidal tails of a given simulated stream, since it is unlikely that such extended tracers would be deemed stream members with statistical confidence. Furthermore, the quality of our derivative estimation will suffer in the limit of having only a few ( 10) tracers to infer dynamical changes along the stream.
These cuts are fiducial, and we find that our results are robust to small variations around these analysis choices. The typical number of stream tracers is ∼ 7000 after these cuts. However, our method is applicable to streams with a significantly smaller number of measured tracers, and can interpolate over missing information. We demonstrate this advantage of our method in §4, where we estimate the acceleration field along an intermediate segment of a given stream without measuring its tracers.

Sampling the Galactic Acceleration Field
We now apply our method introduced in §2 to six stellar streams generated in a triaxial logarithmic potential. We utilize the full 6d phase-space coordinates of stream tracers here, though we comment on the applicability of our approach to streams with missing phase-space dimensions in §6.3.
In Fig. 3, we illustrate our fit to the stream track, x θ (γ), for a single mock stream. The fitting procedure is described in §2.2. The top row is a galactocentric view of stream tracers (black points), with the fitted track overplotted and color-coded by the γ value for each segment of the stream. In the bottom row we plot the cartesian components of the tangent vector to the stream track, Eq. 9, in red. The black points in the bottom panels correspond to the cartesian unit-vector trajectories of individual stream tracers. Our method fits the stream track in both the space of 3d positions (top row) and 3d unit-vector trajectories (bottom row), such that the best fit track is not necessarily required to pass through the local centroid of the stream. This behavior can be seen clearly in the top row, leftmost panel around (x, y) ≈ (0, −6), and the top row middle panel around (y, z) ≈ (5, −4); in these regions, the best fit stream track is offset from the highest density of points in position space. This behavior is due to our fitting procedure, since we fit for both the position-space component of the track and its trajectory simultaneously in Eq. 14. Consequently, the best fit path is not required to pass through the centroid of the stream; rather, the best fit track is the one that simultaneously fits the stream in the space of tracer positions (top row) and trajectories (bottom row). This fitting approach has the benefit of self-regulation, since the best fit track is the path for which our assumptions are most completely satisfied for a given stream (see §2.2 for further discussion).
In order to infer a local acceleration vector along a stream, Eq. 6 requires a differentiable estimate of the local speed of stream tracers as a function of position along the system. As discussed in §2.2, we fit for the speed independently. In Fig. 4, we plot the speed of stream tracers in black as a function of the scalar parameter γ, and our fit to these points, v ϕ , in red. Fig. 4 corresponds to the same stream plotted in Fig. 3. Combined, Figs. 3 and 4 demonstrate that a basic, fullyconnected neural network is able to capture non-trivial bulk motions of stellar tracers in a given stream, while maintaining smoothness and differentiability. Furthermore, we note that training a neural network to generate Figs. 3, 4 is relatively fast and can be performed locally on a personal computer.
Once a differentiable parametrization for each stream is estimated, calculating the acceleration vector, a(γ), along the stream path is achieved by simply differentiating the neural networks x θ (γ) and v ϕ (γ) with respect to γ through Eq. 6. We emphasize that this process is performed independently for each stream. The global acceleration field can be stitched together by fitting a gravitational potential that reproduces the local acceleration field in the vicinity of each stream. In §5, we combine localized constraints to infer the global galactic potential. A static potential is not strictly required for the method presented in this work, provided that the analyzed stream is populated by stars with a locally similar distribution of orbits distributed over slightly different phase angles along the track of the stream. This is not the case for streams with a significant velocity component perpendicular to the track, which can occur in time-dependent potentials (Erkal et al. 2019). We provide a more detailed discussion of time-dependence in §6.2.
Constraints on the galactic acceleration field along six mock stellar streams are illustrated in Fig. 5. Each column corresponds to a different mock-stream generated in the triaxial logarithmic potential, Eq. 15. From top to bottom row, we plot the three cartesian acceleration components a x , a y , and a z . The black points depict the ground-truth acceleration, while the red points depict the inferred accelerations along the stream using our method. The missing gaps correspond to the removed progenitor (discussed in § 3.1). Across the six stellar streams in Fig. 5, we find a generally strong agreement between the reconstructed acceleration components and the ground-truth acceleration field. The typical fractional error is at the few to sub-percent level, demon- strating that our method generally provides an unbiased estimate of the local acceleration field along a given stream. In §5, we demonstrate that our constraints can reproduce the true parameters of the potential without bias (Fig. 8).
Typically, fractional errors are largest near the ends of the streams; this is because near the boundaries of our analysis, it becomes difficult to estimate an accurate derivative for the stream track and speed in these regions. Furthermore, our neural network fits depend on many parameters θ and ϕ that correspond to weights and biases for the stream track and track-speed, respectively. Consequently, the fits in Fig. 5 correspond to the maximum a posteriori estimates of the model parameters θ, ϕ. In a future work, we intend to consider a fully Bayesian treatment of our analysis using variational inference (see Blei et al. 2016 for a review), which enables us to sample from an approximate form of the posterior P (θ, ϕ|D). The result is a probabilistic reconstructed acceleration field, from which confidence limits can be set on galactic accelerations along stellar streams. We have carried out a preliminary test of this approach using Monte-Carlo dropout to estimate the posterior over the model parameters (Gal & Ghahramani 2015), and find that typical deviations from the ground-truth accel-eration field in Fig. 5 fall within the estimated model uncertainty. We postpone a detailed exploration of model uncertainty using variational methods to a future work.
We emphasize that the reconstructed acceleration field in Fig. 5 is the direct result of our analysis-our method does not require a model for the galactic potential to be specified as an intermediate step. This is distinct from previous analyses, which do not have this flexibility. Furthermore, we highlight that the constraints in Fig. 5 are not based on generative modeling; we do not rely on simplified models for the stream distribution function to constrain the acceleration field. Rather, we model the observed stream directly in a data-driven manner. We expand upon these points in §6.

COMPARISON WITH ORBIT FITTING
In this section we compare our method to constraints derived from orbit fitting of stellar streams. Orbit fitting has been used extensively to constrain the galactic potential or other dynamical properties from measurements of stellar streams (see, e.g., Fardal et al. 2006;Koposov et al. 2010;Varghese et al. 2011;Grillmair 2017;Sanderson et al. 2017 Figure 6. Top Row: Mock stellar stream generated in a triaxial logarithmic potential using the "particle-spray" method introduced in Fardal et al. (2015). φ1 and φ2 are angular coordinates which are rotated with respect to the standard equatorial frame. Stream particles are binned by their angular positions, and color-coded by the mean energy relative to the progenitor (around φ1 ≈ −75 deg). For a small band of width δφ1, tracers have a locally similar mixture of energies at slightly different phase-angles. Middle Row: We fit an orbit to the stream illustrated in the top row using a standard χ 2 fitting routine. The binned stream measurements are illustrated by the black points with 1σ error bars, while the red curve is the best fit orbit in the logarithmic potential from Eq. 15. Bottom Row: The black points with error bars and the gray curve are the same as the points and red curve in the middle panel (orbit fitting). We overplot our stream parametrization, color-coded by the energy of the stream track in the ground-truth potential, relative to the progenitor. The color-bar for the bottom row is the same as in the top row. Our method does not require that the stream traces a path of constant energy in any potential, and captures the slowly evolving energy gradient seen among the stream in the top row.
model parameters Γ, and specifying an initial phasespace coordinate (x 0 , v 0 ). Trial orbits are then integrated forward from this initialization in the potential Φ(x|Γ), and the parameters Γ are adjusted between subsequent integrations until the trial orbit passes through the observed stream in phase-space. While appealing in its simplicity, orbit fitting assumes that streams delineate an isoenergy curve in phasespace. In the standard context of stream formation, this cannot be the case since stars become unbound from the stream progenitor with a spread in energies (Helmi & White 1999;Bovy 2014;Johnston 2016). Indeed, Fig. 2 illustrates that stream tracers exhibit a global spread in energy and angular momentum. This can lead to a "misalignment" between the stream track and progenitor orbit, which has been quantified in Sanders & Binney (2013b).
As discussed in §2, our approach does not require a given stream to delineate an isoenergy curve in phasespace under the correct potential, or any potential whatsoever. The phase-space coordinates of the progenitor are also not required. Additionally, our method is model independent, in that it does not require Φ(x|Γ) to be specified. This differs from orbit-fitting algorithms, which require a potential model to be specified beforehand. Our approach also differs from generative methods, which forward model streams or a select number of tracers in a specified potential until a strong agreement is achieved between observational data and simulation (see e.g. Fardal et al. 2013;Gibbons et al. 2014;Bowden et al. 2015;Bonaca et al. 2014;Price-Whelan et al. 2014;Dai et al. 2018;Erkal et al. 2019;Shipp et al. 2021). In addition to requiring a parametric model for the potential, generative methods often require a reliable, efficient prescription for stream formation. Our  Figure 7. The inferred acceleration field in the neighborhood of the stream in Fig. 6 for two methods: orbit fitting (red and cyan), our approach (green), while the ground-truth acceleration a(γ) is plotted in black. The red curve, derived from orbit fitting, assumes the correct model for the galactic potential. The cyan curve assumes an incorrect axisymmetric functional form for the potential. We find orbit fitting can yield a biased estimate of the acceleration field along the stream, since it incorrectly assumes that the stream is characterized by an isoenergy curve in phase-space. Adopting an imperfect model for the galactic potential is also found to yield biased estimates of the acceleration field. Our method is less restrictive, and captures the slowly evolving energy gradient among stream particles as illustrated in Fig. 6 (bottom row). We also do not require an intermediate model for the galactic potential.
method does not require these assumptions or modeling capabilities. Instead, the approach introduced in this work assumes that contiguous segments of the stream are composed of stars with a similar mixture of energies. Hence the change in energy along the stream is small, though not necessarily zero. We next illustrate the bias induced by assuming that stellar streams trace orbits in the galactic potential, and how this bias is alleviated with our approach. We also illustrate bias in estimated accelerations as a result of adopting a slightly incorrect functional form for the galactic potential. To accomplish these experiments, we use the stream illustrated in the top row of Fig 6. This stream is typical in our simulations, and its progenitor is clearly still prominent. In real data, such a stream would be a strong candidate for orbit fitting, since the present day phase-space coordinates of the progenitor can be integrated forward and backward in a specified potential, adjusting model parameters until the orbit passes well through the stream. We carry out this exercise on the synthetic stream in Fig. 6 (top row), using the current phase-space coordinates of the progenitor and 25 other stars at the trailing tail of the stream as trial initial conditions.
Given that we generated the streams used in this analysis, we are at a distinct advantage of choosing an accurate model for the galactic potential. Of course, on real data it is unclear what parametrization should be adopted. We therefore integrate trial orbits in two potentials. The first is the correct model for the potential, Eq. 15, with model parameters Γ = (v c , r h , q 1 , q 2 , q 3 ). The second is an axisymmetric version of the same po-tential, with fixed parameters q 1 = q 2 = 1. We allow q 3 to remain a free parameter, characterizing the z-axis flattening. This is a common choice for the potential in previous works, which have applied orbit fitting to real stellar streams (e.g., Koposov et al. 2010). Trial orbits are integrated in the two potential models separately, and compared to binned stream measurements in a rotated ICRS frame. Model parameters are adjusted until the χ 2 between the observed stream and integrated orbit is minimized for each initial condition using a Scipy (Virtanen et al. 2020) optimizer. This approach is similar to Koposov et al. (2010), where orbit fitting is used to constrain a potential model from the GD-1 stream. The best fit orbit for the stream in Fig. 6 (top panel) using the correct model for the galactic potential is illustrated in the middle row of Fig. 6 in red, where the orbit is plotted in a rotated ICRS frame. The quantities µ φ1 , µ φ2 , µ r , and dist. indicate the two proper motion components, radial velocity, and distance along the orbit, respectively. Binned stream measurements are shown by the black points with error-bars, where the error-bars are derived from the intrinsic scatter of stream observables in each φ 1 bin.
We find that the 5-parameter potential, Eq. 14, is flexible enough to generate an orbit which passes through the binned stream measurements. However, because the stream does not actually delineate an isoenergy curve, the resulting best fit orbit and corresponding potential yield estimates of the acceleration field which are not necessarily accurate. We illustrate this bias in Fig. 7, where we compare our inference of the acceleration field along the stream (green) to orbit fitting (red) against the ground-truth acceleration (dashed black curve) for the cartesian components a = (a x , a y , a z ). The cyan curve is the result of integrating trial orbits in an incorrect axisymmetric model for the galactic potential as previously described. We include this curve to illustrate bias induced by adopting an incorrect functional form for the potential. Both the red and cyan orbit-fitting curves are taken from the negative spatial gradient of the best-fit potential corresponding to the best-fit trial orbits (e.g., Fig. 6 middle row). Estimates of enclosed mass at the progenitor's location for the best fit potential derived from orbit fitting is biased from the truth by roughly 10% when using the correct potential model, and up to 25% when imposing axisymmetry. For the method presented in this work, we find fractional errors on enclosed mass at the sub-percent (∼ 0.9%) level for the dashed black curve in Fig 7. Contrasted with orbit fitting, our method allows for the possibility that globally, streams depart from a single orbit. Indeed, for most regions of the stream in Fig. 7, the true acceleration components and the components inferred from our method are similar. Furthermore, our method does not require priors on the functional form of the potential: this enables increased flexibility and significantly reduces systematic errors from adopting an incorrect model for the galactic potential. Fluctuations in the inferred accelerations are attributed to model uncertainty in our neural network based-fits, which we intend to quantify in a future work using variational inference. We also note that for some streams analyzed in this work, orbit fitting can provide unbiased estimates of the acceleration field along the given stream under the correct model for the potential. However, all streams yield biased constraints on local accelerations under the incorrect functional form for the galactic potential. Substantially increasing the flexibility of the potential by adding more model parameters is one possible solution, though this is also found to yield noisy estimates of the local acceleration field when using orbit-fitting techniques.
In Fig. 6, bottom row, we show our parametrization of the stream track in the same rotated ICRS frame, with the best fit orbit from the top row plotted in grey, and the neural-network-derived stream track illustrated by the multi-colored curve. The stream track is colorcoded by its energy in the correct underlying potential (relative to the progenitor). We find an energy gradient, which closely follows the mean relative energy of tracers in the top row of Fig. 6 (note that the color-scale in the top row of Fig. 6 and bottom row are the same). Taken jointly, Fig. 6 and Fig. 7 illustrate that our approach is significantly less restrictive than orbit-fitting, since we recover the local energy gradient of stream trac-ers (Fig. 6, bottom row) while also recovering the correct acceleration field along the stream (Fig. 7, green curve). Both of these results are achieved without having to model the gravitational potential directly. This is advantageous, since inaccurate models for the potential can lead to bias, even for less restrictive potential reconstruction methods based on generative modeling.
We also use Fig. 7 to illustrate that our analysis can be applied to streams with missing observations. For all streams analyzed using our method in this work, we remove the progenitor system (see §2.2 for discussion). When estimating the acceleration field in Fig. 7, we do not rely on any tracers that fall within the grey band of each panel since this region is roughly the location of the progenitor cluster. This corresponds to the tracers centered around φ 1 ≈ −75 deg in the top row of Fig. 6. We find that our method is capable of flexibly interpolating over this missing information. This is useful for real stellar streams, since the observed tracers are generated from both the stream distribution function and the observational selection function; our method does not require these distributions to be disentangled.

POTENTIAL RECONSTRUCTION
From a Fisher-information analysis of stellar streams, Bonaca & Hogg (2018) found that streams provide a highly localized measurement of the galactic halo. Because of this, they suggest fitting streams independently, and combining constraints on the galactic potential at the global level as a separate step. Our method provides the framework to accomplish this, since we estimate the local acceleration field along each stream rather than global properties of the halo. In this section, we explore how our local estimates of the galactic acceleration field can be combined to constrain a global gravitational potential.
We consider a static potential, though will explore time-dependent potentials in a future work. In principle, a static potential is not strictly required in our analysis. The central requirement of our work is that the stream forms a coherent phase-space object, with the property that adjacent segments of the stream are populated by stars with a locally similar distribution of orbits. For streams with a strong misalignment between the stream-track and stellar trajectories, this assumption is broken. A track-velocity misalignment has been shown to occur in time-dependent potentials with (e.g.) a Large Magellanic Cloud component (Erkal et al. 2019). Still, depending on the magnitude and evolution of the misalignment across the stream our analysis could still be applicable. We provide further discussion of a time-dependent potential in §6.3 and discuss extensions to the present work.
The gravitational potential is related to the acceleration field through Eq. 7. This relation can be used to find the potential, Φ(x), which gives rise to the local acceleration field in the neighborhood of independent streams. Fitting a model for the potential in this framework will typically be computationally inexpensive, provided that we have already estimated the acceleration field along independent streams using the method introduced in this work. This is unique and distinct from previous analyses, which constrain the galactic potential or other dynamical properties of the Milky Way using generative methods and several forward simulations (see, e.g., Fardal et al. 2013;Bonaca et al. 2014;Price-Whelan et al. 2014;Gibbons et al. 2014;Bowden et al. 2015;Pearson et al. 2015;Dai et al. 2018;Erkal et al. 2019;Shipp et al. 2021). These works typically fix an analytic model for the potential with parameters free to vary. However, choosing a model in the first place imposes a functional prior on the potential. That is, a given potential model is never guaranteed to describe the actual observed galaxy. Our method alleviates the need to specify a potential model beforehand, letting the data drive our inference. However, priors on the functional form of the potential can still be adopted (e.g. by assuming that a can only depend on galactocentric distance r).
Using the methods outlined in this work, a range of analytic models for the potential can be constrained quickly using standard optimization algorithms, enabling model comparison techniques which would be too expensive under the generative framework. Alternatively, a highly-flexible representation for the potential can be constrained using (e.g.) basis function expansions or a flexible neural network to represent Φ. This has the advantage of preserving the detailed features of the acceleration field around each measured stream, while interpolating between independent systems. Provided that the data have enough constraining power, flexible models minimize the need to construct an accurate analytic form for the potential beforehand.
We therefore consider two pathways for potential reconstruction. The first fits a standard analytic model for the potential to our estimate of the acceleration field along the six stellar streams analyzed in §3.2. This enables one to derive constraints on the parameters of the model for the Galactic potential. The second uses a neural network to describe the potential that gives rise to the estimated acceleration field in the vicinity of each stream, enabling a flexible stream-based analysis of the galactic halo that is fully data-driven.
The output of our main analysis in §3.2 are acceleration components sampled across several streams. In the subsequent sections, we denote the set of these estimates with {a(x i )} i , where x i is the i th position at which the acceleration field is evaluated.

Constraining an Analytic Potential
In this section we fit an analytic model for the potential to our estimate of the acceleration field along the six independent streams in §3.2. A wide range of parametric models have been used to describe the Galactic potential. Choosing the wrong model can lead to bias (Bonaca et al. 2014, and our Fig. 7), making model selection an important step in constraining the potential. For real data, it is not obvious what model should be adopted beforehand. On simulated data we are provided with a distinct advantage: we know the correct analytic form of the potential. Therefore, in this section we adopt the correct analytic form of the ground-truth potential, Eq. 15, and allow its parameters to vary. However, we emphasize that our analysis enables the use of model comparison techniques without the added computational cost of simulation based inference or generative modeling. A detailed statistical model comparison of the Milky Way potential using stellar streams and e.g. the Bayesian evidence framework has not yet been carried out. This is likely due to the large computational cost of generating a new set of simulations for each model and its corresponding parameter grid. Because our analysis constrains accelerations rather than the potential, it enables us to compare analytic potential fits rapidly. We intend to explore this aspect of our method in a future work with data from real stellar streams.
We constrain the logarithmic potential in Eq. 15 with the five parameters Γ = (v c , r h , q 1 , q 2 , q 3 ) using a Scipy (Virtanen et al. 2020) non-linear least-squares optimizer. In particular, we minimize the χ 2 over the parameters Γ, where σ 2 i are the estimated uncertainties on the inferred acceleration field at each evaluation point x i . In §3.2, we discuss how variational inference can be used to obtain an estimate of this uncertainty, though we set all σ i ≡ 1 for now. We estimate errors on the model parameters Γ by bootstrap resampling the inferred acceleration field along each stream. In particular, we minimize Eq. 16 by sampling the acceleration field at 100 random positions across all streams. We repeat this many times, minimizing Eq. 16 for each realization to estimate errors on the best fit model parameters. We impose the con- straints v c ∈ [50, 400], r h ∈ [.2, 10], q i ∈ [0.1, 5] for i = 1, 2, 3. These constraints are not strictly necessarily, though aid in optimization since the logarithmic potential in Eq. 15 is degenerate in r h and the flattening parameters q 1 , q 2 , q 3 . That is, if one replaces r h with α −1 r h and (q 1 , q 2 , q 3 ) → (αq 1 , αq 2 , αq 3 ) in the potential (Eq. 15) for any constant scale factor α > 0, then Φ changes by an additive constant. However, the underlying accelerations remain unchanged. The resulting distribution of best-fit parameters is illustrated in Fig. 8, where grey lines indicate the true parameter values and contours correspond to regions of 68 and 95% confidence. From Fig. 8, we find that our method for estimating the local acceleration field along independent streams provides combined parameter constraints which fall within the 68% region of the true potential parameters. Adopting different analytic models for the potential Φ in Eq. 16 is straightforward to implement, since our estimate of the acceleration vector, a, along each stream is independent of any particular parametrization of the potential.

Constraining a Flexible Potential
The true Galactic potential and its dark halo are not necessarily well captured by simplified analytic (e.g. spherical) models (Allgood et al. 2006;Law et al. 2009;Vera-Ciro et al. 2011;Rojas-Niño et al. 2015;Alexander et al. 2020;Bonnet et al. 2022;Garavito-Camargo et al. 2021). Therefore, model selection becomes an important step in constraining the galactic and halo potential from stellar streams. Provided that the data have strong constraining power, flexible models for the gravitational potential have the ability to capture a diverse range of features across the galaxy while reducing the bias induced by adopting the wrong parametric model. Basisfunction expansions provide one pathway to parametrizing a flexible potential, and can be fit with standard least-squares optimization routines or MCMC.
We consider an alternative approach, which uses a flexible neural network to describe the Galactic potential. This has the advantage of increased flexibility over basis-function expansions. In particular, a basis and truncation condition do not need to be specified. However, we emphasize that basis expansions can equally be applied to constrain a potential from our estimates of the acceleration field. Following the approach of Green & Ting (2020); Green et al. (2022), we represent the potential Φ with a flexible fully-connected neural network Φ β (x), where β are the parameters of the neural network. Using estimates of the acceleration field localized to several streams, Φ β can be trained through its spatial gradient to reproduce the local acceleration field of stellar streams while interpolating over the unsampled regions. Poisson's equation connects the gravitational potential to the matter density, where G is the gravitational constant and ρ(x) is the matter density of the Galaxy. In order to ensure that the inferred potential is physical, we would like to encourage the condition to promote a non-negative mass density. To accomplish this, we use a loss function that penalizes negative mass densities during training. The loss function is where x * i is sampled within the region of interest (e.g. a 3-dimensional box within which all streams are contained). The λ term in Eq. 19 works to penalize randomly sampled regions with a negative mass density. In practice, we also sample positions along each stream. To encourage smoothness, weight decay is also used to penalize large neural-network parameters. We use a .0225°0 .0200°0 .0175°0 .0150°0 .0125°0 .0100°0 .0250°0 .0225°0 .0200°0 .0175°0 .0150°0 .0125°0 .0100

Residuals°0
.030°0 .025°0 .020°0 .015°0 .010°0 .030°0 .025°0 .020°0 .015°0 .010 weight-decay value of ∼ 10 −5 . We apply this method to the six stellar streams whose local acceleration fields are inferred in Fig. 5, with λ = [0.05, 0.5, 1] after ∼ 200, 500, and 800 training epochs, respectively. We illustrate the reconstructed potential in a Galactocentric view in Fig. 9 for two different slices: the top row shows a constant z slice through our mock galaxy, while the bottom row shows a slice with constant y. The first column depicts the ground truth potential (Eq. 15), the middle column is the inferred potential, Φ β , and the right column depicts fractional errors. We note that the color maps in the first two columns have the same range, allowing for a visual comparison between the true potential and our inference. Across a range of scales, we recover an accurate depiction of the potential with fractional errors typically at the few-percent to sub-percent level. The only imposed constraint in our inference of the potential in this section is the non-negativity of mass. The neural network Φ β knows nothing of the triaxility of the potential beforehand: it discovers this from the inferred acceleration field directly. Using the acceleration field sampled from only six stellar streams (Fig. 5), our flexible neural-network-based reconstruction reproduces the features of the global triaxial potential across a range of scales, even in intermediate regions where mock-streams do not reside. We find the most significant bias at z = 0 near the galactic center, since streams were not generated in this vicinity; they would quickly disrupt anyway due to extreme tidal forces.
The residuals in Fig. 9 are quite structured, especially in the xy-plane. Fractional errors tend to increase away from the six streams, leading to the stream-like structure in the residuals of Fig. 9. This is in agreement with Bonaca & Hogg (2018), who found that fitting streams under flexible potential models constrains the galactic acceleration field localized to the stream(s) of interest. The pericenters of the six streams analyzed fall between [4.5, 8] kpc, with apocenters between [7.5, 16] kpc. Beyond ∼ 16 kpc fractional errors increase.
It appears somewhat surprising that the neural network-based potential can accurately interpolate over unsampled regions of the galactic acceleration field. We aim to explain this in the following way. When fitting an over-parametrized linear regression model using gradient descent, one implicitly optimizes for the smoothest possible solution that is also compatible with the data (Arora et al. 2019). In particular, the optimal parameter set is the pseudoinverse solution, which has the smallest L 2 norm over the model parameters. As an analogous system, we expect that the same is true for deep neural networks trained with gradient descent. Because the L 2 norm over neural-network parameters is directly correlated with the Lipschitz constant of the network (defined as sup x [f (x)]; Gouk et al. 2018), the optimal neural-network parameters also implicitly minimize the upper bound of this constant. In addition, weight decay explicitly minimizes the L 2 norm over the parameters. The result is that neural-network-based regression yields the smoothest possible solution for which the upper bound of the Lipschitz constant is minimized while also maintaining compatibility with the data. For this reason, along with the enforcement of a potential with a non-negative mass density, it is actually not unusual that the best-fit neural-network potential smoothly interpolates over regions without measured accelerations. This is expected for a neural network that characterizes the measured portions of the acceleration field, while also optimizing for smoothness.
The advantage of the neural-network-based potential presented in this section is that one does not need to specify an analytic model for the potential. In this way, a neural-network description enables potential reconstruction from a fully data-driven standpoint. Bonaca & Hogg (2018) demonstrate that flexible models for the potential characterize the local acceleration field in the vicinity of a given stream. A neural-network-based model is flexible enough to fit a potential that precisely reproduces the local acceleration field inferred from independent streams, while also characterizing the potential at a global scale. This is enabled by increased flexibility, and the ability to minimize non-physical artifacts through Eq. 18. Flexible potential models come at the cost of reduced model interpretability, though we intend to explore symbolic regression (e.g., Cranmer et al. 2020) applied to our flexible inference of the potential in a future work. 6. DISCUSSION

Summary
We have introduced a new method to constrain the galactic acceleration field from measurements of stellar streams. Our method is agnostic to the functional form of the underlying potential and does not require streams to delineate orbits. Rather, our central assumption is that streams are populated by an ensemble of stars characterized by a local mixture of orbits that changes gradually along the stream. Therefore, our approach does not require streams to delineate isoenergy curves in phase-space. Rather than relying on generative models and simulation based inference, we fit the dynamics of the bulk motion of the observed stream directly as a function of position along the stream track. Using a flexible neural network parametrization of the stream track, derivatives of dynamical quantities (e.g. position and velocity) can be calculated, providing a direct link to the local acceleration field in which the stream resides. Our method does not require approximate transforma-tions to action-angle coordinates, but has the benefit of operating in phase-space.
Our method for sampling the galactic acceleration field with stellar streams is methodologically different from previous works, which typically rely on orbit fitting, generative models, numerous orbit integrations, or action-angle coordinate transformations to constrain an analytic potential model using stellar streams. We have compared our approach to orbit fitting in §4, where we demonstrated that our approach alleviates inaccuracies in the inferred accelerations induced by assuming that streams delineate isoenergy curves in phasespace. Because our approach estimates accelerations directly without the intermediate requirement of specifying a model for the potential, we also alleviate bias induced by adopting an incorrect model for the potential (Fig. 7). This enables a description of the acceleration field which is independent of any user-defined functional form for the potential, allowing us to accurately probe halo morphologies with complex or unexpected geometries. Furthermore, we do not require a generative prescription for stream formation, nor do we require knowledge of e.g. the present day phase-space coordinates of the stream progenitor. Instead, we assume that small neighboring segments of the stream are populated by stars with a similar distribution of orbits. This makes our analysis applicable to e.g. Gaia-detected stellar streams for which the progenitor has completely disrupted, or streams which are not well modeled in current generative frameworks. Bonaca & Hogg (2018) suggested that stellar streams probe highly localized properties of the underlying galactic potential; our analysis exploits this as a feature of stellar streams to constrain the galactic acceleration field on small scales. As a post-processing step, we demonstrate that highly localized measurements can be combined across several streams to place constraints on a global gravitational potential. Indeed, we find that our method is able to reproduce the correct parameters of a global ground truth static potential without bias ( §5.1). We have also demonstrated that a highly flexible potential model can be constrained, enabling potential reconstruction from a fully data-driven standpoint ( §5.2). For more complicated potentials beyond a triaxial model, the analysis presented in this work has the benefit of flexibility and does not require prior knowledge of the functional form for the gravitational potential. We intend to explore this benefit of our analysis in further detail, through an application to cosmological simulations and real data of Milky Way streams.

Limitations
The Milky Way is currently undergoing a major merger with the Large Magellanic Cloud (LMC; e.g. Besla et al. 2007;Peñarrubia et al. 2016;Erkal et al. 2020Erkal et al. , 2021Conroy et al. 2021). For streams in the Galactic Southern hemisphere in particular, the LMC is likely to impart time-dependent perturbations on the morphology and kinematics of these systems (Erkal & Belokurov 2015;Shipp et al. 2021). Indeed, recent work has explored the dynamical influence of the LMC on the Tucana, Orphan and Sagittarius streams, since these systems may have suffered a direct interaction or close passage with the LMC (Erkal et al. 2018;Koposov et al. 2019;Vasiliev et al. 2021;Lilleengen et al. 2022). Direct interactions with the LMC can lead to a significant misalignment between the stream track and proper motions, which has been shown to occur only in the presence of time-dependent perturbations (Erkal et al. 2019). This misalignment has been used to measure properties of the LMC such as its mass (Erkal et al. 2019;Shipp et al. 2021), but poses a potential challenge in applying our analysis to these systems. The method introduced in this paper requires the local motion of stream stars to roughly fall along the stream track. When this assumption is broken, our method in its current formulation is no longer directly applicable. However, because our analysis has the benefit of self-regulation (see §2), it will be clear when a given stream is a poor candidate for our approach. In a future work, we intend to consider a more general formulation of our methodology which will be applicable to streams which have been influenced by time-dependent perturbations, and therefore have a misalignment between proper motions and the stream track. We discuss one possible formulation in §6.3. Still, there are a number of streams which do not appear to have suffered from time-dependent perturbations (e.g. see the "golden sample" used in Malhan et al. 2020). We therefore anticipate that several streams can be explored using the methods presented in this work, particularly in the galactic Northern hemisphere.
Recent work has detected a relative velocity difference between the outer Milky Way halo (r ∈ [40, 120] kpc) and the galactic disk, thought to be a result of the Milky Way's dynamical response to the LMC (Garavito-Camargo et al. 2019;Erkal et al. 2020;Petersen & Peñarrubia 2020;Erkal et al. 2021). This could lead to a velocity misalignment between the stream track and local motion, as a result of the non-inertial reference frame of the disk. Malhan et al. (2020) has measured the misalignment between the proper motions of stream stars and stream track for streams within ∼ 25 kpc of the galactic center, and found that the misalignment was in agreement with the solar reflex velocity. This indi-cates that the inner halo may have a negligible relative velocity compared to the galactic disk. We therefore anticipate that our analysis will be less sensitive to systematic errors induced by the Milky Way's reflex motion for streams within a ∼ 25 kpc volume.
Our method will be most applicable to static halo potentials, though this condition is not strictly required. In particular, our core assumption is that streams contain stars with a locally similar mixture of orbits at slightly different phase angles (Fig. 2). Provided that the potential is evolving slowly, integrals of motion (e.g. actions) will remain adiabatically invariant so that this assumption will not necessarily be broken. Even for more strongly time-dependent halo potentials in which orbital deviations are non-adiabatic, orbital actions may remain locally comparable so that we can still treat streams as locally similar mixtures of stellar orbits (Buist & Helmi 2015;Vasiliev et al. 2021). Because our method for estimating galactic accelerations works by measuring changes in local dynamical properties along a stream, streams with a high number density of stars and ample phase-space observations will provide more accurate measurements of the galactic acceleration field over recent timescales. That is, if we view the mean stream track as an ensemble of several distinct orbital segments, then the length of each measurable segment becomes smaller as the number of measurable stream stars increases. Provided that the potential is changing on timescales longer than the time it takes for a characteristic star to traverse its orbital segment, our method could provide an accurate "snapshot" of the present day potential. For a more sparsely populated stream, the measured orbital segments become larger, and our analysis becomes more sensitive to the evolving potential. In this regime, we might estimate a mean time-dependent potential, since as stars continue along their orbits the stream could deform. Indeed, the Milky Way halo may be evolving as a result of (e.g.) the LMC and its dark matter wake (Garavito-Camargo et al. 2019Conroy et al. 2021;Lilleengen et al. 2022) or dark matter figure rotation predicted by ΛCDM (Dubinski 1992;Bailin & Steinmetz 2004;Bryan & Cress 2007;Valluri et al. 2021). These effects may lead to detectable morphological deviations in stellar streams, though further work is still needed to explore the impact of these timedependent effects on the local dynamics of streams. In a future work, we intend to apply our method to simulated streams generated in a representative time-dependent potential.
Dark matter subhalos can create gaps in stellar streams as they orbit throughout the galaxy. The GD-1 stream in particular shows signs of such an encounter, due to a gap and spur of stars which bifurcate the stream (Bonaca et al. 2019b;de Boer et al. 2020). Subhalo interactions which lead to gaps along a stream are not particularly challenging to handle in our analysis, and do not need to be modeled separately. Rather, we have demonstrated in §4 that our analysis can interpolate over missing segments of a given stream. Spur-like features that lead to distinct arms can be modeled as independent streams with their own track.
Our main analysis has considered measurements of the full 6d phase-space coordinates of stream tracers. This information is already available for several streams and will be measured with increasing precision across a larger number of systems with the advent of (e.g) Gaia DR3. However, we have also implemented a version of our analysis which constrains the galactic acceleration field using incomplete phase-space measurements with missing dimensions. This makes our approach applicable to the large volume of Gaia-detected streams with missing line-of-sight velocities . We expand upon the case of missing phase-space dimensions in §6.3.

Future Directions
In this paper we have established a new avenue for constraining the galactic potential using measurements of stellar streams. We intend to apply our method to streams with measured 6D phase-space dimensions in a future work. This will enable highly localized estimates of the galactic acceleration field, independent of any particular functional form for the gravitational potential.
A large number of streams have been measured with Gaia (see e.g., Ibata et al. 2019;Borsato et al. 2020), though most stream stars have 5D astrometric solutions (two angular coordinates, 2 proper motions, and a parallax angle). Even then, parallax measurements are not necessarily reliable for dim stars in the stellar halo. We have explored a modification to our analysis which is able to constrain components of the galactic acceleration field along a given stream in the absence of certain phase-space dimensions. In particular, the modeling introduced in §2 is based on vector differential equations (e.g., Eq. 6). These vector expressions can be decomposed into components, allowing us to model streams with incomplete phase-space observations. For instance, if line-of-sight velocities are unavailable then Eq. 6 can be re-written in terms of 3D positions and 2D perpendicular velocity components in the plane of the sky. We can then infer the component of the galactic acceleration vector in the plane of the sky. Without distances we can infer angular accelerations, however, one still needs to estimate the heliocentric distance to the source in order to subtract the solar reflex motion. This can be done using the Galactic Parallax method introduced in Eyre & Binney (2009). Alternatively, photometric distances can be estimated (see, e.g., Shipp et al. 2018). In this way, our analysis is widely applicable to the large number of stellar streams which have 4D astrometric solutions from Gaia . Furthermore, the upcoming Gaia DR3 will include low-resolution spectra for millions of stars, which will be utilized to estimate isochronal distances. These measurements will make our analysis applicable to a large number of streams, providing small-scale constraints on the shape and mass of the Milky Way dark halo. For simplicity, we do not analyze mock-streams with incomplete phase-space observations in the present work. However, this is the subject of future work.
Observations of external galaxies have revealed tidal debris, thought to be remnants of accreted satellite galaxies (see, e.g., Kado-Fong et al. 2018). For M31 in particular, there is mounting evidence for accreted substructure including the Giant stellar stream and a shelffeature also indicative of tidal disruption (Ibata et al. 2001;Ferguson & Mackey 2016;Conn et al. 2016). Radial velocities have been measured for a number of stars in the giant stellar stream, and photometric measurements have constrained a distance gradient along the stream (Ibata et al. 2014;Conn et al. 2016). Our analysis can be applied to these measurements, to obtain a direct estimate of the line-of-sight accelerations along the Giant stellar stream in M31. The Nancy Grace Roman Space Telescope (Spergel et al. 2015) is also expected to detect a large number of streams within M31 and other galaxies in the local volume (Pearson et al. 2019(Pearson et al. , 2022. If radial velocities are obtained for these systems and a distance gradient can be estimated along the stream (using e.g. photometry; Conn et al. 2016), our method will provide an exciting avenue to measure the small-scale structure of the dark halo in external galaxies, placing the Milky Way in a broader context.
The methods developed in this paper are philosophically similar to recent work that discusses measuring the galactic acceleration field using e.g. high-precision radial velocity measurements, pulsar timing arguments, or eclipsing binary star timing arguments (Chakrabarti et al. 2020(Chakrabarti et al. , 2021(Chakrabarti et al. , 2022. In particular, we rely on changes in dynamical properties along a stream to estimate accelerations directly. Analogous to these works, our approach provides highly localized measurements of the galactic acceleration field in which the observed stream resides. Such localized measurements of the dark matter distribution in the stellar halo could prove useful for testing the standard cosmological model, ΛCDM, on small-scales  in a way that does not rely on generative prescriptions for stream formation, stringent energy requirements, and analytic userdefined models for the galactic potential. Comparing our analysis with other small-scale measurements of the dark matter distribution (e.g. gravitational lensing) will enable a test of systematic errors which may bias such highly localized constraints, since different methods are not necessarily subject to the same systematic errors. Our method therefore provides a potentially useful addition to techniques that measure galactic accelerations directly.
In §6.2 we have discussed the limitations of our analysis in a time-dependent potential. Time dependence can lead to a significant misalignment between the local morphology of a stream and the trajectories of its tracers (Erkal et al. 2019). This poses a potential challenge for our analysis. We have begun a preliminary exploration of a modification to our approach that can estimate accelerations even in the presence of such misalignments, by allowing for the possibility that the stream track is not directly coupled to stellar velocities. As discussed in §6.2, we do anticipate that a large number of streams are strong candidates for the analysis presented in this paper already, so we intend to explore further generalizations in a future work.
The aim of this work is to introduce and demonstrate a new method for constraining galactic accelerations from measurements of stellar streams. As such, we have not incorporated a detailed discussion of uncertainty quantification using our method. From a preliminary analysis, we have implemented variational inference (Ganguly & Earp 2021;Gal & Ghahramani 2015) as a technique to provide estimates of model uncertainty on the inferred acceleration field along a stream. Characterizing model uncertainty is especially important on real data, where the ground-truth potential is unknown and the measurements themselves are uncertain. We therefore postpone detailed modeling of uncertainties to a future work, though we note that variational methods make uncertainty quantification using neural networks feasible. We have also already laid the groundwork to incorporate measurement errors in §2.2.
Finally, we conclude by considering a future application of our approach to systems beyond stellar streams. In particular, our method for estimating galactic accelerations requires the existence of a coherent structure in phase-space, with the property that, on average, nearby tracers have similar energies at slightly different phaseangles. In the galactic disk, level-sets of constant element abundances have been shown to roughly delineate stellar orbits . This is expected for a well-mixed equilibrium population, where abundances are analogous to orbital actions and an isoabundance contour is distributed over a range of conjugate angles. Indeed, this formulation has been applied to estimate the mass of the Milky Way disk in Price-Whelan et al. (2021). In addition to this measurement, if an isoabundance curve can be estimated in phase-space using our approach, it would provide a new avenue to constrain the galactic acceleration field in the Milky Way disk without imposing functional priors on the form of the disk potential. We intend to explore this exciting application in a future work. DETERMINING γ The stream track described in §2.2 is a function of the scalar parameter γ, which encodes position along the track. For streams which do not wrap onto themselves, an on-sky position angle can be used to assign an ordering to stream stars. We denote the ordered set of N stream tracers with {T 1 , T 2 , ..., T N }, where T i represents the i th tracer. We choose the convention that increasing i corresponds to motion along the trajectory of the stream's orbit in phase-space (rather than against this trajectory). Once an ordering has been determined, a scalar parameter γ i can be assigned to each tracer, such that for the ordered set {T 1 , T 2 , ..., T N }, γ 1 < γ 2 < ... < γ N .
(A1) Figure 10. Schematic outline of stream unwrapping procedure. Top Left: Input phase-space data, where subsequent points in the dataset are connected by blue lines. Top Middle: We implement a nearest neighbors algorithm to connect tracers which have a small euclidean distance in phase-space. There is an added momentum condition, which encourages the nearest neighbor graph edges in position space to fall roughly tangent to the local direction of motion of the stream. The stream is color-coded by the initial ordering, γinit. Top Right: The momentum condition inevitably leaves out a number of stream members from the nearest neighbors graph, requiring a flexible interpolation to assign an ordering (i.e., γ-value) to tracers which were initially left out of the graph. Interpolation is performed with a neural network, which outputs a membership probability for tracers in addition to the ordering parameter, γ. Bottom Right: Output of the interpolation neural network, color-coded by the interpolated γ-values across the stream. Blue points indicate tracers with the lowest probability membership, though all stars have membership probability > 0.99 since they are stream members. Bottom Middle: The interpolation neural network (top right) is connected to a new neural network called Param-Net, which takes the ordering parameter γ as an input, and outputs the mean position along the stream. γ-values are further refined using the same momentum condition as before: namely, that increasing γ corresponds to motion along the stream from trailing to leading arm. Bottom Left: Example output of the final unwrapped stream, color-coded by the final γ-values.
In this way, γ is analogous to a phase-angle along the ensemble of orbits that streams delineate. Because only the rank-ordering of γ i is meaningful, the numerical value is ambiguous. We therefore adopt the simple convention γ i ∈ [−1, 1]. In practice, we divide this interval by the number of stream tracers, such that γ i+1 − γ i is a positive constant. Eq. A1 implies that once the ordering {T 1 , ..., T N } is fixed, it is straightforward to assign a γ i to each tracer. Determining a suitable ordering, however, is more challenging. To determine an ordering for stream tracers and thus their corresponding γ-values, we adopt a multi-step process illustrated in Fig. 10. We describe the panels of this figure sequentially to highlight our approach.

A.1. Nearest Neighbors with Momentum
An unordered set of phase-space measurements is illustrated in the top left panel of Fig. 10, where subsequent tracers in the dataset are connected by a blue line. The result is a graph with tracers as nodes, though subsequent connections are randomized. We assume that all streams in our analysis are initialized in this way (i.e. with random ordering). For this particular stream a suitable ordering cannot be prescribed using an onsky position angle alone, since the stream wraps onto itself. This differs from the stream in Fig. 2, which can be parametrized in terms of φ 1 . Instead, we adopt a modified nearest neighbors algorithm which connects subsequent tracers in phase-space, with the momentum condition that graph-edges which connect neighboring particles roughly align with the local direction of mo-tion of stream tracers in position space. We provide a pseudocode version of the algorithm in Algorithm 1, where the momentum condition is enforced through the CosSimDist penalty term. We find that this term is redundant for simple stream morphologies. However, for streams which loop back onto themselves one or more times, the momentum condition helps break degeneracies in the nearest neighbors search. By adding the CosSimDist term in Algorithm 1, nearby tracers with velocity vectors pointing along the local graph edge will be preferred; this helps ensure that the ordering induced by the algorithm does not jump from one stream arm to another at an intersection. Additionally, the momentum condition helps ensure that the graph will not be dominated by a few anomalous stars with atypical velocities.
We illustrate the result of the nearest neighbors with momentum algorithm in the top-middle panel of Fig. 10. This panel contains the same stream as in the top left panel. Now that an ordering has been induced for a number of particles which populate the graph, we are able to assign a γ value to these tracers (Eq. A1). We denote these starting values with γ init , since they will momentarily be updated and smoothed over. The leading arm is illustrated in the blue-range of the color-bar, while the trailing arm is in the red-range. The dashedline connects neighboring graph nodes, similar to the blue lines in the top left panel, though now ordered. From this panel, it can be seen that a number of tracer particles were left out of the graph. This is due to the momentum condition, which keeps the nearest neighbors algorithm moving forward along the local trajectory of stream tracers. In order to assign a γ-value to the remaining tracers which were left out of the graph, we utilize a neural network to perform a flexible interpolation.

A.2. Autoencoder Neural Network
From the ordered, incomplete set of stream tracers, we now discuss an interpolation method to assign a γ-value to stars initially left out of the graph. To determine a γ-value for these tracers, we utilize a neural network with an autoencoder architecture, where the bottleneck of the autoencoder is the γ-parameter.
At a conceptual level, the autoencoder is split into two parts. The first we call the Interpolation network ( Fig. 10; top right panel), while the second we call Param-Net ( Fig. 10; bottom row, middle panel). The interpolation network takes 6D phase-space coordinates as an input, and outputs a γ-value and a new scalar p ∈ [0, 1], which indicates stream membership probability. While including velocity components in the interpolation neural network is not strictly required, we find that it helps break degeneracies when unwrapping streams which cross through the same spatial position several times.
We adopt a two-step training process for fitting the parameters of the autoencoder. First, we train the interpolation neural network as follows. From the nearest neighbors with momentum algorithm, we have already determined an initial set of γ-values for a number of tracers; namely, γ init . We adopt these values as labels, and utilize the interpolation network to predict new labels for tracers which were left out of the graph. For tracers which were connected in the initial graph, we assign them each a p i ≡ 1, and train the neural network on these labels as well. We randomly sample points in a phase-space box which encloses the observed stream, assigning these randomized points p ≡ 0. The loss function-dependent on the neural network parameters, θ-is constructed as follows, where (x i , v i ) are the phase-space coordinates of tracers which populate the graph from the nearest neighbors algorithm, and (x i,rand , v i,rand ) are randomly drawn coordinates within a 6D phase-space box that encloses the stream. Training through this loss-function provides a mapping from phase-space to the scalar parameter γ. This enables us to assign a γ-value to tracers that were initially left out of the graph. Furthermore, the probability membership, p θ , is trained such that p θ ≈ 1 will typically correspond to phase-space coordinates which are similar to those used in the training set (e.g. from the nearest neighbors analysis). Alternatively, p θ < 1 indicates that a particular tracer may not be assigned an accurate γ at high-confidence, since it is dissimilar (i.e. out of distribution) compared to the training set. We achieve this behavior from the random sampling of (x i,rand , v i,rand ) in Eq. A2. Because stream members will typically occupy a narrow region in 6D phase-space, randomly sampled coordinates within a phase-space box are unlikely to frequently fall directly along the stream simultaneously in all 6 phase-space dimensions for a finite number of random samples. Additionally, stream members are already trained with p i ≡ 1 over a narrow phase-space volume, thereby reducing contamination from the random sampling scheme. An example output of the interpolation neural network is illustrated in the bottom right panel of Fig. 10, where we have assigned a γ-value to all stream tracers using the interpolation, even for those which were originally left out of the nearest neighbors graph. All stream tracers are found to have p > .99, indicating high probability membership. We plot the top ∼ 20 tracers with the lowest probability membership in the cyan points. Stars in the extended tidal tails have the lowest probability membership, though stars more centrally located along the stream have similar, higher membership probabilities as expected.
We next consider the second half of the autoencoder, Param-Net (Fig 10; bottom row, middle panel). Param-Net characterizes the mapping from γ to position space, (x, y, z). Param-net provides the initial seed for the stream track neural network discussed in §2.2. The stream track neural network differs from Param-Net in that its γ-values are treated as fixed. We again utilize a fully-connected multi-layer perceptron neural network for Param-Net. If we treat the γ-values assigned by the Interpolation neural network in the top right panel of Fig. 10 as fixed values, learning the mapping from γ − → x is straightforward, since each stream tracer already has its own x i and estimated γ i from the Interpolation network. However, a careful inspection of the particles in the bottom-right panel of Fig. 10 reveals that some tracers have been misordered. For instance, there is a lone orange tracer at (x, y) ≈ (15, −4) which should have been assigned a γ value with the opposite sign (e.g. γ ≈ 1). This is an artifact of our incomplete training set for γ init , which is constructed from the nearest neighbors search and inevitably leaves a population of tracers disconnected from the graph.
We address these misclassifications by reinforcing the same momentum condition used in the nearest neighbors search, though in the context of neural network training. That is, we smooth out the ordering of stream tracers inferred by the Interpolation network by encouraging the momentum condition that dx(γ i )/dγ is parallel to the direction of motion of tracers, v i / v i . This is implemented by connecting the Interpolation network to Param-Net through the shared γ-parameter ( Fig. 10; diagonal arrow), reducing two independent neural networks down to one through an autoencoder structure. γ-values are fine tuned by the simultaneous training of these neural networks, through the loss function where T n and T θ are tangent vectors, defined in §2.2, and λ is a hyper-parameter which we set to ∼ 10 2 . Adjusting the γ values inferred by the Interpolation neural network through this loss function, which simultaneously tunes the parameters of Param-Net, leads to a more accurate, smooth ordering of stream tracers illustrated in the bottom left panel of Fig. 10. In this panel, the misplaced tracer at (x, y) ≈ (15, −4) in the bottom right panel has been assigned a reasonable γ-value.
Once a smooth ordering of tracers has been achieved (bottom left panel of Fig. 10), γ-values are fixed in the analyis and no further adjustments to the ordering of stream tracers are made. The stream track and track speed neural networks are trained on these fixed γ-values.

B. NEURAL NETWORK ARCHITECTURES
In this section we provide a description of the neural network architectures utilized in this work.

B.1. Stream Track
The stream track neural network is a fully connected multilayer perceptron with a single input node (γ) and a 3-dimensional output representing the cartesian (x, y, z) coordinates of the stream track. We use 3 hidden layers consisting of 100 nodes each, with tanh activation functions. The output node is a linear combination of weights and biases, and is not mapped through a nonlinear activation function.
Rather than training on the raw (x, y, z) coordinates of stream stars, we standardize each dimension of the training data by its mean µ and standard deviation σ. That is, for each spatial dimension x i we scale the training data as follows: We apply the inverse transformation of Eq. B4 (e.g. x scaled i σ i + µ i ) to infer the unscaled stream track. This neural network has an identical structure to Param-Net, previously described in Appendix A.

B.2. Track Speed
The track-speed neural network is a fully connected multilayer perceptron with a single input node (γ) and a scalar output node representing the average local speed of stream stars, v ≡ v . We use 3 hidden layers consisting of 36 nodes each, with tanh activation functions. The output node is a linear combination of weights and biases, and is mapped through an exponential to yield a positive estimate of the mean speed.

B.3. Interpolation Neural Network
The interpolation neural network is utilized to estimate the parameter γ as a function of position x. It is described in Appendix A.2. The interpolation network is a fully connected multilayer perceptron with 6 input nodes representing the phase-space coordinates (x, v), and two output nodes representing the γ-parameter and a stream membership probability, p. We use 3 hidden layers consisting of 100 nodes each, with tanh activation functions. We adopt the convention that γ ∈ [−1, 1], by applying a tanh activation function to the output node corresponding to the γ-parameter. For the probability parameter p ∈ [0, 1], we use a sigmoid activation function.

B.4. Potential Neural Network
The neural network used to estimate the gravitational potential in §5.1 is a fully connected multilayer perceptron with 3 input nodes representing position x ∈ R 3 relative to the galactic center, and one output node representing the gravitational potential. We use 3 hidden layers, consisting of 36 nodes each with tanh activation functions. The output node is also mapped through a tanh activation function, which is then multiplied by a positive amplitude parameter A. The amplitude is also tuned during neural network training. One subtlety of inferring the gravitational potential from accelerations is that the potential Φ and Φ + const. are both equally valid. We therefore impose the condition that at the origin, x = 0, the potential is zero. We implement this by passing the point x = 0 through the potential neural network for each evaluation point, multiplying by the amplitude A, and subtracting the resulting value from the neural network output corresponding to the arbitrary position x. This ensures that at the origin, the estimated potential will always be zero since all outputs are relative to the potential at the origin. In practice, using a tanh activation function scaled by amplitude A for the output node is found to yield faster convergence to the estimated accelerations and more stable training, likely because mapping the output node through the function A tanh (·) reduces the need for particularly large weights and biases to reproduce the estimated accelerations. Similar results can be achieved without mapping the output node through the tanh function, though fractional errors are found to be slightly higher after the same training period compared to the case with the tanh output.

C. EXISTENCE OF A PHYSICAL STREAM TRACK
In this work we estimate a one-dimensional track that characterizes the mean position and dynamics of a stream. The output of our analysis is an acceleration vector sampled along the stream track, namely, a(γ). From this acceleration vector we can compute the derivative da/dγ. However, because we infer what is effectively a "slice" through the full acceleration field, calculating derivatives like ∂a/∂x i in general (for cartesian coordinate x i ) is not possible. It is therefore interesting to consider under what conditions the inferred accelerations are physical. In particular, for a positive mass density we expect ∇ · a ≤ 0. For a conservative gravitational force, we also expect ∇ × a = 0. Indeed, da/dγ (which we can compute directly) is related to the derivatives involved in computing the divergence and cross product above. Namely, for the i th cartesian component of the acceleration a = (a 1 , a 2 , a 3 ). We define the unit vector which is tangent to the stream track, x(γ). Dividing Eq. C5 by dx/dγ , one obtains Because the left-hand-side of Eq. C7 can be calculated explicitly using our method, we have access to the directional derivative of each acceleration component along the stream track. We now consider working in an orthogonal coordinate system for which one axis is parallel to the stream track for each evaluation point γ. A common choice is the Frenet-Serret frame or TNB frame (Stoker 1969), which is defined by the tangent unit vector (T ; defined in Eq. C6), normal unit vector (N ≡ T (γ)/ T (γ) ), and binormal unit vector (B ≡ T × N ) along a parametrized curve in 3dimensions. Collectively, these unit vectors construct an orthogonal right-handed coordinate system. In our context, each unit vector is a function of the scalar parameter γ, such that the orientation of the coordinate system changes along the stream. For the remainder of this section, we assume that ∇ is the usual gradient operator though expressed in the TNB frame. That is, ∇ ≡ T ∂ ∂T +N ∂ ∂N +B ∂ ∂B . In this frame, Eq. C7 reduces to where i can be T, N or B. We therefore cannot compute derivatives with respect to any quantity other than T . This means that for a given stream track, it will always be possible to construct an acceleration field that satisfies This is because we only constrain the term in the divergence containing the variable T . The other two derivatives are not constrained by our analysis, and can always be chosen so as to satisfy Eq. C9.
For the curl of a, in the TNB frame we have Because we only constrain derivatives of the acceleration components with respect to T , one is guaranteed that the remaining derivatives can be chosen to satisfy ∇ × a = 0, provided that the stream track does not form a closed loop (for a closed loop track, Stokes' theorem can be used to enforce the curl-free condition. However, no such stream is expected to exist). We note that from the curl-free assumption, we can compute the derivatives ∂a T /∂B and ∂a N /∂T . However, these derivatives are not necessary in order to compute Eq. C9. We emphasize that when we reconstruct the gravitational potential ( §5) from the estimated accelerations a(γ) for several streams, we are able to impose a positive mass density as a physically motivated prior. That is, we impose ∇ 2 Φ(x) ≥ 0. We are able to calculate such terms since we model the full scalar potential, rather than a curve parameterized through it.