Numerical code for multi-component galaxies: from N-body to chemistry and magnetic fields

We present a numerical code for multi-component simulation of the galactic evolution. Our code includes the following parts: $N$-body is used to evolve dark matter, stellar dynamics and dust grains, gas dynamics is based on TVD-MUSCL scheme with the extra modules for thermal processes, star formation, magnetic fields, chemical kinetics and multi-species advection. We describe our code in brief, but we give more details for the magneto-gas dynamics. We present several tests for our code and show that our code have passed the tests with a reasonable accuracy. Our code is parallelized using the MPI library. We apply our code to study the large scale dynamics of galactic discs.


Introduction
An explosive growth of observational data on dynamics and chemical composition of various components (subsystems) of galaxies is seen in last decade.To interpret the huge data cubes we need reasonable models included more physical processes like magnetic fields, star formation, feedback, chemical kinetics, dust grain dynamics and so on.Spatial scales of the processes in a galaxy vary from sub-parsecs for turbulent structures in the warm interstellar medium to several kiloparsec size for galactic disc and winds.Formation of molecules and distribution of dust grains strongly affect on the efficiency of star formation, which, in turn, changes the thermal state of a gas through stellar winds and supernova explosions [1].Then, to study the interaction between various components numerically someone needs a code with a lot of physical modules.
In recent years several open source codes were developed for astrophysical applications, such as ZEUS 1 , RAMSES2 , GADGET3 , Athena4 , FLASH 5 and many others.These codes are generalpurpose ones, so that each of them has own advantages and disadvantages in the application to the galactic evolution [2,3,4,5].However galaxies are multi-component systems, and its simulation requires a robust model included special physical processes.So the main goal of this project is to develop the numerical code for simulation of the evolution of disc galaxies including generation of spiral structure, physics of interstellar medium, formation of clouds and stars [6,7].
This paper is organized as follows.In Section 2 the main numerical methods are described.Section 3 presents the results of tests.Section 4 describes the application of our code to the evolution of stellar-gaseous galactic disc.In Section 5 we summarize our results.

Description of the code 2.1. Gas dynamics
The gas dynamics equations for an ideal magnetized gas can be written in the conserved variables U in Cartesian coordinates U = [ρ, ρv x , ρv y , ρv z , E, B x , B y , B z ] T . (1) Then, the conservation laws can be written in a compact form where F, G, and H are vectors of fluxes in x, y, and z-directions, respectively and S is the source vector.The components of these vectors are where P * = P + B 2 2 and total gas energy . For pure gas dynamics we set B ≡ 0.
In multi-dimensional magneto-gas dynamics to obtain zero-value for the magnetic field divergence we apply the Constrained Transport technique for magnetic field transport through computational domain [8,9].In this approach magnetic field strength is defined at faces of a cell, while other gas dynamical variables are defined at the center of a cell.
We solve the equations (2) in Cartesian coordinates (x, y, z).The center of each cell is located at the position (x i , y j , z k ).For each cell the faces being normal to the x direction have coordinates x i±1/2 and sizes δx, δy and δz.The gas dynamic variables (density, momentum, energy) are volume-averaged and defined at the center of a cell.For example, the density value can be written as follow While the magnetic field components are surface-averaged and found at faces of a cell: To calculate the momentum and energy fluxes we apply the second-order cell-centered averages for magnetic field, e.g. for B x : We use a finite-volume discretization for gas dynamics equations and finite-area discretization for magnetic field transport equations.In this case gas-dynamical variables at time t + δt can be found as follow where the vectors i,j,k+1/2 are the time-and area-averaged fluxes through the x, y and z-faces of a cell, correspondingly.
Using the Stokes' law for the magnetic field transport equation we obtain time-evolution of magnetic field equation: where E x , E y and E z are the components of electric field (or electro-magnetic force): The approximation of the electro-magnetic force can be written as follow: where the derivative of E x for each grid cell face is computed by selecting the upwind direction according to the contact mode, namely, Note that in 3D the expressions similar to the above are required to convert the x-and ycomponents of the electric field to the appropriate cell corners.These expressions might be obtained directly cyclic permutation of the (x, y, z) and (i, j, k).
Godunov-type methods are widely used for solving of hyperbolic equations.An idea is to use the Riemann solution for the decay of an arbitrary discontinuity.In this case it is assumed that the solution might be non-smooth and discontinuities might be in every computational cell.That is why, this kind of numerical schemes is a universal tool for simulation of gasdynamical problems of evolution shocks and contact discontinuities.The main problem of Godynov-type methods is a necessity of solving system of nonlinear equations.Usually this leads to iterative procedure, which requires additional computational resources.To resolve this issue there are a lot of approximate solvers of Riemann problem.Our code includes two types of solvers for pure gas dynamics: Harten-Lax-van Leer (HLL) or Harten-Lax-van Leer-Contact (HLLC), and one for magneto-gas dynamics, namely, Harten-Lax-van Leer-Discontinuities (HLLD).These methods are described by [10] in detail.
An exact or approximate Riemann solver requires values for the conserved variables at left and right interfaces of each cell: These values can be reconstructed from cell centered values with second-order or third-order accuracy, on a choice.The interpolation of the values is carried out for primitive variables w ≡ {ρ, P, v, B} as follow and then and minmod is a limiter function [11].The parameter κ can be equal to −1, 0, or 1/3, which represent second-order fully upwind, second-order upwind biased, and third-order upwind biased cases, respectively.Below we use κ = 1/3.

