Gaia DR3 data consistent with a short bar connected to a spiral arm

We use numerical simulations to model Gaia DR3 data with the aim of constraining the Milky Way bar and spiral structure parameters. We show that both the morphology and the velocity field in Milky Way-like galactic disc models are strong functions of time, changing dramatically over a few tens of Myr. This suggests that by finding a good match to the observed radial velocity field, v_R(x,y), we can constrain the bar-spiral orientation. Incorporating uncertainties into our models is necessary to match the data; most importantly, a heliocentric distance uncertainty above 10-15% distorts the bar's shape and v_R quadrupole pattern morphology, and decreases its apparent angle with respect to the Sun-Galactocentric line. An excellent match to the Gaia DR3 v_R(x,y) field is found for a simulation with a bar length R_b\approx3.6 kpc. We argue that the data are consistent with a MW bar as short as ~3 kpc, for moderate strength inner disc spiral structure (A_2/A_0\approx0.25) or, alternatively, with a bar length up to ~5.2 kpc, provided that spiral arms are quite weak (A_2/A_0\approx0.1), and is most likely in the process of disconnecting from a spiral arm. We demonstrate that the bar angle and distance uncertainty can similarly affect the match between our models and the data - a smaller bar angle (20 deg instead of 30 deg) requires smaller distance uncertainty (20% instead of 30%) to explain the observations. Fourier components of the face-on density distribution of our models suggest that the MW does not have strong m=1 and/or m=3 spirals near the solar radius.


