© Ž. Chrobáková et al. 2022

1. Introduction

Although we have an extensive knowledge of our Galaxy, there are still many aspects for which our understanding of the Milky Way is incomplete. Wide discussions about features such as warp, flare, or cut-off show that there is much more to be learned about the Milky Way. To this end, the Gaia Early Data Release 3 Gaia Collaboration (2021b), EDR3) presents an opportunity to study the Galaxy in greater detail than ever before. With its precise positional, proper motions, radial velocity measurements, and distance determinations for millions of stars, it offers the most accurate information about our Galaxy to date; ideal for making advances in all branches of Galactic astrophysics.

The Galactic warp is a well-known feature of the Galactic disc, but its shape is constrained only roughly and there is no consensus on the mechanism causing it. Some of the theories include accretion of intergalactic matter onto the disc (López-Corredoira et al. 2002a), a misaligned halo (Ostriker & Binney 1989), interaction with satellites (García-Ruiz et al. 2002), or an intergalactic magnetic field (Battaner et al. 1990). Currently the kinematical information about the warp is not sufficient to constrain the formation models. Chrobáková et al. (2020) compared the warp of the whole stellar population with the warp of Cepheids and suggest that warp is dependent on the age of the studied population, which would indicate that warp is caused by a non-gravitational mechanism. A similar conclusion was reached by Wang et al. (2020) using LAMOST DR4. In this paper, we recalculate the Galactic warp using the most recent Gaia EDR3 data and compare this with the warp of supergiants, to test this hypothesis.

The flare is an increase in the scale height of the Galactic disc with galactocentric radius, detected in both the gaseous and stellar components. Grabelsky et al. (1987) and May et al. (1997) confirmed flare in the outer disc by tracing molecular clouds, Sánchez-Salcedo et al. (2008) modelled the flare applying modified Newtonian dynamics (MOND), and Narayan & Jog (2002) treated HI, H₂, and stars as gravitationally coupled components of the disc, calculating scale heights for all three components and giving predictions for the flare that matched observations very well. Stellar flare was studied by several authors as well (Alard 2000; López-Corredoira et al. 2002b; Momany et al. 2006). López-Corredoira & Molgó (2014) studied 3D stellar distribution using Sloan Extension for Galactic Understanding and Exploration (SEGUE) data, finding that flare is a prominent feature for Galactocentric distance of R ≳ 15 kpc. Yusifov (2004) reached a similar conclusion studying pulsars. Bovy et al. (2016b) studied the structure of the Galactic disc, distinguishing various stellar populations using APOGEE survey data covering Galactocentric distances of R < 15 kpc. They did not find stellar flaring in the high-[α/Fe] mono-abundance population, while the low-[α/Fe] mono-abundance populations exhibit Galactic flare.

To explain the dependence of the disc thickness on azimuth, Kalberla et al. (2007) modelled it with a ring of dark matter embedded in the disc, and Saha et al. (2009) applied a lopsided dark matter halo. López-Corredoira & Betancort-Rijo (2009) proposed accretion of intergalactic matter onto the disc as a possible mechanism to explain both the flare and its dependence on azimuth.

More recently, the flare was studied with Gaia DR2 using OB stars (Li et al. 2019) or with LAMOST (Wang et al. 2018) and Cepheids (Feast et al. 2014). Possible existence of the flare was also explored in the thick disc (López-Corredoira & Molgó 2014; Mateu & Vivas 2018; Wang et al. 2018). Yu et al. (2021a) studied the warp and the flare traced by OB stars using LAMOST DR5 data and Yu et al. (2021b) investigated the Galactic disc using LAMOST and Gaia Red Clump Sample VII.

In this paper, we use the Gaia EDR3 data to study the warp and flare. We are especially interested in the properties of these features for various stellar populations, and therefore we analyse separately the population of supergiants. The paper is structured as follows. In Sect. 2 we describe our data selection and the extinction map used. In Sect. 3 we present the method used to calculate the density distribution and explain how we chose the sample of supergiants. Section 4 is dedicated to the analysis of the warp and Sect. 5 deals with the analysis of the flare. In Sect. 6 we conclude the paper.

2. Data Selection

