Building an Acceleration Ladder with Tidal Streams and Pulsar Timing

We analyze stellar streams in action-angle coordinates combined with recent local direct acceleration measurements to provide joint constraints on the potential of our Galaxy. Our stream analysis uses the Kullback-Leibler divergence with a likelihood analysis based on the two-point correlation function. We provide joint constraints from pulsar accelerations and stellar streams for local and global parameters that describe the potential of the Milky Way (MW). Our goal is to build an ``acceleration ladder", where direct acceleration measurements that are currently limited in dynamic range are combined with indirect techniques that can access a much larger volume of the MW. To constrain the MW potential with stellar streams, we consider the Palomar 5, Orphan, Nyx, Helmi and GD1 streams. Of the potential models that we have considered here, the preferred potential for the streams is a two-component Staeckel potential. We also compare the vertical accelerations from stellar streams and pulsar timing, defining a function $f(z) = \alpha_{1pulsar}z - \frac{\partial\Phi}{\partial z}$, where $\Phi$ is the MW potential determined from stellar streams, and $\alpha_{1~\rm pulsar}z$ is the vertical acceleration determined from pulsar timing observations. Our analysis indicates that the Oort limit determined from streams is consistently (regardless of the choice of potential) lower than that determined from pulsar timing observations. The calibration we have derived here may be used to correct the estimate of the acceleration from stellar streams.

