Forming disc galaxies in major mergers II. The central mass concentration problem and a comparison of GADGET3 with GIZMO

2.1 MB·

by  S. A. Rodionov

Astronomy&Astrophysicsmanuscriptno.a2015_2 (cid:13)cESO2017
January11,2017
Forming disc galaxies in major mergers II. The central mass
concentration problem and a comparison of GADGET3 with GIZMO
S.A.Rodionov,E.AthanassoulaandN.Peschken
AixMarseilleUniv,CNRS,LAM,Laboratoired’AstrophysiquedeMarseille,Marseille,France
e-mail:[email protected]
ReceivedXXXX,XXX;acceptedXXXX,XXXX
7 ABSTRACT
1
Context.Inaseriesofpapers,westudythemajormergeroftwodiskgalaxiesinordertoestablishwhetherornotsuchamergercan
0
produceadiscgalaxy.
2
Aims.Ouraimhereistodescribeindetailthetechnicalaspectsofournumericalexperiments.

n Methods.Wediscusstheinitialconditionsofourmajormerger,whichconsistoftwoprotogalaxiesonacollisionorbit.Weshowthat
a suchmergersimulationscanproduceanon-realisticcentralmassconcentration,andweproposesimple,parametric,activeGalactic
J nuclei(AGN)-likefeedbackasasolutiontothisproblem.OurAGN-likefeedbackalgorithmisverysimple:ateachtime-stepwetake

