Open Access
Issue
A&A
Volume 643, November 2020
Article Number A75
Number of page(s) 15
Section Galactic structure, stellar clusters and populations
DOI https://doi.org/10.1051/0004-6361/202038535
Published online 05 November 2020

© J.-B. Salomon et al. 2020

Licence Creative Commons
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

The structure and dynamics of our Galaxy have long been studied, revealing, in the process, an ever more complex picture. A great deal can be learned about these properties from the Galaxy’s vertical structure. This knowledge can, in principle, allow us to recover the vertical restoring force above the Galactic plane and teach us about the external perturbations that have affected the disc.

The first studies of the vertical dynamics of Milky Way stars date back to the early works of Kapteyn & van Rhijn (1920), Kapteyn (1922), and Oort (1926, 1932), who noted that the motions of nearby stars perpendicular to the Galactic plane were mostly uncorrelated with horizontal motions, at least for low vertical amplitude. For a stationary stellar population, the collisionless Boltzmann equation or the Jeans equations make it possible to link variations perpendicular to the Galactic plane of velocities and density to the vertical variations of the gravitational potential. Close to 25 works were published up through the 1980s on the measurement of the Galactic potential in the solar neighbourhood. Due to the small size and inhomogeneity of the stellar samples used at the time, it was not possible to obtain an accurate measurement of this potential. The first large homogeneous sample of a few hundreds K dwarfs used to measure both their vertical density variation and kinematics allowed for one of the first accurate determinations of the vertical potential (Kuijken & Gilmore 1989). Moreover, Bienaymé et al. (1987), in comparing a Galactic model of stellar populations (Robin & Crézé 1986) and Galactic star counts, as well as the age-velocity dispersion relation of stellar populations, was also able to accurately measure the vertical variation of the potential at the Galactic position of the Sun. The first accurate measurements of the Oort’s limit, that is, the dynamical determination of the total local mass density, were obtained using HIPPARCOS observations (Crézé et al. 1998; Holmberg & Flynn 2000), which provided the first complete samples of thousands of stars with truly accurate distances within 120 pc around the Sun. The 2000s saw the progressive use of samples of stars more distant from the Galactic plane ranging from a few hundred (Siebert et al. 2003) to a few thousands of stars (Bienaymé et al. 2014); see also Read (2014) and references therein. The advent of the Gaia space mission (Gaia Collaboration 2016) with its successive data releases, has led to a few new analyses (e.g. Hagen & Helmi 2018; Widmark & Monari 2019; Widmark 2019; Buch et al. 2019; Nitschai et al. 2020). For these new samples, the main sources of error were no longer statistical limitations or measurement accuracy, but systematic biases and errors related to the approximate modelling and to the underlying assumptions. Here, we carefully re-examine the problem and, especially, how the north-south asymmetries systematically affect the results.

The Galaxy is not in a stationary state and at present, this appears to be the main limitation in accurately determining the gravity field outside the Galactic plane. Non-stationarity is directly visible, for example, by examining the differences in star counts as well as velocity fields towards the northern and southern galactic poles (Widrow et al. 2012; Yanny & Gardner 2013; Williams et al. 2013; Xu et al. 2015). These north-south asymmetries can be classified into bending (odd parity in density, even in vertical velocity) and breathing (even parity in density, odd in vertical velocity) modes (Widrow et al. 2014), and the Milky Way disc is experiencing both. Faure et al. (2014), Debattista (2014), and Monari et al. (2016) showed that breathing modes can be caused by internal disturbances, such as spiral arms. Bending modes, on the other hand, are likely caused by external perturbations (Gómez et al. 2013; Widrow et al. 2014; Laporte et al. 2018; Chequers et al. 2018), even though, to a certain extent, an internal buckling bar instability can also cause such modes (e.g., Khoperskov et al. 2019). More recently, the extremely precise measurements of the second Gaia data release (Gaia DR2; Gaia Collaboration 2018a) have made it possible to highlight strong correlations between horizontal and vertical velocities up to z = ±1 kpc (Antoja et al. 2018; Binney & Schönrich 2018; Laporte et al. 2019; Bennett & Bovy 2019).

Therefore, the equilibrium of the Galactic disc is no longer a fully credible assumption. Banik et al. (2017) and Haines et al. (2019) have investigated in simulations how departures from equilibrium could affect the determination of the surface density of the disc and of the local dark matter (DM) density. They both concluded that non-stationarity of the disc causes a systematic overestimate of both the vertical force (by a factor that can be as large as 1.5 close to the plane in under-dense regions) and local density (up to 20%) when wrongly assuming equilibrium. They concluded that these biases, especially in under-dense regions, call for the development of non-equilibrium methods to estimate the dynamical matter density.

In this work, we do not propose a non-equilibrium method for the estimate of the dark matter density but, rather, we follow-up on the above-cited works by investigating, for the first time and based on actual data, how making an equilibrium assumption separately for data in the northern and southern Galactic hemispheres affect the estimates. To do so, we carefully selected red clump stars from the Gaia DR2 catalogue and, incidentally, revisited in detail the north-south asymmetries with regard to their density and kinematics. We show, in particular, how the southern kinematics is significantly less perturbed than the northern one.

The article is organised as follows. In Sect. 2, we present the sample used in our study and describe the different cuts applied to target red clump stars from the Gaia DR2 catalogue. We also check the completeness of the sample. In Sect. 3, we specify how we partition the space to obtain both density profiles and, later, velocity dispersion gradients. We then derive the density profiles separately for the north and south and we discuss in detail their north-south asymmetries at large heights. We then analytically approximate them by way of several functions to derive scale heights and scale length. In a similar approach, velocity dispersion profiles and gradients are derived in Sect. 4. Radial and vertical forces are derived in the north and south under the assumption of equilibrium and used to estimate different surface mass density and dark matter density in Sect. 5. An alternative estimate of the dark matter density via the direct modelling of the stellar phase-space distribution is given for comparison in Sect. 6. Discussions and conclusions about the degree of north-south asymmetry of our Galactic suburbs as well as the local dark matter density are given in Sect. 7.

2. Sample selection

2.1. Selection of clump stars

The Gaia DR2 catalogue provides accurate astrometry and photometry for more than 1 billion stars. Detailed information on the data processing, validation, and catalogue content is given in Lindegren et al. (2018) for the astrometry, and Riello et al. (2018) and Evans et al. (2018) for the photometry. The global validation of the catalogue can be found in Arenou et al. (2018) and the description of the Gaia archive1 in Salgado et al. (2017).

Our stellar sample was extracted from Gaia DR2 catalogue in the aim of use stars with the best precision in the six-dimensional space (positions and velocities) to maximise the volume probed and the number of stars and, finally, to achieve the best possible completeness. In order to avoid extinction problems within the plane and contamination by dwarf stars, which are dominant at lower galactic latitudes, we first selected all Gaia DR2 star members with an absolute galactic latitude larger than 45°, which was then reduced to include only the red clump stars. The sample was further restricted to stars with a magnitude in the G-band (G) smaller than 15. The final sample is made up of 2 539 093 stars.

In this sample, the red clump was identified in a colour magnitude diagram as the over-density near the red giant branch (see Fig. 1). The proxy for the absolute Gaia magnitude in the G band was computed for individual stars, using MG = G − 5log10(1000/ϖ)+5, where ϖ is the parallax in milli-arcsec. For the completeness motivation (see Sect. 2.2), only stars with G <  15 were considered.

thumbnail Fig. 1.

Part of the colour-magnitude diagram of stars from the Gaia DR2 catalogue with |b|> 45° and G <  15. The selection of our sample of red clump stars from these data is delimited with the dashed lines.

Open with DEXTER

The vertical cuts in colour GBP − GRP are 1.107 and 1.291. The superior limit is enforced because at larger GBP − GRP, the sample has incomplete velocity information. We need radial velocities in our study to later derive the different velocity gradients in every direction as well as velocity dispersions. The horizontal cuts in absolute magnitude, MG, are 0.185 and 0.83 (see Fig. 1). These strict limits leave us with 43 589 stars. For the measurement of velocity dispersions (Sect. 4), further restrictions based on stars belonging to the Gaia radial velocity spectrometer catalogue (RVS) are set.

To compute the spatial positions of stars, we used simple distance determinations from ϖ. However, this approximation is only valid when σϖ/ϖ≲ 20% (Luri et al. 2018). We check that more than 95% of the selected sample (41 457 stars) comply with this criterion, which justifies our choice. Nevertheless, a distance derived from a parallax with 20% errors can still be biased at the level of a few percents – a caveat to consider. An alternative could have been to work directly in parallax space for the dynamical density estimates, while a joint spectroscopic analysis of the red clump stars could allow for a Bayesian refinement of the distances based on their atmospheric parameters and distance modulus estimates in combination with their parallax. Here, we take the most direct route by simply inverting the parallax, but the above caveats should be kept in mind for potential future improvements.

