© The Authors 2025

1 Introduction

The star formation history (SFH) of the Milky Way disc is one of the fundamental parameters that characterise the evolution of our galaxy, accounting for the internal and external processes responsible for the variation in the star formation rate (SFR). Entangled with the SFH (Haywood et al. 1997; Aumer & Binney 2009; Kroupa & Jerabkova 2021), the initial mass function (IMF; Salpeter 1955) plays a crucial role in the chemical evolution of the Milky Way (and stellar systems in general; e.g. Burbidge et al. 1957; Tinsley 1980; Pagel 2009). Therefore, since the middle of the 20th century, the Galactic research community has focused enormous efforts on the precise determination of both the SFH (e.g. Cignoni et al. 2007; Snaith et al. 2015; Bernard 2018; Haywood et al. 2018; Ruiz-Lara et al. 2020; Sysoliatina & Just 2021; Cukanovaite et al. 2023; Mazzi et al. 2024; Gallart et al. 2024; Spitoni et al. 2024) and the IMF (e.g. Salpeter 1955; Schmidt 1963; Scalo 1998; Kroupa 2001; Chabrier 2003; Sollima 2019; Haghi et al. 2020; Li et al. 2023; Dickson et al. 2023; Kirkpatrick et al. 2024).

A wide range of tools and techniques have been implemented to determine these fundamental stellar-population parameters of the Galaxy. Most of them rely in some way on star counts as a function of magnitude (e.g. Bahcall & Soneira 1980) or as a function of magnitude and colour (Robin et al. 2003; Girardi et al. 2005). With the advent of HIPPARCOS (European Space Agency 1997; Perryman et al. 1995) and Gaia (Gaia Collaboration 2016), an additional valuable piece of information could be added: their parallax measurements allowed for the study of absolute rather than apparent colour–magnitude-diagrams (CMDs; e.g. Vergely et al. 2002) and for the construction of volume-limited samples (e.g. Kirkpatrick et al. 2024). However, none of the previous studies have simultaneously derived both the SFH and the IMF.

Within the different Galactic models, the Besançon Galaxy model (BGM; Robin et al. 2003) is a holistic population synthesis approach that has been widely used in recent decades for the statistical study of the formation, evolution, and present content of the Milky Way. From its new strategy (Czekaj et al. 2014), it became possible to directly use the IMF and the SFH to generate a full-sky mock catalogue to be compared with data from Tycho-2, constituting the starting point of the Galactic parameters inference with BGM. Later, Lagarde et al. (2017) presented a new version of BGM incorporating the STAREVOL stellar evolutionary tracks.

At that point, the execution of a full-sky catalogue of simulated stars using the standard BGM (hereafter BGM Std) implied a computational cost that could not be assumed for the inference of Galactic parameters using iterative methods. This process demanded the use of statistically robust and physically adequate approximations. To respond to these requirements, Mor et al. (2018) presented the BGM fast approximate simulations (BGM FASt), a computationally cheap approximation of BGM that lets us compute fast BGM pseudo-simulations (PSs) weighting the stars of a BGM Std simulation. Using this new tool and the approximate Bayesian computation (ABC) technique, Mor et al. (2019) made a first derivation of the IMF and the SFH of the Galactic thin disc (among other parameters) using Gaia DR2 data up to G = 12.

In parallel, the astrometry outcomes of the Gaia second data release (Gaia Collaboration 2018) opened the window for the improvement of the robustness of the BGM modelling. The problem of dynamical self-consistency was tackled in Robin et al. (2022) by fitting the gravitational potential of the Milky Way to the stellar kinematics and densities of (mostly solar neighbourhood) Gaia DR2 stars down to magnitude G = 17. At that point, it became clear that it is necessary to re-determine the SFH and IMF of the Milky Way in the solar neighbourhood.

We present in this work a new implementation of BGM FASt and, for the first time, its execution considering the full-sky Gaia DR3 data up to magnitude G = 13. The high computational demands of the BGM FASt code have now been overcome. The optimised code presented in this work enables future executions at fainter limiting magnitudes.

This article is structured as follows. We present in Sect. 2 a basic background on the BGM FASt framework and its capabilities, which is needed to understand the upgrades in the method presented in Sect. 3. Then, in Sect. 4 we show in detail the results of this new implementation, which are contextualised and related to the literature in Sect. 5. Finally, we present the conclusions of this work and some future steps to take in Sect. 6.

2 The BGM FASt framework for Galaxy modelling

The BGM FASt framework lets us infer and evaluate Galactic parameters by iteratively fitting our model to observed data. It performs computationally cheap PSs of the BGM, weighting the stars of a pre-sampled BGM Std simulation (the so-called mother simulation, labelled MS; Mor et al. 2018). BGM FASt is ≈ 10⁵ times faster than BGM Std, which opens a window for evaluating and deriving Galaxy modelling ingredients using Bayesian tools. For illustration, the BGM FASt framework is summarised in Fig. 1, while the details of each of its components are given in the following subsections. Following the scheme illustrated in this figure, a BGM FASt PS requires two main inputs: a MS, detailed in Sect. 2.1, and a set of parameters to infer $({\overline{\theta }}_k)$, whose practical implementation for the derivation of the IMF and the SFH is explained in Sect. 3.2.

Fig. 1 Flux diagram of the process involving BGM Std, BGM FASt, and the CMD fitting technique. BGM Std is summarised in Sect. 2.1; we present BGM FASt in Sect. 2.2, where the meaning of the subscripts and superscripts i (indicating the ith Galactic component, e.g. thin disc, thick disc, halo, and bar) and j (associated with the jth interval of the reduced parameter space) are explained in detail; regarding the fitting process, see in Sect. 2.3 the different components of the process as well as an extended explanation of the subscripts and superscripts k (referred to kth iteration within an ABC step), n (indicating the nth bin of the CMD), and t (corresponding to the tth ABC step). See Mor et al. (2018) for a detailed explanation of the equations in the central panel.

2.1 The Besançon Galaxy model

The MS is the result of applying the population synthesis model BGM Std (Robin et al. 2022). This model operates as follows: first, it uses fundamental functions such as the space density distributions, the IMF, the SFH, and the age-metallicity relation of each Galactic component (e.g. thin disc, thick disc, halo, and bar) to generate stars at a given distance with given masses, M, ages, τ, and metallicities, Z, in each volume element (δV in Fig. 1). This process has carried on maintaining dynamical self-consistency thanks to the last update of BGM Std (Robin et al. 2022), explained below. For binary stars, from the available gas at a given position, BGM Std generates a secondary star following a given proportion, binary fraction, semi-major axis, eccentricity, and distribution of the mass ratio dependent on the mass of the primary (Arenou 2011).

Then, BGM Std determines the evolutionary stage of each star and, if the age is shorter than the theoretical lifetime, it computes the position of the star in the Hertzsprung–Russell diagram, from which its astrophysical parameters (effective temperature, T_eff , gravity, log g, and bolometric magnitude, M_Bol) are derived by interpolating in the set of evolutionary tracks. For stars with ages larger than the theoretical lifetime at a given mass, BGM Std computes the mass subsisting as a remnant and the mass released into the interstellar medium, which is instantaneously recycled.

Finally, as is shown in Fig. 1, BGM Std converts the astro- physical parameters of the stars into the observed ones following two steps: first, it applies atmosphere models to get the absolute magnitude and the intrinsic colours; and then, from distance and position it includes the effects of reddening of the interstellar medium by means of extinction maps to obtain the apparent magnitude and the observed colours, which constitute the observables that can be compared with the observed catalogue data after adapting the magnitudes to the corresponding filters. In addition, BGM Std sets an angular separation limit and merges the different components of the system when they are not resolved. In this work, this value is set to 400 mas to mimic Gaia DR3 data (Fabricius et al. 2021).

As was stated before, the star generation in BGM Std uses fundamental functions that depend on each Galactic component. In the inference process made by BGM FASt (see Sect. 2.2), there are two critical global uncertainties of the current best BGM Std model that must be taken into account. On the one hand, the impact of using different stellar evolutionary models (SEMs) on deriving Galactic parameters; on the other hand, how to maintain the dynamical self-consistency of the model with the new set of inferred parameters.

Regarding the choice of stellar models, at present two different SEMs are implemented in BGM Std: PARSEC v1.2 (Bressan et al. 2012; Chen et al. 2014, 2015; Fu et al. 2018) and STAREVOL (Lagarde et al. 2017, 2019) stellar tracks¹. Both stellar models have their specificities (such as the list of nuclear reaction rates or the equations of state) that the reader can find in the relevant references. A summary of them is presented in Appendix A. To compute the Gaia photometry, we make use of the BTSettl models as provided by Allard (2013) available from IVOA².

The spatial distribution of the stars and their kinematics are set using a dynamical self-consistent solution as described in Robin et al. (2022). It is achieved by iteratively fitting an axisymmetric potential and the stellar distribution functions of the model to the Gaia EDR3³ proper motions and parallaxes (see Fig. 2 in Robin et al. 2022 for a better comprehension of the scheme behind the method and Bienaymé et al. 2018 to understand the physical and mathematical basis of the construction of the distribution functions and potential). Consequently, the dynamical self-consistency of a given simulation is dependent on the input parameters of BGM Std. That is why the modification of BGM Std model ingredients such as the IMF, the SFH, or the SEMs will need from a redetermination of the dynamical self-consistency (see Sect. 3.2).

2.2 BGM fast approximate simulations

The complexity of the functions in BGM Std, which vary between Galactic components, requires further exploration to improve the model. BGM FASt addresses this need by defining an N-dimensional space of the parameters we want to infer and a tool to perform fast BGM Std simulations to make it possible (see Mor et al. 2018 for details). We explain here the method for the thin disc, but it can be applied following the same steps to any other Galactic component (e.g. thick disc, halo, and bar). This process is generalised in Fig. 1, where the i^th Galactic component (in this case, the thin disc) is treated as an additional index.

