© ESO 2018

1. Introduction

Recently, the astrophysics community has successfully carried out important ground-based and space missions generating very large data sets. Large sky surveys with photometric and astrometric data, such as Gaia data release 2 (Gaia Collaboration 2018), among others, represent a challenge for Galaxy modelling in terms of both the new types of data and the large amounts of data created. At the same time, new statistical techniques are having a significant effect on modern astronomy. The use of Bayesian statistics for the exploration of large parameter spaces, together with Monte Carlo Markov chains (MCMC) or approximate Bayesian computation (ABC), among others, are in rapid development.

Several attempts have been done to generate fast Milky Way simulations (e.g. Girardi et al. 2005; Jurić et al. 2008; Sharma et al. 2011 or Pasetto et al. 2016). It is demonstrated that the Galaxy models from Sharma et al. (2011) and Pasetto et al. (2016) can be used to explore large parameters spaces under machine-learning algorithms, MCMC, and ABC (e.g. Rybizki & Just 2015; Pasetto et al. 2016). The Galaxia code is able to work in two different modes: it can simulate the Milky Way from a Galaxy model based on Robin et al. (2003), or from N-body simulations. The fast performance of Galaxia relies on its sampling technique and its clever strategy adopted to avoid the simulation of unnecessary stars. The strategy of Pasetto et al. (2016) to perform fast simulations of the Milky Way is based on the use of distribution functions. It constructs colour-magnitude diagrams (CMDs) of single stellar populations from N-body simulations, with low computational cost (Pasetto et al. 2012).

The Besançon Galaxy Model (BGM; Robin et al. 2003) is also a stellar population synthesis model for the Milky Way. It is a very powerful and versatile tool for the statistical analysis of the structure and evolution of the Milky Way. Moreover, it is a valuable tool for the preparation and validation of catalogues for space- and ground-based observational instruments and surveys. Recently, BGM was used to study the kinematics of the local disc from the RAVE survey and the Gaia first data release (DR1; Robin et al. 2017), to evaluate the evolution of the Milky Way’s disc shape over time (Amôres et al. 2017), to constrain Galactic and stellar physics (Lagarde et al. 2017), to constrain the local initial-mass function (IMF) using Galactic Cepheids and Tycho-2 data (Mor et al. 2017), and to study microlensing events in the Galactic bulge (Awiphan et al. 2016). Furthermore, BGM has also been useful, together with other Milky Way models, to study the bulge bar (Simion et al. 2017) and for the validation of the Gaia DR1 (Arenou et al. 2017). Nowadays, a BGM standard (BGM Std) simulation (e.g. Czekaj et al. 2014) has a computational cost that is not adapted to exploring large parameter spaces using modern Bayesian iterative methods that require a very large number of simulations. To overcome this handicap we have developed a theoretical framework to generate very fast Milky Way approximate simulations based on BGM. This framework allows us to explore, among others, the parameter spaces of the IMF, the star formation history (SFH), and the density laws using ABC. The flexibility of the strategy presented here allows for the generation of fast approximate simulations for different Milky Way components, such as thin disc, thick disc, halo, and bulge. Our full strategy is codified to run on Apache Spark¹ (Zaharia et al. 2012) and Apache Hadoop², which are engines coming from business science suited to deal with large surveys. Thanks to its codification, BGM FASt is implemented in the big data infrastructure known as Gaia Data Analytics Framework (GDAF, e.g. Tapiador et al. 2017). As a first application of this complex strategy we use ABC algorithms, BGM, and Tycho-2 data to constrain the IMF, the local SFH, the local stellar mass density, and the thin-disc density laws.

In Sect. 2 we describe the BGM, and in Sect. 3 we present the framework to generate the BGM fast approximate simulations (BGM FASt). In Sect. 4 we describe the treatment of the local dynamical statistical equilibrium in the context of BGM FASt. In Sect. 5 we briefly describe the approximate Bayesian computation technique applied to explore the parameter space of the fundamental functions of the Milky Way. In Sect. 6 we present an evaluation of the BGM FASt performance in the solar neighbourhood. Results are presented in Sect. 7 while a discussion and conclusions are presented in Sects. 8 and 9.

2. The Besançon Galaxy Model

In the present paper we use the following versions of the Galactic components of BGM chosen from a compromise between recent and stable updates.

For the stellar halo component we use the model from Robin et al. (2014) and for the bulge-bar region we use the model described in Robin et al. (2012). For the thick disc component we use the model from the best fit obtained in Robin et al. (2014) which is a thick disc with two main star-formation episodes at 10 and 12 Gyr. For the thin disc component we use the model described in Czekaj et al. (2014) with the updates on the parameters introduced in Mor et al. (2017). The local dynamical statistical equilibrium of BGM is ensured by dynamical constraints based on Bienaymé et al. (1987). The last dynamics and kinematics updates from Bienaymé et al. (2015) and Robin et al. (2017) are not considered in the present paper and will be incorporated in the near future.

2.1. BGM star-generation strategy

The BGM has two main working modes to compute the generation of the stars in the Galaxy. The traditional approach relies on using a precomputed Hess diagram (Robin et al. 2003), while the more updated approach is able to generate the stars from a given set of fundamental functions (e.g. IMF, SFH, age-metallicity among others), making them evolve using a desired set of stellar evolutionary models (Czekaj et al. 2014). For each Galactic component we can choose whether we want to simulate it using the Hess diagram or the updated strategy. Alternatively, from the updated star-generation strategy we can build a Hess diagram from a given set of fundamental functions, and ingest it into BGM code afterwards to be used in a traditional way.

In this section we summarise the stellar generation strategy described in Czekaj et al. (2014), henceforth referred to as our standard strategy. Initially this strategy was developed for the thin-disc component but nowadays it can be used for other Galactic components.

In the BGM Std strategy, to generate stars born τ years ago for a given Galactic i-component (e.g. thin disc, thick disc, halo and bulge-bar), we start from a given total surface mass density at the position of the Sun $({\mathrm{\Sigma }}_{\odot }^i)$. We then use the SFH $({\psi }_{\odot }^i(\tau ))$ to distribute the surface mass density along τ as follows:

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\odot }^i(\tau )\approx {\mathrm{\Sigma }}_{\odot }^i·{\psi }_{\odot }^i(\tau ).\end{array}$$(1)

For simplicity, the current version neglects the radial migration. We setup the model so that the stars are born in the plane. We then redistribute them in the process of secular evolution by using the surface-to-volume mass density ratio at the position of the Sun $(\mathcal{H}(\tau ,{\overline{x}}_{\odot })={\mathrm{\Sigma }}_{\odot }(\tau )/{\rho }_{\odot }(\tau ))$ to compute the volume stellar mass density from ${\mathrm{\Sigma }}_{\odot }^i$, as follows

$$\begin{array}{c}\hfill {\rho }_i(\tau ,x_{\odot },y_{\odot },z_{\odot })=\frac{{\mathrm{\Sigma }}_{\odot }^i(\tau )}{{\mathcal{H}}_i(\tau ,x_{\odot },y_{\odot },z_{\odot })}=\frac{{\mathrm{\Sigma }}_{\odot }^i·{\psi }_{\odot }^i(\tau )}{{\mathcal{H}}_i(\tau ,x_{\odot },y_{\odot },z_{\odot })},\end{array}$$(2)

where we have expressed the position $\overline{x}$ in Cartesian Galactic coordinates as (x, y, z). The volume mass density is distributed throughout the Galaxy as

$$\begin{array}{c}\hfill {\rho }_i(\tau ,x,y,z)={\rho }_i(\tau ,x_{\odot },y_{\odot },z_{\odot })·{\mathcal{R}}_i(\tau ,x,y,z),\end{array}$$(3)

where ℛ_i(τ, x, y, z) are the density laws for the given i-component and ℛ_i(τ, x_⊙, y_⊙, z_⊙)=1. We can then write ℋ_i(τ, x_⊙, y_⊙, z_⊙) as the integral of the density law along the vertical direction at the position of the Sun:

$$\begin{array}{c}\hfill {\mathcal{H}}_i(\tau ,x_{\odot },y_{\odot },z_{\odot })=\frac{{\mathrm{\Sigma }}_{\odot }^i(\tau )}{{\rho }_i(\tau ,x_{\odot },y_{\odot },z_{\odot })}={\displaystyle \underset{\forall z}{\int }{\mathcal{R}}_i(\tau ,x_{\odot },y_{\odot },z)\mathrm{d}z.}\end{array}$$(4)

Finally, from Eqs. (2)–(4), we can write the distribution of the volume mass density along position and age as follows

$$\begin{array}{c}\hfill {\rho }_i(\tau ,\overline{x})=\frac{{\mathrm{\Sigma }}_{\odot }^i·{\psi }_{\odot }^i(\tau )}{{\mathcal{H}}_i(\tau )}·{\mathcal{R}}_i(\tau ,\overline{x}),\end{array}$$(5)

where for simplicity we call ℋ_i(τ) to ℋ_i(τ, x_⊙, y_⊙, z_⊙). As explained in Czekaj et al. (2014), the IMF distributes this mass density in three mass ranges. The star generation process goes through all the volume elements in the Galaxy. First, in a given volume element, for a given age sub-population the age of the star is drawn uniformly within the age limits. Afterwards, the mass of the star is drawn from the IMF. Next, the metallicity is assigned depending on the age and position of the stars. In the most updated versions, the α-elements-to-iron abundance ([α/Fe]) is assigned to each star with a given probability (Lagarde et al. 2017). The evolutionary stage is then assigned to the star by interpolating the stellar evolutionary tracks. In a following step, the process assigns to the generated star a given probability to be the primary component of a stellar multiple system. This probability is assigned following the guidelines of Arenou (2011). Finally, if the star is flagged as a primary component of a stellar multiple system, the standard strategy generates a secondary star with a mass drawn from the probability distributions described in Arenou (2011).

For consideration in the following sections we define a BGM Std simulation to be one that works using the standard stellar generation strategy or a fixed Hess diagram built from the standard stellar generation strategy.

2.2. Standard thin disc component