We then used the classical Galactic cylindrical coordinates (R,ϕ,z), where ϕ equal 0 corresponds to the line connecting the Sun and the Galactic centre and positive towards Galactic rotation. We adopted the distance of the Sun to the Galactic of R0 = 8.34 kpc (Reid et al. 2014). The vertical distance of the Sun to the Galactic plane is 14 pc according to Binney et al. (1997). After correcting the vertical distances of the sample by the height of the Sun, there are 21 065 stars remaining in the north (z >  0) and 22 524 in the south (z <  0).

2.2. Completeness

The sample has been built in order to be complete up to at least 3 kpc from the Galactic plane. At this distance, the stars at our limit magnitude, G = 15, have an absolute magnitude, MG = 2.3, which is 1.5 mag fainter than the clump stars. This ensures that all the clump stars are in our sample, assuming that Gaia is complete within our colour and apparent magnitude selection.

To check this assertion, we made a comparison with the 2 Micron All Sky Survey (2MASS) Point Source Catalogue, which is complete to J = 15.9 and Ks = 14.3 in unconfused regions (Skrutskie et al. 2006). We first used the cross-match of Gaia DR2 with 2MASS provided in the Gaia archive (Marrese et al. 2019) to locate the colour of our clump stars sample in the 2MASS photometric system. The colour of the sample mainly ranges from J − Ks = 0.55 to 0.7, with J <  13.3. Next we selected all stars in 2MASS within this colour range and tried to retrieve them in Gaia DR2. Figure 2 shows the J magnitude distribution of the stars in 2MASS, superimposed with those that are not found in Gaia DR2. It is quite clear from this plot that many stars with J >  13 were not retrieved in Gaia DR2 due to our magnitude cut at G = 15. If we further restrict our sample to the stars brighter than J = 13.3, that is, the faintest stars in our clump sample, we retrieve 98% of the stars.

thumbnail Fig. 2.

Number of stars with respect to their J magnitude. Black line shows stars selected in 2MASS with 0.55 <  J − Ks <  0.7. Grey steps show stars from the previous selection that are not found in Gaia DR2, with a magnitude cut at G = 15.

Open with DEXTER

3. Density profiles

3.1. Sub-samples

From the sample of red clump stars defined above, several sub-samples were compiled in order to make it possible to simultaneously characterise the vertical density, the scale height, and the scale length. This requires the initial volume to be described both in the radial and vertical directions. Hence, the volume of the two cones is divided in sectors (or partial hollow cylinders) in order to respect the cylindrical geometry (to the first order) of the Galaxy. As a summary, Fig. 3 represents selected stars that are colour-coded with respect to their inclusion among the different sub-samples. For simple explanations of how the conical volume is split, we focus on the northern sector centred on the solar position. We briefly note the other sub-samples later in the paper.

thumbnail Fig. 3.

Sub-samples of red clump stars selected from the initial conical selection described in Sect. 2, in Galactic Cartesian coordinates, centred at the solar radius, (x, y) plane (top panel), and (x, z) plane (bottom panel). The five sub-samples delimited in the northern hemisphere (z >  0) are colour-coded. Top panel: cylindrical grid is superimposed in black solid lines. The colours given in the legends are those taking the same notation as described in the text (Sect. 3.1). Bottom panel: distribution of all sub-samples used in this study, 10 volumes, (colour-coded in the north, blue in the south) with respect to the Galactic plane (z = 0). The dashed lines illustrate the intersection of the conical selection and the plane y = 0.

Open with DEXTER

The central sector (called N in the north and S in the south) is characterised by the Galactocentric cylindrical coordinates (R, ϕ, z). It is defined in R, in kpc, such that Rmin ≤ R ≤ Rmax with Rmin = R0 − 0.4 and Rmax = R0 + 0.4. It spans the range in ϕ in radians such that ϕmin = −0.4/R0 and ϕmax = 0.4/R0. The surface of a slice of this sector, perpendicular to the z-axis, is given by .

However, this partial hollow cylinder has to be fully contained within our selection, which is done within galactic cones. Since the cones have an angle of 45°, this involves shifting the sector vertically by the largest distance to the Sun reached in the (x, y) plane. It is reached when R = Rmax and ϕ = ϕmax. This distance is given by kpc. Since our cone summit is not in the plane but at the Sun position, an additional offset of z must be applied to ensure that the north-central sector is within the cone selection. Hence, our central sector is limited in z with zmin = z + d = 0.586 kpc. In order to keep our north and south sub-samples symmetric to the plane, this offset is also applied in the southern-central sector.

To compute the density gradient along the Galactic radius, we defined two sub-samples bordering the central one shifted toward the Galactic centre (Nc and Sc), and one shifted toward the Galactic anti-centre (Na and Sa).

The last parts of our space-cuts are those neighbouring our central sector at the same Galactic radius. We defined two sectors alongside the central sector, one in the same direction of the Galaxy rotation (Nr and Sr), and one in the opposite direction of the Galaxy rotation (No and So).

The geometrical ranges of the sectors and the resulting number of stars in each sub-sample are given in Table 1. We point out that the size of the sectors was designed to conserve a minimum of about 400 stars in order to be able to perform reliable statistics up to high Galactic altitudes at a later stage.

Table 1.

Properties of the sectors.

3.2. Vertical density profiles

To construct vertical profiles, each sub-sample is binned along the z axis. These bins always contain a fixed number of 50 stars. The number of objects was settled after several tests in order to be sensitive enough to the actual density variations but not to statistical fluctuations. A consequence of keeping a fixed star number is that the volume of the bins is not constant since the individual height of slices hi has to be adapted.

If stars are sorted in |z| increasing, from j = 1 to j = Nstars, in each bin i, hi is calculated as the absolute distance between the jth star and the j + 49th star. This provides us with an adequate assessment of density variations (see shallow solid lines in Figs. 46). The volume, Vi, of a bin is defined by the minimum and maximum radii Rmin and Rmax, the minimum and maximum angles ϕmin and ϕmax of the sector to which it belongs, and hi, such that .

thumbnail Fig. 4.

Central sectors vertical density profiles. Top panel: density profiles for the northern and southern central sectors (in stars per kpc3) along the |z| axis (in kpc). Light-coloured solid lines follow sliding bins of a constant number of 50 stars while dots are the values of the density in independent bins of fixed size (see text for details). Bottom panel: normalised density differences between the north and south density profiles as a function of |z| (see Eq. (1)), red when the north is denser, and blue otherwise. Vertical error bars illustrate the normalised uncertainties added in quadrature.

Open with DEXTER

thumbnail Fig. 5.

Same as Fig. 4 but for the sectors shifted in the counter rotation (No and So) on the left and for sectors shifted in the direction of the Galaxy rotation (Nr and Sr) on the right.

Open with DEXTER

thumbnail Fig. 6.

Same as Fig. 4 but for the sectors shifted in the direction of the Galactic centre (Nc and Sc) on the left and for sectors shifted toward the anti-centre (Na and Sa) on the right.

Open with DEXTER

In a second stage (also to get reliable Poisson-like statistics), we derive the densities again, but this time keeping a constant width for bins (denoted ⟨b⟩). This fixed width is given by the first method as the average width on a considered range. Given the large range covered by the density, we have to define four regions to derive ⟨b⟩ in order to get a number of stars (n*) per fixed bin that is neither too high nor too low (actually between ∼20 and ∼100 with a median at 50 stars): |z| ≤ 1.3 kpc; 1.3 < |z| ≤ 1.9 kpc; 1.9 ≤ |z| < 2.5 kpc; 2.5 kpc ≤ |z|. Densities derived with this second method are represented as dots in Figs. 46.

Top panels in Figs. 46 present the vertical density profiles, in stars per kpc3, up to z = 3.5 kpc in each of the sectors, that is, sub-samples in the north and south: (i) for the central sectors; (ii) for the sectors shifted in the opposite direction of the Galaxy rotation and in the direction of the rotation; and (iii) for the sectors shifted towards the Galactic centre and anti-centre. The lower limits on each of these plots are related to the geometry of the selection described in Sect. 3.1 and the upper limits are due to the poor statistics at overly large heights above the plane. Errors bars are illustrating one standard deviation of Poisson noise uncertainties, and our typical tolerance of ∼20% for the distance determinations (see Sect. 2), , and .

