© ESO 2021

1. Introduction

Bars are common structures in spiral galaxies in the local Universe (Eskridge 2000; Menéndez-Delmestre et al. 2007), and they are able to redistribute angular momentum (Lynden-Bell & Kalnajs 1972) from the inner regions of the disc to the outer disc and the dark matter halo. The distribution functions of the disc and halo determine the amount of material that is available to ‘emit’ and ‘absorb’ this angular momentum (Athanassoula 2003). The density distribution of dark matter in the inner regions of haloes is therefore one of the main factors driving the formation and evolution of bars themselves (Debattista & Sellwood 2000); other important factors include the velocity dispersion of the disc and halo, the central mass concentration, and the gas fraction in the disc (Ostriker & Peebles 1973; Combes & Sanders 1981; Athanassoula et al. 2013).

The inner regions of stellar discs lose angular momentum via the bar, by trapping stars onto elongated orbits – and therefore growing in mass and length over time (Athanassoula 2003) – and through dynamical friction with the dark matter halo (Tremaine & Weinberg 1984b; Weinberg 1985). Dynamical friction can act to reduce the angular rotation speed of the bar, commonly referred to as the bar pattern speed, Ω_p, without increasing the length of the bar, R_bar (Debattista & Sellwood 2000). When bars slow down, the corotation radius, R_CR – the radius at which stars on near circular orbits move around the galaxy with the same angular frequency as the bar – increases. As a result, the corotation radius normalised by the bar length, ℛ = R_CR/R_bar, can shed light on how much dynamical friction the halo exerts on the bar, and therefore on the amount of dark matter in the inner galaxy. As is common in the literature, in what follows we refer to bars with ℛ < 1.4 and ℛ > 1.4 as ‘fast’ and ‘slow’, respectively (Debattista & Sellwood 2000).

Numerous observational efforts have been made to measure the pattern speed of bars, finding that bars tend to rotate fast (i.e., with ℛ < 1.4; Corsini 2011; Aguerri et al. 2015; Guo et al. 2019). This suggests that dark matter haloes do not exert much dynamical friction on bars and are therefore likely subdominant in the central regions of galaxies. However, uncertainties in obtaining the mass-to-light ratio (M/L) of stellar discs make determining the baryon-to-dark matter ratio in the inner regions of galaxies difficult. Dynamical studies of massive spiral galaxies in the local Universe tend to find that they are baryon-dominated in their central regions (Weiner et al. 2001; Kranz et al. 2003; Bovy & Rix 2013; Fragkoudi et al. 2017). This is also found by studies that use the less uncertain M/L obtained from stellar population models in the near-infrared (Lelli et al. 2016). On the other hand, studies such as the DiskMass Survey point to contradictory results, suggesting that dark matter haloes dominate in the central regions of spiral galaxies (Bershady et al. 2011). Therefore, whether or not spiral galaxies are baryon-dominated in their central regions is still under debate.

On the other hand, advances in numerical and physical implementations have led to drastic improvements in hydrodynamical cosmological simulations, which are now able to routinely form spiral galaxies with extended discs (Agertz et al. 2011; Vogelsberger et al. 2014; Schaye et al. 2015; Grand et al. 2017). However, the relative contribution of baryonic and dark matter in the inner regions of galaxies in the Lambda cold dark matter (ΛCDM) framework is an ongoing topic of debate as baryonic processes play an important role in reshaping the inner profiles of dark matter haloes and in setting the baryon-to-dark matter ratio in the central regions (Pontzen & Governato 2012; Lovell et al. 2018). Studying the dynamics of barred galaxies in cosmological simulations provides a powerful tool for constraining the relative amount of baryons and dark matter in the inner regions.

To date, a handful of studies have explored the properties of bars in the full ΛCDM cosmological context (e.g., Kraljic et al. 2012; Scannapieco & Athanassoula 2012; Zana et al. 2018; Rosas-Guevara et al. 2020; Blázquez-Calero et al. 2020). However, few of these have explored their pattern speeds, and only two explore this in a large sample of galaxies (Algorry et al. 2017; Peschken & Łokas 2019) using the EAGLE (Schaye et al. 2015) and Illustris (Vogelsberger et al. 2014) simulations. These studies found that by z = 0, bars have slowed down excessively, making them incompatible with observations. This has resulted in a growing tension between ΛCDM cosmological simulations and fast bars, raising the question of whether fast bars are incompatible with the ΛCDM framework.

Here we revisit the question of the slowing of bars in the ΛCDM context using the Auriga simulations, a suite of state-of-the-art high-resolution magneto-hydrodynamical cosmological zoom-in simulations of the formation of Milky Way-mass galaxies. In Sect. 2 we show that bars in Auriga remain fast until z = 0, which is compatible with observations. In Sect. 3 we compare our results to previous findings in the literature and discuss some of the differences in the simulations that can give rise to the different dynamical behaviour. In Sect. 4 we discuss some of the implications of our findings and conclude.

2. Fast-rotating bars in the Auriga simulations

The Auriga simulations (Grand et al. 2017) are a suite of 30 magneto-hydrodynamical cosmological zoom-in simulations of haloes with masses in the range of 1 × 10¹² − 2 × 10¹² M_⊙. By z = 0, the simulations form star-forming disc galaxies with flat rotation curves that reproduce a range of observed scaling relations, such as the Tully-Fisher relation (Grand et al. 2017) and the size-mass relation of HI gas discs (Marinacci et al. 2017). They also form structures, such as bars and bulges, that have properties compatible with those of observed galaxies (Blázquez-Calero et al. 2020; Fragkoudi et al. 2020; Gargiulo et al. 2019). For this study our barred sample includes all the Auriga galaxies with bar strength > 0.2 at z = 0, where bar strength is defined as the maximum amplitude of the m = 2 Fourier mode of the surface density. We excluded galaxies that are undergoing an interaction at z = 0. More details on the sample of Auriga simulations and how the bar strength is derived are provided in Appendices A and B.

