© The Authors 2024

1. Introduction

The determination of the gravitational potential perpendicular to the Galactic plane has an extensive history covering more than a century, beginning with the work of Kapteyn (1922). He noted the analogy between the Boltzmann equation 1D1V for the vertical motion of stars and the barometric equation for a plane-parallel atmosphere. This work was immediately followed by a first determination of the force perpendicular to the Galactic plane, measured as K_z, by Oort (1932). Numerous studies have followed this thread and continued to this day. A few noteworthy publications that have made it possible to trace many of these prior works include: Kerr & Lynden-Bell (1986), Kuijken (1995), Holmberg & Flynn (2000), Bienaymé et al. (2014), Read (2014), Kramer & Randall (2016), and de Salas & Widmark (2021).

Several important steps have led to significant progress and improvements with respect to the measurement of K_z. The use of the same sample of tracer stars to simultaneously measure its density distribution and kinematics perpendicular to the Galactic plane by Kuijken & Gilmore (1989) has done away with the uncertainties related to the lack of homogeneity of combined samples from previous studies.

The development of a wide range of methods has improved the analysis of stellar samples. The distribution functions, which are the exact solutions of the 1D1V Boltzmann equation developed by Kapteyn (1922) and Oort (1932) taking the form of stationary (also known as isothermal) solutions, remain a satisfactory method that allows for smoothing. A separate approach is based on 1D Jeans equations (Xia et al. 2016; Sánchez-Salcedo et al. 2016) which, in principle, avoids putting an a priori on the solutions and also avoids having to provide an explicit form of the distribution functions. In practice, due to the noise of measurements, a smoothing of the data (density and kinematics as a function of distance from the Galactic plane) proves necessary. Thus, it should be noted that the smoothed curves do not necessarily correspond to a stationary solution.

Another method is the direct inversion, which is undoubtedly very elegant, but which relies on an inverse Laplace transform that is known to be an ill conditioned numerical problem. In practice, just a small amount of noise in the measurements of density and kinematics introduces large fluctuations in the recovered potential. Some examples of this problem can be found in the proceedings of a workshop on the K_z problem (Philip & Lu 1989, see in particular the cover page of these proceedings).

A global approach based on more complete modelling of the Galaxy, combining stellar counts and their kinematics, has also been used to measure the gravitational potential and in particular the force perpendicular to the Galactic plane; for example, in Robin & Crézé (1986), Bienaymé et al. (1987), and Crézé et al. (1989). It has also become more commonly developed today (Nitschai et al. 2020; Sysoliatina & Just 2022; Binney & Vasiliev 2023). Significant advances have been made in the consideration of the coupling of radial and vertical motions (Kuijken & Gilmore 1989). At present, one of these approaches involves 3D modelling based on developing distribution functions dependent on three integrals of motion (Sanders & Binney 2014; Bienaymé et al. 2015; Robin et al. 2022; Binney & Vasiliev 2023). These methods are necessary to fully describe the 3D velocities of stars and to accurately measure the potential beyond 500 pc outside the Galactic plane.

Nowadays, the Gaia DR3 survey is several orders of magnitude more accurate than any other survey in terms of distances, 3D kinematics, definition of homogeneous samples, and number of stars measured. This survey has been used in recent publications to develop new methods for analysing the K_z forces. A key independent observation was the measurement of R₀, the distance from the Sun to the Galaxy’s central black hole (GRAVITY Collaboration 2018, 2019). This has finally enabled the community to ascertain the shape of the rotation curve at different Galactic radii (McGaugh 2018; Mróz et al. 2019; Eilers et al. 2019; Deason et al. 2021). It is both the knowledge of the slope of the Galactic rotation curve and measurement of the K_z, the vertical force, that make it possible (via Poisson’s equation) to determine the total mass density in the solar neighbourhood. An uncertainty in the slope of the rotation curve V_c ∼ R^α equivalent to an uncertainty of about ±0.01 in the α coefficient leads to an uncertainty of ∓0.003 M_⊙/pc³ on the local mass density; this is a significant fraction of the expected dark matter density. This α coefficient was known to this low precision before the GRAVITY measurements.

The Gaia observations have outlined the non-stationarity of the stellar disc, which is a limitation to all the previously mentioned methods. The identification of a spiral (Antoja et al. 2018) in the z-w position-velocity space highlights the velocity fluctuations in the Galactic plane. Widrow et al. (2012) had already noted asymmetries in counts between the directions of the north and south Galactic poles as tracers of vertical waves in the Galactic disc. However, this has opened up a new approach using the properties of spirals in z-w space to measure K_z independently of other methods (Widmark & Monari 2019; Widmark et al. 2021, 2022; Guo et al. 2024).

All of these approaches to measuring the gravitational potential in the disc only determine the total mass density, baryonic and dark matter, without distinguishing between them. As long as we remain close to the Galactic plane, most of the mass is baryonic (stellar and ISM). At R₀ and z = 0 (in Galactic coordinates), we would expect only about 10 per cent of the mass density to be dark matter, assuming a spherical dark halo to explain the Galactic rotation curve. However, the uncertainty on the local stellar and ISM volume mass density was known with a precision of the order of the amount of dark matter. To our knowledge, since the work of Flynn et al. (2006) as well as McKee et al. (2015), it is now only the work based on the local and distant stellar counts Gaia DR3 carried out by Robin et al. (2022) that has provided a more realistic estimate of the local stellar volume and surface mass density, which we refer to later in this paper.

In this work, we propose an improvement to the measurement of dark matter density based on a determination of the gravity field at large distances from the Galactic plane up to 3.5 kpc, where dark matter is dominant. At these distances above the Galactic plane, the contribution of the baryonic mass becomes negligible and non-stationary effects are also expected to be weaker for stellar populations that have large velocity dispersions. The very high precision of Gaia observations considerably simplifies the measurement of the K_z. It is therefore all the more necessary to take proper account of the systematic biases identified (Lindegren et al. 2021; Groenewegen 2021; Khan et al. 2023). These biases are far from negligible for parallaxes at distances of several kpc. Here, we apply the various corrections proposed by Lindegren et al. (2021), which depend on the magnitudes, colours, and equatorial coordinates.

