© ESO, 2017

1. Introduction

Stellar ages are difficult to measure. They can be studied only in rare cases when accurate, high-resolution spectroscopy, or asteroseismology is available and for limited evolution stages during which the stars quickly change either their luminosity or colours. Alternatively, ages are often indirectly deduced using chemical or kinematical criteria, in which case their determinations are more model dependent. On the other hand, scenarios of Galactic formation and evolution are best investigated when ages are available.

For this reason, Galactic structure parameters, such as scale length, scale heights, warp and flare have mainly been determined for the Milky Way without accounting for time dependence. In a few exceptions, tracers of different ages have been used to measure these structures. This could partly explain why the thin disc scale length in the literature has been claimed to have values in between 1 to 5 kpc. However, it is noticeable that since the year 2000, most results tend towards short values, smaller than 2.6 kpc: Chen et al. (2001): 2.25 kpc, Siegel et al. (2002): 2–2.5 kpc, López-Corredoira et al. (2002): 1.97 $\begin{array}{c}\mathrm{+}\mathrm{0.15}\\ -0.12\end{array}$ kpc, Cabrera-Lavers et al. (2005): 2.1 kpc, Bilir et al. (2006): 1.9 kpc, Karaali et al. (2007): 1.65–2.52 kpc, Jurić et al. (2008): 2.6 ± 0.5 kpc, Yaz & Karaali (2010): 1.–1.68 kpc, Robin et al. (2012): 2.5 kpc. A few studies have given large values, but mostly with large error bars and correlated parameters: Larsen & Humphreys (2003): 3.5 kpc, Chang et al. (2011): 3.7 ± 1.0 kpc, McMillan (2011): 2.90 ± 0.22 kpc, Cheng et al. (2012): 3.4$\begin{array}{c}\mathrm{+}\mathrm{2.8}\\ -0.9\end{array}$ kpc, Bovy et al. (2012): 3.5 ± 0.2 kpc but changing from 2.4 to 4.4 kpc with metallicity, and Cignoni et al. (2008): 2.24–3.00 kpc from open clusters.

To solve this open debate, it is worth investigating further whether these different studies are considering the same populations, with the same age, and what is the accuracy of their distance estimates and the possible biases. It can be more efficient to consider studies, whenever possible, in which the tracers are better identified, and their ages are more or less known or estimated.

The HI disc is known to have a long scale length (Kalberla & Kerp 2009). Using the Leiden/Argentine/Bonn (LAB) survey (Kalberla et al. 2005), Kalberla & Dedes (2008) have analysed the global properties of the HI distribution in the outer Galaxy, determining its mean surface densities, rotation curve, and mass distribution. They obtained a radial exponential scale length of 3.75 kpc in the mid-plane in the Galactocentric distance range 7–35 kpc.

It has been shown that the young object density laws (OB, A stars, Cepheids, open clusters, etc.) follow longer scale lengths than the mean thin disc. Sale et al. (2010) used A type stars to determine the scale length in the outer Galaxy and showed that those stars of mean age 100 Myr have a typical scale length of 3 kpc.

The density profiles and their dependence with age are crucial to understand the scenario of Galactic evolution and formation. If the Galaxy was formed by a process of inside-out formation as proposed by many authors (Larson 1976; Sommer-Larsen et al. 2003; Rahimi et al. 2011; Brook et al. 2012; Haywood et al. 2013), among others, one can expect the scale length to be time dependent. More precisely, the young disc scale length should be larger than the old disc scale length. However, radial migration in the disc can also perturb this simple idea.

Bovy et al. (2012) noticed that the mean scale length of high metallicity thin disc stars (probably younger in the mean) is significantly shorter than the one of lower metallicity stars (probably older). At the first sight, this can be contradictory to the scenario of inside-out formation. However, thin disc metal-poor stars are also more typical of the outer disc and could reach the solar neighbourhood by migration and high metallicity stars are the Sun position can be older if they come from the inner Galaxy. Haywood et al. (2013) suggested that the evolution of the thin outer disc is disconnected from the thin inner disc and the thin disc scale length would vary with time. It is still a question whether the structure of the outer disc is in a steady state or perturbed by active merging that could manifest by the recently discovered sub-structures like the Monoceros ring (Newberg et al. 2002; Rocha-Pinto et al. 2003) or the Canis Major overdensity (Martin et al. 2004).

It is well known that like many large spirals, the Milky Way is warped and flared in its outskirts. The evidence comes from gas tracers such as HI (Henderson et al. 1982; Burton & te Lintel Hekkert 1986; Burton 1988; Diplas & Savage 1991; Nakanishi & Sofue 2003; Levine et al. 2006) or molecular clouds (Grabelsky et al. 1987; Wouterloot et al. 1990; May et al. 1997). Recent analysis of the outer Galaxy either from 2MASS (López-Corredoira et al. 2002; Momany et al. 2006; Reylé et al. 2009) or from SDSS (Hammersley & López-Corredoira 2011) has led to the conclusion that even the stars follow a warped structure. However, it has not been clearly established that the shape of the warp is similar or deviates from the gas structure. It is now established that the Galaxy flare is evenly populated by young stars, as recently discovered by Feast et al. (2014) from the Cepheids (less than 130 million years old) in the outer disc flare at 1–2 kpc from the plane. Kalberla et al. (2014) have compared the HI distribution with stellar distribution (2MASS, SDSS, SDSS-SEGUE, pulsars, Cepheids) from several authors. They argue for a typical flaring of gas and stars in the Milky Way. Abedi et al. (2014) explored the possibility of determining the warp shape from kinematics in the future Gaia catalogue.

The existence of the warp can originate from perturbations of the Galaxy by the Magellanic Clouds (Weinberg & Blitz 2006). Perryman et al. (2014) argue that the tilt of the disc may vary with time, invoking four reasons for this: i) the combination of the infall of misaligned gas (Shen & Sellwood 2006); ii) the interaction of the infalling gas with the halo (Roškar et al. 2010); iii) the effect of the Large Magellanic Cloud (Weinberg & Blitz 2006); iv) the misalignment of the disc with the halo (Perryman et al. 2014).

In this work, we have investigated the shape of the outer Galactic disc, considering its structural parameters, such as scale length, but also the non-axisymmetric part, as warp, flare, and disc truncation, as well their dependence with age. The inner Galaxy and spiral structure will be presented in a forthcoming paper. To determine accurate estimates of Galactic disc structure parameters, we compare colour–magnitude diagrams observed with 2MASS with the predicted ones by a population synthesis model for the Galactic plane towards the second and third quadrants at | b | ≤ 5.5°. In the Besançon Galaxy Model the stellar ages serve as the driving parameters for stellar evolution, metallicities, scale heights and kinematics. Hence, it is the most useful model to investigate the time dependence of the Galactic structure parameters. The parameter fitting is done using a powerful method for global optimisation called genetic algorithms (GA). The GAs have been extensively employed in different scientific fields for a variety of purposes. The main strategy consists of adjusting the parameters of the thin disc population, such as scale length, warp, flare, disc edge, in order to reproduce the star counts observed by 2MASS.

The thin disc region can suffer from crowding in certain regions, from interstellar extinction and from clumpiness which could be difficult to model with simple assumptions. However, it is possible to model the extinction in 3D appropriately (Amôres & Lépine 2005; Marshall et al. 2006, among others), such that this effect is taken out in the analysis of the Galactic plane stellar content.

