Gaia DR2 reveals a star formation burst in the disc 2-3 Gyr ago

TL;DR

Using Gaia DR2 data and advanced modeling, we identify a significant star formation burst in the Galactic disc 2-3 Gyr ago, revealing complex star formation history and IMF characteristics.

Contribution

First non-parametric star formation history analysis for the Galactic thin disc using Gaia DR2 data within the Besancon Galaxy Model framework.

Findings

Detected a star formation burst 2-3 Gyr ago in the thin disc.

Observed a decreasing star formation rate from 9-10 Gyr to 6-7 Gyr ago.

Estimated that about 50% of stellar mass was formed between 5 and 1 Gyr ago.

Abstract

We use Gaia DR2 magnitudes, colours and parallaxes for stars with G<12 to explore a 15-dimensional space that includes simultaneously the initial mass function (IMF) and a non-parametric star formation history (SFH) for the Galactic disc. This inference is performed by combining the Besancon Galaxy Model fast approximate simulations (BGM FASt) and an approximate Bayesian computation algorithm. We find in Gaia DR2 data an imprint of a star formation burst 2-3 Gyr ago, in the Galactic thin disc domain, and a present star formation rate (SFR) of about 1 Msun. Our results show a decreasing trend of the SFR from 9-10 Gyr to 6-7 Gyr ago. This is consistent with the cosmological star formation quenching observed at redshifts z<1.8. This decreasing trend is followed by a SFR enhancement starting at about 5 Gyr ago and continuing until about 1 Gyr ago which is detected with high statistical significance by discarding the null hypothesis of an exponential SFH with a p-value=0.002. We estimate, from our best fit model, that about the 50% of the mass used to generate stars, along the thin disc life, was expended in the period from 5 to 1 Gyr ago. The timescale and the amount of stellar mass generated during the SFR enhancement event lead us to hypothesise that its origin, currently under investigation, is not intrinsic to the disc. Thus, an external perturbation is needed for its explanation. Additionally we find, for the thin disc, a slope of the IMF of $\alpha_3 \approx 2$ for masses M>1.53 Msun and $\alpha_2 \approx 1.3$ for the mass range between 0,5 and 1,53 Msun. This is the first time that we consider a non-parametric SFH for the thin disc in the Besancon Galaxy Model. This new step, together with the capabilities of the Gaia DR2 parallaxes to break degeneracies between different stellar populations, allow us to better constrain the SFH and the IMF.

Tables

Table 1. Table 1: Age intervals, priors, and posterior PDFs (see Sects. 4 and 5 ) for the 15 inferred parameters for our fiducial case (see text).

parameter	units	age (Gyr)	${\mu }_S$	${\sigma }_S$	posterior
${\Sigma }_{\odot }^1$	$M_{\odot }p{c}^{-2}$	0-0.1	0.17	0.5	${0.16}_{-0.04}^{+0.07}$
${\Sigma }_{\odot }^2$	$M_{\odot }p{c}^{-2}$	0.1-1	1.62	3.0	${2.2}_{-0.5}^{+0.4}$
${\Sigma }_{\odot }^3$	$M_{\odot }p{c}^{-2}$	1-2	2.07	3.0	${6.5}_{-1.4}^{+1.4}$
${\Sigma }_{\odot }^4$	$M_{\odot }p{c}^{-2}$	2-3	2.39	3.0	${8.7}_{-2.1}^{+2.7}$
${\Sigma }_{\odot }^5$	$M_{\odot }p{c}^{-2}$	3-5	5.92	6.0	${12.1}_{-4.5}^{+4.1}$
${\Sigma }_{\odot }^6$	$M_{\odot }p{c}^{-2}$	5-7	7.86	8.0	${7.7}_{-4.1}^{+5.7}$
${\Sigma }_{\odot }^7$	$M_{\odot }p{c}^{-2}$	7-8	5.62	6.0	${5.8}_{-2.9}^{+7.2}$
${\Sigma }_{\odot }^8$	$M_{\odot }p{c}^{-2}$	8-9	5.62	6.0	${7.4}_{-4.4}^{+5.7}$
${\Sigma }_{\odot }^9$	$M_{\odot }p{c}^{-2}$	9-10	5.62	6.0	${9.8}_{-5.3}^{+5.9}$
${\rho }_{\odot }^{young}\cdot {10}^{-3}$	$M_{\odot }p{c}^{-3}$	$\approx 10$	$3.6$	$3.6$	${2.6}_{-0.4}^{+0.7}$
${\rho }_{\odot }^{old}\cdot {10}^{-3}$	$M_{\odot }p{c}^{-3}$	$\approx 12$	$0.5$	0.5	${0.6}_{-0.3}^{+0.7}$
${\alpha }_1$	-	all	0.5	0.5	$-{0.5}_{-0.5}^{+0.8}$
${\alpha }_2$	-	all	2.1	0.5	${1.3}_{-0.3}^{+0.3}$
${\alpha }_3$	-	all	2.9	0.5	${1.9}_{-0.1}^{+0.2}$
$h_R$	pc	0.10-10	2151	274	${1943}_{-370}^{+190}$

