Drivers of disc tilting I: Correlations and possible drivers for Milky Way analogues

TL;DR

This study uses cosmological simulations to investigate the causes of disc tilting in galaxies, finding that cold gas misalignment, rather than dark matter properties, primarily drives the tilting process.

Contribution

It demonstrates that cold gas misalignment is the main factor influencing disc tilting, challenging the assumed role of dark matter in this process.

Findings

No correlation between dark halo properties and disc tilting.

Strong correlation between gas misalignment and disc tilting.

Gas cooling onto the disc is identified as the principal driver.

Abstract

The direction of the spin vectors of disk galaxies change over time. We present the tilting rate of a sample of galaxies in the NIHAO suite of cosmological hydrodynamical simulations. Galaxies in our sample have been selected to be isolated and to have well determined spins. We compare the tilting rates to the predicted observing limit of Gaia, finding that our entire sample lies above the limit, in agreement with previous work. To test the role of dark matter and of gas we compare the weighted Pearson's correlation coefficients between the tilting rates and various properties. We find no correlation between the dark halo's tilting rate, shape, or misalignment with respect to the disc, and the tilting rate of the stellar disc. Therefore, we argue that, in the presence of gas, the dark halo plays a negligible role in the tilting of the stellar disc. On the other hand, we find a strong correlation between the tilting rate of the stellar disc and the misalignment of the cold gas warp. Adding the stellar mass fraction improves the correlation, while none of the dark matter's properties together with the cold gas misalignment improves the correlation to any significant extent. This implies that the gas cooling onto the disc is the principal driver of disc tilting.

Tables

Table 1. Table 1: All error-weighted Pearson’s correlation coefficients reported in this paper, between the predictor listed and the tilting rate of the stellar disc ( Δ ​ θ ∗ / Δ ​ t Δ subscript 𝜃 Δ 𝑡 \Delta\theta_{*}/\Delta t ). Bold p values highlight correlations with p > 0.7 p 0.7 {\rm p}>0.7 . When logarithmic values are given or the values are plotted in log-space the correlation coefficient is calculated after taking the logarithm and assuming a linear relationship.

Baryonic predictor	units	p(${ϵ}_r$)	p(${ϵ}_{\tau }$)
$M_{\ast }$	M⊙	-0.16	-0.19
$\Delta M_{\max ,\ast }/⟨M_{\ast }⟩\Delta t$	Gyr-1	-0.28	-0.11
$\Delta M_{\max ,b}/⟨M_b⟩\Delta t$	Gyr-1	0.05	0.14
$\Delta {\theta }_{\ast }/\Delta t$	∘ Gyr-1	-0.00	-0.01
$\theta (\ast ,\mathrm{corona})$	∘	0.03	-0.02
$\Delta {\theta }_{\mathrm{corona}}/\Delta t$	∘ Gyr-1	0.12	0.14
$\theta (\ast ,\mathrm{cold}\mathrm{gas})$ ($z=0.3$)	∘	0.82	0.82
$\theta (\ast ,\mathrm{cold}\mathrm{gas})$ ($z=0.15$)	∘	0.86	0.89
$\theta (\ast ,\mathrm{cold}\mathrm{gas})$ ($z=0$)	∘	0.62	0.63
log10 sSFR	Gyr-1	0.08	0.25
log10 $⟨$sSFR$⟩$	Gyr-1	-0.17	-0.01
log10 sSFR($z=0.3$)	Gyr-1	-0.42	-0.34
log10 max(sSFR)	Gyr-1	-0.39	-0.35
$M_{\ast }/M_{\ast ,200}$		0.18	0.09
$M_{\ast }/M_{\mathrm{gas},200}$		0.26	0.18
$M_{\mathrm{gas},\mathrm{disc}}/M_{\ast }$		-0.19	-0.09
$M_{\mathrm{gas},\mathrm{disc}}/M_{\mathrm{gas},\mathrm{total}}$		-0.06	0.07
$\left|L_{\ast }\right|\left(t_0\right)$	$M_{\odot }$ kpc km s-1	-0.44	-0.36
$\Delta \left|L_{\ast }\right|/⟨\left|L_{\ast }\right|⟩$		0.63	0.70
$\left|L_{\mathrm{cold}}\right|\left(t_0\right)$	$M_{\odot }$ kpc km s-1	0.11	0.23
$\Delta \left|L_{\mathrm{cold}}\right|/⟨\left|L_{\mathrm{cold}}\right|⟩$		-0.16	-0.10
Dark matter predictor	units	p(${ϵ}_r$)	p(${ϵ}_{\tau }$)
$b/a$		0.11	0.15
$c/a$		0.08	0.07
$T$		-0.26	-0.37
$|$cos $\theta (\ast ,\mathrm{major})|$	∘	-0.35	-0.27
$|$cos $\theta (\ast ,\mathrm{intermediate})|$	∘	-0.23	-0.52
$|$cos $\theta (\ast ,\mathrm{minor})|$	∘	-0.27	-0.18
$\rho /{\rho }_{\mathrm{crit}}$(R < 4 Mpc)		-0.16	-0.08
$\rho /{\rho }_{\mathrm{crit}}$(R < 6 Mpc)		-0.05	0.21
$\rho /{\rho }_{\mathrm{crit}}$(R < 8 Mpc)		0.38	0.50
	∘ Gyr-1	0.89	0.95
	∘ Gyr-1	0.84	0.87
	∘ Gyr-1	0.34	0.64
	∘	-0.17	-0.36
	∘	-0.10	-0.21
	∘	-0.23	-0.24

Equations

