© The Authors 2025

1. Introduction

It is well known that the Milky Way (MW) harbours a stellar bar in the central region (e.g. Liszt & Burton 1980; Binney et al. 1991; Weinberg 1992; Binney et al. 1997; Blitz & Spergel 1991; Hammersley et al. 2000; Wegg & Gerhard 2013). However, even after dedicated efforts in the field, the properties of the MW’s bar still remain poorly constrained. The European Space Agency’s Gaia mission has provided an unprecedented, holistic view of the MW by providing the 6D position-velocity measurements and chemistry of ∼33 million stars in the solar neighbourhood and beyond (Gaia Collaboration 2018, 2023). However, the extinction due to dust has obscured our view towards the Galactic centre, especially at lower latitudes, namely, closer to the Galactic mid-plane (see e.g. Nataf et al. 2013). This, in turn, prevents us from identifying the exact spatial extent of the bar from the stellar density distributions (unlike the external barred galaxies) as well as measuring the properties of the MW’s bar, solely based on the Gaia observations.

The salient properties of the MW’s bar, which are of importance in the field of Galactic dynamics, are the strength, length, and the pattern speed (Ω_bar). The importance of accurately measuring these properties extends beyond the usual notion of comparing MW’s bar with other external barred galaxies. Bars are known to play a pivotal role in driving the secular evolution of disc galaxies (e.g. see Vera et al. 2016, and references therein). Bars play a role in the redistribution of stars and reshaping metallicity distributions via radial migration (e.g. Di Matteo et al. 2013; Kubryk et al. 2013; Halle et al. 2015; Khoperskov et al. 2020b; Haywood et al. 2024), exciting dark gaps along the bar minor axis (e.g. see Kim et al. 2016; Ghosh et al. 2024b), driving vertical breathing (Monari et al. 2015; Khachaturyants et al. 2022) and bending motions (e.g. Khoperskov et al. 2019), as well as exciting ridge-like features in the phase-space (e.g. see Dehnen 2000; Fragkoudi et al. 2019; Trick et al. 2021), impacting the structure of stellar streams in the halo (e.g. see Price-Whelan et al. 2016; Erkal et al. 2017; Bonaca et al. 2020), funnelling gas in the inner region of galaxies thus facilitating in starbursts and formation of nuclear discs (e.g. Shlosman et al. 1990; Sheth et al. 2005), producing large-scale streaming motions in both stars and gas (e.g. see Sellwood & Wilkinson 1993; Athanassoula 1992a,b). Furthermore, a recent study by Ghosh et al. (2023b) showed that the accuracy of recovering the underlying (axisymmetric) potential and the distribution function (DF) in a MW-like barred galaxy critically depends on the location of the survey volume with respect to the bar (also see Khoperskov et al. 2025). Therefore, an accurate measurement of the properties of the bar in the MW is the need of the hour to quantify the bar-driven secular evolution in the MW over cosmic time.

While the jury is still out on the exact values of the bar properties in the MW, multiple efforts have been made in the past towards measuring the pattern speed, length, and orientation of the MW’s bar with respect to the Sun. Earlier bar pattern speed measurements have used techniques ranging from applying the Tremaine-Weinberg method (Tremaine & Weinberg 1984) to the stellar velocity field to matching bar resonance features in the stellar velocity field to matching gas dynamics in presence of a stellar bar. These approaches have favoured a fast-short bar scenario with a somewhat larger value of the pattern speed, ranging from ∼50 to 60 km s⁻¹ kpc⁻¹ (e.g. see Fux 1999; Dehnen 2000; Debattista et al. 2002; Bissantz et al. 2003; Antoja et al. 2014). However, more recent observational measurements have favoured a relatively lower value for the bar pattern speed, mostly converging towards ∼40 km s⁻¹ kpc⁻¹ (see e.g. Sormani et al. 2015; Li et al. 2016; Portail et al. 2017; Bovy et al. 2019; Sanders et al. 2019; Clarke & Gerhard 2022; Li et al. 2022; Lucey et al. 2023). As for the bar length of the MW, an earlier study by Hammersley et al. (1994), using the star counts from the Two-Micron Galactic Survey, estimated the bar length to be ∼4 kpc. In addition, Wegg et al. (2015) used red clump giant (RCG) stars to measure the MW’s bar length to be ∼5 kpc. However, recent study by Lucey et al. (2023), employing a technique based on the maximal extent of trapped bar orbits as an estimate of bar length, determined the bar length to be ∼3.5 kpc. While this latter method relies on robust dynamical arguments, the resulting bar length measurement critically depends on the assumed underlying potential of the Galaxy. Furthermore, the bar length, computed from the stellar density field, can be overestimated if the bar is connected to the spiral structure (e.g. see Hilmi et al. 2020; Ghosh & Di Matteo 2024; Vislosky et al. 2024). As for the orientation of the MW’s bar, Wegg et al. (2015) estimated that the bar is at ∼28° −33° with respect to the Sun. The bar strength and shape are even less well known (however, see the models to compare the strength of the bar resonances in the solar neighbourhood from the Gaia data in Monari et al. 2019).

From the Gaia Data Release 3 (hereafter Gaia DR3), Gaia Collaboration (2023) showed the existence of a ‘quadrupole’ or butterfly-like pattern in the stellar mean radial velocity (⟨V_R⟩) field of red giant branch (RGB) stars (see their Fig. 16). This feature has been reported previously from the APOGEE line-of-sight (LoS) velocities and the Gaia DR2 astrometry (Bovy et al. 2019) and from the APOGEE LoS velocities and the Gaia EDR3 astrometry (Queiroz et al. 2021). Initial theoretical studies, using test particle simulation as well as by using self-consistent N-body simulation of a barred galaxy, showed that a quadrupole feature in the stellar mean radial velocity field (⟨V_R⟩) is excited by a central stellar bar (see e.g. Bovy et al. 2019; Gaia Collaboration 2023). Furthermore, a recent study by Vislosky et al. (2024), carried out by empirically matching the maps of ⟨V_R⟩ from the Gaia DR3 and the barred-spiral galaxies from the TNG50 cosmological simulation, proposed that the MW’s stellar velocity field is consistent with a short bar (with a bar length ∼3.6 kpc) connected to a spiral arm. While it is theoretically understood that a bar can excite a quadrupole feature in the stellar mean radial velocity field, the question it is possible to use the quadrupole feature as a robust kinematic diagnostic to constrain (some of) the properties of the bar is still beyond our grasp. Another directly related question would consider whether we can use the quadrupole feature, as seen from the Gaia data, to put stringent constraints on the bar properties of the MW from a purely dynamical argument. If this is indeed feasible, it could potentially mitigate the existing conundrum of measuring the properties of the bar in the MW from the stellar density field (which is severely affected by the dust attenuation). We aim to pursue this possibility in this work.

In this work, we carried out a systematic study to test the reliability and robustness of the quadrupole pattern in the stellar mean radial velocity field (⟨V_R⟩) as a kinematic diagnostic to put constraints on the bar properties. To this end, we made use of a suite of collisionless N-body models (having both thin and the thick discs), which forms a prominent bar and a boxy/peanut bulge. In addition, we used a sample of barred galaxies, selected from the TNG50 cosmological simulation to augment this study. Within the scope of this paper, we quantify the properties of the quadrupole feature and then we systematically investigate how robustly the properties of the quadrupole trace the properties of the bar. In addition, we investigate how the biases and the errors in measuring parallax, proper motion, and radial velocity of stars, similar to the Gaia DR3, could influence the measured properties of the bar as inferred from the properties of the quadrupole feature. Lastly, we measure the properties of the quadrupole feature in the MW, while using the full 6D phase-space information from the Gaia DR3.