Recently, extreme-precision time-series observations to directly measure Galactic accelerations and the Galactic potential have become feasible. Chakrabarti et al. (2021a) analyzed compiled pulsar timing observations to measure the Galactic acceleration and derive fundamental Galactic parameters, including the Oort limit (the mid-plane density), the local dark matter density, and the shape of the Galactic potential. In the near future, we can expect direct acceleration measurements from extreme-precision radial velocity observations ) and from eclipse timing (Chakrabarti et al. 2021b). Although direct acceleration measurements have provided local constraints on the potential (Chakrabarti et al. 2021a), they do not yet provide constraints for the MW mass or other global parameters for which we require information over a large volume. While stellar streams provide indirect constraints, they add valuable information about the shape and extent of the MW potential at larger distances. Here, we explore the idea that combining complementary information from local acceleration measurements and tidal streams can be used to build an "acceleration ladder" to derive constraints on the Galactic potential.
Several methods have been used to extract information about the Galactic potential from tidal streams. One is to fit an orbit to a stream based on the observed positions and velocities along the stream, as done in Newberg et al. (2010) and Koposov et al. (2010). Agreement between the stream orbit and the data is potential dependent, allowing a fitting procedure to estimate parameters for the MW potential. There are known biases in this methodology that result from differences between the orbits of stream stars and the progenitor (Sanders & Binney 2013). An alternative method is to back-integrate the orbits for the stream member stars in the MW potential (Price-Whelan et al. 2014). If these orbits are calculated using an accurate MW potential, then the member stars should become bound to the progenitor at some point in their orbital history. Another common technique is a forward modelling method, such as the one described in Bonaca et al. (2014), which has several related variations (Varghese et al. 2011;Sanders 2014;Fardal et al. 2015). This method uses a Monte Carlo Markov Chain algorithm that compares simulated streams to observed streams in 6D phase space.
We focus here on stellar streams rather than gaseous streams to take advantage of the vast infrastructure that has been developed for the collisionless components of galaxies (Binney & Tremaine 2008), and the better distance estimates that are available for stars in stellar streams relative to the kinematic distance estimates that are typically available for the gas (Levine et al. 2006). We employ the action space clustering method developed in Sanderson et al. (2015), and later applied in Reino et al. (2021). In the true potential, stream stars that were once part of a disrupted dwarf galaxy are expected to be on similar orbits, and thus have similar positions in action space. We apply this by searching for the potentials that lead to the strongest clustering in the actions of stream stars. Since actions are integrals over one period of the orbit, this method is similar to the "rewinder" algorithm Price-Whelan & Johnston (2013). The analysis of stellar streams by Sanderson et al. (2017) employs a "consensus fit" whereby multi-ple streams provide information on various regions of the MW potential, allowing the potential model to be valid over larger parts of the galaxy; another advantage of this approach is that the consideration of multiple streams leads to a smaller bias than one stream alone (Bonaca et al. 2014) (i.e., due to the disrupted dwarf galaxy being near pericenter, which can bias the results (Reino et al. 2022))It has been demonstrated based on the Aquarius simulations that these methods are effective at recovering the present day properties of the host galaxy potential (Sanderson et al. 2017). A limitation of these techniques is that they do not extend well to models for the potential that are not in equilibrium or are time dependent, so these effects cannot be captured. Our goal is to recover the present day properties of the potential, as this is what the direct accelerations probe, so this method works well for our purposes.
One of the challenges to using tidal streams for measuring the Galactic potential is obtaining 6D phase space information for the member stars. High quality positions and proper motions can typically be obtained from Gaia data (Gaia Collaboration et al. 2016, significantly expanding the sample of stars with available information (Bonaca et al. 2014). To get to full 6D information we also require distances and radial velocities. For some streams, missing measurements can be estimated by fitting the distance or radial velocity along the stream track, then interpolating based on the position in the stream. The same method is used in Reino et al. (2021) for the streams considered here, and for Sagittarius in Vasiliev et al. (2021).
One of our goals is to use recent direct acceleration measurements from pulsar timing to calibrate indirect measurements for the potential of the MW. This would in effect serve to build an "acceleration ladder", similar to the distance ladder that is based on a calibration of indirect distance measures to the basic geometric measure of distance from parallax measurements. Direct acceleration measurements provide the first rung of the ladder or a "pivot point" for the Galactic potential, which can be used to calibrate stream-based potential measurements. In Reino et al. (2021), only models that matched the circular velocity at the Solar circle to within uncertainties were considered, but the measurements of the circular velocity are indirect, while our pulsar accelerations provide a direct measurement of the potential. The paper is organized as follows. We present our methodology in §2 and discuss our data selection in §3. We present our main results in §4, and conclude in §5.

METHODS
We use the python package galpy (Bovy 2015) to calculate the actions for our streams. They are calculated analytically for the Isochrone potential (Binney & Tremaine 2008), and numerically in all other potentials considered. In spherically symmetric potentials, like the Hernquist potential, we use galpy's spherical approximation method, providing accurate and efficient action calculations. Staeckel potentials allow us to use the Staeckel approximations, providing accurate actions requiring only a single numerical integral for each action (Bonaca et al. 2014). For all other potentials, we use the Isochrone approximation method in galpy. This method uses an auxiliary Isochrone potential with integrated orbits of the stars in the desired potential to compute the actions. This requires orbit integration for each star, and is therefore relatively computationally intensive. However, this provides more accurate approximations to the actions in these potentials compared to other methods. In particular, the assumptions made for the other methods break down when the radial or vertical actions are on the same order as L z (Bovy 2014).
For each potential model, we tested the accuracy of the action approximations by calculating the actions of all our stream members along their orbits for a variety of different points in the parameter space. We use a leap-frog algorithm to compute the orbit for at least 1 orbital period, and then calculate the actions at regular intervals along the orbit. The actions should remain constant, and we find that the action variations are consistently less than 1% in J R and J Z . When we use data from multiple streams, we shift the actions in L z for each individual stream such that there is limited overlap between streams. This will not impact the clustering in action space, but helps to avoid error modes with large masses and low scale lengths (Reino et al. 2021).
As we intend to build upon the results from direct acceleration measurements from pulsar timing, we consider the potentials used in Chakrabarti et al. (2021a). Specifically, we consider the α 1 potential and the αγ potential. These are defined in Equations 1 and 2 respectively, where α 1 is the inverse square of the frequency of low-amplitude vertical oscillations, V LSR is the local standard of rest velocity, and γ describes the shape of the potential. The α 1 potential provided the best fit to the pulsar accelerations, and the αγ potential is the "cross-term" model discussed in Chakrabarti et al. (2021a). These models provide a good match to the measured pulsar accelerations, however they do not produce strong clustering in action space. This is expected because these models were designed to produce reasonable orbits for stars in the Solar neighborhood with low eccentricities and inclinations, which is not gen-erally satisfied by our stream members (Quillen et al. 2020).

Kullback-Leibler Divergence
We use the Kullback-Leibler Divergence (KLD) to measure the amount of clustering in action space, as in (Sanderson et al. 2015). The KLD compares two PDFs to each other. We use it to measure the amount of action space clustering by comparing the action distribution with a uniform distribution over the maximum range for each action across all the datasets (potentials and streams). Larger KLD values indicate a larger difference between the distributions, implying that there is more clustering in action space. Equation 3 defines the KLD for a continuous random variable from p(x) to q(x), where in our case p(x) will be the action distribution and q(x) will be our uniform distribution.
We use the kernel density estimator enbid (Sharma & Steinmetz 2011) to estimate the PDF for the actions. We define the resulting distribution for a potential with parameters a as f a (J), and then we can calculate the KLD from equation 4. This integral can be re-written as in (Reino et al. 2021) using a Monte Carlo approximation to sum across points drawn from the distribution. This is a natural choice in our case where we have a discreet set of actions corresponding to our member stars. The exact form of the calculation used corresponds to the wKLD1 given in Reino et al. (2021), which includes a weight term such that the streams contribute equally in joint fits.
For our data set, we use an optimization algorithm to find the maximum value of the KLD within reasonable bounds. For 2 parameter potentials, we apply a simple grid search across the parameter space. The same grid can be used for the error analysis technique described in Reino et al. (2021) and Sanderson et al. (2015), which measures the KLD between the best fitting action distribution and the distributions in other potentials. For potentials with more than 2 parameters, we use a differential evolution optimizer implemented in scipy (Virtanen et al. 2020) to maximize the KLD.

Likelihoods and Error Analysis
We use the two-point correlation function to compute the likelihood from the action distributions in order to compare with the pulsar sample. We use the likelihoods from both methods and combine them to get a joint likelihood function. The likelihoods for the streams are calculated using the methods given in (Yang et al. 2020). This approach is different from that used in Reino et al. (2021), chosen here because it enables the combination of our methods.
We can use the likelihoods to estimate the uncertainty by examining the relative likelihood (i.e., the likelihood divided by the maximum likelihood value) for the parameters that we are interested in. A confidence region in that parameter space can be defined by setting a threshold on the relative likelihoods corresponding to the desired confidence interval. We generally focus on the total mass and Oort limit for the various potentials, but also compute a surface across all parameters to measure their uncertainties. We considered using Fischer matrix analysis to compute the uncertainties. However, this analysis underestimates the uncertainties for the streams. The likelihood surfaces can have multiple local maxima near the maximum likelihood, corresponding to the best fits of individual streams. Since the Fisher matrix analysis is based on derivatives of the likelihood function, it is only sensitive to regions near the maximum likelihood, leading to an underestimation of the uncertainties.
The combined sample (i.e., for pulsars and streams) has a likelihood surface defined by the product of the two individual likelihoods, which gives us the relative likelihoods for the combined sample. In this way we can easily obtain the uncertainties for both the individual and combined samples using the same mechanism.
For all the potentials considered except for the α 1 potential, we examine the total mass and Oort limit for the models. In cases with more than 2 parameters we marginalize over the remaining parameters to produce a two-dimensional likelihood surface. We expect the streams to produce better constraints on the total mass, while the pulsars should provide stronger constraints on the Oort limit. Thus, we expect the combined sample to provide tighter constraints than either individual sample.
Note that we do not use a maximum likelihood analysis to find the best fitting parameters, instead maximizing the KLD as described above. surfaces we have confirmed that the best fitting KLD parameters closely agree with the maximum value of our likelihood surface, so we retain the best-fit identified using the method of (Reino et al. 2021).
The equations from Yang et al. (2020) used to calculate the likelihood are given in Equations 5 and 6. Here D is a distance between particles in action space normalized by the standard deviations of all the actions. P(ln D) is a probability distribution describing the distance between stars in action space, which is computed using equation 5 and the two-point correlation function. We set a maximum value of D max , as at large distances the behavior will change due to an insufficient number of pairs. We have considered five different tidal streams -GD1, Orphan, Palomar 5, Nyx (Necib et al. 2020) and the Helmi stream (Helmi et al. 1999). We find that for some potentials the Helmi stream produces large uncertainties, and is not included in our primary results. Recent work finds that the Helmi stream provides a constraint for the shape of the inner dark matter halo of the Milky Way (Dodd et al. 2022), noting that this stream may be in a resonance. However, if it is truly in a resonance, calculating actions for typical potentials is not valid, which is consistent with our findings here. Meanwhile, Zucker et al. (2021) find that the abundances of likely Nyx stream members are consistent with the thick disk, and may not have an extragalactic origin. It does, however, remain plausible that this is the result of an early minor dwarf galaxy merge (Wang et al. 2022). Therefore our results are calculated using the Nyx, GD1, Palomar 5 and Orphan Streams. However, to avoid contamination if the Nyx stream is actually not extragalactic in origin, all calculations are also performed without these stars included. The rest of our data set is similar to that in (Reino et al. 2021), to reproduce the stream sample and analysis in Reino et al. (2021) as closely as possible. We have produced results in our primary potential model that includes the Nyx stream as well. In this potential, the Nyx Stream does not make a large impact on the derived potential parameters, but does allow for phase space overlap between samples. This overlap allows for better calibrations between the streams and the pulsar accelerations. The Galactocentric R and Z values of the stream member stars can be seen in Figure 1.
The Pal5 stream here is split into two separate samples, a main stream group and a set of 11 stars that appear to be slightly off of the main stream track. To check for any potential biases induced by the inclusion of these sources, we have tracked them in the action calculations, and found that they are not more likely to be outliers in action space than the other stars. This can be seen from the action distributions shown in Figure 2. We have also considered the two separate branches that appear in the Orphan stream data set. All of the analysis that follows has been performed using each branch separately in addition to the full sample. Changing the sample does slightly alter the resulting potential parameters, but at a much smaller level than the estimated uncertainties. These changes do not impact the main results of this paper, so we have retained the full sample here.

RESULTS
Our main results are based on a combination of multiple streams. A single stream will often produce unrealistic potential parameters, while the constraints derived from a combination of streams generates a much more realistic potential (Bonaca et al. 2014 Table 1. The right panel shows a distribution that differs from the best-fit potential by roughly 2σ, based on our error analysis. As before, the red diamonds indicate the locations of the members of Pal5 that appear to be off the main stream. Interestingly, these sources are not the main outliers from the central Pal5 clump in action space. 2018; Reino et al. 2022). We have calculated the best fitting potential for each of our streams individually, as well as with a combination of streams. This is repeated across a range of different potential models, as seen in Table 1. From these results we consider the potential model giving the largest KLD value to be our best fitting potential, in this case a two-component Staeckel potential. This can then be compared with the properties of the potentials that were found to be a good fit with the pulsar sample. Shown in Figure 2 are the radial and azimuthal actions of the Palomar 5, GD1, Nyx and Orphan streams in a two-component Staeckel potential, in both the best-fit potential and a poorly fitting potential.
The results for our combined sample can be seen in Figure 3, which shows the likelihood surface for the Oort limit and MW mass using streams and pulsar timing. We consider both the α 1 potential, which is the best fit to the pulsars, and the two-component Staeckel potential, which is the preferred potential model for the streams. Notes. KLD optimization results for a number of different potentials from the GD1, Palomar 5 and Orphan streams. A two-component Staeckel model provides most clustered actions with a KLD of 1.51. Most of these results lack the Nyx stream in order to minimize potential contamination. The two-component Staeckel model includes this data set however, as it is our best fit stream result without the Nyx stream, and the results are only slightly altered as a result of the inclusion of this data set. This allows for potential constraints with phase space overlap with the pulsar sample. In this case we can see the general trends in the Oort limit that appear across all of our models, where the streams have a lower value for the Oort limit along with higher uncertainties.
Interestingly, the streams lead to a lower Oort limit than the pulsars in every models considered. This is not necessarily surprising as the streams probe distances out to tens of kpc, while the Oort limit is sensitive to the scale height of the disk. Additionally, Staeckel models do not typically produce a thin disk as the measured pulsar accelerations indicate (Batsleer & Dejonghe 1994). They are capable of producing models with Oort limits that match the pulsar value however. The pulsars primarily constrain the acceleration in the vertical direction, so we compare the accelerations from the two methods in the vertical direction. We can write the difference in the accelerations as the function specified in Equation 7. A plot for the vertical accelerations of our preferred models can be seen in Figure 4, where α 1,pulsar = 4073 +1422 −837 Gyr −2 is the value of α 1 that is derived from measured pulsar accelerations, and Φ is the best-fit potential for the streams. Consistent with the lower Oort limit values, the stream samples yield smaller vertical accelerations.
5. CONCLUSIONS • We analyze clustering in action space to determine the best-fit potential for a set of stellar streams that have also been analyzed by Reino et al. (2021), with the addi-  tion of the Nyx stream. We compare the derived fundamental parameters that describe our Galaxy from stellar streams to recent direct acceleration measurements from pulsar timing (Chakrabarti et al. 2021a).
• We consider several potential models, including Hernquist potentials, Isochrone potentials, Staeckel potentials and the potentials used in analyzing the pulsar sample. We analyze Staeckel potentials with both one and two components. The preferred potential for the streams is a two-component Staeckel model, which is the same model considered in Reino et al. (2021) and produces the most clustered actions for our streams.
• We focus primarily here on the mass of the Galaxy and the Oort limit, where the streams provide stronger constraints on the mass and the pulsars provide optimal constraints on the Oort limit. We find that the best-fit potentials for stellar streams underestimate the Oort limit compared to the value derived from measured accelerations from pulsar timing. For the potential masses in the two-component Staeckel model, we obtain M streams = 2.15 ± 0.75 × 10 12 M , M pulsars = 1.95 ± 1.3 × 10 12 M and M combined = 1.98 ± 0.65 × 10 12 M . Our estimated Oort limits in this potential are ρ streams = 0.022 ± 0.02M pc −3 and ρ pulsars = 0.036 ± 0.011M pc −3 . These values can be compared to the Oort limit of the pulsar's preferred potential determined from direct acceleration measurements, which is ρ combined = 0.08 +0.05 −0.02 M pc −3 . The pulsar sample provides better constraints since these local accelerations are highly dependent on the local density, which drives the vertical component of the acceleration.
• The Oort limit obtained from the joint constraints of streams and pulsars is lower than the Oort limit derived from pulsars alone, with a value of 0.028±0.008M pc −3 . However, pulsars alone do not currently constrain the mass of the potential, since the pulsars are distributed within ∼ 1 kpc of the Sun.
• We provide a fitting formula that may be used to calibrate vertical accelerations of stellar streams to the measured pulsar accelerations.