GalactICS with Gas

TL;DR

This paper introduces an enhanced galaxy modeling code capable of creating self-consistent equilibrium models with multiple components, including gas, for realistic galaxy simulations, and demonstrates its effectiveness through various tests.

Contribution

The new GalactICS version incorporates gas dynamics and live halos, enabling more realistic initial conditions for galaxy evolution simulations.

Findings

Models maintain structural properties over 1.5 Gyr.

Bar formation occurs in all disc components at ~1 Gyr.

The simulated bar length matches the Milky Way's observed bar.

Abstract

We present a new version of the GalactICS code that can generate self-consistent equilibrium galaxy models with a two-component stellar disc and a gas disc as well as a centrally-concentrated bulge and extended dark halo. The models can serve as initial conditions for simulations of isolated galaxies that include both hydrodynamics and collisionless dynamics. We test the code by evolving a pair of simple gas disc-halo models, which differ only in the initial temperature of the gas component. The models are similar to the ones considered in the Wang et al. except that here, the halo is live whereas they included the halo as a fixed potential. We find that the basic structural properties of the models, such as the rotation curve and surface density profiles, are well-preserved over 1.5 Gyr. We also construct a Milky Way model that includes thin and thick stellar disc components, a gas disc, a bulge, and a dark halo. Bar formation occurs in all disc-like components at about 1 Gyr. The bar is strongest in the thin disc while the gas disc contains the most prominent spiral features. The length of the bar in our model is comparable to what has been inferred for the Galactic bar.

Equations

$$\begin{array}{rl}f(E,L_z,E_z)=& f_b(E)+f_h(E)+f_{d,1}(E,L_z,E_z)\\ & +f_{d,2}(E,L_z,E_z)+f_g(E,L_z,E_z),\end{array}$$

$${\nabla }^2\Phi =4\pi \rho (R,z,\Phi )=4\pi \int d^{3}vf\left(E,L_z,E_z\right).$$

$$f_i(ϵ)=\frac{1}{\sqrt{8{\pi }^2}}{\int }_E^0\frac{d^2{\tilde{\rho }}_i}{d{\Phi }^2}\frac{d\Phi }{\sqrt{\Phi -E}},$$

$${\sigma }_z^2\simeq {\sigma }_{z,0}^2{e}^{-R/R_d}$$

$${\sigma }_R^2(R)\simeq {\sigma }^2{e}^{(-R/R_d)}.$$

$${\sigma }_{\varphi }(R)\simeq {\sigma }_R(R)\frac{\kappa }{2\omega },$$

$${\rho }_g(R,z)={\rho }_0(R)\text{exp}\left(-\frac{{\Phi }_z(R,z)}{(\gamma -1)ϵ}\right),$$

$$ϵ=\frac{1}{\gamma -1}\frac{k_{B}T}{\mu m_p},$$

$$c_0=\frac{k_{B}T}{\mu m_p},$$

$$\Sigma (R)={\Sigma }_0{e}^{-R/R_g}={\int }_{-\infty }^{\infty }{\rho }_g(R,z)dz,$$

$${\rho }_0(R)=\frac{{\Sigma }_0{e}^{-R/R_g}}{{\int }_{-\infty }^{\infty }{e}^{[-{\Phi }_z/c_0]}dz}.$$

$$V^2(R,z)=R{\frac{\partial \Phi }{\partial R}|}_{z=0}+c_0{\frac{\partial \ln \rho }{\partial \ln R}|}_{z=0},$$

Full text

††pubyear: 2018††pagerange: GalactICS with Gas–References

1 Introduction

A fundamental problem in galactic astronomy is the construction of self-consistent equilibrium models for individual galaxies. Such models can be used as a template for interpreting observations and inferring the gravitational potential of a galaxy and the structure of its dark halo (Widrow et al., 2008; Taranu et al., 2017). Furthermore, equilibrium models can provide initial conditions (ICs) for N-body simulations, which can in turn be used to study dynamical processes such as the formation of bars and spiral structure (Fux, 1997; Athanassoula & Misiriotis, 2002). Not only are these processes interesting in and of themselves, they also provide additional constraints on the models. For example, an equilibrium model that satisfies observational constraints but is nevertheless unstable to the formation of a strong bar is not a suitable model for \textcolorblacka galaxy that has a weak bar or no bar at all (Sellwood, 1985; Widrow et al., 2008; Randriamampandry et al., 2018).

