Calibrating an updated SPH scheme within GCD+

We adapt a modern scheme of smoothed particle hydrodynamics (SPH) to our tree N-body/SPH galactic chemodynamics code GCD+. The applied scheme includes imple- mentations of the artificial viscosity switch and artificial thermal conductivity pro- posed by Morris&Monaghan (1997), Rosswog&Price (2007) and Price (2008), to model discontinuities and Kelvin-Helmholtz instabilities more accurately. We first present hydrodynamics test simulations and contrast the results to runs undertaken without artificial viscosity switch or thermal conduction. In addition, we also explore the different levels of smoothing by adopting larger or smaller smoothing lengths, i.e. a larger or smaller number of neighbour particles, Nnb. We demonstrate that the new version of GCD+ is capable of modelling Kelvin-Helmholtz instabilities to a simi- lar level as the mesh code, Athena. From the Gresho vortex, point-like explosion and self-similar collapse tests, we conclude that setting the smoothing length to keep the number of neighbour particles as high as Nnb~58 is preferable to adopting smaller smoothing lengths. We present our optimised parameter sets from the hydrodynamics tests.


INTRODUCTION
Since it was introduced by Lucy (1977) and Gingold & Monaghan (1977), the smoothed particle hydrodynamics (SPH) methodology has become a regular tool for the numerical simulation of a wide range of astronomical phenomena.Hernquist & Katz (1989) were the first to suggest that the SPH approach would also prove invaluable in the simulation of galaxy formation and evolution.Since then, a number of SPH codes have been developed to simulate such systems, incorporating various physical processes ranging from radiative cooling to star formation and supernovae (SNe) feedback (e.g.Katz 1992;Navarro & White 1993;Katz et al. 1996;Steinmetz & Muller 1995;Mori et al. 1999;Carraro et al. 1998;Kawata 1999;Sommer-Larsen et al. 1999;Springel et al. 2001;Kobayashi 2004;Governato et al. 2004;Okamoto et al. 2008;Saitoh et al. 2008).
Parallel to the development of such particle-based codes, grid-or mesh-based approaches have been employed for modeling the formation and evolution of galaxies (e.g. ⋆ E-mail: dka@mssl.ucl.ac.ukCen & Ostriker 1992).Algorithmic enhancements to a fixed grid approach, such as adaptive mesh refinement (AMR), has led to a massive improvement in the capability of gridbased codes for simulations which require a large dynamic range, including those of galaxy formation (Teyssier 2002;Kravtsov 2003;Tasker & Bryan 2008;Gibson et al. 2008;Joung et al. 2008).Code comparisons between SPH and AMR (Frenk et al. 1999;Ascasibar et al. 2003;Voit et al. 2005;O'Shea et al. 2005;Gibson et al. 2008;Mitchell et al. 2008;Tasker et al. 2008) demonstrate that the competing approaches lead to generally consistent results.That said, comparing the results of non-radiative simulations of the formation of a galaxy cluster, Frenk et al. (1999) claim that SPH codes lead to lower entropy in the central region of the simulated cluster (see also Ascasibar et al. 2003;Voit et al. 2005;O'Shea et al. 2005;Dolag et al. 2005;Mitchell et al. 2008;Wadsley et al. 2008).They suggest that SPH may underestimate turbulence in the central region.Agertz et al. (2007) carried out a series of experiments in order to compare and contrast SPH and AMR in more of a "controlled" environment.They conclude that there is a "fundamental" discrepancy between these approaches, by demonstrating that SPH, at least in its conventional form, cannot capture Kelvin-Helmholtz instabilities (KHI) as accurately as an AMR approach.They further suggest that this discrepancy is not due to one of resolution, but is a fundamental attribute of the scheme itself (see also Imaeda & Inutsuka 2002;Okamoto et al. 2003).
Recently, Price (2008) described a new scheme to improve the conventional implementation of SPH and demonstrated the successful capture of KHI.In what follow, we apply this scheme to our original galactic chemodynamics code, GCD+ (Kawata & Gibson 2003a).GCD+ is a threedimensional tree N -body/SPH code that incorporates selfgravity, hydrodynamics, radiative cooling, star formation, supernova feedback, and metal enrichment.At its heart, the new scheme differentiates itself from the conventional approach via the manner by which diffusion of thermal energy is introduced; we adopt primarily the formalism described by Rosswog & Price (2007).We demonstrate that this new scheme does indeed advance the abilities of SPH codes in a suite of controlled tests (some of which have been described by Price 2008), in addition to leading to more realistic simulations of galaxy formation and evolution.To our knowledge, our work is the first demonstration of this new scheme within the context of cosmological simulations which follow both gas and collisionless dark matter (DM) dynamics simultaneously.We focus only upon non-radiative simulations in this paper, and will study cases including radiative cooling and star formation in a forthcoming work.
Section 2 describes briefly the implementation of this new scheme within GCD+, hereafter differentiating this version from the conventional one by the moniker "GCD+ MkII".Section 3 compares the performance of GCD+ MkII with the conventional version of GCD+ under a range of test conditions, summarising our results in Section 4.