N -body dynamics
We use the "particle-in-cell" method for the dynamics of both dark matter halo ("live" halo) and stellar discs.To obtain the second-order accuracy on time we apply the "flip-flop" integrator.The main advantage of this approach is that it requires only one calculation of forces per time step for each particle.More details about our N -body solver for galaxy dynamics can be found in [12].

Poisson solver
To take into account self-gravity of a gas we solve the Poisson equation with the gas dynamics conservation laws (2): ∆Ψ(x, y, z) = 4πGρ(x, y, z) .
In our code the Poisson equation can be solved using three different methods: FFT-based, TreeCode [13] and over-relaxation Gauss-Seidel method.To calculate self-gravity in Nbody/gasdynamical problems we use regular Cartesian grid.For interpolation of particle density to a mesh we use the second order finite-volume interpolation, which conserves the total mass of the system with a good accuracy.

Test problems
Below we present several tests for our gasdynamical code in 1D, 2D and 3D.

The Sod shock tube test
This is the most simple test demonstrated the formation of shock wave, contact discontinuity and rarefaction wave.The initial distribution of gasdynamical values is taken as usual: ρ, v x , P = {1, 0, 1} at x < 0; ρ, v x , P = {0.125,0, 0.1} at x > 0. Initially the velocity is equal to zero. Figure 1 presents the result of numerical simulation for different resolution and the exact solution.One can see a good coincidence of the numerical solution with the exact one even for low resolution.
Figure 2 shows the initial density distribution (black line).We study the dependence on spatial resolution: the case for N = 200 is exactly smooth, whereas oscillations appear for N = 500, for further increase of resolution (N = 1500) our result becomes closer to the exact solution [14].

The bow shock simulation
A bow shock usually forms due to supersonic motion of star (or planet) through interstellar medium [15] or a galaxy through the intercluster medium [16], so that this is one of the most important astrophysical problem.To test our code we simulate a supersonic gas flow through potential well.The initial distribution of the parameters are homogeneous: ρ = 1.0,P = 1.0, v x = 1.0, and γ equals 1.4.The external potential well is set in the form Ψ(x) = −Ψ 0 exp(−(x/x 0 ) 2 ), where Ψ 0 = 5 and x 0 = 0.1.The boundary conditions are free at the right and the steady inflow at the left.
At t = 0.051 the shock wave is formed on the far edge of the well (towards to the flow) due to supersonic falling of the gas to the potential well (red line).Note that this configuration is unstable.So the shock wave moves upstream through the potential well and after crossing the minimum of the potential it becomes steady-state at front edge of the potential well (blue line).

The Brio-Wu problem
The Brio-Wu test is one dimensional magneto-gas dynamics problem.The solution of this test consists of the fast rarefaction wave, that moves to the left, the intermediate shock wave, the slow rarefaction wave, the contact discontinuity, slow shock wave and another fast rarefaction wave, that moves to the right (see figure 4).The intermediate shock and slow rarefaction waves  This test (and the next one, the Ryu-Jones problem) is compared with the results obtained by the Athena code [18].For both simulations we set the same spatial resolution N = 200 (figure 4).One can see a good agreement between the results, but the artificial oscillations can be found for the longitudinal velocity component for both solutions.