It is invariably easier to model the evolution of collisionless matter (stars and dark matter) than gas. The dynamics of collisionless matter is determined solely by the gravitational field, \textcolorblackwhich can be calculated efficiently using particle mesh or tree codes. On the other hand, the dynamics of collisional matter depends on both the gravitational field and additional physics including hydrodynamical forces, star formation, feedback, and turbulence. Furthermore, many of these processes occur at scales below the resolution limit of simulations.

The complications associated with gas components extend to setting up ICs, especially, if one aims to follow a system from some equilibrium (but possibly unstable) state. ICs for collisionless systems in general require a self-consistent model for the phase space distribution function (DF) and the gravitational potential. Collisional systems also require models for the temperature of the gas as a function of position as well as the equation of state.

Springel et al. (2005) described a method for initializing a system with gas and stellar \textcolorblackdiscs as well as a dark halo. The DF for the collisionless components are found by solving the Jeans equations (see (Hernquist, 1993) and discussion below). For the gas \textcolorblackdisc, they assumed a simple equation of state $P=P(\rho)$. The vertical structure is then determined by solving the equations for hydrostatic equilibrium \textcolorblackin the direction normal to the \textcolorblackdisc plane.

This method was refined Wang et al. (2010), who were particularly interested in initializing \textcolorblackdiscs in their adaptive mesh refinement (AMR) simulations. In brief, they use an iterative approach to solve a set of equations describing the potential, surface density, and scale height of an isothermal, approximately exponential \textcolorblackdisc in the presence of an external potential. However, their simulations were limited to a gas \textcolorblackdisc in a fixed potential. In this paper we describe how to combine the GalactICS (Galaxy Initial ConditionS) code (Kuijken & Dubinski, 1995; Widrow, & Dubinski, 2005; Widrow et al., 2008), which generates collisionless bulge-\textcolorblackdisc-halo systems, with the Wang et al. (2010) method of generating gas \textcolorblackdiscs111The new version of GalactICS is available upon request..

The problem of generating equilibrium ICs for collisionless systems, though more straightforward than that for collisional ones, is still a non-trivial task. The problem amounts to finding the DF for the collisionless components that satisfy (at least approximately) the coupled collisionless Boltzmann and Poisson equations. In the Hernquist method (Hernquist, 1993), the velocity distribution of the stellar and dark matter components are assumed to be \textcolorblackGaussian with a dispersion found by solving the Jeans equations. Since solutions to the Jeans equations are not true solutions to the collisionless Boltzmann equation, these models typically relax to a state different from the one specified in the ICs (Kazantzidis et al., 2004). Other methods \textcolorblackinclude mkgalaxy (McMillan & Dehnen, 2007), which uses a guided relaxation algorithm, GalIC (Yurin & Springel, 2014), whose algorithm has elements of the Schwarzchild method (Schwarzchild, 1979) of building orbit libraries, and the iterative method presented in Rodionov et al. (2009); Rodionov & Athanassoula (2011).

In contrast with the Hernquist method, the GalactICS DFs are based on the Jeans theorem, which states that functions of the integrals of motion are equilibrium solutions to the collisionless Boltzmann equation (CBE). In particular, the DFs for the bulge and halo are assumed to be functions of the energy, which are constructed to yield the \textcolorblacktarget density profiles via the Eddington inversion formula Binney & Tremaine (2008). The DF for the \textcolorblackdisc is a function of the energy, angular momentum about the symmetry axis, and the vertical energy. Since the latter is only approximately conserved, the DFs are not exact solutions to the CBE. However, for thin \textcolorblackdiscs, the approximation is excellent, and the model ICs for even thick \textcolorblackdiscs are are relatively stable. Armed with DFs in terms of the integrals of motion, GalactICS solves the Poisson equation through an iterative algorithm, \textcolorblackwhich meshes well with the Wang et al. (2010) prescription for building an equilibrium gas \textcolorblackdisc. The code has been used to study bar and spiral structure formation in the Milky Way (MW) (Widrow et al., 2008; Fujii et al., 2019), bending waves in \textcolorblackdiscs (Chequers & Widrow, 2017), how radial dispersions and the Galactic bar influence stellar populations (Debattista et al., 2017), merging dwarf galaxies (Łokas et al., 2014), and more.

Recently, a new method based on angle-action variables called AGAMA (Vasiliev, 2019) has become available. It also involves an iterative scheme to solve Poisson’s equation with the additional step that for each iteration, one must construct that action variables in the new potential. One advantage of the method over GalactICS is that it avoids the approximate "third integral", $E_{z}$ and therefore should be able to do a better job of modelling warm \textcolorblackdiscs. On the other hand, the connection between angle-action variables and the usual phase space coordinates of less transparent. Thus, it can be more difficult to build models with specific structural properties, such as a stellar \textcolorblackdisc with constant scale height (see, for example, (Chequers et al., 2018)). In principle, it should be straight forward to add a gas \textcolorblackdisc to an AGAMA model using a method similar to the one described in this paper.