The rest of the paper is organised as follows. Section 2 provides a brief description of the isolated N-body simulations as well the barred galaxies from the TNG50 suite of cosmological simulation, used for this study. Section 3 provides the details of the quantification of the properties of the quadrupole feature as well as their correlation with the properties of the bar. Section 4 contains the details of the influence of Gaia-like uncertainties on inferring the bar properties from the properties of the quadrupole. It also presents results pertaining to the quantification of quadrupole’s properties from Gaia DR3. Section 5 summarises the main findings of this work.

2. Simulated barred galaxies

Here, we briefly describe the initial equilibrium configurations and the structural properties of the suite of the isolated, collisionless N-body models as well as the sample of the barred galaxies, chosen from the TNG50 cosmological simulation.

2.1. Isolated barred models

A total of 14 isolated, collisionless N-body models are used for this work. Below, we briefly mention their structural properties and the initial equilibrium set-up.

2.1.1. Thin+thick models:

A total of 13 thin+thick models (with varying disc geometry and thick-disc mass fraction), taken from Ghosh et al. (2023a), are considered in this work. The initial equilibrium configuration of each of these models consists of a thin and a thick stellar disc, embedded in a live dark matter halo.

Each of the thin and thick discs was modelled with a Miyamoto-Nagai profile (Miyamoto & Nagai 1975), having R_d, z_d, and M_d as the characteristic disc scale length, the scale height, and the total mass of the disc, respectively. The scale heights of the thick and thin discs were fixed to 0.9 kpc and 0.3 kpc, respectively. The total stellar mass is fixed to to 1 × 10¹¹ M_⊙ across the suite of simulations, while f_thick (denoting the mass fraction in the thick disc) was varied from 0 to 0.7 in different models. The dark matter halo was modelled by a Plummer sphere (Plummer 1911), with R_H( = 10 kpc) and M_dm( = 1.6 × 10¹¹ M_⊙) as the characteristic scale length and the total halo mass, respectively. The values of the key structural parameters for the thin and thick discs are noted in Table 1. A total of 1 × 10⁶ particles were used to model the stellar (thin+thick) disc, while a total of 5 × 10⁵ particles were used to model the dark matter halo.

Following the iterative method algorithm by Rodionov et al. (2009), the initial conditions of the discs were obtained, while keeping the velocity dispersion (along the radial and vertical directions) fixed and letting the density to vary until the desired equilibrium confirmation is achieved. The simulations are run using a TreeSPH code by Semelin & Combes (2002), which employs a hierarchical tree method (Barnes & Hut 1986) with an opening angle of θ = 0.7 to compute the gravitational forces. In addition, a Plummer potential was employed for softening the gravitational forces with a softening length ϵ = 150 pc. We evolved all the models for a total time of 9 Gyr. For further details, we refer to Fragkoudi et al. (2017) and Ghosh et al. (2023a). Each of the models forms a prominent stellar bar which subsequently undergoes a vertical buckling instability to form a boxy/peanut (hereafter b/p) bulge (for details, see Ghosh et al. 2023a, 2024a). We mention that in rthickE models, R_d, thick = R_d, thin; in rthickS models, R_d, thick < R_d, thin; and in rthickG models, R_d, thick > R_d, thin, where R_d, thin and R_d, thick denote the scale length for the thin and thick disc, respectively. Furthermore, following the convention used in Ghosh et al. (2023a), any thin+thick model is referred as a unique string ‘[model configuration][thick disc fraction]’ where [model configuration] denotes the corresponding thin-to-thick disc scale length configuration while [thick disc fraction] denotes the value of f_thick.

2.1.2. sim6 model:

This is a higher resolution, collisionless N-body model (i.e. as compared to other thin+thick models used here), which also forms a prominent bar and subsequently undergoes a vertical buckling instability to form a prominent b/p structure. In addition, this model has been extensively used in studying the bar-spiral driven chemo-dynamical evolution of MW-like galaxies (for details, see Fragkoudi et al. 2018, 2019; Khoperskov et al. 2020a,b). The initial equilibrium configuration consists of a thin disc, an intermediate disc, and a thick disc (each modelled with a Miyamoto-Nagai profile with R_d, z_d, and M_d as the characteristic disc scale length, the scale height, and the total mass of the disc) and the stellar discs are embedded in a concentric live dark matter halo. The dark matter halo was modelled by a Plummer sphere, with R_H( = 21 kpc) and M_dm( = 3.7 × 10¹¹ M_⊙) as the characteristic scale length and the total halo mass, respectively. The scale lengths of the thin, the intermediate, and the thick disc were set to 4.8 kpc, 2 kpc, and 2 kpc, respectively. The total stellar mass (thin+intermediate+thick) of the model was fixed to ∼8.7 × 10¹⁰ M_⊙, where the thin disc contributes 50 per cent of the total stellar mass, while the intermediate and the thick disc constitutes 30 per cent and 20 per cent of the total stellar mass, respectively. The scale heights of the thin, intermediate, and thick disc are set to 0.15 kpc, 0.3 kpc, and 0.6 kpc, respectively. For further details, we refer to Fragkoudi et al. (2019).

The initial equilibrium configuration was achieved using the same iterative method algorithm applied in the work by Rodionov et al. (2009). The simulations were run using a parallel MPI tree-code (Khoperskov et al. 2014), which takes into account the adaptive spatial decomposition of particle space between nodes, along with an opening angle of θ = 0.7 to compute the gravitational forces. A Plummer potential was employed for softening the gravitational forces with a softening length of ϵ = 50 pc. A total of 1.5 × 10⁷ particles were used in the model with 1 × 10⁷ used for the disc and 5 × 10⁶ for the dark matter halo. The model was evolved for a total time of 7 Gyr (for further details, see Fragkoudi et al. 2019).

2.2. Barred galaxies from the TNG50 simulations

In this work, we also analysed a sample of the MW and M31 analogues (Pillepich et al. 2024) selected from the TNG50 simulation (Nelson et al. 2019b,a; Pillepich et al. 2019)¹. TNG50 is a magneto-hydrodynamical simulation of the formation and evolution of galaxies in a 51.7 comoving Mpc cube from redshift ≈127 to redshift 0. It was run with the moving-mesh code AREPO (Springel 2010) and uses the fiducial TNG galaxy formation model (Weinberger et al. 2017; Pillepich et al. 2018), with a mass resolution of m_baryon = 8.5 × 10⁴ M_⊙, m_DM = 4.5 × 10⁵ M_⊙; along with a spatial resolution of star-forming gas of ∼150 pc (Pillepich et al. 2024).

The selection criteria of MW/M31 analogues at z = 0 include the following: (i) the galaxy stellar mass is in the following range: M_*(< 30 kpc) = 10^{10.5 − 11.2} M_⊙; (ii) a disc-like stellar morphology; (iii) no other galaxy with stellar mass > 10^10.5 M_⊙ is within 500 kpc distance; and (iv) the total mass of the halo host is smaller than that typical of massive groups < 10¹³ M_⊙ (for further details, see Pillepich et al. 2024). Since we are interested in barred galaxies here, from the parental sample of 198 MW and M31 analogues, we selected 55 galaxies with the most prominent bars.

3. Quantifying bar properties via quadrupole signature in the mean radial velocity

Figure 1 shows one example of the distribution of the stellar density, mean radial velocity (⟨V_R⟩), and the radial velocity dispersion (σ_R) in the (x − y)-plane, calculated at t = 9 Gyr for the rthick0.0 model. We used a galactocentric cylindrical coordinate system (R, ϕ, z) to compute the distribution of ⟨V_R⟩ and σ_R. The model harbours a prominent stellar bar in the central region and the associated ⟨V_R⟩ map displays a prominent quadrupole feature. The model gets heated preferentially along the 2D extent of the bar, similar to the findings of Ghosh et al. (2023a). In Appendix A, we show the distribution of ⟨V_R⟩ in the (x − y)-plane, calculated at the end of the simulation (t = 9 Gyr), for all the thin+thick models considered here (see Fig. A.1 there). In addition, we show the corresponding ⟨V_R⟩ maps in the (x − y)-plane for a sample of TNG50 barred galaxies (with varying bar morphology; see Fig. A.2 in Appendix A). A prominent bar is always associated with a clear quadrupole feature in the ⟨V_R⟩ map, thereby demonstrating that the quadrupole feature in the ⟨V_R⟩ map is a (kinematic) part and parcel of the bar in the density distribution.