We define the parameter space of the thin disc as the N-dimensional space containing all the parameters involved in the distribution function of the generated stars of this Galactic component. Mathematically, the parameter space is defined as $\P =\tau \times M\times Z\times \overline{x}\times \overline{v}\times \alpha \times \overline{p^{\prime }}$, where $\tau ,M,Z,\overline{x},\overline{v}$, and α are the age, the initial mass, the metallicity, the position, the velocity, and the [α/Fe] content of the stars of the thin disc, respectively, and ${\overline{p}}^{\prime }$ accounts for other possible set of independent parameters⁴. From the definition of ¶, we consider a general distribution function of the generated stars $𝒟(\tau ,M,Z,\overline{x},\overline{v},\alpha ,\overline{p^{\prime }})$. By marginalising the general distribution function over ${\overline{p}}^{\prime }$, we obtain the distribution function for the generated thin disc stars in the reduced space 𝒟_r, denoted as $𝒢(\tau ,M,Z,\overline{x},\overline{v},\alpha )$

BGM FASt is based on the idea that the number of stars generated in the j^th interval of the reduced parameter space ¶_r, $\Delta {\P }_r^j$, is proportional to the mass dedicated to generating stars in that interval (Mor et al. 2018). Therefore, it can be computed as the first moment of the mass:

$N_j=K_j{M}_j=K_j{\displaystyle {\int }_{\Delta {\P }_r^j}𝒢}(\tau ,M,Z,\overline{x},\overline{v},\alpha )\cdot Md{\P }_{\text{r}},$(1)

where K_j is the particular coefficient of the relation for the $\Delta {\P }_r^j$ interval. Under the assumption that each star corresponds to the smallest possible interval of the parameter space, we associate the j^th interval with each individual star. The number of stars in the PS is then calculated by weighting each star by the proportion between the amount of mass dedicated to generating stars in that interval in the PS and the MS:

$\frac{N_j^{PS}}{N_j^{MS}}=\frac{M_j^{PS}}{M_j^{MS}}\Rightarrow N_j^{PS}=\frac{M_j^{PS}}{M_j^{MS}}{N}_j^{MS}={\omega }_j{N}_j^{MS},$(2)

where the weight, ω _j, can be expressed in terms of the distribution function as

${\omega }_j=\frac{{\displaystyle {\int }_{\Delta {\P }_r^j}{𝒢}^{PS}}(\tau ,M,Z,\overline{x},\overline{v},\alpha )\cdot Md{\P }_{\text{r}}}{{\displaystyle {\int }_{\Delta {\P }_r^j}{𝒢}^{MS}}(\tau ,M,Z,\overline{x},\overline{v},\alpha )\cdot Md{\P }_{\text{r}}}.$(3)

As was previously mentioned, Eq. (3), which forms the core of BGM FASt, can be applied to stars from different Galactic components. In Sect. 3.2, we focus our efforts on applying this equation first to thin disc stars and later to thin and thick disc stars to determine their SFHs and IMF.

The rigorous mathematical development of the concept as well as the practical implementation shown in Fig. 1 (how to convert $𝒢(\tau ,M,Z,\overline{x},\overline{v},\alpha )$ into the magnitudes we want to infer, in this case the SFH and the IMF) are explained in detail in Mor et al. (2018). We note that in Fig. 1 there are additional subscripts and superscripts in the magnitudes presented in this section that are associated with the iterative process explained in the next section.

BGM FASt concludes when the PS is generated by applying the weights from Eq. (3) to the stars in the MS. At this point, the fitting process begins, as is shown in Fig. 1, using a summary statistics and a distance metric to quantify the similarity between the PS and the observed data.

For the first time, we have introduced in this work a new computational strategy that reduces the time required to compute stellar weights by a factor of ≈ 5–7, significantly improving the efficiency of the inference process compared to previous BGM FASt publications (Mor et al. 2018, 2019). Taking into account the new implementation of BGM FASt described in Sect. 3.2 and summarised in Eq. (5), the weight of a star only depends on its mass and the interval of ages in which it falls. Given the large number of stars in the Gaia DR3 catalogue (over seven million with G < 13, and a similar number in the MS), many of them will have similar ages and masses. Therefore, we can build a two-dimensional grid of masses and ages, create an artificial star with the mean properties of mass and age in each cell, and assign the weight of that star to all the stars that fall in the same cell. Following this process we reduced the number of weights to compute from millions to hundreds of thousands.

Considering steps in mass and age sufficiently small, the influence of this approximation on the final result is negligible. We tested the new method with 100 random combinations of the parameters (IMF and SFH) and found that the mean relative difference between the observed catalogue-PS distance (see Sect. 2.3) computed with the original version of BGM FASt and the new computationally cheap approximation is below 0.03%, with the maximum relative difference being of 0.1%.

2.3 BGM FASt fitting strategy

The fitting process is based on fine-tuning the Galactic parameters to make the final PS as similar as possible to the observed catalogue data. Thanks to the large number of observed and physical parameters present in the BGM FASt PS, this can be done using different magnitudes from the observable space (e.g. apparent magnitudes, colours, parallaxes…) or derived physical stellar parameters (e.g. masses or ages if we have asteroseis- mologic observations). Under any of the previous options, there are several elements to define: 1) the summary statistics, 2) the distance metric, and 3) a tool to perform the inference.

Summary statistics

The chosen summary statistics must contain as much information as possible to characterise the parameters we want to infer. When the summary statistics is statistically sufficient to describe a system, then the posterior probability distribution function (PDF) of the parameters under the summarised data is equivalent to the posterior PDF under the full data and the summary statistics is called sufficient statistics (Mor et al. 2018). Since this is practically impossible, we used a summary statistics containing a reduced sample of the data available in the BGM FASt PS and the observed catalogue including relevant information for the inference of our parameters.

We took as a summary statistics a binned modified colourmagnitude diagram (CMD) $M_G^{\prime }$ versus Gaia colour (from now on CMD), where $M_G^{\prime }=G+5{\log }_{10}(\varpi /1000)+5$ is the modifiedabsolute magnitude in which we take into account the broad- band G filter and the Gaia parallax ϖ (in mas). We note that the computed $M_G^{\prime }$ differs from the intrinsic M_G because the former is not corrected for the extinction of the interstellar medium. The assumption behind the use of $M_G^{\prime }$ is that the extinction map in the MS is good enough to characterise the reddening suffered by the light captured by Gaia. Furthermore, to improve the statistical resemblance between Gaia and the MS, we applied errors to the astrometric and photometric magnitudes of each star generated with BGM Std using the indications of the Gaia- DPAC consortium for Gaia DR3 (see the Gaia instrument model webpage).

We want to emphasise that the BGM FASt framework gives us the flexibility to fit any type of data using different summary statistics, from full-sky Gaia samples up to a given limiting magnitude (as has been done in this work, see Sect. 3), to particular regions of the sky or the CMD, always taking into account the role of completeness.

Distance metric

Following the initial idea of Kendall & Stuart (1973) and Bienayme et al. (1987), already applied in the context of BGM FASt (Mor et al. 2017, 2018), we used a Poissonian distance, δ_P, to quantify how similar the PS (S_sim) and the observed data (S_obs) are:

${\delta }_P(S_{obs},S_{sim})=\left|{\displaystyle \sum _{n=1}^{N_{\text{bins}\ge 100\text{stars}}}{q}_n}[1-R_n+\ln (R_n)]\right|,$(4)

where R_n = f_n/q_n is defined as the quotient between the number of stars in the n^th bin of the CMD in the model ( f_n) and the observed catalogue (q_n). When the PS fits perfectly the data, R_n = 1 ∀n and δ_P(𝒮 _obs, 𝒮 _sim) = 0. Unlike Mor et al. (2018), we increased the robustness of the fitting process considering only in the computation of the distance those bins with ≥ 100 stars in the catalogue, which allowed us to ensure the statistical significance and work with CMDs cleaned from outliers.

Approximate Bayesian computation

We used a Python implementation of a sequential Monte Carlo ABC algorithm (SMC- ABC; Jennings & Madigan 2017) to carry out the iterative process responsible for fitting the model to the observed data to estimate the approximate posterior PDF of the parameters ($\overline{\rho }(\overline{\theta })$ in Fig. 1). The fitting process proceeded as follows. Once we had a PS generated with a given set of the parameters to infer, ABC compared the distance δ^k between the PS and the observed catalogue with the threshold of the t^th ABC step. If it was lower, the set of proposed parameters ${\overline{\theta }}_k$ was accepted and became part of the prior PDF of the next step, ${\overline{\rho }}_{t+1}(\overline{\theta })$. Otherwise, it was rejected.

As is shown in Fig. 1, the process was carried out in two loops. On the one hand, the distance threshold started at a determined upper limit that diminished at each step, t, until the established lower limit was achieved or the imposed maximum number of steps was reached. On the other hand, at each step the incorporation of sets of parameters $\overline{\theta }$ to ${\overline{\rho }}_{t+1}(\overline{\theta })$ was repeated until reaching the number of accepted simulations per step set by the user to ensure a statistically robust PDF. While the distance threshold of the step was above the lower limit or the maximum number of steps had not been reached, the derived ${\overline{\rho }}_{t+1}(\overline{\theta })$ became part of the prior PDF of the next step. Once the lower limit was achieved or the maximum number of steps was reached, we obtained the approximate posterior PDF of the BGM FASt parameters.

2.4 Evaluating capabilities of BGM FASt

The parameters inferred with the BGM FASt framework are conditioned to the MS used to derive them and, therefore, to the Galactic ingredients of its receipt. That is why we can consider it as an additional capability of BGM FASt to obtain the approximate posterior PDF of a set of parameters under different combinations of MS ingredients. This lets us evaluate the individual components of the BGM Std modelling. Additionally, the BGM FASt evaluations can be done using particular regions of the CMD more sensitive to the explored Galactic ingredient, or considering alternative summary statistics or distance metrics to the one presented in Eq. (4). Sect. 3.2 contains a practical application of this concept focusing on different SEMs.