\textcolor

blackYet another approach can be found in Rodionov & Athanassoula (2011) who introduced an alternative iterative approach that allows one to build non-equilibrium galaxy models, such as \textcolorblackdiscs embedded in triaxial halos. In their scheme, the system is evolved on a short time scale via the standard dynamical equations. Structure properties or parameters of the system, such as the surface density, are then reset and the system is again allowed to evolve. Ultimately, the sequence produces a system that is close to equilibrium.

The outline of the paper is as follows: In Section 2 We present the details of this new version of GalactICS. We then examine two gas-halo models in Section 3. The models are similar to the ones consideblack in Wang et al. (2010) except in our simulations, the halo is ’live’. In Section 4, we construct a 5-component MW model based on the best-fit model from McMillan (2017). We conclude with a discussion of our results in Section 5.

2 GalactICS Models

The heart of GalactICS is the construction of a DF for each of the components that together yield an axisymmetric model with the desired structural and kinematic properties. \textcolorblackThe total DF is given by

[TABLE]

where $E$ is the energy, $L_{z}$ is the angular momentum about the symmetry axis of the system, and $E_{z}$ is the energy of the vertical motions in the \textcolorblackdiscs. For a time-independent and axisymmetric system, $E$ and $L_{z}$ are conserved while $E_{z}$ is only approximately conserved for \textcolorblackdisc particles on nearly circular orbits.

\textcolor

blackThe density $\rho$ is determined by the integral of $f$ over all velocities and since $f$ is a function of $E$, which depends on $\Phi$, Poisson’s equation is an implicit function of the $\Phi$:

[TABLE]

Note that here and throughout, we set Newton’s constant $G=1$. GalactICS numerically solves Eq. 2.2 using an iterative approach to obtain a self-consistent density-potential pair. First a target density profile is selected for each component. (The gas \textcolorblackdisc is treated somewhat differently. See Sec. 2.2). The corresponding potential is found by solving Poisson’s equation \textcolorblackusing an expansion in Legendre polynomials. The result is used, in turn, to calculate a new density and the process is repeated until the density-potential pair converges.

In detail, most of the algorithm is identical to that given in Widrow et al. (2008). By design, the bulge has a Sérsic density profile, the \textcolorblackdisc components are truncated exponential-sech2 profiles, and the halo is a truncated double-power law. As with Widrow et al. (2008), the bulge and halo use the Abel integral transformation (Binney & Tremaine, 2008) to get their DF,

[TABLE]

where $i$ is either the bulge or halo. In the presence of any of the \textcolorblackdisc components, this equation is solved using a spherical approximation of the \textcolorblackdisc potential for $\Phi_{tot}$. It is important to note that this approximation is not used for the calculation of $\Phi(R,z)$ or $\rho(R,z)$, but is only used for Eq. 2.3. It is also worth noting that this method of evaluating Eq. 2.3 causes a degree of flattening, but the spheroid components remain isotropic and axisymmetric.

2.1 Stellar \textcolorblackdiscs

The flattened nature of the stellar \textcolorblackdiscs demands a different approach than the simple solutions used for the bulge and halo components. GalactICS uses the DF presented in Kuijken & Dubinski (1995), which itself is based on Shu (1969) and Binney (1987).

The potential for this flattened system is calculated using a combination of spherical harmonics and an analytic ’fake’ density-potential pair (Kuijken & Dubinski, 1995). This pair, $(\rho_{fd},\Phi_{fd})$, has the property that $\rho_{d}=\rho_{fd}+\rho_{r}$ and $\Phi_{d}=\Phi_{fd}+\Phi_{r}$ where $\rho_{r}$ and $\Phi_{r}$ are residuals. The ’fake’ components are designed to account for the higher order moments of the total potential, while the lower order moments are calculated by solving Poisson’s equation for the residuals using a small number of $l$ moments. The two are then summed together to give the total potential of the \textcolorblackdisc given some density. In practice, we find that $l_{max}=10$ is sufficient for most models.

The GalactICS DFs for the \textcolorblackdisc components are constructed to yield a vertical structure that is approximately isothermal. Thus, the vertical velocity dispersion is related to the thickness of the \textcolorblackdisc and the surface density. The DF is designed to yield a scale height that is approximately constant across the \textcolorblackdisc and a surface density that is approximately exponential with a scale radius $R_{d}$. Thus, the vertical dispersion profile is given by