Fig. 1. Tracing the bar with the quadrupole feature: distribution of stellar surface density (left panel), mean radial velocity, ⟨VR⟩ (middle panel), and radial velocity dispersion, σR (right panel), calculated in the (x − y)-plane, for the rthick0.0 model, at the end of the simulation run (t = 9 Gyr). Dashed black lines denote the contours of constant surface density. The cyan dashed circle denotes the bar length, Rbar, and the magenta dashed circle denotes the extent of the quadrupole feature, Rquadrupole. The magenta points denote the spatial distribution of the phase-angle of the m = 4 Fourier moment (φ4). The bar excites a prominent quadrupole pattern in the mean radial velocity field, and the orientation of the quadrupole pattern agrees fairly accurately with the orientation of the bar. We used a galactocentric cylindrical coordinate system (R, ϕ, z) to calculate ⟨VR⟩ and σR.

Next, we quantified the strength, length, and the orientation of the quadrupole feature present in the face-on ⟨V_R⟩ maps. This was achieved by means of the Fourier decomposition of the mean radial velocity field. Subsequently, we studied how they evolve over time for the isolated bar models as well the TNG50 galaxies considered here. Sect. 3.1 provides the details of the temporal evolution of the strength and length of the quadrupole features as well as correlation with the corresponding properties of the bar, while Sect. 3.2 provides the details of the orientation of the quadrupole feature and their linkage with the bar orientation.

3.1. Correlation between strength and length of the bar and the quadrupole feature

The strength and the length of the bar in our models were computed from the m = 2 Fourier coefficients of the underlying density distribution, using

$$\begin{array}{c}\hfill A_2/A_0(R)=\frac{{\sum }_j{m}_j{e}^{i2{\varphi }_j}}{{\sum }_j{m}_j}.\end{array}$$(1)

Here, A₂ denotes the coefficient of the m = 2 Fourier moment of the density distribution, m_j is the mass of the jth particle, and ϕ_j is its azimuthal angle. The summation runs over all the particles within the radial annulus [R, R + ΔR], with ΔR = 0.5 kpc. For each timestep (or equivalently, each snapshot), we calculated the radial profiles of the m = 2 Fourier coefficient from R = 0 to R = 6R_d, thin. At a time, t, we define the strength of the bar, S_bar, as the peak value of the m = 2 Fourier coefficient (A₂/A₀). In addition, at time, t, we define the bar length, R_bar as the radial location where A₂/A₀ drops to 70 per cent of its peak value in the central bar region. For a detailed exposition to different methods of measuring the length of a bar, the reader is referred to a recent study by Ghosh & Di Matteo (2024). The values of S_bar and R_bar, at different times, for all the thin+thick models are taken from Ghosh et al. (2023a), while for the sim6 model and the sample of TNG50 galaxies, we computed the corresponding values of S_bar and R_bar using the method described above.

The strength and the length of the quadrupole feature in our models are computed from the m = 4 Fourier coefficients of the underlying mean radial velocity (⟨V_R⟩) field using

$$\begin{array}{c}\hfill A_4/A_0(⟨V_{\mathrm{R}}⟩)=\frac{{\sum }_j{m}_j|V_{\mathrm{R},j}|e^{i4{\varphi }_j}}{{\sum }_j{m}_j|V_{\mathrm{R},j}|},\end{array}$$(2)

where |V_R, j| denotes the absolute value of the mean radial velocity of the jth particle. We used the same radial binning, as used in Eq. (1). Furthermore, the quantity A₄/A₀(⟨V_R⟩) is calculated from R = 0 to R = 6R_d, thin. The corresponding phase angle, φ₄ can be obtained from Eq. (2), using φ₄ = 1/4tan⁻¹(B₄/A₄), where A₄ and B₄ denotes the real and imaginary parts of the summation in Eq. (2). At time, t, we define the strength of the quadrupole, S_quadrupole as the peak value of the m = 4 Fourier coefficient of the underlying stellar mean radial velocity field (A₄/A₀(⟨V_R⟩)). The corresponding temporal evolution of the S_quadrupole for all thin+thick models as well as for the sim6 model is shown in Appendix A (see top panels of Fig. A.3 there). In addition, at time, t, we define the length of the quadrupole, R_quadrupole, as the radial location where A₄/A₀(⟨V_R⟩) drops to the 70 per cent of its peak value in the central bar region. Using the same definition, we also computed the corresponding strength and length of the quadrupole feature for the sample of TNG50 barred galaxies considered here. Furthermore, we recognise that different isolated thin+thick models and the TNG50 barred galaxies used here have different disc scale lengths. Therefore, to carry out a uniform comparison, we need to normalise the values of R_quadrupole by the corresponding disc scale length. The same argument applies to the bar length (R_bar) as well. For the thin+thick models and the sim6 model, we measure the (average) disc scale length ⟨R_d⟩ using (Ghosh et al. 2023a):

$$\begin{array}{c}\hfill ⟨R_{\mathrm{d}}⟩=\frac{M_{\mathrm{d},\mathrm{thin}}{R}_{\mathrm{d},\mathrm{thin}}+M_{\mathrm{d},\mathrm{thick}}{R}_{\mathrm{d},\mathrm{thick}}}{M_{\mathrm{d},\mathrm{thin}}+M_{\mathrm{d},\mathrm{thick}}},\end{array}$$(3)

where M_d, j denotes the stellar mass and R_d, j denotes the scale length of the jth component (j= thin, thick). As for the TNG50 barred galaxies, we fit a single exponential profile of the form Σ(R)∝Σ₀exp[−R/R_d] to the surface density profiles along the bar major axis, while excluding the central bar region².

Next, we investigated whether there is any correlation between the strength and length of the quadrupole feature and the bar. Fig. 2 (left panel) shows the corresponding correlation between the bar strength, S_bar, and the quadrupole strength, S_quadrupole, for all isolated models and the TNG50 galaxies considered here. We mention that in different thin+thick models, bar forms at different times (for further details, see Ghosh et al. 2023a, 2024a). However, we checked that a prominent bar is always associated with a quadrupole feature, regardless of its formation time. Therefore, only the snapshots after the bar forms are considered for all the thin+thick as well as for the sim6 model. In addition, to quantify whether the quantities S_bar and S_quadrupole are correlated or not, we compute the corresponding Pearson correlation coefficient, ρ, and find that indeed these two quantities are strongly correlated (ρ ≥ 0.75). In other words, the bar and the quadrupole feature in the stellar velocity field evolve in tandem, and this trend holds for all the models considered here. Furthermore, we fit a straight line of the form Y = AX + B to all the points in the S_bar-S_quadrupole plane to check if these two quantities follow a linear scaling law. The best-fit parameters (A = 0.96 ± 0.1; B = −0.18 ± 0.001) suggest that the strength of quadrupole are indeed linearly related to the bar strength. Fig. 2 (right panel) shows the corresponding correlation between the bar length, R_bar and the quadrupole length, R_quadrupole for all isolated models and the TNG50 galaxies considered here. The calculated Pearson correlation coefficient, ρ is found to be greater than 0.75; thereby demonstrating that these two quantities are also strongly correlated. We fit a straight line (of the form Y = AX + B) to all the points in the R_bar − R_quadrupole plane (both quantities being normalised by the same average disc scale length, ⟨R_d⟩). As the best-fit parameters (A = 1.06 ± 0.1; B = 0.01 ± 0.01) suggest, the lengths of the bar and the quadrupole remain linearly related to each other.