Section 2 details the samples used with the density and kinematic profiles measured. Section 3 presents the dynamical model for measuring the potential, with the method being a reworking of that developed in Bienaymé et al. (2015; 2018) and Robin et al. (2022). A discussion on the baryonic mass distribution is given in Section 4. A general discussion and our conclusion are given in Sections 5 and 6.

2. Data sampling

In this work, the measurement of the gravitational potential and the vertical force K_z is based on the main assumption that the distribution of stars is in a stationary state. A sample that is homogeneous in magnitude and colour is requisite for this measurement, noting that homogeneous refers to the statistical point of view, namely, that the sample would be drawn from a stationary distribution law of a model of the Galaxy. This results in a sample of stars that all have the same properties in terms of colour and absolute magnitude, while the outline of the spatial domain under analysis is precisely defined. The Gaia DR3 catalogue makes it possible to constitute such a sample with an accuracy that has not been achieved thus far for this type of study, thanks to the high accuracy of its measurements.

To probe distances with sufficiently large sample sizes, we select the giant stars of the ‘red clump’. The only significant source of error in this study concerns a bias in the DR3 Gaia parallaxes for the most distant stars. This bias has been analysed by Lindegren et al. (2021) and Groenewegen (2021), among others, revealing systematic errors that depend on the colour, magnitude, and position of the stars. This translates into a relative systematic error that increases with the distance of the stars. For stars located at 4 kpc, this error on the distance can vary from 0 to 40 micro-arcsec (up to 20 per cent on distances at 4 kpc towards the Galactic poles) depending on the stars. We have applied the parallax corrections (Figure 1) proposed by Lindegren et al. (2021). They stress that a similar study should be undertaken for proper motions, for which there must also be a bias. We can admit that this bias must be (at most) of the same order of magnitude as for the parallaxes, which translates into biases of less than 1 km s⁻¹ at a distance of 4 kpc. There, we consider this bias negligible. The sample is defined by the colour-magnitude interval that covers the red clump giants. We have chosen the infrared magnitudes J and K from the 2MASS survey, which minimises the effects of absorption by the interstellar medium; towards the poles, the absorption in the K band is of the order or smaller than 0.012 (Groenewegen 2008). Thus, it is neglected.

Fig. 1. Stellar samples towards the Galactic poles: distances from Gaia DR3 parallaxes versus distances determined with corrected parallaxes (correction of 7% at 3 kpc and 20% at 4 kpc).

We set the absolute magnitude to K_abs = K_mag + 5 * log₁₀(π_corr)−10, where π_corr is the parallax of Gaia DR3 corrected according to the recommendations of Lindegren et al. (2021). The selected samples consist of stars of colour J − K ranging from 0.55 to 0.75 and magnitude K_abs from −1.4 to −1.76 (see Figure 2). Thanks to the precision of colour and distance (Figures 3 and 4), with relative errors of just a few percent, many sources of selection bias have been considerably reduced. Examples include the Malmquist or Lutz-Kelker bias, which are linked to sample limits.

Fig. 2. J an K magnitudes for stars with |z|< 4 kpc and |b|> 30 degrees. The sample selected by the magnitude limits is coloured yellow.

Fig. 3. Inverse of the relative error on parallaxes, π/ϵπ, versus distances (left). Uncertainties on apparent 2MASS K magnitudes (right).

Fig. 4. Uncertainties on transverse velocities from proper motions (left) and on the radial velocities versus distances (right).

The nearly completeness of our sample can be claimed, since our first selection from Gaia-DR3 in magnitude and colour ensure completeness for these directions and far away from the Galactic plane where the stellar density is low with nearly no overlap of stellar images. Then considering the DR3 sample with |b|>22 degrees, bp_rp within [1,1.5], and phot_g_mean_mag < 14, this sample covers our red clump sample and only 0.37% of these stars have no parallax or radial velocity. Only a tiny fraction (10⁻⁴) of the selected sources have no photometric counterpart in the 2MASS catalogue, which is itself complete in these unconfused regions (Skrutskie et al. 2006). Thus, our sample of red clump stars is nearly complete.

Here, we use the Galactic Cartesian coordinates (x, y, z; U, V, W) computed from Topcat functions (Taylor 2005) and the Galactocentric cylindrical coordinates (R, ϕ, z; V_R, V_ϕ, V_z). The position of the Sun is assumed to be R₀ = 8.1 kpc and z₀ = 19 pc. The velocity of the Sun relative to the local standard of rest (LSR) is chosen to be U_⊙ = 12.9 km s⁻¹, V_⊙ = 12 km s⁻¹, W_⊙ = 7.8 km s⁻¹, values which are obtained a posteriori after fitting the velocity distributions of our samples (see Section 3). The circular rotation velocity at R₀ is assumed to be V_c(R₀) = 236.27 km s⁻¹, value adopted in Robin et al. (2022) (see also Nitschai et al. 2020).

Two samples towards the North and South Galactic poles (NGP, SGP) were used to measure the K_z force and are obtained by a selection on the Galactic latitude |b|> 22deg, |R_Gal − 8.1|< 0.4 kpc, |ϕ|< 0.1234, (so R₀  ×  ϕ = 1 kpc), and |z|< 4 kpc. In addition, only the stars with the largest angular momentum are retained, L_z > 8.1 kpc × 100 km s⁻¹, to eliminate most of the halo stars (Figure 5). We note that this last selection is based on an integral of motion.

Fig. 5. Red clump stars towards the Galactic poles: Histogram of star counts versus |z| (left) and versus Vϕ (right). Red line: stars with the highest Lz > 100  ×  8100 km s−1 pc. Black line: all Lz.

For these samples, Figure 6 shows at different heights |z|, (z = 500, 1000, 1500, 2500 pc), as a function of ϕ, the density, ρ, the velocity dispersions, σ_R, σ_z, and ⟨Vϕ⟩. This illustrates the relative constancy of these quantities, an indication of the degree of stationarity.

Fig. 6. Red clump star properties towards the Galactic poles (NGP+SGP): from top left to bottom right, distribution of ρ, ⟨Vϕ⟩, σR, and σz with respect to ϕ at four values of |z| (0.5, 1, 1.5 and 2.5 kpc) (resp. blue, grey, pink, and cyan lines).

