A centrally heated dark halo for our Galaxy
We construct a new family of models of our Galaxy in which dark matter and disc stars are both represented by distribution functions that are analytic functions of the action integrals of motion. The potential that is self-consistently generated by the dark matter, stars and gas is determined, and parameters in the distribution functions are adjusted until the model is compatible with observational constraints on the circular-speed curve, the vertical density profile of the stellar disc near the Sun, the kinematics of nearly 200000 giant stars within of the Sun, and estimates of the optical depth to microlensing of bulge stars. We find that the data require a dark halo in which the phase-space density is approximately constant for actions . In real space these haloes have core radii .
keywords:dark matter - galaxies: haloes - solar neighbourhood - Galaxy:disc - Galaxy:fundamental parameters - Galaxy: halo
It is now generally accepted that most of the mass of galaxies like ours is contained in a dark halo that is made up of particles that have yet to be discovered. Since these dark-matter (DM) particles have so far only come to our notice through the gravitational field that they generate, the only way to discover how they are distributed is to model that gravitational field using tracer particles. The tracer particles range from photons (via gravitational lensing), through electrons (via thermal X-ray emission) to stars and neutral atoms. Stars and atoms are the tracer particles of choice in our own Galaxy.
We understand the dynamics of our Galaxy at and inwards of the solar radius much better than we do at large radii, in part because determining distances to and tangential velocities of tracers at is hard, and in part because the density of stars and neutral gas is small at , so statistical uncertainties become large.
One constrains the distribution of DM by building dynamical models of our Galaxy that are consistent with relevant data. These data include (i) the circular-speed curve that one extracts from radio-frequency emission lines of interstellar gas, the proper motion of Sgr A* and the kinematics of stellar masers; (ii) the kinematics of stars that lie close enough to the Sun to have useful proper motions; (iii) star counts, which strongly constrain the vertical structure of the stellar disc at ; (iv) measurements of the optical depth to microlensing of bulge stars, since these measurements constrain the amount of stellar as opposed to interstellar or dark mass at .
As a first approximation one usually assumes that the Galaxy is axisymmetric. Then equilibrium dynamical models of the stars and DM are most conveniently formulated in terms of a distribution function (DF) that depends on the three action integrals , and . quantifies a star’s radial excursions, quantifies its oscillations perpendicular to the Galactic plane, and is the component of angular momentum L abut the symmetry axis.
Piffl et al. (2014) (hereafter P14) used models of this type to obtain very tight limits on the mass of DM at . The central idea of this work is that the kinematics of stars in the RAdial Velocity Experiment (RAVE Steinmetz et al., 2006) essentially fix the dependence of on and , and this dependence essentially fixes the run of vertical gravitational force required to produce the stellar density profile determined from star counts. The required is produced by a combination of matter in the disc and the dark halo, and if the shape of each contribution is known, both normalisations can be recovered from . Thus the data strongly constrain the local surface density of stars and the volume density of DM. If one now assumes, as P14 and all earlier authors did, that one knows the functional form of the DM density profile, then the entire structure of the dark halo can be fixed from local data. The scale-length of the disc then follows from the measured circular-speed curve .
The major uncertainty in this modelling is the flattening of the dark halo. P14 found a one-parameter family of successful models with dark haloes that had essential the same mass inside the isodensity surface through the Sun but differed in the axis ratio of those surfaces. The flatter a dark halo was, the more it contributed to both and , and therefore the smaller was the required mass of the disc.
P14 assumed that in the spherical case the density profile of the dark halo is the NFW profile (Navarro, Frenk & White, 1997)
where and are constants. This profile fits dark haloes in cosmological simulations that do not contain baryons. Piffl, Penoyre & Binney (2015) (hereafter P15) observed that we would not expect this profile to fit the Galaxy’s dark halo, since the latter is subject to the disc’s non-trivial gravitational field. Since Blumenthal et al. (1984) it has been argued that in galaxies like ours, in which baryons have accumulated steadily over a long period, the dark halo would have responded to the strengthening gravitational field adiabatically – for recent work see Katz et al. (2014). That is, the DF of the halo particles would be invariant as the disc and bulge accumulated. Consequently, P15 replaced the assumption of an NFW density profile with the assumption that had a certain form, which they chose such that, in the absence of a disc or bulge, the dark halo would coincide with the NFW halo fitted to the data by P14. When the DF of the disc was also set to that found by P14, the resulting model violated the constraints on at small . Thus, when near the Sun the balance between disc and halo mass is that required by the data and one takes into account the tendency of the dark halo to deform elastically as the disc and bulge grow, excessive mass accumulates near the Galactic centre.
The process of searching model space for models that are consistent with the data becomes much more costly computationally when the dark halo is specified by rather than . Hence P15 only computed a single model, that was based on the parameters found by P14. Binney & Piffl (2015) (hereafter BP15) conducted a systematic search for a model that (i) has a dark halo specified by a DF that would in isolation generate an NFW profile, and (ii) is consistent with and the local kinematic and star-count data. They found such a model. It avoided placing too much mass at small radii by assigning the disc a large scale radius . As a consequence, DM dominated the gravitational force down to small radii, and in this case the optical depth to lensing bulge stars falls below observational requirements. The clear conclusion from this exercise is that the DM did not respond adiabatically to the accumulation of baryons, so now would not in isolation generate an NFW profile. This conclusion is actually not unexpected because in the NFW model the phase-space density tends to infinity at the origin of action space. While it is in principle possible that the infinite phase-space density of the CDM initial conditions survives structure formation at the centres of dark haloes, it is improbable that scattering of DM particles has bot flattened the DM density near the centres of haloes.
The purpose of this paper is to obtain a form of that is consistent with all the observational data and our understanding of cosmology. In Section 2 we summarise the observational constraints imposed on models. In Section 3.1 we modify a DF that generates an NFW halo so the phase-space density of DM tends to a constant for the most bound particles. In Sections 3.2 to 3.5 we specify the other components of our model Galaxy, which include DFs for the thin and thick discs, and explain how the parameters of the DFs are adjusted to obtain self-consistent models that are compatible with the observational constraints. Section 4 describes the some successful models. Section 5 examines these models from the perspective of the microlensing data. Section 6 we investigate the symptoms of adopting an excessively large core for the dark halo. In Section 7 we review studies of the formation of cored dark halos. In Section 8 we sum up and look to the future.
2 Observational inputs
We used the same observational inputs as BP15. Here we summarise the inputs; more detail can be found in BP15. Our solar parameters are given in Table 1.
|0.014||Binney, Gerhard & Spergel (1997)|
|(11.1,12.24,7.25)||Schönrich, Binney & Dehnen (2010)|
2.1 Gas terminal velocities
The distribution of Hi and CO emission in the longitude-velocity plane yield a characteristic maximum (“terminal”) velocity for each line of sight (e.g. Binney & Merrifield, 1998, §9.1.1) which are related to the circular speed . We use the terminal velocities from Malhotra (1995). Following Dehnen & Binney (1998) and McMillan (2011) we neglect data at in order not to be influenced by the Galactic bar, and we assume that the ISM has a Gaussian velocity distribution of dispersion 7 km s.
2.2 Maser observations
We use 103 maser observations from Reid et al. (2014) that provide 6D phase-space information. The maser sources, which are associated with young stars, are assumed to be on nearly circular orbits: their velocities are assumed to be Gaussianly distributed about the circular velocity with dispersion (van der Kruit & Shostak, 1984; McMillan & Binney, 2010). For the likelihood computation we neglected all sources at to prevent the Galactic bar giving rise to a bias.
2.3 Proper motion of SgrA*
We adopt from Reid & Brunthaler (2004) the proper motion
of the radio source SgrA* associated with the super-massive black hole that sits in the Galactic Centre, as an estimate for the solar motion with respect to the GC.
2.4 Vertical density profile from SDSS
We assume that the population from which the RAVE sample is drawn is identical to that studied by Jurić et al. (2008) (hereafter J08). We use the data points shown in the middle panel of Figure 15 in J08, which shows results from M dwarf stars in the colour range , this sample should carry only weak biases in metallicity and age. We omitted the correction of the data for the effects of Malmquist bias and binarity as they had a negligible effect on the results of P14 and we decomposed the density profile into contributions from the disc and stellar halo as in BP15.
2.5 Kinematics from RAVE
We use the stellar parameters and distance estimates in the fourth RAVE data release (Kordopatis et al., 2013). We sort the stars into eight spatial bins, four inside the solar cylinder and four outside. We compute the velocity distributions predicted by our DF at the mean positions (barycentre) of the stars in each bin. We have a histogram for each velocity component, so we accumulate from these 24 histograms.
We modify the model distributions to take into account the effect of errors in the velocity and parallax estimates in the data. This procedure is fully described in P14. Our model selection involves optimising the fit between the data and the velocity histograms after the latter have been modified to allow for the impact of errors in the measurements of velocity and distance.
3 Modelling procedure
Most of the components used in our model are the same as those used by BP15, the main difference being the dark halo. We describe our dark halo DF first and then summarise the other components.
3.1 Distribution Functions for heated DM
The central density cusp of the NFW model implies divergence of the phase-space density of DM particles as their action integrals go to zero because in the cusp the velocity dispersion must tend to zero. Quantitatively, Posti et al. (2015) showed that the DF
with a homogeneous function of degree unity, self-consistently generates a system that closely resembles the NFW profile, with the scale action encoding the scale radius around which the slope of the radial density profile shifts from at small radii to far out. From equation (3) it follows that as .
The message from BP15 is that the DF of DM cannot increase as strongly as as an NFW profile predicts, either because the infinite phase-space density of the CDM initial conditions does not survive structure formation even at the centres of dark haloes, or because the baryons did not accumulate entirely adiabatically, and the fluctuating gravitational potential associated with them has upscattered what were the most tightly bound DM particles. This conclusion parallels what has been learnt from studies of dwarf galaxies: their rotation curves rise less steeply near their centres that would be expected if their dark haloes had cental density cusps like that of the NFW profile (Amorisco, Agnello & Evans, 2013, and references therein), again implying that at some stage even the most bound DM particles have been upscattered.
We seek to model only upscattering of the most bound DM particles – away from the centre of a dark halo the DM particles have been abundantly scattered to low phase-space densities by the potential fluctuations associated with structure formation. We need to model scattering, by whatever agent, that is not correctly captured by the NFW model. This line of argument motivates us to propose analytic approximations to the current DF of Galactic DM based on the assumption that in action space upscattering has produced a constant-density core, from which the DM density falls monotonically as increases such that it asymptotes to the DM-only form (3) proposed by Posti et al. (2015).
To obtain a satisfactory DF we multiply the DF (3) by a function that varies as for small and tends to unity for large . At intermediate values of , should exceed unity by a small amount to ensure that the total mass of DM is conserved. That is, we require
A suitable functional form for is
where is an arbitrary constant with the dimensions of action that sets the scale of the almost constant-density core of the final DF , and the constant , which controls the magnitude in the peak of before it asymptotes to unity, is determined by requiring satisfaction of equation (4).
The function is promising for present purposes, because (i) when , the denominator is dominated by the first term, so which is what is required to annul the divergence of as , and (ii) when it is evident that as required. The peak of is associated with the minimum of the quantity in square brackets, which occurs when and . Naturally we require to ensure that the quadratic expression on the bottom of equation (5) has no real roots. The full curve in Fig. 1 shows the DF obtained for and .
The NFW halo extends to infinity, but we require a halo in which the density vanishes sufficiently far from the origin to simplify the computation of the the potential generated by the halo. Since we are exclusively interested in the structure of the halo at radii , there is no reason not to truncate the halo at a large radius. This we do by subtracting from with argument the value of evaluated at and declaring the DF to be zero if . That is our final DF is
The constant can be any a large action that prevents the dark halo extending at low density to infinity. Our choice, , has negligible effect on the halo’s density within the radius at which the halo’s density becomes 200 times the mean cosmic density of matter. Consequently, no reported property of the halo would be changed by increasing to arbitrarily large values.
The homogeneous function controls the flattening and velocity anisotropy of the halo. Following BP15 we adopt
where and are given by equations (6) and (7) of P15 with to ensure radial anisotropy. Appendix C of P15 gives the rationale for this choice of the dark halo’s DF. On account of the appearance of ratios of epicycle frequencies in equation (7), the functional dependence on J of our halo DF depends on the model’s potential. When the disc and bulge are introduced, this circumstance is inconvenient, so after we have determined the potential that a halo DF generates in isolation, we freeze the functional forms , and of the frequencies that appear in the definition of . Consequently, in the final model the frequencies of circular orbits are only approximately given by these functions.
Since is determined by equation (4), the free parameters in are , which sets the NFW scale radius, , which sets the size of the dark halo’s core, and the normalisation .
The dark green curve in Fig. 2 shows the density profile that is generated in the plane by the halo DF of one of our final models when the halo is isolated. Clearly this halo has a homogeneous core that extends to . We show for comparison the density profiles of: (i) in red the NFW halo flattened to axis ration that P14 fitted; (ii) in dark blue the dark halo fitted by BP15; (iii) in magenta the halo that the BP15 halo DF generates in isolation. The difference between the dark blue and magenta curves represents the adiabatic deformation of the dark halo by the disc and bulge. The red curve has the same form as the magenta curve because the BP15 DF was one which in isolation generates an NFW profile, and the red curve lies above the magenta curve because the former is a fit to the density of the halo in the presence of the disc rather than in isolation. The dark green curve lies below the other curves because it shows what the current halo would look like if the disc and bulge were to be slowly dismantled, allowing the halo to expand. Hence it is most comparable to the magenta curve. It lies below this curve because our dark halo is less massive than that of BP15, and thus allows the disc to place more mass at , and in this way provide adequate optical depth for microlensing.
3.2 The stellar disc
As in P14 and BP15, the DF of the disc is superposition of the “quasi-isothermal” form that was introduced by Binney & McMillan (2011), namely
where and are
Here and control the radial and vertical velocity dispersions, while , and are, respectively, the circular, radial and vertical epicycle frequencies of the circular orbit with angular momentum . The function
where is the radius of the circular orbit, determines the surface density contributed by the quasi-isothermal. To keep the disc scale-height roughly independent of , we take , where the constant . Although there is no compelling reason to do so, it is customary to give the same dependence on .
The DF of the thick disc is taken to be a single quasi-isothermal, while the thin disc’s DF is built up out of a quasi-isothermal for each coeval cohort of stars.
The DF of the thin disc is taken to be a superposition of quasi-isothermal DFs, one for the stars of each age , and the velocity-dispersion parameters depend on the age of the cohort in addition to . As in earlier papers we assume
where , is a constant that determines the velocity dispersion of stars at birth, is the present age of the disc, and is the current velocity dispersion of the oldest thin-disc stars near the Sun.
We assume that the star-formation rate in the thin disc has decreased exponentially with time, with characteristic time scale , so the complete thin-disc DF is
We set the normalising constant that appears in equation (11) to be the same for both discs and use for the complete DF
where is a parameter that controls the fraction of stars that belong to the thick disc. The values of the parameters for our final model are given in Table 5.
We followed P14 in imposing a lower limit of on the value of at which the epicycle frequencies and are evaluated for use in the DF.
3.3 DF of the stellar halo
As in P14 and BP15, when computing velocity histograms for comparison with the RAVE data, we add to the stellar DF a small contribution from the stellar halo, which we presume to have no net rotation. If a halo population is not included, the DF of the thick disc is distorted to provide some stars that are counter-rotating to the disc, since the data include a few such stars. Since the mass of the stellar halo is negligible, we do not include the DF of the stellar halo when we integrate over velocities to determine the total stellar mass. Since we only require the stellar halo at points near the Sun, it is adequate to adopt from Posti et al. (2015) the DF that generates a power-law density profile , with index (see e.g. Binney & Merrifield, 1998, §10.5.2). Thus we take as the (un-normalised) DF of the stellar halo to be
Here , , and . With these choices the stellar halo is approximately spherical.
In summary, when computing kinematics, the total stellar DF is taken to be
with chosen so at to be consistent with the J08 data as explained in BP15.
3.4 The bulge/bar and gas disc
Our modelling technique restricts us to axisymmetric models, so we cannot use a sophisticated model of the bulge/bar. Moreover, the data we use are only sensitive to the bulge’s contribution to radial forces. Therefore we do not represent the bulge by a DF but by a fixed axisymmetric mass distribution. We have updated our model from that used by BP15, which followed McMillan (2011) and was thus based on Bissantz & Gerhard (2002). Our bulge model is now consistent with Wegg & Gerhard (2013) and Portail et al. (2015). These authors fitted dynamical models of the bulge/bar to near-IR photometry from the Vista Variables of the Via Lactae survey (Saito et al., 2012) and line-of-sight velocity measurements of red-clump stars from Ness et al. (2013). They found a total mass within the cuboid with corners at () to be . Due to the axisymmetry of our model, we cannot reproduce this distribution exactly but we construct our bulge to have the same mass within a similar volume taking into account the mass of the discs and the dark halo. This gives us a total mass of the bulge in this volume of . We use a scale length for our exponential that is the average of the two scale lengths in the and direction (the geometric mean is very similar).
The density distributions of the bulge is
The gas disc is likewise represented by an axisymmetric distribution of matter that has density
A non-zero value of the parameter creates a central cavity in the disc. The values of the parameters are given in Table 2. The surface density normalisation is adjusted to maintain the ratio between the gas and stellar surface densities at that is given in Flynn et al. (2006). and are the only parameters that are varied: the other parameters are fixed at the values adopted by P14 and earlier investigators. As in previous papers of this series and several other studies (e.g. Bovy & Rix, 2013; Wegg, Gerhard & Portail, 2016) we set . This setting may over-estimate . However, the local surface density of the gas disc is narrowly confined by the ratio given above, and we have set , so varying will not have much impact on the structure of our model at , which is what is of concern here.
3.5 Fitting algorithm
The algorithm we use to fit the DF to the data is essentially that used by BP15 – Table 3 lists its steps.
Since the main differences between our model and that of BP15 are in the inner Galaxy (), we fixed to their value, . As they argue, this scale action yields a realistic break radius and has little bearing on the dark halo’s contribution to forces in the inner few kiloparsecs. This leaves and the only adjustable parameters of the dark halo’s DF, one more than the single free parameter, , available to BP15. So for a grid of values of , we use the BP15 algorithm to determine and the nine free parameters in the disc DF (two masses, four pseudo velocity dispersions, the radial mass scale length and two radial scale lengths for the thick disc’s velocity dispersions). For every chosen value of we obtain a model that satisfies all observational constraints other than the microlensing data. As increases, the inner disc becomes steadily more massive so the optical depth to microlensing increases.
Fig. 3 shows for several values of the of the J08 data with respect to the model vertical density profile computed in step 10 of the fitting algorithm during the search for the optimum value of the halo normalisation . It is evident that for a given value of remarkably few values of have been tried, with the consequence that our accepted value of almost certainly differs from the best possible value. The sparseness of the search over reflects the computational cost of optimising the disc DF for a given value of : nine parameters can be adjusted, and after each adjustment the self-consistent potential must be determined. The potential is found by integrating the DF over three velocity components at each node of a grid in the plane. Because the potential is determined iteratively, these integrations have to be executed four or five times for each choice of disc DF. Hence core-hours are required to optimise for a given value of .
Since for some values of the adopted value of will come closer to the true optimum value than for other values of , the properties of the recovered models do not change entirely systematically with . However, the jitter is small enough to reveal systematic trends that would dominate if we could always locate the true optimum value of .
|1||Choose values of the scale action and the parameter that controls the mass of the dark halo.||
Choose a pair of plausible double-exponential density distributions to represent the stellar discs, and choose a plausible gas disc.
Find the self-consistent equilibrium of the chosen dark halo in the presence of the bulge, the gas disc and the adopted models of the stellar discs.
Adjust and in the formulae for the stellar and gas discs until the constraints on listed in Section 2 are satisfied. Each time the disc parameters are changed, return to step 3 until changes become negligible.
Choose a DF for the stellar disc that has the scale lengths found in Step 4 and adopt plausible values for its velocity-dispersion parameters.
, , , ,
At the nodes of a spatial grid, integrate the disc DF over velocities to obtain the contribution to the density from the stellar disc. After adding in the contributions from the gas, bulge and dark halo, solve for the associated potential and adopt this potential.
At the nodes of a spatial grid, integrate the halo DF over velocities to obtain the contribution to the density from the dark halo. After adding in the contributions from the gas, bulge and stellar disc, solve for the associated potential and adopt this potential.
If the update to the potential at Step 7 is non-negligible, return to Step 6. Otherwise proceed to Step 9.
Adjust the velocity-dispersion parameters in the disc DF to obtain a good fit of the model’s kinematics to the kinematics of the RAVE giants (including the contribution of the stellar halo). Then return to Step 6 and continue until changes to the velocity-dispersion parameters become negligible.
Compute the residuals between the vertical stellar density profile of J08 and that implied by the DFs of the disc and stellar halo.
If these residuals are unsatisfactory, choose a new normalisation for the dark halo and return to Step 2. Otherwise, finish.
|Thin and thick disc|
Tables 4 and 5 present, respectively, the parameters of the dark haloes and the stellar discs of six models that are consistent with all the data, including the microlensing optical depths – these are the models with in the range . The last column in each table describes a model, described in Section 6, that has and is not satisfactory. Table 4 shows that in the six satisfactory models, the radius at which the density of the dark halo falls to half its central value increases with from to . The local DM density varies in the range . The value of does not vary systematically with and we believe the width of the given range reflects the uncertainty in the data and has no physical significance. Table 5 shows that the scale lengths of the stellar discs vary in the range , again without systematic dependence on . Hence all satisfactory models have disc scale lengths that are significantly shorter than the value recovered by BP15 – Table 5 gives all the parameters of that disc in the column for . When is significantly non-zero, there is less DM at , so, given that the surface density at the Sun is essentially fixed by the J08 and RAVE data, the disc needs to have a shorter scalelength to compensate. By the same token, the mass of the stellar disc must be larger than when . Indeed, in the new models, the total stellar disc mass lies in the range with the mass inside the solar circle in the range . When we include the bulge and gas disc, the total baryonic mass lies in the range . Hence, at 7.3 to 8.8 per cent our baryonic fractions are significantly higher than that, 4.2 per cent, of BP15.
In the upper panel of Fig. 4 we plot the circular speed curves of several models: the circular-speed curves of the new models lie within the very narrow shaded area that is bounded by black lines. The yellow and orange shaded areas show the circular speeds that are generated by the baryons and DM, respectively, in the new models. Naturally, these individual curves cover wider range than the total, since from the baryons rises with and this rise is largely compensated by a fall in from the dark halo. Within the uncertainties, the contributions to become equal at .
The blue curve shows in the model of P14 that has axis ratio ; it differs insignificantly from the dark shaded region of the new models. The red curve shows from BP15. It rises more steeply at small radii than the other curves but is consistent with them in the region for which we have constraints. The broken curves show the circular speeds generated by the dark halo (blue) and the baryons (dark green) in the BP15 model. They differ strongly from the corresponding curves of the new models, in particular crossing at rather than at .
In the lower panel of Fig. 4 we plot the fraction of the total radial force that is contributed by the baryons. In the new models (dark shaded region) this is naturally larger than in the models of either P14 or BP15. In particular, at the baryons contribute of the radial force in the new models.
The curves in Fig. 5 show the vertical density profiles of the thin and thick discs and the stellar halo above the Sun in the model with . The data points show the J08 data, and we see that they agree well with the summed density profile for this model. The other models fit the J08 data with comparable precision.
The red lines in Fig. 6 show the fit to the kinematics of the RAVE giants (data points) achieved by the model with in the four spatial bins that lie inside . The left and right columns show that the model provides fits of outstanding quality to the distributions or radial and vertical velocities. The fit to the distributions of components is good, if not quite as good as that obtained by BP15 (green lines). The most significant discrepancies are under-provision of stars with large further than from the plane (top central panel). This deficit suggests that the scale length of the thick disc should be longer. In fact, J08 reported for the thick disc and for the thin disc, whereas we have assumed from the outset that the two discs have the same scalelengths.
However, evidence is accumulating that the disc formed by -high stars has a smaller value of than the disc formed by the -low stars (Hayden et al., 2015), so if one identifies the former disc with the thick disc, one has a conflict with J08. The likely resolution of this conflict is that the disc formed by the -low stars flares, so the scalelength of all stars that lie at is in fact larger than the scalelength of the stars that lie at . This would explain the deficit in our models of stars with large revealed by the central top panel of Fig. 6.
A significant, but less serious, problem in Fig. 6 is under-provision of stars with low in the central top panel. This problem arises also in the model of BP15 and signals weakness in our choice of the DF of the stellar halo.
Our model provides a better overall fit to the data further from the plane (upper two centre panels) but has an excess of non-rotating stars at low values of .
In Fig. 7 the dark shaded region shows the range of values of the vertical force provided by the new models. This lies slightly above the curve from P14 (blue curve) because the new discs have slightly higher stellar surface densities at ( Table 5) than the P14 disc (). Nonetheless, the new models are consistent with the classic determination of Kuijken & Gilmore (1991). They are also consistent with the findings of Bienaymé et al. (2014), indicated by green boundary curves.
4.1 Structure of the dark halo
The green curves in Fig. 8 show the density profiles in the equatorial plane of the model with (i) in the full model (full curve), and (ii) after removal of the baryons (dashed curve). We see that the gravitational field generated by the baryons increases the halo density by an order of magnitude in the region , but changes it very little at . In the full model the density of the dark halo at is almost the same as in the P14 model (red curve). The density of the dark halo falls well below that of the P14 model only at .
At the density profile of our halo rises more steeply than that of an NFW halo because the inward pull of the baryons’ gravitational field is centrally concentrated. Thus at the two green curves Fig. 8 asymptote to one another and lie below but parallel to the red curve of an NFW halo. On account of the inward pull of the baryons, the full green curve almost touches the red curve around because the density of DM at is tightly constrained by the data, and the density everywhere else follows from this density and the structure of our chosen DF.
The third row of Table 4 gives the values of the mass, , of DM inside the radius where the halo’s mean density is 200 times the mean cosmic density of matter. At these values are significantly less than the value, for the P14 model and that, , found by BP15. This reduction in is a straightfoward consequence of the dynamics described in the previous paragraph. However, there are precedents for such light dark haloes: Peñarrubia et al. (2014) estimate the Milky Way’s dark halo mass as and a recent estimate by Huang et al. (2016) gives a halo mass of .
Our haloes do yield local DM densities that are higher than many encountered in the literature (see the review of Read, 2014). In particular, the dark haloes of Wegg, Gerhard & Portail (2016) yield local densities in the range . However, from the fact that Bienaymé et al. (2014) derived a very similar value to us from a completely different analysis of the RAVE data, we conclude that our large local DM density follows from the unprecedented RAVE data and supersedes earlier estimates.
5 Optical depths to microlensing
In this section we investigate which models are consistent with the
microlensing data. The latter were for many years subject to a regrettable
level of uncertainty on account of the phenomenon of ‘blending’: stars that
will be lensed merging with brighter stars with the result that one
under-estimates the number of stars being monitored and thus over-estimates
the optical depth to lensing (Woźniak & Paczyński, 1997). Recognition of this phenomenon
shifted attention to the statistics for lensing of apparently bright stars,
such as red clump stars, which are less affected by blending
(Popowski et al., 2001, 2005). However, it now
seems that blending is well enough understood for valid optical depths to be
determined from the lensing statistics of all stars, not just the brightest
ones (Sumi & Penny, 2016). In particular, Sumi et al., analysing data from the
Wegg, Gerhard & Portail (2016) have recently evaluated the optical depth to microlensing predicted by their Galaxy models. Inside these models comprise an N-body model from Portail et al. (2015), while at the models comprise a double-exponential stellar disc – the latter consists of a thin disc only with local surface density . In some cases this disc has scale height , and in others it flares from at to at . Given that our bulge is an axi-symmetrised version of a bulge/bar from Portail et al. (2015), these models are very similar to ours.
When viewed close to its major axis, a barred Galaxy will have a larger optical depth to microlensing, , than the axisymmetric model with the same circular-speed curve. As Wegg, Gerhard & Portail (2016) point out, the angle, , between the long axis of the bar and the Sun-Centre line is such that we expect an axisymmetric model to under-estimate only slightly. Most lenses will lie in the region where Wegg, Gerhard & Portail (2016) assume axisymmetry, so their optical depths will be under-estimated to a similar extent to ours.
Wegg, Gerhard & Portail (2016) compute for each of the MOA-II fields taking into account the luminosity function of source stars. Together with the apparent magnitude limit of the survey, the latter determines the distribution along the line of sight of the surveyed stars, and therefore the mean optical depth in the given field. Here we perform a much simpler calculation: we compute
for a star that is located at . Here is the distance to the source, is the distance to the deflector and is the mass density contributed by the deflectors. The value of from equation (21) is typical of what the more sophisticated computation yields after averaging over sources fainter than (see Wegg, Gerhard & Portail, 2016, Fig. 1).
The curves in Fig. 9 show the predicted optical depths. As expected, at any given the optical depth increases with . As increases the decline in with becomes less steep, and the steepest curves are most consistent with the data points. The curve for is certainly too shallow.
BP15 argued that their model was inconsistent with the measured optical depth to microlensing of bulge stars because in it the mass of stars inside was less than the considered essential by Binney & Evans (2001). Our bulge model has at , so by the criterion in Binney & Evans, the stellar disc should contain within . In the new models the mass of the stellar disc within ranges from 3.02 to within of the galactic plane and between 2.55 and within of the galactic plane, so they just about satisfy the criterion of Binney and Evans.
6 An over-large core
A minimum value of is set by the requirement for sufficient microlensing optical depth. Now we show what happens when is made too large.
A large value of implies a large core radius of the dark halo. A compact disc is then required to keep the circular speed quite flat given a radially rising contribution from the dark halo. Since the disc’s surface density at is constrained by the RAVE data, a disc with a small value of is a massive disc. The final column of Table 5 quantifies these points by showing that yields and . The final column of Table 4 shows that with with the local DM density is only yet the reference mass is more than 50 per cent higher than in the acceptable models. With a large core in the halo, a sufficiently large value of can only be achieved alongside a large density at . For this reason, is depressed by large .
Fig. 10 compares the stellar density profiles above the Sun in an acceptable model () and the model with over-large . The full curve for the model with over-large lies above that of the acceptable model. This is not a problem because the data points can be shifted up or down by increasing or decreasing the mass-to-light ratio of the stars. What is an issue is that the full curve is steeper at because the model with over-large has the a more massive disc. As a consequence of this change of shape, by sliding the data points vertically, they can be brought into closer agreement with the broken than with the full curve. The black curve in Fig. 11 shows that the massive thin disc of the model with generates values of that are significantly larger than most estimates.
In the CDM paradigm, structure formation commences with DM at an essentially infinite phase-space density. As the DM aggregates gravitationally, the fluctuating gravitational potential that the DM particles experience, scatters particles to lower phase-space densities. Remarkably, simulations of this process revealed that extremely high, possibly infinite phase-space densities of DM survive at the centres of dark haloes. No analytic function provides a perfect representation of the central structures of the haloes that form in DM-only simulations, but the NFW profile and the Einasto profile (Einasto, 1965) both provide accurate fits (e.g. Duffy et al., 2008). Whereas the NFW profile requires an infinite central density, as the phase-space density associated with the Einasto profile rises steeply to finite value. Even in the Einasto model, however, in no part of action space is the phase-space density essentially constant. We consider this a weakness because on general physical grounds, scattering will create such a region around any peak in the phase-space density.
Numerical experiments have shown that dark haloes have cuspy cores regardless of whether they form primarily through violent relaxation of single large structures, or through repeated merging of many small haloes (Huss, Jain & Steinmetz, 1999; Moore et al., 1999). This finding implies that when a small halo is cannibalised by a larger halo, particles released at high phase-space density by the small halo rather precisely replace in action space the particles of the larger halo that have been kicked up to higher actions.
El-Zant, Shlosman & Hoffman (2001) pointed out that when a baryonic structure falls into a dark halo, the scattering of the host halo’s particles to regions of lower phase-space density is not compensated by release of DM particles by the inspiralling structure, and the central cusp of the host halo will be eliminated. Jardel & Sellwood (2009) quantified this proposal using N-body simulations of spherical systems, and reported that the effect is weaker than El-Zant, Shlosman & Hoffman (2001) predicted. They did find, however, that when is distributed amongst 500 clumps that initially orbit in the halo out to , after a core develops that extends out to , which corresponds to in the case of our Galaxy. Reducing the number of clumps from 500 causes the scale of the core to grow. Cole, Dehnen & Wilkinson (2011) investigated clumps on a variety of initial orbits and found that a clump with mass 1 per cent of the mass of the halo could remove about twice its own mass from the inner halo, transforming a cusp into a core or weak cusp. This effect was enhanced if the clump was subsequently removed, as by a galactic wind.
Hernquist & Weinberg (1992) pointed out that a bar rotating within a spheroidal component will experience dynamical friction as a consequence of upscattering the spheroid’s particles. Subsequent N-body simulations of bars in spheroids (Debattista & Sellwood, 2000; Athanassoula, 2003) have shown that this is an important effect: the bar responds to the loss of energy and angular momentum to the spheroid by becoming longer, slower and stronger. The response of the spheroid to the bar is less well documented. Weinberg & Katz (2002) showed that a bar would tend to produce a core in a dark halo. Holley-Bockelmann, Weinberg & Katz (2005) confirmed that a bar will produce at least a small core in its dark halo, but Sellwood (2003) showed that at small radii the DM density nevertheless rises as the disc forms because at most radii compression by the growing gravitational field of the baryons overwhelms the effect of core formation. This finding underlines the important role played by galactic winds in expelling baryons from the central galaxy after then have upscattered DM.
We believe our proposal is consistent with all these studies. In particular, our favoured core radii are much smaller than and, as Fig. 8 illustrates, upscattering has probably not fully compensated for adiabatic compression through most of the star-forming disc.
Binney, Gerhard & Silk (2001) argued that massive outflows from early galaxies indicate that the centres of galaxies have processed significantly more baryons than now reside in visible galaxies. As baryons sank towards the centres of a dark halo, they probably surrendered a considerable amount of angular momentum to the dark halo, causing it to expand. Later stellar feedback pushed the baryons out into the circum/intergalactic medium, where they now reside. The case for ejection of large quantities of baryons is now established but the impact of these baryons on the dark halo remains uncertain. In a similar spirit Nipoti & Binney (2015) argued that early in the life of a halo, baryons are likely to accumulate in the form of gas until their gravitational field has comparable strength to that of the DM. At that point the body of gas in a low-mass dark halo is likely to break up into a small number of blobs, each of which contains a significant fraction of the total mass interior to its orbit. In these circumstances, a blob will dump its orbital energy into DM within a couple of dynamical times. The take-up of this energy will significantly reduce the central phase-space density of DM.
Navarro, Eke & Frenk (1996) investigated the effect of baryon ejection on a dark halo in the axisymmetric limit and found it to be a significant process, but Gnedin & Zhao (2002) concluded from a study of the spherical limit that the effect is too weak to be of interest. In our view angular-momentum transfer, which is suppressed in the axisymmetric case and eliminated in the spherical case, lies at the core of the process, so these studies are of marginal value.
We conclude that some combination of a bar and massive baryonic lumps are likely capable of smoothing the cusp in the phase-space density at the centre of a dark halo out to as our analysis suggest has happened in our Galaxy.
BP15 showed that an NFW DF is inconsistent with data for our Galaxy, presumably because the most tightly bound DM particles have ben scattered to a nearly constant phase-space density by fluctuations in the overall gravitational potential. Hence we have here introduced a three-parameter family of analytic functions as candidates for the DFs of dark haloes in galaxies like ours. One parameter controls the halo’s scale radius , another controls its mass, and the third, , controls the size of the region around the origin of action space in which the phase-space density of DM is almost constant. The existence of this region distinguishes these haloes from standard NFW haloes, in which the phase-space density diverges as , and from Einasto haloes, in which the phase-space density rises to a sharp peak as .
We have searched the family of new halo DFs and our usual family of stellar DFs for sets of parameters that self-consistently generate a Galaxy model that is consistent with constraints on the rotation curve, the vertical structure of the stellar disc in the solar cylinder and the kinematics of giants in RAVE. There is a range of values of for which these constraints are satisfied and optical depths to microlensing of bulge stars are predicted that are consistent with the results of surveys for microlensing events. In the new models the local DM density  and stellar surface density  are very similar to the values determined by P14 and BP15, but with the new discs are more massive than those of BP15 because the new models have shorter disc scale lengths () than the BP15 model. The haloes have core radii . The microlensing data are critical for these fits because they place a lower limit on the contribution made by stars to the Galaxy’s central gravitational field.
A concern we had regarding the earlier work by P14 and BP15, was that these papers assumed that the stars that entered the RAVE survey are not kinematically biased as a result of some bias in the selection of stars for spectroscopy. A particular worry was the possibility of bias towards younger stars, which would be kinematically cooler. A detailed study of the selection function of the RAVE survey by Wojno & the RAVE collaboration (2016) has now established that the kinematics of the RAVE stars used in BP15 and here is remarkably representative of the underlying population. Specifically, Wojno & the RAVE collaboration (2016) used the tool Galaxia (Sharma et al., 2011) to form a mock-RAVE sample by applying the RAVE selection function to a modified Besançon Galaxy model (Robin et al., 2003). Then they compared the kinematics of the selected stars to those of the underlying population and found that both the giants and the main-sequence stars had essentially identical kinematics in the mock-RAVE sample and the underlying population. Wojno & the RAVE collaboration (2016) find that RAVE’s selection function does introduce a small kinematic bias to the stars in the turn-off region of the colour-magnitude diagram. In future work we will extend the present work to dwarf stars, introducing the appropriate bias to the model before comparing with the data. We will also update the input data to RAVE DR5 (Kunder & the RAVE collaboration, 2016) with the parallaxes and proper motions updated to the add recently released data from Gaia (Gaia Collaboration et al., 2016).
We have not obtained the posterior probability distribution in model space. In fact, as Fig. 3 suggests, we have probably not even identified the most probable model. Our searches of model space have been inadequate because the construction of models specified by DFs and have self-consistent gravitational potentials is computationally expensive. The expense arises from the need to compute the density, by means of a triple integral over velocity, throughout the huge volume occupied by our Galaxy. As a result of software improvements since this work commenced, the cost of model construction has been reduced by about an order of magnitude, and in future it will be possible to search model space more thoroughly.
We thank an anonymous referee for comments that enabled us to improve this paper. 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.
Funding for RAVE has been provided by: the Australian Astronomical Observatory; the Leibniz-Institut für Astrophysik Potsdam (AIP); the Australian National University; the Australian Research Council; the French National Research Agency; the German Research Foundation (SPP 1177 and SFB 881); the European Research Council (ERC-StG 240271 Galactica); the Istituto Nazionale di Astrofisica at Padova; The Johns Hopkins University; the National Science Foundation of the USA (AST-0908326); the W. M. Keck foundation; the Macquarie University; the Netherlands Research School for Astronomy; the Natural Sciences and Engineering Research Council of Canada; the Slovenian Research Agency; the Swiss National Science Foundation; the Science & Technology Facilities Council of the UK; Opticon; Strasbourg Observatory; and the Universities of Groningen, Heidelberg and Sydney. The RAVE web site is at http://www.rave-survey.org.
- Amorisco N. C., Agnello A., Evans N. W., 2013, MNRAS, 429, L89
- Athanassoula E., 2003, MNRAS, 341, 1179
- Bienaymé O. et al., 2014, A&A, 571, A92
- Binney J., Gerhard O., Silk J., 2001, MNRAS, 321, 471
- Binney J., Gerhard O., Spergel D., 1997, MNRAS, 288, 365
- Binney J., McMillan P., 2011, MNRAS, 413, 1889
- Binney J., Merrifield M., 1998, Galactic Astronomy. Princeton University Press
- Binney J., Piffl T., 2015, MNRAS, 454, 3653
- Binney J. J., Evans N. W., 2001, MNRAS, 327, L27
- Bissantz N., Gerhard O., 2002, MNRAS, 330, 591
- Blumenthal G. R., Faber S. M., Primack J. R., Rees M. J., 1984, Nature, 311, 517
- Bovy J., Rix H.-W., 2013, ApJ, 779, 115
- Cole D. R., Dehnen W., Wilkinson M. I., 2011, MNRAS, 416, 1118
- Debattista V. P., Sellwood J. A., 2000, ApJ, 543, 704
- Dehnen W., Binney J., 1998, MNRAS, 294, 429
- Duffy A. R., Schaye J., Kay S. T., Dalla Vecchia C., 2008, MNRAS, 390, L64
- Einasto J., 1965, Trudy Astrofizicheskogo Instituta Alma-Ata, 5, 87
- El-Zant A., Shlosman I., Hoffman Y., 2001, ApJ, 560, 636
- Flynn C., Holmberg J., Portinari L., Fuchs B., Jahreiß H., 2006, MNRAS, 372, 1149
- Gaia Collaboration, Brown A. G. A., Vallenari A., Prusti T., de Bruijne J., Mignard F., Drimmel R., co-authors ., 2016, ArXiv e-prints
- Gnedin O. Y., Zhao H., 2002, MNRAS, 333, 299
- Hayden M. R. et al., 2015, ApJ, 808, 132
- Hernquist L., Weinberg M. D., 1992, ApJ, 400, 80
- Holley-Bockelmann K., Weinberg M., Katz N., 2005, MNRAS, 363, 991
- Huang Y. et al., 2016, ArXiv e-prints
- Huss A., Jain B., Steinmetz M., 1999, ApJ, 517, 64
- Jardel J. R., Sellwood J. A., 2009, ApJ, 691, 1300
- Jurić M. et al., 2008, ApJ, 673, 864
- Katz H., McGaugh S. S., Sellwood J. A., de Blok W. J. G., 2014, MNRAS, 439, 1897
- Kordopatis G. et al., 2013, AJ, 146, 134
- Kuijken K., Gilmore G., 1991, ApJL, 367, L9
- Kunder A., the RAVE collaboration, 2016, arXiv, 000, 000
- Malhotra S., 1995, ApJ, 448, 138
- McMillan P. J., 2011, MNRAS, 414, 2446
- McMillan P. J., Binney J. J., 2010, MNRAS, 402, 934
- Moore B., Quinn T., Governato F., Stadel J., Lake G., 1999, MNRAS, 310, 1147
- Navarro J. F., Eke V. R., Frenk C. S., 1996, MNRAS, 283, L72
- Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
- Ness M. et al., 2013, MNRAS, 432, 2092
- Nipoti C., Binney J., 2015, MNRAS, 446, 1820
- Peñarrubia J., Ma Y.-Z., Walker M. G., McConnachie A., 2014, MNRAS, 443, 2204
- Piffl T. et al., 2014, MNRAS, 445, 3133
- Piffl T., Penoyre Z., Binney J., 2015, MNRAS, 451, 639
- Popowski P. et al., 2001, in Astronomical Society of the Pacific Conference Series, Vol. 239, Microlensing 2000: A New Era of Microlensing Astrophysics, Menzies J. W., Sackett P. D., eds., p. 244
- Popowski P. et al., 2005, ApJ, 631, 879
- Portail M., Wegg C., Gerhard O., Martinez-Valpuesta I., 2015, MNRAS, 448, 713
- Posti L., Binney J., Nipoti C., Ciotti L., 2015, MNRAS, 447, 3060
- Read J. I., 2014, Journal of Physics G Nuclear Physics, 41, 063101
- Reid M. J., Brunthaler A., 2004, ApJ, 616, 872
- Reid M. J. et al., 2014, ApJ, 783, 130
- Robin A. C., Reylé C., Derrière S., Picaud S., 2003, A&A, 409, 523
- Saito R. K. et al., 2012, A&A, 537, A107
- Schönrich R., 2012, MNRAS, 427, 274
- Schönrich R., Binney J., Dehnen W., 2010, MNRAS, 403, 1829
- Sellwood J. A., 2003, ApJ, 587, 638
- Sharma S., Bland-Hawthorn J., Johnston K. V., Binney J., 2011, ApJ, 730, 3
- Steinmetz M. et al., 2006, AJ, 132, 1645
- Sumi T., Penny M. T., 2016, ArXiv e-prints
- van der Kruit P. C., Shostak G. S., 1984, A&A, 134, 258
- Wegg C., Gerhard O., 2013, MNRAS, 435, 1874
- Wegg C., Gerhard O., Portail M., 2016, MNRAS
- Weinberg M. D., Katz N., 2002, ApJ, 580, 627
- Wojno J., the RAVE collaboration, 2016, arXiv, 000, 000
- Woźniak P., Paczyński B., 1997, ApJ, 487, 55