Better Galactic Mass Models through Chemistry

With the upcoming release of the Gaia catalog and the many multiplexed spectroscopic surveys on the horizon, we are rapidly moving into a new data-driven era in the study of the Milky Way’s stellar halo. When combined, these data sets will give us a many-dimensional view of stars in accreted structures in the halo that includes both dynamical information about their orbits and chemical information about their formation histories. Using simulated data from the state-of-the-art Latte simulations of Milky-Way-like galaxies, which include hydrodynamics, feedback, and chemical evolution in a cosmological setting, we demonstrate that while dynamical information alone can be used to constrain models of the Galactic mass distribution in the halo, including the extra dimensions provided by chemical abundances can improve these constraints as well as assist in untangling different accreted components.


Introduction
Our knowledge of the Galaxy in which we live, the Milky Way (MW), is poised to undergo a revolution in the next ten years thanks to a new generation of state-of-the-art stellar surveys.These include photometric surveys like PanSTARRS [1] and LSST [2], spectroscopic surveys like 4MOST [3], DESI [4], WEAVE [5], and Subaru-PFS [6], and astrometric surveys starting with Gaia [7] and extending to LSST and WFIRST [8].The Galactic renaissance will also include far greater insight into the chemical abundances of stars thanks to efforts like the Cannon [9], which can translate the abundance patterns generated from smaller high-resolution spectroscopic surveys like APOGEE [10], GALAH [11], and the Gaia-ESO survey [12] into the larger medium-and low-resolution spectroscopic surveys listed above.It is not an exaggeration to say that ten years from now we can expect to have complete phase-space information for stars in the Galaxy nearly to its virial radius, complemented by 10-20 dimensions of chemical abundance information.This new high-dimensional view of the Galaxy demands new approaches to understanding its contents.Here we discuss how this combined phase and abundance space will be uniquely powerful for setting constraints on the MW's dark matter (DM) distribution and untangling the building blocks of its accreted stellar halo, allowing us to use our Galaxy as a time machine to study ancient dwarf galaxies.
Currently our knowledge of the distribution of stars in the MW is limited to the region within about 120 kpc [13], less than half of what we imagine to be the full extent of its dark halo.Because of this limited information, our knowledge of the shape and mass of the MW's DM halo is equally limited: measurements of its total mass disagree by a factor of ∼4 [14] and our understanding of its shape and radial profile is still poorly constrained enough to be controversial [15,16] especially in light of new estimates of the mass of its largest companion, the Large Magellanic Cloud [17,18].This makes it difficult to place our Galaxy in a cosmological context in order to use it as a test of DM theories, since many predictions from cosmological simulations scale with mass and concentration [19].

Methods
Stars in the accreted stellar halo of the Milky Way provide one of the few sources of information about the outer extents of the Galactic dark halo.Most of these stars are part of structures called tidal streams (e.g., [20]) created by the disruption of ancient dwarf satellite galaxies of the MW, and are expected to reach to the MW's virial radius and beyond [21,22].Since the stars in tidal streams all have similar orbits, tidal streams are, in principle, more sensitive to the radial profile and shape of the halo than the orbits of individual bound dwarf galaxies or globular clusters, and so are valuable for constraining the MW's shape and mass.Furthermore the approximately isotropic nature of accretion (although see [23]) samples orbits well out of the plane of the Galactic disk and at much larger distances, giving us a more complete map of the gravitational forces exerted by our Galaxy than we could otherwise achieve.

