2.2.2 Stars in the Galaxy

The stellar population synthesis model of the Galaxy constructed in Besançon is able to simulate the stellar content of the Galaxy by modelling five distinct stellar populations: the thin disc, two thick discs (a young and a old), the bar and the stellar halo. It can be used to generate stellar catalogues for any given direction, and returns characteristics of each simulated star such as magnitude, colour, and distance as well as kinematics and other stellar parameters.

The approach of the Galactic model is semi-empirical as it is based on theoretical considerations (stellar evolution, Galactic evolution and Galactic dynamics) but is constrained by observations (the local luminosity function, the age-velocity dispersion relation, the age-metallicity relation). The Galactic potential is calculated in order to self-consistently constrain the disc scale height, the thin disc being subdivided into 7 isothermal components of ages varying from 0-0.15 Gyr for the youngest to 7-10 Gyr for the oldest. For computing the scale height at the solar position as a function of age (Bienayme et al. 1987b, a), the Boltzmann equation (first moment at the first order with the plane parallel approximation) is used assuming an age-velocity dispersion relation (AVR) obtained by Robin et al. (2017) by fitting it to Gaia DR1 proper motions and RAVE DR4 radial velocities (Kordopatis et al. 2013).

The distribution in the Hess diagram split into several age bins is obtained from an evolutionary model which starts with a mass of gas, generates stars of different masses assuming an Initial Mass Function (IMF) and a star formation rate history (SFH), and makes these stars evolve along evolutionary tracks. The original evolution model is described in Haywood et al. (1997a, b). It has been revised in Czekaj et al. (2014) and in Mor et al. (2018). The set of tracks are described in Czekaj et al. (2014), based on Padova tracks for stars of mass larger than 0.7 ${ℳ}_{\odot }$ (Bertelli et al. 2008), and Lyon’s tracks for less massive stars (Chabrier and Baraffe 1997). Using Gaia DR2 data, Mor et al. (2018) applied a new code for fast model computation (BGM-Fast) and a Bayesian inference scheme to derive new IMF, SFH, self-consistent eccentricities for the thin disc population, and revised scale length for stars older than 0.15 Gyr. A full description of the resulting parameters are given in the paper.

Based on these inputs (tracks, IMF, SFH) the evolutionary model is used to compute the distribution of stars per element volume in the space (M${}_V$, $T_{\mathrm{eff}}$, Age). This distribution is computed from a simulation of the local sphere, selecting only single stars and primaries. Other components of stellar systems are generated according to empirical multiplicity probabilities (see below).

The Mor et al. (2018) fitting on Gaia DR2 data was performed only for the thin disc evolutionary parameters. For other populations, a single burst population is simulated taking a single isochrone. The bar isochrone is from Padova stellar models Bertelli et al. (2008) for a metallicity z=0.017, helium=0.26, log(age) of 9.9 and IMF slope $\alpha =-0.5$. The young thick disc isochrone is from Bergbusch and Vandenberg (1992) for a metallicity of [Fe/H]=-0.47, alpha enhanced with [O/Fe]=0.23, helium proportion of Y=0.245, age of 10 Gyr and IMF slope of -1. The old thick disc isochrone is from Bergbusch and Vandenberg (1992) for a metallicity of [Fe/H]=-0.78, alpha enhanced with [O/Fe]=0.39, helium proportion of Y=0.239, age of 12 Gyr and IMF slope of -1. The halo isochrone is from Bergbusch and Vandenberg (1992) for a metallicity of [Fe/H]=-1.48, alpha enhanced with [O/Fe]=0.6, helium proportion of Y=0.235, age of 14 Gyr and IMF slope $\alpha$ of -0.5.

Characteristics of evolution of each population are given in Table 2.1.

The density law of each population are given in Table 2.2. In this table for simplicity the density formulae do not include the warp and flare, which are added as a modification of the position and thickness of the mid-plane. There is no more a cut-off of the thin disc with distance. The Sun-Galactocentre distance is assumed to be 8 kpc, and the distance from the Sun to the plane 15 pc.