The remainder of this paper is organised as follows. In Sect. 2 we present the sub-sample of 2MASS data used in the present work. In Sect. 3 we describe the properties of the population synthesis model (its parameters adjusted in the current work) and we compare the standard model with 2MASS data. The basic concepts of the GA method and its implementation in our study are described in Sect. 4. Analysis of scale length and its dependence on age, the warp and flare shapes are presented in Sect. 5. The overall discussion of Galactic structure parameters and its parameters are presented in Sect. 6. In Sect. 7 we address the conclusions of this study and give some final remarks.

2. The data

The 2MASS survey (Skrutskie et al. 2006) represents the most complete (regarding spatial coverage) and homogeneous data set in the near infrared (in the J, H and K_s bands) available for the entire Galactic plane. We have performed a selection on the 2MASS data based on the 2MASS Explanatory Supplement using selections proposed and discussed by Cambrésy et al. (2003). The sources to be accepted should satisfy the following criteria: i) the photometry uncertainty σ ≤ 0.25; ii) signal-to-noise ratio is larger than seven in at least one band; iii) contamination of extended sources must be avoided; iv) the field of photometry quality (Qflag) must be different from: X (no possible brightness estimate), U (source not detected in the band or not resolved), F (poor estimate of the photometric error), E (quality of photometric profile is poor); v) the read flag must be different from zero, four, six, or nine. Those values indicate either non-detections or poor quality photometry and positions.

In order to adjust parameters towards the outer Galaxy, we have used fields located at 80° ≤ ℓ ≤ 280° and | b | ≤ 5.5° distributed every 3° and 10° in longitude for | b | ≤ 3.5° and | b | > 3.5°, respectively. In total, there are 2228 fields. We simulated square fields with an area equal to 0.25 × 0.25 square degree in order to account for the changes in extinction at this scale. Then to obtain higher statistics, we grouped four fields in latitude at any given longitude. The final fields have a size of 0.25° in longitude and 1° in latitude, totalizing 557 fields.

In each field, the completeness limit of the observed sample is estimated, and fainter stars are discarded. To obtain the completeness limits, distributions of star counts as a function of magnitude are built for each filter with bin size equal to 0.2 mag. The bin before the peak gives the respective completeness limit. The source also needs to satisfy either the nominal completeness limits of 2MASS in J (15.8 mag) and K_s (14.3 mag) bands or the completeness for a given field and also to be detected in J and K_s bands. We notice that following the criteria above, the large majority of the stars in our sample that have simultaneous detection in the J and K_s bands, are also detected in the H band. Approximately 55% and 20% of the fields have a completeness limit at K_s equal to 14.1 and 14.3 mag, respectively. The total number of 2MASS stars used in the present work is 886 916.

3. The Besançon Galaxy model

To produce star counts and colour–magnitude diagrams, we make use of the population synthesis model called Besançon Galaxy Model, hereafter BGM (Robin & Crézé 1986; Robin et al. 2003, 2012). It provides a realistic description of the stellar content of the Galaxy, assumptions about the star-formation scenario and evolution in different populations, and includes kinematics and dynamics as further constraints on the mass distribution. One of the main differences between this and other Galaxy models resides in the fact that the BGM is dynamically self-consistent locally (Bienaymé et al. 1987). Here, we recall the main parameters relevant for the present study. The BGM is composed of four components: a thin disc, thick disc, halo and bar.

Since its first version, the BGM has been extensively compared with several large surveys in different wavelengths and at different depths. With regards to the on-line version (Robin et al. 2003) here we have used an update which benefits from new results concerning the 3D extinction model (Marshall et al. 2006), the shape of the warp and flare (Reylé et al. 2009), the bar-bulge region (Robin et al. 2012).

Recently, a more flexible version of the model (Czekaj et al. 2014) has been proposed, which allows modification of the initial mass function (IMF) and star-formation history (SFH) of the Galactic thin disc. As it is still undergoing testing, its use will be deferred to future studies. New revisions concerning the thick disc and halo populations have been discussed by Robin et al. (2014). We do not make use here of the new characteristics of the thick disc and halo, but they should affect the present study only very marginally.

The disc population is assumed to have an age, ranging from 0 to 10 Gyr. The initial values for the evolutive parameters of the disc (star-formation rate history, initial mass function) were obtained by Haywood et al. (1997) from the comparisons with observational data. The density laws for each component can be found in Robin et al. (2003) in which the thin disc follows Einasto laws rather than a double exponential. The disc axis ratio of each age population, presented in Table 1, have been computed assuming an age-velocity dispersion relation from Gómez et al. (1997), using the Boltzmann equation, as explained in Bienaymé et al. (1987) and revised in Robin et al. (2003).

In the present work, we have adjusted parameters of the Einasto law (scale length) and the three most important structures towards outer Galaxy, for example, warp, flare and disc truncation. We postpone the analysis of the inner Galaxy and spiral arms to a further study. The scale length, warp and flare are adjusted considering their dependence on age. While in the standard BGM the disc is truncated at R_gal = 14 kpc, here we use simulations without any truncation to be able to determine it during the fit. The standard values for those parameters are presented in Table 2.

3.1. The warp and flare

We adopt the same representation used by Reylé et al. (2009) in which the height z_warp of the warp (for R > R_warp) is defined as the distance between the mid-plane of the disc and the plane defined by b = 0°. It varies as a function of Galactocentric radius (R) as follows. $\begin{array}{ccc}& & {\mathit{z}}_{\mathrm{warp}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}{\mathit{\gamma }}_{\mathrm{warp}}\mathrm{\times }\mathrm{(}\mathit{R}\mathrm{-}{\mathit{R}}_{\mathrm{warp}}\mathrm{)}\mathrm{\times }\sin \mathrm{(}{\mathit{\theta }}_{\mathrm{u}}\mathrm{-}{\mathit{\theta }}_{\mathrm{warp}}\mathrm{)}\mathit{,}\\ & & \mathrm{where}{\mathit{\theta }}_{\mathrm{u}}\mathrm{=}\mathrm{atan}\mathrm{2}\mathrm{\left(}\mathit{y,x}\mathrm{\right)}\mathit{,}\end{array}$(1)and θ_warp is the node angle; γ_warp and R_warp are the slopes of the amplitude and the Galactocentric radius at which the Galactic warp starts, respectively; x and y are Galactocentric coordinates, x positive in the Sun-Galactic centre direction, and y positive towards rotation (ℓ = 90^o).

As shown from the 2MASS star counts in Reylé et al. (2009), the Galactic warp acts differently on different sides (second and third quadrants) of the Galactic disc. We have considered two values for the warp slope for the second (γ_warp+) quadrant, for example, points towards longitudes less than θ_warp and third (γ_warp−) quadrant, for example, points towards longitudes greater than θ_warp.

Concerning the flare, we use the same representation provided by Gyuk et al. (1999) who modelled the flare by increasing the scale heights by a factor k_flare with an amplitude γ_flare, beyond a Galactocentric radius of R_flare: ${\mathit{k}}_{\mathrm{flare}}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ \mathrm{1}\mathrm{+}{\mathit{\gamma }}_{\mathrm{flare}}\mathrm{(}\mathit{R}\mathrm{-}{\mathit{R}}_{\mathrm{flare}}\mathrm{)}\mathit{,}& \mathrm{(}\mathrm{if}\mathit{R}\mathit{>}{\mathit{R}}_{\mathrm{flare}}\mathrm{)}\\ \mathrm{1}& \mathrm{(}\mathrm{if}\mathit{R}\mathit{<}\mathrm{=}{\mathit{R}}_{\mathrm{flare}}\mathrm{)}\mathit{.}\end{array}$(2)

3.2. Disc truncation