A third sample was chosen according to the same selection in colour, magnitude, and L_z, but along the Galactic radius from 4 kpc to 12 kpc and an azimuth extension of ±500 pc. This sample is used to estimate the radial density and kinematic gradients at different z heights (Figure 7).

Fig. 7. Sample distribution (NGP+SGP) of ρ and σz along the Galactic radius. Same colour coding as in Figure 6.

3. Methods and model

3.1. Distribution function

To measure the vertical potential at the solar Galactic position up to a vertical distance of 3.5 kpc, we adapted and modified the century-old method developed by Kapteyn (1922) and Oort (1932). Indeed, their 1D1V method can only be applied when the stellar oscillations through the Galactic plane are smaller than ∼1 kpc, where the vertical motions remain approximately decoupled from the horizontal ones. By building exact 3D3V stationary solutions with Stäckel potentials, Statler (1989) found corrections of the order of 10% at 1 kpc compared with a 1D model.

Thus, for a correct modelling at larger z, we developed stellar population distribution functions (DFs) for an axisymmetric gravitational potential which are 3D generalisations of Shu (1969) DFs. These DFs are inspired by Statler (1989). Our DFs of positions and velocities of each stellar disc are stationary and are modelled with three isolating integrals of motion:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill f(E,L_z,I_3)=g(L_z){\stackrel{\sim }{\rho }}_{0}exp\left(\frac{R_c-R_0}{{\stackrel{\sim }{H}}_{\rho }}\right)exp\left(-\frac{E_{\Vert }}{{\stackrel{\sim }{\sigma }}_R^2}\right)exp\left(-\frac{E_{\perp }}{{\stackrel{\sim }{\sigma }}_z^2}\right),\end{array}\end{array}$$(1)

$$\begin{array}{cc}\hfill \mathrm{with}& {\stackrel{\sim }{\sigma }}_R={\stackrel{\sim }{\sigma }}_{0,R}exp\left(-\frac{R_c-R_0}{{\stackrel{\sim }{H}}_{{\sigma }_R}}\right),\hfill \\ \hfill & {\stackrel{\sim }{\sigma }}_z={\stackrel{\sim }{\sigma }}_{0,z}exp\left(-\frac{R_c-R_0}{{\stackrel{\sim }{H}}_{{\sigma }_z}}\right),\hfill \\ \hfill \mathrm{and}& g(L_z)=\frac{\mathrm{\Omega }}{\kappa }\frac{1}{2\pi {{\sigma }_R}^2}.\hfill \end{array}$$(2)

Here, E_∥ = (E − E_c)(I_3, max − I₃) and E_⊥ = (E − E_c)I₃ are integrals of motion depending on the total energy, E, as well as on E_c, the energy of the circular orbit of angular momentum, L_z, and on I₃, which is an approximate Stäckel integral (details in Bienaymé et al. 2014, 2015, 2018). Also, E_∥ and E_⊥ are linked respectively to the radial and vertical motion of the stars; R₀ is the Galactic radius at solar position and R_c(L_z) is the radius of the circular orbit with angular momentum, L_z. Then, Ω and κ are the circular and epicyclic frequencies expressed as being dependent on R_c(L_z) (see Equation (8) in Bienaymé et al. 2015). With Stäckel potentials, I_3, max = 1 corresponds to the shell orbit and E = E_c + E_∥ + E_⊥. For non Stäckel-potential such as the Galactic potential, I_3, max(Lz) remains close to 1 and must be computed. The distribution function f(E, L_z, I₃) allows us to compute its various moments, which give us the density, ρ(R, z), the rotational velocity, V_ϕ, the velocity dispersions, σ_R, σ_ϕ, and σ_z, and the tilt angle of the velocity ellipsoid. For small velocity dispersions, the density, and the radial and vertical velocity dispersions have a radial exponential decrease. The three first input parameters, ${\stackrel{\sim }{\rho }}_0$, ${\stackrel{\sim }{\sigma }}_{R,0}$, and ${\stackrel{\sim }{\sigma }}_{z,0}$, are directly related (but not exactly equal) to the computed moments of the DFs at the solar position, respectively, the density ρ(R₀, z = 0) and the dispersions σ_R(R₀, z = 0) and σ_z(R₀, z = 0). The other free parameters, ${\stackrel{\sim }{H}}_{\rho }$, ${\stackrel{\sim }{H}}_{{\sigma }_R}$, and ${\stackrel{\sim }{H}}_{{\sigma }_z}$, of the stellar discs are related to the radial scale lengths for the density and for the kinematics.

We note that the DF (Equation (1)) is approximately isothermal, but exactly isothermal in the case of a separable potential in R and z coordinates. It is also similar to the DFs used by Binney & McMillan (2011) in the case of small departures from circular motions where E_∥ ≃ κJ_R and E_⊥ ≃ νJ_z, with κ and ν the epicyclic and vertical frequencies, J_R and J_z the radial and vertical actions.

Finally, we introduce a cut and set to zero the DF (Equation (1)) for orbits with angular momentum larger than L_z, cut = 8.1 kpc × 100 km s⁻¹. Actually, if this DF is effective to model rotating flat stellar discs, it cannot be used to model the stellar halo; neither can it model stars with small angular momentum. Here, we do not try to model the stellar halo with an ad hoc DF. Applying a cut in L_z, that is an integral of motion, allowed us to keep the stationarity of the DF, whereas we did not model stars with small rotational velocities.

3.2. Gravitational potential

The Galactic potential is involved in the distribution functions through the total energy, E and I₃. The contributions to the potential come from the baryonic components, stars, and ISM, and the dark matter. Here, we only constrain the mass density of the dark matter component.

We consider that the stellar component is already well determined by the precise analysis of stellar counts in the immediate solar neighbourhood and at greater distances, obtained recently by fitting the Besançon Model of the Galaxy (hereafter BGM, Robin et al. 2022). On the other hand, the ISM component is less well characterised. We adopt the stellar and ISM contribution to the gravitational potential as that obtained in Robin et al. (2022).

We model the potential of the dark matter halo in such a way that the Galactic circular rotation curve V_c(R) also remains identical to that of the Besançon model of the Galaxy. We set as a free parameter only its flattening, which makes it possible to adjust its local mass density, ρ_dm(R₀, z = 0), without modifying the circular rotation curve. The analytical form of the dark matter halo is therefore slightly different from that of Robin et al. (2022, see their Equation (4)) which comes as follows:

Φ = −930 300.8952 (R_dm² + R² + z²/q_ϕ, h²)^−γ,

with the parameter q_ϕ, h for the flattening of the dark halo. Then, R_dm = 3315 pc, R, and z are expressed in pc and Φ in (km.s⁻¹)². The exponent γ = 0.05 is used to fit the decreasing rotation curve of R from 20 to 50 kpc (Deason et al. 2021; Robin et al. 2022).

We note that if our fitting procedure modifies the density distribution inside the stellar disc, this only concerns stars in the colour-magnitude range of the red clump giants. It would only marginally modify the distribution of the total mass density in the Galactic disc that has been constrained and measured by Robin et al. (2022) and used here to determine the gravitational potential of the Galactic baryonic component.

3.3. Adjustment of density and velocity profiles

The distribution of the observed moments (density, mean Galactic rotational velocity, and the three velocity dispersions) towards the NGP and the SGP as a function of the distance to the Galactic plane, z, are fitted by summing four elementary distribution functions given by Equation (1), each corresponding to stellar discs of different thicknesses (Figure 8). For each of the four elementary DFs, the fitted parameters are the two local velocity dispersions at the solar position, ${\stackrel{\sim }{\sigma }}_{0,R}$ and ${\stackrel{\sim }{\sigma }}_{0,z}$, and its local density, ${\stackrel{\sim }{\rho }}_0$. The fixed parameters are the radial scale lengths for the velocity dispersions and for the density. These fixed parameters are determined by adjusting the corresponding star counts within the Galactic meridian plane (two examples are shown in Figure 7). The Sun velocity relative to the LSR is also a model parameter. It is fitted separately as detailed in a following paragraph.

Fig. 8. Moments of the DFs of the NGP and SGP samples (two blue open circles at each distance) and model (red dot) for the density, ρ, the velocity dispersions, σR, σΦ, σz, and the mean velocity, Vϕ.

We determined the minimum χ² for the parameters of the model using the MINUIT software (James 2004), which allows us to look for possible multiple minima and to obtain a determination of the variance-covariance matrix:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _i}(\frac{{({\rho }_{\mathrm{mod},i}-{\rho }_{\mathrm{obs},i})}^2}{{ϵ}_{1,i}^2}+\frac{{(V_{\varphi ,\mathrm{mod},i}-V_{\varphi ,\mathrm{obs},i})}^2}{{ϵ}_{2,i}^2}\\ \hfill +\frac{{({\sigma }_{R,\mathrm{mod},i}-{\sigma }_{R,\mathrm{obs},i})}^2}{{ϵ}_{3,i}^2}+\frac{{({\sigma }_{\varphi ,\mathrm{mod},i}-{\sigma }_{\varphi ,\mathrm{obs},i})}^2}{{ϵ}_{4,i}^2}\\ \hfill +\frac{{({\sigma }_{z,\mathrm{mod},i}-{\sigma }_{z,\mathrm{obs},i})}^2}{{ϵ}_{5,i}^2})\end{array},\end{array}$$(3)

with ρ_obs, i = N_i, the number of stars in the i-th bin. The uncertainties are for density ${ϵ}_{1,i}=\sqrt{N_i}$; for V_ϕ, ${ϵ}_{2,i}={\sigma }_{\varphi }/\sqrt{N_i}$; for the dispersions σ_R(j = 3), σ_{ϕ(j = 4)} and σ_z(j = 5), ${ϵ}_{j,i}={\sigma }_j/\sqrt{2N_i}$. For each disc, the fitted parameters are related to the values at z = 0 and R₀ = 8.1 kpc of the density and the radial and vertical velocity dispersions, σ_R and σ_z. The flattening q_ϕ, h of the dark halo is also adjusted. The adjustments are performed by least squares method using both the north and south samples simultaneously. We note that the differences between these two distributions are greater than the uncertainty bars of each of these distributions. The steps of the counting intervals are 125 pc (for z ranging from 375 pc to 750 pc), 250 pc (from 750 pc to 3000 pc), and 500 pc (from 3000 pc to 3500 pc).

The fit over the full interval from 375 pc to 3.5 kpc gives the flattening value q_ϕ, h = 0.843 ± 0.009. The other fitted or fixed parameters are given in Table 1. Very small deviations from the fit in density can appear beyond z = 2.5 kpc (Figure 8), which could indicate that the vertical dark matter distribution used (Equation (1)) does not have sufficient degrees of freedom. Otherwise, with the exception of the slight shift at large z of the σ_ϕ distribution (see Figure 8), the fits are very consistent with the data, lying mostly between the north and south distributions, and well within the error bar limits. This validates the distribution functions used. We obtain for the dark matter density ρ(z = 0)_dm = 0.0146 M_⊙ pc⁻³.

We also carried out a fit restricted to the interval 375 pc–2.5 kpc for which we obtained a very close value, q_ϕ, h = 0.839 ± 0.016, but with a higher uncertainty due to the smaller z interval covered. The formal uncertainty, given by the variance-covariance matrix, is very low insofar as the parameter q_ϕ, h is strongly decorrelated from all the other parameters, with the exception of the dispersion σ_z of the thickest disc.

We examine whether this result is sensitive to changes in the model. The first set of fixed model parameters, the scale lengths, are measured using the third stellar sample along the Galactic radius. The radial gradient of σ_R is small and the uncertainty in the kinematic length scale, H_σR does not change the adjusted parameters. On the other hand, varying the values of the H_ρ and H_σz scales within their uncertainties affects the parameters related to the stellar discs, but not the determination of the q_ϕ, h parameter.

By simultaneously decreasing H_ρ and H_σz, we can increase the contribution of stars with higher σ_z dispersions coming from the inner parts (R < R₀). In this case, it has an impact on the measurement of the potential, but it is only significant for extreme changes in H_σz and H_ρ (beyond the uncertainties on these parameters). The largest effect was obtained by increasing the scale length, H_ρ, of the thickest disc from 2.7 kpc to 3.7 kpc. The measurement of q_ϕ, h then rises from 0.84 to 0.85. We therefore consider that this uncertainty indicates that the systematic errors are less than 1% on the measurement of the halo flattening and density. Another effect accurately modelled by the 3D3V model is the inclination of the velocity ellipsoid, which increases with z. This results in an increasing contribution with z of the major axis component of the ellipsoid (equal to σ_R at z = 0) to vertical velocities. This has the effect of increasing σ_z with z in contrast to the 1D1V model, where σ_z remains constant for an isothermal DF. Taking into account the tilt of the ellipsoid requires the velocity dispersions, σ_z and σ_R, to be adjusted simultaneously.