The flare is the increase of the thickness of the disc with Galactocentric distances, while the warp is the displacement of the mean plane depending on azimuth. The warp shape depends on age.

The warp and flare are modelled following Amôres et al. (2017). They are not symmetrical with regards to the Galactic plane and amplitude and positions depend on the age of each thin disc component. However, the warp amplitudes for old disc (age$>$5 Gyr) used by Amôres et al. (2017) are modified to adopt the value at age of 4 Gyr, because those values suffered of high uncertainties (see their figure 7) and produced a too strong warp amplitude.

Warp parameters as a function of age are given in Table 2.3. The displacement $z_{warp}$ of the mid-plane is computed following the equations:

with ${\theta }_w=atan2(x,-y)$ and $\gamma$ being ${\gamma }_{+}$ or ${\gamma }_{-}$ depending on which side of the Galaxy.

It applies when the Galactocentric radii is larger than $R_w$.

The flare is simulated by multiplying the scale height of a population by a factor $k_{flare}$ for Galactocentric radii larger than $R_{flare}$. The amplitude of the flare ${\gamma }_{flare}$ is defined as:

$k_{flare}=1.+g_{flare}\times (R-R_{flare})$

with $g_{flare}$=0.0001 and $R_{flare}=$ 10 kpc.

In the thin disc population a spiral structure has been added with 2 arms, as determined in a preliminary study by Amores & Robin (priv. comm.). The parameters of the arms are given in Table 2.5. It applies only on the young thin disc of age less than 0.15 Gyr.

The thick disc scale height, scale length and normalisation, and halo parameters are based on (Robin et al. 2014), i.e. fitted to 2MASS and SDSS photometric surveys. The thick disc covers 2 sub-populations with different ages and metallicities. The young thick disc is 10 Gyr old with a metallicity of -0.48 dex, while the old thick disc is 12 Gyr old with a mean metallicity of -0.78 dex. The new values of the scale height, scale length of the discs, and eccentricity and shape of the halo are given in Table 2.2.

In order to compute the $N$ number of stars at any point in the Galaxy of a given population the equation of stellar statistics is used, based on density laws and the distributions in Mv, $T_{\mathrm{eff}}$ and age described above, as follows :

where $\rho$ is the density law of the population, $r$ the distance to the sun, (R,$\theta$,z) the Galactocentric cylindrical coordinates, $\omega$ the solid angle, $N$ is the theoretical number of stars in a volume element with the intrinsic parameters M${}_V$, ${\mathrm{T}}_{\mathrm{eff}}$, and Age. In order to simulate catalogues, from this number a random drawing is performed to produce an integer number of stars, which follows a Poisson distribution of expectation (and variance) $N$.

Metallicity [Fe/H] of the thin disc is computed through an empirical age-metallicity relation from Haywood (2008). For each age and population component the metallicity is drawn from a Gaussian, with a mean and a dispersion, as given in Table 2.6. One also accounts for radial metallicity gradient for the thin disc, -0.7 dex/kpc, but no vertical metallicity gradient, as the gradient is naturally coming from the vertical variation of populations. Alpha elements are computed as a function of the population and the metallicity. For the bar, thick disc and halo, the abundance in alpha element is drawn from a random around a mean value which depends on [Fe/H], truncated at 5 sigma. Formulae are given in Table 2.6. For the thin disc, the formula is a fit to the mean location of thin disc stars from Bensby et al. (2010) study. For the thick disc and the bar the formulae were fitted on APOGEE data in SDSS-DR14 (Abolfathi et al. 2018).