Table 2. Table 2: Summary of the considered model variants, the adopted SFH, extinction map, data used for the fitting and Poissonian distance ( δ P subscript 𝛿 𝑃 \delta_{P} ).

model	SFH	extinction	fitted to	${\delta }_P$
MP-S	exponential	Stilism	Tycho-2 V¡11	$7.5\cdot {10}^5$
G12Exp-S	exponential	Stilism	Gaia DR2 G¡12	$7.4\cdot {10}^5$
G12NP-S	non-param.	Stilism	Gaia DR2 G¡12	$5.6\cdot {\mathrm{𝟏𝟎}}^{\mathrm{𝟓}}$
G12NP-D	non-param.	Drimmel	Gaia DR2 G¡12	$5.4\cdot {10}^5$
G12NP-M	non-param.	Marshall	Gaia DR2 G¡12	$6.0\cdot {10}^5$

Full text

1 Introduction

The star formation history (SFH) of the Milky Way disc contains essential information to understand the Galactic structure and evolution, including key information of its merger history (e.g. Gilmore 2001). Recently, Antoja et al. (2018) discovered, using Gaia data, that an external interaction perturbed the Galactic disc in the last billion years. Moreover, Helmi et al. (2018) suggested that a merger led to the formation of the thick disc. Furthermore, from the cosmological simulations in the framework of $\Lambda CDM$, it is known that the probability that a Milky Way-like Galaxy had a minor merger in the last 10 Gyr is high (e.g. Stewart et al. 2008). These mergers can trigger stellar formation that we expect to detect in the observational catalogues when characterising the SFH of the Galactic disc (e.g. Kruijssen et al. 2018). The analysis of the Milky Way SFH cannot be disentangled from the study of the stellar initial mass function (IMF), as discussed in Haywood et al. (1997) and Aumer & Binney (2009), for example. In this context, the unprecedented accuracy of the Gaia data-release 2 (DR2) data (Gaia Collaboration et al., 2016, 2018) represents a great opportunity to search for hints of star formation bursts in the Galaxy using the population synthesis Besançon Galaxy Model (BGM; Robin et al. 2003). Previous studies performed with BGM used catalogues of colours and apparent magnitudes to perform parameter inference (e.g. Robin et al. 2014, Mor et al. 2018), carrying some degeneracies mostly due to the lack of information of the intrinsic luminosity of the stars. Now, for the first time, Gaia parallaxes help us to break some of the degeneracies between different stellar populations for a large stellar sample. Following the approach proposed in Mor et al. (2018) here we compare synthetic versus observed full-sky magnitude-limited stellar samples by using a Bayesian approach to simultaneously explore a non-parametric SFH, a three truncated power-law IMF, and the disc density laws. This is the first time that, using BGM, the SFH of the thin disc is considered non-parametric. In practice this means that we infer the surface star formation rate (SFR) of nine age bins from 0 to 10 Gyr. We summarise our method in Sect. 2. In Sect. 3 we present the observational sample used and in Sect. 4 we discuss the analysis of the data. The resulting SFH and IMF are presented and discussed in Sect. 5. Finally, in Sect. 6 we present the conclusion.

2 BGM FASt for Gaia

