Galaxy modelling through stellar population synthesis

In recent years, models of the Galaxy based on the spirit of population synthesis have considerably helped to improve our knowledge on the Galaxy and its formation scenario. I describe the different approaches used by several teams, how they differ, and present the significant advances which have been made through this approach using two examples: the thick disc and the bulge/bar populations. 1. GALAXY MODELLING, A HISTORICAL VIEW


GALAXY MODELLING, A HISTORICAL VIEW 1.1 Basic approach
First generations of models dedicated to understanding galactic structure and the stellar population distribution were built since the late 70s'.They were dedicated to predict star count distributions in various direction of observation.This approach was aimed at testing the consistency of our overall view of the Galaxy which was obtained from different aspects.At this time our galactic structure knowledge was believed to be sufficient to make this approach possible.
Basic models start with simple assumptions about the luminosity functions and density laws of 2 to 3 populations (a disc, a halo and sometimes a thick disc).They deduce star count predictions using the fundamental equation of stellar statistics.Among them one can cite the significant models of Bahcall & Soneira [1,2], the approach of Gilmore and Reid introducing the thick disc population [3], the first version of the Besancon Galaxy model [4,19], and a few others [6][7][8] etc...
In the near infrared, galaxy models are mostly interested in modelling the galactic plane, the center and the bulge, parts of the Galaxy too difficult to see at optical wavelengths.Near-infrared models [9][10][11][12], described the HR diagram with 20 to 30 components of different spectral types and luminosity classes in order to reproduce star counts in the J and K bands.Such models allow to determine structural parameters of the Galaxy, even if detailed features are not always reproduced.
In order to study more in detail the stellar population of the Galaxy through star counts, more refined models have been built up.Two kinds of such models may be distinguished: Models dedicated to reproduce in detail the stellar distribution, including inhomogeneities, and models to understand the evolutionary history of the Galaxy ("evolutionary models").They are also called "population synthesis models".
The SKY model [13,14] is the most famous model of the first class.It includes a description of the HR diagram in 87 components of different spectral types and luminosity classes and great efforts to include many classes of late type evolved objects (like AGB stars, planetary nebulae, etc...).This make the model much suitable for predicting star counts in the near infrared from 2 to 35 microns using a grid of low resolution spectra.Model parameters have been fitted to star count distributions in many EPJ Web of Conferences directions in this wavelength range.More recent star counts models have been done based on the new large scale surveys like SDSS and 2MASS [15].

Population synthesis models
Instead, evolutionary models use a physical description of galaxy evolution: From a given amount of gas the number and mass distribution of the stars formed are computed assuming an initial mass function and a star formation rate varying with time.The stars evolve on a grid of evolutionary tracks.The model then gives the distribution in the HR diagram as a function of time and density laws are assumed to describe the star distribution all over the Galaxy.Star count and colour distributions are produced by Monte Carlo simulations, including observational errors (Poisson noise and photometric errors) and an extinction law.
Several groups have been developing such approach in the 90's: a group from Padova [16][17][18], a group from Besançon [19][20][21][22].They have in common a halo, a bulge and a thick disc but with slightly different values for the structural parameters.The main difference results from the dynamical point of view.While in Padova model the disc scale heights are assumed and then adjusted to star counts, in Besançon model the scale heights have been computed self consistently with the galactic potential: assuming an age-velocity dispersion relation, the Boltzmann and Poisson equations are used to compute a potential and resulting scale heights at a given age.It results less free parameters to fit to star counts, then less flexibility, but stronger constraints on the other parameters.
Since then, more groups have developed the population synthesis approach, including the dynamical constraints, such as [23] from Padova group, and a revised version of the Besancon model [24].Contrarily to these models, the TRILEGAL model [25] does not include the dynamical consistency, but it is able to deal with more photometric systems, allowing for more flexibility to simulate a diversity of surveys.[26] added the bulge population to the TRILEGAL model.Other models have included local dynamical constraints on the thin disc, specially for the solar neighbourhood [27,28] and for the latter to explain the thick disc by radial migration.
In the following we take the Besancon Galaxy Model (BGM) as a tool to examplify the use of such model for understanding Galactic evolution history of different stellar populations.