We used the Gaia EDR3 data, collected during the first 34 months of observations. We are interested in sources with G-band (330–1050 nm) magnitude. The photometric uncertainties were ∼0.3 mmag for G < 13, 1 mmag at G = 17, and 6 mmag at G = 20 mag. We chose a magnitude up to G = 19, where the catalogue was sufficiently complete (Fabricius et al. 2021). More details on the catalogue validation are available in Fabricius et al. (2021). To ensure the quality of the dataset, we applied several constraints. To select sources with good astrometry, we only chose data with five- and six-parameter solutions, satisfying the following condition on renormalised unit weight error (RUWE):

$$\begin{array}{c}\hfill \mathtt{RUWE}<1.4\end{array}$$(1)

as suggested by Lindegren et al. (2021). In addition, following Gaia Collaboration (2021a), we applied the following constraint to ensure that the sources had good photometry:

$$\begin{array}{cc}& 0.01+0.039(\mathrm{B}P-RP)<{log}_{10}\mathtt{excess\_flux}\hfill \\ \hfill & <0.12+0.039(\mathrm{B}P-RP),\hfill \end{array}$$(2)

where BP − RP denotes the colour index, with the G_BP band covering the range 330–680 nm and the G_RP band covering the range 640–1050 nm, and excess_flux is the BP and RP flux excess, corrected as suggested by Riello et al. (2021). Moreover, we applied:

$$\begin{array}{cc}& \mathtt{phot\_g\_mean\_flux\_over\_error}>50,\hfill \\ \hfill & \mathtt{phot\_rp\_mean\_flux\_over\_error}>20,\hfill \\ \hfill & \mathtt{phot\_bp\_mean\_flux\_over\_error}>20,\hfill \end{array}$$(3)

which removed variable stars (Gaia Collaboration 2018). We followed the approach of Chrobáková et al. (2020) and chose data with a parallax in the interval [0, 2] mas, and apparent magnitude in the G-band between G = 12 and G = 19. We corrected G fluxes for six-point sources as well, as suggested by Riello et al. (2021), using codes listed in the appendix of Gaia Collaboration (2021b). We also added the zero-point correction as found by Lindegren et al. (2021), using the publicly available Python package¹, which calculates the zero-point as a function of ecliptic latitude, magnitude, and colour.

Extinction map. In order to estimate the extinction, we used the three-dimensional, full-sky extinction map from Bovy et al. (2016a), using the Python package mwdust. This extinction map is a combination of maps 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 the coefficients (Hendy 2018; Rybizki et al. 2018):

$$\begin{array}{cc}\hfill A_G/A_v& =0.859,\hfill \\ \hfill R_V& =A_v/E(B-V)=3.1.\hfill \end{array}$$(4)

3. Methods

We followed the approach of Chrobáková et al. (2020), who used the fundamental equation of stellar statistics (Chandresekhar & Münch 1951) to derive the stellar density:

$$\begin{array}{cc}\hfill \rho (1/\pi )& =\frac{N(\pi ){\pi }^4}{\mathrm{\Delta }\pi \omega {\int }_{M_{G,\mathrm{low}\lim }}^{M_{G,\mathrm{low}\lim }+1}\mathrm{d}{M}_G\mathrm{\Phi }(M_G)},\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill M_{G,\mathrm{low}\lim }& =m_{G,\mathrm{low}\lim }-5{log}_{10}(1/\pi )-10\hfill \\ \hfill & -A_G(1/\pi ),\hfill \end{array}$$(6)

where N(π) are the star counts, ω is the covered angular surface, Δπ 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. For the luminosity function, we used the values given in Table 1 of Chrobáková et al. (2020) for the whole population; or we determined it following the same method for some sub-samples.

The density determination required that we measure the star counts as a function of distance. However, it is well known that the parallax error grows with the distance from the observer, preventing the precise determination of distances. Therefore, we could not simply use the observed star counts to calculate the density, as at distances higher that roughly 5 kpc they are biased. In order to recover correct star counts, we applied a statistical deconvolution method developed by López-Corredoira & Sylos-Labini (2019) based on Lucy’s method (Lucy 1974). They express the observed number of stars per parallax $\overline{N}(\pi )$ as a convolution of the real number N(π) of stars with a Gaussian function:

$$\begin{array}{c}\hfill \overline{N}(\pi )={\displaystyle {\int }_0^{\infty }\mathrm{d}{\pi }^{\prime }N({\pi }^{\prime })G_{{\pi }^{\prime }}(\pi -{\pi }^{\prime }),}\end{array}$$(7)

where

$$\begin{array}{c}\hfill G_{\pi }(x)=\frac{1}{\sqrt{2\pi }{\sigma }_{\pi }}{e}^{-\frac{x^2}{2{\sigma }_{\pi }^2}}.\end{array}$$(8)

For the error σ_π, we averaged parallax errors of every bin. More details about the method can be found in Chrobáková et al. (2020).

With the corrected star counts, we reveal density distribution up to 20 kpc, which can be seen in Fig. 1. This result is almost identical with the one obtained with Gaia DR2 data (Chrobáková et al. 2020). Similarly, we can see overdensities above the plane for azimuths between 300° and 360°. As commented in Chrobáková et al. (2020), these structures are most likely a contamination, since they disappear after integrating the density through the whole disc.

Fig. 1. Density maps in Galactocentric coordinates at various azimuths for the whole population (Sample 0).

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 also made bins of Δπ = 0.01 mas in parallax. We did not use negative parallaxes because these affect the distribution of parallaxes and statistical properties. 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 simply did not use this information because it is not necessary with this approach. Further details on this method and tests of its possible biases can be found in Chrobáková et al. (2020).

4. Sample definition

We analysed this density distribution to reveal the warp and the flare in the whole stellar population, which we will refer to as Sample 0. Moreover, we separated supergiants from the dataset in order to analyse their density distribution separately and find differences with the respect to the whole population. We used two different approaches to separate the supergiants, as described below.

4.1. Sample 1

In the first approach, we chose only stars which we can be certain are supergiants. Based on the error of parallax, we calculated the interval of possible magnitudes [M_min, M_max] in the G- band for every star:

$$\begin{array}{cc}& M_{min}=m_G-5{log}_{10}(1/(\pi -\omega ))-10-A_G(1/\pi )\hfill \\ \hfill & M_{max}=m_G-5{log}_{10}(1/(\pi +\omega ))-10-A_G(1/\pi ),\hfill \end{array}$$(9)

where π is parallax and ω is parallax error. Then we chose only the stars with −5 < M_min < −10 and −5 < M_max < −10. In the end we had a dataset of 331 546 stars.

As the distribution was not homogeneous, we divided stars in bins as follows. For l < 120° and l > 260° and |b|< 4°, we had bins with Δl = 5° and Δb = 2°. For the same range of l for |b|> 4°, we had Δl = 40° and Δb = 20°. We binned the stars in 120° < l < 260° in one bin. We also made bins Δm = 1.0 in magnitude and Δπ = 0.01 mas in parallax. In Fig. 2a, we show the distribution of the sources selected with this method.

Fig. 2. Distribution of supergiants in galactic coordinates for |Z|< 4 kpc. Left: Sample 1. Right: Sample 2 (see Sects. 4.1 and 4.2 for details).

4.2. Sample 2

The second approach was less strict; we chose sources with absolute magnitude in the G-band within the interval −5 < M < −10, regardless of the error of the magnitude. The advantage of this method is that we were complete, and therefore we could make statistical analysis of the data to investigate the flare. As this sample was distributed more homogeneously, we divided the stars in bins the same way as for the total population. In Fig. 2 (b), we show the distribution of the sources selected with this method. In the end, this sample contained about 1.9 × 10⁶ stars. In Fig. 3 we show the Hertzsprung-Russell diagrams (HRD) of all three of the samples.

Fig. 3. Hertzsprung-Russell diagrams (HRD) for the three samples defined in the text. As the datasets are significantly large, in order to avoid the saturation of the diagram, we plot a randomly chosen sub-sample containing 3000 stars in each case. Top: Sample 0. Middle: Sample 1. Bottom: Sample 2.

5. Warp analysis

The first feature that we studied is the Galactic warp. We removed data for 90° < ϕ < 270° as in this part we did not have enough data. Following the approach of Chrobáková et al. (2020), we calculated the average elevation of the plane as