GCD+ MKII: ADVANCING GALACTIC CHEMODYNAMICS
We now describe the specific modifications made to the galactic chemodynamics code GCD+, which themselves are patterned closely after the methodology described by Rosswog & Price (2007).As such, we only outline the final formulae adopted, and refer the interested reader to Rosswog & Price (2007) for their formal derivation.
The density of the i-th SPH particle is defined by where W is the SPH smoothing kernel seen in equation (1) of Kawata (1999), rij ≡ |xi − xj |, and hi is the smoothing length of the i-th particle.We note in passing that GCD+ MkII only takes into account the smoothing length of the i-th particle, hi, to derive the density, while the original version of GCD+ (Kawata 1999) used the pair-averaged smoothing length, hij = (hi + hj)/2.The smoothing length is determined by where we adopt η = 2.4.Note that in our definition of the kernel, our smoothing length corresponds to twice that used by Rosswog & Price (2007), who adopt η = 1.2.The solution of equation ( 2) is calculated iteratively until the relative change between two iterations is smaller than 10 −3 (see Price & Monaghan 2007, for more details).Euler's equation is written as The first term of equation ( 3) corresponds to the pressure gradient, where From equation (2), ∂hi/∂ρi = −hi/(3ρi).The second term of equation ( 3) corresponds to the artificial viscosity (AV), where The signal velocity vsig adopted is where cs,i is the sound velocity of the i-th particle.The amount of AV is controlled by a time-dependent parameter, where in order to suppress AV in pure shear flows.The viscous parameter α AV i (t) varies with time.Rosswog & Price (2007) suggested the following function to evolve this viscous parameter (see also Morris 1997): where α0 = 0.1 and τi = hi 0.2cs,i .
Rosswog & Price (2007) adopt the source term, and set the maximum of the viscous parameter to be α AV max = 2.We call the version of our code which adopts this particular source term "MkII+lowα", hereafter.As shown later, we found that this formulation led to insufficient viscosity for the capture of some strong shocks.As such, our fiducial model uses and employs α AV max = 4.We call this fiducial code "MkII+highα".The third term of equation ( 3) corresponds to the gravitational force, and employs the adaptive gravitational force softening suggested in Price & Monaghan (2007), where the softening length is matched to that of the smoothing length.The fourth term of equation ( 3) is the correction term for adaptive softening, where We apply a cubic splice softening, as suggested by Price & Monaghan (2007); the associated formulae for φ ′ and ∂φ/∂h can also be found in their paper.
The energy equation is written as where ui is the thermal energy of the i-th particle, and Qu,ij Otherwise, it is described by where The second term within the parentheses of equation ( 16) corresponds to the artificial thermal conductivity (AC) (Rosswog & Price 2007;Price 2008).The thermal conductivity parameter, α C , evolves following where the source term is and For collisionless N-body particles (such as DM), only the gravitational force term of equation ( 3) is taken into account, as follows: The adaptive softening length of the DM is calculated using equation ( 2) and for the DM we apply ηDM = 4, as the simulations under consideration generally have DM particles of higher mass.Note that for the DM we use the mass density, ρi,DM = Σj mj,DMW (rij, hi), in equation (2).Therefore, this higher ηDM results in a DM softening length which is larger than that for the gas when the number density of DM particles is comparable to that of the gas.We also set the minimum softening length ǫmin for both N-body and SPH particles to be ǫmin = (cǫmp(M ⊙ )) 1/3 pc, where mp is the mass of the particles.For the work described here, we employ cǫ = 3000, which results in similar softening to that we used in our previous work (Kawata & Gibson 2003a,b).
For SPH particles, we apply an individual timestep scheme to integrate equations ( 3) and ( 15), and for collisionless particles we use equation ( 20).We also employ the timestep limiter suggested by Saitoh & Makino (2008).The timestep for SPH particles is based upon dti = min (dtCFL,ij, dtDYN,i), where the Courant-Friedrich-Levy condition is calculated by where On the other hand, CCFL = 0.1 is adopted for MkII+αhigh, as we found a lower CCFL was required in a one-dimensional shock-tube test.The requirement that the force should not change significantly within one timestep is satisfied by We set CDYN = 0.2 for SPH particles.The timestep for collisionless particles is derived solely from equation ( 22), and we adopt CDYN = 0.2.

