Bringing the Galaxy’s dark halo to life
We present a new method to construct fully self-consistent equilibrium models
of multi-component disc galaxies similar to the Milky Way. We define
distribution functions for the stellar disc and dark halo that depend on phase space
position only through action coordinates. We then use an iterative
approach to find the corresponding gravitational potential. We study the
adiabatic response of the initially spherical dark halo to the introduction
of the baryonic component and find that the halo flattens in its inner
regions with final minor-major axis ratios = 0.75 – 0.95. The extent of
the flattening depends on the velocity structure of the halo particles with
radially biased models exhibiting a stronger response. In this latter case,
which is according to cosmological simulations the most likely one, the new density
structure resembles a “dark disc” superimposed on a spherical halo. We
discuss the implications of these results for our recent estimate of the
local dark matter density.
The velocity distribution of the dark-matter particles near the Sun is very non-Gaussian. All three principal velocity dispersions are boosted as the halo contracts, and at low velocities a plateau develops in the distribution of . For models similar to a state-of-the-art Galaxy model we find velocity dispersions around for and the tangential velocity, , and 150 – for the in-plane radial velocity, , depending on the anisotropy of the model.
keywords:galaxies: kinematics and dynamics
It is now generally accepted that most of the rest mass in the Universe is contained in a type of “dark matter” that is incapable of electromagnetic interactions. Dark matter began to cluster gravitationally earlier than baryonic matter because it was decoupled from the cosmic radiation bath, which at redshifts provided abundant pressure. Hence the formation of a galaxy starts with the accumulation of dark matter to form the dark halo. When and how baryons accumulated in the gravitational potential wells of dark halos is ill understood. In the simplest picture, they simply fell in dissipatively and formed stars near each halo’s centre, but this picture is inconsistent with the small fraction of baryons that are now in stars rather than intergalactic gas. Several lines of evidence indicate that the release of nuclear energy by early stars strongly attenuated the infall of baryons, thus increasing the extent to which galaxies are dark-matter (dm) dominated.
The least dm dominated galaxies are ones with luminosities similar to that, , of our Galaxy: dwarf galaxies are very strongly dm dominated and stars play such a subsidiary role in the most massive dm halos, those of galaxy groups and clusters, that we do not even recognise these halos as galaxies. Yet even in galaxies such as ours, it may be argued that the baryons have played a very subordinate role. Measurements of the optical depth to microlensing of stars located less than from the Galactic centre imply that a significant majority of the mass inside is concentrated in stars (Binney & Evans, 2001; Bissantz et al., 2004). But the majority of the gravitational force that holds the Sun in its orbit comes from dark matter (e.g. Piffl et al., 2014, hereafter P14), and the dominance of dark matter increases sharply with increasing distance from the centre.
Our understanding of how baryons accumulated at the centres of dark halos to form visible galaxies is limited. This is because the physics of this process is complex and is known to involve scales that are by far unresolved in feasible simulations of the formation of galaxies in a cosmological context. Early simulations produced galaxies that contained too large a fraction of the baryons in the universe, and were much more dominated by their spheroidal components than is our Galaxy or our near neighbour M33. In recent years it has become clear that energy released during star formation is very efficient at preventing the accumulation of large masses of gas and stars in the centres of dark halos, by expelling gas that has fallen in to a halo and the products of stellar evolution into intergalactic space, where the great majority of the baryons must reside. The high efficiency of such “feedback” from star formation is surprising and it can be reproduced in simulations only by modifying the encoded laws in essentially arbitrary ways.
Since we cannot trust current simulations of the behaviour of baryons during star formation, it is worthwhile to explore simple limiting cases. One such limiting case is that in which the baryons have accumulated in the Galaxy gradually rather than in lumps. This scenario is favoured by four observations. The first is that the bulge of our Galaxy bears all the hallmarks of having formed through a thin disc of stars that formed on nearly circular orbits in the Galactic plane developing a strong bar, which then buckled (Combes & Sanders, 1981; Raha et al., 1991; Ness et al., 2013). The second relevant observation is that in the last the rate of star formation in the disc has declined by no more than a factor (e.g. Aumer & Binney, 2009). A third relevant observation is that the star-formation rate in a disc appears to be proportional to the surface density of cold, molecular gas (Leroy et al., 2008). A fourth observation is that currently of the baryons in the disc are in gas.
From the indication that all but a tiny fraction of the Galaxy’s stars formed from quiescently orbiting gas in the Galactic plane, we infer that our Galaxy has not experienced a significant merger since it started seriously forming stars. From the sustained rate of star formation over the Galaxy’s life together with the dependence of star formation rate on cold-gas density and the relatively small current stock of gas, we infer that the baryonic mass of our Galaxy has accumulated over its entire life, presumably through sustained cooling of gas from the extended corona into the disc. There is even strong circumstantial evidence of this cooling process taking place now (Marasco et al., 2013).
A dark halo will respond to the quiescent accumulation of baryons in a disc by distorting adiabatically: that is, the orbit of each dm particle will evolve in such a way that its action integrals , and will be constant, with the consequence that the density of dm particles in action space, , is unaffected by the accumulation of baryons and the resulting evolution of the galaxy’s gravitational potential (e.g. Binney & Tremaine, 2008, §4.6.1). The constancy of is enormously useful because the gravitational clustering of dark matter in isolation can be simulated with confidence, so we know what dark halos would now populate the Universe if there were no baryons. That is, from dm only simulations we can infer what would be in the absence of baryons, and the adiabatic principle assures us that will be unchanged by the quiescent infall of baryons in the observed quantities. Hence the hypothesis that the baryons accumulated quiescently so the dark matter responded adiabatically to their arrival, allows us to predict the present structures of galaxies before we understand the complex processes that resulted in the accumulation of baryons.
This general idea was appreciated at an early stage in the development of the cold dark matter (CDM) cosmological paradigm (Blumenthal et al., 1986), but until recently we have lacked the technology required to exploit it fully. The key to its exploitation is the ability to compute the action vector as a function of ordinary phase-space coordinates , for then it is possible to compute the dm density
at any location .
In an axisymmetric potential one action, , is simply the component of angular momentum parallel to the potential’s symmetry axis, and in a spherical potential , where is the length of the angular-momentum vector . Blumenthal et al. (1986) estimated the response of a dark halo to the accumulation of baryons by assuming that all dm particles are on circular orbits, so . A more realistic estimate of halo response was obtained by Sellwood & McGaugh (2005), who lifted the restriction to circular orbits but retained the assumption of spherical symmetry. The latter is awkward because the basic picture is of baryons accumulating in a disc, which will immediately break any pre-existing spherical symmetry of the dark halo.
In this paper we lift the assumption of spherical symmetry, so we are able to compute how the quiescent introduction of the disc would have deformed an initially spherical dark halo into a more centrally concentrated, oblate structure. We are, moreover, able to predict the full, three-dimensional distribution of the velocities of dm particles at any point in the Galaxy, in particular at the Sun. These distributions are triaxial and give non-Maxwellians speed distributions.
Our strategy is to adopt dfs for the Galaxy’s stellar disc and dark halo and to solve for the gravitational potential that these self-consistently generate in the presence of given density distributions for the bulge and gas disc. We now describe our algorithm for the determination of , which an extension of the methodology used by (Binney, 2014, hereafter B14) to compute the observables of a flattened isochrone from an adopted df .
2.1 Iterative scheme
B14 started from a guess at the self-consistent potential. Using this guess and equation (1) he determined the density on a grid in the meridional plane, and then solved Poisson’s equation for the corresponding potential . From this he made an improved guess
at the self-consistent potential, where is parameter chosen to maximise the speed of convergence of the sequence of potentials to the self-consistent potential.
The procedure just described involves a grid of points in the plane at which the density is evaluated by integrating the df over . Since the spatial scales of the disc and the dark halo are very discrepant, it is not cost-effective to use the same grid for both components: the grid would have to be fine enough to resolve the structure of the disc and extensive enough to cover the dark halo, which has a significant amount of mass from the Galactic centre. Hence for each component we use a different grid, and on this grid we store estimates of the density of that component and the contribution that it makes to the total gravitational potential111We actually store coefficients of Legendre-polynomial expansions and store also radial derivatives.
Evaluation of the density of a single component involves evaluation of the actions and in this evaluation we use and the “Stäckel Fudge” of Binney (2012a). Having evaluated the density of a component throughout the component’s grid, we solve Poisson’s equation for the corresponding potential . The values of at an arbitrary point is obtained by interpolation on the values stored on the component’s grid.
We start by determining the self-consistent density and potential implied by the adopted dark-halo df when the halo is isolated. Then our first trial potential for the whole Galaxy, , is the sum of this potential and the potentials of the bulge, gas disc and double-exponential stellar disc. We now compute the density of the dark halo in , and solve for the corresponding potential and use this to update . Then we compute the density of the stellar disc in this estimate of and use the corresponding potential to update . Then we compute the density of the dark halo in the updated , and so on until successive estimates of differ negligibly. This occurs after no more than four computations of the density of each component.
2.2 Reference model
P14 modelled our Galaxy by combining a spheroidal mass distribution that has the classic NFW radial profile (Navarro et al., 1996) and a similar spheroidal bulge with a fixed gas disc and a dynamical stellar disc that has a df of the type introduced by Binney (2012b). The parameters of these components were adjusted until the composite model reproduced a variety of observational constraints, such as tangent velocities extracted from gas kinematics at various longitudes, the proper motion of Sgr A* and maser sources in the disc, and the kinematics of stars with spectra taken in the RAdial Velocity Experiment (RAVE Steinmetz et al., 2006; Kordopatis et al., 2013). Crucially, the stellar disc was also required to be consistent with the density profile above the Sun determined by Gilmore & Reid (1983) and Jurić et al. (2008).
P14 found that the most important unconstrained parameter in their models was the axis ratio of the dark halo. The smaller is, the greater the fraction of the force on the Sun that comes from dark matter rather than stars. Since the potentials of the discs cause a dark halo that is spherical when isolated to flatten in its inner region, we have focused on the model in P14 that has a mildly flattened dark halo (minor-major axis ratio ). Moreover, P14 remark that comparison of their results with those of Bienaymé et al. (2014) favours a flattened model with axis ratio . We shall refer to this model as our ‘reference model’. It was described by P14 (cf. their figures 9, 10 and 11), but they did not specify its parameters. Appendix A gives the functional forms used for the density distributions of the model’s components. Table 4 contains the values of the parameters that appear in these forms.
Many previous authors have, like P14, assumed that the dark halo will have the NFW radial profile. But this is the profile of dark halos in dm only simulations. Here we explore the extent to which the dark halo of the reference model would be modified if the discs determined by P14 were adiabatically inserted into it.
3 The DFs
3.1 The stellar disc
The functional form of the stellar disc’s df, which is made up of contributions from the thick disc and each coeval cohort of thin-disc stars, was given by P14 and is reproduced in Appendix B. Table 5 gives the values of the parameters that appear in the disc’s df.
We did make one minor change to the functional form of the df: we imposed a lower limit, , on the circular radius at which the epicycle frequencies are evaluated. Implicitly this introduces an upper limit on the frequencies , etc., that appear in the df. Crucially, after the limit has been imposed these functions remain well-defined functions of the actions, so the legitimacy of the df is not undermined. But on account of the limit, the frequencies employed in the df differ from the epicycle frequencies at low . Imposing the limit has negligible impact on the structure of the disc near the Sun; its material impact is at small radii, where the steep rise in the epicycle frequencies as the Galactic centre is approached would otherwise depress the disc’s density in an unphysical way. Table 5 gives the values of the parameters in the disc df.
3.2 The stellar halo
P14 used a df also for the stellar halo, but we omit this component because it has negligible mass and served only to improve fits to the stellar kinematics at small , which we do not re-examine here.
3.3 The dark halo
Pontzen & Governato (2013) have fitted dfs to dark halos formed in dm only cosmological simulations. Unfortunately, the dfs they fitted are not usable here because they have the form : a df that describes one component of a composite model with a self-consistently generated potential cannot have energy in its argument list for the following reason. When components are added, the central potential becomes deeper, so the energy of the star that rests at a component’s centre decreases. If appears in the argument list, the phase space density around this orbit, , changes in an uncontrolled way. If increases with as it generally does, the component’s density and mass increases, so increases, and without explicit mass renormalisation, the component’s mass is likely to run away to infinity during successive re-evaluations of the potential. When, by contrast, is taken to a function of just the actions, the phase-space density around each orbit never varies and the iterated potentials converge rapidly.
In dm only simulations dark halos have NFW density profiles (Navarro et al., 1996). Posti et al. (2015) found a simple df that self-consistently generates a spherical model that has almost exactly an NFW profile. This model is essentially isotropic near its centre, and becomes mildly radially biased beyond its scale radius. For the dark halo we adopt dfs that are based on the form proposed by Posti et al. (2015) even though these dfs do not generate an NFW profile in the presence of the stellar disc. The normalisation of the dark-halo df (and thus the halo’s total mass) is chosen to match the circular speed, , in the reference model at the solar radius .
We now specify three different dark-halo dfs that cover the possibilities that the velocity distribution of dm particles is isotropic or radially or tangentially biased. We write
where is a normalisation constant, and are constants that define the exponents of inner and outer power-law slopes of the density profile, is a constant that controls the transition between the two. The constant is a small number required to keep the central density finite. Finally,
is an (approximately) homogeneous function222Posti et al. (2015) argued that should be a homogeneous function of , i.e. they required that . of of order unity with
and . In Appendix C we explain why realistic velocity distributions are obtained only when the values and are functions of the actions.
The parameter controls the model’s velocity anisotropy. For we obtain a near-isotropic model, while values below (above) unity yield tangentially (radially) biased models. The quantities and are epicycle and azimuthal frequencies in the self-consistent spherical model evaluated at the radius of a circular orbit with angular momentum
We take the argument of to be to make it an approximate function of energy, so does not become small, and the epicycle frequencies large, for stars on eccentric and/or highly inclined orbits with small .
Following Posti et al. (2015), we set the exponents and to 5/3 and 2.9, respectively, to obtain a good approximation to the NFW density profile. Table 1 lists all the adopted values of the parameters that define the halo’s df. They generate a halo that (a) closely approximates the NFW profile when the halo is isolated, and (b) in the presence of the discs and bulge has a flattening interior to the Sun similar to that of the reference model’s dark halo.
|Parameter||Radially biased||Isotropic||Tangentially biased|
Once we have obtained a self-consistent model of an isolated halo, we freeze the dependence of on its argument, so while we relax onto the potential that is jointly generated by all the Galaxy’s components, the df stays exactly the same function of . It is essential to freeze the function during the introduction of the discs and the bulge if one seeks to learn how the halo is distorted by the gravitational fields of its companions.
3.4 The bulge and gas disc
We could adopt a df for the final Galactic component, the bulge, but we do not do so for two reasons. First, the bulge contains a rotating bar and at present we cannot represent such an object with a df (but see Sanders & Binney, 2015a, for a non-rotating triaxial example). Second, we are primarily interested in the structure of the dark halo in the vicinity of the Sun, which should be well recovered using an axisymmetrised form of the bulge. Hence we represent the bulge by the fixed density distribution specified in Appendix A. We also represent the Galaxy’s gas disc by a fixed density distribution.
4.1 Isolated halos
We consider three halo dfs, which differ only in the value of the parameter : we set to generate a tangentially biased model, for a nearly isotropic model, and for a radially biased model. In all plots we distinguish these models through the colour scheme orange isotropic, green radially biased, and purple tangentially biased. The lower panel of Fig. 1 shows the density profiles generated by the three dfs in the analytic NFW potential (which has no core). The df parameters were optimized to provide a good match between these density profiles and that of the NFW halo fitted by P14. The upper panel of Fig. 1 shows the radial density profiles of these models after the potential has relaxed to self-consistency. (Fitting a self-consistent df-potential pair to the analytic profile directly is conceptually straightforward, but computationally expensive, and is beyond the scope of this paper.)
The thin curves in Fig. 2 show the anisotropy parameter
for the initial isolated halo as a function of radius in the plane. To provide some context, Fig. 3 shows as a function of the slope of the logarithmic density profile for two of these models. The grey band shows the range of such dependency found by Ludlow et al. (2011) from cosmological simulations, while the black line shows the linear dependence that Hansen & Moore (2006) found described cosmologically simulated halos well. We see that our isotropic and radially biased models bracket the empirical results well, except possibly at the smallest values of , which occur at the centre of a halo.
4.2 Halo distortion
The introduction of the baryonic components causes the initially spherical dark halo to contract. The end point of this contraction should coincide perfectly with what we would have found had we slowly grown the baryonic components in a live N-body simulation of the halo. The almost coincident heavy curves in the upper panel of Fig. 4 show the final, relaxed dm density profiles in the disc () plane. Even though the central parts of the initial profiles differ, with the radially-biased model being most cuspy (light full curves in the upper panel), the final profiles are very similar. The full curves in the lower panel of Fig. 4 show the ratio of initial to final density at each radius. We see that in all three cases (radial bias, isotropic, tangential bias), the central density increases by more than a factor 8.
This figure shows (i) that at the baryonic component distorts the dark halo substantially, and (ii) that the variation of the initial density profiles with the anisotropy parameter is just what is required to ensure that the three profiles coincide after the baryons have been added: as Sellwood & McGaugh (2005) showed, the tangentially biased halo is more strongly compressed than the radially-biased halo, so before the barons are added its density needs to be lower than that of the radially-biased halo.
The dotted curves in the upper panel of Fig. 4 show the profiles one obtains when one starts with the three isolated halos and applies to these objects the classic adiabatic contraction formalism of Blumenthal et al. (1986). This formalism is based on the fiction that all dark matter particles are on circular orbits, i.e. an extreme case of tangential bias. From the lower panel of Fig. 4 we see that the analytic prediction is in excellent agreement with our results for , but under-estimates the contraction for smaller radii.
We can for the first time study how the flattened baryon distribution breaks the spherical symmetry of the halo, and pulls it towards the plane – hitherto analytic prescriptions for adiabatic distortion have been restricted to the spherical case. Fig. 5 shows the axis ratio of the final dark halo as a function of elliptical radius. We see that the radially biased halo is most strongly flattened by the baryons – its axis ratio does not rise above until near . The axis ratio of the initially nearly isotropic halo rises more gradually from similar values at to near the Sun. We emphasise that these axis ratios are arising naturally and are not related to our choice of a reference model with a similarly flattened halo. In fact, we found very similar axis ratios in test runs using the P14 model with a spherical halo and then inspired by these results chose to use the flattened model as our reference.
Fig. 6 shows the fractional change in the dark halo’s density that is induced by adiabatic insertion of the baryons. We see that for any velocity anisotropy the effect is strongly concentrated to the inner Galaxy, . In the case of radial velocity anisotropy, the dark halo’s density is most enhanced in the region , and the enhancement could be (mis-)interpreted as the formation of a dark disc around the baryonic disc.
4.3 Rotation curve
Our P14 reference model was fitted to data constraining the Galaxy’s rotation curve using an NFW halo profile. As we have seen above, the adiabatic contraction of the dark halo significantly alters the halo’s central density profile. To counteract these changes a far as possible, we re-scaled the halo’s mass such that we still obtained the right circular speed in the solar cylinder. However, because contraction steepens the halo’s inner density profile, the halo’s contribution to rises inwards faster than in the P14 model, and when the contracted halo is used in conjunction with the P14 disc and bulge the model no longer matches the constraints on our Galaxy’s inner rotation curve.
Comparison of the black and full green curves in Fig. 7 shows the extent of the mismatch. In the left panel, the model circular-speed curve remains flat down to much lower radii than the curve derived from the observations. The right panel shows the extent of the conflict with observations of the terminal velocities of HI gas (Malhotra, 1995). This is a fundamental problem that ties back to the use of an NFW profile in our reference model. This profile fits baryon-free dark halos and will not fit dark halos to which baryons have been adiabatically added. If the addition of baryons is non-adiabatic, the NFW profile might fit real dark halos, but we have no reason to suppose that the addition of baryons was non-adiabatic to the precise extent required for the NFW profile to remain valid.
4.4 DM velocity distributions
We now examine the velocity distribution of the dark matter particles within the Galactic plane of our models. Fig. 2 shows the radial dependence of the classical anisotropy parameter (equation 9) before and after adiabatic contraction. The most significant changes are in the innermost regions where the baryonic discs are dominant. In the radially biased case (green curves), the anisotropy increases slightly at and scarcely changes further out. In the other two cases the anisotropy increases at all radii, but no qualitative changes are observable outside the disc region.
Fig. 8 shows the distributions of the three velocity components (after marginalising over the other two components) of dm particles at the location of the Sun.333We adopt . The distributions at the actual solar position 14 to above the Galactic plane are indistinguishable from these. The non-Gaussian nature of these velocity distributions is evident. In all cases adiabatic contraction of the halo broadens the distributions for the same reason gas that is adiabatically compressed heats up. As the -direction broadens, its top becomes remarkably flat. The dispersions of the distributions are listed in Table 2.
The sensitivity of particle detectors typically depends on the speed of the particles to be detected (e.h. Green, 2012). The full curves in Fig. 9 show for each final model the distribution of the speeds of dm particles with respect to the Sun. The tangentially-biased model has the most slow-moving and least fast-moving particles, while the radially biased model has more particles with . Hence dm particles will be more readily detected in the radially-biased case, which is fortunately the case to which N-body simulations point (e.g. Hansen & Moore, 2006; Ludlow et al., 2011; Pontzen & Governato, 2013).
The dashed curves in Fig. 9 shows the Maxwellian distributions with the same velocity dispersions as the real distributions. These Maxwellian equivalents are quite similar to one another and differ materially from the actual distributions. In any stellar system the actual velocity distribution has to fall below a Maxwellian at large because it has to vanish above the escape speed, and this effect is evident in Fig. 9.
A model distribution of speeds that falls to zero at a finite speed is that proposed by Tsallis (1988):
which tends to a Maxwellian as . Ling et al. (2010) showed that the Tsallis distribution fitted the speeds of particles in a spherical shell of radius in their cosmological hydrodynamical simulations of baryon infall. The full curves in Fig. 10 are identical to those in Fig. 9, while the dashed curves show the best-fitting Tsallis distributions. Table 3 gives the parameters of these fits. The fit to the speed distribution of the isotropic (orange) model is excellent – in fact, the curve for the Tsallis distribution is frequently invisible under the true distribution. The fit to the tangentially biased model is very good, but the Tsallis distribution completely fails to reproduce the distribution of the radially-biased model.
4.5 The stellar disc
Figs. 11 and 12 show the stellar density distribution in the case of the radially biased halo by showing the disc’s density at fixed as a function of and at as a function of , respectively. The black curves show the result of integrating the disc df in the potential of the reference model, while the green curves are obtained by integrating the same df in the model with the contracted halo.. The differences between these profiles are significant only in the outer disc (Fig. 11) where the final halo is rounder than that of the reference halo. Fig. 12 shows that in the solar cylinder the two density profiles agree extremely well, so the final model still provides an excellent match to the star count data of Jurić et al. (2008) to which the reference model was fitted.
5 Relation to other work
The growth of baryonic discs in dark halos has been extensively studied with particle simulations (e.g. Pillepich et al., 2014, and references therein). In the great majority of studies, the manner in which baryons have infiltrated the dark halo has been controlled by “sub-grid physics”, that is, prescriptions governing gas cooling, star formation, and the consequent heating and ejection of gas. These prescriptions are only loosely based in physical principles and are freely adjusted to bring certain statistical properties of the baryonic structures that form into agreement with observations. In these simulations dark halos do not necessarily evolve adiabatically – in fact in most simulations it is likely that scattering of dm particles by bars and massive gas clouds will have significantly lowered the peak phase-space density of the dark halo.
Debattista et al. (2008) adiabatically introduced rigid discs into triaxial dark halos that had previously been formed by colliding spherical dark halos. After their adiabatic introduction, the discs were adiabatically removed by evaporation, and the final distributions of dark-matter particles were compared with the initial distribution. The principal conclusion of this study was that the initial and final distributions were very similar, which indicates that the dark halos responded adiabatically. In the absence of a disc, the halos were mostly strongly prolate, and their outer regions remained strongly prolate at all times, but the discs distorted the inner regions of their halos to near axisymmetric oblate forms. Thus Debattista et al. (2008) underlines the value of being able to compute the response of a dark halo to the adiabatic insertion of axisymmetric baryonic structures.
Another study in which for a period the dark halo must have evolved adiabatically is that of DeBuhr et al. (2012), for between redshift and they adiabatically grew a massive disc. At the disc became a fully fledged N-body disc capable of developing a bar and strict adiabaticity ended. The paper focuses on the axis ratios of the halo and the evolution of the disc. It provides no quantitative information regarding the evolution of the phase-space density of dm particles.
Pillepich et al. (2014) compared the final distribution of dark matter in two simulations of the formation of an object that is slightly less massive than our Galaxy, in one of which baryons were included while the other contained only dark matter. The radial DM profiles for the models with and without baryons in their Figure 1 are very similar to our plots of the DM density before and after baryon addition in the upper panel of Fig. 4. They found that the inclusion of baryons increased dm density at the Sun by per cent, in excellent agreement with our findings. Two processes contribute to enhancement of the dm density in the plane of the simulation with baryons. One is the pinching of the dark halo by the disc’s gravitational field, and the other is the capture of satellites onto highly inclined orbits through dynamical friction on the baryonic disc (Quinn & Goodman, 1986) and their subsequent tidal shredding into a rotating dark disc (Abadi et al., 2003; Read et al., 2008, 2009; Ruchti et al., 2014; Read, 2014). Our models include the first mechanism but not the second. The second process forms a dark disc that rotates at a significant speed in the same sense as the baryonic disc, while pinching of the dark halo does not set the halo rotating. Pillepich et al. (2014) showed that for the rather quiescent merger history expected for our Galaxy, the effect of pinching predominates, accounting for of the overall dark disc.
Ling et al. (2010) and Pillepich et al. (2014) used hydrodynamical simulations of baryon insertion to examine the velocity distribution of dm particles at the Sun’s location. The simulation of Pillepich et al. (2014) involved 13 million DM and 13 million baryonic particles in the high-resolution region, while the simulations of Ling et al. (2010) employ slightly fewer particles. Pillepich et al. (2014) inferred the local velocity distribution of DM particles from the DM particles in the region , , while Ling et al. (2010) studied the particles in the smaller volume , . Despite the enormous computational resources deployed to run these simulations, in neither study is Poisson noise unimportant, and Pillepich et al. (2014) considered it necessary to smooth their velocity distribution to a resolution. Our models require a tiny fraction of the computational resource and provide velocity distributions that relate precisely to the Sun’s location and yet are completely free of Poisson noise. Neither Ling et al. (2010) nor Pillepich et al. (2014) were able to study the dependence of halo kinematics on the extent of radial bias in the initial dark halo.
The speed distributions at the Sun reported by Ling et al. (2010) and Pillepich et al. (2014) are not unlike the one we obtain in the initially isotropic case. A significant difference between their results and ours is that they find significant net rotation of the local dm particles whereas we have none. For example Ling et al. (2010) find . Streaming in the dark halo is a consequence of non-adiabatic interaction between the baryons and the dark halo. In addition to the shredding of dark halos discussed above, a bar will set a dark halo spinning by losing angular momentum to it (Tremaine & Weinberg, 1984).
6 Discussion & Conclusions
We have shown how to construct fully self-consistent multi-component galaxy models with axisymmetric action-based distribution functions . As a worked example we have taken a Galaxy model presented by P14 and made its dark halo a fully dynamical object. In particular, we determined the gravitational potential that is self-consistently generated by the stellar disc and the dark halo in the presence of pre-defined contributions from a gas disc and an axisymmetric bulge.
To make the dark halo a dynamical object, we must select its velocity anisotropy. To this end we have introduced a new family of dfs. We have used the new dfs to explore three options: radially biased, isotropic and tangentially biased dark halos. The first two models bracket the predictions of cosmological dmonly simulations.
The most important change in the final model is the adiabatic compression of the dark halo by the flattened potential of the discs; we can study this process at a tiny fraction of the computational cost of a cosmological simulation. The contraction has two components. The first component is a spherical shrinkage. At radii , i.e., roughly outside one disc scale height, this is quite well approximated by the simple formalism of Blumenthal et al. (1986). At smaller radii we find stronger contraction. This is the reverse of what is usually found with N-body simulations, in which generally less contraction is observed than Blumenthal et al. predict (Abadi et al., 2010; Gnedin et al., 2010).
The second component of compression is a pinching towards the plane at radii , which flattens the halo. The extent to which the halo flattens on introduction of the discs depends of its velocity anisotropy in the sense that radial bias maximises the flattening. Our initially spherical, radially biased model deforms to an axis ratio out to . The tangentially biased halo, by contrast, has at . This has important consequences, since cosmological simulations predict radially biased configurations. Indeed, the left panel of Fig. 6 shows that in the case of a radially biased dark halo, the difference between the actual, compressed halo and the original spherical halo looks very much like a “dark disc”. A “disc” of this type would be non-rotating, in contrast to a dark disc formed by the capture of satellite halos onto low-inclination orbits followed by tidal shredding (Quinn & Goodman, 1986; Abadi et al., 2003; Read et al., 2008; Ruchti et al., 2014; Read, 2014).
Our models yield the detailed velocity distribution of dm particles at the Sun, upon which depends the detectability of dark matter in direct detection experiments. The distributions of speed with respect to the Sun are distinctly non-Gaussian, being more and less sharply peaked in the radially and tangentially biased models, respectively. Given the tendency of particle detectors to have a threshold speed for detection, the radially biased model offers significantly better chances of detecting dark matter.
Recently, P14 modelled the Galaxy by combining dynamical models of the stellar disc with star counts and the kinematics of RAVE stars to obtain estimates the local dm density. They found that the major model uncertainty was the flattening of the dark halo. The best-fitting local dm density increased with the assumed halo axis ratio as (cf. their Fig. 9). Here we find that even if the dark halo were originally spherical, the dark halo must be flattened inside the solar radius , with to 0.9. Even though varies with radius, this variation is modest within , the relevant region from the perspective of P14.
Using a very different approach to that of P14, Bienaymé et al. (2014) determined the local dm density from a subset of the data used by P14. The local dm densities of these two studies agree for halo axis ratio in the range 0.79 to 0.94. It is interesting that this is just the range of axis ratios the present study leads us to expect if the dark halo were spherical before the disc was added.
The dfs of our dark halos were chosen so they generate isolated dark halos that closely follow the NFW profile. When cohabiting with the P14 disc, the halos becomes more centrally concentrated, with the consequence that the final circular-speed curves do not match the data as well as in the reference model, where the dark halo has an NFW profile in the presence of the discs and bulge. The models do, however, have appropriate values for the circular speed at the Sun.
There are two distinct ways in which we could modify the present model so as to restore a satisfactory fit to the circular-speed curve while keeping the dark halo alive. One way is to compensate for the increased central concentration of the halo by increasing the disc scale lengths from the values adopted by P14. The other is to modify the df of the dark halo in such a way that it approximates the NFW profile when its density is evaluated in the presence of the baryons rather than when isolated. The case for modification of the df is that the introduction of the baryons was not completely adiabatic, with the consequence that decreased near the origin of action space and increased away from the origin to leave the total halo mass constant. The modification of the df could be used to introduce a degree of rotational streaming such as will arise from friction against the bar and entrapment of sub-halos on highly inclined orbits. Cosmological simulations could be used to guide the choice of halo df.
While there is a strong case for modifying the halo df from the very simple form adopted here, it is not clear that the modified df should yield an NFW profile in the presence of the baryons. Therefore, the next step should be to see whether a plausible disc df exists which when combined with the present halo df will produce a satisfactory model of our Galaxy.
We hope soon to report new constraints on the structure of our Galaxy’s dark halo obtained by fitting to observational data model Galaxies in which the stellar disc, stellar halo and dark halo are all represented by dfs of the form . The time is also ripe to fit such models to data for external galaxies, which a new generation of integral-field spectrographs are making extremely rich. We anticipate that these models will constrain the dm distributions strongly because the stars will be described by either distinct dfs for each observationally distinguishable range of metallicities, or by an “extended distribution function” (Sanders & Binney, 2015b).
We thank the member of the Oxford Galactic Dynamics group for helpful comments and suggestions.
JB was supported by STFC by grants R22138/GA001 and ST/K00106X/1. 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.
- Abadi et al. (2003) Abadi M. G., Navarro J. F., Steinmetz M., Eke V. R., 2003, ApJ, 597, 21
- Abadi et al. (2010) Abadi M. G., Navarro J. F., Fardal M., Babul A., Steinmetz M., 2010, MNRAS, 407, 435
- Aumer & Binney (2009) Aumer M., Binney J. J., 2009, MNRAS, 397, 1286
- Bienaymé et al. (2014) Bienaymé O. et al., 2014, A&A, 571, A92
- Binney (2012a) Binney J., 2012a, MNRAS, 426, 1324
- Binney (2012b) Binney J., 2012b, MNRAS, 426, 1328
- Binney (2014) Binney J., 2014, MNRAS, 440, 787
- Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Princeton University Press
- Binney & Evans (2001) Binney J. J., Evans N. W., 2001, MNRAS, 327, L27
- Bissantz et al. (2004) Bissantz N., Debattista V. P., Gerhard O., 2004, ApJ, 601, L155
- Blumenthal et al. (1986) Blumenthal G. R., Faber S. M., Flores R., Primack J. R., 1986, ApJ, 301, 27
- Combes & Sanders (1981) Combes F., Sanders R. H., 1981, A&A, 96, 164
- Debattista et al. (2008) Debattista V. P., Moore B., Quinn T., Kazantzidis S., Maas R., Mayer L., Read J., Stadel J., 2008, ApJ, 681, 1076
- DeBuhr et al. (2012) DeBuhr J., Ma C. P., White S. D. M., 2012, MNRAS, 426, 983
- Gilmore & Reid (1983) Gilmore G., Reid N., 1983, MNRAS, 202, 1025
- Gnedin et al. (2010) Gnedin O. Y., Brown W. R., Geller M. J., Kenyon S. J., 2010, ApJ, 720, L108
- Green (2012) Green A. M., 2012, Modern Physics Letters A, 27, 1230004
- Hansen & Moore (2006) Hansen S. H., Moore B., 2006, New Astronomy, 11, 333
- Jurić et al. (2008) Jurić M. et al., 2008, ApJ, 673, 864
- Kordopatis et al. (2013) Kordopatis G. et al., 2013, AJ, 146, 134
- Leroy et al. (2008) Leroy A. K., Walter F., Brinks E., Bigiel F., de Blok W. J. G., Madore B., Thornley M. D., 2008, AJ, 136, 2782
- Ling et al. (2010) Ling F. S., Nezri E., Athanassoula E., Teyssier R., 2010, J. Cosmology Astropart. Phys., 2, 012
- Ludlow et al. (2011) Ludlow A. D., Navarro J. F., White S. D. M., Boylan-Kolchin M., Springel V., Jenkins A., Frenk C. S., 2011, MNRAS, 415, 3895
- Malhotra (1995) Malhotra S., 1995, ApJ, 448, 138
- Marasco et al. (2013) Marasco A., Marinacci F., Fraternali F., 2013, MNRAS, 433, 1634
- Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
- Ness et al. (2013) Ness M. et al., 2013, MNRAS, 432, 2092
- Piffl et al. (2014) Piffl T. et al., 2014, MNRAS, 445, 3133
- Pillepich et al. (2014) Pillepich A., Kuhlen M., Guedes J., Madau P., 2014, ApJ, 784, 161
- Pontzen & Governato (2013) Pontzen A., Governato F., 2013, MNRAS, 430, 121
- Posti et al. (2015) Posti L., Binney J., Nipoti C., Ciotti L., 2015, MNRAS, 447, 3060
- Quinn & Goodman (1986) Quinn P. J., Goodman J., 1986, ApJ, 309, 472
- Raha et al. (1991) Raha N., Sellwood J. A., James R. A., Kahn F. D., 1991, Nature, 352, 411
- Read (2014) Read J. I., 2014, Journal of Physics G Nuclear Physics, 41, 063101
- Read et al. (2008) Read J. I., Lake G., Agertz O., Debattista V. P., 2008, MNRAS, 389, 1041
- Read et al. (2009) Read J. I., Mayer L., Brooks A. M., Governato F., Lake G., 2009, MNRAS, 397, 44
- Ruchti et al. (2014) Ruchti G. R., Read J. I., Feltzing S., Pipino A., Bensby T., 2014, MNRAS, 444, 515
- Sanders & Binney (2015a) Sanders J. L., Binney J., 2015a, MNRAS, 447, 2479
- Sanders & Binney (2015b) Sanders J. L., Binney J., 2015b, MNRAS, 449, 3479
- Sellwood & McGaugh (2005) Sellwood J. A., McGaugh S. S., 2005, ApJ, 634, 70
- Steinmetz et al. (2006) Steinmetz M. et al., 2006, AJ, 132, 1645
- Tremaine & Weinberg (1984) Tremaine S., Weinberg M. D., 1984, MNRAS, 209, 729
- Tsallis (1988) Tsallis C., 1988, Journal of Statistical Physics, 52, 479
- van der Kruit & Searle (1981) van der Kruit P. C., Searle L., 1981, A&A, 95, 105
Appendix A: the mass model
|Gas disc||Thin disc||Thick disc|
Our mass models have five components: a gas disc, a thin disc, a thick disc, a flattened bulge and a dark halo. Since the stellar halo has negligible mass, we do not include it in the mass model. For the density laws of the disc components we have
A non-zero parameter creates a central cavity in the disc. This is used to model the gas disc while for the other two discs is set to zero. The values of the parameters are given in Table 4.
The density distributions of the dark halo and the bulge components are
The parameters of the reference model are given in Table 4.
Appendix B: the disc DF
The disc df is built up out of “quasi-isothermal” components. The df of such a component is
where and are defined to be
Here , and are, respectively, the circular, radial and vertical epicycle frequencies of the circular orbit with angular momentum , while
where is the radius of the circular orbit, determines the surface density of the disc: to a moderate approximation this surface density can be obtained by using for in equation (17) the angular momentum of the circular orbit with radius . The functions and control the radial and vertical velocity dispersions in the disc and are approximately equal to them at . Given that the scale heights of galactic discs do not vary strongly with radius (van der Kruit & Searle, 1981), these quantities must increase inwards. We adopt the dependence on
so the radial scale-lengths on which the velocity dispersions decline are . Our expectation is that .
In equation (15) the factor containing tanh serves to eliminate retrograde stars; the value of controls the radius within which significant numbers of retrograde stars are found, and should be no larger than the circular angular momentum at the half-light radius of the bulge. Provided this condition is satisfied, the results for the extended solar neighbourhood presented here are essentially independent of .
We take the df of the thick disc to be a single pseudo-isothermal. The thin disc is treated as a superposition of the cohorts of stars that have age for ages that vary from zero up to the age of the thin disc. We take the df of each such cohort to be a pseudo-isothermal with velocity-dispersion parameters and that depend on age as well as on . Specifically, from Aumer & Binney (2009) we adopt
Here is the approximate vertical velocity dispersion of local stars at age , sets velocity dispersion at birth, and is an index that determines how the velocity dispersions grow with age. We further assume that the star-formation rate in the thin disc has decreased exponentially with time, with characteristic time scale , so the thin-disc df is
where is a parameter that controls the fraction of stars that belong to the thick disc.444Note, that is the ratio of the total masses of the thick and the thin discs, while the parameter used for the mass model is the ratio of the local mass densities of the two discs. Hence the two parameters are intimately related but not the same. The values of the parameters are given in Table 5.
|Thin disc||Thick disc|
Appendix C: cusps and dimples in distributions
To implement velocity anisotropy in our halo df we did not opt for the simplest possible solution. Here we explain the reasons for this decision. In the df the two factors and (Eqs. 6, 7) control the anisotropy of the system by enhancing or suppressing the population of orbits with certain properties. These factors are constant except for very radial orbits. This variation makes not a strictly homogeneous function any more, a property that was demanded by Posti et al. (2015) when they proposed this class of dfs.
If we remove the dependency by substituting the function with unity we obtain and which is sufficient to create an anisotropic model. However, when we examine the tangential velocity distribution of such a model (Fig. 13) we find that in the radially biased model, the distribution in has a cusp at , while in the tangentially biased model it has a dimple at ; only the isotropic model has a distribution of the expected domed shape. The cusp and dimple become more pronounced after adiabatic contraction.
The cusp and the dimple arise because the energy, upon which the ergodic df of the isotropic model depends, is a quadratic function of , while is a linear function of . So when the df of the isotropic model is written as a function of the actions and then Taylor expanded in at fixed , the linear dependence of upon that arises from the dependence of on has to be cancelled by linear dependencies of and on . If these linear dependencies do not cancel, is non-zero at , and acquires a cusp or a dimple there. Below we show that this delicate cancellation is disturbed by shifting from unity.
The derivative of the distribution plotted in Fig. 8 is
We require derivatives with respect to at a fixed location in the plane, so
where we have used .
We don’t have analytic expressions for and , but the following argument shows that for . Consider a star with low , and therefore a small guiding radius, that is just able to reach , i.e. has its apocentre there. It is on an eccentric orbit with substantial . Increasing moves the guiding centre outwards and thus diminishes the amplitude of the star’s radial oscillations and lowers its value of .
To see that for , consider the case of a spherical potential. Then with in the Galactic plane, so
which is clearly negative for small . In a strongly flattened potential, the vertical motion largely decouples from the planar motion, so becomes numerically small, but it remains negative.
Inserting all these insights into equation (24) we arrive at
Hence for this balance is disturbed with the consequence that acquires a non-zero gradient at , the gradient being negative when and positive when as is evident in Fig. 13.
While cusps and dimples are not strictly unphysical, we think it unlikely that real dm distributions would exhibit them. A requirement to avoid these features amounts to a restriction on permissible forms of the function as, e.g. the one we choose for this work.