THE BESANCON GALAXY MODEL
In the standard BGM model, described in [24], four stellar populations are taken into account: thin disc, thick disc, stellar halo, and bar/bulge (this last population being split into two in [29]).
A standard evolution model is used to produce the disc population, based on a set of evolutionary tracks, a star formation rate constant over 10 Gyr and a two-slope Initial Mass Function.Using this evolution scheme we populate the thin disc dividing it into 7 age components.For each subcomponent the distribution in absolute magnitude and effective temperature is obtained.The preliminary tuning of the disc evolution parameters against relevant observational data was described in [21,22] and further changes are explained in [24].
The thick disc is simulated as a single burst of age 11 Gyr, mean metallicity of −0.5 dex.The spheroid and the bulge are also single bursts, but of age 13 and 10 Gyr resp.Their metallicities are assumed −1.5 for the spheroid and solar for the bulge.The latter value is still under study and might be subject to change using the forthcoming APOGEE data [30].
BGM also computes the kinematics of the stars, using simple assumptions: the kinematics of the stars is computed from velocity ellipsoids (which follow a age-velocity relation), adding up to the proper rotation of the population around the Galaxy and asymmetric drifts.Parameters of the rotation and velocity ellipsoids can be found in [24].
The star counts are computed using the standard equation of stellar statistics.A 3D extinction model is used to correct the photometry and observational errors are finaly added.
One of the great improvements of the BGM in recent years have been the introduction of a 3D map of extinction [31], allowing to simulate with good accuracy the stellar content and observables in fields in the Galactic plane.It allowed to trace the warp and flare in the outer Milky Way [32], and to detect the dust lanes associated with the bar in the Galactic center region [33].

APPLICATION TO THICK DISC ANALYSIS
Since the discovery of the thick disc [3], many studies have tried to characterize this population, and investigated its origin.Its characteristics of age, metallicity, kinematics, and shape are found intermediate between the thin disc and the halo.Hence it has been difficult to determine them, to extract a pure thick disc sample and to estimate its density with regards to the thin disc, as well as its metallicity distribution.[34] find out that thick disc stars can have metallicity up to solar, while [35] detected a metal weak tail.Both tails -high or low metallicity -can be easily polluted by their neighbour populationthin disc and halo populations resp.This population dominates the counts only at medium latitudes and at distance between 1 to 4 kpc from the Galactic plane.Using only high latitudes star counts, it has been shown [15,36] that the scale height and local density of the thick disc are anticorrelated parameters.In order to break the degeneracy, it is necessary to go at low latitudes.But in this case the star counts are highly dominated by the thin disc.The separation between thin and thick disc might be done using abundance data [37], but in this sample the selection biases necessary to get a pure sample impeaches to estimate the relative proportion of the two populations.
Our approach is to use the well constrained model of the thin disc to obtain the thick disc parameters as residuals from the comparison to data.Using SEGUE data [38] we perform an analysis of the parameters of the thick disc population towards 6 directions at latitudes 12< |b| < 25 deg and 50< l < 190 deg, and use a Monte Carlo Markov Chain method (Metropolis Hasting) to investigate the best parameters with a likehood estimator.We conclude that the local density of the thick disc cannot be larger than 2% of the thin disc in the galactic plane.Combining this information with the high latitude constraints, we constrain the thick disc to have a scale height of about 1200 pc, which is a higher value than [15] but still compatible knowing the degeneracy they found between the thick disc and local density.

BULGE AND BAR IN THE CENTRAL AREA OF THE GALAXY
The triaxiality of the main structure in the Galactic central regions was established from 24 microns observations [39], COBE NIR data [40,41] or visible photometry [42].Numerous attempts have been made in order to understand its formation and evolution.Still, it is unclear whether this structure, often called the outer bulge but sometimes named the bar, had its origin from the early formation of the spheroid, via early gravitational collapse, via early dry or wet mergers, or was formed by a bar instability later in the disc.The question of formation history is crucial and necessary to investigate, as a benchmark for understanding the formation of disc galaxies.Kinematical studies are very important to trace the stellar motions and deduce a dynamical history of the populations in the inner regions.New surveys in the near future will help to solve this question, particularly the Gaia mission of the European Space Agency.In the mean time photometric surveys, especially in the infrared, can help to determine the shape and density laws of the stellar populations in place.The question of the shape, orientation and extent is still open, as one can find in the literature values for the angle of the major axis direction with regard to the Sun-centre direction to be between 10 deg and 45 deg, depending on the fields considered, the wavelength and the method used.Among many attempts, significant ones include [41] and [45] analysing COBE-DIRBE data, [42] using OGLE microlensing survey, [43] for its dynamical approach of the bar, [44] identifying an overdensity at longitude 27 deg and proposing the existence of the long thin bar, [46] studying the shape from DENIS near-infrared photometry, [47] investigating the position of the red clump in 5 directions, [48] analyzing Glimpse survey data in the mid-infrared and [49] proposing a new analyse of the OGLE survey.The difficulty in accurately measuring the orientation of the bulge might be aggravated if, actually, there are more than one population present in the region and if these populations have different origins.[50] presented a smart explanation for the disagreement, which they attribute partly to the cone effect and partly to the putative existence of leading or trailing overdensities at the end of the bar.In this context, metallicity distribution can be useful constraints for distinguishing populations.[51], [52] and [53] highlight a bimodal metallicity distribution in the Baade's window and several fields on the bulge minor axis, identifying a mixture of populations.Concerning the kinematics, proper motions are difficult to obtain due to the large distance of the bulge, while radial velocities are expensive to process for large number of stars.[54] studied proper motions from OGLEIII survey and compared with a dynamical model.However the OGLE fields avoid the Galactic plane and the northern part of the bulge, are limited by extinction being at visible wavelength, and the mean proper motion error is large for giving accurate velocities at the bulge distance.[55] have undertaken a wide survey (BRAVA, for Bulge Radial Velocity Assay) of radial velocities at longitudes between −10 deg and 10 deg at latitudes −4 deg and −8 deg.Their first results from the latitude −4 deg data set showed that the slope of the rotation curve flattens considerably at longitude |l| > 5 deg indicating a probable change of the dominant population at this point.Since 2010, there have been several studies showing the existence of double clumps in a number of fields in the bulge [56][57][58].This is also a feature that should be taken into account in an overall view of the bulge region and that the models should explain.