$$\begin{array}{c}\hfill z_w=\frac{{\int }_{z_{min}}^{z_{max}}\rho z\mathrm{d}z}{{\int }_{z_{min}}^{z_{max}}\rho \mathrm{d}z}\end{array}$$(10)

Then, we fitted z_w with Eq. (11) of Chrobáková et al. (2020), which represents a common warp model:

$$\begin{array}{c}\hfill z_w=[C_{w}R{(pc)}^{{ϵ}_w}sin(\varphi -{\varphi }_w)]\mathrm{pc},\end{array}$$(11)

where C_w, ϵ_w and ϕ_w are free parameters characterizing the warp. We did not account for the height of the Sun on the Galactic plane because some recent studies (e.g., Cheng et al. 2020) suggest that warp starts at smaller radii than previously thought and, therefore, the Solar neighbourhood could be slightly warped. For Sample 0, we find values of warp parameters:

$$\begin{array}{cc}& c_w=1.42\pm 0.15\times {10}^{-8}\mathrm{pc},\hfill \\ \hfill & {ϵ}_w=2.43\pm 0.65,\hfill \\ \hfill & {\varphi }_w=-9.77\pm 7.23\stackrel{°}{.},\hfill \end{array}$$(12)

while for Sample 1, we find:

$$\begin{array}{cc}& c_w=1.92\pm 0.08\times {10}^{-4}\mathrm{pc},\hfill \\ \hfill & {ϵ}_w=1.54\pm 0.18,\hfill \\ \hfill & {\varphi }_w=-8.23\pm 2.95\stackrel{°}{.}\hfill \end{array}$$(13)

and for Sample 2:

$$\begin{array}{cc}\hfill c_w& =4.85\pm 0.18\times {10}^{-5}\mathrm{pc},\hfill \\ \hfill {ϵ}_w& =1.66\pm 0.17,\hfill \\ \hfill {\varphi }_w& =-0.73\pm 2.55\stackrel{°}{.}.\hfill \end{array}$$(14)

The error of c_w stands for the error of the amplitude alone, without the variations of ϵ_w and ϕ_w. For the fit, we used the function curve fit from the python SciPy package, which uses non-linear least squares to fit a function to data.

In Fig. 4 we compare the warp amplitude for supergiants for both samples (Sample 1 and Sample 2). The two samples are in very good agreement, yielding warp amplitude with only negligible differences, showing that the Sample 2 is not significantly contaminated. In Fig. 5, we plot the maximum and minimum amplitudes of warp for the whole population (Sample 0), compared with the supergiants (Sample 2). We confirm that the warp amplitude for the whole population is almost identical to what we found with Gaia DR2 data (Chrobáková et al. 2020). We obtain a maximum warp amplitude of z_w = 0.360 kpc and minimum of z_w = −0.375 kpc at a distance of R = [19.5, 20] kpc, which is slightly higher than the result obtained with Gaia DR2 data, with a small asymmetry between the northern and the southern warp. From comparison with the supergiants, it is clear that the warp amplitude of supergiants is significantly larger that of the whole population, reaching an amplitude twice as large at a distance of R = [19.5, 20] kpc. As supergiants are a young population, a few tens of Myr old on average (e.g., Bouret et al. 2012), whereas the whole population is ∼6−7 Gyr old (e.g., Kilic et al. 2017), there is a clear relationship between the warp amplitude and the age of the studied population, thus confirming conclusions of previous studies (Chrobáková et al. 2020; Wang et al. 2020).

Fig. 4. Comparison of fits of minimum and maximum warp amplitudes for supergiants, chosen by two different approaches (see Sects. 4.1 and 4.2 for details).

Fig. 5. Comparison of fits of maximum and minimum warp amplitudes for the whole population (Sample 0) and supergiants (Sample 2).

6. Flare analysis

In order to investigate the Galactic flare, we considered the density distribution of the Galactic disc as consisting of a thick and a thin component. We adopted the model of the flared disc presented by López-Corredoira & Molgó (2014) in the form:

$$\begin{array}{cc}& {\rho }_{\mathrm{disc}}(R,z)={\rho }_{\mathrm{thin}}(R,z)+{\rho }_{\mathrm{thick}}(R,z),\hfill \\ \hfill & {\rho }_{\mathrm{thin}}(R,z)=(1-f){\rho }_{\odot }\exp (\frac{R_{\odot }}{h_r}+\frac{h_{r,\mathrm{hole}}}{R_{\odot }})\hfill \\ \hfill & \times \exp (-\frac{R}{h_r}-\frac{h_{r,\mathrm{hole}}}{R})\exp (-\frac{|z|}{h_{z,\mathrm{thin}}}),\hfill \\ \hfill & {\rho }_{\mathrm{thick}}(R,z)=f{\rho }_{\odot }\exp (\frac{R_{\odot }}{h_r}+\frac{h_{r,\mathrm{hole}}}{R_{\odot }})\hfill \\ \hfill & \times \exp (-\frac{R}{h_r}-\frac{h_{r,\mathrm{hole}}}{R})\exp (-\frac{|z|}{h_{z,\mathrm{thick}}}),\hfill \\ \hfill \end{array}$$(15)

where the Galactocentric cylindrical coordinate system (R, z) is used. This model takes into account the thick and the thin discs, with an exponential decrease in density in the horizontal and the vertical directions. h_r is the scale length of the whole disc. h_z, thin and h_z, thick are the scale heights of the thin and the thick discs, respectively. The deficit of stars in the inner region of the disc is characterised by the h_r, hole parameter (López-Corredoira et al. 2004). Since we focus on the remote regions of the Galaxy, we kept h_r, hole constant at h_r, hole = 3.74 kpc (López-Corredoira & Molgó 2014). The Galactocentric distance of the Sun is R_⊙ = 8.25 kpc and ρ_⊙ is the volume mass density of the Galactic disc in the solar neighbourhood. The f parameter represents the ratio of thick and thin stars in the solar neighbourhood and we kept it at f = 0.09. We experimented with setting f as a free parameter, but it was varying only slightly, with negligible influence on h_z, and therefore we kept it at the local value as done by López-Corredoira & Molgó (2014). The parameters ρ_⊙, h_r, h_z, thin, and h_z, thick are free and their values were then determined from the fitting procedure. The thin and thick discs were divided based on the geometry of the density profiles in the vertical plane.

6.1. Method

To find the horizontal and vertical star distribution in the Galactic disc, we used the density maps (Fig. 1). We carried the fitting procedure out in two steps: (1) investigating the density profile in the Galactic equatorial plane; (2) investigating the density profile in the vertical direction.

6.2. Density profile in the Galactic equatorial plane

First of all, we focused on the Galactic plane in order to fit the scale length of the Galactic disc. We applied the following constrains on the data: we used stars with |z|< 0.2 kpc and Galactocentric distances R ∈ [5,20] kpc; the bin size was 0.4 kpc in z and 1 kpc in R; and we only considered bins with number of stars N ≥ 50. We applied the disc model with exponential decrease in R in the following form:

$$\begin{array}{cc}& \rho (R)={\rho }_{\odot }\exp (\frac{R_{\odot }}{h_r}+\frac{h_{r,\mathrm{hole}}}{R_{\odot }})\hfill \\ \hfill & \times \exp (-\frac{R}{h_r}-\frac{h_{r,\mathrm{hole}}}{R}),\hfill \end{array}$$(16)

with h_r and ρ_⊙ as free parameters to be fitted. We defined the azimuthal angle ϕ to be measured from the centre-Sun-anticentre direction towards the Galactic rotation, going from 0° to 360°. Since we could not distinguish the individual populations, we neglected the contribution of the thick disc stars in the Galactic equatorial plane. We applied the weighted minimum chi-square method on the data to obtain the values of the fitting parameters in the plane of the Galactic disc given by Eq. (16). The best results of the fitting procedure in the equatorial plane are presented in Table 1, where the comparison with other works is also included. The density profiles in the Galactic plane for the whole data sample are plotted in Fig. 6. Comparing the scale length for various azimuths, one can see the slight dependence of the h_r on ϕ, especially for larger Galactic azimuths; 1.71 ± 0.18 kpc and 1.58 ± 0.15 kpc for Φ ∈ [300°,330°] and for Φ ∈ [30°,60°], respectively. On the other hand, variations of the scale length near the centre-Sun-anticentre direction is not present (2.04 ± 0.05 kpc and 2.26 ± 0.07 kpc for Φ ∈ [330°,0°] and for Φ ∈ [0°,30°], respectively), and the average value of the scale length h_r = 2.19 ± 0.08 kpc (for Φ ∈ [330°,30°]) with little dependence on azimuth. The results are in agreement with previous works. For example, Li et al. (2019) present h_r = 2.10 ± 0.01 kpc for OB stars using Gaia DR2 data, and Chrobáková et al. (2020) show h_r = 2.29 ± 0.08 kpc using Gaia DR2 data. On the other hand, Yu et al. (2021a) show h_r = 1.17 ± 0.05 kpc for OB stars.

