# Disc–halo interactions in ΛCDM

Disc–halo interactions in ΛCDM Abstract We present a new method for embedding a stellar disc in a cosmological dark matter halo and provide a worked example from a Λ cold dark matter zoom-in simulation. The disc is inserted into the halo at a redshift z = 3 as a zero-mass rigid body. Its mass and size are then increased adiabatically while its position, velocity, and orientation are determined from rigid-body dynamics. At z = 1, the rigid disc (RD) is replaced by an N-body disc whose particles sample a three-integral distribution function (DF). The simulation then proceeds to z = 0 with live disc (LD) and halo particles. By comparison, other methods assume one or more of the following: the centre of the RD during the growth phase is pinned to the minimum of the halo potential, the orientation of the RD is fixed, or the live N-body disc is constructed from a two rather than three-integral DF. In general, the presence of a disc makes the halo rounder, more centrally concentrated, and smoother, especially in the innermost regions. We find that methods in which the disc is pinned to the minimum of the halo potential tend to overestimate the amount of adiabatic contraction. Additionally, the effect of the disc on the subhalo distribution appears to be rather insensitive to the disc insertion method. The LD in our simulation develops a bar that is consistent with the bars seen in late-type spiral galaxies. In addition, particles from the disc are launched or ‘kicked up’ to high galactic latitudes. methods: numerical, galaxies: formation, galaxies: kinematics and dynamics, cosmology: theory 1 INTRODUCTION The structure and evolution of galaxies are determined by the spectrum of primordial density perturbations, the dynamics of stars and dark matter, and baryonic physics. Over the past two decades, there has been a concerted effort to incorporate the latter into cosmological simulations (e.g. Katz, Weinberg & Hernquist 1996; Springel & Hernquist 2003; Stinson et al. 2006; Roškar et al. 2010; Pakmor & Springel 2013; Gómez et al. 2017). While these simulations have enhanced our understanding of galaxy formation, their computational cost is high. Adding to the challenge is the complex and subgrid nature of star formation, supernova feedback, and other baryonic processes, which require ad hoc parametric models. In this work, we focus on the dynamics of disc galaxies. Our goal is to study the nature of disc–halo interactions where it is advantageous to be able to control properties of the disc such as its mass, size, and internal kinematics. Such control is not possible in ab initio simulations. Simulations of isolated disc galaxies provide an alternative arena to study galactic structure and dynamics. Moreover, many aspects of disc–halo interactions can be understood by considering the collisionless dynamics of stars and dark matter while ignoring gas physics. For example, simulations of stellar disc–bulge systems embedded in dark haloes have proved indispensable in studies of bar and spiral structure formation (see Sellwood 2013 and references therein). These simulations typically begin with systems that are in equilibrium, or nearly so. For this reason, they usually assume axisymmetric initial conditions, which are manifestly artificial. In short, discs do not come into existence as formed, highly symmetric objects but rather build up through the combined effects of gas accretion, star formation, and feedback (Vogelsberger et al. 2013; Schaye et al. 2015). Moreover, the haloes in which the real discs reside are almost certainly triaxial and clumpy (Navarro, Frenk & White 1997; Klypin et al. 1999; Moore et al. 1999). There now exists a long history of attempts to bridge the gap between simulations of isolated disc–bulge–halo systems, with their pristine initial conditions, and cosmological simulations. Kazantzidis et al. (2008), for example, followed the evolution of a Milky Way like disc in its encounter with a series of satellites whose properties were motivated by cosmological simulations. They found that the satellites ‘heated’ the disc and prompted the formation of a bar and spiral structure. Along similar lines, Purcell et al. (2011) modelled the response of the Milky Way to the gravitational effects of the Sagittarius dwarf galaxy (Sgr) by simulating disc–satellite encounters for different choices of the satellite mass. They concluded that Sgr may have triggered the development of the spiral structure seen in the Milky Way today. Continuing in this vein, Laporte et al. (2016) studied the influence of the Large Magellanic Cloud and Sgr on the Milky Way disc and found that they can create similar warps to what has been observed. The effect of a time-dependent triaxial halo was investigated in Hu & Sijacki (2016) where they found it can trigger grand-design spiral arms. Of course, the disc of the Milky Way lives within a population of satellite galaxies and, quite possibly, pure dark matter subhaloes (Klypin et al. 1999; Moore et al. 1999). With this in mind, Font et al. (2001) simulated the evolution of an isolated disc–bulge–halo model where the halo was populated by several hundred subhaloes. They concluded that substructure played only a minor role in heating the disc, a result that would seem at odds with those of Kazantzidis et al. (2008). Numerical simulations by Gauthier, Dubinski & Widrow (2006) and Dubinski et al. (2008) shed some light on this discrepancy. In those simulations, 10 per cent of the halo mass in an isolated disc–bulge–halo system was replaced by subhaloes with a mass distribution motivated by the cosmological studies of Gao et al. (2004). Gauthier et al. (2006) found that a modest amount of disc heating occurred during the first 5 Gyr, at which point satellite interactions prompted the formation of a bar, which in turn heated the disc more significantly. Not surprisingly, the timing of bar formation varied from 1 to 10 Gyr when the experiment was repeated with different initial conditions for the satellites. The aforementioned simulations have several drawbacks. First, most of them do not allow for halo triaxiality. Secondly, the disc is initialized at its present-day mass whereas real discs form over several Gyr. Finally, the subhaloes are inserted into the halo in an ad hoc fashion. Several attempts have been made to grow a stellar disc in a cosmological halo in an effort to address these shortcomings (Berentzen & Shlosman 2006; DeBuhr, Ma & White 2012; Yurin & Springel 2015). The general scheme proceeds in three stages. During the first stage, a cosmological simulation is run with pure dark matter and a suitable halo is selected. In the second, a rigid disc (RD) potential is grown slowly in the desired halo, thus allowing the halo particles to respond adiabatically to the disc’s time-varying potential. In the third stage, the RD is replaced by a live one and the simulation proceeds with live disc (LD) and halo particles. DeBuhr et al. (2012) used such a scheme to introduce stellar discs into dark matter haloes from the Aquarius Project (Springel et al. 2008). They added a RD at a redshift z = 1.3 with a mass parameter for the disc that grew linearly with the scale factor from an initial value of zero to its final value at z = 1. The disc was initially centred on the potential minimum of the halo and oriented so that its symmetry axis pointed along either the minor or major axis of the halo. During the RD phase, the motion of the disc centre of mass was determined from Newton’s 3rd law. To initialize the LD, DeBuhr et al. (2012) approximate the halo potential as a flattened, axisymmetric logarithmic potential and then determine the disc distribution function (DF) by solving the Jeans equations. Yurin & Springel (2015) introduced a number of improvements to this scheme. Most notably, they use galic to initialize the LD (Yurin & Springel 2014). This code is based on an iterative scheme for finding stationary solutions to the collisionless Boltzmann equation. The general idea for iterative codes is to begin with a set of particles that has the desired spatial distribution and some initial guess for the velocity distribution. The velocities are then adjusted so as to achieve stationarity, as measured by evolving the system and computing a certain merit function. In Yurin & Springel (2015), the initial disc was assumed to be axisymmetric with a DF that depended on two integrals of motion, the energy, E, and angular momentum, Lz. One striking, if not puzzling, result from this work is the propensity of the discs to form very strong bars. These bars are especially common in models without bulges even in cases where the disc is submaximal. In this paper, we introduce an improved scheme for inserting a LD in a cosmological halo. In particular, the centre of mass and orientation of the RD are determined by solving the standard equations of rigid body dynamics. Thus, our RD can undergo precession and nutation. The angular and linear velocities of the RD at the end of the growth phase are incorporated into the LD initial conditions. As in Yurin & Springel (2015) we use an axisymmetric approximation for the halo potential when constructing the disc DF. However, our DF is constructed from an analytic function of E, Lz, and the vertical energy Ez, which is an approximate integral of motion used in galactics (Dubinski & Kuijken 1995; Widrow, Pym & Dubinski 2008). By design, the disc DF yields a model whose density has the exponential-sech2 form. And with a three-integral DF, we have sufficient flexibility to model realistic Milky Way like discs. As discussed below, the initial disc DF may be crucial in understanding the formation of the bar. As a demonstration of our method we grow a Milky Way like disc in an approximately 1012 h−1 M⊙ halo from a cosmological zoom-in simulation. We discuss both disc dynamics and the effect our disc has on the population of subhaloes. Discs have been invoked as a means of depleting halo substructure and thus alleviating the Missing Satellite Problem, which refers to the underabundance of observed Milky Way satellites relative to the number of cold dark matter subhaloes seen in simulations (Moore et al. 1999; Klypin et al. 1999). An earlier study by D’Onghia et al. (2010) found that when a disc potential is grown in a Milky Way size cosmological halo, the abundance of substructure in the mass range 107–109 M⊙ was reduced by a factor of 2–3. Similar results were found by Sawala et al. (2017) and Garrison-Kimmel et al. (2017). The organization of the paper is as follows. In Section 2, we outline our method for inserting a LD into a cosmological simulation. We also present results from a test-bed simulation where a disc is inserted into an isolated flattened halo. We next apply our method to a cosmological zoom-in simulation. In Section 3, we focus on disc dynamics and find that the disc develops a bar, spiral structure and a warp. In addition, disc–halo interactions appear to ‘kick’ stars out of the disc and into regions normally associated with the stellar halo. In Section 4, we present our results for the spherically averaged density profile and shape of the dark matter halo as well as the distribution of subhaloes. Particular attention is paid to the sensitivity of these results to the disc insertion scheme. We conclude in Section 5 with a summary and discuss possible applications of this work. 2 INSERTING A STELLAR DISC INTO A COSMOLOGICAL HALO In this section, we detail our method for inserting a live stellar disc into a cosmological simulation. We begin with an overview of our approach and the five main simulations presented in this paper. We then describe some of the more technical aspects of the method. 2.1 Overview of simulation set Our simulations are performed with the N-body component of gadget-3, which is an updated version of gadget-2 (Springel 2005). For the cosmological simulations, 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. For the cosmological parameters, we use the results from Planck 2013 (Planck Collaboration XVI 2014) with h = 0.679, Ωb = 0.0481, Ω0 = 0.306, $$\Omega _\Lambda = 0.694$$, σ8 = 0.827, and ns = 0.962. We begin by simulating a 50 h−1 Mpc box with Nr = 5123 particles, where Nr is the effective resolution, each with a mass of ∼7.9 × 107 h−1 M⊙. We identify a Milky Way like halo in the present-day snapshot, that is, a ∼1012 M⊙ halo which has experienced no major mergers since z = 1 and which has no haloes with 2 h−1 Mpc more than half the mass of the Milky Way analogue. We then run an intermediate zoom-in simulation targeting all particles within 10 virial radii of the low-resolution halo. The initial conditions for this simulation are generated with music (Hahn & Abel 2013), which creates five nested regions from a coarse resolution of Nr = 1283 in the outskirts to an effective Nr = 20483 resolution in the targeted region. After this simulation reaches z = 0, we select all particles within 7.5 virial radii and regenerate initial conditions with one more level of refinement, giving six nested zoom regions, where now, the highest effective resolution is Nr = 40963. Our final halo is composed of approximately 107 particles, each with a mass of 1.54 × 105 h−1 M⊙. The softening lengths were selected using the criteria in Power et al. (2003) with a softening of 719 comoving pc for the highest resolution particles in the final zoom-in simulation. We found that this repeated zoom-in technique results in remarkably little contamination from coarse resolution particles within the targeted halo, giving a clean region of size 1.9 h−1 Mpc at z = 0. The dark matter only (DMO) simulation not only serves as the basis for four simulations with discs but also provides a control ‘experiment’ for our study of the effect discs have on halo properties. In each of our disc simulations, the potential of a RD is introduced at the ‘growing disc’ redshift zg. The mass parameter of the disc is then increased linearly with the scale factor a = 1/(1 + z) from zero to its final value Md at the ‘LD’ redshift zl. As described in Section 2.5, the radial and vertical scale lengths of the disc are also increased between zg and zl to account for the fact that discs grow in scale as well as mass while they are being assembled. In the first of our disc–halo simulations, dubbed MN, we introduce a rigid Miyamoto–Nagai disc (Miyamoto & Nagai 1975), whose gravitational potential is given by   $$\Phi \left(R, \,z\right) = -\frac{G M_{\rm d}}{\left(R^2 + \left(R_{\rm d} + \left(z^2 + z_{\rm d}^2 \right)^{1/2}\right)^{2}\right)^{1/2}}.$$ (1)For this simulation, which was meant to mimic the scheme in D’Onghia et al. (2010), we assume that the centre of the disc tracks the minimum of the halo potential while the orientation of the disc is fixed to be aligned with the z-axis of the simulation box. Note that this is effectively a random direction for the halo. In the remaining three disc simulations, we grow an exponential-sech2 RD potential in our halo with mass and scale length parameters that grow in time. For our fixed-orientation (FO) simulation, the position and velocity of the disc’s centre of mass are determined from Newtonian dynamics while the spin axis of the disc is initially aligned with the minor axis of the halo at z = zg and kept fixed in the simulation box frame thereafter. In this respect, the simulation is similar to the ones presented in DeBuhr et al. (2012) and Yurin & Springel (2015). For the RD simulation the disc’s orientation, which is now a function of time, is determined from Euler’s rigid body equations. In the MN, FO, and RD runs, we continue the simulation to the present epoch (z = 0) with the assumed RD potential where the mass and scale length parameters are held fixed and the position and orientation are calculated as they were during the growth phase. For our final LD simulation, we swap a LD for the RD disc at zl. Thus, the RD and LD simulations are identical prior to zl. All of our discs have a final mass of Md = 7.2 × 1010 M⊙, a final scale radius of 3.7 kpc, and a final scale height of 0.44 kpc. Our simulation parameters can be found summarized in Table 1. Table 1. A summary of the simulation parameters, as discussed in the text. Md is the final disc mass, Rd, 0 is the final disc scale radius, Nd is the number of particles used to simulate the disc, zg and zl are the redshifts when the disc beings to grow and when it becomes live (respectively), Nr is the effective resolution in the zoom-in region, and Lbox is the comoving size of the box.   DMO  MN  FO  RD  LD  Md (M⊙)  –  7.2 × 1010  7.2 × 1010  7.2 × 1010  7.2 × 1010  Rd, 0 (kpc)  –  3.7  3.7  3.7  3.7  Nd  –  –  106  106  106  zg  –  3.0  3.0  3.0  –  zl  –  1.0  1.0  1.0  1.0  Nr  4096  4096  4096  4096  4096  Lbox(Mpc h−1)  50  50  50  50  50    DMO  MN  FO  RD  LD  Md (M⊙)  –  7.2 × 1010  7.2 × 1010  7.2 × 1010  7.2 × 1010  Rd, 0 (kpc)  –  3.7  3.7  3.7  3.7  Nd  –  –  106  106  106  zg  –  3.0  3.0  3.0  –  zl  –  1.0  1.0  1.0  1.0  Nr  4096  4096  4096  4096  4096  Lbox(Mpc h−1)  50  50  50  50  50  View Large 2.2 Summary of live disc insertion scheme The first step in our disc insertion scheme is to calculate an axisymmetric approximation to the gravitational potential of the DMO halo at z = 0. We do so using an expansion in Legendre polynomials as described below. We then generate a particle representation of a stellar disc that is in equilibrium with this potential using the galactics code (Kuijken & Gilmore 1989; Widrow et al. 2008). Though galactics allows one to generate the phase space coordinates of the disc stars, at this stage, we only need the positions of the stars, which we use to represent the ‘RD’. At the zg snapshot, we incorporate the disc particles into the mass distribution of the system with the disc centred on the potential minimum of the halo. We then rerun the simulation from zg to zl with the following provisos. First, the mass of the disc is increased linearly with a from zero to its final value. Secondly, the size of the disc increases with time, which we account for by having the positions of the disc particles, as measured in the disc frame, expand with time to their final values at zl. Finally, the centre of mass and orientation of the disc are determined by integrating the equations of rigid-body dynamics. At redshift zl, the DF of the disc is recalculated assuming the same structural parameters as before but with an axisymmetric approximation to the new halo potential. An N-body disc is generated from this DF and the simulation proceeds with LD particles. In this paper, we choose zg = 3 and zl = 1 so that the growth period lasts from 2.2 to 5.9 Gyr after the big bang. This period in time roughly corresponds to the epoch of peak star formation in Milky Way like galaxies (e.g. van Dokkum et al. 2013). Our simulations during this epoch can be used to study the effect of a disc potential on the evolution of substructure. On the other hand, our simulations of the LD epoch (zl > z > 0) can also be used to study disc dynamics in a cosmologically motivated dark halo. 2.3 Halo potential Our method requires an axisymmetric approximation to the halo potential centred on the disc. This approximation is found using an expansion in spherical harmonics (see Binney & Tremaine 2008) where only the m = 0 terms are included. The potential is then expressed as an expansion in Legendre polynomials. We divide the region that surrounds the disc into spherical shells and calculate the quantities   $$m_{{\rm l},i} = \sum _{n\in S_i} m_n P_{\rm l}(\cos {\theta _n}),$$ (2)where the sum is over the halo particles of mass mn in the ith shell (Si), Pl are the Legendre polynomials, and (r,  θ,  ϕ) are spherical polar coordinates centred on the disc. The axisymmetric approximation to the potential is then   $$\Phi _{\rm h}\left(r, \, \theta \right) = \sum _{l=0}^\infty A_{\rm l}\left(r\right) P_{\rm l}\left(\cos {\theta }\right),$$ (3)where   \begin{eqnarray} A_{\rm l}(r) \!=\! -G\left( \frac{1}{r^l}\int _0^r {\rm d}r^{\prime } r^{\prime l+2} m_{\rm l}(r^{\prime }) \!+\!r^{l+1}\int _0^r {\rm d}r^{\prime } r^{\prime 1-l} m_{\rm l}(r^{\prime })\right), \end{eqnarray} (4)and ml is given by equation (2) for sufficiently small radial bins. 2.4 Disc DFs Armed with an axisymmetric approximation to the halo potential, we construct a self-consistent DF for the disc following the method outlined in Dubinski & Kuijken (1995). This DF is an analytic function of E, Lz, and Ez. By construction, the DF yields a density law for the disc which is, to a good approximation, given by   $$\rho _{\rm d}\left(R, \,z\right)\simeq \frac{M_{\rm d}}{4\pi R_{\rm d}^2 z_{\rm d}} {\rm e}^{-R/R_{\rm d}} {\rm sech}^2\left(z/z_{\rm d}\right) T\left(R_t, \Delta _t\right),$$ (5)where Md, Rd, and zd are the mass, radial scale length, and vertical scale height of the disc and $$R=\sqrt{r^2-z^2}$$. The truncation function T insures that the density falls rapidly to zero at a radius Rt + Δt. The square of the radial velocity dispersion is chosen to be proportional to the surface density, that is, σR = σR0exp ( − R/2Rd). The central velocity dispersion σR0 controls, among other things, the Toomre Q parameter and thus the susceptibility of the disc to instabilities in the disc plane. The azimuthal velocity dispersion is calculated from the radial velocity dispersion through the epicycle approximation (for details see Binney & Tremaine 2008) while the vertical velocity dispersion is adjusted to yield a constant scale height zd. We stress that although the density law is written as a function of R and z, which are not integrals of motion, the underlying DF is a function of E, Lz, and Ez, which are integrals of motion to the extent that the epicycle approximation is valid and that the potential can be approximated by an axisymmetric function. 2.5 Rigid disc dynamics During the disc growth phase, the disc mass is given by   $$M(a) = M_{\rm d} \left( \frac{a - a_{\rm g}}{a_{\rm l} - a_{\rm g}} \right),$$ (6)where ag is the scale factor evaluated at zg. The positions of the disc particles expand self-similarly in disc or body coordinates. That is, the comoving position of a disc particle in body coordinates is given by $${\boldsymbol s}_i(a) = b(a){\boldsymbol s}_i (a_{\rm l})$$ where $${\boldsymbol s}_i(a)$$ is the comoving position of the ith disc particle in the body frame   $$b(a) = b_{\rm g} + \left(1 - b_{\rm g}\right) \left( \frac{a - a_{\rm g}}{a_{\rm l} - a_{\rm g}} \right),$$ (7)where bg = b(ag), and we choose bg = 0.1. The angular velocity of the disc is described by the vector $$\boldsymbol {\omega } = \left(\omega _x, \, \omega _y, \, \omega _z\right)$$ where ωz corresponds to the spin of the disc about its symmetry axis. We assume that   $$\omega _z(a) = \omega _z(a_{\rm l})\left(\frac{M(a)}{M_{\rm d}}\right)^{1/\alpha }\,b(a),$$ (8)which insures that the disc tracks the Tully–Fisher relation, $$M_{\rm d}\propto V_{\rm d}^\alpha \propto \left(\omega _3R_{\rm d}\right)^\alpha$$ (Torres-Flores et al. 2011). In this work we set α = 3.5. The orientation of the disc is described by its Euler angles. We follow the convention of Thornton & Marion (2008) and use ϕ,  θ,  and ψ where the matrix   $${\cal R} = {\cal R}_z(\phi ) {\cal R}_y(\theta ) {\cal R}_z(\psi )$$ (9)describes the transformation from the disc body frame to the simulation frame. Here, $${\cal R}_i(\alpha )$$ is the matrix for a rotation by angle α about the ith axis. Physically, changes in ϕ and θ correspond to precession and nutation, respectively, while ψ is a degenerate rotation about the disc’s symmetry axis. The equations of motion for the Euler angles and angular velocity of the disc, which must account for the time-dependence of the disc’s moment of inertia as well as the fact that gadget-3 uses comoving coordinates, are derived in Appendix A. These equations allow us to solve for the orientation of the disc under the influence of torque due to dark matter. At the end of the growth phase, we initialize the disc with a DF that is recalculated using galactics. As we will see, during the growth phase the motion of the disc involves a mix of precession and nutation. In general, a LD is not able to support the sort of rapid nutation seen in the RD, essentially because different parts of the disc respond to torques from the halo and the disc itself differently. We therefore initialize the LD with an orientation and precessional motion given by a fixed-window moving average of the RD coordinates. 2.6 Test-bed simulation of an isolated galaxy We test our method by growing a stellar disc in an isolated, flattened halo in a non-cosmological simulation. To initialize the flattened halo we first generate a particle realization of a truncated, spherically symmetric NFW halo (Navarro et al. 1997) whose density profile is given by   $$\rho (r) = \frac{v_{\rm h}^2 a_{\rm h}}{4\pi G} \frac{1}{r\left(r + a_{\rm h}\right)^2},$$ (10)with ah = 8 kpc and vh = 400 km s−1. The halo is truncated at a radius much larger than the radius of the disc. The z and vz coordinates are then reduced by a factor of 2 and the system is evolved until it reaches approximate equilibrium. The result is an oblate halo with an axial ratio of ∼0.8 and a symmetry axis that coincides with the z-axis of the simulation box. We next grow a RD over a period of 1 Gyr to a final mass of 4.9 × 1010 M⊙ and final radial and vertical scale lengths of Rd = 2.5 kpc and zd = 200 pc, respectively. The disc is grown at an incline of 30° relative to the symmetry plane of the halo. Doing so allows us to check the rigid body integration scheme for a case when the symmetry axes of the disc and halo are initially misaligned. At t = 1 Gyr we replace the RD with a live one and evolve the system for an additional 1 Gyr. In a separate simulation, we also follow the evolution of the RD over the same time period, allowing us to compare the evolution of the LD and RD. Fig. 1 shows the Euler angles and angular velocity components for the RD and LD as a function of time. The time-dependence of ωx and ωy is characterized by an interference pattern between short ∼125 Myr period nutations and a decaying long-period precessional motion. Note that θ and ϕ for the LD track the corresponding values for the RD for t > 1 Gyr. By initializing the LD with the angular velocity of the RD, we capture the (small) precessional motion of ∼10° Gyr−1 between t = 1 and 2 Gyr. The disc settles into a preferred plane within the first 300 Myr that is intermediate between its initial symmetry plane and the initial symmetry plane of the halo. More precisely, the vector of the new minor axis is $${\boldsymbol c} = [-0.159,0.146,0.976]$$ measured at 20kpc, 12.5° from the original flattening axis. Figure 1. View largeDownload slide Kinematic variables for the RD and LD in an isolated, flattened halo as a function of time. The upper two panels show the Euler angles θ and ϕ for the RD (dashed curves) and LD (solid curves) where the LD is introduced at t = 1 Gyr (red vertical line). The bottom two panels show the x and y components of the angular velocity, as measured in the body coordinate system. In these two panels the solid curves show the δt ∼ 150 My moving average, which is used to initialize the LD. Figure 1. View largeDownload slide Kinematic variables for the RD and LD in an isolated, flattened halo as a function of time. The upper two panels show the Euler angles θ and ϕ for the RD (dashed curves) and LD (solid curves) where the LD is introduced at t = 1 Gyr (red vertical line). The bottom two panels show the x and y components of the angular velocity, as measured in the body coordinate system. In these two panels the solid curves show the δt ∼ 150 My moving average, which is used to initialize the LD. Fig. 2 shows surface density maps for the disc at two snapshots. The disc develops a weak warp due to its interaction with the halo. The development of the warp is also evident in the surface density, vertical velocity dispersion, and scale height profiles shown in Fig. 3. We see that the surface density within ∼15 kpc or 6Rd is essentially unchanged while at larger radii, there are 10–20 per cent time-dependent fluctuations. The scale height 〈z2〉1/2 increases with time and radius. At early times, the increase is most prominent beyond ∼15 kpc while at late times, the scale height increases more smoothly from the centre to the edge of the disc. Figure 2. View largeDownload slide Face-on projections of the particle distribution for two snapshots of a LD in a flattened halo. The solid line for scale is 25 kpc with a centre coincident with the disc. Figure 2. View largeDownload slide Face-on projections of the particle distribution for two snapshots of a LD in a flattened halo. The solid line for scale is 25 kpc with a centre coincident with the disc. Figure 3. View largeDownload slide Surface density, vertical velocity dispersion, and scale height profiles as a function of Galactocentric radius R for 10 snapshots equally spaced in time. The top panel shows the surface density Σ(R) divided by $$\Sigma _0(R) = \left(M_{\rm d}/2\pi R_{\rm d}^2\right)\exp {\left[-\left(R/R_{\rm d}\right) \right]}$$ in order to highlight departures from a pure exponential disc. Likewise, in the middle panel, we show the ratio σz(R)/σz, 0(R) where σz, 0 = exp ( − R/2Rd). The bottom panel shows the RMS z as a function of R. Figure 3. View largeDownload slide Surface density, vertical velocity dispersion, and scale height profiles as a function of Galactocentric radius R for 10 snapshots equally spaced in time. The top panel shows the surface density Σ(R) divided by $$\Sigma _0(R) = \left(M_{\rm d}/2\pi R_{\rm d}^2\right)\exp {\left[-\left(R/R_{\rm d}\right) \right]}$$ in order to highlight departures from a pure exponential disc. Likewise, in the middle panel, we show the ratio σz(R)/σz, 0(R) where σz, 0 = exp ( − R/2Rd). The bottom panel shows the RMS z as a function of R. 3 COSMOLOGICAL SIMULATIONS We now use our method to insert a LD with prescribed structural properties into a cosmological halo. In this section, we focus on disc dynamics while in the next, we consider the effect the disc has on the dark halo. In Fig.  4, we show the kinematic variables for the RD and LD in the RD and LD simulations. The two simulations are identical prior to z = 1 (a = 0.5) when the LD is swapped in for the rigid one. The short-period (300 Myr) oscillations in ω1 and ω2 are nutations. To initialize the LD, we use the fixed-window moving average of ωx and ωy. By and large, the Euler angles of the RD’s and LD’s track one another for z < 1, indicating that the RD is a reasonable model for a live one, at least in terms of the disc’s orientation. Figure 4. View largeDownload slide Kinematic variables for the RD and LD in our cosmological halo as a function of scale factor a. Line types are the same as in Fig. 1. The LD is introduced at a redshift z = 1 when the scale factor is a = 0.5 (red vertical line). The blue line shows the δa ∼ 0.04 moving average calculated by averaging the last 300 points in the disc integration routine. Figure 4. View largeDownload slide Kinematic variables for the RD and LD in our cosmological halo as a function of scale factor a. Line types are the same as in Fig. 1. The LD is introduced at a redshift z = 1 when the scale factor is a = 0.5 (red vertical line). The blue line shows the δa ∼ 0.04 moving average calculated by averaging the last 300 points in the disc integration routine. In Fig. 5, we show the circular speed curves at z = 1 and z = 0 for our four simulations. We see that the disc in our model is submaximal. To be precise, we have Vd/Vc ≃ 0.68 at R = 2.2Rd where Vd is the circular speed due to the disc and Vc is the total circular speed. In short, the contributions from the disc and halo to the centrifugal force are approximately equal at a radius where the disc contribution reaches its peak value. By comparison, a maximal disc is generally defined to have Vd/Vc > 0.85 (Sackett 1997). If we use Vd/Vc at 2.2Rd as the defining characteristic of the model, then our simulations match up with the F-5 simulation of Yurin & Springel (2015), although our discs are slightly warmer, with a Toomre Q parameter of 1.4 as compared with Q ≃ 0.9 for their discs and our discs are thinner (200 pc ]versus 600 pc). We note that in a two-integral disc DF, Q and the disc thickness are linked whereas in a three-integral DF, they can be set independently. The method of Yurin & Springel (2014) can be extended to consider a three-integral DF, but these models were not considered in Yurin & Springel (2015). Moreover, their two-integral model, which imposes σR = σz, violates the epicycle approximation, leading to transient system behaviour at early times when disc bars first form. Figure 5. View largeDownload slide Circular speed curve decompositions at z = 1 (top row) and z = 0 (bottom row) for (from left to right) our MN, FO, RD, and LD simulations. Halo contributions are represented as dot–dashed lines, disc contributions are represented by dashed lines, and the total rotation curve is given by a solid line. For reference, we have included the circular speed curve for the DMO halo (dot–dashed curve). Figure 5. View largeDownload slide Circular speed curve decompositions at z = 1 (top row) and z = 0 (bottom row) for (from left to right) our MN, FO, RD, and LD simulations. Halo contributions are represented as dot–dashed lines, disc contributions are represented by dashed lines, and the total rotation curve is given by a solid line. For reference, we have included the circular speed curve for the DMO halo (dot–dashed curve). The circular speed curves in Fig. 5 show little change exterior to ∼2Rd after zl, thus providing another indication that the LD was close equilibrium when it was swapped in for the rigid one. The formation of a bar is evident in the circular speed, and we can infer the bar contributes substantially inside 2.2Rd ≃ 8 kpc. The halo contribution at R = 2.2Rd is about 20 per cent larger in the four disc runs than in the DMO one due to adiabatic contraction. Interestingly, the halo in the MN run shows somewhat more contraction than in the RD and LD runs. We note that in the MN run, the disc potential tracks the potential minimum of the halo whereas in the RD/LD case, the disc’s position is determined from Newtonian dynamics. In general, the centre of the disc tracks the halo potential minimum so long as the potential changes slowly with time. However, during a major merger (and indeed, just such an event occurs at z = 2) there are rapid changes in the halo potential and the position of the disc, as determined by Newtonian dynamics, can differ significant from the minimum of the halo potential. Evidently, the ad hoc prescription of growing a disc at the halo’s potential minimum may, in some cases, overestimate the effect of adiabatic contraction. 3.1 Bar formation In Fig. 6, we show orthogonal projections of the disc density in our LD simulation at four epochs between z = 1 (lookback time of 7.9 Gyr) and the present epoch. During the first billion years of LD evolution, the disc develops a bar and spiral structure. In addition, there is a warp in the outer disc extending several kiloparsecs above the mid-plane of the inner disc. By the present epoch, the bar has grown in length and intensified and the edge-on view shows the classic X-pattern. Figure 6. View largeDownload slide Projected density along three orthogonal directions for the LD at four epochs between z = 1 and z = 0. The projections are presented in physical units. The solid line for scale is 37 kpc with a centre coincident with the disc. Figure 6. View largeDownload slide Projected density along three orthogonal directions for the LD at four epochs between z = 1 and z = 0. The projections are presented in physical units. The solid line for scale is 37 kpc with a centre coincident with the disc. We consider the usual parameter bar strength A2 = |c2| where   $$c_m = \frac{1}{M_S}\sum _{j\in S} m_j {\rm e}^{im\phi _j}.$$ (11)Here, S is some circularly symmetric region of the disc (e.g. a circular annulus) and the sum is over all particles labelled by j and with mass mj that are inside S. We find that A2 for the inner 2Rd of the disc reaches 0.43 at t = 6.7 Gyr, decreases to 0.36 by t = 9.2 Gyr, presumably because the bar has buckled, and then increases to 0.47 by the present epoch. On the other hand, A2 for the entire disc increases to 0.27, decreases to 0.23, and then increases to 0.28 for the same epochs. Note that the inner 2Rd of the disc contains 60 per cent of the mass. Thus, the fact that $$A_{2,2R_{\rm d}}/A_{2,{\rm tot}}\simeq 0.6$$ implies that most of the bar mass resides within the inner 2Rd. The bars in Yurin & Springel (2015) seem to be stronger then then ones in our simulations – they find A2 ≃ 0.6 but use a non-standard formula for A2. Moreover, their bars appear to extend across most of the disc. In terms of disc dynamics, the main difference between our simulations is the fact that we use a three-integral DF for the disc whereas they use a two-integral DF. In the latter, the velocity dispersion in the radial and vertical directions are the same. Thus, the radial dispersion, which fixes the Toomre Q parameter, also determines the thickness of the disc. We note that their initial discs are a factor of 2 or 3 thicker than ours. We speculate that the bars that develop in these thick discs are less susceptible to buckling and therefore able to grow stronger and longer. These ideas will be investigated in more detail in a future publication. 3.2 Kicked-up stars and disc heating The outer part of the disc suffers considerable disruption and warping presumably through its interaction with the main halo or substructure. The right-most panel of the a = 0.55 snapshot in Fig. 6, for example, shows a classic integral-sign warp. The other snapshots show that a significant number of disc particles have orbits that now take them to high galactic latitudes. The impressions one has from the density projections are borne out in Fig. 7 where we show the surface density and scale height profiles at different times. Bar formation redistributes mass in the disc leaving a deficit (relative to the initial exponential disc) between 5 and 15 kpc. The disc becomes thicker and its vertical velocity dispersion increases though a combination of disc–halo interactions and the effects of the bar and spiral structure (Gauthier et al. 2006; Dubinski et al. 2008; Kazantzidis et al. 2008). Figure 7. View largeDownload slide Surface density, vertical velocity dispersion, and scale height profiles of the LD for 10 snapshots equally spaced in scale factor a between a = 0.5 (z = 1) and a = 1 (present epoch). Panels are the same as in Fig. 3. Figure 7. View largeDownload slide Surface density, vertical velocity dispersion, and scale height profiles of the LD for 10 snapshots equally spaced in scale factor a between a = 0.5 (z = 1) and a = 1 (present epoch). Panels are the same as in Fig. 3. A striking feature of the simulations are the streams of disc stars well above the disc plane. These stars may represent an example of a kicked-up disc, which has been seen in other N-body simulations (Purcell, Bullock & Kazantzidis 2010; McCarthy et al. 2012) and invoked to explain kinematic and spectroscopic observations of M31 (Dorman et al. 2013) and the Monoceros Ring (e.g. Newberg et al. 2002; Ibata et al. 2003). The idea is that interactions between the disc and both satellite galaxies and halo substructure liberate stars from the disc, launching them to regions of the galaxy normally associated with the stellar halo. Our LD simulation corroborates this hypothesis and is in broad agreement with previous numerical work. Finally, we note that the a = 1 panel of Fig. 6 shows a relatively thin, stream-like structure, above the disc which is qualitatively similar to the Anti-centre Stream (ACS; Grillmair 2006). While the ACS is believed to be due to the disruption of a globular cluster (e.g. Grillmair 2006), Fig. 6 suggests that perturbations to the disc can create similar features. Intriguingly, de Boer, Belokurov & Koposov (2018) recently found that the ACS is rotating in the same sense as the Milky Way disc. 4 HALO SUBSTRUCTURE IN THE PRESENCE OF A DISC In this section, we consider the effect of a disc on a halo’s structural properties such as its spherically averaged density profile, its shape, and its subhalo population. An examination of the DMO simulation shows that the halo we have selected builds up through a series of mergers and accretion events, but that by z = 1 it has settled into a relatively relaxed state with an NFW profile that evolves very little between z = 1 and z = 0 within the inner 100 kpc. Our sequence of simulations, (MN, FO, RD, and LD) allow us to tease out the effects of different disc insertion methods. The MN simulation, for example, pins the centre of the disc to the minimum of the halo potential, whereas the other simulations dynamically evolve the position and velocity of the disc potential via Newtonian mechanics. The MN and FO simulations both assume that the orientation of the disc potential during the growth phase is fixed whereas RD and LD solve for the orientation using rigid body dynamics. 4.1 Global properties of the halo In Fig. 8, we show the ratio of the spherically averaged density profile in the four disc runs to that from the DMO run. At z = 1 the haloes in the FO, RD, and MN runs show evidence for adiabatic contraction with the density in the inner ∼30 kpc increasing by a factor of 1.2–2.1. The effect is strongest in the MN simulation, which is to be expected since the halo in that case always sees the disc potential at the minimum of its potential. Of course, this prescription is unphysical. In general, and especially during a major merger, the disc and halo potential minimum will not necessarily coincide. Figure 8. View largeDownload slide The ratio of halo density to the DMO simulation for MN (green), FO (teal), RD (red), and LD (purple) at z = 1 (dashed) and z = 0 (solid). The presence of the disc significantly increases the central concentration of the halo. Figure 8. View largeDownload slide The ratio of halo density to the DMO simulation for MN (green), FO (teal), RD (red), and LD (purple) at z = 1 (dashed) and z = 0 (solid). The presence of the disc significantly increases the central concentration of the halo. Between z = 1 and z = 0, the mass of the disc is constant. Adiabatic contraction ceases but the halo still responds to the time-varying disc potential. Interestingly, at intermediate radii (between ∼10 and 40 kpc) the density profile of the halo settles back to a state close to that found in the DMO run. Perhaps most striking is the fact that the halo in the LD run becomes more centrally concentrated than the halo in any of the other cases. One possible explanation is that dynamical friction from the disc drags dark matter subhaloes towards the centre of the halo where they are tidally disrupted. In Fig. 9, we show the minor-to-major (cr/ar) and intermediate-to-major (br/ar) axial ratios as a function of radius for both the z = 1 and z = 0 snapshots. The axial ratios are calculated by diagonalizing the moment of inertia tensor in linearly spaced radial shells. The DMO halo is triaxial with cr/ar ≃ 0.75 and br/ar ≃ 0.85. Note that the axial ratio profiles are smoother at z = 0 than at z = 1, which supports the observation that the halo has settled into a more relaxed state over the past 7 or so billion years. In general, discs tend to make haloes more spherical, a result that has been known for some time from both collisionless and hydrodynamical simulations (e.g. Dubinski 1994; Zemp et al. 2012). Figure 9. View largeDownload slide Axial ratios as a function of radius. Shown are the minor-to-major axial ratio (top panel) and the intermediate-to-major axial ratio (bottom panel) at z = 1 (dashed curves) and z = 0 (solid curves). Blue corresponds to DMO, green to MN, teal to FO, red to RD, and purple to LD. Figure 9. View largeDownload slide Axial ratios as a function of radius. Shown are the minor-to-major axial ratio (top panel) and the intermediate-to-major axial ratio (bottom panel) at z = 1 (dashed curves) and z = 0 (solid curves). Blue corresponds to DMO, green to MN, teal to FO, red to RD, and purple to LD. Evidently, the MN halo is rounder, especially in the inner part, the FO halo. Recall that the main difference between these two cases is that the MN disc is pinned to the potential minimum of the halo. It is perhaps not surprising then that, as with adiabatic contraction, it has a stronger effect on the halo’s shape. We also note that the axial ratio profiles for the RD and LD simulations are fairly similar. 4.2 Subhalo populations We now turn our attention to halo substructure. We identify subhaloes and determine their positions and masses using rockstar (Behroozi, Wechsler & Wu 2013), which employs a friends-of-friends algorithm in six phase space dimensions. We consider only those subhaloes with mass ms between mmin = 107.5 M⊙ and mmax = 1010.5 M⊙. Subhaloes at the lower end of this range are marginally resolved with ∼100 particles, above which the subhalo mass can be trusted (e.g. Onions et al. 2012), while those at the upper end contain ∼ 3 per cent of the halo’s virial mass. In Fig. 10, we show the cumulative mass in subhaloes as a function of Galactocentric radius:   $$M_{\rm s}\left(<r\right) = \int _0^r {\rm d}r \int _{m_{\rm min}}^{m_{\rm max}} {\rm d}m_{\rm s}\,m_{\rm s}\, \frac{{\rm d}^2 N}{{\rm d}m_{\rm s}\,{\rm d}r}.$$ (12)We also show the cumulative number of subhaloes. In general, the presence of a disc depletes substructure inside about 30 kpc but leaves the outer substructures unaffected. Figure 10. View largeDownload slide Cumulative mass in subhaloes inside a radius r (upper panel) and cumulative number of subhaloes (lower panel). We consider only subhaloes within 500 kpc of the halo centre and with a mass above 107.5 M⊙. The curves are blue (DMO), green (MN), teal (FO), red (RD), and purple (LD). Figure 10. View largeDownload slide Cumulative mass in subhaloes inside a radius r (upper panel) and cumulative number of subhaloes (lower panel). We consider only subhaloes within 500 kpc of the halo centre and with a mass above 107.5 M⊙. The curves are blue (DMO), green (MN), teal (FO), red (RD), and purple (LD). We next consider the differential mass distribution as a function of subhalo mass. Recall that for a pure dark matter halo, $${\rm d}N/{\rm d}\ln (m_{\rm s})\propto m_{\rm s}^{-p}$$ where p ≃ 0.9 (e.g. Gao et al. 2004). In Fig. 11, we therefore show the quantity m0.9 dN/dln (ms) in order to enhance the differences between the different disc runs. We see that the halo population between ms ≃ mmin and ms ≃ 109 M⊙ is depleted, but only by about 20–30  per cent. Taken together, Figs 10 and 11 imply that the main depletion of the subhaloes occurs within the inner regions of the parent halo, in agreement with D’Onghia et al. (2010), Sawala et al. (2017), and Garrison-Kimmel et al. (2017). That said, the depletion of subhaloes seems rather insensitive to the disc insertion method,although we caution that only a single halo was used. Our results, being mainly in agreement with previous work, should still be viewed with caution due the consideration of a single host halo. Figure 11. View largeDownload slide Differential mass distribution multiplied by M0.9 for subhaloes above 107.5 M⊙. We make an outer radius cut at 500 kpc. The curves are blue (DMO), green (MN), teal (FO), red (RD), and purple (LD). Figure 11. View largeDownload slide Differential mass distribution multiplied by M0.9 for subhaloes above 107.5 M⊙. We make an outer radius cut at 500 kpc. The curves are blue (DMO), green (MN), teal (FO), red (RD), and purple (LD). 4.3 Case study: a Sagittarius-like dwarf Observations of the Milky Way’s dwarf galaxies and associated tidal streams provide a potentially powerful probe of the Galactic potential and thus the Galaxy’s dark halo. One of the best-studied examples is the Sagittarius dwarf (Ibata, Gilmore & Irwin 1994). Fortunately, our simulation has a satellite galaxy with similar properties, namely, a dark matter mass of 1.8 × 1010 M⊙ at z = 1.0 and an orbit that takes it close to the disc. We identify this object in the five simulations using the rockstar halo catalogues. We then gather a list of IDs for all the bound particles at an early time before the dwarf is disrupted and follow these same particles in later snapshots using a binary search tree look-up scheme. In Fig. 12, we show the evolution of this subhalo between redshift zl = 1 and z = 0. The first row shows the baseline evolution in our DMO simulation. The dwarf develops leading and trailing tidal tails during the first few billion years. By the present epoch, the tidal debris has dispersed throughout the halo. Figure 12. View largeDownload slide X–Y projections for a selected Sagittarius-like dwarf galaxy. The rows from top to bottom are no disc, a fixed Miyamoto–Nagai disc, a RD, and a LD. The scale factors in columns from left to right are 0.5, 0.55, 0.6, 0.65, and 0.7. The frame edges are 295 kpc on each side. Figure 12. View largeDownload slide X–Y projections for a selected Sagittarius-like dwarf galaxy. The rows from top to bottom are no disc, a fixed Miyamoto–Nagai disc, a RD, and a LD. The scale factors in columns from left to right are 0.5, 0.55, 0.6, 0.65, and 0.7. The frame edges are 295 kpc on each side. The next four rows show the same satellite in our disc simulations. Perhaps the most noticeable result is that there are stronger features in the tidal debris at the present epoch once a disc is included. The detailed morphology of the tidal debris is certainly different from one disc simulation to the next. By eye, debris in MN and FO look somewhat similar as does the debris in RD and LD. Perhaps the most noticeable result is that the tidal debris extends to larger galactocentric radii when a disc is included. The detailed morphology of the tidal debris clearly depends on the disc insertion method. By eye, tidal debris appears more isotropic in MN and FO than in RD and LD. The implication is that fixed potentials are more efficient at disrupting massive satellites than a potential which can respond to the satellite’s presence. However, we have only a single example of massive satellite disruption, and we caution that more examples of this behaviour are needed to test this hypothesis. 5 CONCLUSIONS Simulations in which a stellar disc is inserted ‘by hand’ into a cosmological N-body halo provide a compromise between simulations of isolated disc–halo systems and cosmological simulations that include gasdynamics and star formation. Our method builds on the scheme used by Berentzen & Shlosman (2006) and DeBuhr et al. (2012), and refined by Yurin & Springel (2015). The basic idea is to introduce, at a redshift zg, a RD with zero mass into a halo within a cosmological zoom-in simulation. Between zg and zl the disc is treated as an external potential with a mass and size that increase adiabatically to their present-day values. At zl, the RD is replaced by an N-body one and the simulation proceeds to the present epoch with LD and halo particles. Our method improves upon previous ones in two important ways. First, during the growth phase (zg > z > zl) the position and orientation of the disc evolve according to the standard equations of rigid-body dynamics. Thus, the disc in our scheme receives its linear and angular momentum with the halo in a self-consistent fashion and is therefore able to move, precess, and nutate due to torques from the halo. While previous methods introduced aspects of rigid-body dynamics during the growth phase none appear to have implemented the full dynamical equations have done here (D’Onghia et al. 2010; DeBuhr et al. 2012; Yurin & Springel 2015). Our sequence MN, FO, RD, and LD of simulations highlights where the details of the disc insertion scheme are important and where they are not. For example, schemes in which the disc tracks the minimum of the halo potential tend to overestimate the effects of adiabatic contraction. On the other hand, the effect of the depletion of halo substructure seems to be rather insensitive to the details of how the disc is introduced into the simulation. Disc insertion schemes such as the one introduced in this paper, provide an attractive arena for studies of galactic dynamics. In particular, they allow one to study the interaction between a stellar disc and a realistic dark halo with computationally inexpensive simulations while maintaining some level of control over the structural parameters of the disc. We fully intend to leverage these advantages in future work. Acknowledgements LMW and JSB are supported by a Discovery Grant with the Natural Sciences and Engineering Research Council of Canada. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement no. 308024. DE acknowledges financial support from the ERC. REFERENCES Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013, ApJ , 762, 109 https://doi.org/10.1088/0004-637X/762/2/109 CrossRef Search ADS   Berentzen I., Shlosman I., 2006, ApJ , 648, 807 https://doi.org/10.1086/506016 CrossRef Search ADS   Binney J., Tremaine S., 2008, Galactic Dynamics , 2nd edn. Princeton University Press, Princeton, NJ D’Onghia E., Springel V., Hernquist L., Keres D., 2010, ApJ , 709, 1138 https://doi.org/10.1088/0004-637X/709/2/1138 CrossRef Search ADS   de Boer T. J. L., Belokurov V., Koposov S. E., 2018, MNRAS , 473, 647 CrossRef Search ADS   DeBuhr J., Ma C.-P., White S. D. M., 2012, MNRAS , 426, 983 https://doi.org/10.1111/j.1365-2966.2012.21910.x CrossRef Search ADS   Dorman C. E. et al.  , 2013, ApJ , 779, 103 https://doi.org/10.1088/0004-637X/779/2/103 CrossRef Search ADS   Dubinski J., 1994, ApJ , 431, 617 https://doi.org/10.1086/174512 CrossRef Search ADS   Dubinski J., Kuijken K., 1995, ApJ , 442, 492 https://doi.org/10.1086/175456 CrossRef Search ADS   Dubinski J., Gauthier J.-R., Widrow L., Nickerson S., 2008, in Funes J. G., Corsini E. M., eds, ASP Conf. Ser. Vol. 396, Formation and Evolution of Galaxy Disks . Astron. Soc. Pac., San Francisco, p. 321 Font A. S., Navarro J. F., Stadel J., Quinn T., 2001, ApJ , 563, L1 https://doi.org/10.1086/338479 CrossRef Search ADS   Gao L., White S. D. M., Jenkins A., Stoehr F., Springel V., 2004, MNRAS , 355, 819 https://doi.org/10.1111/j.1365-2966.2004.08360.x CrossRef Search ADS   Garrison-Kimmel S. et al.  , 2017, MNRAS , 471, 1709 CrossRef Search ADS   Gauthier J.-R., Dubinski J., Widrow L. M., 2006, ApJ , 653, 1180 https://doi.org/10.1086/508860 CrossRef Search ADS   Gómez F. A., White S. D. M., Grand R. J. J., Marinacci F., Springel V., Pakmor R., 2017, MNRAS , 465, 3446 CrossRef Search ADS   Grillmair C. J., 2006, ApJ , 651, L29 https://doi.org/10.1086/509255 CrossRef Search ADS   Hahn O., Abel T., 2013, Astrophysics Source Code Library , record ascl:1311.011 Hu S., Sijacki D., 2016, MNRAS , 461, 2789 https://doi.org/10.1093/mnras/stw1463 CrossRef Search ADS   Ibata R. A., Gilmore G., Irwin M. J., 1994, Nature , 370, 194 https://doi.org/10.1038/370194a0 CrossRef Search ADS   Ibata R. A., Irwin M. J., Lewis G. F., Ferguson A. M. N., Tanvir N., 2003, MNRAS , 340, L21 https://doi.org/10.1046/j.1365-8711.2003.06545.x CrossRef Search ADS   Katz N., Quinn T., Bertschinger E., Gelb J. M., 1994, MNRAS , 270, L71 https://doi.org/10.1093/mnras/270.1.L71 CrossRef Search ADS   Katz N., Weinberg D. H., Hernquist L., 1996, ApJS , 105, 19 https://doi.org/10.1086/192305 CrossRef Search ADS   Kazantzidis S., Bullock J. S., Zentner A. R., Kravtsov A. V., Moustakas L. A., 2008, ApJ , 688, 254 https://doi.org/10.1086/591958 CrossRef Search ADS   Klypin A., Kravtsov A. V., Valenzuela O., Prada F., 1999, ApJ , 522, 82 https://doi.org/10.1086/307643 CrossRef Search ADS   Kuijken K., Gilmore G., 1989, MNRAS , 239, 571 https://doi.org/10.1093/mnras/239.2.571 CrossRef Search ADS   Laporte C. F. P., Gómez F. A., Besla G., Johnston K. V., Garavito-Camargo N., 2016, MNRAS , preprint (arXiv:1608.04743) McCarthy I. G., Font A. S., Crain R. A., Deason A. J., Schaye J., Theuns T., 2012, MNRAS , 420, 2245 https://doi.org/10.1111/j.1365-2966.2011.20189.x CrossRef Search ADS   Miyamoto M., Nagai R., 1975, PASJ , 27, 533 Moore B., Ghigna S., Governato F., Lake G., Quinn T., Stadel J., Tozzi P., 1999, ApJ , 524, L19 https://doi.org/10.1086/312287 CrossRef Search ADS   Navarro J. F., Frenk C. S., White S. D. M., 1994, MNRAS , 267, L1 https://doi.org/10.1093/mnras/267.1.L1 CrossRef Search ADS   Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ , 490, 493 https://doi.org/10.1086/304888 CrossRef Search ADS   Newberg H. J. et al.  , 2002, ApJ , 569, 245 https://doi.org/10.1086/338983 CrossRef Search ADS   Onions J. et al.  , 2012, MNRAS , 423, 1200 https://doi.org/10.1111/j.1365-2966.2012.20947.x CrossRef Search ADS   Oñorbe J., Garrison-Kimmel S., Maller A. H., Bullock J. S., Rocha M., Hahn O., 2014, MNRAS , 437, 1894 https://doi.org/10.1093/mnras/stt2020 CrossRef Search ADS   Pakmor R., Springel V., 2013, MNRAS , 432, 176 https://doi.org/10.1093/mnras/stt428 CrossRef Search ADS   Planck Collaboration XVI, 2014, A&A , 571, A16 CrossRef Search ADS   Power C., Navarro J. F., Jenkins A., Frenk C. S., White S. D. M., Springel V., Stadel J., Quinn T., 2003, MNRAS , 338, 14 https://doi.org/10.1046/j.1365-8711.2003.05925.x CrossRef Search ADS   Purcell C. W., Bullock J. S., Kazantzidis S., 2010, MNRAS , 404, 1711 Purcell C. W., Bullock J. S., Tollerud E. J., Rocha M., Chakrabarti S., 2011, Nature , 477, 301 https://doi.org/10.1038/nature10417 CrossRef Search ADS PubMed  Quinn T., Katz N., Stadel J., Lake G., 1997, eprint (arXiv:e-prints) Roškar R., Debattista V. P., Brooks A. M., Quinn T. R., Brook C. B., Governato F., Dalcanton J. J., Wadsley J., 2010, MNRAS , 408, 783 https://doi.org/10.1111/j.1365-2966.2010.17178.x CrossRef Search ADS   Sackett P. D., 1997, ApJ , 483, 103 https://doi.org/10.1086/304223 CrossRef Search ADS   Sawala T., Pihajoki P., Johansson P. H., Frenk C. S., Navarro J. F., Oman K. A., White S. D. M., 2017, MNRAS , 467, 4383 https://doi.org/10.1093/mnras/stx360 CrossRef Search ADS   Schaye J. et al.  , 2015, MNRAS , 446, 521 https://doi.org/10.1093/mnras/stu2058 CrossRef Search ADS   Sellwood J. A., 2013, Dynamics of Disks and Warps . Springer Science+Business Media, Dordrecht, p. 923 Springel V., 2005, MNRAS , 364, 1105 https://doi.org/10.1111/j.1365-2966.2005.09655.x CrossRef Search ADS   Springel V., Hernquist L., 2003, MNRAS , 339, 289 https://doi.org/10.1046/j.1365-8711.2003.06206.x CrossRef Search ADS   Springel V. et al.  , 2008, MNRAS , 391, 1685 https://doi.org/10.1111/j.1365-2966.2008.14066.x CrossRef Search ADS   Stinson G., Seth A., Katz N., Wadsley J., Governato F., Quinn T., 2006, MNRAS , 373, 1074 https://doi.org/10.1111/j.1365-2966.2006.11097.x CrossRef Search ADS   Thornton S. T., Marion J. B., 2008, Classical Dynamics of Particles and Systems , 5th edn. Cengage Learning Torres-Flores S., Epinat B., Amram P., Plana H., Mendes de Oliveira C., 2011, MNRAS , 416, 1936 https://doi.org/10.1111/j.1365-2966.2011.19169.x CrossRef Search ADS   van Dokkum P. G. et al.  , 2013, ApJ , 771, L35 https://doi.org/10.1088/2041-8205/771/2/L35 CrossRef Search ADS   Vogelsberger M., Genel S., Sijacki D., Torrey P., Springel V., Hernquist L., 2013, MNRAS , 436, 3031 https://doi.org/10.1093/mnras/stt1789 CrossRef Search ADS   Widrow L. M., Pym B., Dubinski J., 2008, ApJ , 679, 1239 https://doi.org/10.1086/587636 CrossRef Search ADS   Yurin D., Springel V., 2014, Astrophysics Source Code Library , record ascl:1408.008 Yurin D., Springel V., 2015, MNRAS , 452, 2367 https://doi.org/10.1093/mnras/stv1454 CrossRef Search ADS   Zemp M., Gnedin O. Y., Gnedin N. Y., Kravtsov A. V., 2012, ApJ , 748, 54 https://doi.org/10.1088/0004-637X/748/1/54 CrossRef Search ADS   APPENDIX A: EULER’S EQUATIONS IN COMOVING COORDINATES The time-evolution of the angular momentum vector $${\boldsymbol L}$$ of a rigid body acted upon by a torque $$\boldsymbol {\tau }$$ is given by   $$\left(\frac{{{\rm d}} {\boldsymbol L}}{{{\rm d}} t}\right)_f = \left(\frac{{{\rm d}} {\boldsymbol L}}{{{\rm d}} t}\right)_{\rm b} + \boldsymbol \omega \times {\boldsymbol L} = \boldsymbol {\tau },$$ (A1)where the subscripts f and b denote the frame of the simulation box and the body frame, respectively. In physical coordinates, $${\boldsymbol L} = {\boldsymbol r \times p}$$. Alternatively, we can write $${\boldsymbol L} = {\boldsymbol s\times q}$$ where $${\boldsymbol s} = a^{-1}{\boldsymbol r}$$ refers to comoving coordinates and $${\boldsymbol q} = a^2\dot{\boldsymbol s}$$ is the conjugate momentum to $${\boldsymbol s}$$ (see Quinn et al. 1997). For a rigid body, the components of the angular momentum are given by Li = Iijωj where i, j run over x,  y,  z and there is an implied sum over j. Since gadget-3 uses comoving coordinates, we write Iij = a2Jij where J is the moment of inertia tensor written in terms of the comoving coordinates, $${\boldsymbol s}$$, rather than the physical coordinates, $${\boldsymbol r}$$. For convenience, we define a ‘comoving’ angular velocity $$\boldsymbol {\varpi } = a^{-2}\boldsymbol {\omega }$$. We then have Li = Jijϖj. Note that because of the symmetry of our disc, the moment of inertia tensor is diagonal with Jxx = Jyy = Jzz ≡ J/2. The equations of motion for the Euler angles and the disc angular velocity are then given by the standard Euler equations of rigid body dynamics, modified to account for the time-dependence of the disc’s moment of inertia:   $$\frac{{\rm d}\phi }{{\rm d}t} = a^{-2}\sin ^{-1}{\theta } \left(\varpi _x\sin (\psi ) + \varpi _y \cos (\psi )\right),$$ (A2)  $$\frac{{\rm d}\theta }{{\rm d}t} = a^{-2}\left(\varpi _1\cos (\psi ) - \varpi _y \sin (\psi )\right),$$ (A3)  $$J\frac{\varpi _x}{{\rm d}t} + \varpi _x\frac{{\rm d}J}{{\rm d}t} + J\varpi _y\varpi _z = \tau _x,$$ (A4)and   $$\frac{\varpi _y}{{\rm d}t} + \varpi _y\frac{{\rm d}J}{{\rm d}t} - J\varpi _x\varpi _z = \tau _y.$$ (A5)We have omitted the equations for ψ (rotations in the body frame about the symmetry axis) and ϖz since these are determined directly from equation (8). © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Monthly Notices of the Royal Astronomical Society Oxford University Press