INTRODUCTION
The most significant stellar component of the Milky Way (MW) disc is its central bar.Stellar bars are non-axisymmetric elongated features present in roughly two-thirds of disc galaxies in the local universe (Knapen et al. 2000;Menéndez-Delmestre et al. 2007;Marinova & Jogee 2007;Sheth et al. 2008;Masters et al. 2011;Erwin 2018).However, the genuine properties of bars, such as their length, pattern speed, and strength, are difficult to resolve from the inner disc, bulge and spiral arms due to their mutual interconnection.
The MW bar was initially identified in near-infrared data (Blitz & Spergel 1991) and in the study of gas kinematics (Binney et al. 1991).Due to the Sun's position within the Galactic plane, it has been difficult to study the bar directly from observations.At first, the MW bar was found to be quite short, with a half-length   ∼ 2.5 kpc, by studying the peaks in the radial distribution of CO gas emission in the inner galaxy (Blitz & Spergel 1991), and by decomposing the distributions from NIR stellar surveys (Wegg et al. 2015;Portail et al. 2017), comparing MW ℓ-v diagrams with hydrodynamical simulations including the effect of the bar (Sormani et al. 2015;Li et al. 2016), and explaining the Hercules moving group with the bar's corotation (CR) (Portail et al. 2017;Monari et al. 2019) or the 4:1 Outer Lindblad Resonance (OLR) (Hunt & Bovy 2018) of a long slow bar, rather than the 2:1 OLR in the case of a faster bar.
The transformative Gaia DR2 and DR3 (Gaia Collaboration et al. 2018aCollaboration et al. , 2023a) ) datasets revealed arches, ridges, and streams in velocity and action space (Antoja et al. 2018;Kawata et al. 2018;Quillen et al. 2018b;Laporte et al. 2019;Bland-Hawthorn et al. 2019;Brown 2021;Queiroz et al. 2021;Poggio et al. 2021), showing unambiguously that the MW disc was out of equilibrium.This confirmed previous expectations based on incomplete pre-Gaia dataset (e.g., Hipparcos, RAVE, and SDSS), that a lot of disc phase-space structure could be explained as phase wrapping (mostly from the effect of the Sagittarius Dwarf galaxy, hereafter Sgr; Ibata et al. 1994;Laporte et al. 2018;Tepper-García et al. 2022), rather than self-gravity, e.g., the arches in the u-v plane (Minchev et al. 2009;Gómez et al. 2012a,b), clumps in the u-v plane (Quillen et al. 2009), disc asymmetries in the vertical direction resembling bending and breathing modes (de la Vega et al. 2015).Many of these structures, however, have also been found consistent with the effect of a slow, long bar (e.g., Fragkoudi et al. 2019;Khoperskov et al. 2020;Sanders et al. 2019;D'Onghia & L. Aguerri 2020;Kawata et al. 2021;Khoperskov & Gerhard 2022), although it has not been easy to break the degeneracy between the tidal effect of an external perturber (e.g., Sgr) and internal perturbations from disc asymmetries, such as the bar and self-sustained spiral arms (e.g, Carrillo et al. 2018Carrillo et al. , 2019;;Gaia Collaboration et al. 2018b;Laporte et al. 2020;Hunt et al. 2022).This drastic change in the bar length and pattern speed estimates, when measured directly from data in the inner disc, rather than modelling the local velocity field, has been rather puzzling.One way to understand this discrepancy was proposed recently by Hilmi et al. (2020).These authors showed that simulated galactic bars exhibit fluctuations in length, amplitude and pattern speed, due to the periodic bar overlap with the inner spiral structure, as already noted by Quillen et al. (2011).Using the ellipse-fitting ( prof ) and Fourier decomposition ( =2 ) methods (Athanassoula & Misiriotis 2002), along with overdensity contour maps ( cont ) of the central bar of MW-like simulations, Hilmi et al. (2020) showed that, while the bar is temporarily connected to a spiral arm it can appear up to twice its true size (see also Petersen et al. 2023), slowing down at the same time but by a smaller fraction, thus causing the ratio R =    /  to become less than 1.This gives rise to "ultrafast" bars, which have been found in observations (Buta & Zhang 2009;Aguerri et al. 2015) but theoretically deemed unphysical.Unlike the  1 stellar orbits which support the bar, orbits outside the bar's CR are perpendicular to its major axis, and so R < 1 is not allowed (Contopoulos & Papayannopoulos 1980;Contopoulos 1980).The work by Hilmi et al. (2020) could then explain the observed ultrafast bars with an overestimation of their length, if they happened to be connected to spiral arms.Indeed, an investigation by Cuomo et al. (2021) of a set of disc galaxies with ultrafast bars from the CALIFA survey, including those from Aguerri et al. (2015), found that they become regular fast bars when the bar length measurement proposed by Lee et al. (2020) was used, based on the analysis of the maps tracing the transverse-to-radial force ratio   (, ) of the galaxy (Combes & Sanders 1981).Another technique to overcome the biases induced by traditional bar length measurements was recently proposed by Petersen et al. (2023).Dubbed 'dynamical length', this method allows to measure the extent of  1 orbits, which defines the unambiguous bar length.
While the MW bar length has been estimated to be ∼ 5 kpc or even larger (e.g., Wegg et al. 2015;Li et al. 2016;Portail et al. 2017), the more recent work by Lucey et al. (2023) found, via orbit integration, that trapped bar orbits extend out to only ≈ 3.5 kpc, although there is an overdensity of stars at the end of the bar, out to 4.8 kpc, which could be related to an attached spiral arm.Another recent study on exploring the bar pattern speed indirectly from the effect on the tidal stream of the Hyades (Thomas et al. 2023) found Ω  ≈ 55 km/s, which is in stark contrast to the direct Tremaine-Weinberg (TW) method measurements (e.g., Bovy et al. 2019;Sanders et al. 2019).Both of the above results are very much in line with the predictions by Hilmi et al. (2020).
Inspired by Gaia Collaboration et al. (2023a), the aim of this work is to model the Gaia DR3 radial velocity field as a function of disc position and find out what we can learn about the Galactic bar length, as well as its orientation with respect to the spiral structure.Gaia Collaboration et al. (2023a) showed that the kinematic manifestation of the MW bar, namely the quadrupole, or butterfly-like radial velocity pattern, when the disc is viewed face-on, seems to be aligned with the Sun-GC line, implying that the bar angle,   , is close to zero.These authors also pointed out that the apparent orientation of the quadrupole is due to distance uncertainty and, based on comparison to simulations, argued for bar angle of about 20 • .Indeed, the consensus agrees on a tilted bar with respect to the Sun-GC line in the range 20-30 • ahead of the Sun (Bland-Hawthorn & Gerhard 2016).
Here, we also seek to explain this disagreement with a set of diverse hydrodynamic simulations of MW-like galaxies by accounting for uncertainties in observable measurements similar to those present in the Gaia DR3 dataset.

GAIA DR3 DATA SELECTION
We use the third data release from the Gaia (ESA) space observatory to study the velocity map of the inner disc region.The unfiltered Gaia DR3 dataset consists of over 33 million stars with 6-dimensional phase space information (Gaia Collaboration et al. 2023a).We use multiple quality criteria to ensure the reliability of the positions and velocities of the MW stars, needed for our analysis.More specifically, we make quality cuts in the renormalized unit weight error (RUWE) < 1.4 (Gaia Collaboration et al. 2023b), rejection of duplicated sources (determined by the Gaia cross-matching algorithm; Torra et al. 2021), and retention of only five-parameter astrometric solutions (astrometric_params_solved = 31, Gaia Collaboration et al. 2023b).
To calculate positions and velocities in the galactocentric restframe, we assumed an in-plane distance of the Sun from the Galactic centre of 8.2 kpc, a velocity of the local standard of rest (LSR), of 240 km s −1 (Reid et al. 2014), and a peculiar velocity of the Sun with respect to the LSR,  ⊙ = 11.1 km s −1 ,  ⊙ = 12.24 km s −1 , and  ⊙ = 7.25 km s −1 (Schönrich et al. 2010).The heliocentric distances were taken from the (Bailer-Jones et al. 2021) catalogue.
Additional quality cuts include constraints on observable velocities and the heliocentric distance, made to ensure the data are mostly free from observation biases stemming from our location in the MW disc.To find out how these cuts affect the galactocentric mean radial velocity map,   (, ), in Fig. 1 we plot six panels with different combinations of cuts in the uncertainty in proper motion,   , distance,   , and line-of-sight velocity,    .It can be seen in the figure that the parameter which most significantly affects both structure in velocity space and the covered disc area is the distance uncertainty, as any cuts below 20% begin to substantially decrease the amount of data beyond the Galactic center (see left column of Fig. 1).To compare to our simulations, we use two distance uncertainty cut:   < 20% 1 (as in Gaia Collaboration et al. 2023b,a) and   < 10%, which shows a different structure inside  = 4 kpc.Looking at the various proper motion cuts (in both RA and DEC), we notice that   < 5% seems to intensify features in the radial velocity field skewed towards the solar neighbourhood, so we maintain a 10% uncertainty bound.Lastly, decreasing    from 5 to 2 km s −1 makes little difference in the radial velocity field of the data (see also Kordopatis et al. 2023), other than to slightly lessen the intensity of the butterfly pattern in the bar region, while it cuts the data sample in half.We thus maintain a 5 km s −1 error limit in   .After all these quality cuts, we are left with a star count of about 17.4 or 15.3 million for the   < 20% and   < 10% cuts, respectively.The variation of the radial velocity field with the above-described combinations of uncertainty cuts is shown in Fig. 1, as indicated in each panel.

SIMULATIONS
We use three hydrodynamical simulations of barred spiral galaxies (two in the cosmological context) with disc properties similar to those of the MW.We only consider the last ∼ 1.4 Gyr of evolution in these models, as we aim at matching the Galactic disc dynamical state at redshift zero.Since the latter is a very strong function of time, we 1 When referring to uncertainty cuts, we use   < 20%, etc., to mean   < 0.2, etc., interchangeably within the paper.examine closely spaced time outputs (from 4.5 to 10 Myr) in each simulation.
Model1 and Model2, as introduced below, were used under the same names by Hilmi et al. (2020) to study the bar length fluctuations due to the bar-spiral interaction resulting from their different pattern speed.We keep the same names here for consistency, although the scaling we do is slightly different, as described below.Model3 is a preassembled stellar disc simulation, including gas dynamics, star formation and chemical evolution.The face-on and edge-on views of the three models for the last time outputs can be seen in Fig. 2.
The angle the Galactic bar semi-major axis makes with the Sun-GC line is referred to as the bar angle,   .In our simulations, we assume a bar angle of   = 30 • , which is consistent with the upper limit of measurements found using distributions of red clump giants from surveys such as the Via Lactea, (OGLE) III, and 2MASS, e.g., 27 The velocities of our models are scaled so that the rotation curve is 240 km/s, while we scale the distances (and thus the bar lengths) such that the observed radial velocity field,   (, ), is reproduced as best as possible.The latter results in different bar lengths, as indicated below.

Model1
The first simulation (Model1) we use was introduced by Buck et al. (2018) as a higher-resolution of the galaxy g2.79e12 from the NI-HAO project (Wang et al. 2015), with a notable boxy/peanut-shaped bulge similar to that of the MW (Buck et al. 2018(Buck et al. , 2019c)).The simulation was made with a modified version of the smoothed particle hydrodynamics (SPH) code GASOLINE2.0(Wadsley et al. 2017).The updated hydrodynamics adopt a metal diffusion algorithm between particles (Wadsley et al. 2008).Star formation in this model follows that described by Stinson et al. (2006) for dense and cool gas.Two modes of stellar feedback are used in the simulation (Stinson et al. 2013); one modelled from luminous young stars before any supernovae, and the other mode from supernovae after 4 Myr of star formation.
Model1 has a flat rotation curve at   = 324 km/s, which we scaled down to match current estimates for the MW of   = 240 km/s (Bland-Hawthorn & Gerhard 2016).
Unlike in Hilmi et al. (2020), we do not scale the distances down, since the Gaia DR3 radial velocity field is matched well with the original bar length,   ≈ 5.2 kpc, as measured from the minimum of its fluctuation in the last ∼ 1.4 Gyr (see Hilmi et al. 2020 for details).The time outputs from this simulation are every 6.9 Myr.

Model2
The second simulation (Model2) is galaxy g106 from a suite of 33 hydrodynamic simulations by Martig et al. (2012) made by extracting merger and accretion histories of a particular halo from a cosmological simulation, then re-simulating with the Particle-Mesh code (Bournaud & Combes 2002, 2003) at high resolution with a galaxy in place of the halo (zoom-in technique introduced by Martig et al. (2009)).This simulation has a mass resolution of 1.5 × 10 4 M ⊙ for gas particles, 7.5 × 10 4 M ⊙ for stars, and 3 × 10 5 M ⊙ for dark matter particles; the spatial resolution is 150 pc.Within the optical radius of 25 kpc, this simulation has a total stellar mass of ∼ 4.3 × 10 10  ⊙ and a dark matter mass of ∼ 3.4 × 10 11  ⊙ .This simulation has also been studied extensively due to its similarity to the MW (e.g.Martig et al. 2014a,b, Kraljic et al. 2012, Minchev et al. 2013, 2014, 2015, 2017;Carrillo et al. 2019;Hilmi et al. 2020).
Similarly to Hilmi et al. (2020), we scale this simulation so that the original rotation curve of  c ≈ 210 km/s matches the MW at 240 km/s and the bar length is fixed at ∼ 3.2 kpc at the final time, by scaling distances down by a factor of 1.46.In the 1.37 Gyr period considered, the bar length increases from ∼ 2.8 to ∼ 3.2 kpc with a time average of 3 kpc.Time outputs here are separated by 4.5 Myr.
We re-scale Model2 distances to shrink the bar, but do not do this for Model1, in order to show that both these models, with largely differing bar lengths, can reproduce the Gaia DR3 radial velocity field, the reason for this being the weaker spiral structure in Model1.

Model3
Our Model3 is an -body/hydrodynamical simulation of a disc galaxy with a total stellar mass and a rotation curve compatible with those of the MW.The simulation starts from a pre-existing axisymmetric stellar disc, including gas and star formation coupled with chemical evolution.The simulation lasts about 3 Gyr of which we consider the last 1.37 Gyr as in our other two models and assume this represents well the last 1.37 Gyr of MW evolution dynamically.A well-defined buckled bar is formed before the time period we consider, as can be seen in Fig. 2. The detailed description of the Model3 setup is as follows: Initially, stellar particles are redistributed following a Miyamoto-Nagai density profile (Miyamoto & Nagai 1975) that has a characteristic scale length of 4 kpc, vertical thicknesses of 0.2 kpc and mass of 4.5×10 10 M ⊙ .Also included is a live dark matter halo (5×10 6 particles), whose density distribution follows a Plummer sphere (Plummer 1911), with a total mass of 6.2 × 10 11 M ⊙ and a radius of 21 kpc.The choice of parameters leads to a galaxy mass model with a circular velocity of ≈ 220 km/s.The gas component is represented by an exponential disc with a scale length of 5 kpc and a total mass of 1.5 × 10 10 M ⊙ .The initial equilibrium state has been generated using the iterative method from the AGAMA software (Vasiliev 2019).A gaseous cell undergoes star formation if: i) the gas mass is > 2 × 10 5 M ⊙ , (ii) the temperature of the gas is lower than 100 K Figure 3. Bar length variation for Model3, measured using the  cont method of Hilmi et al. (2020).Very similar results are found with the  prof method from, e.g., Athanassoula & Misiriotis (2002) (not shown).After the bar is aligned along the x-axis in a face-on view as in the top panels of Fig. 2, each bar half is measured separately, with the "Left Side" (extending at  < 0) shown in the top panel and the "Right Side" (at  > 0) shown in the bottom.The different curves represent different thresholds from 20% to 80% (drop in overdensity along the bar semi-major axis).The bar true length is ∼ 3.6 kpc, estimated as the minimum of the fluctuations for a 50% threshold as done by Hilmi et al. (2020).Lower thresholds result in longer bar measurements, therefore 50% is the line right in the middle for any given time.Five main peaks are seen in both the left and right bar halves (a period of ∼ 270 Myr), however, a slight offset exists, since the two bar halves do not connect to spiral arms always at the same time.The pink lines indicate the times of time outputs shown in Fig. 4. The blue line indicates the time output for which we find the best match to the Gaia DR3 data (see Fig. 6).and (iii) if the cell is part of a converging flow.The efficiency of star formation is set to 0.05, i.e., 5% of the gas eligible to form a new star particle per dynamical time.We consider the ISM as a mixture of several species (H, He, Si, Mg, O, Fe, and other metals), which is sufficient for modelling the Galactic chemical evolution and the newborn stellar particles inherit both kinematics and elemental abundances of their parent gas cells.No chemical information is used in this work.
Following the chemical evolution models by Snaith et al. (2015) and Snaith et al. (2022), at each time step, for newly formed stars we calculate the amount of gas returned, the mass of the various species of metals, the number of SNII or SNIa for a given initial mass and metallicity, the cumulative yield of various chemical elements, the total metallicity, and the total gas released.Feedback associated with the evolution of massive stars is implemented as an injection of thermal energy in a nearby gas cell proportional to the number of SNII, SNI and AGB stars (see Khoperskov et al. 2021, for more details).The hydrodynamical part also includes gas-metallicity-dependent radiative cooling (see details in Khoperskov et al. 2021).The simulations were evolved with the -body+Total Variation Diminishing hydrodynamical code (Khoperskov et al. 2014).For the -body system integration and gas self-gravity, we used our parallel version of the TREE-GRAPE code (Fukushige et al. 2005) with multithread usage under the SSE and AVX instructions.In recent years we already used and extensively tested our hardware-accelerator-based gravity calculation routine in several galaxy dynamics studies where we obtained accurate results with a good performance (Khoperskov et al. 2018a,b;Saburova et al. 2018;Khoperskov et al. 2019).For the time integration, we used a leapfrog integrator with a fixed step size of 0.1 Myr.In the simulation, we adopted the standard opening angle  = 0.7.The dynamics of the ISM is simulated on a Cartesian grid with static mesh refinement and a minimum cell size of ≈ 10 pc in the Galactic plane.
Similarly to Model1 and Model2, we scale this simulation's rotation curve from 220 km/s to 240 km/s and do not scale the distances.The bar length during the period of time we consider is ∼ 3.6 kpc, as estimated from the minima of the 50% threshold in Fig. 3.As in Model1 and Model2, we chose this bar length as it happens to match well the Gaia DR3 radial velocity field.Time outputs here are separated by 10 Myr.

