Modelling Shock Heating in Cluster Mergers: I. Moving Beyond the Spherical Accretion Model

1.3 MB·English

by  Ian G. McCarthy

Mon.Not.R.Astron.Soc.000,1–30(2006) Printed5February2008 (MNLATEXstylefilev2.2)
Modelling Shock Heating in Cluster Mergers:
I. Moving Beyond the Spherical Accretion Model
I. G. McCarthy1⋆, R. G. Bower1, M. L. Balogh2, G. M. Voit3, F. R. Pearce4
T. Theuns1,5, A. Babul6, C. G. Lacey1, and C. S. Frenk1
7
0 1Department of Physics, University of Durham, South Road, Durham, DH13LE
0 2Department of Physics and Astronomy, University of Waterloo, Waterloo, ON, N2L 3G1, Canada
2 3Department of Physics and Astronomy, Michigan State University,East Lansing, MI 48824, USA
4Physics and Astronomy Department, Universityof Nottingham, Nottingham, NG7 2RD
n 5Department of Physics, University of Antwerp, Campus Groenenborger, Groenenborgerlaan 171, B-2020 Antwerp, Belgium
a 6Department of Physics and Astronomy, University of Victoria, Victoria, BC V8P 1A1, Canada
J

1
1
Accepted XXXX.ReceivedXXXX;inoriginalformXXXX