$$\epsilon_{R}(t)=\Bigg{\langle}\phi\Bigg{(}{{\bm{L}}(R,t),\sum_{R}{\bm{L}}}(R,t)\Bigg{)},\;\;\;R=2,3,4,5\Bigg{\rangle}.$$

$$\epsilon_{\tau}(t)=\Bigg{\langle}\phi\Bigg{(}{{\bm{L}}(\tau,t),\sum_{\tau}{\bm{L}}}(\tau,t)\Bigg{)},\;\;\;\tau=\tau_{0,2},\tau_{2,4},\tau_{4,6},\tau_{6,8}\Bigg{\rangle},$$

$$\text{cov}(x,y,\omega )=\frac{{\Sigma }_i{\omega }_i(x_i-⟨x,\omega ⟩))(y_i-⟨y,\omega ⟩)}{{\Sigma }_i{\omega }_i},$$

$$\text{p}=\frac{\text{cov}(x,y,\omega )}{\sqrt{\text{cov}(x,x,\omega )\text{cov}(y,y,\omega )}}.$$

$$I_{ij}=\sum _{\alpha }{m}_{\alpha }{x}_i^{\alpha }{x}_j^{\alpha }/r_{\alpha }^2$$

$$T=\frac{1-(b/a{)}^2}{1-(c/a{)}^2}.$$

$$R^2=1-\frac{{\sum }_i(y_i-{\hat{y}}_i{)}^2}{{\sum }_i(y_i-⟨y⟩{)}^2},$$

$${\bar{R}}^2=1-(1-R^2)\frac{n-1}{n-p-1},$$

$$\frac{\Delta {\theta }_{\ast }}{\Delta t}=3.8\theta (\ast ,\mathrm{cold}\mathrm{gas})+1.6,$$

$$\frac{\Delta {\theta }_{\ast }}{\Delta t}=5.3\theta (\ast ,\mathrm{cold}\mathrm{gas})-3M_{\ast }/M_{\mathrm{gas},200}+3.2,$$

Full text

††pubyear: 2019††pagerange: Drivers of disc tilting I: Correlations and possible drivers for Milky Way analogues–Drivers of disc tilting I: Correlations and possible drivers for Milky Way analogues

1 Introduction

Over the course of its five-year mission, Gaia will provide accurate distances and proper motions for roughly one per cent of the stars in the Milky Way (MW) (Perryman et al., 2001), as well as positions of some $500,000$ quasars, with redshifts between $z=1.5{\text{-}}2$ (Perryman et al., 2001; Claeskens et al., 2006; Lindegren et al., 2008). These quasars will provide an inertial reference frame accurate to $0.2{\text{-}}0.5\mu$as yr*-1*. As a result, Perryman et al. (2014) predicted an accuracy better than $1\mu$as yr*-1* ($0.27^{\circ}$ Gyr*-1*) should be achieved in all inertial spin components of the Gaia reference frame. Binney & May (1986) were the first to propose that discs slew as misaligned angular momentum is accreted by a galaxy. Ostriker & Binney (1989) attributed the formation of warps to disc slewing. Moreover dark haloes also tumble, which may also drive tilting of the disc. If the tilting rates of the stellar discs are similar to those of dark haloes (Bailin & Steinmetz, 2004; Bryan & Cress, 2007), such a tilt would be detectable in the MW. Earp et al. (2017) (hereafter, E17) presented the tilting rates for galaxies comparable to the Milky Way ($9\times 10^{11}\>{\rm M_{\odot}}\leq M_{200}\leq 1.2\times 10^{12}\>{\rm M_{\odot}}$) in a $\Lambda$ Cold Dark Matter cosmological hydrodynamical simulation. E17 measured the tilting rates between $z=0.3$ and [math], excluding mergers, finding that all galaxies exhibited significant tilting, at an average rate of $4.9^{\circ}\pm 2.7^{\circ}$Gyr*-1*, well above Gaia’s predicted detection limit.

Various processes can cause discs to tilt. Minor mergers and small-scale interactions can result in a change in angular momentum (Ostriker & Tremaine, 1975; Toomre, 1981; Huang & Carlberg, 1997; Sellwood et al., 1998; Benson et al., 2004; Read et al., 2008; Kazantzidis et al., 2009; Frings et al., 2017). As satellites fall into the disc they tilt towards alignment with the disc due to the transfer of angular momentum (Huang & Carlberg, 1997; Read et al., 2008). The most massive interaction of the MW is with the Large Magellanic Cloud (LMC), with mass estimates as high as $M_{\mathrm{LMC}}\sim 2\times 10^{11}\>{\rm M_{\odot}}$ (Kallivayalil et al., 2013; Gómez et al., 2015; Peñarrubia et al., 2016), corresponding to $20$ per cent of the mass of the MW. If the LMC is this massive it would put it above the upper mass estimate for the initial mass of the Sagittarius dwarf galaxy (Jiang & Binney, 2000), meaning it is the most massive interaction in some time. Because of the uncertainty in the LMC’s mass, its importance on the tilting of the MW’s disc cannot yet be determined.

