Kinematics and dynamics of Gaia red clump stars

We analyse the kinematics and dynamics of a homogeneous sample of red clump stars selected from the second Gaia data release catalogue in the direction of the Galactic poles. The level of completeness of the sample at heights between 0.6 and 3.5 kpc is asserted by comparison with the 2 Micron All Sky Survey catalogue. We show that both the density distribution and velocity dispersion are significantly more perturbed in the North than in the South, in all analysed regions of our Galactic neighbourhoods. We provide a detailed assessment of these North-South asymmetries at large heights. We then proceed to evaluate how such asymmetries could affect determinations of the dynamical matter density under equilibrium assumptions. 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 to estimate, separately in the North and South, the local dark matter density derived under equilibrium assumptions. In the presence of vertical perturbations, such values should be considered as an upper limit. This Jeans analysis yields values of the local dark matter density above 2~kpc, $\rho_{\rm DM} \sim 0.013 \, {\rm M}_\odot/{\rm pc}^3$ ($ \sim 0.509 \, {\rm GeV/cm}^3$) in the perturbed Northern hemisphere, and $\rho_{\rm DM} \sim 0.010 \, {\rm M}_\odot/{\rm pc}^3$ ($ \sim 0.374 \, {\rm GeV/cm}^3$) in the much less perturbed South. As a comparison, we determine the local dark matter density by fitting a global phase-space distribution to the data. We end up with a value in the range of $\rho_{\rm DM} \sim 0.011 - 0.014 \, {\rm M}_\odot/{\rm pc}^3$ in global agreement with Jeans analysis.


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), Oort (1926), and Oort (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 mea-sure 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 et al. 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 et al. 2018a) have made it possible to highlight strong correlations between horizontal and vertical velocities up to z = ±1kpc (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 followup 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 Section 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 Section 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 Section 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 Section 5. An alternative estimate of the dark matter density via the direct modelling of the stellar phase-space distribution is given for comparison in Section 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 Section 7.

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 archive 1 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 sixdimensional 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 degrees, 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 Figure 1). The proxy for the absolute Gaia magnitude in the G band was computed for individual stars, using M G = G − 5 log 10 (1000/ ) + 5, where is the parallax in milli-arcsec. For the completeness motivation (see Section 2.2), only stars with G < 15 were considered.
The vertical cuts in colour G BP -G RP are 1.107 and 1.291. The superior limit is enforced because at larger G BP -G RP , 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, M G , are 0.185 and 0.83 (see Figure 1). These strict limits leave us with 43 589 stars. For the measurement of velocity dispersions (Section 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% . 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 R 0 = 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).

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, M G = 2.3, which is 1.5 magnitude 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 K s = 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 − K s = 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.