As mentioned above, while the absolute value of the pattern speed of bars can shed light on the exchange of angular momentum in the galaxy, the parameter used to determine the bar slowdown due to dynamical friction is ℛ = R_CR/R_bar, due to the fact that dynamical friction can act to slow down the bar without a corresponding increase in its length. Calculating the length of bars is non-trivial as it oscillates due to the coupling of different structures in the disc, such as the bar with the spiral arms (e.g., Petersen et al. 2019; Hilmi et al. 2020). For the estimate of R_bar, we employed a method similar to that used in analysing observations, namely the ellipse fitting method in Erwin (2005). We verified, however, that our results on ℛ do not change if we opt for a different bar length estimate (see Appendix B for details on the bar length calculation and various tests). The pattern speed in Auriga is calculated from the temporal evolution of the m = 2 Fourier phase of the surface density or via the Tremaine-Weinberg (TW) method (Tremaine & Weinberg 1984a) for simulations without high cadence outputs (see Appendix B for more details).

We plot R_CR versus R_bar in the top panel of Fig. 1 for the barred Auriga galaxies and compare these to observations (Corsini 2011; Aguerri et al. 2015) in which the bar pattern speed is obtained using the TW method. The upper and lower solid lines in Fig. 1 correspond to ℛ = 1.4 and 1, respectively, which denote the regime in which bars are considered fast. For the five Auriga galaxies for which we have high cadence outputs (see Appendix A), we calculated the corotation radius and bar length at redshifts z = 0, 0.25, 0.5, 0.75, and 1 (provided the bar has already formed by the corresponding redshift), as denoted by the colour-coding of the symbols. We find that when bars in the Auriga galaxies are formed they have ℛ < 1.4 (i.e., they are dynamically fast) and that they remain fast throughout their evolution. We note that 25% of the sample have ℛ slightly higher than 1.4 but are compatible with 1.4 within the errors. We therefore find that the Auriga bars are on average fast at z = 0, which is compatible with observations.

Fig. 1. Bar corotation radius vs. bar length for Auriga galaxies. Top panel: corotation radius as a function of bar length for the barred Auriga galaxies, denoted by the symbols in the figure legend. The colour-coding corresponds to measurements at different redshifts, while for the galaxies for which we use the TW method (Au-TW) we estimate the pattern speed at z = 0. These are compared to data from Corsini (2011) and Aguerri et al. (2015) (grey points); the error bars indicate the 1σ uncertainties. The lower and upper dashed grey lines indicate ℛ = 1 and ℛ = 1.4, respectively. Bars in Auriga are fast across all redshifts. Bottom panel: mean ℛ for the aforementioned observations for barred galaxies in Auriga and for the EAGLE and Illustris barred galaxies from Algorry et al. (2017) and Peschken & Łokas (2019) at z = 0. The error bars indicate the 2σ error on the mean. The dashed line indicates ℛ = 1.4, below which bars are considered fast. Bars in the Auriga simulations are fast, compatible with observations, while bars in EAGLE and Illustris are excessively slow at z = 0.

3. Comparison with previous pattern speed estimations from cosmological simulations

In the bottom panel of Fig. 1 we show the mean and 2σ error on the mean of ℛ for barred galaxies in the Auriga simulations, for the above-mentioned observations, and for the barred galaxies studied previously in ΛCDM cosmological simulations (Algorry et al. 2017; Peschken & Łokas 2019); these studies found that bars in the EAGLE and Illustris simulations tend to have ℛ > 2.5, in tension with observations. Bars in Auriga have ℛ < 1.4, thus demonstrating that bars can remain dynamically fast in cosmological simulations within the ΛCDM paradigm. This therefore lessens, to some extent, the previously-reported tension between observed fast bars and ΛCDM.

A natural question to ask is what gives rise to the different behaviour of ℛ in the different simulations. As mentioned above, ℛ depends on the amount of dynamical friction exerted by the halo on the disc, which in turn depends strongly on the amount of dark matter present in the central region of the halo as compared to the disc. The sub-grid physics modelling employed in cosmological simulations can therefore play a critical role in determining the bar properties and, more crucially, its slowdown rate, as it sets the disc-to-halo mass ratio in the central regions of galaxies. In what follows we investigate the effect of numerical resolution, as well as the properties of the galaxies in the different simulations that might lead to differences in ℛ.

Resolution. Low numerical resolution has been extensively discussed in the literature as a parameter that could affect the evolution of bars and the exchange of angular momentum at resonances, with different studies often reaching contradictory conclusions (Weinberg 1998; Weinberg & Katz 2002; Valenzuela & Klypin 2003; Sellwood 2006, 2008).

In order to test the effect of resolution in our models, and to explore whether it could be the primary reason for why previous studies using the EAGLE and Illustris cosmological simulations found large values of ℛ, we carried out a resolution test. We reran two of our Auriga haloes (Au17 and Au18) with eight times lower mass resolution and two times lower spatial resolution. This led to a mass resolution that was closer to, although still slightly higher than, the EAGLE and Illustris resolution. In the low resolution Auriga runs, the stellar and dark matter particles have a mass of 4 × 10⁵ M_⊙ and 3.2 × 10⁶ M_⊙, respectively, while in EAGLE (Illustris) the stellar mass is 1.8 × 10⁶ M_⊙ (1.3 × 10⁶ M_⊙) and the dark matter particle mass is 9.7 × 10⁶ M_⊙ (6.3 × 10⁶ M_⊙). If resolution is the major contributor to the high values of ℛ in previous studies, we would expect to see a significant effect on ℛ when decreasing the resolution to be close to that of EAGLE and Illustris.

In Fig. 2a we show the ℛ values for the level 4 (high resolution) and level 5 (low resolution) runs for Au17 and Au18. We find that ℛ has a slight increase of ∼15% in the low-resolution runs. However, the increase is within the 2σ error of the mean ℛ values for the high-resolution barred galaxies explored in this study. Therefore, this small increase in ℛ for the low-resolution Auriga runs does not appear to be sufficient to explain the much higher values found in the EAGLE and Illustris simulations¹. These considerations suggest that resolution is likely not the main culprit for the high ℛ values found in previous simulations, such as EAGLE and Illustris.