0 allparticleswhoselocalvolumedensityisaboveagiventhresholdvalueandincreasetheirtemperaturetoapresetvalue.Wealso
1 comparetheGADGET3andGIZMOcodes,byapplyingbothofthemtothesameinitialconditions.
Results.Weshowthattheevolutionofisolatedprotogalaxiesresemblestheevolutionofdiskgalaxies,thusarguingthatourproto-

] galaxiesarewellsuitedforourmergersimulations.Wedemonstratethattheproblemwiththeunphysicalcentralmassconcentration
A inourmergersimulationsisfurtheraggravatedwhenweincreasetheresolution.WeshowthatourAGN-likefeedbackremovesthis
G non-physicalcentralmassconcentration,andthusallowstheformationofrealisticbars.NotethatourAGN-likefeedbackmainly
affectsthecentralregionofamodel,withoutsignificantlymodifyingtherestofthegalaxy.Wedemonstratethat,inthecontextofour
h. kindofsimulation,GADGET3givesresultswhichareverysimilartothoseobtainedwiththePSPH(densityindependentSPH)flavor
p ofGIZMO.Moreover,intheexampleswetried,thedifferencesbetweentheresultsofthetwoflavorsofGIZMO–namelyPSPH,and
- MFM(mesh-lessalgorithm)–aresimilartoand,insomecomparisons,largerthanthedifferencesbetweentheresultsofGADGET3
o andPSPH.
r
t Keywords. galaxies:kinematicsanddynamics–methods:numerical
s
a
[

1. Introduction torsishigh,soifwesimplyadopttheKennicutt-Schmidtlaw,we
1
v findthatthestar-formationefficiencywillbegreat(i.e.,thegas
5 consumptiontime-scalewillbesmall).Moreover,starformation
Anumberofobservationalworkshavearguedthatalargefrac-
8 will be strongly enhanced during the merger event (Larson &
tion of present-day spiral galaxies have experienced a major
6 Tinsley 1978). As a consequence, a significant fraction of the
merger event during the last 8 Gyr (e.g., Hammer et al. 2005,
2 gaswillnotsurvivethemerging,andwillturnintostarsbefore
02009a,b; Puech et al. 2012). This has motivated a number of itsend,whichwillleadtohighB/Tfraction.
.numerical simulations modeling the merging and its remnant.
1
However,mostdynamicalsimulationsofmajormergersbetween
0
two spiral galaxies made so far have initial conditions where InA16,weproposedamorerealisticsetupforthedynami-
7
eachgalaxyhasthe propertiesofnearbydiskgalaxies,with,in cal study of a major merger event between two spiral galaxies,
1
:thebestcases,anenhancedgasfractiontomimicdiskgalaxiesat introducinggasintheformofagaseoushalo,existentineachof
vintermediateredshifts(Barnes2002;Springel&Hernquist2005; theprotogalaxies.Wethuscollidenotacoupleoffullyevolved
i
XCox et al. 2006; Hopkins et al. 2009; Lotz et al. 2010; Wang localspiralgalaxies,butacoupleofprotogalaxies,consistingof
et al. 2012; Hopkins et al. 2013; Borlaff et al. 2014; Querejeta only dark matter and gas halos. Each protogalaxy, after 1 − 2
r
aetal.2015;Hopkinsetal.2009;Governatoetal.2009;Hopkins Gyrofevolutioninisolation,resemblesanintermediate-redshift
et al. 2013; Borlaff et al. 2014; Querejeta et al. 2015). Merger disk galaxy, while after 8 − 10 Gyr of evolution, it resembles
remnants from such simulations have, in general, a B/T (bulge a present-day spiral galaxy (see A16 and Appendix A). In the
tototalstellarmass)ratiothatistoohightoadequatelyrepresent current article, the second of a series, we discuss in detail the
spiralgalaxies.Thiscanbeeasilyexplainedassumingthatadisk technicalaspectsofournumericalexperiments.Insection2,we
canbeformedonlyfromthegasthatsurvivedthemerger(Hop- describetheinitialconditionsforoursimulations,whileinAp-
kinsetal.2009).InordertohaveasmallB/Tratio,itisnecessary pendix A we discuss the evolution of individual protogalaxies
thatalargefractionofgassurvivesthemergerevent.Forexam- inisolation.Insect.3,wedemonstratewhyweneedAGN-like
ple,inthemostfavorablecaseof100%gaseousprogenitors,for feedbackinourmodels,describetheadoptedfeedback,andtest
B/T ∼0.3,atleast70%ofthegasmustsurvivethemerger.But it. In sect. 4, we compare results of our simulations calculated
thisisrelativelydifficulttoaccomplishifallthegasisinitiallyin usingeitherGADGET3orGIZMOcodes.Weconcludeinsec-
thedisk.Thesurfacedensityofthegasinsuchgas-richprogeni- tion5.
Articlenumber,page1of14
A&Aproofs:manuscriptno.a2015_2
2. Simulations Theinternalenergyofeachgaseousparticleiscalculatedso
that it is in hydrostatic equilibrium in the halo + gas potential,
2.1. Initialconditions
takingintoaccountthespin,ifitexists.Inthecasewithoutrota-
2.1.1. Individualprotogalaxies tion,fromtheconditionofhydrostaticequilibriumandassuming
thatthepressureiszeroatinfinity,wehave:
In all our simulations, individual protogalaxies are initially
(cid:90) ∞
spherical and composed of a dark matter (DM) halo and a P(r(cid:48))= Gρ(r)M(r)/r2dr, (3)
gaseous halo, of masses M and M , respectively. Stars do
DM gas r(cid:48)
notexistintheinitialconditions,butformduringtheevolution.
where M(r) is the total mass inside a sphere with radius r, and
Towithinamultiplicativeconstant,thetypeoffunctionalform
ρ(r) and P(r) are the density and pressure, respectively, of the
forthedensityofboththeDMandthegaseoushalois
gas at a given radius. From the pressure, assuming the condi-
tion of ideal gas, we can calculate the specific internal energy
C sech(r/|r|) u = P/(ρ (γ − 1)), where γ is the adiabatic index (in all our
ρ(r)= xγi (xη+1)(γo−tγi)/η, (1) simulations,wesetγ=5/3,whichcorrespondstomono-atomic
gas).
where r is the spherical radius, x = r/r , and r and r are Whenever present, the spin is taken into account using the
s s t
characteristic radii of the halo that can be considered as mea- following simple approximation. We derive the thermal energy
sures of the scale length and the tapering radius, respectively. inthepolarregionsu fromequation(3),andcalculatethether-
p
The constants γ, γ , and η characterize the radial density pro- malenergyintheequatorialregionu fromthepressure:
i o e
fileandCisaglobalmultiplicativeconstanttosetthetotalmass (cid:90) ∞
tγooth=ed3e,sainreddηva=lue1.,FworhitchheDcoMrr,esinpotnhdiss,arutpiclteo,wtheeutrsuenγciat=ion1,, Pe(r(cid:48))= r(cid:48) ρ(r)(GM(r)/r2−v2φ/r)dr. (4)
to Navarro–Frenk–White (NFW) models (Navarro et al. 1996, Weinterpolatethethermalenergyas,u=ku +(1−k)u ,where
1997). For the gas, we have γi = 0, γo = 2, and η = 2 (Beta k = |θ|/(π/2), and θ is the angle between a pgiven directeion and
model for β = 2/3, e.g., Miller & Bregman (2015) and refer-
theequatorialplane.However,weshouldnotethatinallmodels
encestherein).Allparametersforeachmodeldiscussedherecan
consideredhere,thecorrectionforthespinisalmostnegligible
befoundinAppendixB.
anditwouldbeenoughtouseequation(3)withoutanycorrec-
Thehalocomponentwasbuiltasasphericalisotropicsystem
tions.
in equilibrium with the total potential including gas, using the
Models constructed in such a way are in equilibrium if the
distributionfunctiontechnique(Binney&Tremaine2008,§4.3).
gasisadiabatic(i.e.,withoutcoolingandfeedback).Duringthe
For this, we used the program mkhalo written by W. Dehnen
runs, however, the gas will cool radiatively and, getting out of
(McMillan & Dehnen 2007), which is part of the NEMO tool
equilibrium,fallinwards.InAppendixAwediscusstheevolu-
box(Teuben1995).Wealsohavethepossibilitytoaddaspinto
tionofsuchprotogalaxiesindetail.
thehalo.Thefractionofparticleswithapositivesenseofrota-
tionisgivenbytheparameter f .Incaseswithnonetrotation
pos
f = 0.5,whileifallparticlesrotateinthepositive(negative) 2.1.2. Collisions
pos
sensethisparameterisequalto1(-1).
In all cases considered here, we have mergers between two
Thegaseouscomponentwasconstructedasfollows.Foreach
equal-massgalaxies.However,caseswithunequalgalaxiescan
gasparticle,allvelocitiesexceptthetangentialvelocityv were
φ beeasilycreatedinourframework,andtheywillbediscussedin
settozero.Thetangentialvelocitywassettothevalueofv¯ ,
H,φ followingarticles.
themeantangentialvelocityinthehaloatthelocationofthegas
EachofthetwoprotogalaxiesisinitiallycreatedwithZasro-
particle.Thisvaluecanbecalculateddirectlyfromthedistribu-
tationaxis.Afterthat,theycanbearbitraryorientated,which,in
tionfunctionofthehalo.First,atagivensphericalradiusr,we
general,givesustwoEuleranglesperprotogalaxy.Lettheplane
calculatev¯(r),themeanofthemoduleofthevelocityvector.
ofthecollisionalorbitbetheX−Y plane.Wedescribethisorbit
byaninitialseparationinphasespace x ,y ,v ,v .Insome
0 0 x,0 y,0
(cid:82)vmaxvp (v)dv cases,wechoosetheinitialseparationfromatwo-bodyproblem
v¯(r)= (cid:82)00vmax pvv(v)dv , (2) (inastswumoipnogintthsa)t,aselltttihnegmthaessdeosfirthedeporrobtiotgeaclcaexnitersiciistyc,oinnciteinatlrasetepd-
aration,anddistanceinpericenter.Seealltheparametersofthe
√
wherev = 2Ψisthemaximalvelocityatagivenradiusr, modelsdiscussedhereinAppendixB.
max
Ψ is the negative of the total gravitational potential at a given
radius, p (v)∼v2DF(Ψ(r)−v2/2)istheprobabilitydensityofv
v 2.2. Code
atagivenradiusandDFisthedistributionfunction(Binney&
Tremaine2008,§4.3).Itiseasytoprovethatin3Dspace,ifwe We use a version of GADGET3 including gas and its physics
startfromanisotropicsystemandaddrotationbychangingthe (Springel 2005). The stars and the dark matter are followed
fractionofparticleswithpositiveandnegativesenseofrotation by N-body particles and gravity is calculated with a tree code.
(by introducing f ), the mean of the tangential velocity v¯ = ThecodeusesanimprovedSPHmethod(Springel&Hernquist
pos φ
(f −0.5)v¯. 2002) and sub-grid physics (Springel & Hernquist 2003). We
pos
It would have been possible to add spin to the system in a usethesameparametersforthesub-gridphysicsasSpringel&
different way. For example f in the halo does not have to be Hernquist (2003), which were calibrated with Kennicutt’s law
pos
constant and could be set as a function f (R,z), where R is (Kennicutt1998)).
pos
the cylindrical radius. Also, rotation in the gas could be set as Except for test simulations, we use a softening length
completelyindependentoftherotationinthehalo. of 25 pc for the gas and stars and of 50 pc for the halo,
Articlenumber,page2of14
S.A.Rodionov,E.AthanassoulaandN.Peschken:Formingdiscgalaxiesinmajormergers.
a 0.005 relative accuracy of the force, and the GADGET & Tremaine 2008). In order to demonstrate this, we consider
opening criterion 1, that is, a relative criterion that tries to two low-resolution test models, each with ten times less parti-
limit the absolute truncation error of the multiple expan- cles than mdf018 (see Appendix B). One model (mdf225) has
sion for every particle-cell interaction (see GADGET man- the same gravitational softening as mdf018: 50 pc for the halo,
ual, http://www.mpa-garching.mpg.de/gadget/users-guide.pdf). and25pcforgasandstars.Theothermodel(mdf214)hasdou-
We also use the GADGET system of units; that is, the unit of blethissoftening:100pcforthehalo,and50pcforgasandstars.
length is 1 kpc, the unit of mass is 1010 M , and the unit of Allparameters,exceptthenumberofparticlesandthesoftening
(cid:12)
velocityis1km/sec.Wecontinuedallsimulationsupto10Gyr. lengths,areidenticaltothoseofmdf018.Low-resolutionmodels
with the same gravitational softening (mdf225) at t = 2.5 Gyr
have a central peak on v (R) with almost the same maximum
c
3. Centralmassconcentrationinmergermodels
asinthehigh-resolutionmodelmdf018(seeFig.1,greenline).
Nevertheless, during the evolution, the central bump decreases
3.1. ModelswithoutAGN-likefeedback
much faster than in mdf018, and by t = 10 Gyr, the low reso-
3.1.1. Thecentralmassconcentrationproblem lution model mdf225 has no central bump on v (R) at all. This
c
fits well with the fact that the relaxation time is smaller for a
AsshowninA16andthefollowingpapers,oursimulationscan
smallernumberofparticles.Low-resolutionmodelswithasoft-
makegalaxieswhosepropertiesaregenerallyingoodagreement
eningtwiceasbig(mdf214),att = 2.5Gyr,haveacentralpeak
with those of real galaxies. Such comparisons include face-on
significantly less pronounced on v (R), which is expected due
and edge-on morphologies, radial density profiles, finding type c
to the decreased spacial resolution. During the evolution, how-
II profiles with inner and outer disc scale lengths, and break
ever, the central bump on v (R) decreases slower than in low-
radiiingoodagreementwiththoseobserved,etc..Wealsofind c
resolutionsimulationswithasmallersoftening(mdf225),which
good agreement with observations for vertical density profiles
isalsoexpected,becausetherelaxationtimeisbiggerforabig-
and thick disk properties, as well as for a number of kinematic
gersoftening(Athanassoulaetal.2001;Rodionov&Sotnikova
properties.Inparticular,theirrotationcurvesareflat,asobserved
2005).Thisanalysisarguesthattheloweringofthecentralpeak
byBosma(1981).
on v (R) in the high-resolution model mdf018 is due, at least
However, our models without AGN feedback have one no- c
partly,tothetwo-bodyrelaxation.
tableproblem,concerningtheverycentralpartofthedisk.Dur-
Asaresult,wecanexpectthathigher-resolutionsimulations,
ingthecollision,themodelsformanextremelydense,and,aswe
namely with a greater number of particles and smaller soften-
discuss below, non-realistic central mass concentration (CMC).
ing, will make the problem with CMC even worse, because a
Thiscentralmassconcentrationcanbebetterseenonthecircular
smallersofteningwillmakethepeakonv (R)higher,andabig-
velocityprofilev (R)1. c
c ger number of particles will prevent its attenuation with time.
LetusconsideratypicalexampleofamodelwithoutAGN:
Ontheotherhand,smallerresolutionspuriouslydiminishesthe
mdf018. All parameters of this model can be found in Ap-
problem(seeFig.1,fortwolow-resolutionmodels,mdf214and
pendixB.Thismodelisthenon-AGNanalogofmdf732,which
mdf225).Thismeansthat,insufficientlylow-resolutionmodels,
was thoroughly analyzed in A16. In Fig. 1 (red line), we show
onewillnothaveproblemswithanon-realisticCMC.Ofcourse,
thecentralpartofthev (R)profilesatdifferenttimesformdf018.
c decreasing the resolution cannot be considered as a solution to
Att=2.5Gyr,thatis,approximately1Gyrafterthemerger,the
thehighcentralconcentrationproblem.
centralbumponv (R)isveryhigh,andreachesalmost360km/s
c
(seefig.1).
Thestrongcentralpeakonvc(R)isunrealisticbyitself.The 3.2. AGN-likefeedback.Howwewillproceed.
circular velocity curves of the models can be compared with
OnecouldexpectthataddingAGNfeedbackwouldsolve,orat
observed HI or Hα rotation curves. Central bumps on rotation
leastalleviate,theproblemwiththecentralpeakinthecircular
curvesmayappearinrealgalaxies,butareconsiderablysmaller
velocityprofile.However,anAGNisnottrivialtomodel.There
(e.g.,Sofueetal.1999).Anotherproblemwiththestrongcentral
areinherentdifficultiesinunderstandingboththeaccretioninto
peaksonv (R)concernsbars.Oursimulationsshowthatsucha
c
thecentralblackhole(hereafterBH)(Hopkins&Quataert2010)
strong central peak stops, or delays beyond 10 Gyr, the forma-
and the subsequent energetic feedback (King 2003; Wurster &
tion of a bar (see also A16), while approximately two thirds of
Thacker2013).Moreover,weobviouslydonothaveenoughres-
diskgalaxieshavebars(Butaetal.2015).Bumpsofsimilarsize
olutiontomodelthesetwoprocessesdirectly.
did also appear in cosmological simulations until relatively re-
Nevertheless, there are several AGN feedback algorithms
cently, but have been lately strongly diminished by increasing
(see, for example, a comparative study of Wurster & Thacker
thefeedback(seeBrooks&Christensen2016,forareview).We
(2013),andpapersbyDuboisetal.(2010);Blechaetal.(2013);
will follow a similar route (see sect. 3.2), but before doing so,
Volonterietal.(2015);Gaboretal.(2016)).Usually,suchalgo-
we would like to continue the discussion of CMCs in a model
rithms include the explicit calculation of the accretion onto the
withoutanAGN.
black hole (thus the evolution of the mass of the black hole),
its movement (advection), and the merging of black holes in
3.1.2. DemiseoftheCMCinlow-resolutionmodels. cases of collisional simulations. But here we propose a much
simplerandlessambitiousapproach.Wedonotaimtocalculate
InFig.1,wecanseethatduringtheevolution,thecentralbump
theevolutionofthesuper-massiveblackholes;weonlywantto
onthev (R)profileinmdf018becomeslessandlessprominent
c removenon-physicalcentralmassconcentrationsinourmodels
withtime,andatt=10Gyrithasonlyamaximumof260km/s.
by adding physically motivated feedback. One way to proceed
Onecansuggestthatthisisduetotwo-bodyrelaxation(Binney
insuchasituationistointroducesomekindofsubgrid-physics
1 In this paper we calculate the circular velocity fro√m the total modeland“calibrate”itwithobservations.
massdistribution,assumingsphericalsymmetry:v (R)= GM(R)/R, Ourintentwastoinclude,inourmodel,asimple,parametric,
c
whereM(R)isthetotalmassinsideasphereofradiusR. AGN-likefeedbackwiththefollowingproperties:
Articlenumber,page3of14
A&Aproofs:manuscriptno.a2015_2
Fig.1.Centralpartofthecircularvelocityprofilesformodelsmdf018(highresolutionmodel),mdf225(lowresolutionmodel),andmdf214(low
resolutionmodelswithtwotimesbiggersoftening)attimes2.5,5,7.5,and10Gyr.Thisiscalculatedfromthetotalmassdistribution,assuming
sphericalsymmetry.
– Beingsimple. thresholdischosensoastoensurethattheparticlesarelocated
– Beingabletosolvetheproblemwiththecentralpeakinthe inthecenter-mostregion.
circularvelocityprofile. Inouralgorithm,wedonothaveasingleparticlerepresent-
– Allowingtheformationofabar. ingtheblackhole,sowedonotdirectlyfollowthemassofthe
– Being physically plausible, thus injecting a reasonable BH.Wecan,however,calculatetheBHmassfromthefeedback
amountofenergyinthecentralregionofthegalaxy. energy by, in some sense, inverting the formalism described in
Springeletal.(2005).WetakethesamevalueasSpringeletal.
– Influencingmainlythecentralregionofthemodel.
(2005) for the radiative efficiency of the BH, (cid:15) = 0.1, and the
r
Thelastpointrequiressomeexplanation.Itiscommonlybe- fractionofradiativeluminosity(cid:15)f =0.05,whichcancouplether-
lieved that AGN feedback could be quite strong and, presum- mallytothesurroundinggas.Consequently:
ably, be able to significantly quench star formation in the en- M (t)= E (t)/(c2(cid:15) (cid:15) )+M (0), (5)
tire galaxy (Springel et al. 2005; Di Matteo et al. 2005; Bundy BH feed f r BH
et al. 2008), thereby influencing more than just the central re- where M (t) is the BH mass at a given time t, and E (t) is
BH feed
gionofthegalaxy.However,aswementionedbefore,therearea the cumulated feedback energy up to time t. We stress that in
lotoffundamentaldifficultiesinpropermodelingofAGNfeed- thebasicversionofourfeedbackalgorithm(withoutEddington
back.So,ifweincludeanAGN-feedbackthatmodifiestheen- limit),wedonotneedtoknowM inthealgorithmitself.
BH
tiregalaxyinourmodel,wewillneedtocalibrateitwithsome The Eddington limit can be added as an option. It requires
observations, which is not trivial. Because of this, we have de- threenewfreeparameters:(cid:15) ,(cid:15) ,andM (0).KnowingM (t),
r f BH BH
cided that it would be reasonable, as a first step, to include an wecancalculatetheEddingtonluminosityas:
AGN-likefeedbackthatmainlyinfluencesthecentralregionof
4πGM (t)m c
themodel,andbarelyaffectstheouterparts. L = BH p . (6)
Edd σ
T
From this, we can calculate the energy available for the feed-
3.3. DescriptionofourAGN-likefeedback
back for the current time-step dt as E = L · dt · (cid:15) . We
Edd Edd f
Our AGN-like feedback is based on two parameters: a volume thencalculateE ,theenergyrequiredtoheatallparticleswith
req
density threshold ρ , and a temperature T . More specif- ρ > ρ uptoatemperatureofT .If E > E ,weheat
AGN AGN AGN AGN req Edd
ically, at every time step, we give internal energy to the gas particleswithaprobability p = E /E .Thisway,westatis-
Edd req
particles whose local volume density is larger than the thresh- tically obey the Eddington limit. For further discussion of the
old ρ , by increasing their temperature to T . The density EddingtonlimitseeSect.3.5.
AGN AGN
Articlenumber,page4of14
S.A.Rodionov,E.AthanassoulaandN.Peschken:Formingdiscgalaxiesinmajormergers.
3.4. ModelswithandwithoutAGN-likefeedback Our last model to compare is mdf791 (T = 1×107 K,
AGN
ρ = 2 M /pc3), which is similar to mdf732 (T =
AGN (cid:12) AGN
Herewecomparefivemodels;allofthemhaveidenticalparam- 1 × 107 K, ρ = 1 M /pc3), but has a ρ twice as high
AGN (cid:12) AGN
eters,exceptfortheAGN-likefeedbackparameters;
astheothers.Wecanseethatbyincreasingρ ,wedecrease
AGN
the efficiency of the AGN and make this model very similar to
– mdf018:ModelwithoutAGN-likefeedback mdf789(T =5×106K,ρ =1M /pc3,seeFigs.3and4).
– mdf789:TAGN =5×106 K,ρAGN =1M(cid:12)/pc3 ThiscanbeAGeaNsilyunderstoodAGbeNcausean(cid:12)increaseinρAGNmeans
– mdf732:TAGN =1×107 K,ρAGN =1M(cid:12)/pc3 thatfewerparticleswillhavetheirthermalenergyincreasedand
– mdf788:TAGN =2×107 K,ρAGN =1M(cid:12)/pc3 therebythetotalfeedbackwillbesmaller.
– mdf791:T =1×107 K,ρ =2M /pc3 The vertical thickness is a further noticeable difference be-
AGN AGN (cid:12)
tweenthemodelswithandwithoutAGN(seefig.2).Inthevery
TheremainingparametersaregiveninAppendixB. centralregion(lessthan2kpc),thedifferenceinthethicknessis
At t=10 Gyr, all these models are very similar. Their ra- duetothepresenceofabar,andparticularlytotheboxy/peanut
dial surface density profiles are virtually identical, except for bulgeassociatedtoit.Butthedifferencesoutsidethebarregion
the innermost region where models with no or weak feedback aremoredifficulttoexplain.ThemodelwithaveryweakAGN
(mdf018, mdf789, and mdf791) have a steeper cusp than the (mdf789), beyond the bar region, is approximately in between
others (fig. 2). It is the extra density in their cusps that is re- themodelwithoutAGNandtherestofthemodels.Presumably
sponsibleforthesteepmaximuminthecircularvelocitycurve, itissomehowconnectedtotheverydensecentralmassconcen-
alreadydiscussedinSect.3.1.1.Kinematicprofilesincludingthe trationbuttheexactmechanismwhichcausessuchadifference
meanazimuthalvelocityandtheradialvelocitydispersionpro- isnotclearforusyet.
filesarealsoverysimilar(Fig.2).However,comparingface-on Ingeneral,wecansaythatoursimpleAGN-likefeedbackis
views, we see differences in morphology (fig. 3). Despite the exactlywhatwewantedittobe.Itintroducesthedesiredconsid-
fact that azimuthally averaged profiles are very similar, in the erablechangesinthecentralpartofthemodel,asaimedfor,and
outerpartofthemodels,wecanseerelativelystrongdifferences onlyrelativelysmallonesintherestofthegalaxy.Morespecifi-
intheshapeofthespiralstructure.Thesedifferencesarepartly cally,itsolvestheproblemwiththecentralpeakofthev profile,
c
duetothetransientnatureofthespirals,butwecannotexclude andallowstheformationofthebar.
that they are also partly due to differences in the central part.
Anyway, the spiral structure has often been referred to as “the
3.5. ModelswithandwithoutEddingtonlimit
frostingonthecake”.Below,wewilldiscussthecentralpartof
themodels. The main reason for including the Eddington limit in an AGN
First,letusanalyzethemodelwithoutAGN(mdf018).Asal- feedback model is to make sure that the amount of feedback
readydiscussedinSect.3.1.1,duringthecollision,anextremely energyisreasonable.Thereis,however,bothobservationaland
dense and presumably non-realistic CMC is created. This can theoreticalevidenceforthepossibilityofsupercriticalaccretion
bebetterseenonthecircularvelocityprofile(vc(R))andisvery (King2003;Sa¸dowskietal.2014).Soevenfromageneralpoint
prominentatt = 5Gyr(seeFig.4forthemodelwithoutAGN, ofview,includingEddingtonlimitsmaynotbeobligatoryinor-
mdf018).Att=10Gyr,thispeakislessprominent,whichcould dertobephysicallyreasonable.
partly be due to the numerical two-body relaxation, but a part Herewewillshowthatinoursimple,parametric,AGN-like
of it, albeit small, could perhaps also reflect true two-body re- feedback model, the Eddington limit can be almost fully com-
laxation, as in globular clusters. Because of this strong central pensatedbyvaryingT .
AGN
mass concentration, we do not have a bar in this model (model Let us compare three models. These models have identical
mdf018,withoutAGN,seefig.3). parameters,exceptforT andtheEddingtonlimit:
AGN
Letusnowconsiderasequenceofmodelswithfixedρ ,
AGN
and increasing T : mdf789 (T = 5 × 106 K), mdf732 – mdf732,withEddingtonlimit,T =1×107 K.
AGN AGN AGN
(T =1×107K),andmdf788(T =2×107K).Thiscanbe – mdf751,withEddingtonlimit,T =1.5×107 K.
AGN AGN AGN
consideredasasequencewithincreasingfeedbackstrength.The – mdf726,withoutEddingtonlimit,T =1×107 K.
AGN
weakAGNinmdf789(T =5×106K)managestoonlypartly
AGN
remove the very dense CMC (see Fig. 4, t = 5 Gyr), which, The remaining parameters for these models are given in Ap-
however, is still relatively prominent, and at t = 10 Gyr, the v pendixB.
c
profileissimilartothatofmdf018(Fig.4,t = 10Gyr).Atthis All three models are very similar, although there is a small
time, we do have a bar, but it is quite small, with a semi-major difference in the mass distribution in the central part (see up-
axisofapproximately1.5kpcat10Gyr(seeFig.3formdf789). perpanelsoffig.5).Ifwecomparemdf732andmdf726,which
In mdf732, we increase T to 107 K, and in this model, the differ only in the presence, or absence of the Eddington limit,
AGN
central bump on the v profile is much less prominent, which we can see that in a model without Eddington limit (mdf732),
c
allows the formation of a bigger bar with a semi-major axis of theAGNissomewhatlessefficientandthecentralbumponthe
approximately3kpc(seefig.4andfig3formdf732).Afurther circular velocity profile is more prominent, which is expected
increaseofT to2×107 inmdf788fullyremovesthecentral becauseoftheEddingtonlimit.However,whenwemakeasim-
AGN
bump on v , but the bar size in this model is similar to mdf732 ulationwithEddingtonlimit,butincreasethevalueofT by
c AGN
(T = 1×107 K).Wemustnotethatonecanexpecttherela- 50%, we fully compensate for the presence of the Eddington
AGN
tionbetweenthebarandtheAGNtobemorecomplicatedthan limit(seemodelmdf751onfig.5).
asimplelinearcorrelation,andthatastrongAGN,byshuffling Letusconsideranothersetofthreemodels:
materialinthecentralregion,could,insomecases,delayoreven
preventtheformationofabar(theanalysisofthissubjectisbe- – mdf737,withEddingtonlimit,TAGN =2×107 K.
yondthescopeofthisarticle).Wenotethatmdf732isoneofour – mdf780,withEddingtonlimit,T =2.5×107 K.
AGN
referencemodelsandwasthoroughlyanalyzedinA16. – mdf730,withoutEddingtonlimit,T =2×107 K.
AGN
Articlenumber,page5of14
A&Aproofs:manuscriptno.a2015_2
Fig.2.ComparisonofvariousradialprofilesforthestellarcomponentoffivemodelswithdifferentAGNfeedbackparameters:mdf018(noAGN),
mdf789(T =5×106 K,ρ =1M /pc3),mdf732(T =1×107 K,ρ =1M /pc3),mdf788(T =2×107 K,ρ =1M /pc3),
AGN AGN (cid:12) AGN AGN (cid:12) AGN AGN (cid:12)
andmdf791(T =1×107 K,ρ =2M /pc3),allatt=10Gyr.Fromlefttorightandtoptobottom:surfacedensity,medianoftheabsolute
AGN AGN (cid:12)
valueofz,whichisagoodapproximationforthickness(Sotnikova&Rodionov2006),meanvalueoftheazimuthalvelocityandradialvelocity
dispersion.Theinlayintheupperleftpanelshowsthesurfacedensityintheinnermostregion(within3kpc).
The remaining parameters for these models are given in Ap- IfwecomparemodelswithandwithouttheEddingtonlimit,but
pendixB. withthesameT ,thenweobservethattheformerhasalarger
AGN
central bump (compare mdf737, which has an Eddington limit,
Thesemodelshaveamoreeccentricorbitandalargerinitial tomdf730,whichdoesnot,inthelowerpanelsofFig.5).Itis,
separationthanthepreviousones.Thisleadstoagreateramount however, possible to compensate for this by simply increasing
of low-angular-momentum gas in the central part, immediately T (seemdf780withEddingtonlimitandT =2.5×107 K
AGN AGN
after the collision, and consequently requires a more energetic vs. mdf730 without Eddington limit and T = 2×107 K in
AGN
AGN (bigger TAGN) in order to remove the central peak in the lower panels of Fig. 5). Thus, it is possible to compensate for
circular velocity profile and make the bar formation possible.
Nevertheless,thisdoesnotchangetheconclusionofthissection.
Articlenumber,page6of14
S.A.Rodionov,E.AthanassoulaandN.Peschken:Formingdiscgalaxiesinmajormergers.
Fig. 3. Face-on (upper row) and edge-on (lower row) views of the stars component for five models with different AGN feedback parameters:
mdf018(noAGN),mdf789(T =5×106 K,ρ =1M /pc3),mdf732(T =1×107 K,ρ =1M /pc3),mdf788(T =2×107 K,
AGN AGN (cid:12) AGN AGN (cid:12) AGN
ρ =1M /pc3),andmdf791(T =1×107K,ρ =2M /pc3),fromlefttoright,respectively,andallatt=10Gyr.Thesizeofeachsquare
AGN (cid:12) AGN AGN (cid:12)
boxcorrespondsto50kpc.
Fig. 4. Circular velocity profiles for five models with different AGN feedback parameters: mdf018 (no AGN), mdf789 (T = 5×106 K,
AGN
ρ =1M /pc3),mdf732(T =1×107K,ρ =1M /pc3),mdf788(T =2×107K,ρ =1M /pc3)andmdf791(T =1×107K,
AGN (cid:12) AGN AGN (cid:12) AGN AGN (cid:12) AGN
ρ =2M /pc3),at5(left)and10(rightpanel)Gyr.
AGN (cid:12)
the effect of the Eddington limit on the circular velocity curve cialviscosityterm(Cullen&Dehnen2010),asdescribedin
byvaryingtheAGN-likefeedbackparameters. Hopkins(2013,sect.3.1).
– MFM (Meshless Finite-Mass): Lagrangian formulation of
mesh-freealgorithmthatconservesparticlemasses(Hopkins
4. ComparisonwithGIZMO
2015).
4.1. Introductoryremarks
Phase boundaries, such as, for example, those between the
The GIZMO simulation code (Hopkins 2013, 2014, 2015) is
hotgasinthehaloandthemuchcoldergasinthedisk,orthose
a fork of GADGET3 (hereafter G3), using the same MPI par-
around cold dense clumps of gas in the hot halo can be better
allelization, domain decomposition, gravity solvers, and so on.
handled by PSPH or MFM, than by TSPH or G3 (see Hopkins
However,incontrasttoG3,GIZMOoffersseveralhydrodynam-
2015,andreferencestherein).
icalsolversfortheusertochoosefrom.Inthisarticle,wecom-
Since our aim here is to see how results will change with
pare the three solvers, which are briefly described below (see
differenthydro-solvers,wekeepthesamesubgridphysics.
Hopkins(2015,sect.4.1)formoreinformation).
ThemaindifferencewefoundbetweenG3andPSPH/MFM
– TSPH (“Traditional” SPH): this is the same, entropy con- is in the hot gaseous halo. This has small gaseous clumps in
serving,density-drivenformulationofSPHastheoneinG3. the simulation run with G3, which are either absent, or much
(Springel&Hernquist2002). lessprominentinthePSPHorMFMruns(seealsoTorreyetal.
– PSPH: this is a density independent version of SPH (Hop- 2012andHaywardetal.2014).However,inourproject,weare
kins 2013). It also includes a better treatment of the artifi- mostlyinterestedinthepropertiesoftheremnant.Therefore,in
Articlenumber,page7of14
A&Aproofs:manuscriptno.a2015_2
Fig.5.Comparisonofcircularvelocityprofiles(at5and10Gyr)formodelswithandwithoutEddingtonlimits.Theupperrowshowstwomodels
withEddingtonlimit(mdf732T =1×107Kandmdf751T =1.5×107K)andonemodelwithoutEddigtonlimit(mdf726T =1×107K).
AGN AGN AGN
Thelowerpanelsalsodisplaythreemodels:twowithEddigtionlimit(mdf737T =2×107K,andmdf780T =2.5×107K)andonemodel
AGN AGN
without(mdf730T =2×107K).
AGN
this comparative study, we will focus on the properties of the Ingeneral,theG3andthePSPHrunsgiveverysimilarresults,in
finalgalaxy. mostcasesdifferingbyapproximatelynomorethantheplotting
accuracy.Ontheotherhand,MFMgivessomewhathigherval-
ues,correspondingtoasomewhathighermass.Wenotethatthe
4.2. SimulationswithoutAGN-likefeedback
CMCfoundintheseruns,whichhavenoAGN-likefeedback,is
sufficienttoprohibitbarformation,oratleasttodelayitbeyond
We first compared G3 to TSPH and, as expected, we found no
differences.WewillthusnotdiscussTSPHfurtherhere. the10Gyrsoverwhichwemakethecomparisons.
Let us now compare three simulations of major mergers
without AGN. All three simulations were run with the same IntheupperleftpanelofFig.7,wecomparetheradialpro-
parameters, those of mdf018 (see Appendix B), except for the jectedsurfacedensityprofiles,asobtainedatt = 10Gyrforthe
softening, which is 50 pc for both dark matter and baryons, threesimulations.Itisclearthattheinnerdiscscalelength(i.e.,
but with a different hydrodynamical solver. We also used the thescalelengthofthepartofthediscwhichiswithinthebreak
same subgrid physics in all three cases. Thus, the only differ- intheprojectedsurfacedensity,aroundR = 17kpc)isapproxi-
enceisthehydro-solvers.Thesethreerunsaremdf624(runwith matelythesameforallthreemodels.Theouterdiscscalelength
G3), mdi101 (run with GIZMO/PSPH), and mdh101 (run with is approximately the same for models G3 and PSPH, while for
GIZMO/MFM). MFMitissomewhatsmaller,leadingtoafasterdropinthesur-
ComparisonsaremadeintheupperpanelsofFigs.6(forthe face densities in the outer disc region (i.e., the part of the disc
circular velocity curves) and 7 (for the projected surface den- beyondthebreak).Themiddlepanelcomparestheradialprofile
sity, the vertical thickness, and the mean tangential velocity of of a measure of the stellar disk vertical thickness, namely the
thestars).Intheformer,weseethatatt=2.5Gyr,allthreeruns azimuthalaveragedmedianoftheabsolutevalueofthezcoordi-
haveaverystrongCMC;thatoftheMFMrunbeingsomewhat natesofallstellarparticles(Sotnikova&Rodionov2006).The
strongerthantheothertwo.Nevertheless,inallthreecases,the threemethodsgivesimilarresults,and,althoughthedifferences
peak is considerably higher than 300 km/sec, for a velocity of arelargerthanwhatwefoundforthecircularvelocitiesandthe
theflatpartofthecircularvelocitycurveofapproximately200 projected surface densities, they are still relatively small, espe-
km/sec.Byt=10Gyr,thisCMChasconsiderablydecreased,as ciallytakingintoaccountthatdistthicknessisverysensitiveto
expected (see sect. 3.1.1), but still has a maximum of approxi- numerical effects (Rodionov & Sotnikova 2013). We also note
mately250km/sec,stillgivinganunrealisticshapeofthecircu- that in this plot, the MFM does not stand out as being further
larvelocitycurve,althoughmuchlesssothanthatatt=2.5Gyr. apartfromtheothertworuns.
Articlenumber,page8of14
S.A.Rodionov,E.AthanassoulaandN.Peschken:Formingdiscgalaxiesinmajormergers.
Fig.6. ComparisonofcircularvelocityprofilesforG3,PSPH,andMFMmodelsatt=2.5Gyr(leftpanels)andt=10Gyr(rightpanels).Theupper
rowcomparesmodelswithoutAGNandthelowerrowcomparesmodelswithAGN.
Finally the rightmost panel compares the radial profiles of Let us now compare three other runs, identical to the three
the three mean stellar tangential velocity radial profiles. Here, comparedintheprevioussubsections,butnowincludingthepro-
again,theG3andthePSPHareverysimilarandtheMFMdiffers posedAGN-likefeedback.
somewhatmore,butneverbytoomuch.
AfurthercomparisonisgiveninFig.8,wherewepresentthe – mdf627 : run with GADGET3. T = 1×107 K, ρ =
AGN AGN
evolutionofthebaryonicmass(gas+stars)insideasphereofra- 1M(cid:12)/pc3,withEddigtonlimit.
dius30kpcwithtime.WecanseethatintheMFMmodel,gas – mdh115:runwithGIZMO/PSPH.T =1×107K,ρ =
accretes on the disk somewhat faster than in the PSPH and G3 AGN AGN
0.75M /pc3,withoutEddigtonlimit.
models.Asaresult,att=10Gyr,thediskintheMFMmodelsis (cid:12)
approximately12%moremassive.Thereisalsosomedifference – mdi115:runwithGIZMO/MFM.TAGN =1×107K,ρAGN =
betweentheG3andthePSPHmodels,butthisismuchsmaller 0.75M(cid:12)/pc3,withoutEddigtonlimit.
than the difference between the MFM and the other two mod-
els. Indeed, at t=10 Gyr, the difference between G3 and PSPH
These simulations were run with 50 pc softening for both dark
amountstoonlyapproximately4%ofthebaryonicmass.These
matter and baryons. The remaining parameters are those of
differencesinaccretionratecouldexplainthecorrespondingdif-
mdf018(seeAppendixB)
ferences between the three circular velocity curves discussed
TheresultsareagaincomparedinFigs.6and7,butnowin
above.
the lower panels. We note that in all three cases this feedback
leadstoarealisticcircularvelocitycurvewithnostrongCMC,
contrarytowhatwefoundinrunswithnoAGN-likefeedback,
4.3. SimulationswithAGN-likefeedback
which all had strong CMCs. Furthermore, we witnessed that,
with the demise of the CMC, bars form in all three cases un-
Here, we again first compared G3 with GIZMO in the case of der consideration. In fact, it is due to these bars that there is
the same hydro solver (G3 vs. TSPH). We found that, with the a maximum of the disk thickness in the central region for all
sameparameters,ourAGN-likefeedbackisslightlylessefficient three models (lower middle panels of Fig. 7). Considering all
inTSPHrunsthaninG3runs.Itshouldbenotedthatevenwith thecomparisonsglobally,weseethatinthecasewithAGN-like
theTSPHhydrosolver,GIZMOisnotidenticaltoG3(Hopkins feedback,theMFMdiffersfromtheothertwomorethaninthe
2015,sect.F1).Nevertheless,thisdifferencecanbeeasilycom- cases without AGN-like feedback. The results of runs with G3
pensatedforbytuningthefeedbackparameters,andwedosoin andPSPHareagainveryclosetoeachother,similarlytocases
thefollowing. withoutAGN-likefeedback.
Articlenumber,page9of14
A&Aproofs:manuscriptno.a2015_2
Fig.7. ComparisonofvariousradialprofilesforG3,PSPH,andMFMmodels.TheupperrowshowsmodelswithoutAGNandthelowerrow
modelswithAGN.Ineachrow,fromlefttoright:surfacedensity,medianoftheabsolutevalueofz,whichisagoodapproximationforthickness
(Sotnikova&Rodionov2006)andmeanvalueoftheazimuthalvelocity
Fig.8. Comparisonofthebaryonicmass(gas+stars),withinasphereof30kpcradiusasafunctionoftimeforG3,PSPH,andMFMmodels.
TheleftpanelisforsimulationswithoutAGN-likefeedbackandtherightonewithit.
5. SummaryandConclusions out any additional feedback, our major merger models have a
very compact and massive CMC (see 3.1). This leads to an
In this paper, we discuss, in detail, the technical aspects of our unphysically high central peak of the circular velocity profile.
wet major merger simulations, where we collide two idealized Moreover, this compact and massive CMC prevents the forma-
protogalaxies. In section 2.1, we describe the initial conditions tionofabar,whileapproximatelytwothirdsofrealdiskgalax-
of our simulations. Each protogalaxy is initially spherical and ies have bars. We demonstrate that low resolution can, in some
consistsofadarkmatterhaloandagaseoushalo.Spinisadded cases, artificially mask this problem. In models with a small
to the system. In Appendix A, we demonstrate that in isolation number of particles, a very compact CMC will be relatively
theirevolutionresemblestheevolutionofdiskgalaxies;after1- quickly attenuated by two-body relaxation, and bigger gravita-
2Gyrtheytendtobesimilartointermediateredshiftdisksand tional softening will make CMC less compact from the begin-
after10Gyrtheyresemblepresentdayspirals. ning. Of course, we cannot consider decreasing the resolution
For the simulations, we use the GADGET3 code (Springel as a “solution” to the problem. We note that higher resolution,
2005) with star-formation and feedback modeled by sub-grid namelyalargernumberofparticlesandasmallersoftening,will
physics(Springel&Hernquist2003).Wedemonstratethat,with- maketheCMCproblemevenmoreacute.
Articlenumber,page10of14