Mon.Not.R.Astron.Soc.000,1–??(2011) Printed11December 2013 (MNLATEXstylefilev2.2) A hybrid SPH/N-body method for star cluster simulations D. A. Hubber1,2, R. J. Allison1,3, R. Smith4, S. P. Goodwin1 1Department of Physicsand Astronomy, University of Sheffield, HicksBuilding, Hounsfield Road, Sheffield, S3 7RH, UK 3 2School of Physics and Astronomy, Universityof Leeds, Leeds, LS2 9JT, UK 1 3Zentrumfu¨r Astronomie derUniversita¨t Heidelberg, Institut fu¨r Theoretische Astrophysik, Albert-Ueberle-Str.2, 61920 Heidelberg, Germany 0 4Departamento de Astronomia, Universidad de Concepcion, Casilla 160-C, Concepcion, Chile 2 n a August15th,2011 J 7 ABSTRACT ] M We present a new hybrid Smoothed Particle Hydrodynamics (SPH)/N-body method formodellingthecollisionalstellardynamicsofyoungclustersinalivegasbackground. I . By derivingthe equations ofmotionfromLagrangianmechanics weobtaina formally h conservative combined SPH-N-body scheme. The SPH gas particles are integrated p with a 2nd order Leapfrog, and the stars with a 4th order Hermite scheme. Our new - o approachisintendedtobridgethedivide betweenthedetailed,butexpensive,fullhy- r t drodynamical simulations of star formation, and pure N-body simulations of gas-free s star clusters. We have implemented this hybrid approach in the SPH code SEREN a [ (Hubber et al. 2011) and perform a series of simple tests to demonstrate the fidelity of the algorithm and its conservation properties. We investigate and present resolu- 1 tion criteria to adequately resolve the density field and to prevent strong numerical v 2 scatteringeffects.Future developmentswillinclude amore sophisticatedtreatmentof 9 binaries. 2 1 Key words: methods: numerical, N-body simulations - hydrodynamics - stellar dy- . namics 1 0 3 1 1 INTRODUCTION : v i Theformation anddissolution ofyoungstellarclustersisanimportant,butcomplexproblemthatrequirescomputersimula- X tionstoexploreindetail.Starsformrapidlyfromturbulentmoleculargas,mostofteninclustersoftenstothousandsofstars r (e.g. Elmegreen 2000; Lada & Lada 2003; McKee & Ostriker 2007; Klessen et al. 2009). Therefore the early phases of the a dynamical evolution of star clusters involves young stars moving through a significant and dynamic gaseous background. To understandthecompleteprocessofclusterformation andevolutionrequiresacombinationofself-gravitating hydrodynamics for the gas from which stars and planets form, and the gravitational N-body dynamics of (multiple) stars and planets once they haveformed. Previously, detailed hydrodynamics and accurate N-body dynamics have been separated. Hydrodynamical simulations havebeenusedtosimulatetheturbulentgasdynamicsleadingtofragmentationandstarformation,whileN-bodysimulations tendtofollow thelategas-freestagesofstarclusterevolution.However,thereisaverysignifcant andimportantphaseinthe life of a star cluster in which stellar dynamics within a gas background is vitally important. In the ‘gas-rich’ phase, which occursaround1–5Myr1 thestarsareinteractingdynamically inalivegasbackground.Thestarformation processtendsto produce binary and multiple systems in complex hiearchical structures. Dynamical interactions between single and multiple systems during the subsequent few Myr changes the binary properties of the stars as well as the structure and dynamics of the whole cluster (see Allison et al. 2009; Goodwin 2010, and references therein). Therefore, the binary properties of stars released into the field after gas expulsion will depend on stellar dynamics during the gas-rich phase (see also Kroupa 1995). Inaddition,theearlystagesofplanetformation will occurduringthisgas-rich phaseandinteractionsmayseriously alterthe architectureandpropertiesofplanetarysystems(Parker & Quanz2011).Accurateobservationsofthebinaryanddynamical 1 Wenotethatstarformationisnotaninstantaneousprocessandthatstarsarestillformingwhilstothersareinteractinganddynamically evolving. However, we make a rough first approximation that most stars form in the first Myr, they then evolve in a dynamical gas backgroundwhichisexpelledat5Myr. 2 Hubber, Allison, Smith & Goodwin properties of clusters are usually only available once gas is expelled (especially those to be provided by Gaia), which means they will havebeen altered by dynamical evolution in the gas-rich phase. In hydrodynamical simulations, we generally replace the dense collapse phase of gas into stars with sink particles (see Bate et al.1995forSPH;Krumholz et al.2004forAMRimplementations,seealsoFederrath et al.2010).Thesesinkparticles can represent individual stars if their sizes are . 1 AU (e.g. Bate et al. 2003; Goodwin et al. 2004; Bate 2009; Offneret al. 2009)orlargerregionsperhapscontainingprimordialmultiplesystemswhichcannotberesolvediftheirsizesare&10AU(e.g. Bonnell et al.2004;Smith et al.2009;Jappsen et al.2005).Sinkshavethehugeadvantageofallowingdense,computationally expensive,andunresolvableregionstobe‘compressed’intoaparticlewhichcaninteractwiththesurroundinggasandaccrete from it. However, sink particles are not point-like N-body particles as, even if each sink represents a single star, (a) their gravity is softened, and (b) sinks accrete from thesurrounding gas. PureN-body simulations of stellar systems havea longhistory (see Aarseth 2003; Heggie & Hut2003),but most ignore theearly gas-rich phasesofastarcluster’s life. Theusualway toincludegas andmodelthegas-rich phaseistointroducean externalpotential,whichisoftenasimplePlummerorKingmodel(e.g.Lada et al.1984;Goodwin1997;Baumgardt & Kroupa 2007;Moeckel & Clarke2011;Smith et al.2011),althoughseeGeyer & Burkert(2001)usedsoftened ‘starparticles’ inalive gas background. Generally, the external potential is allowed to vary with time, e.g. to model the expulsion of gas from a cluster. However,theuseof asimple analytic externalpotential tomodel thegas (which is often themajority of themass in an embedded cluster) is clearly a vast over-simplification. In this paper we introduce a new hybrid N-body/Smoothed Particle Hydrodynamics (SPH) algorithm that has been implemented in the SPH code SEREN (Hubberet al. 2011). The stellar dynamics are computed with a 4th-order integrator allowing the details of N-body interactions between stars to be followed. SPH gas particles are used to represent a live backgroundgaspotentialinwhichthestellardynamicsismodelled.Weemphasisethatthishybridmethodisnotareplacement forfullyself-consistent,high-resolutionstarformationsimulationsordetailedN-bodysimulations.Rather,itrepresentsafast way of exploring stellar dynamics in a live background potential which can be used to perform large suites of simulations to explore large parameter space, or to inform theinitial conditions of pureN-body simulations of thepost-gas phase. In Section 2, we introduce the hydrodynamical and N-body methods used and how they are combined algorithmically. In Section 3, we present a number of simple tests to demonstrate the accuracy and robustness of our method. In Section 4, we discuss various important caveats of our method, in particular understanding resolution effects, and also discuss possible astrophysical problems that can beexplored with this code. 2 NUMERICAL METHOD Self-gravitating hydrodynamical simulations in astrophysics are usually modelled using either a Lagrangian, particle-based approach such as Smoothed Particle Hydrodynamics (Lucy 1977; Gingold & Monaghan 1977), or an Eulerian, grid-based approach such as Adaptive Mesh Refinement Hydrodynamics (Berger & Colella 1989). Whereas SPH derives interaction terms by computing particle-particle force terms and integrating the motion of each particle individually, grid codes operate bycomputingfluxesacrossneighbouringgridcells.SinceN-bodycodesalsoworkbycomputingforcesandintegratingpositions and velocities, SPH is the most natural hydrodynamical method to merge directly with N-body dynamics as the particle- natureof the gas and stars are easily compatible making it straight-forward to derivethe coupling force terms and to merge theirindividual integration schemes. Weuseaconservativeself-gravitatingSPHformulation(Price & Monaghan2007)tomodelthegasdynamicsandinclude the star particles within the SPH formulation as a special type of SPH particle, rather than external N-body particles. Following most modernconservativeSPHschemes(e.g.Springel & Hernquist2002;Price & Monaghan 2007),thesmoothing length of a gas particle i is set by therelation, m 1/3 h = η i , (1) i ρ (cid:18) i (cid:19) and theSPH gas density is given by N ρ = m W(r ,h ). (2) i j ij i j=1 X where r , h , m , ρ are the position, smoothing length, mass and density of particle i respectively, r ≡ r −r , W is the i i i i ij i j SPH smoothing kernel and η is a dimensionless number that controls the mean number of neighbours (usually set to 1.2 to have ∼ 60 neighbours). Since h and ρ depend on each other, we must iterate between Equations 1 and 2 in order to i i reachaconsistentsolution.Incontrast,thestarparticleshaveaconstantsmoothinglengthwhichrepresentsthegravitational softening length to prevent violent 2-body collisions with other stars, in place of using more complicated algorithms such as regularisation (See Aarseth 2003, for a description of common N-body techniques). In order to reduce ‘scattering’ during A hybrid SPH/N-body method for star cluster simulations 3 star-gas interactions, we use the mean-smoothing length approach (See Appendix A of Price & Monaghan 2007), to keep star-gas interactions as smooth as possible. We can now formulate the Lagrangian of the system containing all interaction termsand thenderivetheequations of motion via theEuler-Lagrange Equations. This simple approach allows us todevelop a conservative scheme which in principle can be integrated to arbitrary accuracy (i.e. if direct summation of gravitational forces and a constant, global timestep is used). Due to the larger energy errors often produced by N-body encounters, we use a higher-order Hermite integration scheme (Makino & Aarseth 1992) to integrate star particles, and a simpler 2nd-order Leapfrog kick-drift-kickscheme to integrate thegas particles motion. 2.1 Gravitational force softening in SPH The gravitational force softening between SPH particles can be derived in a number of ways (e.g. Plummer softening, See Dehnen2001).However,ithasbeensuggestedbyBate & Burkert(1997)thatitissafesttousetheSPHkernelitselftoderive thesofteningtermstopreventartificialgravitationalfragmentation.TheyshowedthatforgascondensationswheretheJeans lengthwasoforderthesmoothinglengthorgreaterthan,thenethydrodynamicalforceisstrongerthanthenetgravitational force from all neighbouring particles, thereby suppressing or even reversing the collapse of the condensation and preventing fragmentation. We therefore derive the softening terms from the SPH kernel following the method and nomenclature of Price & Monaghan (2007). First, we consider the case of uniform smoothing length. The gravitational potential at the position of particle i due to a distribution of SPH particles is given by N Φ = G m φ(r ,h), (3) a b ab Xb=1 where φ is the gravitational softening kernel (Price & Monaghan 2007) and h is the smoothing length of all SPH particles. We note that Eqn. 3 requires the softening kernel to be a negative quantity. The potential is related to the density field by Poisson’s Equation, ∇2Φ(r) = 4πGρ(r) (4) where the SPH density defined at r is given by Eqn. 2. This allows us to directly relate the softening kernel to the SPH a smoothing kernel for a consistent formulation of self-gravity in SPH. Substituting Equations 2 & 3 into Poisson’s Equation, we obtain 1 ∂ ∂φ W(r,h) = r2 (r,h) . (5) 4πr2∂r ∂r (cid:18) (cid:19) BydirectintegrationofEqn.5withtheappropriatelimits,wecanobtainthegravitationalsofteningkernelviathegravitational force kernel,φ′ where r ∂φ 4π φ′(r,h) ≡ (r,h)= W(r′,h)r′2dr′, (6) ∂r r2 Z0 and r r Rh 1 φ(r,h) = 4π − W(r′,h)r′2dr′+ W(r′,h)r′dr′− W(r′,h)r′dr′ . (7)  r  Z0 Z0 Z0   where R is thecompact support of thekernel (e.g. R=2 for theM4-kernel). When using variable smoothing lengths, we have two choices for symmetrising the gravitational interaction; (a) use the average of thetwo softening kernels, or (b) use themean smoothing length in thesoftening kernel,i.e. N φ(r ,h )+φ(r ,h ) N Φ =G m ab a ab b , or Φ =G m φ(r ,h ), (8) a b 2 a b ab ab Xb=1 Xb=1 where h ≡ 1(h +h ). Price & Monaghan (2007) advocate using the average softening kernel approach since it requires ab 2 a b less loops over all particles than the mean smoothing length approach. When using only SPH particles where the smoothing length is determined by Eqn. 1, there is little difference in the results since the two methods give similar potentials and forces. However, when we include star particles which can havean arbitrary small smoothing length, therecan be very large discrepancies between the two methods. For example, consider a ‘collision’ between a gas particle and an SPH particle, i.e. where the two particle lie at almost the same position in space. The mean softening kernel approach has two terms, one which will be quite small due to the smoothing length of the gas particle, and the second term using the star smoothing length which can become very large producing a corresponding large scattering force. The mean smoothing length approach 4 Hubber, Allison, Smith & Goodwin however, can never produce a large scattering force, even if the smoothing length of the star becomes zero since the average ofthetwosmoothinglengthscan neverbelessthan 1h or 1h .Furthermore,themeansmoothing lengthmethod allows us 2 i 2 j to use star particles with zero smoothing length, permitting the study of truly collisional stellar dynamics with unsoftened star-star forces, but softened star-gas interactions. Therefore, for the case of our hybrid formulation, we advocate using the mean smoothing length approach and we derive all subsequent equations of motion using this method. 2.2 Conservative SPH The SPH equations of motion can be derived from Lagrangian mechanics, resulting in a set of equations that automati- cally conserve linear momentum, to machine precision, and angular momentum and energy, both to integration error (See Springel & Hernquist 2002; Price & Monaghan 2007). Price & Monaghan (2007) derived the equations of motion for self- gravitating systems with variable smoothing lengths using Lagrangian mechanics. Following their method, we derive the equationsof motion for a set of SPH particles (with variable smoothing length given byEqns.1 & 2) and stars (with a fixed smoothing length). If we have N gas particles with labels b = 1,2,...,N and N star particles with labels i = 1,2,...,N , g g s s then byinserting all terms into theLagrangian, we obtain 1 Ng 1 Ni Ng L = m v2+ m v2− m u +L (9) 2 b b 2 i i b b GRAV Xb=1 Xi=1 Xb=1 where L is the gravitational contribution tothe Lagrangian given by GRAV G Ng Ng G Ns Ns Ng Ns L = − m m φ (h )− m m φ (h )−G m m φ (h ). (10) GRAV 2 b c bc bc 2 i j ij ij b i bi bi Xb=1 Xc=1 Xi=1 Xj=1 Xb=1 Xi=1 Wenotethatthroughoutthispaper,summationsoverallSPHgasparticlesaregivenbytheindicesb,canddandsummations over all star particles by i, j and k. The equations of motion can be obtained by inserting this Lagrangian into the Euler- Lagrange equations, d ∂L ∂L − = 0. (11) dt ∂v ∂r (cid:18) a(cid:19) a After inserting the correct terms and evaluating the algebra (See Appendix A for a full derivation), we obtain the following equation of motion for SPH gas particles Ng P ∂W P ∂W a = − m a ab(h )+ b ab(h ) a b ρ2Ω ∂r a ρ2Ω ∂r b Xb=1 (cid:20) a a a b b a (cid:21) −G Ng m φ′ (h )ˆr −G Ns m φ′ (h )ˆr − G Ng m ζ¯a+χ¯a ∂Wab(h )+ ζ¯b+χ¯b ∂Wab(h ) (12) Xb=1 b ab ab ab Xi=1 i ai ai ai 2 Xb=1 b "(cid:0) Ωa (cid:1) ∂ra a (cid:0) Ωb (cid:1) ∂ra b # where P =(γ−1)ρ u is thethermal pressure of particle a,and Ω , ζ¯ , and χ¯ are defined by a a a a a a N ∂h ∂W Ω =1− a m ab(h ). (13) a ∂ρ b ∂h a a Xb=1 N ∂h ∂φ ζ¯ = a m ab(h ), (14) a ∂ρa Xb=1 b∂hab ab N ∂h ∂φ χ¯ = a m ai(h ). (15) a ∂ρa i=1 i∂hai ai X The Ω term is the familiar ‘grad-h’ correction term that appears in conservative SPH with varying smoothing length (e.g. Springel & Hernquist2002;Price & Monaghan2007).Theζ¯termisthecorrectiontermderivedbyPrice & Monaghan(2007) for gravitational interactions between gas particles in conservative SPH. We obtain an analogous correction term for the star-gasinteraction,χ¯,whichisasummationoverallneighbouringstarparticles.However,χ¯isstillincludedinasummation over all neighbouring gas particles since it is the variation in the smoothing length (which is determined by neighbouring particlepositions) thatgivesrisetothecorrection terms.Wehavesomechoiceinhowtoevolvethethermalpropertiesofthe gas particles. Wechose to evolvethe specific internal energy equation, i.e. du P Ng ∂W a = a m v · ab (16) dt ρ2Ω b ab ∂r a a Xb=1 a A hybrid SPH/N-body method for star cluster simulations 5 For thestar particles, we obtain thefollowing equations of motion for star s, Ns Ng a = −G m φ′ (h )ˆr −G m φ′ (h )ˆr . (17) s i st si si b sb sb sb Xi=1 Xb=1 We note that since only the smoothing lengths of the gas particles are allowed to vary, all correction terms derived via the Lagrangian appear in theequations of motion for thegas particles. 2.3 Coupled-integration scheme The equations of motion for both gas and star particles can be integrated either with a single integration scheme, or with two independent integration schemes in parallel. Current SPH codes typically use 2nd-order schemes, such as the Leapfrog or the Runge-Kutta-Fehlberg, whereas N-body codes use at least 4th-order schemes such as the Hermite scheme. In our implementation,wechosetousea2nd-orderLeapfrogkick-drift-kickscheme(Springel2005)tointegratetheSPHgasparticles coupledwitha4th-orderHermiteintegrationscheme(Makino & Aarseth1992)tointegratethestarparticles.Oneimportant reason for this choice is that a 4th-order Hermite scheme can be considered as the higher-order equivalent of the leapfrog scheme (See Hut et al. 1995), where the force, prediction and correction steps are all computed at the same points in the timestep for both schemes. We discuss the details of our implementation of both integration schemes, and in particular we discuss themodifications to the4th-order Hermitescheme to include SPH smoothing. 2.3.1 2nd-order SPH leapfrog integration scheme We integrate the SPH particles using a 2nd-order Leapfrog kick-drift-kick integration scheme. A traditional leapfrog works byadvancing thepositions and velocities half-a-step apart, i.e. vn+1/2 = vn−1/2+an∆t, (18) a a a rn+1 = rn+vn+1/2∆t. (19) a a a It is possible to transform thetraditional leapfrog equationsinto a form where thepositions and velocities are both updated at theend of the step, 1 rn+1 = rn+vn∆t+ an∆t2, (20) a a a 2 a 1 vn+1 = vn+ an+an+1 ∆t. (21) a a 2 a a (cid:0) (cid:1) This form of the leapfrog (also known as the Velocity-Verlet integration scheme) has 3rd-order accuracy in integrating the positionsand2nd-orderaccuracyinintegratingthevelocities,andtherefore2nd-orderoverall.AswewillseeinSection2.3.2, it also has the useful property of having its acceleration and ‘correction’ terms calculated in sync with the corresponding termsfor the4th orderHermite schemeused for integrating thestar particles. Otherintegration schemes(e.g. the2nd-order Leapfrog drift-kick-drift) donot necessarily share this property. Thetimestepsaredeterminedbytakingtheminimumofthreeseparateconditions,theSPHCourant-Friedrichs-Lewycon- dition(Courant et al.1928),theaccelerationcondition,andtheenergyconditionwhencoolingisemployed(SeeHubberet al. 2011), i.e. γ h γ h 1/2 γ u ∆t = MIN COUR a , ACCEL a , ENERGY a (22) a ((1+1.2α)ca+(1+1.2β)ha|∇·va| (cid:18) |aa|+ǫ (cid:19) |u˙ ·a|+ǫ ) wheretheγ termsaredimensionlesstimestepmultipliersandǫ≪1isasmallnumberusedtopreventadivide-by-zerointhe case of |a |=0 or |u˙ |=0. a a 2.3.2 4th-order N-body Hermite integration scheme We use a standard 4th-order Hermite integrator (Makino & Aarseth 1992) modified by including the same SPH softening schemeasthegasparticles(Hubberet al.2011)tointegratethemotionofthestarparticles.IntheHermitescheme,thefirst time derivativeof the acceleration (often referred to as the‘jerk’) must be calculated explicitly from Equation 17. By taking 6 Hubber, Allison, Smith & Goodwin 100 Force, Monopole 10-1 ForcJee,r kQ, umaodnrouppoollee 10-2 εor, 10-3 al err 10-4 n o acti 10-5 MS fr 10-6 R 10-7 10-8 10-9 0.01 0.1 1 θ max Figure1.TheRMSfractionalerrorusingtreegravityof(a)gravitionalforcewithmonopoletermsonly(bluediamonds),(b)gravitational forceincludingthequadrupoleterms(solidblackcircles),and(c)gravitational‘jerk’withmonopoletermsonly(redtriangles)usingthe geometricMACasafunctionofopeningangle,θ . MAX thetime derivativeof Equation 17 and using Equation 5, we obtain for thejerk term, N m φ′(r ,h ) N m (r ·v )φ′(r ,h ) N m (r ·v )W(r ,h ) a˙n = −G i si si v + 3G i si si si si r − 4πG i si si si si r s |r | si |r |3 si |r |2 si i=1 si i=1 si i=1 si X X X N m φ′(r ,h ) N m (r ·v )φ′(r ,h ) N m (r ·v )W(r ,h ) −G b sb sb v + 3G b sb sb sb sb r − 4πG b sb sb si sb r . (23) |r | sb |r |3 sb |r |2 sb Xb=1 sb Xb=1 sb Xb=1 sb Sincethejerkcanbecomputedfromasinglesumoverallparticles,wecancomputeitexplicitlyatthesametimeascomputing the accelerations. Once both terms are computed, we can calculate the predicted positions and velocities of the stars at the end of the steps, i.e. 1 1 rn+1 = rn+vn∆t+ an∆t2+ a˙n∆t3, (24) s s s 2 s 6 s 1 vn+1 = vn+an∆t+ a˙n∆t2. (25) s s s 2 s The acceleration and jerk are then recomputed at the end of the step, an+1 and a˙n+1. This then allows us to compute the higher ordertime derivatives, 2 −3(an−an+1)−(2a˙n+a˙n+1)∆t a¨n = s s s s , (26) s ∆t2 (cid:0) (cid:1) .a..n = 6 2(ans −ans+1)+(a˙ns +a˙ns+1)∆t . (27) s ∆t3 (cid:0) (cid:1) Finally, we apply the correction step where thehigher-order terms are added to theposition and velocity vectors, i.e. 1 1 ... rn+1 = rn+1+ a¨n∆t4+ an∆t5 (28) s s 24 s 120 s vn+1 = vn+1+ 1a¨n∆t3+ 1 .a..n∆t4. (29) s s 6 s 24 s Wecompute theN-body timesteps using theAarseth timestep criterion (Aarseth 2003), |a ||a¨ |+|a˙ |2 ∆ts = γss|a˙ss||.a..ss|+|a¨ss|2 . (30) where γ is the timestep multiplier for stars. For the very first timestep, we must compute the 2nd and 3rd time derivatives s explicitly (Aarseth 2003) since we do not yet have information on the 2nd and 3rd derivatives (since they are only first computed at theend of thefirst timestep). Hereafter, we use Eqns. 26 & 27 for computing these derivatives. 2.4 Calculating gravitational terms SEREN (Hubberet al. 2011) uses a Barnes-Hut gravity tree (Barnes & Hut 1986) for computing gravitational forces for all self-gravitating gas particles. We use the same tree for computing the gravitational forces due to all gas particles for both A hybrid SPH/N-body method for star cluster simulations 7 SPH gas particles, and star particles. SERENhas a varietyof tree-openingcriteria that can be selected at compilation-time. Formost simulations in this paper,we usetheEigenvalue multipole-acceptance criterion (Hubberet al. 2011) becauseit has better force error control, and therefore ultimately better energy error control, than the standard geometric opening-angle criterion often used. For nearby SPH particles, which require direct computation of the gravitational acceleration, we also computethejerktermwhencomputingtheacceleration ofstarparticles.Fortreecells,wecomputethejerkcontributiondue tothecentreofmassofthecell;however,weignorethecontributionduetothequadrupolemomenttermsforsimplicity.For forces due to star particles, we sum all the contributions for the gravitational acceleration (and jerk) directly without using a gravity tree. In Figure 1, we plot the RMS force (monopole and quadrupole) and jerk (monopole only) errors using the geometric opening-angle criterion (in order to plot both monopole and quadrupole errors) for stars in a star-gas Plummer sphere (as discussed in Section 3.2). We can see that the jerk error scales in the same way as the monopole-only force error. The force errorusingquadrupolemomentsscalesmuchbetterthanbothmonopoleerrors,asexpected.Inordertoallowhigh-accuracyin thecalculationofthejerkwhilestillusingthetree,weusetwodifferenttreeopeningcriteria;oneforSPHgasparticles(which do not need to calculate the jerk), and a stricter one for N-body particles that use the tree. It should be noted that Figure 1 represents an upper limit to the expected jerk error. The dominant contribution to the jerk will be from close encounters with otherstars, which is computed exactly. 2.5 Block timestepping SEREN(Hubberet al.2011)usesastandardblock-timesteppingschemeusedinmanyN-bodyandSPHcodes.Thetimesteps ofallgasandstarparticles arerestrictedtobeing∆t=∆t /2n wheren=0,1,2,3,4,..isapositiveinteger.Allparticles MAX andstarsthereforeareallexactlysynchronisedonthelongesttimestep,whenthetimesteplevelstructureisrecomputed.Inthe standard block-timestepping scheme, particle timesteps are only computed once their current timestep has been completed. At the end of the step, particles are allowed to move to any lower timestep (higher n), or are allowed to move up one level (lower n) provided the new higher level is synchronised with the old level. This approach means particles can rapidly and immediatelydroptoshorttimestepswhenrequired,butareonlyallowedtorisebackupslowlytopreventtimestepsoscillating up and down too frequently. SERENalsocontainstheneighbour-timestepmonitoringprocedureofSaitoh & Makino(2009)topreventlargetimestep differences between SPH neighbours generating large energy, momentum and angular momentum errors. Although this algo- rithmisnotnecessarilyneededforthetestspresentedinthispaper,itwillalmostcertainlyberequiredforfutureapplications where feedback processes can suddenly generate large discontinuities in density and temperature, resulting in large timestep disparities. 3 TESTS WepresentasmallsuiteofnumericaltestsdemonstratingtheaccuracyofourhybridSPH/N-bodyapproach.Testsusingthe SPH and N-body components independently were presented by Hubberet al. (2011). The SPH component was tested using avarietyofshock-tube,Kelvin-Helmholtzinstability andgravitytests,aswellastestsofthetreeandtheerror-scaling ofthe code. The N-body component was tested with some 3-body examples that had known solutions. In this paper, we present tests of the combined scheme to demonstrate that the conservative equations of motion derived in Section 2.2 are correct, andthattheschemedoesnotexhibitanyunexpectednumericaleffects.Weperformasimpletestinvestigatingthescattering between gas and star particles. Using star-gas Plummer spheres, we test theexpected error scalings of theSPH and N-body componentsand thecombined scheme,as well as thestability of thestar-gas Plummer spheres.Finally, weperform asimple testof collidingstar-gas Plummerspheres.Inalltestsweusean adiabaticequationof stateinwhichheatingandcooling are onlyduetoPdV work byexpansion orcontraction,thusenablingustotesttheenergyconservation ofthecode.Wework in dimensionless unitsthroughout, where G=1. 3.1 Star-gas particle scattering OurhybridSPH/N-body method enables investigation of stellar dynamics within a gas potential that may be time-evolving and irregular. For this goal, we must clearly understand the origin and impact of any numerical effects on the results of our simulations. One important difference between gas-only interactions and star-gas interactions in this scheme is that star particles are allowed to ‘pass through’ SPH gas particles, whereas artificial viscosity will prevent gas particle penetration by formingashock.Eventhoughthegravitationalinteractionsbetweenstarsandgasaresmoothed,thestarsarestillinteracting with agravitational fielddefinedbydiscretepointsand thuscan bedeflectedbythosepoints.Therefore, gas particles can in 8 Hubber, Allison, Smith & Goodwin principlescatterstarparticlessignificantlyiftheresolutionistoocoarse.Thistestisdesignedtoinvestigatehowsignificantly gas particles can scatter star particles and to help determine resolution criteria to prevent significant numerical scattering. For point particles obeying Newton’s gravitational law, the scattering angle of a star of mass m due to a hyperbolic s encounterwith a gas particle of mass m in thecentre-of-mass frame is given by g 2m G(m +m ) 2Gm ∆θD = m +gm tan−1 bgv2 s ≈ bv2g (31) g s (cid:18) (cid:19) (e.g. Binney & Tremaine 2008) where b is the impact parameter, v is the relative tangential velocity at infinity, and the approximationisforsmalldeflections.Forsmoothed gravity,wewouldexpectthatthenetdeflectionanglewould bereduced bysmoothing,andforthedependenceonbtobefundamentallyalteredsincethegravitationalforceforneighbouringparticles tendsto zero as the distance becomes zero. For theM4-kernel, thegravitational force reaches a maximum at around |∆r|∼ 0.8h (see Figure 1 of Price & Monaghan 2007) and then decreases to zero (instead of ∝1/r2). Therefore, an approximation to the maximum possible scattering angle can be obtained by setting b ∼ h in Equation 312. Since the smoothing length is dependent solely on the gas particle positions, the strength of the star-gas interaction becomes solely a function of relative velocity, v. From Equation 31, we can define a critical relative velocity, vcrit, where the interaction results in a significant deflection angle, where 2Gm 1/2 vcrit = g . (32) h (cid:18) (cid:19) In this form, the deflection angle is simply ∆θD = vc2rit/v2. Therefore, a simple resolution condition that ensures star-gas scattering is negligible, i.e. ∆θD ≪ 1, is given by v2 ≫ vc2rit. This application of this criterion to various astrophysical scenarios will be discussed later in thepaper. In order to test the validity of this assertion, we perform simulations of a single star that is moving with velocity v through a periodic gas cube of side-length L and uniform density ρ where thegas velocity is fixed to zero everywhere. If the gas densityfieldisperfectly uniform (i.e. inthecontinuumlimit whereN →∞),thenthenet gravitational force duetothe g gas will be zero and the star will simply move with constant velocity without any deviation or deceleration. If the density field is not perfectly uniform, as is the case when represented by a discrete set of particles, then small deflections will alter the path of the star. In our test, the gas particles are first relaxed to a glass (See Hubberet al. 2011, for details) which is the most uniform arrangement of particles used in typical simulations. We use periodic boundary conditions combined with Ewald gravity (Hernquistet al. 1991) to produce a net zero gravitational field in the gas, with the exception of the small deviations due to the smoothed particle distribution. Fellhauer et al. (2000) performed a similar scattering test to quantify theeffects of numerical relaxation, although in the context of large-scale galaxy simulations. The initial velocity v of the star particle is varied between simulations. We test v = 14vcrit, 21vcrit, vcrit, 2vcrit, 4vcrit, 8vcrit,and16vcrit.Foreachcase,weplacetheparticleatarandompositioninthegascubeinordertoremovethesystematic effectsofusingthesameglassarrangement.Wesimulate6realisations foreachoftheselected velocitiesandtaketheaverage and standard deviation of the results. Each star particle moves along its trajectory until t = L/v, i.e. the crossing time of the cube. This ensures, in so far as is possible, that during a simulation a star encounters approximately equal number of gas particles for all choices of initial velocity. To quantify the effect of the scattering on the star particle, we measure the deflection angle, ∆θD, i.e. the angle of the particle’s trajectory relative to the initial velocity along the x-axis.The results of thesetests can be seen in Fig. 2. For v > vcrit, increasing the star particle velocity decreases the net effect of scattering due to the gas particles. In this regime, the scattering angle, ∆θD, falls as v−2, the same dependency as suggested by Eqn. 31 with b ∼ h. We notice that the net average scattering angle is lower than expected (Figure 2; dashed line) due to the effects of smoothing. We note the net deflection is the accumulation of several (∼ 10) deflections, not necessarily in the same direction hence it will ‘random- walk’ with each deflection. The error bars are due to the combination of the random-walk errors and the impact parameter dependence(which is reduced, but not eliminated). However,forv6vcrit,extremelystrongscatteringoccursanddominatesthedynamicsofthestar.Atthesevelocities,the power-law relation between the scattering angle and velocity is broken since the interaction is now effectively a parabolic or elliptical interaction instead of a hyperbolic encounter. The final deflection angle is almost random due to the strong nature of the perturbations. The kinetic energy of the star is smaller than the gravitational potential energy, even accounting for smoothing, and therefore a star can in principle become trapped and bound to individual SPH particles. We note that this wouldnotnecessarilyhappeninamorerealisticastronomicalsimulationbecausethegasparticlesinthistestarestatic(sothe starmovesasatestparticle),andthereforecannotbescatteredoffthestarthemselves.Atsuchvelocities,thecoarsenessofthe 2 Amorerigorousderivationinvolvingtheformof theforcekernel wouldrevealthe exact dependence ofthescattering angleonband v. However, setting b= Ch gives the same qualitative behaviour averaged over all impact parameters where C is some multiplicative constantoforderunity. A hybrid SPH/N-body method for star cluster simulations 9 Figure 2. Deflection angle ∆θD versus star particle velocity (normalised by the critical velocity). As v increases above vcrit, the star particleisincreasinglylessscattered.Thescatteringangleisroughlyproportionaltov−2.However,forv6vcrit,thereisaclearchange inbehaviour.Scatteringisstrong,andthescatteringanglenolongerfollowsthesamepower-lawtrendwithparticlevelocity. SPH particle distribution clearly introduces potentially serious numerical effects which could corrupt any hybrid simulation. Infuturesectionsinthispaper,wediscusspossiblegasresolutioncriteriaasameansofavoidingunphysicalscatteringeffects. 3.2 Stability of star-gas Plummer spheres ThePlummersphere(Plummer1911)iscommonly usedinstellar dynamicssimulationsasitisdescribedbysimple,analytic formula. For a Plummer sphereof mass, M, and characteristic Plummer radius, a,the density distribution, ρ(r),is given by 3M r2 −5/2 ρ(r) = 1+ , (33) 4πa3 a2 (cid:18) (cid:19) and the1-D velocity dispersion, σ(r),is given by GM r2 −1/2 σ2(r) = 1+ . (34) 6a a2 (cid:18) (cid:19) We simulate Plummer spheres that are purely stars, purely gas, or a mixture of stars and gas3. For setting up the stellar component of ourclusters, we use themethod outlined by Aarseth et al. (1974). For the gas distributions, the Plummer model corresponds to a n = 5 polytrope. A polytrope is a self-gravitating gas whose equation of state obeys the form P =Kρ1+n1 and whose density structure is a solution of the Lane-Emden equation (SeeChandrasekhar1939,foran in-depthdescription of polytropes).Forthen=5polytrope,theradial densitydistribution is of the same form as Equation 33. Instead of setting a velocity field complimentary to the density field to support against collapse, a polytrope is supported by a thermal pressure gradient. The thermal energy of the gas is related to the velocity dispersionofthegasbyequatingittothesoundspeedandthenconvertingtospecificinternalenergybyu(r)=σ2(r)/(γ−1) whereγ istheratioofspecificheatsofthegas.Wenotethatthegasitselfinoursimulationdoesnotneedtoobeyapolytropic equation of state, only that the thermal energy distribution of the gas is set-up to mimic the pressure distribution of the equilibriumpolytropeandtherefore remain inhydrostaticbalance. Thegas respondsadiabatically andtherefore canheat by contractionorcool byexpansionasitsettlesorismovedaroundbythepotentialofthestars.Theinitialpositionsareset-up usingthemethodofAarseth et al.(1974),butthethermalenergiesareset usingtheaboveequationandtheinitial velocities are set to zero. Wesimulatetheevolutionof(a)agas-onlyn=5polytrope,(b)astar-onlyPlummersphere,and(c)a50-50mixture(by mass)ofastar-Plummersphereandan=5gaspolytrope.Gas-onlysimulationsareconductedwithN =5,000gasparticles g oftotalmassM ,andstar-only simulationswith N =500 starparticles oftotalmass M .Foreach case,M =M +M =1 g s s s g and a = 1 using dimensionless units (where G = 1). Mixed star-gas simulations have either either N = 100 or 500 star s 3 Note that Plummer spheres are formally infinite in extent. In practice we truncate our Plummer spheres at a radius of 20a. This meansthattheyarenotinexactequilibrium;however thishasanegligibleeffectontheirevolution 10 Hubber, Allison, Smith & Goodwin 10-4 (a) SPHO g(∆ast2 p) ascrtaiclilnegs 10-5 (b) N-bodyO s(t∆atr4 )p ascrtaiclilnegs 10-6 10-7 10-5 10-8 E E ∆E/ ∆E/ 10-9 10-6 10-10 10-11 10-12 10-7 10-13 10-2 10-1 10-2 10-1 τ τ 100 (c) Hybrid (χ correction term) Hybrid (no χ correction) 10-2 O(∆t4) scaling 10-4 E E/ 10-6 ∆ 10-8 10-10 10-1 τ Figure3.Thefractionalglobalenergyerror,∆E/E,asafunctionoftimestepmultiplier,τ for(a)SPHgasparticlesinan=5polytrope integratedwitha2nd-orderLeapfrog,(b)N-bodystarparticlesinaPlummersphereintegratedwitha4th-orderHermite,and(c)SPH gasandstarparticlesinaPlummer/Polytropecombinationwherethegasisintegrated witha2nd-orderLeapfrogandthestarswitha 4th-order Hermite, with and without the χ¯ correction term derived in Section 2.2 and Appendix A. Also plotted for guidance are the 2nd-and4th-orderscalingexpected fortheSPH,N-body andhybridsimulations. particlesandeither10×or100×asmanygasparticles.Inallcasesweuseequal-massstarparticles(m =M /N )andequal- s s s mass gas particles (m =M /N )4. The smoothing length of the stars in all simulations is h=0.0001a. Each simulation is g g g run for 40 crossing times, where we define the crossing time here to be t = a/σ(r = 0) = 6a3/GM = 2.45 code time CR units. p 3.2.1 Resolution criteria Following theideasdiscussedinSection3.1,wecandeterminetheresolution requirementsofanequilibriumPlummersphere tosignificantlyreducetheeffectsofunphysicalstar-gasscattering.LetusassumethatN gasparticlesaccount forafraction g f ofthetotalmass,i.e.Mg =fM =Ngmg whereeachgasparticlehasmassmg.Thecentralgasdensityisρ0 =3Mg/4πa3, andthecentralvelocitydispersion isσ2=GM/6a.Wecanthencalculate thecritical resolution velocityat thecentreofthe 0 Plummer sphereas 2Gm 1/2 2fGM 1/2 3 1/3 v = g = N−1/3 (35) CRIT h ηa 4π (cid:18) (cid:19) (cid:18) (cid:19) (cid:18) (cid:19) where we havesubstituted Equation 1 for h and used theabove expressions for ρ0. In order to avoid catastrophic numerical scattering of star particles off gas particles, star particles must be moving at velocities significantly larger than the critical velocity, i.e. σ0 ≫ vc. This leads to the following resolution criterion for the total gas particle number,N, in the Plummer sphereas 12f 3/2 3 1/2 N ≫ ≈15f3/2 (36) η 4π (cid:18) (cid:19) (cid:18) (cid:19) 4 Equal-massstarparticles isobviouslyasimplificationforthe purposes ofour tests. However equal-mass gas particlesshouldbeused toreduceSPHnoiseandpreventparticleclumping.