3 BGM FASt G13 implementation

In this section we present the ingredients used in the current implementation of BGM FASt for the best fitting to the full-sky Gaia data up to limiting magnitude G = 13.

3.1 Mother simulation ingredients

Most of the MS ingredients used in this work are described in detail in Robin et al. (2022). We used, for the thin disc, the non-parametric SFH implemented there. In the case of the thick disc SFH, we considered two different possibilities, one flat and another Gaussian centred on 10 Gyr ago and with a standard deviation of 2 Gyr. In both cases, the thick disc SFHs span from 8 to 12 Gyr ago, are truncated outside this range, and their total surface density is the same. The reason why we considered two different shapes for the thick disc SFH is that we aimed to run different BGM FASt executions, some fitting the thick disc SFH (for which an MS with a flat thick disc SFH is a better option to avoid influencing the final results), and some others fixing it (in these cases an MS with a Gaussian SFH is more realistic). The values of the SFHs for thin and thick discs in the MSs are found in Table 2.

For the IMF, we considered a three-truncated power-law, ξ (M) ∝ M^−α, with slopes (α₁, α₂, α₃) = (1.0, 1.7, 2.4) corresponding to the mass ranges shown in Table 2. The slopes are very similar to the values of (α₁, α₂, α₃) = (0.85, 1.76, 2.5) used in Robin et al. (2022). Regarding the SEMs, as was noted in Sect. 2.1, we used PARSEC and STAREVOL tracks.

3.2 Parameter space for IMF and SFH fitting

In Mor et al. (2019), we chose a 15-dimensional space including the three slopes of the IMF, nine parameters of a non-parametric SFH covering the range of ages of the thin disc, its radial scale length, and the volume stellar mass density of the young and old thick discs. In this new work, two different frameworks were implemented: one in which we considered a 13-dimensional parameter space composed of two IMF slopes (α₂ and α₃) and 11 parameters of a non-parametric SFH for the thin disc; and another in which the thick disc SFH was unfrozen including in the fit four additional parameters describing it, resulting in a 17-dimensional parameter space. Adding the thick disc SFH parameters in the fit appeared as a possibility thanks to the increase in limiting magnitude. At G < 12, the thick disc represented ≈ 7% of the sample, while this number rose to ≈ 10% at G < 13. This contribution is expected to grow in future works using fainter magnitudes. On the other hand, at limiting magnitude G = 13 we could not fit the low-mass IMF (too few stars), so we fixed α₁ = 1.0, in line with the value it has in the MSs. In addition, the slopes derived in this work correspond to a time- and space-averaged IMF, neglecting the possible variations in the IMF with Galactic epochs and from cluster to cluster (e.g. Jeřábková et al. 2018; Dib & Basu 2018), whose modelling and implementation can be included in BGM FASt future executions in a consistent manner.

We assumed for the current runs that the MSs were dynamically self-consistent (see Sect. 2.1), which allowed us to exclude the radial scale length from the inference parameter space. The limitation of this process is that we broke this dynamical self-consistency throughout the inference process. That is why the whole process requires iterating between BGM Std and BGM FASt to converge into a final solution with the best IMF and SFHs in a dynamically self-consistent context (see Sect. 5.2).

Table 2 displays the mass and age ranges of each parameter. Table 1 shows the set of different configurations considered in this work and named G13(P/S)-(13/17), depending on the evolutionary model that was used to construct the MS (‘P’ stands for PARSEC and ‘S’ for STAREVOL) and the number of BGM FASt fitted parameters (‘13’ or ‘17’ if the thick disc SFH was frozen or not, respectively).

Finally, we centred the Gaussian priors on the MSs values of the parameters, which are found in Table 2. This responds to the fact that we considered the MSs the departing point for BGM FASt and, therefore, the values used for their generation constituted our best knowledge before applying the inference method.

This new implementation of BGM FASt implied some modifications of the original equations (Mor et al. 2018). For instance, the computations regarding the density laws with the Einasto profiles were no longer in use. Under this consideration, we could suppress the spatial dependencies, keeping only the IMF and the

SFH as the BGM FASt parameters to infer. The weight equation, therefore, was simplified significantly to

${\omega }_j=\frac{{\displaystyle {\int }_{\Delta {\tau }_j}{\displaystyle \sum _{\odot }^{PS}(\tau )}}d\tau {\displaystyle {\int }_{\Delta M_j}{\xi }^{PS}}(M)MdM}{{\displaystyle {\int }_{\Delta {\tau }_j}{\displaystyle \sum _{\odot }^{MS}(\tau )}}d\tau {\displaystyle {\int }_{\Delta M_j}{\xi }^{MS}}(M)MdM},$(5)

Where Σ _⊙ (τ) is the star formation rate of the thin or thick disc at the solar neighbourhood at time τ. ξ (M) ∝ = M^−α is the value of the IMF at M. We considered dτ = Δτ _j to be equal to the age intervals given in Table 2⁵, and we took ΔM_j = [M_j, M_j + 0.025 M_⊙ ], M_j being the mass of the j^th star. The value of ω_j obtained from Eq. (5) was then corrected, taking into account stellar multiple systems arising from the imposed IMF and SFH (see Sect. 3.3 from Mor et al. 2018 for more details).

By construction (see Sect. 2.4), the derivation of the IMF and SFH parameters in BGM FASt depends on the ingredients used to generate the MSs and not fitted in the process. This fact confers to the BGM FASt framework the great advantage of evaluating them. In this case, we took advantage of this intrinsic feature by running BGM FASt on MSs generated using different SEMs⁶. The comparison of the obtained best fit to Gaia data allowed us to analyse the behaviour of each of them. We describe this process in Fig. 2: 1) we considered MSs made up of the same Galactic ingredients except for the SEMs, for which we generated a MS with PARSEC and another with STAREVOL; 2) we ran the BGM FASt fitting on the MSs. This can be understood as the computation, at each step, of the distance between the PS and Gaia assuming a given SEM, IMF, and SFH; 3) we obtained as an output the approximate posterior PDF of the IMF and SFH parameters conditioned to each of the MSs, and therefore to each of the SEMs used to build them. The results of this evaluation can be found in Sect. 4.2.

3.3 CMD description

As was mentioned (see Sect. 2.3), our summary statistics comprised a two-dimensional CMD containing $M_G^{\prime }$ in terms of a Gaia colour. In this case, we considered the CMD fitting within the limits $-5\le M_G^{\prime }\le 8.5$. We went up to absolute magnitude $M_G^{\prime }=8.5$ to be able to fit the entire mass range of α₂ considering the limiting apparent magnitude G = 13 of the samples⁷. We note that this choice differs from the previous implementations of BGM FASt (Mor et al. 2018, 2019), in which they fit the CMD in the range $-1\le M_G^{\prime }<5$ and the luminosity function (one-dimensional) for $5\le M_G^{\prime }<8.5$. This change was possible because we have doubled the number of stars at G < 13 with respect to G < 12 and we ensure now the statistical robustness in the computation of the distance (see Sect. 2.3).

For the Gaia colour, we used G − G_RP instead of G_BP − G_RP. This choice is supported by the following arguments: 1) G − G_RP is less affected by extinction than G_BP − G_RP, which directly implies for the latter the loss of more than 300 000 stars due to the limitations of the photometric transformations of considering the resulting most probable value of the distributions, as well as the 16^th and 84^th percentiles of each fitted parameter.

Evans et al. (2018) ( − 0.47 ≤ G_BP − G_RP ≤ 2.73) and the lack of photometric observations in the G_BP passband for some stars; 2) for the same reason, the Gaia colour G G_RP is less dependent on the selected extinction map; 3) using G − G_RP instead of G_BP − G_RP we are able to consider colder stars in the working area of our CMDs, which is crucial to fit the old SFH and the low-intermediate-mass range of the IMF, and 4) the G_BP passband in Gaia has been redefined several times, and the BTSettl grid used in the last version of BGM Std did not include the last modifications. The only drawback we found against using G − G_RP instead of G_BP − G_RP is that the former’s CMDs appear to be slightly more concentrated than the latter’s (due to the smaller wavelength coverage and the shift towards the blue of the G_BP passband), which may imply a dilution of the physical features along the colour axis. Counteracting this drawback, the error of the Gaia colour G_BP − G_RP grows much faster with respect to the error in G − G_RP when increasing the apparent magnitude. All in all, taking into account that our aim in the future is to apply our strategy to samples up to fainter limiting, we consider that the advantages of using G − G_RP overcome its drawbacks.

Regarding the binning, we considered steps of Δ(G − G_RP) = 0.05 mag and $\Delta M_G^{\prime }=0.25$ mag, which gave us a reasonable relation between the conservation of the stars’ information and a statistically sufficient number of stars per bin. Finally, we divided the fit into three different CMDs characterising different Galactic latitude ranges: |b| < 10, 10 < |b| < 30, and |b| > 30. By doing that, we avoided losing valuable information on the effects of a given combination of IMF and SFH in different regions or components (e.g. thin and thick discs, and halo). The catalogue data used for the fitting was the entire Gaia DR3 G < 13 sample of stars with positive parallaxes and available G − G_RP colours (Gaia Collaboration 2023), sampling inhomogeneously a volume of 1 (55% of the stars) to 4 (95% of the stars) kpc.

Fig. 2 Conditional inference scheme. This figure shows how the output parameters describing the IMF and the SFH were derived conditioned to the ingredients of the MSs and, more particularly, to the two SEMs considered in this work. Read Sect. 3.2 for a more detailed description of the process. Besides considering two different SEMs, we also ran BGM FASt on MSs built with two different thick disc SFHs. Considering all combinations (two SEMs and two thick disc SFHs), we have in total four MSs over which we performed the BGM FASt executions (see Tables 1 and 2).