A rough general view of the vertical profiles in the three figures exposes an identical exponentially decreasing trend for all volumes. Sub-samples at the same Galactic radius tend to overlap. As expected, there exists a radial gradient in density, visible in Figs. 6 and 7, since volumes at a smaller radius (Nc and Sc) are denser than their counterparts at a larger radius (Na and Sa). The associated scale heights and scale length are studied in Sect. 3.3.

thumbnail Fig. 7.

Vertical density profiles binned in |z| for the northern (top panel) and southern (bottom panel) sub-samples, superimposed with the best-fit profile at the three R values (upper line: R0 − 0.8 kpc, middle line: R0, lower line: R0 + 0.8 kpc). The colour code is the same as in Fig. 3.

Open with DEXTER

Nevertheless, the density distributions in the north and south are not really identical. It has already been shown, for instance, by Bennett & Bovy (2019) that density differences are significant at z <  1 kpc (under-density of 40% in the north compared to the south at z ∼ 500 pc, which is not analysed here given our selection function). These differences are directly related to the spiral shape visible in phase space (z, Vz) as identified by Antoja et al. (2018). This clearly highlights the non-stationarity near the galactic plane.

When closely inspecting the north and south volumes, several characteristic areas can be identified with deviations from stationarity and symmetry at larger heights. To better visualise the north-south asymmetries, the bottom panels of Figs. 46 present the normalised density differences in opposite hemispheres (Widrow et al. 2012; Bennett & Bovy 2019):

(1)

Firstly, a global mild northern over-density, already hinted at by Bennett & Bovy (2019), ranging from z ∼ 0.6 to ∼1 kpc, can be seen in the top panel of Fig. 4, corresponding to the central sector. Given our selection function, this northern over-density cannot be followed in the other sectors. Furthermore, up to 3.5 kpc, every profile is perturbed in every sector, but generally more so in the northern hemisphere than in the south. Southern profiles tend to generally be smoother, with the notable exception of an interesting bump in the south at z ≃ −2.1 kpc, in the sector corresponding to the direction of Galactic rotation. Northern profiles are rather bumpy in all sectors. At large heights in the central sectors, the north ends up dominating. Interestingly, in the direction of Galactic rotation and the direction opposite to Galactic rotation, it is rather the southern density that tends to dominate at large heights, indicating a global azimuthal wobble, with a northern peak at the Solar azimuth, preceded and followed by southern peaks in the directions of rotation and counter-rotation. A similar wobble can be noted in the radial direction, with a global northern domination at large heights in the inner and central sectors, followed by a southern domination in the anti-centre sector. These indicate clearly the wave-like nature of the general north-south asymmetries, with a slightly shorter wavelength in the azimuthal direction.

The more localised wobbles and bumps in the vertical density can be classified into three categories: In the first, some bumps and wiggles are jointly occurring in the north and south, and, therefore, they are not visible on the bottom panels (representing the asymmetry). A noticeable example is a correlated bump around |z|≃1 kpc in the central (N, S) sector. One can again find such small correlated bumps at |z|≃1.8 kpc in the direction of rotation (Nr, Sr), at |z|≃1.95 kpc in (No, So), at |z|≃1.65 kpc in (Na, Sa), or at |z|≃2.05 kpc and 2.2 kpc in (Nc, Sc).

The second kind of disturbances are revealed by some bumps and wiggles which are clearly north-south antisymmetric, meaning that they represent anti-correlated oscillation of the density. This is clearly visible in the central sectors (N, S) at |z| ∼ [1.15, 1.25], [1.25, 1.35], and [2, 2.2] kpc. These anti-correlated features are creating large asymmetries in the bottom panels of the figures. In the anti-centre (Na, Sa), we find them at [1.4, 1.7], [1.7, 2], and [2.25, 2.95] kpc. In the sectors toward the Galactic centre (Nc, Sc), we find one large anti-correlated feature at [2.3, 2.9] kpc. In the sectors opposite to the Galactic rotation (No, So), we find them at [1.35, 1.55], [1.75, 2.05], and [2.2, 2.45] kpc. Finally in the direction of rotation, (Nr, Sr), we find small anti-correlated features at 1.65 kpc and 2.55 kpc. Obviously, such features are present in all sectors between 2 and 3 kpc, but vastly more so in the radial directions (centre and anti-centre), with only a tiny asymmetry in the three sectors at R = R0.

The last noticeable features are some bumps which occur only in one hemisphere. For instance, we can note in the central sectors (N, S) that the southern density has a plateau which globally dominates over the north at [1.4−1.8] kpc, while the north goes through a bump at 1.6 kpc, and then another one at 1.8 kpc to catch up with the south. These bumps have no counterpart in the south. The most striking of these bumps in only one hemisphere is the one at 2.1 kpc in Sr. A similar, less prominent, bump is present in the north at 1.7 kpc in Nc.

3.3. Scale heights and scale length

After having obtained the detailed density profiles in different sectors, we now want to describe their global variations. The goal here is to describe the global density scale height and scale length, which will be later used in Jeans equations. Given the north-south asymmetries noted previously, the two Galactic hemispheres will be treated separately. We want to insist that the purpose here is not to give physically fully reliable parameters for those scales but, rather, to derive approximate analytical expressions both vertically and radially, which will be useful for a Jeans analysis. In other words, we aim to parametrise, in the best possible way, the general trend of the density evolution of red clump stars with |z| between 0.6 and 3.5 kpc. At the same time, the radial variation will be also parametrized but only with three bins, namely, the central value of each radial sector (R0 − 0.8 kpc, R0, R0 + 0.8 kpc).

Several analytical expressions based on exponential or secant square functions are tested to simultaneously fit radial and vertical profiles (using from three up to nine free parameters, taken from multiple components with exponential and secant square profiles, with or without scale-length). We use the Markov chains Monte Carlo (MCMC) technique to adjust our data with analytical expressions. The best likelihood (ℒ) stemming from MCMC chains and the number of free parameters (k) in each analytical expression make the application of the Akaike information criterion (AIC, Akaike 1974) possible: AIC = 2k − 2ln(ℒ). The following expression using an exponential law in radius and height plus one in height has been retained:

(2)

The MCMC technique is employed with the five free parameters with flat priors: two normalisation factors αν, βν, two scale heights hza and hzb, and the scale length hR.

When maximising the likelihood, we assume the data follow Poisson statistics. Accuracies on the height and length scales are constrained by the amplitude of the density deviations from the adjusted model. Accuracy is not related to the model, which only partially adjusts the counts, but is dominated by variations between the different sub-samples. The measurement of this accuracy is important since it is the dominant source of uncertainty for measuring the vertical and horizontal force fields. A key advantage of using a MCMC method is the production of a probability density function for each output parameter, providing access to realistic uncertainties (for the posterior distribution correlation matrix, see Appendix A).

Profiles are simultaneously adjusted from |z| equal 0.6 kpc up to 3.5 kpc, and for R = [7.54,8.34,9.14] kpc for each hemisphere. Maximum likelihood parameters and their respective uncertainties are given in Table 2 for the north and the south. The vertical density profiles are shown in Fig. 7 for the north (top panel) and south (bottom panel) sub-samples. The best-fit density profile is superimposed with black lines for the three values of R.

Table 2.

Best-fitting parameters of the density variations.

As noted in Sect. 3.2, density distributions differ by Galactic hemisphere. The scale length in the north (2.162 kpc) is shorter than in the south (3.113 kpc). Two well-defined scale heights have been found, namely, 0.219 and 0.638 kpc for the south and 0.262 and 0.827 kpc for the north. The difference between these two scale heights, hza and hzb, is larger for the north (565 pc) than for the south (419 pc), which may indicate that the density transition in the north between the two vertical density components is sharper or more pronounced than in the south. Although our goal is not to make a detailed thin-thick disc decomposition, which would not be useful without additional chemical information, these two components may, nonetheless, be broadly seen as the geometrical thin and thick discs. It is, thus, reassuring that the values are in broad agreement with the literature both for the scale heights and the scale length (see e.g. Jurić et al. 2008; Robin et al. 2014). Because of the way we built the sample, we do not have information about the radial density variation below ∼1.3 kpc in |z| (see Table 1): therefore, we cannot constraint the scale length of the thin disc and the scale-length we obtained is more representative of the thick (or thicker) disc component. Extrapolating our fits towards the plane, we have a local “thick-to-thin” disc density normalisation (βν/αν) of about 4% in the north and 11% in the south.