Fig. 2. Tracing the bar properties with the quadrupole feature: correlation between the bar strength, Sbar and the strength of the quadrupole, Squadrupole (left panel), and correlation between the bar length, Rbar and the length of the quadrupole, Rquadrupole (right panel), computed using all isolated thin+thick models and the TNG50 barred galaxies (see the legend). The black dash line denotes the best-fit straight line (of the form Y = AX + B) while the grey shaded region denotes the 5-σ scatter around the best-fit line. The properties of the bar (strength and extent) remain strongly correlated with the properties of the quadrupole structure (Pearson correlation coefficient, ρ > 0.75).

Lastly, we investigated whether the choice of kinematically cold thin-disc or kinematically hot thick-disc stars can affect the earlier found correlations between the properties of the bar and the quadrupole feature. We mention that in the sim6 and all the thin+thick models, we can identify and separate, by construction, which stars are members of the thin-disc component at initial time (t = 0) and which stars are members of the thick-disc component at t = 0; in addition, we can track them as the system evolves self-consistently. This, in turn, allows us to test the dependence (if any) of the correlations found between the properties of the bar and the quadrupole feature on the choice of thin-disc or thick-disc stars. To achieve this, we recalculated the strength and the length of both the bar and the quadrupole feature, using the thin- and thick-disc particles separately. The corresponding correlations, as a function of stellar populations with varied degree of velocity dispersion (i.e. kinematically colder thin disc and kinematically hotter thick disc)³ are shown in Fig. 3. As seen clearly from Fig. 3, both the length and the strength of bar and the quadrupole feature remain strongly correlated (Pearson correlation coefficient, ρ > 0.75), regardless of the choice of thin-disc or thick-disc stars. As for the linear scaling relation, the choice of thin-disc or thick-disc stars does not change appreciably (less than 10 per cent) the value of the best-fit slope of straight line when fitted to R_bar and R_quadrupole (see bottom panels in Fig. 3). However, for the quantities S_bar and S_quadrupole, the best-fit straight line is shallower when computed using only the thick-disc stars, as compared to when they were calculated using only the thin-disc stars (see the corresponding best-fit values in top panels of Fig. 3). As shown in Ghosh et al. (2023a), the thin-disc stars constitute a stronger bar and the thin-disc stars (being the kinematically colder component) are more perturbed kinematically by the non-axisymmetric structures (e.g. m = 2 bar here, also see Debattista et al. (2017)). Therefore, a steeper best-fit straight line for the thin-disc stars in the S_bar − S_quadrupole plane is likely to be a combined result of the two above-mentioned dynamical effects.

Fig. 3. Dependence on the choice of the thin- or thick-disc stars: correlation between the bar strength, Sbar and the strength of the quadrupole, Squadrupole (top row), and correlation between the bar length, Rbar and the length of the quadrupole, Rquadrupole (bottom row), computed using thin-disc particles (left panels) and thick-disc particles (right panels), for all thin+thick models and the sim6 model (see the legend). The black dash line denotes the best-fit straight line (of the form Y = AX + B), while the grey shaded region denotes the 5-σ scatter around the best-fit line. Regardless of the choice of thin or thick-disc stellar particles, the length and strength of the quadrupole remain strongly correlated with the length and the strength of the bar.

3.2. Orientation of quadrupole and the bar

In the earlier section, we demonstrated that the strength and the extent of the quadrupole feature in the ⟨V_R⟩ distribution is an excellent proxy for the strength and the length of a bar. We note that the Fourier decomposition of the mean radial velocity distribution (see Eq. (2)) not only provides the amplitude of the m = 4 quadrupole feature, but the corresponding phase angle (φ₄) also provides the information about the orientation of the quadrupole feature. In Fig. 1, we show the spatial variation in the phase angle of the m = 4 Fourier moment (φ₄) of the ⟨V_R⟩ distribution, within and beyond the bar region. Even a basic visual inspection reveals that the phase angle, φ₄, remains constant within the bar region and the distribution of the φ₄ follows the orientation of the bar. We checked that this trends hold true for all models considered here; thereby implying that the orientation of the quadrupole feature in the ⟨V_R⟩ can be used as a proxy for the bar orientation in the density field. However, in the MW, we almost do not see the other side of the Galactic centre; thus, it is difficult to reconstruct the full quadrupole feature in the corresponding mean radial velocity field. The question remains: in this case (and for arbitrary orientations of the bar), how robustly the orientation of the quadrupole feature (in ⟨V_R⟩ distribution) can trace or recover the bar orientation (in the stellar density field). Here, we test this approach in detail.

Figure 4 shows an example of how the distribution of ⟨V_R⟩ in the (x − y)-plane would appear to a hypothetical observer at a solar-like position (R = −8 kpc, ϕ = 0^o, z = 0) when the bar is placed at different orientations. In addition, to mimick a MW-like scenario, we have only considered the negative-half (x ≤ 0) of the distribution of ⟨V_R⟩ in the (x − y)-plane. We then repeated the Fourier decomposition of the distribution of ⟨V_R⟩ (taking only the negative half and bar placed at different orientations) to recompute the strength and the extent of the quadrupole. In Fig. 5, we show the corresponding strength and orientation of the quadrupole, while placing the bar at different angles. As seen from Fig. 5, the quadrupole strength is recovered within 10 per cent (relative) errors for all assumed bar orientations and considering the negative-half (x ≤ 0). The quadrupole orientation (φ₄) also recovers the bar orientations, albeit with lesser accuracy. Furthermore, the φ₄ values in the central regions, show a departure from the ‘true’ values (however, they remain within ±4°). We notice that for larger bar orientation angles (ϕ_bar > 10°), the ⟨V_R⟩ distribution no longer remains bi-symmetric (w.r.t. x = 0 line) as the part of the lobe from the other-half (i.e. x > 0) starts appearing to the hypothetical observer at a solar-like position (compare cases for ϕ_bar = 0° and ϕ_bar = 40° in Fig. 4). This, in turn, introduces a systematic fluctuation (around the true value of ϕ_bar) in inferring the bar orientation angle from the measured φ₄ values (see the bottom panel of Fig. 5).

Fig. 4. Distribution of the mean radial velocity (⟨VR⟩) in the (x − y)-plane, calculated at t = 9 Gyr for the model rthick0.0, with the bar placed at different viewing angles with respect to a hypothetical observer (shown in diamond) at a solar-like position (R = −8 kpc, ϕ = 0°, z = 0). The cyan circles in each sub-panel denote the variation in the phase-angle of the m = 4 Fourier coefficient (φ4), computed from the distribution of ⟨VR⟩. The dashed straight line denotes the ‘true’ bar angle, computed from the intrinsic density distribution of the stellar particles.

Fig. 5. Top panel: radial variation in the amplitude of the m = 4 Fourier moment of the mean radial velocity, ⟨VR⟩, calculated at t = 9 Gyr for the model rthick0.0, while putting the bar at different orientations (see the legend). Bottom panel: radial variation in the corresponding phase-angle (φ4) of m = 4 Fourier moment of the mean radial velocity, ⟨VR⟩. The points denoting the radial variation computed from the mean radial velocity, ⟨VR⟩, while the horizontal lines denote the corresponding true bar orientation, and the shaded region around the true value denotes a 8° scatter. ‘Half-extent’ refers to the scenario where stars falling only in the negative half (x ≤ 0) are considered while computing the Fourier moments. For further details, see the text. The vertical dash-dotted black line denotes the length of the bar (Rbar = 5.75 kpc).

3.3. Influence of spiral arms in measuring the quadrupole length