Fig. 3 Maximum, minimum, and mean distance evolution over ABC steps for the BGM FASt execution G13P-13. The maximum distance can be considered the threshold of the distance at a given step. The pink-shaded region covers the last 70 steps, which were considered for the computation of the posterior approximate PDFs of the parameters to avoid the priors’ influence.

3.4 Best ABC parameters

The parameters to be set by the ABC include the distance thresholds, the number of accepted simulations per step, and the maximum number of steps. For the former, the goal is to establish intelligent upper and lower limits to obtain a statistically robust posterior approximate PDF of the parameters representing a scientific improvement with respect to previous studies. That is why we chose for the upper limit, δ_max, the distance between the MSs and the Gaia DR3 data. Behind this selection, we accepted that the MSs were currently the best modelling of the Milky Way within the BGM Std framework. For the lower limit, δ_min, we set an arbitrarily small value that let us explore the possibilities of the inference process without constraints.

For each step, we ran as many simulations as were needed to collect 500 accepted sets of parameters. This value was decided after checking that the final solution was stable and equivalent within the statistical fluctuations to the one derived by doubling the number of accepted simulations per step. As is seen in Fig. 3, we fixed the maximum number of steps to 100. In Fig. B.1, we show the convergence of each of the fitted parameters. As can be seen, the chosen value of 100 steps is suitable for achieving the convergence in the parameters with enough information (see Appendix B for details). No better results were observed by increasing the number of steps to 200. For the computation of the final posterior of the parameters, we discarded the first 30 steps to avoid the influence of the priors and keep the information of the equivalent solutions (degeneracies).

4 The Gaia G13 – BGM FASt scientific outputs

In Sect. 4.1, we present and analyse the approximate posterior PDFs of the inferred parameters describing the IMF and the SFHs. The best CMDs resulting from the BGM FASt fits are examined and compared in Sect. 4.2, where we evaluate PARSEC and STAREVOL evolutionary models against Gaia data.

4.1 Posterior PDFs of the inferred IMF and SFH parameters

In Table 2 and Figs. 4 and 5, we present the most probable values of the SFHs and IMF with their 16^th and 84^th percentiles. In Fig. B.1, we show the evolution of the parameters for the execution G13P-13 over the 100 processed ABC steps (see Sect. 3.4), which allows us to evaluate the degree of convergence (or divergence) of each of them. In Figs. C.1 and C.2, we provide the corner plots with the approximate posterior PDFs of the two executions with the frozen thick disc SFH, which contain information on the correlations among the parameters. In Fig. 6, we zoom in on the most relevant features of these plots for those parameters directly related to the joint inference of the IMF and SFH processed here. Detailed comments on Figs. B.1, C.1, and C.2 can be found in Appendices B and C, respectively.

Fig. 4 Posterior SFHs inferred with BGM FASt in the different executions. Those performed fixing the thick disc are found on the left, while the executions including the fit of the thick disc SFH are set on the right. For completeness, in the left plot, we also show in grey the fixed thick disc SFH. We distinguish in blue and green the SEM used to build the MS of each execution, PARSEC and STAREVOL, respectively. The marginalised posterior distributions for each of the weights of each age bin are represented as vertical shaded areas (the so-called violin plot).

Fig. 5 Posterior IMF inferred with BGM FASt in the executions G13P- 13 (blue) and G13S-13 (green). In dashed grey lines are shown the limits of the three power-law IMF ranges, spanning 0.015–0.5 M⊙, 0.5– 1.53 M⊙, and 1.53–120 M⊙. We recall that the value of the low-mass slope of the IMF is fixed to α1 = 1.0.

4.1.1 Initial mass function

Fig. 5 shows that the IMF slopes α₂ and α₃ inferred in the G13P-13 execution lead to a smooth and continuous evolution across the entire mass range (blue lines and shaded areas). The inference of these parameters in G13S-13 (green lines), however, produced a non-gradual behaviour around 1.5 M_⊙ (α₂ > α₃). If real, this would imply a notable break in the IMF that has not been observed in earlier work and is hard to reconcile with current star-formation models.

We highlight that both inferences used the same input parameters in the MSs – except for the SEM, as well as the same set of priors. Furthermore, no significant correlations are found between α₂ and α₃ in either G13P-13 or G13S-13 (see Fig. 6, top panel). More importantly, as shown in Fig. B.1, all the IMF parameters achieve a fast convergence during the ABC fit, always before the first 50 steps. Thus, the markedly different final most probable values obtained in G13P-13 and G13S-13 (see Table 2) confirm the crucial role of the SEMs in the inference process. Regarding the IMFs derived in executions G13P-17 and G13S- 17, almost the same behaviours are observed with respect to G13P-13 and G13S-13, respectively.

4.1.2 The thin and thick discs’ SFHs

In Fig. 4, we compare the SFHs resulting from the different BGM FASt executions: fixing the thick disc SFH on the left and fitting it on the right. All solutions have in common a tiny SFH at early ages (<1 Gyr ago) followed by a wide enhancement in the star formation between 1–2 and 5–6 Gyr ago and a hiatus around ≈ 5–7 Gyr ago. Each of these three main features of the thin disc evolution deserves special discussion.

Although young stars (ages <1 Gyr) in our up to G=13 sample are distributed along the full range of masses, they are more concentrated at the upper main sequence, the region of the CMD where the inference of stellar masses and ages is most reliable. That explains why our strategy allowed us to achieve very narrow distributions (small error bars) for the inferred parameters of the young thin disc SFH at ages between 0 and 2 Gyr $({\Sigma }_{\odot }^0,{\Sigma }_{\odot }^1,\text{and}{\Sigma }_{\odot }^2)$. We can state that all solutions confirm an abrupt decrease in star formation in the solar neighbourhood approximately 1–1.5 Gyr ago.

As was mentioned, all solutions confirm the presence of a wide enhancement of star formation at intermediate ages. However, while it spans from ≈ 2 to ≈ 5 Gyr ago in PARSEC solutions (G13P-13 and G13P-17), it is clearly shifted in about ∼ 1 Gyr in STAREVOL executions (G13S-13 and G13S-17). We note that this shift supports the statement in Haywood et al. (2016) that the current uncertainty in the SEMs at intermediate ages can trigger differences in the derivation of the SFH of, at least, ∼1–1.5 Gyr.

All solutions reproduce the hiatus in the SFH of the thin disc widely discussed in recent literature (see Snaith et al. 2015, subsequent studies and Sect. 5.1). Furthermore, it appears more pronounced in the case of PARSEC solutions. We find this hiatus located at a lookback time of ≈ 6.5 Gyr for PARSEC and ≈ 5.5 Gyr for STAREVOL. The correlation among points adjacent to the hiatus, which can be found in the mid panel of Fig. 6, is very low. The interpretation of this hiatus in terms of Galactic disc evolution is discussed in Sect. 5.

At lookback times >6 Gyr, PARSEC solutions present a pronounced enhancement of star formation. The STAREVOL run, although less pronounced, also detects this feature. However, at these older ages, this enhancement has to be discussed in the full scenario of the thin plus thick disc evolution. Our strategy, when considering only stars up to G=13, as done in this first paper, suffers from the degeneracy between the old stars in the CMD, especially the thin and the thick discs. In Fig. 4, we have added to the thin disc SFH (${\Sigma }_{\odot }^0\text{to}{\Sigma }_{\odot }^{10}$, ages from 0–10 Gyr) the most probable values we derive for the four thick disc SFH parameters (${\Sigma }_{\odot }^{9T}\text{to}{\Sigma }_{\odot }^{12T}$, ages from 8–12 Gyr) in the right panel and the values used in the MSs in the left panel. As can be seen in Table 2, the derived values for executions G13P-17 and G13S-17 remain very close to the initial priors, demonstrating that we cannot resolve the thick disc SFH at this limiting magnitude. This is justified by the lack of old stars at G < 13 and the limitations of using CMDs as summary statistics without taking into account chemistry and kinematics (see Sect. 5). Due to these limitations, in this work, we consider executions G13P-13 and G13S-13, obtained by fixing the thick disc SFH, more reliable than the ones derived fitting it, and we shall limit our analysis to them from now on in this article.

4.1.3 Identification and impact of parameter correlations

The consistency of these results requires a thorough and rigorous analysis of the corner plots (Figs. C.1 and C.2, and zoom images in Fig. 6) and of the figures showing the evolution of convergence of the fitted parameters (Fig. B.1). The corner plots show different zones with non-negligible values of the Pearson coefficient (R > 0.4), which quantifies the linear correlation between two parameters. First, we must highlight both for G13P-13 and G13S- 13 the important degeneracy between the third slope of the IMF and the recent SFH from 1 to 3 Gyr ago. We see, for instance, $R({\alpha }_3,{\Sigma }_{\odot }^2)=-0.800$ and $R({\alpha }_3,{\Sigma }_{\odot }^2)=-0.898$ for G13P-13 and G13S-13, respectively. This is an important correlation given that the young (τ < 3 Gyr) and massive (M > 1.53 M_⊙ ) stars are well described at limiting magnitude G = 13, representing 40% of the sample in the MSs. In Sect. 5.2, we provide insights on how to overcome this challenging degeneracy in upcoming works.

Secondly, we observe also in both solutions (even though more pronounced for G13S-13), correlations in the zone along the diagonal describing the posterior PDFs of consecutive–or second consecutive–parameters in the young SFH (τ < 4 Gyr) For example, we see correlations between ${\Sigma }_{\odot }^1-{\Sigma }_{\odot }^3$ or ${\Sigma }_{\odot }^2-{\Sigma }_{\odot }^4$,with R = 0.404 and R = 0.441 for G13P-13 and R= 0.605 and R = 0.630 for G13S-13. The smaller values of R in this region for G13P-13 than for G13S-13 make us think that the PARSEC’s solution for the young SFH is more robust. Although this type of correlation could be considered as something expected as consecutive parameters represent stars with similar features in the CMD, these results will require further analysis in future executions.

