Characterising stellar halo populations I: An extended distribution function for halo K giants
We fit an Extended Distribution Function (EDF) to K giants in the Sloan Extension for Galactic Understanding and Exploration (SEGUE) survey. These stars are detected to radii and span a wide range in [Fe/H]. Our EDF, which depends on [Fe/H] in addition to actions, encodes the entanglement of metallicity with dynamics within the Galaxy’s stellar halo. Our maximum-likelihood fit of the EDF to the data allows us to model the survey’s selection function.
The density profile of the K giants steepens with radius from a slope to at large radii. The halo’s axis ratio increases with radius from 0.7 to almost unity. The metal-rich stars are more tightly confined in action space than the metal-poor stars and form a more flattened structure. A weak metallicity gradient dex/kpc, a small gradient in the dispersion in [Fe/H] of dex/kpc, and a higher degree of radial anistropy in metal-richer stars result. Lognormal components with peaks at and are required to capture the overall metallicity distribution, suggestive of the existence of two populations of K giants.
The spherical anisotropy parameter varies between 0.3 in the inner halo to isotropic in the outer halo. If the Sagittarius stream is included, a very similar model is found but with a stronger degree of radial anisotropy throughout.
keywords:Galaxy: halo - Galaxy: kinematics and dynamics - Galaxy: stellar content - methods: data analysis
About 1 of the stellar mass in the Galaxy comprises the stellar halo, primarily in the form of old and metal-poor stars. Some of these stars may have formed in-situ, but others come from disrupted satellite galaxies and globular clusters (ibata+95; majewski+03; belokurov+07; bullock+05; abadi+06; cooper+11). An essential step towards distinguishing between processes for halo formation is to gain a clear characterisation of the halo’s present chemodynamical structure.
Since the Galaxy’s potential is in an approximately steady state, we can construct steady-state models of the halo by adopting distribution functions (DFs) that depend only on the constants of stellar motion (jeans+16). There are clear advantages in taking these constants of motion to be the actions. First, action coordinates can be complemented by canonically conjugate variables, the angles, to obtain a complete coordinate system for phase space. Second, it is straightforward to add distinct DFs for the thin and thick discs, the bulge and the dark halo to the DF of the stellar halo to build up a complete Galaxy model (binney_piffl2015). Third, actions are adiabatic invariants and therefore can be used to examine phenomena such as adiabatic contraction (e.g. piffl+15). Finally, the actions , and have simple physical interpretations in that they approximately describe excursions of an orbit in the radial, azimuthal and vertical directions.
A DF maps the halo into a distribution of stars in an easily imagined three-dimensional action space. In general we expect the dynamics of a population of stars to depend on their chemistry, so we assume that is a function of [Fe/H] in addition to . That is, we assume that the halo can be well represented by an Extended Distribution Function (EDF) of the type introduced by sanders+15 for the Galaxy’s disc(s).
From an EDF one can predict the velocity distribution of stars of any metallicity at any location, so an EDF provides the power to pull together results from several surveys into a single coherent picture.
bell+08 established that much of the halo is contributed by substructures that are thought to be relics of disrupted satellites and globular clusters. An EDF of the form can represent such structures only to the extent that they are old enough to have become well phase mixed. Moreover, each disruption event will have created a clump of stars in action space, so the halo’s EDF is likely not a smooth function of . Nonetheless, in this paper we investigate the extent to which existing data for the halo can be represented by an EDF which is a smooth function of . If this proves possible, it will suggest that the data do not yet provide sufficient phase-space resolution to demonstrate the existence of clumps in action space. Alternatively, an unsuccessful attempt to fit a smooth EDF to the data will bring into focus action-space structures arising from accretion events.
The stellar halo can be traced through several stellar populations. Blue horizontal-branch (BHB) stars are luminous, primarily metal-poor, and their distances can be related simply to their apparent magnitudes, colours, and metallicities (deason+11; fermani+13a). RR Lyrae stars are also primarily metal-poor and are associated with a period-luminosity function that can be used to determine their distances. Metal-poor K giants, which have high luminosities, provide useful tracers of the outer halo.
Early efforts to determine the halo’s density profile were limited to Galactocentric radii . The halo proved to be oblate with an axis ratio , and its density declined as with (e.g. preston+91; robin+00; yanny+00; juric+08). More recent work, reaching to greater distances, has found evidence for a break in the stellar density profile at kpc. Star counts of main-sequence turn-off stars (bell+08; sesar+11), RR Lyrae stars (watkins+09; sesar+13), BHB stars (deason+11), and K giants (kafle+13) all exhibit a break in the number density around with at and further out. An examination of data for BHBs and blue stragglers from the SDSS by deason+14 found evidence that beyond and at still larger radii. More recently however, xue+15 found that if the halo’s flattening is permitted to decrease with radius in the case of K giants observed in SEGUE, the best-fitting model has at all radii, but a break in the flattening at from axis ratio at to much further out.
gould+03 and smith+09 find that the velocity ellipsoid is generally well aligned with spherical polar coordinates, as it would be if the Galaxy’s potential were spherical.
It is unclear whether the halo rotates (fermani+13a), and if so, whether at the largest radii it rotates in the opposite sense to the disc (gould+03; beers+12; schonrich+11; schonrich+14). Dynamical models in a spherical potential have generally found radial velocity dispersions to be bigger than the tangential dispersions (deason+12a; kafle+14; williams+15b). Some authors have, however, found the opposite in the region kpc (kafle+12; deason+13). If there really is a radial range in which tangential velocity dispersion is dominant, this may reflect the presence of substructures in the halo, especially the Sagittarius stream (chiba+00; king+15).
The existence of a metallicity gradient in the halo is controversial. Early work such as that of zinn+85 found no gradient outside the solar radius, but carollo+07 claim a negative gradient, with the outer halo significantly more metal poor than the local halo. This was later confirmed by several studies (e.g. dejong+10; kafle+13; allende+14; chen+14), with some also finding metal-poorer stars in retrograde motion and metal-richer stars in prograde motion (carollo+07; kafle+13). schonrich+14 argue that metal-poorer stars can be seen at greater distances than metal-richer stars, and if this effect is not correctly included in the adopted selection function, a metallicity gradient can be erroneously inferred. However work by peng+12 examining main-sequence stars found little evidence for a gradient, although their analysis is limited to metallicities above -2.0. The latest study on the topic by xue+15, who take the selection function into account in a sample of K giant stars, does find a modest, but significant metallicity gradient within the outer halo.
To obtain a coherent picture of the current state of the halo, a model needs to be created that simultaneously considers the various constraints on the density, the kinematics, and the chemistry, with a detailed consideration of the various selection effects in play. A forward-modelling approach based on an EDF is best because then one can take full account of the constraints imposed by dynamics and one can model the impact of the selection function on the data. In this paper, we develop an EDF for K giants, which span a range of metallicities.
2 A model stellar halo
An extended distribution function (EDF) gives the probability density of stars in the space spanned by the phase-space coordinates and the variables that characterise stars, such as mass , age , and chemistry ([Fe/H], ). Given the success of the hypothesis of a universal initial mass function (IMF) in astronomy, we assume that the halo’s EDF is simply proportional to this IMF below the mass at which the stellar lifetime becomes equal to , and zero at higher masses. We assume that all halo stars have age (jofre+11), and adopt the IMF of kroupa+93
We use the -enhanced B.A.S.T.I. isochrones of pietrin+06. With these simplifying assumptions, the only explicit dependencies of the EDF are on phase-space coordinates and [Fe/H], so the EDF can be considered either a function or a function .
We take the gravitational potential to be axisymmetric, so the natural choice of phase-space coordinates are cylindrical polar coordinates . For actions we use and . We use the Stäckel Fudge (binney+12) to convert to .
We adopt a left-handed coordinate system in which positive is away from the Galactic centre and positive is in the direction of Galactic rotation. To convert from the Galactocentric coordinates to the heliocentric coordinates we assume that the Sun is located at (schonrich+12; binney+97), that the Local Standard of Rest (LSR) has an azimuthal velocity of , and that the velocity of the Sun relative to the LSR is (schonrich+12).
2.1 An extended distribution function (EDF)
The formulae below are significantly simplified by introducing as a measure of metallicity the variable
is an increasing function of [Fe/H] with a vertical asymptote at . The upper panel of Fig. 1 shows the relationship between and [Fe/H] when .
Our EDF is
where is a normalisation constant enforcing a total probability of one. For , which is approximately the overall metallicity DF, we adopt
This function, which is plotted in the lower panel of Fig. 1 for and , provides a well defined peak at and a tail towards lower metallicities. Although is a Gaussian in , is not a Gaussian in either or [Fe/H]. Several of these components could be superimposed to capture the overal metallicity DF; we start with just one.
The part of the EDF (3) that involves , , is a generalisation of the DF proposed by posti+15 for spheroidal objects. Specifically,
For , is dominated by , which is a homogeneous function of the actions, and for , is dominated by , which would also be a homogeneous function of but for a dependence on the core action, . This dependence ensures that tends to a finite value as (the orbit that involves sitting at the centre of the Galaxy). Outside the tiny core region, ensures that the phase-space density declines as for , while ensures that the phase-space density declines as at much larger actions.
posti+15 showed that DFs of the form (5) generate self-gravitating systems in which the density profile approximates a double power law in radius, so , where the inner and outer values of are related in a simple way to , and , respectively. Thus by making and functions of metallicity we can arrange for the metal-rich stars (for example) to form a more compact body than the metal-poor stars, with the consequence that there is a metallicity gradient within the halo.
Recasting the metallicity dependence of in terms of rather than , gives
Since is an even function of , it does not endow the halo with rotation.
The parameters control the shape of the density distribution and the velocity ellipsoids (posti+15). Rescaling all the and by the same factor will have no effect on the model if it is accompanied by rescaling of and . We eliminate this degeneracy by imposing the conditions .
Reducing or increases the velocity dispersion in the th coordinate direction in the inner/outer halo. Increasing reduces the flattening of the density distribution. The flattening can also be reduced by lowering . Thus the and control both the kinematics and the physical shape of the halo. We could have made them functions of [Fe/H], but in this first attempt at an EDF for the halo we do not do so.
2.2 The gravitational potential
Our potential is based on the functional forms proposed by dehnen+98. It is generated by thin and thick stellar discs, a gas disc, and two spheroids representing the bulge and the dark halo. The densities of the discs are given by
where is the scale length, , is the scale height, and controls the size of the hole at the centre of the disc, which is only non-zero for the gas disc. The densities of the bulge and dark halo are given by
Here sets the density scale, is a scale radius, and the parameter is the axis ratio of the isodensity surfaces. The exponents and control the inner and outer slopes of the radial density profile, and is a truncation radius.
The adopted parameter values are taken from piffl+14 and given in Table 1. The dark halo is assumed to be a spherical NFW halo, and not truncated (). The contribution of the stellar halo, which has negligible mass, can be considered included in the contributions of the bulge and dark halo.
3 Observational constraints
Here, we introduce the observations that we will use to constrain parameters of the EDF. An ideal catalogue would provide Galactic coordinates , parallaxes , proper motions , apparent magnitudes, and colours, from which a photometric metallicity can be derived . If spectra are also available, we have additionally line-of-sight velocities , spectroscopic metallicities , -abundances, effective temperatures , and surface gravities . From isochrones we can relate apparent magnitudes to a heliocentric distance111We shorten ‘heliocentric distance’ to distance for the remainder of the paper. .
3.1 SEGUE data
We take values for the seven components of the vector
from Data Release 9 of the SDSS. These data form part of the Sloan Extension for Galactic Understanding and Exploration survey II (SEGUE II). We do not use SEGUE-I data because the selection function of these data is hard to determine. We use only the K-giant sample of xue+14. We supplement the catalogue values with proper motions downloaded from SkyServer’s CasJobs222http://skyserver.sdss.org/CasJobs/, by cross matching the data sets with the requirement that right ascensions and declinations match within 15 arcsec. We remove stars on cluster and test plates and apply the cut to reduce the contamination from disc stars (schonrich+12). We compile two data sets, one including the Sagittarius stream and one that does not include plates that intersect the two polygons given by fermani+13a as containing the stream. The final data sets contain 1689 and 1096 stars.
|Star type||Proper motions||Apparent||Colour||Metallicity|
3.2 Selection criteria
Our selection functions are intersections of the original SEGUE-II targeting criteria, and further criteria imposed by xue+14 and by us.
The selection on sky positions is given by the coverage of the SEGUE-II plates, and on each plate by the completeness of the spectroscopic sample, i.e., the fraction of available targets that were observed. We assume that the completeness is independent of apparent magnitude, an assumption which we checked for a handful of plates. The remaining selection criteria, relating to proper motions, apparent magnitudes, colours, and metallicity are summarised in Table 2. Grouping by selection criteria results in three subsamples for the K giants. The SEGUE-II I-colour K giants comprises the largest sample, so we model these stars. The phase-space and metallicity distributions of the sample are shown in Figs 2 and 3, with and without the Sagittarius stream. Removing the Sagittarius stream affects the distributions of the stars in sky positions, results in a slightly steeper observed density radial profile, and shifts the distribution of towards more negative values because the Sagittarius stream has a mean (ibata+97).
4 Fitting the data
We determine the best-fit parameters of the EDF using the likelihood method of mcmillan+13, which sanders+15 used to fit an EDF to the Galaxy’s discs. We now discuss the form of the likelihood and the methods used to determine the contributing terms.
4.1 The likelihood from Bayes’ law
The likelihood of a model is given by the product over all stars of the probabilities . The are the probabilities of measuring the star’s catalogued coordinates given the model and that it is in the survey . By Baye’s law this is
where is the number of stars. is the probability that the star is in the survey given the observables i.e. the ‘selection function’. is the EDF convolved with the error distribution of the observables. is the probability that a randomly chosen star in the model enters the catalogue. We maximise the total log-likelihood
We assume a small core action , so is defined by the parameters ().
4.2 Selection function
There is no selection on . A star is selected when its total proper motion lies within the limits specified in Table 2 and is excluded otherwise. Hence the selection function is separable as
4.2.1 Selection on sky positions
The selection on and depends on the coordinates of the SEGUE-II plates (which sample the sky sparsely), the completeness fraction per plate, and whether the Sagittarius stream is excluded. The completeness of a plate depends strongly on because close to the plane available targets are numerous so an individual star has a low probability of acquiring a spectrum.
For coordinates within of the centre of a plate that is not excluded because it intersects the Sagittarius stream, the selection function equals the completeness fraction for that plate, which as stated earlier, we assume to be independent of apparent magnitude. The completeness fractions are given by xue+15. In Fig. 4 the colour scale shows the selection function in sky coordinates. Near the poles the selection function is large but there are not many stars.
4.2.2 Selection on distance and metallicity
We use isochrones to tabulate on a grid in and [Fe/H] the probability that a star randomly chosen at birth passes the sample’s magnitude and colour cuts. Given our assumption that all halo stars have age , this involves integrating the IMF over values of the mass that correspond to the appropriate absolute magnitudes
We use the -enhanced B.A.S.T.I. isochrones for mass-loss parameter (pietrin+06). We store each isochrone in one-dimensional interpolants that predict absolute magnitude and colours given the mass. For given values of and [Fe/H] we randomly choose 5000 masses from a uniform distribution that covers the range of masses that correspond to stars luminous enough to enter the sample, and if a star passes the magnitude and colour cuts, we add its IMF value to .
The colour scale in Fig. 5 shows . At any given value of [Fe/H], first rises with distance beyond the bright apparent magnitude cut-off and then declines sharply when the limiting apparent magnitude is reached. At any value of , declines when [Fe/H] is reduced because the colours predicted by the isochrones start falling out of the allowed colour boxes.
The symbols in Fig. 5 mark the locations of the observed stars, which are largely, but not entirely, confined to regions of non-negligible . The density of stars should, of course, be proportional to the product of the EDF and the selection function integrated through a volume.
4.3 Convolution of the EDF with the error distribution
We neglect errors in sky coordinates and adopt (independent) Gaussian error distributions for , , , [Fe/H], and . Thus the multi-variate error distribution is
Here true values of observables are indicated by primes and
The convolution of the EDF with the error distribution of a given star is
The Jacobian determinant here is proportional to (mcmillan+12). The integral is calculated using a fixed Monte Carlo sample of 2000 points to eliminate the impact of Poisson noise on the likelihood (mcmillan+13). The calculation of the integral for all the observed stars is parallelised using OpenMP over as many cores as are available.
4.4 The normalisation factor
The normalisation of is given by
The selection function could also be convolved with the observed error distribution as some stars will scatter into the observational volume, while others scatter out by observational errors. We do not do this on account of computational constraints. We approximate the integrals over sky coordinates by a sum of the values taken by the remaining five-dimensional integral at the centre of each included SEGUE plate.
Since is multiplied by in equation (11), it must be calculated to a high degree of accuracy (mcmillan+12). We calculate it to an accuracy of 0.1% using the cubature method implemented in the Cuhre routine (hahn05). The package is highly versatile, with a range of both Monte Carlo and deterministic algorithms, and automatic core parallelization.
4.5 Maximising the likelihood
To derive the probability distribution of the model parameters would require a Markov Chain Monte Carlo (MCMC) approach. Standard MCMC provides an estimate of the propagation of random errors and is useful for examining correlations between parameters. Here systematics arising because the true EDF is not smooth, our selection function and potential are imperfect, and the EDF is parameterized, likely dominate the uncertainties. Hence, we restrict ourselves to finding the EDF that maximises the likelihood of the data. We do this using the Nelder-Mead ‘amoeba’ algorithm (nelder+65), subject to the following priors. We impose and , where the angled brackets denote a mean. This ensures a density profile that decreases outwards for all metallicities, has an outer slope that is either the same as or steeper than the inner slope, and corresponds to a finite total mass. We require to drive the fitting procedure towards models in which all three actions play a role in setting the value of the EDF.
As stated earlier, we fix to . On account of the large number of parameters to fit, we seek the maximum-likelihood model in three stages. We start by fitting the base power-law indices ( and ) and the scale action () with the other parameters set to sensible values. Specifically, we set the metallicity indices ( and ) to zero and the lognormal parameters ( and ) to the values obtained by fitting the global metallicity distribution. We set the action weights (, , , and ) to the values given in Col. 1 of Table 3. These parameters produce a flattened density ellipsoid with an orbital structure that becomes more radially anisotropic outwards. The initial amoeba for this first fit is kept relatively large to allow a thorough scan of the three-dimensional parameter space. In the second step, is fixed to the value found in the first step and the other parameters are varied with the exception of the action weights, , , , and . For this step the initial amoeba is narrower in and than in the other parameters. In the final step, we allow the action weights to vary, again with the initial amoeba narrower in parameters already fitted and wide in the action weights. The outputs of each stage are given in Table 3.
To validate the likelihood methodology, we used a sampling-rejection algorithm to draw 1000 actions from the EDF with parameters representing an approximately isotropic, round model (Col. 1 of Table 3). We complement these actions with angles uniformly distributed on and used torus mapping (binney_mcmillan15) to convert to the first six coordinates in equation (9). We then added 5% Gaussian errors to these mock observations, and used them to fit an EDF using the likelihood methodology described above. We repeated this procedure ten times to study how the fitted parameters scattered around the true values. The true values were captured in the range of best-fit parameters derived from the mock observations.
Here we discuss the parameters of the best-fitting EDF, the quality of the fit, and the structure of this EDF in action and phase space.
5.1 The best-fit parameters
In Table 3 we give the best-fit parameters without (Cols. 2–4, 6) and with (Col. 5) the Sagittarius stream. The parameters are similar in the two cases. The halo’s inner region has velocity ellipsoids that are strongly elongated in the radial direction () and such that the azimuthal velocity dispersion is larger than the dispersion perpendicular to the Galactic plane (). The inner halo is flattened as a consequence. The outer halo is a nearly spherical system with almost isotropic velocity ellipsoids (). The halo’s density profile always steepens around the radius corresponding to (). The principal effects of removing the stream are (i) to strengthen the dependence upon [Fe/H] of the inner density profile ( versus ), and (ii) to reduce the radial anisotropy of the outer region ( versus 1.14).
5.2 Fits to the observables
We assess the fit of the model to the observables when the stream is excluded by generating a mock catalogue, using a sampling-rejection method and taking into account the selection function and errors. The sampling-rejection method involves constructing a proposal density that is a function of the coordinates to be sampled. Our proposal density is the product of Gaussian distributions in each of , , , , , , and . We pick a point in the seven-dimensional space by sampling the proposal density, and evaluate the EDF there. We compute the point’s heliocentric coordinates and Gaussian noise is added to them and to . The selection function is evaluated at the point just derived and multiplied by the previously derived value of the EDF. If the ratio of this product to the proposal density is greater than another uniformly generated random number, the mock observables are added to the sample.
The red points joined by lines in Figs. 6 and 7 show histograms of the mock catalogue. The region covered by analogous histograms of 2000 resamplings of the error distributions of the measured observables are shown by the grey regions. In plots for observables with small errors, such as and and , the grey regions form fairly well defined curves. In plots for observables with large errors, namely , the grey regions fill out a region of significant width. In every case but one, the red curve lies near the centre of this swathe, indicating that the EDF and selection function are together doing a good job in the sense that the data would be reasonably likely to occur if the EDF and selection function accurately evaluated the probabilities. The exception to this statement is the metallicity. Although the general shape of the observed metallicity DF is captured, the fit is poor. So we generate a new mock catalogue of by adding a second lognormal component to
where and are as before, , and is associated with a new and . A non-rigorous exploration of the space finds a much improved fit to the observed metallicities for parameters given in Col. 6 of Table 3 and illustrated in Fig. 7. These parameters generate peaks at and . The exercise of fitting the revised EDF with a second Gaussian in ab initio is left for a future paper, in which substructures are sought in chemodynamical space.
Fig. 8 compares histograms for the joint distribution of pairs of phase-space observables , , , etc., when the Sagittarius stream is excluded and the EDF is defined in terms of a single Gaussian in . The colour-filled contours show distributions of mock observables, while the black contours show the distributions of measured observables. In general there is good agreement between the mock and observational distributions.
5.3 The distribution of stars in action space and metallicity
Fig. 9 shows for metallicities , , and plots of the best-fitting along the line in action space given by , thus probing similar halo orbits at different energies. The values of have been scaled to agree at . We see that declines more steeply with for metal-richer stars, implying that these stars are tightly confined near the origin of action space.
Each panel of Fig. 10 shows for a different metallicity the value of in the action-space plane , again in the case the stream is excluded. Contours of constant are approximately straight lines because they are similar to the lines on which either or are constant. At small values of , the slopes of the contours are , and at large they tend to . Comparing the top panel for with the bottom panel for , the closer confinement of metal-rich stars to the origin of action space is evident.
5.4 Moments of the best-fitting DF
Fig. 11 shows the density predicted by the best-fitting EDF when the Sagittarius stream is excluded. The flattening of the halo can be discerned. By fitting ellipses to the isodensity curves in this figure we obtain the radial density profiles shown in the top panel of Fig. 12. The full black line shows the profile obtained by excluding the Sagittarius stream, while the dashed black line shows the profile obtained by including the stream. The difference between these two profiles is barely discernible. The red dotted line shows the radial profiles of the most metal-rich stars (), while the green dotted curve shows that of the most metal-poor stars (). Again the increase in the compactness of the halo as [Fe/H] is increased is evident. The bottom panel of Fig. 12 shows the axis ratios of the entire halo (black curves) and its mono-abundance components (coloured dotted curves). The inner halo is more flattened () than the outer halo (), and flattening decreases as one moves from the most metal-rich to the most metal-poor component. The central panel of Fig. 12 shows plots of logarithmic density gradients. The gradient of the overall profile (full curve) steepens from in the inner halo to in the outer halo.
We characterise the mixed velocity moments by calculating the tilt in the velocity ellipsoid using (smith+09)
where and is the angle between the major axis of the projection of the velocity ellipsoid in the coordinate plane and a coordinate direction. We find the tilt to be close to zero for all projections, which means the velocity ellispoids all point towards the centre.
Fig. 13 shows that contours of constant are mildly elongated, and those of constant are strongly elongated along the axis. By contrast, contours of constant are flattened. The fourth panel of Fig. 13 plots the spherical velocity-dispersion anisotropy parameter
This is generally positive, implying radial bias, with the highest values attained at on the equatorial plane. The anisotropy varies between moderate radial bias in the inner halo () to near isotropy () in the outer halo. To compare anisotropy profiles when the Sagittarius stream is excluded and included, we calculated a one-dimensional spherical anisotropy profile by averaging over bins in ellipsoidal radii. Fig. 14 shows the resulting profiles of for the whole halo (black lines) and for the halo binned by metallicity (dotted coloured lines). Anisotropy tends to decrease outwards. Including the Sagittarius stream (broken black line) changes to roughly constant between 0.28 and 0.35 throughout. Since the metal-rich stars are concentrated at low actions where the anisotropy is largest, at any given radius anisotropy increases with metallicity, being everywhere below at .
The top panel of Fig. 15 shows the mean metallicity as a function of . Surfaces of constant are moderately flattened, and declines outwards from near the centre to at the largest radii, the typical metallicity gradient being dex/kpc. The lower panel of Fig. 15 shows the dispersion in , which increases outwards from to . Hence at any radius the range of metallicities represented spans the range covered by the mean metallicity over a full .
6.1 Comparison with previous work
6.1.1 Radial density profile
The density of stars decreases with radius with a power-law slope of in the inner halo, steepening to in the outer halo. The halo’s axis ratio increases from in the inner halo to in the outer halo. The findings do not depend on whether we include or omit the Sagittarius stream, and are consistent with the profile previously determined for K giants (xue+15). BHBs also have a density profile that becomes steeper with radius but the logarithmic slope is claimed to vary much more rapidly (deason+11), with potentially a third steeper power-law index required to describe the density profile beyond (deason+14). We cannot comment on the sharpness of the transition from the inner to the outer slope in the K giants because this is fixed by our adopted form of EDF. By adding an additional parameter to the EDF, it would be straightforward to vary the sharpness of the transition analogously to how zhao+96 generalised the nfw+96 profile, but we have not explored this avenue.
carollo+07 inferred the global structure of the halo from the kinematics of stars now located within of the Sun. They argued that the sample divided into “inner” and “outer” populations. The MDF of the inner population peaks at , while that of the outer population peaks at . From local kinematics of a larger sample carollo+10 inferred that the inner population has a steeper, more flattened density profile (, ) than the outer population (, ). Our EDF is consistent with these findings as regards the peak metallicities and flattenings of the two populations. However, both the samples and the analysis method of carollo+07 and carollo+10 differ sharply from ours. Their dynamics was restricted to a spherical potential, and their samples were local and contained a mix of stellar types. schonrich+14 argue that their evidence for dichotomy in the halo is an artefact arising from their failure to model their sample’s selection function.
6.1.2 The velocity ellipsoid
We find the velocity ellipsoids to be aligned with spherical polar coordinates, which may be the consequence of assuming a spherical dark matter halo. Previous studies have found spherical alignment (e.g. gould+03; smith+09), but not in the context of a fully dynamical model.
The velocity ellipsoid is radially biased in the inner halo, but when the Sagittarius stream is excluded, we find it to be nearly isotropic in the outer halo. Studies of BHBs have generally found the velocity ellipsoid to be more radially biased (deason+12a; williams+15b), although from the same sample of BHBs hattori+13 reached a similar conclusion to our finding for the K giants. carollo+10 also find the outer-halo stars to be on less eccentric orbits.
Simulations predict a higher degree of radial anisotropy for stars that have been accreted by the halo through minor mergers (bullock+05; abadi+06). The simulations may be unreliable, but it is also true that the anisotropy in the outer halo is more uncertain than in the inner halo because stars in the outer halo are less well represented in the data and their tangential components of velocities are very uncertain.
6.1.3 Metallicity components
Two lognormal components are required to capture adequately the metallicity DF, even when the dependence of the selection function on metallicity is taken into account. An ad-hoc exploration of the parameters of a DF defined by two Gaussians finds a metal-richer component peaking at , with an upper limit at , and a metal-poorer component peaking at , with an upper limit at . carollo+07, chen+14, an+15, and xue+15 reached similar conclusions, based on very different samples and technqiues. Of these, only an+15 and xue+15 consider selection effects.
6.1.4 Dependence of phase-space structure on metallicity
The metal-rich stars are more tightly confined to the region around the origin of action space than the metal-poor stars. As a consequence, the radial density profile of the metal-rich stars is steeper than that of the metal-poor stars, leading to a weak metallicity gradient of dex/kpc. This result is consistent with the analysis of K giants by xue+15, but takes the work further by demonstrating dynamical consistency between the density and dispersion profiles in a realistic Galactic potential. The literature regarding determination of gradients in other stellar populations (e.g. carollo+07; peng+12; chen+14; allende+14; schonrich+14; an+15) shows no consensus on the existence of a gradient. The lack of consensus may arise from improper treatment of the selection function in some cases, with metal-poorer stars being preferentially detected at larger distances (schonrich+14). In other cases, it may reflect genuine differences between stellar populations, again emphasising the need to employ EDFs, and to extend them to cover other characteristics of stars.
A gradient arises naturally from metallicity gradients in the satellites digested by the Galaxy: the outer, more metal-poor stars are stripped at large radii, and the more metal-rich core at smaller radii. Direct evidence for this effect is provided by the metallicity gradient along the Sagittarius stream (chou+07). Moreover, dynamical friction drags the most massive and metal-rich satellites inwards, while the less massive ones experience less friction and are tidally shredded at large radii.
The dispersion in [Fe/H] increases with radius, indicating that a wider range of metallicities is contributing to the outer halo. Since metal-rich stars can be found only in more massive satellites, this finding suggests the outer halo is produced by stripping a more heterogeneous group of satellites than produced the inner halo. However, we cannot discount the possibility that some of the metal-richer stars now found at large radii have been scattered onto energetic orbits by the Galaxy’s bar.
At all metallicities, the radial bias of the velocity ellipsoids decreases outwards. Since the metal-richer stars reside deeper in action space, they have a more radially biased velocity ellipsoid than the metal-poorer stars. This result mirrors conclusions drawn by carollo+10 from a sample of mixed stellar types, and by hattori+13 from BHBs.
6.2 Uncertainties in the analysis
6.2.1 Impact of substructure
Removing the Sagittarius stream diminishes the radial bias of the velocities in the halo. We do not know how substructures that we have not masked impact our assessment of the halo’s structure. However, the ability to reproduce adequately the phase-space observations after excluding the stream implies that the current data for K giants is too sparse in phase space to resolve the halo’s substructure. Indeed, the latter has only been detected in counts of dwarf stars, which are much more numerous than K giants. It seems that our EDF gives a satisfactory account of the structure of the halo on the medium to large scales that can currently be probed with K giants.
6.2.2 Stellar population assumptions
Our evaluation of the metallicity-distance selection function depended on the assumption of a single (old) age, and on relations from isochrones between the age, mass, and metallicity of a star and its luminosity in various wavebands. Systematic errors arising from faulty isochrones are difficult to assess. However, the assumption of a single age for the K giants is likely to have the following impact. Isochrones for older stars result in very similar selection functions. However if there is a significant proportion of younger stars, the true selection function would be higher than the one we used at the metal-richer and shorter-distance end. Using our selection function will lead to an EDF that over-estimates the number of metal-poor and distant stars. The ability of our EDF to reproduce the density profile for K giants found by xue+15 suggests that our selection function is not materially in error.
6.2.3 Parameterised form of the EDF
The imposition of a functional form on the EDF can bias the results by restricting the set of possible solutions. We believe the phase-space distribution function is sufficiently general, allowing a range of density and anisotropy profiles, but the dependence of the EDF on metallicity is rather limited in that only the power-law indices and depend on [Fe/H], not the scale action or the parameters and that determine the shape of the velocity ellipsoid and the halo’s flattening. It is clear however that a metallicity gradient is not forced by the EDF; if none were needed by the data, the indices would have been found to be independent of metallicity.
6.2.4 Adoption of a fixed potential
With sufficiently perfect data for a stellar population, one should be able to determine the potential, because an EDF can perfectly model the data only in the true potential. We chose a specific potential at the outset and sought an EDF that reproduced the observations in this potential, and imperfections in our choice of potential will both bias our EDF away from the true EDF and give rise to discrepancies between our best model and the data. Our success in reproducing the data suggests that our chosen potential is not seriously in error. The main weakness of our work probably lies in the small number of very distant stars in the sample and the large uncertainties in their tangential velocities.
7 Conclusions and further work
We used a sample of I-color K giants from the SEGUE-II survey to probe the chemodynamical structure of our galaxy’s stellar halo. To these data we fitted an EDF that provides for a density profile that gradually steepens with radius from a shallow asymptotic slope at small radii to a steeper slope at large radii. These slopes can vary with [Fe/H]. Using an EDF enables one to build a dynamically consistent model that allows the phase-space distribution of stars to vary with metallicity, whilst incorporating a realistic selection function that takes into account restrictions on sky positions, proper motions, apparent magnitudes, and colours. We find that our models reproduce the observations well, presumably because the data for K giants are not yet rich enough to resolve most halo substructures.
We find that the power-law slope of the density profile of the K giants steepens from to in the outer halo. The halo is flattened to axis ratio on the inside but becomes almost spherical on the outside. The overall metallicity DF is best captured with two lognormal distributions, peaking at and . The metal-rich stars are more tightly confined in action space than the metal-poor stars, with the consequence that the halo has a weak metallicity gradient. The dispersion in metallicity at any point increases outwards. The metal-rich stars form a more flattened body than the metal-poor ones. The extent of velocity anisotropy depends on whether the Sagittarius stream is included or not. With the stream included, the velocity ellipsoid is everywhere moderately radially biased, while when it is excluded, moderate radial bias in the inner halo gives way to isotropy at large radii.
There are several possible directions for further work. The EDF could be applied to detect substructures in a denser sample of K giants in the phase-space-metallicity domain, thus complementing the work of janesh+16, who look for substructures in SEGUE K giants using sky positions, heliocentric distance, and radial velocity. The EDF could be elaborated to include rotation and to allow for explicit dependence on [Fe/H] of the parameters , that jointly control the system’s flattening and velocity anisotropy. The EDF could be changed to make the transition between the inner and outer asymptotic slopes of the density profile sharper. The EDF could be elaborated to include dependence on . The EDF could be fitted to a larger sample of K giants and to other halo tracers, such as blue horizontal branch stars – we hope to report on these exercises shortly.
The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement no. 321067. The work was also supported through grant ST/K00106X/1 by the UK’s Science and Technology Facilities Council.
We thank GitHub for providing free private repositories for educational use, thus enabling seamless version control and collaboration. PD thanks Jason Sanders, Wyn Evans, Gus Williams, and Eugene Vasilyev for useful conversations, and Jason Sanders, Paul McMillan and Til Piffl for routines. Finally, PD would like to thank Xiang-Xiang Xue for providing the K giant catalogues and completeness fractions for the SEGUE-II plates.