The thin disc component is described in Czekaj et al. (2014). The stars are generated as described in Sect. 2.1. The thin disc population is divided in seven age sub-populations. Usually the chosen age intervals are those described in Bienaymé et al. (1987) but for the two youngest populations we use the age limits described in Mor et al. (2017). The density distribution of each sub-population of the thin disc is assumed to follow an Einasto density profile as described in Robin et al. (2012), except for the youngest sub-population which follows the expression described in Robin et al. (2003). These profiles are characterised by the eccentricities of the ellipsoid (i.e. the axis ratio), the radial scale length of the disc (h_R), and the radial scale length of the disc hole (h_Rh). A velocity dispersion as a function of age is adopted and the dynamical statistical equilibrium is ensured by using the strategy described in Bienaymé et al. (1987). Stellar evolutionary tracks and model atmosphere, combined with an age-metallicity relation, allow us to go from masses, ages, and metallicities to the space of the observables. In this process a three-dimensional (3D) interstellar extinction map is adopted.

3. Framework for the Besançon Galaxy Model Fast Approximate Simulation

A BGM Std simulation has a computational cost of ∼432 h of CPU time for a simulation of 10⁶ stars, excluding the use of iterative methods like ABC or MCMC to explore large parameter spaces. Hence, we have developed a new method, called BGM FASt, which is able to robustly simulate the Galaxy with a computational cost of ∼240 s of CPU time for a simulation of 10⁶ stars. Thanks to the use of Apache Hadoop and Apache Spark environments (Zaharia et al. 2012) the computational cost should not scale with the number of stars as is the case in standard environments (Julbe, priv. comm.).

3.1. The BGM FASt concept

The BGM FASt is a population-synthesis simulation of the Milky Way, obtained from a clever modification of a BGM Std simulation. The BGM FASt development is based on the distribution function of the generated stars (𝒟_i). The 𝒟_i carries on the information about the generation of the stars in the i-component of the Galaxy (e.g. thin disc, thick disc, halo, bulge-bar) throughout the life of the given component up to the present day. This distribution function contains the chemo-dynamical information that is classically expressed by fundamental functions such as the IMF, the SFH, density distribution, the age-metallicity relation, and the radial metallicity gradient, among others. The 𝒟_i is defined in a N dimensional space (${\mathrm{\P }}^i$) for each of the i-components of the Galaxy. This N dimensional space contains all the parameters that can be involved in a distribution function of the generated stars in the Galaxy. Let us introduce the parameter space as follows:

$$\begin{array}{c}\hfill {\mathrm{\P }}^i\equiv \tau \times M\times Z\times \overline{x}\times \overline{\upsilon }\times \overline{p},\end{array}$$(6)

where τ is the present age of the stellar object, M and Z are its initial mass and metallicity, and $\overline{x}$ and $\overline{\upsilon }$ are position and velocity, respectively. p̄ accounts for other independent parameters that, for some specific purposes, would be interesting to have explicitly introduced in the distribution function. The α-elements-to-iron abundance ratio ([α/Fe]) is an example of one of the possible p̄ parameters and we treat this in the following section.

The strategy to generate a BGM FASt begins with the choice of a specific Mother Simulation with an imposed set of fundamental functions. We use Mother Simulation to refer to a BGM Std simulation used as a seed to generate a BGM FASt simulation. This Mother Simulation is used as a main constituent to generate one or several BGM FASt simulations with different assumptions for the fundamental functions. The idea behind the BGM FASt strategy is that the number of stars generated in a given interval of the parameter space (${\mathcal{N}}_i(\mathrm{\Delta }\mathrm{\P })$) is proportional to the mass dedicated to generate stars for that given interval (${\mathcal{M}}_i(\mathrm{\Delta }\mathrm{\P })$):

$$\begin{array}{c}\hfill {\mathcal{N}}_i(\mathrm{\Delta }\mathrm{\P })\propto {\mathcal{M}}_i(\mathrm{\Delta }\mathrm{\P }),\end{array}$$(7)

where $\mathrm{\Delta }\mathrm{\P }\equiv (\mathrm{\Delta }\tau ,\mathrm{\Delta }M,\mathrm{\Delta }Z,\mathrm{\Delta }\overline{x},\mathrm{\Delta }\overline{\upsilon },\mathrm{\Delta }\overline{p})$. Equation (7) is valid for both the Mother Simulation and the BGM FASt simulation. If $\mathrm{\Delta }\mathrm{\P }$ is small enough, we can write a proportion between the number of stars and the masses relating both the Mother Simulation and the BGM FASt simulation:

$$\begin{array}{c}\hfill \frac{{\mathcal{N}}_i^{\mathrm{FASt}}(\mathrm{\Delta }\mathrm{\P })}{{\mathcal{N}}_i^{\mathrm{MSt}}(\mathrm{\Delta }\mathrm{\P })}\propto \frac{{\mathcal{M}}_i^{\mathrm{FASt}}(\mathrm{\Delta }\mathrm{\P })}{{\mathcal{M}}_i^{\mathrm{MSt}}(\mathrm{\Delta }\mathrm{\P })}·\end{array}$$(8)

Then we can approximate the number of stars for a given interval for a BGM FASt simulation as follows:

$$\begin{array}{c}\hfill {\mathcal{N}}_i^{\mathrm{FASt}}(\mathrm{\Delta }\mathrm{\P })\approx \frac{{\mathcal{M}}_i^{\mathrm{FASt}}(\mathrm{\Delta }\mathrm{\P })}{{\mathcal{M}}_i^{\mathrm{MSt}}(\mathrm{\Delta }\mathrm{\P })}·{\mathcal{N}}_i^{\mathrm{MSt}}(\mathrm{\Delta }\mathrm{\P }).\end{array}$$(9)

Let us call weight to the mass ratio of Eq. (9):

$$\begin{array}{c}\hfill w_i=\frac{{\mathcal{M}}_i^{\mathrm{FASt}}(\mathrm{\Delta }\tau ,\mathrm{\Delta }M,\mathrm{\Delta }Z,\mathrm{\Delta }\overline{x},\mathrm{\Delta }\overline{\upsilon },\mathrm{\Delta }\overline{p})}{{\mathcal{M}}_i^{\mathrm{MSt}}(\mathrm{\Delta }\tau ,\mathrm{\Delta }M,\mathrm{\Delta }Z,\mathrm{\Delta }\overline{x},\mathrm{\Delta }\overline{\upsilon },\mathrm{\Delta }\overline{p})},\end{array}$$(10)

where we compute the mass dedicated to generate stars, in a given interval, from the distribution function of the generated stars that we present in following sections.

In practice, we generate a BGM FASt simulation by applying a weight to each star of the Mother Simulation. This is done according to the parameters of the star, such as mass, age, position, and distance, among others. Thus, the resulting simulation is an approximation of the BGM Std simulation that would be obtained with the standard BGM star generation process.

We present the theoretical framework and the practical implementation for the generation of a BGM FASt as follows.

First, in Sect. 3.2, we describe the distribution function of the generated stars 𝒟_i in the most generic context and its relation with the masses involved in Eq. (10). This allows us to introduce the classical fundamental functions, such as the IMF and the SFH, and functions describing more complex scenarios. In a following step we consider a set of assumptions and approximations to reach a 𝒟_i function compatible with the 𝒟_i implicitly involved in a BGM Std simulation.

Next, in Sect. 3.3, we discuss the treatment of stellar multiple systems as modelled in a BGM. We introduce the probability of obtaining a binary system at birth in our approximated 𝒟_i. The obtained expressions are useful for both the process ensuring the local dynamical statistical equilibrium and the computation of the surface mass stellar density at the position of the Sun. Finally after describing a generalizable weight expression, in Sect. 3.4, we constrain it to the BGM context including stellar multiple systems.

3.2. The distribution function of the generated stars

3.2.1. Generic context and fundamental functions

In this section we describe the distribution function of the generated stars in its most generic context. Under the given definition of the parameter space (${\mathrm{\P }}^i$) we can precisely define some of the parameters belonging to p̄. It is convenient for future purposes to write the distribution function ${\mathcal{D}}_i(\tau ,M,Z,\overline{x},\overline{\upsilon },\overline{p})$ accounting explicitly for the ratio [α/Fe]. Subsequently, [α/Fe] is considered as one of the p̄ parameters and we can write ${\mathcal{D}}_i(\tau ,M,Z,\overline{x},\overline{\upsilon },[\alpha /\mathrm{Fe}],\overline{p\prime })$ in the N dimensional space:

$$\begin{array}{c}\hfill {\mathrm{\P }}^i\equiv \tau \times M\times Z\times \overline{x}\times \overline{\upsilon }\times \overline{p}=\tau \times M\times Z\times \overline{x}\times \overline{\upsilon }\times [\alpha /\mathrm{Fe}]\times \overline{p^{\prime }}.\end{array}$$(11)

In our line of action, for the moment, we are interested in explicitly representing age, mass, metallicity, position, velocity, and the ratio [α/Fe] in the distribution function. Subsequently we marginalize ${\mathcal{D}}_i(\tau ,M,Z,\overline{x},\overline{\upsilon },[\alpha /\mathrm{Fe}],\overline{p\prime })$ over the rest of the p̄′ parameters:

$$\begin{array}{c}\hfill {\mathcal{G}}_i(\tau ,M,Z,\overline{x},\overline{\upsilon },[\alpha /\mathrm{Fe}])={\displaystyle \underset{\forall \overline{p^{\prime }}\in {\mathrm{\P }}^i,}{\int }{\mathcal{D}}_i(\tau ,M,Z,\overline{x},\overline{\upsilon },[\alpha /\mathrm{Fe}],\overline{p^{\prime }})\mathrm{d}\overline{p^{\prime }},}\end{array}$$(12)

where 𝒢_i is the distribution function of the generated stars for the i-component in the reduced space ${\mathrm{\P }}_r^i$:

$$\begin{array}{c}\hfill {\mathrm{\P }}_r^i\equiv \tau \times M\times Z\times \overline{x}\times \overline{\upsilon }\times [\alpha /\mathrm{Fe}].\end{array}$$(13)

For simplicity let us henceforth refer to [α/Fe] using only α.

The ${\mathcal{G}}_i(\tau ,M,Z,\overline{x},\overline{\upsilon },\alpha )$ distribution function is such that the integral over all the parameters that belong to the parameter space is the total number of generated stars in the i-component of the Galaxy:

$$\begin{array}{c}\hfill {\displaystyle \underset{\forall \tau ,M,\overline{Z},\overline{x},\overline{\upsilon },\alpha \in {\mathrm{\P }}_r^i}{\int }{\mathcal{G}}_i(\tau ,M,Z,\overline{x},\overline{\upsilon },\alpha )\mathrm{d}\tau \mathrm{d}M\mathrm{d}Z\mathrm{d}\overline{x}\mathrm{d}\overline{\upsilon }\mathrm{d}\alpha =N_i,}\end{array}$$(14)