The Action Space of Stellar Orbits
One approach to mapping the Milky Way's dark halo with tidal streams is to work in the space of their actions J, which are defined for stars on bound orbits in a system described by a Hamiltonian H as a canonical transformation of their positions x and momenta p, such that for each component x i , where the integral is over one cycle of the coordinate x i .When defined in this way, a star's actions J are adiabatic invariants (i.e., integrals of the motion) and the angles θ evolve linearly in time with frequencies Ω: For more background on action-angle variables we refer the reader to the discussions in Chapter 3 of [24] or Chapter 10 of [25]; for more detailed information on the techniques described in this section see [26,27].
For a tidal stream, in which all the stars began life within the small (compared to the MW) phase-space volume of a dwarf galaxy, their actions are persistently similar unless the host galaxy (in this case the MW) undergoes a major merger that non-adiabatically transforms the gravitational potential.The actions thus conveniently label stars in each tidal stream by their orbit family, and the small cluster in J representing the stream can be treated using linear transformations [28].Furthermore, a stellar halo made of accreted small galaxies, as predicted by theories of cold DM, will exhibit a high degree of clustering in the space of its actions even for structures that are spread across the Galaxy (Figure 1, panels A and B).This space is thus ideal for separating different accreted structures from one another [29][30][31].The same mock stellar halo shown in a projection of action space, calculated using the best-fit spherical NFW approximation for the gravitational potential of the dark halo.The stars in each progenitor are now clustered together around the parent orbits of the accreted galaxies they formed from.Panel C: The same mock stellar halo with actions calculated using a potential that is roughly 1σ away from the best-fit approximation.The clustering is notably less pronounced in this incorrect action space.Figure adapted from [27].
Studying tidal streams in action space requires assuming a global gravitational potential, but the structure of action space itself gives insight into the correct potential to use.As shown by comparing panels B and C of Figure 1, only when the assumed potential is a close representation of the true mass distribution will the action space of the stellar demonstrate the expected degree of clustering [26,33,34].A statistical measure of the degree of clustering can thus be used as the figure of merit to determine the best-fit potential [26,27,35].We will describe the quantification of clustering in Section 2.2.The advantage of this method of determining the Galactic gravitational potential, especially in the context of the many upcoming stellar surveys, is twofold.First, it does not require that individual streams are identified or stars assigned to membership in a particular stream, merely that the stars in the fitting sample come mainly from the accreted stellar halo so that the assumption of an intrinsically clustered distribution is satisfied.Second, the method is agnostic as to which labels are used as dimensions, so long as the stars in each stream demonstrate some correlation with each other in those labels and at least some of them depend on the gravitational potential.This means that if chemical abundance information is also available for the stars in the fitting sample, these labels can be added as extra dimensions.The chemical abundances of stars from a given progenitor galaxy should be correlated through their shared star formation history, so this extra information should improve the ability to constrain the gravitational potential by, for example, differentiating stars from progenitors with different masses that accreted on similar orbits.

Measuring Clustering in Actions
To quantify the degree of clustering in the distribution of actions J and abundances Z, given a guess for the gravitational potential Φ, we use the mutual information (MI) statistic [36], defined as the relative entropy between the target distribution p Φ (J, Z) and a "shuffled" version pΦ (J, Z): The distribution pΦ is constructed by randomizing the J and Z coordinates of each star relative to one another.This transformation preserves the one-dimensional marginal distributions along each coordinate, and is equivalent to the product of the marginals: To calculate the MI for a given set of stars, we first calculate the stellar actions assuming the potential Φ and use the density estimation code EnLink [37] to calculate the densities p Φ and pΦ in (J, Z) space at the location of each star.We then use Monte Carlo integration to estimate the MI by summation over all stars in the sample: As shown in Figure 2, the more clustered the distribution is, the greater the difference with the product of its marginals, and hence the larger the value of MI(Φ).The maximum value of MI(Φ) identifies the most clustered distribution and hence the best-fit gravitational potential.
Figure 2. Measuring the degree of clustering using the mutual information.All panels show the distribution of energies and angular momenta for a region of a simulated stellar halo from one of the Latte galaxies [38].The colors denote log projected density, with dark red denoting low density regions and light yellow regions of high density.The side-plots show the one-dimensional projected density along each dimension.The top row shows a case with a high degree of mutual information (MI, Equation ( 4)), where the distribution obtained by randomizing the coordinates of the stars while retaining the one-dimensional projected distributions (right) differs greatly from the original correlated distribution (left).The bottom row shows a case with a low degree of mutual information, where the distribution is similar whether or not the coordinates of the stars are shuffled.The mutual information can thus be used to measure the degree of clustering in the space of actions or constants of motion for a set of stars.