In this section, we derive density distributions and their gradients by adjusting simple functions over extended intervals and defining the most likely values for the scale heights and scale length. It should be noted that deviations from the mean values are greater than expected from random statistical fluctuations and serve as indicators of non-stationarity and of the differences between the northern and southern sectors. We also note that due to this state of affairs, the formal errors on the parameters in Table 2 should only be taken as lower limits, since the reduced χ2 values are larger than 1 (1.12 and 1.90 in the north and south). This may also partly explain the differences in measured galactic vertical forces obtained in works prior to the publication of the Gaia DR2 catalogue. We go on to explore how velocity dispersions differ between the north and south.

4. Velocity dispersion gradients

To be able to derive gravitational forces under equilibrium assumptions from Jeans equations, velocity dispersion variations are needed. Of course, only stars with line-of-sight velocity data can be used for this. The selection will thus contain only objects with six dimensional phase space information from the Gaia RVS. About 76% of the sample that we previously used now remains, consisting of 33 005 red clump giants (15 702 in the north, 17 303 in the south). Actually, the radial velocity condition to build the new sample results in an additional cut in magnitude. This effect is clearly shown in Fig. 8, which represents the number of stars per bin of magnitude. The first sample of this study used to derive densities in grey is limited to G = 15, as imposed in Sect. 2. The second sample in green, where only stars with radial velocity information are kept, is limited to a magnitude of 14 and strongly reduced up to a magnitude of 13. It is obvious that the sample used in this step is not complete (whereas the selection built for the density analysis is complete). Nevertheless, this is not of primary importance as the assumption that our partial sample follows the same dynamics as the complete sample is reasonable. Moreover, the number of stars is still large enough to get good statistics across the studied range.

thumbnail Fig. 8.

Distribution of red clump samples designed for this study with respect to their GaiaG magnitude. The grey distribution is the main sample used to derive the density profiles (Sect. 3). Distribution in green is a selection from the main sample where only stars with radial velocity measurements are considered in order to obtain velocity dispersion profiles (Sect. 4).

Open with DEXTER

We then apply the same space partition in sectors inside the conic selection as in Sect. 3.1. Similarly to the way the vertical density was derived, data are binned along vertical direction keeping a constant number of 50 stars per bin in a first step. Then, the average width ⟨b⟩ is calculated in the four ranges in z as defined in Sect. 3.2 in order to derive the density in bins of fixed width (see dots in Figs. 9 and 10). We adopt the conventional Galactic cylindrical coordinates (R, ϕ, z) defined in Sect. 2. Velocities are converted in this frame (as VR, Vϕ, Vz) from Gaia DR2 proper motions and radial velocity (μα, μδ, Vlos) taking the same solar velocity employed in Gaia Collaboration (2018b): 240 km s−1 for the circular velocity (Reid et al. 2014) and (U, V, W) = (11.1,12.24,7.25) km s−1 for the Sun’s peculiar velocity with respect to the local standard of rest (Schönrich et al. 2010).

thumbnail Fig. 9.

Vertical squared velocity dispersion on radial direction profiles binned in |z|, for the five sub-samples in the north (upper panel) and for the five in the south (bottom panel). The colour code per sector is the same as in Fig. 3. Black lines represent the best fit following Eq. (5) and the parameters given in Table 3.

Open with DEXTER

To calculate the velocity profiles along z in each direction (vR(z), vϕ(z), vz(z)), a σ clipping method is applied. For each bin with n*, the median and the mean velocity and the corresponding velocity dispersion (σj) are calculated. Values lying beyond 4σ from the median are rejected. Then the mean, the median, and the dispersion are re-calculated. The sequence is repeated until no values are rejected. We carefully applied this process in our study to eventually eliminate strong outlier values coming from extreme velocity stars, but in the end, less than 0.5% of stars were excluded.

Velocities per bin, on each axis j, are then calculated as:

(3)

where nex is the number of stars excluded by the σ clipping. Their corresponding velocity dispersions are derived with the following equation:

(4)

Figures 9 and 10 show the distribution of the squared velocity dispersions and in R and z directions, respectively, with respect to the Galactic height. The colour code is the same as in Fig. 3. Error bars are representing the classical one σ standard deviation, on 20% of the distance along z () and on the squared dispersion ( and ).

thumbnail Fig. 10.

Vertical squared velocity dispersion on vertical direction profiles binned in |z|, for the five sub-samples in the north (upper panel), and in the south (bottom panel). The colour code per sector is the same as in Fig. 3. Black curves represent the best fit following Eq. (6) and the parameters given in Table 4, plotted at the three R values (7.54, 8.34, and 9.14 kpc from the highest to the lowest).

Open with DEXTER

Once vertical dispersion profiles are built, their variations have to be characterised. A similar MCMC machinery as described in Sect. 3.3 is used to approximate the different profiles. In the same way that we searched for the best possible simple function to fit density profiles using the Akaike information criterion values, we now test several functions to simultaneously fit the dispersion profiles of each hemisphere based on combinations of exponential and linear functions for |z| between 0.6 and 3.5 kpc. To express the squared vertical dispersion in R, we ultimately use a simple linear equation:

(5)

where aR is the slope in km2 s−2 kpc−1 and bR is the normalisation in km2 s−2. The best-fitting parameters are reported in Table 3.

Table 3.

Best-fitting parameters of the vertical profile for the radial velocity dispersion squared.

We notice there is no need to introduce a radial dependency since the likelihood is not significantly improved when taking into account the radial information. But we may have to concede that the number of bins is more limited than for the density (as the sample is 25% smaller). Consequently, a potential radial gradient would be detected only if it is obvious. On top of this, the scale length of the velocity dispersion is expected to be large and it is, therefore, not unsurprising that we do not detect it. We note, however, that incorporating the uncertainty associated with a possible radial gradient as a nuisance parameter would have led to larger error estimates.

A vertical gradient exists in , which is stronger in the north. Moreover, in this hemisphere, the dispersion around the fit is much larger than in the south with strong bumps and wiggles from one bin to another. It suggests once again that the Galaxy is not in a stationary state and is far from being north-south symmetric.

To fit the squared vertical velocity dispersion in z, the best choice turns out to be the equation:

(6)

where az and bz are normalisation factors in km2 s−2. Parameters and are respectively the scale length and the scale height of in kpc. We note that since the bz normalisation factors are typically negative, should not be considered along the lines of “classical” scale heights. Indeed, the vertical velocity dispersion increases with height. Unlike variations of , there is a need for a radial gradient to best fit, , in the different sectors. This is represented by the black lines in Fig. 10 for the three different radii considered in this study. The parameters of maximum likelihood are reported in Table 4 along with their uncertainties. Following the same trend as the radial dispersion, the gradient of the vertical dispersion is smoother for the south. It accounts for the fact the Galaxy might be locally more disturbed in the north. We note, however, that in the direction of Galactic rotation, the southern velocity dispersions display an interesting bump at large heights, both in terms of radial and vertical velocity dispersions, which is probably related to the density bump noted in the previous section.

Table 4.

Best-fitting parameters of the profile for the vertical velocity dispersion squared.

5. Vertical force and dark matter density using Jeans equations

Here, we investigate how making an equilibrium assumption separately for data in the northern and the southern Galactic hemispheres affects the estimates of the force field and gravitational potential with the Jeans equations. We mainly use the notation of Sánchez-Salcedo et al. (2016) and concentrate on the local (“central”) sector of the previous sections. The old method proposed by Oort and Kapteyn (Kapteyn & van Rhijn 1920; Oort 1932) assumed the decorrelation between vertical and horizontal motions. Here, with samples that are distant from the Galactic plane, we have to take into account horizontal and vertical couplings. We nevertheless simplify the problem as most past studies did, that is, by assuming an axisymmetric Galaxy and considering that the inclination of the ellipsoid is spherically aligned. The goal is to explore the different dynamical and dark matter density estimates that one can get in the north and south when ignoring disequilibrium.

5.1. Vertical forces

We define two global scale parameters. They can be directly calculated from the density profile derived in Sect. 3 (see Eq. (2)):

(7)

(8)

Vertical forces Kz are calculated with the following equation:

(9)

The velocity dispersion covariance term is given by:

(10)

Based on Eqs. (5), (6), and (10), we get:

(11)

To latter derive the vertical forces, we independently draw a solution in the last 20% of the three Markov chains obtained from the fits of ν(R, z), and . It lets us with a combination of 11 parameters (αν, hza, βν, hzb, hR, aR, bR, az, , bz, ) reflecting a possible parametrization of the data.

5.2. Rotation curve corrections

A correction term describing the radial forces has to be taken into account to accurately derive the surface mass density at large distance from the plane from the vertical force.

We calculate this radial force correction term from the rotation curve as:

(12)

such that the dynamical surface mass density (Σdyn = Σ* + ΣDM), equal to the stellar plus the dark matter mass surface density, is

(13)

Here, Vc(R, z) is the circular velocity related to the radial force KR, through . In the Galactic plane, this velocity has been set up to Vc(R0, 0) = 240 ± 8 km s−1 (see Sect. 4). In order to encompass the majority of recent values from the literature (for example 229.0 ± 0.2 km s−1 from Eilers et al. 20192) we extend the standard deviation to 1.5σ, meaning Vc(R0, 0) = 240 ± 12 km s−1.

As a consequence, we have to do the same for the solar radius by taking R = 8.34 ± 0.24 kpc which allows us to include in the uncertainty values from GRAVITY Collaboration (2019).

The variation of Vc with |z| is constrained by the Besançon Galactic Model (BGM) potential (Bienaymé et al. 1987). We enforce the circular velocity to vary at the same rate as it evolves in the BGM. It approximately corresponds at the solar radius to a linear decrease such that Vc(R0, z) ≈ Vc(R0, 0) − 9.5 × z(kpc).

For the radial gradient of the circular velocity, we chose to take an average value from the literature in the plane at the solar radius (Kawata et al. 2019; Mróz et al. 2019; Eilers et al. 2019), which is constant with radius and height such that ∂Vc(R, z)/∂R = ∂Vc(R0, 0)/∂R = −2.18 ± 1.82 km s−1 kpc−1.

Following the dispersions of the quantities described above, we drew a set of three parameters (R0, Vc(R0, 0), ∂Vc(R0, 0)/∂R). With the addition of the 11 parameters obtained earlier (see Sect. 5.1), it is possible to derive both vertical (Eq. (9)) and radial force correction (Eq. (12)) with respect to galactic radius and height. We then obtain the dynamical surface mass density Σdyn = Σ* + ΣDM, from Eq. (13).

By reiterating this drawing sequence, dispersions of the forces calculation and mass surface density are natural outcomes. Figure 11 is depicting the evolution of the Σdyn with respect to the height at the Solar radius for the north and south hemispheres. The width of the two curves represents the one σ standard deviation. It stems from 10 000 draws in the Markov chains for the adjusted parameters and in the standard deviation for fixed quantities.

thumbnail Fig. 11.

Surface mass density (Σdyn) in M pc−2 at the Solar position (Eq. (13)) along the direction perpendicular to the Galactic plane (|z| in kpc) for the north (red) and the south (blue). There are 10 000 superimposed curves per hemisphere, which arise from the same number of draws in parameter uncertainties. The shallow red and blue shapes encompass one σ standard deviation of solutions. The two solid lines in red and blue (approximately in the middle of each shape) represent median curves.

Open with DEXTER

At first glance, the north and south curves appear to be following the same trend since the two values for Σdyn are increasing with z on a similar order of magnitude. Nevertheless, the derived southern mass surface density is clearly much more regular, expressively reflecting the fact that disequilibrium are stronger in the north. Interestingly, we can immediately note that this southern surface density is smaller at low z, which, despite neglecting disequilibrium, does seem to reflect the slight stellar over-density in the north noted in Sect. 3.2. Then the clear slope break around 1.2 kpc for both north and south hemispheres indicates the transition between an environment largely dominated by baryons to an environment with a lower baryonic density in the thick disc, where dark matter becomes progressively dominant. Between ∼1.2 and 2.2 kpc, Σdyn presents a kind of oscillation for the north which, despite averaging trends with simple functions, still reflects the strong northern disequilibrium noted in previous sections. Bearing in mind that unresolved systematic errors discussed in the previous sections would inevitably increase the error bars, we can nevertheless note that with our formal errors, the two deduced dynamical surface densities do not overlap at one σ between 1 and 1.5 kpc.

Haines et al. (2019) investigated, using simulations, how departures from equilibrium could affect the determination of the surface density of the disc. Their conclusion was that the derived surface densities from a perturbed disc can give similar orders of magnitude in the north and south separately while still overestimating it by a factor that can be as large as 1.5 close to the plane. Therefore those surface densities should be considered as upper limits.

5.3. Dark matter density

We assume here the dark matter content is dominating the baryonic part of the Galaxy at a large height from the plane, above two kpc. Hence, it gives dΣdyn/dz = dΣDM/dz for |z|> 2 kpc. Then it allows us to derive the dark matter mass density at these heights:

(14)

The density so determined is only slightly smaller (few percents at 2 kpc) than the one in the plane (ρDM(R0)) under the assumption of a spherical halo. Obviously, for a flattened oblate DM halo, the local dark matter density in the plane should then be even larger. Hence, the value derived here is a lower limit for the plane dark matter density (e.g. Piffl et al. 2014).

Above 2 kpc, the variation of Σdyn is well-approximated with a linear variation, especially in the south (see Fig. 11). Thus, for each of the 10 000 solutions, we derive the average gradient of Σdyn for each hemisphere. It leads us to obtain the dark matter mass densities with their associate uncertainties:

(15)

(16)

The choice of 2 kpc as the lower limit is somewhat arbitrary, so it is interesting to check what we obtain for 1.5 kpc and 2.5 kpc, respectively. Interestingly, in the south, the density is fully unaffected, 0.0098 ± 0.0015 M pc−3 and 0.0098 ± 0.0043 M pc−3, respectively. In the northern hemisphere, however, we get different values: above 1.5 kpc, we have 0.0099 ± 0.0017 M pc−3 and above 2.5 kpc, we have 0.0161 ± 0.0044 M pc−3.

Based on our finding that the northern hemisphere is much more perturbed than the southern one, the wobbling dynamical surface density yields dark matter densities that are much more varying, and a priori less reliable than in the south. With these caveats in mind, we can also proceed to carry out a direct modelling of the stellar phase-space distribution function, as a matter of comparing it with the results obtained from the Jeans analysis.

6. An alternative modelling of the vertical force and dark matter density

With the caveats of the previous sections in mind, it is instructive to proceed, as a means of comparing, to a global classical determination of Kz and of the local dark matter density by modelling stellar density distributions and velocity dispersions using isothermal components, and slightly redefining the sample selection geometry. This allows us to check the robustness of the previous results under equilibrium assumptions, meaning that at least the orders of magnitudes of the obtained values are similar, and to quantify more precisely the role of the local value of the circular velocity and its gradient in the systematic uncertainties.

We adopted, in particular, the method developed by Bienaymé et al. (2014), which uses pseudo-isothermal stellar phase-space distribution function components dependent on three integrals of motion. This extends the validity of the classical method of Oort (1932) well beyond 1 kpc above the Galactic plane.

For the star counts, we adopt the fits obtained for the “central” sample of stars towards the Galactic poles and, thus, we simply use the observed density distributions shown in Fig. 7 (central north and south) reproduced in Fig. 12. For the vertical velocity dispersion distributions towards the poles, instead of a cylindrical sector volume, we use a conical volume directed towards the Galactic poles with a 20 × 20° aperture (extending to 2.8 kpc). We compute the velocity dispersions in 200 pc z-wide bins (see Fig. 13). We note that in the study of the Kz of Bienaymé et al. (2014) as well as in that of Zhang et al. (2013), the stellar samples are subdivided into several sub-samples of different metallicities in order to reinforce the constraints on the Kz measurement. As metallicities are not available for all stars in our study, we did not subdivide the samples. The velocity dispersions σz are determined per 200 pc height interval with σ clipping.

thumbnail Fig. 12.

Vertical density profiles of red clump stars towards the north (left) and the south (right) Galactic poles. Star counts and (continuous lines) best-fit model.

Open with DEXTER

6.1. Fitting a self-consistent model

For the determination of the vertical force Kz, we use the Galactic mass distribution model of Bienaymé et al. (2014). This model incorporates components that are Stäckel potentials (e.g. Famaey & Dejonghe 2003) including two Kuzmin-Kutuzov spheroids: one of them, more flattened, represents the galactic disc whose thickness is set at 300 pc, which is typical for the old galactic stellar disc. The other component is the dark halo. The adjusted parameters are the masses of these two components and the flattening of the dark halo. A third component with a highly concentrated mass distribution adds an additional degree of freedom to the baryonic mass distribution of the disc and thus allows us to model the contribution of the disc to the Kz in a more general and realistic way.

These three-components models allow to represent the gravitational potential of the Galaxy (Famaey & Dejonghe 2003). The density and velocity dispersions of our northern and southern samples are modelled using distribution functions that are generalisations of the Shu (1969) distributions (see Bienaymé et al. 2015). They are defined using three integrals of motion E, Lz, and I3, which are precise with regard to Stäckel potentials. The associated distribution functions allow us to then adjust the gravitational potential with an excellent level of accuracy (see Bienaymé et al. 2014 for a detailed description).