Sustaining star formation in galaxies requires ongoing gas accretion onto the disc. As gas falls into the dark matter’s potential well it is shock heated to the halo virial temperature $T_{\rm vir}=10^{6}(\nu_{\rm circ}/167\text{km s}^{-1})^{2}$K, forming a hot gas corona (Spitzer, 1956; Rees & Ostriker, 1977; Silk, 1977; White & Rees, 1978; Savage & de Boer, 1979; White & Frenk, 1991; Dahlem, 1997; Wang et al., 2001; Fukugita & Peebles, 2006). This hot gas cools and eventually settles into the disc (Fall & Efstathiou, 1980; Brook et al., 2004; Kereš et al., 2005; Robertson et al., 2006; Brooks et al., 2009). Gas cooling from the corona contributes angular momentum to the disc. Cosmological simulations have shown that the coronae of MW progenitor galaxies form early and thereafter their angular momentum evolves differently from that of the dark matter (e.g. Obreja et al., 2019; Jiang et al., 2019). The angular momentum of the corona is usually misaligned with that of their stellar disc (van den Bosch et al., 2002; Roškar et al., 2010; Velliscig et al., 2015; Stevens et al., 2017). This addition of misaligned angular momentum results in the disc tilting. Debattista et al. (2015) (hereafter, D15) showed that under these circumstances, the orientation of the disc’s angular momentum is determined by a balance between the torques from the triaxial dark halo, and the net inflow of angular momentum via cooling gas. As a result, star-forming galaxies where gas is continually cooling onto the disc are generally misaligned with the principal axes of their dark haloes, as has been found in large surveys (Sales & Lambas, 2004; Brainerd, 2005; Agustsson & Brainerd, 2006; Yang et al., 2006; Azzaro et al., 2007; Faltenbacher et al., 2007; Wang et al., 2008; Wang et al., 2010; Nierenberg et al., 2011; Li et al., 2013).

Dubinski (1992) measured the tilting rates of dark haloes in the mass range $(1-2)\times 10^{12}\>{\rm M_{\odot}}$; he found that they rotated uniformly with rotation rates in the range $6^{\circ}-96^{\circ}$Gyr*-1*. Similarly Bailin & Steinmetz (2004) found smooth figure rotation with an average tilting rate of $6.2^{\circ}$Gyr*-1*. Bryan & Cress (2007) found that $63$ haloes exhibited an average pattern speed of $13.8^{\circ}\;h\;$Gyr*-1*. To investigate coupling between the orientation of discs and dark haloes, Yurin & Springel (2015) inserted live stellar discs into eight, MW-sized, high-resolution dark haloes from the aquarius simulation. They found tilting rates of $5^{\circ}-6^{\circ}$Gyr*-1*, comparable to the halo tilting rates measured in pure N-body simulations.

Bett & Frenk (2012) examined the consequences of minor mergers and flybys on the spin of dark haloes. They measured the angular momentum of haloes with mass $12.0\leq\log_{10}(M/\>{\rm M_{\odot}})h^{-1}\leq 12.5$ at $z=0$. They found that such events only caused small changes to the angular momentum of the entire halo. However, the inner halo ($R<0.25R_{200}$), which Binney et al. (1998) show is very tightly coupled gravitationally to the stellar disc, was not so stationary, with 47 per cent experiencing a change in their angular momentum orientation of at least $45^{\circ}$ during their lifetimes.

Therefore, there is a good theoretical basis to expect the MW’s disc to be tilting, with several possible driving mechanisms. In this paper, we investigate the role of the dark matter and of gas on the tilting rates of stellar discs from the NIHAO suite of cosmological hydrodynamical zoom-in simulations. The paper is organized as follows. Section 2 presents the suite of simulations used in this paper. Section 3 details the selection of a subsample of galaxies that have well-determined tilting rates. Section 4 compares the tilting rates to various possible predictors. Section 5 presents our conclusions.

2 Numerical Simulation

For this paper we used the NIHAO (Numerical Investigation of Hundred Astrophysical Objects)111’nihao’ is the Chinese word for ’hello’ simulations suite. NIHAO is a sample of $\approx 100$ hydrodynamical cosmological zoom-in simulations performed using the SPH code gasoline2 (Wadsley et al., 2017). The code includes gas heating via ultraviolet (UV) heating and ionization and cooling due to hydrogen, helium and metals (Shen et al., 2010). The star formation and feedback modelling follows what was used in the MaGICC simulations (Stinson et al., 2013), adopting a threshold for star formation of $n_{th}>10.3$ cm*−3*. Stars can feed energy back into the interstellar medium via SN feedback (Stinson et al., 2006) and via ionizing radiation from massive stars (early stellar feedback) before they turn into SN (Stinson et al., 2013). We refer the reader to Wang et al. (2015) for a more detailed description of the code and the simulations.

The NIHAO simulations are the largest set of cosmological zoom-in simulations covering the halo mass range $10^{10}$ - $10^{12}\>{\rm M_{\odot}}$, they combine very high spatial and mass resolution with a statistical sample of haloes. NIHAO has proven very successful in reproducing several key properties of observed galaxies including their cold gas masses and sizes (Stinson et al., 2015; Macciò et al., 2016), the stellar and baryonic Tully-Fisher relations (Dutton et al., 2017) and stellar disc kinematics (Obreja et al., 2016). A key property of the NIHAO galaxies is that they lie on the stellar mass-halo mass relation across their full mass range. This is shown for the galaxies used in this work in Figure 1, where the virial mass versus the stellar mass within 10 per cent of $R_{200}$ for NIHAO galaxies is compared to the observed relationship from Kravtsov et al. (2018), derived using halo abundance matching (black line). The green star on Figure 1 denotes the observed MW values, with virial mass $M_{200}\sim 1.1\times 10^{12}\>{\rm M_{\odot}}$ and stellar mass $M_{*}\sim 5\times 10^{10}\>{\rm M_{\odot}}$ (Bland-Hawthorn & Gerhard, 2016). Having galaxies with realistic stellar content at all masses allows us to extend the analysis performed in E17, and hence to consider the tilting rates of lower mass galaxies.