Fig. 2. Comparing properties of discs in Auriga to those in EAGLE and Illustris. (a) Resolution test for ℛ: the value of ℛ for the two lower resolution runs (Au17-5 and Au18-5) compared to the standard resolution runs (Au17-4 and Au18-4). For reference, we also show the mean and 2σ error of ℛ for the whole bar sample of Auriga, EAGLE, and Illustris. (b) Disc thickness: the flattening, ϵ = 1 − c/a, of disc galaxies in Auriga, Illustris, and EAGLE. The average values of the distributions (vertical dashed lines) are 0.42, 0.36, and 0.46, respectively. (c) Gas fraction vs. radius: the average gas density over stellar density in the Auriga galaxies as compared to disc galaxies in EAGLE and Illustris. (d) Gas fraction vs. normalised radius: as in panel c, with the radius normalised by the average bar length in each of the simulations.

Disc thickness. The thickness of discs, which is related to their vertical velocity dispersion, has been reported in the literature as playing a role in the slowing of bars (Misiriotis & Athanassoula 2000; Klypin et al. 2009). While the exact mechanism for this is still unclear, it has been suggested that bars are slower in thicker discs (Klypin et al. 2009). In Fig. 2b we investigate the differences in disc thickness between Auriga, Illustris, and EAGLE by exploring the flattening of their mass distributions, defined as ϵ = 1 − c/a, where c (a) corresponds to the smallest (largest) axis (i.e., ϵ = 0 corresponds to a spheroidal distribution). We obtained the axes of the ellipsoid by calculating the reduced inertia mass tensor of the galaxies (for more details, see Appendix C). We find that the disc galaxies in Auriga, Illustris, and EAGLE have similar distributions of flattening, with Illustris having a slightly lower average value (0.36) than Auriga (0.42) and EAGLE (0.46). While a slight difference in disc thickness might be contributing to the slowing of bars in Illustris, the flattening of discs in Auriga and EAGLE is very similar, and it is therefore unlikely that disc thickness alone can explain the bar slowdown.

Gas fraction. The fraction of gas within discs, and in particular within the bar radius, can play a role in both the formation and strength of bars. It has also been suggested that gas can play a role in maintaining fast bars (Athanassoula 2014), although the topic has given rise to considerable debate (e.g., Bournaud et al. 2005; Berentzen et al. 2007; Athanassoula 2014; Sellwood & Debattista 2014 and references therein). In the two lower panels of Fig. 2 we explore the gas content of disc galaxies in Auriga, EAGLE, and Illustris. Auriga and Illustris have similar gas fractions within 10 kpc, while EAGLE has lower gas fractions, which could possibly contribute to the slower bars in EAGLE as compared to Auriga. However, when the radius is normalised by the average bar length in each simulation, we find that, within the bar radius, Auriga galaxies have gas fractions similar to those of EAGLE and lower than those of Illustris. Therefore, while the gas fraction might contribute to the differences in pattern speed between the different simulations, it is likely not the main culprit for allowing Auriga to maintain fast bars.

Baryon dominance. We therefore turn our attention to the relative contribution of the stellar component compared to dark matter in the Auriga simulations. We explore this in Fig. 3 and compare Auriga to the EAGLE and Illustris cosmological simulations. In the left panel we plot the average and 1σ relation for the stellar mass as a function of the dark matter halo mass (M_⋆ vs. M_h) for the Auriga galaxies (red), as well as for the Illustris and EAGLE discs (green and purple, respectively). We compare these to the commonly used abundance matching relation (Moster et al. 2018; dashed blue line), which is obtained by matching the observed stellar mass function to the halo mass function in cosmological simulations. Auriga galaxies are offset from the relation and lie above both the EAGLE and Illustris discs (i.e., they are globally more baryon-dominated).

Fig. 3. Global and local stellar-to-dark matter ratios in the Auriga galaxies as compared to EAGLE and Illustris. Left: mean (lines) and 1σ (shaded region) values of the M⋆ vs. Mh of Auriga, EAGLE, and Illustris disc galaxies at z = 0 compared to the abundance matching relation from Moster et al. (2018). Middle: efficiency with which gas is turned into stars, i.e.,, $f_{\star }=\frac{M_{\star }}{f_b{M}_{\mathrm{h}}}$ as a function of stellar mass for the Auriga, EAGLE and Illustris galaxies, compared to the relation from abundance matching. The mean and 1σ values for disc (barred) galaxies are shown with the circles/thick lines (squares/thin lines) and shaded region respectively. For comparison, we also show the high-mass disc galaxies from the SPARC sample explored in Marasco et al. (2020), where the error bars denote 16th–84th percentile uncertainties. Right: baryon dominance, V⋆/Vtot, as a function of radius of Auriga, EAGLE, and Illustris disc galaxies in the mass range 3 × 1010 < M⋆/M⊙ < 1011.

In the middle panel of Fig. 3 we plot the efficiency with which – given a universal baryon-to-dark matter ratio, f_b = Ω_b/Ω_c – galaxies convert their gas into stars (i.e.,, f_⋆ = M_⋆/(f_bM_h)) as a function of stellar mass. This is shown for the Auriga, EAGLE, and Illustris simulations, for both the entire disc sample (circles) and for the barred galaxies (squares). The Auriga galaxies are offset from the relation predicted by abundance matching – by more than 2σ at M_⋆ = 9 × 10¹⁰ M_⊙. They also lie above both the EAGLE and Illustris galaxies and thus have higher global stellar-to-dark matter ratios for a given stellar mass at z = 0. Recent work by Posti et al. (2019) explored galaxies in the SPARC sample (Lelli et al. 2016), finding that massive spirals lie above the relation for f_⋆ derived from abundance matching, which implies that massive spirals are overly efficient at converting gas to stars (the so-called ‘failed feedback problem’). A subsequent study by Marasco et al. (2020) found that massive spirals in SPARC are more baryon-dominated than spiral galaxies in cosmological simulations, such as EAGLE and IllustrisTNG100. Interestingly, the Auriga galaxies follow a similar trend as the high-mass spiral galaxies in the SPARC sample, which are denoted by the grey symbols.