We highlight a third region of the corner plots, which is the one described by α₃ and the very old SFH (τ > 9 Gyr in G13P- 13 and 6 < τ/Gyr < 8 in G13S-13). For G13S-13, this region of correlations can be extended for the whole combinations of ${\Sigma }_{\odot }^7$ and ${\Sigma }_{\odot }^8$ with the parameters describing the young SFH, mainly ${\Sigma }_{\odot }^1$ and ${\Sigma }_{\odot }^2$ . If we compare the position of high-mass stars in the CMD with that of old stars, we find that a non-negligible number of both falls in the red clump (more than 900 000 stars with M > 1.53 M_⊙ and around 700 000 stars with τ > 8 Gyr in G13P- 13). Young–and massive–stars and old stars clearly constitute two different populations in the Galaxy. The observed correlations in this region demands for including additional sources of information in the future, such as chemistry or kinematics, in order to break the current degeneracies in the CMD.

Fig. 6 Correlations among some of the BGM FASt inferred parameters with their corresponding projected approximate posteriors PDF (with a Gaussian fit) from our fiducial execution G13P-13. The resulting most probable value and 16th and 84th percentiles are also marked with solid and dashed lines, respectively. It is also indicated with a dashed line the median of the distributions. In black, at the top right of each plot, we show Pearson’s correlation coefficient. High values of this coefficient are highlighted with an intense blue frame in the plot. Finally, the parameters adopted by the MS are marked with a magenta cross and the final values with an orange star. The full corner plot is presented in Appendix C.

4.2 BGM CMDs versus Gaia data: An evaluation of the stellar evolutionary models

In Appendix D (Figs. D.1 and D.2), we present the complete set of CMDs describing the comparison between the MSs and the final BGM FASt PSs (top and bottom panels of each figure, respectively) with the Gaia data for executions G13P-13 and G13S-13. To provide a more focused analysis, Fig. 7 highlights the most relevant subplots of the entire set of CMDs. These plots allow us to analyse the regions of the CMD in which the new IMF and SFH significantly improve the resemblance with the Gaia catalogue. This process also allows us to assess the caveats of this method and the still most significant discrepancies between Gaia and the model.

The evolution of the distance metric defined in Eq. (4) over the inference process (similarity between catalogue – initial MS and catalogue – final PS) is presented in Table D.1. For G13P-13, we start with an initial value of 1 772 695 and it is reduced to a residual distance of 922 755, which implies an improvement of almost a factor of 2. The outcome of BGM FASt G13S-13 yields a less favourable solution. It begins with a higher initial distance to Gaia for the MS (1 825 293) that is reduced by only a factor of 1.3 after applying BGM FASt (1 362 537), far from the clear improvement in G13P-13.

Regarding the analysis of the CMDs, despite the overall agreement, several regions show a discrepancy between the model and the data. The discrepancies can be seen in the third column of Figs. D.1 and D.2, where we present the distance for each given bin (higher the distance, larger the discrepancy), and in the fourth column, where the colour represents the simple difference between the model and the data star counts. Focusing first on the CMDs of execution G13P-13 (Fig. D.1 with a zoom in Fig. 7), it is interesting to compare the discrepancies between the initial MS and Gaia (top) with the final discrepancies after the BGM FASt execution (bottom):

• Upper main sequence: this region initially showed a significant excess of stars in the MS, especially visible at low and intermediate latitudes (as expected for massive stars). After the fit, the discrepancy was almost completely resolved thanks to the new IMF and SFH, responsible for reducing the number of stars from 8 572 604 in the MS to 7 113 842 in the final PS, much closer to the value of 7 266 095 stars present in the Gaia G < 13 sample. However, a slight shift of the entire main sequence in colour is still present. This region is where BGM FASt has been most efficient.

• Asymptotic giant branch (AGB): in this region, the distance δ_p showed large values, mostly at intermediate and low latitudes, and did not improve much after the fit. Comparing the first and second columns of Fig. D.1, we see that the upper AGB is absent in the MS. This is due to the tracks, which do not cover this region well, so the fit cannot improve the resemblance between the simulations and the catalogue in it. Furthermore, in the fourth column, we observe that the absolute differences in star counts are low compared to the significant number of stars present in the main sequence and the red giant branch. Therefore, this region has no impact in the fit.

• Red clump: the MS produced a red clump that is roughly positioned correctly but appears less compact (more extended in colour) than in the catalogue. The new IMF and SFH fit improves this region, particularly at high and mid latitudes. However, the discrepancies persist at lower latitudes, which could be partially due to some remaining inaccuracies in the extinction map used in the MS or/and to nearby star-forming regions.

• Subgiants and lower main sequence: these regions initially exhibited a slight lack of stars in the MS, which was reduced after the fit. Although the lack is still present, we suspect it could be completely removed if the colour shift is corrected.

In the case of using STAREVOL (Fig. D.2), similar patterns are seen, but the discrepancies are generally stronger. Notably, one of the major issues with these tracks is the position of the red giant branch, which is redder in the simulation than in the data. This discrepancy could not be corrected with our process, likely contributing to the lower accuracy of the fit obtained with STAREVOL.

Complementing the physical information presented in the CMDs, the differences in the mass-age distribution resulting from the G13P-13 and G13S-13 executions are shown in Fig. 8. As can be seen, the discrepancies between the two SEMs show a non-continuous trend. We have to highlight that this behaviour is reflected in all processes of IMF and SFH parameter inference and it is the consequence of two effects: 1) the differences in the tracks of the PARSEC and STAREVOL evolutionary models and 2) the lack of correspondence between the IMF and SFH derived in each case (see Table 2 and Fig. 4). While the effect of the particular tracks is difficult to understand in the mass-age distribution, it is relatively easy to observe the consequences of using different IMFs or SFHs in this space. For example, we clearly see the shift of 1–1.5 Gyr commented in Sect. 4.1.2.

Fig. 7 Low-latitude ( |b| < 10) CMDs illustrating the distance per bin, as defined in Eq. (4), for both the initial MSs (left) and the final solutions (right) of the BGM FASt executions G13P-13 (top) and G13S-13 (bottom).

Fig. 8 Bin-wise differences in stellar number density in mass-age space between the two final solutions obtained considering a fixed thick disc SFH, G13P-13 and G13S-13. For reference, we show in the mass axis the limits of the three power-law IMF ranges shown in Table 2 (vertical dashed grey lines).

5 Discussion

In Sect. 5.1, we compare our findings on the IMF and SFH to recent studies in the literature. Section 5.2 presents a deep analysis of the caveats and limitations of the present work and advances new strategies under development.

5.1 IMF and SFH of the solar vicinity

The IMF and its critical role in many fields of stellar, Galactic, and extragalactic astrophysics are still a matter of debate within the community. The results presented in this work (see Sect. 4.1.2) illustrates how the derivation of this fundamental function critically depends on the accuracy of the stellar evolution models being used. Our IMF derivation, based on a full-sky sample up to magnitude G = 13 where ∼60% of the stars fall in the range of masses 0.5–1.53 M_⊙, suggests that, apart of the dependency on the SEMs, the values for α₂ derived here are statistically robust when compared with other recent estimates. Furthermore, our derivations for α₃ have a high added value thanks to the high number of massive stars in our limiting magnitude sample compared to other solutions (these stars are intrinsically bright so we are considering here all the massive Gaia sources up to approximately 1–2 kpc from the Sun). Nonetheless, as discussed in Sect. 4.1.3, the high correlation between α₃ and ${\Sigma }_{\odot }^2$ demands further developments (see Sect. 5.2).

In the bottom right panel of Fig. 9, we compare our G13P-13 and G13S-13 IMFs with widely used values from the literature and recent estimates. We find that, in all the works shown in the figure except for the STAREVOL solution G13S-13, the slopes of the IMF progressively grow with mass (α₁ < α₂ < α₃). As mentioned in Sect. 4.1, the fact that α₃ < α₂ in G13S-13 could be one of the drawbacks of this execution. The value of the second IMF slope we obtained in G13P-13 for the mass range 0.5–1.53 M_⊙ has small error bars (see Table 2) and it is not correlated with other parameters (see Fig. C.1). This value is slightly below the commonly used slopes from Kroupa (2001) and Salpeter (1955) (see Fig. 9). For masses larger than 1.53 M_⊙ , the value of the IMF slope derived using PARSEC (G13P-13) is consistent with the Salpeter’s IMF and with recent values derived from completely independent methods such as the one obtained by Dickson et al. (2023) fitting a multi-mass model to globular clusters data.

We next analyse results for the SFH of the solar vicinity, whose derivation has seen a tremendous revival in recent times following the publication of large astrometric, spectroscopic, and photometric surveys. Intending to provide a synthetic view of the local SFH, in Fig. 9 we compare the results obtained using different methods. We believe that the representative compilation shown in these figures present many of the features that are currently under debate. The figures include the results of three different approaches to the problem. One of the techniques pursue fitting Galactic chemical evolution models to the abundances derived from spectroscopic data (e.g. Snaith et al. 2015; Haywood et al. 2016; Spitoni et al. 2023; Palla et al. 2024). Another consists on mimicking the CMDs obtained from photometric observations (in recent years, mainly Gaia data) by considering population synthesis models of the Galaxy (e.g. Mor et al. 2019; Mazzi et al. 2024; Gallart et al. 2024, and the present approach). All studies based on Gaia CMDs provide estimates of a dynamically evolved SFH (an SFH convolved with the effects of stellar dynamical mixing) rather than the true SFH, as is explained in, for example, Gallart et al. (2024). A third group aims to derive the SFH using specific tracers, such as the white dwarf population. As detailed in the recent review by Tremblay et al. (2024), local volumes of white dwarfs are essential for studies of Galactic SFH.