3 Sample selection

3.1 Tilting rates

To measure the angular momentum of our galaxies, it is vital that their centres are correctly determined. We use two different methods to determine the centre of each galaxy: the centre of mass within a shrinking sphere, and the lowest potential dark matter particle. The shrinking sphere method follows the iterative technique of Power et al. (2003): at each step, the centre of the sphere is determined as the centre of mass of the previous step, and the radius of the sphere is reduced by 2.5 per cent. This is iterated down to the smallest sphere containing at least 100 particles, at which point the centre of mass of this final sphere is returned. For the second technique, we simply look for the dark matter particle with the lowest potential energy and use its position as the centre. As in E17, we find that the lowest potential dark matter particle gives more reliable centres for the sample, due to some galaxies having high-density substructure away from the true centre. After removing galaxies where the centres were still poorly determined, the sample is reduced to 85 galaxies from 91. We are predominantly interested in MW mass galaxies, therefore, we impose a lower limit on the stellar mass of $M_{*}>10^{9}\>{\rm M_{\odot}}$, leaving us with 30 galaxies.

The angular momentum of the stellar disc is calculated using only stars with ages $\tau<1$ Gyr within a range of radii, from 0.5 per cent of $R_{200}$, out to 10 per cent of $R_{200}$, in 0.5 per cent increments. We adopt the angular momentum of the young stars within 5 per cent of $R_{200}$ ($\sim 10$ kpc for a MW sized galaxy) for the stellar disc. This is motived by warps often being present at radii, $R>0.05R_{200}$. We determine the tilting rates by measuring the angle, $\Delta\theta$, between the angular momentum vectors of the stellar disc at $z=0.3$ and at [math], then dividing by the time difference ($\sim 3.7$ Gyr).

We have verified, by visual inspection, that the evolution of the tilt is uniform over 20 steps between $z=0.3$ and $z=0$, indicating that the tilting is coherent. Thus the integrated tilt is not a spurious effect of random noise. We will present examples of this evolution in a followup paper.

3.2 Isolation criterion

Following E17, we are primarily interested in galaxies evolving in relative isolation. In such galaxies, the change in stellar mass should be driven primarily by in-situ star formation. Therefore, we compare the maximum fractional change in stellar mass $\Delta M_{*}(t_{\rm peak})/(\langle M_{*}\rangle\Delta t)$ to the specific star formation rate (sSFR), given by the ratio of the star formation rate to the stellar disc mass, at the same time step, $t_{\rm peak}$. We further refine our high mass galaxies by excluding any galaxy that has gained more than 50 per cent of its stellar mass via accretion at any time step; leaving us with 26 galaxies. We then divide these 26 galaxies into two mass subsamples, high mass (6 galaxies) and low mass (20 galaxies), with $M_{200}=10^{12}\>{\rm M_{\odot}}$ being the separator. Figure 2 shows the resulting distribution in the space of mass growth versus sSFR. Galaxies in the shaded region will have gained a majority of stellar mass through accretion and are excluded. The (black) diamonds denote the low mass galaxies and the (red) squares show the high mass galaxies. We estimate the mass lost due to supernova feedback by multiplying the rate for type-II supernovae in the MW ($1.9\pm 1.1$ per century Diehl et al. (2006)) by the integral of the Kroupa IMF (Kroupa, 2002) for high mass stars ($M_{*}=8-40\>{\rm M_{\odot}}$) derived by the integral for all mass. Combining this feedback rate with the Galaxy’s sSFR rate obtained from the current star formation rate $1.65\pm 0.19\>{\rm M_{\odot}}$ yr*-1* and the current stellar mass, $M_{*}\sim 5\times 10^{10}\>{\rm M_{\odot}}$ (Bland-Hawthorn & Gerhard, 2016), we estimate the stellar mass change of the MW. This implies that $\sim 30$ per cent of the mass gained from star formation is lost due to supernovae feedback. The (green) dashed line in Fig. 2 applies this offset to the one-to-one relationship.

3.3 Error measurements

The spin direction errors were obtained using two different methods. First, we calculate the angular momentum vectors of the young stellar disc within 2, 3, 4 and 5 per cent of $R_{200}$. We then measure the angle between each of these four vectors, $\phi$, and their sum. The mean of these values is then taken as our error,

[TABLE]

We take the root mean squared sum of $\epsilon_{R}(t)$ at $z=0.3$ and [math], resulting in our first error $\epsilon_{R}$. We use only stars with ages $\tau<1$ Gyr for this calculation, as these are the least contaminated by the bulge and stellar halo; because of this we must ensure that the young stars are closely aligned to the rest of the disc. Therefore, we employ a second method by calculating the angular momentum vectors of star particles with ages in 2 Gyr bins from 0 to 8 Gyr. As before, we measure the angle between each of these four vectors, $\phi$, and their sum. We then take the mean value,

[TABLE]

where $\tau_{n,m}$ are the stars with ages between $n$ and $m$ Gyr. We take the root mean squared sum of $\epsilon_{\tau}(t)$ at $z=0.3$ and [math], resulting in $\epsilon_{\tau}$. We construct our final subsample of galaxies by imposing an upper limit on both $\epsilon_{R}$ and $\epsilon_{\tau}$ of $\epsilon\leq 5^{\circ}$. The upper limit on $\epsilon_{R}$ reduces our subsample down to 21 galaxies and the upper limit on $\epsilon_{\tau}$ reduces our subsample down to 12 galaxies (7 low mass and 5 high mass, which we refer to as the high mass and low mass samples, respectively, hereafter) that have well-determined tilting rates. We have verified by visual inspection that the galaxies removed by these cuts are poor analogues of the Milky Way, generally because they are either elliptical or interacting.

