EUNHA: A NEW COSMOLOGICAL HYDRODYNAMIC SIMULATION CODE

Abstract

We develop a parallel cosmological hydrodynamic simulation code designed for the study of formation and evolution of cosmological structures. The gravitational force is calculated using the TreePM method and the hydrodynamics is implemented based on the smoothed particle hydrodynamics. The initial displacement and velocity of simulation particles are calculated according to second-order Lagrangian perturbation theory using the power spectra of dark matter and baryonic matter. The initial background temperature is given by Recfast and the temperature uctuations at the initial particle position are assigned according to the adiabatic model. We use a time-limiter scheme over the individual time steps to capture shock-fronts and to ease the time-step tension between the shock and preshock particles. We also include the astrophysical gas processes of radiative heating/cooling, star formation, metal enrichment, and supernova feedback. We test the code in several standard cases such as one-dimensional Riemann problems, Kelvin-Helmholtz, and Sedov blast wave instability. Star formation on the galactic disk is investigated to check whether the Schmidt-Kennicutt relation is properly recovered. We also study global star formation history at different simulation resolutions and compare them with observations.

1. INTRODUCTION

Over the last few decades, the cosmological gravitational N-body simulation (Kim et al. 2009, 2011; Rasera et al. 2013; Angulo et al. 2012; Kuhlen et al. 2012) has proven to be a powerful tool for study of the Large-scale Structures (LSS) of the Universe. The observed distribution of galaxies is well reproduced with the simulations in terms of the two-point correlations (Masaki et al. 2013) or power spectrum of galaxies (Tegmark et al. 2006), cosmic topology (Choi et al. 2010, 2013), and the abundance of groups (Nurmi et al. 2013) or the largest structure in the universe (Park et al. 2012). These statistical studies enable us to refine the current cosmological model by constraining model parameters at the percent level of accuracy.

Below the galaxy scale, however, hydrodynamic force begins to overtake gravity. In the potential well of a dark matter halo, baryonic matter decouples from the surrounding dark matter inflow and forms distinct and complex inner structures such as galactic disks, bulges, spiral arms, and so on. Therefore, it would be indispensable to include the hydrodynamical effects if one wants to study galaxy formation and evolution by numerical simulations.

The Two Degree Field Galaxy Redshift Survey (2dF-GRS; Colless et al. 2001) and a series of the Sloan Digital Sky Survey (SDSS; Aihara et al. 2011), for example, have enabled us to study physical properties of and the environmental effects on individual galaxies. In the current paradigm of hierarchical clustering, it is believed that a galaxy is a result of a series of mergers of smaller objects. Therefore, the secular evolution of a galaxy and interaction between galaxies began to be investigated in great detail in the cosmological context. Park & Hwang (2009) reported that the morphology of a galaxy is correlated with that of the nearest neighbors. This correlation maybe reasonably attributed to hydrodynamic interactions within the galaxy and its neighbors.

An observed galaxy is an outcome of various physical processes accumulated over a long period of time. Galaxies have been shaped not only by gravity and hydrodynamic force but also by other astrophysical processes like cooling, heating, star formation, and supernova feedback. Baryonic matter is believed to settle down in the deep potential wells of dark matter and stars form therein, in high density environments. After consuming all available fuel by nuclear fusion, massive stars at last explode as supernovae dispersing energy and metal-enriched material into the interstellar medium. This explosion may trigger another star formation, and so the baryonic component starts another life cycle of metal enrichment. Consequently, the new generation of stars has higher metallicity than before. Therefore, a galaxy may have a wide spectrum of stellar populations and metallicities.

There are many cosmological hydrodynamic simulation codes available today (Kravtsov 1999; Fryxell et al. 2000; Teyssier 2002; O’Shea et al. 2004; Wadsley, Stadel, & Quinn 2004; Springel 2005; Wetzstein et al. 2009; Springel 2010, among others). Largely, they can be grouped in two categories; Eulerian and Lagrangian codes. The Eulerian scheme use regular grids to compute the hydrodynamic interaction between grid points. However, the Lagrangian code is based on particles bearing hydrodynamic properties. Mutual interactions between particles are computed using the SPH (smoothed particle hydrodynamics).