We use an approximate Bayesian computation (ABC) algorithm (Jennings & Madigan, 2017) together with the Besançon Galaxy model fast approximate simulations (BGM FASt, Mor et al. 2018) to infer a parameter space with 15 dimensions. BGM FASt is an analytical framework to perform very fast Milky Way simulations based on BGM. The theory of BGM FASt and the basis of the parameter inference strategy that we use in this work is extensively described in Mor et al. (2018). Summarising, our iterative parameter inference strategy works as follows. First we sample a set of 15 parameters from the prior probability distribution functions (PDFs); we choose these to be wide Gaussians centred on the results of Mor et al. (2018) (see their Figs. 6 and 7). Subsequently, we perform a new BGM FASt simulation using the sampled parameters as inputs. We then define $M_{\varpi}$ as a combination of Gaia observables: $M_{\varpi}=G+5\cdot log_{10}(\varpi/1000)+5$, where $\varpi$ is the parallax of the star. If the parallax accuracy and the interstellar absorption go to zero ($\sigma_{\varpi}\to 0$ and $A_{G}\to 0$), the $M_{\varpi}$ becomes the absolute magnitude of the star. We then use the Poissonian distance metric111We know from Kendall & Stuart (1973) and Bienaymé et al. (1987) that the Poissonian distance is a good choice for the comparison. For simplicity it can be understood as a goodness-of-fit: the shorter the Poissonian distance, the closer the simulation to the observations. ($\delta_{P}$; Eq. 58 from Mor et al. 2018) to compare synthetic versus Gaia DR2 $M_{\varpi}$-colour ($G_{Bp}-G_{Rp}$) distributions for the whole sky divided into three latitude ranges ($|b|<10$, $10<|b|<30$ and $30<|b|<90$). If the resulting Poissonian distance is smaller than an imposed threshold, the given set of 15 parameters is accepted as part of the posterior PDF. Otherwise it is rejected. We set the threshold to be small enough to ensure that we discard all the combinations of parameters that give worse results than the best model in Mor et al. (2018).

We infer the 15 parameters listed in Table 1 which are the following: the thin disc radial scale length ($h_{R}$) for populations older than 0.10 Gyr; the three slopes of a three truncated power-law IMF, $\alpha_{1}$ (for $0.09M_{\odot}<M<0.5M_{\odot}$), $\alpha_{2}$ (for $0.5M_{\odot}<M<1.53M_{\odot}$), and $\alpha_{3}$ (for $1.53M_{\odot}<M<120M_{\odot}$); the present volume stellar mass density of the thick disc at the position of the Sun for the BGM young ($\rho^{young}_{\odot}$) and old ($\rho^{old}_{\odot}$) components of the thick disc (Robin et al., 2014); and the surface stellar mass density at the position of the Sun ($\Sigma^{j}_{\odot}$) of the generated stars along the life of the thin disc, for nine intervals of age. These nine values of $\Sigma^{j}_{\odot}$ divided by the interval of age become the mean surface SFR per age bin ($M_{\odot}Gyr^{-1}pc^{-2}$). All of them together constitute a non-parametric SFH. Each complete and robust inference of the full set of parameters requires $2\cdot 10^{4}$ CPU hours in the Spark environment of the Big Data platform at the University of Barcelona. We used more than $10^{5}$ hours of CPU.

The fixed model ingredients are described in Mor et al. (2018) following Robin et al. (2003), Robin et al. (2012), and Czekaj et al. (2014). We adopt the photometric transformation of Evans et al. (2018) to transform the simulated data from Johnson to Gaia bands. The error modelling of astrometric and photometric data and the angular resolution of the stellar multiple systems (0.04 arcsec) are chosen accordingly to Gaia Collaboration et al. (2018). We define as our fiducial case the one that uses a non-parametric SFH and the Stilism extinction map (Lallement et al., 2018). This is the most recent extinction map specifically developed to be used in BGM. The prior PDFs adopted for our fiducial case are shown in Table 1.

3 Gaia DR2 observational sample