4 Drivers of tilting

We find error-weighted averages for the tilting rate of the entire low$+$high mass subsample of $3.8\pm 2.3^{\circ}$ Gyr*-1* and $3.6\pm 2.4^{\circ}$ Gyr*-1* for $\epsilon_{R}$ and $\epsilon_{\tau}$, respectively. For just the high mass galaxies we find error-weighted averages of $3.8\pm 2.8^{\circ}$ Gyr*-1* and $3.7\pm 2.7^{\circ}$ Gyr*-1* for $\epsilon_{R}$ and $\epsilon_{\tau}$, respectively.

4.1 Pearson’s correlation coefficient

Throughout this section we report the error-weighted Pearson correlation coefficient, $\text{p}(x,y,\omega)$ between the tilting rate of the stellar disc and each possible predictor. The weights are defined as $\omega=\epsilon^{-2}$, where $\epsilon$ is the error on the tilting rate. With this set of weights we determine error-weighted covariance as

[TABLE]

where $\langle x,\omega\rangle$ denotes the weighted mean of $x$. The error weighted Pearson correlation coefficient is given by

[TABLE]

We report p values for the entire low mass$+$high mass subsample, due to the small number of galaxies. Correlation coefficients of $|{\rm p}|<0.4$ will be treated as null results, values between $0.4<|{\rm p}|<0.7$ will be referred to as weak correlations and values of $|{\rm p}|>0.7$ will be referred to as strong correlations. Table 1 provides a full list of all correlation coefficients calculated using both error methods, $\epsilon_{R}$ and $\epsilon_{\tau}$, with strong correlations indicated in bold; however, in the text we will refer only to p values calculated using $\epsilon_{R}$.

4.2 The minimal role of dark matter

In this section we compare the properties of the dark halo to the tilting rates of the stellar discs they harbour, to investigate the extent of the halo’s influence.

Figure 3 presents the tilting rates of the stellar disc and the dark matter’s angular momentum. The stellar disc’s angular momentum is measured within $0.05R_{200}$ for all three panels, whereas, the dark matter’s angular momentum is measured within $0.05R_{200}$ (left), $0.1R_{200}$ (middle) and $R_{200}$ (right). The distribution roughly follows the one-to-one relationship (black dashed line) within $R<0.1R_{200}$, with a very strong correlation for both the left and middle panels. For $R_{200}$ (right panel), this correlation vanishes, suggesting that the tilting of the entire dark halo is not related to that of the disc. Instead, just the inner part of the halo, which is directly affected by the disc, follows the disc closely.

As the shape of the halo becomes less spherical the torques imposed will increase, all other things being equal. Therefore, if the dark halo is the primary driver of the tilting stellar disc, one might expect a correlation between the shape of the halo and the tilting rate. To measure the shape of the halo, we define the principal axes of the dark halo with the condition that $a>b>c$ and, following Kazantzidis et al. (2004), the principal axis ratios are given by $s=b/a$ and $q=c/a$. Using an iterative process, the shape of the dark halo is calculated, starting with a spherical ellipsoid. The modified inertia tensor $I_{ij}$ is defined as:

[TABLE]

where $x_{i}^{\alpha}$ is the $i$th coordinate of the $\alpha$th particle and $r_{\alpha}$ is the elliptical radius defined as $r_{\alpha}^{2}=x^{2}_{\alpha}+y^{2}_{\alpha}/s^{2}+z_{\alpha}^{2}/q^{2}$. The eigenvalues of the modified inertia tensor are used as the new values of $s$ and $q$. These iterations continue until the values of $s$ and $q$ converge to a fractional difference less than $10^{-2}$. Following Franx et al. (1991) we measure the triaxiality as

[TABLE]

$T$ values of $1$, [math] and $\sim 0.5$, correspond to prolate, oblate and triaxial spheroids, respectively. Figure 4 shows the tilting rate of the disc with respect to the resulting intermediate to major axis ratios ($b/a$, left), minor to major axes ratios ($c/a$, middle) and the triaxiality parameter ($T$, right), measured at $z=0$. Although previous authors have claimed that they are able to cause stellar discs to tilt purely by the gravitational dynamical interactions between the disc and halo, in the case of these hydrodynamical simulations we find no correlation between the shape of the halo and the tilting rate of the disc. We have also verified this with the halo shape measured at $z=0.3$, again finding no correlations.

D15 showed that red galaxies tended to be aligned such that their disc angular momentum was parallel to the minor axis of their parent dark halo, whereas blue galaxies tended to have random orientations. To test if stellar discs tilt towards alignment with one of the principal axes of the dark halo, Figure 5 shows the distance between the angular momentum vector of the stellar disc and the major- (left), intermediate- (middle) and minor-axis (right); we find no correlations for any of the alignments, neither measured at $z=0.3$ nor at $z=0$.

Figure 6 shows the distribution of alignments between stellar discs and their host haloes. The right panel shows that the majority of the galaxies are not closely aligned to any of the dark halo’s principal axes. Almost all the galaxies in the low mass subsample are star-forming blue galaxies, with sSFRs generally above the level of the MW. This is in line with the result of D15 that blue galaxies generally have random alignments with respect to their dark halo. We also considered the alignment between the hot ($T>50000$K) gas corona and the dark halo’s angular momenta, again finding no tendency to align. Therefore, as the angular momentum of the corona is a product of the ongoing gas flow from the surrounding IGM and feedback, we argue that there is no preference for the angular momenta of the inflow to be aligned with that of the dark halo.