The GOTPM (Dubinski et al. 2004) code is well designed for massive N-body simulations with efficient use of memory space and fast calculation speed. For example, the Horizon-Runs (HRs; Kim et al. 2009, 2011), which were the largest ones at that time, were performed with 70 – 350 billion particles. The GOTPM code adopts non-recursive walks on oct-sibling tree structures which speeds up the code significantly.

On top of the GOTPM code, we placed the SPH algorithms to properly simulate cosmological structures down to the galaxy scale the hydrodynamic force is more dominant than gravity. Among the SPH algorithms, we adopt the entropy-conservation scheme (Springel 2005). We used the CLOUDY 90 package (Ferland et al. 1998) to arrange a table of the heating and cooling rates for reference during the simulation run. We set a global value for the reionization epoch and the heating rated before and after reionization are different. After reionization, Jeans’s mass may increase and small mass objects are prevented from forming due to baryonic pressure resisting gravitational contraction. We include the self-shielding of gas particles from uniform UV background, which may have a larger effect on the formation of the dwarf galaxies in the early universe (Tajiri & Umemura 1998; Sawala et al. 2010).

This paper is organized as follows: In Section 2, we show the basic equation of motion of N-body particles. The hydrodynamic equations are listed in Section 3. Section 4 is devoted to a description of the individual time stepping adopted in the code and Section 5 shows how to implement astrophysical processes in the code. Simulation tests and the results are discussed in Section 6. We summarize our work in Section 7.

 

2. BASIC EQUATIONS OF MOTION

In an expanding background, the comoving distance (x) is related to the physical distance (r) with the scale factor (a) as

where am is the scale factor at z = 0. In the GOTPM code, we use the following variables

where xs is the position of a particle in unit of the mean particle separation (ℓb ≡ Lb/N), Lb is the simulation box size on one side, N is the number of grids in one direction, mr is the physical mass of the particle, ms is the particle mass in simulation unit, and ⟨ρ0⟩ is the mean density at the current epoch (a = am). Then, the physical acceleration on a particle can be expressed as

Henceforth, the subscript r and s mean that the quantity is given in physical units and simulation units, respectively. In comoving space, Equation (4) can be rewritten as

The second term in the parenthesis of the right-hand side of the equation is the acceleration due to the homogeneous background and can be neglected if the universe is isotropic.

The acceleration, ar, can be decomposed into two parts, due to the gravitational and hydrodynamical force as below,

where as is defined as the acceleration in simulation unit, ), is the critical density at z = 0, and the summation index is running over all other particles. After some calculations, it can be shown that Equation (6) has the following form;

Here, the initial matter density, Ωi is given by

In an expanding medium, it is helpful to split the velocity into two parts; the Hubble flow and the peculiar velocity. Then, the Hubble flow between two points of separation, is obtained as

where g1 = H(z)ℓb=(1 + z). The peculiar velocity (vp) is obtained from the simulation velocity (vs) as

where g2 = ag1.

 

3. HYDRODYNAMICS

We adopt the same entropy-conservation scheme of the SPH as used in the Gadget code (Springel 2005). But we use a different method to identify the N-nearest neighbors. Rather than using the predict-correct method adopted by the Gadget code, we apply an improved method, a direct search with the Oct-Sibling Tree, which is fast and reliable in identifying neighbors. We call this cosmological hydrodynamic code the EUNHA (Evolution of the Universe simulated with N-body and Hydrodynamic Algorithms).

3.1. Basic Equations

The pressure (Pr) and specific internal energy (ur) can be measured from the temperature (T) and density (ρr) of ideal gas as

where mH is the mass of hydrogen atom, 𝜇 is the mean molecular weight, kB is the Boltzmann constant, and γ is the adiabatic index. In an adiabatic process, the pressure is related to the gas density as , where Ar is the entropy.

Now we adopt the following relations of conversion between the physical and simulation units,

where ⟨ρz⟩ = ⟨ρ0⟩ (1 + z)3 = Ωm(1 + z)3 and

It should be noted that the simulation entropy (As) may change with time even when Ar is unchanged with time.

3.2. Smoothed Particle Hydrodynamics