[TABLE]

The radial dispersion profile appears as a free function in the \textcolorblackdisc DF. Here, we assume

[TABLE]

\textcolor

blackThis is based on the observations of Bottema (1993). The tangential dispersion is found through the epicycle approximation:

[TABLE]

where $\omega$ is the angular frequency and $\kappa$ is the epicyclic frequency. Note that in the actual DF, $\sigma_{R}$, $\kappa$, and $\omega$ are written as functions of the guiding radius $R_{c}$, which is a function of $L_{z}$ and hence an integral of motion. The approximations in the above equations reflect the fact that $R_{c}$ is only approximately equal to $R$.

2.2 Gas \textcolorblackdisc

Wang et al. (2010) introduced two methods for generating isothermal equilibrium gas \textcolorblackdiscs in a general galactic potential. We have adapted their ’potential’ method. Since the \textcolorblackdisc is isothermal, the scale height increases as a function of the radius. We assume a target exponential surface density for gaseous component, which causes the space density in the midplane to be a decreasing function of radius.

Following Wang et al. (2010), we assume that the density of the gas \textcolorblackdisc is given by

[TABLE]

where $\rho_{0}(R)$ is the mid-plane density, $\gamma$ is the adiabatic index, $\epsilon$ is the specific internal energy, and $\Phi_{z}(R,z)=\Phi(R,z)-\Phi(R,0)$. The specific internal energy depends on the temperature through the equation of state

[TABLE]

where $k_{B}$ is the Boltzmann constant, $\mu$ is the atomic weight of the gas, and $m_{p}$ is the mass of the proton. For simplicity, that gas is assumed to be hydrogen and adiabatic so $\mu=1$, $\gamma=5/3$, and $\epsilon$ only depends on the gas temperature. It is convenient to re-parameterize the gas temperature as\textcolorblack

[TABLE]

which is similar to the specific internal energy and related to the sound speed. Combining the target exponential surface density,

[TABLE]

with $\rho_{g}(R,z)$ gives

[TABLE]

The result is then included in Poisson’s equation (Eq. 2.2). The algorithm for solving this equation proceeds as follows: First an initial surface density profile is used to estimate the gas \textcolorblackdisc potential. That potential is then used to calculate the scale height profile and total density, which are in turn used to get a new surface density and potential. The whole process is repeated until the system converges. Note that the surface density of the final system may deviate slightly from the target exponential surface density. As with the stellar \textcolorblackdisc, we use an analytic fake density-potential pair when solving Poisson’s equation.

The gas \textcolorblackdisc density-potential pair is constructed at the same time as the total density-potential pair. The full GalactICS algorithm for obtaining a self-consistent density potential pair is then given as follows:

Define the target profiles for all components. 2. 2.

Estimate the total potential using Poisson’s equation. 3. 3.

Calculate new densities for the collisionless components using their DFs. 4. 4.

Calculate the gas density using the current scale height profile and total potential. 5. 5.

Calculate new gas surface density and scale height profiles. 6. 6.

Repeat steps ii-v to achieve convergence.

We use a greatly simplified velocity structure for the gas \textcolorblackdisc compared to the stellar \textcolorblackdiscs. The gas temperature is treated as a kinetic temperature that accounts for all turbulent motions. Therefore, the particles are initialized on \textcolorblackpurely rotational orbits where the speed is

[TABLE]

with no random motions. \textcolorblackNear the center it is possible for the gas pressure term to set $V^{2}<0$. While this condition has not occurred in any of the subsequent simulations, GalactICS currently sets $V^{2}=0$ if the gas pressure term yields a negative result. In such cases we set $V^{2}=0$. An alternate method is to use the softening procedure outlined in Hernquist (1993). At this point we have not yet implemented and tested this in the context of GalactICS. In addition, GalactICS sets $v_{z}=0$ for all gas particles as the vertical gravitational force is balanced against the gas pressure.

While this section has focused on the generation of a gas \textcolorblackdisc within the GalactICS methodology, it is worth noting that this same algorithm can be included in the AGAMA method with little in the way of modifications since AGAMA also uses an iterative approach to calculate the total potential and density.\textcolorblack However, since GalIC, and the Hernquist method do not include an iterative adjustment of the total potential, it would be difficult to modify those algorithms to utilize this \textcolorblackparticular method of constructing a gas \textcolorblackdisc. It would also be difficult to implement this method of gas disc construction into mkgalaxy due to differences in how the iterative calculation of the potential is performed. Nonetheless, gas discs could be generated in all these algorithms using alternate methods than the process utilized here.