Disc truncations were first discovered in external galaxies (van der Kruit 1979). The first studies in our Galaxy were done by Habing (1988) from OH/IR stars and by Robin et al. (1992a) from UBV photometry. Whether or not the disc is truncated, and how the truncation scale occurs can play a role in the determination of the thin disc scale length. The truncation can be related to the star-formation threshold in outer galaxies (Kennicutt 1989). Due to the presence of the warp and flare, it might happen not to be circular. Seiden et al. (1984) found evidence for disc truncation from the point of view of a process of stochastic star formation.

For the gaseous component, Wouterloot et al. (1990) found a decrease in the CO density distribution between 18 and 20 kpc (kinematic distance method).

To model the disc truncation, we have considered a radial cut as proposed by Robin et al. (1992b). The thin disc truncation is computed by multiplying the density law by a factor f which is defined by a Gaussian truncation with a scale h_cut: $\mathit{f}\mathrm{=}\{\begin{array}{c}\\ \exp \left(\mathrm{-}{\mathrm{(}{\frac{\mathit{R}\mathrm{-}{\mathit{R}}_{\mathrm{dis}}}{{\mathit{h}}_{\mathrm{cut}}}}^{\mathrm{)}}}^{\mathrm{2}}\right)\mathit{,}& \mathrm{if}\mathit{R}\mathrm{\ge }{\mathit{R}}_{\mathrm{dis}}\\ \mathrm{1.}& \mathrm{if}\mathit{R}\mathit{<}{\mathit{R}}_{\mathrm{dis}}\end{array}$(3)where h_cut is the truncation scale and R_dis is the Galactocentric radius of disc truncation.

3.3. Other simulation parameters

The simulations made by the BGM have to assume realistic photometric errors, and the comparison should be made in a magnitude range where the data are complete. The completeness limit for each field was applied for two filters, J and K_s, using the completeness limits obtained from 2MASS data. The observational errors are assumed to follow an exponential law as a function of magnitude (Bertin 1996) as in Eq. (4). ${\mathit{\sigma }}_{\mathit{M}}\mathrm{=}\mathit{A}\mathrm{+}\exp \mathrm{(}\mathit{C}\mathrm{\times }\mathit{m}\mathrm{-}\mathit{B}\mathrm{)}\mathit{,}$(4)where m is the magnitude in the observed band, and A, B, C are parameters obtained by fitting on 2MASS observational errors. These parameters are computed for each 2MASS field, because of varying observational conditions. The fit was performed on 2MASS data before cutting them at the completeness limit.

Since the fields used in the present work are in the Galactic plane, the interstellar extinction is crucial in the analysis of colour–magnitude diagrams (Amôres & Lépine 2005; Marshall et al. 2006). Hence, it is mandatory to use a model with good spatial resolution and to estimate the reddening accurately. However the parameters which are being studied (warp, flare, scale length, disc truncation) affect only the density, meaning that the amplitudes and not the shape of the histograms (J−K_s), and cannot mimic a reddening. To have the best estimate of the distribution of extinction along every line of sight, we adopt the 3D extinction model proposed by Marshall et al. (2006) with a posterior revision Marshall (2009, priv. comm.) covering the entire Galactic plane including Galactic anti-centre regions.

3.4. Catalogue of simulated stars

We start the process by creating a catalogue of simulated stars from the standard axisymmetric model. A typical simulation using the model provides a catalogue of stars with their properties, such as distance, colours, magnitude, age, luminosity class, effective temperature, gravity and metallicities. We apply the same selection function on the model simulation as on the data. In total, there are 912 417 pseudo-stars produced by the model and 886 916 observed by 2MASS with the same selection.

Figure 1 shows the distribution of simulated stars according to their age bin, see Table 1, where the age class (AC) equal to eight corresponds to thick disc stars.

Next, we distributed the stars in colour bins of size 0.5 in J−K_s to obtain star counts N(J−K_s, ℓ, b). Finally, we merged the bins in colour (J−K_s) in which the number of stars is less than 100 (either for 2MASS or BGM) with the right neighbour bin in J−K_s.

We discarded a few bins (43) in J−K_s which show large discrepancies with their neighbours in the (ℓ, b) space, where we suppose that either the data suffer from errors or the extinction is not well modelled. The final total number (N_bin in Eq. (5)) of bins is 1615. In the catalogue of simulated stars, there are 23 193 and 7292 stars for the range of 15 kpc ≤R< 18 kpc and R ≥ 18 kpc, respectively, which allows us to constrain the distribution of stars at large distances in the Galactic plane.

3.5. Comparison between the standard model and 2MASS

Next, to analyse the differences between the standard model with 2MASS data we computed the relative residuals in each colour bin, as defined below. ${\mathit{ϵ}}_{\mathrm{rel}}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{bin}}}\mathrm{(}{\mathit{N}}_{\mathit{i,}\mathrm{obs}}\mathrm{-}{\mathit{N}}_{\mathit{i,}\mathrm{model}}\mathrm{)}\mathit{/}{\mathit{N}}_{\mathit{i,}\mathrm{obs}}\mathit{,}$(5)in which N_i,model, and N_i,obs are the model (BGM standard) and observed counts in the space (J−K_s, ℓ, b), and N_bin is the number of bins.

The differences in star counts as a function of longitude between the standard model and 2MASS data are presented in Fig. 2 and as a map of relative residuals in Fig. 3.

4. The optimisation procedure

In order to fit model parameters, we have used a genetic algorithm, a method based on evolutionary mechanisms and theories, using selections of models similar to natural selection and genetics. This method presents some characteristics which make this technique more efficient than the usual heuristic methods based on calculus, either random or enumerative processes as pointed out by Holand (1975) and Mitchell (1996), and references therein. The GAs are also related to artificial intelligence, having the ability to learn by experimenting, as used in many computational domains.

In astronomy, different versions of the GA have been employed in several different applications, involving Galaxy modelling, such as Sevenster et al. (1999) and stellar population diagnostics by colour–magnitude diagrams (Ng 1998; Ng et al. 2002). In another instance, Larsen & Humphreys (2003) applied the GA to retrieve eight parameters of the Galactic structure. They performed comparisons between their model counts and the data from the Automated Plate Scanner Catalog for 88 fields.

In the present work, we have used the version of GA called PIKAIA¹. This optimisation subroutine was first presented by Charbonneau (1995) and Charbonneau & Knapp (1996). PIKAIA works with 12 parameters (Charbonneau 1995). The choice for the values of those parameters depends on the application. After some tests, we chose the values of these parameters described in Table 3. Some of them are slightly different from the standard ones; they are appropriate when a large number of parameters are adjusted. They force a higher search in the space of parameters, as for instance, the crossover probability and the maximum mutation rate.

The ngen parameter defines the number of the generations used. In our case, a too-small value would cause a premature solution. On the other hand, we identified that for our problem, there is no significant improvement in the χ² with ngen values of greater than 300.

The nd parameter defines the number of significant digits retained in chromosomal encoding (Charbonneau 1995). For the parameters used in the present work, a value of 5 for nd was found to be the most reasonable value as a compromise between accuracy and performance. Indeed a run with nd = 6 gave similar parameter values and χ². Regarding CPU time, the difference was approximately 10% larger than with nd = 5.

We also have used the elitism technique that consists of storing away the parameters of the best-fit member of the current population and later copying them into the offspring population (Charbonneau 1995). Use of this technique also avoids a lost solution of good gens by either mutation or crossover. As we have used elitism, it is necessary to use the options (irep) 1 or 2 (see Table 3) in the reproduction plan. We have used the option irep = 1 (full generational replacement). We also did a run using irep = 2, but the χ² was similar, as well as the parameters, considering the range of standard deviation.