The basic equations of the smoothed particle hydrodynamics are summarized as follows;

where h is smoothing length defined as the distance to the Nn’th nearest neighbor, W is the smoothing kernel, and fi is defined as

In this study, we used the spline kernel discussed in Monaghan (1992).

We have adopted the typical artificial viscosity effect to the acceleration as

where 𝚷ij is the viscosity factor introduced to capture the shock front. We adopt the form proposed by Monaghan (1997)

where and α is viscosity coefficient. We adopt α = 1 in this study. The simulation sound speed, cs, is defined as

and this is identical to the physical sound speed, cr. Finally, we get the viscous force as

3.3. Neighbor Findings

We use the Oct-Sibling Tree (OST) to find the N-nearest neighbor gas particles. The OST has been ex- tensively used in the Tree force calculation and has proven to be very fast because of the non-recursive nature of the sibling connections (Dubinski et al. 2004; Kim et al. 2011) between Tree nodes and particles. We exploit the OST for identifying the N-nearest neighbors (or smoothing length) with a simple modification of the original Tree-gravity routine.

The advantage of the tree searching is that it identifies the N-nearest neighbors without any assumption on the initial trial value (Thacker et al. 2000; Springel 2005). Although the predict-correct method has been widely used to determine the smoothing length, there is a tension between successive iterations when the smoothing length changes significantly, especially in regions where the number of neighbors dramatically changes with a small change of the searching length.

 

4. INDIVIDUAL TIME STEP

4.1. Subtime step

In a particular time step (Δt), the change in the position of a particle is quite simply

where v is the velocity and a is the acceleration measured at the position of the particle. This raises the question of determining the time step Δt for which the above expression is a valid approximation. In many cases, it is sufficient to constrain the time step by |r′ − r| ≤ 𝜖, which means that the position change should be less than the force resolution (𝜖) times the step size of the simulation. If the velocity of a particle is larger than the acceleration, it is reasonable to use

The time-step size for a given gas particle is determined as (Springel 2005)

where C is the Courant number, hsml is the smoothing length, and Vsig is the maximum signal velocity of the particle. The signal velocity between particle i and j is defined as

where is the normalized mutual displacement vector and vij is the relative velocity. Here, Vsig is the maximum value among ’s. This ensures that for a given time the “signal distance” should be less than the smoothing length scale multiplied by the Courant factor, which is fixed to 0.15 in this paper. It is important to note that Equation (27) includes not only the sound speed but also the relative radial speed between particles.

During the simulation run, particles experience various forces and their velocities change with time. For example, a supernova explosion may expel nearby gas particles in the radial direction and, therefore, the timestep size should be reduced to properly capture the supernova shock. In the original GOTPM code, however, global time stepping was adopted. In underdense regions, particle positions vary slowly and hence evolution over larger periods of times maybe computed in fewer iterations. Thus, adopting global time steps maybe wasteful and lead to bottlenecks in improving simulation performance. Hence, in the EUNHA code we implemented the individual time steps for both the N-body and hydro parts.

4.2. Merging N-body and Hydrodynamic Subtime Steps

For the EUNHA code, we adopt the global time step blocking time stepping, whose subtime step is obtained by dividing the global step size by the integer power of two;

where dt and Δt are the individual and global time step sizes, respectively, and p is an integer. This equation can be changed to p = ceil(log2(Δt=dt)) where ceil() is a round up function. Hereafter, we call p the subtime step power. As we are using two kinds of forces (gravitation and hydrodynamic forces), there are two individual step sizes. Dark matter particles have only N-body subtime step (Equation (26)) while gas particles have an additional hydrodynamic subtime step (Equation (27)). If a gas particle has different step sizes, we adopt the smaller value.

However, the individual time step may not properly capture the shock front. It is because that a particle in a pre-shock region (low T and Cs) may possibly have a significantly larger step size than the shock-passing time scale. In some extreme cases, the pre-shock particle may not even experience the shock front due to its large size of time step. In order to avoid this situation and to relax the tension between neighboring particles which have a large difference in the step size, we adopt the time step limiter (Saitoh & Makino 2009), which propagates the step size into neighbor particles so that the difference of the subtime step power among neighbor particles should not be larger than two. For a full description of the method, see Saitoh & Makino (2009).

 