Simulating the Accreted Stellar Halo
To test the hypothesis that adding chemical abundance information can improve constraints on the Galactic gravitational potential, we used the simulated galaxy m12i from the Latte simulation suite, described in [38].This MW-like simulated galaxy is the product of a high-resolution, cosmological zoom-in simulation that uses the FIRE-2 feedback recipe [39] to follow star formation and the evolution of stellar populations.The simulation tracks the evolution of ten chemical elements produced by supernovae (types I and II) and winds from evolved stars.For the examples shown here we used the second-highest resolution version of the simulation, in which star particles have m p ∼56,500 M .For this preliminary test, we selected star particles in a region between 50 and 100 kpc from the main galaxy's center to avoid the galactic stellar and gas disk.
Using this selection as the fitting sample, we modeled the galactic mass distribution with a spherical Navarro-Frenk-White form [40] with scale radius r s and scale density ρ s , parameterized in terms of the enclosed mass M s at the scale radius, In terms of M s and r s , the potential is where G is Newton's constant.The choice of M s over ρ s as the mass parameter is motivated by the sensitivity of stellar orbits to the enclosed mass rather than the local density, and avoids severe degeneracies between parameters during the fit.This functional form ignores the triaxiality of the halo but can adequately approximate the radial profile of the gravitational potential over a limited radial range.To check the results of the fit using the halo sample, we also fit the same functional form directly to the gravitational potentials of the stars as calculated directly by the simulation (Figure 3).over the range shaded in gray.The density of stars is higher at smaller r within the range being fit, so agreement at smaller r is prioritized (consistent with the behavior of the action-space fitting scheme).

Results
Calculating the MI for a distribution involves estimating the density from a finite number of points, so adding extra dimensions is costly and can introduce noise in the MI value and hence in the best fit.So before adding all ten chemical abundances tracked by the simulation, we first looked at the change in the MI when adding each abundance separately to the actions, while assuming the potential fit directly to the values calculated from the simulation and shown in Figure 3.A larger increase in MI when a given abundance is added signals a larger degree of clustering in that dimension and should correspond to more additional constraining power on the potential when added to the actions.The left panel of Figure 4 shows that indeed some abundances appear to contain more extra information than others, probably as a result of the assumptions about how different elements are produced in the simulation.Guided by this result we selected nitrogen, calcium, and iron to use as additional labels in fitting the gravitational potential.We caution that the result that these three abundances were most informative should not be over-interpreted, given the simplified yield tables used in the simulations; this is more a proof of concept and a demonstration of a technique that could be used on future data sets.Next we carried out a fit of the gravitational potential by varying the parameters M s and r s of the model over an iteratively refined grid in parameter space, first using only the actions and then adding the extra chemical abundance dimensions one at a time.We then looked at the contrast between the maximum value of the MI (at the best-fit potential) and the value far from the best-fit peak, and likewise at the value of the MI for the parameters obtained when the NFW function was fit to the potential values calculated directly from the simulation (i.e., Figure 3).We compare to this baseline MI far from the best fit because the MI is an entropy and therefore an extensive quantity, meaning that its absolute value will increase somewhat when a dimension is added even if that dimension is not informative.However, as is evident from the right panel of Figure 4, the difference between the peak and baseline MI values increases as new abundance dimensions are introduced, suggesting that the height of the best-fit peak is increasing (and therefore the confidence interval around the best fit should shrink).
Once the best-fit potential is obtained, the different building blocks of the accreted stellar halo can be separated from one another by running a group-finding algorithm on the same space of actions and abundances.Here again one expects that adding abundance information to constants-of-motion space can help distinguish the different components.Figure 5 illustrates how groups found in constants of motion space alone using EnLink (left) indeed represent different tidal streams (center) but also have distinct, though overlapping, abundance distributions (right).This bodes well for future attempts to reconstruct the component dwarf galaxies of the Milky Way's stellar halo.[37] in the space of constants of motion for star particles with 50 < d < 300 kpc in a high-resolution simulated Milky Way (m p = 7070 M ) from the Latte suite, including a more realistic model for the turbulent diffusion of chemical elements [39] that results in more plausible abundance distributions.Different groups are plotted in different colors.Center: the same groups shown in Cartesian coordinates (with origin at the center of the main simulated galaxy) show the features of tidal streams.Right: a selection of 5 groups shows the variety of different distributions in iron abundance among the different building blocks of the stellar halo.