In Fig. 9, we provide our fiducial PARSEC SFH solution together with literature values in three different ways: the top left panel shows the comparison with works that derive the SFH in units of surface density (M_⊙ pc⁻² Gyr⁻¹); the top right panel shows comparisons with normalised (relative) SFHs; finally, the bottom left panel shows comparisons in terms of the SFH per local volume density (in M_⊙ pc⁻³ Gyr⁻¹).

Overall, we observe how most of the recent derivations of the local SFH present a first enhancement of star formation at intermediate ages followed by a pronounced hiatus that is succeeded by a complex and discrepant evolution at ages older than 8 Gyr, just when the thick disc begins to play a relevant role. In more detail, as can be seen in Fig. 9, both the width of the enhancement and the location of the hiatus change from one solution to another. In the following paragraphs, we discuss how different models, assumptions, sets of data, and fitting strategies contribute to the complex understanding of these discrepancies.

Snaith et al. (2015) derived an SFH for the inner disc (R < 10 kpc) by fitting a closed-box-like system in this region to the age-[Si/Fe] relation of the HARPS GTO high-resolution (R = 110 000) spectra sample of FGK stars from Adibekyan et al. (2012), with ages determined in Haywood et al. (2013). According to their model, most accretion takes place at early times, when substantial star formation has not yet occurred. Therefore, only a single infall of gas is considered in this case in the very early times of the Galaxy, and the later SFH variations are a consequence of its internal evolution; for example, the formation of the bar or the end of the thick-disc phase (Haywood et al. 2016). According to the authors, the cessation of star formation during a period of ∼ 1–2 Gyr could reproduce the bimodal distribution found in the [Fe/H]-[α/Fe] diagram of APOGEE data, among other spectroscopic surveys (Haywood et al. 2016). Looking at Fig. 9 (top right), our fiducial PARSEC SFH matches their normalised SFH in the interval 1–10 Gyr ago with only a shift of 1–2 Gyr.

Different approaches are taken in Spitoni et al. (2023) and Spitoni et al. (2024). Both works consider Galactic chemical models in which different infalls of metal-poor gas from external accretion could be responsible for the distinguished high- and low-alpha sequences. Spitoni et al. (2023) compare their three-infall model to the Gaia DR3 GSP-Spec chemical abundances for α elements, focusing their efforts on the characterisation of the population of young and massive stars with a deficiency in metallicity [M/H] in the solar neighbourhood (guiding radius between 8.1 and 8.4 kpc). In contrast, Spitoni et al. (2024) adjust a two-infall chemical model to the [Fe/α] versus [α/H] relation of a sample of red giant stars in the solar neighbourhood (galactocentric distances between 7 and 9 kpc and |z| < 2 kpc) observed by APOGEE. It is worth mentioning that none of the previous works included stellar migration, a process that was instead considered in the three-infall model of Palla et al. (2024), leading to similar results as the three-infall model of Spitoni et al. (2023). An interesting feature of Palla et al. (2024) is that they use open clusters to study the variation in the SFR with the Galactic radius. All in all, when comparing our G13P-13 SFH to the solutions arising from the chemical evolution models of Spitoni et al. (2023), Spitoni et al. (2024) and Palla et al. (2024), and also with that of Snaith et al. (2015) discussed above, we can confirm that our approach, different both in modelling and data, stands coherent with the existence of the hiatus, also unveiled in our previous fit to the Gaia DR2 up to G = 12 catalogue (Mor et al. 2019).

The comparison of our results with those from Mazzi et al. (2024) and Gallart et al. (2024), both based on CMD-fitting techniques, is difficult due to the lack of detail of the former at ages older than 2–3 Gyr (they fit a logarithmic SFH instead of a linear one) and the high level of discretisation at recent times for the latter. It is especially interesting to contrast our inference with that of Gallart et al. (2024). Although based on a similar methodology – a CMD-fitting process applied to Gaia data, some of the strategies and basic ingredients used in their work are significantly different than ours. We shall mention some of them here and, more importantly, discuss the scientific outputs these differences may entail: (1) Gallart et al. (2024) used a volume complete sample restricted to the 331 312 stars currently in the 100 pc solar neighbourhood (GCNS), whereas our all-sky catalogue, also based on a well-defined observational cut (apparent limiting magnitude G = 13), includes more stars by a factor of 20 ( ∼ 7 million sources); (2) while Gallart et al. (2024) imposed the Kroupa et al. (1993) IMF, the simultaneous fit of the IMF and the SFH done in this work may be crucial to successfully understand the dependencies and degeneracies among both fundamental quantities. For instance, as seen in the bottom right panel of Fig. 9, the high value of IMF slope at masses >1 M _⊙ used by Gallart et al. (2024) could be behind their large values of the SFH at recent times (τ < 2 Gyr); and 3) as is deeply discussed in Sect. 3.3, we know that systematic differences are also present due to the different SEMs used, PARSEC and STAREVOL in our case and BASTI-IAC in Gallart et al. (2024). Even considering all the effects mentioned so far, further work is needed to justify the significant differences observed in the bottom left panel of Fig. 9.

Our fitting technique has been demonstrated to be statistically robust when deriving the more recent phases of disc evolution. In particular, we observe a sudden decrease in star formation in the last 3 Gyr, and the SFH values describing this behaviour present tiny error bars (see Fig. 4). This trend seems to be, in terms of shape, in agreement with the SFH determined by Mazzi et al. (2024), which is the result of collapsing the different SFHs as a function of the distance from the Galactic plane obtained from fitting the Gaia DR3 CMD within a cylinder of 200 pc radius centred on the Sun and spanning 1.3 kpc above and below the plane. However, the total surface mass density obtained in their work, 118.7 ± 6.2 M_⊙pc⁻², appears much larger than ours $\left({\Sigma }_{\odot }^T={53.71}_{-6.89}^{+5.48}{M}_{\odot }{\text{pc}}^{-2}\text{for G13P-13}\right)$ and other works in the literature (see the area under the curves in top-left panel of Fig. 9). There are several aspects to take into account to understand these significant discrepancies. In our view, the derivation of this parameter requires the full achievement of the dynamical self-consistency in the process, a challenge that we plan to address through the iteration between the BGM FASt solutions for the IMF and SFH, and the fit to the gravitational potential to kinematic and densities from Gaia data (see Fig. 2 from Robin et al. 2022). Until this is achieved, estimates of this parameter can suffer for significant biases coming from the combined use of different volume samples and simplified expressions for the radial or vertical spatial distributions. As an example, as found in Robin et al. (2022), a fit to kinematic Gaia DR2 data indicates that the effective radial scale length of a population varies with the distance from the Galactic plane. Facing this would require a more complex modelling than classical exponential or sech2 vertical scale heights, which could affect the vertical structure of the SFH infered by Mazzi et al. (2024). Additionally, recent tools derived by the GaiaUnlimited team (e.g. Castro-Ginard et al. 2023) shall be considered to correct the survey selection function when using different volume samples.

Regarding the present-day SFH in the solar neighbourhood, the Galactic chemical evolution model proposed by Prantzos et al. (2018) predicts a value of 2–5 M _⊙ pc⁻²Gyr⁻¹ (see mid panel of Fig. 9 in the authors’ article). Independently of the SEM being used, our results indicate an SFR of about 1.3–2.5 M _⊙ pc⁻²Gyr⁻¹ for the last Gyr (see in Table 2 values for ${\Sigma }_{\odot }^0$ and ${\Sigma }_{\odot }^1$ ), very close to the value derived in Mor et al. (2019). Furthermore, our results show a fast and stable convergence of these parameters in all BGM FASt executions (see Fig. B.1). Recent estimations from Mazzi et al. (2024) and Spitoni et al. (2023, 2024) point in the same direction, deriving a present-day of the SFH close to the lower limit of the wide range of 2–5 M _⊙ pc⁻²Gyr⁻¹ from Prantzos et al. (2018).

Fig. 9 Comparison of our inferred SFH and IMFs with the literature. Top left: comparison of our fiducial solution G13P-13 with the results from Mazzi et al. (2024), Spitoni et al. (2024), Spitoni et al. (2023), and Mor et al. (2019), all presented in original units (surface density). The error bars of our SFH and IMFs correspond to the 16th and 84th percentiles of the marginalised posterior PDFs of each parameter, whose values can be found in Table 2. Top right: findings of Snaith et al. (2015) contrasted with our fiducial solution normalised to the maximum value of each SFH. Bottom left: Comparison of the solution of our execution G13P-13 with the SFH from Gallart et al. (2024), presented in terms of volume density. To convert our surface densities to volume densities we employed the scale heights of each population in BGM. Bottom right: slopes of the IMFs derived in executions G13P-13 and G13S-13 alongside those from Salpeter (1955), Kroupa et al. (1993) (assumed in Gallart et al. 2024), Kroupa (2001), and Kirkpatrick et al. (2024). Vertical dashed lines delineate the mass ranges of this work, detailed in Table 2.

5.2 Limits of our analysis

In this section, we summarise the main caveats of our analysis and how can we overcome them in the future. First and foremost, as discussed in Sect. 3.3, the choice of the SEMs introduces a considerable systematic uncertainty, even when rather canonical models are used. This work has been a first approximation of the problem, but we are aware that a more thorough analysis must be done considering alternative tracks such as BASTI-IAC and PARSEC v2.0, which we shall implement in future executions of BGM FASt.