4.3 Local cosmic overdensity

Another possible driver of disc tilting is the proximity of high-density regions. E17 found that galaxies in higher density regions tend to tilt at higher rates and this correlation strengthened for larger volumes. Figure 7 compares the tilting rate to the local density calculated within volumes of radii between 4 and 8 Mpc. We find no correlation between the tilting rate of the stellar disc and the local cosmic overdensity out to any radius, at $z=0.3$ or at $z=0$. This is in contrast to the results of E17. However, the galaxies in high-density regions in E17 were more likely to be in cluster/group environments, which might explain this apparent discrepancy, whereas the NIHAO sample was specifically designed for more isolated galaxies. Indeed the range of cosmic overdensities in E17 spanned 0.7 to 10.1, with an average and median of 2.5 and 1.7, respectively, whereas for our NIHAO subsample the range is much smaller, 0.4 to 1.6, with average and median values of 0.85 and 0.78, respectively. In comparison, the local overdensity for the local group is roughly unity (Klypin et al., 2003; Karachentsev, 2005).

4.4 The impact of gas

In hydrodynamical simulations, the stellar disc and hot gas corona are often misaligned (e.g. van den Bosch et al., 2002; Velliscig et al., 2015; Stevens et al., 2017). This misalignment results in gas cooling from the corona onto the disc with misaligned angular momentum, forming warps (Roškar et al., 2010), and directly contributing misaligned angular momentum to the disc. As in E17, we find that all our discs have angular momentum misaligned with that of their hot gas coronae. Following E17, Figure 8 compares both the angular misalignment between the hot ($T>50000$K) gas corona and the stellar disc (left), as well as the tilting rate of the hot gas corona to the tilting rate of the disc (right). We do not reproduce the weak correlation between the angular misalignment of the hot gas corona and disc found in E17.

Gas falling onto the hot gas corona from the intergalactic medium contributes misaligned angular momentum. As a result, the net angular momentum of the hot gas corona tilts. This gas then cools from the corona and contributes its misaligned angular momentum to the disc. Therefore, the two tilting rates may be correlated. The right panel of Figure 8 compares the two tilting rates; as in E17, we find no correlation.

Roškar et al. (2010) showed that the angular momentum of the warp gives a good indication of the angular momentum of the corona. As such, the misalignment of the warp and the stellar disc should shed light on the angular momentum that is being added to the disc from the corona, or directly through cold flows. As our galaxies vary in disc scale-length and size we measure the angular momentum of all the gas out to $0.1R_{200}$ to determine the angular momentum of the warp. Figure 9 compares the angle between the angular momentum vectors of the cold gas ($T<20000$K) and the stellar disc versus the tilting rate of the stellar disc, at $z=0.3$ (left), $z=0.15$ (middle) and $z=0$ (right). We find a strong correlation between the cold gas misalignment with the tilting rate for $z=0.3$ and $z=0.15$, and a slightly weaker correlation at $z=0$. A comparable analysis for cool gas ($T<50000$K) gives similar results.

4.5 Is the star formation rate a proxy for the tilting rate?

Observationally detecting gas accretion or the alignment of the hot gas corona is extremely difficult. However, as the galaxy requires this ongoing accretion of gas to fuel its star formation, the star formation rate itself is a proxy for the amount of gas reaching the disc. We, therefore, measure the star formation rate (SFR) using the mass of all star particles born between subsequent time steps. Since our subsample spans a wide range of masses, we use the specific star formation rate (sSFR), by dividing the SFR by the stellar mass within $R<0.1R_{200}$, at the subsequent time step. Figure 10 compares the sSFR at $z=0$ and averaged over the time range to the tilting rate of the stellar disc. The green dashed line shows a recent estimate for the MW’s current and averaged sSFR. The left panel takes the value of $1.65\pm 0.19\>{\rm M_{\odot}}$ yr*-1* for the present day MW SFR, with the current stellar mass assumed to be $M_{*}\sim 5\times 10^{10}\>{\rm M_{\odot}}$ (Bland-Hawthorn & Gerhard, 2016). We find no correlation between either the present day sSFR or the average sSFR and the tilting rate. We do find a weak correlation between the tilting rate and the sSFR at $z=0.3$ as well as with the peak sSFR.

4.6 Baryonic mass fractions

The left panel of Figure 11 compares the ratio of stellar mass out to $0.1R_{200}$ and total stellar mass out to $R_{200}$ to the tilting rate of the stellar disc. The green line in this figure denotes an upper limit for the MW, assuming the LMC is the dominant stellar mass outside $0.1R_{200}$ and has a virial mass of $M_{\rm LMC}\sim 2\times 10^{11}\>{\rm M_{\odot}}$ (Gómez et al., 2015; Kallivayalil et al., 2013; Peñarrubia et al., 2016). Determining the stellar mass from $M_{200}$ was done by assuming the LMC follows an observed abundance matching relation (Kravtsov et al., 2018). A value of $M_{*}\sim 5\times 10^{10}\>{\rm M_{\odot}}$ was assumed for the stellar mass of the MW (Bland-Hawthorn & Gerhard, 2016). We find no correlation between the stellar mass ratio and the tilting rate of the stellar disc, indicating that the presence of nearby satellites within $R_{200}$ is not the main driver of the tilting rate; therefore, in the MW it is not likely that the LMC would be responsible for the tilting of the disc. The right panel of Figure 11 compares the ratio of stellar (out to $0.1R_{200}$) to total gas mass (out to $R_{200}$) against the tilting rate of the stellar disc. We find no correlation between either the mass of the hot gas corona or the sSFR and the tilting rate. Moreover, we find no correlations between the tilting and either the stellar mass or the fractional change in baryonic/stellar mass.