In previous sections, we demonstrated that the properties of the quadrupole feature are aptly correlated with the properties of bar. However, the MW harbours other non-axisymmetric features, for example, spirals (e.g. Oort et al. 1958; Georgelin & Georgelin 1976; Gerhard 2002; Churchwell et al. 2009; Reid et al. 2014), which also excite non-zero mean radial velocities in the disc region (e.g. see Siebert et al. 2011, 2012). Interestingly, the non-zero mean radial velocity, excited by spirals, can often overlap or is connected spatially with the non-zero mean radial velocity excited by the bar (i.e. the quadrupole feature, e.g. see Vislosky et al. 2024). This, in turn, can pose a problem in disentangling the dynamical effect of the bar on mean radial velocities, and can result in overestimating the extent of the quadrupole feature. This also resembles the dynamical situation where spirals emerge from the end tip of the bar, causing an overestimation of the bar length (for a detailed discussion, see Hilmi et al. 2020; Ghosh & Di Matteo 2024). In such a dynamical situation, using the constancy of the phase-angle (φ₂) of the m = 2 Fourier moment (of the density distribution) can potentially decrease the overestimation of the bar length due to the presence of spirals, as demonstrated in Ghosh & Di Matteo (2024). Here, we pursue a similar strategy.

In Fig. 6, we show the distribution of ⟨V_R⟩ in the (x − y)-plane at different times for the model rthick0.0 which harbours a bar+spiral feature. As seen clearly from a visual inspection, non-zero ⟨V_R⟩ in the disc region dominated by the spirals, are connected to the quadrupole feature excited by the bar. Figure 6 further shows the distribution the m = 4 Fourier phase-angle (from the ⟨V_R⟩ distribution) in the disc regime where the m = 2 bar is dominant, as well as in the outer disc region where the spirals are dominant. In all such bar+spiral scenarios, the phase-angle, φ₄, remains constant (within 5 − 8°) in the central bar region (also see Fig. 1); meanwhile in the outer disc region (dominated by the spirals), the phase-angle, φ₄ does not remain constant. In other words, phase-angle, φ₄ displays a characteristic modulation as one moves out from the central bar-dominated region to outer spiral dominated region. Ghosh & Di Matteo (2024) showed a similar characteristic change in the m = 2 phase-angle (from the density distribution) for a bar+spiral asymmetry. Therefore, we adopted a new (conservative) definition of the extent of the quadrupole, R_{quadrupole, φ4} as the radial extent within which the m = 4 Fourier phase-angle remains constant (within ∼5 − 8°). Fig. 7 shows the corresponding comparison between the two different length estimators for the quadrupole feature, namely, R_quadrupole and R_{quadrupole, φ4}, at different times for the model rthick0.0. As seen from Fig. 7, the values of R_{quadrupole, φ4} remains systematically lower than the values of R_quadrupole when the model harbours both a stellar bar and (transient) spirals. It is only towards the end phase of the evolution that when the model no longer hosts spirals, these two values match fairly well (see Fig. 7). In other words, the presence of (transient) spirals systematically overestimates the length of the quadrupole feature (with a median relative difference of ∼25 per cent). We checked that this trend holds true for other models, including those displaying a bar+spirals scenario (for the sake of brevity, the models are not shown here). Therefore, our results outline the importance of using R_{quadrupole, φ4} as a more appropriate estimator for measuring the length of the quadrupole (and, in turn, inferring the bar length) in a bar+spirals dynamical scenario.

Fig. 6. Distribution of the mean radial velocity, ⟨VR⟩ in the (x − y)-plane, computed at different times when the model rthick0.0 harbours a strong bar and prominent spirals. Dashed black lines denote the contours of constant surface density. The magenta points denote the spatial distribution of the phase-angle of the m = 4 Fourier moment (φ4). Within the bar region, φ4 remains constant whereas in presence of prominent spirals, φ4 shows a characteristic modulation and this trend holds for all bar+spiral scenarios shown here.

Fig. 7. Correlation between the two methods of measuring the quadrupole extent: Rquadrupole, measured solely from the amplitude of the m = 4 Fourier coefficient of the mean radial velocity ⟨VR⟩, and Rquadrupole, φ4, measured solely from the constancy of the phase-angle, φ4 of the mean radial velocity ⟨VR⟩, as a function of bar age (tbar-age, see the colour bar). Both the quantities, Rquadrupole and Rquadrupoleφ4, are normalised by the mean disc scale length, ⟨Rd⟩. For details, see the text. The dotted straight line denotes the 1:1 relation. The colour bar denotes the age of the bar.

4. Implication on constraining the MW’s bar properties from the quadrupole feature

In previous sections, we have demonstrated that for simulated bars (from a wide variety of numerical simulations), the quadrupole feature in the stellar mean radial velocity distribution serves as an excellent kinematic diagnostic for constraining the bar’s properties. However, we still need to test whether the quadrupole feature still works as a kinematic diagnostic in presence of dust extinction and Gaia-like uncertainties and systematic errors for the MW. We pursue this investigation below. In Sect. 4.1, we describe the sample selection from the Gaia DR3, whereas in Sect. 4.2, we provide the details of creating a Gaia-like mock dataset from the high-resolution sim6 model (for details, see Sect. 2.1). Lastly, our findings are provided in Sect. 4.3.

4.1. Sample selection from Gaia DR3

We selected a sample of RGB stars from the Gaia DR3, broadly following the same procedure as described in Gaia Collaboration (2023), by selecting stars with 3000 K < T_eff < 5500 K and log g < 3.0. We limited ourselves to stars with full 6D kinematics and apparent magnitudes within the brightness limits 5 < G < 15.5 that have high-fidelity astrometric solutions (fidelity_v2 > 0.5 from Rybizki et al. (2021)). The parallax zero points were corrected using the method described in Groenewegen (2021). The sample was augmented with extinction E using the all-sky extinction map from Zhang & Green (2025). The extinction map extends out to 5 kpc (from the solar location), beyond which it is extrapolated using a double-exponential model for the differential dust density. We further applied quality cuts radial_velocity_error <  20 km s⁻¹ and ϖ/σ_ϖ > 5. We calculated the distance using the simple inverse of the measured parallax, d = 1/ϖ. The positions and velocities of the stars are transformed to Galactocentric coordinates using the astropy.coordinates module (Astropy Collaboration 2022). For the transformation, we used z_⊙ = 20.8 pc (Bennett & Bovy 2019), R_⊙ = 8.277 kpc, and v_⊙ = (9.3 km s⁻¹, 251.5 km s⁻¹, 8.59 km s⁻¹) (GRAVITY Collaboration 2022). Thus, we have the Sun located at ( − R_⊙, 0, z_⊙) moving with velocity v_⊙. Our final selected sample contains 5 667 443 stars. The corresponding stellar density and the ⟨V_R⟩ distribution in the Galactocentric (x − y)-plane, computed for our final selected sample from the Gaia DR3 are shown in Fig. 8 (see the right column). A prominent quadrupole feature is present in the ⟨V_R⟩ distribution calculated in the (x − y)-plane, in agreement with Gaia Collaboration (2023).

Fig. 8. Comparison between the Gaia-like mock data and the Gaia DR3 sample. Top panels: ⟨VR⟩ distribution and the bottom panels showing the stellar density, both in the (x − y)-plane, calculated for all the stellar particles from model sim6 (left column), the Gaia-like mock dataset from model sim6 (middle column) and the selected sample of RGB stars from the Gaia Data Release 3. As shown in this work, the mock dataset qualitatively reproduces the behaviour seen in the Gaia sample, notably, observational effects increase the extent of the quadrupole and reduce the angle of its main axis compared to y = 0. The green lines in the left column corresponds to the major and minor axis of the bar, angled at 28°. The black filled circle in the right column denotes the assumed solar location (for further details, see Sect. 4.1).

4.2. Building Gaia-like mock datasets

