skip to main content

gaia data release 3 documentation

2.2 Universe Model

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 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 (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 (MV, Teff, 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 α=-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 α of -0.5.

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

Table 2.1: Age in Gyr, eccentricities (for the thin disc and stellar halo), scale height hz in pc (for the thick disc), initial mass function (IMF), and star formation history (SFH) of the stellar components, obtained from a fit of BGM-Fast on Gaia DR2 data (Mor et al. 2018).
Age ecc IMF SFH
Disc 0-0.15 0.0140
0.15-1 0.0268 dn/dm m-α
1-2 0.0375     α = 0.5 for m <0.5
2-3 0.0555     α = 2.1 for 0.5 <m<1.53 exp.
3-5 0.0703     α = 2.9 for m >1.53 decreasing
5-7 0.0794 γ = 0.14
7-10 0.0799
Stellar halo 14 0.774 dn/dm m-0.5 one burst
Age hz IMF SFH
Thick disc young 10 400. dn/dm m-1 one burst
Thick disc old 12 795. dn/dm m-1 one burst
Bar 10 300. dn/dm 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 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.

Table 2.2: Density laws and associated parameter of the stellar components. a2=R2+z2ϵ2 where R is the galactocentric distance, z is the height above the Galactic plane, and ϵ is the axis ratio. Values of ϵ are given in Table 2.1. For simplicity the disc density law is given here without the warp and flare (see text for their characteristics). d0 are normalisation constants. For the bar, x,y and z are in the bar reference frame and parameter values are given in Table 2.4.
Thin Disc
age0.15 ρ0/d0/kflare×{exp(-(a/hR+)2)-exp(-(a/hR-)2)}
   with hR+ = 5000 pc, hR- = 3000 pc
age >0.15 ρ0/d0×{exp(0.5-(0.25+a2/hR+2)1/2)
   with hR+ = 2011 pc, hR- = 1310 pc
Thick disc ρ0/d0/kflare×exp(-R-RhR)×cosh(z/(2.×hz))
   with hR = 2040 pc, hz = 400 pc young thick disc
   with hR = 2919 pc, hz = 795 pc old thick disc
Spheroid ρ0/d0×a-αh×(rc+a)-δh
   with αh=1.,δh=2.777231 and rc=2175.827881
Bar N×(cosh(rs))-2
   with rPerp=[(x/x0)cPerp+(y/y0)cPerp](1./cPerp)

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 zwarp of the mid-plane is computed following the equations:


with θw=atan2(x,-y) and γ being γ+ or γ- depending on which side of the Galaxy.

It applies when the Galactocentric radii is larger than Rw.

Table 2.3: Warp parameters based on Amôres et al. (2017) study of 2MASS data. Distances are in pc, angles in radian. Rw is the minimum Galactocentric radius at which the warp occurs. θw is the node angle. γ+ and γ- refer to the gradient of the warp on positive and negative sides of the Galactic plane (resp.)
Age Rw θw γ+ γ-
0-0.15 10169.875 1.80685 1.00652 0.06815
0.15-1 10042.375 1.77785 0.84525 0.10915
1-2 9806.5 1.72420 0.61675 0.18500
2-3 9551.5 1.66620 0.47175 0.26700
3-5 9169. 1.57920 0.45300 0.39000
5-7 8659. 1.46320 0.45300 0.39000
7-10 8021.5 1.31820 0.45300 0.39000

The flare is simulated by multiplying the scale height of a population by a factor kflare for Galactocentric radii larger than Rflare. The amplitude of the flare γflare is defined as:


with gflare=0.0001 and Rflare= 10 kpc.

Table 2.4: Bar parameters from Awiphan (priv. comm) fitted on VVV photometric near-infrared survey and HST data, based on sech2 shape from Freudenreich (1998) generalised ellipses.
Parameter Value Unit Comment
x0 1076. pc
y0 490. pc
z0 300. pc
α 0.509 rad bar angle with regards to Sun-galactic centre axis
γ 0 rad bar roll angle
β 0 rad bar tilt angle
D0 0.4508 ×1011 normalisation
RMax 4195.55 pc cut-off radius
cPara 3.007 boxyness in (x,y) plane
cPerp 3.329 boxyness in (x,z) plane
invhEnd 4.e-6 pc-1 sigma of the cut-off
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.

Table 2.5: Provisional parameters of the 2 arms of the spiral structure.
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 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, Teff and age described above, as follows :


where ρ is the density law of the population, r the distance to the sun, (R,θ,z) the Galactocentric cylindrical coordinates, ω the solid angle, N is the theoretical number of stars in a volume element with the intrinsic parameters MV, Teff, 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).

Table 2.6: Abundances in iron and α elements used in the simulations for each stellar population. σ indicated in column 3 and 6 are the cosmic variances used to generate the sample with random Gaussians. The drawn values are however truncated to not exceed 5 sigmas.
Population [Fe/H] σ d[Fe/H]dR [α/Fe] σ
Disc +0.01 to -0.12 0.10 to 0.18 -0.07 -0.121×[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-exp(1.19375×[Fe/H]-1.6) 0.05
Thick disc old -0.8 0.30 0.00 0.32-exp(1.19375×[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-exp(1.19375×[Fe/H]-1.6) 0.05


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 kms-1/pc on σU and σV, and no gradient on σW are applied.

The rotation curve is taken from Sofue (2015), giving a VLSR of 230.6 kms-1 at R=8 kpc. The Sun peculiar velocities are assumed to be (12.75, 0.93, 7.10) in kms-1 following the fit to TGAS+RAVE data (Robin et al. 2017).

Table 2.7: Velocity dispersions (σU,σV,σW) for each sub-component.
Age σU σV σW
(Gyr) (km s)-1 (km s)-1 (km s)-1
Disc 0-0.15 12.691 07.214 05.87
0.15-1 14.528 08.258 06.72
1-2 19.851 11.283 09.182
2-3 24.368 13.851 11.272
3-5 30.361 17.257 14.044
5-7 36.889 20.968 17.063
7-10 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.

For thin disc stars the asymmetric drift is computed from tabulated (R,z,age) distributions. However this table covers the Galactocentric distance range 2R16 kpc. Outside this range the asymmetric drift needs to be extrapolated but the extrapolation is not accurate in the present version. At R<2 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

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

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.

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 O7-B4 stars, 20% of B5-B7 stars and 3% of B8-B9 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 Section 2.2.5). This parameter is the ratio between the envelope radius and the stellar radius (R_env/R_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 vsini< 50 kms-1 and 23% have 50<vsini< 100 kms-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 : vsini< 120 kms-1.

Wolf Rayet stars

Wolf Rayet (WR) densities are computed following the observed statistics from the VIIth catalogue of Galactic Wolf-Rayet stars of van der Hucht (2001). The local column density of WR stars is 2.9×10-6 pc-2 or volume density 2.37×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  in the mean, an absolute V magnitude of -4, a gravity of -0.5, and an effective temperature of 50,000K.