The Ryu-Jones problem
The Ryu-Jones problem tests the rotation of the magnetic field components.The solution consists of two fast shock waves with velocities 1.22 and 1.28 Mach numbers, two slow shock  waves with 1.09 and 1.07 Mach numbers, two rotational and one contact discontinuities [19].
The initial distribution is: ρ = 1.08,P = 0.95, v x = 1.2, v y = 0.01, v z = 0.5, One can see a good agreement between results obtained both in our code and in the Athena code.However the peaks of v y and B y in our simulation are smoother than these obtained in the Athena code [18].

The rotor problem
This test demonstrates fast rotation of the cylinder in the non-moving medium with homogeneous magnetic field.Initially there is a disc in the center of computational domain: ρ = 10, v x = −v 0 y/r 0 and v y = v 0 x/r 0 at r < r 0 = 0.1; ρ = 1 + 9f (r), v x = −v 0 f (r)y/r 0 , v y = v 0 f (r)x/r 0 at r 0 < r < r 1 = 0.2, and ρ = 1 at r > r 1 ; P = 1 and B x = 5/ √ 4π are constant in the whole computational domain.Here the function f (r) = (r 1 − r)/(r 1 − r 0 ) is the linear interpolation of the gas density and velocity.This configuration is strongly unstable, because the centrifugal forces are not balanced.A rotating gas is redistributed gradually within the computational domain capturing the stationary gas.Under these conditions the magnetic Figure 7.The Sedov-Taylor blast-wave.The 2D density map for pure gas dynamics (left panel) and the radial pressure profile at t = 0.1 for our simulations (red dots in left middle panel) compared with the exact solution (blue line).The 2D density and magnetic pressure maps at t = 0.1 for magneto-gas dynamics (two right panels).field keeps the rotating material in the flattened form (Figure 6).

The Sedov-Taylor blast-wave
One of the most well-known self-similar solution describes a strong shock wave originated from the point explosion.For zeroed magnetic field a numerical solution can be compared with the exact solution.In 2D we simulate a strong explosion with the following parameters: ρ = 1, P = 10 −5 are set in the whole computational domain, P = π 2 r 2 0 /(γ − 1) at r < r 0 .We set r = 0.01 and take a computational domain [−0.5, 0.5] × [−0.75, 0.75] and number of cells 512 × 768. Figure 7 shows the 2D density map (left panel) and the radial pressure profile for our simulations (red dots in left middle panel) compared with the exact solution (blue line).One can see a good agreement between numerical data and analytic curve.
We consider a blast wave in a magnetized medium.At t = 0 we set B x = B y = 1/ √ 8π.Despite of the symmetric initial distribution the shock front expands larger along the lines of the magnetic field.

The Kelvin-Helmholtz instability
Shear instabilities are very common in astrophysical objects, so that numerical algorithms should reproduce such inabilities well.Here we simulate the evolution of two gas layers moving with different velocities and separated initially by a contact discontinuity: ρ = 1, v x = −0.5 at |y| < 0.25 and ρ = 2, v x = 0.5 otherwise, and P = 2.5 everywhere.We add a random perturbation of the velocity with the amplitude 0.001.The grid size is 512 × 512.
Figure 8 shows the evolution of gas density for pure gas dynamics (upper row of panels) and magneto-gas dynamics (bottom row of panels).One can note the turbulent flow formation due to shear instability at early times, which produces large scale vortexes at further nonlinear stage.The presence of magnetic field strongly changes the picture: a small scale perturbations disappear for initial longitudinal magnetic field B x = 0.5 (see bottom line in Figure 8), whereas a large scale oscillation of a gas along magnetic field can be seen.

Spherical collapse of gravitating gas
To test a gravity solver we use a standard 3D cosmological problem -the collapse of a gaseous sphere with initial density profile ρ(r) = 1/r.We set the cubic computational domain (x, y, z) ∈ [−1, 1] with cell number 100 3 .Figure 9 shows the evolution of the density profile.