While the M_⋆ – M_h relation tells us about the global ratio of baryons to dark matter, the local distribution of baryons to dark matter is more relevant for the dynamics of disc galaxies and for the bar instability. We show this in the right panel of Fig. 3, where we plot V_⋆/V_tot (i.e., the ratio of the stellar component to the total rotation curve) as a function of radius for disc galaxies in Auriga, Illustris, and EAGLE (in the mass range 3 × 10¹⁰ < M_⋆/M_⊙ < 10¹¹). This reveals how ‘maximal’ or baryon-dominated a galaxy is², which in turn will partially determine whether the disc is able to form a bar that remains fast. We find that the Auriga barred galaxies are overall more baryon-dominated than those in EAGLE and Illustris at all radii. In Fig. D.1 we also explore this trend at redshift z = 0.5 and find that Auriga galaxies are already more baryon-dominated at higher redshifts.

What causes Auriga to be more baryon-dominated than the EAGLE and Illustris simulations is the complex interplay of the sub-grid implementations of the baryonic physics in the simulations. A full exploration of this is well beyond the scope of this Letter, but here we outline a few differences that are likely important and refer the reader to the papers describing the simulations for more details (Schaye et al. 2015; Vogelsberger et al. 2014; Grand et al. 2017). In terms of the stellar feedback, the wind model in Auriga has significant differences from its predecessor, Illustris, in both its parametrisation and its implementation, for example, the winds are isotropic in Auriga versus bipolar in Illustris. This might lead to the winds being more effective at higher redshifts in Auriga where galaxies are irregular. On the other hand, EAGLE employs a thermal feedback prescription which is ‘burstier’, heating particles stochastically to sufficiently high temperatures to avoid catastrophic cooling, which is likely more effective at removing baryons from the central regions of galaxies. In terms of the active galactic nucleus (AGN) feedback, the prescription in Auriga is ‘smoother’ than that in Illustris and EAGLE. In Auriga the ‘bubble’ radio-mode feedback provides a gentler and more distributed heating of the circumgalactic medium than in Illustris, in which a smaller number of very energetic bubbles were able to blow out all the gas from the halo; in EAGLE there is only one mode of AGN feedback in which thermal energy is injected stochastically by heating particles around the black hole to very high temperatures. Furthermore, Auriga includes magnetic fields, in contrast to both EAGLE and Illustris. The differences in these physical modelling assumptions, which combine non-linearly, give rise to discs that are more massive in Auriga.

4. Conclusions

We find that bars in the Auriga cosmological simulations are fast, with ℛ < 1.4, in agreement with observations and contrary to the results from previous cosmological simulations. The considerations of the previous section provide us with a likely explanation for why bars in other cosmological simulations tend to slow down so dramatically by z = 0: their host galaxies are embedded in dominant dark matter haloes, both locally and globally. In other words, they have a higher stellar-to-dark matter ratio, and thus their bars suffer more dynamical friction. Disc galaxies in the Auriga simulations, on the other hand, are more baryon-dominated and are thus able to remain fast for over half the age of the Universe. This alleviates the previously reported tension between fast bars and the ΛCDM cosmological paradigm. However, our findings suggest that in order to reproduce the dynamics of barred galaxies, massive spirals must be baryon-dominated and should lie above the abundance matching relation. This is consistent with the recently reported findings (e.g., Posti et al. 2019; Marasco et al. 2020), which highlight the tension between the stellar-to-dark matter ratio of high-mass spirals obtained dynamically, and via methods such as abundance matching.

Our results imply that for cosmological simulations of galaxy formation and evolution to reproduce the dynamical properties of observed barred spirals, they need to build up more stellar mass in their discs than what is predicted by abundance matching. This suggests that care is needed to not make the stellar and AGN feedback prescriptions overly efficient as this can lead to lower stellar-to-dark matter ratios in high-mass spirals, which will prevent them from maintaining fast bars. As bars are present in a significant fraction (≥50%) of the population of spiral galaxies in the local Universe (Eskridge 2000; Menéndez-Delmestre et al. 2007), this effect is likely to considerably increase the expected scatter around the abundance matching relation, particularly for spirals at the high-mass end, where bars are more numerous (Masters et al. 2012; Gavazzi et al. 2015). Whether these results have repercussions on the expected value of the abundance matching relation or simply on its scatter will depend on the relation followed by the global population of galaxies at this mass range, including early type galaxies. In terms of the implications for lower-mass galaxies, there is still considerable debate in the observational community with regards to the frequency of bars in low-mass spirals (Nair & Abraham 2010; Erwin 2018), and, as such, an exploration of the properties of lower-mass barred spirals is timely. Our results highlight the importance of taking the dynamics of barred spiral galaxies into account when constraining models of galaxy formation and evolution.

Acknowledgments

We thank the anonymous referee for a constructive report. FF thanks L. Posti, B. Moster and C. Damiani for useful discussions and N. Peschken for providing the halo IDs of barred galaxies in Illustris. FM acknowledges support through the Program “Rita Levi Montalcini” of the Italian MUR.

References

Appendix A: Barred galaxies in the Auriga simulations