5. ASTROPHYSICAL PROCESSES

In addition to the basic hydrodynamic algorithms, we implement the following astrophysical processes in the EUNHA code: (1) non-adiabatic evolution of the gas particles through radiative cooling, (2) global reionization heating, (3) star formation, and (4) energy and metallicity feedbacks by supernova type II (SNII) explosion.

5.1. Radiative Heating and Cooling

Using the CLOUDY package (version 10.10; Ferland et al. 1998), we calculated the heating and cooling rates and tabulated them in four terms of the gas density (ρ), temperature (T), metallicity (Z), and redshift (z) combining the effects of Compton heating/cooling, inverse Compton cooling, atomic/molecular cooling, and background UV heating. We assume that the whole simulation box is instantaneously full of the UV photons emitted by massive Pop-III stars at zre = 8.9 (Haardt & Madau 1996). With this step-function like global reionization process, we calculate the collisional ionization before zre and the photoionization after zre. Also to reflect the self-shielding effect when a dense gas cloud is optically thick against the background UV radiation, we adopt the critical hydrogen number density, , above which the UV background radiation is effectively blocked. In this study, we set following Tajiri & Umemura (1998) and Sawala et al. (2010). Figure 1 shows the heating/cooling rates of UV exposed/shielded gas. The UV-irradiated gas is gradually heated up to the equilibrium temperature Teq, at which the UV heating and molecular cooling are balanced with each other. However, the UV-shielded gas may continuously be cooled down.

Figure 1.Radiative heating/cooling rates of UV exposed/shielded gas as a function of temperature. The black lines denote the radiative heating/cooling curves when gas are irradiated by the UV background radiation, while the red line denotes the collisional cooling curve where the UV is shielded. Except for the existence of the UV radiation, the other conditions are the same as nH = 0.03 cm −3, Z=1 Z☉, and z = 8.

5.2. Star Formation

During the simulation run, we transform a gas particle into a star particle when the gas particle satisfies all the following star formation criteria (Katz 1992): (1) nH ＞ 0.1 cm−3, (2) cosmic virialization condition or ρg ＞ 57.5 ⟨ρg⟩, where ⟨ρg⟩ is the global value of gas density ρg at the redshift, (3) T ＜ 104 K, and (4) ∇·v ＜ 0 for the convergent flow.

The star formation rate is calculated according to the Schmidt law as

where ρ* is the density of newly born stars, c* is the characteristic star formation efficiency, and tdyn is the dynamical timescale of the gas particle. In this study, we adopt the standard star formation efficiency (c* = 0:0333) as adopted by Abadi et al. (2003). The dynamical time scale of a gas particle is defined as

The probability of a gas particle to be a star particle is given by the exponential law as

and for each time step we determine whether a gas particle satisfying all the star formation criteria would be converted to a star particle by generating a random number and applying the probability function in Equation (32). As the mass resolution of cosmological simulations may not reach individual stellar mass, each star particle actually represents a star cluster whose individual stars follow the stellar mass function of Kroupa (2001) with the mass range of 0.1–100 M☉.

5.3. Supernova Feedback

Stars follow different evolutionary tracks for different stellar mass and metallicity. As the SNII explosion releases a large amount of energy and metals, and affects the phase state of the interstellar medium, it would be important to include the SNII explosions in a high resolution cosmological simulation. For each star particle, we calculate the number of massive stars which may eventually end up as SNII explosions using the stellar evolutionary model of Hurly, Pols & Tout (2000). Similar to the star formation recipe, we consider the SNII feedback using the probability distribution (Okamoto, Nemmen, & Bower 2008) as

where rSN(t,Z) is the SNII explosion rate as a function of time and the metallicity (Z). tmax(Z) is the maximum life time of the massive stars. Typically, tmax(Z) is known to be less than ten million years.

We assume that the supernova explosion affects nearby gas particles through heating and metal enrichment. For this purpose the SPH-like scheme is adopted to distribute the energy and metals to the nearby gas particles. The total amount of energy released by a supernova explosion is fixed to 1051 erg.