3 Gas Only Models

In this section we describe testbed simulations that involve a gas \textcolorblackdisc and dark halo. The ICs roughly correspond to the Gas0 and Gas4 models of Wang et al. (2010). They have the same structural properties for the two components and differ only in the gas temperature. In particular, the dark halo has an NFW profile with scale length $r_{h}=17.8\,{\rm kpc}$ and velocity scale parameter (in GalactICS parameterization) $\sigma_{h}=389\,{\rm km\,s}^{-1}$. This velocity scale parameter implies an NFW density parameter of $\rho_{0}=\sigma_{h}^{2}\left(2\pi r_{h}\right)^{-2}\simeq 1.77\times 10^{7}\,M_{\odot}\,{\rm kpc}^{-3}$. The gas \textcolorblackdisc has a mass of $M_{d}=10^{10}\,\mathrm{M}_{\odot}$ and an exponential scale length $R_{g}=3.5$ kpc. The resultant rotation curve rises slowly to about $180\,{\rm km\,s}^{-1}$ at a radius of $30$ kpc. The gas contribution to the rotational speed is always subdominant indicating that the \textcolorblackdisc will be stable to global perturbations. Following Wang et al. (2010) we set the gas temperature in the GalactICS-Gas0 model to $T=4\times 10^{4}$ K and the temperature in the GalactICS-Gas4 model to $T=8\times 10^{3}$ K. With these values, we expect the GalactICS-Gas0 model to be stable to local perturbations and the colder GalactICS-Gas4 model to be unstable (see Figure 4 of Wang et al. (2010)).

The DFs for the models are sampled with $1$M gas particles and $5$M halo particles and evolved for 1.6 Gyr using the Gadget 2 N-body code (Springel, 2005). \textcolorblackWe use an adaptive softening length with an initial length set to 0.005 kpc and a maximal softening lengths of 0.5 kpc. This is smaller than the recommended softening lengths from Rodionov & Sotnikova (2005). However, the gas \textcolorblackdiscs for both models are quite thin in the central kpc. As such, we feel that this small scale length is justified. In the simulation, we use an adaptive time step with a maximal step of 0.01 Gyr (which also matches the snapshot frequency). The Courant factor in the Gadget 2 simulation is set equal to 0.25.

We stress that in our simulations, the halo is ’live’ whereas the halo in the Wang et al. (2010) simulations is treated as a static halo potential. Moreover, Wang et al. (2010) ran their simulations with the adaptive-mesh-refinement code RAMSES (Teyssier, 2002) whereas Gadget 2 models the gas using smooth particle hydrodynamics.

3.1 The GalactICS-Gas0 model

Figure 1 shows the initial and final surface and cross-sectional densities for the GalactICS-Gas0 model. It is clear that there has been some evolution of the system. In particular, the central surface density of the \textcolorblackdisc has increased (upper panels), as has its scale height (lower panels).

The changes in the gas surface density and thickness are highlighted in Figs. 2 and 3 respectively. In Fig. 2 we see that the surface density has increased significantly in the central regions. In addition, the sharp cut-off at the edge of the \textcolorblackdisc has been smoothed out. On the other hand, at intermediate radii ($1<R<15\,{\rm kpc}$), the evolved surface density is consistent with the initial one. The \textcolorblackdisc thickness increases at all radii, with the largest increase occurring in the central $0.5$ kpc.

Further clues as to the evolution of the system can be found in Fig. 4, which shows that the azimuthally averaged rotation curve \textcolorblack(calculated from the gravitational potential) as well as a scatter plot of gas particle azimuthal velocities. \textcolorblackAlthough the gas \textcolorblackdisc is initialized with purely rotational motions, it is clear from the figure that the gas particles gain random azimuthal (and presumably radial and vertical) motions. These random motions represent a transfer of rotational energy to the random kinetic energy of the gas particles. It is therefore not surprising that the \textcolorblackdisc "puffs up" a bit and also compresses slightly in the radial direction, thereby increasing the central surface density. Note that the \textcolorblackgravitational rotation curve stays very nearly constant over the duration of the simulation.

To investigate the transients in this simulation we show, in Fig. 5, the face-on and edge-on projections of the density after just 50 Myr. A ring of particles corresponding to a axisymmetric density wave, is clearly visible.

