2.2.2 Stars in the Galaxy
Author(s): Annie Robin
A reliable model of the stellar distribution over the sky up to magnitude 21 is most important, with reliable spectral distributions and motions for these objects. The current version of a Galaxy model has been implemented in java, based on the Besançon Galaxy Model (Robin et al. (2003), here after BGM) but with noticeable updates.
Galactic structure model inputs
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 semiempirical as it is based on theoretical considerations (stellar evolution, Galactic evolution and Galactic dynamics) but is constrained by observations (the local luminosity function, the agevelocity dispersion relation, the agemetallicity relation). The Galactic potential is calculated in order to selfconsistently constrain the disc scale height, the thin disc being subdivided into 7 isothermal components of ages varying from 00.15 Gyr for the youngest to 710 Gyr for the oldest. For computing the scale height at the solar position as a function of age (Bienayme et al. 1987, 1987), the Boltzmann equation (first moment at the first order with the plane parallel approximation) is used assuming an agevelocity 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 ${\mathcal{M}}_{\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 (BGMFast) and a bayesian inference scheme to derive new IMF, SFH, selfconsistent 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}$, ${\mathrm{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.
Age  ecc  IMF  SFH  
Disc  00.15  0.0140  
0.151  0.0268  dn/dm $\propto $ m${}^{\alpha}$  
12  0.0375  $\alpha $ = 0.5 for m $$ ${\mathcal{M}}_{\odot}$  
23  0.0555  $\alpha $ = 2.1 for 0.5 $$ ${\mathcal{M}}_{\odot}$  exp.  
35  0.0703  $\alpha $ = 2.9 for m $>1.53$ ${\mathcal{M}}_{\odot}$  decreasing  
57  0.0794  $\gamma $ = 0.14  
710  0.0799  
Stellar halo  14  0.774  dn/dm $\propto $ m${}^{0.5}$  one burst 
Age  ${h}_{z}$  IMF  SFH  
Thick disc young  10  400.  dn/dm $\propto $ m${}^{1}$  one burst 
Thick disc old  12  795.  dn/dm $\propto $ m${}^{1}$  one burst 
Bar  10  300.  dn/dm $\propto $ m${}^{0.5}$  one burst 
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 midplane. There is no more a cutoff of the thin disc with distance. The SunGalactocentre distance is assumed to be 8 kpc, and the distance from the Sun to the plane 15 pc.
Thin Disc  
age$\le 0.15$  ${\rho}_{0}/{d}_{0}/{k}_{flare}\times \{\mathrm{exp}({(a/{h}_{{R}_{+}})}^{2})\mathrm{exp}({(a/{h}_{{R}_{}})}^{2})\}$  
with ${h}_{{R}_{+}}$ = 5000 pc, ${h}_{{R}_{}}$ = 3000 pc  
age $>0.15$  ${\rho}_{0}/{d}_{0}\times \{\mathrm{exp}(0.5{(0.25+{a}^{2}/{h}_{{R}_{+}}^{2})}^{1/2})$  
$\mathrm{exp}(0.5{(0.25+{a}^{2}/{h}_{{R}_{}}^{2})}^{1/2})\}$  
with ${h}_{{R}_{+}}$ = 2011 pc, ${h}_{{R}_{}}$ = 1310 pc  
Thick disc  ${\rho}_{0}/{d}_{0}/{k}_{flare}\times \mathrm{exp}(\frac{R{R}_{\odot}}{{h}_{R}})\times \mathrm{cosh}(z/(2.\times {h}_{z}))$  
with ${h}_{R}$ = 2040 pc, ${h}_{z}$ = 400 pc  young thick disc  
with ${h}_{R}$ = 2919 pc, ${h}_{z}$ = 795 pc  old thick disc  
Spheroid  ${\rho}_{0}/{d}_{0}\times {a}^{{\alpha}_{h}}\times {(rc+a)}^{{\delta}_{h}}$  
with ${\alpha}_{h}=1.,{\delta}_{h}=2.777231$ and $rc=2175.827881$  
Bar  $N\times \mathrm{cosh}(r{s}^{2})$  
with $rPerp=[{(x/x0)}^{cPerp})+{(y/y0)}^{cPerp})]{}^{(1./cPerp)}$  
$rs={[(rPer{p}^{cPara})+{(z/z{p}_{0})}^{cPara}]}^{1./cPara}$  
$z{p}_{0}=z0\times (1.+dzbone\times \mathrm{sin}(x/xbone)$ 
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 midplane is computed following the equations:
$${z}_{warp}=(R{R}_{w})\times \mathrm{cos}(\theta {\theta}_{w})\times \gamma $$ 
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}$.
Age  ${R}_{w}$  ${\theta}_{w}$  ${\gamma}_{+}$  ${\gamma}_{}$ 
00.15  10169.875  1.80685  1.00652  0.06815 
0.151  10042.375  1.77785  0.84525  0.10915 
12  9806.5  1.72420  0.61675  0.18500 
23  9551.5  1.66620  0.47175  0.26700 
35  9169.  1.57920  0.45300  0.39000 
57  8659.  1.46320  0.45300  0.39000 
710  8021.5  1.31820  0.45300  0.39000 
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.
Parameter  Value  Unit  Comment 
x0  1076.  pc  
y0  490.  pc  
z0  300.  pc  
$\alpha $  0.509  rad  bar angle with regards to Sungalactic center axis 
$\gamma $  0  rad  bar roll angle 
$\beta $  0  rad  bar tilt angle 
D0  0.4508 $\times {10}^{11}$  normalisation  
${R}_{Max}$  4195.55  pc  cutoff radius 
cPara  3.007  boxyness in (x,y) plane  
cPerp  3.329  boxyness in (x,z) plane  
invhEnd  4.e6  pc${}^{1}$  sigma of the cutoff 
dzbone  0.534  peanut shape amplitude  
xbone  870.  pc  distance of the peak of the peanut to z axis 
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.
Parameter  1st arm  2nd arm 
Internal radius in kpc  3.426  3.426 
Pitch angle in radian  4.027  3.426 
Phase angle of start in radian  0.188  2.677 
Amplitude  1.823  2.013 
Thickness in the plane  4.804  4.964 
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 subpopulations 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 excentricity 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, Teff and age described above, as follows :
$$N={\int}_{0}^{{R}_{max}}\rho (R,\theta ,z)\times \mathrm{\Phi}({\mathrm{M}}_{V},{\mathrm{T}}_{\mathrm{eff}},\mathrm{Age})\omega {r}^{2}dr$$ 
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$.
Abundances
Metallicity [Fe/H] of the thin disc is computed through an empirical agemetallicity 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 SDSSDR14 (Abolfathi et al. 2018).
Population  $$ [Fe/H] $>$  $\sigma $  $\frac{d[\mathrm{Fe}/\mathrm{H}]}{dR}$  $$ [$\alpha $/Fe] $>$  $\sigma $ 
Disc  +0.01 to 0.12  0.10 to 0.18  0.07  0.121$\times $[Fe/H] + 0.0259 if [Fe/H]$$0.1  0.02 
0. if [Fe/H]$>$0.1  0.02  
Thick disc young  $$0.5  0.30  0.00  0.32$\mathrm{exp}(1.19375\times $[Fe/H]1.6)  0.05 
Thick disc old  $$0.8  0.30  0.00  0.32$\mathrm{exp}(1.19375\times $[Fe/H]1.6)  0.05 
Stellar halo  $$1.5  0.50  0.00  0.3  0.05 
Bar  0.20  0.20  0.00  0.32$\mathrm{exp}(1.19375\times $[Fe/H]1.6)  0.05 
Kinematics
The kinematics is based on a dynamically selfconsistent 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 subcomponent of the thin and thick discs. Finally, the agevelocity 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.16E03 km 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 km s${}^{1}$ at R=8 kpc. The Sun peculiar velocities are assumed to be (12.75, 0.93, 7.10) in km s${}^{1}$ following the fit to TGAS+RAVE data (Robin et al. 2017).
Age  ${\sigma}_{U}$  ${\sigma}_{V}$  ${\sigma}_{W}$  
(Gyr)  (km s${}^{1})$  (km s${}^{1})$  (km s${}^{1})$  
Disc  00.15  12.691  7.214  5.87 
0.151  14.528  8.258  6.72  
12  19.851  11.283  9.182  
23  24.368  13.851  11.272  
35  30.361  17.257  14.044  
57  36.889  20.968  17.063  
710  42.699  24.270  19.751  
Thick disc young  10  40.018  31.855  27.889 
Thick disc old  12  75.645  55.406  66.426 
Spheroid  14  131  106  85 
Bar  8  150  115  100 
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.
White dwarfs
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.49e3 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 VI. The computation of colours in the Gaia bands is performed using the suitable stellar librairies (Sect. 2.2.5).
White dwarf density laws are the same as for other spectral types of a given age.
Premain sequence stars
Premain sequence stars are generated using evolutionary tracks from VallsGabaud, 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.
Rare objects
For the needs of the simulator, some rare objects have been added to the BGM Hess diagram.
Be stars
Be stars are simulated as a proportion of 29% of O7B4 stars, 20% of B5B7 stars and 3% of B8B9 stars (Jaschek et al. 1988; Zorec and Briot 1997). They are attributed a parameter linked with the line strength in order to attribute them a spectrum from the spectral library (see Sect. 2.2.5). This parameter is the ratio between the envelope radius and the stellar radius ($R\mathrm{\_}env/R\mathrm{\_}star$). In order to simulate the variation of the emission lines with time this parameter is set randomly between 1.2 and 6.9 at any time.
Peculiar metallicity stars
Two types of peculiar metallicity stars are simulated, following Kurtz (1982); Kochukhov (2007); Stift and Alecian (2009).