However the cooling time scale of the shock-heated gas is usually less than the hydro subtime step. Then, even though there is a temperature balance between radiative cooling and background heating, the simulated gas would not attain the equilibrium. This is the time-resolution problem in the gas cooling. To overcome this resolution problem, we precalculated the temperature evolutions with much finer time steps as

where Γ is the cooling rate, Λ is the heating rate, and the temperature time step size is numerically set as dtT = dt/100. At every global time step, we tabulate the temperature evolution as a function of the gas density, temperature, metallicity, and time interval (dt). The metal enrichment is measured based on the tables given in Woosley & Weaver (1995).

 

6. CODE TESTS

To test the SPH algorithms and the implemented astrophysical processes, we perform five test simulations which are the most prominent topics that have been studied : (1) one-dimensional Riemann problems, (2) Kelvin-Helmholtz instability, (3) three-dimensional blast shock wave, (4) star formation on the isolated galactic disk, and (5) global star formation history in the cosmological context. Since the first three tests are designed to verify the hydrodynamics, we turn off gravity and astrophysical processes. We assume static backgrounds except in the last case. In the last two tests, we include gravity and astrophysical processes. Cosmic expansion is included only in the last case.

6.1. One-Dimensional Riemann Problems

In this subsection, we consider three cases of the onedimensional Riemann problem in various shock conditions: a weak shock (Problem 1), a strong rarefaction shock (Problem 2), and an extremely strong shock (Problem 3). We adopt the initial conditions given in Springel (2010) and briefly describe them below. Periodic boundary conditions are imposed with 0 ≤ x ＜ 1. The initial conditions simulating a shock front are chosen in the form of discontinuity in density (ρ), pressure (P), and particle velocity (v). Here, the subscript L denotes the region 1 (x ≤ 0.5) and subscript R is used for region 2 (x ＞ 0.5). For Problem 1, the initial conditions are (ρ, P, v) = (1, 1, 0)L and (0.125, 0.1, 0)R for Region 1 and 2, respectively. For Problem 2, we set (1, 0.4, −2)L and (1, 0.4, 2)R. And for Problem 3, we arrange (1, 1000, 0)L and (1, 0.01, 0)R. We used γ = 1.4 in these tests and the number of neighbors is fixed to 7.

We compare the SPH results with the analytical so- lutions of Sod (1978) at t = 0:05 in Figure 2. The den- sity, particle velocity, and pressure are shown for Problem 1 (left column), Problem 2 (middle), and Problem 3(right), respectively. In the weak and strong rarefaction shock model, the simulated profiles well match the analytic solutions. However, there are scatters at the contact discontinuity and the shock front is not sharp due to the smoothing features of the SPH. Also, in the extremely-strong shock model we observe that the simulated profiles of the density, velocity, and pressure show larger deviations from the analytic solutions between the rarefaction and the shock front. However, these features are also seen in the standard SPH simulation as shown in other papers.

Figure 2.Three comparisons of the one-dimensional Riemann shock tube problem between analytic solution (red line) and the simulations (blue dots).

6.2. Kelvin-Helmholtz Instability

Agertz et al. (2007) have argued that the standard SPH algorithm cannot correctly model contact discontinuities like the Kelvin-Helmholtz (KH) instability. There have been several attempts to properly handle the instability by modifying the SPH algorithms (Price 2008; Wadsley et al. 2008; Read et al. 2010). We examine the issue of contact discontinuities in the case of the shear flows.

The initial condition of the test is laid out to satisfy the periodic boundary conditions of the cubic box on a side length of Lbox = 0.001 Mpc, in which we distribute 1283 gas particles. We divide the simulation box into two regions, the central region (hereafter Region 1) of | y/Lbox − 0.5 |＜ 0.25 and the outer region (Region 2). Region 1 has a density ratio of ρ/ρc = 2, where ρc = 7.22 × 107 M☉/kpc3, an x-directional velocity of vx = 0.5 Vs, where Vs = 8.8 km/s. The outer region is set to have ρ/ρc = 1 and vx = −0.5Vs. Region 1 has a temperature 50000 K and the temperature of Region 2 is doubled to maintain pressure equilibrium at the contact discontinuity.