# Disc–halo interactions in ΛCDM

, Volume 476 (1) – May 1, 2018
12 pages

/lp/ou_press/disc-halo-interactions-in-cdm-lrJTgwD0ri
Publisher
The Royal Astronomical Society
ISSN
0035-8711
eISSN
1365-2966
D.O.I.
10.1093/mnras/sty120
Publisher site
See Article on Publisher Site

### Abstract

Abstract We present a new method for embedding a stellar disc in a cosmological dark matter halo and provide a worked example from a Λ cold dark matter zoom-in simulation. The disc is inserted into the halo at a redshift z = 3 as a zero-mass rigid body. Its mass and size are then increased adiabatically while its position, velocity, and orientation are determined from rigid-body dynamics. At z = 1, the rigid disc (RD) is replaced by an N-body disc whose particles sample a three-integral distribution function (DF). The simulation then proceeds to z = 0 with live disc (LD) and halo particles. By comparison, other methods assume one or more of the following: the centre of the RD during the growth phase is pinned to the minimum of the halo potential, the orientation of the RD is fixed, or the live N-body disc is constructed from a two rather than three-integral DF. In general, the presence of a disc makes the halo rounder, more centrally concentrated, and smoother, especially in the innermost regions. We find that methods in which the disc is pinned to the minimum of the halo potential tend to overestimate the amount of adiabatic contraction. Additionally, the effect of the disc on the subhalo distribution appears to be rather insensitive to the disc insertion method. The LD in our simulation develops a bar that is consistent with the bars seen in late-type spiral galaxies. In addition, particles from the disc are launched or ‘kicked up’ to high galactic latitudes. methods: numerical, galaxies: formation, galaxies: kinematics and dynamics, cosmology: theory 1 INTRODUCTION The structure and evolution of galaxies are determined by the spectrum of primordial density perturbations, the dynamics of stars and dark matter, and baryonic physics. Over the past two decades, there has been a concerted effort to incorporate the latter into cosmological simulations (e.g. Katz, Weinberg & Hernquist 1996; Springel & Hernquist 2003; Stinson et al. 2006; Roškar et al. 2010; Pakmor & Springel 2013; Gómez et al. 2017). While these simulations have enhanced our understanding of galaxy formation, their computational cost is high. Adding to the challenge is the complex and subgrid nature of star formation, supernova feedback, and other baryonic processes, which require ad hoc parametric models. In this work, we focus on the dynamics of disc galaxies. Our goal is to study the nature of disc–halo interactions where it is advantageous to be able to control properties of the disc such as its mass, size, and internal kinematics. Such control is not possible in ab initio simulations. Simulations of isolated disc galaxies provide an alternative arena to study galactic structure and dynamics. Moreover, many aspects of disc–halo interactions can be understood by considering the collisionless dynamics of stars and dark matter while ignoring gas physics. For example, simulations of stellar disc–bulge systems embedded in dark haloes have proved indispensable in studies of bar and spiral structure formation (see Sellwood 2013 and references therein). These simulations typically begin with systems that are in equilibrium, or nearly so. For this reason, they usually assume axisymmetric initial conditions, which are manifestly artificial. In short, discs do not come into existence as formed, highly symmetric objects but rather build up through the combined effects of gas accretion, star formation, and feedback (Vogelsberger et al. 2013; Schaye et al. 2015). Moreover, the haloes in which the real discs reside are almost certainly triaxial and clumpy (Navarro, Frenk & White 1997; Klypin et al. 1999; Moore et al. 1999). There now exists a long history of attempts to bridge the gap between simulations of isolated disc–bulge–halo systems, with their pristine initial conditions, and cosmological simulations. Kazantzidis et al. (2008), for example, followed the evolution of a Milky Way like disc in its encounter with a series of satellites whose properties were motivated by cosmological simulations. They found that the satellites ‘heated’ the disc and prompted the formation of a bar and spiral structure. Along similar lines, Purcell et al. (2011) modelled the response of the Milky Way to the gravitational effects of the Sagittarius dwarf galaxy (Sgr) by simulating disc–satellite encounters for different choices of the satellite mass. They concluded that Sgr may have triggered the development of the spiral structure seen in the Milky Way today. Continuing in this vein, Laporte et al. (2016) studied the influence of the Large Magellanic Cloud and Sgr on the Milky Way disc and found that they can create similar warps to what has been observed. The effect of a time-dependent triaxial halo was investigated in Hu & Sijacki (2016) where they found it can trigger grand-design spiral arms. Of course, the disc of the Milky Way lives within a population of satellite galaxies and, quite possibly, pure dark matter subhaloes (Klypin et al. 1999; Moore et al. 1999). With this in mind, Font et al. (2001) simulated the evolution of an isolated disc–bulge–halo model where the halo was populated by several hundred subhaloes. They concluded that substructure played only a minor role in heating the disc, a result that would seem at odds with those of Kazantzidis et al. (2008). Numerical simulations by Gauthier, Dubinski & Widrow (2006) and Dubinski et al. (2008) shed some light on this discrepancy. In those simulations, 10 per cent of the halo mass in an isolated disc–bulge–halo system was replaced by subhaloes with a mass distribution motivated by the cosmological studies of Gao et al. (2004). Gauthier et al. (2006) found that a modest amount of disc heating occurred during the first 5 Gyr, at which point satellite interactions prompted the formation of a bar, which in turn heated the disc more significantly. Not surprisingly, the timing of bar formation varied from 1 to 10 Gyr when the experiment was repeated with different initial conditions for the satellites. The aforementioned simulations have several drawbacks. First, most of them do not allow for halo triaxiality. Secondly, the disc is initialized at its present-day mass whereas real discs form over several Gyr. Finally, the subhaloes are inserted into the halo in an ad hoc fashion. Several attempts have been made to grow a stellar disc in a cosmological halo in an effort to address these shortcomings (Berentzen & Shlosman 2006; DeBuhr, Ma & White 2012; Yurin & Springel 2015). The general scheme proceeds in three stages. During the first stage, a cosmological simulation is run with pure dark matter and a suitable halo is selected. In the second, a rigid disc (RD) potential is grown slowly in the desired halo, thus allowing the halo particles to respond adiabatically to the disc’s time-varying potential. In the third stage, the RD is replaced by a live one and the simulation proceeds with live disc (LD) and halo particles. DeBuhr et al. (2012) used such a scheme to introduce stellar discs into dark matter haloes from the Aquarius Project (Springel et al. 2008). They added a RD at a redshift z = 1.3 with a mass parameter for the disc that grew linearly with the scale factor from an initial value of zero to its final value at z = 1. The disc was initially centred on the potential minimum of the halo and oriented so that its symmetry axis pointed along either the minor or major axis of the halo. During the RD phase, the motion of the disc centre of mass was determined from Newton’s 3rd law. To initialize the LD, DeBuhr et al. (2012) approximate the halo potential as a flattened, axisymmetric logarithmic potential and then determine the disc distribution function (DF) by solving the Jeans equations. Yurin & Springel (2015) introduced a number of improvements to this scheme. Most notably, they use galic to initialize the LD (Yurin & Springel 2014). This code is based on an iterative scheme for finding stationary solutions to the collisionless Boltzmann equation. The general idea for iterative codes is to begin with a set of particles that has the desired spatial distribution and some initial guess for the velocity distribution. The velocities are then adjusted so as to achieve stationarity, as measured by evolving the system and computing a certain merit function. In Yurin & Springel (2015), the initial disc was assumed to be axisymmetric with a DF that depended on two integrals of motion, the energy, E, and angular momentum, Lz. One striking, if not puzzling, result from this work is the propensity of the discs to form very strong bars. These bars are especially common in models without bulges even in cases where the disc is submaximal. In this paper, we introduce an improved scheme for inserting a LD in a cosmological halo. In particular, the centre of mass and orientation of the RD are determined by solving the standard equations of rigid body dynamics. Thus, our RD can undergo precession and nutation. The angular and linear velocities of the RD at the end of the growth phase are incorporated into the LD initial conditions. As in Yurin & Springel (2015) we use an axisymmetric approximation for the halo potential when constructing the disc DF. However, our DF is constructed from an analytic function of E, Lz, and the vertical energy Ez, which is an approximate integral of motion used in galactics (Dubinski & Kuijken 1995; Widrow, Pym & Dubinski 2008). By design, the disc DF yields a model whose density has the exponential-sech2 form. And with a three-integral DF, we have sufficient flexibility to model realistic Milky Way like discs. As discussed below, the initial disc DF may be crucial in understanding the formation of the bar. As a demonstration of our method we grow a Milky Way like disc in an approximately 1012 h−1 M⊙ halo from a cosmological zoom-in simulation. We discuss both disc dynamics and the effect our disc has on the population of subhaloes. Discs have been invoked as a means of depleting halo substructure and thus alleviating the Missing Satellite Problem, which refers to the underabundance of observed Milky Way satellites relative to the number of cold dark matter subhaloes seen in simulations (Moore et al. 1999; Klypin et al. 1999). An earlier study by D’Onghia et al. (2010) found that when a disc potential is grown in a Milky Way size cosmological halo, the abundance of substructure in the mass range 107–109 M⊙ was reduced by a factor of 2–3. Similar results were found by Sawala et al. (2017) and Garrison-Kimmel et al. (2017). The organization of the paper is as follows. In Section 2, we outline our method for inserting a LD into a cosmological simulation. We also present results from a test-bed simulation where a disc is inserted into an isolated flattened halo. We next apply our method to a cosmological zoom-in simulation. In Section 3, we focus on disc dynamics and find that the disc develops a bar, spiral structure and a warp. In addition, disc–halo interactions appear to ‘kick’ stars out of the disc and into regions normally associated with the stellar halo. In Section 4, we present our results for the spherically averaged density profile and shape of the dark matter halo as well as the distribution of subhaloes. Particular attention is paid to the sensitivity of these results to the disc insertion scheme. We conclude in Section 5 with a summary and discuss possible applications of this work. 2 INSERTING A STELLAR DISC INTO A COSMOLOGICAL HALO In this section, we detail our method for inserting a live stellar disc into a cosmological simulation. We begin with an overview of our approach and the five main simulations presented in this paper. We then describe some of the more technical aspects of the method. 2.1 Overview of simulation set Our simulations are performed with the N-body component of gadget-3, which is an updated version of gadget-2 (Springel 2005). For the cosmological simulations, 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. For the cosmological parameters, we use the results from Planck 2013 (Planck Collaboration XVI 2014) with h = 0.679, Ωb = 0.0481, Ω0 = 0.306, $$\Omega _\Lambda = 0.694$$, σ8 = 0.827, and ns = 0.962. We begin by simulating a 50 h−1 Mpc box with Nr = 5123 particles, where Nr is the effective resolution, each with a mass of ∼7.9 × 107 h−1 M⊙. We identify a Milky Way like halo in the present-day snapshot, that is, a ∼1012 M⊙ halo which has experienced no major mergers since z = 1 and which has no haloes with 2 h−1 Mpc more than half the mass of the Milky Way analogue. We then run an intermediate zoom-in simulation targeting all particles within 10 virial radii of the low-resolution halo. The initial conditions for this simulation are generated with music (Hahn & Abel 2013), which creates five nested regions from a coarse resolution of Nr = 1283 in the outskirts to an effective Nr = 20483 resolution in the targeted region. After this simulation reaches z = 0, we select all particles within 7.5 virial radii and regenerate initial conditions with one more level of refinement, giving six nested zoom regions, where now, the highest effective resolution is Nr = 40963. Our final halo is composed of approximately 107 particles, each with a mass of 1.54 × 105 h−1 M⊙. The softening lengths were selected using the criteria in Power et al. (2003) with a softening of 719 comoving pc for the highest resolution particles in the final zoom-in simulation. We found that this repeated zoom-in technique results in remarkably little contamination from coarse resolution particles within the targeted halo, giving a clean region of size 1.9 h−1 Mpc at z = 0. The dark matter only (DMO) simulation not only serves as the basis for four simulations with discs but also provides a control ‘experiment’ for our study of the effect discs have on halo properties. In each of our disc simulations, the potential of a RD is introduced at the ‘growing disc’ redshift zg. The mass parameter of the disc is then increased linearly with the scale factor a = 1/(1 + z) from zero to its final value Md at the ‘LD’ redshift zl. As described in Section 2.5, the radial and vertical scale lengths of the disc are also increased between zg and zl to account for the fact that discs grow in scale as well as mass while they are being assembled. In the first of our disc–halo simulations, dubbed MN, we introduce a rigid Miyamoto–Nagai disc (Miyamoto & Nagai 1975), whose gravitational potential is given by   $$\Phi \left(R, \,z\right) = -\frac{G M_{\rm d}}{\left(R^2 + \left(R_{\rm d} + \left(z^2 + z_{\rm d}^2 \right)^{1/2}\right)^{2}\right)^{1/2}}.$$ (1)For this simulation, which was meant to mimic the scheme in D’Onghia et al. (2010), we assume that the centre of the disc tracks the minimum of the halo potential while the orientation of the disc is fixed to be aligned with the z-axis of the simulation box. Note that this is effectively a random direction for the halo. In the remaining three disc simulations, we grow an exponential-sech2 RD potential in our halo with mass and scale length parameters that grow in time. For our fixed-orientation (FO) simulation, the position and velocity of the disc’s centre of mass are determined from Newtonian dynamics while the spin axis of the disc is initially aligned with the minor axis of the halo at z = zg and kept fixed in the simulation box frame thereafter. In this respect, the simulation is similar to the ones presented in DeBuhr et al. (2012) and Yurin & Springel (2015). For the RD simulation the disc’s orientation, which is now a function of time, is determined from Euler’s rigid body equations. In the MN, FO, and RD runs, we continue the simulation to the present epoch (z = 0) with the assumed RD potential where the mass and scale length parameters are held fixed and the position and orientation are calculated as they were during the growth phase. For our final LD simulation, we swap a LD for the RD disc at zl. Thus, the RD and LD simulations are identical prior to zl. All of our discs have a final mass of Md = 7.2 × 1010 M⊙, a final scale radius of 3.7 kpc, and a final scale height of 0.44 kpc. Our simulation parameters can be found summarized in Table 1. Table 1. A summary of the simulation parameters, as discussed in the text. Md is the final disc mass, Rd, 0 is the final disc scale radius, Nd is the number of particles used to simulate the disc, zg and zl are the redshifts when the disc beings to grow and when it becomes live (respectively), Nr is the effective resolution in the zoom-in region, and Lbox is the comoving size of the box.   DMO  MN  FO  RD  LD  Md (M⊙)  –  7.2 × 1010  7.2 × 1010  7.2 × 1010  7.2 × 1010  Rd, 0 (kpc)  –  3.7  3.7  3.7  3.7  Nd  –  –  106  106  106  zg  –  3.0  3.0  3.0  –  zl  –  1.0  1.0  1.0  1.0  Nr  4096  4096  4096  4096  4096  Lbox(Mpc h−1)  50  50  50  50  50    DMO  MN  FO  RD  LD  Md (M⊙)  –  7.2 × 1010  7.2 × 1010  7.2 × 1010  7.2 × 1010  Rd, 0 (kpc)  –  3.7  3.7  3.7  3.7  Nd  –  –  106  106  106  zg  –  3.0  3.0  3.0  –  zl  –  1.0  1.0  1.0  1.0  Nr  4096  4096  4096  4096  4096  Lbox(Mpc h−1)  50  50  50  50  50  View Large 2.2 Summary of live disc insertion scheme The first step in our disc insertion scheme is to calculate an axisymmetric approximation to the gravitational potential of the DMO halo at z = 0. We do so using an expansion in Legendre polynomials as described below. We then generate a particle representation of a stellar disc that is in equilibrium with this potential using the galactics code (Kuijken & Gilmore 1989; Widrow et al. 2008). Though galactics allows one to generate the phase space coordinates of the disc stars, at this stage, we only need the positions of the stars, which we use to represent the ‘RD’. At the zg snapshot, we incorporate the disc particles into the mass distribution of the system with the disc centred on the potential minimum of the halo. We then rerun the simulation from zg to zl with the following provisos. First, the mass of the disc is increased linearly with a from zero to its final value. Secondly, the size of the disc increases with time, which we account for by having the positions of the disc particles, as measured in the disc frame, expand with time to their final values at zl. Finally, the centre of mass and orientation of the disc are determined by integrating the equations of rigid-body dynamics. At redshift zl, the DF of the disc is recalculated assuming the same structural parameters as before but with an axisymmetric approximation to the new halo potential. An N-body disc is generated from this DF and the simulation proceeds with LD particles. In this paper, we choose zg = 3 and zl = 1 so that the growth period lasts from 2.2 to 5.9 Gyr after the big bang. This period in time roughly corresponds to the epoch of peak star formation in Milky Way like galaxies (e.g. van Dokkum et al. 2013). Our simulations during this epoch can be used to study the effect of a disc potential on the evolution of substructure. On the other hand, our simulations of the LD epoch (zl > z > 0) can also be used to study disc dynamics in a cosmologically motivated dark halo. 2.3 Halo potential Our method requires an axisymmetric approximation to the halo potential centred on the disc. This approximation is found using an expansion in spherical harmonics (see Binney & Tremaine 2008) where only the m = 0 terms are included. The potential is then expressed as an expansion in Legendre polynomials. We divide the region that surrounds the disc into spherical shells and calculate the quantities   $$m_{{\rm l},i} = \sum _{n\in S_i} m_n P_{\rm l}(\cos {\theta _n}),$$ (2)where the sum is over the halo particles of mass mn in the ith shell (Si), Pl are the Legendre polynomials, and (r,  θ,  ϕ) are spherical polar coordinates centred on the disc. The axisymmetric approximation to the potential is then   $$\Phi _{\rm h}\left(r, \, \theta \right) = \sum _{l=0}^\infty A_{\rm l}\left(r\right) P_{\rm l}\left(\cos {\theta }\right),$$ (3)where   \begin{eqnarray} A_{\rm l}(r) \!=\! -G\left( \frac{1}{r^l}\int _0^r {\rm d}r^{\prime } r^{\prime l+2} m_{\rm l}(r^{\prime }) \!+\!r^{l+1}\int _0^r {\rm d}r^{\prime } r^{\prime 1-l} m_{\rm l}(r^{\prime })\right), \end{eqnarray} (4)and ml is given by equation (2) for sufficiently small radial bins. 2.4 Disc DFs Armed with an axisymmetric approximation to the halo potential, we construct a self-consistent DF for the disc following the method outlined in Dubinski & Kuijken (1995). This DF is an analytic function of E, Lz, and Ez. By construction, the DF yields a density law for the disc which is, to a good approximation, given by   $$\rho _{\rm d}\left(R, \,z\right)\simeq \frac{M_{\rm d}}{4\pi R_{\rm d}^2 z_{\rm d}} {\rm e}^{-R/R_{\rm d}} {\rm sech}^2\left(z/z_{\rm d}\right) T\left(R_t, \Delta _t\right),$$ (5)where Md, Rd, and zd are the mass, radial scale length, and vertical scale height of the disc and $$R=\sqrt{r^2-z^2}$$. The truncation function T insures that the density falls rapidly to zero at a radius Rt + Δt. The square of the radial velocity dispersion is chosen to be proportional to the surface density, that is, σR = σR0exp ( − R/2Rd). The central velocity dispersion σR0 controls, among other things, the Toomre Q parameter and thus the susceptibility of the disc to instabilities in the disc plane. The azimuthal velocity dispersion is calculated from the radial velocity dispersion through the epicycle approximation (for details see Binney & Tremaine 2008) while the vertical velocity dispersion is adjusted to yield a constant scale height zd. We stress that although the density law is written as a function of R and z, which are not integrals of motion, the underlying DF is a function of E, Lz, and Ez, which are integrals of motion to the extent that the epicycle approximation is valid and that the potential can be approximated by an axisymmetric function. 2.5 Rigid disc dynamics During the disc growth phase, the disc mass is given by   $$M(a) = M_{\rm d} \left( \frac{a - a_{\rm g}}{a_{\rm l} - a_{\rm g}} \right),$$ (6)where ag is the scale factor evaluated at zg. The positions of the disc particles expand self-similarly in disc or body coordinates. That is, the comoving position of a disc particle in body coordinates is given by $${\boldsymbol s}_i(a) = b(a){\boldsymbol s}_i (a_{\rm l})$$ where $${\boldsymbol s}_i(a)$$ is the comoving position of the ith disc particle in the body frame   $$b(a) = b_{\rm g} + \left(1 - b_{\rm g}\right) \left( \frac{a - a_{\rm g}}{a_{\rm l} - a_{\rm g}} \right),$$ (7)where bg = b(ag), and we choose bg = 0.1. The angular velocity of the disc is described by the vector $$\boldsymbol {\omega } = \left(\omega _x, \, \omega _y, \, \omega _z\right)$$ where ωz corresponds to the spin of the disc about its symmetry axis. We assume that   $$\omega _z(a) = \omega _z(a_{\rm l})\left(\frac{M(a)}{M_{\rm d}}\right)^{1/\alpha }\,b(a),$$ (8)which insures that the disc tracks the Tully–Fisher relation, $$M_{\rm d}\propto V_{\rm d}^\alpha \propto \left(\omega _3R_{\rm d}\right)^\alpha$$ (Torres-Flores et al. 2011). In this work we set α = 3.5. The orientation of the disc is described by its Euler angles. We follow the convention of Thornton & Marion (2008) and use ϕ,  θ,  and ψ where the matrix   $${\cal R} = {\cal R}_z(\phi ) {\cal R}_y(\theta ) {\cal R}_z(\psi )$$ (9)describes the transformation from the disc body frame to the simulation frame. Here, $${\cal R}_i(\alpha )$$ is the matrix for a rotation by angle α about the ith axis. Physically, changes in ϕ and θ correspond to precession and nutation, respectively, while ψ is a degenerate rotation about the disc’s symmetry axis. The equations of motion for the Euler angles and angular velocity of the disc, which must account for the time-dependence of the disc’s moment of inertia as well as the fact that gadget-3 uses comoving coordinates, are derived in Appendix A. These equations allow us to solve for the orientation of the disc under the influence of torque due to dark matter. At the end of the growth phase, we initialize the disc with a DF that is recalculated using galactics. As we will see, during the growth phase the motion of the disc involves a mix of precession and nutation. In general, a LD is not able to support the sort of rapid nutation seen in the RD, essentially because different parts of the disc respond to torques from the halo and the disc itself differently. We therefore initialize the LD with an orientation and precessional motion given by a fixed-window moving average of the RD coordinates. 2.6 Test-bed simulation of an isolated galaxy We test our method by growing a stellar disc in an isolated, flattened halo in a non-cosmological simulation. To initialize the flattened halo we first generate a particle realization of a truncated, spherically symmetric NFW halo (Navarro et al. 1997) whose density profile is given by   $$\rho (r) = \frac{v_{\rm h}^2 a_{\rm h}}{4\pi G} \frac{1}{r\left(r + a_{\rm h}\right)^2},$$ (10)with ah = 8 kpc and vh = 400 km s−1. The halo is truncated at a radius much larger than the radius of the disc. The z and vz coordinates are then reduced by a factor of 2 and the system is evolved until it reaches approximate equilibrium. The result is an oblate halo with an axial ratio of ∼0.8 and a symmetry axis that coincides with the z-axis of the simulation box. We next grow a RD over a period of 1 Gyr to a final mass of 4.9 × 1010 M⊙ and final radial and vertical scale lengths of Rd = 2.5 kpc and zd = 200 pc, respectively. The disc is grown at an incline of 30° relative to the symmetry plane of the halo. Doing so allows us to check the rigid body integration scheme for a case when the symmetry axes of the disc and halo are initially misaligned. At t = 1 Gyr we replace the RD with a live one and evolve the system for an additional 1 Gyr. In a separate simulation, we also follow the evolution of the RD over the same time period, allowing us to compare the evolution of the LD and RD. Fig. 1 shows the Euler angles and angular velocity components for the RD and LD as a function of time. The time-dependence of ωx and ωy is characterized by an interference pattern between short ∼125 Myr period nutations and a decaying long-period precessional motion. Note that θ and ϕ for the LD track the corresponding values for the RD for t > 1 Gyr. By initializing the LD with the angular velocity of the RD, we capture the (small) precessional motion of ∼10° Gyr−1 between t = 1 and 2 Gyr. The disc settles into a preferred plane within the first 300 Myr that is intermediate between its initial symmetry plane and the initial symmetry plane of the halo. More precisely, the vector of the new minor axis is $${\boldsymbol c} = [-0.159,0.146,0.976]$$ measured at 20kpc, 12.5° from the original flattening axis. Figure 1. View largeDownload slide Kinematic variables for the RD and LD in an isolated, flattened halo as a function of time. The upper two panels show the Euler angles θ and ϕ for the RD (dashed curves) and LD (solid curves) where the LD is introduced at t = 1 Gyr (red vertical line). The bottom two panels show the x and y components of the angular velocity, as measured in the body coordinate system. In these two panels the solid curves show the δt ∼ 150 My moving average, which is used to initialize the LD. Figure 1. View largeDownload slide Kinematic variables for the RD and LD in an isolated, flattened halo as a function of time. The upper two panels show the Euler angles θ and ϕ for the RD (dashed curves) and LD (solid curves) where the LD is introduced at t = 1 Gyr (red vertical line). The bottom two panels show the x and y components of the angular velocity, as measured in the body coordinate system. In these two panels the solid curves show the δt ∼ 150 My moving average, which is used to initialize the LD. Fig. 2 shows surface density maps for the disc at two snapshots. The disc develops a weak warp due to its interaction with the halo. The development of the warp is also evident in the surface density, vertical velocity dispersion, and scale height profiles shown in Fig. 3. We see that the surface density within ∼15 kpc or 6Rd is essentially unchanged while at larger radii, there are 10–20 per cent time-dependent fluctuations. The scale height 〈z2〉1/2 increases with time and radius. At early times, the increase is most prominent beyond ∼15 kpc while at late times, the scale height increases more smoothly from the centre to the edge of the disc. Figure 2. View largeDownload slide Face-on projections of the particle distribution for two snapshots of a LD in a flattened halo. The solid line for scale is 25 kpc with a centre coincident with the disc. Figure 2. View largeDownload slide Face-on projections of the particle distribution for two snapshots of a LD in a flattened halo. The solid line for scale is 25 kpc with a centre coincident with the disc. Figure 3. View largeDownload slide Surface density, vertical velocity dispersion, and scale height profiles as a function of Galactocentric radius R for 10 snapshots equally spaced in time. The top panel shows the surface density Σ(R) divided by $$\Sigma _0(R) = \left(M_{\rm d}/2\pi R_{\rm d}^2\right)\exp {\left[-\left(R/R_{\rm d}\right) \right]}$$ in order to highlight departures from a pure exponential disc. Likewise, in the middle panel, we show the ratio σz(R)/σz, 0(R) where σz, 0 = exp ( − R/2Rd). The bottom panel shows the RMS z as a function of R. Figure 3. View largeDownload slide Surface density, vertical velocity dispersion, and scale height profiles as a function of Galactocentric radius R for 10 snapshots equally spaced in time. The top panel shows the surface density Σ(R) divided by $$\Sigma _0(R) = \left(M_{\rm d}/2\pi R_{\rm d}^2\right)\exp {\left[-\left(R/R_{\rm d}\right) \right]}$$ in order to highlight departures from a pure exponential disc. Likewise, in the middle panel, we show the ratio σz(R)/σz, 0(R) where σz, 0 = exp ( − R/2Rd). The bottom panel shows the RMS z as a function of R. 3 COSMOLOGICAL SIMULATIONS We now use our method to insert a LD with prescribed structural properties into a cosmological halo. In this section, we focus on disc dynamics while in the next, we consider the effect the disc has on the dark halo. In Fig.  4, we show the kinematic variables for the RD and LD in the RD and LD simulations. The two simulations are identical prior to z = 1 (a = 0.5) when the LD is swapped in for the rigid one. The short-period (300 Myr) oscillations in ω1 and ω2 are nutations. To initialize the LD, we use the fixed-window moving average of ωx and ωy. By and large, the Euler angles of the RD’s and LD’s track one another for z < 1, indicating that the RD is a reasonable model for a live one, at least in terms of the disc’s orientation. Figure 4. View largeDownload slide Kinematic variables for the RD and LD in our cosmological halo as a function of scale factor a. Line types are the same as in Fig. 1. The LD is introduced at a redshift z = 1 when the scale factor is a = 0.5 (red vertical line). The blue line shows the δa ∼ 0.04 moving average calculated by averaging the last 300 points in the disc integration routine. Figure 4. View largeDownload slide Kinematic variables for the RD and LD in our cosmological halo as a function of scale factor a. Line types are the same as in Fig. 1. The LD is introduced at a redshift z = 1 when the scale factor is a = 0.5 (red vertical line). The blue line shows the δa ∼ 0.04 moving average calculated by averaging the last 300 points in the disc integration routine. In Fig. 5, we show the circular speed curves at z = 1 and z = 0 for our four simulations. We see that the disc in our model is submaximal. To be precise, we have Vd/Vc ≃ 0.68 at R = 2.2Rd where Vd is the circular speed due to the disc and Vc is the total circular speed. In short, the contributions from the disc and halo to the centrifugal force are approximately equal at a radius where the disc contribution reaches its peak value. By comparison, a maximal disc is generally defined to have Vd/Vc > 0.85 (Sackett 1997). If we use Vd/Vc at 2.2Rd as the defining characteristic of the model, then our simulations match up with the F-5 simulation of Yurin & Springel (2015), although our discs are slightly warmer, with a Toomre Q parameter of 1.4 as compared with Q ≃ 0.9 for their discs and our discs are thinner (200 pc ]versus 600 pc). We note that in a two-integral disc DF, Q and the disc thickness are linked whereas in a three-integral DF, they can be set independently. The method of Yurin & Springel (2014) can be extended to consider a three-integral DF, but these models were not considered in Yurin & Springel (2015). Moreover, their two-integral model, which imposes σR = σz, violates the epicycle approximation, leading to transient system behaviour at early times when disc bars first form. Figure 5. View largeDownload slide Circular speed curve decompositions at z = 1 (top row) and z = 0 (bottom row) for (from left to right) our MN, FO, RD, and LD simulations. Halo contributions are represented as dot–dashed lines, disc contributions are represented by dashed lines, and the total rotation curve is given by a solid line. For reference, we have included the circular speed curve for the DMO halo (dot–dashed curve). Figure 5. View largeDownload slide Circular speed curve decompositions at z = 1 (top row) and z = 0 (bottom row) for (from left to right) our MN, FO, RD, and LD simulations. Halo contributions are represented as dot–dashed lines, disc contributions are represented by dashed lines, and the total rotation curve is given by a solid line. For reference, we have included the circular speed curve for the DMO halo (dot–dashed curve). The circular speed curves in Fig. 5 show little change exterior to ∼2Rd after zl, thus providing another indication that the LD was close equilibrium when it was swapped in for the rigid one. The formation of a bar is evident in the circular speed, and we can infer the bar contributes substantially inside 2.2Rd ≃ 8 kpc. The halo contribution at R = 2.2Rd is about 20 per cent larger in the four disc runs than in the DMO one due to adiabatic contraction. Interestingly, the halo in the MN run shows somewhat more contraction than in the RD and LD runs. We note that in the MN run, the disc potential tracks the potential minimum of the halo whereas in the RD/LD case, the disc’s position is determined from Newtonian dynamics. In general, the centre of the disc tracks the halo potential minimum so long as the potential changes slowly with time. However, during a major merger (and indeed, just such an event occurs at z = 2) there are rapid changes in the halo potential and the position of the disc, as determined by Newtonian dynamics, can differ significant from the minimum of the halo potential. Evidently, the ad hoc prescription of growing a disc at the halo’s potential minimum may, in some cases, overestimate the effect of adiabatic contraction. 3.1 Bar formation In Fig. 6, we show orthogonal projections of the disc density in our LD simulation at four epochs between z = 1 (lookback time of 7.9 Gyr) and the present epoch. During the first billion years of LD evolution, the disc develops a bar and spiral structure. In addition, there is a warp in the outer disc extending several kiloparsecs above the mid-plane of the inner disc. By the present epoch, the bar has grown in length and intensified and the edge-on view shows the classic X-pattern. Figure 6. View largeDownload slide Projected density along three orthogonal directions for the LD at four epochs between z = 1 and z = 0. The projections are presented in physical units. The solid line for scale is 37 kpc with a centre coincident with the disc. Figure 6. View largeDownload slide Projected density along three orthogonal directions for the LD at four epochs between z = 1 and z = 0. The projections are presented in physical units. The solid line for scale is 37 kpc with a centre coincident with the disc. We consider the usual parameter bar strength A2 = |c2| where   $$c_m = \frac{1}{M_S}\sum _{j\in S} m_j {\rm e}^{im\phi _j}.$$ (11)Here, S is some circularly symmetric region of the disc (e.g. a circular annulus) and the sum is over all particles labelled by j and with mass mj that are inside S. We find that A2 for the inner 2Rd of the disc reaches 0.43 at t = 6.7 Gyr, decreases to 0.36 by t = 9.2 Gyr, presumably because the bar has buckled, and then increases to 0.47 by the present epoch. On the other hand, A2 for the entire disc increases to 0.27, decreases to 0.23, and then increases to 0.28 for the same epochs. Note that the inner 2Rd of the disc contains 60 per cent of the mass. Thus, the fact that $$A_{2,2R_{\rm d}}/A_{2,{\rm tot}}\simeq 0.6$$ implies that most of the bar mass resides within the inner 2Rd. The bars in Yurin & Springel (2015) seem to be stronger then then ones in our simulations – they find A2 ≃ 0.6 but use a non-standard formula for A2. Moreover, their bars appear to extend across most of the disc. In terms of disc dynamics, the main difference between our simulations is the fact that we use a three-integral DF for the disc whereas they use a two-integral DF. In the latter, the velocity dispersion in the radial and vertical directions are the same. Thus, the radial dispersion, which fixes the Toomre Q parameter, also determines the thickness of the disc. We note that their initial discs are a factor of 2 or 3 thicker than ours. We speculate that the bars that develop in these thick discs are less susceptible to buckling and therefore able to grow stronger and longer. These ideas will be investigated in more detail in a future publication. 3.2 Kicked-up stars and disc heating The outer part of the disc suffers considerable disruption and warping presumably through its interaction with the main halo or substructure. The right-most panel of the a = 0.55 snapshot in Fig. 6, for example, shows a classic integral-sign warp. The other snapshots show that a significant number of disc particles have orbits that now take them to high galactic latitudes. The impressions one has from the density projections are borne out in Fig. 7 where we show the surface density and scale height profiles at different times. Bar formation redistributes mass in the disc leaving a deficit (relative to the initial exponential disc) between 5 and 15 kpc. The disc becomes thicker and its vertical velocity dispersion increases though a combination of disc–halo interactions and the effects of the bar and spiral structure (Gauthier et al. 2006; Dubinski et al. 2008; Kazantzidis et al. 2008). Figure 7. View largeDownload slide Surface density, vertical velocity dispersion, and scale height profiles of the LD for 10 snapshots equally spaced in scale factor a between a = 0.5 (z = 1) and a = 1 (present epoch). Panels are the same as in Fig. 3. Figure 7. View largeDownload slide Surface density, vertical velocity dispersion, and scale height profiles of the LD for 10 snapshots equally spaced in scale factor a between a = 0.5 (z = 1) and a = 1 (present epoch). Panels are the same as in Fig. 3. A striking feature of the simulations are the streams of disc stars well above the disc plane. These stars may represent an example of a kicked-up disc, which has been seen in other N-body simulations (Purcell, Bullock & Kazantzidis 2010; McCarthy et al. 2012) and invoked to explain kinematic and spectroscopic observations of M31 (Dorman et al. 2013) and the Monoceros Ring (e.g. Newberg et al. 2002; Ibata et al. 2003). The idea is that interactions between the disc and both satellite galaxies and halo substructure liberate stars from the disc, launching them to regions of the galaxy normally associated with the stellar halo. Our LD simulation corroborates this hypothesis and is in broad agreement with previous numerical work. Finally, we note that the a = 1 panel of Fig. 6 shows a relatively thin, stream-like structure, above the disc which is qualitatively similar to the Anti-centre Stream (ACS; Grillmair 2006). While the ACS is believed to be due to the disruption of a globular cluster (e.g. Grillmair 2006), Fig. 6 suggests that perturbations to the disc can create similar features. Intriguingly, de Boer, Belokurov & Koposov (2018) recently found that the ACS is rotating in the same sense as the Milky Way disc. 4 HALO SUBSTRUCTURE IN THE PRESENCE OF A DISC In this section, we consider the effect of a disc on a halo’s structural properties such as its spherically averaged density profile, its shape, and its subhalo population. An examination of the DMO simulation shows that the halo we have selected builds up through a series of mergers and accretion events, but that by z = 1 it has settled into a relatively relaxed state with an NFW profile that evolves very little between z = 1 and z = 0 within the inner 100 kpc. Our sequence of simulations, (MN, FO, RD, and LD) allow us to tease out the effects of different disc insertion methods. The MN simulation, for example, pins the centre of the disc to the minimum of the halo potential, whereas the other simulations dynamically evolve the position and velocity of the disc potential via Newtonian mechanics. The MN and FO simulations both assume that the orientation of the disc potential during the growth phase is fixed whereas RD and LD solve for the orientation using rigid body dynamics. 4.1 Global properties of the halo In Fig. 8, we show the ratio of the spherically averaged density profile in the four disc runs to that from the DMO run. At z = 1 the haloes in the FO, RD, and MN runs show evidence for adiabatic contraction with the density in the inner ∼30 kpc increasing by a factor of 1.2–2.1. The effect is strongest in the MN simulation, which is to be expected since the halo in that case always sees the disc potential at the minimum of its potential. Of course, this prescription is unphysical. In general, and especially during a major merger, the disc and halo potential minimum will not necessarily coincide. Figure 8. View largeDownload slide The ratio of halo density to the DMO simulation for MN (green), FO (teal), RD (red), and LD (purple) at z = 1 (dashed) and z = 0 (solid). The presence of the disc significantly increases the central concentration of the halo. Figure 8. View largeDownload slide The ratio of halo density to the DMO simulation for MN (green), FO (teal), RD (red), and LD (purple) at z = 1 (dashed) and z = 0 (solid). The presence of the disc significantly increases the central concentration of the halo. Between z = 1 and z = 0, the mass of the disc is constant. Adiabatic contraction ceases but the halo still responds to the time-varying disc potential. Interestingly, at intermediate radii (between ∼10 and 40 kpc) the density profile of the halo settles back to a state close to that found in the DMO run. Perhaps most striking is the fact that the halo in the LD run becomes more centrally concentrated than the halo in any of the other cases. One possible explanation is that dynamical friction from the disc drags dark matter subhaloes towards the centre of the halo where they are tidally disrupted. In Fig. 9, we show the minor-to-major (cr/ar) and intermediate-to-major (br/ar) axial ratios as a function of radius for both the z = 1 and z = 0 snapshots. The axial ratios are calculated by diagonalizing the moment of inertia tensor in linearly spaced radial shells. The DMO halo is triaxial with cr/ar ≃ 0.75 and br/ar ≃ 0.85. Note that the axial ratio profiles are smoother at z = 0 than at z = 1, which supports the observation that the halo has settled into a more relaxed state over the past 7 or so billion years. In general, discs tend to make haloes more spherical, a result that has been known for some time from both collisionless and hydrodynamical simulations (e.g. Dubinski 1994; Zemp et al. 2012). Figure 9. View largeDownload slide Axial ratios as a function of radius. Shown are the minor-to-major axial ratio (top panel) and the intermediate-to-major axial ratio (bottom panel) at z = 1 (dashed curves) and z = 0 (solid curves). Blue corresponds to DMO, green to MN, teal to FO, red to RD, and purple to LD. Figure 9. View largeDownload slide Axial ratios as a function of radius. Shown are the minor-to-major axial ratio (top panel) and the intermediate-to-major axial ratio (bottom panel) at z = 1 (dashed curves) and z = 0 (solid curves). Blue corresponds to DMO, green to MN, teal to FO, red to RD, and purple to LD. Evidently, the MN halo is rounder, especially in the inner part, the FO halo. Recall that the main difference between these two cases is that the MN disc is pinned to the potential minimum of the halo. It is perhaps not surprising then that, as with adiabatic contraction, it has a stronger effect on the halo’s shape. We also note that the axial ratio profiles for the RD and LD simulations are fairly similar. 4.2 Subhalo populations We now turn our attention to halo substructure. We identify subhaloes and determine their positions and masses using rockstar (Behroozi, Wechsler & Wu 2013), which employs a friends-of-friends algorithm in six phase space dimensions. We consider only those subhaloes with mass ms between mmin = 107.5 M⊙ and mmax = 1010.5 M⊙. Subhaloes at the lower end of this range are marginally resolved with ∼100 particles, above which the subhalo mass can be trusted (e.g. Onions et al. 2012), while those at the upper end contain ∼ 3 per cent of the halo’s virial mass. In Fig. 10, we show the cumulative mass in subhaloes as a function of Galactocentric radius:   $$M_{\rm s}\left(<r\right) = \int _0^r {\rm d}r \int _{m_{\rm min}}^{m_{\rm max}} {\rm d}m_{\rm s}\,m_{\rm s}\, \frac{{\rm d}^2 N}{{\rm d}m_{\rm s}\,{\rm d}r}.$$ (12)We also show the cumulative number of subhaloes. In general, the presence of a disc depletes substructure inside about 30 kpc but leaves the outer substructures unaffected. Figure 10. View largeDownload slide Cumulative mass in subhaloes inside a radius r (upper panel) and cumulative number of subhaloes (lower panel). We consider only subhaloes within 500 kpc of the halo centre and with a mass above 107.5 M⊙. The curves are blue (DMO), green (MN), teal (FO), red (RD), and purple (LD). Figure 10. View largeDownload slide Cumulative mass in subhaloes inside a radius r (upper panel) and cumulative number of subhaloes (lower panel). We consider only subhaloes within 500 kpc of the halo centre and with a mass above 107.5 M⊙. The curves are blue (DMO), green (MN), teal (FO), red (RD), and purple (LD). We next consider the differential mass distribution as a function of subhalo mass. Recall that for a pure dark matter halo, $${\rm d}N/{\rm d}\ln (m_{\rm s})\propto m_{\rm s}^{-p}$$ where p ≃ 0.9 (e.g. Gao et al. 2004). In Fig. 11, we therefore show the quantity m0.9 dN/dln (ms) in order to enhance the differences between the different disc runs. We see that the halo population between ms ≃ mmin and ms ≃ 109 M⊙ is depleted, but only by about 20–30  per cent. Taken together, Figs 10 and 11 imply that the main depletion of the subhaloes occurs within the inner regions of the parent halo, in agreement with D’Onghia et al. (2010), Sawala et al. (2017), and Garrison-Kimmel et al. (2017). That said, the depletion of subhaloes seems rather insensitive to the disc insertion method,although we caution that only a single halo was used. Our results, being mainly in agreement with previous work, should still be viewed with caution due the consideration of a single host halo. Figure 11. View largeDownload slide Differential mass distribution multiplied by M0.9 for subhaloes above 107.5 M⊙. We make an outer radius cut at 500 kpc. The curves are blue (DMO), green (MN), teal (FO), red (RD), and purple (LD). Figure 11. View largeDownload slide Differential mass distribution multiplied by M0.9 for subhaloes above 107.5 M⊙. We make an outer radius cut at 500 kpc. The curves are blue (DMO), green (MN), teal (FO), red (RD), and purple (LD). 4.3 Case study: a Sagittarius-like dwarf Observations of the Milky Way’s dwarf galaxies and associated tidal streams provide a potentially powerful probe of the Galactic potential and thus the Galaxy’s dark halo. One of the best-studied examples is the Sagittarius dwarf (Ibata, Gilmore & Irwin 1994). Fortunately, our simulation has a satellite galaxy with similar properties, namely, a dark matter mass of 1.8 × 1010 M⊙ at z = 1.0 and an orbit that takes it close to the disc. We identify this object in the five simulations using the rockstar halo catalogues. We then gather a list of IDs for all the bound particles at an early time before the dwarf is disrupted and follow these same particles in later snapshots using a binary search tree look-up scheme. In Fig. 12, we show the evolution of this subhalo between redshift zl = 1 and z = 0. The first row shows the baseline evolution in our DMO simulation. The dwarf develops leading and trailing tidal tails during the first few billion years. By the present epoch, the tidal debris has dispersed throughout the halo. Figure 12. View largeDownload slide X–Y projections for a selected Sagittarius-like dwarf galaxy. The rows from top to bottom are no disc, a fixed Miyamoto–Nagai disc, a RD, and a LD. The scale factors in columns from left to right are 0.5, 0.55, 0.6, 0.65, and 0.7. The frame edges are 295 kpc on each side. Figure 12. View largeDownload slide X–Y projections for a selected Sagittarius-like dwarf galaxy. The rows from top to bottom are no disc, a fixed Miyamoto–Nagai disc, a RD, and a LD. The scale factors in columns from left to right are 0.5, 0.55, 0.6, 0.65, and 0.7. The frame edges are 295 kpc on each side. The next four rows show the same satellite in our disc simulations. Perhaps the most noticeable result is that there are stronger features in the tidal debris at the present epoch once a disc is included. The detailed morphology of the tidal debris is certainly different from one disc simulation to the next. By eye, debris in MN and FO look somewhat similar as does the debris in RD and LD. Perhaps the most noticeable result is that the tidal debris extends to larger galactocentric radii when a disc is included. The detailed morphology of the tidal debris clearly depends on the disc insertion method. By eye, tidal debris appears more isotropic in MN and FO than in RD and LD. The implication is that fixed potentials are more efficient at disrupting massive satellites than a potential which can respond to the satellite’s presence. However, we have only a single example of massive satellite disruption, and we caution that more examples of this behaviour are needed to test this hypothesis. 5 CONCLUSIONS Simulations in which a stellar disc is inserted ‘by hand’ into a cosmological N-body halo provide a compromise between simulations of isolated disc–halo systems and cosmological simulations that include gasdynamics and star formation. Our method builds on the scheme used by Berentzen & Shlosman (2006) and DeBuhr et al. (2012), and refined by Yurin & Springel (2015). The basic idea is to introduce, at a redshift zg, a RD with zero mass into a halo within a cosmological zoom-in simulation. Between zg and zl the disc is treated as an external potential with a mass and size that increase adiabatically to their present-day values. At zl, the RD is replaced by an N-body one and the simulation proceeds to the present epoch with LD and halo particles. Our method improves upon previous ones in two important ways. First, during the growth phase (zg > z > zl) the position and orientation of the disc evolve according to the standard equations of rigid-body dynamics. Thus, the disc in our scheme receives its linear and angular momentum with the halo in a self-consistent fashion and is therefore able to move, precess, and nutate due to torques from the halo. While previous methods introduced aspects of rigid-body dynamics during the growth phase none appear to have implemented the full dynamical equations have done here (D’Onghia et al. 2010; DeBuhr et al. 2012; Yurin & Springel 2015). Our sequence MN, FO, RD, and LD of simulations highlights where the details of the disc insertion scheme are important and where they are not. For example, schemes in which the disc tracks the minimum of the halo potential tend to overestimate the effects of adiabatic contraction. On the other hand, the effect of the depletion of halo substructure seems to be rather insensitive to the details of how the disc is introduced into the simulation. Disc insertion schemes such as the one introduced in this paper, provide an attractive arena for studies of galactic dynamics. In particular, they allow one to study the interaction between a stellar disc and a realistic dark halo with computationally inexpensive simulations while maintaining some level of control over the structural parameters of the disc. We fully intend to leverage these advantages in future work. Acknowledgements LMW and JSB are supported by a Discovery Grant with the Natural Sciences and Engineering Research Council of Canada. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement no. 308024. DE acknowledges financial support from the ERC. REFERENCES Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013, ApJ , 762, 109 https://doi.org/10.1088/0004-637X/762/2/109 CrossRef Search ADS   Berentzen I., Shlosman I., 2006, ApJ , 648, 807 https://doi.org/10.1086/506016 CrossRef Search ADS   Binney J., Tremaine S., 2008, Galactic Dynamics , 2nd edn. Princeton University Press, Princeton, NJ D’Onghia E., Springel V., Hernquist L., Keres D., 2010, ApJ , 709, 1138 https://doi.org/10.1088/0004-637X/709/2/1138 CrossRef Search ADS   de Boer T. J. L., Belokurov V., Koposov S. E., 2018, MNRAS , 473, 647 CrossRef Search ADS   DeBuhr J., Ma C.-P., White S. D. M., 2012, MNRAS , 426, 983 https://doi.org/10.1111/j.1365-2966.2012.21910.x CrossRef Search ADS   Dorman C. E. et al.  , 2013, ApJ , 779, 103 https://doi.org/10.1088/0004-637X/779/2/103 CrossRef Search ADS   Dubinski J., 1994, ApJ , 431, 617 https://doi.org/10.1086/174512 CrossRef Search ADS   Dubinski J., Kuijken K., 1995, ApJ , 442, 492 https://doi.org/10.1086/175456 CrossRef Search ADS   Dubinski J., Gauthier J.-R., Widrow L., Nickerson S., 2008, in Funes J. G., Corsini E. M., eds, ASP Conf. Ser. Vol. 396, Formation and Evolution of Galaxy Disks . Astron. Soc. Pac., San Francisco, p. 321 Font A. S., Navarro J. F., Stadel J., Quinn T., 2001, ApJ , 563, L1 https://doi.org/10.1086/338479 CrossRef Search ADS   Gao L., White S. D. M., Jenkins A., Stoehr F., Springel V., 2004, MNRAS , 355, 819 https://doi.org/10.1111/j.1365-2966.2004.08360.x CrossRef Search ADS   Garrison-Kimmel S. et al.  , 2017, MNRAS , 471, 1709 CrossRef Search ADS   Gauthier J.-R., Dubinski J., Widrow L. M., 2006, ApJ , 653, 1180 https://doi.org/10.1086/508860 CrossRef Search ADS   Gómez F. A., White S. D. M., Grand R. J. J., Marinacci F., Springel V., Pakmor R., 2017, MNRAS , 465, 3446 CrossRef Search ADS   Grillmair C. J., 2006, ApJ , 651, L29 https://doi.org/10.1086/509255 CrossRef Search ADS   Hahn O., Abel T., 2013, Astrophysics Source Code Library , record ascl:1311.011 Hu S., Sijacki D., 2016, MNRAS , 461, 2789 https://doi.org/10.1093/mnras/stw1463 CrossRef Search ADS   Ibata R. A., Gilmore G., Irwin M. J., 1994, Nature , 370, 194 https://doi.org/10.1038/370194a0 CrossRef Search ADS   Ibata R. A., Irwin M. J., Lewis G. F., Ferguson A. M. N., Tanvir N., 2003, MNRAS , 340, L21 https://doi.org/10.1046/j.1365-8711.2003.06545.x CrossRef Search ADS   Katz N., Quinn T., Bertschinger E., Gelb J. M., 1994, MNRAS , 270, L71 https://doi.org/10.1093/mnras/270.1.L71 CrossRef Search ADS   Katz N., Weinberg D. H., Hernquist L., 1996, ApJS , 105, 19 https://doi.org/10.1086/192305 CrossRef Search ADS   Kazantzidis S., Bullock J. S., Zentner A. R., Kravtsov A. V., Moustakas L. A., 2008, ApJ , 688, 254 https://doi.org/10.1086/591958 CrossRef Search ADS   Klypin A., Kravtsov A. V., Valenzuela O., Prada F., 1999, ApJ , 522, 82 https://doi.org/10.1086/307643 CrossRef Search ADS   Kuijken K., Gilmore G., 1989, MNRAS , 239, 571 https://doi.org/10.1093/mnras/239.2.571 CrossRef Search ADS   Laporte C. F. P., Gómez F. A., Besla G., Johnston K. V., Garavito-Camargo N., 2016, MNRAS , preprint (arXiv:1608.04743) McCarthy I. G., Font A. S., Crain R. A., Deason A. J., Schaye J., Theuns T., 2012, MNRAS , 420, 2245 https://doi.org/10.1111/j.1365-2966.2011.20189.x CrossRef Search ADS   Miyamoto M., Nagai R., 1975, PASJ , 27, 533 Moore B., Ghigna S., Governato F., Lake G., Quinn T., Stadel J., Tozzi P., 1999, ApJ , 524, L19 https://doi.org/10.1086/312287 CrossRef Search ADS   Navarro J. F., Frenk C. S., White S. D. M., 1994, MNRAS , 267, L1 https://doi.org/10.1093/mnras/267.1.L1 CrossRef Search ADS   Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ , 490, 493 https://doi.org/10.1086/304888 CrossRef Search ADS   Newberg H. J. et al.  , 2002, ApJ , 569, 245 https://doi.org/10.1086/338983 CrossRef Search ADS   Onions J. et al.  , 2012, MNRAS , 423, 1200 https://doi.org/10.1111/j.1365-2966.2012.20947.x CrossRef Search ADS   Oñorbe J., Garrison-Kimmel S., Maller A. H., Bullock J. S., Rocha M., Hahn O., 2014, MNRAS , 437, 1894 https://doi.org/10.1093/mnras/stt2020 CrossRef Search ADS   Pakmor R., Springel V., 2013, MNRAS , 432, 176 https://doi.org/10.1093/mnras/stt428 CrossRef Search ADS   Planck Collaboration XVI, 2014, A&A , 571, A16 CrossRef Search ADS   Power C., Navarro J. F., Jenkins A., Frenk C. S., White S. D. M., Springel V., Stadel J., Quinn T., 2003, MNRAS , 338, 14 https://doi.org/10.1046/j.1365-8711.2003.05925.x CrossRef Search ADS   Purcell C. W., Bullock J. S., Kazantzidis S., 2010, MNRAS , 404, 1711 Purcell C. W., Bullock J. S., Tollerud E. J., Rocha M., Chakrabarti S., 2011, Nature , 477, 301 https://doi.org/10.1038/nature10417 CrossRef Search ADS PubMed  Quinn T., Katz N., Stadel J., Lake G., 1997, eprint (arXiv:e-prints) Roškar R., Debattista V. P., Brooks A. M., Quinn T. R., Brook C. B., Governato F., Dalcanton J. J., Wadsley J., 2010, MNRAS , 408, 783 https://doi.org/10.1111/j.1365-2966.2010.17178.x CrossRef Search ADS   Sackett P. D., 1997, ApJ , 483, 103 https://doi.org/10.1086/304223 CrossRef Search ADS   Sawala T., Pihajoki P., Johansson P. H., Frenk C. S., Navarro J. F., Oman K. A., White S. D. M., 2017, MNRAS , 467, 4383 https://doi.org/10.1093/mnras/stx360 CrossRef Search ADS   Schaye J. et al.  , 2015, MNRAS , 446, 521 https://doi.org/10.1093/mnras/stu2058 CrossRef Search ADS   Sellwood J. A., 2013, Dynamics of Disks and Warps . Springer Science+Business Media, Dordrecht, p. 923 Springel V., 2005, MNRAS , 364, 1105 https://doi.org/10.1111/j.1365-2966.2005.09655.x CrossRef Search ADS   Springel V., Hernquist L., 2003, MNRAS , 339, 289 https://doi.org/10.1046/j.1365-8711.2003.06206.x CrossRef Search ADS   Springel V. et al.  , 2008, MNRAS , 391, 1685 https://doi.org/10.1111/j.1365-2966.2008.14066.x CrossRef Search ADS   Stinson G., Seth A., Katz N., Wadsley J., Governato F., Quinn T., 2006, MNRAS , 373, 1074 https://doi.org/10.1111/j.1365-2966.2006.11097.x CrossRef Search ADS   Thornton S. T., Marion J. B., 2008, Classical Dynamics of Particles and Systems , 5th edn. Cengage Learning Torres-Flores S., Epinat B., Amram P., Plana H., Mendes de Oliveira C., 2011, MNRAS , 416, 1936 https://doi.org/10.1111/j.1365-2966.2011.19169.x CrossRef Search ADS   van Dokkum P. G. et al.  , 2013, ApJ , 771, L35 https://doi.org/10.1088/2041-8205/771/2/L35 CrossRef Search ADS   Vogelsberger M., Genel S., Sijacki D., Torrey P., Springel V., Hernquist L., 2013, MNRAS , 436, 3031 https://doi.org/10.1093/mnras/stt1789 CrossRef Search ADS   Widrow L. M., Pym B., Dubinski J., 2008, ApJ , 679, 1239 https://doi.org/10.1086/587636 CrossRef Search ADS   Yurin D., Springel V., 2014, Astrophysics Source Code Library , record ascl:1408.008 Yurin D., Springel V., 2015, MNRAS , 452, 2367 https://doi.org/10.1093/mnras/stv1454 CrossRef Search ADS   Zemp M., Gnedin O. Y., Gnedin N. Y., Kravtsov A. V., 2012, ApJ , 748, 54 https://doi.org/10.1088/0004-637X/748/1/54 CrossRef Search ADS   APPENDIX A: EULER’S EQUATIONS IN COMOVING COORDINATES The time-evolution of the angular momentum vector $${\boldsymbol L}$$ of a rigid body acted upon by a torque $$\boldsymbol {\tau }$$ is given by   $$\left(\frac{{{\rm d}} {\boldsymbol L}}{{{\rm d}} t}\right)_f = \left(\frac{{{\rm d}} {\boldsymbol L}}{{{\rm d}} t}\right)_{\rm b} + \boldsymbol \omega \times {\boldsymbol L} = \boldsymbol {\tau },$$ (A1)where the subscripts f and b denote the frame of the simulation box and the body frame, respectively. In physical coordinates, $${\boldsymbol L} = {\boldsymbol r \times p}$$. Alternatively, we can write $${\boldsymbol L} = {\boldsymbol s\times q}$$ where $${\boldsymbol s} = a^{-1}{\boldsymbol r}$$ refers to comoving coordinates and $${\boldsymbol q} = a^2\dot{\boldsymbol s}$$ is the conjugate momentum to $${\boldsymbol s}$$ (see Quinn et al. 1997). For a rigid body, the components of the angular momentum are given by Li = Iijωj where i, j run over x,  y,  z and there is an implied sum over j. Since gadget-3 uses comoving coordinates, we write Iij = a2Jij where J is the moment of inertia tensor written in terms of the comoving coordinates, $${\boldsymbol s}$$, rather than the physical coordinates, $${\boldsymbol r}$$. For convenience, we define a ‘comoving’ angular velocity $$\boldsymbol {\varpi } = a^{-2}\boldsymbol {\omega }$$. We then have Li = Jijϖj. Note that because of the symmetry of our disc, the moment of inertia tensor is diagonal with Jxx = Jyy = Jzz ≡ J/2. The equations of motion for the Euler angles and the disc angular velocity are then given by the standard Euler equations of rigid body dynamics, modified to account for the time-dependence of the disc’s moment of inertia:   $$\frac{{\rm d}\phi }{{\rm d}t} = a^{-2}\sin ^{-1}{\theta } \left(\varpi _x\sin (\psi ) + \varpi _y \cos (\psi )\right),$$ (A2)  $$\frac{{\rm d}\theta }{{\rm d}t} = a^{-2}\left(\varpi _1\cos (\psi ) - \varpi _y \sin (\psi )\right),$$ (A3)  $$J\frac{\varpi _x}{{\rm d}t} + \varpi _x\frac{{\rm d}J}{{\rm d}t} + J\varpi _y\varpi _z = \tau _x,$$ (A4)and   $$\frac{\varpi _y}{{\rm d}t} + \varpi _y\frac{{\rm d}J}{{\rm d}t} - J\varpi _x\varpi _z = \tau _y.$$ (A5)We have omitted the equations for ψ (rotations in the body frame about the symmetry axis) and ϖz since these are determined directly from equation (8). © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society

### Journal

Monthly Notices of the Royal Astronomical SocietyOxford University Press

Published: May 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations