Abstract
We investigate the connection between the vertical structure of stellar discs and the formation of bars using high-resolution simulations of galaxies in isolation and in the cosmological context. In particular, we simulate a suite of isolated galaxy models that have the same Toomre Q parameter and swing amplification parameter but that differ in the vertical scale height and velocity dispersion. We find that the onset of bar formation occurs more slowly in models with thicker discs. Moreover, thicker discs and discs in simulations with larger force softening are found to be more resilient to buckling, which acts to regulate the length and strength of bars. We simulate disc-halo systems in the cosmological environment using a disc-insertion technique developed in a previous paper. In this case, bar formation is driven by the stochastic effects of a triaxial halo and subhalo-disc interactions and the initial growth of bars appears to be relatively insensitive to the thickness of the disc. On the other hand, thin discs in cosmological haloes do appear to be more susceptible to buckling than thick ones and therefore bar strength correlates with disc thickness as in the isolated case. Thus, one can form discs in cosmological simulations with relatively weak bars or no bars at all provided the discs as thin as the discs we observe and the softening length is smaller than the disc scale height.
1 INTRODUCTION
The problem of bar formation in disc galaxies tests our understanding of cosmological structure formation and galactic dynamics. In principle, theories of galaxy formation should yield predictions for the fractional distribution of bars in terms of their strength, length, and pattern speed. While it is often difficult to make precise, quantitative statements about bars from observations, general properties of their distribution have emerged (see Sellwood & Wilkinson 1993; Sellwood 2013; Binney & Tremaine 2008 for reviews). Roughly 25–30 per cent of disc galaxies exhibit strong bars (that is bars that dominate the disc luminosity) and another 20 per cent or more have relatively weak bars (Barazza, Jogee & Marinova 2008; Aguerri, Méndez-Abreu & Corsini 2009; Masters et al. 2010; Sellwood 2013). The bar fraction appears to increase with time. Approximately one-tenth of disc galaxies between 0.5 ≤ z ≤ 2 have visually identifiable strong bars, which is a factor of 3–4 smaller than the fraction in the local Universe (see Simmons et al. 2014 and references therein).
The bar fraction also varies with galaxy type. Masters et al. (2010) find that 70 ± 5 per cent of the so-called passive red spirals have bars as compared to a 25 ± 5 per cent bar fraction for blue spirals. Since the red spirals are interpreted as old galaxies that have used up their star-forming gas, this result is consistent with a bar fraction that increases with time.
Further evidence of a colour correlation with galaxy bar fractions was presented in Masters et al. (2011). They also found an increased bar fraction in galaxies with dominant bulges. Other correlations were explored by Nair & Abraham (2010), who found correlations between bar fractions galaxies of different morphological types and characteristic masses. They also found that central mass concentrations were correlated with different strong bar fractions depending on the morphological type and mass of the galaxies considered. In the same line of thought, Masters et al. (2012) found a correlation between the gas content of a galaxy and the likelihood of a strong bar being present.
Putting all of these observations together, disc galaxies in the local Universe divide into three roughly equal bins: those with strong bars, those with weak bars, and those with no detectable bar. These observations suggest that bars are capable of forming at a wide range of cosmic times, but once formed are difficult to destroy. Though bars can, in theory, be destroyed through an increased central mass concentration, the requisite mass exceeds typical masses for supermassive black holes (Athanassoula, Lambert & Dehnen 2005).
Intuition tells us that properties of a bar should depend on the properties of its host galaxy and the environment in which that galaxy lives. Theoretical arguments indicate that cold, thin discs are susceptible to local ‘Toomre’ instabilities (Binney & Tremaine 2008). Furthermore, discs that are strongly self-gravitating, that is discs that contribute the bulk of the gravitational force required to maintain their rotational motion, are susceptible to global instabilities. Thus, one can construct initially axisymmetric galaxy models that form bars with vastly different properties (or no bars at all) by changing the internal disc dynamics or trading off disc mass for bulge and halo mass. The implication is that the structure of bars provides an indirect means for inferring a disc’s kinematics and mass-to-light ratio as well as the distribution of mass in a galaxy’s dynamically ‘hot’ components, namely its bulge and dark matter halo. However, it is difficult to infer the mass distribution of a galaxy from bar strengths alone since it depends on the halo velocity dispersion, which cannot be observed directly (Athanassoula 2003).
A galaxy’s ability to resist local instabilities is typically expressed in terms of the (kinetic) Toomre Q-parameter (see equation 1) (Toomre 1964) while its ability to resist global instabilities is encapsulated in the swing-amplification X-parameter (see equation 2) (Goldreich & Tremaine 1978, 1979). Both parameters are defined so that large values imply a more stable disc. The hypothesis that they determine a galaxy’s susceptibility to bar formation has been tested by simulations of isolated, idealized galaxy models (Ostriker & Peebles 1973; Zang & Hohl 1978; Combes & Sanders 1981; Sellwood 1981). Typically, the initial galaxy is represented by an N-body (Monte Carlo) realization of an equilibrium solution to the collisionless Boltzmann equation comprising a disc, dark matter halo, and often, a bulge. Equilibrium does not imply stability, and a galaxy can develop spiral structure and a bar through instabilities that are seeded by shot noise (Efstathiou, Lake & Negroponte 1982) and amplified by feedback loops such as swing amplification (Sellwood 2013). A common way to suppress the mechanisms which give rise to these effects is by increasing either Q or X. For example, in dynamically warm discs (that is, discs with high Q) perturbations are randomized on time-scales short enough to prevent the feedback loop from starting (Athanassoula & Sellwood 1986). Likewise, submaximal discs, that is, discs with high X, avoid the bar instability presumably because the disc lacks the self-gravity to drive the bar instability (Efstathiou et al. 1982; Christodoulou, Shlosman & Tohline 1995; Sellwood 2013).
As one might imagine, the parameters Q and X do not uniquely describe a galaxy’s susceptibility to bar formation. Widrow, Pym & Dubinski (2008) present a grid of models in the Q − X plane that all satisfy observational constraints for the Milky Way. These simulations confirm the basic notion that susceptibility to instabilities increases with decreasing Q and X. However, a careful study of bar strength and length as a function of time across these simulations suggests a more complicated picture. In particular, the bar strength is not a perfectly monotonic function of X at fixed Q or vice versa. The implication is that additional parameters are required to fully predict how bar formation will proceed from some prescribed initial conditions. In short, bar formation may proceed very differently within a family of models that have the same Q and X but vary in other ways.
One property of a disc not captured by either Q or X is its thickness, or alternatively, its vertical velocity dispersion. (The Toomre parameter depends only on the radial velocity dispersion.) Klypin et al. (2009) use a suite of simulations to demonstrate that the thickness of the disc plays a profound role in the development of a bar. In particular, their thick disc model forms a stronger and more slowly rotating bar as compared with the bar that forms in a thin disc model with the same initial radial dispersion profile and rotation curve decomposition. Moreover, simulation parameters such as mass resolution and time-step also influence the growth of the bar instability and slowdown of the bar due to angular momentum transfer with the dark halo (Dubinski, Berentzen & Shlosman 2009).
Of course, galaxies are neither isolated nor born as axisymmetric, equilibrium systems. In N-body realizations of these idealized systems, instabilities are seeded from the shot noise of the particle distribution. On the other hand, interactions in real galaxies can be driven by interactions with the given galaxy’s environment. As such, the process of bar formation may be very different in idealized galaxies as compared with galaxies in a cosmological setting. For example, halo substructure in the form of satellite galaxies and dark matter subhaloes can pass through and perturb the disc. Gauthier, Dubinski & Widrow (2006), Kazantzidis et al. (2008), and Dubinski et al. (2009) showed that an apparently stable disc galaxy model can develop a bar when a fraction of the ‘smooth’ halo is replaced by substructure in the form of subhaloes. The effect is stochastic; subhalo-triggered bar formation seems to require subhaloes whose orbits take them into the central regions of the disc in a prograde sense. More recently, Purcell et al. (2011) showed that Sagittarius dwarf alone could have been responsible for the Milky Way’s spiral structure and bar.
Cosmological haloes also possess large-scale time-dependent tidal fields, which impart torques to the disc and cause it to precess, nutate, and warp (Dubinski & Kuijken 1995a; Binney, Jiang & Dutta 1998; Dubinski & Chakrabarty 2009; Bauer, Widrow & Erkal 2018). In turn, stellar discs can reshape the inner parts of the dark matter haloes via adiabatic contraction and dynamical friction (Blumenthal et al. 1986; Ryden & Gunn 1987; Dubinski 1994; Dubinski & Kuijken 1995b; DeBuhr, Ma & White 2012; Yurin & Springel 2015; Bauer et al. 2018). In principle, one can turn to ab initio hydrodynamic cosmological simulations to capture the effects of both triaxiality and substructure. Indeed, galaxy formation studies in the cosmological environment that include the formation of supermassive black holes, stars, and the feedback from these objects on galaxy formation have had remarkable success over the last decade (see for example Vogelsberger et al. 2013; Schaye et al. 2015). Unfortunately, feedback techniques are extremely computationally expensive. Moreover, the simulator cannot control the properties of the disc such as its mass and radial scale length that form within a particular haloes. This restriction makes it difficult to explore the ‘nature versus nurture’ question: Do bars reflect the structure of their host galaxy, substructure interactions of the disc’s lifetime, or large-scale tidal fields from the dark halo?
Techniques developed by Berentzen & Shlosman (2006), DeBuhr et al. (2012), Yurin & Springel (2015), Bauer et al. (2018), and others allow one to insert a collisionless disc into a dark matter halo. These techniques provide a compromise between fully cosmological simulations and simulations of isolated galaxies. In general, the first step is to run a pure dark matter simulation of a cosmological volume and select a halo suitable for the galaxy one intends to study. The simulation is then rerun with a rigid disc potential that is adiabatically grown from some early epoch (say redshift z = 3) and an intermediate one (say z = 1). Doing so allows the halo to respond to ‘disc formation’. At the intermediate redshift, a suitable N-body disc that is approximately in equilibrium is swapped in for the rigid disc potential and the simulation continues, now with live disc and halo particles.
Perhaps the most striking and perplexing result to come from recent disc-insertion simulations is the prevalence of strong bars. Yurin & Springel (2015), for example, find that all of the discs in their bulgeless simulations for strong bars even though some of these discs are decidedly submaximal. In addition, over half of the discs in simulations with classical bulges form strong bars. Their results suggest that discs in the cosmological setting are more prone to forming strong bars and that bulges play an essential role in explaining the presence of disc galaxies with weak bars or no bars at all.
Though the Yurin & Springel (2015) models vary in Q and X they share two important properties. First, the ratio of the radial and vertical velocity dispersion is set to unity throughout the disc. Secondly, the ratio of the vertical and radial scale lengths is fixed at 0.2, which is roughly a factor of two larger than that of the Milky Way’s thin disc. In addition, the gravitational softening length in their simulations is fixed at |$680\, {\rm pc}$|. Thus, the discs in their simulations are thicker and (vertically) warmer than what one might expect for Milky Way-like galaxies. The results from Klypin et al. (2009) suggest that these properties rather than (or together with) halo substructure and triaxiality are responsible for the preponderance of strong bars in the Yurin & Springel (2015) models.
In this paper, we attempt to isolate the different effects that determine whether a galaxy forms a strong bar, a weak one, or no bar at all. The core of the paper is a a sequence of N-body simulations that include simulations of isolating disc galaxies and galaxies in cosmological haloes. Our choice of models is meant to complement those of Yurin & Springel (2015). In particular, we choose models that have the same Q and X but that vary in their vertical structure.
We begin in section Section 2 by discussing the dimensionless parameters that arise when one constructs equilibrium disc galaxy models. These include the aforementioned Q and X parameters as well as the ratio of the vertical and radial velocity dispersion and the ratio of the vertical and radial scale lengths. We also discuss possible effects of gravitational softening. In Section 3, we describe our sequence of simulations and discuss how they fit in with previous work. We discuss results for our isolated galaxy simulations in Section 3.1 and our cosmological ones in Section 5. We conclude with a discussion of the implications of our results in Section 6.
2 THEORETICAL CONSIDERATIONS
In this section, we consider the structural properties of disc-halo models with an eye towards understanding the formation of bars in these systems. We begin with the Q and X parameters and then discuss the vertical structure of the disc as defined by its scale height, vertical velocity dispersion, and surface density. Finally, we consider the effect softening might have on bar formation.
2.1 Q and X
2.2 Vertical structure of stellar discs
As discussed in Klypin et al. (2009), the vertical structure of a stellar disc plays a key role in determining the properties of any bar that forms. In general, the vertical structure is characterized by the vertical velocity dispersion σz, surface density Σ, and scale height. For a self-gravitating plane-symmetric isothermal disc these quantities are connected through the relation |$\sigma _z^2 = \sqrt{12}G\Sigma z_{\rm rms}$|, where zrms is the root mean square distance of ‘stars’ from the mid-plane, and zd is the disc scale height (see Section 4 and Spitzer 1942; Camm 1950).
The dimensionless ratio |$\sigma _{\rm z}^2/\pi G\Sigma z_{\rm d}$| as a function of ρh/ρ0 for the models considered in this paper (stars), the disc-halo models from Yurin & Springel (2015) (filled squares) and the model from Gauthier et al. (2006) (filled triangle). The straight line is the function f = 1 + (2π/3)1/2ρh/ρ0 discussed in the text.
Summary of parameters for the simulations considered in this paper, the disc-halo simulations considered in Yurin & Springel (2015) (labelled YS15) and the Gauthier et al. (2006) (G06). Md is the final disc mass in units of |$10^{10}\, {\rm M}_\odot$|, Rd is the disc scale radius in units of kpc, and Vc and σR are the circular speed and radial velocity dispersion in units of |${\rm km\, s^{-1}}$| and evaluated at Rp = 2.2Rd. For the disc aspect ratio, we quote zd/Rd, where zd is the sech2-scale length. The velocity dispersion ratio σR/σz, the X and Q parameters, the ratio of the halo density in the mid-plane to that of the disc, and the logarithmic derivative of the circular speed are also measured at Rp. Finally, the softening length ϵ is given in units of kpc. Simulations A.III and B.III are the same as A.I and B.I except that they are run with vertical motions isotropized so as to shut off the buckling instability.
| Md | Rd | Vc | σR | zd/Rd | σR/σ|$\rm z$| | X | Q | ρh/ρ0 | ϵ | Description | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| A.I/A.III | 3.49 | 2.50 | 216 | 25.3 | 0.10 | 1.27 | 2.34 | 1.00 | 0.14 | 0.15 | Thin, fiducial/thin, isotropized |
| A.II | 3.49 | 2.50 | 216 | 25.3 | 0.10 | 1.27 | 2.34 | 1.00 | 0.14 | 0.50 | Thin, high softening |
| B.I/B.III | 3.49 | 2.50 | 213 | 25.3 | 0.20 | 0.97 | 2.34 | 1.00 | 0.28 | 0.15 | Mid, fiducial/mid, isotropized |
| B.II | 3.49 | 2.50 | 213 | 25.3 | 0.20 | 0.97 | 2.34 | 1.00 | 0.28 | 0.50 | Mid, high softening |
| C.I | 3.49 | 2.50 | 208 | 25.3 | 0.40 | 0.77 | 2.34 | 1.00 | 0.50 | 0.15 | Thick, fiducial |
| C.I.Ag | 3.49 | 2.50 | 216 | 25.3 | 0.40 | 0.77 | 2.34 | 1.00 | 0.50 | 0.15 | Thick, AGAMA ICs |
| D.I | 5.82 | 3.70 | 245 | 25.2 | 0.10 | 1.14 | 2.45 | 1.00 | 0.29 | 0.18 | Thin, cosmological |
| E.II | 5.82 | 3.70 | 245 | 27.4 | 0.25 | 0.72 | 2.45 | 1.00 | 0.73 | 0.74 | Thick, high softening, cosmological |
| YS15.A5 | 5.00 | 3.00 | 263 | 30.7 | 0.2 | 1.00 | 3.22 | 1.38 | 0.21 | 0.68 | – |
| YS15.B5 | 5.00 | 3.00 | 211 | 26.6 | 0.2 | 1.00 | 2.06 | 0.96 | 0.11 | 0.68 | – |
| YS15.C5 | 5.00 | 3.00 | 270 | 30.3 | 0.2 | 1.00 | 3.31 | 1.42 | 0.23 | 0.68 | – |
| YS15.D5 | 5.00 | 3.00 | 236 | 26.6 | 0.2 | 1.00 | 2.58 | 1.12 | 0.16 | 0.68 | – |
| YS15.E5 | 5.00 | 3.00 | 233 | 27.1 | 0.2 | 1.00 | 2.58 | 1.11 | 0.15 | 0.68 | – |
| YS15.F5 | 5.00 | 3.00 | 219 | 27.0 | 0.2 | 1.00 | 2.22 | 1.02 | 0.11 | 0.68 | – |
| YS15.G5 | 5.00 | 3.00 | 227 | 28.2 | 0.2 | 1.00 | 2.45 | 1.09 | 0.13 | 0.68 | – |
| YS15.H5 | 5.00 | 3.00 | 244 | 28.6 | 0.2 | 1.00 | 2.85 | 1.21 | 0.16 | 0.68 | – |
| G06 | 7.77 | 5.57 | 226 | 17.1 | 0.06 | 1.80 | 2.78 | 1.43 | 0.10 | 0.15 | – |
| Md | Rd | Vc | σR | zd/Rd | σR/σ|$\rm z$| | X | Q | ρh/ρ0 | ϵ | Description | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| A.I/A.III | 3.49 | 2.50 | 216 | 25.3 | 0.10 | 1.27 | 2.34 | 1.00 | 0.14 | 0.15 | Thin, fiducial/thin, isotropized |
| A.II | 3.49 | 2.50 | 216 | 25.3 | 0.10 | 1.27 | 2.34 | 1.00 | 0.14 | 0.50 | Thin, high softening |
| B.I/B.III | 3.49 | 2.50 | 213 | 25.3 | 0.20 | 0.97 | 2.34 | 1.00 | 0.28 | 0.15 | Mid, fiducial/mid, isotropized |
| B.II | 3.49 | 2.50 | 213 | 25.3 | 0.20 | 0.97 | 2.34 | 1.00 | 0.28 | 0.50 | Mid, high softening |
| C.I | 3.49 | 2.50 | 208 | 25.3 | 0.40 | 0.77 | 2.34 | 1.00 | 0.50 | 0.15 | Thick, fiducial |
| C.I.Ag | 3.49 | 2.50 | 216 | 25.3 | 0.40 | 0.77 | 2.34 | 1.00 | 0.50 | 0.15 | Thick, AGAMA ICs |
| D.I | 5.82 | 3.70 | 245 | 25.2 | 0.10 | 1.14 | 2.45 | 1.00 | 0.29 | 0.18 | Thin, cosmological |
| E.II | 5.82 | 3.70 | 245 | 27.4 | 0.25 | 0.72 | 2.45 | 1.00 | 0.73 | 0.74 | Thick, high softening, cosmological |
| YS15.A5 | 5.00 | 3.00 | 263 | 30.7 | 0.2 | 1.00 | 3.22 | 1.38 | 0.21 | 0.68 | – |
| YS15.B5 | 5.00 | 3.00 | 211 | 26.6 | 0.2 | 1.00 | 2.06 | 0.96 | 0.11 | 0.68 | – |
| YS15.C5 | 5.00 | 3.00 | 270 | 30.3 | 0.2 | 1.00 | 3.31 | 1.42 | 0.23 | 0.68 | – |
| YS15.D5 | 5.00 | 3.00 | 236 | 26.6 | 0.2 | 1.00 | 2.58 | 1.12 | 0.16 | 0.68 | – |
| YS15.E5 | 5.00 | 3.00 | 233 | 27.1 | 0.2 | 1.00 | 2.58 | 1.11 | 0.15 | 0.68 | – |
| YS15.F5 | 5.00 | 3.00 | 219 | 27.0 | 0.2 | 1.00 | 2.22 | 1.02 | 0.11 | 0.68 | – |
| YS15.G5 | 5.00 | 3.00 | 227 | 28.2 | 0.2 | 1.00 | 2.45 | 1.09 | 0.13 | 0.68 | – |
| YS15.H5 | 5.00 | 3.00 | 244 | 28.6 | 0.2 | 1.00 | 2.85 | 1.21 | 0.16 | 0.68 | – |
| G06 | 7.77 | 5.57 | 226 | 17.1 | 0.06 | 1.80 | 2.78 | 1.43 | 0.10 | 0.15 | – |
Summary of parameters for the simulations considered in this paper, the disc-halo simulations considered in Yurin & Springel (2015) (labelled YS15) and the Gauthier et al. (2006) (G06). Md is the final disc mass in units of |$10^{10}\, {\rm M}_\odot$|, Rd is the disc scale radius in units of kpc, and Vc and σR are the circular speed and radial velocity dispersion in units of |${\rm km\, s^{-1}}$| and evaluated at Rp = 2.2Rd. For the disc aspect ratio, we quote zd/Rd, where zd is the sech2-scale length. The velocity dispersion ratio σR/σz, the X and Q parameters, the ratio of the halo density in the mid-plane to that of the disc, and the logarithmic derivative of the circular speed are also measured at Rp. Finally, the softening length ϵ is given in units of kpc. Simulations A.III and B.III are the same as A.I and B.I except that they are run with vertical motions isotropized so as to shut off the buckling instability.
| Md | Rd | Vc | σR | zd/Rd | σR/σ|$\rm z$| | X | Q | ρh/ρ0 | ϵ | Description | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| A.I/A.III | 3.49 | 2.50 | 216 | 25.3 | 0.10 | 1.27 | 2.34 | 1.00 | 0.14 | 0.15 | Thin, fiducial/thin, isotropized |
| A.II | 3.49 | 2.50 | 216 | 25.3 | 0.10 | 1.27 | 2.34 | 1.00 | 0.14 | 0.50 | Thin, high softening |
| B.I/B.III | 3.49 | 2.50 | 213 | 25.3 | 0.20 | 0.97 | 2.34 | 1.00 | 0.28 | 0.15 | Mid, fiducial/mid, isotropized |
| B.II | 3.49 | 2.50 | 213 | 25.3 | 0.20 | 0.97 | 2.34 | 1.00 | 0.28 | 0.50 | Mid, high softening |
| C.I | 3.49 | 2.50 | 208 | 25.3 | 0.40 | 0.77 | 2.34 | 1.00 | 0.50 | 0.15 | Thick, fiducial |
| C.I.Ag | 3.49 | 2.50 | 216 | 25.3 | 0.40 | 0.77 | 2.34 | 1.00 | 0.50 | 0.15 | Thick, AGAMA ICs |
| D.I | 5.82 | 3.70 | 245 | 25.2 | 0.10 | 1.14 | 2.45 | 1.00 | 0.29 | 0.18 | Thin, cosmological |
| E.II | 5.82 | 3.70 | 245 | 27.4 | 0.25 | 0.72 | 2.45 | 1.00 | 0.73 | 0.74 | Thick, high softening, cosmological |
| YS15.A5 | 5.00 | 3.00 | 263 | 30.7 | 0.2 | 1.00 | 3.22 | 1.38 | 0.21 | 0.68 | – |
| YS15.B5 | 5.00 | 3.00 | 211 | 26.6 | 0.2 | 1.00 | 2.06 | 0.96 | 0.11 | 0.68 | – |
| YS15.C5 | 5.00 | 3.00 | 270 | 30.3 | 0.2 | 1.00 | 3.31 | 1.42 | 0.23 | 0.68 | – |
| YS15.D5 | 5.00 | 3.00 | 236 | 26.6 | 0.2 | 1.00 | 2.58 | 1.12 | 0.16 | 0.68 | – |
| YS15.E5 | 5.00 | 3.00 | 233 | 27.1 | 0.2 | 1.00 | 2.58 | 1.11 | 0.15 | 0.68 | – |
| YS15.F5 | 5.00 | 3.00 | 219 | 27.0 | 0.2 | 1.00 | 2.22 | 1.02 | 0.11 | 0.68 | – |
| YS15.G5 | 5.00 | 3.00 | 227 | 28.2 | 0.2 | 1.00 | 2.45 | 1.09 | 0.13 | 0.68 | – |
| YS15.H5 | 5.00 | 3.00 | 244 | 28.6 | 0.2 | 1.00 | 2.85 | 1.21 | 0.16 | 0.68 | – |
| G06 | 7.77 | 5.57 | 226 | 17.1 | 0.06 | 1.80 | 2.78 | 1.43 | 0.10 | 0.15 | – |
| Md | Rd | Vc | σR | zd/Rd | σR/σ|$\rm z$| | X | Q | ρh/ρ0 | ϵ | Description | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| A.I/A.III | 3.49 | 2.50 | 216 | 25.3 | 0.10 | 1.27 | 2.34 | 1.00 | 0.14 | 0.15 | Thin, fiducial/thin, isotropized |
| A.II | 3.49 | 2.50 | 216 | 25.3 | 0.10 | 1.27 | 2.34 | 1.00 | 0.14 | 0.50 | Thin, high softening |
| B.I/B.III | 3.49 | 2.50 | 213 | 25.3 | 0.20 | 0.97 | 2.34 | 1.00 | 0.28 | 0.15 | Mid, fiducial/mid, isotropized |
| B.II | 3.49 | 2.50 | 213 | 25.3 | 0.20 | 0.97 | 2.34 | 1.00 | 0.28 | 0.50 | Mid, high softening |
| C.I | 3.49 | 2.50 | 208 | 25.3 | 0.40 | 0.77 | 2.34 | 1.00 | 0.50 | 0.15 | Thick, fiducial |
| C.I.Ag | 3.49 | 2.50 | 216 | 25.3 | 0.40 | 0.77 | 2.34 | 1.00 | 0.50 | 0.15 | Thick, AGAMA ICs |
| D.I | 5.82 | 3.70 | 245 | 25.2 | 0.10 | 1.14 | 2.45 | 1.00 | 0.29 | 0.18 | Thin, cosmological |
| E.II | 5.82 | 3.70 | 245 | 27.4 | 0.25 | 0.72 | 2.45 | 1.00 | 0.73 | 0.74 | Thick, high softening, cosmological |
| YS15.A5 | 5.00 | 3.00 | 263 | 30.7 | 0.2 | 1.00 | 3.22 | 1.38 | 0.21 | 0.68 | – |
| YS15.B5 | 5.00 | 3.00 | 211 | 26.6 | 0.2 | 1.00 | 2.06 | 0.96 | 0.11 | 0.68 | – |
| YS15.C5 | 5.00 | 3.00 | 270 | 30.3 | 0.2 | 1.00 | 3.31 | 1.42 | 0.23 | 0.68 | – |
| YS15.D5 | 5.00 | 3.00 | 236 | 26.6 | 0.2 | 1.00 | 2.58 | 1.12 | 0.16 | 0.68 | – |
| YS15.E5 | 5.00 | 3.00 | 233 | 27.1 | 0.2 | 1.00 | 2.58 | 1.11 | 0.15 | 0.68 | – |
| YS15.F5 | 5.00 | 3.00 | 219 | 27.0 | 0.2 | 1.00 | 2.22 | 1.02 | 0.11 | 0.68 | – |
| YS15.G5 | 5.00 | 3.00 | 227 | 28.2 | 0.2 | 1.00 | 2.45 | 1.09 | 0.13 | 0.68 | – |
| YS15.H5 | 5.00 | 3.00 | 244 | 28.6 | 0.2 | 1.00 | 2.85 | 1.21 | 0.16 | 0.68 | – |
| G06 | 7.77 | 5.57 | 226 | 17.1 | 0.06 | 1.80 | 2.78 | 1.43 | 0.10 | 0.15 | – |
As discussed in the next section, equation (8) holds at the 10 per cent level for our equilibrium models. Departures from equation (8) might come from radial gradients and the rotation of the disc. (see e.g. Read 2014).
2.3 Effect of gravitational softening
Numerical effects can significantly alter the development of bars in simulated galaxies. For example, in simulations of an isolated galaxy that is initially in equilibrium, the onset of bar formation is delayed when mass resolution is increased (Dubinski et al. 2009) essentially because the bar instability is seeded by shot noise. The importance of mass resolution as well as force resolution and time-stepping have been discussed in Klypin et al. (2009).
In this section, we focus on the effects of force softening. Gravitational softening is necessary to avoid large (possibly infinite) forces that arise from Newtonian two-body interactions. Equilibrium models, such as the ones used as initial conditions in isolated galaxy simulations, satisfy the collisionless Boltzmann and Poisson equations. When evolved with force softening, they will begin slightly out of equilibrium. This effect should be most noticeable when the softening length is comparable to or larger than the thickness of the disc. Detailed discussions of the effect softened gravity has on the secular evolution of N-body galaxies can be found in Romeo (1994), Athanassoula & Sellwood (1986), Merritt (1996), Weinberg (1996), and Romeo (1997, 1998).
Softening may have other effects on the development of the bar. In principle, softening should suppress the Toomre instability on small scales. However, this instability develops on scales comparable to or larger than the Jeans length, which is typically much larger than the thickness of the disc and hence larger than the softening length for most simulations. On the other hand, softening may suppress buckling, a bending instability, which is strongest on small scales. As discussed below, buckling appears to be responsible for regulating the growth of bars.
3 MODELS AND SIMULATIONS
3.1 Initial conditions for isolated galaxy simulations
We follow the evolution of isolated disc-halo systems using the N-body code gadget-3 (Springel 2005). The initial conditions for most of our isolated galaxy simulations are generated with GalactICS (Kuijken & Dubinski 1995; Widrow et al. 2008), which allows users to build multicomponent, axisymmetric equilibrium systems with prescribed structural and kinematic properties. Disc particles are sampled from a distribution function (DF) that is a semi-analytic function of the total energy E, the angular momentum about the disc symmetry axis Lz, and the vertical energy |$E_{\rm z} = \Phi (R,0) - \Phi (R,z) + \frac{1}{2} v_{\rm z}^2$|, where Φ is the gravitational potential and v|$\rm z$| is an orbit’s vertical velocity. By design, the disc DF yields a density law in cylindrical |$(R,\, \phi ,\, z)$| coordinates given, to a good approximation, by |$\rho (R,z) = \Sigma (R)\, {\rm sech}^2(z/z_{\rm d})$|. Here, Σ(R) is exponential surface density profile (equation 4) and zd is the scale height. Note that |$z_{\rm rms}= \pi /\sqrt{12} z_{\rm d}$| while the ‘half-mass’ scale height used in Yurin & Springel (2015) is given by z1/2 ≃ 0.549zd ≃ 0.605zrms. The disc DF is also constructed to yield a radial velocity dispersion profile that is exponential in R with scale length 2Rd. The halo DF is designed to yield a truncated NFW profile (Navarro, Frenk & White 1997) as described in Widrow et al. (2008).
While Ez used in the GalactICS disc DF is conserved to a good approximation for nearly circular orbits it varies considerable for stars that make large excursions in R and z. Thus, the initial conditions for ‘thick’ or ‘warm’ discs will not represent true equilibrium solutions to the dynamical equations. To test whether non-conservation of vertical energy affects our results, we compare a thick disc model with GalactICS initial conditions with a similar one where the initial conditions are generated with agama (Vasiliev 2018), an action-based code that does not rely on the epicycle approximation. agama initial conditions should be closer to the true equilibrium than GalactICS models, especially for the thick discs. In principle, this action-based code should yield initial conditions that are closer to a true equilibrium system than ones based on Ez especially for thick discs.
3.2 Description of simulations
In this section, we describe a suite of simulations where Q and X are fixed and where the velocity dispersion and scale length ratios are allowed to vary. Our aim is to test the hypothesis that scale height plays a key role in the development of bars. The parameters for our simulations are summarized in Table 1. Our suite of isolated galaxy simulations form a sequence A, B, C in increasing thickness. The models have the same rotation curve decomposition, which is shown in the top panel of Fig. 2. By design, the contribution to the rotation curve from the disc is slightly below that of the halo at Rp. Therefore our models have X slightly greater than 2 and should be susceptible to global instabilities.
Rotation curve decomposition for our models. Total rotation curves are shown as solid lines while the separate contributions from the disc and halo are shown as dashed and dot–dashed curves, respectively. Blue curves in the top panel are for the isolated galaxy simulations with GalactICS initial conditions while the green curves are for the simulations C.I.Ag run with agama initial conditions. Bottom panel shows initial rotation curve decomposition for the runs D.I and E.II.
The fiducial simulations are run with a softening length of |$184,\, {\rm pc}$|, which is about two-thirds of the scale height of our thinnest model (A.I). The simulations A.II and B.II use a softening length of |$736\, {\rm pc}$|, which is close to the value assumed in Yurin & Springel (2015). The simulation C.I.Ag is similar to C.I (large scale height) but run with agama initial conditions. A comparison of its rotation curve decomposition with that for model C.I is shown in the top panel of Fig. 2. The contributions from the discs in the two models are nearly the same and the contributions from the haloes differ significantly only beyond |$10\, {\rm kpc}$|. The simulations A.III and B.III use a scheme to isotropize vertical motions and effectively shut off buckling and are discussed in Section 4.4.
In addition to these isolated galaxy simulations, we run two cosmological simulations using the disc insertion scheme described in Bauer et al. (2018). The initial conditions for these models, labelled D.I and E.II, are identical except for the vertical scale height and softening length, which are larger in E.II. Thus, these models are cosmological analogues to A.I and B.II. The rotation curves for these models are shown in the bottom panel of Fig. 2. The models themselves are discussed in Section 5.
3.3 Comparison with previous work
While the parameters Q and X allow one to predict the rapidity and vigour with which instabilities develop in disc galaxies that are actually imperfect predictors of the strength and length of bars at late times. The point is illustrated in Widrow et al. (2008) where results for a suite of 25 simulations that explore the Q − X plane are presented. By design, the initial conditions for the models satisfy observational constraints for the Milky Way such as the rotation curve, the local vertical force, and the velocity dispersion toward the bulge. (See Hartmann et al. 2014 for a further analysis of these simulations.) As expected, the onset of the bar instability is delayed in models with large initial values for Q and/or X. However, the dependence on these parameters of the bar strength and length is more complicated. In Fig. 3, we show the distribution of models in the Q − X and zd/Rd − σR/σz planes from a few papers which consider Milky Way like galaxies. It is common for bar formation studies to consider only a small subset of the total parameter space. As we will see, models which are close together in Q and X may exhibit very different bar formation behaviour.
Distribution of simulations considered in this paper in the Q − X and the zd/Rd − σR/σz planes. Stars are simulations run for this paper (A–E); filled squares denote the series of simulations described in Yurin & Springel (2015); the filled triangle denotes the simulation of M31 run in Gauthier et al. (2006); open pentagons denote the simulations described in Widrow et al. (2008).
In Fig. 4, we show the bar strength parameter A2 and length of the bar across the models considered in Widrow et al. (2008). Evidently, the models that form the strongest and longest bars have intermediate values of Q and X. The implication is that models where the instabilities grow too quickly lead to weaker and somewhat shorter bars. Bar formation appears to be a self-regulating process.
Strength and length of bars for the simulations considered in Widrow et al. (2008). The 25 models span the Q–X plane. Top panel shows the A2 parameter while the bottom panel shows the bar length, lb. Both are measured at |$5\, {\rm Gyr}$| (the final snapshot of the simulations). The bar length is the full bar length.
Table 1 gives the relevant parameters for the eight disc-halo models from Yurin & Springel (2015) as well as the disc-bulge-halo model for M31 from Gauthier et al. (2006). In the Yurin & Springel (2015) simulations discs are inserted into dark matter haloes from the cosmological Aquarius simulation. In this respect, they are similar to the disc-insertion simulations described in Section 5. The initial discs in these models all have a scale height to scale length ratio of 0.2 and a radial to vertical velocity dispersion ratio of 1. As discussed above, these choices mean that their discs were chosen from a one-parameter family of models within the Q−X parameter space.
4 ISOLATED GALAXY SIMULATIONS
4.1 Morphology of bar forming galaxies
Face on surface density maps for models A.I, A.II, B.I, B.II, C.I, and C.I.Ag are shown in Fig. 5. All discs form bars by the end of the simulation (|$t=10\, {\rm Gyr}$|). However, bar formation appears to be delayed in models B.I and C.I relative to that in model A.I while the final bar in A.I is shorter than those in B.I and C.I. Other m = 2 features are also evident. These include two-armed spiral structure, most clearly seen in A.I and B.III and elliptical rings, as, for example in C.I.
Surface density maps for isolated galaxy simulations at select times. Time proceeds from 0 to 10 Gyr, left to right, and the models span top to bottom in order of their appearance in Table 1. The overlaid red circles have radii |$R_{\rm p}=5.5\, {\rm kpc}$| and 20 kpc. Recall RP = 2.2Rd, the radius of peak circular velocity.
Evidently, the dominant mode for in-plane perturbations is m = 2 Nevertheless, there are strong m = 3 structures in the |$1.5\, {\rm Gyr}$| snapshot of the A.I and A.II simulations and hints of a weak m = 3 structure in the same snapshot of B.II.
A larger softening length seems to lead to stronger bars at intermediate times. We see this in the comparison of A.I and A.II or B.I and B.II in the 3.0 and |$4.5\, {\rm Gyr}$| snapshots.
4.2 Bar strength parameter A2
Fig. 6 shows a plot of the mean A2 inside the radius Rp as a function of time for the fiducial simulations, the two simulations with high softening, and the thick disc simulation with initial conditions from agama. Consider first the fiducial (low-softening) simulations. Initially, A2 grows roughly exponentially with a growth rate that decreases with increasing thickness. In simulations A.I and B.I, the end of exponential growth is followed by a decrease in A2 after which A2 again increases, now, approximately linearly with time. In the thick disc case (C.I) exponential growth transitions directly to linear growth. The trend is for exponential growth to end at later times as one goes to thicker discs. It is worth noting that the value of A2 at |$10\, {\rm Gyr}$| is similar in the three low-softening simulations.
Mean bar strength parameter inside a cylindrical radius Rp, A2(< Rp), as a function of time. Curves are smoothed in time with a top-hat moving window of width |$1\, {\rm Gyr}$|. Line colours are blue, red, and orange for models A, B, and C, respectively. Results for the fiducial runs A.I, B.I, and C.I are shown as solid curves while the results for the runs with high softening length, A.II and B.II, are shown as dashed curves. The agama model C.I.Ag is shown as a dot–dashed curve.
In the thin disc case, an increase in softening appears to delay the onset of exponential growth as well as the time at which exponential growth ends. Furthermore, the drop in A2 is less severe. Though the value of A2 at the end of the simulation is approximately the same in the low and high softening cases, the bar strength, as measured by A2 is larger in the high-softening case at intermediate times between 4 and |$8\, {\rm Gyr}$|. For the intermediate thickness case (B.I and B.II) softening has little effect on the initial growth rate of A2. But as in the thin disc case, softening allows exponential growth to continue to later times and the final bar is about 20 per cent stronger as compared with the low-softening case. Once again we see that the effect of high softening is to produce stronger bars at intermediate times.
The evolution of A2 for the thick disc runs with GalactICS and agama initial conditions are fairly similar. In particular, the initial growth rate is almost identical as are the final values.
Fig. 6 encapsulates bar strength into a single number, the mean m = 2 Fourier amplitude inside 2.2 disc scale lengths, or about |$5.5\, {\rm kpc}$|. A more complete picture of bar strength is presented in Fig. 7, where we plot A2 as a function of R and t. The figure is constructed by calculating c2 (equation 11) for cylindrical rings of radius |$200\, {\rm pc}$|. Also shown is the corotation radius, which is determined from the pattern speed Ωp and rotation curve. The former is given by a numerical estimate of dϕ2/dt; corotation is found by determining the radius at which Ωp = Vc/R. Thus, since our galaxy models have roughly flat rotation curves beyond |$5\, {\rm kpc}$|, the corotation essentially gives the inverse patter speed or pattern period.
Bar strength parameter A2 as a function of radius and time. The trajectory of corotation is shown by the dashed red line.
From Fig. 7, we see that the corotation radius tends to grow with time and provides an envelope for the bar and other m = 2 structures such as two-armed spirals and elliptical rings. The bar pattern speed is therefore decreasing with time, presumably due to dynamical friction between the bar and both the disc and dark halo (Debattista & Sellwood 1998, 2000b). It is worth noting that the corotation radius increases more rapidly in simulations with high softening. The naive interpretation is that softening somehow increases the frictional coupling between the bar and disc or halo particles. A more likely explanation is that with a high softening length comes stronger bars. Since the acceleration on the bar due to dynamical friction scales as the mass of the bar, stronger bars should spin-down more rapidly.
As in Fig. 6, we see that bar formation is delayed in models with thicker discs. Bar formation is well underway by |$2\, {\rm Gyr}$| in A.I but does not really take hold until |$4\, {\rm Gyr}$| in C.I. Moreover, the first hints of m = 2 power in C.I arise further out at radii closer to |$5\, {\rm kpc}$|.
The dip in bar strength is clearly seen between 2.5 and 3 Gyr in A.I and between 3.5 and 4 Gyr in B.I. As discussed below, we attribute this dip to buckling.
4.3 Vertical structure and velocity dispersion
Fig. 8 shows the zrms radial profiles for a sequence of snapshots in various models. The evolution of zrms is studied during the initial |$500\, {\rm Myr}$| of the simulation. The top panels show the zrms profiles for simulations A.I and A.II and illustrate the effect softening has on the evolution from ‘equilibrium’ initial conditions. As discussed in Section 2, a system that is initialized to be in equilibrium under the assumption of Newtonian gravity will be out of equilibrium if evolved with softened gravity. In particular, the mean potential energy will be systematically low and the system will puff up. For our thin disc model, |$z_{\rm rms}\simeq 230\, {\rm pc}$|. In the high softening case, |$\epsilon = 736\, {\rm pc} \simeq 2.2z_{\rm rms}$|, we estimate the virial ratio for the vertical structure to be 2T/W ≃ 1.7. Of course, the excess kinetic energy will redistribute itself into both kinetic and potential energy. The upshot is that the system quickly settles into a new state with a thickness somewhat larger than the initial one as seen in the right-hand panel. One might be tempted to suggest that a high softening implies that a disc behaves like a thicker disc. This picture is not quite correct, since a thicker disc will not substantially redistribute energy if it is initialized in equilibrium.
Root mean square height zrms as a function of cylindrical radius R for 10 snapshots equally spaced over the first |$500\, {\rm Myr}$|. Panels are for simulations A.I (upper left), A.II (upper right), C.I (lower left), and C.I.Ag (lower right).
The bottom panels in Fig. 8 provide a comparison of zrms profiles for the thick disc simulations with GalactICS and agama initial conditions. We first note that zrms is approximately constant in the C.I but varies by about |$200\, {\rm pc}$| in C.I.Ag. This difference is simply a reflection in how the initial conditions are set-up. In both cases, the scale height depends implicitly on the functional form of the DFs, which are written in terms of either |$E,\, E_{\rm z},\,$| and Lz or the actions. The GalactICS case does exhibit transient wavelike perturbations with a peak to trough amplitude of |$100\, {\rm pc}$| at radii |$R\gt 4\, {\rm kpc}$|. A plausible explanation for these oscillations is that they are due to the fact that Ez is not a true constant of motion. In any case, the system quickly settles to to a new equilibrium state not too different from the initial one.
Fig. 9 shows profiles for zrms and the diagonal components of the velocity dispersion tensor, σz, σR, and σϕ. All models exhibit significant in-plane heating across the disc and throughout the simulations. While the initial radial velocity dispersion profiles are exponential in R the final profiles are relatively flat within the central 5 kpc. Thus, the greatest increase in radial velocity dispersion is at about this radius, which corresponds to the end of the bar. Likewise the σϕ profile develops a ‘bump’ around R ≃ 5–10 kpc.
Diagonal components of the velocity dispersion tensor and zrms as a function of R for different snapshots between 0 and |$10\, {\rm Gyr}$|. Shown, from to to bottom, are profiles for zrms, σz, σR, and σϕ for the same size models included in Fig. 5.
Vertical heating and thickening also appear to be connected with bar formation. The σz profile increases in the A.I/A.II and B.I/B.II simulations such that by 10 Gyr, the central value is 100 km s−1. This value is roughly equal to the initial value in the C.I and C.I.Ag (thick disc) simulations. Note as well that vertical heating in the inner disc occurs rapidly beginning around 3-4 Gyr, roughly when the bar is forming. By contrast, vertical heating of the outer disc occurs gradually over the entire simulation.
The bottom row in Fig. 9 shows the zrms profiles for the different simulations. Again, we see that the inner disc in the A.I/A.II and B.I/B.II simulations begins to thicken around 3–4 Gyr. Moreover, the thickness of the inner disc increases with R reaching a maximum at around R ∼ 5–7 kpc. Curiously enough, the innermost region of the discs in C.I and C.I.Ag decrease with time. As with σz, the thickness of the discs for R ≃ 0 at late times are all approximately the same. Evidently, bar formation tends to drive the innermost regions of the host disc to a common vertical structure.
Fig. 10 provides a different perspective on these results. Here, we plot σR, σϕ, σz, and zrms as a function of time for three representative radii. From the upper three rows, we see that the increase in velocity dispersion is most abrupt in the inner disc, implying that heating there is driven by bar formation. Furthermore, the vertical velocity dispersion in the thick disc models is very nearly time-independent. As in Fig. 9, we see that the discs in all simulations are thickest for |$R \simeq 2.2\, R_{\rm d} \simeq 5{-}6$| kpc. Finally, we again observe that in the thick disc simulations, the inner disc becomes thinner with time.
The velocity dispersions, σR, σϕ, and σz, and zrms as a function of time for our six fiducial models. The lines are averaged in the circular region of radius 0.25 kpc (solid), an annulus of width 0.25 kpc centred near 2.2 Rd (dashed), and a 0.25 kpc width annulus centred at around 5 Rd (dot–dashed).
4.4 Simulations where buckling is suppressed
Buckling is a well-known phenomena often seen in simulations of bar-forming galaxies where the bar bends in and out of the disc plane. Eventually, these coherent oscillations are converted to random vertical motions (Binney & Tremaine 2008). Buckling typically leads to shorter and weaker bars (Martinez-Valpuesta & Shlosman 2004; Debattista et al. 2006)
To isolate the effects of buckling, we implement a simple scheme that prevents the instability from taking hold. Essentially, at each time-step, we reverse the vertical components of the position, velocity, and acceleration for a fraction p of disc particles. In practice, we choose p = 0.25 though the results are insensitive to the exact value.
Fig. 11 shows the effect suppressing buckling has on the disc evolution. In the thin disc case, the bar instability develops a bit faster when buckling is suppressed. More importantly, the drop in A2 seen in simulation A.I is not as strong, thus confirming the notion that buckling regulates the strength of bars. Buckling has a similar effect on our intermediate thickness runs. Furthermore, the effect of suppressing buckling is similar, in some respects, to the effect of increasing softening as can be seen by noting similarities between A.II and B.III. Finally, we note that buckling does not appear to occur in our thick disc simulations.
Mean bar strength parameter inside the cylindrical radius Rp, A2(< Rp), as a function of time. The figure is essentially the same as Fig. 6 though this time we include simulations A.III, B.III, and C.III where buckling is suppressed.
5 COSMOLOGICAL SIMULATIONS
Disc galaxies simulated from axisymmetric, equilibrium initial conditions, as was done in the previous section, form bars at rates and with strengths that depend on their intrinsic scale height of the disc and on the force resolution of the simulation. In this section, we investigate the extent to which these results hold in a cosmological environment. In particular we follow the evolution of a thin disc with moderate softening and a thick disc with high softening that are embedded in identical cosmological haloes.
5.1 Simulation set-up; inserting discs into cosmological haloes
We model a stellar disc in a cosmological halo using the disc insertion scheme described in Bauer et al. (2018). This scheme, which builds on the methods developed by Berentzen & Shlosman (2006), DeBuhr et al. (2012), and Yurin & Springel (2015) uses an iterative procedure to initialize the disc. The first step is to run a pure dark matter simulation and identify a suitable halo. The system is then rerun from redshift zg to zl, this time with a disc potential that grows slowly in mass and radius. Doing so allows the halo particles to respond to the gravitational field of the would-be disc. At zl, the rigid disc is replaced by an N-body system and the ‘live’ disc-halo system is evolved to the present epoch.
For our pure dark matter simulation, we implement the zoom-in technique of Katz et al. (1994) and Navarro, Frenk & White (1994), broadly following the recommendations of Oñorbe et al. (2014), which allows us to achieve very high spatial and mass resolution for a single halo while still accounting for the effects of large-scale tidal fields. We choose cosmological parameters based on the results from Planck 2013 (Planck Collaboration et al. 2014) with |$H_0=67.9\, {\rm km\, s}^{-1}\, {\rm kpc}^{-1}$|, Ωb = 0.0481, Ω0 = 0.306, |$\Omega _\Lambda = 0.694$|, σ8 = 0.827, and ns = 0.962. N-body initial conditions for the dark matter particles are generated with the music code (Hahn & Abel 2013). We select a suitably sized halo for a Milky Way-like galaxy, namely one with a z = 0 mass of |$1.23\times 10^6\, h^{-1} {\rm M}_\odot$| that comprises 106.
During its growth phase from zg = 3 to zl = 1, the disc is treated as a rigid body whose orientation and centre-of-mass position evolve according to the standard equations of rigid body dynamics. At zl, we swap a live disc for the rigid one using the GalactICS code (Kuijken & Gilmore 1989; Widrow et al. 2008), which generates a three-integral DF disc in the best axisymmetric approximation to the halo Bauer et al. (2018).
We run two simulations, D.I, which assumes a thin disc with a softening length of |$184\, {\rm pc}$| and E.II, which assumes a thick disc with a softening length of |$736\, {\rm pc}$|. The softening length chosen for D.I is in accord with the criteria outlined in Power et al. (2003). The simulations D.I and E.II roughly correspond to A.I and B.II, respectively. As well, E.II is similar to the discs considered by DeBuhr et al. (2012), Yurin & Springel (2015) and Bauer et al. (2018), whereas D.I is more consistent with typical discs considered in isolated galaxy suites like Widrow et al. (2008).
5.2 Results
Results from our two cosmological simulations are displayed in Figs 12 and 13. The former shows projections of the mass density at three epochs while the latter gives A2(< Rp) as a function of time. Evidently, the discs in both cases roughly follow the same evolutionary sequence that was seen in the isolated galaxy simulations: rapid growth of the bar strength followed by a period where the bar strength decreases, presumably due to buckling, and finally steady strengthening of the bar. The three epochs chosen in Fig. 12 correspond to the initial growth phase of the bar (|$a=0.6,~t=6.3\, {\rm Gyr}$|), an epoch after buckling (|$a=0.7,~t = 9.2\, {\rm Gyr}$|), and the present epoch at |$t=13.8\, {\rm Gyr}$|. Visually, the bar appears to be stronger and longer in the E.II run than D.I one at each of these epochs but perhaps most notably in the final one. Indeed, the disc in E.II looks very similar to those seen in the simulations of DeBuhr et al. (2012), Yurin & Springel (2015), and Bauer et al. (2018). The fact that the bar in D.I is weaker than the one in E.II is consistent with the results from our isolated galaxy simulations that thicker discs produce stronger bars (see Fig. 5).
Projections for the D.I (left three columns) and E.II (right three columns). The three columns for each simulation correspond,from left to right, to 2.2, 5.9, and |$13.7\, {\rm Gyr}$| after the big bang. The overlaid red circles have radii Rp and 20 |$h^{-1} \,$|kpc.
A2(<Rp) as a function of the age of the Universe for simulations D.I (solid curve) and E.II (dashed curve).
The most significant difference between bar formation in the cosmological setting and bar formation in isolated galaxies concerns the initial growth of the bar. For isolated galaxies, Fig. 6 clearly shows that the onset of bar formation is delayed for thicker discs. Conversely, in the cosmological case, A2 rapidly grows to a value of 0.17 within the first few hundred Myr after the disc ‘goes live’ regardless of the disc thickness. At this point, the bar in the thin disc model decreases in strength with A2 dropping to 0.11 before resuming its growth. By contrast, the bar in the thick disc model continues to grow monotonically. As in the isolated galaxy simulations, self-regulating processes such as buckling are more efficient in the thin disc case and so A2 in simulation D.I lags behind that of E.II. We note that in both cases, A2 drops significantly at around |$t=7.5\, {\rm Gyr}$| and grows steadily thereafter.
Our interpretation of these results is as follows: In isolation, where discs start from axisymmetric initial conditions, the only source of the m = 2 perturbations that drive bar formation is shot noise from the N-body distribution. Evidently, making a disc thicker slows the growth of these perturbations. On the other hand, m = 2 perturbations abound in the cosmological environment where haloes are clumpy and triaxial. The initial growth of the bar may, in fact, be relatively insensitive to the thickness of the disc, once discs are placed in a cosmological setting. On the other hand, disc thickness does effect the resilience of the bar to self-regulating processes, such that buckling and therefore thick discs tend to have stronger bars.
Finally, we note that in both D.I and E.II, a significant number of particles are found at high galactic latitudes. These particles represent stars ‘kicked-up’ from the disc presumably by the large-scale tidal fields of the halo and interactions between the disc and halo substructure. Kicked-up stars have been seen in cosmological simulations by Purcell, Bullock & Kazantzidis (2010), McCarthy et al. (2012) and Tissera et al. (2013). Their existence was inferred in a combined analysis of kinematic and photometric data for the Andromeda galaxy (Dorman et al. 2013). Furthermore, the idea of kicked-up stars has been invoked by (Price-Whelan et al. 2015) to explain the Triangulum–Andromeda stellar clouds (Rocha-Pinto et al. 2003; Martin et al. 2014) and by (; Laporte et al. 2018; Sheffield et al. 2018) to explain the Monoceros Ring (Yanny et al. 2000; Newberg et al. 2002) and associated A13 stellar overdensity (Sharma et al. 2010).
In Fig. 14, we show the fraction of stars with |z| > d for different regions of the discs in our two cosmological simulations. The results are strikingly similar for the two simulations as is already evident from a visual inspection of Fig. 12. The implication is that the processes by which stellar orbits are perturbed out of the disc plane are relatively insensitive to the vertical structure of the disc. We see that very few of the stars with cylindrical radius |$R\lt 15\, {\rm kpc}$| and only |$1-2\, {{\ \rm per\ cent}}$| of the stars between 15 and |$20\, {\rm kpc}$| are kicked-up to distances greater than |$3\, {\rm kpc}$| though some stars from the |$15-20\, {\rm kpc}$| region do end up with |$|z|\gt 10\, {\rm kpc}$|. On the other hand, |$20\, {{\ \rm per\ cent}}$| of the stars from the region beyond |$20\, {\rm kpc}$| end up with |$|z|\gt 3\, {\rm kpc}$| from the mid-plane and |$10\, {\rm kpc}$|. Of course, the actual number of stars is certainly larger since a fraction of the kicked-up stars will be passing through the disc with large vertical velocities.
Fraction of particles with distance from the mid-plane greater than some distance d as a function of d. The difference colours correspond to different bins in cylindrical radius R: |$0\lt R\lt 5\, {\rm kpc}$| – black; |$5\, {\rm kpc} \lt R \lt 10\, {\rm kpc}$| – blue; |$10\, {\rm kpc} \lt R \lt 15\, {\rm kpc}$| – red; |$15\, {\rm kpc} \lt R \lt 20\, {\rm kpc}$| – green; |$20\, {\rm kpc} \lt R\, {\rm kpc}$| – magenta.
6 CONCLUSIONS
We briefly summarize our findings as:
Bar strength in a cosmological setting depends on X and Q, as well as additional properties such as thickness and vertical dispersion.
In isolation, thick discs appear more resilient to buckling.
In a cosmological setting, environment appears to be a much larger driver of bar formation than internal disc dynamics.
Disc scale height has little-to-no impact on the number of stars being kicked out of a stellar disc.
The seminal work of Ostriker & Peebles (1973) introduced the notion that disc dynamics provides a powerful constraint on the structure of discs and the haloes in which they reside. In short, discs that are dynamically cold and that account for a substantial fraction of the gravitational force that keeps their stars on nearly circular orbits are unstable to the formation of strong bars and spiral structure. The existence of galaxies with weak bars or no bars at all tells us that at least some discs are relatively low in mass (i.e. submaximal) and/or dynamically warm.
The theoretical analysis presented in Section 2 showed with a few simple assumptions (e.g. exponential surface density profile) one can derive a relation among the structural parameters of a disc in approximate equilibrium and thus a constraint on initial conditions that one might choose for simulations. For example, if one fixes zd/Rd and σR/σz, as was done in Yurin & Springel (2015), then there is an approximately one-to-one relationship between Q and X. Likewise, fixing Q and X implies a relationship between zd/Rd and σR/σz. These results have important implications for applying disc dynamics as a constraint on models of galaxy formation. In particular, inconsistencies between bar demographics in a galaxy formation model and in observational surveys may reflect differences in the scale height and vertical velocity dispersion of model and real galaxies.
One lesson from our work and the work of others is that the relation between structural parameters of galaxies and bar strength and length is often rather complicated. This observation is no doubt due, at least in part, to the self-regulating nature of bar formation. When bars develop rapidly, they tend to buckle, which leads to weaker and shorter bars (Martinez-Valpuesta & Shlosman 2004). Thick discs appear to be more resilient to buckling, which may explain why bars in these models often end up stronger and longer than bars in thin-disc models (Klypin et al. 2009). For similar reasons, gravitational softening can affect the development and ultimate strength of bars.
In simulations of isolated galaxies from ‘pristine’ equilibrium initial conditions, bar formation is seeded by the shot noise of the N-body distribution. On the other hand, bars in a cosmological environment are subjected to large perturbations including the m = 2 ones that drive bar formation. Thus, the fact that bar formation is delayed in thick disc models of isolated galaxies may be purely academic – bar formation in the cosmological environment will be initiated by a variety of stochastic effects regardless of the thickness of the disc. On the other hand, the resilience of thick discs to buckling is relevant in the cosmological setting and may explain why thick discs tend to form strong bars. The upshot is that a proper understanding the distribution of bars in cosmological models must go hand-in-hand with a proper understanding of the vertical structure of discs.
Clearly, a more exhaustive exploration of the model parameter space is in order. One might, for example, include galaxy scaling relations to further constrain the space of models. In addition, it would be of interest to insert different discs (and for that matter, nearly identical ones) into different haloes in order to explore the random nature of disc–halo interactions. Ultimately, improvements in observations together with a more complete survey of models via simulations should allow us to fully exploit bars in discs as a means of testing and constraining theories of structure formation.
ACKNOWLEDGEMENTS
LMW and JSB are supported by a Discovery Grant with the Natural Sciences and Engineering Research Council of Canada. JSB acknowledges the assistance of Matthew Chequers and Keir Darling in understanding the agama program interface. The authors acknowledge the contributions of the Python scientific open source community, particularly the SciPy project (Jones et al. 2001), NumPy (Oliphant 2015), Pandas (McKinney 2010), and matplotlib (Hunter 2007).