4.7 Angular Momentum

Figure 12 shows the tilting rate of the stellar disc versus the angular momentum of the stellar disc (top) at $z=0.3$ (left), and its difference between $z=0.3$ and $z=0$ normalized by the mean angular momentum (right). A rough estimate for the angular momentum of the Milky Way’s stellar disc is $|L|=M_{d}R_{d}V_{c}=3.1\times 10^{13}\>{\rm M_{\odot}}$ km s*-1* kpc (green dashed line), obtained by assuming a disc stellar mass to be $5\times 10^{10}\>{\rm M_{\odot}}$ (Bland-Hawthorn & Gerhard, 2016), all at the scale radius $R=2.6$ kpc (Bland-Hawthorn & Gerhard, 2016), and with a circular velocity of $240$ km s*-1* Schönrich (2012); Bland-Hawthorn & Gerhard (2016). The bottom panels show the angular momentum of the cold gas disc ($T<20000$K) at $z=0.3$ (left), and its change normalized by the average angular momentum over the same time period (right), versus the tilting rate of the stellar disc. We find a weak anti-correlation between the amount of angular momentum in the stellar disc and the tilting rate, meaning galaxies with higher angular momentum tilt more slowly. Moreover, we find a weak correlation between the normalized change in angular momentum of the stellar disc and its tilting rate. This hints at a connection between the amount of angular momentum being gained by the disc and its tilting rate. Therefore, the amount of angular momentum able to reach the disc might be an indicator of how fast the disc can tilt. In the remaining two panels, we find no correlation.

4.8 Multiple regression: double variable models

To measure the statistical significance of including multiple variables we use weighted least squares (WLS) linear regression to determine the adjusted coefficient of determination ($\bar{R}^{2}$). The coefficient of determination is given by,

[TABLE]

where $y_{i}$ is the observed value for each galaxy and $\hat{y}_{i}$ is the predicted value for each galaxy. A value of $R^{2}=1$ would mean the model can explain all variability in the data, whereas, $R^{2}=0$ would mean the model fails to explain any of the variability. When adding more parameters the value of $R^{2}$ can only remain the same or increase, therefore, we turn to the adjusted coefficient of determination, given as

[TABLE]

where $p$ is the number of predictors, and $n$ is the number of observations. If we include an additional variable that does not improve the model, $R^{2}$ will stay the same, but $\bar{R}^{2}$ will decrease. The addition of a statistically significant variable will increase both $R^{2}$ and $\bar{R}^{2}$. As $p$ becomes bigger, the change in $R^{2}$ needed for the variable to be significant will also increase. The WLS linear model is constructed with the statsmodels package in Python, using $1/\epsilon_{R}^{2}$ for our weights. To find the maximum value for $\bar{R}^{2}$ we use the forward selection form of stepwise regression. For this method, we measure $\bar{R}^{2}$ for a model with just a constant term, then measure $\bar{R}^{2}$ after the addition of each variable. The variable that statistically improves the model the most is chosen as our first parameter. We then measure $\bar{R}^{2}$ for a model with our first parameter and the inclusion of a second parameter from the remaining variables. This process is repeated until the model no longer statistically improves or until we reach a specified number of parameters. We consider all the baryonic predictors listed in Table 1 finding a best fitting, single parameter, linear model for the tilting rate of the stellar disc given by

[TABLE]

where $\theta(*,{\rm cold\;gas})$ is the angle between the stellar disc and cold gas warp measured at $z=0.15$, with $R^{2}=0.74$ and $\bar{R}^{2}=0.71$. This result is unsurprising, as the correlation between $\theta(*,{\rm cold\;gas})$ and the tilting rate exhibits the largest p-value. We compare this best fitting model to the best fitting model with two parameters,

[TABLE]

where $M_{*}/M_{{\rm gas,200}}$ is the ratio of disc stellar mass to virial gas mass. Equation 2 has a larger $R^{2}$ of $0.897$ and is a statistically significant improvement over Equation 1 with $\bar{R}^{2}=0.874$. Therefore, a better fitting linear model uses two parameters $\theta(*,{\rm cold\;gas})$ and $M_{*}/M_{{\rm gas,200}}$.

Figure 13 shows the residuals between the expected values for the tilting rates using Equation 1 and the measured tilting rates against the ratio of disc stellar mass to virial gas mass. We see a negative trend which, due to the increase in $\bar{R}^{2}$ between Equations 1 and 2, is statistically significant.

We used the same method on the dark matter predictors listed in Table 1. We find no model that can predict the tilting rate better than $R^{2}=0.47$ and $\bar{R}^{2}=0.35$. Therefore, no model using dark matter predictors is able to match the tilt rate as well as models using gas predictors.

5 Discussion and Conclusions

We have presented the tilting rates for galaxies in a suite of hydrodynamical cosmological simulations and tested for correlations with various properties. As in E17, we find that all our galaxies tilt with rates that would be detectable by Gaia, assuming a detection limit of $0.27^{\circ}$ Gyr*-1* (Perryman et al., 2014). We find error-weighted average tilting rates for our subsample of $3.8\pm 2.3^{\circ}$ Gyr*-1*. If we restrict our subsample to high mass galaxies, we find error-weighted average tilting rates of $3.8\pm 2.8^{\circ}$ Gyr*-1*.