Finally, the result obtained above calls for a correction linked to the slope of the Galactic rotation curve. It is well established in the solar neighbourhood and is dV_c/dR ∼ −1.7 km s⁻¹ kpc⁻¹ (Eilers et al. 2019), namely, d log V_c/d log R = −0.058. The rotation curve of the BGM used here fits the observed curves precisely (see Figure 1 in Robin et al. 2022), but its derivative in R₀ is almost zero. If we had used the correct slope of the rotation curve, this would not change the results concerning the vertical variation of the gravitational potential. However, for the calculation of the local dark matter density, the term (1/R) d[R dΦ/dR)]/dR of the Poisson equation depends on this slope. This introduces a correction term ρ_dm, corr which here is −0.0018 M_⊙ pc⁻³. Taking this correction into account, we finally obtain for the local dark matter density: ρ_dm(0) = 0.0146 − 0.0018 = 0.0128 ± 0.0002 M_⊙ pc⁻³.

Thus, we get the following values for the total local density, local density of baryons and dark matter (Figure 9): ρ_tot, 0 = 0.0844 M_⊙ pc⁻³, ρ_bar, 0 = 0.0716 M_⊙ pc⁻³, and ρ_dm, 0 = 0.0128 M_⊙ pc⁻³.

Fig. 9. Volume mass density at the solar radius position as a function of Galactic height z for baryons (green), dark matter (blue), and total (red).

We obtained for the force, K_z/(2πG), as seen in Figure 10, at different heights, 0.5, 1.1, 2, and 3 kpc, respectively: 45.4, 63.6, 81.5, and 96.4 M_⊙ pc⁻².

Fig. 10. Surface mass density at the solar radius position as a function of Galactic height z for the baryons (green), dark matter (blue), and total (red). The vertical force Kz is superimposed in black.

3.4. Sun velocity relative to the LSR

For the velocity components of the Sun relative to the LSR, we obtained U_⊙ = 12.9 km s⁻¹ and W_⊙ = 7.8 km s⁻¹, deduced from the mean velocity of the samples. On the other side, the V_⊙ component is obtained with the best-fitting of the distribution functions. Then, W_⊙ is quite stable from z = 500 pc to 3 kpc, while the U_⊙ component varies by 3 km s⁻¹. The velocity component V_⊙ = 12 km s⁻¹ is derived by correcting the asymmetric drift, which changes with z and by adjusting V_⊙ so that the observed and modelled distributions are superimposed with by posing $V_{\varphi ,\mathrm{mod}}(z)=V_c(R_0)+V_{\odot }+\overline{V_{\mathrm{obs}}(z)}$ (Figure 8). Then, $\overline{V_{\mathrm{obs}}(z)}$ is the mean observed velocity at different height, z. This is close to the values given by Wang et al. (2021) for the mean solar motion (11.69 ± 0.68, 10.16 ± 0.51,7.67 ± 0.10) km s⁻¹, and Schönrich et al. (2010) who found (11.1 ± 0.69, 12.24 ± 0.47, 7.25 ± 0.37) km s⁻¹. It is also comparable to the values obtained in Robin et al. (2022) (U_⊙ = 10.79 km s⁻¹, V_⊙ = 11.06 km s⁻¹, and W_⊙ = 7.66 km s⁻¹). We note that our method to determine V_⊙ from data at large z is different from these based on the extrapolation of the asymmetric drift velocity to σ_R relation to null velocity dispersion. Such methods are more affected by non-stationary perturbations, whose relative amplitudes are high for small velocity dispersions and at low z values.

4. Discussion on the baryonic surface mass density

Measuring the K_z gives access to the mass density. However, this does not allow us to distinguish between baryonic and dark matter. The slope of the K_z(z) decreases sharply at 600 pc when the density of baryonic matter drops rapidly and it is only around z = 1.4 kpc that the surface density of baryonic matter becomes less than that of the dark matter density, ρ_dm. As a result, determinations of the dark matter density based on the K_z measurements below ∼1.4 kpc are strongly correlated with the value of the integrated baryonic surface mass density, Σ_bar. We refer to Figure 8 in Read (2014) and Figure 13 in Guo et al. (2024). Hence, an error of 10% in the estimated surface mass density of baryons (stars and gas) can produce a systematic error up to 35% in the deduced density of dark matter.

More directly, Caldwell & Ostriker (1981) determined the dark matter density by developing a global mass model of the Galaxy that includes stellar discs and an assumed spherical dark matter halo. The mass of the stellar discs is constrained by stellar counts (e.g. Bahcall & Soneira 1980). Caldwell & Ostriker (1981) obtained a local dark matter density of 0.011 M_⊙ pc⁻³. This was probably the first determination of the local dark matter density after Bosma (1978, 1981a,b, 2023) definitively demonstrated its existence in disc galaxies. A simultaneous combination of the two approaches is the dynamically coherent Galactic model of global stellar population synthesis (Robin & Crézé 1986; Bienaymé et al. 1987; Robin et al. 2022), which also makes it possible to constrain ρ_dm.

4.1. Local baryonic mass density and dark matter halo flattening:

To correctly evaluate the flattening of the dark halo in addition to measuring the dark matter density, it is therefore necessary to have an accurate characterisation of the distribution of the baryonic mass of the Galaxy. To this end, we rely here on the determinations obtained by fitting the BGM (Robin et al. 2022) to the most recent Gaia DR3 observations of positions, parallaxes, proper motions, and radial velocities in a sphere up to about 3 kpc from the Sun. The BGM is a model of stellar population synthesis of the Galaxy, bringing together information on the evolutionary stages and masses of stars as a function of their ages.