Next, to investigate any plausible effect of the Gaia-like uncertainties and systematic errors on the robustness of the quadrupole feature as a kinematic diagnostic, we first generated Gaia-like mock datasets, using the snapshots from the sim6 model at different times (and thus for different bar strengths; see Fig. A.3). For each snapshot, we placed a hypothetical observer at a solar-like position at ( − R_⊙, 0, z_⊙) and at a bar viewing angle of ϕ_bar = 28° to imitate the qualitative behaviour of a sample of RGB stars in Gaia DR3.

Next, we modelled the errors in the measured parallax, radial velocity, proper motion, and apparent magnitude. The observational errors for each simulated star were taken from randomly chosen stars in the Gaia DR3 sample with apparent magnitude G > 15. The apparent magnitude threshold is chosen because most of the stars in Gaia DR3 for which the effect of observational errors is significant are further away and have their apparent magnitudes in the corresponding range. The observed values of parallax, radial velocity, and proper motion are then simulated by adding normally distributed errors to the true values, with the standard deviations as reported above: X_observed ∼ 𝒩(μ = X_true, σ² = X_error²) for an observable X. The above process is repeated ten times, in order to end up with a comparable number of particles as is in the Gaia dataset.

To account for the selection function, we adopted the same procedure described in Zhang et al. (2024). This was done by matching each real star in the Gaia DR3 dataset with its closest counterpart in the simulation. Simulated particles that failed to obtain a match or that were more than 0.1 kpc away from the matched real star are discarded. The final mock dataset has a spatial distribution that closely mimics that of the Gaia DR3 sample (see Fig. B.1 in Appendix B). This approach has the benefit of implicitly accounting for complicated selection effects such as crowding, dust, and the Gaia scanning law. However, it can fail to properly account for the origin of the particles. For example, in Gaia, a star hidden behind a thick dust cloud will never be observable due to the dust cloud blocking the light, but in the mock dataset, it could radially scatter into the visible volume. The ⟨V_R⟩ distribution for one of the mocks, calculated at t = 3 Gyr for the sim6 model, is shown in Fig. 8 (see middle column).

4.3. Quadrupoles in Gaia-like mocks: results and limitations

First, we chose the snapshot at t = 3 Gyr from the sim6 model to investigate the properties of the quadrupole in ⟨V_R⟩ in Gaia-like mock datasets. We note that there is nothing special about the choice of t = 3 Gyr, except the fact that it harbours a prominent stellar bar in the central region. Later, we have tested the same robustness of quadrupole in ⟨V_R⟩ for other snapshots as well. Looking at the distribution of ⟨V_R⟩ in the (x − y)-plane for the MW RGB subsample and the mock dataset highlighted in Fig. 8, we observe largely the same behaviour that would be exhibited by a prominent and extended quadrupole signature near the galactic centre, which tapers off around R = 4 kpc. For the simulated model, even a visual inspection reveals that the bar angle in the mock data is lying closer to the y = 0 axis (as derived from the orientation of the quadrupole as φ₄ ∼ 5°); this is in sharp contrast with the true bar angle of ∼28°. In addition, using the methods described in Sect. 3.1, we calculated the strength of the quadrupole from the mock dataset to be S_quadrupole = 0.37 and the length, R_quadrupole = 3.8 kpc. When compared with the ‘true’ values (i.e. directly calculated from the simulation at t = 3 Gyr), we find that the S_quadrupole is overestimated by ∼40 per cent and the R_quadrupole is overestimated by ∼35 per cent. We repeated this procedure for a few other snapshots that harbour prominent bars. We find that the overestimation obtained in inferring S_quadrupole varies in the range 40 − 45 per cent (of the ‘true’ value). For R_quadrupole, the corresponding overestimation varies in the range 35 − 45 per cent (of the ‘true’ value). Therefore, the values of the quadrupole’s properties are not well recovered from the Gaia-like mock data. The large discrepancy between the values derived from the Gaia-like mock dataset and the true values (directly obtained from the simulation) is quite puzzling and warrants further investigation. We pursue this issue next.

In Fig. 9, we demonstrate how the distribution of ⟨V_R⟩ in the (x − y)-plane, changes in the mock dataset as we sequentially incorporate the uncertainties from proper motion, radial velocity, and parallax measurements. As seen from Fig. 9 (top-middle panel), when the mock data does not include any uncertainties, but is completeness-limited, it still shows a quadrupole feature (albeit fainter). The weakening of the quadrupole is due to the effect of dust (implicit in the applied magnitude cut), which essentially blocks the stars towards the Galactic centre and closer to the mid-plane. Since the effect of bar decreases as one moves away from the mid-plane; therefore, the quadrupole structure imprinted in the kinematics of the stars (far away from the mid-plane) becomes weaker. The introduction of uncertainties in radial velocity and proper motion does not drastically change the overall distribution of ⟨V_R⟩ (see the top-right and bottom-left panels in Fig. 9). We find that the bulk of the effect from observational errors comes from the uncertainty in parallax. A signal-to-noise ratio (S/N) of 10 causes an uncertainty of ∼1 kpc along the LoS near the Galactic centre. Since $V_{\mathrm{R}}=\mathit{V}·\widehat{R}$ is the projection of a star’s velocity on the cylindrical radial unit vector, which has a singularity at the Galactic centre, V_R gets contaminated by other velocity components as results of the shift of its observed position around the Galactic Centre. In Fig. 9 (bottom right panel), we further demonstrate how a scenario where the effect of the intrinsic V_R of the particles is eliminated by setting all the particles to artificially have zero V_R prior to the application of observational effects reproduces the central observed quadrupole. Thus, the ‘observed’ quadrupole present in the final mock dataset (bottom middle panel of Fig. 9) is essentially a manifestation of the intrinsic fainter quadrupole which is greatly enhanced by the uncertainties in parallax measurement. This clearly outlines the pitfall of inferring MW’s bar properties when using the stellar kinematic information from the Gaia DR3 if proper care is not taken when considering the observational biases of the Gaia survey.

Fig. 9. Investigating the effects of different uncertainties: Distribution of ⟨VR⟩ in the (x − y)-plane for the Gaia-like mock dataset at t = 3 Gyr, showcasing how the observed quadrupole feature is affected by errors in proper motion, radial velocity (VR), and parallax. Top-left panel: stars in the snapshot (no uncertainties and no magnitude cut). Top-middle panel: magnitude-limited G ≤ 15.5 mock dataset with no errors. The top-right and the bottom-left panels show the mock datasets with only proper motion (PM) uncertainty and both the PM and radial velocity uncertainties, respectively. The bottom-middle panel shows the mock dataset resampled according to the errors in the Gaia DR3 RGB subsample (see Section 4.2). The bottom-right panel shows the mock dataset where all the particles have had their VR set to zero prior to applying observational effects. For further details, see Sect. 4.3.

Lastly, we measured the strength and the length of the quadrupole (for definitions, see Sect. 3.1) for the MW, using our selected sample of RGB stars (for details of the sample selection, see Sect. 4.1). The strength of the quadrupole, S_quadrupole = 0.37, and the extent, R_quadrupole = 4.75 kpc, for the selected sample of the RGB stars in the MW. Now, if we assume the empirical relations for the S_bar − S_quadrupole and the R_bar − R_quadrupole, as obtained in Sect. 3.1 (see Fig. 2) from a wide variety of simulated bar models, the inferred S_bar and R_bar would be ∼0.57 and ∼4.5 kpc, respectively. Furthermore, the computed φ₄ for the selected sample of RGB stars turns out to be ∼5°, thereby implying the MW’s bar will be at an angle of ∼5° with respect to the Sun (for details, see Sect. 3.2). At this point, it is tempting to estimate the ‘true’ bar strength and length for the MW, given how the Gaia-like mock data overestimate these values by 35–45 per cent, as shown previously. However, for a rigorous inference, first we need to check the universality and robustness of the Gaia-like mocks produced from different high-resolution simulated MW-like galaxies. This is beyond the scope of this work and will be taken up in a future work.