We use, from Gaia DR2, the G mean magnitude, the colours $(G_{Bp}-G_{Rp}),$ and the parallaxes ($\varpi$) for a full-sky sample limited to stars with magnitude G¡12. The completeness of the sample is estimated using the pre-cross-match of Gaia DR2 with Tycho-2 catalogue from the Gaia archive. First, we take all stars in this cross-match with $V_{T}<11$, where Tycho-2 is $99\%$ complete (Høg et al., 2000). We then compare the obtained number of stars with the number of stars in Tycho-2 with the same magnitude limit. The results obtained show that the cross-match has $\sim 2\%$ less stars than Tycho-2. Additionally, Gaia DR2 is known to be complete from G=12 to G=17 (Gaia Collaboration et al., 2018). Therefore, it is plausible to assume that, for G¡12, the catalogue is more complete closer to G=12. As a consequence, we expect that for $G<12$ the catalogue is at least as complete as for $V_{T}<11$. Therefore, we estimate that the Gaia DR2 catalogue is about $97\%$ complete up to $G=12$. Additionally, we feel it necessary to mention that about $1\%$ of the data have either no colours or have no parallax. We also limit our model-versus-data comparison in the colour range where the photometric transformation of Evans et al. (2018) is valid; this is the range of colour ($G_{Bp}-G_{Rp}$) from $-0.47$ to $2.73$. To avoid the white and brown dwarfs, which for the moment are treated independently of the thin disc in BGM, we consider only stars with $M_{\varpi}<10$. The total number of stars in the Gaia DR2 subsample used is 2890208.

4 Analysis of the data

In Table 2 we present the model variants considered here and the resulting Poissonian distance for each one. The MP-S was obtained from a fit to Tycho-2 photometry and its main parameters are reported in Fig. 7 of Mor et al. (2018). The remaining model variants in the present work result from fitting Gaia DR2 data using photometric data and parallaxes. As shown in this table, the Poissonian distance of the models that use a non-parametric SFH is smaller than that of the models with an exponential SFH. We want to emphasise the fact that the difference in the Poissonian distance between G12NP-S (our fiducial case; see Sect. 2) and G12NP-D is not large enough to settle on which extinction map is better.

In Fig. 1 we show a density map, $M_{\varpi},$ as a function of the Gaia colour $(G_{Bp}-G_{Rp})$ for the three latitude ranges considered. For the stars with $10<M_{\varpi}<-1$ we set the bin size to 0.05 mag in colour and 0.25 mag in $M_{\varpi}$. For stars with $M_{\varpi}<-1$ the bin size is enlarged to 0.5 mag in colour and 1 mag in $M_{\varpi}$ to allow a robust statistical analysis; these stars, which represent $7\%$ of the sample, are not shown in the figure. Even though $M_{\varpi}$ is not strictly the absolute magnitude, we refer to this density map as if it were a true Hertzsprung-Russel diagram. The first column shows the Gaia DR2 data, the second column shows the MP-S model variant, and the third column shows the best-fit model obtained in this work for our fiducial case (G12NP-S). The fourth and fifth columns show the differences, in star counts per bin, between our Gaia DR2 sub-sample and both the old MP-S and the new G12NP-S. With this new fit we improve the agreement of the model with the data for the three latitude ranges. Focusing on the bright end ($M_{\varpi}<-1$), the improvement is significant mostly in the Galactic plane.

The first feature of Fig. 1 that we want to comment on is the following. In the Gaia DR2 data, we can see a blob of stars (mostly in blue) below the main sequence (MS) that is not reproduced in the simulations. These are stars that are flagged with a bad colour-excess (Evans et al., 2018) in the Gaia catalogue and represent only 0.1% of our sample. By comparing the fourth and fifth columns of Fig. 1, we can see how the agreement with Gaia data is better when using the non-parametric SFH (G12NP-S). Both the excess of stars detected in the MP-S around the MS region and the deficit of stars in the region with $0.5<(G_{Bp}-G_{Rp})<1.0$ and $1.5<M_{\varpi}<3$ are clearly diminished in the new G12NP-S.