The Auriga simulations (Grand et al. 2017) are a suite of 30 magneto-hydrodynamical cosmological zoom-in simulations of haloes with masses in the range of 1 × 10¹² − 2 × 10¹² M_⊙ that run from redshift z = 127 to z = 0 with cosmological parameters Ω_m = 0.307 and Ω_b = 0.048 and a Hubble constant of H₀ = 67.77 km s⁻¹ Mpc⁻¹ (Planck Collaboration XVI 2014); unless otherwise stated, we use these parameters throughout the Letter. The simulations were performed with the magneto-hydrodynamic code AREPO (Springel 2010; Pakmor et al. 2016) with a comprehensive galaxy formation model (Vogelsberger et al. 2013; Marinacci et al. 2014; Grand et al. 2017), and they form star-forming disc galaxies with flat rotation curves that reproduce a range of observed scaling relations, such as the Tully-Fisher relation (Grand et al. 2017) and the size-mass relation of HI gas discs (Marinacci et al. 2017). They form structures, such as bars and boxy-peanuts, that have properties compatible with those of observed bars (Blázquez-Calero et al. 2020; Fragkoudi et al. 2020) and mainly consist of so-called pseudo bulges (Gargiulo et al. 2019). For more details we refer the reader to the aforementioned papers and references therein.

In this study our ‘barred sample’ includes all the Auriga galaxies with bar strength A₂ > 0.2 at z = 0. This excludes five Auriga galaxies that are undergoing an interaction at z = 0 (Au1, Au11, Au20, Au25, and Au30). For these galaxies the pattern speed cannot be reliably measured using the TW method, for which the continuity equation must hold. We therefore have 16 Auriga galaxies in our barred sample. Five of these haloes (Au9, Au17, Au18, Au26, and Au28) are reruns of the original Auriga haloes with higher cadence outputs (i.e., every 10 Myr), which we included in order to be able to determine the pattern speed from the temporal evolution of the bar and to test our implementation of the TW method (due to the high computational cost, we did not rerun all the Auriga galaxies). For the reruns, the initial conditions of the haloes and the physics implementations were the same as those of the original Auriga haloes. We refer the reader to the following section, where we describe in detail how we obtain the bar strength, length, and pattern speed of our simulated galaxies.

Appendix B: Bar strength, length, and corotation radius

Bar strength. We define the bar strength as the maximum amplitude of the normalised m = 2 Fourier mode of the surface density. We selected the stellar particles in the disc (i.e., with |z|< 1 kpc) in annuli of width 0.5 kpc to calculate

$$\begin{array}{c}\hfill a_m(R)={\displaystyle \sum _{i=0}^N{m}_{i}cos(m{\theta }_i),m=0,1,2,\dots ,}\end{array}$$(B.1)

$$\begin{array}{c}\hfill b_m(R)={\displaystyle \sum _{i=0}^N{m}_{i}sin(m{\theta }_i),m=0,1,2,\dots ,}\end{array}$$(B.2)

where m_i is the mass of particle i, R is the cylindrical radius, N is the total number of particles in that annulus, and θ is the azimuthal angle. To obtain a single value for the bar strength at each time-step, we took the maximum value of the normalised m = 2 component,

$$\begin{array}{c}\hfill A_2=\max \left(\frac{\sqrt{a_2^2(R)+b_2^2(R)}}{a_0(R)}\right),\end{array}$$(B.3)

in radius.

We also carried out visual inspections of the snapshots in order to ensure that the m = 2 mode is due to the bar and not some other spurious short-lived effect. For example, at high redshifts – when mergers are frequent – short-lived m = 2 modes can appear due to the centre of mass of the merging system not being at the centre of the disc (e.g., in the final stages of a merger). This would therefore lead our algorithm to detect a spurious and short-lived m = 2 mode that is unrelated to a bar mode in the disc. However, we note that such spurious and short-lived spikes in m = 2 are essentially only found at high redshifts (see e.g., the large peak in m = 2 for Au26 at ∼9 Gyr in Fig. B.1).

Fig. B.1. Bar properties as a function of time. Top panel: bar strength for the five Auriga galaxies for which we have high cadence outputs as a function of lookback time. The horizontal dashed line indicates A2 = 0.2, above which we consider the bar to have formed. Bottom panel: bar pattern speed for these galaxies once the bar has formed as a function of lookback time.

Pattern speed. For the five galaxies for which we have high cadence outputs, the bar pattern speed was obtained by examining the temporal evolution of the bar phase, that is, by calculating the m = 2 phase from the aforementioned Fourier decompositions in each snapshot and as a function of radius as

$$\begin{array}{c}\hfill \theta (R)=0.5arctan\left(\frac{b_2(R)}{a_2(R)}\right).\end{array}$$(B.4)

The bar pattern speed, Ω_p, was then calculated as an annular average within the bar radius as

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{p}}=\frac{\mathrm{\Delta }⟨\theta ⟩}{\mathrm{\Delta }t}.\end{array}$$(B.5)

The uncertainty on this value mainly derives from the oscillations in pattern speed at different snapshots (see for example the bottom panel of Fig. B.1), which are typically around ±3–5 km⁻¹s⁻¹ kpc⁻¹. We explored the uncertainty this implies in terms of the corotation radius and found that it translates to ∼ ± 0.5 kpc.

Changes in the pattern speed as a function of time can occur due to coupling with spiral arms or odd modes, such as m = 1. To check whether m = 1 can have a significant effect on our estimates of the pattern speed, in the top panel of Fig. B.2 we show the variation in the m = 2 and m = 1 modes as well as the pattern speed as a function of lookback time. We find that at higher redshifts, when gas accretion and minor mergers are more frequent, m = 1 is non-negligible, and this induces some noise in the pattern speed of the bar. However, by low redshifts the contribution from m = 1 is negligible, of the order of a few percent, and as such does not seem to affect our estimates of the bar pattern speed – especially since we excluded the galaxies that are undergoing interactions at z = 0 from our sample.

Fig. B.2. Bar strength, pattern speed, and length variation for Au18 with time. Top panel: pattern speed as estimated from the phase of the m = 2 Fourier mode (right axis – blue), the m = 2 Fourier mode (bar strength) of the surface density, and the m = 1 Fourier mode (dashed red line) as a function of lookback time. Middle panel: various measures of bar length as a function of lookback time (see text). Bottom panel: variation in ℛ as a function of time using different estimates for the bar length.