Fig. 6. Dependence of the density on the Galactocentric distance in the Galactic equatorial plane for the azimuth ϕ ∈ [330°,30°]. The data points were obtained as weighted mean in bins of size 1 kpc in R and 0.4 kpc in |z|, and were fitted with the model defined in Eq. (16).

6.3. The scale height

In order to investigate the flaring of the Galactic disc, we fitted the vertical density profiles of the whole data sample, which are presented in Fig. 1, with the model of flared thick and thin discs described by Eq. (15). We applied the weighted minimum chi-square method to obtain the values of the scale height of the thick and thin discs, while the scale length calculated in Sect. 6.2 remained fixed. We divided the data in bins with size ΔR = 1 kpc and Δz = 0.2 kpc. The vertical density profiles for various values of Galactocentric distances (R = 11 kpc, R = 14 kpc and R = 18 kpc) are plotted in Fig. 7. The values of the scale height are presented in Table 2 and shown in Fig. 8. The h_z fitting function is a polynomial of the second order. Here, we stress again that we divide the Galactic disc into thin and thick discs, based solely on geometric properties.

Fig. 7. Dependence of the density on |z| for various values of Galactocentric distance. The Galactic azimuth is ϕ ∈ [330°,30°]. The data points were obtained as weighted mean in bins of size 1 kpc in R and 0.2 kpc in |z|.

Fig. 8. Dependence of the scale height of the thick and the thin discs on the Galactocentric distance. The Galactic azimuth is ϕ ∈ [330°,30°]. The dashed line is the second-order polynomial fit to the data points.

The Galactic flare is significant in Sample 0, which represents the whole population. The scale height of the thin disc rises from 0.26 ± 0.01 kpc in the solar neighbourhood (R = 8 kpc) to 0.77 ± 0.08 kpc in the remote regions of the Galactic disc (R = 20 kpc). The flaring of the thick disc is even more significant, where the scale height increases from 0.75 ± 0.01 kpc in R = 8 kpc to reach 3.35 ± 0.50 kpc in R = 19 kpc. The strong flare in the thick disc is present despite large error bars of h_z for R > 16 kpc. Comparison of the results with other works is presented in Fig. 9. The flaring of the thin and the thick discs is presented in the results of all authors, but our data exhibits stronger flaring of the thick disc in the remote regions for R > 14 kpc compared to Yu et al. (2021a) and López-Corredoira & Molgó (2014).

Fig. 9. Comparison of the flare for the whole population (Sample 0) with other works. Our work is represented by the polynomial fits to the scale height data points (for more details and data points with error bars, see Table 2 and Fig. 8).

6.4. The northern and southern flare

We also focused on the differences between the northern and the southern flare. In Fig. 10 we plot the comparison of the northern, southern and northern+southern flares. There is no significant difference in the dependence of the scale height with the Galactocentric distance for R < 15 kpc. However, for larger distances, the flaring of the thick disc is asymmetric. The value of the southern scale height is approxomately 2 kpc higher than the h_z of the northern flare in R = 17 kpc, although the error bars due to the lack of robust datasets in this region have to be taken into account. The difference decreases for R > 17 kpc and the scale height error bars of the northern and southern flares overlap.

Fig. 10. Dependence of the scale height of the thick and the thin discs on the Galactocentric distance. The Galactic azimuth is ϕ ∈ [330°,30°]. The northern, the southern, and the northern+southern flares are compared.

6.5. The azimuthal dependence of the scale height