Thermodynamical module
To study the evolution of the interstellar medium in galaxies we develop a module of thermal processes.In this module we can switch between tabulated cooling/heating rates and calculation of rates for main chemical species in the interstellar medium.For the former we can use tables of cooling rates in the temperature range 10 K< T < 10 8 K calculated by [20,21,22].In the latter we calculate the cooling and heating rates produced by specific emission processes and in the temperature range T < 2 × 10 4 K, that is of a great importance for studying thermal state of the interstellar medium (see [23] for instance).We can investigate the thermal evolution of the interstellar medium for different abundances of heavy elements (metallicities).Our thermal block includes following cooling processes: cooling due to recombination and collisional excitation and free-free emission of hydrogen [24], molecular hydrogen cooling [25], cooling in the fine structure and metastable transitions of carbon, oxygen and silicon [26], energy transfer in collisions with the dust particles [27] and recombination cooling on the dust [28].The heating rate takes into account photoelectric heating on the dust particles [28,27], heating due to H 2 formation on dust and the H 2 photodissociation [29] and the ionization heating by cosmic rays [30].For multi-level atoms the cooling rates are obtained from the level population equation assuming the optically thin steady-state regime (see e.g.[26]).
Figure 10 presents cooling rates in the low temperature (T < 2 × 10 4 K) range adopted for our model of the interstellar medium.H , for solar metallicity.The total cooling rate is depicted by black line.

H 2 kinetics
Molecular hydrogen (H 2 ) regulates star formation in galaxies [23], then to study the evolution of the interstellar medium we need to take into consideration the H 2 kinetics.Because of significant difference between dynamic and chemical timescales the H 2 kinetics is nonequilibrium and to get a source term in equation ( 2) we solve a system of ordinary differential equations for the following chemical species: H, H + , H 2 .Molecular hydrogen is formed on the surface of dust grains and dissociated by ultraviolet Lyman-Werner photons and cosmic rays.Because of strong absorption and scattering of ultraviolet photons in the neutral hydrogen and dust to mimic radiative transfer and H 2 self shielding we use a simplified approach introduced by [31] for calculation of neutral and molecular hydrogen column densities.

Star formation implementation
Stellar feedback effects can significantly influence on the galactic evolution.However, there is no consensus about the implementation of the star formation into gas dynamics, that is clearly demonstrated in [32].Their results of the galaxy evolution simulations within the ΛCDM paradigm strongly depend on the star formation and feedback recipes: stellar mass, size, morphology and gas content of the galaxy at z = 0 vary significantly due to the different implementations of star formation and feedback.Despite this problem we include the star formation effects using an intuitive approach, that partly based on the well-known method offered in [33].
To form a star or in general stellar particle we need to find cells in a computational domain, which satisfy to our conditions for star formation.These conditions can be more or less sophisticated and sometimes may have intricate nature.Usually the following criteria for the stellar particle creation in a given cell are considered: the surface density should be greater than the threshold for star formation Σ t and simultaneously the temperature should be less than some critical level T t .In the simulations presented below we have assumed Σ t > 10 4 cm −3 and T t ≤ 10 2 K.If these conditions are realised in a given cell, then we create a test stellar particle.The Salpeter initial mass function is assumed for each particle, which usually consists of 10 3 − 10 4 stars.The initial velocity of this test particle is taken to be the same as its host cell.The kinematics of test particles are followed by the second-order integrator mentioned in Section 2.2.We take into account several feedback effects: stellar winds from massive stars, radiative pressure, supernova explosions, and stellar mass loss by low-mass stars [34].Metals ejected by supernovae and lost by low-mass stellar population can strongly change chemical composition and thermal evolution of a gas.Because of local character of enrichment process the re-distribution (mixing) of metals is of great importance for further star formation in galaxies [35,36,37,38].