RESULTS
Having outlined the improvements made to GCD+, we now compare its performance with that of the original version.
As noted in Section 2, we explore the impact of two different forms for the AV source term: MkII+lowα (using equation 12) and MkII+highα (using equation 13).

One-dimensional shock-tube test
Not surprisingly, our first experiments involve a version of the classical one-dimensional shock-tube test.The initial conditions are set by assuming the simulation region spans from x = −0.5 to x = 0.5; the region for which x < 0 is set to (ρ1, P1, v1) = (1, 1, 0), using 960 particles, and the region for which x > 0 is set to (ρ2, P2, v2) = (0.125, 0.1, 0), using 120 particles, adopting a γ = 5/3 throughout.Figs. 1, 2, and 3 show the results for the original version of GCD+, MkII+lowα, and MkII+highα, respectively.One can see a clear jump in thermal energy and pressure at contact discontinuities in the original version, while in the MkII+lowα version, the contact discontinuity is resolved, and the pressure and thermal energy distribution is much smoother.It is also remarkable that the number of particles employed to resolve the shock front in this case is so low.Although less smooth compared with that of the MkII+lowα version, the MkII+highα implementation also resolves the contact discontinuity better    than that found using the original GCD+.However, because we assume a higher α AV in this case, a greater number of particles is required to resolve the shock front.A smooth pressure distribution at the contact discontinuity is the key to simulate accurately KHI (Price 2008); as such, it would appear that both the MkII+lowα and MkII+highα versions of GCD+ are promising tools for modeling KHI within an SPH framework.