The high quality of the new information given by Gaia data reveals new discrepancies between model and data. Most of these discrepancies could come from several assumptions on the fixed ingredients of the BGM model that, as largely discussed in Sect. 7.3 of Mor et al. (2018), can impact our parameter inference. Here we discuss the discrepancies that we see in the fifth column of Fig. 1, grouped into three main areas: the red giant branch (RGB), the MS, and the region of stars with $3<M_{\varpi}<4$ and $0.5<(G_{Bp}-G_{Rp})<1,$ hereafter referred to as the square region. In the RGB region, we notice that in the simulations the position of the RGB clump is shifted and that it is less extended than in the observations. In this region, the asymptotic giant branch bump is less dense in the simulations; this effect is mostly seen at intermediate latitudes. In the MS region, we detect a clear sequence where the simulation has an excess of stars and, immediately above, a clear sequence with a deficit of stars. We also notice that in the square region the simulation has a deficit of stars. This last effect is stronger at high and intermediate latitudes. As expected, these differences are caused by a mixture of factors. From a first analysis of our simulations we conclude that the thick disc modelling and the stellar evolutionary models are the two main ingredients causing the discrepancies. There are other ingredients that can contribute to these discrepancies however: the assumed colour transformation; the radial metallicity distribution; the age-metallicity relation; and the atmosphere models, which mostly have their impact in the MS and RGB; the rate of mass loss assumed in the stellar evolutionary models; the assumption of the extinction map, mostly affecting the RGB; and finally, the assumed resolution of the stellar multiple systems mostly affecting the MS. Work is in progress to more thoroughly analyse these discrepancies, to confirm their causes, and to improve the BGM model accordingly.

5 The resulting SFH and IMF

In the last column of Table 1 we show the results of the 15 inferred parameters for our fiducial case. In this section we focus on the discussion of the resulting IMF 222As largely discussed in Mor et al. (2017), the IMF considered in BGM is a composite IMF (or Integrated Galactic IMF; IGIMF). and the non-parametric SFH of the thin disc. In Fig. 2 we present the nine values of the local mean surface SFR as a function of age that constitute the SFH. In this figure we show the results for our fiducial case (G12NP-S). Additionally, to evaluate the impact of the choice of the extinction model in our results, we present the inferred SFH when using the Drimmel & Spergel (2001) (G12NP-D) and Marshall et al. (2006) 333The Marshall et al. (2006) extinction map covers the longitude ranges $-100<l<100$ and the latitude ranges $|b|<10$. Therefore, Drimmel map is used for the rest of the sky. (G12NP-M) extinction maps (see Table 2). We notice that the differences in the SFR among them are not larger than $1.5\sigma$. Regardless of the choice of the extinction map we can see a general decreasing trend from 9-10 Gyr to 6-7 Gyr ago followed by a SFR enhancement event beginning about 5 Gyr ago. This enhancement event is of about 4 Gyr in duration with a maximum at about 2-3 Gyr ago and with a final decreasing trend until the present time. We would like to point out that our results do not rule out an earlier beginning for this SFR enhancement event, nor a constant (or slightly decreasing) SFR with a value of about $7M_{\odot}\cdot Gyr^{-1}pc^{-2}$ from 10 Gyr ago until 1 Gyr ago, with a very sharp and fast drop in the last 1 Gyr. We estimate from our best-fit model that about $50\%$ of the mass used to generate stars throughout the life of the thin disc was expended in the period from 5 to 1 Gyr ago.

To evaluate the statistical significance of the enhancement event we compare the values of the non-parametric SFH of our fiducial case (G12NP-S) with both: (1) the values of the SFH obtained when imposing an exponential SFH in BGM and performing the fit with Gaia data (G12Exp-S results) and (2) the values of an exponential shape fitted to the G12NP-S results using the least squares method (grey dashed line in Fig. 2). This fitted exponential shape is purely mathematical and does not necessary make physical sense. The statistical significance of the points in the SFR enhancement event for the first and second cases are the following: $2.8\sigma$ and $2.5\sigma$ for the point at 2.5 Gyr ago (the relative maximum); $3\sigma$ and $1.5\sigma$ for the point at 1.5 Gyr ago, and finally $1.5\sigma$ and $0.8\sigma$ for the point at 4 Gyr ago. To evaluate the significance of the event as a whole we perform two tests with the following two null hypotheses: (1) the SFH is exponential and follows the results of the G12Exp-S variant and (2) the SFH follows the fitted mathematical exponential shape (grey dashed line in Fig. 2). In a subsequent step, for both tests, we assume that the distribution of the SFR at each age bin follows a Gaussian centred in the values given by the null hypothesis and with $\sigma$ estimated from the obtained posterior PDF. Afterwards, for both null hypotheses, we compute the p-value for the points at 1.5, 2.5, and 4 Gyr ago. We finally obtain a p-value of the global event by combining the individual p-values using Fisher’s Method. The results give p-values of ¡0.001 and =0.002, for hypotheses (1) and (2), respectively; we therefore reject both null hypotheses.