We perturb the equilibrium by adding a y-directional velocity vy as

where X ≡ x/Lbox, Y ≡ y/Lbox, ω0 = 0.1, and following Springel (2010).

Figure 3 shows the temporal evolution of density field at three different epochs, t = 0.25, 0.75 and 1.25 Ts, where Ts = 55.52 Myr. As the vertical perturbation grows, the horizontal (x-direction) shear ows generate vortex structures around the contact plane. The expected whirlpool-like structure gets prominent with time. It is well known that the standard SPH algorithm has a problem in reproducing the shear vortex of the KH instability (McNally et al. 2012; Hubber et al. 2013; Read et al. 2010). Many authors have suggested that a modification to the SPH algorithms would be required to produce KH vortex, for example, the artificial thermal conductivity (Price 2008; Agertz et al. 2007), a suitable smoothing kernel (Valcke et al. 2010), the Gudunov-SPH formalism (Cha et al. 2010), a diffusion term (Wadsley et al. 2008), the pressure-entropy formulation (Hopkins 2013), or the moving mesh (Heß & Springel 2010). Our findings are similar to that of Price (2008) who argued that particle noise at the discontinuity may suppress the growth of instability. The typical method to identify neighbor particles is to use the predict-correct algorithm, which is sometimes erroneous leading to noise. However, it is too far-fetched to tell whether our implementation of the neighbor searching method would solve the KH instability problem. The only thing we can say is that our neighbor finding may be one of the “partial” solutions to the KH instability problem. Further tests would be required to quantitatively compare EUNHA results with others.

Figure 3.Density maps of ows at t = 0.25, 0.75, and 1.25 Ts. The red color denotes high density region moving right, while the blue color denote low density region moving to the left. The vortex-like structures become more prominent with time.

6.3. Three-Dimensional Blast Wave

We also investigated the three-dimensional Sedov blast wave test. A nearly homogeneous glass-like distribution of 2563 particles is created in the 100 kpc simulation box with total particle mass 7.22 × 107M☉. The temperature of a gas particle is uniformly set to 10 K. We assign a huge temperature T = 106 K to the particle located at the center of the box to simulate the supernova explosion. This generates a shock wave, which propagates outwards and sweeps away the surrounding cool and less-dense gas.

Figure 4 shows the comparison of density, temperature, and pressure profiles around the ground zero of explosion between the simulation (points) and semi-analytic solution (line; Sedov 1959) at t = 120.5 Myr. Good agreement is seen both in the location of the shock front and the maximum density, with a substantial scatter around the analytic solution. The simulated upstream has a scatter of density with a finite slope because of the finite smoothing scale of the SPH and the anisotropic initial particle distribution. The regular lattice point of initial particle positions was adopted by Springel & Hernquist (2002) and Merlin et al. (2010) while a new setup is adopted to reduce anisotropic and inhomogeneous distribution in the shock propagation (Rosswog & Price 2007).

Figure 4.Density (top), temperature (middle) and pressure (bottom panel) profiles of three-dimensional Sedov blast wave. Here, r is the distance from the central particle. The red dots denote the simulation results and the black line shows the semi-analytic solution of Sedov (1959). The x-axis is scaled to the box size and the y-axis is normalized to the initial density, temperature and pressure of the test simulation.

6.4. Isolated Galaxy

Now, we investigate the stability of galaxy disk and the star formation rate on the disk plane using the isolated galaxy model. The initial conditions of a multi-component galaxy were generated according to the isolated Milky-Way model proposed by Hernquist (1993). An exponential gas disk is arranged with a characteristic mass Mg = 4.125 × 109M☉, disk scale height z = 0.3 kpc, and disk scale length of h = 3.33 kpc. An exponential stellar disk is built with Ms = 3.715× 1010 M☉, z = 0.3 kpc, and h = 3.33 kpc. The bulge component follows Hernquist profile with Mb = 1.375 × 1010 M☉ and scale length of a = 0.8 kpc. The halo component follows the Hernquist model with parameters Mh = 2.2 × 1011 M☉ and a = 10 kpc. The three-dimensional locations and velocities of particles are calculated by ZENO software package developed by Joshua E. Barnes1. The gas and stellar disks are composed of 983,040 particles with particle mass of Mp = 4.196 × 104 M☉ while the bulge and halo struc- tures are treated as fixed external potentials. We set the initial temperature of gas particles to be 104 K.