The kinematics is based on a dynamically self-consistent model, following Bienaymé et al. (2018), where the third integral of motion is obtained by approximating the BGM Galactic potential with a Staeckel potential. The velocity distributions in (R,z) are then tabulated to be further interpolated when computing the star kinematics. The asymmetric drift is also provided in a table, as a function of R and z and for each sub-component of the thin and thick discs. Finally, the age-velocity dispersion relations have been established by comparing simulations with Gaia DR1 (TGAS) proper motions and RAVE radial velocities (Robin et al. 2017). Values are given in Table 2.7 at the solar position. An exponential radial gradient of -0.16E-03 $\mathrm{km}{\mathrm{s}}^{-1}$/pc on ${\sigma }_U$ and ${\sigma }_V$, and no gradient on ${\sigma }_W$ are applied.

The rotation curve is taken from Sofue (2015), giving a $V_{LSR}$ of 230.6 $\mathrm{km}{\mathrm{s}}^{-1}$ at R=8 kpc. The Sun peculiar velocities are assumed to be (12.75, 0.93, 7.10) in $\mathrm{km}{\mathrm{s}}^{-1}$ following the fit to TGAS+RAVE data (Robin et al. 2017).

Due to a mistake in the implementation, very young stars (age less than 0.15 Gyr) have wrong circular velocity (the circular motion from the rotation curve was not applied on them). Hence their circular velocity is null in average instead of being equal to the mean rotation curve. Their velocity corresponds only to their peculiar motion (assumed Gaussian with sigma according to the velocity ellipsoid). The mistake concerns about one million very young stars in the Milky Way. For the user it is possible to correct this bug by adding to the V velocity the mean rotation curve of the Galaxy for each star with age smaller than 0.15 Gyr. Notice also that these stars which are expected to be related to spiral arms or in the Gould Belt should also have other systematic motions and deviations from the circular motion, which are not either taken into account in the present simulation.

For thin disc stars the asymmetric drift is computed from tabulated (R,z,age) distributions. However this table covers the Galactocentric distance range 2$\le$R$\le$16 kpc. Outside this range the asymmetric drift needs to be extrapolated but the extrapolation is not accurate in the present version. At R2 kpc there are only a few thin disc stars, mainly because the major contribution comes from bulge and thick disc stars. At R$>$16 kpc, there are still a few thin disc stars that can show overestimated circular velocity due to this wrong asymmetric drift.

An error has also been identified on velocities of stars which are on the other side of the Galaxy (at Xgal$>$0, beyond the Galactic centre). The code used a wrong formula to compute the projected velocities to the heliocentric (u,v,w) frame. The bug concerns only stars with heliocentric distances larger than 8 kpc, a small subsample of Gaia stars.

White dwarfs are generated from tables as in the BGM. For the thin disc population, the theoretical luminosity function from Wood (1992) is used, normalised to the local density of white dwarfs estimated from Gaia DR2 by Hollands et al. (2018) in a volume of 20 pc from the Sun. The local value of ${4.4910}^{-3}$ per cubic parsec is then corrected for binarity assuming that 12% are secondaries of binary systems to have an estimate of the primary white dwarfs. In the simulator, white dwarfs are generated in secondaries as well. Only DA and DB type white dwarfs are generated in this version. For the thick disc population white dwarfs are generated similarly assuming an age of 7 to 10 Gyr, while halo white dwarfs follow Chabrier luminosity function for an age of 14 Gyr and a Salpeter IMF (Chabrier 1999, and private communication).

Colours of white dwarfs are taken from synthetic photometry provided by models of Holberg and Bergeron (2006), which allows to compute V and V-I. The computation of colours in the Gaia bands is performed using the suitable stellar libraries (Section 2.2.5).

White dwarf density laws are the same as for other spectral types of a given age.

Pre-main sequence stars are generated using evolutionary tracks from Valls-Gabaud, Pichon, et al (in prep.). They have been produced using the code of stellar evolution CESAM for masses between 0.7 and 3 solar masses. The volume density along the tracks have been computed with a constant star formation rate.

For the needs of the simulator, some rare objects have been added to the BGM Hess diagram.