and if we multiply by the mass before the integration, we have the total mass of the generated stars for the i-component of the Galaxy:

$$\begin{array}{c}\hfill {\displaystyle \underset{\forall \tau ,M,\overline{Z},\overline{x},\overline{\upsilon },\alpha \in {\mathrm{\P }}_r^i}{\int }{\mathcal{G}}_i(\tau ,M,Z,\overline{x},\overline{\upsilon },\alpha )·M\mathrm{d}\tau \mathrm{d}M\mathrm{d}Z\mathrm{d}\overline{x}\mathrm{d}\overline{\upsilon }\mathrm{d}\alpha =M_i.}\end{array}$$(15)

The mass of Eq. (7) can then be expressed as the integral of the distribution function for a given interval in the parameter space $(\mathrm{\Delta }{\mathrm{\P }}_r^i)$:

$$\begin{array}{cc}& {\mathcal{M}}_i(\mathrm{\Delta }\tau ,\mathrm{\Delta }M,\mathrm{\Delta }Z,\mathrm{\Delta }\overline{x},\mathrm{\Delta }\overline{\upsilon },\mathrm{\Delta }\alpha )\hfill \\ \hfill & ={\displaystyle \underset{\forall \tau ,M,\overline{Z},\overline{x},\overline{\upsilon },\alpha \in \mathrm{\Delta }{\mathrm{\P }}_r^i}{\int }{\mathcal{G}}_i(\tau ,M,Z,\overline{x},\overline{\upsilon },\alpha )·M\mathrm{d}\tau \mathrm{d}M\mathrm{d}Z\mathrm{d}\overline{x}\mathrm{d}\overline{\upsilon }\mathrm{d}\alpha ,}\hfill \end{array}$$(16)

and we can write Eq. (10) as

$$\begin{array}{c}\hfill w_i=\frac{{\displaystyle \underset{\mathrm{\Delta }{\mathrm{\P }}_r^i}{\int }{\mathcal{G}}_i^{\mathrm{FASt}}(\tau ,M,Z,\overline{x},\overline{\upsilon },\alpha )·M\mathrm{d}{\mathrm{\P }}_r}}{{\displaystyle \underset{\mathrm{\Delta }{\mathrm{\P }}_r^i}{\int }{\mathcal{G}}_i^{\mathrm{MSt}}(\tau ,M,Z,\overline{x},\overline{\upsilon },\alpha )·M\mathrm{d}{\mathrm{\P }}_r}},\end{array}$$(17)

where $\mathrm{d}{\mathrm{\P }}_r\equiv \mathrm{d}\tau \mathrm{d}M\mathrm{d}Z\mathrm{d}\overline{x}\mathrm{d}\overline{\upsilon }\mathrm{d}\alpha$.

The true ${\mathcal{G}}_i(\tau ,M,Z,\overline{x},\overline{\upsilon },\alpha )$ distribution function is unknown, but its marginalization over combinations of parameters results in deeply studied functions such as the IMF, the SFH, the age-metallicity relation, the radial metallicity gradient, and also functions carrying information about the density distribution of the Galaxy or dynamical and chemo-dynamical information. Let us exemplify mathematically how some of these fundamental functions can be treated related to the distribution function.

The marginalization over the parameters $\tau ,Z,\overline{x},\overline{\upsilon }$ and α within the values that belong to the ${\mathrm{\P }}_r^i$ space can be written as:

$$\begin{array}{c}\hfill {\displaystyle \underset{\forall \tau ,Z,\overline{x},\overline{\upsilon },\alpha \in {\mathrm{\P }}_r^i}{\int }{\mathcal{G}}_i(\tau ,M,Z,\overline{x},\overline{\upsilon },\alpha )\mathrm{d}\tau \mathrm{d}Z\mathrm{d}\overline{x}\mathrm{d}\overline{\upsilon }\mathrm{d}\alpha ={\xi }_i(M),}\end{array}$$(18)

where ξ_i(M) is the composite IMF for each one of the i-components. Marginalizing now the 𝒢_i over $M,Z,\overline{x},\overline{\upsilon },\alpha$ within the values that belong to the ${\mathrm{\P }}_r^i$ space we have the τ distribution of the generated stars that can be interpreted as the SFH of the whole i-component:

$$\begin{array}{c}\hfill {\displaystyle \underset{\forall M,Z,\overline{x},\overline{\upsilon },\alpha \in {\mathrm{\P }}_r^i}{\int }{\mathcal{G}}_i(\tau ,M,Z,\overline{x},\overline{\upsilon },\alpha )\mathrm{d}M\mathrm{d}Z\mathrm{d}\overline{x}\mathrm{d}\overline{\upsilon }\mathrm{d}\alpha ={\mathrm{\Psi }}_i(\tau ).}\end{array}$$(19)

If we are interested in studying how the τ distribution depends on the position, we can then perform a marginalization over M, Z, $\overline{\upsilon }$ and α, obtaining

$$\begin{array}{c}\hfill {\displaystyle \underset{\forall M,Z,\overline{\upsilon },\alpha \in {\mathrm{\P }}_r^i}{\int }{\mathcal{G}}_i(\tau ,M,Z,\overline{x},\overline{\upsilon },\alpha )\mathrm{d}M\mathrm{d}Z\mathrm{d}\overline{\upsilon }\mathrm{d}\alpha =f_i(\tau ,\overline{x}).}\end{array}$$(20)

If in Eq. (20) we set $\overline{x}={\overline{x}}_{\odot }$, the position of the Sun, then the resulting function can be interpreted as the SFH at the position of the Sun.

The functions involving metallicity, such as the age-metallicity relation or the radial metallicity gradient, could also be considered in BGM FASt. If we marginalize ${\mathcal{G}}_i(\tau ,M,Z,\overline{x},\overline{\upsilon },\alpha )$ over mass, α, and phase-space we get the Z distribution of stars formed τ years ago:

$$\begin{array}{c}\hfill {\displaystyle \underset{\forall M,\overline{x},\overline{\upsilon },\alpha \in {\mathrm{\P }}_r^i}{\int }{\mathcal{G}}_i(\tau ,M,Z,\overline{x},\overline{\upsilon },\alpha )\mathrm{d}M\mathrm{d}\overline{x}\mathrm{d}\overline{\upsilon }\mathrm{d}\alpha ={\chi }_i(\tau ,Z).}\end{array}$$(21)

The radial metallicity gradient can be deduced from a more complex expression obtained marginalizing G_i over age, mass, velocity and α:

$$\begin{array}{c}\hfill {\displaystyle \underset{\forall \tau ,M,\overline{\upsilon },\alpha \in {\mathrm{\P }}_r^i}{\int }{\mathcal{G}}_i(\tau ,M,Z,\overline{x},\overline{\upsilon },\alpha )\mathrm{d}\tau \mathrm{d}M\mathrm{d}\overline{\upsilon }\mathrm{d}\alpha ={\eta }_i(Z,\overline{x}).}\end{array}$$(22)

This latter expression is the position and metallicity distribution of the generated stars throughout the life of the i-component.

Finally, information about chemo-dynamics and kinematics can be introduced with the following two equations.

$$\begin{array}{cc}\hfill {\displaystyle \underset{\forall M,\overline{x},\alpha \in {\mathrm{\P }}_r^i}{\int }{\mathcal{G}}_i(\tau ,M,Z,\overline{x},\overline{\upsilon },\alpha )\mathrm{d}M\mathrm{d}\overline{x}\mathrm{d}\alpha }& ={\mathcal{Q}}_i(\tau ,Z,\overline{\upsilon }),\hfill \end{array}$$(23)

$$\begin{array}{cc}\hfill {\displaystyle \underset{\forall M,Z,\alpha \in {\mathrm{\P }}_r^i}{\int }{\mathcal{G}}_i(\tau ,M,Z,\overline{x},\overline{\upsilon },\alpha )\mathrm{d}M\mathrm{d}Z\mathrm{d}\alpha }& ={\mathcal{K}}_i(\tau ,\overline{x},\overline{\upsilon }).\hfill \end{array}$$(24)

The spatial distribution of the volume mass density (${\rho }_g^i(\overline{x})$) that has been dedicated to generating stars for the Galactic i-component throughout its life can be written as follows

$$\begin{array}{c}\hfill {\rho }_g^i(\overline{x})={\displaystyle \underset{\forall \tau ,M,\overline{Z},\overline{\upsilon },\alpha \in {\mathrm{\P }}_r^i}{\int }{\mathcal{G}}_i(\tau ,M,Z,\overline{x},\overline{\upsilon },\alpha )·M\mathrm{d}\tau \mathrm{d}M\mathrm{d}Z\mathrm{d}\overline{\upsilon }\mathrm{d}\alpha .}\end{array}$$(25)

In the following steps it is useful to have the equation of the mass density dedicated to generating stars born τ years ago in a position $\overline{x}$:

$$\begin{array}{c}\hfill {\rho }_g^i(\tau ,\overline{x})={\displaystyle \underset{\forall M,\overline{Z},\overline{\upsilon }\in {\mathrm{\P }}_r^i}{\int }{\mathcal{G}}_i(\tau ,M,Z,\overline{x},\overline{\upsilon },\alpha )·M\mathrm{d}M\mathrm{d}\overline{Z}\mathrm{d}\overline{\upsilon }\mathrm{d}\alpha .}\end{array}$$(26)

Until now, we have described the distribution function of the generated stars in the Galaxy in a generic context. We have emphasised that our strategy can be generalizable and any of the fundamental equations described above (Eqs. (18)–(24)) can be used if we are able to write an analytical expression for them.

3.2.2. The approximate solution

In this section we find an approximation to the distribution function of the generated stars compatible with BGM. At the same time this approximation aims to be extensible to other models of the Galaxy that use similar star-generation strategies. As the exploration of the velocity spaces is not included in the present paper, for simplicity we do not consider the kinematic part here. Our goal is therefore to find an approximate solution to the integral ${\int }_{\forall \overline{\upsilon }\in {\mathrm{\P }}_r^i}{\mathcal{G}}_i(\tau ,M,Z,\overline{x},\overline{\upsilon },\alpha )\mathrm{d}\overline{\upsilon }$.

In this context the first assumption comes from a traditional strategy (e.g. Tinsley 1980), assuming that mass and age distributions are separated. Splitting the mass function from the function of τ and $\overline{x}$ as follows