Figure 5 shows the density-weighted temperature map of the gas disk at t = 0.5 Gyr. Gas particles of T ＜ 104 K are distributed along spiral structures, where density is so high that the gas is able to cool down to lower temperature through efficient radiative cooling. Density-temperature distribution of gas particles under radiative heating and cooling process are shown in Figure 6. We can see that gas particles in the low density regions keep the initial temperature of T ~ 104 K, while the gas particles in the higher density regions may cool down to lower temperature toward the equilibrium state. Stars form in the high-density and low-temperature regions along the spirals. Tens of million years after star formation, the SNII explosions of the massive stars redistribute thermal energy into the surrounding interstellar medium. The red clumps in Figure 5 represent the highly heated gas particles near the supernova explosion sites.

Figure 5.Density-weighted temperature map of gas disk. The dense spiral arms are cooled down more efficiently than the less-dense regions through the more efficient radiative cooling. Red clumps along the spiral structures represent the shock-heated gas by the SNII explosions among the recently-forming massive stars.

Figure 6.nH–T diagram of gas particles at T = 0.5 Gyr. All gas particles are irradiated by UV backgrounds because we assume that there is no self-shielding in this test. Thus, most of simulated particles are distributed along the equilibrium temperature line. The scatters around the equilibrium temperature is mainly due to adiabatic expansion and contraction. Points in the upper-right corner are the shock-heated gas particles due to the nearby SNII explosions. After a transient period, their density will drop due to expansion and points will move to the upper-left corner of the box.

Finally, we compare the simulated star formation rate with the observed value in Figure 7. The projected star formation rates on the gas disk ∑SFR for a given column density ∑gas well reproduce the observed Schmidt-Kennicutt relation (Kennicutt 1998).

Figure 7.Projected star formation rate ∑SFR versus the column density ∑gas on the simulated disk at 0.5 Gyr. The solid line enclosed by the gray-shaded region represents the observed Kennicutt relation with errors (Kennicutt 1998), and squares are the simulation results.

6.5. Cosmological Model