Simulation of the multicomponent galactic disc
Here we present an example of the gravitationally unstable stellar-gaseous galactic disc evolution.In this simulation we take into account magnetic field and radiation processes.The stellar disc is simulated using N -body method with number of particles N = 10 7 .The number of cells for gas dynamics equals to 1024 × 1024 × 128.We start our simulation from the equilibrium state of the galactic gaseous disc within the fixed potential of the dark matter halo for Toomre's parameter Q T = 1.2.A construction of initial equilibrium configuration for stellar-gaseous disc is described by [39] in detail.We assume that the galactic disc consists of stellar component and gas with temperature 10 4 K in the middle plane.The magnetic field is set as the toroidal configuration with B = 6 µG.
Figure 11 presents the result of this simulation.The distribution of stellar particles (top row of three left panels) and gas surface density (bottom row of three left panels) are shown at time moments t = 0; 300; 600 Myr.One can see how the gravitational instability drives the formation of spiral pattern, which is clearly seen in both stellar and gaseous components.At t = 300 Myr it is easily found narrow and dense galactic shocks (bottom row of panels).Because of relatively small initial value Q T a flocculent structure is rapidly developed in the spiral arms [40].The spatial structure of a gas becomes very complex due to numerous shear flows, thermal and gravitational instabilities, magnetic field pressure influence as well.
One can see the formation of the large number of dense gaseous complexes with number density ∼ 300 − 1000 cm −3 .These complexes are mostly located in galactic spiral arms.This result is almost independent on numerical resolution and slightly dependent on the initial magnetic field.The increase of numerical resolution leads to formation of clouds on smaller scales.Such clouds are expected to give birth of stars.
Figure 11(upper right panel) shows the structure of magnetic field at galactic disc at t = 600 Myr.The amplitude of magnetic field decreases with radius from |B| = 8µG at 4 kpc down to |B| = 4µG at the outer part of the disc.This is in a good agreement as with observations [41] and with cosmological simulations [42].The vector map (right bottom panel in Figure 11) represents the chaotic and regular velocity components.The regular field corresponds to the rotation of a gas and large-scale spiral structure in the disc.

Conclusion
In this paper we have described our three-dimensional numerical code for multi-component simulation of the galactic evolution.This code has been mainly developed to study the evolution of disc galaxies taking into account generation of spiral structure, physics of interstellar medium, formation of clouds and stars.Our code includes the following ingredients: N -body dynamics, ideal magneto-gas dynamics, self gravity of gaseous and stellar components, cooling and heating processes, star formation, chemical kinetics and multi-species gaseous and particle (for dust grains) advection.We present several tests for our code and show that our code have passed the tests with a resonable accuracy.It should be noted that our code is parallelized using the MPI library.We apply our code to study the large scale dynamics of galactic discs in context of formation and evolution of galactic spiral structure [39], molecular clouds formation [43] and disc-to-halo interactions [44].

Figure 1 .
Figure 1.The Sod shock tube.The exact solution is shown by black line, the other lines represent the distribution at t = 0.5 with grid resolution N = 200 (green), N = 500 (red), N = 1500 (blue).One can find five regions in the solution:the regions 1 and 5 correspond to the initial unperturbed state of a gas, the rarefaction wave is found in the region 2, the regions 3 and 4 are separated by contact discontinuity, and the shock wave is located between regions 4 and 5.

Figure 2 .
Figure 2. The Shu-Osher test.The initial density distribution is shown by black line, the other lines represent the distribution at t = 0.18 with grid resolution N = 200 (green), N = 500 (red), N = 1500 (blue).

Figure 4 .
Figure 4.The Brio-Wu problem.Our result is depicted by red dots, the result obtained in the Athena code is shown by blue line.

Figure 5 .
Figure 5.The Ryu-Jones problem.Our result is depicted by red dots, the result obtained in the Athena code is shown by blue line.

Figure 6 .
Figure 6.The rotor problem.The maps of density, pressure, magnetic pressure and Mach number are shown from left to right correspondingly.

Figure 8 .
Figure 8.The Kelvin-Helmholtz instability.The gas density distributions at t = 0, 0.1, 4, 12 (from left to right) in the pure gas dynamics (top row of panels) and the same, but in the presence of magnetic field (bottom row of panels).

Figure 9 .
Figure 9.The density profile for t = 0 and t = 1 is depicted by cyan and black lines, respectively, the pressure profile for t = 1 is shown by red line (left panel).The velocity profile for t = 1 is shown at right panel.

Figure 10 .
Figure 10.Cooling functions, Λ/n 2 H , for solar metallicity.The total cooling rate is depicted by black line.

Figure 11 .
Figure 11.The Galactic disc evolution.First three columns: stellar (top) and gaseous (bottom) distributions at t = 0, 300, 600 Myr.Right pictures: radial variations of the module of magnetic field (top) and vectors of magnetic field in the disc plane.