The main procedure consisted of comparing the counts obtained from drawn parameters with the observed data, and making the parameters evolve to improve the figure of merit. In order to optimise the process and to avoid recomputing all simulations each time the model parameters are changed, we attribute to each star (s) of coordinates (x, y, z) a weight (w_s) which is the ratio between the new density (ρ_new(x,y,z)) obtained by parameter fitting and the standard one (ρ_std(x,y,z)) obtained with BGM standard parameters, as described below: ${\mathit{w}}_{\mathrm{s}}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{new}}\mathrm{\left(}\mathit{x,y,z}\mathrm{\right)}}{{\mathit{\rho }}_{\mathrm{std}}\mathrm{\left(}\mathit{x,y,z}\mathrm{\right)}}\mathrm{·}$(6)This was only applied to the thin disc population. Other populations were unchanged in the process. Positions are given as a function of Galactocentric cartesian coordinates (x, y, z).

The total number of stars (N_i,new) modelled (for the given set of parameters) for each bin (i) is given by Eq. (7): ${\mathit{N}}_{\mathit{i,}\mathrm{new}}\mathrm{=}\sum _{\mathit{s}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathit{i,}\mathrm{std}}}{\mathit{w}}_{\mathrm{s}}\mathit{.}$(7)N_i,std is the number of stars in the bin (i) in the standard model.

The merit function is presented in Eq. (8): ${\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{bin}}}\mathrm{(}{\mathit{N}}_{\mathit{i,}\mathrm{obs}}\mathrm{-}{\mathit{N}}_{\mathit{i,}\mathrm{new}}{\mathrm{)}}^{\mathrm{2}}\mathit{/}\mathrm{(}{\mathit{N}}_{\mathit{i,}\mathrm{obs}}\mathrm{+}{\mathit{N}}_{\mathit{i,}\mathrm{new}}{\mathrm{)}}^{\mathrm{2}}\mathit{,}$(8)N_i,obs is the number of stars in bin (i) observed by 2MASS, N_bin is the number of bins (1615) and N_i,new defined in Eq. (7).

We preferred to use this relation instead of traditional χ² in order to have a reasonable weight for bins with a large number of stars. We note that a relative χ² avoids overweighting the contribution of latitude bins with high star counts but without necessarily high contrast. Larsen & Humphreys (2003) used a similar one. As a reference, the total χ² of the standard model is 33.32 for 1615 bins, distributed in 14.70 towards second quadrant (784 bins) and 18.62 towards third quadrant (831 bins). Table 4 shows the range of parameters involved in those simulations.

5. Results

5.1. Fitting the parameters of the standard model

The standard BGM works with two scale lengths for disc stars, first one, kp₁, for stars with age < 0.15 Gyr and the second one, kp_2...7, for stars ranging from 0.15 to 10.0 Gyr (see also Table 1). We performed a set of fitting procedures (sim-1) with 100 independent runs of 300 generations each, considering these two scale lengths for disc stars, and warp and flare parameters. Table 5 summarises the fitted parameters in sim-1, as well as the median and standard deviation for each parameter.

The scale length for stars with age larger than 0.15 Gyr is best fitted by a scale length of 2.48 kpc. This value is similar to the scale length used in the standard BGM, and also in agreement with many recent studies (see discussion in Sect. 6). As shown in Table 5, the warp is found asymmetrical, with a stronger slope in the second quadrant than in the third quadrant. The flare is also very significant. The starting radius of the flare is close to the starting radius of the warp.

Figure 2 shows longitudinal profiles for sim-1 versus BGM standard and 2MASS data for three ranges of Galactic latitudes. There is a net improvement with the new fits at nearly all longitudes.

Depending on latitudes, the sim-1 model performs better or similarly to the standard BGM. For example, the agreement is better at | b | < 0.5°, apart from the longitude ℓ< 100° where it is similar. Though, at the same longitudes but higher latitudes (b> 3.5°) sim-1 is much better. Overall, sim-1 performs much better than the standard model, although it is not perfect everywhere. We expect that a more complex model might give a better solution (see Sect. 5.2).

Indeed, part of the differences at ℓ ~ 80° could be attributed to the Local Arm (Drimmel 2000) that is not modelled in the present work. He also pointed out that the counterpart of the Local arm is centred at ℓ ~ 260°, where we also see some disagreement between our model and 2MASS data. The fitting of spiral arms using BGM with 2MASS data is outside of the scope of the present paper. It will be analysed in detail in a forthcoming paper where we attend to fit a spiral model including data from the inner Galaxy.

The irregular profile can be due to the interval of bins in longitude and because each field has its extinction, photometric errors and completeness limits. Another source of error resides on the fact that the map of Marshall (2009, priv. comm.) has a lower resolution for outer Galaxy fields than for inner part. For instance in some fields, we got extinction determination for a distance of approximately 0.5 deg of the centre of our field.

Figures 3a and b show the counts observed by 2MASS and modelled by BGM Standard, respectively. Figures 3c and d show the (ℓ, b) map of the relative residuals for the standard model and sim-1, respectively. In Fig. 3c, some regions in orange, red, yellow show where the standard model predicts more counts than observed. An interesting pattern is that the position of the mid-plane slightly differs in the model from the data. This is particularly the case in the second quadrant between ℓ ~ 120° and 180°. In the third quadrant, the mid-plane seems to be almost correctly placed but, the model overestimates the counts at latitudes ~0°. A similar feature, due to dissymmetry of the disc in the warp and probably in disc truncation, has been pointed out by Reylé et al. (2009). They show that the warp slope was different in the second and third quadrants. While the differences towards ℓ ~ 90° and 270° can be attributed mainly to the warp, the overdensities found in Galactic anti-centre are probably due to the disc scale length, disc edge and flare effects, as it is investigated below.

In Fig. 3d, in the comparison of sim-1 with 2MASS data, we can identify that most of the fields for both second and third quadrants have a relative difference in the range of 0− 20%. A minuscule area can be seen at (ℓ, b) ~ (210°, –2°). We note that in analysing the Marshall (2009, priv. comm.) extinction map, this region has a high extinction with A_V of about 10 to 12 mag, while in the neighbouring fields the extinction drops to four magnitudes in A_V. This explains the reduced counts not only in the model counts but also in 2MASS data as can be seen in Fig. 3a. This high extinction area is not well taken into account in the 3D extinction model used.

It can also be seen that both second and third quadrants present a significant improvement with regards to the standard model. We note that it is still better in the third quadrant, apart from the region of the local arm. Some residuals appear larger than 10% at ℓ ~ 140°, which could be related to the outer arm, not included here. After the optimisation procedure, most of the fields with green colour, for example, relative differences around 30–40%, decreases to residuals smaller than 20%. Figure 4 (upper panel) shows the histogram of the differences N_i,model−N_i,obs per bin for the standard model and sim-1. As can be seen in Fig. 4 (lower panel), there is a significant increase of bins with a relative difference of less than 20% in sim-1.

In order to investigate whether the discrepancies between simulations and data can be attributed to the interstellar extinction, we have analysed eight colour–magnitude diagrams (K_s, J−K_s) towards the regions with substantial differences in the modelled counts. Figure 5 shows those diagrams for 2MASS data and simulations with the standard model. It is clear that the extinction cannot be invoked to interpret the difference in star counts seen between 2MASS data and the model, the colours being in good agreement. We note that both standard BGM simulations and the fitted models named sim-1 and sim-2 (see next section) use the same extinction map.

To improve the model, we have explored whether dependencies of parameters with age can better describe the observed data. This is presented in the next section.

5.2. Dependence of the warp, flare and scale length with age