Point-like explosion test
We next consider the Sedov-Taylor-type spherical explosion test.Following Springel & Hernquist (2002), we set a periodic boundary box with low-temperature and homogeneous density (ρ = 1), using 64 3 SPH particles.At t = 0, we deposit E = 1 energy on the central particle and simulate the evolution thereafter.The analytic solution can then be derived via the adoption of Sedov-Taylor self-similarity.Fig. 4 shows the result at t = 0.04 using the original version of GCD+, while the grey line represents the analytic solution; it is clear that the original version of the code was incapable of capturing a strong shock such as this.The MkII+lowα version show much better performance in this regard (Fig. 5), although even still, particles do "penetrate" the shock front.
Conversely, the MkII+highα version does capture successfully the associated strong shock, as seen in Fig. 6.As such strong shocks can be an important byproduct of the processes governing galaxy formation and evolution, we henceforth adopt the MkII+highα version as our fiducial code.Fig. 7 demonstrates that the energy conservation form of    velocity distribution (Fig. 9), with the shock front significantly more apparent.The MkII+highα version possesses even less scatter in the velocity distribution, due to the higher AV, although the "sharpness" of the shock front appears comparable to that seen with MkII+lowα.It would appear that in three-dimensional simulations, the high α AV scheme may not change the resolution of the shock significantly.et al. (2007) introduced a straightforward test which allows for a given code (particle-or grid-based) to be assessed in terms of its ability to resolve KHI (see also Price 2008).In this section, we demonstrate that our fiducial model ( MkII+highα) is now successful in resolving such instabilities.Following Agertz et al. (2007), we consider a periodic boundary box with x = −0.5, 0.5, y = −0.5, 0.5, and z = −1/64, 1/64.The region within |y| < 0.25 is set to be the high-density region (with ρ h = 9.6), while the rest is the low-density region (with ρ l = 1).Equal mass particles are adopted in both regions, and 544 (256) particles are used to cover the x-axis for the high-density (lowdensity) region.The two regions are in pressure equilibrium and we assume P h = P l = 2.5.The high-density (lowdensity) region has velocity vx = −0.5 (vx = 0.5).We also added sinusoidal perturbations on the vertical velocity, using vy(x) = δvy sin(λ2πx), setting δvy = 0.025 and λ = 1/6.These perturbations are only imposed to the particles within |y| = 0.225, 0.275.As before, we assume γ = 5/3; this condition leads to a timescale for KHI of τKHI = 0.57 Fig. 11 shows the density and pressure distribution at t = τKHI for the MkII+αhigh version; this can be compared directly with Fig. 13 of Agertz et al. (2007).The latter demonstrates that instabilities do not appear in their SPH simulations.On the other hand, simulations with the grid codes employed capture the development of KHI within the expected KHI timescale, leading Agertz et al. (2007) to argue that this is a fundamental problem for conventional SPH codes.Okamoto et al. ( 2003) also show in their Fig. 9 that    although SPH can develop weak KHI, they quickly disappear owing to numerical momentum transfer.Fig. 11 demonstrates that the new version of our SPH code is capable of capturing KHI, and leads to similar results to those found using the AMR codes in Agertz et al. (2007).We confirm that the new SPH scheme proposed by Rosswog & Price (2007) appears to solve this apparent fundamental problem of SPH (see also Price 2008, who demonstrates it with two dimensional test).

Blob test
In the previous section, we demonstrated that the new version of our SPH code -i.e.MkII+highα -is capable of capturing KHI in a simple shear flow simulation.Agertz et al. ( 2007) also introduced a so-called "blob test", a more appropriate analogue for simulating galaxy formation and evolution.Following Agertz et al. ( 2007) (Oscar Agertz, priv  comm), we set a high-density spherical "blob", with temperature, T cl = 10 6 K, density, ρ cl = 0.07281 M ⊙ kpc −3 , and radius 193 kpc, within a low-density field of temperature Text = 10 7 K and density ρext = 0.007281 M ⊙ kpc −3 (note: the two regions are in pressure equilibrium).We consider a periodic box with {Lx, Ly, Lz} = {1, 1, 3 Mpc}.The blob is set at {0, 0, −0.75 Mpc} and initially assumed to be static, but the low density field is assumed to have a velocity of vz = 1000 km s −1 .The particles with the same mass are distributed randomly so that the above condition is satisfied.
Fig. 12 shows the results at t ∼ 2τKHI for the simulations with the new MkII+highα version, adopting 10 6 particles within the blob.One can see that the blob is disrupted significantly due to the instability.We should note that disruption is reduced relative to the results seen with AMR codes (see Fig. 4 of Agertz et al. 2007); as such, it may be the case that our new code still underestimates the effects of KHI.Fig. 13 shows the same simulation, but with lower resolution (where the blob consists of 10 5 particles).Comparing with Fig. 12, the low-resolution simulation significantly underestimates the effects of KHI.Note that the blob test in Fig. 4 of Agertz et al. (2007) has the same resolution as our low resolution simulation.We also carried out a simulation with the same initial conditions (kindly provided by Oscar Agertz) as Agertz et al. (2007), however the results were comparable to their published SPH results.We conclude that even with our new version of GCD+, in order to accurately capture KHI, high-resolution is required in both the low-and high-density regions.One must be ever-vigilant when it comes to the appropriate resolution when simulating systems where KHI become important.
For reference, we also show the results employing high-(Fig.14) and low-(Fig.15) resolution simulations, but with the original version of GCD+.Surprisingly, the original version show significantly more disruption than the new version.However, as shown in Sections 3.2 and 3.3, the original version employs an artificial viscosity which is too small, which is the agent responsible for allowing the instability to become more developed.