5.1 The lack of a dark matter influence

Using cosmological pure N-body simulations Yurin & Springel (2015) reported tilting rates for their stellar discs greater than the tilting rates presented here. In such a scenario the tilting of the stellar disc is driven by the torques, dynamical friction and interactions imposed by the tilting host DM halo and surrounding substructure. To test whether torquing by the dark halo drives the tilting of the stellar disc, we compare the tilting rate of the stellar disc to the axis ratios, $b/a$ and $c/a$ as well as the triaxiality of the halo and find no correlations. We also compare the tilting rate of the dark halo’s angular momentum to that of the stellar disc, finding strong correlations for the tilting rates within $0.1R_{200}$. However, at such low radii, the halo’s angular momentum is dominated by the stellar disc. Therefore, it is more plausible that the tilting of the dark matter is driven by the tilting of the stellar disc, not vice versa. Moreover, we compared the tilting rate of the disc to the misalignment with respect to each axis of dark halo, finding no correlations.

We investigate the correlation found in E17 that galaxies in denser regions tilted with higher rates. We find no such correlation between the local overdensity and the tilting rate at any radius between 4 and 8 Mpc. This difference is most likely due to selection criteria imposed on the NIHAO sample, which were absent in the galaxies presented in E17. E17 did not choose galaxies to be preferentially in the field, resulting in some galaxies residing in cluster/group environments with overdensities as large as $\rho/\rho_{\rm crit}\sim 10$, whereas the largest values for NIHAO are $\sim 2$. Dense environments could greatly affect the tilting of galaxies due to the higher frequency of tidal interactions.

Finally, we tried to combine the dark matter predictors into a single model using multiple regression, but were unable to create a model that can predict the tilting rate with an $\bar{R}^{2}>0.35$, far worse than models using gas predictors. Combining these five results, it becomes hard to argue that the torques imposed by the dark matter are the primary drivers for disc tilting in the presence of gas.

Although the majority of galaxies in NIHAO are star forming, one of the galaxies in our sample has sSFR lower than the MW. This galaxy also has the lowest tilting rate, $0.68\pm 0.23^{\circ}$ Gyr*-1*, and its stellar disc is aligned to within $8^{\circ}$ of the minor axis of its dark halo at $z=0$. This makes its disc the most aligned with the minor axis of any galaxy in our sample, with the lowest sSFR rate of $0.022$ Gyr*-1*. Therefore, in the absence of gas, the dark halo can begin to drive disc tilting, as D15 predicted.

5.2 The strong influence of gas

The angle between the angular momentum vectors of the gas warp and the stellar disc, at both $z=0.3$ and $z=0.15$, correlates very strongly with the tilting rate of the disc. The angular difference at $z=0$ provides a weaker correlation with the tilting rate. Since we calculate the tilting rate between $z=0.3$ and $z=0$, the misalignment of the warp at $z=0$ should not be a good indicator of the tilting rate since $z=0.3$. Of all the correlations explored between the gas and the stellar disc, these are the strongest, suggesting that the warp may be a good indicator of the tilting rate. We improved on this correlation by building a linear multiple regression model. We found that the combination of angular misalignment of the cold gas warp at $z=0.15$ and the stellar disc to total gas ratio were the best predictors finding $R^{2}=0.897$ and $\bar{R}^{2}=0.87$.

We find a weak anti-correlation between the amount of angular momentum in the stellar disc and its tilting rate. This is not unexpected, as discs with greater angular momentum will be harder to tilt. When we take the change in the magnitude of the disc’s angular momentum and normalize by its average over the same time interval, we also find a weak correlation. Therefore, galaxies that have gained fractionally more angular momentum exhibit faster tilting.

We compared the tilting rate of the disc to the sSFR measured at $z=0$, and at $z=0.3$, finding only a weak correlation for sSFR at $z=0.3$. We also compared the average sSFR since $z=0.3$ and the peak star formation rate, still finding no correlation.

To determine the effect of the hot ($T>50000$K) gas corona, we test for correlations between the misalignment and tilting rate of the angular momentum of the hot gas corona and the stellar disc’s tilting rate. We find no correlations between the tilting rates of the two different components and the angular misalignment. We also tested this result using a higher temperature for the cutoff ($T>100000$K), finding similar results.

D15 reported that redder galaxies tended to be aligned with the minor axis of their dark halo. We measured the alignment between the stellar disc and the principal axes of the dark halo for the predominantly blue galaxies within the NIHAO sample. We find no preferential alignments, in agreement with D15.

We have presented the tilting rates of galaxies in cosmological simulations and a possible link to the warp and the role of misaligned gas accretion. As a next step, we will demonstrate, directly, the role of misaligned gas on the tilting of stellar discs.

Acknowledgements

SWFE would like to thank the Max-Planck-Institut für Astronomie, Heidelberg for their hospitality, as well as Aaron Dutton, Benjamin MacFarlane and Min Du for useful conversations. VPD is supported by STFC Consolidated grant # ST/R000786/1. TB acknowledges support from the Sonderforschungsbereich SFB 881 ‘The Milky Way System’ (subproject A2) of the German Research Foundation (DFG). Simulations were carried out on the High Performance Computing resources at New York University Abu Dhabi. This work made use of the pynbody package (Pontzen et al., 2013) to analyse simulations and the python package matplotlib (Hunter, 2007) to generate all the figures for this work. The data analysis for this work was carried out using the python packages scipy, numpy, ipython and jupyter (Jones et al., 2001; Van Der Walt et al., 2011; Pérez & Granger, 2007; Kluyver et al., 2016).