Characterizing stellar halo populations II: The age gradient in blue horizontal-branch stars
MMNRAS , 1–18 (2015) Preprint 11 November 2018 Compiled using MNRAS L A TEX style file v3.0
Characterizing stellar halo populations II: The agegradient in blue horizontal-branch stars
Payel Das (cid:63) , Angus Williams , and James Binney Rudolf Peierls Centre for Theoretical Physics, University of Oxford, OX1 3NP, UK Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK
11 November 2018
ABSTRACT
The distribution of Milky Way halo blue horizontal-branch (BHB) stars is examinedusing action-based extended distribution functions (EDFs) that describe the locationsof stars in phase space, metallicity, and age.The parameters of the EDFs are fitted using stars observed in the Sloan Exten-sion for Galactic Understanding and Exploration-II (SEGUE-II) survey that trace thephase-space kinematics and chemistry out to ∼
70 kpc. A maximum a posteriori proba-bility (MAP) estimate method and a Markov Chain Monte Carlo method are applied,taking into account the selection function in positions, distance, and metallicity for thesurvey. The best-fit EDF declines with actions less steeply at actions characteristic ofthe inner halo than at the larger actions characteristic of the outer halo, and older agesare found at smaller actions than at larger actions. In real space, the radial densityprofile steepens smoothly from − ∼ − ∼ . β s varyingfrom isotropic to between ∼ . ∼ . − .
03 Gyr kpc − , with someuncertainty in the distribution of the larger ages. These results are consistent with ascenario in which older, larger systems contribute to the inner halo, whilst the outerhalo is primarily comprised of younger, smaller systems. Key words:
Galaxy: halo - Galaxy: kinematics and dynamics - Galaxy: stellar content- methods: data analysis
If we assume that the stellar halo is in an approximatelysteady state, we can characterise it with distribution func-tions (DFs) f ( J ) that depend only on the constants of stel-lar motion J i (Jeans 1916). Using actions as the constantsof motion has several clear advantages. First, action coordi-nates can be complemented by canonically conjugate vari-ables, the angles, to obtain a complete coordinate systemfor phase space. Second, it is straightforward to add DFs forthe thin and thick discs, the bulge and the dark halo to theDF of the stellar halo to build up a complete Galaxy model(Piffl et al. 2015). Third, actions are adiabatic invariants andtherefore can be used to examine phenomena such as adia-batic contraction (e.g. Piffl et al. 2015). Finally, the actions J r , J φ and J z quantify excursions of an orbit in the radial,azimuthal and vertical directions, and thus are a natural setof labels for categorizing orbits. (cid:63) E-mail:[email protected]
Bell et al. (2008) established that much of the halo iscomprised of substructures, thought to be relics of disruptedsatellites and globular clusters. These clumps disperse in po-sitions and velocities but their ages and metallicities remainbunched (Bell et al. 2010). We expect each clump of halostars sharing the same metallicity and age to have its ownDF f ( J , θ, [Fe / H] , τ ). Over time as the clump phase mixes,its DF becomes a function of actions, metallicity, and ageonly, i.e. an extended distribution function (EDF) of thetype introduced by Sanders & Binney (2015) for the MilkyWay disc(s), and developed by Das & Binney (2016) forthe Milky Way halo K giants. The EDF of the entire halois simply the sum of these individual EDFs. The ages andmetallicities of stars are thus associated with separationsin action space that may manifest as gradients in real andvelocity space.Several methods have been explored in the literaturefor determining the distribution of ages of halo stars. Theseinclude estimating the main-sequence turn-off temperatureand combining it with metallicities and isochrones, finding c (cid:13) a r X i v : . [ a s t r o - ph . GA ] A ug Payel Das . ± . . ± . B − V increases outwardsto 40 kpc. Interpreting this as an age gradient amounts toa spread of roughly 2–2 . An extended distribution function (EDF) gives the probabil-ity density of stars in the space specified by the phase-spacecoordinates ( x , v ) and the variables that characterise stellarproperties, such as mass m , age τ , and chemistry ([Fe/H],[ α/ Fe] , . . . ). Below the mass at which the stellar lifetime be-comes equal to τ , we assume that the halo’s EDF is propor-tional to the Kroupa IMF (Kroupa et al. 1993), (cid:15) ( m ) = . m − . if 0 . (cid:54) m < . . m − . if 0 . (cid:54) m < . . m − . if m (cid:62) . . (1)where m is in solar masses. The EDF vanishes at highermasses. Thus the only explicit dependencies of the EDF areon phase-space coordinates, [Fe/H], and age, and the EDFcan be considered either a function f ( x , v , [Fe / H] , τ ) or afunction f ( J , [Fe / H] , τ ). For actions we use J r , J φ ≡ L z and J z . Where required, we use the St¨ackel Fudge (Binney 2012)to convert ( x , v ) to J . We take the gravitational potentialto be axisymmetric, and thus cylindrical polar coordinates( R, φ, z, v R , v φ , v z ) are a natural choice.Below we introduce the EDFs and the gravitational po-tential. The former includes two ‘pseudo’ EDFs that dependonly on positions and metallicities. We include these both as an intermediate step to constructing the full phase-spaceEDFs and as a means to comparing with past works thatonly fit density profiles. They are equivalent to full phase-space EDFs integrated over velocities. The EDF is specified asd N d x d[Fe / H] = f ( x , G ) = Af bpl ( R, z ) f m ( G ) , (2)where A is a normalization constant enforcing a total prob-ability of one and G = − ln([Fe / H] max − [Fe / H]) . (3)The function f bpl is the density profile and is specified by abroken power law and a constant axis ratio (Deason et al.2012) f bpl ( R, z ) = (cid:16) R + z /q r (cid:17) − α in / , if R + z /q (cid:54) r (cid:16) R + z /q r (cid:17) − α out / , otherwise (4)We consider the metallicity to be a lognormal in [Fe / H] (Das& Binney 2016) f m ( G ) = e G e − G σ σ √ π , (5)where σ = − ln([Fe / H] max − [Fe / H] peak ) . (6)Thus, the distribution is specified by a maximum metallicity,and a metallicity at which the distribution peaks. The differ-ence between these metallicities cannot be greater than 1. Aglance at the distribution of metallicities in the observations(Fig. 1) suggests this is suitable.Model 1 is thus specified by the parameter set M ( q, α in , α out , r b , [Fe / H] max , [Fe / H] peak ) (7) The EDF is specified asd N d x d[Fe / H] = f ( x , G ) = Af spl ( R, z ) f m ( G ) , (8)where A is a normalization constant enforcing a total prob-ability of one. G and f m are defined by Equations (3) and(5), respectively. f spl is the density profile and is specifiedby a single power law with a variable axis ratio (Xue et al.2015) f spl ( R, z ) = ( R + z /q ( r ) ) − α/ , (9)where r = ( R + z ) − / q ( r ) = q ∞ − ( q ∞ − q ) exp (cid:32) − (cid:112) r + r r (cid:33) . (10)Model 2 is thus specified by the parameter set M ( q ∞ , q , α, r , [Fe / H] max , [Fe / H] peak ) (11) MNRAS000
Bell et al. (2008) established that much of the halo iscomprised of substructures, thought to be relics of disruptedsatellites and globular clusters. These clumps disperse in po-sitions and velocities but their ages and metallicities remainbunched (Bell et al. 2010). We expect each clump of halostars sharing the same metallicity and age to have its ownDF f ( J , θ, [Fe / H] , τ ). Over time as the clump phase mixes,its DF becomes a function of actions, metallicity, and ageonly, i.e. an extended distribution function (EDF) of thetype introduced by Sanders & Binney (2015) for the MilkyWay disc(s), and developed by Das & Binney (2016) forthe Milky Way halo K giants. The EDF of the entire halois simply the sum of these individual EDFs. The ages andmetallicities of stars are thus associated with separationsin action space that may manifest as gradients in real andvelocity space.Several methods have been explored in the literaturefor determining the distribution of ages of halo stars. Theseinclude estimating the main-sequence turn-off temperatureand combining it with metallicities and isochrones, finding c (cid:13) a r X i v : . [ a s t r o - ph . GA ] A ug Payel Das . ± . . ± . B − V increases outwardsto 40 kpc. Interpreting this as an age gradient amounts toa spread of roughly 2–2 . An extended distribution function (EDF) gives the probabil-ity density of stars in the space specified by the phase-spacecoordinates ( x , v ) and the variables that characterise stellarproperties, such as mass m , age τ , and chemistry ([Fe/H],[ α/ Fe] , . . . ). Below the mass at which the stellar lifetime be-comes equal to τ , we assume that the halo’s EDF is propor-tional to the Kroupa IMF (Kroupa et al. 1993), (cid:15) ( m ) = . m − . if 0 . (cid:54) m < . . m − . if 0 . (cid:54) m < . . m − . if m (cid:62) . . (1)where m is in solar masses. The EDF vanishes at highermasses. Thus the only explicit dependencies of the EDF areon phase-space coordinates, [Fe/H], and age, and the EDFcan be considered either a function f ( x , v , [Fe / H] , τ ) or afunction f ( J , [Fe / H] , τ ). For actions we use J r , J φ ≡ L z and J z . Where required, we use the St¨ackel Fudge (Binney 2012)to convert ( x , v ) to J . We take the gravitational potentialto be axisymmetric, and thus cylindrical polar coordinates( R, φ, z, v R , v φ , v z ) are a natural choice.Below we introduce the EDFs and the gravitational po-tential. The former includes two ‘pseudo’ EDFs that dependonly on positions and metallicities. We include these both as an intermediate step to constructing the full phase-spaceEDFs and as a means to comparing with past works thatonly fit density profiles. They are equivalent to full phase-space EDFs integrated over velocities. The EDF is specified asd N d x d[Fe / H] = f ( x , G ) = Af bpl ( R, z ) f m ( G ) , (2)where A is a normalization constant enforcing a total prob-ability of one and G = − ln([Fe / H] max − [Fe / H]) . (3)The function f bpl is the density profile and is specified by abroken power law and a constant axis ratio (Deason et al.2012) f bpl ( R, z ) = (cid:16) R + z /q r (cid:17) − α in / , if R + z /q (cid:54) r (cid:16) R + z /q r (cid:17) − α out / , otherwise (4)We consider the metallicity to be a lognormal in [Fe / H] (Das& Binney 2016) f m ( G ) = e G e − G σ σ √ π , (5)where σ = − ln([Fe / H] max − [Fe / H] peak ) . (6)Thus, the distribution is specified by a maximum metallicity,and a metallicity at which the distribution peaks. The differ-ence between these metallicities cannot be greater than 1. Aglance at the distribution of metallicities in the observations(Fig. 1) suggests this is suitable.Model 1 is thus specified by the parameter set M ( q, α in , α out , r b , [Fe / H] max , [Fe / H] peak ) (7) The EDF is specified asd N d x d[Fe / H] = f ( x , G ) = Af spl ( R, z ) f m ( G ) , (8)where A is a normalization constant enforcing a total prob-ability of one. G and f m are defined by Equations (3) and(5), respectively. f spl is the density profile and is specifiedby a single power law with a variable axis ratio (Xue et al.2015) f spl ( R, z ) = ( R + z /q ( r ) ) − α/ , (9)where r = ( R + z ) − / q ( r ) = q ∞ − ( q ∞ − q ) exp (cid:32) − (cid:112) r + r r (cid:33) . (10)Model 2 is thus specified by the parameter set M ( q ∞ , q , α, r , [Fe / H] max , [Fe / H] peak ) (11) MNRAS000 , 1–18 (2015) he age gradient in halo BHBs The EDF is specified asd N d x d v d[Fe / H] = f ( J , G ) = Af ps ( J ) f m ( G ) , (12)where A is a normalization constant enforcing a total prob-ability of one. f m is given by Equation (5), and f ps is fromPosti et al. (2015) f ps ( J ) = [1 + J /h ( J )] β in [1 + g ( J ) /J ] β out h ( J ) = a r J r + a φ | J φ | + a z J z g ( J ) = b r J r + b φ | J φ | + b z J z . (13)For | J | (cid:29) J , f ps is dominated by g ( J ), and for | J | (cid:28) J , f ps is dominated by h ( J ). Both g ( J ) and h ( J ) are homogeneousfunctions of the actions of degree one.The parameters ( a r , a φ , a z , b r , b φ , b z ) control the shapeof the density and velocity ellipsoids. Rescaling the a i and b i by the same factor has no effect on the model if accompa-nied by a rescaling of J . This degeneracy is eliminated byimposing the conditions (cid:80) i a i = (cid:80) i b i = 3.Model 3 is thus specified by M ( β in , β out , J , a r , a φ , b r , b φ , [Fe / H] max , [Fe / H] peak ) . (14) The EDF is specified asd N d x d v d[Fe / H] d τ = f ( J , G, τ )= Af psr ( J ) f m ( G ) f psa ( J , τ ) , (15)where A and f m are as above. f psr is given by f psr ( J ) = R ( J φ ) f ps ( J ) R ( J φ ) = 1 + x tanh (cid:18) J φ J (cid:19) , (16)where f ps is given by Equation (13). The prefactor R ( J φ )splits the phase-space DF into even and odd components, in-troducing the possibility for rotation. x governs the strengthof the rotation. f psa is given by f psa = δ (cid:18) τ − (cid:20) a τ + b τ ln J t J (cid:21)(cid:19) , (17)where J t is the total action J t = (cid:113) J r + J φ + J z . (18)This implies a single age at each total action, which is aguide to apocentric radius. b τ encodes the dependence onactions, and a τ is the age for which the total action is equalto the transition action J . Increasing | b τ | increases the agegradient within the halo, with b τ < a τ makes the halo olderat every radius.Model 4 is specified by M ( β in , β out , J , a r , a φ , b r , b φ , x, [Fe / H] max , [Fe / H] peak , a τ , b τ ) . (19) Table 1.
Parameters of the Galactic potential.Component Parameter ValueThin R d (kpc) 2.682 z d (kpc) 0.196Σ d ( M (cid:12) kpc − ) 5.707 × Thick R d (kpc) 2.682 z d (kpc) 0.701Σ d ( M (cid:12) kpc − ) 2.510 × Gas R d (kpc) 5.365 z d (kpc) 0.040Σ d ( M (cid:12) kpc − ) 9.451 × R hole (kpc) 4.000Bulge ρ ( M (cid:12) kpc − ) 9.490 × q γ δ r (kpc) 0.075 r t (kpc) 2.100Dark halo ρ ( M (cid:12) kpc − ) 1.815 × q γ δ r (kpc) 14.434 r t (kpc) ∞ As in previous works, we use the composite potential pro-posed by Dehnen & Binney (1998), generated by thin andthick stellar discs, a gas disc, and two spheroids representingthe bulge and the dark halo. The densities of the discs aregiven by ρ d ( R, z ) = Σ z d exp (cid:20) − (cid:18) RR d + | z | z d + R hole R (cid:19)(cid:21) , (20)where R d is the scale length, z d , is the scale height, and R hole controls the size of the hole at the centre of the disc,which is only non-zero for the gas disc. The densities of thebulge and dark halo are given by ρ ( R, z ) = ρ (1 + m ) ( γ − δ ) m γ exp (cid:2) − ( mr /r t ) (cid:3) , (21)where m ( R, z ) = (cid:112) ( R/r ) + ( z/qr ) . (22) ρ sets the density scale, r is a scale radius, and the pa-rameter q is the axis ratio of the isodensity surfaces. Theexponents γ and δ control the inner and outer slopes of theradial density profile, and r t is a truncation radius.The adopted parameter values are taken from Piffl et al.(2014) and given in Table 1. They specify a spherical NFWhalo that is not truncated ( r t = ∞ ). The stellar halo con-tributes only negligible mass, and thus can be consideredincluded in the contributions of the bulge and dark halo. Here, we introduce the observations that we will use to con-strain parameters of the EDFs. We adopt a left-handed coor-
MNRAS , 1–18 (2015)
Payel Das l ( o ) P r o b a b ili t y d e n s i t y −50 0 50 b ( o ) P r o b a b ili t y d e n s i t y
13 14 15 16 17 18 19 20 21 g (mag) P r o b a b ili t y d e n s i t y
13 14 15 16 17 18 19 20 21 r (mag) P r o b a b ili t y d e n s i t y −400 −200 0 200 400 v || (km/s) P r o b a b ili t y d e n s i t y −15 −10 −5 0 5 10 15 µ ∗l (mas/year) P r o b a b ili t y d e n s i t y −15 −10 −5 0 5 10 15 µ b (mas/year) P r o b a b ili t y d e n s i t y −3.0−2.8−2.6−2.4−2.2−2.0−1.8−1.6−1.4 [Fe/H] (dex) P r o b a b ili t y d e n s i t y Figure 1.
One-dimensional distributions for the SEGUE-II BHBs in Galactic coordinates, for sky positions, apparent magnitudes,line-of-sight velocities, proper motions, and metallicities. dinate system in which positive v R is away from the Galacticcentre and positive v φ is in the direction of Galactic rota-tion. To convert from Galactocentric coordinates to helio-centric coordinates we assume that the Sun is located at( R , z ) = (8 . , . − , and that the velocity of the Sun rel-ative to the LSR is ( v R , v φ , v z ) = ( − . , . , .
25) km s − (Sch¨onrich 2012). We use the BHB sample of Xue et al. (2011), which is basedon Data Release 10 of the Sloan Digital Sky Survey (SDSS)and in particular the Sloan Extension for Galactic Under-standing and Exploration surveys (SEGUE) within it. Thesample consists of equatorial coordinates ( α, δ ), apparentmagnitudes, colours, line-of-sight velocities v || , and spec-troscopic metallicities [Fe / H]. We supplement the cataloguevalues with proper motions ( µ ∗ l = ˙ l cos b and µ b = ˙ b ) down-loaded from SkyServer’s CasJobs , by cross matching within15 arcsec. We include proper motion measurements to useall available data, but note that the uncertainties are of-ten > / H] (cid:54) − . http://skyserver.sdss.org/CasJobs/ colours, and metallicity to a heliocentric distance s for theBHBs. Thus our observables are defined by the vector u = ( l, b, s, v || , µ ∗ l , µ b , [Fe / H]) . (23) Our selection function is the overlap between the originalspectroscopic targeting criteria, criteria imposed by Xueet al. (2011), and further criteria applied by us. The selectionon sky positions is given by the coverage of the spectroscopicplates on the sky, and on each plate by the completeness ofthe spectroscopic sample, i.e. the number of stars observedcompared to the potential number of targets. The remainingcombined selection criteria relate to apparent magnitudesand various colour indices. These are summarised in Ta-ble 2. The phase-space, colour, and metallicity distributionsfor the sample are shown in Fig. 1.
We construct the likelihood of models by the method ofMcMillan & Binney (2013) that Das & Binney (2016) usedto fit an EDF to halo K giants. The form and contributingterms of the likelihood are described in detail here.
The total likelihood L of a model M is given by the productover all stars i of the individual likelihoods L i = P ( u i | SM )of measuring the star’s catalogued coordinates u i given themodel M and that it is in the survey S . By Bayes’ law thisis L i = P ( S | u i ) P ( u i | M ) P ( S | M ) . (24) We shorten ‘heliocentric distance’ to distance for the remainderof the paper. MNRAS000
The total likelihood L of a model M is given by the productover all stars i of the individual likelihoods L i = P ( u i | SM )of measuring the star’s catalogued coordinates u i given themodel M and that it is in the survey S . By Bayes’ law thisis L i = P ( S | u i ) P ( u i | M ) P ( S | M ) . (24) We shorten ‘heliocentric distance’ to distance for the remainderof the paper. MNRAS000 , 1–18 (2015) he age gradient in halo BHBs Table 2.
Combined selection criteria in terms of apparent magnitude, colour indices, and metallicity. n ∗ is the number of stars afterapplying the combined selection criteria, and u , g , r , and i refer to apparent magnitudes in Sloan’s ugriz colour-magnitude system.Programme n ∗ Apparent Colour MetallicitymagnitudeSEGUE II 701 15 . < g < . . < ( u − g ) < . < − . r > . − . < ( g − r ) < . l( o ) −50050 b ( o ) Figure 2.
Selection function as a function of sky positions, p ( S | l, b ). P ( S | u i ) is the probability that the star is in the survey giventhe observables i.e. the ‘selection function’ (Section 4.2). P ( u i | M ) is the EDF convolved with the error distributionof the observables (Section 4.3). P ( S | M ) is the probabilitythat a randomly chosen star in the model enters the cata-logue (Section 4.4). The total log-likelihood islog L = n ∗ (cid:88) i = k log L i = n ∗ (cid:88) i = k log (cid:16) P ( S | u i ) (cid:17) + n ∗ (cid:88) i = k log (cid:16) P ( u i | M ) (cid:17) − n ∗ log ( P ( S | M )) , (25)where n ∗ is the number of stars. There is no selection on line-of-sight velocities or proper mo-tions. Moreover, we assume that the selection function isseparable as P ( S | u i ) = p ( S | l, b, s, [Fe / H])= p ( S | l, b ) p ( S | s, [Fe / H]) . (26)The selections on sky positions p ( S | l, b ) and dis-tance/metallicity p ( S | s, [Fe / H]) are described below.
The selection on l and b , p ( S | l, b ), depends on the coordi-nates of the SEGUE-II plates and the completeness frac-tion, i.e. the fraction of photometrically identified targetsfor which spectra were obtained. The completeness fraction depends strongly on | b | because close to the plane availabletargets are numerous so an individual star has a low proba-bility of being allocated a fibre. For coordinates within 1 . ◦ of the centre of a plate, the selection function equals thecompleteness fraction for that plate. Thus P ( S | l, b ) = N spec , plate N phot , plate , (27)where N spec , plate is the number of BHBs on the plate thatmake the spectroscopic sample and N phot , plate is the numberof BHBs in the photometric sample within the same patch inthe sky. We evaluate these fractions by searching for BHBsin the spectroscopic and photometric samples in the regionscovered by each of the plates, using SkyServer’s CasJobs. InFig. 2 the colour scale shows p ( S | l, b ). The dependence of p ( S | l, b ) on b is evident. The selection probability in terms of distance and metal-licity depends on the assumed IMF, the survey selectionwith respect to apparent magnitudes and colours, and theisochrones used to relate apparent magnitudes and coloursto intrinsic properties. We assume the survey selection oncolour to be uniform over the ranges in Table 2. The signal-to-noise ratio however decays with apparent magnitude, andtherefore the survey selection is not uniform within the im-posed apparent magnitude range. Fig. 3 shows the g -bandapparent magnitude distribution for our sample of SEGUE-II stars and for stars in the SDSS photometric sample foundin Section 4.2.1 to lie within the plate dimensions and theselection criteria of Table 2. The left panel shows the distri-butions normalized to unity for g -band apparent magnitudesbetween 15 . .
0, just within the range imposed by theSEGUE-II targetting criteria for BHBs (15 . < g < . g -band apparent magnitude of19 .
0. Therefore we take the selection on apparent magnitudeto be uniform within the ranges specified in Table 2 and zerootherwise.We use the α -enhanced B.A.S.T.I. isochrones for amass-loss parameter η = 0 . s and [Fe/H] the prob-ability that a star randomly chosen at birth passes the mag-nitude and colour cuts p ( S | τ, s, [Fe / H]) = (cid:90) d m (cid:15) ( m ) p ( S | m, τ, s, [Fe / H]) . (28)For a given value of s , we integrate each isochrone (speci-fied by [Fe/H] and τ ) over the mass distribution, by adding a MNRAS , 1–18 (2015)
Payel Das (a) (b)
Figure 3.
Probability (a) and cumulative (b) g -band apparent magnitude distributions for the photometric and spectroscopic samples. star’s IMF value to p ( S | s, [Fe / H] , τ ) only if it passes the mag-nitude and colour cuts. We repeat this for a grid of distances,and create an interpolant of p ( S | s, [Fe / H] , τ ) that gives theselection function probability for any distance, metallicity,and age within the specified domain.The colour scale in Fig. 4 shows p ( S | s, [Fe / H] , τ ) for τ =9 ,
11, and 13 Gyr. The top plot shows that the BHBs cannotbe described by a single age of 9 Gyr, as stars of metallicitieshigher than − . p ( S | s, [Fe / H] , τ = 11) is approximately constant within well-defined boundaries. Therefore to avoid the time-consumingengagement with the isochrones, we assume a constant prob-ability within a box, so p ( S | s, [Fe / H]) = (cid:54) s (cid:54)
50 kpcand − . (cid:54) [Fe / H] (cid:54) − .
40 otherwise . (29)We adopt this age distribution for the models specified inSections 2.1, 2.2, and 2.3 ( M , M , and M ). For the modeldescribed in Section 2.4 ( M ), we assume the age to dependon actions within the same distance and metallicity limitsas specified in Equation (29). We neglect errors in sky coordinates and adopt independentGaussian error distributions for v (cid:107) , µ ∗ l , µ b , [Fe/H], and log s .Thus the multi-variate error distribution is C (7) ( u i , u (cid:48) , p i ) ≡ (cid:89) l =1 δ ( u il − u (cid:48) l ) (cid:89) l =3 C ( u il , u (cid:48) l , p il ) , (30) where true values of observables are indicated by primes and C ( u i , u (cid:48) , p i ) ≡ exp (cid:2) − ( u i − u (cid:48) ) / (2 p i ) (cid:3) √ πp i . (31)The convolution of the EDF with the error distribution of agiven star is P ( u i | M ) = (cid:90) d u (cid:48) C ( u i , u (cid:48) , p i ) f ( x , v , [Fe / H])) (cid:12)(cid:12)(cid:12)(cid:12) ∂ ( x , v , [Fe / H]) ∂ ( u (cid:48) ) (cid:12)(cid:12)(cid:12)(cid:12) . (32)The Jacobian determinant here is proportional to s cos b (McMillan & Binney 2012). The integral is calculated using afixed Monte Carlo sample of 2000 points per star to eliminatePoisson noise in likelihood evaluations (McMillan & Binney2013). The normalization of L is given by P ( S | M ) = (cid:90) d u (cid:48) P ( S | u (cid:48) ) f ( x , v , [Fe / H]) (cid:12)(cid:12)(cid:12)(cid:12) ∂ ( x , v , [Fe / H] ∂ ( u (cid:48) ) (cid:12)(cid:12)(cid:12)(cid:12) . (33)We approximate the integrals over sky coordinates by sumsof the remaining five-dimensional integrals evaluated at thecentre of each SEGUE-II plate (Das & Binney 2016). Since P ( S | u (cid:48) ) is multiplied by n ∗ in Equation (25), the integralsmust be calculated to a high degree of accuracy (McMil-lan & Binney 2012). We calculate them to an accuracy of0.1% using a Python wrapper for the cubature method . Theprocedure is parallelised over 16 cores using the distributedmemory tool in the multiprocessing package of Python. This wrapper can be downloaded from https://github.com/saullocastro/cubature . MNRAS000
40 otherwise . (29)We adopt this age distribution for the models specified inSections 2.1, 2.2, and 2.3 ( M , M , and M ). For the modeldescribed in Section 2.4 ( M ), we assume the age to dependon actions within the same distance and metallicity limitsas specified in Equation (29). We neglect errors in sky coordinates and adopt independentGaussian error distributions for v (cid:107) , µ ∗ l , µ b , [Fe/H], and log s .Thus the multi-variate error distribution is C (7) ( u i , u (cid:48) , p i ) ≡ (cid:89) l =1 δ ( u il − u (cid:48) l ) (cid:89) l =3 C ( u il , u (cid:48) l , p il ) , (30) where true values of observables are indicated by primes and C ( u i , u (cid:48) , p i ) ≡ exp (cid:2) − ( u i − u (cid:48) ) / (2 p i ) (cid:3) √ πp i . (31)The convolution of the EDF with the error distribution of agiven star is P ( u i | M ) = (cid:90) d u (cid:48) C ( u i , u (cid:48) , p i ) f ( x , v , [Fe / H])) (cid:12)(cid:12)(cid:12)(cid:12) ∂ ( x , v , [Fe / H]) ∂ ( u (cid:48) ) (cid:12)(cid:12)(cid:12)(cid:12) . (32)The Jacobian determinant here is proportional to s cos b (McMillan & Binney 2012). The integral is calculated using afixed Monte Carlo sample of 2000 points per star to eliminatePoisson noise in likelihood evaluations (McMillan & Binney2013). The normalization of L is given by P ( S | M ) = (cid:90) d u (cid:48) P ( S | u (cid:48) ) f ( x , v , [Fe / H]) (cid:12)(cid:12)(cid:12)(cid:12) ∂ ( x , v , [Fe / H] ∂ ( u (cid:48) ) (cid:12)(cid:12)(cid:12)(cid:12) . (33)We approximate the integrals over sky coordinates by sumsof the remaining five-dimensional integrals evaluated at thecentre of each SEGUE-II plate (Das & Binney 2016). Since P ( S | u (cid:48) ) is multiplied by n ∗ in Equation (25), the integralsmust be calculated to a high degree of accuracy (McMil-lan & Binney 2012). We calculate them to an accuracy of0.1% using a Python wrapper for the cubature method . Theprocedure is parallelised over 16 cores using the distributedmemory tool in the multiprocessing package of Python. This wrapper can be downloaded from https://github.com/saullocastro/cubature . MNRAS000 , 1–18 (2015) he age gradient in halo BHBs q s (kpc) −3.5−3.0−2.5−2.0−1.5 [ F e / H ] s (kpc) −3.5−3.0−2.5−2.0−1.5 [ F e / H ] s (kpc) −3.5−3.0−2.5−2.0−1.5 [ F e / H ] Figure 4.
Selection function as a function of distance and metal-licity, assuming a single age of 9 (top), 11 (middle), and 13 (bot-tom) Gyr, with locations of observed stars superimposed.
We explore the posterior distribution using two algorithms.The first method is the Nelder-Mead amoeba method (Nelder& Mead 1965), a downhill-simplex optimization routine forlocating an extremum of a multi-dimensional function whenthe derivatives of that function are not known. This algo-rithm has been shown to work well in many non-linear opti-mization problems, although it suffers from the common is-sue that one may only discover a local maximum, as opposedto a global one (if it exists). The second algorithm we useis the affine-invariant Monte Carlo Markov Chain approachimplemented in the emcee package for Python (Foreman-Mackey et al. 2013). The approach uses an interacting en-semble of ‘walkers’ and has been shown to be more effec-tive at dealing with narrow degeneracies than simpler ap-
Table 3.
Median and 68% confidence intervals for M , M , and M , and MAP estimates for M . A ‘*’ indicates that the param-eter is either set prior to the runs, or fixed by other parameters.Parameter M M M M q . +0 . − . - - - α in . +0 . − . - - - α out . +0 . − . - - - r b /kpc 29 . . − . - - - q - 0 . +0 . − . - - q ∞ - 0 . +0 . − . - - r /kpc - 7 . +1 . − . - - α - 4 . +0 . − . - - β in - - 2 .
05 2 . +0 . − . β out - - 4 .
72 4 . +0 . − . J - - 1635 1635 ∗ a r - - 1 .
31 0 . +0 . − . a φ - - 0 .
93 0 . +0 . − . a z - - 0 . ∗ . ∗ b r - - 1 .
11 1 . +0 . − . b φ - - 0 .
60 0 . +0 . − . b z - - 1 . ∗ . ∗ x - - 0 .
34 0 . +0 . − . [Fe / H] max − . +0 . − . − . +0 . − . − . +0 . − . − . +0 . − . [Fe / H] peak − . +0 . − . − . +0 . − . − . +0 . − . − . +0 . − . a τ - - - 12 . +0 . − . b τ - - - − . +0 . − . proaches, such as basic Metropolis-Hastings sampling. Wenow outline our priors and method of sampling for each ofour four models: • M : Spatial DF and metallicity DF with constant axisratio (Equation 2) and box-uniform distance-metallicity se-lection function (Equation 29). We fit all parameters using emcee with 40 walkers and 10,000 steps each. The priors areset as uniform within the following ranges:(i) 0 < q < < α in < α out < < r b <
40 kpc(iv) 0 < [Fe / H] max − [Fe / H] peak < • M : Spatial DF and metallicity DF with variable axisratio (Equation 8) and box-uniform distance-metallicity se-lection function (Equation 29). We fit the spatial DF pa-rameters only using emcee with 40 walkers and 10,000 stepseach. The metallicity DF parameters are independent andtherefore the same as found for M . The priors are uniformwithin the following ranges:(i) q > q ∞ > < r <
50 kpc(iv) 0 < α < • M : Action-based DF and metallicity DF (Equation12) with box-uniform distance-metallicity selection function(Equation 29). We first fit [ β in , β out , J ] using the amoeba algorithm, with the flattening/anisotropy and rotation pa- MNRAS , 1–18 (2015)
Payel Das
Figure 5. emcee results for fitting the DF M . The plots along the diagonal illustrate 1-D probability distributions for each of the modelparameters. The remaining plots show joint probability distributions between all pairs of parameters. The contours show the 1- σ and2- σ confidence levels. rameters fixed at those for an isotropic and round halo, andthe metallicity parameters fixed at those found for M . Wethen fit [ a r , a φ , b r , b φ ] using the amoeba algorithm, with thedensity parameters fixed from the first stage, and metallic-ity parameters fixed from M . The priors are uniform withinthe following ranges:(i) 0 < β in < β out > a r > . a φ > . a r + a φ < . b r > . b φ > . b r + b φ < . • M : Action-based DF, metallicity DF, and age DF(Equation 15) with age-varying distance-metallicity selec-tion function. We first fit [ x , a r , a φ , b r , b φ ] using the amoeba algorithm, with the density and metallicity parameters fixedat those found for M . Then we fit [ a τ , b τ ] using the amoeba algorithm, with the density, flattening/anisotropy, rotation,and metallicity parameters fixed at those found for the first stage. We then fit all parameters using the amoeba algorithm,and again using emcee with 22 walkers and 500 steps each.We require fewer steps than for M and M as we start veryclose to the MAP estimate. However, it cannot be guaran-teed that the chains have converged - the limit is simplya function of computational resources (22 walkers and 500steps take about three weeks to run across 16 cores). Thepriors are as for M with the following set to be uniform inthe given ranges:(i) a τ < . AGAMA ).This is a galaxy-modelling library in C++ consisting of sev-eral layers that together provide a complete package forconstructing galaxy models. The innermost layer providesa range of mathematical tools that include integration, in-terpolation, multi-dimensional samplers, and units. The cen- AGAMA can be downloaded from https://github.com/GalacticDynamics-Oxford/AGAMA
MNRAS000
MNRAS000 , 1–18 (2015) he age gradient in halo BHBs Figure 6. emcee results for fitting the DF M . The plots along the diagonal illustrate 1-D probability distributions for each of the modelparameters. The remaining plots show joint probability distributions between all pairs of parameters. The contours show the 1- σ and2- σ confidence levels. tral layers provide classes for gravitational potentials, actionfinders converting phase-space coordinates to actions, andDFs. The outermost layer comprises Python wrappers forseveral of the C++ classes, and an additional suite of Pythonroutines for fitting DFs, self-consistent modelling, generat-ing mock catalogues, isochrone interpolation, and selectionfunctions. A more detailed description of the library will begiven elsewhere (Vasiliev et al. in prep). Here we discuss the favoured parameters determined for ourfour models, and their uncertainties where derived. We as-sess the quality of the fit of the models to the observablesand investigate the moments associated with M . Table 3 gives the median of the recovered parameters for M , M , and M and the MAP estimates for the parameters afterthe multi-stage fit for M . From M , an axis ratio q ∼ . α in ∼ − . α out ∼ − . r b ∼ M , the axis ratio varies from q ∼ . q ∞ ∼ .
8, transitioning at a radius r ∼ α ∼ − .
7. Both densityprofiles would yield a divergent mass if extended towards thecentre, and are therefore not physical there. The metallicityDF has a cut-off at [Fe / H] max = − .
83 dex and peaks at[Fe / H] peak = − .
77 dex.For models M and M , the slope in actions steepenssmoothly from β in ∼ − . − . ∼ − to a slope of β out ∼ − . M results in a minorchange in the inner and outer halo power-law indices. Sinceany rotation within the halo is weak ( x ∼ J φ probably has little effecton the optimal values of the weights on the actions a r , b r and so on. However running an emcee chain, which ensuresthe parameter space is fully explored, finds different actionweights between M and M . In the latter, the isodensityellipsoids are flattened at low and high actions ( a φ (cid:28) a z and b φ (cid:28) b z ) and the velocity ellipsoids are elongated in MNRAS , 1–18 (2015) Payel Das − . − . − . [Fe/H] cutoff − . − . − . [ F e / H ] p e a k − . − . − . [Fe/H] peak Figure 7. emcee results for the metallicity DF of M , M , and M . The plots along the diagonal illustrate 1-D probability dis-tributions for each of the model parameters. The remaining plotshows the joint probability distribution between the two param-eters of the metallicity DF. The contours show the 1- σ and 2- σ confidence levels. the radial direction in the inner halo ( a r (cid:28) a z ). The agemodel predicts a mean age of ∼ . Table 3 gives the 68% confidence intervals for parametersof M , M , M (metallicity DF parameters only) and M .Figs. 5, 6, 7, and 8 present the 1-D and joint probabilitydistributions from the emcee runs. The metallicity DF isindependent of the spatial DF in M , M and of the phase-space DF in M . Its parameters are thus shown separately.The uncertainties on the parameters of M and M areof the order of ∼ M , whichare more uncertain ( ∼ M , and the outer axis ratio, flattening transition radius,and slope of M .For models M , M , and M , the uncertainties on themaximum and peak metallicities are smaller, at the levelof ∼ M vary greatly.Those on the parameters of the metallicity DF are greaterthan the uncertainties in M and M because of the flex-ibility introduced with the age model, which modifies thedistance-metallicity selection function. The uncertainties onthe action weights vary between ∼
15 to 35%. The uncer- tainty on the rotation parameter is high and translates to a68% confidence interval ∼ [ − ,
30] kms − for the rotationspeed.The mean age in the case of no dependence on ac-tions has an uncertainty ∼ .
33 Gyr, primarily towardshigher ages, because the distance-metallicity selection func-tion varies much less there. The dependence on age is quiteuncertain, again partly because of the insensitivity of thedistance-metallicity selection function to the oldest ages.The 1-D marginalized probability distributions in Fig. 8are unimodal, but often skewed. Asymmetrical distributionsinclude the mean age ( a τ ) in the case of no dependence onactions, and the dependence on actions ( b τ ), both of whichhave an extended tail towards higher values. The inner halopower-law index ( β in ) has an asymmetrical distribution withan extended tail towards lower values. Several correlationsexist between parameters. There is a weak, negative corre-lation between the power-law indices of the inner and outerhalo ( β in and β out ), i.e. a lower value of the index in theinner halo can be partly compensated by a larger value ofthe index in the outer halo. There is also a correlation be-tween the weight on the radial action in the outer halo ( b r )and the weight on the angular action there ( b φ ), and to alesser extent in the inner halo. This highlights a connectionbetween the flattening and anisotropy, i.e. more flattenedsystems tend to have a higher degree of radial anisotropy.There is also the correlation between the parameters of themetallicity DF. Other correlations may also exist but areobscured by the lack of sufficient resolution in the emcee runs. We assess the quality of the model fits by generating mockcatalogues at the measured sky positions, using an adap-tive sampling-rejection method from
AGAMA . The productof the SF and EDF is sampled in log( s/ kpc), log [Fe / H],and for M and M in tanh( v r / km s − ), tanh( µ ∗ l / mas), andtanh( µ b / mas). In these coordinates the value of the EDFvaries less strongly. Each proposed sample has noise addedto it according to the errors measured on the observablesat those sky positions. The parameters used to generate themock samples in each case correspond to the MAP esti-mates.The coloured points (cyan - M , orange - M , red - M , and green - M ) joined by lines in Figs. 9 and 10 showhistograms of the mock catalogues. The region covered byanalogous histograms of 2000 resamplings of the error dis-tributions of the measured observables are shown by thegrey regions. In plots for observables with small errors, suchas l , b , s , and v (cid:107) , the grey regions form fairly well definedcurves. In plots for observables with large errors, namely( µ ∗ l , µ b , [Fe / H]), the grey regions fill out a region of signifi-cant width. The coloured curves generally overlap with thisregion, indicating that the EDF and SF are together doinga good job at describing the location of the observables.Fig. 11 compares histograms for the joint distributionof pairs of phase-space observables ( l, b ), ( s, b ), ([Fe / H] , b ),etc. in the case of parameters corresponding to the MAPestimate obtained for M . The colour-filled contours showdistributions of mock observables, while the black contoursshow the distributions of measured observables. In general MNRAS000
AGAMA . The productof the SF and EDF is sampled in log( s/ kpc), log [Fe / H],and for M and M in tanh( v r / km s − ), tanh( µ ∗ l / mas), andtanh( µ b / mas). In these coordinates the value of the EDFvaries less strongly. Each proposed sample has noise addedto it according to the errors measured on the observablesat those sky positions. The parameters used to generate themock samples in each case correspond to the MAP esti-mates.The coloured points (cyan - M , orange - M , red - M , and green - M ) joined by lines in Figs. 9 and 10 showhistograms of the mock catalogues. The region covered byanalogous histograms of 2000 resamplings of the error dis-tributions of the measured observables are shown by thegrey regions. In plots for observables with small errors, suchas l , b , s , and v (cid:107) , the grey regions form fairly well definedcurves. In plots for observables with large errors, namely( µ ∗ l , µ b , [Fe / H]), the grey regions fill out a region of signifi-cant width. The coloured curves generally overlap with thisregion, indicating that the EDF and SF are together doinga good job at describing the location of the observables.Fig. 11 compares histograms for the joint distributionof pairs of phase-space observables ( l, b ), ( s, b ), ([Fe / H] , b ),etc. in the case of parameters corresponding to the MAPestimate obtained for M . The colour-filled contours showdistributions of mock observables, while the black contoursshow the distributions of measured observables. In general MNRAS000 , 1–18 (2015) he age gradient in halo BHBs . . . β in − . b τ a τ − . − . − . − . [ F e / H ] p e a k − . − . − . [ F e / H ] m a x − . . x . . b φ . . . b r . . . a φ . . a r . . β o u t . . β out . . a r . . . a φ . . . b r . .
75 1 b φ − . . x − . − . − . [Fe/H] max − . − . − . − . [Fe/H] peak
11 12 13 a τ − − . b τ Figure 8. emcee results for fitting the EDF M . The plots along the diagonal are 1-D probability distributions for each of the modelparameters. The remaining plots show joint probability distributions between all pairs of parameters. The contours show the 1- σ and2- σ confidence levels. there is good agreement between the mock and observationaldistributions. M We now describe the model M with the parameters corre-sponding to the MAP estimate. Fig. 12 shows the shape of the density distribution. Thecolour scale in the top panel shows ρ ( R, z ). Flattening ofthe contours is evident. By fitting ellipses to the isoden-sity curves, we obtain the radial density profile shown bythe green curve in the second panel. The steepening of thedensity profile with increasing radius is shown by the greencurve in the third panel, which gives the logarithmic radialdensity gradient. The slope steepens smoothly from − . ∼ ∼ − q (cid:39) . MNRAS , 1–18 (2015) Payel Das
10 20 30 40 50 60 s (kpc) N u m b e r −400 −200 0 200 400 v || (km/s) N u m b e r −15 −10 −5 0 5 10 15 µ ∗l (mas/year) N u m b e r −15 −10 −5 0 5 10 15 µ b (mas/year) N u m b e r Figure 9.
In colour: one-dimensional distributions of mock phase-space observables (cyan - M , orange - M red - M , green - M ).The error bars show Poisson errors. Grey regions: 2000 Monte Carlo resamplings of the data from the error distributions. Progressivelylighter shades of grey indicate 1 σ to 3 σ and 100% regions for the resampled observables. −3.0 −2.8 −2.6 −2.4 −2.2 −2.0 −1.8 −1.6 −1.4 [Fe/H] (dex) N u m b e r Figure 10.
In colour and grey: as Fig. 9 but for metallicities. .
8) throughout. Radial profiles of the logarithmic density,logarithmic density gradient, and axis ratio are shown alsofor M (cyan) and M (orange), each plotted against ellipti-cal radius. The density profile of M is steeper in the innerand outermost parts than that of M , while that of M issteeper throughout. The axis ratio of M is very similar tothat of M . The axis ratio of M is significantly lower thanthat of M and M in the inner halo and moderately higherin the outer halo. Fig. 13 shows the velocity dispersions of M . σ r gen-erally dominates at high z throughout, implying radialanisotropy there. At low z , σ φ dominates, implying tangen-tial anisotropy. The contours of constant σ r and σ θ are elon-gated in the z direction, while σ φ contours are elongated inthe R direction. Fig. 14 shows the spherical anisotropy pa-rameter β s = 1 − σ θ + σ φ σ r . (34)The degree of radial anisotropy increases from a tangentialbias in the equatorial plane to ∼ . z -axis. The lower panel of Fig. 14 shows radialprofiles of the spherical anisotropy parameter against spher-ical radius along a range of polar angles, θ , where θ = 0 isalong the z axis and θ = 90 o is in the equatorial plane. Inthe equatorial plane, the orbits vary from isotropic to mildlytangential. Nearer the z axis, the profiles become more ra-dially anisotropic. The top panel of Fig. 15 shows the distribution of meanages. Contours of constant mean age are flattened, approx-imately as the contours of constant density. The bottompanel shows the age map inferred from the actions of the
MNRAS000
MNRAS000 , 1–18 (2015) he age gradient in halo BHBs −60−3003060 b ( o ) s ( k p c ) −450−300−1500150300 v || ( k m / s ) −10−50510 µ ∗ l ( m a s / y e a r ) −16−80816 µ b ( m a s / y e a r )
60 120 180 240 l ( o ) −2.50−2.25−2.00−1.75−1.50 [ F e / H ] ( d e x ) −40 0 40 80 b ( o )
10 20 30 40 50 s (kpc) −400 −200 0 200 v || (km/s) −12 −6 0 6 µ ∗l (mas/year) −10 0 10 20 µ b (mas/year) Figure 11.
Colour-filled contours: two-dimensional distributions of measured observables. Black contours: distributions of the mockobservables for M . stars in the SEGUE-II sample, convolved with their errordistributions. In general stars at lower R and z have higherages, but a high velocity implies large actions and thereforea relatively young age. Therefore the gradient in age withactions manifests as a small negative age gradient with ra-dius ∼ − .
03 Gyr kpc − . The plot also suggests that we arebiased towards observing younger stars. Here we focus our discussion on the results of fitting the fullphase-space EDFs, how they compare to the literature, andhighlight uncertainties that may impact these results.
Traditionally metallicity and age gradients have been ex-amined as a function of radius. However, stars are bettercharacterised by their actions, which do not vary along or-bits. A clear separation in action space becomes ‘smeared’in real space, and a glance at the picture in action space canbe enlightening. We find that the EDF declines more rapidlywith actions in the outer halo (slope ∼ −
5) than in the in-ner halo (slope ∼ − z increases. The part of the EDF odd in J φ MNRAS , 1–18 (2015) Payel Das
10 20 30 40 50 60 70
R (kpc) z ( k p c ) ln(ρ) −18.5−16.0−13.5−11.0−8.5−6.010 Elliptical radius (kpc) −25−20−15−10−5 l n ( ρ ) Elliptical radius (kpc) −5.0−4.5−4.0−3.5−3.0−2.5−2.0 d l n ρ / d l n r Elliptical radius (kpc) q Figure 12.
Real-space density: (1) contour map of M (toppanel), (2) radial density profiles of models M (cyan), M (or-ange), and M (green) (second panel), (3) logarithmic radial den-sity gradient (third panel), and (4) axis ratio for models M , M ,and M with the same colour coding. is negligible. We find the ages of the stars to be well pre-dicted by a log-linear dependence on the total action, withhigher ages at smaller actions. The gradient is thus negative( − .
69 Gyr dex − ). A single-age model is, however, also ableto reproduce the observations. The slopes in action space translate to a density profile inreal space that steepens with radius from a slope of ∼ − ∼ ∼ − ∼ − .
03 Gyr kpc − ,subject to a significant degree of uncertainty.The weights on actions determine the shape of the den-sity and velocity ellipsoids. The halo’s axis ratio is roughlyconstant at ∼ .
7, very similar to what is implied bythe broken power-law model. The orbital structure variesfrom mildly tangential to moderately radially anisotropicthroughout, becoming more radial as you move from theequatorial plane towards the z axis.The negligible part of the EDF that is odd in J φ gener-ates a very small level of rotation with a large uncertainty(our 68% confidence interval extends between −
10 to 30 kms − ).Finally, we did not probe the distribution of metallici-ties in action space directly. It is however linked to the EDFthrough the selection function in distance and metallicity.The metallicity DF is well described by a single lognormaldistribution with a peak at ∼ − . ∼ − . There is compelling evidence for differences between the in-ner and outer halo in action space, primarily due to thedifference in the inner and outer slopes in actions. This man-ifests in real space as a steepening of the density profile withradius, a non-negligible negative age gradient with radius,and a variation in the anisotropy with radius.There have been several claims of two populations inthe halo (e.g. Carollo et al. 2007; Beers et al. 2012; Deasonet al. 2013; Hattori et al. 2013). Although our sample istoo small to rule out such a dichotomy, our models, whichpredict smooth transitions between the inner and outer halo,are sufficient to reproduce the current data.
Several authors have determined density profiles for haloBHBs, finding a power-law index of ∼ − . − . − . − . . . . . MNRAS000
Several authors have determined density profiles for haloBHBs, finding a power-law index of ∼ − . − . − . − . . . . . MNRAS000 , 1–18 (2015) he age gradient in halo BHBs
10 20 30 40 50 60 70
R (kpc) z ( k p c ) σ r (km/s) R (kpc) z ( k p c ) σ θ (km/s) R (kpc) z ( k p c ) σ φ (km/s) Figure 13.
Velocity dispersions predicted by the best-fitting EDF. Going from the left to right: spherical radial velocity dispersion,spherical angular velocity dispersion, and spherical azimuthal velocity dispersion.
10 20 30 40 50 60 70
R (kpc) z ( k p c ) β s −0.35−0.30−0.25−0.20−0.15−0.10−0.050.000.050.100.150.200.250.3010 20 30 40 50 60 70 r (kpc) −0.4−0.20.00.20.40.6 β s θ=0 o θ=30 o θ=50 o θ=85 o Figure 14.
Spherical anisotropy parameter predicted by the best-fitting EDF plotted as a map in R and z (top panel) and againstspherical radius for a range of colatitudes θ (bottom panel). Deason et al. (2011a) and Hattori et al. (2013) found theBHB stars in the Milky Way halo to exhibit a dichotomy be-tween a prograde-rotating, comparatively metal-rich compo-nent ([Fe / H] >
2) and a retrograde-rotating, comparativelymetal-poor ([Fe / H] <
2) component. Deason et al. (2011a)attribute the prograde metal-rich population to the accre-tion of a massive satellite ( ∼ M (cid:12) ) and the metal-poorpopulation to the primordial stellar halo. The net retrograderotation might then reflect an underestimate in the adoptedLSR circular velocity. Fermani & Sch¨onrich (2013b) however remeasured the rotation of the Milky Way stellar halo on twosamples of BHB halo stars from SDSS with four differentmethods, and found a weakly prograde or non-rotating haloin all cases. They attributed the rotation gradient acrossmetallicity measured by Deason et al. (2011a) on a similarsample of BHB stars to the inclusion of regions in the appar-ent magnitude-surface gravity plane known to be contami-nated by substructures. Sirko et al. (2004) did not find anyrotation in their sample of BHBs and Smith et al. (2009a)do not find any rotation in their sample of subdwarfs.Several authors (Deason et al. 2012; Kafle et al. 2012;Williams & Evans 2015; Cunningham et al. 2016) found a ra-dially biased velocity ellipsoid overall but some find a regionaround ∼
20 kpc with tangential anisotropy. From a similarsample of BHBs, Sirko et al. (2004) found isotropy and Hat-tori et al. (2013) found the metal-rich component to exhibitmild radial anisotropy, and the metal-poor component toexhibit tangential anisotropy. Analysis of other halo tracerpopulations have also arrived at a mixture of conclusions.Dehnen et al. (2006) found radial anisotropy in a mixture ofglobular clusters, horizontal-branch and red-giant stars, anddwarf spheroidal satellites. Smith et al. (2009b) found thevelocity ellipsoid of SDSS subdwarfs to be radially biased.Carollo et al. (2010) found inner-halo metal-richer stars onradially anisotropic orbits, and outer-halo stars to be on lesseccentric orbits.The diversity in conclusions regarding halo anisotropymay be a result of several things. Spectroscopic surveys canhave a strong selection bias on metallicity (not so much inBHBs, but in K giants, see Das & Binney (2016)), and ifthere is a correlation between metallicity and dynamics, thena lack of treatment of such a bias can lead to erroneous con-clusions about anisotropy. It is also true that the proper mo-tions currently available are highly uncertain. Furthermore,it should be emphasised that our models are designed to fitthe smooth, phase-mixed halo, and this may be why we donot reproduce the ‘dip’ in anisotropy found by several au-thors at 20 kpc. On a related note, we have attempted to re-move substructures where possible - the exact samples usedby various authors differ slightly in their observed velocitydistributions and thus derive different anisotropy profiles.Our model qualitatively agrees with simulations at high z (i.e. radial bias increases outward), though with a lesserdegree of radial bias throughout. In such simulated haloes,the primary mechanism for growth since z ∼ MNRAS , 1–18 (2015) Payel Das
10 20 30 40 50 60 70
R (kpc) z ( k p c ) Age (Gyr)
R (kpc) z ( k p c ) Figure 15.
Age map predicted from the total action moment(top panel) and for the stars in the SEGUE-II sample (bottompanel). be accretion onto the halo through minor mergers (Bullock& Johnston 2005; Abadi et al. 2006), and accreting objectshave rather radial orbits.
Carollo et al. (2007) and Beers et al. (2012) found evidencefor two metallicity components in a mixed sample of stellartypes. An et al. (2015) analysed a sample of main-sequencehalo stars and found two components peaking at [Fe/H] ∼− . ∼ − .
3. Xue et al. (2015) and Das & Binney (2016)reached similar conclusions about SDSS K giants.We find one lognormal component is sufficient to de-scribe the metallicity DF of the BHBs. The discrepancy mayarise from a difference between the types of systems thoughtto contribute to halo stars. There is an age-metallicity bi-modality in the Milky Way globular cluster system (e.g. Lea-man et al. 2013). Fiorentino et al. (2015) analysed the peri-ods and luminosity amplitudes of field RR Lyrae stars andfound that dwarf spheroidals lacked high-amplitude short-period variable stars; whereas these are found in globularclusters and massive dwarf irregulars such as the Sagittar-ius stream.Preston et al. (1991) detected a colour gradient in BHBsout to ∼
12 kpc, which Santucci et al. (2015) consolidatedwith a larger sample, extending out to ∼
40 kpc. Theyclaimed that the gradient is independent of metallicity andtherefore indicate a gradient in age. More massive systemspenetrate deeper into the gravitational potential. We alsoexpect these systems to have the oldest components becausethey would have grown over a longer timescale. From a simi-lar argument however, we would also expect chemical evolu- tion to have occurred more rapidly in these systems, there-fore producing a metallicity gradient. It is unclear why thedifference in the make-up of the inner and outer halos wouldmanifest solely as an age gradient rather than a metallicitygradient. The metallicity gradient could just be small.
The errors quoted in this work represent statistical uncer-tainties only, rather than systematic errors, which are moredifficult to characterise. We discuss possible sources of sys-tematic errors below.
We have assumed a fully-integrable potential in which thenumber of isolating integrals of motion is equal to the num-ber of degrees of freedom. In such cases a transformationfrom phase-space coordinates to angle-action coordinatescan be done globally. Real potentials however permit fami-lies of resonantly trapped orbits. These resonant orbits pos-sess actions, but require a different transformation. In non-integrable potentials, some fraction of the orbits will bechaotic. Chaotic orbits are not bound to the surface of atorus and instead fill the spaces between tori. Chaotic or-bits have no orbital actions. Binney (2016) concluded thatthe net impact of resonant trapping on the dynamics of halostars is likely to be small.
An EDF can perfectly model the data only in the true poten-tial. Therefore any imperfections in our choice of potentialwill both bias our EDF away from the true EDF and giverise to discrepancies between our best model and the data.Our success in reproducing the data suggests that our cho-sen potential is not seriously in error. The supposition of aparticular functional form for the EDF can bias the resultsby restricting the set of possible solutions, despite allowinga range of density and anisotropy profiles. Our ansatz re-garding the dependence of stellar ages on actions representsjust one, physically-motivated, possibility that is simple tocalculate. An age gradient is not forced by the EDF how-ever; if none were needed by the data, age would have beenfound to be independent of actions.
Our evaluation of the distance-metallicity selection func-tion depended on relations from isochrones between the age,mass, and metallicity of a star and its luminosity in variouswavebands. Systematic errors arising from faulty isochronesare difficult to assess. The ability of our EDF to produce asimilar density profile for the BHBs to that in the literaturesuggests that our selection function is not significantly inerror.
We have masked the Sagittarius stream in this analysis, butwe do not know how other, unmasked substructures may im-
MNRAS000
MNRAS000 , 1–18 (2015) he age gradient in halo BHBs pact our assessment of the halo’s structure. Our ability tosufficiently reproduce the phase-space observables after ex-cluding the stream implies that either we are predominantlyprobing the smooth stellar halo with the data or the currentdata are too sparse to resolve the halo’s substructure. We probed the chemodynamical structure of Milky Way haloBHBs by combining spatial and action-based EDFs that de-scribe the locations of stars in phase space, metallicity, andage. The analysis allows a more natural description of theages of BHBs in action space in which their separation isclearer than in real and velocity space. The specificationof an EDF enables the incorporation of a realistic selectionfunction that takes into account restrictions on sky posi-tions, apparent magnitudes, and colours. In general, ourmodels reproduce the observations well. This may be anargument that there is enough phase-mixed debris for ac-tion space to be smoothly populated (at least for relativelytightly bound orbits). i.e., there may be a part of the innerhalo that will always be well represented by smooth models.Alternatively it may be because the data for BHBs are notyet rich enough to resolve most halo substructures.The EDF of the BHBs is steeper at larger actions thanat smaller actions. Older stars are found at smaller actionsand younger stars at larger actions. The spatial distributionof the stars is similarly well reproduced by a broken powerlaw with a constant axis ratio, a single power law with a vari-able axis ratio, and a gradually steepening power law with avariable axis ratio. Fitting positions and velocities simulta-neously yields a density profile that steepens smoothly from ∼ − ∼ ∼ − ∼ . ∼ − . ∼ − . −
10 kms − to 30 km s − at most but the median result favours norotation. The stellar velocity ellipsoid varies from tangentialbias in the equatorial plane to radially elongated at high z .Allowing a dependence of stellar ages on actions leads toan age gradient ∼ − .
03 kpc − , with moderate uncertainty.However, an EDF assuming approximately a single age of11 Gyr is also able to fit the observables well.There are several possible directions for further work.The EDF could be applied to detect substructures in a richersample of halo stars in the phase-space-metallicity domain.The EDF could be changed to make the transition betweenthe inner and outer asymptotic slopes of the density profilesharper. The EDF could be further elaborated to include adependence on [ α /Fe]. ACKNOWLEDGEMENTS
The research leading to these results has received fund-ing from the European Research Council under the Euro-pean Union’s Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement no. 321067. PD thanks GitHubfor providing free private repositories for educational use. AW acknowledges the support of STFC. PD is also gratefulfor fruitful discussions with members of the Oxford GalacticDynamics group.
REFERENCES
Abadi M. G., Navarro J. F., Steinmetz M., 2006, MNRAS, 365,747An D., Beers T. C., Santucci R. M., Carollo D., Placco V. M.,Lee Y. S., Rossi S., 2015, ApJ, 813, L28Beers T. C., et al., 2012, ApJ, 746, 34Bell E. F., et al., 2008, ApJ, 680, 295Bell E. F., Xue X. X., Rix H.-W., Ruhland C., Hogg D. W., 2010,AJ, 140, 1850Binney J., 2012, MNRAS, 426, 1324Binney J., 2016, MNRAS,Binney J., Gerhard O., Spergel D., 1997, MNRAS, 288, 365Bullock J. S., Johnston K. V., 2005, ApJ, 635, 931Carollo D., et al., 2007, Nature, 450, 1020Carollo D., et al., 2010, The Astrophysical Journal, 712, 692Cunningham E. C., et al., 2016, ApJ, 820, 18Das P., Binney J., 2016, MNRASDe Propris R., Harrison C. D., Mares P. J., 2010, ApJ, 719, 1582Deason A. J., Belokurov V., Evans N. W., 2011a, MNRAS, 411,1480Deason A. J., Belokurov V., Evans N. W., 2011b, MNRAS, 416,2903Deason A. J., Belokurov V., Evans N. W., An J., 2012, MNRAS,424, L44Deason A. J., Van der Marel R. P., Guhathakurta P., Sohn S. T.,Brown T. M., 2013, ApJ, 766, 24Dehnen W., Binney J., 1998, MNRAS, 294, 429Dehnen W., McLaughlin D. E., Sachania J., 2006, MNRAS, 369,1688Fermani F., Sch¨onrich R., 2013a, MNRAS, 430, 1294Fermani F., Sch¨onrich R., 2013b, MNRAS, 432, 2402Fiorentino G., et al., 2015, ApJ, 798, L12Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013,PASP, 125, 306Guo J.-C., Liu C., Liu J.-F., 2016, Research in Astronomy andAstrophysics, 16, 008Hattori K., Yoshii Y., Beers T. C., Carollo D., Lee Y. S., 2013,The Astrophysical Journal Letters, 763, L17Hawkins K., Jofr´e P., Gilmore G., Masseron T., 2014, MNRAS,445, 2575Jeans J. H., 1916, MNRAS, 76, 567Jofr´e P., Weiss A., 2011, A&A, 533, A59Kafle P. R., Sharma S., Lewis G. F., Bland-Hawthorn J., 2012,ApJ, 761, 98Kalirai J. S., 2012, Nature, 486, 90Kinman T. D., Suntzeff N. B., Kraft R. P., 1994, AJ, 108, 1722Kroupa P., Tout C. A., Gilmore G., 1993, MNRAS, 262, 545Leaman R., VandenBerg D. A., Mendel J. T., 2013, MNRAS, 436,122Marquez A., Schuster W. J., 1994, A&AS, 108McMillan P. J., Binney J., 2012, MNRAS, 419, 2251McMillan P. J., Binney J. J., 2013, MNRAS, 433, 1411Nelder J. A., Mead R., 1965, The Computer Journal, 7, 308Pietrinferni A., Cassisi S., Salaris M., Castelli F., 2006, ApJ, 642,797Piffl T., et al., 2014, MNRAS, 445, 3133Piffl T., Penoyre Z., Binney J., 2015, MNRAS, 451, 639Posti L., Binney J., Nipoti C., Ciotti L., 2015, MNRAS, 447, 3060Preston G. W., Shectman S. A., Beers T. C., 1991, ApJ, 375, 121Sanders J. L., Binney J., 2015, MNRAS, 449, 3479Santucci R. M., et al., 2015, ApJ, 813, L16Sch¨onrich R., 2012, MNRAS, 427, 274MNRAS , 1–18 (2015) Payel Das
Sesar B., et al., 2013, AJ, 146, 21Sirko E., et al., 2004, AJ, 127, 914Sluis A. P. N., Arnold R. A., 1998, MNRAS, 297, 732Smith M. C., et al., 2009a, MNRAS, 399, 1223Smith M. C., Wyn Evans N., An J. H., 2009b, ApJ, 698, 1110Watkins L. L., et al., 2009, MNRAS, 398, 1757Williams A. A., Evans N. W., 2015, MNRAS, 454, 698Xue X.-X., et al., 2011, ApJ, 738, 79Xue X.-X., Rix H.-W., Ma Z., Morrison H., Bovy J., Sesar B.,Janesh W., 2015, preprint ( arXiv:1506.06144 ) MNRAS000