Structure of the outer Galactic disc with Gaia-DR2
AAstronomy & Astrophysics manuscript no. gaia_maps_final c (cid:13)
ESO 2020April 8, 2020
Structure of the outer Galactic disc with
Gaia
DR2
Ž. Chrobáková , , R. Nagy , M. López-Corredoira , Instituto de Astrofísica de Canarias, E-38205 La Laguna, Tenerife, Spain Departamento de Astrofísica, Universidad de La Laguna, E-38206 La Laguna, Tenerife, Spain Faculty of Mathematics, Physics, and Informatics, Comenius University, Mlynská dolina, 842 48 Bratislava, SlovakiaReceived xxxx; accepted xxxx
ABSTRACT
Context.
The structure of outer disc of our Galaxy is still not well described, and many features need to be betterunderstood. The second Gaia data release (DR2) provides data in unprecedented quality that can be analysed to shedsome 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=20kpc with a deconvolution technique for the parallax errors. Then we analyse the density in order to study the structureof the outer Galactic disc, mainly the warp.
Methods.
In order to carry out the deconvolution, we used the Lucy inversion technique for recovering the correctedstar counts. We also used the Gaia luminosity function of stars with M G < to extract the stellar density from thestar counts. Results.
The stellar density maps can be fitted by an exponential disc in the radial direction h r = 2 . ± . kpc, with a weak dependence on the azimuth, extended up to 20 kpc without any cut-off. The flare and warpare clearly visible. The best fit of a symmetrical S-shaped warp gives z w ≈ z (cid:12) + (37 ± . stat. ) − . syst. )) pc · ( R/R (cid:12) ) . ± . stat. )+0 . syst. ) sin ( φ + 9 . ◦ ± . ◦ ( stat. ) + 4 . ◦ ( syst. )) for the whole population. When we analysethe northern and southern warps separately, we obtain an asymmetry of an ∼ larger amplitude in the north. Thisresult may be influenced by extinction because the Gaia G band is quite prone to extinction biases. However, we testedthe 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 ofthe map, which may influence the final result.The analysis was also carried out for very luminous stars alone ( M G < − ), which on average represents a youngerpopulation. We obtain similar scale-length values, while the maximum amplitude of the warp is − larger thanwith the whole population. The north-south asymmetry is maintained. Key words.
Galaxy:disc – Galaxy: structure
1. Introduction
Studying the Galactic structure is crucial for our under-standing of the Milky Way. Star counts are widely used forthis purpose (Paul 1993), and the importance of this toolhas increased in the past decades with the appearance ofwide-area surveys (Bahcall 1986; Majewski 1993), whichmade it possible to obtain reliable measurements of theGalactic thin- and thick-disc and halo (Chen et al. 2001; Ju-rić et al. 2008; Bovy et al. 2012; Robin et al. 2012). It is com-mon to simplify the Galactic disc as an exponential or hy-perbolic secant form, but there are many asymmetries suchas 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 Wayusing the LAMOST (The Large Sky Area Multi-Object Fi-bre Spectroscopic Telescope) RGB (red-giant branch) stars;Skowron et al. (2019a), who constructed a map of the MilkyWay from classical Cepheids; or Anders et al. (2019), whoused the second Gaia data release (DR2).The warp was first detected in the Galactic gaseousdisc in 21 cm HI observations (Kerr 1957; Oort et al.1958). Since then, the warp has also been discovered in thestellar disc (Carney & Seitzer 1993; López-Corredoira et al. 2002b; Reylé et al. 2009; Amôres et al. 2017; Chen et al.2019), and the kinematics of the warp has been studied aswell (Dehnen 1998; Drimmel et al. 2003; López-Corredoiraet al. 2014; Schönrich & Dehnen 2018).Vertical kinematics in particular can reveal much aboutthe mechanism behind the formation of warp. Poggio et al.(2018) found a gradient of − km/s in the verticalvelocities of upper main-sequence stars and giants locatedfrom 8 to 14 kpc in Galactic radius using Gaia DR2data, revealing the kinematic signature of the warp. Theirfindings suggest that the warp is principally a gravitationalphenomenon. Skowron et al. (2019b) also found a stronggradient in vertical velocities using classical Cepheidssupplemented by the OGLE (Optical Gravitational Lens-ing Experiment) survey. Lopez-Corredoira et al. (2020)investigated the dynamical effects produced by differentmechanisms that can explain the radial and vertical com-ponents of extended kinematic maps of López-Corredoira& Sylos-Labini (2019), who used Lucy’s deconvolutionmethod (see Sect. 4.1) to produce kinematical maps upto a Galactocentric radius of 20 kpc. Lopez-Corredoiraet al. (2020) found that vertical motions might be dom- Article number, page 1 of 19 a r X i v : . [ a s t r o - ph . GA ] A p r &A proofs: manuscript no. gaia_maps_final inated by external perturbations or mergers, althoughwith a minor component due to a warp whose amplitudeis evolving with time. However, the kinematic signatureof the warp is not enough to explain the observed velocities.To date, the shape of the warp has been constrained onlyroughly, and the kinematical information is not satisfyingenough to reach consensus about the mechanism causingthe warp. Theories include accretion of intergalactic mat-ter onto the disc (López-Corredoira et al. 2002a), interac-tion with other satellites (Kim et al. 2014), the intergalacticmagnetic field (Battaner et al. 1990), a misaligned rotatinghalo (Debattista & Sellwood 1999), and others.We now have a new opportunity to improve our knowl-edge about the Milky Way significantly through the Gaiamission of the European Space Agency (Gaia Collaborationet al. 2016). Gaia data provide unprecedented positionaland radial velocity measurements and an accurate distancedetermination, although the error of the parallax measure-ment increases with distance from us. It brings us the mostaccurate data about the Galaxy so far, ideal to advance inall branches of Galactic astrophysics and study our Galaxyin greater detail than ever before. Gaia DR2 has been usedby Anders et al. (2019), who provided photo-astrometricdistances, extinctions, and astrophysical parameters up tomagnitude G=18, making use of the Bayesian parameterestimation code StarHorse . After introducing the obser-vational data and a number of priors, their code finds theBayesian stellar parameters, distances, and extinctions. Theauthors also present density maps, which we compare withour results in Section 4.3. Gaia data have also been usedto study the structure of outer Galactic disc, especially thewarp and the flare. The first Gaia data release brought someevidence of the warp (Schönrich & Dehnen 2018), but themore extensive second data release provides a better op-portunity to study the warp attributes. Poggio et al. (2018)combined Gaia DR2 astrometry with 2MASS (Two MicronAll-Sky Survey) photometry and revealed the kinematic sig-nature 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 onlyOB stars of the Gaia DR2. In this work, we make use ofGaia DR2 data as described in Section 2 and use star countsto obtain the stellar density by applying Lucy’s inversiontechnique. Then we analyse the density maps to determinethe warp.The paper is structured as follows: in Section 2 we de-scribe the Gaia data and extinction maps that we used, inSection 3 we present the luminosity function used in ourcalculations, in Section 4 we explain the methods for ob-taining our density maps, and in Section 5 we discuss theresults. In Section 5.4 we present the exponential fits of thedensity, in Section 5.5 we study the warp, and in Section 5.6we repeat the previous analysis of the young population.
2. Data selection
We used data of the second Gaia data release (Gaia Col-laboration et al. 2018) here, which were collected duringfirst 22 months of observation. We are interested in starswith known five-parameter astrometric solution: more than1.3 billion sources. G magnitudes, collected by astromet-ric instrument in the white-light G-band of Gaia (330–1050nm) are known for all sources, with precisions varying fromaround 1 millimag at the bright (G<13) end to around 20 millimag at G=20. For the details on the astrometric dataprocessing and validation of these results, see Lindegrenet al. (2018). We chose stars with apparent magnitude up toG=19, where the catalogue is complete up to 90% (Arenouet al. 2018). We chose data with a parallax in the interval[0,2] mas.In our analysis, we did not consider any zero-point biasin the parallaxes of Gaia DR2, as found by some authors(Lindegren et al. 2018; Arenou et al. 2018; Stassun & Torres2018; Zinn et al. 2019), except in Sect. 4.3, where we repeatour main calculation including a non-zero value of the zero-point to prove that this effect is negligible in our results.
Extinction maps
We used two different extinction maps. For the luminosityfunction (Sect. 3), we used the extinction map of Greenet al. (2018) through its Python package dustmaps , choos-ing the
Bayestar17 version. This map covers 75% of thesky (declinations of δ (cid:38) − ◦ ) and provides reddening insimilar units as Schlegel et al. (1998, SFD).To calculate the density (Sect. 4), we need to cover thewhole sky, therefore we used the three-dimensional less ac-curate but full-sky extinction map of Bovy et al. (2016)through its Python package mwdust . This map combinesthe results of Marshall et al. (2006), Green et al. (2015),and Drimmel et al. (2003) and provides reddening as de-fined in Schlegel et al. (1998).In order to convert the interstellar reddening of these mapsinto E ( B − V ) , we used coefficients (Hendy 2018; Rybizkiet al. 2018) A G /A v = 0 . ,R V = A v /E ( B − V ) = 3 . . (1)
3. Luminosity function
To construct the luminosity function, we chose all stars withheliocentric distance d < . kpc (distances determined as /π , where π is the parallax). We did not find many brightstars ( M G < − ) in this area, therefore we also chose aspecific region with Galactic height | z | < kpc and Galac-tocentric distance R < kpc, in which we only selectedstars with absolute magnitude M G < − . We normalisedthe counts of stars with high magnitude and then joinedthese two parts to create the luminosity function.In the range of distance that we used for the luminos-ity function, the star counts are complete for the absolutemagnitude that we are calculating, except perhaps for thepossible loss of the brightest stars through saturation at M G < − . Moreover, the error in the parallax for thesestars is negligible, so that the calculation of the absolutemagnitude from the apparent magnitude is quite accurate.We did not take into account the variations of the lumi-nosity function throughout the Galactic disc. We assumedthat 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 = [ − , because the values outside this interval areunreliable, and we used the extrapolation of the spline func- Article number, page 2 of 19. Chrobáková, R. Nagy, M. López-Corredoira: Structure of the outer Galactic disc with
Gaia
DR2 - 1 0 - 8 - 6 - 4 - 2 0 2 4 6 8 10M1e-81e-61e-40.01 N Fig. 1: Luminosity function.
15 10 5 0 5 10 15M10 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 N interpolationBahcall & Soneira (1980)data Fig. 2: Interpolation of the luminosity function with a splinecompared with the luminosity function of Bahcall & Soneira(1980). These two functions are not directly comparable be-cause Bahcall & Soneira (1980) used a slightly different fil-ter in the visible, but it shows that our luminosity functionis reasonable.tion to lower magnitudes. The values of the luminosity func-tion are listed in Table 1.
4. Density maps
To calculate the stellar density, we need to measure starcounts as a function of distance. However, the error of par-allax increases with distance from us, which means thatour analysis would be correct only within roughly kpcfrom the Sun. To be able to reach higher distances, we cor-rected for this effect using the method developed by López-Corredoira & Sylos-Labini (2019), who used Lucy’s decon-volution method (Lucy 1974; see Appendix A) to obtain anaccurate distance measurement up to R = 20 kpc. They M G N-10 . · − -9 . · − -8 . · − -7 . · − -6 . · − -5 . · − -4 . · − -3 . · − -2 . · − -1 . · − . · − . · − . · − . · − . · − . · − . · − N ( π ) as a convolution of the real number N ( π ) of stars with aGaussian function N ( π ) = (cid:90) ∞ d π (cid:48) N ( π (cid:48) ) G π (cid:48) ( π − π (cid:48) ) , (2)where G π ( x ) = 1 √ πσ π e − x σ π . (3)For the error σ π we averaged errors of every bin, which wecalculated from values given by Gaia DR2.We only used the parallax between [0,2] mas. For the upperlimit the relative error of parallax is very small and doesnot produce any bias. For the lower limit, the truncationavoiding the negative parallaxes affects the distribution ofparallaxes and statistical properties (average, median, etc.)(Luri et al. 2018, Section 3.3). However, in our method wedo not calculate the average distance from the average par-allax. We used Lucy’s method, which iterates the countsof the stars with positive parallaxes until we obtained thefinal solution. This does not mean that we truncated thestar counts with negative parallaxes. We used only the starswith positive parallaxes as is required by our method, ex-plained in the Appendix A. N ( π ) for negative values of π can also be calculated and fitted, but they are not usedin our calculation. In other words, we did not assume thatthe number of the stars with negative parallaxes is zero, wesimply did not use this information because it is not neces-sary. The fact that this method does not produce any biasis tested in section 4.2. Article number, page 3 of 19 &A proofs: manuscript no. gaia_maps_final
In order to test the reliability of the inversion method, weperformed Monte Carlo simulations to determine whetherwe can recover the original function after deconvolution. Wecreated datasets with randomly distributed particles. Thenwe convolved this distribution with a Gaussian. We appliedLucy’s deconvolution method to the dataset to determinewhether we can recreate the original distribution. The re-sults are shown in Fig. 3. We conclude that regardless of theoriginal distribution, we can accurately recover the originaldata up to kpc or more, which is satisfying to study theMilky Way. We also studied the dependence of the methodon the parallax error. We used various values of the aver-age parallax error in Eq. 3 from the interval [0.05,0.4] mas,which are the most common values for the average parallaxerror in our data. In Fig. 4 we plot the result, which showsthat even though the precision of the method depends onthe parallax error, we obtain a satisfying result up to 20kpc even in the worst case with the highest parallax error. We divided the data into bins of Galactic longitude (cid:96) ,Galactic latitude b, and apparent magnitude m . For thevalues of b we made bins of length ◦ and corresponding (cid:96) in bins of ◦ /cos ( b ) . We divided each of the lines of sightin magnitude, binned with size ∆ m = 1 . between G=12and G=19. We obtained 29 206 different areas in which wecalculated the density independently.We made use of the fundamental equation of stellarstatistics, where the number of stars N ( m ) of apparentmagnitude m is expressed per unit solid angle and per unitmagnitude interval (Chandresekhar & Münch 1951), N ( m ) = (cid:90) ∞ ρ ( r )Φ( M ) r d r , (4)where we substitute r ( m ) = (1 /π ) = 10 ( m − M +5 − A G (1 /π )) / , (5)which yields for the density ρ (1 /π ) = N ( π ) π ∆ πω (cid:82) M G,low lim +1 M G,low lim d M G Φ( M G ) , (6) M G,low lim = m G,low lim − log (1 /π ) − − A G (1 /π ) , (7)where ω is the covered angular surface (
10 degrees inour case), ∆ π is the parallax interval (0.01 mas in our case),which must be added in the equation because we did notuse the unit parallax, Φ( M G ) is the luminosity function inthe G filter, m G,low lim is the limiting maximum apparentmagnitude, and A G ( r ) is the extinction, as a function ofdistance.After this, we calculated the weighted mean densityfor all seven ranges of magnitude in each line of sight. Then we transformed this into cylindrical coordinates andmade bins of Galactocentric radius R of length 0.5 kpc,in Galactic height z of 0.1 kpc and in azimuth of ◦ . Wedefine the azimuthal angle φ to be measured from thecentre-Sun-anticentre direction towards the Galactic rota-tion, going from ◦ to ◦ . We interpolated the missingbins with NearestNDInterpolator from the python
SciPy package, which uses nearest-neighbour interpolation in Ndimensions. We plot the resulting density maps in Fig. 5-6.In Fig. 5 we plot the density in cylindrical coordinates as afunction of Galactic radius R for different azimuths. We donot plot the results for azimuths ◦ < φ < ◦ because inthis area the extinction is significant and we observe starsfarther than the Galactic centre for which the errors aretoo large, therefore we cannot see any structure in density.However, we can see even by eye that a northern warp ispresent in the azimuths ◦ < φ < ◦ and a southernwarp in the azimuths ◦ < φ < ◦ . Another structurethat can be seen from the plots is the flaring of the disc.We analyse these structures below. In Fig. 6. (a)-(c) weplot the density map in Cartesian coordinates, and inFig 6. (d) we plot the density in cylindrical coordinates,integrated through all ranges of azimuths, except for theareas that were excluded from the analysis. The Cartesiancoordinates are defined such that X (cid:12) = 8 . kpc. Inthese plots we note a flat disc with some fluctuations indensity, but no apparent features. However, some slightsoverdensities both above and below the Galactic plane arevisible. The features above the plane are present only inFig. 6. (b)-(c), but not in Fig. 6. (d), which suggests that itmight be a contamination. The feature below the Galacticplane is present in all the three plots. As the direction ofthese overdensities is towards the Magellanic Clouds, itmight be an effect of the Milky Way pulling stars out ofMagellanic Clouds, as suggested by Anders et al. (2019).Another possible explanation for these overdensities is thefinger of God artefact, which is caused by the foregrounddust clouds and causes elongated overdensities that pointto the Sun. This artefact has previously been seen in Gaiadata, as shown in the Gaia DR2 documentation . So far, we did not consider any zero-point bias in paral-laxes. Lindegren et al. (2018) found a global mean offset of − . mas, meaning that Gaia DR2 parallaxes are lowerthan the true value. We repeated our calculations with thiscorrection and present the results in Fig. 7, where we chosesome of the lines of sight to show the comparison. We findthat these results are very similar to our original results,and this correction brings a negligible effect. We also trieda value of − . mas, found by Riess et al. (2018). In Fig.7 we show that the difference between the different zero-point values is very small, therefore we only use the valueof -0.029 mas in the further calculations.For the analysis of the warp in Sections 5.5 and 5.6., werepeated the analysis of Section 4 with the value of paral-lax corrected for the zero-point. We find that this brings a https://gea.esac.esa.int/archive/documentation/GDR2/Data_analysis/chap_cu8par/sec_cu8par_validation/ssec_cu8par_validation_additional-validation.html Article number, page 4 of 19. Chrobáková, R. Nagy, M. López-Corredoira: Structure of the outer Galactic disc with
Gaia
DR2 -1 d [kpc] -5 -4 -3 -2 -1 N convolutionoriginal distributiondeconvolution (a) Gamma distribution -1 d [kpc] -4 -3 -2 -1 N convolutionoriginal distributiondeconvolution (b) Exponential distribution -1 d [kpc] -1 N convolutionoriginal distributiondeconvolution (c) Geometric distribution -1 d [kpc] -2 -1 N convolutionoriginal distributiondeconvolution (d) Logarithmic distribution Fig. 3: Monte Carlo simulation of deconvolution. We recover random distributions, convolved with a Gaussian.small correction to the warp parameters, which we state asthe systematic error in the results.
To test how accurate the extinction map is, we analysedthe map of Green et al. (2015) using the function query ,which returns the standard deviation σ G for a given line ofsight. We calculate a new extinction as A ∗ G ( r ) = A G ( r ) + f ∗ σ G ( r ) , (8)where A G is the extinction given by the map, r is the dis-tance, and f is a factor chosen randomly from a Gaussiandistribution with µ = 0 and σ = 1 .In Fig. 8 we show the relative error of the density δ = ( ρ ( A G ) − ρ ( A ∗ G )) /ρ ( A G ) . For all lines of sight that we tested, the difference is negligible, except for the areain the centre of the Galaxy, which we know is problematic.However, in the outer disc, where we carried out our analy-sis, the extinction is determined quite accurately. We mustof course take into account that we used the map of Bovyet al. (2016), which combines different maps and is less ac-curate and therefore can give different results than Greenet al. (2015) in some areas. Moreover, we estimated onlythe statistical error of the extinction, but we recall that wedo not have information about the systematic error of theextinction map. However, for our purposes, the extinctionmap gives satisfying results in the area we analysed. Thestellar warp has been studied using star counts by manyother authors (López-Corredoira et al. 2002b; Reylé et al.2009; Amôres et al. 2017, and others), therefore this methodis most likely not especially flawed. Article number, page 5 of 19 &A proofs: manuscript no. gaia_maps_final -1 d [kpc] -3 -2 -1 N convolutionoriginal distributiondeconvolution (a) σ π = 0 . mas -1 d [kpc] -3 -2 -1 N convolutionoriginal distributiondeconvolution (b) σ π = 0 . mas -1 d [kpc] -3 -2 -1 N convolutionoriginal distributiondeconvolution (c) σ π = 0 . mas -1 d [kpc] -4 -3 -2 -1 N convolutionoriginal distributiondeconvolution (d) σ π = 0 . mas 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.
In the previous analysis we considered only the thin-discpopulation because the luminosity function presented inSection 3 is calculated in thin-disc regions. However, wecan also analyse high Galactic heights, where the influenceof the thick disc is significant. To test the importance ofthe change in luminosity function, we tested the densitycalculations with a tentative thick-disc luminosity functionthat reduces the number of bright stars. We used the sourcetable of Wainscoat et al. (1992), who give the ratio of allthe components of the Galaxy for all stellar classes. Basedon this comparison, we altered our luminosity function toconstruct a theoretical thick-disc luminosity function, asdepicted in Fig. 9. Then we repeated our calculation withthis new luminosity function. In Fig. 10 we show the re-sult for some lines of sight. In the area where we carriedout the analysis, the difference between the two approaches is clearly visible starting at ∼ kpc. Our density analy-sis is made in the area below 20 kpc, where the differencebetween the two densities is negligible. We note that thisdifference changes with line of sight, which is caused bythe extinction. In the areas where the extinction is signif-icant, the difference between densities derived from thin-and thick-disc luminosity functions is more important, butthese areas are removed from our analysis. Therefore ourmaps are also valid for thick-disc areas.
5. Analysis of the density maps
Recently, similar maps were created by Anders et al. (2019).In their analysis, they used the code
StarHorse , originallydeveloped to determine stellar parameters and distances forspectroscopic surveys (Queiroz et al. 2018). This code com-
Article number, page 6 of 19. Chrobáková, R. Nagy, M. López-Corredoira: Structure of the outer Galactic disc with
Gaia
DR2 Z [ k p c ] [ s t a r s / k p c ] (a) ◦ < φ < ◦ Z [ k p c ] [ s t a r s / k p c ] (b) ◦ < φ < ◦ Z [ k p c ] [ s t a r s / k p c ] (c) ◦ < φ < ◦ Z [ k p c ] [ s t a r s / k p c ] (d) ◦ < φ < ◦ Z [ k p c ] [ s t a r s / k p c ] (e) ◦ < φ < ◦ Z [ k p c ] [ s t a r s / k p c ] (f) ◦ < φ < ◦ Fig. 5: Density maps for various azimuths between ◦ an ◦ . Article number, page 7 of 19 &A proofs: manuscript no. gaia_maps_final
15 10 5 0 5 10 15 20X [kpc]1510505101520 Y [ k p c ] [ s t a r s / k p c ] (a)
15 10 5 0 5 10 15 20Y [kpc]7.55.02.50.02.55.07.510.0 Z [ k p c ] [ s t a r s / k p c ] (b)
15 10 5 0 5 10 15 20X [kpc]7.55.02.50.02.55.07.510.0 Z [ k p c ] [ s t a r s / k p c ] (c) Z [ k p c ] [ s t a r s / k p c ] (d) Fig. 6: Density maps.pares observed quantities to a number of stellar evolution-ary models. It finds the posterior probability over a gridof stellar models, distances, and extinctions. To do this,it needs many priors, including stellar initial mass func-tion, density laws for main Milky Way components, andthe broad metallicity and age of those components. After-wards, the authors applied various criteria on their sampleto choose only accurate results.When we compare our results, we can observe similar struc-tures, except in the area of the Galactic bulge, where ourdata are not reliable and the data of Anders et al. (2019)are much more accurate. However, because data with higherrors in parallax were removed, Anders et al. (2019) wereunable to reach such high distances, which are necessaryto study features of the outer disc such as the flare or thewarp. Another advantage of our method is that we did notassume any priors about the Milky Way. Furthermore, ourdensity maps are a representation of the complete numberof stars per unit volume up to some given absolute magni-tude, taking into account the luminosity function, whereas Anders et al. (2019) gave the stars observed by Gaia, amuch larger number in the solar neighbourhood, thus notuseful to quantify absolute trends in the density distribu-tion. Nevertheless, we consider the results of Anders et al.(2019) very useful because they improve the accuracy of thedata significantly and can be used to study parts of MilkyWay where our data fail.
There has been some discussion about the cut-off in theMilky Way. Some authors have reported to find a cut-offstarting at about 14 kpc from the Galactic centre (Robinet al. 1992, 2003; Minniti et al. 2011). However, Carraro(2014) argued that these finding are erroneous eitherbecause the dataset is biased or because the warp and flareis confused with the cut-off. The absence of the cut-off hasbeen confirmed by several studies (López-Corredoira &Molgó 2014; Sale et al. 2010; Brand & Wouterloot 2007).Our results show that there is no cut-off in the Galactic
Article number, page 8 of 19. Chrobáková, R. Nagy, M. López-Corredoira: Structure of the outer Galactic disc with
Gaia
DR2 mas ]10 [ s t a r s / k p c ] l=[0.0, 5.007] b=[2, 4] original= 0.029 mas= 0.046 mas (a) mas ]10 [ s t a r s / k p c ] l=[304.053, 310.014] b=[-34, -32] original= 0.029 mas= 0.046 mas (b) mas ]10 [ s t a r s / k p c ] l=[63.925, 95.887] b=[80, 82] original= 0.029 mas= 0.046 mas (c) mas ]10 [ s t a r s / k p c ] l=[180.027, 185.028] b=[0, 2] original= 0.029 mas= 0.046 mas (d) mas ]10 [ s t a r s / k p c ] l=[129.41, 134.586] b=[14, 16] original= 0.029 mas= 0.046 mas (e) mas ]10 [ s t a r s / k p c ] l=[348.944, 357.253] b=[-54, -52] original= 0.029 mas= 0.046 mas (f) Fig. 7: Comparison of densities for various lines of sight. Orange and green curves represent the density including thezero-point correction of the parallax, and the blue curve shows the density without this correction.
Article number, page 9 of 19 &A proofs: manuscript no. gaia_maps_final mas ]0.000.010.020.030.040.050.060.070.08 l=[180.027, 185.028] b=[0, 2] (a) mas ]0.000.050.100.150.200.250.300.35 l=[90.344, 95.363] b=[4, 6] (b) mas ]0.00.20.40.60.81.01.21.4 l=[0.0, 5.007] b=[2, 4] (c) mas ]0.0000.0020.0040.0060.0080.010 l=[249.586, 254.679] b=[10, 12] (d) mas ]0.000.010.020.030.040.05 l=[65.488, 70.526] b=[-8, -6] (e) mas ]0.0000.0050.0100.0150.0200.0250.0300.0350.040 l=[120.255, 125.483] b=[16, 18] (f)
Fig. 8: Relative error δ = ( ρ ( A G ) − ρ ( A ∗ G )) /ρ ( A G ) of the density calculated including the standard deviation of theextinction. Article number, page 10 of 19. Chrobáková, R. Nagy, M. López-Corredoira: Structure of the outer Galactic disc with
Gaia
DR2 -7 -6 -5 -4 -3 -2 -1 N interpolation thick diskinterpolation thin diskdata thick diskdata thin disk Fig. 9: Comparison of luminosity functions for the thin andthick disc.disc, at least up to 20 kpc.
We define the solar neighbourhood as the area where 7.5kpc < R < 8.5 kpc and | z | < . kpc and calculate theaverage density in this area. We find ρ (cid:12) = 0 . stars/pc ,which is close to other values in the literature, for example . stars/pc obtained by Chang et al. (2011), who useda three-component model to fit data from 2MASS. Eatonet al. (1984) found ρ (cid:12) = 0 . stars/pc , which is lowerthan our result, but this value is influenced by the rangeof the luminosity function, which is where the differencebetween the values stems from. In our case, we measuredstars with M G < . To describe the radial volume mass density distribution inthe Galactic equatorial plane, we used a modified exponen-tial disc with a deficit of stars in the inner in-plane regionadopted from López-Corredoira et al. (2004) in the follow-ing form: ρ ( R ) = ρ × exp (cid:18) R (cid:12) h r + h r,hole R (cid:12) (cid:19) × exp (cid:18) − Rh r − h r,hole R (cid:19) , (9)where h r is the scale length, h r,hole = 3 . kpc is thescale of the hole, R (cid:12) is the Galactocentric distance ofthe Sun, and R is the Galactocentric distance. We ne-glected the contribution of the thick disc and analysedonly the thin disc. We divided the Galactic equatorialplane into three regions according to the Galactic az-imuth [ − ◦ , − ◦ ] , [ − ◦ , ◦ ] , [15 ◦ , ◦ ] . We focused onthe Galactic equatorial plane, therefore we considered starsin the close vicinity of the plane with a vertical distance | z | < . kpc and R>6 kpc. We fitted the density forvarious azimuths with the corresponding exponential fitsbased on Eq. (9). The scale length slightly depends on theGalactic azimuth; it reaches the highest value for the Sun- anticentre direction and φ = +30 ◦ , h r = 2 . ± . kpc,and h r = 2 . ± . kpc. On the other hand, the low-est value of the scale length is h r = 1 . ± . kpc for φ = − ◦ . This results in an average of h r = 2 . ± . kpc, with small dependence on the azimuth. We can com-pare the results with published papers. López-Corredoira &Molgó (2014) used SDSS-SEGUE (Sloan Digital Sky Survey- Sloan Extension for Galactic Understanding and Explo-ration) data to investigate the density distribution in theGalactic disc. They obtained the scale length for the thickand for the thin disc, h r,thin = 2 . kpc and h r,thick = 2 . kpc for the azimuth φ ≤ ◦ , which is consistent with ourresults. Li et al. (2019) studied OB stars using Gaia DR2data and the derived scale length of the Galactic disc, andfound h r = 2 . ± . kpc, which is in accordance with ourresults.We also plot the dependence of the density in the Galac-tic equatorial plane on azimuth for various values of Galac-tocentric distance in Fig. 11. The density is slightly depen-dent on the Galactic azimuth for all radii, but this depen-dence is very small. An analysis of the scale height and itscorresponding flare will be given in a forthcoming paper(Nagy et al. 2020, in preparation). The density maps (Fig. 5) directly show a northern warpin azimuth ◦ and southern warp in azimuth ◦ . Here,we analyse these structures in greater detail. We removedthe azimuths ◦ < φ < ◦ and radii R<6 kpc from ouranalysis because these data have low quality and influencethe results negatively.We calculated the average elevation above the plane z w as z w = (cid:82) z max z min ρz d z (cid:82) z max z min ρ d z (10)and fit this quantity with models of the warp. In our firstapproach, we used the model by López-Corredoira et al.(2002b, Eq. 20), z w = [ C w R ( pc ) (cid:15) w sin ( φ − φ w ) + 17] pc . (11)The 17 pc term compensates for the elevation of theSun above the plane (Karim & Mamajek 2017). C w , (cid:15) w ,and φ w are free parameters of the model, which were fittedto our data. An asymmetry is observed between the north-ern and southern warp for the gas (Voskes & Butler Burton2006) and for the young population (Amôres et al. 2017),therefore we also explore the northern and southern warpseparately here. The fit of our data yields maximum am-plitudes z w = 0 . kpc for the northern and z w = − . kpc for the southern warp, both at a distance R=[19.5,20]kpc, revealing a small asymmetry between the north andsouth. For the fit, we used the function curve fit from thepython SciPy package, which uses non-linear least squaresto fit a function to data. The parameters of the best fit forthis model for the whole dataset are
Article number, page 11 of 19 &A proofs: manuscript no. gaia_maps_final d [ kpc ]10 [ s t a r s / k p c ] l=[170.821, 175.997] b=[14, 16] thin diskthick disk (a) d [ kpc ]10 [ s t a r s / k p c ] l=[184.589, 191.426] b=[42, 44] thin diskthick disk (b) d [ kpc ]10 [ s t a r s / k p c ] l=[180.247, 185.254] b=[2, 4] thin diskthick disk (c) d [ kpc ]10 [ s t a r s / k p c ] l=[183.116, 189.22] b=[-36, -34] thin diskthick disk (d) d [ kpc ]10 [ s t a r s / k p c ] l=[160.024, 165.025] b=[-2, 0] thin diskthick disk (e) d [ kpc ]10 [ s t a r s / k p c ] l=[195.953, 206.267] b=[60, 62] thin diskthick disk (f) Fig. 10: Densities for different lines of sight. Blue curves are densities calculated with the thin-disc luminosity function.Orange curves are densities calculated with the thick-disc luminosity function.
Article number, page 12 of 19. Chrobáková, R. Nagy, M. López-Corredoira: Structure of the outer Galactic disc with
Gaia
DR2 R [ kpc ]5.56.06.57.07.58.08.5 l o g () [ s t a r s / k p c ] fit - [315.0 , 345.0 ]fit - [345.0 , 15.0 ]fit - [15.0 , 45.0 ]data - [315.0 , 345.0 ]data - [345.0 , 15.0 ]data - [15.0 , 45.0 ] Fig. 11: Dependence of the density on azimuth nearthe centre-Sun-anticentre direction for various values ofGalactocentric distance. The data points are obtained asweighted mean in bins of size 1 kpc in R and 0.4 kpc in | z | . Only bins with a number of points N ≥ points areplotted. C w = 1 . · − pc ± . · − pc( stat. ) − . · − pc( syst. ) ,(cid:15) w = 2 . ± . stat. ) + 0 . syst. ) , (12) φ w = − . ◦ ± . ◦ ( stat. ) + 4 . ◦ ( syst. ) . Here, the error of c w stands for the error of the amplitudealone, without the variations of (cid:15) w and φ w . The plot of theresults is shown in Fig. 12, where we show the comparisonof minimum and maximum value of z w ( R ) . The averageelevation of the plane is highest for azimuths [60 ◦ , ◦ ] and [90 ◦ , ◦ ] in most of the cases, whereas the minimum isreached for azimuths [240 ◦ , ◦ ] in most of the cases. Aslight asymmetry between the northern and southern warpis also clearly visible.Another approach that we used is based on the workof Levine et al. (2006), who studied the vertical structureof the outer disc of the Milky Way by tracing neutral hy-drogen gas. They analysed the Galactic warp using a Lombperiodogram analysis. They concluded that the first twoFourier modes are the strongest modes. We use the expres-sion derived by Levine et al. (2006) in the following form: z w = z + z · sin ( φ − φ ) + z · sin (2 φ − φ ) , (13)where z w is the average elevation above the plane, z i for i ∈ (0 , , are the amplitudes of the warp, φ i for i ∈ (1 , are the phases. The dependence of the amplitudes of thewarp on the Galactocentric distances is z i = k + k · ( R − R k ) + k · ( R − R k ) for i = 0 , , , (14)where k i and R k are free parameters of the fit. We fitted ourdata with Eqs. (13) and (14) for various values of Galacto-centric distances R < kpc. We plot the data and the fitsfor R ∈ (13 . , . , . kpc in Fig. 13. Fig. 14 showsthe azimuth of the maximum and minimum of the Galac-tic warp as a function of the Galactocentric distance. In R [ kpc ]0.40.20.00.20.4 z w [ k p c ] z w maximum z w minimum [ d e g ] [ d e g ] Fig. 12: Minimum and maximum average elevation of theplane as a function of radius. The warp fit is based on Eq.(11), and the error bars represent the uncertainty in thedistance in the Lucy method.our analysis, we excluded data for the Galactic azimuths φ ∈ (120 ◦ , ◦ ) because of the high error values in ourdata. We used a ◦ binning in azimuth. Fig. 13 shows thatthe data for ◦ < φ < ◦ are somewhat noisy, whichcan be caused by problems with extinction or with the Lucymethod in a particular line of sight. Therefore we tested afit without these points, which turned out to produce aninsignificant difference. For instance, the minimum ampli-tude obtained without these points changed by 10% in theworse case, and the maximum amplitude changed by 2%.Figs. 13 and 14 clearly show that the warp is presentin our analysis. The azimuth of the maximum of the warp(the northern warp) is an increasing function of the Galac-tocentric distance ( ◦ < φ < ◦ ). On the other hand, theazimuth of the minimum of the warp is in ◦ < φ < ◦ and corresponds to the southern warp. The strongest de-viation of the average elevation of the Galactic plane fromthe Galactic equatorial plane rises with Galactocentric dis-tance. The highest amplitude of the northern and southernwarp is z w = 0 . kpc and z w = − . kpc, respectively.An asymmetrical warp is clearly present.The value of the line of nodes from the fit is φ = − . ◦ . We plot the changes in amplitude of the Galac-tic warp fit [Eq.(14)] with Galactocentric distance in Fig.15. Similar results were obtained by Li et al. (2019), whoused OB stars of Gaia DR2 to measure the warp. They fittheir data with a sinusoidal function similar to ours andobtained a warp with a mean magnitude up to z = 0 . kpc. However, they did not account for the asymmetry ofthe warp, therefore they found the same result for the northand south. Chen et al. (2019) used Cepheids from the WISE(Wide-field Infrared Survey Explorer) catalogue and tracedthe warp up to R=20 kpc. Their results show a warp ex-tended up to | z | = 1 . kpc, which we cannot confirm usingthe whole population. Poggio et al. (2018) studied the kine-matics of the Milky Way using Gaia DR2 and found thewarp up to 7 kpc from the Sun. This agrees with our re-sults, but we show that the warp extends to a higher radius,at least up to 20 kpc. In Fig. 16, we compare the maximumamplitudes of our model with other works. We obtain a very Article number, page 13 of 19 &A proofs: manuscript no. gaia_maps_final deg ]0.60.40.20.00.20.40.6 z w [ k p c ]
300 350[ deg ] The Galactic Warp - fitThe Galactic Warp - data (a) R = . kpc deg ]0.60.40.20.00.20.40.6 z w [ k p c ]
300 350[ deg ] The Galactic Warp - fitThe Galactic Warp - data (b) R = . kpc deg ]0.60.40.20.00.20.40.6 z w [ k p c ]
300 350[ deg ] The Galactic Warp - fitThe Galactic Warp - data (c) R = . kpc Fig. 13: Average elevation of the plane as a function of the Galactic azimuth for various values of the Galactocentricdistance. Red markers represent values of binned data, and the blue line represents a fit to the data.
10 12 14 16 18 R [ kpc ]0.40.20.00.20.40.6 z w [ k p c ]
324 321 320 318 317 316 315 314 313 313 31256 55 55 54 54 53 53 53 53 53 52 [ d e g ] [ d e g ] Fig. 14: Minimum and maximum of the average elevationof the plane as a function of Galactocentric distance. Thewarp fit is based on Eq. 13. The colours code the azimuth ofthe minimum and maximum of the warp fit, and the errorbars represent the uncertainty on the distance in the Lucymethod.low amplitude, especially in comparison with Cepheids. Onthe other hand, the closest result is that of Chen et al.(2019), who used OB stars from Gaia DR2. This significantdifference between the amplitude of various populations isin favour of the formation of the warp through accretiononto the disc (López-Corredoira et al. 2002a), which causesthe gas and young stars to warp more strongly than theremaining population.Momany et al. (2006) studied the stellar warp using 2MASSred clump and red giant stars, selected at fixed heliocentricdistances of 3, 7, and 17 kpc. They found a rather sym-metric warp and argued that a symmetric warp can be ob-served as asymmetric for two reasons. First, the Sun is notlocated at the line of nodes, and second, the northern warpis located behind the Norma-Cygnus arm, which can causevariation in extinction that can produce an apparent asym-metric warp. As for the first point, the position of Sun onthe line of nodes is a problem when we observe the warp
10 12 14 16 18 R [ kpc ]0.20.10.00.10.20.3 z w [ k p c ] z z z Fig. 15: Changes of the amplitudes of the Galactic warp fitaccording to Eqs. 13 and 14.
10 12 14 16 18 20R [kpc]0.00.20.40.60.81.01.21.4 z [ k p c ] our work - whole populationour work - young populationLópez-Corredoira et al. (2002b)Yusifov (2004), Pulsars Reylé et al. (2009)Chen et al. (2019), CepheidsSkowron et al. (2019a), CepheidsLi et al. (2019), OB stars Fig. 16: Comparison of maximum amplitudes of our model(based on Eq. (13)) with other works.
Article number, page 14 of 19. Chrobáková, R. Nagy, M. López-Corredoira: Structure of the outer Galactic disc with
Gaia
DR2
10 8 6 4 2M10 -5 -4 -3 -2 -1 N Fig. 17: Luminosity function used in the Eq. 6 for the anal-ysis of the young population.at a fixed distance. However, we have a 3D distribution,which ensures that the position from which we look doesnot influence how we perceive the warp. As for the secondremark, as we showed in Section 2 that the extinction is de-termined quite accurately by the extinction map of Greenet al. (2015). However, some variations might influence thefinal shape of the warp and may not have been taken intoaccount, therefore we need to keep that in mind when weinterpret our results.
In this section, we apply the previous analysis to theyoung population. To do so, we only chose stars brighterthan an absolute magnitude M G = − (see the lumi-nosity function in Fig. 17) and repeated all the stepsas described in Section 4.3. Then we produced den-sity maps and analysed the scale length and the warpof this population using methods from Sections 5.4 and 5.5.The exponential fits of the density for the young popula-tion yield h r = 2 . ± . kpc for φ = 30 ◦ , h r = 1 . ± . kpc for φ = 0 ◦ , and h r = 2 . ± . kpc for φ = 30 ◦ ,which is similar to the whole population. This results in h r = 2 . ± . kpc on average. The variation with azimuthis still insignificant, as in the case of the entire population.Fig. 18 shows that the variation of density with azimuth isalso negligible in the case of the young population.For the warp, as previously, we removed the azimuths ◦ < φ < ◦ from the analysis. The fit of Eq. (11)to the young population yields C w = 4 . · − pc ± . · − pc( stat. )+ 5 . · − pc( syst. ) ,(cid:15) w = 3 . ± . stat. ) − . syst. ) , (15) φ w = − . ◦ ± . ◦ ( stat. ) − . ◦ ( syst. ) . We also repeated the analysis with the approach us-ing Eq. (13). Fig. 19 presents the Galactic warp of the R [ kpc ]2.02.53.03.54.04.55.05.5 l o g () [ s t a r s / k p c ] fit - [315.0 , 345.0 ]fit - [345.0 , 15.0 ]fit - [15.0 , 45.0 ]data - [315.0 , 345.0 ]data - [345.0 , 15.0 ]data - [15.0 , 45.0 ] Fig. 18: Same as Fig. 11, but the young population alone isconsidered.young stellar population for various Galactocentric dis-tances, and Fig. 20 shows the amplitudes of the fits ofthe Galactic warp and the azimuth of the maximum andminimum. In this case, the warp of the young stellar pop-ulation is stronger than the case considering all stars inour dataset. The azimuth of the maximum of the northernwarp is an increasing function of the Galactocentric dis-tance ( ◦ < φ < ◦ ), and the azimuth of the minimum ofthe warp is in ◦ < φ < ◦ . The highest amplitude ofthe northern and the southern warp is z w = 0 . kpc and z w = − . kpc, respectively. For the line of nodes, we find φ = − . ◦ , which agrees with the whole population.Chen et al. (2019) used Cepheids from the WISE surveyand a number of optical surveys to measure the warp, andSkowron et al. (2019a) used Cepheids from the OGLE cat-alogue supplemented by other surveys. Chen et al. (2019)obtained a rather symmetric warp with an amplitude ofabout 1.5 kpc in R=20 kpc. Skowron et al. (2019a) ob-tained a similar result with an amplitude 0.74 kpc in R=15kpc. These values are much higher than our findings, whichis probably due to differences in the population: our youngpopulation is older than the Cepheids. In Fig. 21 we plot thevariation of the line of nodes with radius for the whole andthe young population, compared with other works. We usetwo different methods to plot the line of nodes for our work.First, we plot the angle φ w for Eq. (12). Another methodis to use the Eq. (11) to find the value of the angle φ when z w = 0 . We would expect that our young population liesbetween the total population and the young Cepheids. Thisis true only for R>12 kpc. At shorter distances, the warpis not very strong and is more difficult to detect, thereforethe error bars are larger in this area. Moreover, the errorbars of the young populations are very large because of thelower number of stars in the sample combined with possibleproblems in determining extinction. For these reasons, thevalue of the line of nodes for R<12 kpc is rather unreliable.
6. Conclusions
We produced density maps from Gaia DR2 data andanalysed them to study the Galactic warp. The density
Article number, page 15 of 19 &A proofs: manuscript no. gaia_maps_final deg ]0.60.40.20.00.20.40.6 z w [ k p c ]
300 350[ deg ] The Galactic Warp - fitThe Galactic Warp - data (a) R = . kpc deg ]0.60.40.20.00.20.40.6 z w [ k p c ]
300 350[ deg ] The Galactic Warp - fitThe Galactic Warp - data (b) R = . kpc deg ]0.60.40.20.00.20.40.6 z w [ k p c ]
300 350[ deg ] The Galactic Warp - fitThe Galactic Warp - data (c) R = . kpc Fig. 19: Dataset containing a young population of stars. The average elevation of the plane as a function of the Galacticazimuth for various values of the Galactocentric distance. Red markers represent values of binned data, and the blue linerepresents a fit to the data.
10 12 14 16 18 R [ kpc ]0.60.40.20.00.20.40.6 z w [ k p c ]
265 294 301 305 308 310 311 313 314 314 31550 51 51 52 53 53 53 53 54 54 54 [ d e g ] [ d e g ] Fig. 20: Minimum and maximum of the average elevationof the plane as a function of the Galactocentric distance.The warp fit is based on Eq. (13). The colours code theazimuth of the minimum and the maximum of the warpfit. The dataset containing a young population of stars isconsidered, and the error bars represent the uncertainty onthe distance in the Lucy method.maps directly show a northern warp in the azimuths ◦ < φ < ◦ and a southern warp in the azimuths ◦ < φ < ◦ . Our maps reach a Galactocentric radiusof 20 kpc, and we note that up to this distance, thedensity decreases exponentially and we do not observe acut-off. Another feature in the density maps is a Galacticflare, that is, an increase in scale height towards theouter Galaxy. The analysis of the flare will be given ina forthcoming paper (Nagy et al. 2020, in preparation).We used the maps to calculate the scale length, where wefind h r = 2 . ± . kpc, with a small dependence of h r from the Galactic azimuth. The lowest value of h r that wefound is . ± . kpc for φ ≈ ±− ◦ and the highestvalue is . ± . kpc for the Sun-anticentre directionand . ± . φ ≈ ± ◦ .
11 12 13 14 15 16 17 18 19 R [ kpc ]402002040 [ d e g ] our work - all stars - for z w =0our work - young population - for z w =0our work - all stars - w (Eq.12)our work - young population - w (Eq.15) López-Corredoira et al. (2002)Yusifov (2004)Skowron et al. (2019)Chen et al. (2019) Fig. 21: Comparison of line of nodes for our model (basedon Eq. 11) with other works. We use two different methodsto plot the line of nodes for our work. First, we plot theangle φ w for Eq. (12). Another method is to use the Eq.(11) to find the value of the angle φ when z w = 0 .From our maps, we calculated the average elevation ofthe plane and fitted it with different warp models. We fit-ted the northern and southern warp separately with a sim-ple sinusoidal model, and we found a small asymmetry:the northern warp reaches an amplitude . kpc for theazimuth ◦ < φ < ◦ and the southern warp reaches − . kpc for the azimuth ◦ < φ < ◦ , both atR=[19.5,20.0] kpc. Then we fitted the warp with a modelcombining two sinusoids to detect the asymmetry withoutassuming its existence, and we found values of amplitude ∼ . for the northern and ∼ − . for the southern warpboth at R=[19.5,20.0] kpc, revealing the asymmetry foundwith the previous approach. The azimuths of the warp max-imum and minimum for this model are ◦ < φ < ◦ and ◦ < φ < ◦ , respectively. In terms of Galactocen-tric radius, we find that warp starts to manifest itself fromabout 12 kpc and extends at least up to 20 kpc. We repeatedthis analysis on the young population, where we find that Article number, page 16 of 19. Chrobáková, R. Nagy, M. López-Corredoira: Structure of the outer Galactic disc with
Gaia
DR2 it follows the result for the whole population, but reaches ahigher amplitude of warp and similar values of scale height.The comparison of our amplitude of the warp with otherworks showed that we obtain a significantly lower amplitudethan an analysis carried out with very young stars such asCepheids. This supports the formation of the warp throughaccretion onto the disc (López-Corredoira et al. 2002a).A future analysis of the next Gaia data release combinedwith the deconvolution method based on Lucy’s method ofinversion, as described in Section 4.1, will allow us to ex-plore distances larger than 20 kpc. The future data releasewill provide a much deeper magnitude limit and much lowerparallax errors, which will allow us to extend the range ofGalactocentric distances and study the morphology of thedisc and of the stellar halo at very large distances.
Acknowledgements.
We thank the anonymous referee for helpful com-ments, which improved this paper, and Astrid Peter (language edi-tor of A&A) for proof-reading of the text. ZC and MLC were sup-ported by the grant PGC-2018-102249-B-100 of the Spanish Min-istry of Economy and Competitiveness (MINECO). RN was sup-ported by the Scientific Grant Agency VEGA No. 1/0911/17. Thiswork made use of the IAC Supercomputing facility HTCondor(http://research.cs.wisc.edu/htcondor/), partly financed by the Min-istry of Economy and Competitiveness with FEDER funds, codeIACA13-3E-2493. This work has made use of data from the Euro-pean Space Agency (ESA) mission
Gaia ( ), processed by the Gaia
Data Processing and AnalysisConsortium (DPAC, ). Funding for the DPAC has been provided by nationalinstitutions, in particular the institutions participating in the
Gaia
Multilateral Agreement. The reduced catalogue of Gaia with m G < was produced by Pedro Alonso Palicio. Article number, page 17 of 19 &A proofs: manuscript no. gaia_maps_final
References
Amôres, E. B., Robin, A. C., & Reylé, C. 2017, A&A, 602, A67Anders, F., Khalatyan, A., Chiappini, C., et al. 2019, A&A, 628, A94Arenou, F., Luri, X., Babusiaux, C., et al. 2018, A&A, 616, A17Bahcall, J. N. 1986, ARA&A, 24, 577Bahcall, J. N. & Soneira, R. M. 1980, ApJS, 44, 73Balázs, L. G. 1995, Inverse Problems, 11, 731Battaner, E., Florido, E., & Sanchez-Saavedra, M. L. 1990, A&A, 236,1Bovy, J., Rix, H.-W., Green, G. M., Schlafly, E. F., & Finkbeiner,D. P. 2016, ApJ, 818, 130Bovy, J., Rix, H.-W., Liu, C., et al. 2012, ApJ, 753, 148Brand, J. & Wouterloot, J. G. A. 2007, A&A, 464, 909Carney, B. W. & Seitzer, P. 1993, AJ, 105, 2127Carraro, G. 2014, in IAU Symposium, Vol. 298, Setting the scene forGaia and LAMOST, ed. S. Feltzing, G. Zhao, N. A. Walton, &P. Whitelock, 7–16Chandresekhar, S. & Münch, G. 1951, ApJ, 113, 150Chang, C.-K., Ko, C.-M., & Peng, T.-H. 2011, ApJ, 740, 34Chen, B., Stoughton, C., Smith, J. A., et al. 2001, ApJ, 553, 184Chen, X., Wang, S., Deng, L., et al. 2019, Nature Astronomy, 3, 320Craig, I. J. D. & Brown, J. C. 1986, Inverse problems in astronomy.A guide to inversion strategies for remotely sensed data (AdamHilger, Bristol, UK)Debattista, V. P. & Sellwood, J. A. 1999, ApJ, 513, L107Dehnen, W. 1998, AJ, 115, 2384Drimmel, R., Cabrera-Lavers, A., & López-Corredoira, M. 2003, A&A,409, 205Eaton, N., Adams, D. J., & Giles, A. B. 1984, MNRAS, 208, 241Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A,616, A1Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A,595, A1Green, G. M., Schlafly, E. F., Finkbeiner, D., et al. 2018, MNRAS,478, 651Green, G. M., Schlafly, E. F., Finkbeiner, D. P., et al. 2015, ApJ, 810,25Hendy, Y. H. M. 2018, NRIAG Journal of Astronomy and Geophysics,7, 180Jurić, M., Ivezić, Ž., Brooks, A., et al. 2008, ApJ, 673, 864Karim, M. T. & Mamajek, E. E. 2017, MNRAS, 465, 472Kerr, F. J. 1957, AJ, 62, 93Kim, J. H., Peirani, S., Kim, S., et al. 2014, ApJ, 789, 90Levine, E. S., Blitz, L., & Heiles, C. 2006, ApJ, 643, 881Li, C., Zhao, G., Jia, Y., et al. 2019, ApJ, 871, 208Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, A&A, 616,A2Liu, C., Xu, Y., Wan, J.-C., et al. 2017, Research in Astronomy andAstrophysics, 17, 096López-Corredoira, M., Abedi, H., Garzón, F., & Figueras, F. 2014,A&A, 572, A101López-Corredoira, M., Betancort-Rijo, J., & Beckman, J. E. 2002a,A&A, 386, 169López-Corredoira, M., Cabrera-Lavers, A., Garzón, F., & Hammers-ley, P. L. 2002b, A&A, 394, 883López-Corredoira, M., Cabrera-Lavers, A., Gerhard, O. E., & Garzón,F. 2004, A&A, 421, 953Lopez-Corredoira, M., Garzon, F., Wang, H. F., et al. 2020, arXive-prints, arXiv:2001.05455López-Corredoira, M., Hammersley, P. L., Garzón, F., Simonneau, E.,& Mahoney, T. J. 2000, MNRAS, 313, 392López-Corredoira, M. & Molgó, J. 2014, A&A, 567, A106López-Corredoira, M. & Sylos-Labini, F. 2019, A&A, 621, A48Lucy, L. B. 1974, ApJ, 79, 745Luri, X., Brown, A. G. A., Sarro, L. M., et al. 2018, A&A, 616, A9Majewski, S. R. 1993, ARA&A, 31, 575Marshall, D. J., Robin, A. C., Reylé, C., Schultheis, M., & Picaud, S.2006, A&A, 453, 635Minniti, D., Saito, R. K., Alonso-García, J., Lucas, P. W., & Hempel,M. 2011, ApJ, 733, L43Momany, Y., Zaggia, S., Gilmore, G., et al. 2006, A&A, 451, 515–538Nagy, R. et al. 2020, in preparationOort, J. H., Kerr, F. J., & Westerhout, G. 1958, MNRAS, 118, 379Paul, E. R. 1993, The Milky Way Galaxy and Statistical Cosmology,1890-1924Poggio, E., Drimmel, R., Lattanzi, M. G., et al. 2018, MNRAS, 481,L21Queiroz, A. B. A., Anders, F., Santiago, B. X., et al. 2018, MNRAS,476, 2556 Reylé, C., Marshall, D. J., Robin, A. C., & Schultheis, M. 2009, A&A,495, 819Riess, A. G., Casertano, S., Yuan, W., et al. 2018, ApJ, 861, 126Robin, A. C., Creze, M., & Mohan, V. 1992, ApJ, 400, L25Robin, A. C., Marshall, D. J., Schultheis, M., & Reylé, C. 2012, A&A,538, A106Robin, A. C., Reylé, C., Derrière, S., & Picaud, S. 2003, A&A, 409,523Rybizki, J., Demleitner, M., Fouesneau, M., et al. 2018, PASP, 130Sale, S. E., Drew, J. E., Knigge, C., et al. 2010, MNRAS, 402, 713Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525Schönrich, R. & Dehnen, W. 2018, MNRAS, 478, 3809Skowron, D. M., Skowron, J., Mróz, P., et al. 2019a, Science, 365, 478Skowron, D. M., Skowron, J., Mróz, P., et al. 2019b, Acta Astron., 69,305Stassun, K. G. & Torres, G. 2018, ApJ, 862, 61Turchin, V. F., Kozlov, V. P., & Malkevich, M. S. 1971, Soviet PhysicsUspekhi, 13, 681Voskes, T. & Butler Burton, W. 2006, arXiv e-prints, astroVozoff, K. & Jupp, D. L. B. 1975, Geophysical Journal of the RoyalAstronomical Society, 42, 977Wainscoat, R. J., Cohen, M., Volk, K., Walker, H. J., & Schwartz,D. E. 1992, APJS, 83, 111Yusifov, I. 2004, in The Magnetized Interstellar Medium, ed.B. Uyaniker, W. Reich, & R. Wielebinski, 165–169Zinn, J. C., Pinsonneault, M. H., Huber, D., & Stello, D. 2019, ApJ,878, 136
Article number, page 18 of 19. Chrobáková, R. Nagy, M. López-Corredoira: Structure of the outer Galactic disc with
Gaia
DR2
Appendix A: Lucy’s method for the inversion ofFredholm integral equations of the first kind
The inversion of Fredholm integral equations of the firstkind such as Eq. (2) is ill-conditioned. Typical analyticalmethods for solving these equations (Balázs 1995) cannotachieve a good solution because the kernel is sensitive tothe noise of the star counts (Craig & Brown 1986, chap-ter 5). Because the functions in these equations have astochastic rather than analytical interpretation, it is to beexpected that statistical inversion algorithms are more ro-bust (Turchin et al. 1971; Vozoff & Jupp 1975; Balázs 1995).These statistical methods include the iterative method ofLucy’s algorithm (Lucy 1974; Turchin et al. 1971; Balázs1995; López-Corredoira et al. 2000), which is appropriatehere. Its key feature is the interpretation of the kernel as aconditioned probability and the application of Bayes’ the-orem.In Eq. (2), N ( π ) is the unknown function, and the kernelis G ( x ) , whose difference x is conditioned to the parallax π (cid:48) . The inversion is carried out as N ( π ) = lim n →∞ N n ( π ) , (A.1) N n +1 ( π ) = N n ( π ) (cid:82) ∞ N ( π (cid:48) ) N n ( π (cid:48) ) G π (cid:48) ( π − π (cid:48) ) dπ (cid:48) (cid:82) ∞ G π (cid:48) ( π − π (cid:48) ) dπ (cid:48) , (A.2) N n ( π ) = (cid:90) ∞ N n ( π (cid:48) ) G π (cid:48) ( π − π (cid:48) ) dπ (cid:48) . (A.3)The iteration converges when = N n ( π ) ≈ N ( π ) ∀ π , thatis, when N n ( π ) ≈ N ( π ) ∀ π . The first iterations produce aresult that is close to the final answer, with the subsequentiterations giving only small corrections. In our calculation,we set as initial function of the iteration N ( π ) = N ( π ) , and we carry out a number of iterations until the Pearson χ test N p − N p − (cid:88) j =2 [ N n ( π j ) − N ( π j )] N n ( π j ) , (A.4)reaches the minimum value. Further iterations would enterwithin the noise.This algorithm has a number of beneficial properties(Lucy 1974, 1994): all the functions are defined as beingpositive, the likelihood increases with the number of itera-tions, the method is insensitive to high-frequency noise in N ( π ) , and so on. We note, however, that precisely becausethis method only works when N are positive functions, itdoes not work with negative ones.are positive functions, itdoes not work with negative ones.