We perform a least-squares fit, separately in the north and south directions, of the vertical density and velocity dispersion distributions (see Figs. 1213). The model parameters are three masses and two thicknesses for the Galaxy’s gravitational potential, as well as (for each of the three distribution functions) the local stellar density in the Galactic plane at the position of the Sun and the local vertical velocity dispersion.

thumbnail Fig. 13.

Vertical velocity dispersion of red clump stars towards the north (left) and the south (right) Galactic poles. Measured dispersions and the (continuous line) best-fit model.

Open with DEXTER

The best fits show significantly different mass distributions between north and south that can reproduce the observed ν(z) and σz(z) in the distance intervals considered (400 pc to 3 kpc). Distinct solutions are found and the parameters obtained for the mass distribution of the different components are highly correlated. The dark matter density of the dark halo is not constrained and a wide range of values of ρDM seems to be compatible with the observed spatial and kinematic distributions of the red clump giants. However these different models differ from each other by their circular rotation curves, producing different values of the circular rotation speed Vc(R0) at the Solar Galactic radius as well as by different gradients of the rotation curve at the position of the Sun ∂Vc(R0, 0)/∂R. These models all have an almost spherical dark halo.

6.2. Adding rotation curve constraints

To distinguish among our different models that equally fit star counts and velocity dispersion, one can add two supplementary constraints: the circular rotation speed at R0 and the slope of the rotation curve at R0.

For instance, fixing those values to Vc(R0) = 233 km s−1 (McGaugh 2018) and ∂Vc(R0, 0)/∂R = −1.7 km s−1 kpc−1 (Eilers et al. 2019), we find ρDM = 0.0108 ± 0.004 M pc−3 for the combined sample. Figure 14 (left) shows the vertical force for the north and south taken together, and also the contributions of the dark halo component and the baryonic disc to the Kz, while Fig. 14 (right) shows the corresponding integrated surface density Σ(z). The formal error on the estimate of ρDM is on the order of a few percent. We note that the reduced χ2 of the best fit is, again, larger than 1 (namely, 1.7). However, we note that with this method, the recovered value of the local dark matter density is identical in the north and south. The formal error on this value is much smaller than the systematic error coming from the different possible values of the two local rotation curve parameters.

thumbnail Fig. 14.

Left: vertical force (in M pc−2 units) perpendicular to the galactic plane versus the vertical distance. Red shows the total vertical force. Blue shows the disc contribution to the vertical force. Black shows the dark matter contribution to the vertical force. Right: surface density in M pc−2 integrated from −z to z. Red shows the total surface density. Blue shows the disc contribution. Black shows the dark matter contribution to the integrated surface density. The differences between the left and right panels can be understood from Eqs. (12) and (13) for each component. The rotation curve of the disc component is decreasing and that of the halo is increasing, leaving the total rotation curve almost flat (albeit slightly decreasing). Therefore, the surface density of the disc is smaller than its corresponding vertical force (blue curves) and the opposite is true for the halo (black curves). The total surface density is slightly smaller than the corresponding vertical force up to z ∼ 2 kpc (red curves), but the trend is then reversed due to the details of the mass distribution at larger heights.

Open with DEXTER

A systematic correction can be applied to the local dark matter density based on the adopted value of these two parameters:

This means that, in adopting the same values as in the rest of the paper, namely, Vc(R0) = 240 km s−1 and ∂Vc(R0, 0)/∂R = −2.18 km s−1 kpc−1, we find ρDM(R0)M pc−3 = 0.0144, which is slightly larger than the results from the Jeans analysis.

The formalism used in this section would also benefit in the future from a full MCMC treatment to further explore the correlation between rotation curve parameters and associated uncertainties and the local dark matter density. Including the stellar samples away from the Sun in a single simultaneous model fit could also prove extremely interesting, bearing in mind the caveats about non-stationarity. This will be the subject of future works.

7. Discussion and conclusion

In this paper, we analyse in detail the kinematics and dynamics above the Galactic plane of a homogeneous and 98% complete sample of 43 589 red clump giants with G ≤ 13.3 from Gaia DR2. The spatial distribution of these stars is limited to the Galactic poles, in a cone of 45°. This was done for completeness purposes, which were validated against 2MASS catalogue up to at least |z| = 3.5 kpc.

We have sliced our sample in ten sectors, five in the northern hemisphere and five in the south, separated by 0.8 kpc, covering both the inner disc and outer disc regions, as well as the regions in the direction of Galactic rotation and counter-rotation. We then analyse in detail the north-south asymmetries reflecting the non-stationarity of the Galactic disc. In particular, a global azimuthal wobble has been noted with a northern density peak at the Solar azimuth, preceded and followed by southern density peaks in the directions of rotation and counter-rotation. A similar wobble was found in the radial direction, with a global northern domination at large heights in the inner and central sectors, followed by a southern domination in the anti-centre sector. These indicate clearly the wave-like nature of the general north-south asymmetries, with a slightly shorter wavelength in the azimuthal direction. The more localised wobbles and bumps in the vertical density can be classified into three categories, with: (i) some bumps and wiggles jointly happening in the north and south; (ii) some being clearly north-south antisymmetric, meaning that they represent anti-correlated oscillation of the density, present in all sectors between 2 and 3 kpc, but vastly more so in the radial directions (centre and anti-centre); and (iii) some bumps happening only in one hemisphere, such as a striking southern bump at z = −2.1 kpc in the direction of rotation. Notwithstanding, one of our main conclusions is that in our Galactic suburbs, the northern hemisphere is generally more perturbed than the southern one. It will be interesting to understand how this can inform our knowledge of the interaction of the Galactic disc with external perturbers, such as the Sagittarius dwarf (e.g. Laporte et al. 2019).

Averaging over the bumps and wiggles, we derive the scale-length and a double scale-height of our sample separately in the north and south. Although, in the absence of chemical information, our goal was not to make a detailed thin-thick disc decomposition, these two scale-heights roughly correspond to the contribution of the geometrical thin and thick discs. We note that the “thick disc” scale height is significantly larger in the north (∼827 pc) than in the south (∼638 pc). Due to our selection function, the radial scale-length is essentially constrained by the thick disc and it is significantly shorter in the north (∼2.2 kpc) than in the south (∼3.1 kpc). It will be particularly interesting to conduct in the future a detailed chemical analysis of our sample of red clump giants (e.g. with WEAVE and 4MOST) to better understand these spatial structures of the thick disc in the northern and southern Galactic hemispheres.

Restricting the sample to Gaia RVS stars, which leaves us with 33 005 stars, we analyse the velocity dispersion gradients in the northern and southern hemispheres and, again, we find a much more perturbed northern Galactic hemisphere at heights above z ∼ 1 kpc. However, we also noted an interesting bump at large heights in the southern velocity dispersions, restricted to the sector in the direction of Galactic rotation.

After this analysis of the north-south asymmetries, we investigated how making an equilibrium assumption separately for data in the northern and southern Galactic hemispheres affects the estimates of the dynamical matter density and dark matter density. We note that many systematics could enlarge our error estimates, meaning that the formal errors quoted in this study are to be taken as lower limits. For this task, we concentrated on the local sample (“central” sector) with most of the data. Given the findings above, this assumption of equilibrium is likely to hold better in the south than in the north. We find that a Jeans analysis delivers relatively similar vertical forces and integrated dynamical surface densities at large heights above the plane in both hemispheres. At these heights, the densities of stars and gas are very low and the surface density is largely dominated by dark matter, which allows us to estimate, separately in the north and in the south, the local dark mater density derived under equilibrium assumptions. We find a local dark matter density (estimated above 2 kpc) ρDM ∼ 0.013 M pc−3 (∼0.509 GeV cm−3) in the perturbed northern hemisphere and ρDM ∼ 0.010 M pc−3 (∼0.374 GeV cm−3) in the much less perturbed south.

Finally, with all these caveats in mind, we obtained a value of ρDM ∼ 0.011−0.014 M pc−3 (depending on the local value of the circular velocity and its gradient) through a global phase-space distribution function fit to the data under global equilibrium assumptions. While in global agreement, we note that with the exact same rotation curve constraints, this distribution function analysis yields a slightly larger value of the local dark matter density (0.014 M pc−3) than the Jeans analysis. We conclude that at this stage, and once the rotation curve is fully fixed, only future developments of non-equilibrium methods would allow for a full evaluation of potential systematics due to non-stationary effects and, thus, allowing for more solid estimates of the dynamical matter density in the Galactic disc to be obtained.