To mathematically characterise the SFR enhancement event we fit a bounded exponential plus a Gaussian function to the results (black dashed line in Fig. 2), obtaining for the Gaussian component $\mu=2.57Gyr$ and $\sigma=1.25Gyr$. In Fig. 2 we additionally show the exponential part of this last fit (grey solid line) where we see how the SFH of our fiducial case follows an exponential shape from 10 Gyr until 6-7 Gyr ago. From all the performed tests we conclude that the SFR enhancement that we find is statistically significant.

Our findings that the thin disc SFH does not follow a simple decreasing shape until the present are in good agreement with Snaith et al. (2015), Haywood et al. (2016), and Haywood et al. (2018) who found, using data with metallicities and assuming a fixed IMF, the existence of an SFR quenching followed by a reactivation. Kroupa (2002b), using stellar kinematics, found the SFH to behave similarly. The relative maximum of the SFR that we find at 2-3 Gyr ago is compatible with the results of Vergely et al. (2002) and Cignoni et al. (2006) that, using Hipparcos data in a sphere of 80 pc around the Sun and assuming a fixed IMF, found maximum peaks at 1.75-2 Gyr ago and 2-3 Gyr ago, respectively. Recently, Bernard (2018), in a tentative work using TGAS data, pointed towards the existence of a relative maximum also located 2-3 Gyr ago.

In Fig. 3 we show the resulting slopes $(\alpha)$ of the inferred IMF as a function of stellar mass and a compilation of results in the literature. For the mass range between $0.5M_{\odot}$ and $1.53M_{\odot}$ we find $\alpha_{2}=1.3\pm 0.3$, in very good agreement with Rybizki & Just (2015) who found $\alpha=1.49\pm 0.08$ (in the range $0.5M_{\odot}<M<1.4M_{\odot}$). For masses larger than $1.53M_{\odot}$ we find $\alpha_{3}=1.9_{-0.1}^{+0.2}$, which is flatter than the $\alpha_{3}$ obtained by Salpeter (1955) and Kroupa (2002a). For the low-mass range ($0.09M_{\odot}<M<0.5M_{\odot}$) we obtain values between $\alpha_{1}=-1$ and $\alpha_{1}=0.5$. We must keep in mind our estimation that about $99.6\%$ of the stars in our sample have masses between $0.5M_{\odot}$ and $10M_{\odot}$, with only $0.1\%$ of the stars in our sample belonging to the lowest mass range, and that of the order of $10^{4}$ stars have $M>10M_{\odot}$. We also want to compare our results with two works that consider a non-universal IMF. These are the recent works of Dib & Basu (2018) and Jeřábková et al. (2018) (see Appendix A). We note that, as in our results, their IMFs have a shallower shape than the values of Kroupa or Salpeter. The information from Gaia parallaxes, when imposing an exponential SFH, brings the resulting $\alpha_{3}$ to be $\approx 2.5$. This is flatter than in our previous works (Mor et al., 2017, 2018) but compatible with the $\alpha_{3}\approx 2.7$ of the IGIMF (e.g. Kroupa et al. 2013). When we adopt a non-parametric SFH we find $\alpha_{3}\approx 2$, more in the direction of Zonoozi et al. (2019). We know from Mor et al. (2018) (e.g. their Fig. 7) that when characterising the IGIMF from star counts the correlation of $\alpha_{3}$ with the SFH is high. In our previous works, the fact that we were imposing an exponential SFH resulted in a steeper $\alpha_{3}$. These correlations between the $\alpha_{3}$ and the SFH are also observed in our present work. We find that the $\alpha_{3}$ is anti-correlated with the four $\Sigma^{j}_{\odot}$ values for the age bins from 0.1 to 5 Gyr, with a Pearson’s correlation coefficient going from about $-0.8$ to about $-0.5$. We also note that these four surface densities ($\Sigma^{j}_{\odot}$) are also correlated among them, with coefficients from about 0.3 to 0.5. We want to emphasise that the effects of these correlations in the results (Table 1, Figs. 2 and 3) are already taken into account when we provide the 0.16 and 0.84 quantiles of the posterior PDFs. From the surface SFR of the youngest population ($1.6_{-0.4}^{+0.7}M_{\odot}Gyr^{-1}pc^{-2}$) we find the present SFR in the disc to be about $1M_{\odot}/yr$. This result is very sensitive to the disc scaling of the youngest population. We also find a radial scale length of $h_{R}=1943^{+190}_{-370}$ compatible with Robin et al. (2012). For robustness we repeated our analysis by adding to the Gaia DR2 parallaxes the offset of 0.029 mas reported in Lindegren et al. (2018), concluding that the impact on the derived IMF and SFH is much smaller than the impact of the choice of the extinction map.

