Gas stripping and mixing in galaxy clusters: A numerical comparison study
The ambient hot intra-halo gas in clusters of galaxies is constantly fed and stirred by in-falling galaxies, a process that can be studied in detail with cosmological hydrodynamical simulations. However, different numerical methods yield discrepant predictions for crucial hydrodynamical processes, leading for example to different entropy profiles in clusters of galaxies. In particular, the widely used Lagrangian smoothed particle hydrodynamics (SPH) scheme is suspected to strongly damp fluid instabilities and turbulence, which are both crucial to establish the thermodynamic structure of clusters. In this study, we test to which extent our recently developed Voronoi particle hydrodynamics (VPH) scheme yields different results for the stripping of gas out of infalling galaxies, and for the bulk gas properties of cluster. We consider both the evolution of isolated galaxy models that are exposed to a stream of intra cluster medium or are dropped into cluster models, as well as non-radiative cosmological simulations of cluster formation. We also compare our particle-based method with results obtained with a fundamentally different discretisation approach as implemented in the moving-mesh code AREPO. We find that VPH leads to noticeably faster stripping of gas out of galaxies than SPH, in better agreement with the mesh-code than with SPH. We show that despite the fact that VPH in its present form is not as accurate as the moving mesh code in our investigated cases, its improved accuracy of gradient estimates makes VPH an attractive alternative to SPH.
keywords:methods: numerical – hydrodynamics
In the hierarchical structure formation model (White & Frenk, 1991), interactions or mergers of galaxies with other galaxies are common processes, and depending on the environment, a diverse set of physical processes can become important. For example, when a galaxy falls into a cluster of galaxies, gravity may affect its integrity through gravitational tidal forces. The closer the galaxy gets to the cluster the more it is exposed to the hot intra-cluster medium (ICM). This is experienced as a headwind by the galaxy that compresses and removes the gas of the galaxy (Gunn & Gott, 1972). This may lead to the formation of a bow shock in front of the galaxy if it supersonically ploughs through the ICM.
Whereas the gas at the front of the galaxy is primarily compressed by ram pressure in this situation, the sides are exposed to strong shear flows which may trigger fluid instabilities, leading to stripping and mixing of the gas with the ICM. This in turn crucially determines how rapidly star formation is turned off in an infalling galaxy, and how metals are mixed into the group or cluster ICM (see e.g. Larson et al., 1980; Balogh et al., 2000; Domainko et al., 2006). The hydrodynamical interaction processes between galaxies and groups/clusters may thus influence important observational effects, such as the density morphology relation (Dressler, 1980; van der Wel et al., 2010; Boselli & Gavazzi, 2006), or the existence of a tight red cluster galaxy sequence (Baldry et al., 2004; Cortese & Hughes, 2009). Another interesting aspect besides the stripping is the generation of turbulence by galaxies infalling into clusters. Numerical simulations (Dolag et al., 2005; Vazza et al., 2006; Iapichino & Niemeyer, 2008; Dolag et al., 2009) indicate that turbulence can be generated in the wake of an infalling galaxy. This interesting phenomenon may play an important role in determining the dynamics of the cluster gas and could in principle be detected through X-ray observations of the broadening of sharp metal lines (Sunyaev et al., 2003).
Modelling these processes reliably goes back to Gunn & Gott (1972) and is a considerable challenge, as direct hydrodynamic simulations of the relevant processes are very difficult to model due to the large dynamic range, and the uncertainties associated with the modelling of star formation and feedback processes. Furthermore, it is not clear whether the hydrodynamical techniques presently in use are strongly affected by systematic errors in this regime (e.g. Agertz et al., 2007). Indeed, previous studies of galaxy-cluster interactions with different numerical schemes disagreed strongly about how fast a galaxy loses its gas to the cluster. For example, a grid-based simulation of Quilis et al. (2000) predicted a much higher stripping rate than SPH-based work by Abadi et al. (1999).
In this paper, we concentrate on the numerical aspects of the stripping of a galaxy’s gas as it interacts with the ICM. Agertz et al. (2007) has recently shown that the popular smoothed particle hydrodynamics (SPH) technique has problems to properly account for fluid instabilities, prompting significant concerns about a possible unphysical suppression of stripping processes (see also Price, 2008; Wadsley et al., 2008; Read et al., 2010; Cha et al., 2010; Valcke et al., 2010; Junk et al., 2010; Abel, 2011). In a recent study (Heß & Springel, 2010) we have therefore proposed a new ‘Voronoi particle hydrodynamics’ (VPH) method that improves on the widely used SPH technique (Lucy, 1977; Gingold & Monaghan, 1977; Larson, 1978) in several respects. By employing a Voronoi tessellation for the local density estimate, a consistent decomposition of the simulation volume is achieved, contact discontinuities can be resolved much more sharply, and a ‘surface tension’ effect across them is avoided. Our preliminary tests of VPH based on the ‘blob-test’ of Agertz et al. (2007) already suggested that stripping is more efficient in VPH compared with SPH. Here we shall investigate this in more detail using more realistic set-ups that mimic galaxy evolution processes. In addition to the two particle-based Lagrangian schemes for hydrodynamics, SPH and VPH, we will also carry out comparison simulations with the moving-mesh code AREPO (Springel, 2010a). Whereas AREPO uses a Voronoi tessellation as well, this code employs an entirely different methodology for fluid dynamics, based on a finite volume Godunov scheme that calculates hydrodynamical fluxes with a Riemann solver across mesh boundaries. The comparison of this diverse set of three numerical methods is useful to understand and quantify the systematic uncertainties of the different methods.
In our tests simulations, we ideally want to investigate as realistic conditions as possible, including a full treatment of gravity and dark matter, as well as star formation and cooling. This in principle calls for simulations in which a well-resolved galaxy model is dropped into an equally well represented cluster of galaxies. Because the substantial computational cost of such a set-up severely limits the resolution that can be achieved, we base part of our study on a number of more idealised simulations, for example by placing galaxy models into a ‘wind tunnel’ that mimics the impinging flow of gas onto a galaxy in orbit in a cluster. Finally, we also carry out cosmological simulations of the formation of the ‘Santa Barbara cluster’ (Frenk et al., 1999), primarily to study how well VPH performs in non-radiative cosmological simulations of cluster formation compared to the other techniques. Already in past studies, the Santa Barbara cluster comparison project has led to important insights into how numerical effects impact the thermodynamic structure of simulated clusters, and for this test problem, there is already a considerable body of results in the literature (e.g. Wadsley et al., 2004; Springel, 2005; Thacker & Couchman, 2006).
This paper is structured as follows. In Section 2, we introduce the numerical methods we will use in our numerical comparison study. In Section 3, we check how well the different numerical schemes represent the evolution of an isolated galaxy, and whether there are already significant differences at this level. We then consider in Section 4 wind-tunnel experiments in which we expose the isolated galaxy models to a supersonic wind, allowing a detailed examination of the stripping process. In Section 5, we follow up on these experiments by studying the behaviour of a galaxy model that falls into an isolated spherical cluster, comparing the VPH results with those obtained with SPH and AREPO. Finally, in Section 6, we investigate the performance of VPH in cosmological simulations of cluster formation, using re-simulations of rich clusters identified in the Millennium Simulation as well as the Santa Barbara cluster initial conditions. We give a discussion of our results in the context of a vortex test problem in Section 7, combined with a summary of our conclusions.
We begin by introducing the different codes and computational methods we use in this study. In the interests of brevity, we shall only describe the most important characteristics of the codes; full details can be found elsewhere (Springel, 2005; Heß & Springel, 2010; Springel, 2010a).
2.1 Smoothed particle hydrodynamics
Smoothed particle hydrodynamics (SPH) is a particle based approach to fluid dynamics that has seen widespread use in astronomy since it was first introduced more than three decades ago (Lucy, 1977; Gingold & Monaghan, 1977). SPH discretises the mass of the fluid in terms of particles, whose dynamics is governed by a fluid Lagrangian of the form
Here describes the thermal energy per unit mass, which for an ideal gas depends only on density and specific entropy . An adaptive, spherically symmetric smoothing kernel is employed to estimate the density based on the spatial distribution of an approximately fixed number of nearest neighbours. (As a standard value within this work we set the number of neighbours to .) The Lagrangian in Eq. (1) then uniquely determines the equations of motion in a form that simultaneously conserves energy and entropy (Springel & Hernquist, 2002; Price & Monaghan, 2007; Rosswog, 2009; Springel, 2010b). However, an artificial viscosity needs to be added to allow for dissipative shock capturing and to produce a damping of post-shock oscillations. If not explicitly denoted otherwise we use the viscosity parametrisation of the GADGET code with a viscosity strength of , and we apply the viscosity limiter in the presence of shear introduced by Balsara (1995).
Among the primary advantages of SPH are its very good conservation properties, the manifest Galilean invariance at the discretised level of the equations, and the absence of preferred spatial directions. However, it has recently become clear that contact discontinuities, in particular, are challenging for the method. For example, they are associated with a spurious surface tension effect in SPH (Springel, 2010b), and fluid instabilities such as the Kelvin-Helmholtz instability are suppressed across them (Agertz et al., 2007).
2.2 Voronoi particle hydrodynamics
Because of the accuracy problems of SPH in certain situations, we have recently developed another particle-based hydrodynamic method that differs from SPH in important ways, while still sharing a number of similarities. In this ‘Voronoi Particle Hydrodynamics’ (VPH) method (Heß & Springel, 2010), the fluid is also discretised in terms of mass elements , and the same fluid Lagrangian as in Eq. (1) is used. Also identically to SPH we assume a polytropic equation of state with an adiabatic coefficient in all our tests. However, the kernel estimation technique of SPH is replaced by a very different approach to calculate the gas density. This is done by constructing a Voronoi tessellation for the point set, in which to each of the particles the volume of all the space which is closer to this point than to any other point is assigned. The density can then be obtained as the particle mass divided by the volume of the associated Voronoi cell. We note that such a tessellation technique allows a close to optimum exploitation of the density information contained in the particle distribution (Pelupessy et al., 2003), and in particular, very sharp discontinuities can be resolved over the distance of a single cell.
With the density computed, and given specified entropies for every particle, we can use the ideal gas Lagrangian to derive unique equations of motion for VPH. The key quantity that is needed in doing this is the derivative of the particle volume with respect to the coordinates of the surrounding particles, but thanks to the mathematical properties of the Voronoi tessellation, the corresponding geometric factors can be computed relatively easily (Serrano & Español, 2001; Heß & Springel, 2010) and have a clear geometric meaning. For example, the primary force between two VPH particles is their pressure difference times the area of their common tessellation face. Because the dynamics in VPH is derived from the Lagrangian given in Eqn. (1), energy, momentum and entropy are conserved exactly, just like in SPH. However, there is no surface tension effect, and the VPH particle volumina add up to the total simulated volume exactly, something that is not the case in SPH in general. One can hence hope that VPH provides an interesting improvement over SPH. Indeed, in Heß & Springel (2010) we have shown that this is the case for several test problems, such as the growth of Kelvin-Helmholtz instabilities at contact discontinuities with a large density jump. It is one of the primary goals of this paper to see whether these improvements are also relevant in cosmological structure formation problems.
In VPH, we need to rely on an artificial viscosity to treat shocks, similarly to SPH. We here follow the standard SPH approach by Gingold & Monaghan (1977) and Balsara (1995) in introducing an extra friction force that reduces the kinetic energy and transforms it into heat in regions of rapid compression. Due to its similarities to its SPH counterpart we usually use a value of for our artificial viscosity in VPH as well, but we also vary explicitly to demonstrate its implications. In addition, in VPH we can also add some weak, non-dissipative “shape forces” that are designed to maintain a reasonably regular mesh geometry, which tends to improve the overall accuracy of the VPH dynamics. As discussed extensively in Heß & Springel (2010), there are different possibilities to introduce such extra forces. We derive them by adding extra terms to the Lagrangian which penalise highly distorted mesh cells. For the shape correction scheme we usually choose and . However, the optimum choice for the strength of these terms is unfortunately not (yet) clear. For this reason, we have typically carried out all our tests with several different settings for the strength of the artificial viscosity and the shape forces, finding in most cases that our results are quite insensitive to the detailed choices over a broad range of parameter settings.
Both the SPH and VPH simulations we study here have have been carried out with GADGET-3, an improved version of the publicly available cosmological code GADGET-2 (last described in Springel, 2005). The code is parallelized for distributed memory machines and uses a hierarchical tree algorithm for calculation the self-gravity of the gas. A collision-less stellar or dark matter component can be optionally included as well.
2.3 Hydrodynamical moving-mesh code
Finally, the third numerical technique for hydrodynamics we examine is implemented in the moving-mesh code AREPO (Springel, 2010a). It uses an entirely different approach in which the volume is discretised. The hydrodynamical equations are solved in terms of a Godunov scheme where at each mesh interface a Riemann problem is calculated and the resulting fluxes are used for an exchange of conserved fluid quantities between the cells. For higher spatial accuracy, this is combined with a piece-wise linear reconstruction of the gas distribution and a second-order time integration scheme, yielding overall second-order integration accuracy in smooth parts of the flow.
The particular mesh used by AREPO is given by the Voronoi tessellation of a set of mesh-generating points. The use of this particular type of mesh allows the method to be formulated such that the mesh can move along with the gas, continuously adjusting its topology to the fluid motion. The method exploits the property of a Voronoi tessellation that the mesh changes continuously when the generating points are moved, without showing problematic mesh-tangling effects. If the points are moved with the local flow velocity (which is the standard way to operate the code), the method then has a similar Lagrangian character as SPH/VPH. However, unlike in SPH and VPH, the masses of AREPO’s resolution elements (the cells) are not required to remain strictly constant. In fact, if the mesh-generating points are being fixed on a static Eulerian grid, then AREPO is equivalent to the MUSCL-Hancock scheme, a popular method for second-order Eulerian gas dynamics on structured grids. AREPO should hence be viewed as being closely related to ordinary Eulerian mesh-based hydrodynamics, except that its ability to move the mesh along with the flow considerably reduces advection errors.
3 Isolated galaxies and their evolution
Before we analyse the interaction of individual galaxies within a cluster environment, we here study galaxy models in isolation to see how well the different numerical schemes treat them and whether our initial conditions are sufficiently stable. This is also intended to identify possible subtle systematic differences between the numerical schemes that may be difficult to detect in more complicated set-ups. We hence use identical initial conditions for the three codes such that differences introduced by varying start-up procedure are minimised. We we will see however that even for identical particle positions, masses, and temperatures, there can be variations in density (and hence specific entropy since we start with a prescribed temperature) directly after the start of a simulation, as a result of the different density estimation methods. We are in fact especially interested in the question whether these small differences have any bearing on the further evolution of isolated galaxy models.
To set-up compound galaxies in isolation, we construct dark matter halo models and stellar disks in dynamical equilibrium through an approximate solution of the Jeans equations (Hernquist, 1993; Springel & White, 1999; Springel et al., 2005). The primary assumption made is that the local velocity dispersion can be approximated reasonably well with a triaxial Gaussian. For the models considered here, a dark matter halo with a Hernquist profile is used, but alternatively a truncated NFW profile could be used as well, with no difference to our conclusions. The baryonic matter is represented with a central stellar bulge, an exponential stellar disk, and an exponential gaseous disk. In addition, we adopt a gas distribution of very low density in the halo, taken to be in hydrostatic equilibrium in the total potential such that this gas is close to the virial temperature. This component helps to avoid problems in VPH and AREPO due to “empty space” – unlike SPH, these codes can not easily deal with vacuum boundary conditions as they tessellate space. By having the low density background gas in the halo, which we extend to a box region around the galaxy, we can avoid this problem.
Following the formalism and nomenclature of Mo et al. (1998) and Springel et al. (2005), we have chosen a virial velocity of , yielding a total mass of for the system, of which a fraction is assumed to be in the disk, a fraction of the mass is in the bulge and the rest is distributed in a halo with a concentration of and a spin parameter of . The Hubble constant is set to . Since we want to concentrate on the gas dynamics, we make the disk very gas-rich by giving it a gas fraction of . The gas component in the halo is given a mass fraction , considerably less than the amount of gas in the disk. This additional gas in the halo is represented with particles/cells of mass equal to the gas particles in the disk, and its spatial distribution is clipped outside a box of side-length centred on the galaxy, allowing us to impose periodic boundary conditions for the gas there.
We include radiative cooling (for a primordial composition), and model star formation and supernova feedback with the simple multi-phase model for the ISM of Springel & Hernquist (2003). For the supernova feedback model we adopt an elevated effective temperature of and an evaporation efficiency of . This is done in order to shift the density threshold for the onset of star-formation to a 10 times higher values than in the default model of Springel & Hernquist (2003), while keeping the star-formation timescale at and maintaining a good agreement with the Kennicutt (1998) relation. With this change, the star formation threshold corresponds to a hydrogen number density of . Such an elevated star formation density threshold has recently been advocated to facilitate the formation of more realistic spiral galaxies (Guedes et al., 2011), and it will sharpen any potential issues brought about by the density contrast between disk and gaseous corona, which is what we are specifically interested in in this work. If not denoted otherwise we employ this model throughout this paper whenever star formation and/or cooling is required.
In our default intermediate resolution (hereafter called ‘R4’), the gas in the disk is represented by particles and an equal number of collision-less particles for the stellar disk. The bulge is represented by particles, and the dark matter halo by heavier particles that are given a larger gravitational softening length of compared to for the softening length of the baryonic particles. The gas in the halo out to the boundaries of the simulation box is represented with an additional particles. Besides this default resolution we have also simulated realisations of the same galaxy model both at lower and higher resolutions, in order to get a good sense of numerical resolution effects. To this end we lowered the resolution in three steps by factors of 2 in all particle components, and we increased it in three steps by factors of 2 as well. We hence arrived at 7 different resolutions R1-R7, with , , , , , , and resolution elements in the gaseous disk, respectively (see Table 1). The gravitational softening was varied in proportion to the cube root of the mass resolution, , in these simulations. All our simulations with SPH, VPH and the AREPO code were started from the same identical initial conditions files in order to yield a direct comparison that also allows an identification of different start-up systematics. For the simulations with AREPO, we have activated on-the-fly refinements and derefinements similar to Vogelsberger et al. (2011) in order to keep the mass resolution always close to the initial value.
|run name||gas particles||disk + bulge particles||softening gas/stars||DM particles||softening DM|
3.1 Gas density maps and structure of the disk
As we have remarked above, identical particle positions and masses do not necessarily imply equal densities in our different schemes. To examine this point, we consider in Figure 1 each method’s estimate of the initial gas density profile in the -direction, perpendicular to the disc. The solid curves represent averages of the densities estimated for the individual particles/cells as a function of distance from the disk plane, while the dashed line is the actual average mass density distribution obtained by binning the particles directly. Since both VPH and AREPO use a Voronoi tessellation for the density estimate, they agree on the densities at the initial time, but the kernel estimates of SPH (based on smoothing neighbours) show an interesting difference. In particular, the particles close to the dense disk gas (for ) show an elevated density estimate because they “see” some of the dense particles within their smoothing radius. This effect is not unexpected given that sharp discontinuities will always be smoothed out to some extent in SPH, by construction. The density estimate of VPH has a somewhat better resolving power, but here a different systematic effect becomes noticeable. Particles in the surface layer of the dense disk sometimes extend their volume quite far into the region above the disk until they “see” a particle of the background corona, thereby underestimating the density in this region. Finally, there is an interesting systematic difference between SPH and VPH at the very centre of the disk, where SPH shows a small positive bias in the density estimate. This comes about in part due to the Poisson sampling present in the initial conditions, which in SPH causes larger positive excursions in the density estimates. For a system in pressure equilibrium, this effect is weaker, but it not necessarily vanishes because also here the SPH density may be biased and the sum of the volumes associated with each particle is not guaranteed to add up to the total volume.
In Figure 2, we analyse the vertical structure of the gas after the isolated galaxy has been evolved for a time of . In the top panel of the figure, we again compare measurements of the mean density estimated for particles/cells as a function a distance above the disk mid-plane with the actual mean mass profile present in the simulation. We give results for all three numerical methods (solid lines), and since the actual gas distribution can have evolved differently in the three methods, this is given individually as well (dashed lines). We see that the two measures of the density stratification track each other well both in VPH and AREPO. However, the moving-mesh code shows higher density just above the disk compared with VPH. If anything, SPH shows still higher densities in this region, although its actual mass distribution is very similar to that of VPH. This situation arises due to a significant difference SPH exhibits in the estimated density versus the real mass distribution in the regime of the ‘shoulder’ of the disk’s density distribution. Also, there is actually a ‘gap’ in the run of SPH’s density estimates, arising simply because we find no SPH tracer particles for measuring the density there. This gap becomes more explicit in the bottom panel of Figure 2, where we show counts of fluid tracer particles in logarithmic bins in the same region above the disk. The sampling gap at in SPH is evident in this figure, but it does not occur in VPH. We note that Agertz et al. (2007) pointed out that such a gap is seen generically at strong density jumps in SPH, and that it may be a primary cause for the suppression of fluid instabilities across such contact discontinuities. Recently, a number of authors (Price, 2008; Read & Hayfield, 2011) have suggested that SPH should be outfitted with additional dissipation schemes that smooth out such structures to improve on this behaviour.
In Figure 3, we show face-on maps of the galaxy’s projected gas density distribution, after a period of of evolution. It is reassuring but probably not too surprising that the overall morphology is extremely similar, as the gravity from the dark matter halo and the stellar disk are primary driving forces of the disc dynamics. However, upon close inspection, one can identify some interesting systematic differences between the techniques that manifest themselves in the gas density maps. Compared with SPH, the higher effective spatial resolution of VPH produces crispier but marginally noisier looking features such as spiral arms. Perhaps the most prominent difference is however the lower density of the gas found in VPH compared to SPH in some regions between spiral arms. The gas distribution of the mesh-based AREPO looks a bit sharper and less noisy than both particle-based schemes.
Nevertheless, there is reassuring similarity of the gas surface density distribution between the different techniques, an impression that is corroborated by maps of the projected stellar mass density in the disks, which we show in Figure 4. Here we show only the newly formed stars in the simulations, so that any structure in this component directly reflects the underlying gas dynamics. We find good agreement for these collisionless particles, which can also be interpreted as evidence that the gravitational dynamics of collisionless particles is followed with equal quality in all three codes. This verifies that differences in the evolution of our galaxy models are expected to arise only from the different treatment of the hydrodynamics, and not from differences in the way the collisionless dynamics of stars and dark matter is followed.
Despite the systematics we identified in the density estimates and in the projected gas distributions, we overall find that all three codes evolve the galaxy model in a similar and stable fashion. This is further confirmed by radial surface density profiles of the gas and the newly formed stars shown in Figure 5. All three methods are clearly able to retain the initial structure of the galaxy model to a similar degree, yielding only very minor differences after some time.
3.2 Star formation rate evolution
Despite the good agreement between the different methods noted above, some minor but interesting discrepancies are revealed when looking in detail at the time evolution of the star formation rate at different resolutions, as done in Figure 6. In the top panel, we compare the results for SPH, VPH and AREPO at our default resolution of R4, which lies in the middle of our extended set of models used in our resolution study. There is in general good agreement between the runs at late times, where SPH and AREPO agree extremely well and VPH lies at most a few per cent lower. However, at earlier times some large differences are noticeable. First of all, at , the star formation rate for SPH is about per cent higher than that of VPH and AREPO, a difference entirely caused by different density estimation methods (see below). But after an initial transient, when the disk settles from the approximate equilibrium realised in the initial conditions to a proper self-consistent equilibrium, somewhat larger differences are present for an extended period of time. VPH tends to show a slightly smaller star formation rate than SPH, with AREPO lying somewhere in the middle.
Some clues to the origin of these differences are provided by the lower three panels of Figure 6, which give the evolution of the SFR for all seven numerical resolutions we considered for our isolated galaxy model, and for each numerical method. We see that all three techniques show some residual drift in the SFR at low resolution, but they tend to converge reasonably well towards higher resolution. Given the fact that the mass resolution is varied by a factor of 64 here, a scatter in the star formation rate of order per cent can probably be considered satisfactory. Nevertheless, there are clearly interesting differences in the strength of the resolution-dependence of the numerical techniques. In particular, VPH shows a somewhat stronger reduction of the SFR at the lowest resolution compared to SPH, whereas the latter is remarkably resilient to resolution changes. Perhaps counter-intuitively, this arises despite the relatively strong biases in SPH’s density estimate at the interface between disk and corona.
Another look at this issue is provided by Figure 7, in which we show how much gas mass is estimated by the codes to lie at a certain density value. In the top panel of Fig. 7, this is given for the initial time. The fact that SPH systematically estimates more gas to lie at a density value above the star formation threshold than the Voronoi-tessellation schemes (despite an identical point distribution in this case) explains the offset in the SFRs at . However, also at later times, as seen in the bottom panel of Fig. 7, such a systematic difference persists. The variations in the star formation rates calculated by the different codes appear hence to be primarily driven by the way the sharp edge of the star-forming disk is represented. These differences are however quite minor overall and can be largely ignored for our subsequent investigations.
4 Galaxy in a wind tunnel
Simulating the interaction of a galaxy with the sparse intra-cluster gas of a galaxy cluster at high-resolution is computationally challenging (e.g. Abadi et al., 1999; Tittley et al., 2001; Schulz & Struck, 2001; Roediger & Brüggen, 2006, 2007, 2008), especially with the particle-based methods SPH and VPH, which do not easily lend themselves to adaptive refinement techniques. In fact, they basically require particles of equal or very similar mass to work well (Ott & Schnetter, 2003). Hence, if one simply wants to let a galaxy model fall into a massive cluster, the latter needs to be represented with very high particle number (because it is so massive), such that one ends up spending only a tiny fraction of the computational effort onto the galaxy and the surrounding gas, where the actual interaction takes place. A further complication is that some of the stripped gas mass may be distributed over a large region across the cluster, corroborating the need to simultaneously account for the whole cluster at high resolution.
We will nevertheless consider such simulations later on in Section 5 to check the validity of our approach, but here we first investigate an alternative, more idealised set-up, which makes it much easier to reach an adequate resolution. In this approach, we effectively put a galaxy model into a ‘wind tunnel’, i.e. a rectangular box in which we let gas stream onto the galaxy with prescribed density, velocity and temperature, matched to what we expect during a cluster passage. The much smaller volume that needs to be simulated around the galaxy in this situation drastically lowers the computational expense, and even more importantly, the simplification brought about by such a controlled set-up makes it much easier to distinguish between different numerical effects.
4.1 Initial conditions
For definiteness, we put our model galaxy into a rectangular box of side-length . At one side of the box, gas is constantly injected, which we realise in the case of SPH and VPH by creating new particles in a suitable fashion, mimicking an infinitely extended grid of particles that moves into the simulation domain. On the opposite side of the box, we effectively implement outflow boundary conditions by removing particles once they start moving out of the box. Since we here focus on supersonic winds, the removal of particles at the outflow side does not lead to the propagation of perturbations upstream to the inflowing gas. In the case of the AREPO code, the inflow and outflow regions are realised in an equivalent, yet technically different fashion. We here use the on-the-fly refinement and derefinement features of the AREPO code to produce new mesh-generating points in the inflow region, and to remove them near the outflow side. The net result is identical to the SPH and VPH cases, i.e. the galaxy at the centre of the box is hit by a wind of particles/cells with prescribed density and velocity. This impinging wind has the same gas mass resolution as our target galaxy model.
For all the boundaries of the box that form the enclosing sides of our wind tunnel, we adopt periodic boundary conditions for the gas as far as the hydrodynamics is concerned. However, periodic boundary conditions are not imposed for self-gravity, which we fully include (note that our galaxy is held together by gravity). Since the mass of the gas in the wind is quite small, most of it is unbound to the galaxy, and our restricted treatment of self-gravity to the region of the box is still a good approximation. We note that a number of similar wind tunnel experiments have previously been carried out (e.g. Agertz et al., 2007; Iapichino & Niemeyer, 2008), but without a treatment of self gravity. While our setup is not capable of accounting for the tidal forces that arise when an extended galaxy travels through the gravitational potential of a cluster, these forces can be expected to be sub-dominant compared to the hydrodynamical forces in the outer parts of a cluster, which is the region we are most interested in. However, gravitational forces between dark matter, stars and gas inside the galaxy are very important, and our approach is the first that allows a detailed study of the response of the gravitationally coupled disk-halo system to the ram-pressure of the impinging wind. In particular, this allows a realistic measurement of the acceleration of the whole galaxy due to the ram-pressure force exerted by the wind.
In our default set-up, we have chosen a constant density of and a velocity of for the wind, letting it flow onto the galaxy in a face-on orientation. A density and velocity of this magnitude would be typical close to apocentre for a galaxy that has fallen into a large cluster. Of course, in reality the strength of the wind will be a function of the orbital phase and of parameters such as cluster size and angular momentum of the galaxy orbit. Changing these values should however not change our results at a qualitative level. For the galaxy model studied in our simulations, we have adopted structural properties identical to those described in Section 3. We place the galaxy suddenly into the supersonic wind, refraining from any attempt to impose a more gradual start-up in which the wind would somehow be slowly turned on.
Hydrodynamically, the galaxy represents an obstacle in the supersonic wind, causing a bow shock in the upstream region ahead of the galaxy. Inside the bow shock and ahead of the front side of the galaxy, the compressed wind is slowed down and exerts significant ram-pressure onto the galaxy (Gunn & Gott, 1972), whereas around the sides, a region of strong shear flow is produced as the wind streams around the galaxy. In this region, the formation of Kelvin-Helmholtz instabilities is expected, which accelerate the stripping of gas from the galaxy’s disk and mix it with the downstream flow.
Agertz et al. (2007) has previously investigated the simpler case of a spherical, non-self-gravitating gaseous obstacle (the ‘blob-test’) and found substantial discrepancies between different hydrodynamical methods for the rate of gas stripping and the time until eventual complete disruption of the gas blob. Here, we are particularly interested to see whether such differences also occur between SPH, VPH and AREPO when the more realistic disk-wind setup we are investigating is considered.
In Figure 8, we show projected gas density maps of the galaxy in the wind tunnel in an edge-on orientation, with time evolving from bottom to top. While the morphology of the stripping process looks similar at early times in all three codes considered, significant differences develop over time. Both of the particle-based methods loose small clumps of gas and exhibit string-like features of quite dense gas. This effect is particularly strong in SPH. In contrast, the AREPO disk stays comparatively coherent, with all of the stripped gas moving to lower densities quickly. Hence there remains no dense debris in the downstream flow. At the latest time shown in the figure (top row), the residual gas disks of AREPO and VPH are visually significantly smaller than in SPH, suggesting that more gas has been stripped.
This impression is born out by a quantitative measurement of the gas mass that is still in the ISM phase of the galaxy model. We consider a gas particle (or cell) to be part of the ISM of the galaxy when it is sufficiently dense and cold. Specifically, we require the conditions
to be met, where and are the density and temperature of the wind. Furthermore, we impose an additional condition on the maximum allowed separation of a particle from the centre-of-mass of all gas that fulfils the criterion of Eq. (2). Only for the gas particle is considered to be part of the ISM of the primary galaxy. This condition is introduced in order to discard cold clumps of gas that have been stripped out of the galaxy but have still stayed cold and dense. Indeed, it turns out that this secondary criterion is quite important in our SPH simulations, where the galactic disk tends to shed dense clumps under the action of the ram pressure. These fragments can still fulfil Eqn. (2) despite being stripped, but as they separate from the main galaxy, they eventually violate the distance condition and are then counted as stripped gas.
In Figure 9 we show the ratio of the remaining ISM gas mass to the initial mass a function of time, comparing wind tunnel simulations carried out with the three different numerical techniques. Clearly, the mass loss in SPH is smallest overall, with about 10% of the initial mass remaining after , whereas at this time the ISM of the VPH and AREPO runs is almost stripped completely. Note that this measure of the mass loss counts most of the dense gas blobs that stay coherent in SPH as lost gas. If we were to measure the amount of “non-dispersed” gas instead, the difference would be considerably larger because the gas that is stripped in SPH from the disk is largely not dispersed, unlike in the other two approaches (see below). The relative similarity of the overall shape of the gas mass evolution in Figure 9 is therefore not simply a time delay in SPH relative to the other methods, even though at intermediate times the mass loss rates are similar.
It is interesting to note that the stripped dense gas in the SPH simulation, and to a lesser extent in the VPH run, can continue to form stars at some level. This is particularly evident in Figure 10, which shows maps of the star formation rate density in an edge-on projection, at time . This highlights that in the case of the SPH simulations the stripped gas clumps are not only very concentrated, they also continue to form stars. This phenomenon is still present in VPH, albeit at a much weaker level and with decreasing star formation rate over distance. In contrast, AREPO does not show such a behaviour.
In SPH, there are at least two effects that can help to explain the enhanced star formation in the stripped clumps. First, it is known that spurious surface tension forces in SPH exist that will slightly compress the gas of the clumps and keep them coherent. Second, there is the issue of a potential artificial suppression of gas mixing in the turbulent wake behind the galaxy. Both SPH and VPH, by construction, prevent that the low entropy gas of the ICM is mixed at the particle level with the high entropy gas of the impacting wind. This imposed Lagrangian behaviour ignores any small-scale mixing processes that could happen on scales below the spatial resolution limit. In contrast, the mesh-based approach of AREPO implicitly allows for such processes. Here, the mixing manifests itself through mass exchanges between the cells, giving rise to a much more diffuse wake behind the galaxy (left column in Figure 8) where no gas remains that is still dense enough to support star formation. In the mesh-based approach, one may have in fact an opposite problem compared to the particle schemes. Here the need to advect the fluid over mesh interfaces can easily lead to excessive numerical mixing. However, the moving-mesh technique of AREPO should reduce such errors significantly compared to more traditional Eulerian methods.
Another view on the different behaviour with respect to stripping is provided by Figure 11, which shows the time evolution of the sum of the star forming gas mass (which is just the total gas mass above the star formation threshold) plus the stellar mass formed since the beginning of the simulation. Interestingly, this quantity is nearly constant in time in the case of SPH, showing that essentially none of the initially dense and star forming gas is dragged to low enough density to stop star formation. This is consistent with our earlier findings, but also makes it quantitatively evident that essentially no mixing of this low entropy gas with higher entropy material occurs. In contrast, the sum of star forming gas mass plus stellar mass declines significantly in AREPO, by more than a factor of 2 over the timescale of that is shown here. The results for VPH are intermediate, showing some evidence for gas moving to higher entropy or lower densities during the stripping process and thus reducing the amount of total gas available for star formation.
Interestingly, the discrepancies in the stripping rates also lead to different accelerations of the whole galaxies, because they experience unequal ram pressure forces as a result of different effective areas of the remaining gas disks. Since the gas in the disk is gravitationally coupled to the stellar disk and the dark matter halo, it is important to note that the wind not only strips the gas component of the galaxy, it also accelerates the whole galaxy. In Figure 12, we show maps of the projected stellar densities of the stellar disks in SPH and VPH at an identical time, corresponding to . Clearly, the whole SPH galaxy has been pushed more into the downstream direction, demonstrating that it has experienced a higher effective ram pressure, as a consequence of the reduced stripping rate of its disk. In real galaxy clusters, this will in turn lead to a larger dynamical friction force and a somewhat faster decay of the orbit of the galaxy. An additional effect is that in the SPH case we see small features emanating in the downstream direction from the stellar disk. Figure 12 shows two small humps in the stellar disk, as well as a faint stellar blob at some distance further downstream. These are in fact stars of the original stellar disk that have been gravitationally pulled out of the disk by dense gas. Both effects are absent in the equivalent VPH and AREPO calculations, emphasising again how much more (artificial) “coherence” the dense gas phase exhibits in the case of SPH, where comparatively massive gas clumps are removed from the disk that stay coherent afterwards.
4.3 Dependence on resolution and artificial viscosity
|Code||resolutions||soft. baryons||softening DM|
|VPH||4k, 20k, 80k, 320k|
|SPH||4k, 20k, 80k, 320k|
|AREPO||4k, 20k, 80k, 320k|
In order to examine the quantitative robustness of our wind tunnel results with respect to numerical resolution, we repeated our simulations with SPH, VPH and AREPO using both a lower and a higher resolution than in our default runs discussed thus far, where the galaxy is resolved with gas particles in the initial gas disk. With respect to our default resolution, our low resolution simulations have 4 times fewer particles in each component, whereas the high resolution runs have 4 times mores particles, as denoted in Table 2. We here use gravitational softening lengths independent of resolution.
In Figure 13, we show results of our resolution study for the gas stripping out of the ISM as a function of time, and for the sum of remaining star forming gas and newly formed stellar mass. We find that the general trends remain very similar at all three resolution levels, i.e., gas is stripped the fastest in AREPO and slowest in SPH, with VPH taking on an intermediate role. Interestingly, in the two particle codes, the stripping proceeds more slowly when the resolution is poorer, whereas in AREPO the opposite trend is observed, here the stripping tends to accelerate slightly for better resolution. We hence find that the results become more similar for better resolution. In particular, in the high resolution case, VPH is in fact quite close to AREPO. When the sum of the star forming gas mass and the newly formed stars is considered, we see that this quantity is very robust for SPH and VPH as a function of resolution, with the offset probably simply reflecting different systematics related to the density estimation at the interface between gas disk and surrounding medium, as in Figure 7. On the other hand, the relative constancy of this quantity at late times reflects the inherent inability of these particle schemes to mix low entropy gas with gas of much higher entropy.
Both SPH and VPH rely on an artificial viscosity to capture shock waves and to damp out small scale numerical noise. The choice of the artificial viscosity parametrisation and the associated strength of the viscous forces are crucial for the overall robustness and accuracy of the scheme. Using a large viscosity leads to more robust shock capturing, but on the other hand it may produce substantial viscous effects in regions where the gas really ought to flow without dissipation (Cullen & Dehnen, 2010). In particular, a too high viscosity may suppress the growth of fluid instabilities that are important in gas stripping processes. We have therefore checked how our wind tunnel results for VPH depend on the viscosity setting. To this end we have reduced the viscosity parameter from our default value of (high viscosity) to (intermediate) and (low artificial viscosity). In Figure 14, we compare the corresponding results for the stripping as a function of time in our low-resolution wind tunnel set-up. We see that for reduced viscosity the gas is actually stripped a bit faster in VPH than for the high viscosity setting. This is consistent with our expectation that the numerical viscosity will tend to damp small-scale fluid instabilities and turbulence, which as a side effect reduces the stripping rate.
Another numerical nuisance parameter in VPH lies in the strength of the shape correction forces that can be introduced into the technique in order to regularise the mesh and encourage a quasi-regular particle distribution. As explained in Heß & Springel (2010), this is accomplished by adding a small term to the Lagrangian that penalises large aspect-ratio distortions of cells and large offsets between the generating point of a cell and its centre-of-mass. While not strictly necessary for VPH to work, we have found the technique to be considerably less noisy in certain conditions when such weak shape correction forces are used. In our standard implementation, their strength is controlled by two dimensionless parameters and (see Heß & Springel, 2010, for the exact definition of these parameters). In Figure 15, we show the dependence of the stripping results for VPH when instead stronger shape correction forces described by or even are invoked. In both cases, this has only a very weak influence on the results, as desired.
5 Stripping of a galaxy during cluster in-fall
In order to complement our wind-tunnel experiments carried out in the previous section with more realistic setups, we here want to conduct a few simulations where galaxy models are in-falling into live models of galaxy clusters. This approach includes a full treatment of gravity as well as a modelling of cooling and star formation. In particular, it accounts correctly for the tidal effects of a galaxy travelling within the cluster potential. The latter inevitably changes the structure of the galaxy, most prominently by reducing its dark matter mass through tidal truncation, which in turn changes the conditions under which the hydrodynamical processes occur.
Besides looking at the stripping of the gas, we will also compare how the star formation rates in the in-falling galaxy decline in our different simulations, since we expect that different numerical ram pressure will strongly affect this quantity (e.g. Kronberger et al., 2008; Kapferer et al., 2009). We note that a number of recent studies both with grid-based codes (Roediger & Brüggen, 2007; Iapichino & Niemeyer, 2008) and SPH codes (Jáchym et al., 2007; McCarthy et al., 2008) have begun to investigate this question in detail. Our focus here will be much more limited though, and only be concerned with a comparison of our new VPH scheme relative to SPH and AREPO.
5.1 Setup of galaxy-cluster interaction simulations
For definiteness, we consider a parabolic encounter of a disk galaxy with a small galaxy cluster. The galaxy is constructed as a compound system as described in Section 3 but scaled down to fewer particles than in our earlier simulations. This was necessary to reduce the computational cost to an acceptable level, given that in this setup we need to represent the much larger cluster with the same mass resolution as the galaxy in order to avoid numerical problems in SPH (Ott & Schnetter, 2003; Read & Hayfield, 2011). In our default set-up, we have therefore chosen gas particles for the ISM of the in-falling galaxy, and gas particles of identical mass for the IGM of the galaxy cluster. This low resolution will of course limit our ability to properly resolve hydrodynamical instabilities. We note however that similar and often still lower resolution is routinely used in cosmological simulations of galaxy formation, so our results are directly representative of such studies. They in any case allow us the extend our previous tests in Section 4 to lower resolution and investigate how reliably they can be assumed to carry over to more complicated situations.
|Code||softening baryons||softening DM||regularisation|
For the cluster, we adopt a gas-fraction of and a NFW concentration of . A detailed list of main numerical parameters adopted for our simulations in this section is given in Table 3. The trajectory of the encounter starts at a distance of , just outside of the virial radius, and follows a parabolic orbit where the minimal distance of the two objects would be if they were point masses. Relative to the orbital plane, the galaxy’s disk is tilted by and then turned by in the azimuth (see Duc et al., 2000, for a sketch of the orbital geometry). We note that we can expect to be able to draw quite general conclusions from a single choice of inclination angles since the orbital geometry has been shown to have no significant effect on the gas stripping (Roediger & Brüggen, 2006).
5.2 Properties of the head wind
In Figure 16, we show the properties of the wind encountered by the in-falling galaxy as a function of time. We give the time evolution of the galaxy’s environment as found in a sector with opening angle of a spherical shell in the radial range centred around the direction the galaxy is heading to. The velocity difference shown in the fist panel of Figure 16 is computed with respect to the centre of mass of the galaxy. The figure shows that density, pressure and temperature of the head wind reach their peak at the pericentre, as expected.
We include results for VPH, SPH and AREPO in Fig. 16, and as the comparison shows, there is generally very good agreement between the three different simulation techniques. All three runs show a strong rise of density and pressure as the galaxy approaches and passes pericentre, while the temperature is roughly constant, reflecting the nearly isothermal conditions in the cluster gas. There are no significant deviations in the orbit of the galaxy between the different methods, hence any difference in the evolution of the galaxies can only arise from differences in the hydrodynamical treatment of the interaction of galaxy and cluster gas.
5.3 Gas stripping and star formation truncation
In order to define whether a gas particle or cell still belongs to the galaxy we use the condition:
Furthermore, we additionally require that the separation of particles from the centre-of-mass of the remaining gravitationally bound dark matter of the galaxy is less than . We note that with this definition the galaxy can in principle also accrete new gas from the cluster which may then cool onto the ISM and help to provide fuel for star formation.
In Figure 17, we show a visual comparison of the gas stripping in VPH and SPH simulations, where it is readily apparent that the gas mass loss proceeds much slower in SPH than in VPH. This is confirmed by quantitative measurements shown in Figure 18, which include results for gas stripping in three different simulations of the encounter of a galaxy with the cluster ICM. Interestingly, VPH looses gas here even somewhat faster than AREPO, but both methods yield a substantially faster loss of gas than SPH. While at time the VPH simulation has lost all of the gas, the galaxy in SPH still retains half of it.
This substantial difference is corroborated by the behaviour of the star formation rate in the galaxies, shown in cumulative form in Figure 19. Here SPH shows the slowest decline overall, while both VPH and AREPO lead to a rapid termination of star formation, which implies a quick reddening of the galaxies. It is a well known problem in galaxy formation to understand the colours of cluster galaxies in detail. Semi-analytic models of galaxy formation have commonly predicted a rather quick truncation of star formation upon cluster in-fall, yielding cluster populations that are actually too red (e.g. Weinmann et al., 2006). It can be expected, due to the faster stripping, that our VPH and AREPO simulations suffer from the same problem, potentially making it even more severe. This may indicate that the stripping efficiencies in all of the simulations are substantially too large, probably because the ISM is still under-resolved and appears as a homogeneous dense phase instead of being resolved into a true multi-phase medium. The latter would make the medium more porous, allowing very dense clouds of gas to resist the ICM wind for a longer time and to stay in the in-falling galaxy.
6 Cosmological cluster simulations
We now turn to results for fully cosmological simulations of cluster formation. We first carry out a comparison of the gas content of satellite systems in ‘zoom’ simulations of the formation of rich galaxy clusters, carried out with VPH, SPH and AREPO. We here primarily want to check whether the differences we have observed in the stripping of dense ISM gas out of galaxies in our earlier more idealised simulations manifest themselves also in non-radiative simulations without star formation, where the density contrast is much smaller. Secondly, we study the well known Santa Barbara cluster of Frenk et al. (1999) in order to investigate whether VPH produces a higher entropy in the cluster centre compared to SPH, which would then make it closer to the results of mesh codes that have been applied to this problem. We also use this cluster in order to assess the numerical convergence of VPH for the properties of the intra-cluster gas.
6.1 Gas stripping in non-radiative zoom simulations of galaxy clusters
In order to simulate the formation of rich galaxy clusters in the CDM cosmology, we extract a massive halo from the Millennium Simulation (Springel et al., 2005), and re-simulate this cluster with the addition of gas. To this end we trace the particles that make up the cluster back to the original initial conditions at , thereby finding the Lagrangian region out of which the cluster has formed. This region is then populated both with dark matter and gas particles, which are perturbed with the original displacement field. We can also increase the resolution compared to that in the original simulation if desired, in which case additional small-scale fluctuation power is added in the region between the old and new Nyquist frequencies. Further away from the region that holds the cluster material and its immediate surroundings, the resolution is reduced by combining particles into progressively heavier ‘boundary’ particles. In this way, the resolution gradually declines with distance from the cluster while the gravitational tidal field that influences its formation is still accurately determined. With this standard ‘zoom’ technique, the computational cost is concentrated in the small region of interest, allowing high resolution simulations of individual objects in comparably short time.
For definiteness, we have picked a cluster of virial mass , which we re-simulate with a baryon fraction . The other cosmological parameters are the same as in in the original Millennium run and are given by , , , , and . We have initialised the re-simulations at and evolved them with VPH, SPH and AREPO to the present epoch. We employed an identical gravitational softening length of for all particle types and hydrodynamical schemes, and treated the gas as a non-radiative mix of hydrogen and helium.
|Code||resolution||softening baryons||softening DM|
We identify halos in the simulations by applying the FOF algorithm to the high-resolution dark matter particles. Each gas particle is assigned to the group in which its closest dark matter particle lies. We then apply the SUBFIND algorithm (Springel et al., 2001) to the particle groups we found in this way in order to decompose them into gravitationally bound (sub)groups. Using the IDs attached to each particle we can track individual halos/subhalos as a function of time, and, in particular, study how the gas content of subhalos declines as a halo falls into the cluster. In order to reduce numerical noise in our simulation comparison, we however restrict our substructure selection to a sample where more than of the DM particles can be found in every studied simulation run.
As an example, we show in Figure 20 the evolution of the gas fraction of two different substructures as a function of their distance to the cluster centre (the corresponding times are indicated by the redshift labels on the upper axis), comparing results for VPH, SPH and AREPO simulations. For both substructures, we find good agreement in the early in-fall phase, for distances larger than about . At smaller separations, the influence of the cluster is however very noticeable already, and the substructures loose gas quickly. In both of these examples, we identify a considerably larger gas loss in the VPH run than in SPH at the last time when we can still find the substructures before they are disrupted. This is consistent with our earlier findings for the stripping of the dense ISM gas.
In Figure 21, we now consider all of the substructures around the cluster at the final time, and simply compare the gas mass fraction in substructures as a function of cluster-centric radius. Even though the results are a bit noisy, we find a clear statistical trend of a smaller gas mass fraction in the VPH run relative to SPH, especially in the region , where the gas fraction declines rapidly. Interestingly, a simulation with AREPO for the same cluster initial conditions shows a bound gas fraction that is still smaller in-falling dark matter halos. We hence conclude that even at the level of non-radiative simulations there are already significant differences in the stripping of gas out of in-falling substructure, which also implies that the mixing of gas in massive halos can be expected to differ in the three examined numerical techniques.
6.2 The Santa Barbara Cluster
The ‘Santa Barbara cluster comparison project’ of Frenk et al. (1999) analysed the results of a large number of different cosmological hydrodynamical codes for the formation of a rich cluster of galaxies with non-radiative gas. The comparison involved both SPH codes and hydrodynamical mesh codes, and focused, in particular, on the resulting thermodynamic properties of the intra-cluster gas. An important result that emerged from the study was that the different methods systematically disagree in the amount of entropy predicted for the central regions of the cluster, with the mesh-based approaches yielding consistently higher central entropy and correspondingly lower density than the SPH codes.
The Santa Barbara (SB) cluster has become a standard test problem for cosmological hydrodynamic codes, with results reported in numerous studies. Recently, some studies have suggested that the difference seen between the various techniques is primarily associated with differences in the treatment of mixing (Mitchell et al., 2009; Vazza, 2011), which is suppressed in SPH by construction, and this may artificially lower the central entropy. Specifically, it has been suggested that this problem may be related to a suppression of Raleigh-Taylor fluid instability in SPH (Wadsley et al., 2008). If the difference is caused by the lack of mixing at the particle level (Tasker et al., 2008; Wadsley et al., 2008) the VPH results should in principle agree with SPH.
In the spirit of the original project of Frenk et al. (1999), we here re-run the SB cluster with our new hydrodynamical VPH method. We are especially interested in the question whether VPH differs in its predictions for the central cluster region compared with SPH, which perhaps could arise due to the different stripping efficiency of this technique. We use the same initial conditions that have been used in the original SB cluster comparison project, where an Einstein-de-Sitter cosmological model with mean density , Hubble constant , and baryon fraction of was used. The initial conditions were constructed as a constrained realisation where a rich cluster corresponding to a peak was imposed to form in the centre of a cubic box with a side-length of , in an otherwise random realisation.
We have simulated the SB cluster at a variety of different resolutions using VPH, ranging from to particles with parameters as denoted in Table 4, in order to see how well the method converges for the primary cluster properties. In Figure 22, we show the results of this convergence test, in terms of radial profiles for gas density, dark matter density, enclosed gas fraction, temperature, entropy, and specific X-ray emissivity. We find that the convergence is in general good for the outer profiles and the dark matter properties. Only in the very centre some interesting differences can be observed. For the most part they can be simply understood as an imprint of the varying spatial and mass resolution across the sequence, where shortly before the resolution limit is reached the lower resolution runs peel away from the converged result seen in the higher resolution simulations. This is for example the case for the entropy and gas density profiles.
It is clear however that VPH does not predict the formation of a large constant entropy core region, unlike what is typically seen in mesh codes and consistent with the falling entropy profiles that are typically observed with SPH all the way to the centre of halos. In fact, in general our results obtained with VPH for the Santa Barbara cluster are in good agreement with SPH results reported in the literature (e.g. Frenk et al., 1999; Ascasibar et al., 2003; Springel, 2005), suggesting that it makes comparatively little difference for the bulk thermodynamic structures which particle-based method is used in non-radiative cosmological simulations. In particular, the different stripping efficiency alone and the slightly higher spatial resolving power of the Voronoi-based density estimate do not yield an entropy core as generally found in hydrodynamical mesh-codes.
7 Discussion and Conclusions
In this study, we have carried out a systematic comparison of the properties of our new Voronoi particle hydrodynamics (VPH) method (Heß & Springel, 2010) with respect to standard SPH and moving-mesh hydrodynamics as implemented in the AREPO code. We have focused on stripping processes in galaxies and halos upon in-fall into galaxy clusters. Here, it is expected that the outer parts of gaseous disks are quickly removed due to ram pressure stripping (Gunn & Gott, 1972), but the subsequent more gradual gas loss sensitively depends on the ability of hydrodynamical codes to capture fluid instabilities occurring in shear flows around the galaxies. The recent findings that SPH appears to exhibit severe inaccuracies in this regime has prompted us to develop the alternative VPH method. In the present paper, we have studied how this new technique compares with traditional SPH and the new moving-mesh code technique in a number of basic problems relevant for galaxy formation.
To this end, we first compared results for isolated compound galaxy models, both in isolation and in wind tunnels where they were exposed to a supersonic head wind. This set-up allowed relatively high-resolution simulations of wind–ISM interactions. Our simulations have revealed non-negligible differences in the rate at which dense ISM gas is stripped and dispersed, and in the appearance of this gas in the downstream part of the flow. SPH showed a lower stripping rate than both VPH and the mesh-code AREPO. As a result, the SPH galaxy also experienced the largest displacement due to the ram pressure of the impinging wind. Also, we were able to show that essentially none of the ISM gas in SPH and VPH could ever be transferred to much lower density. Instead, if gas was stripped, it stayed in coherent dense blobs, where even star formation could continue. Furthermore due to elevated pressure caused by the surface tension effect in SPH, the star formation in these stripped blobs remained at an unphysically high level. AREPO, in contrast, showed a rapid loss of gas out of the disk, which was furthermore efficiently mixed with other gas, so that lower densities were reached quickly by the stripped gas and star formation was stopped.
We followed up these simulations with numerical experiments where we dropped galaxy models in live cluster models. Even though here the resolution was substantially lower, we obtained results in good qualitative agreement with our wind tunnel runs. Likewise, in non-radiative cosmological simulations of galaxy cluster formation, we followed individual subhalos as they fell into the forming cluster, finding again that the gas content of satellite systems declined slowest in SPH, while the stripping in VPH and especially in AREPO proceeded noticeably faster.
Finally, we considered simulations of the Santa Barbara cluster, which has become an important test problem for evaluating cosmological hydrodynamical codes. While the VPH runs revealed a slightly elevated entropy compared to SPH at the smallest radii, they in general agreed quite well with SPH and in particular did not provide evidence at the highest resolution that an entropy core similar to those found in mesh-codes such as AREPO is formed. This is perhaps to be expected if the entropy core indeed primarily arises from mixing processes (Mitchell et al., 2009) that are largely absent in SPH and VPH, by construction.
Overall, it thus appears that VPH offers some improvements over SPH without however changing its fundamental character. We argue that the most important origin of these differences lies in an improved gradient estimate in VPH compared to SPH. VPH is second-order accurate in the gradient estimates, i.e. a linear gradient is always reproduced exactly, independent of the particle distribution (Springel, 2010a). In contrast, SPH has a so-called zero-th order error in its gradient estimate (e.g. Read et al., 2010). This in particular means that even for equal pressures for all particles the pressure force not necessarily vanishes (Abel, 2011), and furthermore, the absolute size of the gradient error grows with the pressure itself. The gradient errors in SPH have also been linked to the creation of small-scale velocity noise in studies of subsonic turbulence (Bauer & Springel, 2011).
A good test problem for appreciating this difference in the gradient accuracy is the flow of a stable vortex as suggested by Gresho & Chan (1990). In this ‘triangular vortex problem’ a fluid is set up with an azimuthal velocity profile (see Appendix A for details), such that the pressure gradients counter the centrifugal force, and the vortex evolves in a time-independent, stable fashion. Figure 23 shows the radial velocity profile after the vortex has been evolved for a time with the VPH, SPH and AREPO codes in 2D, using a Cartesian grid for the initial conditions in the domain . It can be seen clearly that SPH shows a much larger velocity scatter than the other two codes, and its solution has already degraded quite noticeably relative to VPH and AREPO. Especially in the inner domain, where the fluid rotates like a solid body, SPH deviates systematically from the analytic solution. We note that part of this degradation can be influenced by the artificial viscosity setting (Springel, 2010b; Read & Hayfield, 2011), but if a higher viscosity is used to suppress the velocity noise it typically also leads to a faster viscous erosion of the vortex. In VPH, some velocity scatter is seen as well, but it is appreciably smaller than in SPH, which we interpret as a consequence of the more accurate gradient estimates in VPH.
Arguably one of the best ways to quantify the accuracy of the different numerical techniques for this analytic test problem is to look at an objective error measure as a function of resolution. In Figure 24, we compare the L1-error norm for the azimuthal streaming velocity of the Gresho test as a function of resolution for all three techniques. It is evident that ordinary SPH converges only very slowly, whereas VPH shows a considerably improved convergence rate. This directly demonstrates an important improvement brought about by VPH, which we can directly trace to better gradient estimates. The latter appears also as the primary reason for the better results for stripping we obtained with VPH compared to SPH, which is due to the more accurate treatment of contact discontinuities and the avoidance of the ‘gap’ phenomenon of SPH.
Note that the difference in convergence rate also means that VPH will outperform SPH in terms of computational cost if very high accuracy needs to be achieved, provided its computational cost differs only by a constant factor of order unity which is indeed the case in our present implementation. However, the Voronoi mesh construction adds substantial computational cost compared to ordinary “standard” SPH, making the VPH code about a factor 3-4 slower for pure hydrodynamics when the same number of particles is used. This is mitigated to some extent if self-gravity is included (which is typically about as expensive or slightly more expensive than SPH-based hydrodynamics), reducing the difference to less than a factor of 2. We note that some alternative suggestions to improve standard SPH, such as the scheme by Read et al. (2010) which involves a dramatic increase of the number of smoothing neighbours, also require an increased computational cost per resolution element. Which of these schemes is ultimately the most efficient one (in the sense of reaching a given accuracy goal with the smallest computational cost) is difficult to answer in general, and is in any case a problem-dependent and implementation-dependent question.
According to Fig. 24, VPH still falls short of the better convergence rate obtained with the moving-mesh code AREPO (and similarly with fixed grid mesh codes such as ATHENA, see Stone et al., 2008; Springel, 2010a). The above discussion suggests that higher order density estimates combined with at least equally accurate gradient estimates are needed to improve on SPH and VPH in this respect. Some suggestions in this direction have recently been made (e.g. Read & Hayfield, 2011; Maron et al., 2011; McNally et al., 2011), and it will be interesting to see whether they can successfully yield significant accuracy improvements in cosmological applications such as those discussed here.
The authors thank the anonymous referee for insightful comments about the paper. VS acknowledges support by the DFG Research Centre SFB-881 ‘The Milky Way System’ through project A1.
- Abadi et al. (1999) Abadi M. G., Moore B., Bower R. G., 1999, MNRAS, 308, 947
- Abel (2011) Abel T., 2011, MNRAS, 413, 271
- Agertz et al. (2007) Agertz O., Moore B., Stadel J., Potter D., Miniati F., Read J., Mayer L., Gawryszczak A., Kravtsov A., Nordlund Å., Pearce F., Quilis V., Rudd D., Springel V., Stone J., Tasker E., Teyssier R., Wadsley J., Walder R., 2007, MNRAS, 380, 963
- Ascasibar et al. (2003) Ascasibar Y., Yepes G., Müller V., Gottlöber S., 2003, MNRAS, 346, 731
- Baldry et al. (2004) Baldry I. K., Glazebrook K., Brinkmann J., Ivezić Ž., Lupton R. H., Nichol R. C., Szalay A. S., 2004, ApJ, 600, 681
- Balogh et al. (2000) Balogh M. L., Navarro J. F., Morris S. L., 2000, ApJ, 540, 113
- Balsara (1995) Balsara D. S., 1995, J. Comp. Phys., 121, 357
- Bauer & Springel (2011) Bauer A., Springel V., 2011, ArXiv e-prints, 1109.4413
- Boselli & Gavazzi (2006) Boselli A., Gavazzi G., 2006, Publications of the Astronomical Society of the Pacific, 118, 517
- Cha et al. (2010) Cha S.-H., Inutsuka S.-I., Nayakshin S., 2010, MNRAS, 403, 1165
- Cortese & Hughes (2009) Cortese L., Hughes T. M., 2009, MNRAS, 400, 1225
- Cullen & Dehnen (2010) Cullen L., Dehnen W., 2010, MNRAS, 408, 669
- Dolag et al. (2009) Dolag K., Borgani S., Murante G., Springel V., 2009, MNRAS, 399, 497
- Dolag et al. (2005) Dolag K., Vazza F., Brunetti G., Tormen G., 2005, MNRAS, 364, 753
- Domainko et al. (2006) Domainko W., Mair M., Kapferer W., van Kampen E., Kronberger T., Schindler S., Kimeswenger S., Ruffert M., Mangete O. E., 2006, A&A, 452, 795
- Dressler (1980) Dressler A., 1980, ApJ, 236, 351
- Duc et al. (2000) Duc P., Brinks E., Springel V., Pichardo B., Weilbacher P., Mirabel I. F., 2000, AJ, 120, 1238
- Frenk et al. (1999) Frenk C. S., White S. D. M., Bode P., Bond J. R., Bryan G. L., Cen R., Couchman H. M. P., Evrard A. E., Gnedin N., Jenkins A., Khokhlov A. M., Klypin A., Navarro J. F., Norman M. L., Ostriker J. P., Owen J. M., Pearce F. R., Pen U.-L., Steinmetz M., et al. T., 1999, ApJ, 525, 554
- Gingold & Monaghan (1977) Gingold R. A., Monaghan J. J., 1977, MNRAS, 181, 375
- Gresho & Chan (1990) Gresho P. M., Chan S. T., 1990, International Journal for Numerical Methods in Fluids, 11, 621
- Guedes et al. (2011) Guedes J., Callegari S., Madau P., Mayer L., 2011, ApJ, 742, 76
- Gunn & Gott (1972) Gunn J. E., Gott III J. R., 1972, Astrophys. J., 176, 1
- Gunn & Gott (1972) Gunn J. E., Gott III J. R., 1972, ApJ, 176, 1
- Hernquist (1993) Hernquist L., 1993, ApJS, 86, 389
- Heß & Springel (2010) Heß S., Springel V., 2010, MNRAS, 406, 2289
- Iapichino & Niemeyer (2008) Iapichino L., Niemeyer J. C., 2008, MNRAS, 388, 1089
- Jáchym et al. (2007) Jáchym P., Palouš J., Köppen J., Combes F., 2007, A&A, 472, 5
- Junk et al. (2010) Junk V., Walch S., Heitsch F., Burkert A., Wetzstein M., Schartmann M., Price D., 2010, MNRAS, 407, 1933
- Kapferer et al. (2009) Kapferer W., Sluka C., Schindler S., Ferrari C., Ziegler B., 2009, A&A, 499, 87
- Kennicutt (1998) Kennicutt Jr. R. C., 1998, ApJ, 498, 541
- Kronberger et al. (2008) Kronberger T., Kapferer W., Ferrari C., Unterguggenberger S., Schindler S., 2008, A&A, 481, 337
- Larson (1978) Larson R. B., 1978, Journal of Computational Physics, 27, 397
- Larson et al. (1980) Larson R. B., Tinsley B. M., Caldwell C. N., 1980, ApJ, 237, 692
- Lucy (1977) Lucy L. B., 1977, AJ, 82, 1013
- Maron et al. (2011) Maron J. L., McNally C. P., Mac Low M.-M., 2011, ArXiv e-prints, 1110.0835
- McCarthy et al. (2008) McCarthy I. G., Frenk C. S., Font A. S., Lacey C. G., Bower R. G., Mitchell N. L., Balogh M. L., Theuns T., 2008, MNRAS, 383, 593
- McNally et al. (2011) McNally C. P., Maron J. L., Mac Low M.-M., 2011, ArXiv e-prints, 1110.0836
- Mitchell et al. (2009) Mitchell N. L., McCarthy I. G., Bower R. G., Theuns T., Crain R. A., 2009, MNRAS, 395, 180
- Mo et al. (1998) Mo H. J., Mao S., White S. D. M., 1998, MNRAS, 295, 319
- Ott & Schnetter (2003) Ott F., Schnetter E., 2003, ArXiv Physics e-prints, 0303112
- Pelupessy et al. (2003) Pelupessy F. I., Schaap W. E., van de Weygaert R., 2003, A&A, 403, 389
- Price (2008) Price D. J., 2008, Journal of Computational Physics, 227, 10040
- Price & Monaghan (2007) Price D. J., Monaghan J. J., 2007, MNRAS, 374, 1347
- Quilis et al. (2000) Quilis V., Moore B., Bower R., 2000, Science, 288, 1617
- Read & Hayfield (2011) Read J. I., Hayfield T., 2011, ArXiv e-prints, 1111.6985
- Read et al. (2010) Read J. I., Hayfield T., Agertz O., 2010, MNRAS, 405, 1513
- Roediger & Brüggen (2006) Roediger E., Brüggen M., 2006, MNRAS, 369, 567
- Roediger & Brüggen (2007) Roediger E., Brüggen M., 2007, MNRAS, 380, 1399
- Roediger & Brüggen (2008) Roediger E., Brüggen M., 2008, MNRAS, 388, 465
- Rosswog (2009) Rosswog S., 2009, New Astronomy Reviews, 53, 78
- Schulz & Struck (2001) Schulz S., Struck C., 2001, MNRAS, 328, 185
- Serrano & Español (2001) Serrano M., Español P., 2001, Phys Rev E, 64, 046115
- Springel (2005) Springel V., 2005, MNRAS, 364, 1105
- Springel (2010a) Springel V., 2010a, MNRAS, 401, 791
- Springel (2010b) Springel V., 2010b, ARA&A, 48, 391
- Springel et al. (2005) Springel V., Di Matteo T., Hernquist L., 2005, MNRAS, 361, 776
- Springel & Hernquist (2002) Springel V., Hernquist L., 2002, MNRAS, 333, 649
- Springel & Hernquist (2003) Springel V., Hernquist L., 2003, MNRAS, 339, 289
- Springel & White (1999) Springel V., White S. D. M., 1999, MNRAS, 307, 162
- Springel et al. (2005) Springel V., White S. D. M., Jenkins A., Frenk C. S., Yoshida N., Gao L., Navarro J., Thacker R., Croton D., Helly J., Peacock J. A., Cole S., Thomas P., Couchman H., Evrard A., Colberg J., Pearce F., 2005, Nature, 435, 629
- Springel et al. (2001) Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726
- Stone et al. (2008) Stone J. M., Gardiner T. A., Teuben P., Hawley J. F., Simon J. B., 2008, ApJS, 178, 137
- Sunyaev et al. (2003) Sunyaev R. A., Norman M. L., Bryan G. L., 2003, Astronomy Letters, 29, 783
- Tasker et al. (2008) Tasker E. J., Brunino R., Mitchell N. L., Michielsen D., Hopton S., Pearce F. R., Bryan G. L., Theuns T., 2008, MNRAS, 390, 1267
- Thacker & Couchman (2006) Thacker R. J., Couchman H. M. P., 2006, Computer Physics Communications, 174, 540
- Tittley et al. (2001) Tittley E. R., Pearce F. R., Couchman H. M. P., 2001, ApJ, 561, 69
- Valcke et al. (2010) Valcke S., de Rijcke S., Rödiger E., Dejonghe H., 2010, MNRAS, 408, 71
- van der Wel et al. (2010) van der Wel A., Bell E. F., Holden B. P., Skibba R. A., Rix H.-W., 2010, ApJ, 714, 1779
- Vazza (2011) Vazza F., 2011, MNRAS, 410, 461
- Vazza et al. (2006) Vazza F., Tormen G., Cassano R., Brunetti G., Dolag K., 2006, MNRAS, 369, L14
- Vogelsberger et al. (2011) Vogelsberger M., Sijacki D., Keres D., Springel V., Hernquist L., 2011, ArXiv e-prints, 1109.1281
- Wadsley et al. (2004) Wadsley J. W., Stadel J., Quinn T., 2004, New Astronomy, 9, 137
- Wadsley et al. (2008) Wadsley J. W., Veeravalli G., Couchman H. M. P., 2008, MNRAS, 387, 427
- Weinmann et al. (2006) Weinmann S. M., van den Bosch F. C., Yang X., Mo H. J., Croton D. J., Moore B., 2006, MNRAS, 372, 1161
- White & Frenk (1991) White S. D. M., Frenk C. S., 1991, ApJ, 379, 52
Appendix A Gresho Vortex Test
In the test problem of a stable vortex as suggested by Gresho & Chan (1990) the azimuthal velocity profile has the form
in a gas of constant density equal to and adiabatic index of . The pressure is chosen as
so that the pressure gradients balance the centrifugal force.
In Figure 23, we showed the radial velocity profile after the vortex has been evolved for a time with the VPH, SPH and AREPO codes in 2D, using a Cartesian grid for the initial conditions in the domain . Interestingly, one can also clearly see that the velocity noise in VPH is still lower in the solid body rotation part of the vortex, for . Here the initial pressure gradients are correct in VPH, and they will in principle stay that way during the motion of the solid body part. However, the strongly shearing part between will eventually disturb the evolution – and this can in some sense be blamed on the density estimate.
To see this, let us imagine for the moment that the fluid particles would all move under the exact analytic pressure gradient, then they would move on circular orbits with constant angular velocity. The resulting time evolution of the point set is shown in the top four panels of Figure 25, together with the resulting density estimates at time for VPH and SPH in the lower panel. The shearing motion creates substantial irregularities in the particle distribution, manifesting itself in significant fluctuations of the density estimates around the desired background density of . These fluctuations tend to be roughly of comparable magnitude for the kernel-based density estimates in SPH and the Voronoi-based density estimates in VPH. In any case, the density fluctuations of course cause pressure fluctuation that necessarily prevent that the particles move on exactly circular orbits. Instead, additional motions are created that locally restore the pressure equilibrium. Motions that do not agree with the analytic characteristics are here needed to mitigate pressure fluctuations and to prevent that density errors as large as implied by the analytic trajectories occur in the first place. In contrast to SPH and VPH, in AREPO, the density values can stay correct because this code calculates mass exchanges between cells in such a way that the density of cells stays (nearly exactly) constant, despite the strong shearing of the tessellation. Combined with the accurate gradient operator in this code (which is the same as in VPH), this avoids spurious contributions to the pressure gradients, leading to higher accuracy in the final solution.
A further layer of complexity is added when the ambiguous role of the shape correction forces in VPH is considered. On one hand they induce deviations of particle motions from circular orbits in order to make cells more regular, which can degrade the accuracy with which the analytic solution is followed. On the other hand they maintain regularity of the cell configuration, which improves the robustness of the scheme and helps to ensure accurate density and gradient estimates (Heß & Springel, 2010). It is well possible that an improved cell regularisation scheme is able to improve VPH’s performance for the Gresho problem, something we leave for future investigations.