09003-p.3 EPJ Web of Conferences
The complete understanding of the bulge structure and formation would require the study of the abundances and kinematics in various fields at different longitudes and latitudes, as for example the APOGEE project [30].In the mean time we explore the bulge structure based on a re-analysis of existing photometric data.We intend to check whether a detailed modelling of the stellar populations and interstellar extinction present in the central region might enlighten the problem.We use the Two Micron All Sky Survey [59] photometry and a model fitting procedure based on the Besancon Galaxy Model [24] and explore a wide region −20 deg < l < 20 deg and −10 deg < b < 10 deg to fit a multidimensionnal space of parameters defining the density laws of the bulge and disc populations.The analysis proceeds in several steps.First we attempt to fit the whole area with a Monte Carlo scheme to explore the parameter space using the differential counts in bins of K and J-K in a set of selected windows.We further check the 3D extinction distribution and compare the resulting model with integrated counts for having a complete map of residuals and to visually check the realism of the resulting model.In a second step we compare simulated CMD in a sample of directions to check the detailed structure.In a final step we compare model simulations of metallicity distribution function with

09003-p.4
Assembling the Puzzle of the Milky Way data over the minor axis of the bar, allowing to put constraints on the mean metallicity of the composite populations in the bulge region.
From this method we find that the data are not well fitted by a unique bulge population simulated by a Ferrer's ellipsoid.But allowing to fit the sum of two ellipsoids gives a very good fit to the overall star counts and to the CMDs. Figure 1 shows the comparison of 2MASS data with star counts simulated by several models: first a model assuming a single bulge/bar population, second a model with two distinct populations.We show [29] that the model with two populations (a bulge and a bar) reproduces best 2MASS data.With these 2 populations simulated CMDs are well in agreement with the data, as shown for example in figure 2 in Baade's window.Moreover, this model reproduces also the bimodality which is seen in metallicity distribution in fields along the minor axis from [51].
The populations which give the best fit are: a main bar population modeled by a triaxial sech 2 Ferrer ellipsoid with a major axis pointing at an angle of 13 deg with the Sun-Center direction, and scale lengths of 1.46/0.49/0.39kpc, and a total mass of 6.1 × 10 9 M .The age is old, 8 Gyr or more, and the mean metallicity is about solar.The second component is modeled by a triaxial exponential Ferrer ellipsoid less massive (2.6 × 10 8 ) M , with scale lengths of 4.44/1.31/0.80kpc, also old and with a poorer metallicity of −0.35 dex.We emphasize that this second population can be either a classical bulge or a inner counterpart of the thick disc.This model with two components explains well the observed metallicity distribution function and its gradient along the minor axis by a transition between the two components dominating the counts [29].

CONCLUSION AND PERSPECTIVES
The synthetic approach is a very useful tool for testing scenarios and hypothesis about the Galaxy formation and evolution.It allows to better use the informations from kinematics, abundances and density distributions of the stellar populations in order to interpret the data in terms of evolution.It is also useful to simulate the data and for this purpose is used in the preparation of the Gaia mission [60].It is generally a challenge to analyse multivariate data.New efficient methods will need to be developed in order to get profit from such a large and complex survey.To this purpose we also prepare ourself to use the BGM model for bayesian classification of the Gaia data, and for overall comparisons and data assimilation of the whole survey for the benefit of understanding in much more detail the overall scheme of Galaxy formation and evolution.

Figure 1 .
Figure 1.Star counts up to limiting magnitude K S from 2MASS data (top left) compared with 1-ellipsoid model (top middle) and its residuals (N mod -N obs )/ N obs (bottom middle), 2-ellipsoid model (top right) and its residuals (bottom right).In the residual map, contours are drawn at intervals of 20% model overestimate (black) and 20% model underestimate (white).The 2-ellipsoid model gives a much better fit to the dat counts over a wide area.Near the Galactic centre the nuclear bar population is missing in the model and thus appears in the residuals.

Figure 2 .
Figure 2. Colour-magnitude diagrams for the 2MASS field at l = 0, b = −4.(a) data; (b) best fit model with 1 ellipsoid; (c) best fit model with 2 populations; (d) histograms of data (solid red) and models (1 ellipsoid: long dashed green; 2 populations: dotted blue).The CMD and differential star counts are well reproduced by a 2-population model with a bar and a bulge.2MASS data in this field are incomplete at K > 12.5.