Santa Barbara Cluster
Our final test is based upon the so-called Santa Barbara Cluster (SBC) simulation, a now-standard test for cosmological hydrodynamics codes (Frenk et al. 1999).The initial conditions correspond to a cosmological simulation in which a galaxy cluster forms at the centre of the final volume.We extract a spherical region with a radius of 32 h −1 Mpc which is expanding with the Hubble flow, and assume an open boundary.The gas and DM particle masses in the highest-resolution region are 8.67 × 10 8 and 7.80 × 10 9 M ⊙ , respectively.We run the simulation using both the original version of GCD+ and the new MkII+highα version, and compare the results.The minimum softening (and therefore smoothing) lengths for gas and DM particles are set to 13.8 and 28.6 kpc, respectively, for the MkII+highα version.Conversely, the original version of GCD+ uses fixed Plummer softening lengths of 11.6 and 23.0 kpc, respectively, for gas and DM particles.Frenk et al. (1999) demonstrate that the various codes employed in their analysis lead to generally consistent results.However, they also found a clear difference in the entropy profile between the particle-and grid-based codes, with the latter resulting in systematically higher entropy than the SPH codes, where entropy is defined as s = ln(T /ρ 2/3 ).We compare density, temperature, and entropy profiles at z = 0 in Fig. 16.Comparing with Figs. 12, 17, and 18 of Frenk et al. (1999), the original version of GCD+ shows consistent results with those from the SPH codes used by Frenk et al. (1999).It is very interesting that the MkII+highα version of GCD+ shows significantly different results from the original version of GCD+.The MkII+highα version results in a much lower density, higher temperature, and hence significantly higher entropy at the centre of the cluster.As a result, the profiles of the MkII+highα version of the simulation are similar to those for the grid codes used by Frenk et al. (1999) 1 .In particular, our results show a marked similarity to those of the AMR code labeled therein as "Bryan" (a precursor to what is now referred to as Enzo).
We also ran a simulation with the MkII+highα version of GCD+, but without AC, the results of which are shown as crosses in Fig. 16.Such a test leads to a result similar to that seen when using the original version of GCD+.We conclude that it is AC which ultimately leads to the difference seen in the profiles generated with the original version of GCD+ and the new MkII+highα version.
Figs. 17 and 18 show the density, temperature and line-of-sight velocity distributions, in the central 4 Mpc region of the cluster.We can see that the original version of GCD+ shows more small-scale perturbations, particularly in the temperature distribution.On the other hand, the MkII+highα version results in a much smoother distribution.Finally, we compare the temperature and density distributions for the gas particles within virial radius in Fig. 19.As expected, the MkII+highα version shows a much smoother distribution.Differences in the temperature and density distributions are seen in the simulations of Section 3.3.Therefore, we expect the same results even in the case of spherically-symmetric monolithic collapse -i.e., in the absence of hierarchical clustering.We conclude by noting that we have still lack a clear understanding as to why the new version of GCD+ results in entropy profiles which differ from conventional SPH codes.