In order to investigate whether the warp and flare parameters depend on age (as the scale length), and to constrain the warp formation scenario, four sets of simulations were performed dividing the ages into two groups, as the number of parameters to fit would be too high to estimate the warp shape individually for each Age Class bin. For instance, the first run of optimisation considered one group for the stars with AC = 1 and another group for stars with AC ≥ 2; the second run, one group for stars with AC ≤ 2, and the other group for stars with AC ≥ 3, and so on. Those tests showed that the dependence on age can be modelled linearly for most of the parameters, except for the warp in the second quadrant, which follows roughly a second order polynomial.

Then, we performed new fits assuming age dependencies of parameters, either linearly or as a second order polynomial, as given in Eqs. (9) to (15) below with the units provided in Table 6. $\begin{array}{ccc}& & {\mathit{\gamma }}_{\mathrm{war}{\mathrm{p}}_{\mathrm{+}}}\mathrm{=}{\mathit{a}}_{{\mathit{\gamma }}_{\mathrm{war}{\mathrm{p}}_{\mathrm{+}}}}\mathrm{-}{\mathit{b}}_{{\mathit{\gamma }}_{\mathrm{war}{\mathrm{p}}_{\mathrm{+}}}}\mathrm{\times }\mathrm{⟨}\mathrm{ag}{\mathrm{e}}_{\mathrm{i}}\mathrm{⟩}\mathrm{+}{\mathrm{c}}_{{\mathit{\gamma }}_{\mathrm{war}{\mathrm{p}}_{\mathrm{+}}}}\mathrm{\times }\mathrm{⟨}\mathrm{ag}{\mathrm{e}}_{\mathrm{i}}{\mathrm{⟩}}^{\mathrm{2}}\mathit{,}\\ & & {\mathit{\gamma }}_{\mathrm{war}{\mathrm{p}}_{\mathrm{-}}}\mathrm{=}{\mathit{a}}_{{\mathit{\gamma }}_{\mathrm{war}{\mathrm{p}}_{\mathrm{-}}}}\mathrm{+}{\mathit{b}}_{{\mathit{\gamma }}_{\mathrm{war}{\mathrm{p}}_{\mathrm{-}}}}\mathrm{\times }\mathrm{⟨}\mathrm{ag}{\mathrm{e}}_{\mathrm{i}}\mathrm{⟩}\mathit{,}\\ & & {\mathit{R}}_{\mathrm{warp}}\mathrm{=}{\mathit{a}}_{{\mathit{R}}_{\mathrm{warp}}}\mathrm{-}{\mathit{b}}_{{\mathit{R}}_{\mathrm{warp}}}\mathrm{\times }\mathrm{⟨}\mathrm{ag}{\mathrm{e}}_{\mathrm{i}}\mathrm{⟩}\mathit{,}\\ & & {\mathit{\theta }}_{\mathrm{warp}}\mathrm{=}{\mathit{a}}_{{\mathit{\theta }}_{\mathrm{warp}}}\mathrm{-}{\mathit{b}}_{{\mathit{\theta }}_{\mathrm{warp}}}\mathrm{\times }\mathrm{⟨}\mathrm{ag}{\mathrm{e}}_{\mathrm{i}}\mathrm{⟩}\mathit{,}\\ & & {\mathit{R}}_{\mathrm{flare}}\mathrm{=}{\mathit{a}}_{{\mathit{R}}_{\mathrm{flare}}}\mathrm{+}{\mathit{b}}_{{\mathit{R}}_{\mathrm{flare}}}\mathrm{\times }\mathrm{⟨}\mathrm{ag}{\mathrm{e}}_{\mathrm{i}}\mathrm{⟩}\mathit{,}\\ & & {\mathit{\gamma }}_{\mathrm{flare}}\mathrm{=}{\mathit{a}}_{{\mathit{\gamma }}_{\mathrm{flare}}}\mathrm{+}{\mathit{b}}_{{\mathit{\gamma }}_{\mathrm{flare}}}\mathrm{\times }\mathrm{⟨}\mathrm{ag}{\mathrm{e}}_{\mathrm{i}}\mathrm{⟩}\mathit{.}\end{array}$The scale length also has a tendency to decrease with age in those simulations. Hence, we adopted a dependence as shown in Eq. (15), which adequately follows the variations seen in the tests with two age groups. ${\mathit{h}}_{\mathrm{R}}\mathrm{=}{\mathit{h}}_{\mathrm{ra}}\mathrm{+}{\mathit{h}}_{\mathrm{rb}}\mathrm{\times }\exp \mathrm{(}\mathrm{-}{\mathit{h}}_{\mathrm{rc}}\mathrm{⟨}\mathrm{ag}{\mathrm{e}}_{\mathrm{i}}\mathrm{⟩}\mathrm{)}\mathit{,}$(15)where ⟨ age_i ⟩ is the mean age of Age Class i, as indicated in Table 1.

Using the equations above, with the range of parameters given in Table 6 (second column), we have performed a new set of optimisations (called sim-2) with 100 independent runs, using the ranges of parameters presented in Table 6, together with the best fit values and estimated uncertainties. Equation (15) is applied on ages from 0.15 to 10 Gyr, while the first age bin is fitted with the parameter kp₁, as in sim-1. The standard errors can be high in some cases, such as c_γwarp+ and b_γflare. Although their impact on our estimates of the warp and flare shapes is limited and concerns mainly oldest stars, as will be shown in Figs. 7 and 8.

As can be seen in the comparison of χ² in Tables 5 and 6, there is an improvement of χ² in sim-2 compared to sim-1. To estimate the number of parameters and the effective improvement in the χ², we compute the Bayesian information criterion (BIC; Schwarz et al. 1978). We follow the recipes of Robin et al. (2014) first computing the reduced likelihood (Lr, Eq. (16)) for a binomial statistics (Bienaymé et al. 1987), as described below: $\begin{array}{ccc}\mathit{Lr}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}{\mathit{q}}_{\mathit{i}}\mathrm{\times }\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{R}}_{\mathit{i}}\mathrm{+}\ln \mathrm{(}{\mathit{R}}_{\mathit{i}}\mathrm{\left)}\mathrm{\right)}\mathit{,}& & \end{array}$(16)where f_i and q_i are the number of stars in bin i in the model and the data, and ${\mathit{R}}_{\mathit{i}}\mathrm{=}\frac{{\mathit{f}}_{\mathit{i}}}{{\mathit{q}}_{\mathit{i}}}$.

The BIC is then computed following the formula from Schwarz et al. (1978), BIC = −2 × Lr + k × ln(n). It penalises models with a larger number of parameters and allows to compare the goodness of fit of different models with a different number of fitted parameters. n is the number of bins on which the data are fitted and k the number of fitted parameters. In our case, n = 1615, and k = 10 and 19 for sim-1 and sim-2, respectively.

Then, Lr (Eq. (16)) and the BIC were computed for the three fitting procedures with their values showed in Table 5 (sim-1) and Table 6 (sim-2 and sim-3). The BIC for sim-2 and sim-3 are quite similar despite the fact that sim-3 incorporates extra parameters to model the ripples, but we did not adjust these parameters. Hence the number of fitted parameters k is the same in sim-2 and sim-3. The comparison between sim-1 and sim-2 highlights the fact that sim-2, despite using more parameters than sim-1, is preferred on the criterion of goodness of fit.

The values obtained for each parameter in sim-2 are given in Table 6 (third column). The map of the relative differences between sim-1 and sim-2 (N_sim−2 – N_sim−1)/N_sim−1 is presented in Fig. 6. The average difference is around 10%. There are two main large area where the two fitted models differ. The first one at ℓ ~ 150° and b< −2° where sim-2 produces more counts than sim-1, and the area around ℓ ~ 200° where it produces less. The shape of the residuals seems to indicate that the change in the mean scale length with age in sim-2 affects the counts in the anti-centre, as expected. But also the shape of the warp varying with time has an impact on the star counts at the level of about 10 to 15% in the specific area at ℓ ~ 150° and b< −2°.