5. Summary and future direction

In summary, we investigated the formation and temporal evolution of the quadrupole feature in the stellar mean radial velocity field, excited by an m = 2 stellar bar. In addition, we carried out a thorough study to test whether this quadrupole feature can be used as a robust kinematic diagnostic to put stringent constraints on the bar properties. To achieve this, we used a suite of 14 isolated, collisionless N-body models (having both thin and thick stellar discs) as well as a sample of barred galaxies from the TNG50 cosmological simulation. Finally, we investigated the possibility of constraining the MW’s bar properties from the quadrupole feature in the stellar ⟨V_R⟩ field. Our main findings are listed below.

In the absence of observational errors:

• A prominent stellar bar always excites a ‘butterfly’ or quadrupole structure in the stellar mean radial velocity ⟨V_R⟩ field. In this work, we devised a quantitative method to measure the strength (S_quadrupole) and length (R_quadrupole) via a Fourier decomposition of the ⟨V_R⟩ distribution. The strength and the length of the quadrupole are strongly correlated with the strength and length of the bar, irrespective of the choice of thin- or thick-disc stars. In addition, the orientation of the quadrupole in the stellar kinematics (calculated from the m = 4 Fourier phase angle) robustly captures the bar orientation in the density distribution. These trends hold true for all barred models (isolated or from the TNG50) used in this work.

• The strengths of the bar and the quadrupole tend to follow an empirical linear scaling relation of the form Y = AX + B (A = 0.96 ± 0.01; B = −0.18 ± 0.005). Similarly, the lengths of the bar and the quadrupole tend to follow an empirical linear scaling relation of the form Y = AX + B (A = 1.06 ± 0.01; B = 0.01 ± 0.01). However, the best-fit values of these linear relations change (within ∼10 − 25 per cent), depending on the choice of thin- or thick-disc stars.

• The presence of transient spirals systematically overestimates the length of the quadrupole feature (with a median relative difference of ∼25 per cent). In such a bar+spiral dynamical scenario, as in the case of the MW, the constancy of φ₄ (within ∼5 − 8°) serves as a much robust proxy for measuring the length of the quadrupole and, in turn, inferring the bar length.

In presence of observational errors:

• In Gaia-like mock datasets, constructed from the simulated model while incorporating the dust extinction and the broad trends of Gaia -like uncertainties and systematic errors, the quadrupole properties are overestimated by ∼35 − 45 per cent when compared with their ‘true’ values. We demonstrate that the majority of this effect comes from the uncertainty in parallax measurement from the Gaia survey.

In conclusion, we have demonstrated that the quadrupole feature in the ⟨V_R⟩ distribution is indeed an excellent kinematic diagnostic that enables us to put stringent constraint on the bar properties (strength, length, and orientation), provided there are no significant observational errors involved. It will be worthwhile to check whether the scaling relations derived in this work between the length and strength of the bar and the quadrupole hold true for a diverse bar models. Furthermore, we caution that inferring the MW’s bar properties by using the stellar kinematic information from Gaia DR3 when proper care is not taken with respect to the observational errors (predominantly the uncertainty in parallax measurement) of the Gaia survey can result in misleading assumptions on MW dynamics. We point out that this effect can be reduced by considering stars with more stringent cuts on parallax error. However, this comes at a cost for the number of available stars near the Galactic centre. With future data releases of Gaia DR4, this scenario is expected to improve. Alternatively, combining the Gaia parallax estimates with photometric surveys, for example, those in Zhang & Green (2025) can also improve this scenario.

Acknowledgments

We thank the anonymous referee for useful comments which helped to improve this paper. S.G. acknowledges funding from the IIT-Indore, through a Young Faculty Research Seed Grant (project: ‘INSIGHT’; IITI/YFRSG/2024-25/Phase-VII/02). T.K., S.G., and G.M.M. acknowledge funding from the Alexander von Humboldt Foundation, through a Sofja Kovalevskaja Award. P.D.M., D.K., and M.H. acknowledge the support of the French Agence Nationale de la Recherche (ANR), under grant ANR-13-BS01-0005 (project ANR-20-CE31-0004-01 MWDisc). This work has made use of the computational resources obtained through the DARI grant A0120410154 (P.I.: P. Di Matteo).

References

Appendix A: Temporal evolution of quadrupole properties

Figure A.1 shows the distribution of the stellar mean radial velocity (⟨V_R⟩) in the (x − y)-plane, calculated at t = 9 Gyr, for all 13 thin+thick models with varying f_thick values. A prominent bar is always accompanied by a prominent quadrupole feature in the ⟨V_R⟩ distribution. In each case, the extents of the bar and the quadrupole feature as well as the orientation of the quadrupole feature are indicated. In all cases, the orientation of the quadrupole feature traces accurately the bar orientation, thereby demonstrating that the quadrupole feature is a robust kinematic tracer of bar properties. Similarly, in Fig. A.2, we show the corresponding distribution of the ⟨V_R⟩ in the (x − y)-plane, for a sample of TNG50 MW-like barred galaxies. The quadrupole feature traces accurately the bar orientation for these TNG50 galaxies as well. In addition, Fig. A.3 shows the temporal evolution of the strength and extent of the quadrupole (for definitions, see Sect. 3.1) for all 14 isolated thin+thick bar models.

Fig. A.1. Distribution of the mean radial velocity, ⟨VR⟩ in the (x − y)-plane, calculated at the end of the simulation run (t = 9 Gyr), for all thin+thick models considered here. The dashed black lines denote the contours of constant surface density. The cyan dashed circle denotes the bar length, Rbar, and the magenta dashed circle denotes the extent of the quadrupole feature, Rquadrupole. The magenta points denote the spatial distribution of the phase-angle of the m = 4 Fourier moment (φ4). Top row: rthickS models. Middle and bottom rows: rthickE and rthickG models, respectively. The thick-disc mass fraction (fthick) varies from 0.1 to 0.7 (from left to right panels). A prominent quadrupole feature is present in all thin+thick models considered here.

Fig. A.2. Distribution of the mean radial velocity, ⟨VR⟩ in the (x − y)-plane, for some of the barred galaxies from the TNG50 simulations considered in this work. For details, see the text. The dashed black lines denote the contours of constant surface density. A prominent bar is always accompanied by a prominent quadrupole feature in the ⟨VR⟩ distribution. The colour bar is same as in Fig. A.1.

Fig. A.3. Temporal evolution of the strength of the quadrupole, Squadrupole and the extent of the quadrupole, Rquadrupole (normalised by the mean disc scale length, ⟨Rd⟩), for thin+thick models with different fthick values (see the colour bar), as a function of bar age (tbar − age). Left panels: rthickS models. Middle and right panels: rthickE and rthickG models, respectively. The strength and extent of the quadrupole for the model sim6 are shown in blue lines (see the left panels).

Appendix B: Further comparison of Gaia-like mock datasets

Figure B.1 shows the on-sky distribution of stars in Galactic coordinates for the Gaia-like mock dataset created from the bar model sim6 and our selected sample of Gaia DR3 stars. Even a visual inspection reveals that the final mock dataset has a spatial distribution that closely mimics that of the Gaia DR3 sample once we incorporate the dust extinction, Gaia selection function, and broad trends of Gaia-like uncertainties and systematic errors.

Fig. B.1. On-sky map in Galactic coordinates, for the Gaia-like mock dataset created from the bar model sim6 (top panel) and our selected sample of Gaia DR3 stars (bottom panel).

All Tables

All Figures