The luminosity function, excluding white dwarf stars, was obtained by MCMC fitting to the local sphere (the Gaia catalogue of nearby stars, hereafter GCNS, Gaia Collaboration 2021) and to distant star counts in different directions (see Robin et al. 2022). Assuming null extinction, Figure 11 presents on the top panel the distribution of absolute G magnitudes, taking into account the selection criteria (limiting magnitude G < 17, parallaxes > 10 mas). The bottom panel shows the present-day mass function for stars with mass lower than 0.9 M_⊙, within a sphere of 20 pc radius centred on the Sun, obtained by Kirkpatrick et al. (2024) and from a BGM simulation with the limiting magnitude G < 17 that reduces the number of very low-mass stars.

Fig. 11. Local stellar luminosity function (top). Observations with G < 17 and π > 10 mas are shown in green and BGM model in red. Present-day mass function within 20 pc (in green), see text, and from BGM simulations (in red) with G < 17 (bottom).

In this paper, the adjustment of the counts and kinematics of the red clump giant stars is performed at a a large distance outside the plane at z = 3.5 kpc, in a region where baryonic matter is almost absent. This enables the measurement of the dark matter density to be almost completely decoupled from the value of the surface density of baryonic mass. However, this is a measurement at large z. If we want to correctly determine, at z = 0, the local density of dark matter and the flattening of the dark matter halo, it is still necessary to know precisely the distribution of baryons in the disc. In this case, an uncertainty on the value of Σ_bar also introduces another second uncertainty on the component of the dark halo. In this case, an uncertainty of 10% on the baryonic surface density produces a systematic error of about 10% on the local density of the dark halo. This is because the contributions at (R₀, z = 0) of the baryons and the dark halo to the radial force, K_R, are of the same order of magnitude.

This model combines elements rarely found together in other Galactic models. The dynamical coherence links the gravitational potential to all other Galactic components. The stellar distribution functions, density, and kinematics; hence, it is consistent with the potential. The approach involves less uncertainty than alternative methods previously described for the following reasons: (i) the vertical and horizontal stellar density distributions are dynamically justified and (ii) the stellar luminosity function is constrained by stellar mass luminosity relationships provided by the most recent evolutionary tracks based on Gaia observations. The Starevol stellar evolution models (e.g. Lagarde et al. 2012) cover a wide mass range from 0.6 to 6 solar masses and wide ranges of metallicity and α-abundance, which are incorporated into the BGM by Lagarde et al. (2017). They ensure the reliability of the mass-luminosity relation in this mass range.

The adjustments to this model give the local stellar surface density Σ_* = 33.1 M_⊙ pc⁻². For the ISM contribution, we have assumed Σ_ISM = 11 M_⊙ pc⁻², which, in this case, is the least well known quantity, with an uncertainty of at least 20%, i.e. 2 M_⊙ pc⁻². In the BGM model, the ISM density distribution is a double exponential law with scale length and height: h_R = 700 pc, h_z = 200 pc, and a local density ρ_ISM, 0 = 0.0275 M_⊙ pc⁻³. The value of Σ_ISM adopted here is comparable to that choosen by other authors, who propose a decomposition of the ISM into sub-components. We refer, for instance, to Holmberg et al. (2006) or Table 3 in McKee et al. (2015), but lower than in the case of Read (2014).

Thus, for all baryons, Σ_bar = Σ_* + Σ_ISM = 44.1 ± 2 M_⊙ pc⁻². Our measurement with the sample of red clump giants gives for the local density of the model:

$$\begin{array}{cc}\hfill {\rho }_{\mathrm{tot},0}& =0.0844{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-3},\hfill \\ \hfill {\rho }_{\ast ,0}& =0.0441{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-3},\hfill \\ \hfill {\rho }_{\mathrm{ISM},0}& =0.0275{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-3},\hfill \\ \hfill {\rho }_{\mathrm{bar},0}& =0.0716{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-3},\hfill \\ \hfill \mathrm{and}{\rho }_{\mathrm{dm},0}& =0.0128{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-3}.\hfill \end{array}$$

The stellar density can be compared for instance with the recent determination of ρ_*, 0 = 0.043 ± 0.004 M_⊙ pc⁻³ obtained by McKee et al. (2015). These values are also similar to those recently published, see for example references in Horta et al. (2024). In particular, we obtain values close to Binney & Vasiliev (2023), de Salas & Widmark (2021), McKee et al. (2015), Salomon et al. (2020), and Bland-Hawthorn & Gerhard (2016).

In the light of the discussions developed above, we consider that the analysis of Gaia observations using the Galactic model for the synthesis of stellar populations by Robin et al. (2022) provides a more complete and rigorous approach for the determination of local densities. This is even if the version used here is a simplified, axisymmetric version of the Besançon model.

The formal random error resulting from the adjustments of the moments of the distribution functions of the stars towards the Galactic poles (Figure 8) is very small. On the other hand, the uncertainty of the ISM contribution introduces the largest possible source of systematic error. An uncertainty of 2 M_⊙ pc⁻² in the ISM surface density (i.e. an error of 20%) translates into an uncertainty of 5% in the local surface density of baryonic discs. In order to keep the Galactic rotation curve unchanged, this also means an uncertainty of ±5% on the mass of the dark matter halo. As we determine the dark matter density at z heights between 2 and 3.5 kpc, this implies an uncertainty on the flattening q_ϕ of about 4% and on ρ_dm(z = 0) of 7%. So the flattening of the potential linked to the dark halo is q_ϕ, dm = 0.843 ± 0.035 (or a flattening of the density of the dark halo q_ρ, dm = 0.78 ± 0.06) and ρ_dm(z = 0) = 0.0128 ± 0.0008 M_⊙ pc⁻³.

5. Discussion

With the arrival of the Gaia measurements, considerable progress has been made in our understanding of the Galaxy. The precision and abundance of the Gaia data have clearly revealed many properties and characteristics of the Galactic structure. At the same time, new, more comprehensive and detailed methods of analysis have been rapidly implemented. With regard to the local measurement of the gravity field, we now place our results in the context of recent works.

The most original and unexpected result is certainly the new measurement of the K_z based on the analysis of local non-stationary effects, in particular the spiral structures in the phase space (z, w) visible at z < 1 kpc (Widmark et al. 2021). This method is totally different from the traditional approach of Kapteyn (1922) and Oort (1932), who assumed a stationary state. Using this method Guo et al. (2024) recently obtained Σ_tot(z = 1.1 kpc) = 63 M_⊙ pc⁻², comparable to our results.