2

The main recent literature values are 236 ± 3.8 km s−1 from Mróz et al. (2019) based on the analysis of Cepheids, that of McGaugh (2018), 233.3 ± 1.4 km s−1 based on a new analysis of the gas rotation and the determination of the distance to the Galactic centre (GRAVITY Collaboration 2018) or the analysis of Kawata et al. (2019), 236 ± 3 km s−1, and McMillan (2017), 232.8 ± 3.0 km s−1.

Acknowledgments

The authors would like to thank the International Space Science Institute, Bern, Switzerland, for providing financial support and meeting facilities. This work was supported by the Programme National Cosmology et Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, co-funded by CEA and CNES. JBS acknowledges financial support of the CNES. BF acknowledges funding from the Agence Nationale de la Recherche (ANR project ANR-18-CE31-0006 and ANR-19-CE31-0017) and from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 834148). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

References

  1. Akaike, H. 1974, IEEE Trans. Autom. Control, 19, 716 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  2. Antoja, T., Helmi, A., Romero-Gómez, M., et al. 2018, Nature, 561, 360 [NASA ADS] [CrossRef] [Google Scholar]
  3. Arenou, F., Luri, X., Babusiaux, C., et al. 2018, A&A, 616, A17 [CrossRef] [EDP Sciences] [Google Scholar]
  4. Banik, N., Widrow, L. M., & Dodelson, S. 2017, MNRAS, 464, 3775 [NASA ADS] [CrossRef] [Google Scholar]
  5. Bennett, M., & Bovy, J. 2019, MNRAS, 482, 1417 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bienaymé, O., Robin, A. C., & Crézé, M. 1987, A&A, 180, 94 [Google Scholar]
  7. Bienaymé, O., Famaey, B., Siebert, A., et al. 2014, A&A, 571, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Bienaymé, O., Robin, A. C., & Famaey, B. 2015, A&A, 581, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Binney, J., & Schönrich, R. 2018, MNRAS, 481, 1501 [NASA ADS] [CrossRef] [Google Scholar]
  10. Binney, J., Gerhard, O., & Spergel, D. 1997, MNRAS, 288, 365 [NASA ADS] [CrossRef] [Google Scholar]
  11. Buch, J., Leung, J. S. C., & Fan, J. 2019, JCAP, 2019, 026 [CrossRef] [Google Scholar]
  12. Chequers, M. H., Widrow, L. M., & Darling, K. 2018, MNRAS, 480, 4244 [CrossRef] [Google Scholar]
  13. Crézé, M., Chereul, E., Bienaymé, O., & Pichon, C. 1998, A&A, 329, 920 [Google Scholar]
  14. Debattista, V. P. 2014, MNRAS, 443, L1 [NASA ADS] [CrossRef] [Google Scholar]
  15. Eilers, A.-C., Hogg, D. W., Rix, H.-W., & Ness, M. K. 2019, ApJ, 871, 120 [NASA ADS] [CrossRef] [Google Scholar]
  16. Evans, D. W., Riello, M., De Angeli, F., et al. 2018, A&A, 616, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Famaey, B., & Dejonghe, H. 2003, MNRAS, 340, 752 [NASA ADS] [CrossRef] [Google Scholar]
  18. Faure, C., Siebert, A., & Famaey, B. 2014, MNRAS, 440, 2564 [NASA ADS] [CrossRef] [Google Scholar]
  19. Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Gaia Collaboration (Brown, A. G. A., et al.) 2018a, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Gaia Collaboration (Katz, D., et al.) 2018b, A&A, 616, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Gómez, F. A., Minchev, I., O’Shea, B. W., et al. 2013, MNRAS, 429, 159 [NASA ADS] [CrossRef] [Google Scholar]
  23. GRAVITY Collaboration (Abuter, R., et al.) 2018, A&A, 615, L15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. GRAVITY Collaboration (Abuter, R., et al.) 2019, A&A, 625, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Hagen, J. H. J., & Helmi, A. 2018, A&A, 615, A99 [CrossRef] [EDP Sciences] [Google Scholar]
  26. Haines, T., D’Onghia, E., Famaey, B., Laporte, C., & Hernquist, L. 2019, ApJ, 879, L15 [NASA ADS] [CrossRef] [Google Scholar]
  27. Holmberg, J., & Flynn, C. 2000, MNRAS, 313, 209 [NASA ADS] [CrossRef] [Google Scholar]
  28. Jurić, M., Ivezić, Ž., Brooks, A., et al. 2008, ApJ, 673, 864 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  29. Kapteyn, J. C. 1922, ApJ, 55, 302 [NASA ADS] [CrossRef] [Google Scholar]
  30. Kapteyn, J. C., & van Rhijn, P. J. 1920, ApJ, 52, 23 [NASA ADS] [CrossRef] [Google Scholar]
  31. Kawata, D., Bovy, J., Matsunaga, N., & Baba, J. 2019, MNRAS, 482, 40 [CrossRef] [Google Scholar]
  32. Khoperskov, S., Di Matteo, P., Gerhard, O., et al. 2019, A&A, 622, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Kuijken, K., & Gilmore, G. 1989, MNRAS, 239, 605 [NASA ADS] [CrossRef] [Google Scholar]
  34. Laporte, C. F. P., Johnston, K. V., Gómez, F. A., Garavito-Camargo, N., & Besla, G. 2018, MNRAS, 481, 286 [NASA ADS] [CrossRef] [Google Scholar]
  35. Laporte, C. F. P., Minchev, I., Johnston, K. V., & Gómez, F. A. 2019, MNRAS, 485, 3134 [NASA ADS] [CrossRef] [Google Scholar]
  36. Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, A&A, 616, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Luri, X., Brown, A. G. A., Sarro, L. M., et al. 2018, A&A, 616, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Marrese, P. M., Marinoni, S., Fabrizio, M., & Altavilla, G. 2019, A&A, 621, A144 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. McGaugh, S. S. 2018, Res. Notes Am. Astron. Soc., 2, 156 [NASA ADS] [CrossRef] [Google Scholar]
  40. McMillan, P. J. 2017, MNRAS, 465, 76 [NASA ADS] [CrossRef] [Google Scholar]
  41. Monari, G., Famaey, B., Siebert, A., et al. 2016, MNRAS, 461, 3835 [NASA ADS] [CrossRef] [Google Scholar]
  42. Mróz, P., Udalski, A., Skowron, D. M., et al. 2019, ApJ, 870, L10 [Google Scholar]
  43. Nitschai, M. S., Cappellari, M., & Neumayer, N. 2020, MNRAS, 494, 6001 [CrossRef] [Google Scholar]
  44. Oort, J. H. 1926, Publications of the Kapteyn Astronomical Laboratory Groningen, 40, 1 [Google Scholar]
  45. Oort, J. H. 1932, Bull. Astron. Inst. Neth., 6, 249 [NASA ADS] [Google Scholar]
  46. Piffl, T., Binney, J., McMillan, P. J., et al. 2014, MNRAS, 445, 3133 [NASA ADS] [CrossRef] [Google Scholar]
  47. Read, J. I. 2014, J. Phys. G Nucl. Phys., 41, 063101 [NASA ADS] [CrossRef] [Google Scholar]
  48. Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2014, ApJ, 783, 130 [NASA ADS] [CrossRef] [Google Scholar]
  49. Riello, M., De Angeli, F., Evans, D. W., et al. 2018, A&A, 616, A3 [CrossRef] [EDP Sciences] [Google Scholar]
  50. Robin, A., & Crézé, M. 1986, A&A, 157, 71 [Google Scholar]
  51. Robin, A. C., Reylé, C., Fliri, J., et al. 2014, A&A, 569, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Salgado, J., González-Núñez, J., Gutiérrez-Sánchez, R., et al. 2017, Astron. Comput., 21, 22 [NASA ADS] [CrossRef] [Google Scholar]
  53. Sánchez-Salcedo, F. J., Flynn, C., & de Diego, J. A. 2016, ApJ, 817, 13 [CrossRef] [Google Scholar]
  54. Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829 [NASA ADS] [CrossRef] [Google Scholar]
  55. Shu, F. H. 1969, ApJ, 158, 505 [NASA ADS] [CrossRef] [Google Scholar]
  56. Siebert, A., Bienaymé, O., & Soubiran, C. 2003, A&A, 399, 531 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163 [Google Scholar]
  58. Widmark, A. 2019, A&A, 623, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Widmark, A., & Monari, G. 2019, MNRAS, 482, 262 [NASA ADS] [CrossRef] [Google Scholar]
  60. Widrow, L. M., Gardner, S., Yanny, B., Dodelson, S., & Chen, H.-Y. 2012, ApJ, 750, L41 [NASA ADS] [CrossRef] [Google Scholar]
  61. Widrow, L. M., Barber, J., Chequers, M. H., & Cheng, E. 2014, MNRAS, 440, 1971 [NASA ADS] [CrossRef] [Google Scholar]
  62. Williams, M. E. K., Steinmetz, M., Binney, J., et al. 2013, MNRAS, 436, 101 [NASA ADS] [CrossRef] [Google Scholar]
  63. Xu, Y., Newberg, H. J., Carlin, J. L., et al. 2015, ApJ, 801, 105 [NASA ADS] [CrossRef] [Google Scholar]
  64. Yanny, B., & Gardner, S. 2013, ApJ, 777, 91 [NASA ADS] [CrossRef] [Google Scholar]
  65. Zhang, L., Rix, H.-W., van de Ven, G., et al. 2013, ApJ, 772, 108 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Parameter correlations