$$\begin{array}{c}\hfill {\xi }_i(M)·{\mathcal{F}}_i(\tau ,\overline{x}).\end{array}$$(27)

Moreover, the BGM assumes a metallicity distribution that depends on position and age. We can therefore introduce, in the equations, the probability that a star of a given age in a given position has a metallicity Z: ${\mathcal{P}}_i(Z|\tau ,\overline{x})$. Furthermore, BGM has recently included the possibility to use [α/Fe], a parameter which affects the stellar evolutionary tracks (Lagarde et al. 2017), assuming, from observational surveys, a certain probability that a given star has a given [α/Fe]. In general, for a given i-component, this probability depends on the age, the position, and the metallicity of the star and we denote it as ${\mathcal{P}}_i(\alpha |\tau ,Z,\overline{x})$.

From Eq. (27), and the metallicity and [α/Fe] distributions, we can write:

$$\begin{array}{cc}\hfill {\displaystyle \underset{\forall \overline{\upsilon }\in {\mathrm{\P }}_r^i}{\int }{\mathcal{G}}_i(\tau ,M,Z,\overline{x},}& \overline{\upsilon },\alpha )\mathrm{d}\overline{\upsilon }\hfill \\ \hfill & \approx {\xi }_i(M)·{\mathcal{F}}_i(\tau ,\overline{x})·{\mathcal{P}}_i(Z|\tau ,\overline{x})·{\mathcal{P}}_i(\alpha |\tau ,Z,\overline{x}),\hfill \end{array}$$(28)

where the assumptions and approximations behind the mathematical expression of the functions ${\mathcal{F}}_i(\tau ,\overline{x})$, ${\mathcal{P}}_i(Z|\tau ,\overline{x})$ and $\mathcal{P}(\alpha |\tau ,Z,\overline{x})$ are imposed to be compatible with the BGM. Equation (28) assumes, from the statistical point of view, that the probability to generate a star with mass M and the probability to generate a star τ years ago in a given position are conditionally independent. This means that the IMF is assumed to be independent of time and position.

The standard star-generation strategy described in Sect. 2.1 guides us by using Eq. (5) to approximate the function ${\mathcal{F}}_i(\tau ,\overline{x})$ as follows

$$\begin{array}{c}\hfill {\mathcal{F}}_i(\tau ,\overline{x})\approx \frac{{\mathrm{\Sigma }}_{\odot }^i·{\psi }_{\odot }^i(\tau )}{{\mathcal{H}}_i(\tau )}·{\mathcal{R}}_i(\tau ,\overline{x}),\end{array}$$(29)

where ${\mathrm{\Sigma }}_{\odot }^i$ is the stellar surface mass density ( * /pc²) of the generated stars at the position of the Sun for the Galactic i-component, ℋ_i(τ) is the surface-to-volume-density ratio at the position of the Sun, ${\mathcal{R}}_i(\tau ,\overline{x})$ is the density distribution, and ${\psi }_{\odot }^i$ is the SFH in the solar neighbourhood.

Taking into account the approximations implicitly or explicitly adopted in standard BGM we present a solid solution for the following integral

$$\begin{array}{cc}\hfill {\displaystyle \underset{\forall \overline{\upsilon }\in {\mathrm{\P }}_r^i}{\int }{\mathcal{G}}_i}& (\tau ,M,Z,\overline{x},\overline{\upsilon },\alpha )\mathrm{d}\overline{\upsilon }\hfill \\ \hfill & \approx \frac{{\mathrm{\Sigma }}_{\odot }^i·{\psi }_{\odot }^i(\tau )}{{\mathcal{H}}_i(\tau )}·{\mathcal{R}}_i(\tau ,\overline{x})·{\xi }_i(M)·{\mathcal{P}}_i(Z|\tau ,\overline{x})·{\mathcal{P}}_i(\alpha |\tau ,Z,\overline{x}),\hfill \end{array}$$(30)

where the IMF is normalized as ${\int }_{\forall M\in {\mathrm{\P }}_r^i}\xi (M)·M\mathrm{d}M=1$, the SFH in the solar neighbourhood is normalized as ∫_{∀ τ ∈ Πr}ψ_⊙(τ) dτ = 1, and by definition ${\mathcal{P}}_i(Z|\tau ,\overline{x})$ and ${\mathcal{P}}_i(\alpha |\tau ,Z,\overline{x})$ are normalized to 1. ${\mathcal{R}}_i(\tau ,\overline{x})$ and ℋ_i(τ) are such that the integral over the whole parameter space of ${\mathcal{G}}_i(\tau ,M,Z,\overline{x},\overline{\upsilon },\alpha )$ gives the total number of generated stars in the Galactic i-component.

For simplicity, as BGM Std does, we assume an axi-symmetric structure with no radial migration. Additionally, we assume that for a given volume element, $\frac{\mathrm{d}\rho }{\mathrm{d}t}\approx 0$, where ρ is the stellar volume mass density. These assumptions allow us to use the density distributions (${\mathcal{R}}_i(\tau ,\overline{x})$) described in Robin et al. (2003) and Robin et al. (2014), derived to match the present density distribution of each Galactic i-component.

3.3. Handling of the stellar multiple systems

The probability of obtaining a multiple system at birth in a star formation event is unknown, very complex, and can depend on many parameters, such as the metallicity, the mass of the molecular cloud, the turbulence, mass segregation, mass competition, and mass accretion rate (e.g. Bonnell et al. 2007; Kroupa et al. 2013). For these reasons, the BGM approach (e.g. Czekaj et al. 2014; Robin et al. 2012) is performed adopting an empirical law to match the observed present distribution function of multiple stellar systems (Arenou 2011). This approach assigns a certain probability to a generated star of being the primary component of a multiple system, depending on its mass and its luminosity class. Inside our framework this means a dependence on mass, age, and, through the stellar evolutionary models, Z and α. This probability is therefore P(bin|τ, M, Z, α).

Let us refer to the generated star susceptible to be flagged as either the single star or primary component of the binary system as the “primal star”. Once a primal star is decided as being a primary component of a multiple system, a secondary star is generated with a mass m. The mass m of the secondary is assigned following a probability distribution function (PDF) that depends mainly on the mass M and the luminosity class of the primary component, leading to P(m|τ, M, Z, α).

We introduce the stellar multiple systems in BGM FASt by simplifying their treatment under two assumptions: (1) The masses of primal stars (singles and primaries) are drawn from the IMF while the mass of the secondary follows the empirical laws described above, and (2) the age of the primary star follows the SFH, while the age of the secondary is assumed to be the same as the primary (i.e. both components were born together). With these two assumptions, the multiple stellar systems can be introduced in our analytical approach at low computational cost. We can compute first, at the position of the Sun, the surface mass density of secondary stars (${\mathrm{\Sigma }}_{\odot }^{i,\$}$) as a function of the surface mass density of the primal stars (${\mathrm{\Sigma }}_{\odot }^{i,▵}$) and afterwards compute ${\mathrm{\Sigma }}_{\odot }^{i,▵}$ as a function of ${\mathrm{\Sigma }}_{\odot }^i$. We begin by expressing ${\mathrm{\Sigma }}_{\odot }^i$ as

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\odot }^i={\mathrm{\Sigma }}_{\odot }^{i,▵}+{\mathrm{\Sigma }}_{\odot }^{i,\$}.\end{array}$$(31)

Subsequently, the stellar volume mass density of primary stars at the position of the sun for the i-component is given by

$$\begin{array}{cc}\hfill {\rho }_g^{i,▵}({\overline{x}}_{\odot })& \approx {\displaystyle \underset{\forall \tau ,M,Z,\alpha \in {\mathrm{\P }}_r^i}{\int }\frac{{\mathrm{\Sigma }}_{\odot }^{i,▵}·{\psi }_{\odot }^i(\tau )}{{\mathcal{H}}_i(\tau )}·\mathcal{R}(\tau ,{\overline{x}}_{\odot })·{\xi }_i(M)·M}\hfill \\ \hfill & ·{\mathcal{P}}_i(Z|\tau ,{\overline{x}}_{\odot })·{\mathcal{P}}_i(\alpha |\tau ,Z,{\overline{x}}_{\odot })\mathrm{d}\tau \mathrm{d}M\mathrm{d}Z\mathrm{d}\alpha ,\hfill \end{array}$$(32)

and the stellar volume mass density of secondary stars at the position of the sun for the i-component is given by

$$\begin{array}{cc}\hfill {\rho }_g^{i,\$}({\overline{x}}_{\odot })& \approx {\displaystyle \underset{\forall \tau ,M,Z,\alpha ,\in {\mathrm{\P }}_r^i}{\int }\underset{m_{min}}{\overset{m_{max}}{\int }}\frac{{\mathrm{\Sigma }}_{\odot }^{i,▵}·{\psi }_{\odot }^i(\tau )}{{\mathcal{H}}_i(\tau )}·\mathcal{R}(\tau ,{\overline{x}}_{\odot })·{\xi }_i(M)}\hfill \\ \hfill & ·{\mathcal{P}}_i(Z|\tau ,{\overline{x}}_{\odot })·{\mathcal{P}}_i(\alpha |\tau ,Z,{\overline{x}}_{\odot })·P(\mathrm{bin}|\tau ,M,Z,\alpha )\hfill \\ \hfill & ·P(m|\tau ,M,Z,\alpha )m\mathrm{d}m\mathrm{d}\tau \mathrm{d}M\mathrm{d}Z\mathrm{d}\alpha .\hfill \end{array}$$(33)

We also want to compute the surface mass density; then, if we leave out ℋ_i from the integrand of Eq. (33), we get the surface mass density of secondary stars at the position of the Sun:

$$\begin{array}{cc}\hfill {\mathrm{\Sigma }}_{\odot }^{i,\$}& \approx {\mathrm{\Sigma }}_{\odot }^{i,▵}{\displaystyle \underset{\forall \tau ,M,Z,\alpha \in {\mathrm{\P }}_r^i}{\int }\underset{m_{min}}{\overset{m_{max}}{\int }}{\psi }_{\odot }^i(\tau )·{\xi }_i(M)·{\mathcal{P}}_i(Z|\tau ,{\overline{x}}_{\odot })}\hfill \\ \hfill & ·P_i(\alpha |\tau ,Z,{\overline{x}}_{\odot })·P(\mathrm{bin}|\tau ,M,Z,\alpha )\hfill \\ \hfill & ·P(m|\tau ,M,Z,\alpha )·m\mathrm{d}m\mathrm{d}\tau \mathrm{d}M\mathrm{d}Z\mathrm{d}\alpha .\hfill \end{array}$$(34)