SUMMARY
We implement new AV and AC schemes within an Nbody/SPH galaxy formation code, applying adaptive softening to both the N-body and SPH particles.In this paper, we study how these new schemes work within the context of non-radiative simulations, and focus on the effect of the combination of the new AV and AC.
We demonstrate that this combination helps to "smooth" the thermal energy at the contact discontinuity.Because of this improvement, the new code succeeds in capturing KHI in a manner that conventional SPH treatments do not.In essence, this confirms that the AV and AC scheme proposed by Rosswog & Price (2007) and Price (2008) remedies the fundamental problem of SPH outlined by Agertz et al. (2007).It is also shown that to properly simulate the KHI, resolution remains critical, even for our new implementation.Therefore, one must remain vigilant when it comes to resolution required in situations where KHI may become important.
We also find that to capture strong shocks, like that expected in supernova explosions for example, a relatively high AV parameter, α AV , is required, which may be higher than what was been used previously for conventional galaxy formation simulations.We have applied this higher α AV to the new code, demonstrating that the new code is still able to model KHI, while simultaneously capturing the strong shock.
It is also intriguing that our new code leads to similar results found when employing AMR codes to simulate the formation of a galaxy cluster within a cosmological context.We believe that the new AV and AC schemes significantly help to improve conventional SPH codes, and the resulting code should provide a means to simulate more accurately the formation and evolution of galaxies.
In a forthcoming paper, we will carry out more realistic simulations of galaxy formation and evolution, including radiative cooling, star formation, SNe feedback and chemical evolution, comparing and contrasting its behaviour with conventional SPH approaches.

Figure 1 .
Figure 1.Results of a one-dimensional shock-tube test with the original version of GCD+ at t = 0.2.The grey line represents the analytic solution.

Figure 2 .
Figure 2. As in Fig. 1, but now using the MkII+lowα version of GCD+.

Figure 3 .
Figure 3.As in Fig. 1, but now using the MkII+highα version of GCD+.

Figure 4 .
Figure 4. Radial density (upper) and pressure (lower) distributions at t = 0.04 in the point-like explosion test using the original version of GCD+.The grey line shows the analytic solution.

Figure 5 .
Figure 5.As in Fig. 4, but now using the MkII+lowα version of GCD+.

Figure 6 .
Figure 6.As in Fig. 4, but now using the MkII+highα version of GCD+.

Figure 7 .
Figure 7.Total energy as a function of time for the original version of GCD+ (dotted line) and the new MkII+highα version (solid line), for the point-like explosion test.

Figure 8 .
Figure 8. Normalised velocity (upper), density (middle), and pressure (lower) distribution at an arbitrary time, employing the Bertschinger (1985) test with the original version of GCD+.The grey line represents the analytic solution.

Figure 9 .
Figure 9.As in Fig. 8, but now using the MkII+lowα version of GCD+.

Figure 10 .
Figure 10.As in Fig. 8, but now using the MkII+highα version of GCD+.

Figure 11 .
Figure 11.Density (left) and pressure (right) distributions at t = t KH within the KHI test, using MkII+highα version of GCD+.

Figure 12 .
Figure 12.Density (left) and pressure (right) distribution at t = 2τ KHI for the blob test, using the MkII+highα version of GCD+.The initial blob consists of 10 6 particles.

Figure 13 .
Figure13.As in Fig.12, but for a lower-resolution simulation.The initial blob consists of 10 5 particles.

Figure 14 .
Figure 14.As in Fig. 12, but employing the original version of GCD+.

Figure 15 .
Figure 15.As in Fig. 13, but employing the original version of GCD+.

Figure 17 .
Figure 17.Gas density (left), temperature (middle), and line-of-sight velocity (right) distributions in an arbitrary slice through the centre of the simulated cluster, using the original version of GCD+.

Figure 19 .
Figure 19.Gas temperature vs. density map for the original version of GCD+ (left), the MkII+highα version (middle), and the MkII+highα version without AC (right).