Advanced
EUNHA: A NEW COSMOLOGICAL HYDRODYNAMIC SIMULATION CODE
EUNHA: A NEW COSMOLOGICAL HYDRODYNAMIC SIMULATION CODE
Journal of The Korean Astronomical Society. 2014. May, 47(3): 87-98
Copyright © 2014, null
  • Received : September 13, 2013
  • Accepted : February 11, 2014
  • Published : May 20, 2014
Download
PDF
e-PUB
PubReader
PPT
Export by style
Share
Article
Author
Metrics
Cited by
TagCloud
About the Authors
JIHYE SHIN
School of Space Research, Kyung Hee University, Yongin, Kyeonggi 446-701, Korea
JUHAN KIM
Center for Advanced Computation, Korea Institute for Advanced Study, 85 Hoegiro, Dongdaemun-gu, Seoul 130-722, Korea;kjhan@kias.re.kr
SUNGSOO S. KIM
Department of Astronomy and Space Science, Kyung Hee University, Yongin, Kyeonggi 446-701, Korea
CHANGBOM PARK
School of Physics, Korea Institute for Advanced Study, 85 Hoegiro, Dongdaemun-gu, Seoul 130-722, Korea

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.
Keywords
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
PPT Slide
Lager Image
where am is the scale factor at z = 0. In the GOTPM code, we use the following variables
PPT Slide
Lager Image
PPT Slide
Lager Image
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
PPT Slide
Lager Image
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
PPT Slide
Lager Image
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,
PPT Slide
Lager Image
where a s is defined as the acceleration in simulation unit,
PPT Slide
Lager Image
),
PPT Slide
Lager Image
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;
PPT Slide
Lager Image
PPT Slide
Lager Image
Here, the initial matter density, Ω i is given by
PPT Slide
Lager Image
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,
PPT Slide
Lager Image
is obtained as
PPT Slide
Lager Image
where g 1 = H(z)ℓb =(1 + z ). The peculiar velocity ( v p ) is obtained from the simulation velocity ( v s ) as
PPT Slide
Lager Image
where g 2 = ag 1 .
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 EU N HA (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
PPT Slide
Lager Image
PPT Slide
Lager Image
where mH is the mass of hydrogen atom, 𝜇 is the mean molecular weight, k B is the Boltzmann constant, and γ is the adiabatic index. In an adiabatic process, the pressure is related to the gas density as
PPT Slide
Lager Image
, where Ar is the entropy.
Now we adopt the following relations of conversion between the physical and simulation units,
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
where ⟨ ρz ⟩ = ⟨ ρ 0 ⟩ (1 + z ) 3 = Ω m
PPT Slide
Lager Image
(1 + z ) 3 and
PPT Slide
Lager Image
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;
PPT Slide
Lager Image
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
PPT Slide
Lager Image
In this study, we used the spline kernel discussed in Monaghan (1992).
We have adopted the typical artificial viscosity effect to the acceleration as
PPT Slide
Lager Image
PPT Slide
Lager Image
where 𝚷 ij is the viscosity factor introduced to capture the shock front. We adopt the form proposed by Monaghan (1997)
PPT Slide
Lager Image
where
PPT Slide
Lager Image
and α is viscosity coefficient. We adopt α = 1 in this study. The simulation sound speed, cs , is defined as
PPT Slide
Lager Image
and this is identical to the physical sound speed, cr . Finally, we get the viscous force as
PPT Slide
Lager Image
- 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
PPT Slide
Lager Image
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
PPT Slide
Lager Image
The time-step size for a given gas particle is determined as (Springel 2005)
PPT Slide
Lager Image
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
PPT Slide
Lager Image
where
PPT Slide
Lager Image
is the normalized mutual displacement vector and v ij is the relative velocity. Here, Vsig is the maximum value among
PPT Slide
Lager Image
’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. MergingN-body and Hydrodynamic Subtime Steps
For the EU N HA 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;
PPT Slide
Lager Image
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(log 2 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 EU N HA 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 (SN II ) 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,
PPT Slide
Lager Image
, above which the UV background radiation is effectively blocked. In this study, we set
PPT Slide
Lager Image
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 T eq , at which the UV heating and molecular cooling are balanced with each other. However, the UV-shielded gas may continuously be cooled down.
PPT Slide
Lager Image
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) n H > 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 < 10 4 K, and (4) ∇· v < 0 for the convergent flow.
The star formation rate is calculated according to the Schmidt law as
PPT Slide
Lager Image
where ρ * is the density of newly born stars, c * is the characteristic star formation efficiency, and t dyn 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
PPT Slide
Lager Image
The probability of a gas particle to be a star particle is given by the exponential law as
PPT Slide
Lager Image
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 SN II 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 SN II explosions in a high resolution cosmological simulation. For each star particle, we calculate the number of massive stars which may eventually end up as SN II explosions using the stellar evolutionary model of Hurly, Pols & Tout (2000). Similar to the star formation recipe, we consider the SN II feedback using the probability distribution (Okamoto, Nemmen, & Bower 2008) as
PPT Slide
Lager Image
where r SN ( t,Z ) is the SN II explosion rate as a function of time and the metallicity ( Z ). t max ( Z ) is the maximum life time of the massive stars. Typically, t max ( 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 10 51 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
PPT Slide
Lager Image
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.
PPT Slide
Lager Image
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 L box = 0.001 Mpc, in which we distribute 128 3 gas particles. We divide the simulation box into two regions, the central region (hereafter Region 1) of | y / L box − 0.5 |< 0.25 and the outer region (Region 2). Region 1 has a density ratio of ρ/ρc = 2, where ρc = 7.22 × 10 7 M /kpc 3 , 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.5 Vs . 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
PPT Slide
Lager Image
where X x/Lbox , Y y/Lbox , ω 0 = 0.1, and
PPT Slide
Lager Image
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 EU N HA results with others.
PPT Slide
Lager Image
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 256 3 particles is created in the 100 kpc simulation box with total particle mass 7.22 × 10 7 M . The temperature of a gas particle is uniformly set to 10 K. We assign a huge temperature T = 10 6 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).
PPT Slide
Lager Image
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 × 10 9 M , 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× 10 10 M , z = 0.3 kpc, and h = 3.33 kpc. The bulge component follows Hernquist profile with Mb = 1.375 × 10 10 M and scale length of a = 0.8 kpc. The halo component follows the Hernquist model with parameters Mh = 2.2 × 10 11 M and a = 10 kpc. The three-dimensional locations and velocities of particles are calculated by ZENO software package developed by Joshua E. Barnes 1 . The gas and stellar disks are composed of 983,040 particles with particle mass of Mp = 4.196 × 10 4 M while the bulge and halo struc- tures are treated as fixed external potentials. We set the initial temperature of gas particles to be 10 4 K.
Figure 5 shows the density-weighted temperature map of the gas disk at t = 0.5 Gyr. Gas particles of T < 10 4 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 ~ 10 4 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 SN II 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.
PPT Slide
Lager Image
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.
PPT Slide
Lager Image
nHT 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).
PPT Slide
Lager Image
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 EU N HA 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 L box = 32 h −1 Mpc. 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
PPT Slide
Lager Image
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 n H < 0.014 cm −3 are heated up to Tg ~ 10 4 K, thus the distribution of the gas particles is more diffuse than dark matter particles. Meanwhile, gas particles located in high density regions ( n H > 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.
PPT Slide
Lager Image
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.
PPT Slide
Lager Image
nHT 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.
PPT Slide
Lager Image
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
PPT Slide
Lager Image
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 L box = 4 h −1 Mpc and 512 3 particles show the lowest
PPT Slide
Lager Image
but highest ρ SFR .
7. SUMMARY AND CONCLUSION
We have developed a new cosmological hydrodynamic simulation code (EU N HA) 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 EU N HA 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 EU N HA 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.
http://www.ifa.hawaii.edu/~barnes/zeno/index.html
Acknowledgements
Authors thank the anonymous referee for his/her invaluable comments on the draft, which help us enhance the consistency of the content. This work was supported by the BK21 plus program through the National Research Foundation (NRF) funded by the Ministry of Education of Korea. SSK and JS’s works were supported by Mid-career Research Program (No. 2011-0016898) through the NRF grant funded by the Ministry of Science, ICT and Future Planning of Korea. JS deeply appreciates Jeong-Sun Hwang for her help in making the initial conditions for compound model galaxies.
References
Abadi Mario G. , Navarro Julio F. , Steinmetz Matthias , Eke Vincent R. 2003 Simulations of Galaxy Formation in a Cold Dark Matter Universe. I. Dynamical and Photometric Properties of a Simulated Disk Galaxy ApJ 591 499 -    DOI : 10.1086/375512
Agertz O. , Moore B. , Stadel J. , Potter D. , Miniati F. , Read J. , Mayer L. , Gawryszczak A. , Kravtsov A. , Nordlund A , Pearce F. , Quilis V. , Rudd D. , Springel V. , Stone J. , Tasker E. , Teyssier R. , Wadsley J. , Walder R. 2007 Fundamental Differences between SPH and Grid Methods MNRAS 380 963 -    DOI : 10.1111/j.1365-2966.2007.12183.x
Aihara H. , Allende P. C. , An D. , Anderson S. F. , Aubourg E. , Balbinot E. , Beers T. C. , Berlind A. A. , Bickerton S. J. , Bizyaev D. 2011 The Eighth Data Release of the Sloan Digital Sky Survey: First Data from SDSS-III ApJSS 193 29 -    DOI : 10.1088/0067-0049/193/2/29
Angulo R. E. , Springel V. , White S. D. M. , Jenkins A. , Baugh C. M. , Frenk C. S. 2012 Scaling Relations for Galaxy Clusters in the Millennium-XXL Simulation
Cha S. , Inutsuka S. , Nayakshin S. 2010 Kelvin-Helmholtz Instabilities with Godunov Smoothed Particle Hydrodynamics MNRAS 403 1165 -    DOI : 10.1111/j.1365-2966.2010.16200.x
Choi J.-H. , Nagamine K. 2012 On the Inconsistency be-tween the Estimates of Cosmic Star Formation Rate and Stellar Mass Density of High-Redshift Galaxies MNRAS 419 1280 -    DOI : 10.1111/j.1365-2966.2011.19788.x
Choi Y.-Y. , Kim J. , Rossi G. , Kim S. S. , Lee J.-E. 2013 Topology of Luminous Red Galaxies from the Sloan Digital Sky Survey ApJS 209 19 -    DOI : 10.1088/0067-0049/209/2/19
Choi Y.-Y. , Park C. , Kim J. , Gott J. R. , Weinberg D. H. , Vogeley M. S. , Kim S. S. 2010 Galaxy Clustering Topology in the Sloan Digital Sky Survey Main Galaxy Sample: A Test for Galaxy Formation Models ApJS 190 181 -    DOI : 10.1088/0067-0049/190/1/181
Colless M. , Dalton G. , Maddox S. , Sutherland W. , Norberg P. , Cole S. , Bland-Hawthorn J. , Bridges T. , Can-non R. , Collins C. 2001 The 2dF Galaxy Redshift Survey: Spectra and Redshifts MNRAS 328 1039 -    DOI : 10.1046/j.1365-8711.2001.04902.x
Dubinski J. , Kim J. , Park C. , Humble R. 2004 GOTPM: a Parallel Hybrid Particle-Mesh Treecode New Astronomy 9 111 -    DOI : 10.1016/j.newast.2003.08.002
Ferland G. J. , Korista K. T. , Verner D. A. , Ferguson J. W. , Kingdon J. B. , Verner E. M. 1998 CLOUDY 90: Numerical Simulation of Plasmas and Their Spectra PASP 110 761 -
Fryxell B. , Olson K. , Ricker P. , Timmes F. X. , Zingale M. , Lamb D. Q. , MacNeice P. , Rosner R. , Truran J. W. , Tufo H. 2000 ApJSS 131 273 -    DOI : 10.1086/317361
Haardt F. , Madau P. 1996 Radiative Transfer in a Clumpy Universe. II. The Ultraviolet Extragalactic Background ApJ 461 20 -    DOI : 10.1086/177035
Hernquist L. 1993 N-Body Realizations of Compound Galaxies ApJSS 86 38 -
Hess S. , Springel V. 2010 Particle Hydrodynamics with Tessellation Techniques MNRAS 406 2289 -    DOI : 10.1111/j.1365-2966.2010.16892.x
Hopkins A. M. , Beacom J. F. 2006 On the Normalization of the Cosmic Star Formation History ApJ 651 142 -    DOI : 10.1086/506610
Hopkins P. F. 2013 A General Class of Lagrangian Smoothed Particle Hydrodynamic Methods and Implications for Fluid Mixing Problems MNRAS 428 2840 -    DOI : 10.1093/mnras/sts210
Hubber D. A. , Batty C. P. , McLeod A. , Whitworth A. P. 2011 SEREN-A New SPH Code for Star and Planet Formation Simulations A&A 529 27 -    DOI : 10.1051/0004-6361/201014949
Hubber D. A. , Falle S. A. E. G. , Goodwin S. P. 2013 Convergence of AMR and SPH Simulaitons-I. Hydro-dynamical Resolution and Convergence Tests MNRAS 432 711 -    DOI : 10.1093/mnras/stt509
Hurley J. R. , Pols O. R. , Tout C. A. 2000 Comprehensive Analytic Formulae for Stellar Evolution as a Function of Mass and Metallicity MNRAS 315 543 -    DOI : 10.1046/j.1365-8711.2000.03426.x
Jenkins A. 2010 Second-Order Lagrangian Perturbation Theory Initial Conditions for Resimulations MNRAS 403 1859 -    DOI : 10.1111/j.1365-2966.2010.16259.x
Katz N. 1992 Dissipational Galaxy Formation. II-Effects of Star Formation ApJ 391 502 -    DOI : 10.1086/171366
Kennicutt R. C. 1998 The Global Schmidt Law in Star-Forming Galaxies ApJ 498 541 -    DOI : 10.1086/305588
Kim J. , Park C. , Gott J. R. III. , Dubinski J. 2009 The Horizon Run N-Body Simulation: Baryon Acoustic Oscillations and Topology of Large-scale Structure of the Universe ApJ 701 1547 -    DOI : 10.1088/0004-637X/701/2/1547
Kim J. , Park C. , Rossi G. , Lee S. , Gott J. R. III 2011 The New Horizon Run Cosmological N-Body Simulations JKAS 44 (6) 217 - 234    DOI : 10.5303/JKAS.2011.44.6.217
Kravtsov A. V. 1999 High-Resolution Simulations of Structure Formation in the Universe, PhD thesis New Mexico State Univ.
Kroupa P. 2001 On the Variation of the Initial Mass Function MNRAS 322 231 -    DOI : 10.1046/j.1365-8711.2001.04022.x
Kuhlen M. , Vogelsberger M. , Angulo R. 2012 Numerical Simulations of the Dark Universe: State of the Art and the Next Decade Physics of the Dark Unverse 1 50 -    DOI : 10.1016/j.dark.2012.10.002
Masaki S. , Hikage C. , Takada M. , Spergel D. N. , Sugiyama N. 2013 Understanding the Nature of Luminous Red Galaxies (LRGs): Connecting LRGs to Central and Satellite Subhalos MNRAS 436 2286 -    DOI : 10.1093/mnras/stt1729
McNally C. P. , Lyra W. , Passy J.-C. 2012 AWell-Posed Kelvin-Helmholtz Instibility Test and Comparison ApJS 201 18 -    DOI : 10.1088/0067-0049/201/2/18
Merlin E. , Buonomon U. , Grassi T. , Piovan L. , Chiosi C. 2010 EvoL: the New Padove Tree-SPH Parallel Code for Cosmological Simulations A&A 513 36 -    DOI : 10.1051/0004-6361/200913514
Monaghan J. J. 1992 Smoothed Particle Hydrodynamics ARA&A 30 54 -
Monaghan J. J. 1997 SPH and Riemann Solvers Comput. Phys. 136 298 -    DOI : 10.1006/jcph.1997.5732
Navarro J. F. , Frenk C. S. , White S. D. M. 1996 The Structure of Cold Dark Matter Halos ApJ 463 563 -
Nurmi P. , Heinamaki P. , Sepp T. , Tago E. , Saar E. , Gramann M. , Einasto M. , Tempel E. , Einasto J. 2013 Groups in the Millennium Simulation and in SDSS DR7 MNRAS 436 380 -    DOI : 10.1093/mnras/stt1571
Okamoto T. , Nemmen R. S. , Bower R. G. 2008 The Impact of Radio Feedback from Active Galactic Nuclei in Cosmological Simulations: Formation of Disc Galaxies MNRAS 385 161 -    DOI : 10.1111/j.1365-2966.2008.12883.x
O’Shea B. W. , Bryan G. , Bordner J. , Norman M. L. , Abel T. , Harkness R. , Kritsuk A. 2004 Introducing Enzo, an AMR Cosmology Application
Park C. , Choi Y.-Y. , Kim J. , Gott J. R. , Kim S. S. , Kim K.-S. 2012 The Challenge of the Largest Structures in the Universe to Cosmology ApJ 759 7 -    DOI : 10.1088/0004-637X/759/1/7
Park C. , Hwang H. S. 2009 Interactions of Galaxies in the Galaxy Cluster Environment ApJ 699 1595 -    DOI : 10.1088/0004-637X/699/2/1595
Price D. J. 2008 Modelling Discontinuities and Kelvin Helmholtz Instabilities in SPH JCoPh 227 10040 -
Rasera Y. , Corasaniti P.-S. , Alimi J.-M. , Bouillot V. , Reverdy V. , Balmes I. 2013 Cosmic Variance Lim-ited Baryon Acoustic Oscillaions from the DEUS-FUR ACDM Simulation
Read J. I. , Hayfield T. , Agertz O. 2010 Resolving Mixing in Smoothed Particle Hydrodynamics MNRAS 405 1513 -
Rosswog S. , Price D. 2007 MAGMA: a Three-Dimensional, Lagrangian Magnetohydrodynamics Code for Merger Applications MNRAS 379 915 -    DOI : 10.1111/j.1365-2966.2007.11984.x
Saitoh T. R. , Makino J. 2009 A Necessary Condition for Individual Time Steps in SPH Simulations ApJ 697 99 -    DOI : 10.1088/0004-637X/697/2/L99
Sawala T. , Scannapieco C. , Maio U. , White S. 2010 Formation of Isolated Dwarf Galaxies with Feedback MNRAS 402 1599 -    DOI : 10.1111/j.1365-2966.2009.16035.x
Seager S. , Sasselov D. D. , Scott D. 1999 A New Calculation of the Recombination Epoch ApJ 523 L1 - 1    DOI : 10.1086/312250
Sedov L. I. 1959 Similarity and Dimensional Methods in Mechanics Academic Press New York
Sod G. A. 1978 A Survey of Several Finite Difference Methods for Systems of Nonlinear Hyperbolic Conservation Laws JocoPh 27 1 -
Springel V. , Hernquist L. 2002 Cosmological Smoothed Particle Hydrodynamics Simulations: the Entropy Equation MNRAS 333 649 -    DOI : 10.1046/j.1365-8711.2002.05445.x
Springel V. 2005 The Cosmological Simulation Code GADGET-2 MNRAS 364 1105 -    DOI : 10.1111/j.1365-2966.2005.09655.x
Springel V. 2010 Smoothed Particle Hydrodynamics in Astrophysics ARA&A 48 391 -    DOI : 10.1146/annurev-astro-081309-130914
Tajiri Y. , Umemura M. 1998 A Criterion for Photoionization of Pregalactic Clouds Exposed to Diffuse Ultraviolet Background Radiation ApJ 502 59 -    DOI : 10.1086/305898
Tegmark M. 2006 Cosmological Constraints from the SDSS Luminous Red Galaxies Phy. Rev. D 74 123507 -    DOI : 10.1103/PhysRevD.74.123507
Teyssier R. 2002 Cosmological Hydrodynamics with Adaptive Mesh Refinement. A New High Resolution Code Called RAMSES A&A 385 337 -    DOI : 10.1051/0004-6361:20011817
Thacker R. J. , Tittley E. R. , Pearce F. R. , Couchman H. M. P. , Thomas P. A. 2000 Smoothed Particle Hydrodynamics in Cosmology: a Comparative Study of Im-plementations MNRAS 319 619 -    DOI : 10.1111/j.1365-8711.2000.03927.x
Valcke S. , de Rijcke S. , Rodiger E. , Dejonghe H. 2010 Kelvin-Helmholtz Instabilities in Smoothed Particle Hydrodynamics MNRAS 408 71 -    DOI : 10.1111/j.1365-2966.2010.17127.x
Wadsley J. W. , Veeravalli G. , Couchman H. M. P. 2008 On the Treatment of Entropy Mixing in Numerical Cos-mology MNRAS 387 427 -    DOI : 10.1111/j.1365-2966.2008.13260.x
Wadsley J. W. , Stadel J. , Quinn T. 2004 Gasoline: a Flexible, Parallel Implementation of TreeSPH New Astronomy 9 137 -    DOI : 10.1016/j.newast.2003.08.004
Wetzstein M. , Nelson A. F. , Naab T. , Burkert A. 2009 Vine-A Numerical Code for Simulating Astrophysical Systems Using Particles. I. Description of the Physics and the Numerical Methods ApJS 184 298 -    DOI : 10.1088/0067-0049/184/2/298
Woosley S. E. , Weaver Thomas A. 1995 The Evolution and Explosion of Massive Stars. II. Explosive Hydrodynamics and Nucleosynthesis ApJ 101 181 -    DOI : 10.1086/192237