Fig. 1. Tracing the bar with the quadrupole feature: distribution of stellar surface density (left panel), mean radial velocity, ⟨VR⟩ (middle panel), and radial velocity dispersion, σR (right panel), calculated in the (x − y)-plane, for the rthick0.0 model, at the end of the simulation run (t = 9 Gyr). Dashed black lines denote the contours of constant surface density. The cyan dashed circle denotes the bar length, Rbar, and the magenta dashed circle denotes the extent of the quadrupole feature, Rquadrupole. The magenta points denote the spatial distribution of the phase-angle of the m = 4 Fourier moment (φ4). The bar excites a prominent quadrupole pattern in the mean radial velocity field, and the orientation of the quadrupole pattern agrees fairly accurately with the orientation of the bar. We used a galactocentric cylindrical coordinate system (R, ϕ, z) to calculate ⟨VR⟩ and σR.

In the text

Fig. 2. Tracing the bar properties with the quadrupole feature: correlation between the bar strength, Sbar and the strength of the quadrupole, Squadrupole (left panel), and correlation between the bar length, Rbar and the length of the quadrupole, Rquadrupole (right panel), computed using all isolated thin+thick models and the TNG50 barred galaxies (see the legend). The black dash line denotes the best-fit straight line (of the form Y = AX + B) while the grey shaded region denotes the 5-σ scatter around the best-fit line. The properties of the bar (strength and extent) remain strongly correlated with the properties of the quadrupole structure (Pearson correlation coefficient, ρ > 0.75).

In the text

Fig. 3. Dependence on the choice of the thin- or thick-disc stars: correlation between the bar strength, Sbar and the strength of the quadrupole, Squadrupole (top row), and correlation between the bar length, Rbar and the length of the quadrupole, Rquadrupole (bottom row), computed using thin-disc particles (left panels) and thick-disc particles (right panels), for all thin+thick models and the sim6 model (see the legend). The black dash line denotes the best-fit straight line (of the form Y = AX + B), while the grey shaded region denotes the 5-σ scatter around the best-fit line. Regardless of the choice of thin or thick-disc stellar particles, the length and strength of the quadrupole remain strongly correlated with the length and the strength of the bar.

In the text

Fig. 4. Distribution of the mean radial velocity (⟨VR⟩) in the (x − y)-plane, calculated at t = 9 Gyr for the model rthick0.0, with the bar placed at different viewing angles with respect to a hypothetical observer (shown in diamond) at a solar-like position (R = −8 kpc, ϕ = 0°, z = 0). The cyan circles in each sub-panel denote the variation in the phase-angle of the m = 4 Fourier coefficient (φ4), computed from the distribution of ⟨VR⟩. The dashed straight line denotes the ‘true’ bar angle, computed from the intrinsic density distribution of the stellar particles.

In the text

Fig. 5. Top panel: radial variation in the amplitude of the m = 4 Fourier moment of the mean radial velocity, ⟨VR⟩, calculated at t = 9 Gyr for the model rthick0.0, while putting the bar at different orientations (see the legend). Bottom panel: radial variation in the corresponding phase-angle (φ4) of m = 4 Fourier moment of the mean radial velocity, ⟨VR⟩. The points denoting the radial variation computed from the mean radial velocity, ⟨VR⟩, while the horizontal lines denote the corresponding true bar orientation, and the shaded region around the true value denotes a 8° scatter. ‘Half-extent’ refers to the scenario where stars falling only in the negative half (x ≤ 0) are considered while computing the Fourier moments. For further details, see the text. The vertical dash-dotted black line denotes the length of the bar (Rbar = 5.75 kpc).

In the text

Fig. 6. Distribution of the mean radial velocity, ⟨VR⟩ in the (x − y)-plane, computed at different times when the model rthick0.0 harbours a strong bar and prominent spirals. Dashed black lines denote the contours of constant surface density. The magenta points denote the spatial distribution of the phase-angle of the m = 4 Fourier moment (φ4). Within the bar region, φ4 remains constant whereas in presence of prominent spirals, φ4 shows a characteristic modulation and this trend holds for all bar+spiral scenarios shown here.

In the text

Fig. 7. Correlation between the two methods of measuring the quadrupole extent: Rquadrupole, measured solely from the amplitude of the m = 4 Fourier coefficient of the mean radial velocity ⟨VR⟩, and Rquadrupole, φ4, measured solely from the constancy of the phase-angle, φ4 of the mean radial velocity ⟨VR⟩, as a function of bar age (tbar-age, see the colour bar). Both the quantities, Rquadrupole and Rquadrupoleφ4, are normalised by the mean disc scale length, ⟨Rd⟩. For details, see the text. The dotted straight line denotes the 1:1 relation. The colour bar denotes the age of the bar.

In the text

Fig. 8. Comparison between the Gaia-like mock data and the Gaia DR3 sample. Top panels: ⟨VR⟩ distribution and the bottom panels showing the stellar density, both in the (x − y)-plane, calculated for all the stellar particles from model sim6 (left column), the Gaia-like mock dataset from model sim6 (middle column) and the selected sample of RGB stars from the Gaia Data Release 3. As shown in this work, the mock dataset qualitatively reproduces the behaviour seen in the Gaia sample, notably, observational effects increase the extent of the quadrupole and reduce the angle of its main axis compared to y = 0. The green lines in the left column corresponds to the major and minor axis of the bar, angled at 28°. The black filled circle in the right column denotes the assumed solar location (for further details, see Sect. 4.1).

In the text

Fig. 9. Investigating the effects of different uncertainties: Distribution of ⟨VR⟩ in the (x − y)-plane for the Gaia-like mock dataset at t = 3 Gyr, showcasing how the observed quadrupole feature is affected by errors in proper motion, radial velocity (VR), and parallax. Top-left panel: stars in the snapshot (no uncertainties and no magnitude cut). Top-middle panel: magnitude-limited G ≤ 15.5 mock dataset with no errors. The top-right and the bottom-left panels show the mock datasets with only proper motion (PM) uncertainty and both the PM and radial velocity uncertainties, respectively. The bottom-middle panel shows the mock dataset resampled according to the errors in the Gaia DR3 RGB subsample (see Section 4.2). The bottom-right panel shows the mock dataset where all the particles have had their VR set to zero prior to applying observational effects. For further details, see Sect. 4.3.

In the text

Fig. A.1. Distribution of the mean radial velocity, ⟨VR⟩ in the (x − y)-plane, calculated at the end of the simulation run (t = 9 Gyr), for all thin+thick models considered here. The dashed black lines denote the contours of constant surface density. The cyan dashed circle denotes the bar length, Rbar, and the magenta dashed circle denotes the extent of the quadrupole feature, Rquadrupole. The magenta points denote the spatial distribution of the phase-angle of the m = 4 Fourier moment (φ4). Top row: rthickS models. Middle and bottom rows: rthickE and rthickG models, respectively. The thick-disc mass fraction (fthick) varies from 0.1 to 0.7 (from left to right panels). A prominent quadrupole feature is present in all thin+thick models considered here.

In the text

	Fig. A.2. Distribution of the mean radial velocity, ⟨VR⟩ in the (x − y)-plane, for some of the barred galaxies from the TNG50 simulations considered in this work. For details, see the text. The dashed black lines denote the contours of constant surface density. A prominent bar is always accompanied by a prominent quadrupole feature in the ⟨VR⟩ distribution. The colour bar is same as in Fig. A.1.
In the text

Fig. A.3. Temporal evolution of the strength of the quadrupole, Squadrupole and the extent of the quadrupole, Rquadrupole (normalised by the mean disc scale length, ⟨Rd⟩), for thin+thick models with different fthick values (see the colour bar), as a function of bar age (tbar − age). Left panels: rthickS models. Middle and right panels: rthickE and rthickG models, respectively. The strength and extent of the quadrupole for the model sim6 are shown in blue lines (see the left panels).

In the text

	Fig. B.1. On-sky map in Galactic coordinates, for the Gaia-like mock dataset created from the bar model sim6 (top panel) and our selected sample of Gaia DR3 stars (bottom panel).
In the text