The variation of the warp, flare and scale length with age over the 100 different independent solutions are presented in Figs. 7 to 9, respectively. The box plots were produced by using IDL Coyote’s Library²; the box encloses the InterQuartile Range (IQR), defined as Q₃ – Q₁, in which Q₃ and Q₁ are the upper and lower quartile, respectively. The whiskers extend out to the maximum or minimum value of 100 independent solutions, or to the 1.5 times either the Q₃ or Q₁, if there is data beyond this range. The small circles are outliers, defined as values either greater or lower than 1.5 ×Q₃ or Q₁.

We note that all warp parameters significantly change with age (Fig. 7). The slope of the warp for both second and third quadrants varies as a function of age. The shapes significantly vary from second to third quadrant, even though error bars are larger in the third quadrant and for stars older than 6 Gyr. One common aspect resides on the fact that the amplitude increases for stars with age greater than 4 Gyr for both quadrants.

These figures show that the warp shape changes significantly with time (or with the age of the tracers), and the variations are different in the second and in the third quadrants. The starting radius of the warp changes monotonically with age, as well as the node angle. But for young stars, the slope of the warp is very different between second and third quadrants. This confirms the strong dissymmetry between the region of the plane moved to positive latitudes (first and second quadrants) and the opposite region (third and fourth quadrant), but points out that this dissymmetry is mainly due to young stars, while old stars follow a more symmetrical warp, even though the error bars for old stars are larger.

The asymmetry of the warp was long identified in HI gas (Burton 1988; Levine et al. 2006; Kalberla & Dedes 2008) and even in stars (López-Corredoira et al. 2002; Reylé et al. 2009). But this is the first time that it is shown that the shape varies with the age of the tracers considered.

The flare also exhibits a time dependence (Fig. 8). Despite the fact that error bars of the flare amplitude increase with time, there is a clear trend showing that the amplitude increases with age. The starting radius of the flare also varies with age. We shall consider in Sect. 5.4 whether these features can be due to correlations.

Figure 9 shows the scale length variation with age obtained with sim-2 using Eq. (15). The values of the fitted parameters are given in Table 6. As can be seen in this figure, the scale length clearly decreases with age. This result, which is related to the thin disc formation scenario, will be discussed in Sect. 6.

5.3. Disc truncation

As seen in Table 5 with sim-1, the disc truncation (R_dis) is found at approximately 16.1 ± 1.3 kpc. The scale (h_cut) of cutoff is equal to 0.72 ± 0.27 kpc. It is such that the density drops by 90% within 1 kpc, therefore, the cutoff is sharp. In sim-2, when age dependence for warp, flare and scale length is considered the disc truncation obtained is 19.4 ± 1.4 kpc and considering the ripples (sim-3, see Sect. 5.5) is 18.7 ± 1.6 kpc, still in agreement at the 1σ level with its value in sim-1.

We note that the disc truncation is found at a larger Galactocentric distance than in the standard model, where it was assumed 14 kpc. The determination is here more robust than in Robin et al. (1992a) from which the standard model was based, since only one direction was used for that determination.

5.4. Parameter correlations

In order to investigate whether the parameters are correlated, we computed the correlation coefficients in the optimisation procedure sim-1, shown in Table B.1, and in Fig. A.1. R_flare and γ_flare are fairly correlated because we have considered only low latitude fields in our study. The same problem was mentioned by López-Corredoira et al. (2002) and López-Corredoira & Molgó (2014). A significant correlation also appears between R_warp and the warp amplitude specially in the third quadrant.

Table B.2 shows the correlations among the parameters for sim-2. As can be seen, some parameters are correlated (notably the starting radius and the amplitudes of the warp and flare). This is expected because the photometric distances are not precise enough. Hence, the exact position where the flare starts is not well determined, and it probably starts more smoothly than a linear increase starting at a given point. But the overall measurements of the warp and flare dependencies with age are clear enough and not impacted by those correlations. To avoid these correlations, in the future several complementary studies should be considered. Firstly, looking at higher latitudes, covering latitudes from –30° to +30° in the anti-centre direction would allow disentangling the effect of the flare from the radial truncation. Secondly, the correlation between the starting radius and the amplitude is due to a simplistic linear shape used. Using samples with good distance estimates, such as Gaia data or red clump giants, would allow to determine better the shape of the starting radius and the amplitude of the flare.

5.5. Is there an improvement considering the ripples?

We also considered in another set of simulations (sim-3), the contribution of ripples, an oscillation in the number of star counts, as described by Xu et al. (2015) without fitting any parameter of the ripple itself. They have used SDSS data towards 110°<ℓ ≤ 229° in order to study the structures in the outer Galaxy. We have used the same equations for the northern and southern, as provided by Xu et al. (2015) in Sect. 6 of their paper.

Basically, when computing z_w we add the z_wripple given by Xu et al. (2015). The resulting parameters and the χ² of the solution are given in Table 6, showing no improvements compared with sim-2 but rather a slightly higher BIC value.

6. Discussion

6.1. Thin disc scale length

The disc scale length is the most important unknown in disentangling the contributions from the disc and the dark halo to the mass distribution of the Galaxy (Dehnen & Binney 1998; Bovy & Rix 2013). Recently, the thin disc scale length has been studied from photometric and spectroscopic large scale surveys. For example, Bovy et al. (2012) determine various scale lengths for mono-abundance populations. Assuming that the thin disc is defined by the low α abundance population, they show that the scale length can depend on the metallicity within the thin disc. In their sample high metallicity stars have a shorter scale length than low metallicity ones, with noticeably large error bars, which is in contradiction with our work if metallicity is anti-correlated with age. However the mean orbital Galactic radii of the low metallicity stars are much larger, which implies that their sample of low metallicity stars comes from the outer Milky Way, while the high metallicity sample comes from the inner Milky Way. In the end, the sample from the outer Milky Way has larger scale length, as expected from the inside-out scenario. The scale lengths of the mono-abundance populations of the thin disc range from 2.4 to 4.4 kpc, while we found a very similar scale length range of 2.3 to 3.9 kpc.

Our result of a short scale length for the main thin disc (as stars older than 3 Gyr dominates the counts in the Milky Way) is in general agreement with previous results. Table 7 shows most of the scale length determinations from different authors during the last twenty years as well as the tracer used and field coverage. It is noticeable that most studies using IR data found a short scale length for the thin disc, compared with longer scale length obtained from visible data. This is understandable if the IR tracers are in the mean older than the visible tracers (young stars emit mostly in the visible). However, this could also be due to the extinction which perturbs the analysis in the visible.

López-Corredoira et al. (2002) constrained the outer disc (scale length, warp and flare) from a study of 2MASS data in 12 fields but only four different longitudes (155°, 165°, 180° and 220°), which does not include fields where the warp is stronger (around 90° and 270°). They only use a selection of red giant stars from CMDs and they apply a simple modelling approach assuming a luminosity function and density law for the disc. They find a thin disc scale length of 1.91$\begin{array}{c}\mathrm{+}\mathrm{0.20}\\ -0.16\end{array}$ kpc, which is in excellent agreement with our determination for stars older than 3 Gyr, but smaller than our mean scale length. Their sample is dominated by red giants. Hence it does not include higher mass and younger stars, which might explain the difference.

Yusifov (2004) studied the distribution of pulsars from the Manchester et al. (2005) catalogue. Although Momany et al. (2006) claims that the catalogue is incomplete in the outer Galaxy, their scale length for these young objects is 3.8 ± 0.4 kpc, in close agreement with our young star sample.