In Fig. B.1 we show the evolution of the bar strength and pattern speed as a function of lookback time for five of the Auriga galaxies for which we have high cadence outputs. In the top panel of the figure we show the evolution of the bar strength. Bars form at a range of lookback times in our simulations, with some forming at z > 1 (e.g., Au17 and Au18) and some as late as 1 Gyr ago (Au28). Bars in Auriga are long-lived structures (i.e., once the bar has formed it does not dissolve). In the bottom panel of Fig. B.1 we show the bar pattern speed, Ω_p, for these five barred galaxies as a function of lookback time. Most of the bars have pattern speeds that remain roughly constant or are slightly decreasing (e.g., Au26 after t_lookback ∼ 2 Gyr). There is also the curious case of Au17, in which the bar pattern speed increases from the formation time until t_lookback ∼ 3 Gyr due to a resonant interaction with another nearby massive system (a more detailed exploration of this system will be presented elsewhere). The galaxy with the largest decrease in Ω_p, Au26, shows a prototypical example of bar growth in an isolated disc: when the galaxy evolves in isolation, the bar grows gradually stronger over time and Ω_p decreases. In a number of the other cases shown here, bar formation is triggered after a significant merger, with the disc losing a large fraction of its angular momentum by torques induced during the merger (Fragkoudi et al. 2020) with the bar subsequently not growing much stronger or becoming slower.

For the rest of the Auriga barred sample for which we do not have high cadence outputs, we employed the commonly-used TW method (Tremaine & Weinberg 1984a), which uses the continuity equation to estimate the bar pattern speed from a single snapshot at z = 0. We tested our method on the five simulations for which we have high cadence outputs – and therefore a reliable estimate for the bar pattern speed – and found that we can typically recover the pattern speed to within ±5 km⁻¹s⁻¹kpc⁻¹. This translates to an error on the corotation radius of ∼ ± 1 kpc. These are the uncertainties employed in Fig. 1.

Bar length. There is no standard way of estimating the bar length, and various methods have been used in observations and simulations (e.g., visual estimation, isophote ellipse fitting, measurements based on Fourier analysis, etc.), each with its own advantages and disadvantages (for a discussion of the various methods, see Erwin 2005, Petersen et al. 2019, and Hilmi et al. 2020). Here we calculated the bar length using a method commonly used in observations (such as the observations presented in Fig. 1): we fitted ellipses to isodensity contours of face-on images of the galaxies.

We followed a prescription similar to the one used in Erwin (2005), in which the bar length is defined as the minimum of two following ellipse fitting measures: the radius of the first minimum after the maximum of ellipticity or the radius at which the position angle has an extremum outside the bar. An example of this is shown in the left and middle panels of Fig. B.3. If there is no extremum in position angle, then we instead used the radius at which the position angle changes by 10 degrees. For this technique, we did not implement an automated method but rather inspected the outputs on a snapshot-by-snapshot basis. Therefore, this method is not efficient at obtaining bar lengths for a large number of snapshots, but it does give a rather accurate bar length estimate. Other techniques that are more conducive to being automatised exist, such as defining the bar length as the radius at which A₂ drops below a certain percentage of its maximum value or the radius at which there is a certain change in the phase θ of the m = 2 component. These methods, however, can lead to an over- or underestimation of the bar length due to the complex interplay between various modes in the disc; for example, when the spiral arms and the bar align, this can lead to an overestimate of the bar length (see e.g., Hilmi et al. 2020 for a detailed exploration of various estimates of the bar length and the effects of bar-spiral coupling on these estimates), or an underestimation of the bar length can occur when a secondary mode (such as an inner bar) forms with a different pattern speed.

Fig. B.3. Estimating the bar length and corotation radius. Left panel: bar length estimated by fitting ellipses to isodensity contours. The vertical dashed lines denote the estimated bar length. Middle panel: ellipticity and position angle as a function of the semi-major axis of the ellipses. The bar length is the lesser of the two following values: the radius of the first minimum after the ellipticity maximum (dot-dashed orange line) or the radius at which there is the first extremum in position angle outside the bar (dashed line). Right panel: bar corotation radius (marked by the vertical dashed blue line), the radius at which stars in the disc have the same angular frequency (solid black line) as the bar pattern speed (horizontal dashed blue line). The dashed and dot-dashed lines denote the Ω − κ/2 and Ω + κ/2 curves for the galaxy, respectively.

To explore whether the value of ℛ we obtain for the Auriga galaxies in this work would be significantly altered if we used a different measure of the bar length, we explore, in the middle and bottom panels of Fig. B.2, different bar length estimates. We employed two commonly used methods from the literature, namely when the m = 2 Fourier mode drops below 60% of the maximum (orange) and when the phase of the m = 2 Fourier mode changes by more than Δθ_m = 2 = 7.5 degrees (blue). Employing different values for the A₂ threshold or Δθ_m = 2 does not qualitatively affect these results. The estimates we used based on the ellipse fitting method (which is used in the main part of the paper) are shown with black crosses. In the middle panel of Fig. B.2 we see that our ellipse-fitting method for estimating the bar length coincides with the lower values obtained from the other two methods, in agreement with what is estimated as being closer to the ‘true’ bar length in Hilmi et al. (2020). The bar length we use in the main part of the paper therefore gives values on the low side of that obtained from the fluctuations in bar length due to bar-spiral coupling. This in turn leads to higher values of ℛ, as shown in the bottom panel of Fig. B.2. That is to say, for most of the temporal range our method gives consistent or even slightly higher values of ℛ than the other two methods. Thus, we see that using another method to estimate the bar length would not change the results of our study by any significant amount (i.e., we would not find values of ℛ > 1.4) and would thus not change our conclusion on whether Auriga bars are fast or slow. We note that at z ∼ 0 the algorithm that obtains the bar length from the difference in the m = 2 phase (Δθ_m = 2) gives shorter bar estimates than the ones we obtain from our ellipse-fitting method and the maximum of A₂ method. This is due to the growth of an ‘inner bar’-like structure at late times in the galaxy, which has a different pattern speed than the ‘main’ bar of the galaxy. This leads to abrupt phase changes at low radii. However, upon visual inspection it is clear that the primary bar is longer.