6 Conclusion

For the first time we have considered a non-parametric SFH for the thin disc in the BGM model. This new step, together with the capability of the Gaia DR2 parallaxes to break degeneracies between different stellar populations, allowed us to better constrain the thin disc SFH and IMF. The resulting SFH shows a decreasing trend from 9-10 to 6-7 Gyr ago that is consistent with the quenching observed in a cosmological context for redshifts $z<1.8$ (e.g. Rowan-Robinson et al. 2016) and is also compatible with the evidence of the quenching in the Milky Way reported in Haywood et al. (2018). The quenching that we find could be linked to a previous merger event. Simulations in the framework of $\Lambda CDM$ show that after a merger, there is an enhancement of the star formation followed by a quenching (e.g. Di Matteo et al. 2008, Fig. 4). This would be compatible with the thick-disc formation scenario recently proposed in Helmi et al. (2018) with a merger which occurred more than 10 Gyr ago. As suggested in other works, the quenching that we find could also be partially produced by the presence of a galactic bar (Haywood et al., 2016; Khoperskov et al., 2018). The two quenching mechanisms discussed here are not mutually exclusive but complementary. After the quenching, we detect a 4 Gyr duration SFR-enhancement event starting at about 5 Gyr ago and with a maximum at 2-3 Gyr ago. The large timescale of this recent SFR enhancement event, together with the large amount of mass that we estimate to be involved in it (see Sect. 5), lead us to propose that this recent event is not intrinsic to the disc but is produced by an external perturbation. Furthermore, the slow increase of the star formation process, its duration, as well as the high absolute value of the maximum suggest that this could be produced by a recent merger with a gas-rich satellite galaxy that could have started between about 5 and 7 Gyr ago. However, an analysis of other stellar parameters (e.g. metallicities) would be needed to favour this hypothesis over other possible scenarios.

Work is in progress to more thoroughly analyse the Gaia DR2 data by extending our study to fainter magnitudes, updating the stellar evolutionary models and the thick disc modelling.

Acknowledgements.

We thank the anonymous referee for the constructive report which helped to improve the quality of the work. Special thanks are to N. Lagarde, C. Charbonnel, and C. Reylé for very useful discussions. This work was supported by the MINECO (Spanish Ministry of Economy) - FEDER through grant ESP2014-55996-C2-1-R and MDM-2014-0369 of ICCUB (Unidad de Excelencia ’María de Maeztu’), the French Agence Nationale de la Recherche under contract ANR-2010-BLAN-0508-01OTP and the European Community’s Seventh Framework Programme (FP7/2007-2013) under grant agreement GENIUS FP7 - 606740. We also acknowledge the team of engineers (GaiaUB-ICCUB) in charge to set up the Big Data platform (GDAF) at University of Barcelona. We also acknowledge the International Space Science Institute, Bern, Switzerland for providing financial support and meeting facilities. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Appendix A Treatment of the IMFs from Dib & Basu (2018) and Jeřábková et al. (2018)

To be able to compare the IMFs of Dib & Basu (2018) and Jeřábková et al. (2018) with the IMF that we obtain in this work we need to perform an adequate treatment. For the case of Dib & Basu (2018), to be able to compare the slopes, we fit a three truncated power-law IMF to their results when they assume $a_{\Gamma}=0.5$, $a_{\gamma}$=0.5, and $a_{M_{ch}}=0.5$ (see their Fig. 1). We plot this result in Fig. 3. The case of Jeřábková et al. (2018) is slightly more complex as their IMF depends on the SFR and metallicity. From our results we estimate a mean SFR (in $M_{\odot}/yr$) and a mean metallicity for each one of the age bins considered. Then from Yan et al. (2017) and Jeřábková et al. (2018) we get a resulting IMF for each age bin (using galIMF; https://github.com/Azeret/galIMF). Finally we compute a weighted mean of the IMF depending on the total mass for each age bin. The total integrated galactic IMF is then plotted in Fig. 3.