Spiral structure
The spiral arms of Model1 are more tightly wound and multi-armed (see Fig. 2 and Fig. 1 by Buck et al. 2018), compared to Model2 or Model3, where they are more open and dominated by two or four arms (see also Fig. 1 by Minchev et al. 2013), signifying that they are stronger.We measured the spiral structure overdensity from the ratio of the amplitude of the m=1,2,3 and 4 components to the m=0 Fourier component of the stellar density, as a function of galactic disc radius and that for three radii in Fig. A3.We found that the spiral overdensity is typically 10% for Model1 and about 20-25% Model2 and Model3.One important difference is that Model3 displays the strongest odd modes: m=1 and m=3, which correspond to a one-and three-armed spirals (see Fig. A3).
For the MW we expect spiral-arm overdensity of ∼ 15% from modelling the radial velocity field of RAVE data (Siebert et al. 2012), ∼ 25% from considering the migration rate of open clusters near the Sun (Quillen et al. 2018a), ∼ 20% from matching the radial velocity field of stars from a compilation of data (Eilers et al. 2020).These estimates are larger than the spiral overdensity of Model1, but consistent with our Model2 and Model3.

Matching the selection and uncertainties of Gaia DR3
In order to compare properly to the observations, we need to match the geometry of the Gaia DR3 sample and to introduce synthetic uncertainties into our models, consistently with the data.
We first transform our simulation data from the native galactocentric cylindrical coordinate system to a Galactic spherical coordinate frame centred on the Sun.This is done with the astropy.coordinatesmodule using the SkyCoord.transform_toobject method (see Astropy Collaboration et al. 2022).In this transformation we specify the position of the Sun's barycentre at (, , ) = (−8.2,0, 0) kpc in the galactocentric frame, as this is shifted to the origin of the new Galactic frame.To match the reference frame from which Gaia measures kinematic observables, we perform another simple coordinate transformation of the simulated data, taking it from Galactic spherical to ICRS coordinates.
To match the geometry of our Gaia sample, we picked the Sun position in the disc so that the bar is oriented at 30 • ahead of the Sun-galactocentric line (Bland-Hawthorn & Gerhard 2016), but we test the case of 20 • as well.Moreover, we made a cut in Galactic latitude, || < 1.0 • about the midplane, in order to exclude the dust-obscured regions in the data and match the Gaia footprint.This achieves a similar effect as our quality cuts in the Gaia DR3 data, which preferentially reject stars within the Galactic midplane.
The Gaia DR3 uncertainties for our dataset revealed roughly Gaussian distributions in proper motion and line-of-sight velocity uncertainties and a complex, skewed Gaussian distribution in distance uncertainty which is coupled to distance, as shown in the left panel of A1.This skewing towards larger distances was already shown by Gaia Collaboration et al. (2023a).As for the data, we introduced synthetic uncertainties in the line-of-sight velocity    = 5 km/s, in heliocentric distance   = 0.2 or   = 0.1, and in both proper motions   = 0.1, using Gaussian distributions.To model the distribution of relative uncertainties in distance (  /) in two different ways.First we use a Gaussian distribution of width 20% or 10% centered on zero, which is used for most of the figures.We also fit the data using skewed probability distributions fitting the data in different distance bins.These functions are then interpolated to create a continuous probability density function dependent on distance.The result is shown in the right panel of Fig. A1 and provides a very good match to the data on the left.After convolving the uncertainties, we converted from ICRS back to Galactocentric cylindrical coordinates.Now our models include the biases in the Gaia DR3 observables and can allow for proper comparison to the data.