Finally, ${\mathrm{\Sigma }}_{\odot }^{i,▵}$ is computed from ${\mathrm{\Sigma }}_{\odot }^i$ using Eqs. (31) and (34).

This result is very useful for both determining the weight function when considering multiple stellar systems as implemented in BGM (Sect. 3.4) and computing the local dynamical statistical equilibrium (Sect. 4). In Sect. 4 we also particularize the computation of ${\mathrm{\Sigma }}_{\odot }^{i,▵}$ when specifically using the thin disc component as described in Czekaj et al. (2014).

3.4. The weight

The strategy developed in previous sections allows us to give the analytical expression to compute the weights. These weights are able to transform the distribution of stars in the Mother Simulation, linked to a given set of Galactic fundamental functions, into the distribution of the stars linked to other Galactic fundamental functions adopted for a BGM FASt simulation. From (17) and (30) we obtain an expression for the weights applicable to Galaxy models that sample the stars from a distribution function of the form of Eq. (28), considering the assumptions discussed in Sect. 3:

$$\begin{array}{cc}& w_i(\tau ,M,Z,\overline{x},\alpha )\approx \frac{\underset{\mathrm{\Delta }{\mathrm{\P }}_r^i}{\int }\frac{{\mathrm{\Sigma }}_{\odot }^{i,FASt}·{\psi }_{\odot }^{i,FASt}}{{\mathcal{H}}_i^{\mathrm{FASt}}(\tau )}·{\mathcal{R}}_i^{\mathrm{FASt}}(\tau ,\overline{x})}{\underset{\mathrm{\Delta }{\mathrm{\P }}_r^i}{\int }\frac{{\mathrm{\Sigma }}_{\odot }^{i,MSt}·{\psi }_{\odot }^{i,MSt}}{{\mathcal{H}}_i^{\mathrm{MSt}}(\tau )}·{\mathcal{R}}_i^{\mathrm{MSt}}(\tau ,\overline{x})}\hfill \\ \hfill & ·\frac{{\xi }_i^{\mathrm{FASt}}(M)·M·{\mathcal{P}}_i^{\mathrm{FASt}}(Z|\tau ,\overline{x})·{\mathcal{P}}_i^{\mathrm{FASt}}(\alpha |\tau ,Z,\overline{x})·\mathrm{d}{\mathrm{\P }}_r}{{\xi }_i^{\mathrm{MSt}}(M)·M·{\mathcal{P}}_i^{\mathrm{MSt}}(Z|\tau ,\overline{x})·{\mathcal{P}}_i^{\mathrm{MSt}}(\alpha |\tau ,Z,\overline{x})·\mathrm{d}{\mathrm{\P }}_r}·\hfill \end{array}$$(35)

In the BGM context, including multiple stellar systems and taking into account the scenario described in Sect. 3.3, the weight that we apply in practice to generate a BGM FASt simulation is the following

$$\begin{array}{cc}& w_i(\tau ,M,Z,\overline{x},\alpha )\approx \frac{\underset{\mathrm{\Delta }{\mathrm{\P }}_r^i}{\int }\frac{{\mathrm{\Sigma }}_{\odot }^{i,▵,FASt}·{\psi }_{\odot }^{i,FASt}}{{\mathcal{H}}_i^{\mathrm{FASt}}(\tau )}·{\mathcal{R}}_i^{\mathrm{FASt}}(\tau ,\overline{x})}{\underset{\mathrm{\Delta }{\mathrm{\P }}_r^i}{\int }\frac{{\mathrm{\Sigma }}_{\odot }^{i,▵,MSt}·{\psi }_{\odot }^{i,MSt}}{{\mathcal{H}}_i^{\mathrm{MSt}}(\tau )}·{\mathcal{R}}_i^{\mathrm{MSt}}(\tau ,\overline{x})}\hfill \\ \hfill & ·\frac{{\xi }_i^{\mathrm{FASt}}(M)·M·{\mathcal{P}}_i^{\mathrm{FASt}}(Z|\tau ,\overline{x})·P_i^{\mathrm{FASt}}(\alpha |\tau ,Z,\overline{x})·\mathrm{d}{\mathrm{\P }}_r}{{\xi }_i^{\mathrm{MSt}}(M)·M·{\mathcal{P}}_i^{\mathrm{MSt}}(Z|\tau ,\overline{x})·{\mathcal{P}}_i^{\mathrm{MSt}}(\alpha |\tau ,Z,\overline{x})·\mathrm{d}{\mathrm{\P }}_r},\hfill \end{array}$$(36)

where we have substituted ${\mathrm{\Sigma }}_{\odot }^i$ with ${\mathrm{\Sigma }}_{\odot }^{i,▵}$ from Eq. (35).

In our approach, the stellar evolutionary tracks are the same for both the Mother Simulation and the BGM FASt simulation. For simplicity, the probabilities ${\mathcal{P}}_i(Z|\tau ,\overline{x})$ and ${\mathcal{P}}_i(\alpha |\tau ,Z,\overline{x})$ are imposed and we are not going to explore them; they are equal for both the Mother Simulation and the BGM FASt simulation and we can marginalize the numerator and denominator over Z and α to obtain the final expression for the weights:

$$\begin{array}{cc}\hfill w_i(\tau ,M,\overline{x})& \approx \frac{\underset{\mathrm{\Delta }\tau ,\mathrm{\Delta }M,\mathrm{\Delta }\overline{x}}{\int }\frac{{\mathrm{\Sigma }}_{\odot }^{i,▵,FASt}·{\psi }_{\odot }^{i,FASt}}{{\mathcal{H}}_i^{\mathrm{FASt}}(\tau )}·{\mathcal{R}}_i^{\mathrm{FASt}}(\tau ,\overline{x})}{\underset{\mathrm{\Delta }\tau ,\mathrm{\Delta }M,\mathrm{\Delta }\overline{x}}{\int }\frac{{\mathrm{\Sigma }}_{\odot }^{i,▵,MSt}·{\psi }_{\odot }^{i,MSt}}{{\mathcal{H}}_i^{\mathrm{MSt}}(\tau )}·{\mathcal{R}}_i^{\mathrm{MSt}}(\tau ,\overline{x})}\hfill \\ \hfill & ·\frac{{\xi }_i^{\mathrm{FASt}}(M)·M·\mathrm{d}\tau \mathrm{d}M\mathrm{d}\overline{x}}{{\xi }_i^{\mathrm{MSt}}(M)·M·\mathrm{d}\tau \mathrm{d}M\mathrm{d}\overline{x}}.\hfill \end{array}$$(37)

In practice, the weight in Eq. (37) is applied to each single star and each stellar system. This means that both the primary and the secondary components of the stellar system are weighted with the same value, according to the parameters of the primary component. We choose this way to apply the weights because in BGM Std the mass and age of the secondary star are drawn according to the mass and the age of the primary star. This choice ensures that the mass and age distributions of the primal stars in a BGM FASt simulation follow the pertinent IMF and SFH, and that the mass distribution of the secondary stars follows the empirical distributions described in Arenou (2011).

The choice of the intervals for the integrals must be discussed for each particular case. Generally they must be small enough to preserve Eqs. (8) and (9). For tests and cases presented in the present paper, the choice of the mass interval is set to be very small (0.025 M_⊙). For a given star of age τ we set the limits of the integral to be the age limits of the age sub-populations to which the star belongs. Finally, as the BGM Std generation strategy assigns the density of its centre to the whole volume element, in BGM FASt we do not need to perform a volume element integral.

4. Local dynamical statistical equilibrium (LDSE)

By local dynamical statistical equilibrium we understand that, at the position of the Sun, the mass density distribution and the potential satisfy both the Poisson equation and the first-order moment of the collisionless Boltzmann equation for the vertical direction (Bienaymé et al. 1987). In this approach we consider axi-symmetry to solve both equations. To solve the first-order moment of the collisionless Boltzmann equation in the vertical direction, we also assume steady state, isothermal state, and decoupled radial and vertical motions. One possible methodology to ensure the LDSE is described in Czekaj et al. (2014) and summarised in Sect. 4.1. As it requires a computational time that is not affordable for the methodology presented in this paper, we develop here analytical (Sect. 4.2) and approximate (Sect. 4.4) methods that ensure LSDE, significantly reducing the computational cost.

4.1. Full LDSE

The iterative strategy described in Czekaj et al. (2014) performs, for a given set of Galactic fundamental functions, a local normalization and a sphere simulation around the Sun to compute at the position of the Sun: the surface mass density of the generated stars, Σ_⊙, the volume mass density of the stars generated τ years ago, that is ${\rho }_g(\tau ,{\overline{x}}_{\odot })$, and the volume mass density for stars generated τ years ago which are not stellar remnants at present, that is $\rho (\tau ,{\overline{x}}_{\odot })$ ³. We define the stars which are not remnant at present as those with τ ≤ T_lim(M, Z, α), where T_lim(M, Z, α) is the maximum age that a star of a given mass, metallicity, and [α/Fe] reaches without becoming a remnant object. The stellar mass density at the position of the Sun is fitted with the stars that are not remnants at present and the white dwarfs’ density is added separately. The mass lost by the stars during their evolution and the interaction between components of a stellar multiple system are neglected. The total mass in stars, plus that of the interstellar medium, the dark matter, and the central mass, allows us to compute the radial force. At this stage, the dark matter density distribution and the central mass is adjusted such that the model rotation curve fits the observations, the fit is done using the least-squares method in velocity. Finally, the Poisson and the first-order moment of the collisionless Boltzmann equation in the vertical direction are iteratively solved. The whole strategy is iterated from the beginning until convergence is reached.

4.2. Analytical LDSE

In this analytical approach to ensure LDSE, we follow the process described in Sect. 4.1, but instead of using a simulation of a sphere around the Sun, we analytically derive the surface mass density of the generated stars (Σ_⊙), the volume mass density of the stars generated τ years ago (${\rho }_g(\tau ,{\overline{x}}_{\odot })$), and the volume mass density for stars generated τ years ago which are not remnant at present ($\rho (\tau ,{\overline{x}}_{\odot })$). As this approach concerns only the thin-disc component, from now on we avoid the use of the index i.

Using Eq. (31) we can express Σ_⊙ as the sum of the surface density for primal stars $({\mathrm{\Sigma }}_{\odot }^{▵})$ and the surface density for secondary stars $({\mathrm{\Sigma }}_{\odot }^{\$})$. In an equivalent way, the ${\rho }_g(\tau ,{\overline{x}}_{\odot })$ for the generated stars can be expressed as