Figures A.1 and A.2 represent the triangle plots of the two-dimensional probability density distributions for the five parameters we used to fit density profiles. Most importantly, it appears clear that despite the various (anti-)correlations, chains are converging towards a defined area following the best likelihood, highlighting that there is no degeneracy.

thumbnail Fig. A.1.

Two-dimensional posterior distributions (last steps of the Markov chain) of the five parameters from the adjusted Eq. (2) for the northern hemisphere, represented in a triangle correlation matrix. Colours represent the density in an arbitrary log-scale from purple (high density) to red (low density).

Open with DEXTER

thumbnail Fig. A.2.

Same as Fig. A.1 but for the southern hemisphere.

Open with DEXTER

As expected, the normalised density factor βν (respectively (αν)) is anti-correlated with its scale height hzb (hza). The height of the thick (thin) disk is correlated with its density at R0. In the same spirit, but to a lower level, βν (respectively (αν)) is mildly anti-correlated with its opposite scale height hza (hzb) meaning the height of the thin (thick) disk relies on the density of the thick (thin) disk. The height of the thick disk hzb is correlated with the height of the thin disk hza in the same way the density of the thick disk βν is correlated with the density of the thin disk αν. Lastly, hR presents only very weak correlations with any parameters.

All Tables

Table 1.

Properties of the sectors.

Table 2.

Best-fitting parameters of the density variations.

Table 3.

Best-fitting parameters of the vertical profile for the radial velocity dispersion squared.

Table 4.

Best-fitting parameters of the profile for the vertical velocity dispersion squared.

All Figures

thumbnail Fig. 1.

Part of the colour-magnitude diagram of stars from the Gaia DR2 catalogue with |b|> 45° and G <  15. The selection of our sample of red clump stars from these data is delimited with the dashed lines.

Open with DEXTER
In the text
thumbnail Fig. 2.

Number of stars with respect to their J magnitude. Black line shows stars selected in 2MASS with 0.55 <  J − Ks <  0.7. Grey steps show stars from the previous selection that are not found in Gaia DR2, with a magnitude cut at G = 15.

Open with DEXTER
In the text
thumbnail Fig. 3.

Sub-samples of red clump stars selected from the initial conical selection described in Sect. 2, in Galactic Cartesian coordinates, centred at the solar radius, (x, y) plane (top panel), and (x, z) plane (bottom panel). The five sub-samples delimited in the northern hemisphere (z >  0) are colour-coded. Top panel: cylindrical grid is superimposed in black solid lines. The colours given in the legends are those taking the same notation as described in the text (Sect. 3.1). Bottom panel: distribution of all sub-samples used in this study, 10 volumes, (colour-coded in the north, blue in the south) with respect to the Galactic plane (z = 0). The dashed lines illustrate the intersection of the conical selection and the plane y = 0.

Open with DEXTER
In the text
thumbnail Fig. 4.

Central sectors vertical density profiles. Top panel: density profiles for the northern and southern central sectors (in stars per kpc3) along the |z| axis (in kpc). Light-coloured solid lines follow sliding bins of a constant number of 50 stars while dots are the values of the density in independent bins of fixed size (see text for details). Bottom panel: normalised density differences between the north and south density profiles as a function of |z| (see Eq. (1)), red when the north is denser, and blue otherwise. Vertical error bars illustrate the normalised uncertainties added in quadrature.

Open with DEXTER
In the text
thumbnail Fig. 5.

Same as Fig. 4 but for the sectors shifted in the counter rotation (No and So) on the left and for sectors shifted in the direction of the Galaxy rotation (Nr and Sr) on the right.

Open with DEXTER
In the text
thumbnail Fig. 6.

Same as Fig. 4 but for the sectors shifted in the direction of the Galactic centre (Nc and Sc) on the left and for sectors shifted toward the anti-centre (Na and Sa) on the right.

Open with DEXTER
In the text
thumbnail Fig. 7.

Vertical density profiles binned in |z| for the northern (top panel) and southern (bottom panel) sub-samples, superimposed with the best-fit profile at the three R values (upper line: R0 − 0.8 kpc, middle line: R0, lower line: R0 + 0.8 kpc). The colour code is the same as in Fig. 3.

Open with DEXTER
In the text
thumbnail Fig. 8.

Distribution of red clump samples designed for this study with respect to their GaiaG magnitude. The grey distribution is the main sample used to derive the density profiles (Sect. 3). Distribution in green is a selection from the main sample where only stars with radial velocity measurements are considered in order to obtain velocity dispersion profiles (Sect. 4).

Open with DEXTER
In the text
thumbnail Fig. 9.

Vertical squared velocity dispersion on radial direction profiles binned in |z|, for the five sub-samples in the north (upper panel) and for the five in the south (bottom panel). The colour code per sector is the same as in Fig. 3. Black lines represent the best fit following Eq. (5) and the parameters given in Table 3.

Open with DEXTER
In the text
thumbnail Fig. 10.

Vertical squared velocity dispersion on vertical direction profiles binned in |z|, for the five sub-samples in the north (upper panel), and in the south (bottom panel). The colour code per sector is the same as in Fig. 3. Black curves represent the best fit following Eq. (6) and the parameters given in Table 4, plotted at the three R values (7.54, 8.34, and 9.14 kpc from the highest to the lowest).

Open with DEXTER
In the text
thumbnail Fig. 11.

Surface mass density (Σdyn) in M pc−2 at the Solar position (Eq. (13)) along the direction perpendicular to the Galactic plane (|z| in kpc) for the north (red) and the south (blue). There are 10 000 superimposed curves per hemisphere, which arise from the same number of draws in parameter uncertainties. The shallow red and blue shapes encompass one σ standard deviation of solutions. The two solid lines in red and blue (approximately in the middle of each shape) represent median curves.

Open with DEXTER
In the text
thumbnail Fig. 12.

Vertical density profiles of red clump stars towards the north (left) and the south (right) Galactic poles. Star counts and (continuous lines) best-fit model.

Open with DEXTER
In the text
thumbnail Fig. 13.

Vertical velocity dispersion of red clump stars towards the north (left) and the south (right) Galactic poles. Measured dispersions and the (continuous line) best-fit model.

Open with DEXTER
In the text
thumbnail Fig. 14.

Left: vertical force (in M pc−2 units) perpendicular to the galactic plane versus the vertical distance. Red shows the total vertical force. Blue shows the disc contribution to the vertical force. Black shows the dark matter contribution to the vertical force. Right: surface density in M pc−2 integrated from −z to z. Red shows the total surface density. Blue shows the disc contribution. Black shows the dark matter contribution to the integrated surface density. The differences between the left and right panels can be understood from Eqs. (12) and (13) for each component. The rotation curve of the disc component is decreasing and that of the halo is increasing, leaving the total rotation curve almost flat (albeit slightly decreasing). Therefore, the surface density of the disc is smaller than its corresponding vertical force (blue curves) and the opposite is true for the halo (black curves). The total surface density is slightly smaller than the corresponding vertical force up to z ∼ 2 kpc (red curves), but the trend is then reversed due to the details of the mass distribution at larger heights.

Open with DEXTER
In the text
thumbnail Fig. A.1.

Two-dimensional posterior distributions (last steps of the Markov chain) of the five parameters from the adjusted Eq. (2) for the northern hemisphere, represented in a triangle correlation matrix. Colours represent the density in an arbitrary log-scale from purple (high density) to red (low density).

Open with DEXTER
In the text
thumbnail Fig. A.2.

Same as Fig. A.1 but for the southern hemisphere.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.