As pointed above, Kalberla & Dedes (2008) by fitting HI distribution obtained a radial exponential scale length of 3.75 kpc in the mid-plane. Even if they do not give an estimate of their error, it is in excellent agreement with our very young star scale length of 3.90 ± 0.28 kpc.

For the first time, we report a clear evolution of the scale length with time among field stars. There have been previous evidences of young stellar associations and young open clusters at distances where the old stars seem to be no more present (Robin et al. 1992b). However, it was long known that the HI disc has a longer scale length than the stars in the mean. Hence, it is not too surprising to see that very young stars have a scale length similar to the gas.

6.2. Disc truncation

The disc truncation has been often observed in external galaxies. The existence of this cut-off can be related to a threshold in the star formation efficiency, when the gas density drops under certain value (Kennicutt 1989; van der Kruit 1979).

We found a radius for disc truncation equal to R_dis = 16.1 ± 1.3 kpc in the case of sim-1. This value increases up to 18 kpc in the case of sim-2 when the scale length is assumed to vary with time. We have not tested a possible dependence of the truncation on age, in order to avoid increasing too much the number of free parameters. We shall consider this point in the future by extending the fit to higher latitudes, in order to decrease the degeneracy between the flare and the truncation.

Habing (1988) from OH/IR stars and Robin et al. (1992a) from UBV photometry found the edge of Galactic disc located at 14–15 kpc and 14 kpc, respectively. Ruphy et al. (1996) using DENIS data similarly found 15 ± 2 kpc for the disc truncation distance. These values are in good agreement with our determination. Using photometry for IPHAS stars towards 160° ≤ ℓ ≤ 200° (| b | ≤ 1.0°), Sale et al. (2010) obtained a truncation radius of 13.0 ± 1.1 kpc, slightly smaller than ours, but still within the uncertainties. López-Corredoira et al. (2002) did not find a disc truncation at distances closer than 20 kpc, but indicated that they could not distinguish between a truncation or a flare.

6.3. Warp parameters and origin

Several scenarii have been proposed to explain the Galactic warps, the most cited being a dark halo-disc misalignment and possible interactions with nearby galaxies (or flyby galaxies). Perryman et al. (2014) stated that the misalignment of the disc inside the dark halo might produce changes in the warp angle with time. But in this case, all components of the Galaxy should be perturbed independently of their age. The warp angle dependence that we see could be consistent with this scenario if the time dependence reflects a precession. Alternatively, Weinberg & Blitz (2006) produced dynamical simulations of the flyby of the Magellanic Clouds around the Galaxy which could produce a strongly asymmetric warp varying with time.

However, Reshetnikov et al. (2016) analysed the global structure of 13 edge-on spiral galaxies using SDSS data. They pointed out that the warps found in those galaxies are generally slightly asymmetrical. They studied the relation of the strength and asymmetry of the warps with the dark halo mass. They showed that galaxies with massive halos have weaker and more symmetric warps, concluding that these dark halos play an important role in preventing strong and asymmetric warps. In the case of our Galaxy, the M_tot/M_∗ amount to about 10 (Robin et al. 2003) and the warp angle (computed as seen in edge-on galaxies) is less than 1°. Typical galaxies in Reshetnikov’s sample with this mass ratio have warp angles less than 10° and an asymmetry of less than 5°. Hence, our Galaxy compares well with these edge-on galaxies having a relatively large dynamical to stellar mass ratio and a weak warp, even though it is asymmetrical.

The Galactic warp can be observed in different components, such as dust, gas and stellar. Reylé et al. (2009) using BGM and 2MASS described the warp and flare, comparing their determination with those provided by several authors from different tracers, such as dust, gas and stars. In this work, they considered R_warp = 8.4 kpc, γ_warp = 0.09 pc kpc^-1 for a scale length equal to 2.2 kpc, the same value obtained by Derrière & Robin (2001). Their starting radius was a bit shorter than our mean value of 9.18 kpc.

They found evidence that the warp is asymmetric but they were not able to find a good S-shape for the warp in the third quadrant. In the present work, we have adjusted distinct values for the slope of the warp at second and third quadrants. We found a value larger than the one obtained by Reylé et al. (2009), with a slope in the second quadrant approximately three times larger than in the third quadrant. Analysing dust extinction distributions, Marshall et al. (2006) also determined that the warp in the second quadrant has a larger amplitude than in the third quadrant, γ_warp = 0.14 and 0.11 pc kpc^-1, respectively. They also pointed out that the warp seems to start earlier in the third quadrant. However, the two parameters (R_warp and γ_warp) are not completely independent when the warp model is assumed linear.

The difference can be well explained by the fact that we have varied other parameters which were not considered by these authors. They also did not consider the dependence with age, and in the present study we considered a larger range of longitudes, latitudes and limiting magnitude to constrain the model.

Drimmel & Spergel (2001) modelled the Galactic structure by using COBE/DIRBE data at near and far-infrared assuming a warp with a quadratic function z_warp(R) = 27.4(R−R_warp)²sin(φ), with R_warp = 7 kpc. López-Corredoira et al. (2002) model the warp with a different shape than ours. The amplitude of their warp is 2.1 × 10^-19 × R^5.25. It means that the mid-plane is going up to 1.2 kpc at R = 14 kpc, their node angle is found to be −5°± 5°. They are less sensitive to very young stars than we are because they select mainly red giants in their sample.

Yusifov (2004) models the warp that pulsars follow and find a warp node angle of 14.5°, at 1σ of our angle for youngest stars. Levine et al. (2006) present a very complex warp shape in HI, which extends at much larger distances, difficult to compare with our simple model.

Momany et al. (2006) studied the distribution of RGB stars to model the outer Galaxy and attempt to explain the Canis Major overdensity. They argue for no outer disc truncation, and that the presence of a strong warp and flare at longitudes of about 240° explain well the overdensity. The study of this particular sub-structure as well as the Monoceros ring is postponed to a future paper.

Figure 10 compares the amplitude of the warp from different authors with our result for three age class. A good agreement is obtained with Burton (1988), López-Corredoira et al. (2002) and Yusifov (2004) for stars with mean age equal to 2.5 Gyr in the distance range, 10.5 <R< 13.5 kpc. However, the warp slopes were determined using tracers in different Galactocentric distance ranges. For instance, López-Corredoira et al. (2002) and Yusifov (2004) measurements are claimed by their authors to be valid up to 15 and 18 kpc, respectively. There is also a good agreement between Drimmel & Spergel (2001) and our results for mean age equal to 6.0 Gyr for 9.0 <R< 11.5 kpc. In addition, our estimates of the warp for stars with mean age equal to 0.575 Gyr at negative z_w are in good agreement with Levine et al. (2006), at less than 1.5σ. The agreement is less good at positive values of z_w, as if the warp effect on stars was different from the HI warp.

We conclude that there is a general consistency between our results and previous ones. The difference between our estimates of the warp shape and other authors can be attributed to the different ranges of longitudes and latitudes covered in each study, as well as to the tracers and methods used. But we are first to claim a warp dependence with the age of the tracers, which is a clue for the scenario of warp formation.

6.4. How does the Galaxy flare in the outskirts

It has been long shown that the HI gas is flaring in the outer Galaxy. This is most probably related to the vertical force (Crézé et al. 1998; Holmberg & Flynn 2000; Kerr & Lynden-Bell 1986; Moni Bidin et al. 2012; Siebert et al. 2003; Sánchez-Salcedo et al. 2011), which is dominated in the inner Galaxy by the stellar disc but in the outer galaxy by dark matter. Hence a good characterization of the K_z and the flare in the outer Galaxy would lead to constrains on its dark matter content.