$$\begin{array}{c}\hfill {\rho }_g(\tau ,{\overline{x}}_{\odot })={\rho }_g^{▵}(\tau ,{\overline{x}}_{\odot })+{\rho }_g^{\$}(\tau ,{\overline{x}}_{\odot }),\end{array}$$(38)

and the stellar volume mass density of the stars with τ ≤ T_lim(M, Z, α) can be expressed as

$$\begin{array}{c}\hfill \rho (\tau ,{\overline{x}}_{\odot })={\rho }^{▵}(\tau ,{\overline{x}}_{\odot })+{\rho }^{\$}(\tau ,{\overline{x}}_{\odot }).\end{array}$$(39)

Now we want to derive the stellar volume mass density of the primal stars at the position of the Sun ${\rho }_g^{▵}(\tau ,Z,{\overline{x}}_{\odot },\alpha )$. Therefore, from equation (32) we obtain

$$\begin{array}{cc}\hfill {\rho }_g^{▵}(\tau ,Z,{\overline{x}}_{\odot },\alpha )& \approx {\displaystyle \underset{\forall M\in {\mathrm{\P }}_r}{\int }\frac{{\mathrm{\Sigma }}_{\odot }^{▵}·{\psi }_{\odot }(\tau )}{\mathcal{H}(\tau )}·\xi (M)·M}\hfill \\ \hfill & ·\mathcal{P}(Z|\tau ,{\overline{x}}_{\odot })·\mathcal{P}(\alpha |\tau ,Z,{\overline{x}}_{\odot })\mathrm{d}M.\hfill \end{array}$$(40)

At this stage we assume that at the position of the Sun all the stars have solar metallicity Z = Z_⊙ and for the moment BGM assumes that the α-elements-to-iron abundance for the thin disc is α = α_⊙. The probabilities for Z and α then become