Gaia DR3 radial velocity map
It has been previously shown (e.g., Carrillo et al. 2019;Fragkoudi et al. 2019;Bovy et al. 2019) that a central bar produces a quadruple pattern in the disc radial velocity field when viewed face-on.This signature was first identified for the MW by Bovy et al. (2019) and then more clearly by Queiroz et al. (2021), combining APOGEE (Majewski et al. 2017) spectroscopy with earlier Gaia DRs astrometry, using a few tens of thousands of stars.The advent of Gaia DR3 confirmed the existence of such a kinematic pattern using millions of stars, and extended to spiral arms all the way outside the solar circle.3) with the bar oriented at 30 • ahead of the Sun-galactocentric line.The middle and right panels show the effect of 15% and 30% distance uncertainty, respectively.The measurement errors, most importantly the distance uncertainty, cause the bar angle to shift significantly towards the Sun-galactocentric line, as indicated at the bottom of each panel.Moreover, the central contours are affected more, resulting in a distorted bar shape.The white contour marks the same density level, which changes from the 6th, to the 5th, to the 4th most dense contour for a   = 0, 20%, 30%, respectively.This is caused by the central density spread out due to the distance uncertainty increase.We measured 25 • for   < 15% and 12 • for   < 30% taking all density inside the bar length (3.6 kpc) into account.However, it is clear from the innermost couple of contours that the effect is even larger there due to the more circular initial distribution.
The expected orientation of the bar density (ellipse in Fig. 1) should lie along the line delineating negative from positive   lobes, as shown in Fig. 4.However, we find that the bar semi-major axis is aligned with the blue lobe instead.We will show later that this is the effect of the distance uncertainty.Note that once we cut at   < 10%, the inner ∼ 2 kpc velocity structure aligns better with the bar major axis.

Rapid variations in Galactic disc morphology and radial velocity field
When more than one perturber is present in a galactic disc, such as a central bar and spiral arms moving at different pattern speeds, one expects a strong variation of both the density and the velocity field on short timescales (e.g., Carrillo et al. 2019;Asano et al. 2022).Indeed, using Model1 and Model2, Hilmi et al. (2020) showed that the bar length, amplitude, and pattern speed can all fluctuate on a dynamical timescale consistent with the beat frequency between the bar and inner spiral modes.
In Fig. 3 we show the bar length evolution with time for our Model3.The two half lengths are measured separately with the one near the Sun shown in the top panel (Left Side) and the one farther from the Sun in the bottom (Right Side).About 5 fluctuations are seen from the number of peaks and troughs.The length measurement method used is  cont , tracing the drop in overdensity along the bar major axis, as described by Hilmi et al. (2020), who reasoned that the minimum in the 50% threshold in the density drop was closest to the true bar length, which happens when the spiral is fully disconnected from it.We see that within this definition, the bar fluctuates between ∼ 3.6 to ∼ 7 kpc in length, with a period of about 250-300 Myr (the beat frequency between the bar and the dominant inner spiral mode).
To understand how the bar length fluctuations seen in Fig. 3 affect the inner disc morphology and its radial velocity field, in Fig. 4 we plot the disc face-on view for five snapshots from our Model3.Those are separated by 40-120 Myr and are picked according to the relative orientation between the bar and spiral.The left column shows the stellar overdensity, computed as Σ(, ) = (Σ(, )−Σ 0 ())/Σ 0 (), where Σ(, ) is the density as seen in the top panel of Fig. 2 and Σ 0 is the azimuthally averaged density for radial bins of 0.3 kpc.
In the right column, we plot the galactocentric radial velocity field,   (, ).In all panels, the bar is oriented at 30 • with respect to the Sun-GC line for a Sun position indicated by the black dot at  = 8.2 kpc and  = 0.The two dashed circles show the solar radius and 4 kpc to guide the eye.We can see that in different rows, spirals have different orientations with respect to the fixed bar, due to their lower pattern speed.In the reference frame of the bar, spirals move counterclockwise, although Galactic rotation is clockwise.Note that the Gaia DR3 uncertainties and sample selection, as described in Sec.3.5, are not applied to Fig. 4.
The time outputs shown in Fig. 4 span about 290 Myr, starting and ending with a complete separation between the near bar half and the spiral arms.This corresponds to one period of the bar length fluctuations, seen in Fig. 3 (pink vertical lines).From top to bottom, the bar half nearer the Sun is well separated from the spiral arm and thus at a minimum in Fig. 3, in the process of connecting (second row), fully connected (third row) and thus a maximum in Fig. 3, in the process of disconnecting (fourth row), and again fully disconnected (bottom row).
It is easy to infer from both the overdensity and   plots the times for which the bar is separated from the spiral.In the top and bottom panels, the bar fits well within the 4-kpc dashed circle, while in the third row (fully connected) it is extending well beyond it and the positive   lobe covers roughly four times larger area.
The above-described variations outside the bar region in both density and velocity mean that comparison between the Gaia DR3 data and models should be done carefully, studying the detailed time evolution of the disc.To accomplish this, for all our models we use time outputs between 4.5 and 10 Myr, depending on the simulations (see Sec. 3).