Kalberla & Dedes (2008) studied the HI distribution in the outer Milky Way and showed that the gas is distributed in two populations. The main gas layer goes to R ~ 35 kpc; the second one goes up to 60 kpc with a very shallow scale length of 7.5 kpc. The first and dominant component is flaring and lopsided. They explain the asymmetry by a dark matter wake towards the Magellanic Clouds. This is much shallower than the flare indicated (see below) by López-Corredoira et al. (2002).

Alard (2000) using 2MASS data studied the flare and warp for longitudes located at ℓ ~ 66.0°, ℓ ~ 180.0° and ℓ ~ 240.0° for | b | < 50° founding evidences for an asymmetry associated with the Galactic warp. He also argued that the flaring and warping seen in the stellar disc is very similar to the characteristics observed in HI disc. López-Corredoira et al. (2002) estimated the flare distance scale to be 4.6 ± 0.5 kpc, which gives an increase of the scale height by a factor of two at a Galactocentric distance of about 11 kpc and a factor of 10 at 18 kpc. Instead of an exponential, we are using a linear increase of the scale height with Galactocentric radius, which gives a factor of two at 12.2 kpc close to López-Corredoira et al. (2002), and a factor of 4 at about 20 kpc, a less extreme value.

Figure 11 (upper panel) shows a comparison of our flare factor (Eq. (2)) with the one of the Alard (2000), López-Corredoira et al. (2002), Yusifov (2004), Kalberla et al. (2014), López-Corredoira & Molgó (2014). To compute the flare factor and compare it with our results, we have divided the flare from other authors by h_z(R_Sun) when the flare factor was not given. The values of h_z(R_Sun) were taken from Eq. (4) of Kalberla et al. (2014), Eq. (3) of López-Corredoira & Molgó (López-Corredoira & Molgó (2014)), the equation quoted in the abstract of López-Corredoira et al. (2002) and Eq. (7) of Yusifov (2004). The values of h_z(R_Sun) used by each author are indicated in the caption of Fig. 11.

Concerning the estimations of errors on the flare factor, Yusifov (2004) mentioned an error around 30% on the pulsar distances determination, López-Corredoira et al. (2002) estimated random errors ranging from 5 to 7%, and López-Corredoira & Molgó (2014) estimated the error to be approximately 18%.

At small distances, all studies give very similar flare factors. In young stars (⟨ age ⟩ = 0.075 Gyr) our flare factor appears similar to López-Corredoira et al. (2002) flare up to 13.0 kpc, and the values are slightly larger for the older stars. For intermediate age stars (⟨ age ⟩ = 2.5 Gyr), a good agreement is seen for a wide range of Galactocentric radius (10.5 kpc ≤R< 16.0 kpc) in comparison with López-Corredoira & Molgó (2014). These authors determine the flare from the SEGUE imaging survey. They fitted a model to F8-G5 stars selected by colours in a sample which did not cover low latitudes | b | < 8°, contrarily to our study. They claimed to have found a significant flare but no outer disc truncation. Their flare factor amounts to 3.7 for thin disc stars between the Sun and a Galactocentric distance of 15 kpc, and a factor of 8.3 at a distance of 20 kpc. In comparison, our result points to a factor of 3.8 at 15 kpc and 5.8 at 20 kpc. We assumed a linear slope of the flare to limit the number of parameters. The agreement up to 15 kpc is remarkable. The difference for further distance is most probably due to the disc truncation that we see in-plane.

Yusifov (2004) has a flare scale length of 14 kpc which gives a shallow increase of the scale height with Galactocentric radius, shallower than the HI flare. However, their sample is small (1600 pulsars) and probably incomplete in the outer Galaxy. No estimate of the error bars was given. Hammersley & López-Corredoira (2011) by using five SDSS fields with | b | > 11° found a smaller amplitude for the flare than obtained, for instance, by Yusifov (2004).

Recently, Feast et al. (2014) identified a few Cepheids at a large distance from the Galactic plane. These Cepheids are a clear example that the young disc is also flaring.

Minchev et al. (2015) attempt to study the effect of the flare on the radial gradient of age, that can be seen in the “geometrical” thick disc (as defined by the population located at some distance from the Galactic plane, where one expects to be dominated by thick disc stars). They show from numerical simulations that a radial gradient in age naturally emerges due to the flaring of the populations, if younger populations are more extended than older ones. Martig et al. (2016) by using APOGEE data reinforces the results found by Minchev et al. (2015). In our study, the younger disc populations are indeed more extended than older ones, and they also flare earlier. Hence it is expected to produce such age radial gradient at a few kiloparsec distances from the plane. We have not tested this yet, but it will be considered in a future study.

7. Conclusions

We have studied the outer Galaxy structure by using the Besançon Galaxy Model with 2MASS data in order to constrain external disc parameters, such as warp and flare shape, scale length and disc truncation parameters. After the parameter optimisation, the relative difference between the model and 2MASS is generally 10–20% at most, which includes the cosmic variance of the star counts as well as the effect of patchy extinction.

We show strong evidence that the thin disc scale length, as well as the warp and flare shapes, changes with age, and proposed expressions for these dependencies. These results impact directly our comprehension about, not only the shape of those components, but also their origin.

The warp can come from misalignment of the disc inside the dark halo, which might change with time (precession). Alternatively, it can be due to the interaction with the Magellanic Clouds. In both cases, it can imply that tracers of different ages show different spatial distribution, either because there is precession, or because they react differently to the perturbation. The flare could be caused by the fact that the intergalactic gas (or gas coming from Galactic fountains) is falling into the disc more slowly and on longer time scales in the outskirts than in the inner Galaxy. We show in the present work some evidence for their dependence with age that reinforces the point that the warp has a dynamical origin as also demonstrated by López-Corredoira et al. (2007), among others.

Our results clearly show the variation of scale length with the age. The larger scale length for youngest stars (3.9 kpc) is well in agreement with the values found in the HI gas, while the shorter value (2 kpc) for the oldest thin disc is similar to the one of the thick disc. This also directly affects our comprehension of the history of Galactic formation and evolution, reinforcing the idea of an inside-out process of formation.

We also found a disc truncation with R_dis = 16.1 ± 1.3 kpc. But when we allow the parameters to vary with time, the disc truncation is no more clear and could be at larger distances than 19 kpc, a distance where the stars are scarce anyway, simply due to the exponential fall off. However, this value is not completely independent of other parameters such as the flare and disc scale length.

A further study involving higher latitudes will be considered to strengthen the conclusion about the disc truncation and its dependence with age.

---

Acknowledgments

We thank Misha Haywood, Paola Di Matteo and Vanessa Hill for fruitful discussions about this work. We also thank the anonymous referee for useful, valuable, and detailed comments on the manuscript. Eduardo Amôres would like to thank very much all staff of the Observatoire de Besançon for the hospitality and courtesy during his sejour and collaboration visits. We also thanks Dr. Douglas Marshall for let available the 3D extinction model prior to its publication. BGM simulations were executed on computers from the UTINAM Institute of the Université Bourgogne Franche-Comté, supported by the Région de Franche-Comté and Institut des Sciences de l’Univers (INSU). Financial support for this work was provided by Conselho Nacional para o Desenvolvimento Científico e Tecnológico (CNPq), the Coordenacão de Aperfeiçoamento Pessoal de Nível Superior (CAPES) and the CNRS. Eduardo Brescansin de Amôres had a CNPq fellowship (2003/201113-3) and also thanks INCT−A (Brazil). This publication made use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.

References

Appendix A: Parameters correlations

Appendix B: Additional tables

All Tables

All Figures