1
v
5 ABSTRACT
3
3 The thermal history of the intracluster medium (ICM) is complex. Heat input
1 from cluster mergers, from AGN, and from winds in galaxies offsets and may even
0 preventthe cooling of the ICM. Consequently,the processes that setthe temperature
7 anddensity structureofthe ICMplaya keyroleindetermininghowgalaxiesform.In
0 this paper we focus onthe heating ofthe intraclustermedium during cluster mergers,
/
h with the eventual aim of incorporating this mechanism into semi-analytic models for
p galaxy formation.
- We generate and examine a suite of non-radiative hydrodynamic simulations of
o
mergers in which the initial temperature and density structure of the systems are set
r
t using realistic scaling laws. Our collisions cover a range of mass ratios and impact
s
parameters, and consider both systems composed entirely of gas (these reduce the
a
: physical processes involved), and systems comprising a realistic mixture of gas and
v
darkmatter.WefindthattheheatingoftheICMcanbeunderstoodrelativelysimply
i
X by considering evolution of the gas entropy during the mergers. The increase in this
quantityinoursimulationscloselycorrespondstothatpredictedfromscalingrelations
r
a based on the increase in cluster mass.
Weexaminethephysicalprocessesthatsucceedingeneratingtheentropyinorder
tounderstandwhypreviousanalyticalapproachesfailed.Wefindthat:(1)Theenergy
that is thermalised during the collision greatly exceeds the kinetic energy available
whenthesystemsfirsttouch.Thesmallersystempenetratesdeepintothegravitational
potential before it is disrupted. (2) For systems with a large mass ratio, most of the
energy is thermalised in the massive component.The heating of the smaller systemis
minoranditsgassinkstothecentreofthefinalsystem.Thiscontrastswithspherically
symmetric analytical models in which accreted material is simply added to the outer
radius of the system. (3) The bulk of the entropy generation occurs in two distinct
episodes. The first episode occurs following the collision of the cores, when a large
shock wave is generated that propagates outwards from the centre. This causes the
combinedsystemtoexpandrapidlyandovershoothydrostaticequilibrium.Thesecond
entropy generation episode occurs as this material is shock heated as it re-collapses.
Both heating processes play an important role, contributing approximately equally
to the final entropy. This revised model for entropy generation improves our physical
understanding of cosmologicalgas simulations.
Key words: cosmology:theory — galaxies: clusters: general — galaxies: formation
⋆ E-mail:[email protected](IGM)
2 I. G. McCarthy et al.
1 INTRODUCTION ings expected for their mass and temperature distributions
(Evrard&Henry1991; Kaiser1991; Bower 1997; Balogh et
ThestructureoftheX-rayemittinggasthatisseeningalaxy al. 1999; McCarthy et al. 2003). By knowing the entropy
groupsandclustersisstillnotunderstood.Thisgas(thein- distribution of the gas, we can determine the density and
tracluster medium, hereafter ICM) dominates the baryonic temperatureprofilesoftheICMwithin agivendarkmatter
mass content of clusters and is material which is left over gravitational potential (e.g., Voit et al. 2002). This process
from galaxy formation. As such its properties present an requires that an outer boundary condition is set (e.g. by
important clue to the galaxy formation process. If the pro- computingtheexternalpressureduetoinfallinggas)butthe
cesses that set its temperature and density structure could detailsoftheboundaryconditionhaveonlyaweakeffecton
be understood they should provide valuable constraints on theprofile.Discussingclusters inthelanguage ofentropyis
galaxy formation models. However, linkingtogethermodels extremely powerful, since the gas responds adiabatically to
forgalaxyformation andaccuratenumericalmethodscapa- slowchangesinthegravitationalpotential.Furthermore,the
ble of tracing the hydrodynamic evolution of the ICM has buoyancyoftheICMensuresthatarelaxedsystemwillhave
proveddifficult.Ifcoolingisneglected,theemission proper- an ordered structurewith thelowest entropygas located in
ties of clusters can be calculated in a robust manner (e.g., thedeepestpartofthegravitationalpotential.Onlycooling
Evrard et al. 1996); however such simulations cannot form andshockheatingormixingeventsaltertheentropydistri-
galaxies and thus omit a fundamental component of the butionofthegas:coolinglowersthegasentropy,whileshock
physics.Ontheotherhand,ifcoolingisincludedinthesim- heating and mixing can only raise it. Cooling in clusters is
ulations, the results are not stable to changes in numerical already well understood and can be modelled with simple
resolution unless an effective form of heating is introduced semi-analytictechniques(e.g.,McCarthy etal.2004). How-
tolimitthecoolingrateinsmallorearlyhalos(e.g.,Balogh ever,weneedtounderstandhowentropyisgeneratedduring
et al. 2001; Borgani et al. 2006). Many different “feedback” cluster growth to self-consistently model the cosmological
mechanisms have been tried to oppose the cooling insta- formation of structure in a semi-analytic way.
bility. Examples include galaxy winds and thermal energy
So what sets the entropy distribution in clusters? Pre-
injection (e.g., Springel & Hernquist 2003), delayed cooling
vious papers have looked at how entropy is generated in
(e.g.,Kauffmannetal.1999),preheatingofintergalacticgas
smoothly infalling gas (Cavaliere et al. 1998; Abadi et al.
(e.g., Kaiser 1991; Evrard & Henry 1991), thermal conduc-
2000; Tozzi & Norman 2001; DosSantos& Dor´e2002; Voit
tion(e.g.,Bensonetal.2003;Dolagetal.2004),andheating
et al. 2003). For a spherically symmetric smooth shell, it is
by AGN (e.g., Dalla Vecchia et al. 2004; Sijacki & Springel
possible to determine the infall velocity and density of this
2006). These schemes have all met with various levels of
material at the cluster virial radius. The bulk infall veloc-
success, buthavenotyet achieveda simultaneousmatch to
ityisconvertedintointernalthermalenergyatashocksur-
the observed galaxy luminosity function, stellar mass frac-
roundingthecluster.Theentropygeneratedintheshockcan
tionandtemperatureanddensitystructureoftheICM.One
becomputedfromtheRankine-Hugoniotshockjumpcondi-
factorslowingthedevelopmentofthesetheoreticalmodelsis
tions.Modellingthegrowthofclustersinthiswayproduces
the lack of understanding of the physical processes at work
realistic entropy profiles, however the normalisation of the
in the simulations. Because it is so difficult to quantify the
entropy is somewhat too high, so that the model predicts
likely effect of introducing new processes, the models must
clusterswith average densitiesρgas, thatare lower than ob-
bedeveloped on a largely trial and error basis.
served (see Voit et al. 2003).
An alternative approach is to develop a semi-analytic However, the smooth accretion approximation is not
model of the thermal history of the ICM, in which the valid in a cold dark matter (CDM) dominated universe
complex physicsis encapsulated in a small numberof semi- since most of the mass accreted by clusters has previously
empiricalequations.Thisapproachhasbeenverysuccessful collapsed into smaller virialised mass concentrations (e.g.,
in improving our understanding of galaxy formation (e.g., Lacey&Cole1993;Rowleyetal.2004;Cohn&White2005).
Cole et al 2000; Springel et al. 2005; Bower et al. 2006). This leads to a problem for the picture above since the en-
The semi-analytic approach has been tried in several stud- tropygenerated dropsrapidly if theaccreted material is al-
ies (e.g., Wu et al. 2000; Bower et al. 2001; Benson et al. ready dense and warm prior to the accretion shock. Impor-
2003), but the results have been limited because the tech- tantly, if this deficit is present at every step in a system’s
niquesforsettingthegasdistributioninthehaloshavebeen merger history, the problem becomes greatly compounded
adhoc.Whatisneededistotakeastepback,andtobetter over time. This effect is far larger than the overestimate of
understand the results of simple gas hydrodynamic mod- the cluster entropy obtained in the smooth accretion case.
els. A better physical understanding of how entropy of the Modellingarealisticstructureintheaccretedmaterialleads
ICM is set in these simulations would equip us with the to systems that are much denser and more luminous in X-
techniquesneededtoattackthemuchgreaterphysicalcom- rays than observed systems. This issue is discussed by Voit
plexity of the problem when cooling and galaxy formation et al. (2003).
are included. This is the subject of this paper. Recently, considerable progress has been made in sim-
The starting point for the next generation of semi- ulating the growth of clusters using purely numerical tech-
analytical approaches to the ICM should be to describe niques(e.g., Borgani etal. 2004; Kravtsovetal. 2005). The
the physical state of the gas in terms of its entropy dis- numericalsimulationsclearlyshowthatρgas ≈fbρtot(where
tribution function. Entropy is a powerful concept for un- fb is the cosmic baryon fraction, and ρtot is total baryon
derstanding the density and temperature structure of clus- + dark matter density) at large radii, in agreement with
ters, and was initially introduced as ameans of quantifying observations (e.g., Vikhlinin et al. 2006; McCarthy et al.
the departures of cluster structures from the simple scal- 2007). This result holds for a wide range of simulation pa-
Shock Heating in Cluster Mergers 3
rametersandislargelyinsensitivetonumericaltechniqueor Explorationofidealisedtwo-bodymergerscouldprovidethe
resolution (Frenket al. 1999). Theagreement between sim- key insights necessary to achieve this goal. These idealised
ulations andobservationsindicates thatwe donot yethave mergers should be representative of typical mergers in cos-
an adequate understanding of the lumpy accretion process, mological simulations and should themselves reproduce the
sinceanalyticattemptstomodelthisprocessyieldresultsin propagation of (near) self-similarity. However, this is by no
strongdiscord with theobservations.Ourimmediate aim is means a guaranteed result. First, idealised approaches such
to improveour understandingof these simulations. Clearly, as those adopted below may be too simplistic to mimic
theapproximatemethodswewishtodevelopwillnotreplace closelyenoughthetypicalmergereventincosmologicalsim-
direct simulations, butthey doprovidean essential tool for ulations.Forexample,thedegreetowhichtheself-similarity
understandingtheirresults,forestimatingtheimpactoflim- ofsystemsformedinnon-radiativecosmological simulations
itednumericalresolution,andforexploringthevastparam- issensitivetotheexactinitialconditionsofthesystems(i.e.,
eter space of heating mechanisms used in galaxy formation structuralproperties) ortotheproperties of themerger or-
models.Insubsequentpapers,wewillapplythisunderstand- bits is unclear. There are few studies in the literature to
ingtoimprovethetreatmentofgasheatinginsemi-analytic date that examine the dependence of the properties of the
codesgalaxyformationmodels(e.g.,Coleetal.2000;Bower baryonson thedetailed merger history of a system. On the
et al. 2006). By introducing our models for the shock heat- otherhand, theadopted resolution of ouridealised mergers
ing of the ICM, we will be able to self-consistently model istypicallybetterthanthatofcosmological simulations. So
additional heating from supernovae and AGN. At present, it is possible that we may resolve phenomena that are not
thebesteffortstotracktheimpactofAGNheatinginsemi- yet adequately resolved in cosmological simulations. In ei-
analyticmodels(e.g.,Boweretal.2001)cannotfollowprop- thercase,thereproductionofself-similarity inouridealised
erly thedevelopment of entropyand usean ad hoc descrip- simulations is not an automatic result. Below, we explore
tion based on energy conservation. underwhat conditions self-similarity is achieved in ouride-
In this paper we address the problem of lumpy accre- alised merger simulations.
tion by focusing on idealised simulations of merging sys- The following subsections present a description of our
tems. We consider simple two-body collisions and explore simulation setup and methods. This discussion goes into
how the entropy of the final system is generated. In setting some detail and the reader who is mainly interested in the
up our two-body collisions, we remain as faithful as possi- results of the simulations may wish to skip ahead to§3.
ble to the results of cosmological simulations. We initialise
oursystemswithstructuralpropertiesguidedbytheresults
2.1 Initial conditions
ofrecentcosmological simulationsandexploreawiderange
of representative mass ratios and orbits. These simulations We make use of the public version of the parallel TreeSPH
show that the spherical accretion model previously consid- code GADGET-2 (Springel 2005) for our merger simula-
ered is apoor representation of thetruephysicalprocess in tions. The Lagrangian nature of this code makes it ideal
hierarchical cosmological models. We show that the kinetic for tracking the entropy evolution of specific sets of par-
energy of the infall lumps is much greater than previously ticles (or in principle even on a particle by particle basis)
estimated. Furthermore, we show that, if the mass ratio of over the course of the simulations. By default the code im-
the accretion event is large, only a small fraction of the ki- plementstheentropy-conservingSPHschemeof Springel&
neticenergyisdissipatedinthesmallercomponent.Mostof Hernquist (2002). This is ideal for our purposes since the
the heating effect is felt by the larger system. Our simula- codeexplicitlyguaranteesthattheentropyofagasparticle
tions allow us to present a much improved physical model will be conserved during any adiabatic process. Thus, we
for entropy generation in clusters. are assured that any increase in this quantity is rooted in
The present paper is outlined as follows. In §2, we physics.
present the details of our simulation setup, including a dis- We perform a series of binary merger simulations in-
cussion of the initial structural conditions of our systems, volving systems ranging in mass from 1014M⊙ 6 M200 6
the adopted orbital parameters of the mergers, and other 1015M⊙ with the mass of the primary system (henceforth,
characteristics of the simulations. In §3, we analyse the en- we refer to the ‘primary’ system as the more massive of
tropyhistoryofthesystemsduringthemergingprocessand the two systems) in all collisions set by default to M200 =
qualitatively compare it with the standard spherical accre- 1015M⊙.(M200 isdefinedin§2.1.1.)Thus,wesimulatecolli-
tion model. In §4, we present an analytic model that en- sions characterised by mass ratios ranging from 10:1 to 1:1.
capsulates the essential physics of the merging process in Although we have experimented with a variety of mass ra-
our idealised simulations. Finally, in §5, we summarise and tiosinthisrange,wepresenttheresultsforsimulationswith
discuss our findings. mass ratios of 10:1, 3:1, and 1:1 only, but note that these
yield representative results. Even though we have chosen
to focus on the high end of the mass spectrum, the results
should also be applicable to mergers involving lower mass
2 SIMULATION SETUP AND METHODS
systems (e.g., galaxies) so long as the mass ratio is similar.
The baryons in virialised systems formed in non-radiative Thisisbyvirtueofthefactthatgravity,whichisscale-free,
cosmological simulationsapproximatelytracethedarkmat- is completely driving the evolution of the systems. One of
terwith ρgas ≈fbρtot.Anyreasonable modelofshock heat- the reasons for choosing to focus on clusters is that grav-
ing ought to reproduce this basic result. However, as out- itational processes dominate their final properties (X-ray
linedin§1,ithassofarprovedquitedifficulttoconstructa observations demonstrate that non-gravitational processes
physicalanalytic model that can successfully pass thistest. such as cooling and AGN feedback influence the very cen-
4 I. G. McCarthy et al.
tral regions only of massive clusters; see, e.g., McCarthy et dP(r) GM(r)
al. 2007). So we expect that our idealised merger simula- dr =− r2 ρgas(r) (3)
tions will be directly applicable to the bulk of the baryons The entropy1, K, is then deduced through the equa-
in real clusters, although we leavea comparison with X-ray tion of state P = Kρ5/3. A boundary condition must be
gas
(and/orSZeffect)observationsforfuturework.Asoutlined
suppliedbefore equation(3) can besolved andwe specify a
in §1, our primary goal is to develop a physical description
value for the pressure of the gas at r200, where we truncate
of entropy generation in mergers.
thegasprofiles(unlessstatedotherwise).Thereissomefree-
dominourchoiceofthepressureattheedgeofthesystem.
2.1.1 Gas-only models However,forphysicallyreasonablevaluesofP(r200)thereis
in fact very little difference between the resulting profiles.
We have set up and run collisions of systems composed For the bulk of the gas, we find that the requirement of
purelyofgas.Ofcourse,observationsdemonstratethatdark hydrostatic equilibrium forces K(Mgas) into a near power-
matter dominates by mass the baryons in galaxies, groups, law distribution with an index of approximately 1.3. This
andclustersandisoffundamentalimportanceintheprocess can be understood by noting that at intermediate radii the
ofstructureformation.Inthisrespect,itmightbeexpected NFW profile is approximated well by an isothermal profile
that the gas-only merger simulations will have limited ap- (i.e., ρ ∝ r−2), which has an entropy profile characterised
plicability.Ontheotherhand,weanticipatethattheseide- by dlogK/dlogMgas =4/3. To ensure that our initial sys-
alisedcollisionswillbeconsiderablymorestraightforwardto tems are structural copies of one another and to ease the
interpret than the gas + dark matter (hereafter, gas+DM) comparison between the initial and final systems later on,
simulations (described below), given that the former have avalueofP(r200)thatestablishes thisnearpurepower-law
only a single phase to consider. If a physical understand- entropy distribution all the way out to the system’s edge is
ing of the gas-only runs can be achieved then this could be selected. In real systems, the hot diffuse baryons are sup-
quite helpful for understanding the role of dark matter in posedly confined by the ram pressure of infalling material.
thecombinedruns.Infact,itisdemonstratedin§3thatthe Experimentingwith aphysicalmodelforthisram pressure,
gas+DM collisions exhibit what may be regarded as fairly Voit et al. (2002) have demonstrated that groups and clus-
minordeviationsfromtheresultsofthegas-onlyruns.Thus, ters are indeed expected to have near power-law entropy
wefindthegas-onlysimulationsarepotentiallyanexcellent distributions out to large radii.
tool forelucidating thekeyphysicalprocesses in mergers of Theanalyticprofilesmustbediscretisedintoindividual
massive systems. particles for inputintotheGADGET-2 code. Inparticular,
The systems are constructed to initially be structural theinitialparticlepositions,velocities,andinternalenergies
copies of one another. In particular, the gas is assumed to (perunit mass) must be specified.
havea NFW density profile (Navarro et al. 1997): For the initial particle positions, we start by generat-
ρ ing a glass using the “MAKEGLASS” compiler option of
ρ(r)= (r/r )(1+s r/r )2 (1) GADGET-2.Basically, aPoisson distributionofparticles is
s s
generated and is then run with GADGET-2 using −G in
where ρs=Ms/(4πrs3) and place of G for Newton’s constant. We allow the simulation
to run for ≈1000 time steps to ensure the glass achieves a
Ms = ln(1+r200/rs)−(Mr220000/rs)/(1+r200/rs). (2) uniformdistribution.Thisuniformdistributionisthenmor-
phed into the mass profile that corresponds to the density
In the above, r200 is the radius within which the mean profilegiveninequation(1).Thisisdonebyselectingapoint
density is 200 times the critical density, ρcrit, and M200 ≡ inside the distribution of particles, which we have taken to
M(r200)=(4/3)πr2300×200ρcrit. bethecentreof mass, ranking all of theparticles according
To define the density profile for a halo of mass M200, to the distance from this point, and then moving each par-
we must set the scale radius, rs. The scale radius is often ticle radially such that the desired mass profile is achieved.
expressed in terms of a concentration parameter, c200 ≡ For the particle velocities, the baryons are initially as-
r200/rs.Therearenumerousstudiesthathaveexaminedthe sumed to be at rest in their hydrostaticconfiguration.
trend between concentration parameter and mass in pure The specific internal energy, I, defined as I =
darkmattercosmological simulations(e.g.,Ekeetal.2001). (3/2)P/ρgas, of each particle is assigned by using theparti-
The fact that there is a trend at all implies that the dark cle’s distance from the centre of the halo and interpolating
matterhalosinthesesimulationsarenotstrictlyself-similar. within theanalytic profiles derived above.
Butthemassdependenceoftheconcentration parameteris Lastly, we surround the systems with a low density
weak, varying by only a factor of 2 or so from individual medium with a pressure set to P(r200). This medium is
galaxiestogalaxyclusters(withlowermasshalosbeingmore dynamically-negligible and issimply putinplace toconfine
concentrated).Afixedconcentrationparameterof c200 =4, the gas particles near the system’s edge. Without a confin-
avaluetypicalofgalaxyclusterssimulatedinaΛCDMcon- ing medium as much as 20% of the system’s mass can leak
cordance cosmology, is adopted for all of oursystems. beyondr200.
Inordertocompletely specify theproperties ofthegas
we must choose an entropy profile and an outer boundary
condition. As outlined in §1, the gas entropy is probably 1 We adopt this common re-definition of entropy, which is re-
the most useful quantity to track during the simulations. lated to the standard thermodynamic specific entropy via a log-
In order to specify the entropy profiles of the systems, it is arithm and an additive term that depends only on fundamental
assumed that thegas is initially in hydrostaticequilibrium: constants.
Shock Heating in Cluster Mergers 5
σr(r25)=rG3Mr2525 (5)
We have also experimented with the boundary condi-
1 tion that σ2(r)ρ(r)→0 as r →∞ and obtained very simi-
r
lar results within r25. Henceforth, results are presented for
the runsthat make use of the boundary condition given by
equation (5).
Eachofthethreevelocitycomponentsforadarkmatter
particleareassignedarandomvaluepickedfromaGaussian
distribution with width equal to the 1-D velocity disper-
sion,σ ,atthatradius.However,itisknownthatthislocal
r
Maxwellian approximation does not result in a system that
0.1 isinitially inequilibrium(thoughitisnotfarremovedfrom
one).Ofmost concernisthatasharptruncationofthesys-
tem at some finite radius3 results in a significant amount
of the mass near this radius seeping out to much larger
distance. Thus, the resulting equilibrium mass profile can
be significantly different from that which was originally in-
tended.Infact,thishappensatallradiibutinthesystem’s
0.1 1 interior the flux of particles moving to larger radii is offset
byparticlesmovinginwards(whichwereoriginallyatlarger
radii).Thisobviouslycannothappennearthesystem’sedge
Figure 1. Testing the local Maxwellian approximation for the where there are no particles originally at larger radii (and
primary system. The dashed blue curve shows the initial mass moving inwards) to replace those moving away from the
profileextendedouttor25.Thesolidredcurveshowstheresulting system’s centre. One way to overcome this problem is to
equilibrium mass profile after running the dark matter halo in
apply a smooth tapering of the density profile beyond the
isolationfor50Gyr(i.e.,manydynamicaltimes).Overtherange
system’s edge and compute the distribution function of the
the 0.1 6r/r200 61.0 the initial and final mass profiles follow
eachotherclosely. particles as opposed to making the Maxwellian assumption
(e.g.,Kazantzidisetal.2004;Pooleetal.2006).However,we
adoptamoresimplistic,butstilleffective,methodforover-
2.1.2 Gas + dark matter (gas+DM) models comingthisproblem.Inparticular,asmentionedabove,the
For the systems containing both gas and dark matter a initial darkmatterprofileis extendedwell beyondr200,out
slightlydifferentprocedureisusedtoconstructthesystems. tor25.Theaimistoensurethatwithinr200,ourregionofin-
terest,thereisalargeenoughinfluxofparticlestomaintain
First, we set up systems composed entirely of dark matter.
thedesired mass profile.
The density profiles of the systems are given by equation
Thepuredarkmatterhalosareruninisolationformany
(1). However, for reasons described below, we extend the
dynamical times. In Figure 1, we show that by extending
dark matter well beyond r200, to an overdensity of 25 (for
our adopted concentration r25 ≈ 2.44r200). As in the gas- the dark matter halo out to r25, the resulting equilibrium
mass profile for the primary system traces the initial (in-
only setup,a glass particle distribution ismorphed intothe
desiredmassprofile.Unlikethegas,however,thedarkmat- tended) profile within r200 quite well. A slight deviation
terparticleshavenothermalpressureandmustbeassigned is apparent for r < 0.1r200, which is due to limited mass
resolution. However, we explicitly demonstrate in the Ap-
appropriatevelocitiesinordertomaintainthismassprofile.
pendix that this deviation has negligible consequences for
Toachievethis,wesolvethetheJeansequation(seeBinney
our equal mass mergers. In the case of our unequal mass
& Tremaine 1987):
mergers,thesecondarysystemsareresolvedwithfewerpar-
d[σ2(r)ρ(r)] GM(r) ticles and therefore the deviation at small radii is slightly
r =− ρ(r) (4)
dr r2 largerinthesesystems.However,aswedemonstratein§4.2
anddiscussfurtherintheAppendix,thevastbulkoftheen-
for the velocity dispersion profile, σ (r), of the systems.
r
ergyisthermalisedinthemoremassiveprimarysystem.So,
Equation (4) implicitly assumes that the orbits of the dark
matter particles are purely isotropic2. aslongastheprimarysystemiswellresolvedtheproperties
of the final merged system should be robust. This is what
To solve equation (4) it is assumed that virial equilib-
likely accounts for the fact that the properties of massive
rium holds at the edge of the dark matter halo (see, e.g.,
virialisedsystemsformedinnon-radiativecosmologicalsim-
Ricker& Sarazin 2001); i.e.,
ulations are robust to relatively large changes in numerical
resolution (see, e.g., Frenk et al. 1999).
2 Cosmologicalsimulationsdemonstratethatpureisotropyisvi- A drawback of extending the halo out to much larger
olated in the outer regions of clusters. However, given that our radiiisthatasignificantfractionofthemassislocatedout-
(hydrostatic) gas-only simulations yieldresults that are remark- sidetheregionofinterest(thus,wepotentiallywasteagood
ablyconsistentwithourgas+DMsimulations(see§3.2),thisim-
plies that the resulting structural properties of our systems are
not sensitive to how the initial (pre-merger) mass profiles are 3 Note that some sort of truncation is necessary since the mass
maintained. profilecorrespondingtoequation(1)divergesatlargeradii.
6 I. G. McCarthy et al.
dealofcomputationaleffort).Inordertomitigatethisissue,
Table 1.MergerOrbitalParameters.
wetaketheequilibriumconfigurationafterrunningthedark
halos in isolation and simply clip particles that lie beyond
Massratio Sim.type d0 vr vt
r50.Wehaveruntheclippedhalosforafurther10Gyrand (kpc) (km/s) (km/s)
verify that the mass profiles within r200 is virtually unaf-
1:1 gas-only 4126 1444 0
fected.This clippedequilibrium configuration isusedtoset
3:1 gas-only 3494 1444 0
up our initial gas+DM systems. We keep the current posi-
10:1 gas-only 3021 1444 0
tionsandvelocitiesforourinitialvaluesforthedarkmatter
particlesinthecombinedruns.Forthegasparticleswealso 1:1 gas+DM 4342 1444 0
use the positions of the dark matter particles (that is, for 1:1 gas+DM 4342 1400.9 350.2
thoseparticlesthatliewithinr200)butreflectthemthrough 1:1 gas+DM 4342 1291.6 645.8
thecentreofmasssothatthegasparticlesdonotliedirectly 3:1 gas+DM 3658 1444 0
3:1 gas+DM 3658 1400.9 350.2
on top of the dark matter particles. The gas particle veloc-
3:1 gas+DM 3658 1291.6 645.8
ities are set to zero and the internal energy densities are
10:1 gas+DM 3170 1444 0
determined by placing the gas in hydrostatic equilibrium
10:1 gas+DM 3170 1400.9 350.2
within the total gravitational potential. As in the gas-only
10:1 gas+DM 3170 1291.6 645.8
models described above, a pressure boundary condition is
selected such that K(Mgas) is nearly a pure power-law out
to r200. thepeak also dependson thetotal mass of systems(see his
The above procedure implies that we will have equal
Fig. 4). In particular, the orbits of mergers involving mas-
numbers of gas and dark matter particles in our combined
sive systems are more radial and less tangential than those
runs within r200 (i.e., double the number of particles used of mergers involving low mass halos, presumably because
intheisolateddarkmatterrunwithinthatradius).Inorder
moremassivestructurestendtoformattheintersectionsof
toconservemass, sothat theequilibrium stateforthedark
large-scalefilaments.Unfortunately,thesamplewasn’tlarge
matter is still valid, we re-assign the dark matter particle
enough to fully quantify this mass dependence.
masses to (1−f ) times that used in theisolated run while
b Forthepurposesofthepresentstudy,weadoptthefol-
the gas particles are assigned a mass of f times the dark
b lowing set of orbital parameters, which are representative
matterparticlemassintheisolatedrun.Heref istheratio
b of the results of Benson (2005). Unless stated otherwise,
of gas mass to total mass within r200, which is set to the a fixed relative total velocity corresponding to the circular
universal valueof Ω /Ω =0.02h−2/0.3=0.136.
b m velocityoftheprimaryatr200,≈1444km/s,isadoptedini-
The end result of the above procedure is that the gas tially for all of our simulations4. This is slightly lower than
traces the dark matter, both phases are initially in equilib-
the mean value found by Benson (2005), but consistent at
rium, and within r200 both phases have a form that is very the1-sigmalevel.Forthegas+DMsimulations,weexamine
nearly the intended NFW distribution. As in the gas-only
three different orbits for each mass ratio, corresponding to
models, we surround the gas+DM systems with a low den-
v = 0, v = v /4, and v = v /2. We refer to these as the
t t r t r
sity pressure-confining gaseous medium.
“headon”,“smallimpactparameter”and“largeimpactpa-
rameter” runs, respectively. In practical terms, this implies
2.2 Merger orbits and other simulation vt ≈ 0.243vc,p(r200) and vt ≈ 0.447vc,p(r200) for the latter
characteristics two. This choice spans the lower half of the vt/vr distribu-
tion for massive halos seen in Fig. 4 of Benson (2005). For
We follow the approach of Poole et al. (2006) and use the thegas-only simulations, which we run simply to help coax
results of cosmological simulations to help specify the or- out the key physics during mergers, we simulate head on
bital properties of our idealised merger simulations. In par- collisions only [i.e., vr =vc,p(r200) and vt =0].
ticular, we turn to the study of Benson (2005; but see also Each ofthesimulations areinitialised with thegaseous
Tormen1997;Vitvitskaetal.2002;Wangetal.2005;Khoch- components of the primary and secondary systems just
far & Burkert 2006). Benson (2005) used a large collec- barely touching. Thus, the initial separation, d0, is just the
tionofN-bodysimulations carried outbytheVIRGOCon- summation of r200 for the primary and r200 for the sec-
sortium to study the distribution of orbital parameters of ondary.Forthegas+DMsimulations,thisimpliesthatthere
substructure falling onto massive (primary) halos. Particu- isinitiallysomeoverlapofdarkmatterhalosoftheprimary
lar emphasis was given to the 2-D distribution of the rela- and secondary systems, which extends beyond the gaseous
tive radial and tangential velocities of the substructure as component.Acompletesummaryoftheadoptedorbitalpa-
they passed through a spherical shell with a radius set to rameters of all of our simulations is presented in Table 1.
r = (1±0.2)rvir, where rvir is the virial radius of the pri- Note that for a fixed mass ratio d0 differs slightly between
maryhalo. Asanticipated, thepeakof thedistributioncor- the gas+DM and the gas-only simulations. This difference
responds to approximately the circular velocity of the pri- arises owing to the slight deviation of the mass profiles in
maryhaloatthevirialradius(seehisFig.2).Inparticular,
v≡(vr2+vt2)1/2 ≈1.1vc(rvir),wherevr and vt aretherela-
tiveradialandtangentialvelocities,respectively.Forthein-
4 Wenotethatthecircularvelocityatr200differsbyonlyasmall
dividualcomponents,thepeakofthedistributionisapprox-
amount from the circular velocity at the virial radius, which is
imately centred on vr ≈ 0.9vc(rvir) and vt ≈ 0.65vc(rvir) whatBensonactuallycompareshisvelocitiesto.Thisjustreflects
andit isclear thatthereisacorrelation betweenvr andvt. the fact that the NFW mass profile gives rise to a nearly flat
However, as pointed out by Benson (2005), the centring of circularvelocityprofileatlargeradii.
Shock Heating in Cluster Mergers 7
the gas+DM simulations from the intended NFW distribu-
tion.Finally,thesimulationshavebeensetupsuchthatthe 2
physical and centre-of-mass (of the entire system) coordi-
1.8
nates and velocities are identical, where we have approxi-
mated theinitial configuration as two point masses. 1.6
We have carried out a mass resolution study of one of
1.4
our simulations (see the Appendix). Based on this study
we adopt the following conditions. The primary systems 1.2
in our default runs have 50,000 gas particles within r200. 1
This is the case for both the gas-only and gas+DM runs
and implies the gas particle mass is 2×1010M⊙ in the for- 0.8
mer and fb(2×1010M⊙) = 0.272×1010M⊙ in the latter. 1.8
In the gas+DM runs, the primary systems have ≈76,500
darkmatterparticleswithinr50and50,000withinr200 with 1.6
a dark matter particle mass of (1 − fb)(2 × 1010M⊙) = 1.4
1.728 ×1010M⊙. These particle masses are held fixed for
thesecondarysystemsaswell(thus,theycontainfewerpar- 1.2
ticles than theprimary).
1
Forthegravitationalsofteninglength,weadoptafixed
physical size of 10 kpc for both the gas and dark matter 0.8
0 5 10 0 5 10
particles.(Wehavealsoexperimentedwithsofteninglengths
of5kpcand20kpcandfoundthedifferencesintheentropy
evolution to be negligible.) We use a fairly standard set of
Figure2.Entropyevolutioninthe1:1gas-onlymerger.Thered
SPH parameters, including a viscosity parameter, αvisc, of and blue lines/hatched regions represent the primary and sec-
0.8, a Courant coefficient of 0.1, and the number of SPH ondarysystems,respectively.Theupperandlowerboundsofthe
smoothing neighbours, Nsph, is set to 50±2. hatched regions represent the 75th and 25th percentiles, respec-
Finally,eachsimulationisrunforadurationof13Gyr, tively, while the central (thick) lines represent the median. The
i.e., approximately a Hubble time. Particle data is written panelsshowtheevolutionoftheentropyofparticlesseparatedac-
out in regularly spaced intervals of 0.1 Gyr. cordingtotheirinitialpositioninthecumulativegasmassprofile.
Sincethecumulative gas massprofileisamonotonic function of
radius,the panels alsoseparate the particles according to initial
distancefromthecentreoftheirrespectivehalos.
3 SIMULATION RESULTS
Belowwepresentadetaileddiscussionoftheentropyevolu-
Some discussion of how the particles in each system
tionofoursimulations.Thereaderwhoismainlyinterested
actually get theirentropyiswarranted. Tohelp aidthedis-
in a general understanding of this progression may wish to
cussion we plot in Fig. 3 a sequence of snapshots of the 1:1
skip ahead to §3.3, which presents a summary of our find-
merger with particles colour coded according to the frac-
ings.
tional change in entropy since the previous simulation out-
put.Thegeneralprogressionofthemergercanbedescribed
as follows. Since the two systems are initially just barely
3.1 Gas-only simulations
touching and have a relative velocity of vc,p(r200), they be-
3.1.1 Entropy evolution gininteractingassoonasthesimulationstarts.Inparticular,
twolargeshockfrontsarequicklyestablishedasthesystems
We start by considering the evolution of the entropy in the
approach oneanother(seethepanelcorresponding tot=1
gas-only merger simulations. Plotted in Figure 2istheevo-
Gyr in Fig. 3). However, there is very little increase in the
lution for the 1:1 merger subdivided into spherical shells.
entropy of the particles in any of the four regions plotted
Immediately apparent is the high degree of symmetry in
in Fig. 2 until after the cores of the two systems collide at
the 1:1 merger. In all four regions the 25th, 50th, and 75th
t ≈ 2 Gyr into the simulation. The implication is that the
percentilesfollow each otherquiteclosely. Inaddition,with
initial shocks are actually quite weak (as the colour coding
the exception of the outermost ring, it is evident that the
in Fig. 3 would also indicate). A plausible explanation for
particlesinbothsystemshaveachievedanearlyconvergent
this behaviour comes when one compares the sound speed
state by the end of the simulation. Of course, it is possi-
of thegas, c , definedas
bletocontinuerunningthesimulation todetermineexactly s
what the convergent state of the outermost regions will be.
∂P γP 5 k T
However, since we have already run the simulation for a cs ≡r∂ρ =r ρ =r3µmb (6)
p
Hubble time this implies that the outer regions of systems
formedfromsimilarmergersincosmologicalsimulationswill withtheinitialrelativevelocityofthemerger.Sincetheini-
also not have had sufficient time to completely virialise by tial relative velocity of the merger is set to the circular ve-
the present day. Since one of our aims is to understand the locity,we are effectively comparing thesound crossing time
resultsofsuchsimulations,welimitthedurationofoursim- ofasystemwithitsdynamicaltime.Theconditionofhydro-
ulations to13 Gyr.Below we will compare thisfinalconfig- static equilibrium ensures that the gas temperature will be
uration with a scaled up copy of theinitial systems. such that these two timescales are comparable. In this case
8 I. G. McCarthy et al.
Figure3.Fractionalchangeinparticleentropyasafunctionoftimeforthe1:1gas-onlymerger.Shownareparticlesinaslicethrough
the centre of thickness 250 kpc (i.e., |z|<125 kpc). Particles are colour coded according to the fractional change in entropy since the
lastsimulation output (0.1 Gyr ago). Black points arefor afractional change of less than 2%, blue are fora 2%-10% change, cyan are
fora10%-20% change, greenarefora20%-30%change, yellow arefora30%-40% change, andredarefora >40% change. Forclarity
thesurroundingpressure-confiningmediumisnotdisplayed.Eachpanel is10Mpconaside.
wefindthemeansoundspeedofthesystemsis≈1200km/s. tem, succeeds in reversing the infall of the material. This
Compared to the initial relative velocity (see Table 1), this outflow of gas proceeds nearly adiabatically for a period of
impliesaninitialMachnumberofonlyM≈1.2.Therefore, time that is dictated by the initial distance of the particle
we should expect only mild shock heating to occur during fromitssystem’scentre.Forexample,fortheoutermostring
theearly stages ofthemerger(asobserved inFig. 2). How- plotted in Fig. 2 this period of adiabatic expansion occurs
ever,asweillustratein§4,bythetimethecorescollide the from t≈2.5 Gyrfor a duration of nearly 3.5 Gyr.Gas par-
relative velocity can approach or exceed nearly 3 times the ticles initially located close to the centre of their respective
initial value, effectively bringing themerger into the strong systemdonotspendsuchalargetimeinthisoutflowphase,
shock regime. as they end up being more deeply embedded within in the
merged system’s potential well.
Following the collision of the cores, a large shock wave
isgenerated.Thisshockquicklypropagatesoutwards,heat- The outflowing gas eventually halts and begins to re-
ing the gas that was initially predominantly on the far side accreteontothecoreofthemergedsystem.Althoughthegas
of each system and has not yet finished falling in (see the accretes from virtually all directions, the symmetry of the
panels corresponding to t=2 and 3 Gyrin Fig. 3). In fact, 1:1merger issuchthattheaccretion happenspreferentially
thisshockwave,combinedwiththeexpansionoftheshocked alongandperpendiculartotheoriginalcollisionalaxis.This
high pressure/density gas near the core of the merged sys- occursforasignificantperiodoftime,untilthemergedcore
Shock Heating in Cluster Mergers 9
2
1.8
2
1.6
1.4
1.5
1.2
1
1
0.8
1.8
2
1.6
1.4
1.5
1.2
1
1
0.8
0 5 10 0 5 10 0 5 10 0 5 10
Figure 4. Entropy evolution in the 3:1 (left) and 10:1 (right) gas only mergers. The red and blue lines/hatched regions represent the
primaryandsecondarysystems,respectively.Theupperandlowerboundsofthehatchedregionsrepresentthe75thand25thpercentiles,
respectively, whilethe central (thick) curves represent the median. The panels show the evolution of theentropy of particles separated
according to their initial position in the cumulative gas mass profile. Since the cumulative gas mass profile is a monotonic function of
radius,thepanelsalsoseparatetheparticlesaccordingto initialdistancefromthecentreoftheirrespectivehalos.
hasachieved anearspherical symmetry.Following this,gas 21/2 overtheentiresystem.Therefore,toalargedegree,one
continuestoaccrete at theoutskirtsof thesystem butdoes expects(andthesimulationsbearthisout)thatmostofthe
so in a more or less spherical fashion. As evidenced from particles will be shock heated to a similar degree. As a re-
Figs. 2 and 3, the gas is slowly being shock heated during sult,convectivemixingisminimal.Thisagreeswellwiththe
this period of re-accretion. The character of this episode of idealised cluster merger simulations of Poole et al. (2006).
entropy production therefore differs significantly from the Theentropyevolutionofthe3:1and10:1gas-onlymerg-
firstabruptepisodefollowingcorecollision.However,asFig. ers is plotted in Figures 4a and 4b, respectively. In a qual-
2 indicates, both episodes are of comparable importance in itative sense, both the primary and secondary systems in
terms of setting the final state of the gas. Finally, we note thesemergers behavein asimilar mannertothesystemsin
that the late time (nearly spherical) accretion of material the1:1case.Forexample,theybothgeneratetworelatively
is what gives rise physically to the continued increase of weak shock fronts early on, with the secondary driving a
entropyintheoutermostregionsofthesystemuntiltheend bow shock into the primary and vice-versa. As in the 1:1
ofthesimulations(e.g.,asinthebottom-rightpanelofFig. case, they exhibit two main periods of entropy production,
2). thefirstbeingassociatedwithastrongshockapproximately
An interesting question is whether ornot muchmixing whentheircorescollideandthesecondbeingassociatedwith
occurs as a result of the merging process. We have tested an extendedperiod of re-accretion shocks. Furthermore,we
this by tracking not only the entropy evolution of individ- alsoseelittleevidenceformixinginthe3:1and10:1mergers.
ual particles but also their spatial distributions. Interest- Ofcoursetherearesomedifferencesbetweenthethreesetsof
ingly, even though the shock heating of the systems occurs simulationsindetail.Forexample,becausethesoundspeed
in an “inside-out” fashion, we see very little evidence for of the gas in the secondary in the 10:1 collision is smaller
mixing. For example, particles that were initially located than that in the 1:1 case, the initial Mach number for the
nearthecentresofthetwo(pre-merger)halos(i.e.,particles secondary is higher. Consequently,the initial bout of shock
thatinitially haverelativelylowentropies) endupbeinglo- heatingbeforethecorescollideisrelativelymoreimportant
cated near the centre of the final merged system, whereas for thesecondary in high mass ratio mergers.
thelarge-radii(highentropy)particlesendupoccupyingthe It is clear from a comparison of Figs. 2 and 4a,b that
outerregionsofthefinalsystem.Aplausibleexplanationfor the higher the mass ratio of the merger the more strongly
this behaviour is that the initial Mach numberdistribution heated the secondary is relative to the primary. This again
of the particles varies relatively weakly with radius. This maybetiedtothefactthatasonedecreasesthemassofthe
can be understood by considering the initial temperature secondaryitscharacteristic temperaturedecreases resulting
profiles. The condition of hydrostatic equilibrium results in in a decreased sound speed and an increased Mach num-
initialtemperatureprofilesthatvarybylessthanafactorof ber. It is also expected if the preservation of self-similarity
2from theclustercentretoitsperiphery.Thisthenimplies is achieved by heating both the secondary and the primary
theMachnumberdistributionvariesbylessthanafactorof up to same level. To explain, K = P/ρ5/3 ∝ T/ρ2/3, but
gas gas
10 I. G. McCarthy et al.
2
initial primary
scaled primary
final (total) system 1.5
1
0.5
total
primary
secondary
0
10 100 1000 0.2 0.4 0.6 0.8 1
Figure5.Comparisonoftheresultingentropydistributionsfromthe1:1gas-onlymergerwithself-similarexpectations.Left:Comparison
ofradialprofiles.Theshort-dashedredcurverepresentstheinitialentropyprofileoftheprimary,thelong-dashedmagentacurverepresents
aself-similarlyscaledupversionofthisprofile,andthesolidgreencurverepresentsthefinalentropyprofileofthetotalmergedsystem.
Right: Comparison of Lagrangian gas mass profiles. The short-dashed red, long-dashed blue, and solidgreen curves represent the final
K(Mgas) distributions of the primary, secondary, and total merged systems, respectively. The horizontal dotted line represents the
self-similarscaling.
self-similarity means that all systems have the same inter- However,deviationsareclearly visible,particularly at large
nal mass structure and therefore the characteristic entropy radii. These are discussed in more detail below.
depends only on the temperature of the system. The virial Theentropydistributionsofobservedsystemsaremost
theorem relates the temperature of a system with its mass commonlypresentedinradialcoordinates(e.g.,Piffarettiet
via M ∝T3/2, and therefore thecharacteristic entropy of a al. 2005; Donahue et al. 2006; Pratt et al. 2006), and so
system scales simply as it is useful to present theoretical models in such units if
the goal is to explain the properties of observed systems.
K ∝M2/3 (7)
However, our immediate goal in this paper is to develop
Initially, therefore, the ratio of the secondary’s char- a physical shock heating model. For this purpose, the en-
acteristic entropy to that of the primary is K /K = tropy distributions are best presented in Lagrangian (gas
s p
(M /M )2/3. So the secondary requires extra heating to mass) coordinates. In Figure 5b, therefore, we plot the fi-
s p
overcome its initial deficit compared to the primary if self- nal K(Mgas) distributions5 for the merged system and for
similarity is achieved by heating both systems to the same theprimaryandsecondary individuallyforthe1:1 gas-only
level. merger. The gas mass coordinate has been normalised by
the total gas mass in the systems (in this case 1015M⊙ for
the primary and secondary and 2×1015M⊙ for the final
3.1.2 Final configurations mergedsystem).Theentropycoordinatehasbeenscaled to
Howdothefinalentropydistributionsofthegas-onlymerg- theself-similar expectations;i.e.,theinitial K(Mgas)distri-
erscomparewithself-similarexpectations?Figure5presents butionoftheprimaryscaledupbyafactorof(Mtot/Mp)2/3.
(Equivalently, we could use the initial entropy distribution
this comparison for the 1:1 gas-only merger both in tra-
ditional radial coordinates (left panel) and in physical gas of the secondary scaled up by a factor of [Mtot/Ms]2/3.)
Perhaps the simplest way to achieve this state is by hav-
mass coordinates (right panel). In Figure 5a, we compare
ing the primary and secondary systems individually obey
theinitial entropyradial profileoftheprimary systemwith
self-similarity following the merger. In this case, we would
the final radial profile of the total merged system. We also
expectthefinaldistributionoftheparticlesbelongingtothe
showaself-similarlyscaledupversionoftheprimary’sinitial
entropy profile, achieved by scaling the entropy coordinate primary to be (Mtot/Mp)2/3 times the initial distribution,
u(Mptboty/M(Mp)t1o/t3/M=p2)21//33.=Fig2.2/53aashnodwtshtehraatdtihael cscoaolreddinuaptevebry- w(Mhitloet/tMhes)s2e/c3onladragreyr’sthfiannalitdsisitnriitbiaultidoinstwriobuultdiobne. afactorof
sion of theprimary’s initial entropy profile follows thefinal
entropyprofileofthemergedsystemoverapproximatelytwo 5 The physical meaning of the K(Mgas) distribution is most
decadesinradius. Thisindicatesthat,byandlarge, the1:1 straightforwardly thought of in terms of its inverse, Mgas(K),
gas-onlymergerhasscaledself-similarlythroughthemerger. whichisthetotal massofgaswithentropylowerthan K.