Spiral arm in stellar mass near bar end due to overlap of multiple modes
To first order, when the bar is in the process of connecting to, or disconnecting from, a spiral arm (second and fourth rows in Fig. 4, respectively), a leading or trailing arm, respectively, appears attached to the bar.Overall, the velocity field seen in the right panels shows similar morphology in the positive   nearby lobe.It should be kept in mind, however, that multiple spiral modes with different multiplicity (typically m=1,2,3, and 4) and different pattern speeds, are always present just outside the bar, as seen in both numerical simulations (e.g., Sellwood & Sparke 1988;Quillen et al. 2011;Minchev et al. 2012;Hilmi et al. 2020) and observations (e.g., Elmegreen et al. 1992;Rix & Rieke 1993;Henry et al. 2003;Meidt et al. 2009).This has been shown to be the case also for Model1 and Model2 by Hilmi et al. (2020), by constructing power spectrograms.Therefore, the spiral overdensity seen in the mass in the left column of Fig. 4 is due to the overlap of all these modes and not caused by a single spiral pattern (although dominated by the strongest mode).It is thus possible that while the bar appears to be connecting to, or disconnecting from this apparent spiral (which, in fact, is an overdensity associated with the overlap of the multiple modes), signatures of both disconnecting and connecting spiral arms are present in the   field in the same snapshot.Indeed, in the second row of Fig. 4 we can see a trailing arm in positive   (right column) just outside the 4-kpc dashed circle, in addition to the leading arm seen in the density we reported above, although we see no spiral overdensity in the left panel.
We estimated the m=1,2,3, and 4 Fourier components from the face-on density of each model as functions of time with the results shown in Fig. A3.It can be seen there that as stated before, Model1 has significantly weaker spiral structure than the other models.While Model2 and Model3 have similar 2-and 4-armed modes, Model3 has stronger m=1 and m=3 components overall.We later argue that these odd modes are not expected to be very strong for the MW.

Effect of Gaia DR3 uncertainties on the bar shape and orientation
We now study the effect of Gaia DR3-like uncertainties on the disc morphology of our Model3, but the results are very similar for our other two Models.Fig. 5 shows the face-on stellar density when no uncertainties are considered (left), for an adopted distance uncertainty of   < 15% (middle), and   < 30% (right).
It is immediately obvious that the bar angle is strongly decreased from its true 30 • when synthetic uncertainties are included in the simulation, as already expected from Gaia mock catalogues (Romero- Gómez et al. 2015) and the work by Gaia Collaboration et al. (2023a).We measured ∼ 25 • for   < 15% and 12 • for   < 30% taking all density inside the bar length (3.6 kpc) into account.However, it  6), but with different combinations of distance uncertainty and bar angle, as indicated.When the implemented distance uncertainty is 20% and the bar angle is 30 • (leftmost column), as in the data (rightmost column), an upward kink inside the oval representing the bar is seen for the positive   lobe, which is not present in the data.To achieve the flatness of the transition between negative and positive   , we propose two solutions: decreasing the bar angle from 30 • to 20 • , while setting the distance error at 20% (second column), or using an uncertainty of 30% and a bar angle of 30 • (third column, as in Fig. 6).The latter solution appears to give a better match to the data (rightmost column), suggesting that Gaia DR3 distance uncertainties are underestimated.
is clear from the innermost couple of contours that the effect is even larger there (close to zero in the right panel), due to the more circular initial distribution.The white contour in each panel marks the same density level, which shifts from the 6th, to the 5th, to the 4th most dense contour with an increase in uncertainty:   = 0, 20%, 30%, respectively.This indicates that the central density spreads out as the distance uncertainty increase.This stretch in the initially circular central density contour can be also seen in the stellar velocity dispersion, as shown by Hey et al. (2023).
In addition to the decrease in bar angle, the uncertainties cause the bar to appear less centrally concentrated and offset from the galactic centre toward the direction of the Sun, as seen in the figure (especially for the 30% error).It is notable that the contours that encapsulate the bar are affected differently by   .The highest density contour is almost aligned with the Sun-GC line for both distance uncertainty cuts, which can be linked to its originally nearly circular shape.
As it can be already expected, we show in the next section that this apparent decrease in the bar orientation strongly affects the observed orientation of the central radial velocity field,   (, ), as well.

Matching the Gaia DR3 radial velocity field
Fig. 6 presents a snapshot at  = 12.4 Gyr from our Model3, exploring the effect of Gaia DR3 uncertainties and providing a match to the disc radial velocity field.The left column shows the model   (, ) (top) and the stellar overdensity (bottom), with no uncertainties included.Since this is a snapshot when the near half of the bar is connected to a spiral, the bar appears much longer than its true length, which at this particular time is   ≈ 3.2 (see minimum at  ≈ 12.37 Gyr in Fig. 3).Using the ellipse fitting method,  prof , we measured   ≈ 5.5 kpc apparent bar length, however, due to the gap present in the stellar overdensty along the bar semi-major axis (bottom left panel of Fig. 6), the  cont method introduced by Hilmi et al. ( 2020), measured a lower value, more consistent with the real length.It should be kept in mind that observationally  cont cannot be applied to the MW,2 thus a bar at this configuration will be likely miss-measured by a factor of ∼ 1.7.This also results in a larger positive radial velocity lobe, compared to when the bar is disconnected from the spiral structure, as it was already illustrated in Fig. 4.
The middle column of Fig. 6 shows the radial velocity field as in top left, but including 30% (top) and 15% (bottom) synthetic error, in addition to the uncertainties in heliocentric radial velocity,    , and proper motion,   (see Sec. 2).As expected from the decrease in bar angle caused by the uncertainty that we saw in Fig. 5, the butterfly pattern in the centre is rotated counterclockwise so that the bar semi-major axis (dashed line) passes through the negative velocity lobe (blue), instead of the interface between positive and negative (as in top left).We note that our results of the bar angle decreasing with distance error are in agreement with conclusions by Gaia Collaboration et al. (2023a) and Hey et al. (2023).The more important effect of the relative bar-spiral orientation, however, has not been discussed before.
Finally, in the right column of Fig. 6 we show the radial velocity field of Gaia DR3 data with   < 0.2 (top) and   < 0.1 (bottom).The middle column provides an excellent match to the entire Gaia DR3 dataset, especially inside the solar circle in the following: • the sizes and shapes of the negative and positive radial velocity lobes associated with the bar's near half; • the upward extending positive   arm attached to the positive   lobe; The second row shows Model2, which has a true size of ∼ 3 kpc (varying from ∼ 2.8 to ∼ 3.2 kpc over the examined time period).The bottom row shows Model3, with a true bar length of ∼ 3.6 kpc, which provided the best match the the Gaia DR3 data in Fig. 6 (we do not repeat that time output here).For all snapshots we used the same synthetic uncertainties as in Fig. 6.It is remarkable that this large range of bar lengths (from ∼ 3 to ∼ 5.2 kpc) among our three models can match relatively well the structure inside the solar circle.We attribute this mostly to the relative overdensity between the bar and spiral arms -the stronger the spiral arms, the smaller the bar that can reproduce the observations.
• the innermost 2-3 kpc for the uncertainty in data of 10% (15% in model, see Sec. 4.6), where the semi-major axis of the oval (3.5 kpc long) cleanly separates positive and negative velocities before it flattens sharply at ∼ 2 kpc from the Galactic centre (see red arrows); • the sharp decrease in positive   area below the positive   lobe; • the wide positive   arm along the left side of the plot, with a bifurcation in the upper left quadrant.
The reason in Fig. 6 we compared data with 10% and 20% distance uncertainty cut to 15% and 30% in the model is because increasing the distance uncertainty was the only way to achieve the flatness of the transition between negative and positive   lobes in Fig. 6 (see Sec.4.6 and Fig. 7 for the effect of 10% and 20% uncertainties).This suggests that the data distance uncertainties are underestimated, but see the next section for more discussion on this.