Corotation radius. To obtain the corotation radius of Auriga barred galaxies, we first needed to calculate the circular velocity of the galaxy, $V_{\mathrm{c}}^2(r)=r\frac{\partial \mathrm{\Phi }}{\partial r}$. We made the simplifying assumption of spherical symmetry, which reduces the circular velocity to $V_{\mathrm{c}}=\sqrt{GM(r)/r}$, where M(r) is the enclosed mass within a sphere of a radius, r. This simplification will tend to under-predict the true circular velocity, which will be higher for flatter systems. We tested this by assuming a Miyamoto-Nagai disc (Binney & Tremaine 2008) with realistic flattening of the disc, for which we can calculate V_c analytically. We then calculated the approximate value of V_c using the spherical symmetry assumption. We find that the error in V_c is of the order of ∼5−10% depending on the flattening of the disc, with the error increasing for thinner discs. For galaxies in our sample that also have a spherical component due to the dark matter halo, this error will be smaller, roughly of the order of 5%, which is acceptable for the purposes of this study. The angular frequency of stars was then obtained according to Ω = V_c/r, and the epicyclic frequency was calculated as $\kappa =\sqrt{R\frac{\mathrm{d}{\mathrm{\Omega }}^2}{\mathrm{d}R}+4{\mathrm{\Omega }}^2}$ (Binney & Tremaine 2008). The corotation radius is the radius at which stars on nearly circular orbits in the disc have the same angular frequency as the bar. We therefore calculated this as the intersection of the Ω curve and the bar pattern speed, Ω_p (right panel of Fig. B.3).

We note that for the ℛ values of barred galaxies in EAGLE and Illustris shown in the bottom panel of Fig. 1, we used the values reported in Algorry et al. (2017) and Peschken & Łokas (2019), respectively. For the Illustris ℛ, we divided the value reported in Peschken & Łokas (2019) by two since the authors used a proxy for the bar length that, according to them, underestimates the true bar length by a factor of ∼2.

Appendix C: Estimating the galaxy flattening

Galaxies with thinner discs (i.e., with a smaller scaleheight to scalelength ratio) have been found to host fast bars, while models with thicker discs appear to have slower bars (Misiriotis & Athanassoula 2000; Klypin et al. 2009). We therefore explored whether disc thickness could be the main culprit for making bars in Auriga faster as compared to Illustris or EAGLE. We calculated the thickness of the discs in Auriga and compared our results to publicly available values obtained for Illustris and EAGLE (see Genel et al. 2015 and Thob et al. 2019). For Illustris and EAGLE we took only discs galaxies in the mass range 3 × 10¹⁰ < M_⋆/M_⊙ < 10¹¹, the mass range covered by the Auriga galaxies.

To calculate the disc thickness in Auriga, we employed the same method as that in Illustris and EAGLE. In other words, we modelled the galaxies as three-dimensional ellipsoids and calculated their reduced inertial mass tensor,

$$\begin{array}{c}\hfill I_{\mathit{ij}}=\frac{1}{M}{\displaystyle \sum _n{m}_n\frac{r_{n,i}{r}_{n,j}}{r_n^2},}\end{array}$$(C.1)

where i, j = {0, 1, 2} denote the three principal axes, r_n, i is the ith component of the coordinate vector of particle n, m_n is the particle’s mass, and M is the total mass of the galaxy within the radius we are considering. To calculate the reduced mass tensor, we used only stellar particles inside 2×r₅₀, where r₅₀ is the half-mass radius of the galaxy. The axis lengths of the galaxy are then defined as the square root of the eigenvalues λ_i of this matrix. We next compared the flatness of the discs in the different simulations, defined as ϵ = 1 − c/a, where c denotes the smallest axis and a the largest axis of the ellipsoid. The results of this comparison can be seen in Fig. 2b.

Appendix D: The abundance matching relation, M_⋆ and M_h, and maximality

In Figs. 3 and D.1 we use the abundance matching relation from Moster et al. (2018), employing the values for ‘all centrals’ from their Table 8 and using their derived relation for the scatter (their Eq. (25); Moster et al. 2018). The values of M_⋆ and M_h for Auriga, Illustris, and EAGLE used in the aforementioned figures are obtained as follows: for Auriga, we calculated M_⋆ by summing the mass of stellar particles inside 10% of the virial radius of the halo; M_h is the total mass of dark matter particles within the virial radius of the halo. For Illustris, we extracted the M_⋆ within 20 kpc and the total dark matter mass inside the virial radius, M_h, for disc galaxies using the publicly available Illustris data (Nelson et al. 2015). We used the definition of disc galaxies employed by Peschken & Łokas (2019), that is, selecting galaxies that have more than 20% of their stellar mass with a circularity parameter above 0.7 and with flattened distributions, that is, they have a flatness ratio (defined by the ratio eigenvalues of the stellar mass tensor) < 0.7. For the barred galaxies in Illustris, we took those in the ‘total bar’ sample from Peschken & Łokas (2019), where barred galaxies are defined as those disc galaxies with A₂ > 0.15 with some additional criteria (see Sect. 2.2 of their paper for more details). For EAGLE, we used the publicly available data from McAlpine et al. (2016) and, where relevant, the subsequent particle data release (The EAGLE team 2017) to obtain M_⋆ (stellar mass within 20 kpc) and M_h (dark matter mass within the virial radius) values for disc galaxies in EAGLE, which we define as having stellar velocity rotation-to-dispersion ratios (RotToDispRatio) and a stellar disc-to-total ratio from counter rotation (DiscToTotal) larger than 1.7 and 0.7, respectively. For the barred galaxies, we use the IDs of the galaxies identified as barred in Algorry et al. (2017).