Second, as is shown in Sect. 4.1.2, some parameters, such as the SFH of the old thin disc, are highly degenerated with the parameters of the thick disc SFH, and these correlations are impossible to break with the magnitude-limited sample G < 13 used here. Furthermore, parameters related to old stars in general (τ > 8 Gyr) suffer from a lack of convergence throughout the ABC process (see Fig. B.1) and, as demonstrated in Fig. 6 and Sect. 4.1.3, some of them are clearly correlated. The effects of degeneracy, correlation, and convergence are usually neglected in similar works, resulting in underestimated uncertainties for the derived SFH. In conclusion, the current executions of BGM FASt up to G = 13 does not allow us to break the degeneracies and correlations found in the earliest phases of disc evolution (stars with ages >8 Gyr). We shall tackle this problem in the future by considering samples up to further limiting magnitudes (including more old stars) and additional information such as chemistry or kinematics (better characterising the populations). These approximations will also be useful to tackle the current degeneracies between α₃ and the young stars SFH (see Sect. 4.1.3).

Third, as was mentioned in Sect. 3.2, throughout the BGM FASt inference we break the dynamical self-consistency of the MSs, which is clearly seen in Table 2 comparing the total surface densities of the MSs and the final solutions. We are working in two different possibilities to overcome this problem: 1) the iteration between BGM Std and BGM FASt until reaching a dynamical self-consistent solution with the best IMF and SFH; 2) the inclusion of new constraints in BGM FASt to maintain the dynamical self-consistency in the process.

6 Conclusions and future work

We present the joint inference of the SFH and the IMF of the Galactic disc in the solar neighbourhood obtained by iteratively fitting our dynamically self-consistent BGM to the Gaia DR3 CMD up to G = 13. The process was carried out using an improved version of the BGM FASt tool developed within our team. Given the nature of the process, the inference of these two Galactic properties – the SFH and the IMF – was done by deeply analysing the critical role of the adopted SEMs in our population synthesis model. We consider the derivation with PARSEC SEM to be our fiducial solution due to its significantly higher likelihood compared to the STAREVOL solution.

From Gaia DR3 and BGM FASt up to G = 13, we report a new derivation of the slopes of the IMF in the range of [0.5, 120] M _⊙ . Although statistically robust, they appear to be discrepant when using different SEMs. In our fiducial execution with PARSEC, for the range 0.5–1.53 M_⊙ the slope takes a value of ${\alpha }_2={1.45}_{-0.12}^{+0.19}$, while for masses larger than 1.53 M_⊙ we obtain ${\alpha }_3={1.98}_{-0.05}^{+0.13}$. Using STAREVOL, the inferred values are ${\alpha }_2={2.48}_{-0.11}^{+0.09}$ and ${\alpha }_3={1.64}_{-0.02}^{+0.15}$. These parameters were derived with a remarkable level of convergence and high statistics due to the nature of the G < 13 sample. For the SFH of the Galactic disc at the solar neighbourhood, all our solutions have in common a clear hiatus around ≈ 5–7 Gyr ago with a shift of ≈ 1 Gyr depending on the evolutionary model. The presence of a secular mechanism is supported by the wide plateau in the bump between ≈ 6 and ≈ 2 Gyr ago, and the high convergence achieved for the young SFH confirms an abrupt decrease in star formation in the solar neighbourhood approximately 1–1.5 Gyr ago. The fact that the statistically robust SFH derived here has several features in common with recent determinations in the literature (the hiatus, the burst, and the abrupt recent decrease) indicates that our SFH is ready for the complex astrophysical interpretation of the disc evolution.

In this work, we have identified the areas in the CMD where PARSEC and STAREVOL match Gaia or need improvement. In all BGM FASt executions, we tackled the initial discrepancies at the upper main sequence and in the region of the subgiants, while issues in the asymptotic giant branch and the red clump seem to be related to the stellar tracks and are therefore not solvable when fitting the IMF and the SFH. PARSEC SEMs produce smaller final distance metrics with respect to Gaia, whereas non- continuous features are observed as residuals of the fit in the case of STAREVOL.

The future BGM FASt executions being prepared pursue the goal of achieving a robust fitting of the SFH at old ages. For that, the current implementation of the code is ready to significantly increase the limiting apparent magnitude of the sample and, similarly to other analyses, to take into account the selection function of the Gaia survey (Cantat-Gaudin et al. 2023; Castro-Ginard et al. 2023). The code is also ready to evaluate and even derive the stellar binary fraction, another key ingredient that plays a crucial role in the dynamical evolution of any stellar system. In parallel, our new BGM FASt executions will require reintroducing the BGM dynamical self-consistency. This will be tackled by incorporating spatial density and kinematics constraints. Last but not least, future BGM FASt executions must incorporate the use of chemical and kinematical data to address the complex evolution of the old thin and thick discs. All this work will be developed in parallel with the use of BGM FASt as a tool for testing, not only new SFHs and IMFs coming from literature but also other key ingredients of the population synthesis models arising in the near future as by-products of new Gaia Data Releases.

Acknowledgements

We thank the anonymous referee for a constructive report that helped to revise and improve the quality of the manuscript. 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). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This work was (partially) supported by the Spanish MICIN/AEI/10.13039/501100011033, by “ERDF A way of making Europe” by the “European Union” through grant PID2021-122842OB-C21, by the Institute of Cosmos Sciences University of Barcelona (ICCUB, Unidad de Excelencia ‘María de Maeztu’) through grant CEX2019-000918-M, by the OCRE awarded project Galactic Research in Cloud Services (Galactic RainCloudS, funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement no. 824079) and by MCIN with funding from European Union NextGenerationEU(PRTR-C17.I1) and by Generalitat de Catalunya. All BGM simulations and part of BGM FASt fits were executed on computers from the Utinam Institute of the Université de Franche-Comté, supported by the Région de Franche- Comté and the Institut des Sciences de l’Univers (INSU). FA acknowledges financial support from MCIN/AEI/10.13039/501100011033 and European Union NextGenerationEU/PRTR through grant RYC2021-031638-I.

Appendix A Main characteristics of the PARSEC and STAREVOL stellar models

Comparing PARSEC and STAREVOL evolutionary models, one notices small differences in the treatment of the convection, both using the Mixing Length Theory, but with different α_MLT of 1.74 in PARSEC and 1.6264 in STAREVOL. Overshooting is taken into account in PARSEC from the convective core, not in STAREVOL. None of the two sets of stellar models implemented in BGM includes the rotation. The mass loss coefficient, η, (Reimers 1975) is 0.5 in STAREVOL and 0.2 in PARSEC. The solar abundances differ slightly, being from Caffau et al. (2011) in the case of PARSEC and from Asplund et al. (2009) in STAREVOL. For the initial helium abundances, the values are slightly different as well, PARSEC using the formula Y = 0.2485 + 1.78 × Z, while STAREVOL assumes ΔY/ΔZ = 1.29.

Appendix B Evolution and convergence of the BGM FASt fitted parameters

We show in Fig. B.1 the evolution of the 13 parameters inferred in executions G13P-13 (blue) and G13S-13 (green) over the 100 ABC steps, as well as their final most probable values (MPV) with 16^th and 84^th percentiles as error bars. In the case of G13P- 13 solution (MS built using PARSEC stellar tracks and thick disc SFH not considered in the inference), we observe a strong convergence of the two slopes of the IMF (α₁ and α₂), as well as the parameters describing the young SFH up to 3 Gyr ago. We do not find strong trends for ${\Sigma }_{\odot }^4,{\Sigma }_{\odot }^5$, even though we cannot con-sider these parameters effectively converged. In the case of ${\Sigma }_{\odot }^6$ , the stability of its MPV contrasts with its non-negligible error bars. The most important problems appear for the SFH parameters describing stars with τ > 6 Gyr, for which we find a lack of convergence and important trends. As explained in Sect. 5.2, we propose this problem to be related with the lack of information we have for old stars at G < 13 and using only photometric data.

Regarding the convergence of the parameters for G13S-13 (same as G13P-13 but using STAREVOL stellar tracks to build the MS), we find a similar picture as in G13P-13 for the two slopes of the IMF, the very young SFH (τ < 1 Gyr) and the parameters describing the old SFH (6 > τ/Gyr > 8 Gyr). However, in this case we observe a worrying steep evolution of SFH parameters for stars with 1 < τ/Gyr < 4. On the other hand, more stable values are found for ${\Sigma }_{\odot }^9$ and ${\Sigma }_{\odot }^{10}$, even though their important uncertainties are evidence of a lack of characterisation of stars at those ages.

Appendix C Corner plots of the IMF and SFH inferred parameters

We show in Figs. C.1 and C.2 the posterior distributions and the correlations of the BGM FASt inferred parameters for executions G13P-13 (MS built using PARSEC stellar tracks and thick disc SFH not considered in the fit) and G13S-13 (same but using STAREVOL stellar tracks), respectively. The Pearson coefficient R, shown in the top-right corner and highlighted in the border of each subplot of the figures, characterises the linear correlation between two distributions. It can take values from 1 (linearly anticorrelated) to − 1 (linearly correlated), being 0 the representation of two linearly independent distributions. A description of the most relevant correlations among parameters is given in Sect. 4.1.3.

Appendix D Colour magnitude diagrams of the fiducial executions

Figures D.1 and D.2 expand the information provided in Fig. 7 and show the complete picture of the comparisons between Gaia and our model before and after the BGM FASt inference for executions G13P-13 (MS built using PARSEC stellar tracks and thick disc SFH not considered in the fit) and G13S-13 (same but using STAREVOL stellar tracks), respectively. In Table D.1 we present the quantification of the comparisons shown in the aforementioned figures, with the evolution of the distance and the difference in number of stars between Gaia and the model in terms of the latitude bin. We complement here Sect. 4.2 with additional comments analysing this table.