Interplay between distance uncertainty and bar angle
In Fig. 6 we showed a time output from our Model3 simulation, which provided an excellent match to the Gaia DR3 data.However, the implemented distance uncertainty in the simulation was 15% and 30%, rather than the 10% and 20% cuts in the data.To justify this, in the top left panel of Fig. 7 we show the same Model3 snapshot but with 10% and 20%.It can be seen that the flatness of the transition between negative and positive   cannot be achieved with the 20% error, although the difference between 10 and 15% is not so dramatic.Therefore, an uncertainty of 30% needs to be used, as we did in Fig. 6, to match the data cut of 20%.
We also considered the possibility that the bar angle is smaller than 30 • , which would then require smaller uncertainty in the simulation to align the bar with the Sun-GC line.In Fig. 7 we explore how well this particular Model3 snapshot (at 12.4 Gyr) matches the data when the implemented error is 20% as in the data, but changing the bar angle from 30 to 20 • , which can be seen as a lower limit (Bland-Hawthorn & Gerhard 2016;Gaia Collaboration et al. 2023a).
Indeed, we can see that, even though   = 20% combined with 30 • (top left panel) does not provide a good match to the data (top rightmost panel), when the angle goes down to 20 • the comparison with Gaia DR3 is much better, though arguably not as good as in the top-third panel (same as the top-middle panel of Fig. 6. Finally, using the Gaussian distance uncertainty modeled as in the data (Fig. A1), we show in Fig. A2 that we also require a bar angle of 20 • to match well the Gaia data.More work is needed to explore the interplay between bar angle and distance uncertainties.

Matches to data, considering 𝑣 𝑅 (𝑥, 𝑦) inside solar circle
Focusing mostly on the   structure inside the solar circle, we identified good matches to the Gaia DR3   field, by visually inspecting all snapshots from out three simulations in the studied period of 1.37 Gyr.A representative sample of those is displayed in Fig. 8, with the Gaia DR3 data shown in the lower-leftmost panel.It is remarkable that Model1, with a bar size   ≈ 5.2 kpc, gives similar   morphology as a bar as short as ∼ 3 kpc (Model2), or Model3's ∼ 3.6 kpc bar.This is, however, not surprising since we know that the fluctuation in bar parameters with time is strongly dependent on the strength of spiral structure outside the bar region (Hilmi et al. 2020).We scaled the distances in each model3 so that we could obtain good matches to the observed   (, ) field.This has naturally resulted in an arrangement, such that smaller bars are accompanied by stronger spiral structure (Model2 and Model3), as reported in Sec.3.4 and seen in Fig. A3 (see dominant m=2 mode).
In Fig. 9 we show all the times at which good matches to the data are obtained for our three models (some of which were displayed in Fig. 8), in order to understand if they always correspond to a certain bar-spiral orientation.The three blocks of two rows each, show the time evolution of the bar half-length for Model1 (top), Model2 (middle), and Model3 (bottom).The vertical lines indicate the times when good matches to the Gaia DR3 data take place.The higher frequency of oscillations in the top panels reflects the fast bar of Model1, causing it to meet the spiral arms more often.From the number of wave packets seen in Model2 and Model3 it is clear that both of these bars are slower than that of Model1, while Model2's bar is the slowest.
It can be seen from Fig. 9 that, in the same period of time, many more matches are found for Model1 than for the other models, as it should be expected if its bar encounters the spiral arms more often, as discussed above.For all models, the number of good matches is roughly equal to the number of peaks, i.e., to first order, it is expected that a good match occurs when the bar-spiral orientation is the same.But how do we determine the bar-spiral orientation?
As discussed in Sec.4.3, the spiral seen in the stellar mass (Figures 4 and 6) results from the constructive interference of all spiral modes of different multiplicities overlapping at a given time just outside the bar.We can see that in most matching time outputs from Model1 (20 out of 35) and Model2 (12 out of 17), the bar is to the right of the nearest peak, indicating it is in the process of separating from the spiral.
However, only 5 out of the 14 matches from Model3 are for a bar disconnecting from a spiral, according to the above criterion.This includes our best match (see Fig. 6), shown by the blue vertical line in the bottom panel of Fig. 9. Upon another inspection of the bottom panels of Fig. 8, we can see that the Model3 matches (mostly focused on the upward positive   leading arm, stemming from the positive   lobe), have a common flaw outside the solar radius: a trailing positive   arm extends also downward along the solar circle, which is not seen in the Gaia DR3 data, nor in the other two models.This may be due to a spiral mode present in this simulation, that either does not exist in the MW, or is simply not as strong as the one causing the upward   arm.In other words, the bar in Model3 is connecting to one spiral mode while disconnecting from another, as discussed in Sec.4.3 and seen in the bottom left panel of Fig. 6.While this would also happen for the other two models, the difference is that these two modes in Model3 are of similar strength, judging from the similar response seen in the radial velocity field.
This expectation is confirmed in Fig. A3, which shows the Fourier components estimated from the face-on disc density for modes m=1-m=4 as functions of time.We can see that in the range 7.2-9.3kpc, which is where the positive   arm that stretches downward in most matching snapshots of Model3 (bottom row of Fig. 8) is located, and which is not seen in the data, results from these odd modes.Indeed, Model3 shows the strongest m=1 and m=3 modes, corresponding to one-armed and three-armed spiral structure.The matches to the data appear to happen near m=1 and/or m=3 maxima, including the best match shown in Fig. 6 (blue vertical).This suggests that the MW lacks such strong m=1 or m=3 modes in the radial range shown.
Although we found the best match to the data in Model3 considering the overall radial velocity field, Model1 and Model2 are consistently showing better matches to both the upper and lower right quadrants of the   (, ) plane simultaneously.It should be noted that the spiral structure responsible for the radial velocity outside the solar circle is expected to be due to yet slower moving modes, different from the ones reaching the bar, which further complicates the problem.
We conclude that most likely the MW near bar side is in the process of disconnecting from a spiral arm, even though our best match is for a connecting one according to the bar length fluctuation with time (Fig. 9), but a disconnecting one according to the overdensity seen in the lower left panel of Fig. 6.Again, this complication is due to the presence of multiple modes and their interference as a function of time.More work is needed to understand better this behaviour.

DISCUSSION AND CONCLUSIONS
In this work we used three MW-like simulations of galactic discs to study the Gaia DR3 radial velocity field,   (, ).For all models we examined the last 1.37 Gyr of evolution, using frequent time outputs, from 4.5 to 10 Myr depending on the simulation.This allowed to resolve well the   (, ) time variation caused by the interaction of multiple patterns in the disc, most importantly for this project -the bar-spiral periodic overlap.Our models' true bar lengths, resulting when the bar is separated from the inner spiral structure, are about 5.2, 3, and 3.6 kpc, for Model1, Model2, and Model3, respectively.
Our results can be summarized as follows: • We showed that the Galactic disc radial velocity field,   (, ), is a strong function of time, due to the relative orientation between the bar and spiral structure.The butterfly pattern in the bar region can thus vary dramatically both in shape and size over periods of a few tens of Myr (see Fig. 4).
• Because of the above, the morphology of the Gaia DR3   (, ) field can be used to constrain the relative orientation between the bar and the inner spiral structure, although this is not straightforward.We found a very good match to the observations for a snapshot from our Model3, for a bar in the process of connecting to a spiral arm.However, identifying the times of all possible matches to the inner disc radial velocity field morphology, for all three models, we concluded that most likely the MW bar is in the process of disconnecting from a spiral, likely the Scutum-Centaurus (see discussion in Sec.4.7).
• The dominating factor distorting the bar's shape and decreasing its position angle with respect to the Sun-GC line, is the heliocentric distance uncertainty (Fig. 5).While this affects the   (, ) morphology, the bar-spiral orientation produces more important variations in both the apparent bar length and the   butterfly pattern (see Fig. 4).
• We require a distance uncertainty of   < 30% in the models to match well the Gaia DR3 data with   < 20%, in order to reproduce the flatness of the transition between negative and positive   , which cannot be achieved with the 20% error in the simulations (see Fig. 7).This may suggest that the data distance uncertainties are underestimated or that the bar angle is 20 • , rather than the nominal value of 30 • .
• We also considered the possibility that the bar angle is smaller than 30 • , which would require smaller distance uncertainty in the simulation to align the transition in the bar butterfly pattern with the Sun-galactocentric line.We found that a bar at a 20 • angle and   < 20% uncertainty can produce a good match to the data, although not as good as the 30 • angle and   < 30% uncertainty (see Fig. 7).More work is needed to explore the interplay between bar angle and distance uncertainties.
• We showed that a range in bar lengths can reproduce the Gaia DR3 radial velocity field (focusing on structure inside the solar circle; see Fig. 8), provided smaller bars are accompanied by stronger spiral structure.Our simulations have bars with genuine lengths of about 5.2, 3.6, and 3 kpc and corresponding spiral structure overdensity of about 10% for Model1 and 20-25% for Model2 and Model3.Considering the MW spirals are expected to have ∼ 20% overdensity, it is tempting to conclude that the MW bar length is consistent with a bar of size below or around 4 kpc.
• We calculated the Fourier components for our three models and found that Model3 has the strongest m=1 and m=3 spirals.This likely results in an additional feature in the   (, ) field (a positive   arm stretching downward from the positive   lobe), which does not exist in the data.This suggests that the MW lacks strong m=1 and/or m=3 modes.
Our conclusion that the MW bar's length is affected by an attached spiral is supported by the work by Rezaei Kh. et al. (2018), which presented an extinction map using red clump and giant stars from the APOGEE survey, showing that the location of the Scutum-Centaurus spiral arm is likely connected to the bar's near side (as shown in their fig.4).Another piece of evidence is the recent work of Lucey et al. (2023), who measured the maximal extent of trapped bar orbits in APOGEE DR17 to extend to ∼ 3.5 kpc, very much consistent with our best-match model (∼ 3.6 kpc).
One should be particularly careful in the interpretation of the velocity field when features are found along a line-of-sight from the solar position.It can be seen in the top-left panel of Fig. 6 that the leading positive   arm, stemming from the positive   lobe associated with the bar, which our simulation so well reproduces with the 30% distance error, is in fact broken at (, ) ≈ (−2, 5) when no uncertainty is added.The 30% distance error, however, causes the gap to disappear and match the data.A hint of this gap is found for distance uncertainty of 20% in the top-left and middleleft panels of Fig. 7, but we need probably less than 10-15% error to identify it unambiguously, as in the bottom panels of the figure.Note that the 10% distance uncertainty cut in the Gaia DR3 data (bottom-rightmost panel) indeed seems to suggest the arm is broken at (, ) ≈ (−4, 5) kpc.
While we found a very good match to the radial velocity field, there are other constraints that should be considered in future work.An obvious one is the tangential component of the velocity,   (, ), which would require to assume a rotation curve in order to subtract the Galactic disc rotation and exhibit the residuals.One can also consider the velocity dispersion, or velocity moments (e.g., Mühlbauer & Dehnen 2003;Hey et al. 2023).
The effect of the beat frequency between the bar and spiral structure needs to be explored, by keeping all other parameters the same.It is feasible that a lower beat frequency (i.e., a slow bar) would result in a stronger effect on the central   (, ) morphology, since fast bars will not have sufficient time for interaction with the spiral.This may be another reason why our Model1, which hosts a bar at the allowed lower limit in terms of the fraction of CR radius to bar length, R =   /  ≈ 1.02 (to be compared to R ≈ 1.75 for Model2's slow bar, see Hilmi et al. 2020), is not producing much variations in the size of the velocity field butterfly pattern, in addition to its weak spiral structure.9. Model3 shows the strongest m=1 and m=3 modes, corresponding to one-armed and three-armed spiral structure.The matches to the data appear to happen near m=1 and/or m=3 maxima, including the best match shown in Fig. 6 (blue vertical).we can conclude that the positive   arm, which stretches downwards in the matching snapshots of Model3 (bottom row of Fig. 8), and which is not seen in the data, results from these odd modes.This would suggest that the MW lacks such strong m=1 or m=3 modes.

Figure 1 .
Figure 1.Gaia DR3 radial velocity field,   ( , ) for different uncertainty cuts, as indicated.Top:   < 10% and    < 5 km/s (top left),   < 20% and    < 5 km s −1 (top middle), and   < 20% and    < 2 km/s (top right).Proper motion uncertainty is   < 10% in all top panels.Bottom: As in top, but with   < 5% (left and middle) and   < 20% (right).Significant difference in the disc area covered and the velocity map morphology are seen only when the distance measurement uncertainty changes.For comparison with our models we use the top left and middle panels.An ellipse with a semi-major axis of 3.5 kpc depicts the bar, oriented at 30 • ahead of the Sun-GC line.The panels with   < 20% reproduce well the top-left panel of Fig. 16 by Gaia Collaboration et al. 2023a, but show a larger range in   and cover a larger disc area.

Figure 2 .
Figure 2. Face-on (top) and edge-on (bottom) stellar density plots for our Model1 (left), Model2 (middle) and Model3 (right panel).X-shaped side-on bar is seen for Model1 and Model3.which have gone through a buckling phase, in contrast to Model2, which has not.

Fig. 1
presents the Gaia DR3 radial velocity field,   (, ) for different uncertainty cuts, as indicated at the top of each panel.The panels with   < 20% (all except for the top left one) reproduce well the top-left panel of Fig. 16 by Gaia Collaboration et al. (2023a), but show a larger range in   (colour bar) and cover a larger disc area.This allows us to see better the positive and negative velocity lobes on the other side of the GC, the emergence of an additional armlike feature in negative velocity (blue) at the upper-left quadrant of the plots, and an area of negative radial velocity in the lower-right corner.An ellipse with a semi-major axis of 3.5 kpc depicts the MW bar in the figure, assumed to be 30 • ahead of the Sun-GC line (e.g., Bland-Hawthorn & Gerhard 2016).

Figure 4 .
Figure 4.Illustrating the changes to the Model3 disc morphology (left column) and its radial velocity field (right column) when the bar is disconnected from the inner spiral structure (top row), in the process of connecting (second row), fully connected and at a maximum measurable length (third row), in the process of disconnecting from spiral (fourth row), and disconnected once again (bottom row).Drastic differences are seen among different panels over these very short time intervals.In all the plots the bar is oriented at 30 • with respect to the Sun-galactocentric line.The Sun's location is indicated by the black dot at  = 8.2 kpc and  = 0.The two dotted circles show the solar radius and 4 kpc.

Figure 5 .
Figure 5.The effect of Gaia DR3-like uncertainties on the bar shape and orientation.The left panel shows the face-on stellar density of Model3 at  = 12.4 Gyr (bar length near a peak, see Fig.3) with the bar oriented at 30 • ahead of the Sun-galactocentric line.The middle and right panels show the effect of 15% and 30% distance uncertainty, respectively.The measurement errors, most importantly the distance uncertainty, cause the bar angle to shift significantly towards the Sun-galactocentric line, as indicated at the bottom of each panel.Moreover, the central contours are affected more, resulting in a distorted bar shape.The white contour marks the same density level, which changes from the 6th, to the 5th, to the 4th most dense contour for a   = 0, 20%, 30%, respectively.This is caused by the central density spread out due to the distance uncertainty increase.We measured 25 • for   < 15% and 12 • for   < 30% taking all density inside the bar length (3.6 kpc) into account.However, it is clear from the innermost couple of contours that the effect is even larger there due to the more circular initial distribution.

Figure 6 .
Figure 6.A snapshot from our Model3, exploring the effect of Gaia DR3 uncertainty and providing a match to the Gaia DR3 radial velocity field.Left column: Model3 radial velocity field (top) and stellar overdensity (bottom).Middle column: As top left, but including 30% (top) and 15% (bottom) synthetic uncertainty in heliocentric distance,   .Right column: Radial velocity field of Gaia DR3 data with   < 20% (top) and   < 10% (bottom).The middle column provides an excellent match to the data.Examining Fig. 3, the time 12.4 Gyr (blue vertical line) corresponds to an increase in the bar's length, thus, the bar can be thought of as being in the process of connecting to a spiral arm.Looking at the morphology of the overdensity plot (bottom-left panel), however, it appears that the bar is in the process of disconnecting from an arm (see Sec. 4.3 for discussion on this conundrum in terms of material arms versus spiral density wave modes).The red arrows in the bottom middle and right panels indicate a feature, where the semi-major axis of the oval cleanly separates positive and negative velocities before it flattens sharply at ∼ 2 kpc from the Galactic centre in both model and data.

Figure 7 .
Figure 7. Exploring the interplay between distance uncertainty and bar angle.The left three columns show the Model3 radial velocity field,   ( , ), for the matching snapshot at 12.4 Gyr (see Fig.6), but with different combinations of distance uncertainty and bar angle, as indicated.When the implemented distance uncertainty is 20% and the bar angle is 30 • (leftmost column), as in the data (rightmost column), an upward kink inside the oval representing the bar is seen for the positive   lobe, which is not present in the data.To achieve the flatness of the transition between negative and positive   , we propose two solutions: decreasing the bar angle from 30 • to 20 • , while setting the distance error at 20% (second column), or using an uncertainty of 30% and a bar angle of 30 • (third column, as in Fig.6).The latter solution appears to give a better match to the data (rightmost column), suggesting that Gaia DR3 distance uncertainties are underestimated.

Figure 8 .
Figure 8.A selection of matches to the Gaia DR3 radial velocity field (shown in bottom left panel) from our three different simulations, by visual inspection, considering different time outputs and focusing mostly on structure inside the solar circle.The top row shows Model1 with a true bar length of ∼ 5.2 kpc, but weaker spiral arms.The second row shows Model2, which has a true size of ∼ 3 kpc (varying from ∼ 2.8 to ∼ 3.2 kpc over the examined time period).The bottom row shows Model3, with a true bar length of ∼ 3.6 kpc, which provided the best match the the Gaia DR3 data in Fig.6(we do not repeat that time output here).For all snapshots we used the same synthetic uncertainties as in Fig.6.It is remarkable that this large range of bar lengths (from ∼ 3 to ∼ 5.2 kpc) among our three models can match relatively well the structure inside the solar circle.We attribute this mostly to the relative overdensity between the bar and spiral arms -the stronger the spiral arms, the smaller the bar that can reproduce the observations.

Figure 9 .
Figure 9. Variation in bar length with time, as in Fig. 3.The three blocks of two rows show Model1 (top), Model2 (middle), and Model3 (bottom).Each block shows the Left Side (top) and Right Side (bottom) of the bar, as indicated.The vertical lines show the times when good matches to the Gaia DR3 data are achieved, some of which were presented in Fig. 8.The higher frequency of oscillations in the top panels reflects the fast bar of Model1, causing it to meet the spiral arms more often.From the number of wave packets seen in Model2 and Model3, it is clear that these bars are both slower than that of Model1, while the Model2 bar is the slowest.

Figure A1 .
Figure A1.Left: Distribution of relative heliocentric distance uncertainties (  /) versus distance, , for the data.Right: Fit to the data on the left, using skewed Gaussian distributions for different distance bins.These are then interpolated to create a continuous probability density function dependent on distance as in the data.Yellow indicates higher density of stars.

Figure A2 .
Figure A2.Match from Model3 to Gaia data, but using the Gaussian distance uncertainties as in the data, shown in Fig. A1.The data are still reproduced well, suggesting that the upward positive   arm may have a break as in the model.In order to reproduce the flatness of the transition between the negative and positive   lobes, a 20 • bar angle was required.

Figure A3 .
Figure A3.Fourier components estimated from the face-on disc density for modes m=1-m=4, as functions of time (see e.g.,Athanassoula & Misiriotis 2002).Curves of different color show three radii in the range 7.2-9.3kpc, as indicated.The vertical lines indicate the times when good matches to the Gaia DR3 data are achieved for the Right side of the bar -same as in the bottom panels of the three blocks found in Fig.9.Model3 shows the strongest m=1 and m=3 modes, corresponding to one-armed and three-armed spiral structure.The matches to the data appear to happen near m=1 and/or m=3 maxima, including the best match shown in Fig.6(blue vertical).we can conclude that the positive   arm, which stretches downwards in the matching snapshots of Model3 (bottom row of Fig.8), and which is not seen in the data, results from these odd modes.This would suggest that the MW lacks such strong m=1 or m=3 modes.