Now, we turn our attention to the application of EUNHA to cosmological star formation. We solve the equation of motion of the gas and dark matter particles in a cubic simulation box of a side length Lbox = 32 h−1Mpc. The total number of gas and dark matter particles is 5123. The cosmological model adopted in the simulation is the WMAP-5 year cosmology with parameters of Ωm = 0:26, ΩΛ = 0.74, Ωb = 0.044, σ8 = 0.76, ns = 0.96, and h = 72 km s−1 Mpc−1. The simulation starts with glass-like pre-initial condition, which has a virtually null net gravitational force. To perturb pre-initial particle distribution, we applied the second-order linear perturbation method described in Jenkins (2010). The initial power spectra for baryonic and dark matter at redshift z = 150 are calculated by the CAMB package (http://camb.info/sources). The initial background temperature of gas is calculated using the RecFast package (Seager, Sasselov, & Scott 1999). In generating initial conditions, the density fluctuations of the gas particles may create temperature fluctuations. We determine the temperature fluctuation for each gas particle using the adiabatic model as

Figure 8 is an example of the gas density at z = 6. The distribution of gas particles by and large follow those of dark matter particles such as clusters, filamentary structures, and cosmic voids. After the predefined cosmic reionization epoch at zre = 8.9, most of gas particles in low density regions of nH ＜ 0.014 cm−3 are heated up to Tg ~ 104 K, thus the distribution of the gas particles is more diffuse than dark matter particles. Meanwhile, gas particles located in high density regions (nH ＞ 0.014 cm−3), which are optically thick to the cosmic UV background, may cool down to lower temperature and, thus, settle down to the inner region of dark matter halos. The detailed density-temperature relation of gas particles are shown in Figure 9. Due to the gas inflow to the halo center, typically halo regions have ∇ · v ＜ 0. And as the inner regions of the halo have a high hydrogen density and low temperature, gas particles are able to cool and finally form stars therein. The red dots in Figure 8 represent the newly generated stars in clustered regions.

Figure 8.Projected density map of gas particles in a z-directional slab with a width equal to three times of the mean particle separation at the epoch of z = 6. Red dots represent newly formed star particles, and they are located in cluster regions.

Figure 9.nH–T diagram of gas particles at z = 6. The solid line marks the UV-background shielding density of nH = 0.014 cm −3.

The observed global star formation rate (ρSFR) of the Universe is known to be a function of redshift with a peak located among z = 2−3. In Figure 10, we compare the star formation history of simulations (lines) with observations (symbols, Hopkins & Beacom 2006). The overall star formation history agrees with each other.

Figure 10.Comparison of star formation histories of the universe between the simulation results (three different lines) and the observations (shaded symbols with error bars). The simulation results are differed by resolutions, different box size and particle number. The observed data is compiled by Hopkins & Beacom (2006).

However, there are differences in the gradient of star formation rate between three different simulations, and it is because of the different resolution. Choi & Nagamine (2012) reported that the lower mass galaxies are preferred sites for star formation in the higher redshift, thus slope of ρSFR(z) for the finer resolution is shallower than others. Among cosmological simulations, the highest resolution simulation with Lbox = 4 h−1Mpc and 5123 particles show the lowest but highest ρSFR.

 

7. SUMMARY AND CONCLUSION

We have developed a new cosmological hydrodynamic simulation code (EUNHA) by combining the pre-existing gravitational N-body code of GOTPM (Dubinski et al. 2004) with the standard smoothed particle hydrodynamics (SPH) algorithm. This code fully exploits the advantage of the Oct-Sibling Tree (OST) of the GOTPM to identify the N nearest neighbors.

For astrophysical evolution of gas particles, we also implemented the processes of (1) non-adiabatic evolution of the gas particles through radiative heating/cooling, (2) global reionization, (3) star formation, and (4) energy and metallicity feedbacks by supernova type II explosion. To demonstrate our new implementations of the SPH and the astrophysical processes, we made five test simulations: (1) one-dimensional Riemann problems, (2) Kelvin-Helmholtz instability, (3) three-dimensional blast shock wave, (4) star formation on the isolated galactic disk, and (5) global star formation history in the cosmological context.

It is interesting to see the growth of the shear vortex in the Kelvin-Helmholtz instability test because the only difference in EUNHA from other standard SPH methods is the neighbor searching. Even though it is premature to conclude in this simple test, we may get a hint from Price (2008), who showed that the particle noise on the shear contact plane may be one of the possible causes to the suppressed KH instabilities. Our improved neighbor searching method may reduce this kind of neighboring noise in the predict-correct method. The SPH measures hydrodynamic quantities of a gas particle using nearby interacting neighbors contained by a finite smoothing length. Therefore, the gas density, pressure, and corresponding acceleration may change for different smoothing lengths. It means that during a simulation run a sudden change in the number of interacting neighbors may also have a sudden impact on the gas motion. Then, small perturbations may be buried in the numerical noise and may be suppressed from developing on the contact layer. However, it is too hasty to jump to the conclusion with this simple result because some authors already showed that with different initial conditions the KH vortex can grow even in the standard-SPH test (Hopkins 2013) but with much less surface mixing than shown in other grid-based tests. Also Price (2008) showed that a standard SPH test may yield a small KH vortex by adjusting the test parameters. So, it is unclear whether our implementation of the advanced neighbor searching solves the KH instability problem in the Lagrangian code. Our tentative conclusion is that our code may develope the KH vortex but the surface mixing on the contact plane seems to be not so strong as other improved versions of the SPH or Eulerian codes. More rigorous tests will be necessary to draw any conclusion on this issue.

The EUNHA code is originally intended to serve as the basis for the further development of the cosmological hydrodynamic simulation code. Most hydrodynamic routines are built for fast and efficient communication between the parallel ranks, and the SPH functions is carefully organized to be isolated (or modularized) for exibility in case of future updates. Therefore, we may be able to easily exchange the SPH routines with those of another type of the particle-based hydrodynamic algorithms.