•
Am: 12% of A stars on the main sequence in the Teff range 7400 K  10200 K are set to be Am stars. 70% of these Am stars belong to a binary system with period from 2.5 to 100 days. They are forced to be slow rotators : 67% have $$ 50 km s${}^{1}$and 23% have $$ 100 km s${}^{1}$.

•
Ap/Bp: 8% of A and B stars on the main sequence in the Teff range 8000 K  15000 K are set to be Ap or Bp stars. They are forced to have a smaller rotation : $$ 120 km s${}^{1}$.
Wolf Rayet stars
Wolf Rayet (WR) densities are computed following the observed statistics from the VIIth catalogue of Galactic WolfRayet stars of van der Hucht (2001). The local column density of WR stars is $2.9\times {10}^{6}$ pc${}^{2}$ or volume density 2.37$\times {10}^{8}$ pc${}^{3}$. Among them 50% are WN, 46% are WC and 3.6% are WO. The absolute magnitude, colours, effective temperature, gravity and mass have been estimated from the literature. The masses and effective temperature vary considerably from one author to another. As a conservative value, it is assumed that the WR stars have masses of 10 ${\mathcal{M}}_{\odot}$ in the mean, an absolute V magnitude of 4, a gravity of 0.5, and an effective temperature of 50,000K.