Apart from the scattering of \textcolorblackdisc particles, a second possible cause of the transients may be in the use of the ’fake’ gas \textcolorblackdisc density-potential pair. The gas \textcolorblackdisc is extremely thin near its centre with a scale height of less than 100 pc. Moreover, the thickness changes rapidly with radius there. Thus, it may be difficult for our fake-\textcolorblackdisc/Legendre polynomial expansion scheme to full capture the potential. Small errors in the potential solver can lead to an incorrect scale height and an over/under pressure, thereby causing a transient wave, as is seen in Fig. 5.

To be clear the approximation that causes the transient only fails in regions where the value of $\Phi_{z}$ is dominated by the gas \textcolorblackdisc contribution. At moderate radii the halo dominates $\Phi_{z}$ and the small errors from the gas potential approximation are inconsequential. To illustrate this point, we generated a Gas0-stars model where the gas \textcolorblackdisc is replaced by a two-component \textcolorblackdisc of gas and stars with a mass ratio of $1:3$. The velocity dispersion for the stars is assumed to be exponential dispersion as in Eq. 2.10 with a central dispersion of $\sigma_{r,0}=100\mathrm{km~{}s}^{-1}$. The right panel of Fig. 5 shows the gas \textcolorblackdisc of this system at 50 Myr. In contrast with the GalactICS-Gas0 model, there is no transient wave even though the system has the same surface density profile. It is worth noting that the thickness of the stellar \textcolorblackdisc also thickens the gas \textcolorblackdisc in the central regions while decreasing the degree of flaring in the outer regions.

3.2 The GalactICS-Gas4 Model

We next consider the GalactICS-Gas4 model, which is unstable to local perturbations. In the Wang et al. (2010) realization of this model, the gas \textcolorblackdisc rapidly fragments (see the lower left panel of their Figure 5). \textcolorblackOur realization of this model is also unstable, but it evolves slightly differently.

\textcolor

black Fig. 6 shows our gas \textcolorblackdisc at 50 Myr and 1 Gyr. The initially very thin \textcolorblackdisc thickens significantly by 1 Gyr. At 50 Myr, the \textcolorblackdisc has begun fragmenting into a strong spiral structure. This fragmentation is less symmetric than that seen in the Wang et al. (2010) model due to both the random sampling procedure of GalactICS and the live halo. By 1 Gyr this instability has evolved in such a way to produce thick central region with poorly defined spiral arms.

\textcolor

black

\textcolor

black Figure 7 shows the initial \textcolorblackdisc thickness and the thickeness at $T=1.0$ Gyr. Comparing this to Fig. 3 it is clear that the Gas4 gas \textcolorblackdisc is initially about half the thickness of the GalactICS-Gas0 \textcolorblackdisc. The height in the central region increases to about the same as the Gas0 model in the final timesteps, but rather than staying at some constant thickness before flaring in the outer radii, the Gas4 thickness decreases to a minimum near 3 kpc before flaring to the outer radii.

\textcolor

black

\textcolor

black Figure 8 shows the evolution of the azimuthally averaged surface density profile. The changes to the averaged profile are similar in scale and origin to the changes that occur for the GalactICS-Gas0 model.

\textcolor

black It is important to note that both the GalactICS-Ga0 and Gas4 models are fairly different from real galaxies. Nonetheless, these two models demonstrate the ability of GalactICS to generate both stable and unstable models. These models are relatively straightforward, consisting of only a halo and a thin (or very thin) gas disc.

\textcolor

black

4 A Milky Way Model

In this section, we demonstrate the utility of our new GalactICS code by generating and evolving a realistic five-component model for a MW-like galaxy. Our model is based on the mass models from McMillan (2017). In that paper observations of Maser line-of-sight velocities, the proper motion of Sgr A*, the terminal velocity curve, the local vertical force, and the total mass within 50 kpc were used to constrain the parameters of an analytic model for the mass and gravitational potential of the Milky Way. The McMillan (2017) model comprised an NFW halo, thin and thick stellar \textcolorblackdiscs, a bulge based on Bissantz & Gerhard (2002), and both a HI and \textcolorblackH2 gas \textcolorblackdiscs.