In order to study the azimuthal dependence of the scale height, we divided the Galactic disc into sectors. We used the following azimuth intervals ϕ ∈ [300°,330°], ϕ ∈ [330°,0°], ϕ ∈ [0°,30°], and ϕ ∈ [30°,60°]. The results for various Galactocentric distances R ∈ [13, 15, 17] kpc are plotted in Fig. 11.

Fig. 11. Dependence of the scale height on the Galactic azimuth ϕ for various Galactocentric distances: R = 13 kpc (red lines); R = 15 kpc (blue lines); R = 17 kpc (green lines). Dotted lines represent the scale height of the thick disc and solid lines represent the scale height of the thin disc. Azimuth is binned with size Δϕ = 30°.

The dependence of the scale height is not significant in our data for R < 17 kpc. The h_z value of the thick disc in R = 17 kpc for the azimuth Φ ∈ [300° ,330° ] is significantly lower than for the rest of azimuth intervals. However, considering the error bars of h_z, the azimuthal asymmetry almost vanishes.

6.6. Supergiants

We also investigated the flaring of supergiants. We only used Sample 2, defined in Sect. 4.2, as Sample 1 was incomplete and therefore could not be used to model the flare. Since we only put constraints on the absolute magnitude of the stars (−5 < M < −10), contamination is expected in the dataset, although as mentioned is Sect. 5, it is not significant. We followed the flare analysis procedure described in Sects. 6.1–6.3. We considered the thin disc and we did not put any constrains on the azimuth. Due to the significantly smaller size of the dataset, we used larger bins: ΔR = 2 kpc and Δz = 0.3 kpc. We calculated the scale length h_r = 1.99 ± 0.13 kpc of the disc. The flaring of Sample 2 stars is presented in Table 3. One can see the subtle increase in h_z, from values h_z ≈ 0.2 kpc in the solar neighbourhood (R ≈ 8 kpc) to h_z ≈ 0.8 kpc in the remote regions of the Milky Way (R ≈ 18 kpc). A significant increase in the scale height appears for R ≥ 13 kpc, where h_z ∈ [0.48, 0.84] kpc. Although the error bars are huge in this interval of the Galactocentric distances (from ±0.2 kpc to ±0.7 kpc), the higher value of the scale height for R ∈ [13, 15] kpc (hz ≈ 0.5 kpc) might be a real feature in the studied data sample. The comparison with other works using OB stars is plotted in Fig. 12.

Fig. 12. Comparison of the thin disc scale heights of the supergiants (Sample 2) with other works.

7. Conclusions

We used Gaia EDR3 to study the outer Galactic disc using supergiants and compared them with the whole population. We concentrated on the Galactic warp and flare. The warp of the whole population is similar to the results from previous works (e.g., Chrobáková et al. 2020), reaching a maximum amplitude of z_w = 0.360 kpc and a minimum amplitude of z_w = −0.375 kpc at a distance R = [19.5, 20] kpc, revealing a small asymmetry between the northern and the southern warp. The warp of the supergiants, which are notably younger than the whole population, reaches a much larger maximum amplitude of z_w = 0.658 kpc and a minimum of z_w = −0.717 kpc at a distance R = [19.5, 20] kpc, with the north-south asymmetry maintained. The difference between the warps of the two populations is significant, confirming a significant relationship between the age of the studied population and the warp amplitude. This result suggests that the warp is induced by a non-gravitational mechanism, such as accretion of intergalactic matter onto the disc or an intergalactic magnetic field.

We find a significant flare of the whole population, especially in the thick disc. The scale height increases from h_z, thick = 0.75 ± 0.01 kpc and h_z, thin = 0.26 ± 0.01 kpc for R = 8 kpc, to h_z, thick = 2.63 ± 0.34 kpc and h_z, thin = 0.64 ± 0.07 kpc at the Galactocentric distance R ≈ 18 kpc, if the azimuth Φ ∈ [330° ,30° ] is considered. We also investigated the dependence of the scale height on the azimuth, which is not present for R < 17 kpc. For R = 17 kpc and Φ ∈ [300° ,330° ],  the changes of the h_z are visible. However, considering the error bars of the h_z, the azimuthal asymmetry almost vanishes. On the other hand, we find a small north-south asymmetry, especially in the thick disc. The asymmetry appears for R > 15 kpc and the value of h_z of the southern flare is approximately 1 kpc higher than for the northern flare. However, the error bars of the h_z for the northern and the southern flares overlap for R > 17 kpc. A subtle flare is present in the supergiants population, in comparison with the whole population of the disc. We find a significant increase in the h_z for R ≥ 13, but considering error bars, this rise of h_z remains inconclusive. The h_z error bars for R ≥ 13 reach higher values due to a significant dispersion of densities in vertical profiles and asymmetry between northern and southern flares. It is therefore clear that the population of supergiants is warped more significantly than the whole population, while the flare for this population is less prominent.

