DraftversionFebruary1,2008 PreprinttypesetusingLATEXstyleemulateapjv.04/03/99 IMPLEMENTING FEEDBACK IN SIMULATIONS OF GALAXY FORMATION: A SURVEY OF METHODS R. J. Thacker TheoreticalPhysicsInstitute, DepartmentofPhysics, UniversityofAlberta,Edmonton, Alberta,T6G2J1,Canada and DepartmentofPhysicsandAstronomy, UniversityofWesternOntario,London,Ontario,N6A3K7,Canada. Currentaddress DepartmentofAstronomy, UniversityofCaliforniaatBerkeley,Berkeley,CA,94720. 0 and 0 0 H. M. P. Couchman 2 DepartmentofPhysicsandAstronomy, UniversityofWesternOntario,London,Ontario,N6A3K7,Canada. n Currentaddress a DepartmentofPhysicsandAstronomy, J McMasterUniversity,1280MainSt. West, Hamilton,Ontario,L8S4M1,Canada. Draft version February 1, 2008 5 1 ABSTRACT 1 Wepresentadetailedinvestigationofanumberofdifferentapproachestomodellingfeedbackinsimu- v lations of galaxy formation. Gas-dynamic forces are evaluated using Smoothed Particle Hydrodynamics 6 (SPH). Starformationandsupernovafeedbackareincluded usingathree parametermodelwhichdeter- 7 minesthestarformationrate(SFR)normalization,feedbackenergyandlifetimeoffeedbackregions. The 2 star formation rate is calculated for all gas particles which fall within prescribed temperature, density 1 and convergentflow criteria,and for cosmologicalsimulations we also include a self-gravitycriterion for 0 gas particles to prevent star formation at high redshifts. A LagrangianSchmidt law is used to calculate 0 0 thestarformationratefromtheSPHdensity. Conversionofgastostarsisperformedwhenthestarmass / fora gasparticleexceedsacertainlimit, typicallyhalfthatofthe gasparticle. Feedbackis incorporated h by returning a precalculated amount of energy to the ISM as thermal heating. We compare the effects p ofdistributingthis energyoverthe smoothingscaleordepositingitonasingleparticle. Radiativelosses - o are prevented from heated particles by adjusting the density used in radiative cooling so that the decay r of energyoccursovera sethalf-life, orby turning off coolingcompletely andallowingfeedbackregionsa t s brief period of adiabatic expansion. We test the models on the formation of galaxies from cosmological a initial conditions and also on isolated disk galaxies. For isolated prototypes of the Milky Way and the : v dwarf galaxy NGC 6503 we find feedback has a significant effect, with some algorithms being capable i of unbinding gas from the dark matter halo (‘blow-away’). As expected feedback has a stronger effect X on the dwarf galaxy, producing significant disk evaporation and also larger feedback ‘bubbles’ for the r same parameters. In the critical-density CDM cosmological simulations, evolved to a redshift z =1, we a find the reverse to be true. Further, feedback only manages to produce a disk with a specific angular momentum value approximately twice that of the run with no feedback, the disk thus has an specific angular momentum value that is characteristic of observed elliptical galaxies. We argue that this is a result of the extreme central concentration of the dark halos in the standard CDM model and the pervasiveness of the core-halo angular momentum transport mechanism (even in light of feedback). A simulation with extremely violent feedback, relative to our fiducial models, leads to a disk that resem- bles the other simulations at z = 1 and has a specific angular momentum value that is more typical of observed disk galaxies. At later times, z = 0.5, a large amount of halo gas which does not suffer an angular momentum deficit is present, however the cooling time is too long to accrete on to the disk. We further point out that the disks formed in hierarchicalsimulations are partially a numerical artifact producedby the minimum mass scaleofthe simulationacting asa highly efficient ‘support’mechanism. Disk formationis stronglyaffectedby the treatmentofdenseregionsinthe SPHmethod. The problems inherentinthetreatmentofhighdensityregionsinSPH,inconcertwiththedifficultyofrepresentingthe hierarchicial formation process, means that realistic simulations of galaxy formation require far higher particle resolution than currently used. 1. INTRODUCTION clustering of matter is determined almost solely by gravi- tational forces, a large number of physical processes con- A detailed understanding of galaxy formation in hier- tribute to galaxy formation. Further complication is evi- archically clustering universes remains one of the primary dentinthatasignificantnumberoftheseprocessescannot goals of modern cosmology. Whereas on large scales the be modelled fromfirst principles in a simulationof galaxy 1 2 formation, the main example of this being star formation. Notably,arecentletterbyDominguez-Tenreiro,Tissera& Analytic andsemi-analytictheoriesofgalaxyformation Saiz (1998), based upon analytic work by Christodoulou, arewell-developed(seeWhite 1994,foranoverview). The Shlosman & Tohline (1995) and van den Bosch (1998), theoretical framework for studying the condensation of has shown that the inclusion of star formation may go baryons in dark matter halos was laid out by White & some way in helping to resolve the angular momentum Rees (1978) following the foundational work on hierarchi- problem. They suggestthat a secondmechanism, bar for- cal clustering by Peebles (1980, and references therein). mation,also contributesto the loss ofangularmomentum White andReesilluminatedthefactthatforbaryonscon- withinthedisk. Includingaphotoionizingbackground(Ef- densing in dark matter halos of sub-galactic (> 106 M ) stathiou 1992) appears to make the problem slightly less ⊙ size the cooling time of the gas is always sh∼orter than severe (Sommer-Larsen, Gelato & Vedel 1998) but this the free-fall collapse time. In a related paper, Fall & Ef- is the subject of debate; for a conflicting viewpoint see stathiou (1980) demonstrated that disk galaxies could be Navarro & Steinmetz (1997). formed in a hierarchical clustering model, provided that Individualsupernovaecannotbemodelledingalaxyfor- the gas maintains its angular momentum. Later work by mationsimulations. Consequently,thereisnoapriorithe- White & Frenk (1991) further developed this hierarchi- ory for how feedback energy should be distributed in the cal clustering model andidentified the cooling catastrophe simulation. A decision must be made whether feedback for CDM cosmologies. The cooling catastrophe is caused energy should be passed to the thermal or kinetic sector by the large mass fraction in small high-density halos at of the simulation. As a direct result of this ‘freedom’, a early times. The gas resident in these halos cools on the numberofdifferentalgorithmsforthedistributionoffeed- time-scaleofaMyr,thusprecipitatingmassivestarforma- back energy have been proposed (Katz 1992, Navarro & tionveryearlyoninthedevelopmentoftheCDMcosmol- White 1993, Mihos & Hernquist 1994 for example). Since ogy (at odds with the observations of our Universe). To the interstellarmediumismulti-phase(McKee&Ostriker circumvent this problem White & Frenk introduced star 1977), the feedback process should ideally be represented formation and the associated feedback from supernovae as an evaporation of molecular clouds. Note that in a andshowedthat, givenplausible assumptions,the cooling single-phasemodel, the analogueofthe cloud evaporation catastrophe can be avoided. The main deficiency in the process is thermal heating. A seminal attempt at repre- semi-analyticprogrammesisthattheycannotdescribethe senting this process has been made by Yepes et al.(1997), geometry of mergers which is exceptionally important in and recently Hultman and Pharasyn(1999) have adapted the assembly of galaxies. theYepesetal.modeltoSPH.Atthemomentitisunclear SmoothedParticleHydrodynamicsimulationsofgalaxy how well motivated these models are since (1) they must formation (Katz 1992, Navarro & White 1994, Evrard, make a number of assumptions about the physics of the Summers & Davis 1994, for example) detail the hierar- ISM and (2) they are not truly multiphase, since the dy- chical merging history, but have been limited in terms of namicsarestilltreatedusingasinglephase. Consequently, resolution. Achieving high resolution is particularly dif- given the uncertainties inherent in multiphase modeling, ficult. Since galaxy formation is affected by long range the investigation presented here explores a single phase tidalforces,thesimulationmustbelargeenoughtoinclude model. these, which in turn enforces a low mass resolution. Two As a first approximation,the star formation algorithms solutions to this problem exist; the first samples the long canbe dividedinto threegroups. Thefirstgroupcontains range fields using lower particle resolution than the main algorithms which rely upon experimental laws to derive simulation region (the multiple mass technique, Porter the star formation rate (Mihos & Hernquist 1994, for ex- 1985) while the second includes the long range fields as a ample). The second group contains those which predict pre-calculated low resolution external field (e.g. Sommer- the SFR from physical criteria (Katz 1992, for example). Larsen, Gelato & Vedel 1998). Simulations performed in The third group is comprised of algorithms which do not thismannerrepresentlengthsscalesfrom50Mpcto1kpc, attempttopredicttheSFR,butinsteadsetadensitycrite- a dynamic range of 5 104. rionforthegassothatwhenthislimitisreachedthegasis × It is well known that in SPH simulations it is compar- convertedinto stars (Gerritsen & Icke 1997,for example). atively easy to form flattened disk structures resembling In this work the first approach is taken. disk galaxies. However, the resulting gaseous structures Giventhatfeedbackoccursonsub-resolutionscales,itis are deficient in angularmomentum (see Navarro,Frenk & difficulttodecideuponthescaleoverwhichenergyshould White1995,foracomparisonofanumberofsimulationsto bereturned. However,SPHincorporatesaminimumscale observedgalaxies). The lossofangularmomentumoccurs automatically—the smoothing scale. Katz (1992)was the during the merger process as dense gas cores lose angular first to show that simply returning a specified amount of momentumto the darkmatterhalos(see Barnes1992,for energy to the ISM is ineffective: the characteristic time- an explanation of the mechanism). This is a very signifi- scaleofradiativecoolingathighdensityisfarshorterthan cant problem since a fundamental requirement of the disk the simulationdynamicaltime-scale. This is, atleastpar- formationmodel,presentedinFall&Efstathiou(1980),is tially, anotherdrawbackoftrying to model sub-resolution that the gas must maintain a similar specific angular mo- physics at the edge of the resolution scale. Further, the mentum to the dark matter to form a disk. The solution strong disparity between time-scales means that not only to this problem is widely believed to be the inclusion of are length scales sub-resolution, so is the characteristic star formation and feedback from supernovae and stellar temporal evolution. This point has been made several winds. By including these effects the gas should be kept times in relation to simulations with cooling (e.g. Katz in a more diffuse state which consequently does not suffer 1992) but it seems even more important in the context of fromthecore-haloangularmomentumtransportproblem. models with feedback. Thus, in an attempt to circumvent 3 this problem, a new thermal feedback model is presented Iftheminimumgravitationalresolutionis2ǫandthe in which the radiative losses are reduced by changing the minimum hydrodynamic resolution 4h then the min density used in the radiative cooling equation. The den- two of them match with this definition of h . min sity used is predicted from the ideal gas equation of state and then integrated forward, decaying back to the SPH Once the smoothing length of a particle reaches • density (in a prescribed period). The effect of preventing hmin, all particles within 2hmin are smoothed over. radiative losses using a brief adiabatic period of evolution Thus at this scale the code changes from an adap- for the feedback region(Gerritsen 1997)is also examined. tivetononadaptivescheme. Thisavoidsmismatched An alternative method of returning thermal energy is to gravitationalandhydrodynamicforcesatscalesclose heat an individual SPH particle (Gerritsen 1997). This to that of the resolution (Bate & Burkert 1997). method has been shown to be effective in simulations of isolateddwarfgalaxies. Thefinalmechanismconsideredis Radiative cooling is implemented using an assumed 5% onethatattemptstoincreasetheenergyinputfromSNto Z⊙ metallicity. The precise cooling table is interpolated accountforthe highradiativelosses. Mechanicalfeedback from Sutherland and Dopita (1993). Radiative cooling is boosts are not considered for the following reasons; (1) calculated in the same fashion as discussed in Thomas & parameters that are physically motivated (e.g. feedback Couchman (1992). efficiency of 10%) produce too much velocity dispersion, and (2) the method makes no account for force softening. 3. STARFORMATIONPRESCRIPTION The response of the (simulated) ISM to a single feed- 3.1. Implementation details back event is examined to gain an understanding of the Thestarformationalgorithmisbasedonthatpresented qualitativeandquantitativeperformanceofeachfeedback inMihosandHernquist(1994). Kennicutt(1998)haspre- algorithm. To determine whether the parameters of the sentedastrongargumentthattheSchmidtLaw,withstar star formation are more important than the dynamics of formation index α = 1.4 0.15, is an excellent model. It the simulation, an exploration of the parameter space of ± characterizes star formation over five decades of gas den- one of the algorithms is undertaken. Since feedback is ex- sity, although it does exhibit significant scatter. pected to have a more significant effect on dwarf systems Forcomputationalefficiency,aLagrangianversionofthe (due to the lower escape velocity) the effect of the feed- SchmidtLawisused,thatcorrespondstoastarformation back algorithms on a Milky Way prototype is contrasted index α=1.5, againstamodelofthe dwarfgalaxyNGC6503. Following this investigation, a series of cosmological simulations is dM conducted to examine whether the conclusions from the ∗ =C ρ1/2 M , (1) dt sfr g g isolated simulations hold in a cosmological environment. Particularattentionis paidto the rotationcurvesand an- where C is the star formation rate normalization, the gular momentum transport between the galaxy and halo. sfr g subscripts denotes gas and the subscript stars. As- Thestructureofthepaperisasfollows: Insection2,im- ∗ suming approximately constant volume over a time-step, portant features of the numerical technique are reviewed. both sides may divided by the volume and the standard Insections3-4,thestarformationprescriptionispresented Schmidt Law with index α = 1.5 is recovered. A value and a detailed analysis of its performance on isolated test of α = 1.5 is preferred since it leads to the square root objectspresentedconducted. Resultsfromsimulationsare of the gas density, which is numerically efficient and the reviewed in section 4, and a summary of the findings is value is within the error bounds. As written the constant given in section 6. C hasunitsofρ−1/2t,butcanbemadedimensionlessby sfr 2. NUMERICALMETHOD multiplying by (4πG)1/2. Hence, equation 1 may be writ- ten with dimensionless constants, leading to the following An explicit account of the numerical method, including form, the equationofmotionandartificialviscosityused,ispre- dM c∗M sented in Thacker et al.(1998). The code is a significant ∗ =√4πGc∗ ρ1/2 M = g, (2) developmentofthepublicallyavailableHYDRAalgorithm dt g g tff (Couchman, Thomas & Pearce 1995). The main features where c∗ is the dimensionless star formation rate (Katz of the method are summarized for clarity. Gravitational 1992). The rangeforthis parameteris reviewedinsection forces are evaluated using the adaptive Particle-Particle, 3.4. Limiting the star formation rate by the local cooling Particle-Mesh algorithm (AP3M, Couchman 1991). Hy- time-scale was not considered. The reason for this is that drodynamic evolution is calculated using SPH. Notable typically the cold dense cores which form stars are at the features of algorithm relevant to the hierarchical simula- effective temperature minimum of the simulation for non- tions follow, voidregions,namely 104 K, which corresponds to the end The neighbour smoothing is set to attempt to of the radiative cooling curve. • smooth over 52 neighbour particles, usually leading Starformationisallowedtoproceedinregionsthatsat- to a particle having between 30 and 80 neighbours. isfy the following criteria, This number of neighbors assures a stable integra- tionandreducesconcernovertherelativefluctuation 1. the gasexceedsthe densitylimitof2 10−26gcm−3 × in neighbor counts. 2. the flow is convergent,( .v<0) The minimum hydrodynamic resolution scale is set ∇ • byh =ǫ/2whereǫisthegravitationalsoftening. 3. the gas temperature is less than 3 104 K min × 4 4. the gas is partially self-gravitating, ρ >0.4ρ cloudcomplexesinanygalaxyhaveanassociatedstarfor- gas dm mation history. The first criterion associates star formation with dense There are two notable drawbacks to this algorithm. regions (regardless of the underlying dark matter struc- Firstly, since star particles are only formed when the star ture). The second criterion is included to link star forma- mass exceeds a certain threshold there is a delay in form- tion with regions that are collapsing. The third prevents ing in stars. As a consequence, prior to the star particle star formation from occurring in regions where the aver- being formed, the SPH density used in the calculation of age gas temperature is too high for star formation. The the SFRwillbe overestimated. By selectingthestarmass final criterionis particularly relevant to cosmologicalsim- tobeonehalfthatofthegasmassthisproblemisreduced, ulationssince itlimits starformationto regionswherethe but it is not removed. Secondly, until the first star parti- dynamics are at least partially determined by the baryon cle is spawned, the trajectory of the stellar component is density. Sincethemassscalesprobedbythesimulationare entirely determined by the gas dynamics. considerablylargerthanthemassscalesofself-gravitating Once a star particle is created the associated feedback cores in which star formation would occur, we use a pre- must be evaluated. A simple prescription is utilized, factor of 0.4 to help compensate for this disparity. namely that for every 100M of stars formed there is one ⊙ Representingthegrowthofthestellarcomponentofthe supernovae which contributes 1051 erg to the ISM. This simulation requires compromises. It is clearly impossible value is used in Sommer-Larsen et al.(1998), and feeds to add particles with masses dM at each time-step since back 5 1015 erg g−1 of star particle to the surrounding ∗ × this would lead to many millions of small particles being ISM. Since this value is a specific energy, a temperature addedto the simulation. Alternatively,a gasparticle may can be associated with feedback regions (see section 3.2). be viewed as having a fractional stellar mass component, Navarro&White(1993),usingaSalpeterIMF,withpower i.e. the particle is a gas-star hybrid (in the terminology law slope 1.5 and mass cut-offs at 0.1 and 40 M , derive ⊙ of Mihos & Hernquist 1994). Thus, as gas is turned into that 2 1015 erg g−1 of star particle created is fed back × stars, the stellar mass increases while the gas mass de- to the surrounding gas. The actualvalue is subject to the creases. Mihos and Hernquist take this idea to its limit IMF, but scaling between values can be achieved using a by only calculatinggas forcesusing the gas mass ofa par- singleparameterwhichwelabel e∗. Fore∗ =1 the energy ticle (gravitational forces are unaffected because of mass returnis5 1015 ergg−1,whilee∗ =0.4givestheNavarro × conservation). A drawback of this method is that it still & White value. Only brief attention is paid to changing enforces a collisional trajectory on the collisionless stellar thisparametersinceitismoreconstrainedthantheothers content of a particle. While this is a good description of inthemodel. Variationsinmetallicitycausedbythefeed- real physics, since stars take a number of galactic rota- back process (Martel & Shapiro 1998) are not considered, tions todepartfromthe gascloudinwhichthey areborn, as this involves a complicated feedback loop involving the it is in strong disagreement with what would happen to gas density and cooling rate. collisional and collisionless particles in a simulation. To 3.2. Energy feedback decide when to form stars, Katz (1992) and Steinmetz & Muller (1994, 1995) use a star formation efficiency. Once The following sections describe each of the feedback al- the star mass of the gas particle reaches a set (efficiency) gorithms considered. percentage of the gas particle mass then a star particle 3.2.1. Energy smoothing (ES) is created with that mass value. This also leads to the potential spawning of very many gas particles. As a com- The firstofthe methods is comparativelystandard: the promise, the scheme used in this work is as follows: Once totalfeedbackenergyisreturnedtothelocalgasparticles. M˙ has been evaluated, the associated mass increase over For simplicity in the following argument, it is assumed ∗ the time-step is added to the ‘star mass’ of the hybrid a tophat kernel is used to feedback the energy over the gas-star particle. Note, in this model the ‘star mass’ is a nearest neighbour particles. The number of neighbours is ‘ghost’ variable and does not affect dynamics in any way, constrainedto be between30to 80whichmaydivorcethe i.e. until the star particle is createdthe particle is treated feedbackscalefromtheminimumSPHsmoothingscaleset as entirely gaseous when calculating hydrodynamic quan- by h . This occurs in the densest regionswhere the av- min tities. Twostarparticles(ofequalmass)canbecreatedfor erage interparticle spacing becomes significantly less than everygasparticle. Thisensuresthatfeedbackeventsoccur h . This actually renders feedback more effective in min frequently since we do not have to wait for the entire gas these regions than it would be if all the energy were re- contentofaregiontobeconsumedbeforefeedbackoccurs, turnedtotheparticleswithin2h region,buttheeffect min and at the same time prevents the spawning of too many isnotsignificant. Giventhese assumptionsforstarforma- star particles. The creation of a star particle occurs when tion, and workingin units of internalenergy,it is possible thestarmassofagas-starparticlereachesonehalfthatof to evaluate the temperature increase, ∆T of a feedback the mass of the initial gas particles. The gas-star particle region as follows, mass is then decremented accordingly. The second star 2µm E M particle is created when the star mass reaches80% of half ∆T = p SN ∗, (3) 3 k N M the initial gas mass. SPH forces are calculated using the s g total mass of the hybrid particles. Note, in all the follow- which given N = 52, E = 5 1015 erg g−1, and s SN ing sections,“gasparticle” shouldbe interpretedmeaning × M /M =1/2 yields, ∗ g “gas-starhybridparticle”. Theseassumptionsyieldastar formation algorithm where each gas particle has a star M ∆T 2.4 107 ∗ 2.3 105 K. (4) formation history. This assumption is motivated because ≃ × N M ≃ × s g 5 Clearlythis boostmaybe increasedordecreasedby alter- 1 SPH density inganyoneofthevariablesESN,NsandtheratioM∗/Mg. 0.9 Cooling density KeepingE constant,butreducingN to32andincreas- SN s 0.8 ingM /M to1wouldyield∆T 7 105 K,demonstrat- ∗ g ingthesensitivityoffeedbacktot≃heS×PHsmoothingscale. 0.7 y nsit 0.6 3.2.2. Single particle feedback (SP) e de 0.5 v Alternatively,alloftheenergymaybereturnedtoasin- elati 0.4 gle SPH particle. Gerritsen (1997) has shown that this is R 0.3 aneffectiveprescriptioninsimulationsofevolvedgalaxies, yieldingaccuratemorphologyandphysicalparameters. In 0.2 this case the temperature boost is trivially seen to be 0.1 ∆T 2.4 107K. (5) 00 2 4 6 8 10 12 ≃ × Time (Myr) as N = 1 and M /M = 1. There is one minor problem Fig.1.—EvolutionofthecoolingdensityandtheSPHdensity s ∗ g in that when the gas supply is exhausted, the mechanism following asinglefeedback eventin theESnascheme. t =5 1/2 has no way of returning the energy (unless one continues Myr, and clearly the values converge within 2t . Note the 1/2 to make star particles of smaller and smaller masses). As small initial drop in the SPH density in response to the feed- a compromise the nearest SPH particle is found and the backenergy,followedbyaslowerexpansion. Theinitialcooling energy is given to this particle. densityisapproximately5%oftheSPHdensity,consistentwith thefeedbacktemperatureof2×105 Kinanambient10,000 K 3.2.3. Temperature smoothing (TS) region. The final feedback mechanism considered is one that Tocalculatethedecayrate(i.e.topredictthecoolingden- accounts for the fast radiative losses by increasing the en- sity at the time-step n+1 from that at time-step n) the ergy input. The first step is to calculate the temperature following function is used asingleparticlewouldhaveifalltheenergywerereturned to it (as in the SP model). Then this value is smoothed ρncool+dt×∆ρn(t2f/t31/2)e−0.33(tf/t1/2)3, tf ≤t1/2; overthelocalparticlesusingtheSPHkernel. Thismethod lrteheapedrEsesSteonfteasedvaabcsaatclskye)ho.figFehoxertrrteehmneeeirsfgeoyeladitnbepaduckstim(tehusaslnaenttithoianelslo,ytth5he2erstcimoaonelds- ρnco+o1l =ρρnSnS+PP1HH,−∆ρne−0.693dt/t1/2, t31t/12/2≤≤tftf<,3t1/2; (7) ing mechanism(see below)wasnotadjustedinthis model wheretistimesincethefeedbackeventoccurred,∆ρn = f since the feedbackregionsin the disk have coolingtime of ρn ρn , and dt is the time-step increment. Once a several time-steps (at least 10). In the cosmological sim- SPH − cool region passes beyond 3t , the cooling density is forced 1/2 ulations, the alternative cooling mechanisms were consid- to be equivalent to the SPH density, although usually the ered since the cooling time is only a few time-steps (Katz valuesconvergewithin2t . This comparativelycomplex 1/2 1992, Sommer-Larsen, Gelato & Vedel 1998). functionwaschosenbecauseinasimpleexponentialdecay model,thecoolingdensityincreasesbythelargestamount 3.2.4. Preventing immediate radiative energy losses immediatelyfollowingthefeedbackevent. Tohaveanyef- As was previously discussed, Katz (1992) was the first fect the cooling density must be allowed to persist at its to show that feedback energyreturned to the ISM is radi- low value for a reasonable period of time. In figure 1 the ated awayextremely quickly inhigh density regions. This two densities calculated after a single feedback event are isaresultofthecharacterisiticdynamicaltimeofthesim- compared. This cooling mechanism is denoted by a suf- ulation being far longer than that of cooling. fix na (for non-adiabatic) on the energy input acronym. Thefirstmethodforpreventingtheimmediateradiative Springel and White (in prep) have considered a similar loss of the feedback energy is to alter the density value model (Pearce,private communication). used in the radiative cooling mechanism. This change is Gerritsenallowedhis SPHparticlesto remainadiabatic motivated by Gerritsen’s (1997) tests on turning off ra- for 3 107 years, approximately the lifetime of a 8M⊙ × diative cooling in regions undergoing feedback. Assuming star. It is difficult to argue what this value should be and pressure equilibrium between the ISM phases (which is hence the t parameterspace wasexplored. Note that if 1/2 not true after a SN shell explodes but is a good starting during the decay period another feedback event occurs in point) one may derive the estimated density that the re- the localregion,the density value usedis the minimum of gion would have after the SN shell has exploded. If the the current decaying one and the new calculated value. local gasenergy is increasedby E , then the perfect gas As asecondmethodforpreventingradiativelossesGer- SN equation of state yields ritsen’s idea was utilized: make the feedback region adia- batic. Thisiseasilyachievedinourcodebyusingthesame E ρ mechanism that calculates the estimated density value. i i ρ = . (6) est E +E Providedthe estimateddensityvalueis lessthanhalfthat i SN of the local SPH value then the particle is treated as adi- Followingafeedbackeventtheestimateddensityisallowed abatic. Above this value the estimated density is used in to decay back to it’s local SPH value with a half-life t . the radiative cooling equation. This cooling mechanismis 1/2 6 Total Energy Smoothing: t*=5 Myr, no adiabatic period Total Energy Smoothing: t*=5 Myr, with adiabatic period 7 7 K) K) ure (6 ure (6 mperat5 mperat5 g Te 0 g Te 0 Lo4 100 Lo4 100 200 200 3 300 3 300 0 400 0 400 1 2 500 1 2 500 3 4 600 3 4 600 5 6 7 8 800700 Radius (pc) 5 6 7 8 800700 Radius (pc) Time (Myr) Time (Myr) Single Particle: t*=5 Myr, no adiabatic period Temperature Smoothing: t*=5 Myr, no adiabatic period 7 7 K) K) ure (6 ure (6 mperat5 mperat5 g Te 0 g Te 0 Lo4 100 Lo4 100 200 200 3 300 3 300 0 400 0 400 1 2 500 1 2 500 3 4 600 3 4 600 5 6 7 8 800700 Radius (pc) 5 6 7 8 800700 Radius (pc) Time (Myr) Time (Myr) Fig. 2.—Evolution of the cooling density and the SPH density followingasinglefeedback event inthe ESnascheme. t =5Myr, and 1/2 clearly the values converge within 2t . Note the small initial drop in the SPH density in response to the feedback energy, followed by a 1/2 slowerexpansion. Theinitialcoolingdensityisapproximately5%oftheSPHdensity,consistentwiththefeedbacktemperatureof2×105 K inanambient10,000Kregion. denoted with a suffix a. tion, so that the estimated density is significantly lower, does increase the cooling time sufficiently. However, since 3.3. Comparison of methods introducing the adiabatic phase allows the feedback en- To test each of the feedback methods and gain insight ergy to induce expansion, we prefer this method, rather into their effect on the local ISM, a single feedback event than trying to adjust the estimated density method. For was set up within a prototype isolated Milky Way galaxy. theSPrunthedensityreductionismuchhigher(sincethe The evolution of the particles within 3h of the feedback total energy, ESN, is applied to the single particle) and event were followed. The time evolution of the tempera- hence the mechanism does allow the particle to remain ture versus radius for each scheme is shown in figure 2. hot. There is little perceptible difference between the SPa and SPna profiles. 3.3.1. Qualitative discussion Both SP feedback, TS, and ESa induce noticeable ex- pansion of the feedback region. Thus after the heat input The adjusted cooling mechanism has little effect on the has been radiated away, the continued expansion intro- ESrunbecausetheestimateddensityisnotlowenoughto duces adiabatic cooling (since the temperature of the re- increase the cooling time significantly beyond the length gionfallsbelowthe 10,000Kcutoffofthe coolingcurve). of the time-step. Including a prefactor of 0.1 in the equa- 7 It is particular noticeable in the TS plot where the low Run c∗ ta e∗ Nb 1/2 step temperature plateau continues to widen quite drastically 5001 0.001 0 0.4 1999 andthisismanifestinthesimulationwiththe appearance 5002 0.003 0 0.4 2001 of a large bubble in the disk. Caution should be exer- 5003 0.01 0 0.4 1924 cised in interpreting these bubbles in any physical man- 5004 0.03 0 0.4 1977 ner, since their size is set solely by the resolution scale of 5005 0.1 0 0.4 1989 the SPH.ESaalsoproducesbubbles, butdue to the lower 5006 0.3 0 0.4 2140 temperature there is less expansion. Single particle feed- 5007 1.0 0 0.4 2087 back produces the smallest bubbles and often ejects the 7001 0.033 1. 0.4 1989 hotparticleverticallyfromthedisk. Thisoccurswithreg- 7002 0.033 5. 0.4 2000 ularity since the smoothing scale is typically larger than 7003 0.033 10. 1.0 1947 the disk thickness: if the hot particle resides close to the 7004 0.033 1. 1.0 1938 edge of the disk the pressure forces from the surrounding 7005 0.033 10. 0.4 1940 particles will be asymmetric resulting in a strong ‘push’ out of the disk. Note that this is phenomenologicallysim- Table 1: Summary of star formation parameter space simu- ilar to the mechanism by which SN gas is ejected from lations. The simulations were of a rotating cloud collapse (see disks (McKee & Ostriker 1977), although this should not Thacker et al 1998). a t = 0 denotes that feedback was re- 1/2 be over-interpreted. moved from the simulation. b The number of steps are given The only scheme which stands out in this investigation to t=1.13, thefinal point of theparameter space plot. is ESna: the cooling mechanism fails to prevent drastic radiativelosses. The remainderofthe algorithmsproduce 3.4.1. Simple collapse test an effect on both the thermal properties of the ISM and its physical distribution. To gain an understanding of the algorithms in a sim- ple collapse model (that also may be run in a short wall- clock time) the rotating cloud collapse model of Navarro and White (1993, also see Thacker et al 1998) was uti- 3.3.2. Effect of methods on time-step criterion lized. Such models actually bear little resemblance to the Since our code does not have multiple time-steps, it is hierarchical formation picture, but they do allow a fast important to discern whether one method allows longer exploration of the parameter space. time-steps than another. This is a desirable feature since For this test the self-gravity requirement was removed. it reduces the wall-clock time for simulations. Of course The reason for this is that in cosmological simulations it analgorithmwhichhasafastwall-clocktimebutproduces is virtuallyguaranteedthatthe gasin acompactdisk will poor results would never be chosen, howeverfor two algo- be self-gravitating. This is due to the low number of dark rithms with similar results this criterion provides a useful matterparticlesinthecoreofthehalorelativetothenum- parameter to choose one over the other. A comparison of ber of gas particles. SPna,TS,ESnaandthenofeedback(NF)runisdisplayed The most important parameter in the star formation in figure 3. model is the c∗ parameter since it governs the SFR nor- The simulation time versus the number of time-steps malization. Therefore, an ensemble of models was run wascomparedfordatafromtheMilkyWayprototyperuns with c∗ [0.001,1](and e∗ =0.4). The secondaryparam- ∈ (see section 4.1). The SP methods produce the shortest eter in the model, t1/2, is expected to have comparatively time-step, requiring almost double the number of time- little effect on the star formationrate (due to the low vol- steps than the NF model. The ES variants require only ume factor of regions undergoing feedback). Hence only 10% more time-steps than the run without feedback. In a range of plausible alternatives were considered, namely SPfeedbackthe accelerationfeltbythehotparticlelimits t = 1,5,10 Myr, in the ESa model. The simulation 1/2 the time-step significantly. TS also requires more time- parameters are detailed in table 1. steps than ES due to the rapid expansion of feedback 3.4.2. Results regions. Typically though, the number of time-steps re- quiredissomewhatless(20%)thanthatforsingleparticle Figure4displaysaplotofthec∗ parameterspace,show- feedback. ing SFR and c∗ versus time. Feedback was effectively turnedoffinthissimulationbyreducingtheenergyreturn efficiency to 10−7. Although a severe amount of smooth- inghashadtobeapplied(arunningaverageover40time- 3.4. Explored parameter space of the algorithm steps,followedbylinearinterpolationontothegrid)there All models exhibit dependencies on the free parameters are a number of interesting results. c∗ ande∗,correspondingtotheSFRnormalizationandef- The time at which the peak SFR occurs is almost con- ficiency of the feedback energy return. The models which stant across all values of c∗. This is an encouraging result useamodifiedcoolingformalismalsoexhibitadependence since it indicates that the time at which star formation upon t1/2, the approximate half-life of feedback regions. peaks is dictated by dynamics and not by the parame- To simplify matters, anensemble with e∗ =0.4 (the value ters of the model (at least without feedback). In fact the used in Navarro & White 1993) was run, allowing us to SFR peak time correspondsto the time when the collapse concentrate on the effect of the c∗ and t . To determine model reaches its highest density, following this moment 1/2 the effect of varying e∗, two more simulations with e∗ =1 asignificantamountofrelaxationoccursandthegashasa were run. 8 1000 900 NF SPna 800 TS ESna 700 yr) 600 M me ( 500 Ti 400 300 200 100 0 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 Time-step Fig.3.—Effectoffeedbackschemeonthetime-stepselection Fig. 5.—Dependence of the SFR on the t1/2 and e∗ param- in thesimulation. Datafrom theMilky Wayprototyperunsis eters (for the ESa algorithm). Comparing the runs with the shown. Twiceasmanytime-stepsarerequiredforSPfeedback same e∗ parameter shows that varying t from 1 to 10 has 1/2 comparedtorunswithoutfeedback. TheESaandSParunsare no noticeable effect. Conversely, changing e∗ from 0.4 to 1.0 notshown,butareapproximately5%lowerthantherespective removes thesecondary and tertiary peaks in the SFR. runs without theadiabatic period. run produces enough heating to provide a significant amountofexpansionrendering .v 0 forthe firstfeed- ∇ ≫ backregion. Notethatsincelessofthegasisusedatearly times, the SFR at later epochs is higher. 250 Insummary,whilethec∗parameterclearlysetstheSFR 200 normalization, it does not change the epoch of peak star −1SFR (M yr)o110500 trfohergemiooavntesiroiannll.tShFFeuRretvhdoeulrve,ettdohestyhtse1t/el2omwp.avroalmumeteerfrhaacstiolinttoleffeeffeedcbtaocnk 50 4. APPLICATIONTOISOLATED‘REALISTIC’MODELS 0 0 This section reports the results of applying the algo- −0.5 rithm to idealized models of mature isolated galaxies. −1 Thesemodelsarecreatedtofittheobservedparametersof suchsystems,i.e. therotationcurveanddiskscalelength. −1.5 Thecharacteristicsofeachmodelarediscussedwithinthe −2 sectionsdevotedtothem. Inthisinvestigationtherelative −2.5 1.15 Log c* 1.1 gas to stellar fraction is low, compared to the primordial 10 −3 1.05 1 ratio, and thus star formation is not as rapid as would be Time (1010 yr) Fig. 4.—Dependence of the SFR on the c∗ parameter in a expected in the early stages of the cosmological simula- tions(section5). BothmodelsweresuppliedbyDr. Fabio model with no feedback. The data for seven runs was linearly Governato. We note that they both have a sufficiently interpolated toform theplottedsurface. Thetimeofthepeak SFR moves very little with changing c∗, and almost all the high particle number (104 SPH particles)to representthe gas dynamic forces with reasonable accuracy (Steinmetz modelscanbefittedtoexponentialdecaymodelsfollowingthe & Muller 1993, Thacker et al 1998). A summary of the peak SFR epoch. simulations is presented in table 2. 4.1. Milky Way prototype lower average density. Note, however, that this is an ide- alized model with no feedback and a uniform collapse. The first prototype model is an idealized one of the Figure 5 shows the dependence of the SFR on the t1/2 Galaxy. It is desired that the model should reproduce ande∗parameters. TodetecttrendsintheSFR,arunning the measuredSFR 1M yr−1 andalsothe velocitydis- ⊙ average is shown, calculated over 40 time-steps. Clearly, persionin the disk.∼Evolvedgalaxieshavea lower relative there is little distinction between the runs with t1/2 = 1 gas content than protogalaxies. Further, because of hi- and10. This canbe attributedto the lowvolume fraction erarchical clustering, they are significantly more massive. of regions undergoing feedback. It is interesting to note Hencefeedbackisexpectedtohavelesseffectonthismodel that the SFR is more sensitive to the amount of energy thanonprotogalaxiesformedinsimulationsofhierarchical returned, dictated by the e∗ parameter, than the lifetime merging. for which this energy is allowed to persist. The line cor- responding to e∗ =1 (the standardenergy returnvalue of 4.1.1. Model Parameters 5 1015 ergg−1)doesnothavethesecondaryandtertiary The Milky Way prototype contains stars, gas and dark pe×aks in the SFR exhibited by the e∗ = 0.4 runs. This is matter. The total mass of each sector is 5 1010 M , ⊙ probably attributable to the .v criterion: the e∗ =1 × ∇ 9 SPa run (see section 3.3.1). TS produces the most sig- Run Simulation objecta feedbackb Nc N nificant disturbance in the disk, which is to be expected step SPH 1001 NGC 6503 none 3010 10240 given that it injects more energy into the ISM than the 1002 NGC 6503 SPna 3453 10240 other methods. For TS feedback, 7% of the disk gas had 1003 NGC 6503 SPa 3986 10240 been ejected (falls outside an arbitrary horizontal 6 kpc 1004 NGC 6503 TS 5345 10240 band) by t=323 Myr rising to 14% by t=506 Myr. Note 1005 NGC 6503 ESna 3453 10240 thattheamountofejectedgasiscalculatedrelativetothe 1006 NGC 6503 ESa 4535 10240 total gas in the simulation at the time of measurement. 2001 Milky Way none 387 10240 This is a decreasing function of time, but is similar for all 2002 Milky Way SPna 952 10240 simulations since the SFRs differ little. Particle ejection 2003 Milky Way SPa 943 10240 velocities, vz, were close to 300 km s−1, although some 2004 Milky Way TS 723 10240 did achieve escape velocity ( 500 km s−1 at solar ra- ≃ 2005 Milky Way ESna 424 10240 dius). Hence, while TS can project particles out of the 2006 Milky Way ESa 453 10240 halo (‘blow-away’) it preferentially leaves them bound in 6001 Cosmological SPa 4792 17165 the halo (‘blow-out’). SP feedback (both SPa and SPna) 6002 Cosmological SPna 4710 17165 has a similar evolution, but only tends to eject the single 6003 Cosmological ESa 4319 17165 heated particle during each feedback event, thus leading 6004 Cosmological ESna 4319 17165 to a lower mass-loss rate (1% of the disk gas had been 6005 Cosmological TSna 4335 17165 ejected by t=323 Myr for both versions). Particles were 6006 Cosmological TSa 4322 17165 often ejected with vz > 600 km s−1, which is larger than 6007 Cosmological NF 4341 17165 the escape velocity, leading to a proportionally stronger 6008 Cosmological NSF 4314 17165 tendencyforblow-away. ESaalsoejectsparticlesfromthe 6010 Cosmological TS 4342 17165 disk(0.4%ejectedbyt=323Myr),althoughingeneralthe particleshavelowervelocities(v 200kms−1)thanthe z ≃ particles ejected by either TS or SP. Hence, almost all of Table 2: Summary of the main simulations using realistic theejectedgasremainsboundtothesystem,i.e. ESleads models. a‘Cosmological’ refers to the object formed being de- only to blow-out. ESna does not eject particles since the rived from a cosmological simulation. bSPa=Single particle feedback regions cool sufficiently fast (0.01% ejected by adiabaticperiod,SPna=singleparticlenoadiabaticperiodbut t=323 Myr). adjusted cooling density, ESna=Total energy smoothing with TheSFRsforthreeofthesimulations(NF,TSandESa) adjusted cooling density but no adiabatic period, ESa=Total are plotted in figure 7. Most noticeable is the reduction energy smoothing with adiabatic period, TS=Temperature in the SFR produced by TS and ESa (TS is 35% lower smoothing (normal cooling), NF=no feedback, NSF=no star than no feedback at t=500 Myr, ESa is 10% lower). This formation. cFor the cosmological simulations the number of is due to three factors; (a) the ejection of matter from time-steps to z = 1.0 is given. Since the isolated simulations the disk depletes the cold gas available for star formation were not run to the same final time, the average number of (seefigure6),(b)smoothingfeedbackenergyleadstospa- time-steps per 100 Myris shown. tially extended ‘puffy’ feedback regions in the disk, which in turn leads to a lower average density, and hence lower 9 109 M⊙,and3 1011 M⊙ respectively. 11980particles SFR,(c)particlesinthe feedback regionswilltypicallybe × × were used to represent the stars, 10240 to represent the abovethe temperature threshold, which prevents star for- gas and 10240 to represent the stars. The individual par- mation which further reduces the SFR. For single particle ticlemasseswere4 106 M⊙,9 105 M⊙,and3 107 M⊙ feedback the lower mass loss rate leads to a higher SFR × × × respectively. The (stellar) radial scale length was 3.5 kpc than for the TS or ESa runs. Of the versions not plotted, and the scale height 0.6 kpc. Density and velocities were ESna resembles the no feedback run (SFR approximately assigned using the method described in Hernquist (1993). 3-4%loweronaverage),sincemostoftheenergyisrapidly The maximal radius of the dark matter halo was 85 kpc. radiatedaway. SPnaresemblestheSPasincethefeedback Theartificialviscositywasnotshear-correctedinthissim- events produce very similar effects (see section 3.3.1). ulation since the simulation was integrated through only To calculate radial profiles of the disks, an arbitrary slightly more than two rotations and the particle resolu- plane of thickness 6 kpc was centered on the disk. This tion in the object of interest is quite high. thickness ensured that the stellar bulge was contained A comparatively large softening length of 0.5 kpc was within the band. Radial binning was then performed on used,renderingtheverticalstructureofthediskpoorlyre- this data set using cylindrical bins. solved. However, this is in line with the softening lengths In figure 8, gas rotation curves are compared for the that are typically used in cosmologicalsimulations (of or- simulations at t=323 Myr. To provide a fairly accurate der2kpc). Shortersofteninglengthsallowhigherdensities depiction of the rotation curve that would be measured in the SPH, whichin turn leads to higher SFRs. The self- the rotation curve was calculated by radial averaging gravity requirement was again removed. r v/r rather than by calculating the circular velocity | × | | | from GM(<R)/R. The main drawback of this method 4.1.2. Results isthatinthe coreregions,wheretherearefew particlesin In figure 6, the gas particle distributions are shown for thebipns,themeasurementcanbecome‘noisy’. Clearly,fig- theNF,SPa,ESaandTSruns. Oftheversionsnotshown, ure 8 showsthat there is little difference between schemes ESna has a smooth disk since the feedback regions do not (a maximum of 9% at 4 scale lengths—ignoring the persist as long and the SPna disk resembles that from the 10 Fig. 6.—MorphologyofMilkyWaysimulationsatt=323Myr. z-andx-projectionsareshowntodetaildiskandgashalostructure. Both SPafeedback andTSleadtosignificantejectionofmatter fromthedisk. ESa‘inflates’thediskbutdoesnotejectasmuchmatterasTS. dispersion. Unfortunately it is difficult to relate the mea- surements made here to those of molecular clouds (Mal- hotra 1995), primarily because the mass scales are signifi- cantly different. Nonetheless, it is interesting to compare each of the separate feedback schemes. In figure 8 the radial velocity dispersion, σ , is plotted for three of the r simulations at t=323 Myr (again only considering matter withinthe6kpcband). Temperaturesmoothingproduces alargeamountofdispersionduetotheexcessiveenergyin- put(interiorto8scalelengthsitvariesbetweenbeing40% to300%higherthanothervalues). Notethatthereisadi- rect correlation between a higher velocity dispersion and a lowered measured rotation curve. This is asymmetric drift: feedbackeventsproducevelocitydispersionwhichin turnincreasestherelativedriftspeed,v ,betweenthe gas a and the local circularvelocity (Binney & Tremaine 1987). Theremainingalgorithms(SPa,SPna,ESa,ESna)exhibit similar velocity dispersions. Thus, for all but the TS al- gorithm, the bulk dynamics remain the most important Fig. 7.—SFRs for the Milky Way prototype (time-averaged factor in determining the velocity dispersion. Although over 160 time-steps to show trend). The SP variants are not the velocity dispersionpresented here is not directly com- shown sincetheirevolutionis similar tothatoftheES version patible with that of the value for local molecular clouds, plotted. TS producesasignificant (35% at t=500 Myr) reduc- it is interesting to note that measurements for the Milky tion in the SFR as compared to no feedback. ES also reduces Way suggest σ =9 1 km s−1(Malhotra 1995). vcloud ± the SFR, but has a less significant effect (10% reduction at Duetothefinitesizeofthecomputationalgrid,thecode t=500 Myrversus nofeedback) than TS. wasunabletofollowallofthesimulationstothedesiredfi- nalepoch (500Myr). This limitation wasmost noticeable in the SP simulations where ejected gas particles reached under-sampled central values). At large radii the curves the edge of the computational domain within 320 Myr. It match precisely since there are no feedback events in the is possible to remove these particles from the simulation low-density outer regions of the disk, except for the TS since they are comparatively unimportant to the remain- run where a feedback event has ejected particles to the der of the simulation. However maintaining an accurate outerregions. Comparingtotheinitialrotationcurve(not integration was considered to be more important. shown),thediskhasclearlyrelaxed,extendingbothinthe tailandtowardthecenter. Thecurveswerealsoexamined at t=507 Myr (for those simulations integratedthat long) and similar results were found with maximum differences being in the 10% range. A more telling characteristic is the gas radial velocity