Discussion
These results, obtained using the star particle data from a realistic simulation of a MW-like galaxy, indicate that studying tidal streams in the combined space of orbits and abundances is a promising direction both for constraining the Galactic gravitational potential and reconstructing the destroyed galaxies of the halo.Most interestingly, these results were obtained with a simulation in which the diffusion of chemical elements into star-forming gas was not implemented, which produces broader abundance distributions compared to those shown in Figure 5, and broader than the known abundance distributions of bound MW dwarf galaxies (Escala et al., in prep).Thus these results can be understood as a conservative estimate for the extra information contributed by elemental abundances.
On the other hand, these results assumed perfect data without observational uncertainties; furthermore, those uncertainties depend strongly on the type of stellar tracer for both phase-space coordinates and abundances.For example, variable-star standard candles like RR Lyra stars will have extremely precise distances [13,41] but pose difficulties when obtaining precise radial velocities or elemental abundances thanks to the same physics that makes them standard candles.Conversely, bright, cool giant stars will be visible to great distances and have well-constrained abundances and RVs, but less precise distances.Even more fundamentally, different stellar tracers have different specific frequencies depending on the underlying metal content of the population; for example, RR Lyr are more common in metal-poor populations and M giants are more common in metal-rich populations.As a result, they probe different regions of the Galaxy, and different regions in the stellar halo in particular, as illustrated in the left panel of Figure 6.
To better account for both the variety of stellar tracers and the range of expected observational errors, we are now updating the synthetic survey code Galaxia [42] to self-consistently resample the high-resolution Latte simulations with synthetic stars drawn from PARSEC isochrones [43].A preliminary example is shown in the right panel of Figure 6.When complete, this machinery will allow us to make robust, detailed comparisons between the Milky Way and simulated galaxies, and accurately forecast the performance of various methods to constrain the Galactic gravitational potential and untangle the accretion history of the stellar halo.

Figure 1 .
Figure 1.Panel A: tidal streams in an accreted mock stellar halo from[32] shown distributed on the sky in Galactic coordinates.Each different color represents stars from a different progenitor dwarf galaxy.Panel B: The same mock stellar halo shown in a projection of action space, calculated using the best-fit spherical NFW approximation for the gravitational potential of the dark halo.The stars in each progenitor are now clustered together around the parent orbits of the accreted galaxies they formed from.Panel C: The same mock stellar halo with actions calculated using a potential that is roughly 1σ away from the best-fit approximation.The clustering is notably less pronounced in this incorrect action space.Figure adapted from[27].

Figure 3 .
Figure 3. Gravitational potential Φ * of star particles in the simulated galaxy calculated directly from the simulation (blue) as a function of radius r, with the best direct fit of the NFW functional form (red) over the range shaded in gray.The density of stars is higher at smaller r within the range being fit, so agreement at smaller r is prioritized (consistent with the behavior of the action-space fitting scheme).

Figure 4 .
Figure 4. Left: Mutual information of stars in the Latte halo sample, in the space of constants of motion plus one additional label as shown on the horizontal axis.The solid black line shows the value of the mutual information (MI) calculated for constants of motion alone.The elements nitrogen, calcium, and iron appear to be most strongly correlated with the orbits of stars in the tidal streams.Right: MI values at the best-fit (blue), direct-fit (green), and far from the best fit (red) parameters for the gravitational potential, for the Latte stellar halo sample.When additional abundances (Ca, Fe, N) are included as extra dimensions, the contrast between the peak and the background increases relative to the fit when only actions are used ("dyn only").

Figure 5 .
Figure5.Left: groups found by EnLink group finder[37] in the space of constants of motion for star particles with 50 < d < 300 kpc in a high-resolution simulated Milky Way (m p = 7070 M ) from the Latte suite, including a more realistic model for the turbulent diffusion of chemical elements[39] that results in more plausible abundance distributions.Different groups are plotted in different colors.Center: the same groups shown in Cartesian coordinates (with origin at the center of the main simulated galaxy) show the features of tidal streams.Right: a selection of 5 groups shows the variety of different distributions in iron abundance among the different building blocks of the stellar halo.

Figure 6 .
Figure 6.Left: Cumulative distance distributions for two stellar tracers, M giants (red) and RR Lyr (cyan), generated from the simulated Latte halo used in this work using a modification of the code Galaxia [42].While the RR Lyr closely trace the distribution of star particles from the simulation (black), the M giants have significantly different, far more centrally concentrated distribution.Right: Sky distribution of the synthetic RR Lyr survey generated from the Latte halo.The color scale indicates Galactocentric distance.