In the first column of the table and for each dataset (G13P- 13 and G13S-13 before and after BGM FASt, and Gaia G < 13 data), we show the total number of stars in the sample and the contribution of each latitude bin. The initial guess in the MSs has in all cases too many stars, and the distribution shows an excess of them at low latitudes compared to Gaia. This is completely solved by the new SFH and IMF. Especially for G13P-13, the discrepancy in the total number of stars with respect to Gaia data improves from an initial relative contribution of 17% to a residual value of 2%, with an almost perfect distribution in the different latitude bins. This is probably mainly due to the significant reduction of stars in the upper main sequence at |b| < 10 thanks to the new SFH and IMF pointed out in Sect. 4.2. Simi-lar conclusions can be made regarding G13S-13. However, while the distribution between latitude bins clearly improves, the total number of stars in the final PS still falls far from the value in Gaia, with a relative discrepancy of ∼ 8% mainly concentrated in mid-latitudes (10 < |b| < 30).

In the third column of each dataset in Table D.1 we find the total distance between the model and Gaia before and after the inference and the contribution to this value of the stars in the different latitude bins. Taking the results for G13P-13, we see that the contribution of stars at |b| < 10 diminishes 10 points throughout the BGM FASt process, from 48.9% to 38.7%, while the contribution to the total discrepancy at mid and especially high latitudes increases significantly. This is in concordance with the caveats of this work, explained in detail in Sects. 5.2 and 6. Our fit up to magnitude G = 13 lets us mainly tackle issues with young (τ < 2 − 3 Gyr) and intermediate-mass stars, while the characterisation of old and low-mass stars remains unsolved. The former population has typically more circular orbits and less vertical heating, which locates their stars in the very thin disc, extensively observed at low latitudes. On the other hand, the old star population is found further away from the Galactic plane and being part of both the thin and the thick discs. That is the reason why their presence is more noticeable at high latitudes. Therefore, the limitations of this work and the nature of the star populations explain why the contribution to the total distance at high latitudes increases significantly with the inference process.

Fig. B.1 Evolution of the median at each step of the 13 BGM FASt inferred parameters across the 100 ABC steps (refer to the derivation process in Sect. 2.3) for executions G13P-13 (blue) and G13S-13 (green). At the end of each plot, we present the final most probable values (MPVs) of these parameters, along with their 16th and 84th percentiles as determined in Sect. 3.4.

Fig. C.1 Corner plot containing the parameters derived from our fiducial execution G13P-13, including two slopes of the IMF and 11 parameters for the SFH of the thin disc. At the top of the columns, the projected approximate posterior PDF of each parameter with its Gaussian fit is shown, as well as the resulting most probable value and 16th and 84th percentiles, which are also marked with solid and dashed lines, respectively. It is indicated with a dashed line the median of the distributions. In black at the top right of each plot we show Pearson’s correlation coefficient. High values of this coefficient are highlighted with an intense blue frame in the plot. Finally, the parameters adopted by the MS are marked with a magenta cross and the final values with an orange star.

Fig. C.2 Same as in Fig. C.1 but for execution G13S-13.

Fig. D.1 Top: from left to right-hand side, columns show CMDs of the PARSEC MS used in execution G13P-13, Gaia DR3 G < 13 data, distance per pixel between them (each of the elements of the summation in Eq. (4)), and their absolute difference in number of stars. From the top to the bottom rows, we find the CMDs for the three considered latitude ranges 30 < |b| < 90, 10 < |b| < 30, 0 < |b| < 10. Bottom: same as on top but for the final BGM FASt PS obtained from the application of the fitting process described in Sect. 2.1 on PARSEC MS.

Fig. D.2 Same as in Fig. D.1 but for the STAREVOL MS used in execution G13S-13.

References

All Tables

All Figures

Fig. 1 Flux diagram of the process involving BGM Std, BGM FASt, and the CMD fitting technique. BGM Std is summarised in Sect. 2.1; we present BGM FASt in Sect. 2.2, where the meaning of the subscripts and superscripts i (indicating the ith Galactic component, e.g. thin disc, thick disc, halo, and bar) and j (associated with the jth interval of the reduced parameter space) are explained in detail; regarding the fitting process, see in Sect. 2.3 the different components of the process as well as an extended explanation of the subscripts and superscripts k (referred to kth iteration within an ABC step), n (indicating the nth bin of the CMD), and t (corresponding to the tth ABC step). See Mor et al. (2018) for a detailed explanation of the equations in the central panel.

In the text

Fig. 2 Conditional inference scheme. This figure shows how the output parameters describing the IMF and the SFH were derived conditioned to the ingredients of the MSs and, more particularly, to the two SEMs considered in this work. Read Sect. 3.2 for a more detailed description of the process. Besides considering two different SEMs, we also ran BGM FASt on MSs built with two different thick disc SFHs. Considering all combinations (two SEMs and two thick disc SFHs), we have in total four MSs over which we performed the BGM FASt executions (see Tables 1 and 2).

In the text

	Fig. 3 Maximum, minimum, and mean distance evolution over ABC steps for the BGM FASt execution G13P-13. The maximum distance can be considered the threshold of the distance at a given step. The pink-shaded region covers the last 70 steps, which were considered for the computation of the posterior approximate PDFs of the parameters to avoid the priors’ influence.
In the text

Fig. 4 Posterior SFHs inferred with BGM FASt in the different executions. Those performed fixing the thick disc are found on the left, while the executions including the fit of the thick disc SFH are set on the right. For completeness, in the left plot, we also show in grey the fixed thick disc SFH. We distinguish in blue and green the SEM used to build the MS of each execution, PARSEC and STAREVOL, respectively. The marginalised posterior distributions for each of the weights of each age bin are represented as vertical shaded areas (the so-called violin plot).

In the text

	Fig. 5 Posterior IMF inferred with BGM FASt in the executions G13P- 13 (blue) and G13S-13 (green). In dashed grey lines are shown the limits of the three power-law IMF ranges, spanning 0.015–0.5 M⊙, 0.5– 1.53 M⊙, and 1.53–120 M⊙. We recall that the value of the low-mass slope of the IMF is fixed to α1 = 1.0.
In the text

Fig. 6 Correlations among some of the BGM FASt inferred parameters with their corresponding projected approximate posteriors PDF (with a Gaussian fit) from our fiducial execution G13P-13. The resulting most probable value and 16th and 84th percentiles are also marked with solid and dashed lines, respectively. It is also indicated with a dashed line the median of the distributions. In black, at the top right of each plot, we show Pearson’s correlation coefficient. High values of this coefficient are highlighted with an intense blue frame in the plot. Finally, the parameters adopted by the MS are marked with a magenta cross and the final values with an orange star. The full corner plot is presented in Appendix C.

In the text

	Fig. 7 Low-latitude ( |b| < 10) CMDs illustrating the distance per bin, as defined in Eq. (4), for both the initial MSs (left) and the final solutions (right) of the BGM FASt executions G13P-13 (top) and G13S-13 (bottom).
In the text

	Fig. 8 Bin-wise differences in stellar number density in mass-age space between the two final solutions obtained considering a fixed thick disc SFH, G13P-13 and G13S-13. For reference, we show in the mass axis the limits of the three power-law IMF ranges shown in Table 2 (vertical dashed grey lines).
In the text

Fig. 9 Comparison of our inferred SFH and IMFs with the literature. Top left: comparison of our fiducial solution G13P-13 with the results from Mazzi et al. (2024), Spitoni et al. (2024), Spitoni et al. (2023), and Mor et al. (2019), all presented in original units (surface density). The error bars of our SFH and IMFs correspond to the 16th and 84th percentiles of the marginalised posterior PDFs of each parameter, whose values can be found in Table 2. Top right: findings of Snaith et al. (2015) contrasted with our fiducial solution normalised to the maximum value of each SFH. Bottom left: Comparison of the solution of our execution G13P-13 with the SFH from Gallart et al. (2024), presented in terms of volume density. To convert our surface densities to volume densities we employed the scale heights of each population in BGM. Bottom right: slopes of the IMFs derived in executions G13P-13 and G13S-13 alongside those from Salpeter (1955), Kroupa et al. (1993) (assumed in Gallart et al. 2024), Kroupa (2001), and Kirkpatrick et al. (2024). Vertical dashed lines delineate the mass ranges of this work, detailed in Table 2.

In the text

	Fig. B.1 Evolution of the median at each step of the 13 BGM FASt inferred parameters across the 100 ABC steps (refer to the derivation process in Sect. 2.3) for executions G13P-13 (blue) and G13S-13 (green). At the end of each plot, we present the final most probable values (MPVs) of these parameters, along with their 16th and 84th percentiles as determined in Sect. 3.4.
In the text

Fig. C.1 Corner plot containing the parameters derived from our fiducial execution G13P-13, including two slopes of the IMF and 11 parameters for the SFH of the thin disc. At the top of the columns, the projected approximate posterior PDF of each parameter with its Gaussian fit is shown, as well as the resulting most probable value and 16th and 84th percentiles, which are also marked with solid and dashed lines, respectively. It is indicated with a dashed line the median of the distributions. In black at the top right of each plot we show Pearson’s correlation coefficient. High values of this coefficient are highlighted with an intense blue frame in the plot. Finally, the parameters adopted by the MS are marked with a magenta cross and the final values with an orange star.

In the text

	Fig. C.2 Same as in Fig. C.1 but for execution G13S-13.
In the text

Fig. D.1 Top: from left to right-hand side, columns show CMDs of the PARSEC MS used in execution G13P-13, Gaia DR3 G < 13 data, distance per pixel between them (each of the elements of the summation in Eq. (4)), and their absolute difference in number of stars. From the top to the bottom rows, we find the CMDs for the three considered latitude ranges 30 < |b| < 90, 10 < |b| < 30, 0 < |b| < 10. Bottom: same as on top but for the final BGM FASt PS obtained from the application of the fitting process described in Sect. 2.1 on PARSEC MS.

In the text

	Fig. D.2 Same as in Fig. D.1 but for the STAREVOL MS used in execution G13S-13.
In the text