Investigation into the warping and flaring mechanisms (e.g., various mechanisms of disc heating, mergers, magnetic field, etc.) is beyond the scope of this study. However, the forthcoming Gaia data releases promise significant improvement in positional and kinematic data, which will pave the way for a better understanding of the dynamics forming these structural features in the remote regions of the Milky Way.

Acknowledgments

We thank the anonymous referee for helpful comments, which improved this paper. Z.C. and M.L.C. were supported by the grant PGC-2018-102249-B-100 of the Spanish Ministry of Economy and Competitiveness (MINECO). R.N. was supported by the VEGA – the Slovak Grant Agency for Science, grant No. 1/0761/21, by the Slovak Research and Development Agency under the contract No. APVV 18-0103 and by the Erasmus+ programe of the European Union under grant No. 2020-1-CZ01-KA203-078200. 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 IACA13-3E-2493. 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.

References

All Tables

All Figures

	Fig. 1. Density maps in Galactocentric coordinates at various azimuths for the whole population (Sample 0).
In the text

	Fig. 2. Distribution of supergiants in galactic coordinates for |Z|< 4 kpc. Left: Sample 1. Right: Sample 2 (see Sects. 4.1 and 4.2 for details).
In the text

	Fig. 3. Hertzsprung-Russell diagrams (HRD) for the three samples defined in the text. As the datasets are significantly large, in order to avoid the saturation of the diagram, we plot a randomly chosen sub-sample containing 3000 stars in each case. Top: Sample 0. Middle: Sample 1. Bottom: Sample 2.
In the text

	Fig. 4. Comparison of fits of minimum and maximum warp amplitudes for supergiants, chosen by two different approaches (see Sects. 4.1 and 4.2 for details).
In the text

	Fig. 5. Comparison of fits of maximum and minimum warp amplitudes for the whole population (Sample 0) and supergiants (Sample 2).
In the text

	Fig. 6. Dependence of the density on the Galactocentric distance in the Galactic equatorial plane for the azimuth ϕ ∈ [330°,30°]. The data points were obtained as weighted mean in bins of size 1 kpc in R and 0.4 kpc in |z|, and were fitted with the model defined in Eq. (16).
In the text

	Fig. 7. Dependence of the density on |z| for various values of Galactocentric distance. The Galactic azimuth is ϕ ∈ [330°,30°]. The data points were obtained as weighted mean in bins of size 1 kpc in R and 0.2 kpc in |z|.
In the text

	Fig. 8. Dependence of the scale height of the thick and the thin discs on the Galactocentric distance. The Galactic azimuth is ϕ ∈ [330°,30°]. The dashed line is the second-order polynomial fit to the data points.
In the text

	Fig. 9. Comparison of the flare for the whole population (Sample 0) with other works. Our work is represented by the polynomial fits to the scale height data points (for more details and data points with error bars, see Table 2 and Fig. 8).
In the text

	Fig. 10. Dependence of the scale height of the thick and the thin discs on the Galactocentric distance. The Galactic azimuth is ϕ ∈ [330°,30°]. The northern, the southern, and the northern+southern flares are compared.
In the text

	Fig. 11. Dependence of the scale height on the Galactic azimuth ϕ for various Galactocentric distances: R = 13 kpc (red lines); R = 15 kpc (blue lines); R = 17 kpc (green lines). Dotted lines represent the scale height of the thick disc and solid lines represent the scale height of the thin disc. Azimuth is binned with size Δϕ = 30°.
In the text

	Fig. 12. Comparison of the thin disc scale heights of the supergiants (Sample 2) with other works.
In the text