Sub-samples
From the sample of red clump stars defined above, several subsamples 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, Figure 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. 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 R min ≤ R ≤ R max with R min = R 0 − 0.4 and R max = R 0 + 0.4. It spans the range in φ in radians such that φ min = −0.4/R 0 and φ max = 0.4/R 0 . The surface of a slice of this sector, perpendicular to the z-axis, is given by φ max (R 2 max − R 2 min ) = 0.64 kpc 2 . 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 degrees, this involves shifting the sector vertically by the largest distance to the Sun reached in the (x, y) plane. It is reached when R = R max and φ = φ max . This distance is given by d = 2R 0 (1 − cos φ max )(R 0 + 0.4) + 0.4 2 = 0.572 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 z min = 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.  Notes. Geometrical limits of the sectors, as defined in the text. The resulting number of stars in each sector is also given. . Sub-samples of red clump stars selected from the initial conical selection described in Section 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. On the top panel, the 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 (Subsection 3.1). The bottom panel illustrates the 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.

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 h i has to be adapted. If stars are sorted in |z| increasing, from j = 1 to j = N stars , in each bin i, h i is calculated as the absolute distance between the j th star and the j + 49 th star. This provides us with an adequate assessment of density variations (see shallow solid lines in Figures 4, 5, and 6). The volume, V i , of a bin is defined by the minimum and maximum radii R min and R max , the minimum and maximum angles φ min and φ max of the sector to which it belongs, 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 Figures 4, 5, and 6.
Top panels in Figures 4, 5, and 6 present the vertical density profiles, in stars per kpc 3 , 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 Section 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 Section 2), z = 0.2|z|/ √ n * , 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 Figures 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 Section 3.3.
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, V z ) 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 Figures 4, 5, and 6 present the normalised density differences in opposite hemispheres (Widrow et al. 2012;Bennett & Bovy 2019): 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 Figure 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)  Finally in the direction of rotation, (Nr, Sr), we find small anticorrelated 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.

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 Article number, page 5 of 16 A&A proofs: manuscript no. 38535corr_reply2 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 (R 0 − 0.8kpc, R 0 , R 0 + 0.8kpc).
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 (L) 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 − 2 ln(L). The following expression using an exponential law in radius and height plus one in height has been retained: The MCMC technique is employed with the five free parameters with flat priors: two normalisation factors α ν , β ν , two scale heights h za and h zb , and the scale length h R . 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 Figure 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.
As noted in Section 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, h za and h zb , 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  Notes. Parameters from Equation 2 derived from the MCMC method for two sets of volumes: the five in the north and the five in the south. The parameters α ν and β ν are the normalisation factors such that the density of stars at zero height and solar radius is ν(R 0 , 0) = α ν + β ν . The parameters h za and h zb are the scale heights and h R is the scale length.
∼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 'thickto-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.

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 (15702 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 Figure 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 Section 2. The second sample in green, where only Article number, page 7 of 16 A&A proofs: manuscript no. 38535corr_reply2  Figure 3 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.
We then apply the same space partition in sectors inside the conic selection as in Section 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 Section 3.2 in order to derive the density in bins of fixed width (see dots in Figures 9 and 10.) We adopt the conventional Galactic cylindrical coordinates (R, φ, z) defined in Section 2. Velocities are converted in this frame (as V R , V φ , V z ) from Gaia DR2 proper motions and radial velocity (µ α , µ δ , V los ) taking the same solar velocity employed in Gaia Collaboration et al. (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).
To calculate the velocity profiles along z in each direction (v R (z),v φ (z),v z (z)), a σ clipping method is applied. For each bin  Fig. 8. Distribution of red clump samples designed for this study with respect to their Gaia G magnitude. The grey distribution is the main sample used to derive the density profiles (Section 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 (Section 4).
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: where n ex is the number of stars excluded by the σ clipping. Their corresponding velocity dispersions are derived with the following equation: Figures 9 and 10 show the distribution of the squared velocity dispersions σ 2 R and σ 2 z in R and z directions, respectively, with respect to the Galactic height. The colour code is the same as in Figure 3. Error bars are representing the classical one σ standard deviation, on 20% of the distance along z ( d = 0.2|z|/ √ n * ) and on the squared dispersion ( σ 2 R = σ 2 R √ 2/(n * − 1) and σ 2 z = σ 2 z √ 2/(n * − 1)). Once vertical dispersion profiles are built, their variations have to be characterised. A similar MCMC machinery as described in Section 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:  Table 3. where a R is the slope in km 2 s −2 kpc −1 and b R is the normalisation in km 2 s −2 . The best-fitting parameters are reported in Table 3. Notes. Parameters obtained from Equation 5 for the northern and southern sub-samples. The parameter a R is the slope of the gradient, b R is the normalisation factor.
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 σ 2 R (z), 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.  Figure 3. Black curves represent the best fit following Equation 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).
To fit the squared vertical velocity dispersion in z, the best choice turns out to be the equation: where a z and b z are normalisation factors in km 2 s −2 . Parameters h σ 2 za and h σ 2 zb are respectively the scale length and the scale height of σ 2 z in kpc. We note that since the b z normalisation factors are typically negative, h σ 2 zb should not be considered along the lines of 'classical' scale heights. Indeed, the vertical velocity dispersion increases with height. Unlike variations of σ 2 R , there is a need for a radial gradient to best fit, σ 2 z , in the different sectors. This is represented by the black lines in Figure 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.

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.

Vertical forces
We define two global scale parameters. They can be directly calculated from the density profile derived in Section 3 (see Equation 2): Vertical forces K z are calculated with the following equation: The velocity dispersion covariance term is given by: Based on Equations 5, 6, and 10, we get: 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), σ 2 R (z) and σ 2 z (R, z). It lets us with a combination of 11 parameters (α ν , h za , β ν , h zb , h R , a R , b R , a z , h σ 2 za , b z , h σ 2 zb ) reflecting a possible parametrization of the data.

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: such that the dynamical surface mass density (Σ dyn = Σ * + Σ DM ), equal to the stellar plus the dark matter mass surface density, is Here, V c (R, z) is the circular velocity related to the radial force K R , through V 2 c = −RK R . In the Galactic plane, this velocity has been set up to V c (R 0 , 0) = 240±8 km.s −1 (See Section 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. 2019 2 ) we extend the standard deviation to 1.5 σ, meaning V c (R 0 , 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 et al. (2019).
The variation of V c 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 V c (R 0 , z) ≈ V c (R 0 , 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 ∂V c (R, z)/∂R = ∂V c (R 0 , 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 (R 0 , V c (R 0 , 0), ∂V c (R 0 , 0)/∂R).
With the addition of the 11 parameters obtained earlier (see Section 5.1), it is possible to derive both vertical (Equation 9) and radial force correction (Equation 12) with respect to galactic radius and height. We then obtain the dynamical surface mass density Σ dyn = Σ * + Σ DM , from Equation 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 10000 draws in the Markov chains for the adjusted parameters and in the standard deviation for fixed quantities. ) along the direction perpendicular to the Galactic plane (|z| in kpc) for the north (red) and the south (blue). There are 10000 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.
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 Section 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.

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 : The density so determined is only slightly smaller (few percents at 2 kpc) than the one in the plane (ρ DM (R 0 )) 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 Figure 11). Thus, for each of the 10000 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: 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.0015M .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.

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 K z 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 phasespace 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 Figure 7 (central north and south) reproduced in Figure 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 degree aperture (extending to 2.8 kpc). We compute the velocity dispersions in 200 pc z-wide bins (see Figure 13). We note that in the study of the K z 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 K z 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.

Fitting a self-consistent model
For the determination of the vertical force K z , 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 K z 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, L z , and I 3 , 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 Figures 12-13). 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.
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 V c (R 0 ) at the Solar Galactic radius as well as by different gradients of the rotation curve at the position of the Sun ∂V c (R 0 , 0)/∂R. These models all have an almost spherical dark halo.

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 R 0 and the slope of the rotation curve at R 0 .
For instance, fixing those values to V c (R 0 ) = 233 km.s −1 (McGaugh 2018) and ∂V c (R 0 , 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 K z , while Figure 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.
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, V c (R 0 ) = 240 km.s −1 and ∂V c (R 0 , 0)/∂R = −2.18 km.s −1 .kpc −1 , we find ρ DM (R 0 ) 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.

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 degrees. 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 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 Equations 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. 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 scalelength 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 phasespace 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 nonstationary effects and, thus, allowing for more solid estimates of the dynamical matter density in the Galactic disc to be obtained. 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).