In what follows we report published results of gravity field measurements obtained at large heights, z, for the whole Galactic disc. Nitschai et al. (2020) carry out a complete modelling of the Galaxy kinematics over a range from 5 to 12 kpc in radius and up to |z| = 2 kpc based on the older Gaia (DR2) measurements. They derive the density ρ_dm = 0.0115 ± 0.0020 M_⊙ pc⁻³, very close to our values. However, they obtain a very low value of ρ_tot = 0.0640 ± 0.0043 M_⊙ pc⁻³ and an inaccurate measurement of the flattening of the dark halo q = 1.14 ± 0.21. These different results could be explained by the less precise measurements in the Gaia DR2 catalogue.

With Gaia DR3 and APOGEE data, Horta et al. (2024) modelled a large area of the Galaxy, with R ranging from 8 to 11 kpc and z up to ∼1.5 kpc. They got similar results to ours for the local measurements: Σ_⊙(z = 1.1 kpc) = 72${}_{-9}^{+6}$ M_⊙ pc⁻², ρ_tot(z = 0) = 0.081${}_{-0.009}^{+0.015}$M_⊙ pc⁻³, using the most constraining abundance ratio, [Mg/Fe]. This corresponds to a dark matter contribution to the surface density of Σ_dm(z = 1.1 kpc) = 24 ± 4 M_⊙ pc⁻²

Their approach is innovative in the sense that they are able to split the samples into several still highly populated subsamples with various element abundances and also kinematics. These numerous subsamples allow us to put constraints at a range of heights z (unlike other works which use only two or three subsamples, constraining the K_z force at a limited number of z heights). Their second ingenious approach is that they do not fit the moments of the distribution function, density, and vertical velocity dispersion. But instead, they directly adjust the (z, w) distribution identifying the isolevels in this phase space, resulting in stronger constraints on the K_z force. However, they applied a 1D1V dynamical model that we consider is not appropriate to analyse the K_z force at z higher than 1 kpc. They intend to improve this point in a future work.

Cheng et al. (2024) address the problem of measuring the Galactic gravity field for a wide range of Galactic radii and for vertical heights up to 3 kpc. They separate the thin and thick disc populations on the basis of chemical abundances using spectroscopic measurements from the APOGEE survey. They obtain significantly different results for the K_z with the thin or the thick disc samples. They conclude that traditional K_z analyses are ‘challenging’ and highlight the problem of stationarity. However, if we look at their graphs, we can see that the deviations from north-south symmetry and the stationarity problems appear essentially close to the plane, with abrupt changes in the velocity dispersions below z = 500 pc.

In particular, their measurement of the K_z is close to our determination for R between 8 and 9 kpc and for their analysis of the thick disc at large z up to 4 kpc (see their Figures 6 and 7, for R = 8.5 kpc, and their measurement of Σ(z)). We note, however, that their distances based on GSP-phot photometric distances deduced from Gaia BP/RP spectra and parallaxes (all from Gaia DR3) do not seem to have been corrected for the biases mentioned by Lindegren et al. (2021).

Binney & Vasiliev (2023) also carried out a global model of the Galaxy by fitting 3D3V Gaia DR3 measurements within a neighbourhood of 3 kpc from the Sun. The model is based on stellar distribution functions that depend on integrals of motion. Their dark halo is so constructed that it is dynamically consistent with the rest of the components of their Galactic model. They fit 1D projected velocity distributions instead of moments. As the 1D projections of the velocities are strongly non-Gaussian, this prevents valuable information from being discarded. Their results are extremely close to what we obtained:

$$\begin{array}{cc}\hfill {\rho }_{\mathrm{dm},0}& =0.0121{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-3},\hfill \\ \hfill {\mathrm{\Sigma }}_{\mathrm{tot}}(1.1\mathrm{kpc})& =61.2{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-2},\hfill \\ \hfill K_z(1.1\mathrm{kpc})/(2\pi G)& =64.1{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-2}.\hfill \end{array}$$

However, if we consider their Figure 13 (in Binney & Vasiliev 2023), the flattening of their dark matter halo seems to be close to 0.95, nearly spherical. A provisional limitation of the model is that the adjustment of the vertical density of the disc is made using old stellar counts. While justifying this temporary choice, they have left the use of counts taken directly from Gaia observations for a future work.

Ibata et al. (2024) has given an estimate of the halo gravitational potential from the trajectories of stellar streams. They include the largest combination of streams probing the halo and disc gravitational potential leading to extremely small formal errors of the gravitational forces. This method of investigating the potential, now by far the most accurate of all, requires orders of magnitude fewer stars than conventional methods based on the Jeans equations. This translates into extreme precision, because the gravity field is determined at nearly every position of the stars within stellar streams. These authors found values close to ours with ρ_dm, 0 = 0.0114 ± 0007 M_⊙ pc⁻³ and a halo density flattening q_ρ, halo = 0.75 ± 0.03, similar to our value. In fact we notice that they set the vertical force at z = 1.1 kpc as K_1.1/(2πG) = 71 ± 6 M_⊙ pc⁻². This quite certainly stiffens the solution for determining the local dark matter density. Ibata et al. (2024) recognised that tight constraints derived on the local dark matter are clearly (in part) ‘an artefact of the rigidity of the analytic function’. We note that the flattening of the dark halo they obtain is based on a global modelling of the Galactic potential, while ours is more local.

5.1. Stationarity:

All these results remain subject to the stationarity assumption, which is only partially satisfied and whose impact on the study of the K_z has been studied by Read (2014). Widmark & Naik (2024) show that this impact is moderate. They find fluctuations of only about 20%, by tracing the local spiral arms in the Galactic plane through dynamical measurements of the gravity field close to the plane. At high |z|, however, the North-South symmetry is better, and the deviations between the north and south directions remain small for our samples. We have therefore ignored the disequilibria of the Milky Way, considering that these effects are weak (see also the comments in Horta et al. 2024, paragraph 6.4).

5.2. Other models:

The mild flattening of the dark halo and the small variation of the dark matter density between z = 0 and 3.5 kpc do not seem to favour the hypothesis of late accretion of satellite galaxies, which would create a very thick disc or slowly rotating or a flattened spheroidal component of dark matter as suggested by the cosmological simulations of Read et al. (2009) or Pillepich et al. (2014); nor does this slight flattening seem to be very favourable to MOND’s predictions of the presence of a thick disc of phantom matter (Milgrom 2001; Bienaymé et al. 2009). However, truly firm conclusions could only be reached by adjusting the observations of the kinematic-density profiles within the framework of an accurate modelling of the gravitational potentials predicted by these two hypotheses.

6. Conclusion

In this work, we re-examine the gravity field measurements and the mass distribution in the solar neighbourhood towards the Galactic poles. As noted by Cheng et al. (2024): the measured surface density is ‘highly dependent on the assumptions made in its calculation’. This is certainly also true for other measured quantities such as the vertical gravitational forces or the baryonic and dark matter densities, even if the most recent works often publish values that are quite close to each other. Here, we have carried out a study to minimise the dependence of the results on the underlying assumptions.

To this end, we measured the gravity field sufficiently far outside the Galactic plane, where the dark matter mass density largely dominates that of baryonic matter and where the variation of the gravity field depends almost solely on the dark matter density. Thus, we have analysed the gravity field up to z heights of 3.5 kpc, which is significantly higher than in previous studies.

On the other hand, for the accurate modelling of the near-plane gravity field, we assess the mass of stellar discs based directly on the most recent stellar counts and mass-luminosity relations. This provides a more accurate measure of the surface density of baryonic matter, which is also necessary for a more robust determination of the density of dark matter within the Galactic plane. This point is important since it is well established that the local measurement of the dark matter density is correlated with the value adopted for the surface density of baryonic matter (Read 2014; Guo et al. 2024). For that reason, our determination of the distribution of baryonic matter in the solar neighbourhood is based on recent adjustments of the stellar counts and their kinematics, the local luminosity function, and the IMF determined by Robin et al. (2022) and Lagarde et al. (2017).

Finally, for this study, the sample of tracer stars used are red clump giants towards the Galactic poles. By limiting ourselves to bright, ‘nearby’ stars with distances of less than 4 kpc, we have obtained a sample with high measurement accuracy, which profoundly reduces the sources of observational bias. The moments of the observed distributions of the stars in the sample (density, asymmetric drift and velocity dispersions σ_R, σ_ϕ, and σ_z) were fitted with analytical stationary distribution functions. These functions depend on three integrals of motion to properly model the correlations between radial and vertical motions, which is necessary for z heights greater than 1 kpc.

Our results include the measurement of the gravity field from z = 0 to 3.5 kpc, the variation of the dark matter density with z, and the flattening of the dark matter halo:

$$\begin{array}{cc}\hfill {\rho }_{\ast ,0}& =0.0441{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-3},\hfill \\ \hfill {\rho }_{\mathrm{ISM},0}& =0.0275{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-3},\hfill \\ \hfill {\rho }_{\mathrm{bar},0}& =0.0716{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-3},\hfill \\ \hfill {\rho }_{\mathrm{dm}}(z=0)& =0.0128\pm 0.0008{\mathrm{M}}_{\odot }{\mathrm{pc}}^{-3}=0.486\pm 0.030\mathrm{Gev}{\mathrm{cm}}^{-3},\hfill \\ \hfill q_{\varphi ,\mathrm{dm}}& =0.843\pm 0.0035\mathrm{and}{q}_{\rho ,\mathrm{dm}}=0.781\pm 0.0055.\hfill \end{array}$$

We present a model that provides a good fit to the observations. The spatial distribution and kinematics of our stellar sample can be described to recover characteristics of our Galaxy. Further progress can be expected; for instance, instead of simply adjusting the moments of the stellar distribution functions, an adjustment of the 1D projections of the velocities would offer better constraints on the modelling.

Acknowledgments

During the analysis, we have made extensive use of the astronomical java software TOPCAT and STILTS (Taylor 2005). This work has made use of data from the European Space Agency (ESA) mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, http://www.cosmos.esa.int/ web/gaia/dpac/consortium). This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.

References

All Tables

All Figures

	Fig. 1. Stellar samples towards the Galactic poles: distances from Gaia DR3 parallaxes versus distances determined with corrected parallaxes (correction of 7% at 3 kpc and 20% at 4 kpc).
In the text

	Fig. 2. J an K magnitudes for stars with |z|< 4 kpc and |b|> 30 degrees. The sample selected by the magnitude limits is coloured yellow.
In the text

	Fig. 3. Inverse of the relative error on parallaxes, π/ϵπ, versus distances (left). Uncertainties on apparent 2MASS K magnitudes (right).
In the text

	Fig. 4. Uncertainties on transverse velocities from proper motions (left) and on the radial velocities versus distances (right).
In the text

	Fig. 5. Red clump stars towards the Galactic poles: Histogram of star counts versus |z| (left) and versus Vϕ (right). Red line: stars with the highest Lz > 100  ×  8100 km s−1 pc. Black line: all Lz.
In the text

	Fig. 6. Red clump star properties towards the Galactic poles (NGP+SGP): from top left to bottom right, distribution of ρ, ⟨Vϕ⟩, σR, and σz with respect to ϕ at four values of |z| (0.5, 1, 1.5 and 2.5 kpc) (resp. blue, grey, pink, and cyan lines).
In the text

	Fig. 7. Sample distribution (NGP+SGP) of ρ and σz along the Galactic radius. Same colour coding as in Figure 6.
In the text

	Fig. 8. Moments of the DFs of the NGP and SGP samples (two blue open circles at each distance) and model (red dot) for the density, ρ, the velocity dispersions, σR, σΦ, σz, and the mean velocity, Vϕ.
In the text

	Fig. 9. Volume mass density at the solar radius position as a function of Galactic height z for baryons (green), dark matter (blue), and total (red).
In the text

	Fig. 10. Surface mass density at the solar radius position as a function of Galactic height z for the baryons (green), dark matter (blue), and total (red). The vertical force Kz is superimposed in black.
In the text

	Fig. 11. Local stellar luminosity function (top). Observations with G < 17 and π > 10 mas are shown in green and BGM model in red. Present-day mass function within 20 pc (in green), see text, and from BGM simulations (in red) with G < 17 (bottom).
In the text