EUNHA: A NEW COSMOLOGICAL HYDRODYNAMIC SIMULATION CODE

Journal of The Korean Astronomical Society.
2014.
Jun,
47(3):
87-98

- Received : September 13, 2013
- Accepted : February 11, 2014
- Published : June 30, 2014

Download

PDF

e-PUB

PubReader

PPT

Export by style

Article

Metrics

Cited by

TagCloud

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.
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.
x
) is related to the physical distance (
r
) with the scale factor (
a
) as
where
a_{m}
is the scale factor at
z
= 0. In the GOTPM code, we use the following variables
where
x_{s}
is the position of a particle in unit of the mean particle separation (
ℓ_{b}
≡
L_{b}/N
),
L_{b}
is the simulation box size on one side,
N
is the number of grids in one direction,
m_{r}
is the physical mass of the particle,
m_{s}
is the particle mass in simulation unit, and ⟨
ρ
_{0}
⟩ is the mean density at the current epoch (
a = a_{m}
). 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
a
_{s}
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
g
_{1}
=
H(z)ℓ_{b}
=(1 +
z
). The peculiar velocity (
v
_{p}
) is obtained from the simulation velocity (
v
_{s}
) as
where
g
_{2}
=
ag
_{1}
.
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).
P_{r}
) and specific internal energy (
u_{r}
) can be measured from the temperature (
T
) and density (
ρ_{r}
) of ideal gas as
where
m_{H}
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
, where
A_{r}
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 (
A_{s}
) may change with time even when
A_{r}
is unchanged with time.
where
h
is smoothing length defined as the distance to the
N_{n}
’th nearest neighbor,
W
is the smoothing kernel, and
f_{i}
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,
c_{s}
, is defined as
and this is identical to the physical sound speed,
c_{r}
. Finally, we get the viscous force as
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.
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,
h_{sml}
is the smoothing length, and
V_{sig}
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 v
_{ij}
is the relative velocity. Here,
V_{sig}
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.
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;
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
C_{s}
) 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).
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.
ρ
), 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
z_{re}
= 8.9 (Haardt & Madau 1996). With this step-function like global reionization process, we calculate the collisional ionization before
z_{re}
and the photoionization after
z_{re}
. 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
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.
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 n _{H} = 0.03 cm ^{−3}, Z=1 Z_{☉}, and z = 8.
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
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
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
_{☉}
.
_{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
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
where Γ is the cooling rate, Λ is the heating rate, and the temperature time step size is numerically set as
dt_{T}
=
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).
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.
Three comparisons of the one-dimensional Riemann shock tube problem between analytic solution (red line ) and the simulations (blue dots ).
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
v_{x}
= 0.5
V_{s}
, where
V_{s}
= 8.8 km/s. The outer region is set to have
ρ/ρ_{c}
= 1 and
v_{x}
= −0.5
V_{s}
. 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
v_{y}
as
where
X
≡
x/L_{box}
,
Y
≡
y/L_{box}
, ω
_{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
T_{s}
, where
T_{s}
= 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.
Density maps of ows at t = 0.25, 0.75, and 1.25 T_{s} . 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.
^{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).
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.
M_{g}
= 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
M_{s}
= 3.715× 10
^{10}
M
_{☉}
,
z
= 0.3 kpc, and
h
= 3.33 kpc. The bulge component follows Hernquist profile with
M_{b}
= 1.375 × 10
^{10}
M
_{☉}
and scale length of
a
= 0.8 kpc. The halo component follows the Hernquist model with parameters
M_{h}
= 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
M_{p}
= 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.
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 SN_{II} explosions among the recently-forming massive stars.
n _{H}–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 SN_{II} 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).
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.
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,
n_{s}
= 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
z_{re}
= 8.9, most of gas particles in low density regions of
n
_{H}
＜ 0.014 cm
^{−3}
are heated up to
T_{g}
~ 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.
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.
n _{H}–T diagram of gas particles at z = 6. The solid line marks the UV-background shielding density of n _{H} = 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.
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
L
_{box}
= 4
h
^{−1}
Mpc and 512
^{3}
particles show the lowest
but highest
ρ
_{SFR}
.
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.
1 http://www.ifa.hawaii.edu/~barnes/zeno/index.html

1. INTRODUCTION

Over the last few decades, the cosmological gravitational
2. BASIC EQUATIONS OF MOTION

In an expanding background, the comoving distance (
PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

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
- 3.1. Basic Equations

The pressure (
PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

- 3.2. Smoothed Particle Hydrodynamics

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

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

- 3.3. Neighbor Findings

We use the Oct-Sibling Tree (OST) to find the
4. INDIVIDUAL TIME STEP

- 4.1. Subtime step

In a particular time step (Δ
PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

- 4.2. MergingN-body and Hydrodynamic Subtime Steps

For the EU
PPT Slide

Lager Image

5. ASTROPHYSICAL PROCESSES

In addition to the basic hydrodynamic algorithms, we implement the following astrophysical processes in the EU
- 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 (
PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

- 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)
PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

- 5.3. Supernova Feedback

Stars follow different evolutionary tracks for different stellar mass and metallicity. As the SN
PPT Slide

Lager Image

PPT Slide

Lager Image

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 ≤
PPT Slide

Lager Image

- 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
PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

- 6.3. Three-Dimensional Blast Wave

We also investigated the three-dimensional Sedov blast wave test. A nearly homogeneous glass-like distribution of 256
PPT Slide

Lager Image

- 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
PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

- 6.5. Cosmological Model

Now, we turn our attention to the application of EU
PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

7. SUMMARY AND CONCLUSION

We have developed a new cosmological hydrodynamic simulation code (EU
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.

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**

Citing 'EUNHA: A NEW COSMOLOGICAL HYDRODYNAMIC SIMULATION CODE
'

@article{ CMHHBA_2014_v47n3_87}
,title={EUNHA: A NEW COSMOLOGICAL HYDRODYNAMIC SIMULATION CODE}
,volume={3}
, url={http://dx.doi.org/10.5303/JKAS.2014.47.3.87}, DOI={10.5303/JKAS.2014.47.3.87}
, number= {3}
, journal={Journal of The Korean Astronomical Society}
, publisher={The Korean Astronomical Society}
, author={SHIN, JIHYE
and
KIM, JUHAN
and
KIM, SUNGSOO S.
and
PARK, CHANGBOM}
, year={2014}
, month={Jun}