Fig. D.1. Stellar-to-dark matter ratios at higher redshifts. Top: f⋆ vs. M⋆ for Auriga, EAGLE, and Illustris disc galaxies at z = 0.5. The dotted blue line denotes the relation determined using the Moster et al. (2018) abundance matching relation, and the shaded region denotes the 1σ scatter around the relation. Bottom: V⋆/Vtot for the galaxies shown in the top panel. Auriga is already locally more baryon-dominated than EAGLE and Illustris at higher redshifts.

To explore the local baryon dominance of barred galaxies in Fig. 3, we followed the definition of disc maximality used in previous works (i.e., we employed the ratio of the circular velocity due to the stellar component as compared to the total circular velocity). We obtained the stellar and total circular velocities as described in the previous section.

All Figures

Fig. 1. Bar corotation radius vs. bar length for Auriga galaxies. Top panel: corotation radius as a function of bar length for the barred Auriga galaxies, denoted by the symbols in the figure legend. The colour-coding corresponds to measurements at different redshifts, while for the galaxies for which we use the TW method (Au-TW) we estimate the pattern speed at z = 0. These are compared to data from Corsini (2011) and Aguerri et al. (2015) (grey points); the error bars indicate the 1σ uncertainties. The lower and upper dashed grey lines indicate ℛ = 1 and ℛ = 1.4, respectively. Bars in Auriga are fast across all redshifts. Bottom panel: mean ℛ for the aforementioned observations for barred galaxies in Auriga and for the EAGLE and Illustris barred galaxies from Algorry et al. (2017) and Peschken & Łokas (2019) at z = 0. The error bars indicate the 2σ error on the mean. The dashed line indicates ℛ = 1.4, below which bars are considered fast. Bars in the Auriga simulations are fast, compatible with observations, while bars in EAGLE and Illustris are excessively slow at z = 0.

In the text

Fig. 2. Comparing properties of discs in Auriga to those in EAGLE and Illustris. (a) Resolution test for ℛ: the value of ℛ for the two lower resolution runs (Au17-5 and Au18-5) compared to the standard resolution runs (Au17-4 and Au18-4). For reference, we also show the mean and 2σ error of ℛ for the whole bar sample of Auriga, EAGLE, and Illustris. (b) Disc thickness: the flattening, ϵ = 1 − c/a, of disc galaxies in Auriga, Illustris, and EAGLE. The average values of the distributions (vertical dashed lines) are 0.42, 0.36, and 0.46, respectively. (c) Gas fraction vs. radius: the average gas density over stellar density in the Auriga galaxies as compared to disc galaxies in EAGLE and Illustris. (d) Gas fraction vs. normalised radius: as in panel c, with the radius normalised by the average bar length in each of the simulations.

In the text

Fig. 3. Global and local stellar-to-dark matter ratios in the Auriga galaxies as compared to EAGLE and Illustris. Left: mean (lines) and 1σ (shaded region) values of the M⋆ vs. Mh of Auriga, EAGLE, and Illustris disc galaxies at z = 0 compared to the abundance matching relation from Moster et al. (2018). Middle: efficiency with which gas is turned into stars, i.e.,, $f_{\star }=\frac{M_{\star }}{f_b{M}_{\mathrm{h}}}$ as a function of stellar mass for the Auriga, EAGLE and Illustris galaxies, compared to the relation from abundance matching. The mean and 1σ values for disc (barred) galaxies are shown with the circles/thick lines (squares/thin lines) and shaded region respectively. For comparison, we also show the high-mass disc galaxies from the SPARC sample explored in Marasco et al. (2020), where the error bars denote 16th–84th percentile uncertainties. Right: baryon dominance, V⋆/Vtot, as a function of radius of Auriga, EAGLE, and Illustris disc galaxies in the mass range 3 × 1010 < M⋆/M⊙ < 1011.

In the text

	Fig. B.1. Bar properties as a function of time. Top panel: bar strength for the five Auriga galaxies for which we have high cadence outputs as a function of lookback time. The horizontal dashed line indicates A2 = 0.2, above which we consider the bar to have formed. Bottom panel: bar pattern speed for these galaxies once the bar has formed as a function of lookback time.
In the text

Fig. B.2. Bar strength, pattern speed, and length variation for Au18 with time. Top panel: pattern speed as estimated from the phase of the m = 2 Fourier mode (right axis – blue), the m = 2 Fourier mode (bar strength) of the surface density, and the m = 1 Fourier mode (dashed red line) as a function of lookback time. Middle panel: various measures of bar length as a function of lookback time (see text). Bottom panel: variation in ℛ as a function of time using different estimates for the bar length.

In the text

Fig. B.3. Estimating the bar length and corotation radius. Left panel: bar length estimated by fitting ellipses to isodensity contours. The vertical dashed lines denote the estimated bar length. Middle panel: ellipticity and position angle as a function of the semi-major axis of the ellipses. The bar length is the lesser of the two following values: the radius of the first minimum after the ellipticity maximum (dot-dashed orange line) or the radius at which there is the first extremum in position angle outside the bar (dashed line). Right panel: bar corotation radius (marked by the vertical dashed blue line), the radius at which stars in the disc have the same angular frequency (solid black line) as the bar pattern speed (horizontal dashed blue line). The dashed and dot-dashed lines denote the Ω − κ/2 and Ω + κ/2 curves for the galaxy, respectively.

In the text

Fig. D.1. Stellar-to-dark matter ratios at higher redshifts. Top: f⋆ vs. M⋆ for Auriga, EAGLE, and Illustris disc galaxies at z = 0.5. The dotted blue line denotes the relation determined using the Moster et al. (2018) abundance matching relation, and the shaded region denotes the 1σ scatter around the relation. Bottom: V⋆/Vtot for the galaxies shown in the top panel. Auriga is already locally more baryon-dominated than EAGLE and Illustris at higher redshifts.

In the text