In this paper, we build a GalactICS model that uses the "most likely" parameters from McMillan (2017). We assume an NFW halo with $\sigma_{h}=298~{}\mathrm{km~{}s}^{-1}$ and $r_{h}=19.6~{}\textrm{kpc}$ ($\rho_{0}=8.55\times 10^{6}\,M_{\odot}\,{\rm kpc}^{-3}$). For the thin \textcolorblackdisc we assume $M_{d,{\rm thin}}=3.52\times 10^{10}~{}M_{\odot}$, $R_{d,{\rm thin}}=2.5~{}\textrm{kpc}$ and $z_{d,{\rm thin}}=0.28~{}\textrm{kpc}$ while for the thick \textcolorblackdisc we assume $M_{d,{\rm thick}}=7.2\times 10^{9}~{}M_{\odot}$, $R_{d,{\rm thick}}=3.02~{}\textrm{kpc}$, $z_{d,{\rm thick}}=0.83~{}\textrm{kpc}$. Note that the scale height parameters are somewhat smaller than those used in McMillan (2017) to account the fact that the vertical structure in those models is exponential in $|z|$ whereas GalactICS \textcolorblackdiscs have an approximately ${\rm sech(z)}^{2}$ form. The latter form naturally arises from the assumption that the \textcolorblackdiscs are vertically isothermal.

By design GalactICS bulges have a surface density profile that is given, to a good approximation by the Sérsic profile. To convert the Bissantz & Gerhard (2002) bulge from McMillan (2017) to a Sérsic one we performed a simple parameter search using the inner 2 kpc of the density profile. This search yielded $n\simeq 2$ for the Sérsic index, $R_{b}\simeq 0.64$ kpc for the radial scale length, and and $\sigma_{b}=304~{}\mathrm{km~{}s}^{-1}$ for the velocity scale (see Widrow et al. (2008)).

Finally, we consider the gas components from McMillan (2017). Since our current version of GalactICS has just a single exponential gas \textcolorblackdisc we fit this model to the total gas surface density from the combined HI and \textcolorblackH2 gas \textcolorblackdiscs of McMillan (2017). Fig. 9 shows the comparison of the GalactICS \textcolorblackdisc to the two \textcolorblackdiscs in McMillan (2017). Our \textcolorblackdisc has $M_{g}=2.39\times 10^{10}~{}M_{\odot}$ and $R_{g}=13.1~{}\textrm{kpc}$. Note that the McMillan (2017) \textcolorblackdiscs have a hole in the centre so the structure of the gas components in the two models is rather different though the total mass in gas is very similar. Our model is constructed with a gas temperature of $10^{4}~{}K$.

As with the models in the previous section, the system is evolved using Gadget-2. Our N-body model comprises $10^{6}$ gas particles, $5\times 10^{5}$ bulge particles, $2\times 10^{6}$ thin \textcolorblackdisc particles, $10^{6}$ thick \textcolorblackdisc particles, and $5\times 10^{6}$ halo particles. The system is evolved for a total of 2 Gyr.

Fig. 10 shows circular speed decomposition of our model at $T=0$ and $T=1$ Gyr. We see that the \textcolorblackdisc components (primarily the thin \textcolorblackdisc) dominate the radial force at radii of about two \textcolorblackdisc scale lengths. Thus, we expect the model to form a bar and spiral structure, which indeed it does. The surface density profiles of the different components are shown in Fig. 11. The thin \textcolorblackdisc dominates the baryon mass budget at radii $R<13\,{\rm kpc}$ except in the innermost region where the bulge makes a comparable contribution to the surface density. On the other hand, the gas \textcolorblackdisc, which in this model has a larger radial scale length, dominates the outer \textcolorblackdisc.

The evolution of the circular speed and azimuthally averaged surface density profiles is due to the formation of a bar and spiral arms. From the surface density profiles we see that mass in the \textcolorblackdisc components moves outward from $R\sim 3-5$ kpc to $R\sim 9-11$ kpc. This mass redistribution leads to a change in the shape of the rotation curve.

Figure 12 shows a snapshot of the model at $T=1$ Gyr. As noted above, the \textcolorblackdisc components have all developed bar and spiral structure, although the bar strength varies considerably between the various \textcolorblackdisc components, \textcolorblack. The difference between the components can be seen more clearly in Fig. 13, which shows the evolution of $A_{2}$ for the three \textcolorblackdisc components as a function of radius and time where $A_{2}$ is the amplitude of the $m=2$ Fourier component for the surface density where $m$ is the usual azimuthal mode number. At all times and radii, the thin and thick \textcolorblackdiscs show similar structure, although the $A_{2}$ moment is always weaker in the thick \textcolorblackdisc. Both stellar \textcolorblackdiscs show a clear bar component with peaks in the $A_{2}$ moment at $R\approx 2.5$ kpc for \textcolorblack$T\geq 1.0$ Gyr. The gas \textcolorblackdisc has a much more transitory structure. While there is some structure in the bar region, the dominant $m=2$ structures are at larger radii and clearly correspond to the two-armed spiral structure seen in Fig. 12.

\textcolor

