Issue 
A&A
Volume 637, May 2020



Article Number  A96  
Number of page(s)  17  
Section  Galactic structure, stellar clusters and populations  
DOI  https://doi.org/10.1051/00046361/201937289  
Published online  27 May 2020 
Structure of the outer Galactic disc with Gaia DR2
^{1}
Instituto de Astrofísica de Canarias, 38205 La Laguna, Tenerife, Spain
email: zofiachrobakova@gmail.com
^{2}
Departamento de Astrofísica, Universidad de La Laguna, 38206 La Laguna, Tenerife, Spain
^{3}
Faculty of Mathematics, Physics, and Informatics, Comenius University, Mlynská dolina, 842 48 Bratislava, Slovakia
Received:
9
December
2019
Accepted:
31
March
2020
Context. The structure of outer disc of our Galaxy is still not well described, and many features need to be better understood. The second Gaia data release (DR2) provides data in unprecedented quality that can be analysed to shed some light on the outermost parts of the Milky Way.
Aims. We calculate the stellar density using star counts obtained from Gaia DR2 up to a Galactocentric distance R = 20 kpc with a deconvolution technique for the parallax errors. Then we analyse the density in order to study the structure of the outer Galactic disc, mainly the warp.
Methods. In order to carry out the deconvolution, we used the Lucy inversion technique for recovering the corrected star counts. We also used the Gaia luminosity function of stars with M_{G} < 10 to extract the stellar density from the star counts.
Results. The stellar density maps can be fitted by an exponential disc in the radial direction h_{r} = 2.07 ± 0.07 kpc, with a weak dependence on the azimuth, extended up to 20 kpc without any cutoff. The flare and warp are clearly visible. The best fit of a symmetrical Sshaped warp gives z_{w} ≈ z_{⊙} + (37 ± 4.2(stat.) − 0.91(syst.))pc ⋅ (R/R_{⊙})^{2.42 ± 0.76(stat.) + 0.129(syst.)}sin(ϕ + 9.3° ±7.37° (stat.) + 4.48° (syst.)) for the whole population. When we analyse the northern and southern warps separately, we obtain an asymmetry of an ∼25% larger amplitude in the north. This result may be influenced by extinction because the GaiaG band is quite prone to extinction biases. However, we tested the accuracy of the extinction map we used, which shows that the extinction is determined very well in the outer disc. Nevertheless, we recall that we do not know the full extinction error, and neither do we know the systematic error of the map, which may influence the final result. The analysis was also carried out for very luminous stars alone (M_{G} < −2), which on average represents a younger population. We obtain similar scalelength values, while the maximum amplitude of the warp is 20 − 30% larger than with the whole population. The northsouth asymmetry is maintained.
Key words: Galaxy: disk / Galaxy: structure
© ESO 2020
1. Introduction
Studying the Galactic structure is crucial for our understanding of the Milky Way. Star counts are widely used for this purpose (Paul 1993), and the importance of this tool has increased in the past decades with the appearance of widearea surveys (Bahcall 1986; Majewski 1993), which made it possible to obtain reliable measurements of the Galactic thin and thickdisc and halo (Chen et al. 2001; Jurić et al. 2008; Bovy et al. 2012; Robin et al. 2012). It is common to simplify the Galactic disc as an exponential or hyperbolic secant form, but there are many asymmetries such as the flare and warp that need to be taken into account. These structures can be seen from 3D distribution of stars, as shown by Liu et al. (2017), who mapped the Milky Way using the LAMOST (The Large Sky Area MultiObject Fibre Spectroscopic Telescope) RGB (redgiant branch) stars; Skowron et al. (2019a), who constructed a map of the Milky Way from classical Cepheids; or Anders et al. (2019), who used the second Gaia data release (DR2).
The warp was first detected in the Galactic gaseous disc in 21 cm HI observations (Kerr 1957; Oort et al. 1958). Since then, the warp has also been discovered in the stellar disc (Carney & Seitzer 1993; LópezCorredoira et al. 2002a; Reylé et al. 2009; Amôres et al. 2017; Chen et al. 2019), and the kinematics of the warp has been studied as well (Dehnen 1998; Drimmel et al. 2003; LópezCorredoira et al. 2014; Schönrich & Dehnen 2018).
Vertical kinematics in particular can reveal much about the mechanism behind the formation of warp. Poggio et al. (2018) found a gradient of 5 − 6 km s^{−1} in the vertical velocities of upper mainsequence stars and giants located from 8 to 14 kpc in Galactic radius using Gaia DR2 data, revealing the kinematic signature of the warp. Their findings suggest that the warp is principally a gravitational phenomenon. Skowron et al. (2019b) also found a strong gradient in vertical velocities using classical Cepheids supplemented by the OGLE (Optical Gravitational Lensing Experiment) survey. LópezCorredoira et al. (2020) investigated the dynamical effects produced by different mechanisms that can explain the radial and vertical components of extended kinematic maps of LópezCorredoira & SylosLabini (2019), who used Lucy’s deconvolution method (see Sect. 4.1) to produce kinematical maps up to a Galactocentric radius of 20 kpc. LópezCorredoira et al. (2020) found that vertical motions might be dominated by external perturbations or mergers, although with a minor component due to a warp whose amplitude is evolving with time. However, the kinematic signature of the warp is not enough to explain the observed velocities.
To date, the shape of the warp has been constrained only roughly, and the kinematical information is not satisfying enough to reach consensus about the mechanism causing the warp. Theories include accretion of intergalactic matter onto the disc (LópezCorredoira et al. 2002b), interaction with other satellites (Kim et al. 2014), the intergalactic magnetic field (Battaner et al. 1990), a misaligned rotating halo (Debattista & Sellwood 1999), and others.
We now have a new opportunity to improve our knowledge about the Milky Way significantly through the Gaia mission of the European Space Agency (Gaia Collaboration et al. 2016). Gaia data provide unprecedented positional and radial velocity measurements and an accurate distance determination, although the error of the parallax measurement increases with distance from us. It brings us the most accurate data about the Galaxy so far, ideal to advance in all branches of Galactic astrophysics and study our Galaxy in greater detail than ever before. Gaia DR2 has been used by Anders et al. (2019), who provided photoastrometric distances, extinctions, and astrophysical parameters up to magnitude G = 18, making use of the Bayesian parameter estimation code StarHorse. After introducing the observational data and a number of priors, their code finds the Bayesian stellar parameters, distances, and extinctions. The authors also present density maps, which we compare with our results in Sect. 4.3. Gaia data have also been used to study the structure of outer Galactic disc, especially the warp and the flare. The first Gaia data release brought some evidence of the warp (Schönrich & Dehnen 2018), but the more extensive second data release provides a better opportunity to study the warp attributes. Poggio et al. (2018) combined Gaia DR2 astrometry with 2MASS (Two Micron AllSky Survey) photometry and revealed the kinematic signature of the warp up to 7 kpc from the Sun. Li et al. (2019) found the flare and the warp in the Milky Way, using only OB stars of the Gaia DR2. In this work, we make use of Gaia DR2 data as described in Sect. 2 and use star counts to obtain the stellar density by applying Lucy’s inversion technique. Then we analyse the density maps to determine the warp.
The paper is structured as follows: in Sect. 2 we describe the Gaia data and extinction maps that we used, in Sect. 3 we present the luminosity function used in our calculations, in Sect. 4 we explain the methods for obtaining our density maps, and in Sect. 5 we discuss the results. In Sect. 5.4 we present the exponential fits of the density, in Sect. 5.5 we study the warp, and in Sect. 5.6 we repeat the previous analysis of the young population.
2. Data selection
We used data of the second Gaia data release (Gaia Collaboration et al. 2018) here, which were collected during first 22 months of observation. We are interested in stars with known fiveparameter astrometric solution: more than 1.3 billion sources. G magnitudes, collected by astrometric instrument in the whitelight Gband of Gaia (330–1050 nm) are known for all sources, with precisions varying from around 1 millimag at the bright (G < 13) end to around 20 millimag at G = 20. For the details on the astrometric data processing and validation of these results, see Lindegren et al. (2018). We chose stars with apparent magnitude up to G = 19, where the catalogue is complete up to 90% (Arenou et al. 2018). We chose data with a parallax in the interval Anders et al. (2019) mas.
In our analysis, we did not consider any zeropoint bias in the parallaxes of Gaia DR2, as found by some authors (Lindegren et al. 2018; Arenou et al. 2018; Stassun & Torres 2018; Zinn et al. 2019), except in Sect. 4.3, where we repeat our main calculation including a nonzero value of the zeropoint to prove that this effect is negligible in our results.
We used two different extinction maps. For the luminosity function (Sect. 3), we used the extinction map of Green et al. (2018) through its Python package dustmaps, choosing the Bayestar17 version. This map covers 75% of the sky (declinations of δ ≳ −30°) and provides reddening in similar units as Schlegel et al. (1998, SFD).
To calculate the density (Sect. 4), we need to cover the whole sky, therefore we used the threedimensional less accurate but fullsky extinction map of Bovy et al. (2016) through its Python package mwdust. This map combines the results of Marshall et al. (2006), Green et al. (2015), and Drimmel et al. (2003) and provides reddening as defined in Schlegel et al. (1998).
In order to convert the interstellar reddening of these maps into E(B − V), we used coefficients (Hendy 2018; Rybizki et al. 2018)
3. Luminosity function
To construct the luminosity function, we chose all stars with heliocentric distance d < 0.5 kpc (distances determined as 1/π, where π is the parallax). We did not find many bright stars (M_{G} < −5) in this area, therefore we also chose a specific region with Galactic height z < 1 kpc and Galactocentric distance R < 5 kpc, in which we only selected stars with absolute magnitude M_{G} < −5. We normalised the counts of stars with high magnitude and then joined these two parts to create the luminosity function.
In the range of distance that we used for the luminosity function, the star counts are complete for the absolute magnitude that we are calculating, except perhaps for the possible loss of the brightest stars through saturation at M_{G} < −5. Moreover, the error in the parallax for these stars is negligible, so that the calculation of the absolute magnitude from the apparent magnitude is quite accurate. We did not take into account the variations of the luminosity function throughout the Galactic disc. We assumed that it does not change.
The luminosity function we obtained is shown in Fig. 1. We interpolated the luminosity function with a spline N = spl(M) of the first degree. The result is shown in Fig. 2. For the interpolation, we used values between magnitudes M = [ − 5, 10] because the values outside this interval are unreliable, and we used the extrapolation of the spline function to lower magnitudes. The values of the luminosity function are listed in Table 1.
Fig. 1. Luminosity function. 
Fig. 2. Interpolation of the luminosity function with a spline compared with the luminosity function of Bahcall & Soneira (1980). These two functions are not directly comparable because Bahcall & Soneira (1980) used a slightly different filter in the visible, but it shows that our luminosity function is reasonable. 
Values of the luminosity function.
4. Density maps
4.1. Deconvolution of star counts
To calculate the stellar density, we need to measure star counts as a function of distance. However, the error of parallax increases with distance from us, which means that our analysis would be correct only within roughly 5 kpc from the Sun. To be able to reach higher distances, we corrected for this effect using the method developed by LópezCorredoira & SylosLabini (2019), who used Lucy’s deconvolution method (Lucy 1974; see Appendix A) to obtain an accurate distance measurement up to R = 20 kpc. They expressed the observed number of stars per parallax as a convolution of the real number N(π) of stars with a Gaussian function
where
For the error σ_{π} we averaged errors of every bin, which we calculated from values given by Gaia DR2.
We only used the parallax between Anders et al. (2019) mas. For the upper limit the relative error of parallax is very small and does not produce any bias. For the lower limit, the truncation avoiding the negative parallaxes affects the distribution of parallaxes and statistical properties (average, median, etc.) (Luri et al. 2018, Sect. 3.3). However, in our method we do not calculate the average distance from the average parallax. We used Lucy’s method, which iterates the counts of the stars with positive parallaxes until we obtained the final solution. This does not mean that we truncated the star counts with negative parallaxes. We used only the stars with positive parallaxes as is required by our method, explained in the Appendix A. N(π) for negative values of π can also be calculated and fitted, but they are not used in our calculation. In other words, we did not assume that the number of the stars with negative parallaxes is zero, we simply did not use this information because it is not necessary. The fact that this method does not produce any bias is tested in Sect. 4.3.
4.2. Monte Carlo simulation to test the Lucy inversion method
In order to test the reliability of the inversion method, we performed Monte Carlo simulations to determine whether we can recover the original function after deconvolution. We created datasets with randomly distributed particles. Then we convolved this distribution with a Gaussian. We applied Lucy’s deconvolution method to the dataset to determine whether we can recreate the original distribution. The results are shown in Fig. 3. We conclude that regardless of the original distribution, we can accurately recover the original data up to 50 kpc or more, which is satisfying to study the Milky Way. We also studied the dependence of the method on the parallax error. We used various values of the average parallax error in Eq. (3) from the interval [0.05, 0.4] mas, which are the most common values for the average parallax error in our data. In Fig. 4 we plot the result, which shows that even though the precision of the method depends on the parallax error, we obtain a satisfying result up to 20 kpc even in the worst case with the highest parallax error.
Fig. 3. Monte Carlo simulation of deconvolution. We recover random distributions, convolved with a Gaussian. (a) Gamma distribution; (b) exponential distribution; (c) geometric distribution; (d) logarithmic distribution. 
Fig. 4. Monte Carlo simulation of deconvolution. In all cases we recover a gamma distribution convolved with a Gaussian, varying the average error of parallax. (a) σ_{π} = 0.05 mas; (b) σ_{π} = 0.1 mas; (c) σ_{π} = 0.25 mas; (d) σ_{π} = 0.4 mas. 
4.3. Application to fullsky GaiaDR2 data
We divided the data into bins of Galactic longitude ℓ, Galactic latitude b, and apparent magnitude m. For the values of b we made bins of length 2° and corresponding ℓ in bins of 5° /cos(b). We divided each of the lines of sight in magnitude, binned with size Δm = 1.0 between G = 12 and G = 19. We obtained 29 206 different areas in which we calculated the density independently.
We made use of the fundamental equation of stellar statistics, where the number of stars N(m) of apparent magnitude m is expressed per unit solid angle and per unit magnitude interval (Chandresekhar & Münch 1951),
where we substitute
which yields for the density
where ω is the covered angular surface (10 degrees^{2} in our case), Δπ is the parallax interval (0.01 mas in our case), which must be added in the equation because we did not use the unit parallax, Φ(M_{G}) is the luminosity function in the G filter, m_{G, low lim} is the limiting maximum apparent magnitude, and A_{G}(r) is the extinction, as a function of distance.
After this, we calculated the weighted mean density for all seven ranges of magnitude in each line of sight. Then we transformed this into cylindrical coordinates and made bins of Galactocentric radius R of length 0.5 kpc, in Galactic height z of 0.1 kpc and in azimuth of 30°. We define the azimuthal angle ϕ to be measured from the centreSunanticentre direction towards the Galactic rotation, going from 0° to 360°. We interpolated the missing bins with NearestNDInterpolator from the python SciPy package, which uses nearestneighbour interpolation in N dimensions. We plot the resulting density maps in Figs. 5 and 6. In Fig. 5 we plot the density in cylindrical coordinates as a function of Galactic radius R for different azimuths. We do not plot the results for azimuths 90° < ϕ < 270° because in this area the extinction is significant and we observe stars farther than the Galactic centre for which the errors are too large, therefore we cannot see any structure in density. However, we can see even by eye that a northern warp is present in the azimuths 60° < ϕ < 90° and a southern warp in the azimuths 270° < ϕ < 300°. Another structure that can be seen from the plots is the flaring of the disc. We analyse these structures below. In Figs. 6a–c we plot the density map in Cartesian coordinates, and in Fig. 6d we plot the density in cylindrical coordinates, integrated through all ranges of azimuths, except for the areas that were excluded from the analysis. The Cartesian coordinates are defined such that X_{⊙} = 8.4 kpc. In these plots we note a flat disc with some fluctuations in density, but no apparent features. However, some slights overdensities both above and below the Galactic plane are visible. The features above the plane are present only in Figs. 6b–c, but not in Fig. 6d, which suggests that it might be a contamination. The feature below the Galactic plane is present in all the three plots. As the direction of these overdensities is towards the Magellanic Clouds, it might be an effect of the Milky Way pulling stars out of Magellanic Clouds, as suggested by Anders et al. (2019). Another possible explanation for these overdensities is the finger of God artefact, which is caused by the foreground dust clouds and causes elongated overdensities that point to the Sun. This artefact has previously been seen in Gaia data, as shown in the Gaia DR2 documentation^{1}.
Fig. 5. Density maps for various azimuths between 0° and 360°. (a) 0° < ϕ < 30°; (b) 30° < ϕ < 60°; (c) 60° < ϕ < 90°; (d) 270° < ϕ < 300°; (e) 300° < ϕ < 330°; (f) 330° < ϕ < 360°. 
Fig. 6. Density maps. 
4.4. Zeropoint correction in parallaxes
So far, we did not consider any zeropoint bias in parallaxes. Lindegren et al. (2018) found a global mean offset of −0.029 mas, meaning that Gaia DR2 parallaxes are lower than the true value. We repeated our calculations with this correction and present the results in Fig. 7, where we chose some of the lines of sight to show the comparison. We find that these results are very similar to our original results, and this correction brings a negligible effect. We also tried a value of −0.046 mas, found by Riess et al. (2018). In Fig. 7 we show that the difference between the different zeropoint values is very small, therefore we only use the value of −0.029 mas in the further calculations.
Fig. 7. Comparison of densities for various lines of sight. Orange and green curves represent the density including the zeropoint correction of the parallax, and the blue curve shows the density without this correction. 
For the analysis of the warp in Sects. 5.5 and 5.6, we repeated the analysis of Sect. 4 with the value of parallax corrected for the zeropoint. We find that this brings a small correction to the warp parameters, which we state as the systematic error in the results.
4.5. Error of the extinction
To test how accurate the extinction map is, we analysed the map of Green et al. (2015) using the function query, which returns the standard deviation σ_{G} for a given line of sight. We calculate a new extinction as
where A_{G} is the extinction given by the map, r is the distance, and f is a factor chosen randomly from a Gaussian distribution with μ = 0 and σ = 1.
In Fig. 8 we show the relative error of the density . For all lines of sight that we tested, the difference is negligible, except for the area in the centre of the Galaxy, which we know is problematic. However, in the outer disc, where we carried out our analysis, the extinction is determined quite accurately. We must of course take into account that we used the map of Bovy et al. (2016), which combines different maps and is less accurate and therefore can give different results than Green et al. (2015) in some areas. Moreover, we estimated only the statistical error of the extinction, but we recall that we do not have information about the systematic error of the extinction map. However, for our purposes, the extinction map gives satisfying results in the area we analysed. The stellar warp has been studied using star counts by many other authors (LópezCorredoira et al. 2002a; Reylé et al. 2009; Amôres et al. 2017, and others), therefore this method is most likely not especially flawed.
Fig. 8. Relative error of the density calculated including the standard deviation of the extinction. 
4.6. Thickdisc areas
In the previous analysis we considered only the thindisc population because the luminosity function presented in Sect. 3 is calculated in thindisc regions. However, we can also analyse high Galactic heights, where the influence of the thick disc is significant. To test the importance of the change in luminosity function, we tested the density calculations with a tentative thickdisc luminosity function that reduces the number of bright stars. We used the source table of Wainscoat et al. (1992), who give the ratio of all the components of the Galaxy for all stellar classes. Based on this comparison, we altered our luminosity function to construct a theoretical thickdisc luminosity function, as depicted in Fig. 9. Then we repeated our calculation with this new luminosity function. In Fig. 10 we show the result for some lines of sight. In the area where we carried out the analysis, the difference between the two approaches is clearly visible starting at ∼20 kpc. Our density analysis is made in the area below 20 kpc, where the difference between the two densities is negligible. We note that this difference changes with line of sight, which is caused by the extinction. In the areas where the extinction is significant, the difference between densities derived from thin and thickdisc luminosity functions is more important, but these areas are removed from our analysis. Therefore our maps are also valid for thickdisc areas.
Fig. 9. Comparison of luminosity functions for the thin and thick disc. 
Fig. 10. Densities for different lines of sight. Blue curves are densities calculated with the thindisc luminosity function. Orange curves are densities calculated with the thickdisc luminosity function. 
5. Analysis of the density maps
5.1. Comparison with the maps of Anders et al. (2019)
Recently, similar maps were created by Anders et al. (2019). In their analysis, they used the code StarHorse, originally developed to determine stellar parameters and distances for spectroscopic surveys (Queiroz et al. 2018). This code compares observed quantities to a number of stellar evolutionary models. It finds the posterior probability over a grid of stellar models, distances, and extinctions. To do this, it needs many priors, including stellar initial mass function, density laws for main Milky Way components, and the broad metallicity and age of those components. Afterwards, the authors applied various criteria on their sample to choose only accurate results.
When we compare our results, we can observe similar structures, except in the area of the Galactic bulge, where our data are not reliable and the data of Anders et al. (2019) are much more accurate. However, because data with high errors in parallax were removed, Anders et al. (2019) were unable to reach such high distances, which are necessary to study features of the outer disc such as the flare or the warp. Another advantage of our method is that we did not assume any priors about the Milky Way. Furthermore, our density maps are a representation of the complete number of stars per unit volume up to some given absolute magnitude, taking into account the luminosity function, whereas Anders et al. (2019) gave the stars observed by Gaia, a much larger number in the solar neighbourhood, thus not useful to quantify absolute trends in the density distribution. Nevertheless, we consider the results of Anders et al. (2019) very useful because they improve the accuracy of the data significantly and can be used to study parts of Milky Way where our data fail.
5.2. Cutoff in the Milky Way
There has been some discussion about the cutoff in the Milky Way. Some authors have reported to find a cutoff starting at about 14 kpc from the Galactic centre (Robin et al. 1992, 2003; Minniti et al. 2011). However, Carraro et al. (2014) argued that these finding are erroneous either because the dataset is biased or because the warp and flare is confused with the cutoff. The absence of the cutoff has been confirmed by several studies (LópezCorredoira & Molgó 2014; Sale et al. 2010; Brand & Wouterloot 2007). Our results show that there is no cutoff in the Galactic disc, at least up to 20 kpc.
5.3. Stellar density in the solar neighbourhood
We define the solar neighbourhood as the area where 7.5 kpc < R < 8.5 kpc and z < 0.05 kpc and calculate the average density in this area. We find ρ_{⊙} = 0.064 stars pc^{−3}, which is close to other values in the literature, for example 0.03 stars pc^{−3} obtained by Chang et al. (2011), who used a threecomponent model to fit data from 2MASS. Eaton et al. (1984) found ρ_{⊙} = 0.056 stars pc^{−3}, which is lower than our result, but this value is influenced by the range of the luminosity function, which is where the difference between the values stems from. In our case, we measured stars with M_{G} < 10.
5.4. Exponential fits of the density
To describe the radial volume mass density distribution in the Galactic equatorial plane, we used a modified exponential disc with a deficit of stars in the inner inplane region adopted from LópezCorredoira et al. (2004) in the following form:
where h_{r} is the scale length, h_{r, hole} = 3.74 kpc is the scale of the hole, R_{⊙} is the Galactocentric distance of the Sun, and R is the Galactocentric distance. We neglected the contribution of the thick disc and analysed only the thin disc. We divided the Galactic equatorial plane into three regions according to the Galactic azimuth [−45°,−15°], [−15°,15°], [15°,45°]. We focused on the Galactic equatorial plane, therefore we considered stars in the close vicinity of the plane with a vertical distance z < 0.2 kpc and R > 6 kpc. We fitted the density for various azimuths with the corresponding exponential fits based on Eq. (9). The scale length slightly depends on the Galactic azimuth; it reaches the highest value for the Sunanticentre direction and ϕ = +30°, h_{r} = 2.78 ± 0.13 kpc, and h_{r} = 2.29 ± 0.21 kpc. On the other hand, the lowest value of the scale length is h_{r} = 1.88 ± 0.12 kpc for ϕ = −30°. This results in an average of h_{r} = 2.29 ± 0.08 kpc, with small dependence on the azimuth. We can compare the results with published papers. LópezCorredoira & Molgó (2014) used SDSSSEGUE (Sloan Digital Sky Survey – Sloan Extension for Galactic Understanding and Exploration) data to investigate the density distribution in the Galactic disc. They obtained the scale length for the thick and for the thin disc, h_{r, thin} = 2.1 kpc and h_{r, thick} = 2.5 kpc for the azimuth ϕ ≤ 30°, which is consistent with our results. Li et al. (2019) studied OB stars using Gaia DR2 data and the derived scale length of the Galactic disc, and found h_{r} = 2.10 ± 0.1 kpc, which is in accordance with our results.
We also plot the dependence of the density in the Galactic equatorial plane on azimuth for various values of Galactocentric distance in Fig. 11. The density is slightly dependent on the Galactic azimuth for all radii, but this dependence is very small. An analysis of the scale height and its corresponding flare will be given in a forthcoming paper (Nagy et al., in prep.).
Fig. 11. Dependence of the density on azimuth near the centreSunanticentre direction for various values of Galactocentric distance. The data points are obtained as weighted mean in bins of size 1 kpc in R and 0.4 kpc in z. Only bins with a number of points N ≥ 50 points are plotted. 
5.5. Warp
The density maps (Fig. 5) directly show a northern warp in azimuth 90° and southern warp in azimuth 270°. Here, we analyse these structures in greater detail. We removed the azimuths 150° < ϕ < 240° and radii R < 6 kpc from our analysis because these data have low quality and influence the results negatively.
We calculated the average elevation above the plane z_{w} as
and fit this quantity with models of the warp. In our first approach, we used the model by LópezCorredoira et al. (2002a, Eq. (20)),
The 17 pc term compensates for the elevation of the Sun above the plane (Karim & Mamajek 2017). C_{w}, ϵ_{w}, and ϕ_{w} are free parameters of the model, which were fitted to our data. An asymmetry is observed between the northern and southern warp for the gas (Voskes & Butler Burton 2006) and for the young population (Amôres et al. 2017), therefore we also explore the northern and southern warp separately here. The fit of our data yields maximum amplitudes z_{w} = 0.317 kpc for the northern and z_{w} = −0.287 kpc for the southern warp, both at a distance R = [19.5, 20] kpc, revealing a small asymmetry between the north and south. For the fit, we used the function curve fit from the python SciPy package, which uses nonlinear least squares to fit a function to data. The parameters of the best fit for this model for the whole dataset are
Here, the error of c_{w} stands for the error of the amplitude alone, without the variations of ϵ_{w} and ϕ_{w}. The plot of the results is shown in Fig. 12, where we show the comparison of minimum and maximum value of z_{w}(R). The average elevation of the plane is highest for azimuths [60° ,90° ] and [90° ,120° ] in most of the cases, whereas the minimum is reached for azimuths [240° ,270° ] in most of the cases. A slight asymmetry between the northern and southern warp is also clearly visible.
Fig. 12. Minimum and maximum average elevation of the plane as a function of radius. The warp fit is based on Eq. (11), and the error bars represent the uncertainty in the distance in the Lucy method. 
Another approach that we used is based on the work of Levine et al. (2006), who studied the vertical structure of the outer disc of the Milky Way by tracing neutral hydrogen gas. They analysed the Galactic warp using a Lomb periodogram analysis. They concluded that the first two Fourier modes are the strongest modes. We use the expression derived by Levine et al. (2006) in the following form:
where z_{w} is the average elevation above the plane, z_{i} for i ∈ (0,1,2) are the amplitudes of the warp, ϕ_{i} for i ∈ (1,2) are the phases. The dependence of the amplitudes of the warp on the Galactocentric distances is
where k_{i} and R_{k} are free parameters of the fit. We fitted our data with Eqs. (13) and (14) for various values of Galactocentric distances R < 20 kpc. We plot the data and the fits for R ∈ (13.25,16.25,19.25) kpc in Fig. 13. Figure 14 shows the azimuth of the maximum and minimum of the Galactic warp as a function of the Galactocentric distance. In our analysis, we excluded data for the Galactic azimuths ϕ ∈ (120°,240°) because of the high error values in our data. We used a 10° binning in azimuth. Figure 13 shows that the data for 250° < ϕ < 270° are somewhat noisy, which can be caused by problems with extinction or with the Lucy method in a particular line of sight. Therefore we tested a fit without these points, which turned out to produce an insignificant difference. For instance, the minimum amplitude obtained without these points changed by 10% in the worse case, and the maximum amplitude changed by 2%.
Fig. 13. Average elevation of the plane as a function of the Galactic azimuth for various values of the Galactocentric distance. Red markers represent values of binned data, and the blue line represents a fit to the data. (a) R = 13.25 kpc; (b) R = 16.25 kpc; (c) R = 19.25 kpc. 
Figures 13 and 14 clearly show that the warp is present in our analysis. The azimuth of the maximum of the warp (the northern warp) is an increasing function of the Galactocentric distance (52° < ϕ < 56°). On the other hand, the azimuth of the minimum of the warp is in 312° < ϕ < 324° and corresponds to the southern warp. The strongest deviation of the average elevation of the Galactic plane from the Galactic equatorial plane rises with Galactocentric distance. The highest amplitude of the northern and southern warp is z_{w} = 0.48 kpc and z_{w} = −0.38 kpc, respectively. An asymmetrical warp is clearly present.
Fig. 14. Minimum and maximum of the average elevation of the plane as a function of Galactocentric distance. The warp fit is based on Eq. (13). The colours code the azimuth of the minimum and maximum of the warp fit, and the error bars represent the uncertainty on the distance in the Lucy method. 
The value of the line of nodes from the fit is ϕ_{0} = −1.18°. We plot the changes in amplitude of the Galactic warp fit [Eq. (14)] with Galactocentric distance in Fig. 15.
Similar results were obtained by Li et al. (2019), who used OB stars of Gaia DR2 to measure the warp. They fit their data with a sinusoidal function similar to ours and obtained a warp with a mean magnitude up to z = 0.5 kpc. However, they did not account for the asymmetry of the warp, therefore they found the same result for the north and south. Chen et al. (2019) used Cepheids from the WISE (Widefield Infrared Survey Explorer) catalogue and traced the warp up to R = 20 kpc. Their results show a warp extended up to z = 1.5 kpc, which we cannot confirm using the whole population. Poggio et al. (2018) studied the kinematics of the Milky Way using Gaia DR2 and found the warp up to 7 kpc from the Sun. This agrees with our results, but we show that the warp extends to a higher radius, at least up to 20 kpc. In Fig. 16, we compare the maximum amplitudes of our model with other works. We obtain a very low amplitude, especially in comparison with Cepheids. On the other hand, the closest result is that of Chen et al. (2019), who used OB stars from Gaia DR2. This significant difference between the amplitude of various populations is in favour of the formation of the warp through accretion onto the disc (LópezCorredoira et al. 2002b), which causes the gas and young stars to warp more strongly than the remaining population.
Momany et al. (2006) studied the stellar warp using 2MASS red clump and red giant stars, selected at fixed heliocentric distances of 3, 7, and 17 kpc. They found a rather symmetric warp and argued that a symmetric warp can be observed as asymmetric for two reasons. First, the Sun is not located at the line of nodes, and second, the northern warp is located behind the NormaCygnus arm, which can cause variation in extinction that can produce an apparent asymmetric warp. As for the first point, the position of Sun on the line of nodes is a problem when we observe the warp at a fixed distance. However, we have a 3D distribution, which ensures that the position from which we look does not influence how we perceive the warp. As for the second remark, as we showed in Sect. 2 that the extinction is determined quite accurately by the extinction map of Green et al. (2015). However, some variations might influence the final shape of the warp and may not have been taken into account, therefore we need to keep that in mind when we interpret our results.
5.6. Young population
In this section, we apply the previous analysis to the young population. To do so, we only chose stars brighter than an absolute magnitude M_{G} = −2 (see the luminosity function in Fig. 17) and repeated all the steps as described in Sect. 4.3. Then we produced density maps and analysed the scale length and the warp of this population using methods from Sects. 5.4 and 5.5.
The exponential fits of the density for the young population yield h_{r} = 2.5 ± 0.22 kpc for ϕ = 30°, h_{r} = 1.92 ± 0.15 kpc for ϕ = 0°, and h_{r} = 2.04 ± 0.15 kpc for ϕ = 30°, which is similar to the whole population. This results in h_{r} = 2.09 ± 0.09 kpc on average. The variation with azimuth is still insignificant, as in the case of the entire population. Figure 18 shows that the variation of density with azimuth is also negligible in the case of the young population.
For the warp, as previously, we removed the azimuths 150° < ϕ < 240° from the analysis. The fit of Eq. (11) to the young population yields
We also repeated the analysis with the approach using Eq. (13). Figure 19 presents the Galactic warp of the young stellar population for various Galactocentric distances, and Fig. 20 shows the amplitudes of the fits of the Galactic warp and the azimuth of the maximum and minimum. In this case, the warp of the young stellar population is stronger than the case considering all stars in our dataset. The azimuth of the maximum of the northern warp is an increasing function of the Galactocentric distance (50° < ϕ < 54°), and the azimuth of the minimum of the warp is in 265° < ϕ < 315°. The highest amplitude of the northern and the southern warp is z_{w} = 0.57 kpc and z_{w} = −0.5 kpc, respectively. For the line of nodes, we find ϕ_{0} = −6.56°, which agrees with the whole population.
Fig. 19. Dataset containing a young population of stars. The average elevation of the plane as a function of the Galactic azimuth for various values of the Galactocentric distance. Red markers represent values of binned data, and the blue line represents a fit to the data. (a) R = 13.25 kpc; (b) R = 16.25 kpc; (c) R = 19.25 kpc. 
Fig. 20. Minimum and maximum of the average elevation of the plane as a function of the Galactocentric distance. The warp fit is based on Eq. (13). The colours code the azimuth of the minimum and the maximum of the warp fit. The dataset containing a young population of stars is considered, and the error bars represent the uncertainty on the distance in the Lucy method. 
Chen et al. (2019) used Cepheids from the WISE survey and a number of optical surveys to measure the warp, and Skowron et al. (2019a) used Cepheids from the OGLE catalogue supplemented by other surveys. Chen et al. (2019) obtained a rather symmetric warp with an amplitude of about 1.5 kpc in R = 20 kpc. Skowron et al. (2019a) obtained a similar result with an amplitude 0.74 kpc in R = 15 kpc. These values are much higher than our findings, which is probably due to differences in the population: our young population is older than the Cepheids. In Fig. 21 we plot the variation of the line of nodes with radius for the whole and the young population, compared with other works. We use two different methods to plot the line of nodes for our work. First, we plot the angle ϕ_{w} for Eq. (12). Another method is to use the Eq. (11) to find the value of the angle ϕ when z_{w} = 0. We would expect that our young population lies between the total population and the young Cepheids. This is true only for R > 12 kpc. At shorter distances, the warp is not very strong and is more difficult to detect, therefore the error bars are larger in this area. Moreover, the error bars of the young populations are very large because of the lower number of stars in the sample combined with possible problems in determining extinction. For these reasons, the value of the line of nodes for R < 12 kpc is rather unreliable.
Fig. 21. Comparison of line of nodes for our model (based on Eq. (11)) with other works. We use two different methods to plot the line of nodes for our work. First, we plot the angle ϕ_{w} for Eq. (12). Another method is to use the Eq. (11) to find the value of the angle ϕ when z_{w} = 0. 
6. Conclusions
We produced density maps from Gaia DR2 data and analysed them to study the Galactic warp. The density maps directly show a northern warp in the azimuths 60° < ϕ < 90° and a southern warp in the azimuths 270° < ϕ < 300°. Our maps reach a Galactocentric radius of 20 kpc, and we note that up to this distance, the density decreases exponentially and we do not observe a cutoff. Another feature in the density maps is a Galactic flare, that is, an increase in scale height towards the outer Galaxy. The analysis of the flare will be given in a forthcoming paper (Nagy et al., in prep.). We used the maps to calculate the scale length, where we find h_{r} = 2.29 ± 0.08 kpc, with a small dependence of h_{r} from the Galactic azimuth. The lowest value of h_{r} that we found is 1.88 ± 0.12 kpc for ϕ ≈ ±−30° and the highest value is 2.78 ± 0.13 kpc for the Sunanticentre direction and 2.29 ± 0.21 ϕ ≈ ±30°.
From our maps, we calculated the average elevation of the plane and fitted it with different warp models. We fitted the northern and southern warp separately with a simple sinusoidal model, and we found a small asymmetry: the northern warp reaches an amplitude 0.317 kpc for the azimuth 60° < ϕ < 90° and the southern warp reaches −0.287 kpc for the azimuth 240° < ϕ < 270°, both at R = [19.5, 20.0] kpc. Then we fitted the warp with a model combining two sinusoids to detect the asymmetry without assuming its existence, and we found values of amplitude ∼0.5 for the northern and ∼ − 0.4 for the southern warp both at R = [19.5, 20.0] kpc, revealing the asymmetry found with the previous approach. The azimuths of the warp maximum and minimum for this model are 52° < ϕ < 56° and 312° < ϕ < 324°, respectively. In terms of Galactocentric radius, we find that warp starts to manifest itself from about 12 kpc and extends at least up to 20 kpc. We repeated this analysis on the young population, where we find that it follows the result for the whole population, but reaches a higher amplitude of warp and similar values of scale height. The comparison of our amplitude of the warp with other works showed that we obtain a significantly lower amplitude than an analysis carried out with very young stars such as Cepheids. This supports the formation of the warp through accretion onto the disc (LópezCorredoira et al. 2002b).
A future analysis of the next Gaia data release combined with the deconvolution method based on Lucy’s method of inversion, as described in Sect. 4.1, will allow us to explore distances larger than 20 kpc. The future data release will provide a much deeper magnitude limit and much lower parallax errors, which will allow us to extend the range of Galactocentric distances and study the morphology of the disc and of the stellar halo at very large distances.
Acknowledgments
We thank the anonymous referee for helpful comments, which improved this paper, and Astrid Peter (language editor of A&A) for proofreading of the text. ZC and MLC were supported by the Grant PGC2018102249B100 of the Spanish Ministry of Economy and Competitiveness (MINECO). RN was supported by the Scientific Grant Agency VEGA No. 1/0911/17. This work made use of the IAC Supercomputing facility HTCondor (http://research.cs.wisc.edu/htcondor/), partly financed by the Ministry of Economy and Competitiveness with FEDER funds, code IACA133E2493. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. The reduced catalogue of Gaia with m_{G} < 19 was produced by Pedro Alonso Palicio.
References
 Amôres, E. B., Robin, A. C., & Reylé, C. 2017, A&A, 602, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Anders, F., Khalatyan, A., Chiappini, C., et al. 2019, A&A, 628, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Arenou, F., Luri, X., Babusiaux, C., et al. 2018, A&A, 616, A17 [Google Scholar]
 Bahcall, J. N. 1986, ARA&A, 24, 577 [NASA ADS] [CrossRef] [Google Scholar]
 Bahcall, J. N., & Soneira, R. M. 1980, ApJS, 44, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Balázs, L. G. 1995, Inverse Prob., 11, 731 [NASA ADS] [CrossRef] [Google Scholar]
 Battaner, E., Florido, E., & SanchezSaavedra, M. L. 1990, A&A, 236, 1 [NASA ADS] [Google Scholar]
 Bovy, J., Rix, H.W., Liu, C., et al. 2012, ApJ, 753, 148 [NASA ADS] [CrossRef] [Google Scholar]
 Bovy, J., Rix, H.W., Green, G. M., Schlafly, E. F., & Finkbeiner, D. P. 2016, ApJ, 818, 130 [NASA ADS] [CrossRef] [Google Scholar]
 Brand, J., & Wouterloot, J. G. A. 2007, A&A, 464, 909 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Carney, B. W., & Seitzer, P. 1993, AJ, 105, 2127 [NASA ADS] [CrossRef] [Google Scholar]
 Carraro, G. 2014, in Setting the Scene for Gaia and LAMOST, eds. S. Feltzing, G. Zhao, N. A. Walton, & P. Whitelock, IAU Symp., 298, 7 [NASA ADS] [Google Scholar]
 Chandresekhar, S., & Münch, G. 1951, ApJ, 113, 150 [NASA ADS] [CrossRef] [Google Scholar]
 Chang, C.K., Ko, C.M., & Peng, T.H. 2011, ApJ, 740, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, B., Stoughton, C., Smith, J. A., et al. 2001, ApJ, 553, 184 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, X., Wang, S., Deng, L., et al. 2019, Nat. Astron., 3, 320 [NASA ADS] [CrossRef] [Google Scholar]
 Craig, I. J. D., & Brown, J. C. 1986, Inverse Problems in Astronomy. A Guide to Inversion Strategies for Remotely Sensed Data (Bristol, UK: Adam Hilger) [Google Scholar]
 Debattista, V. P., & Sellwood, J. A. 1999, ApJ, 513, L107 [NASA ADS] [CrossRef] [Google Scholar]
 Dehnen, W. 1998, AJ, 115, 2384 [NASA ADS] [CrossRef] [Google Scholar]
 Drimmel, R., CabreraLavers, A., & LópezCorredoira, M. 2003, A&A, 409, 205 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Eaton, N., Adams, D. J., & Giles, A. B. 1984, MNRAS, 208, 241 [NASA ADS] [CrossRef] [Google Scholar]
 Gaia Collaboration (Prusti, T., et al.) 2016, A&A, 595, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Green, G. M., Schlafly, E. F., Finkbeiner, D. P., et al. 2015, ApJ, 810, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Green, G. M., Schlafly, E. F., Finkbeiner, D., et al. 2018, MNRAS, 478, 651 [NASA ADS] [CrossRef] [Google Scholar]
 Hendy, Y. H. M. 2018, NRIAG J. Astron. Geophys., 7, 180 [NASA ADS] [CrossRef] [Google Scholar]
 Jurić, M., Ivezić, Ž., Brooks, A., et al. 2008, ApJ, 673, 864 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Karim, M. T., & Mamajek, E. E. 2017, MNRAS, 465, 472 [NASA ADS] [CrossRef] [Google Scholar]
 Kerr, F. J. 1957, AJ, 62, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Kim, J. H., Peirani, S., Kim, S., et al. 2014, ApJ, 789, 90 [NASA ADS] [CrossRef] [Google Scholar]
 Levine, E. S., Blitz, L., & Heiles, C. 2006, ApJ, 643, 881 [NASA ADS] [CrossRef] [Google Scholar]
 Li, C., Zhao, G., Jia, Y., et al. 2019, ApJ, 871, 208 [NASA ADS] [CrossRef] [Google Scholar]
 Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, A&A, 616, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Liu, C., Xu, Y., Wan, J.C., et al. 2017, Res. Astron. Astrophys., 17, 096 [NASA ADS] [CrossRef] [Google Scholar]
 LópezCorredoira, M., & Molgó, J. 2014, A&A, 567, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 LópezCorredoira, M., & SylosLabini, F. 2019, A&A, 621, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 LópezCorredoira, M., Hammersley, P. L., Garzón, F., Simonneau, E., & Mahoney, T. J. 2000, MNRAS, 313, 392 [NASA ADS] [CrossRef] [Google Scholar]
 LópezCorredoira, M., CabreraLavers, A., Garzón, F., & Hammersley, P. L. 2002a, A&A, 394, 883 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 LópezCorredoira, M., BetancortRijo, J., & Beckman, J. E. 2002b, A&A, 386, 169 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 LópezCorredoira, M., CabreraLavers, A., Gerhard, O. E., & Garzón, F. 2004, A&A, 421, 953 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 LópezCorredoira, M., Abedi, H., Garzón, F., & Figueras, F. 2014, A&A, 572, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 LópezCorredoira, M., Garzon, F., Wang, H. F., et al. 2020, A&A, 634, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lucy, L. B. 1974, ApJ, 79, 745 [Google Scholar]
 Lucy, L. B. 1994, A&A, 289, 983 [NASA ADS] [Google Scholar]
 Luri, X., Brown, A. G. A., Sarro, L. M., et al. 2018, A&A, 616, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Majewski, S. R. 1993, ARA&A, 31, 575 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Marshall, D. J., Robin, A. C., Reylé, C., Schultheis, M., & Picaud, S. 2006, A&A, 453, 635 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Minniti, D., Saito, R. K., AlonsoGarcía, J., Lucas, P. W., & Hempel, M. 2011, ApJ, 733, L43 [NASA ADS] [CrossRef] [Google Scholar]
 Momany, Y., Zaggia, S., Gilmore, G., et al. 2006, A&A, 451, 515 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Oort, J. H., Kerr, F. J., & Westerhout, G. 1958, MNRAS, 118, 379 [NASA ADS] [CrossRef] [Google Scholar]
 Paul, E. R. 1993, The Milky Way Galaxy and Statistical Cosmology, 1890 [Google Scholar]
 Poggio, E., Drimmel, R., Lattanzi, M. G., et al. 2018, MNRAS, 481, L21 [NASA ADS] [CrossRef] [Google Scholar]
 Queiroz, A. B. A., Anders, F., Santiago, B. X., et al. 2018, MNRAS, 476, 2556 [NASA ADS] [CrossRef] [Google Scholar]
 Reylé, C., Marshall, D. J., Robin, A. C., & Schultheis, M. 2009, A&A, 495, 819 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Riess, A. G., Casertano, S., Yuan, W., et al. 2018, ApJ, 861, 126 [Google Scholar]
 Robin, A. C., Creze, M., & Mohan, V. 1992, ApJ, 400, L25 [NASA ADS] [CrossRef] [Google Scholar]
 Robin, A. C., Reylé, C., Derrière, S., & Picaud, S. 2003, A&A, 409, 523 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Robin, A. C., Marshall, D. J., Schultheis, M., & Reylé, C. 2012, A&A, 538, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rybizki, J., Demleitner,, M., Fouesneau, M., et al. 2018, PASP, 130 [Google Scholar]
 Sale, S. E., Drew, J. E., Knigge, C., et al. 2010, MNRAS, 402, 713 [NASA ADS] [CrossRef] [Google Scholar]
 Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525 [NASA ADS] [CrossRef] [Google Scholar]
 Schönrich, R., & Dehnen, W. 2018, MNRAS, 478, 3809 [NASA ADS] [CrossRef] [Google Scholar]
 Skowron, D. M., Skowron, J., Mróz, P., et al. 2019a, Science, 365, 478 [NASA ADS] [CrossRef] [Google Scholar]
 Skowron, D. M., Skowron, J., Mróz, P., et al. 2019b, Acta Astron., 69, 305 [NASA ADS] [Google Scholar]
 Stassun, K. G., & Torres, G. 2018, ApJ, 862, 61 [NASA ADS] [CrossRef] [Google Scholar]
 Turchin, V. F., Kozlov, V. P., & Malkevich, M. S. 1971, Sov. Phys. Usp., 13, 681 [Google Scholar]
 Voskes, T., & Butler Burton, W. 2006, ArXiv eprints [arXiv:astroph/0601653] [Google Scholar]
 Vozoff, K., & Jupp, D. L. B. 1975, Geophys. J. R. Astron. Soc., 42, 977 [CrossRef] [Google Scholar]
 Wainscoat, R. J., Cohen, M., Volk, K., Walker, H. J., & Schwartz, D. E. 1992, ApJS, 83, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Yusifov, I. 2004, in The Magnetized Interstellar Medium, eds. B. Uyaniker, W. Reich, & R. Wielebinski, 165 [Google Scholar]
 Zinn, J. C., Pinsonneault, M. H., Huber, D., & Stello, D. 2019, ApJ, 878, 136 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Lucy’s method for the inversion of Fredholm integral equations of the first kind
The inversion of Fredholm integral equations of the first kind such as Eq. (2) is illconditioned. Typical analytical methods for solving these equations (Balázs 1995) cannot achieve a good solution because the kernel is sensitive to the noise of the star counts (Craig & Brown 1986, chapter 5). Because the functions in these equations have a stochastic rather than analytical interpretation, it is to be expected that statistical inversion algorithms are more robust (Turchin et al. 1971; Vozoff & Jupp 1975; Balázs 1995). These statistical methods include the iterative method of Lucy’s algorithm (Lucy 1974; Turchin et al. 1971; Balázs 1995; LópezCorredoira et al. 2000), which is appropriate here. Its key feature is the interpretation of the kernel as a conditioned probability and the application of Bayes’ theorem.
In Eq. (2), N(π) is the unknown function, and the kernel is G(x), whose difference x is conditioned to the parallax π′. The inversion is carried out as
The iteration converges when ∀π, that is, when N_{n}(π) ≈ N(π) ∀π. The first iterations produce a result that is close to the final answer, with the subsequent iterations giving only small corrections. In our calculation, we set as initial function of the iteration and we carry out a number of iterations until the Pearson χ^{2} test
reaches the minimum value. Further iterations would enter within the noise.
This algorithm has a number of beneficial properties (Lucy 1974, 1994): all the functions are defined as being positive, the likelihood increases with the number of iterations, the method is insensitive to highfrequency noise in , and so on. We note, however, that precisely because this method only works when N are positive functions, it does not work with negative ones.
All Tables
All Figures
Fig. 1. Luminosity function. 

In the text 
Fig. 2. Interpolation of the luminosity function with a spline compared with the luminosity function of Bahcall & Soneira (1980). These two functions are not directly comparable because Bahcall & Soneira (1980) used a slightly different filter in the visible, but it shows that our luminosity function is reasonable. 

In the text 
Fig. 3. Monte Carlo simulation of deconvolution. We recover random distributions, convolved with a Gaussian. (a) Gamma distribution; (b) exponential distribution; (c) geometric distribution; (d) logarithmic distribution. 

In the text 
Fig. 4. Monte Carlo simulation of deconvolution. In all cases we recover a gamma distribution convolved with a Gaussian, varying the average error of parallax. (a) σ_{π} = 0.05 mas; (b) σ_{π} = 0.1 mas; (c) σ_{π} = 0.25 mas; (d) σ_{π} = 0.4 mas. 

In the text 
Fig. 5. Density maps for various azimuths between 0° and 360°. (a) 0° < ϕ < 30°; (b) 30° < ϕ < 60°; (c) 60° < ϕ < 90°; (d) 270° < ϕ < 300°; (e) 300° < ϕ < 330°; (f) 330° < ϕ < 360°. 

In the text 
Fig. 6. Density maps. 

In the text 
Fig. 7. Comparison of densities for various lines of sight. Orange and green curves represent the density including the zeropoint correction of the parallax, and the blue curve shows the density without this correction. 

In the text 
Fig. 8. Relative error of the density calculated including the standard deviation of the extinction. 

In the text 
Fig. 9. Comparison of luminosity functions for the thin and thick disc. 

In the text 
Fig. 10. Densities for different lines of sight. Blue curves are densities calculated with the thindisc luminosity function. Orange curves are densities calculated with the thickdisc luminosity function. 

In the text 
Fig. 11. Dependence of the density on azimuth near the centreSunanticentre direction for various values of Galactocentric distance. The data points are obtained as weighted mean in bins of size 1 kpc in R and 0.4 kpc in z. Only bins with a number of points N ≥ 50 points are plotted. 

In the text 
Fig. 12. Minimum and maximum average elevation of the plane as a function of radius. The warp fit is based on Eq. (11), and the error bars represent the uncertainty in the distance in the Lucy method. 

In the text 
Fig. 13. Average elevation of the plane as a function of the Galactic azimuth for various values of the Galactocentric distance. Red markers represent values of binned data, and the blue line represents a fit to the data. (a) R = 13.25 kpc; (b) R = 16.25 kpc; (c) R = 19.25 kpc. 

In the text 
Fig. 14. Minimum and maximum of the average elevation of the plane as a function of Galactocentric distance. The warp fit is based on Eq. (13). The colours code the azimuth of the minimum and maximum of the warp fit, and the error bars represent the uncertainty on the distance in the Lucy method. 

In the text 
Fig. 15. Changes of the amplitudes of the Galactic warp fit according to Eqs. (13) and (14). 

In the text 
Fig. 16. Comparison of maximum amplitudes of our model (based on Eq. (13)) with other works. 

In the text 
Fig. 17. Luminosity function used in the Eq. (6) for the analysis of the young population. 

In the text 
Fig. 18. Same as Fig. 11, but the young population alone is considered. 

In the text 
Fig. 19. Dataset containing a young population of stars. The average elevation of the plane as a function of the Galactic azimuth for various values of the Galactocentric distance. Red markers represent values of binned data, and the blue line represents a fit to the data. (a) R = 13.25 kpc; (b) R = 16.25 kpc; (c) R = 19.25 kpc. 

In the text 
Fig. 20. Minimum and maximum of the average elevation of the plane as a function of the Galactocentric distance. The warp fit is based on Eq. (13). The colours code the azimuth of the minimum and the maximum of the warp fit. The dataset containing a young population of stars is considered, and the error bars represent the uncertainty on the distance in the Lucy method. 

In the text 
Fig. 21. Comparison of line of nodes for our model (based on Eq. (11)) with other works. We use two different methods to plot the line of nodes for our work. First, we plot the angle ϕ_{w} for Eq. (12). Another method is to use the Eq. (11) to find the value of the angle ϕ when z_{w} = 0. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.