$$\begin{array}{c}\hfill \mathcal{P}(Z|\tau ,{\overline{x}}_{\odot })\{\begin{array}{cc}1,\hfill & \hspace{0.17em}\text{if}Z=Z_{\odot }\hfill \\ 0,\hfill & \text{otherwise}\hfill \end{array},\end{array}$$(41)

and

$$\begin{array}{c}\hfill \mathcal{P}(\alpha |\tau ,Z,{\overline{x}}_{\odot })\{\begin{array}{cc}1,\hfill & \hspace{0.17em}\text{if}\alpha ={\alpha }_{\odot }\hfill \\ 0,\hfill & \text{otherwise}\hfill \end{array}.\end{array}$$(42)

Combining Eqs. (40)–(42) we can approximate ${\rho }_g^{▵}(\tau ,{\overline{x}}_{\odot })$ by

$$\begin{array}{c}\hfill {\rho }_g^{▵}(\tau ,{\overline{x}}_{\odot })\approx {\displaystyle \underset{\forall M\in {\mathrm{\P }}_r}{\int }\frac{{\mathrm{\Sigma }}_{\odot }^{▵}·{\psi }_{\odot }(\tau )}{\mathcal{H}(\tau )}·\xi (M)·M\mathrm{d}M.}\end{array}$$(43)

As mentioned above, the fit with the observational value of the volume stellar mass density at the position of the Sun is done with the stars with τ ≤ T_lim(M, Z, α). Therefore, as a next step, we need to derive an analytical expression for the density of non-remnant primal stars at the position of the Sun. This is given by

$$\begin{array}{cc}\hfill {\rho }^{▵}(\tau ,Z,{\overline{x}}_{\odot },\alpha )& \approx {\displaystyle \underset{\forall M\in {\mathrm{\P }}_r}{\int }\mathcal{O}(\tau ,M,Z,\alpha )·\frac{{\mathrm{\Sigma }}_{\odot }^{▵}·{\psi }_{\odot }(\tau )}{\mathcal{H}(\tau )}·\xi (M)·M}\hfill \\ \hfill & ·\mathcal{P}(Z|\tau ,{\overline{x}}_{\odot })·\mathcal{P}(\alpha |\tau ,Z,{\overline{x}}_{\odot })\mathrm{d}M,\hfill \end{array}$$(44)

where the function 𝒪(τ, M, Z, α) discussed in Sect. 4.3 is defined as

$$\begin{array}{c}\hfill \mathcal{O}(\tau ,M,Z,\alpha )=\{\begin{array}{cc}0,\hfill & \tau >T_{\lim }(M,Z,\alpha )\hfill \\ 1,\hfill & \tau \le T_{\lim }(M,Z,\alpha )\hfill \end{array}.\end{array}$$(45)

As before, we assume all the stars at the position of the Sun have Z = Z_⊙ and α = α_⊙ (Eqs. (41) and (42)), thus

$$\begin{array}{c}\hfill {\rho }^{▵}(\tau ,{\overline{x}}_{\odot })\approx {\displaystyle \underset{\forall M\in {\mathrm{\P }}_r}{\int }\mathcal{O}(\tau ,M,Z_{\odot },{\alpha }_{\odot })·\frac{{\mathrm{\Sigma }}_{\odot }^{▵}·{\psi }_{\odot }(\tau )}{\mathcal{H}(\tau )}·\xi (M)·M·\mathrm{d}M.}\end{array}$$(46)

Up to now we have derived analytical expressions for the volume densities of the primal stars ((43) and (46)). To complete Eqs. (38) and (39) we need to derive the analytical expressions for the mass density of secondary components. Considering Eq. (33) we can write ${\rho }_g^{\$}(\tau ,Z,{\overline{x}}_{\odot },\alpha )$:

$$\begin{array}{cc}\hfill {\rho }_g^{\$}(\tau ,Z,{\overline{x}}_{\odot },\alpha )& \approx {\displaystyle \underset{\forall M\in {\mathrm{\P }}_r^i}{\int }\underset{m_{min}}{\overset{m_{max}}{\int }}\frac{{\mathrm{\Sigma }}_{\odot }^{▵}·{\psi }_{\odot }(\tau )}{\mathcal{H}(\tau )}·\xi (M)·\mathcal{P}(Z|\tau ,{\overline{x}}_{\odot })}\hfill \\ \hfill & ·\mathcal{P}(\alpha |\tau ,Z,{\overline{x}}_{\odot })·P(\mathrm{bin}|\tau ,M,Z,\alpha )\hfill \\ \hfill & ·P(m|\tau ,M,Z,\alpha )·m\mathrm{d}m\mathrm{d}M.\hfill \end{array}$$(47)

Next, to continue reducing computational time, we assume that the mass of the secondary is only dependent on the mass of the primary component P(m|τ, M, Z, α)≈P(m|M) and we consider this probability is uniform and given by

$$\begin{array}{c}\hfill p(m|M)\approx \{\begin{array}{cc}\frac{1}{M-M_{min}},\hfill & \hspace{0.17em}\text{if}{M}_{min}<m\le M\hfill \\ 0,\hfill & \text{otherwise}\hfill \end{array},\end{array}$$(48)

where M_min is the minimum mass needed to generate a star. This expression states that the mass of the secondary star is always equal or less massive than the corresponding primary component.

We introduce a second approximation assuming that for a given generated primal star the probability of being the primary component of a multiple system is given by

$$\begin{array}{c}\hfill P(\mathrm{bin}|\tau ,M,Z,\alpha )\approx \{\begin{array}{cc}0,\hfill & \tau >T_{\lim }\hfill \\ P(\mathrm{bin}|M),\hfill & 0<\tau \le T_{\lim }\hfill \end{array}.\end{array}$$(49)

This probability, as defined, is not null only for the non-remnant stars. Therefore, using Eq. (45), it can also be expressed as p(bin|τ, M, Z, α)=𝒪(τ, M, Z, α)⋅p(bin|M). In other words, this approach assumes that the probability of being a primary component of a binary system is null for those stars which are remnant at present. Introducing Eqs. (48) and (49) in (47), we have

$$\begin{array}{cc}\hfill {\rho }_g^{\$}(\tau ,Z,{\overline{x}}_{\odot },\alpha )& \approx {\displaystyle \underset{\forall M\in {\mathrm{\P }}_r^i}{\int }\underset{m_{min}}{\overset{m_{max}}{\int }}\mathcal{O}(\tau ,M,Z_{\odot },{\alpha }_{\odot })·\frac{{\mathrm{\Sigma }}_{\odot }^{▵}·{\psi }_{\odot }(\tau )}{\mathcal{H}(\tau )}·\xi (M)}\hfill \\ \hfill & ·\mathcal{P}(Z|\tau ,{\overline{x}}_{\odot })·\mathcal{P}(\alpha |\tau ,Z,{\overline{x}}_{\odot })\hfill \\ \hfill & ·P(\mathrm{bin}|M)·P(m|M)·m\mathrm{d}m\mathrm{d}M.\hfill \end{array}$$(50)

The above assumptions imply that secondary stars are never remnants, thus ${\rho }^{\$}(\tau ,Z,{\overline{x}}_{\odot },\alpha )={\rho }_g^{\$}(\tau ,Z,{\overline{x}}_{\odot },\alpha )$.

As done for primal stars, we assume that all the stars at the position of the Sun have solar metallicity Z = Z_⊙ and α = α_⊙. Therefore,

$$\begin{array}{cc}\hfill {\rho }^{\$}(\tau ,{\overline{x}}_{\odot })& ={\rho }_g^{\$}(\tau ,{\overline{x}}_{\odot })\approx {\displaystyle \underset{M=M_{min}}{\overset{M=M_{max}}{\int }}\underset{m=M_{min}}{\overset{m=M_{max}}{\int }}\mathcal{O}(\tau ,M,Z_{\odot },{\alpha }_{\odot })}\hfill \\ \hfill & ·\frac{{\mathrm{\Sigma }}_{\odot }^{▵}·{\psi }_{\odot }(\tau )}{\mathcal{H}(\tau )}·\xi (M)·P(\mathrm{bin}|M)·P(m|M)·m\mathrm{d}m\mathrm{d}M.\hfill \end{array}$$(51)

As done for primal stars, we leave out ℋ(τ) from (51) to reach the expression for ${\mathrm{\Sigma }}_{\odot }^{\$}$:

$$\begin{array}{cc}\hfill {\mathrm{\Sigma }}_{\odot }^{\$}& \approx {\displaystyle \underset{\tau =0}{\overset{\tau ={\tau }_{\mathrm{end}}}{\int }}\underset{M=M_{min}}{\overset{M=M_{max}}{\int }}\underset{m=M_{min}}{\overset{m=M_{max}}{\int }}\mathcal{O}(\tau ,M,Z_{\odot },{\alpha }_{\odot })·{\mathrm{\Sigma }}_{\odot }^{▵}·{\psi }_{\odot }(\tau )·\xi (M)}\hfill \\ \hfill & ·P(\mathrm{bin}|M)·P(m|M)·m\mathrm{d}m\mathrm{d}M\mathrm{d}\tau .\hfill \end{array}$$(52)

At this point we have on hand all the analytical expressions to compute the Σ_⊙ that fits the observed stellar mass density at the position of the Sun (${\rho }_{\text{obs}}^{\odot }$) for the thin disc stars with τ ≤ T_lim(M, Z, α). To do so we use Eqs. (51) and (46) in (39) and set $\rho (\tau ,{\overline{x}}_{\odot })={\rho }_{\mathrm{obs}}^{\odot }$, with ${\rho }_{\text{obs}}^{\odot }$ being the observed stellar mass volume density at the position of the Sun for stars which are not remnant at present. We then solve the resulting equation for ${\sum }_{\odot }^{\Delta }$ and compute Σ_⊙ from (31) and (52).

For the practical implementation of the analytical approach developed in this section we adopt M_min = 0.09 M_⊙ and M_max = 120 M_⊙, the mass range of the evolutionary models we are using at present (Chabrier & Baraffe 1997; Bertelli et al. 2008, 2009). For p(bin|M) in Eq. (49) we follow the expression for main sequence stars in Arenou (2011).

4.3. The non-remnant fraction (Ω)

The 𝒪 function (45) is directly related with the stellar evolutionary tracks considered through the expressions of T_lim(M, Z, α). We define T_lim(M, Z, α) as the maximum age for which a star of a given mass, metallicity, and α-elements-to-iron abundance is still not a stellar remnant. As discussed in the previous section, for the 𝒪 function we consider solar metallicity Z = Z_⊙ and solar α = α_⊙. We use, for the moment, the stellar evolutionary tracks Bertelli et al. (2008) and Chabrier & Baraffe (1997) considered in Czekaj et al. (2014) which do not consider the α-element abundance; therefore T_lim(M, Z, α)=T_lim(M, Z). From the mentioned stellar evolutionary tracks we derive T_lim(M, Z_⊙) fitting three truncated logarithmic expressions to the two-dimensional (2D) grid of mass and age limit for Z = Z_⊙. We end up with the following expressions

for M/M_⊙ ≥ 7.

$$\begin{array}{c}\hfill T_{\lim }(M,Z_{\odot })\approx \exp (-1.6·\ln (M)+20.8)\end{array}$$(53)

for 2.2 < M/M_⊙ < 7.0,

$$\begin{array}{c}\hfill T_{\lim }(M,Z_{\odot })\approx \exp (-2.7·\ln (M)+23.0)\end{array}$$(54)

for 2.0 < M/M_⊙ < 2.2,

$$\begin{array}{c}\hfill T_{\lim }(M,Z_{\odot })\approx \exp (-2.7·\ln (2.2)+23.0)\end{array}$$(55)

for M/M_⊙ ≤ 2.2,

$$\begin{array}{c}\hfill T_{\lim }(M,Z_{\odot })\approx \exp (-3.5·\ln (M)+23.3).\end{array}$$(56)

With these latter expressions we have everything we need for the execution of the analytical LDSE as described in Sect. 4.2. At this point, to take into account that the BGM Std strategy samples the age of the stars uniformly inside each age sub-population (see Sect. 2.2) we substitute the 𝒪 function by the Ω function in the equations of Sect. 4.2. The Ω(τ, M, Z_⊙) function is defined as

$$\begin{array}{c}\hfill \mathrm{\Omega }(\tau ,M,Z_{\odot })=\{\begin{array}{cc}1,\hfill & T_{\lim }(M,Z_{\odot })>T_{\mathrm{end}}^j(\tau )\hfill \\ \frac{T_{\lim }(M,Z_{\odot })-T_{\mathrm{ini}}^j(\tau )}{T_{\mathrm{end}}^j(\tau )-T_{\mathrm{ini}}^j(\tau )},\hfill & T_{\mathrm{ini}}^j(\tau )\le T_{\lim }(M,Z_{\odot })\le T_{\mathrm{end}}^j(\tau )\hfill \\ 0,\hfill & T_{\lim }(M,Z_{\odot })<T_{\mathrm{ini}}^j(\tau ),\hfill \end{array}\end{array}$$(57)

where $T_{\text{ini}}^j(\tau )\hspace{0.17em}\text{and}\hspace{0.17em}{T}_{\text{end}}^j(\tau )$ are respectively the low and high age boundaries for each of the seven age sub-populations of the thin disc (j = 1…7). We note that $T_{\text{ini}}^j(\tau )\hspace{0.17em}\text{and}\hspace{0.17em}{T}_{\text{end}}^j(\tau )$ obviously depend on τ, because the age τ of the star establishes the j-sub-population to which the star belongs. We emphasise that the Ω function defined above gives the fraction of stars which are not remnant at present and is constant with τ inside the limits of each age sub-population. Furthermore, when defining 𝒪 and Ω functions we are neglecting both the entangled evolution of stars belonging to the multiple stellar systems and the mass lost by the star during its evolution. If in the future we want to implement the effect of the mass lost by the stars during its evolution, we can do this by introducing a new function, ℒ(τ, M, Z, α), into the equations, accounting for the fraction of mass that a given star has lost during its evolution up to the present time. If at some point we are interested in accounting for the entangled evolution of a stellar multiple system then we need to implement the desired evolutionary model of stellar systems directly in the BGM Std strategy, generate a Mother Simulation, and finally modify T_lim accordingly.

4.4. Approximate LDSE

In Sects. 4.2 and 4.3 we explained how, in the full process to ensure the LDSE (Sect. 4.1), we can substitute the entire simulation of a sphere around the Sun by crafted analytical expressions which are computationally undemanding. Here we complete the construction of an approximate LDSE strategy by complementing Sects. 4.2 and 4.3 with a final assumption to make the process even faster. As described in Sect. 4.1 in the full LDSE process, the central mass and the dark matter density distribution is adjusted such that the model rotation curve fits the observations. We can avoid the computational cost of this adjustment by assuming that the central mass density and the dark matter density distribution are invariant under variations of the Galactic fundamental functions⁴.

We now test the performance of the approximate LDSE to evaluate the assumptions that we made and to constrain its range of validity. In the following section we present a set of tests comparing the results obtained when applying the full LDSE process against the results of our approximate LDSE.

4.5. Validation of the approximate LDSE

We develop three tests to evaluate the goodness of our approximate method (Sects. 4.2–4.4). These tests quantify the differences, when using the approximate instead of full method, in some of the key parameters resulting from the LDSE process. All the tests are performed using the seven model variants described in Mor et al. (2017) that were built with different assumptions of the IMF, the SFH, and the density laws. In Table 1 we present their main parameters. These model variants cover the range of parameters that we want to explore in this paper well, and therefore they are a good set to perform the tests. The parameters for the thin disc density laws which are not listed in Table 1 are adopted to be the same as used in Mor et al. (2017), being the functional forms of the density laws listed in Robin et al. (2003). The rest of the model ingredients which are not specifically indicated in Table 1, for example the atmosphere models, the stellar evolutionary tracks, the age-metallicity relation, and the age-velocity dispersion, are adopted to be the same as in “Model B” listed in Table 5 of Czekaj et al. (2014).

In the first test we show that our approximate method obtains a rotation curve compatible with the one resulting from the full method. For all model variants, we obtain differences of less than 2% in the rotational velocity between a galactocentric radius from 3 kpc to 14 kpc. Moreover these differences are much smaller than discrepancies between the observational values (e.g. Caldwell & Ostriker 1981 and Sofue 2015). In Fig. 1 we present the results for the two model variants where we found the highest discrepancies between the rotation curve obtained from both methods. The discrepancies along the curve from 3 kpc to 14 kpc are always smaller than 5 km s⁻¹. Thus we consider that the approximated LDSE is valid within these galactocentric radii.

Fig. 1. Two examples of the comparison between the rotation curve obtained with both the full process to ensure LDSE and the rotation curve obtained with the approximate LDSE. Top panel: results for the DAV variant from Mor et al. (2017). Bottom panel: results for the HRBV variant from Mor et al. (2017). The green triangles show the data points derived from the Caldwell & Ostriker (1981) rotation curve assuming R0 = 8000 pc and V0 = 230 km s−1 for the Sun. The magenta triangles show the data points of the rotation curve from Sofue (2015). For Sofue data, error bars are provided by the author, while we have estimated the errors of the data from Caldwell & Ostriker (1981) following the expressions and tables in their paper.

The second test compares for each model variant the eccentricities of the Einasto density profiles which are fitted inside the process that ensures the LDSE. The Einasto eccentricities are computed for each one of the seven age sub-populations of the thin disc. In all cases we find differences smaller than 0.7%.

Our third test compares the local volume stellar mass density of all the thin-disc age sub-populations. For each model variant we check our capability to recover the density values of Table 1 of Mor et al. (2017). The discrepancies are within the error bars and are always smaller than 5%.

We want to point out that these tests are the first empirical demonstration that our strategy for the BGM FASt simulation is well founded. Whereas it is true that they validate the equations behind the weight function, they do not evaluate the capabilities of the weight function itself to generate a BGM FASt simulation. To do so we need to directly test distributions of observable parameters obtained with BGM FASt simulations. Tests on age and mass distribution are also needed to deeper validate the BGM FASt strategy. These tests are presented in Sect. 6 and Appendix A.1.

5. ABC to infer Galactic fundamental functions

5.1. The ABC method

For very complex models such as the BGM, to obtain the exact likelihood function is mathematically impossible or computationally prohibitive. In these cases the ABC algorithms allow us to compute an approximate posterior PDF of the explored parameters. In this paper we use a sequential Monte Carlo approximate Bayesian computation algorithm (SMC-ABC) because of its balance between high acceptance rate and independence of the outcomes. The algorithm that we use is further described in Jennings & Madigan (2017), and is available as a python package named astroABC. The SMC-ABC method is an improvement of basic ABC acceptance-rejection sampling (ARS-ABC) to optimise the acceptance rate. The basic ARS-ABC contains the theoretical bases of the SMC-ABC algorithm. This latter first generates a proposed set of parameters θ̄ from the prior PDF. Subsequently, the algorithm generates the simulated data using the proposed set of parameters in a given model, and finally the set of parameters is accepted as part of the posterior PDF if the simulated data (𝒟_simu) are equal to the real data (𝒟); otherwise θ̄ is rejected and we start the process again. This algorithm is very restrictive because it is almost impossible to find a model with a given set of parameters that perfectly reproduces the data. Usually one uses the ARS-ABC, relaxing the condition 𝒟_simu = 𝒟 and approximating the method as follows

• (1)

  Generate θ̄ from the prior PDF.

• (2)

  Simulate data 𝒟_sim from the model ℳ with parameters θ̄

• (3)

  Calculate the distance δ(𝒟, 𝒟_sim) between 𝒟 and 𝒟_sim.

• (4)

  Accept θ̄ if δ is smaller than a given threshold (υ) δ ≤ υ; return to 1.

This approximate algorithm needs to adopt an adequate distance metric δ and a threshold υ. When υ → 0, the algorithm is sampling exactly from the posterior PDF. Very small values of υ diminish the acceptance rate of the algorithm, and therefore the choice of υ must be a compromise between computational cost and accuracy of the posterior PDF. We want to emphasise that ARS-ABC samplers generate independent outcomes; each iteration is independent from the previous.

Other ABC methods, like the so-called free likelihood MCMC samplers (MCMC-ABC), are also built to optimise the acceptance rate. The MCMC-ABC increases the acceptance rate by modifying MCMC sampling algorithms (such as Metropolis-Hastings) to be able to work without the need of the likelihood (Marjoram et al. 2003), but the price paid in this case is that the outcomes are dependent (Marjoram et al. 2003). In a different way, the SMC-ABC improves the acceptance rate with a double entangled optimisation. For each iteration, the algorithm uses kernels to assign a higher sampling probability to the sets of parameters with better results. The resulting new sampling probability is used in the next iteration. In other words, regions of the parameter space with better results are visited more frequently. The application of kernels in the limits of the prior PDF can sometimes produce a posterior PDF slightly wider than the limits stated by the prior PDF. The use of the kernels is complemented with an adaptive threshold υ, where the upper and lower limits have to be set. This allows to go from a more relaxed threshold in the first iterations down to a small threshold for the last ones, optimising the sampling of the posterior PDF in terms of computational time. We emphasise that for a given expression for the distance metric and a given υ, the outcomes are still independent in the SMC-ABC.

Now we need to choose the data in accordance with the parameters that we want to infer. Observations can only provide a small subset of all the potential data defining the Milky Way system. This limitation forces us to search for summary statistics, which are of a lower dimension and are incomplete. If the chosen summary statistics, 𝒮, are statistically sufficient for 𝒟 then the posterior PDF under the summary data P(θ|𝒮) is equivalent to the posterior PDF under the full data P(θ|𝒟). In practice it is very difficult, or impossible, to formally identify rigorous summary statistics sufficient for 𝒟. As suggested in Marjoram et al. (2003) we use a more heuristic approach for our specific problem. In the following section we propose summary statistics 𝒮 that capture information on the θ̄ essential parameters of the Milky Way. Once sufficient statistics are established, 𝒟 is replaced by 𝒮 in the algorithm as done in Marjoram et al. (2003). The theoretical basis for these algorithms can be found in Marin et al. (2011), Beaumont et al. (2009) and Sisson & Fan (2010), for example.

5.2. Bayesian inference in the solar neighbourhood

For the appropriate use of the ABC algorithm described above, we need to define, for each of the specific scientific goals that we want to achieve, our sufficient statistics 𝒮, a distance metric δ, and a threshold υ. For both the evaluation of the BGM FASt performance (Sect. 6) and the science demonstration cases (Sect. 7), we consider as our sufficient statistics 𝒮 the star counts in a binned four-dimensional space of position, apparent magnitude, and observed colour. This means that for the purposes of this paper we define our sufficient statistics as the number of stars in each bin of latitude, longitude, visual Tycho magnitude (V_T), and Tycho-2 colour (B − V)_T. Specifically, our 𝒮 is the colour-magnitude diagram split in three latitude ranges: (|b|< 10), (10 < |b|< 30), and (30 < |b|< 90). We choose the bin size of the colour-magnitude diagram to be large enough to allow for a robust statistical analysis to be performed, but small enough to avoid loosing information. As mentioned in the previous section, our choice of 𝒮 is based on our previous experience when comparing BGM with observational data. We know that the observed colour-magnitude diagrams in different regions of the sky offer very valuable information to constrain the IMF, the SFH, the local stellar mass density, and the density laws (e.g. Mor et al. 2017; Robin et al. 2014, 2012). In some of these papers these colour-magnitude diagrams have already been used as sufficient statistics combined with ABC algorithms (e.g. Robin et al. 2014). Furthermore, in Czekaj et al. (2014) it was demonstrated that the division of the sky into the three above-mentioned latitude ranges is very useful to analyse the Galactic fundamental functions by fitting BGM to Tycho-2 data. In our star counts analysis we always work inside the completeness limits of the observational catalogues and therefore our sufficient statistics when using Tycho-2 data are limited at V_T = 11 where Tycho-2 is complete up to 99%.

Once our sufficient statistics are defined, we need to choose the distance metric to quantify differences between the observed and simulated data, δ(𝒮_obs,𝒮_simu). We use this to compare the simulated 𝒮_simu and the observed 𝒮_obs data. For this paper we use the following expression, which we call Poissonian distance:

$$\begin{array}{c}\hfill {\delta }_{\mathrm{P}}({\mathcal{S}}_{\mathrm{obs}},{\mathcal{S}}_{\mathrm{simu}})=\left|{\displaystyle \sum _{i=1}^N{q}_i·(1-R_i+\ln (R_i))}\right|,\end{array}$$(58)

where R_i is defined as the quotient R_i = f_i/q_i and q_i and f_i are the number of stars in the data and the model, respectively. Additionally, we penalise the cases where there are no stars observed in the bin but the model predicts stars to be present by replacing q_i and f_i in the formula with q_i + 1 and f_i + 1, respectively. The defined distance metric becomes zero when the simulation and the observations have the same number of stars in each bin. The smaller the value of the distance metric, the closer 𝒮_simu is to 𝒮_obs. From Kendall & Stuart (1973) and Bienaymé et al. (1987) we know this formula is a good choice for the comparison between observed and simulated colour-magnitude diagrams in terms of star counts. Furthermore, in Mor et al. (2017) we already introduced the idea that expression (58) can be understood as a distance metric.

Finally, we choose an upper limit and a lower limit for the threshold υ by experimenting with different values to have a good balance between computational cost and accuracy of the posterior PDF.

6. Evaluating BGM FASt at the solar neighbourhood

In this section we evaluate the behaviour of BGM FASt in the solar neighbourhood and we demonstrate that the performance of the BGM FASt strategy is independent of the Mother Simulation. In Sect. 6.1 we use both BGM FASt and BGM Std to simulate the solar neighbourhood. We then analyse the resulting samples limited in apparent magnitude by comparing colour distributions and CMDs. Additionally, we extend the comparison to mass and age distributions for a deeper analysis of the BGM FASt performance. We do these comparisons in the completeness regime of Tycho-2 catalogue, and therefore we use samples limited in Tycho visual apparent magnitude (V_T) up to 11, where Tycho-2 is complete at 99%.

In Sect. 6.2 we present a test for the BGM FASt framework together with the ABC algorithm, consisting in exploring the parameter space of the SFH, demonstrating that we are able to correctly recover an imposed SFH in the solar neighbourhood.

6.1. BGM FASt versus BGM Std

To evaluate the performance of the BGM FASt we compare simulations generated with the same sets of fundamental functions obtained from both the BGM FASt and BGM Std strategy. We have selected five model variants from Mor et al. (2017): DAV, DBV, DCV, HRV and SV. All of them constitute a good framework to analyse the BGM FASt behaviour in the solar neighbourhood as they adequately cover the parameter space that we want to explore with BGM FASt in Sect. 7. For each model variant we perform the tests using both the Drimmel & Spergel (2001) and Marshall et al. (2006) extinction maps⁵. The main set of comparisons aims to evaluate the behaviour of BGM FASt when using as Mother Simulation the best fit variant from Mor et al. (2017) (the DAV variant) obtained using Galactic Cepheids and Tycho-2 data. The DAV variant is our choice for the Mother Simulation when exploring a six-dimensional (6D) parameter space in Sect. 7.2. In the tests we evaluate the effects of changing one or more parameters of the fundamental functions when generating BGM FASt simulations according to the model variants DBV, DCV, HRV and SV.

Additionally we repeat the comparisons but use the DCV variant (constant SFH) as the Mother Simulation. These additional tests help us to evaluate the BGM FASt performance in several aspects. First, they allow us to demonstrate whether or not the performance of the BGM FASt strategy is independent of the Mother Simulation. Second, we are able to analyse possible dependencies on the total number of stars of the Mother Simulation. Third, they allow us to study simultaneous changes of the IMF and the SFH. Finally, they provide us with more data, which can be used to better interpret the obtained results.

In Table 2 we summarise the parameters of the IMF, the SFH, and the density laws for the variants involved in the tests. The functional form for the SFH is assumed to be an exponential law, the IMF is assumed to be a three-times-truncated power law and the density profiles are assumed to be of Einasto shape. In the first column we show the parameters that we use for the Mother Simulation: the DAV and DCV. In the second column we show the parameters for the BGM FASt simulations. The same parameters are also used to perform BGM Std simulations. We mark in bold text the parameters that have been modified from the Mother Simulation to generate the BGM FASt simulation. In the third column (FASt vs. Std) we present a summary of the global results comparing BGM FASt with BGM Std simulations, δ_p is for the Poissonian distance from Eq. (58) and % refers to the discrepancies of the total number of stars between BGM FASt and BGM Std simulations, that is ((# Stars in BGM FASt – # Stars in BGM Std)/# Stars in BGM Std) 100. We note that for all cases except the SV variant the discrepancies in the total number of stars is smaller than 4%. When looking to the Poissonian distance metric we note that all the values are below 2000, except, again, for the SV variant. The Poissonian distance between BGM FASt and BGM Std simulations is one order of magnitude smaller than the difference between the best fit variant from Mor et al. (2017) and Tycho-2 data.