black The difference in the bar strengths between the thin and thick discs is similar what is observed in both nature and simulations. Athanassoula (1983) ran simulations of galaxies with single discs and noted that dynamically colder discs formed thinner bars than hotter discs. Athanassoula (2003) followed this investigation with a theoretical and numerical exploration of the relationship between angular momentum transfer and bar properties. In that work it was found that increased disc velocity dispersions inside the co-rotation radius decreased the bar strength. More recently Athanassoula et al. (2016) ran simulations of merging galaxies. They found that the thin disc also formed a stronger bar component than the thick disc. The bar present in their thick discs have the same length and orientation as in the thin disc bar, but have larger vertical extents and are significantly more oval in nature. Athanassoula (2018) found similar results using a more advanced analysis of these merging simulations clearly show the difference between the thin and thick discs. The simulations of Bekki & Tsujimoto (2011) examined the formation of the bulge from a single or two-disc scenario. They found that a thick disc is less influenced by the presence of a bar, allowing for a steep vertical metallicity gradient. Similarly, Debattista et al. (2017) utilized variety of different simulations to show that populations with lower radial velocity dispersions form stronger bars than those that are dynamically hotter. Fragkoudi et al. (2017) found similar results in their simulations, noting that the thin disc bar was $\sim$50$\%$ stronger than the bars found in the thick disc, as well as being significantly more elongated.

The bar and spiral structure seen in our evolved model are similar to what are observed in the Milky Way. For example, the bar in the MW has a length of about $3-4$ kpc (see Bissantz & Gerhard (2002) and references therein), which is comparable to what we find in our simulation. Likewise, the spiral structure shows a strong two-armed pattern with evidence for weaker four or more armed spiral structure. Our structures do seem to have a tighter winding pattern than what is seen in the Galaxy.

The origin of the bar and spiral structure pattern of the Milky Way is still an open question. For example, Purcell et al. (2011) argued that these structures might be the result of an encounter between the \textcolorblackdisc and the Sagittarius dwarf. Our results suggest that the \textcolorblackdisc can form structures of this type through internal instabilities. \textcolorblack However, this simulation has only been run for a few orbits. Further investigation must be done to examine whether the bar and spiral structure seen here will persist over dozens of dynamical times.

5 Summary and Conclusions

In this paper, we have introduced a new version of the GalactICS code that is able to generate equilibrium models for a galaxy consisting of thin and thick stellar \textcolorblackdiscs, a stellar bulge, a dark halo, and gas \textcolorblackdisc. The models can be used to generate ICs for N-body simulations that, in turn, can be used to study the dynamics of real galaxies and, in particular, the formation of bars and spiral structure.

\textcolor

blackTest bed simulations of models with just a gas and dark halo showed the ability of GalactICS to produce both stable and unstable models. A transient is present in some gas-halo only models. The inclusion of a stellar \textcolorblackdisc suppresses these transients, yielding stable models. Regardless of the presence of a transient, the azimuthally averaged surface densities remain roughly constant. The tangential velocities of the gas particles have an increase in random motions, but the gravitationally calculated circular speed remains constant.

We demonstrated the applicability of the model by constructing a multi-component dynamical model for a MW-like galaxy based on the mass model of McMillan (2017). We found that this model produces a relatively strong bar and spiral structure by \textcolorblack$T=1.0$ Gyr. The bar is present in all of the \textcolorblackdisc-like baryonic components though strongest in the thin \textcolorblackdisc. Spiral arms are also present in all of the \textcolorblackdisc components though noticeably thinner and more prominent in the gas \textcolorblackdisc.

GalactICS provides an excellent tool for dynamical studies of galaxies. In particular, it generates N-body ICs that that can be used to study the dynamics of observationally-motivated galaxy models. The inclusion of a gas \textcolorblackdisc opens up new and exciting possibilities that we are eager to explore in future work. Indeed, the models have already been used to study the rotation curves of barred spiral galaxies (Randriamampandry et al., 2018). Many standard algorithms for galaxy modelling fail for particular bar orientations. A suite of tailored ICs generated with GalactICS allowed us to model galaxies of this type. The simulations shown in the present work suggest that dynamical simulations combined with observations might allow one to constrain the models in a fashion not possible with mass models of the type considered in McMillan (2017).

CC’s work is based upon research supported by the South African Research Chairs Initiative (SARChI) of the Department of Science and Technology (DST), the SKA SA and the National Research Foundation (NRF). ND’s work is supported by a SARChI’s South African SKA Fellowship. LMW is supported by the Natural Sciences and Engineering Research Council of Canada through Discovery Grants. The numerical simulations were performed at the Centre for High Performance Computing.