Accepted, January15,2007 PreprinttypesetusingLATEXstyleemulateapjv.10/09/06 DISSIPATIONLESS COLLAPSES IN MOND Carlo Nipoti Dept. ofAstronomy,UniversityofBologna,I-40127Bologna,Italy Pasquale Londrillo INAF-BolognaAstronomicalObservatory,I-40127Bologna,Italy Luca Ciotti 7 Dept. ofAstronomy,UniversityofBologna,I-40127Bologna,Italy 0 Accepted, January 15, 2007 0 2 ABSTRACT n Dissipationless collapses in Modified Newtonian Dynamics (MOND) are studied by using a new a particle-mesh N-body code based on our numerical MOND potential solver. We found that low J surface-density end-products have shallower inner density profile, flatter radial velocity-dispersion 5 profile, and more radially anisotropic orbital distribution than high surface-density end-products. 1 The projected density profiles of the final virialized systems are well described by Sersic profiles with 1 index m<4, down to m 2 for a deep-MOND collapse. Consistently with observations of elliptical v galaxies,∼the MOND end-∼products, if interpreted in the context of Newtonian gravity, would appear 8 to have little or no dark matter within the effective radius. However, we found impossible (under 1 4 the assumption of constant mass-to-light ratio) to simultaneously place the resulting systems on the 1 observedKormendy,Faber-JacksonandFundamentalPlanerelationsofellipticalgalaxies. Finally,the 0 simulations provide strong evidence that phase mixing is less effective in MOND than in Newtonian 7 gravity. 0 / Subject headings: gravitation — stellar dynamics — galaxies: kinematics and dynamics — galaxies: h p elliptical and lenticular, cD – galaxies: formation — methods: numerical - o r 1. INTRODUCTION throughout this paper we use t s y In Bekenstein & Milgrom’s (1984, hereafter BM) La- a µ(y)= . (5) : grangianformulationofMilgrom’s(1983)ModifiedNew- 1+y2 v Xi tonian Dynamics (MOND), the Poisson equation In the so-called ‘deep pMOND regime’ (hereafter ar ∇2φN =4πGρ (1) adM),OµN(yD)),=dyesacnridbisnogelqowua-aticocnel(e2r)atsiiomnpslyifisteesmtos (k∇φk≪ 0 for the Newtonian gravitationalpotential φ is replaced N ( φ φ)=4πGa ρ. (6) by the field equation 0 ∇· k∇ k∇ φ The source term in equation (2) can be eliminated by µ k∇ k φ =4πGρ, (2) using equation (1), giving ∇· a ∇ (cid:20) (cid:18) 0 (cid:19) (cid:21) where a0 ≃ 1.2× 10−10ms−2 is a characteristic accel- µ k∇aφk ∇φ=∇φN+S, (7) eration, ... is the standard Euclidean norm, φ is the (cid:18) 0 (cid:19) k k MOND gravitational potential produced by the density where S = curlh is a solenoidal field dependent on ρ distribution ρ, and in finite mass systems φ 0 for and in general different from zero. When S = 0 equa- ∇ → x . The MOND gravitationalfield g experienced tion(7)reducestoMilgrom’s(1983)formulationandcan k k→∞ by a test particle is besolvedexplicitly. Suchreductionispossibleforconfig- urations with spherical, cylindrical or planar symmetry, g= φ, (3) −∇ whicharespecialcasesofa moregeneralfamily ofstrat- and the function µ is such that ifications (BM; Brada & Milgrom 1995). Though the solenoidal field S has been shown to be small for some y for y 1, µ(y) ≪ (4) configurations(Brada& Milgrom1995;Ciotti, Londrillo ∼ 1 for y 1; (cid:26) ≫ &Nipoti2006,hereafterCLN),neglectingitwhensimu- Electronicaddress: carlo.nipoti@unibo.it latingtime-dependentdynamicalprocesseshasdramatic 2 Nipoti, Londrillo & Ciotti effects such as non-conservation of total linear momen- tion 4 and discussed in Section 5. tum (e.g. Felten 1984; see also Section 3.1). 2. THEN-BODYCODE Nowadays severalastronomical observationaldata ap- While most N-body codes for simulations in New- pear consistent with the MOND hypothesis (see, e.g., tonian gravity are based on the gridless multipole ex- Milgrom 2002; Sanders & McGaugh 2002). In addition, pansion treecode scheme (Barnes & Hut 1986; see also Bekenstein(2004)recentlyproposedarelativisticversion Dehnen2002),thenon-linearityoftheMONDfieldequa- of MOND (Tensor-Vector-Scalar theory, TeVeS), mak- tion(2)forcesonetoresorttoothermethods,suchasthe ing it an interesting alternative to the cold dark matter particle-meshtechnique(seeHockney&Eastwood1988). paradigm. However,dynamicalprocessesinMONDhave Inthisapproach,particlesaremovedundertheactionof beeninvestigatedverylittlesofar,mainlyduetodifficul- a gravitational field which is computed on a grid, with ties posed by the non-linearity of equation (2). Here we particle-mesh interpolation providing the link between recall the spherically symmetric simulations (in which the two representations. In our MOND particle-mesh S = 0) of gaseous collapse in MOND by Stachniewicz N-body code, we adopt a spherical grid of coordinates &Kutschera(2005)andNusser&Pointecouteau(2006). (r, ϑ, ϕ), made of N N N points, on which the Theonlygenuinethree-dimensionalMONDN-bodysim- r × ϑ × ϕ MONDfieldequationissolvedasinCLN.Particle-mesh ulations(inwhichequation[2]issolvedexactly)arethose interpolations are obtained with a quadratic spline in by Brada & Milgrom (1999, 2000), who studied the sta- eachcoordinate,while time stepping is givenby aclassi- bility of disk galaxies and the external field effect, and cal leap-frog scheme (Hockney & Eastwood 1988). The those of Tiret & Combes (2007). Other attempts to time-step ∆t is the same for all particles and is allowed study MOND dynamical processes have been conducted to vary adaptively in time. In particular, according to usingthree-dimensionalN-body codesbyarbitrarilyset- the stability criterion for the leap-frog time integration, ting S = 0: Christodoulou (1991) investigated disk sta- weadopt∆t=η/ max 2φ,whereη<0.3isadimen- bility, while Nusser (2002) and Knebe & Gibson (2004) sionless parameter. We f|o∇und|that η =∼0.1 assures good explored cosmological structure formation1. p conservation of the total energy in the Newtonian cases In this paper we present results of N-body simula- (see Section 3.1). In the present version of the code, all tions of dissipationless collapse in MOND. The simula- the computations on the particles and the particle-mesh tions were performed with an originalthree-dimensional interpolations can be split among different processors, particle-mesh N-body code, based on the numerical while the computations relative to the potential solver MOND potential solver presented in CLN, which solves are not performed in parallel. The solution of equa- equation (2) exactly. These numerical experiments are tion (2) over the grid is then the bottleneck of the sim- interesting both from a purely dynamical point of view, ulations: however, the iterative procedure on which the allowing for the first time to explore the relaxation pro- potential solver is based (see CLN) allows to adopt as cesses in MOND, and in the context of elliptical galaxy seed solution at each time step the potential previously formation. In fact, the ability of dissipationless collapse determined. at producing systems strikingly similar to real ellipti- TheMONDpotentialsolvercanalsosolvethe Poisson calsisaremarkablesuccessofNewtoniandynamics(e.g., equation(obtained by imposing µ=1 in equation 2), so van Albada 1982; Aguilar & Merritt 1990; Londrillo, Newtonian simulations can be run with the same code. Messina & Stiavelli 1991; Udry 1993; Trenti, Bertin & We exploited this property to test the code by running vanAlbada 2005;Nipoti, Londrillo,& Ciotti 2006,here- severalNewtoniansimulationsofbothequilibriumdistri- after NLC06), while there have been no indications so butions and collapses, comparing the results with those far that MOND can work as well in this respect. Here ofsimulations(startingfromthe sameinitialconditions) we study the structural and kinematical properties of performed with the FVFPS treecode (Londrillo, Nipoti the end-products of MOND simulations, and we com- & Ciotti 2003; Nipoti, Londrillo & Ciotti 2003). One of pare them with the observed scaling relations of ellipti- these tests is described in Section 4.1.2. cal galaxies: the Faber–Jackson (FJ) relation (Faber & We also verified that the code reproduces the New- Jackson 1976), the Kormendy (1977) relation, and the tonian and MOND conservation laws (see Section 3.1): Fundamental Plane (FP) relation (Djorgovski & Davis note that the conservation laws in MOND present some 1987, Dressler et al. 1987). peculiarities with respect to the Newtonian case, so we The paper is organized as follows. The main features give here a brief discussion of the subject. As already ofthenewN-bodycodearepresentedinSection2,while stressed by BM, equation (2) is obtained from a vari- Section 3 describes the set-up and the analysis of the ational principle applied to a Lagrangian with all the numerical simulations. The results are presented in Sec- required symmetries, so energy, linear and angular mo- 1CosmologicalN-bodysimulationsinthecontextofarelativistic mentum are conserved. Unfortunately, as also shown by MONDtheorysuchasTeVeShavenotbeenperformedsofar. BM, the total energy diverges even for finite mass sys- Dissipationless collapses in MOND 3 in practice 2 systems of arbitrarymass and size. Each TABLE 1 ∞ Time, velocity,andenergyunitsforNewtonian and of them is obtained by assigning specific values to the MOND(subscript n),anddMOND (subscriptd) N-body length and mass units, r∗ and M∗, in which the initial simulations. conditions are expressed. Also dMOND gravity is scale free, because a appears only as a multiplicative factor 0 t∗n=r∗3/2(GM∗)−1/2 t∗d=r∗(GM∗a0)−1/4 in equation (6), and so a simulation in dMOND gravity v∗n=(GM∗)1/2r∗−1/2 v∗d=(GM∗a0)1/4 representssystems with arbitrarymass and size (though E∗n=GM∗2r∗−1 E∗d=(Ga0)1/2M∗3/2 inprinciple the results applyonly to systems with accel- tems,thusposingacomputationalchallengetocodeval- erations much smaller than a ). MOND simulations can 0 idation. Wesolvedthisproblembycheckingthevolume- also be rescaled, but, due to the presence of the charac- limited energy balance equation teristicaccelerationa inthe non-linearfunctionµ, each 0 d k+ρφ+ a20 ||∇φ|| d3x= sMim∗uclaantinoontdbeescchriobseesnoinnldyep∞en1desynsttlyemosf,eabcehcaoutsheerr.∗ and dt 8πGF a ZV0(cid:20) (cid:18) 0 (cid:19)(cid:21) Onthebasisoftheabovediscussion,wefixthephysical 1 ∂φ µ < φ,nˆ >da, (8) units as follows (see Appendix C for a detailed descrip- 4πG ∂t ∇ Z∂V0 tion of the scaling procedure). Let the initial density which is derived in Appendix A. In equation (8) V0 is distribution be characterized by a total mass M∗ and a an arbitrary (but fixed) volume enclosing all the system characteristicradiusr∗. Werescalethefieldequationsso mass, k is the kinetic energy per unit volume, and that the dimensionless source term is the same in New- tonian, MOND and dMOND simulations. We also re- y (y) 2 µ(ξ)ξdξ, (9) quire that the Second Law of Dynamics, when cast in F ≡ Zy0 dimensionless form, is independent of the specific force where y0 is an arbitrary constant; note that only finite law considered, and this leads to fix the time unit. As a quantities are involved. Another important relation be- result,NewtonianandMONDsimulationshavethesame tween global quantities for a system at equilibrium (in time unit t∗n = r∗3/2(GM∗)−1/2, while the natural time MOND as in Newtonian gravity) is the virial theorem unit in dMOND simulations is t∗d = r∗(GM∗a0)−1/4. Note that MOND simulations are characterized by the 2K+W =0, (10) dimensionless parameter κ = GM∗/r∗2a0, and scaling of where K is the total kinetic energy and W = TrWij is aspecific simulationisallowedprovidedthe value ofκ is the trace of the Chandrasekhar potential energy tensor maintained constant. So, simulations with lower κ val- ∂φ(x) ues describe lower surface-density, weaker acceleration W ρ(x)x d3x (11) ij i systems; dMOND simulations represent the limit case ≡− ∂x Z j κ 1, while Newtonian ones describe the regime with (e.g., Binney & Tremaine 1987). Note that in MOND ≪ κ 1. With the time units fixed, the corresponding ve- K + W is not the total energy, and is not conserved. ≫ locity and energy units are v∗n r∗/t∗n, v∗d r∗/t∗d, However, W is conserved in the limit of dMOND, be- E∗n =M∗v∗2n, and E∗d =M∗v∗2d≡(see Table 1 fo≡r a sum- ing W = (2/3) Ga M3 for all systems of finite total 0 ∗ mary). − mass M∗ (see Appendix B for the proof). As a con- p sequence, in dMOND the virial theorem writes simply 3.1. Initial conditions and analysis of the simulations σV4 = 4GM∗a0/9, where σV ≡ 2K/M∗ is the system We performed a set of five dissipationless-collapse N- virial velocity dispersion (this relation was proved for p body simulations, starting from the same phase-space dMOND spherical systems by Gerhard & Spergel 1992; configuration: theinitialparticledistributionfollowsthe seealsoMilgrom1984). Inoursimulationswealsotested Plummer (1911) spherically symmetric density distribu- thatequation(10)issatisfiedatequilibrium,andthatW tion is conserved in the dMOND case (see Sections 3.1 and 4). ρ(r)= 3M∗r∗2 , (12) 4π(r2+r2)5/2 ∗ 3. NUMERICALSIMULATIONS where M∗ is the total mass and r∗ a characteristic ra- The choice of appropriate scaling physical units is an dius. The choice of a Plummer sphere as initial con- important aspect of N-body simulations. This is espe- dition is quite artificial, and not necessarily the most cially true in the presentcase,in whichwe wantto com- realistic to reproduce initial conditions in the cosmolog- pareMONDandNewtoniansimulationshavingthesame ical context (e.g., Gunn & Gott 1972). We adopt such a initial conditions. As well known, due to the scale-free distribution to adhere to other papers dealing with col- nature of Newtonian gravity, a Newtonian N-body sim- lisionless collapse (e.g., Londrillo et al. 1991; NLC06; ulation starting from a given initial condition describes see also Section 5, in which we present the results of a 4 Nipoti, Londrillo & Ciotti Fig. 1.— Time evolution of 2K/|W|, K, W, and K+W for simulations D, M1, and N. K, W, and K+W are in units of E∗d (left column),andE∗n (central andrightcolumns). Forclarity,thetimeaxes arezoomed-inbetween 0and10. set of simulations starting from different initial condi- and a grid with N = 64, N = 16 and N = 32) are r ϑ ϕ tions). Theparticlesareatrest,sotheinitialvirialratio evolved up t = 150t . In all cases the modulus of the dyn 2K/W = 0. What is different in each simulation is center of mass position oscillates around zero with r.m.s | | the adopted gravitationalpotential, which is Newtonian <0.1r∗; similarly, the modulus of the total angular mo- in simulation N, dMOND in simulation D, and MOND m∼entum oscillates around zero2 with r.m.s. <0.02, in with acceleration ratio κ in simulations Mκ (κ=1, 2, 4). unitsofr∗M∗v∗n(simulationsMκandN)ando∼fr∗M∗v∗d For each simulation we define the dynamical time t (simulationD).K+W intheNewtoniansimulationand dyn as the time at which the virial ratio 2K/W reaches its W inthedMONDsimulationareconservedtowithin2% | | maximum value. In particular, we find tdyn 2t∗d in and 0.6%, respectively. The volume-limited energy bal- ∼ simulation D, and tdyn 2t∗n in simulations N, M1, M2 anceequation(8)isconservedwithanaccuracyof1%in ∼ andM4. Wenotethattdyn GM∗5/2(2K+W )−3/2 in MONDsimulations,independentlyoftheadoptedV0. To ∼ | | simulation N. estimate possible numerical effects, we reran one of the Following NLC06, the particles are spatially dis- MOND collapse simulations (M1) using N = 2 106, × tributed according to equation (12) and then randomly N = 80, N = 24, and N = 48: we found that the r ϑ ϕ shifted in position (up to r∗/5 in modulus). This arti- end-products of these two simulations do not differ sig- ficial, small-scale ”noise” is introduced to enhance the 2 As an experiment we also ran a simulation, with the same phase mixing at the beginning of the collapse, because initial conditions and parameter κ as M1, in which the force was the numerical noise is small, and the velocity dispersion calculated from equation (7) imposing S = 0. In this simulation is zero (see alsoSection 4.2). As such, these fluctuations thelinearandangular momentum arestronglynotconserved: for are not intended to reproduce any physical clumpiness. instance, the center of mass is already displaced by ∼ 7r∗ after All the simulations (realized with N = 106 particles, ∼30tdyn. Dissipationless collapses in MOND 5 nificantly,asfarasthe propertiesrelevanttothe present generaloverviewofthetimeevolutionofthevirialquan- work are concerned. titiesinoursimulations,postponingtoSection4.2amore The intrinsic and projected properties of the collapse detaileddescriptionofthephase-spaceevolution. Inpar- end-products are determined as in NLC06. In partic- ticular, in Fig. 1 we show the time evolution of 2K/W , | | ular, the position of the center of the system is deter- K,W, and K+W for simulations D, M1, and N. In the mined using the iterative technique described by Power diagramstime isnormalizedto t ,soplots referringto dyn et al. (2003). Following Nipoti et al. (2002),we measure different simulations are directly comparable (the values theaxisratiosc/aandb/aofthe inertiaellipsoid(where of t in time units for the five simulations are given dyn a,bandcarethemajor,intermediateandminoraxis)of in Section 3.1). In simulation N (right column) we find the final density distributions, their angle-averagedpro- thewellknownbehaviorofNewtoniandissipationlesscol- file and half-mass radius r . We fitted the final angle- lapses: 2K/W has a peak, then oscillates, and eventu- h | | averageddensityprofileswiththeγ-model(Dehnen1993; ally converges to the equilibrium value 2K/W =1; the | | Tremaine et al. 1994) total energy K +W is nicely conserved during the collapse, though it presents a secular drift, a well known ρ r4 ρ(r)= 0 c , (13) feature of time integration in N-body codes. The time rγ(r +r)4−γ c evolution of the same quantities is significantly different where the inner slope γ and the break radius r are inadMONDsimulation(leftcolumn). Inparticular,the c free parameters, and the reference density ρ is fixed virial ratio 2K/W quickly becomes close to one, but is 0 | | by the total mass M∗. The fitting radial range is still oscillating at very late times because of the oscilla- 0.06 < r/r < 10. Inordertoestimatetheimportanceof tionsofK,whileW isconstantasexpected. Asweshow h proje∼ctioneff∼ects,foreachend-productweconsiderthree inSection4.2, these oscillationsare relatedto a peculiar orthogonalprojectionsalongtheprincipalaxesofthein- behavior of the system in phase space. Finally, simula- ertia tensor, measuring the ellipticity ǫ = 1 b /a , the tionM1(centralcolumn)representsanintermediatecase e e − circularizedprojecteddensityprofileandthecircularized betweenmodels N andD: the systemstarts as dMOND, effectiveradiusR √a b (wherea andb arethema- but soon its core becomes concentrated enough to enter e e e e e ≡ jorandminorsemi-axisoftheeffectiveisodensityellipse). the Newtonian regime. After the initial phases of the We fit (over the radial range 0.1 < R/R < 10) the cir- collapse, Newtonian gravity acts effectively in damping e cularized projected density profil∼es of the∼end-products the oscillationsofthe virialratio. Overall,it is apparent with the R1/m Sersic (1968) law: how the system is in a “mixed” state, neither Newto- nian (K +W is not conserved) nor dMOND (W is not 1/m R constant). I(R)=I exp b(m) 1 , (14) e (− "(cid:18)Re(cid:19) − #) 4.1. Properties of the collapse end-products whereI I(R )andb(m) 2m 1/3+4/405m(Ciotti e e ≡ ≃ − 4.1.1. Spatial and projected density profiles &Bertin1999). Inthefittingproceduremistheonlyfree parameter, because Re and Ie are determined by their We found that all the simulated systems, once viri- measuredvaluesobtainedbyparticlecount. Inaddition, alized, are not spherically symmetric. However, while we measure the central velocity dispersion σ0, obtained the dMOND collapse end-product is triaxial (c/a 0.2, ∼ by averaging the projected velocity dispersion over the b/a 0.4), MOND and Newtonian end-products are ∼ circularized surface density profile within an aperture of oblate (c/a c/b 0.5). The ellipticity ǫ of the pro- ∼ ∼ Re/8. Some of these structural parameters are reported jected density distributions (measured for each of the in Table 2 for the five simulations described above, as principal projections) is found in the range 0.5 0.8 in − wellasforthreeadditionalsimulations,whichstartfrom D, and 0 0.5 in M1, M2, M4 and N. These values are − different initial conditions (see Section 5). consistentwiththoseobservedinrealellipticals,withthe exception of ǫ in model D (see Table 2), which would b 4. RESULTS correspond - if taken at face value - to an E8 galaxy. In Newtonian gravity, collisionless systems reach viri- These result could be just due to the procedure adopted alization through violent relaxation in few dynamical to measure the ellipticity (see Section 3.1), however we times, as predicted by the theory (Lynden-Bell 1967) find interesting that dMOND gravity could be able to and confirmed by numerical simulations (e.g. van Al- produce some system that would be unstable in New- bada 1982). Onthe other hand, due to the non linearity tonian gravity. We remark that a similar result, in the of the theory and the lack of numerical simulations, the different context of disk stability in MOND, has been detailsofrelaxationprocessesandvirializationinMOND obtained by Brada & Milgrom (1999). are much less known. Thus, before discussing the spe- Inordertodescribethe radialmassdistributionofthe cificpropertiesofthe collapseend-productswepresenta final virialized systems, we fitted their angle-averaged 6 Nipoti, Londrillo & Ciotti TABLE 2 End-product properties. κ c/a b/a γ rc/rh ma mb mc ǫa ǫb ǫc D - 0.21 0.41 0.17+0.39 0.27+0.09 2.87±0.01 2.50±0.03 2.16±0.02 0.48 0.80 0.58 −0.17 −0.02 M1 1 0.47 0.85 1.24+0.40 0.44+0.20 3.20±0.07 3.00±0.09 3.07±0.13 0.42 0.51 0.17 −0.36 −0.12 M2 2 0.48 0.92 1.45+0.26 0.53+0.19 3.38±0.08 3.24±0.08 3.28±0.12 0.49 0.45 0.08 −0.34 −0.13 M4 4 0.47 0.90 1.54+0.26 0.58+0.20 3.55±0.10 3.40±0.11 3.34±0.15 0.51 0.45 0.10 −0.32 −0.15 N - 0.45 0.91 1.69+0.13 0.74+0.14 4.21±0.07 4.35±0.08 3.96±0.13 0.48 0.55 0.12 −0.15 −0.12 D′ - 0.25 0.45 0.72+0.36 0.37+0.12 3.06±0.06 2.90±0.04 2.71±0.08 0.44 0.76 0.56 −0.50 −0.11 M′ 20 0.42 0.83 1.26+0.44 0.47+0.25 3.41±0.09 3.36±0.06 3.20±0.13 0.49 0.57 0.16 −0.40 −0.15 N′ - 0.45 0.93 1.78+0.15 0.78+0.24 4.29±0.10 4.56±0.15 4.19±0.22 0.51 0.55 0.09 −0.18 −0.18 Firstcolumn: nameofthesimulation. κ=GM∗/r∗2a0: acceleration ratio. c/aandb/a: minor-to-majorandintermediate-to-majoraxis ratios. γ,rc: best-fitγ-modelparameters. ma,mb,mc andǫa,ǫb,ǫc: best-fitSersicindicesandellipticitiesforprojectionsalongthe principalaxes. Fig. 2.—Angle-averageddensity,radialvelocity-dispersionandanisotropy-parameterprofiles(frombottomtotop)fortheend-products of simulations D,M1, and N. Dotted lines inthe bottom panels represent ρ∝r−1 profiles, whichareshown for reference. Empty circles intherightcolumnshowthecorrespondingprofilesobtained withtheFVFPStreecodefromthesameinitialconditions. Dissipationless collapses in MOND 7 Fig. 3.— Line-of-sightvelocity-dispersionprofiles(top),circularizedprojecteddensityprofilesandresidualsoftheSersicfit(bottom)of theend-products ofsimulationsD,M1,andN(squares;1-σ errorbarsarealwayssmallerthanthesymbolsize). Thedotted linesarethe best-fittingSersicmodels. 8 Nipoti, Londrillo & Ciotti density profiles with the γ-model (13) over the radial anisotropic for r>r . For each model projection we h range 0.06 < r/r < 10. The best-fit γ and r for the computed the line∼-of-sight velocity dispersion σ , con- h c los final distrib∼ution of∼each simulation are reported in Ta- sidering particles in a strip of width R /4 centered on e ble 2 together with their 1σ uncertainties (calculated the semi-major axis of the isophotal ellipse. The line-of- from ∆χ2 = 2.30 contours in the space γ r ). As sightvelocity-dispersionprofiles (for the major-axispro- c − also apparent from Fig. 2 (bottom), the Newtonian col- jection)areplottedinthetoppanelsofFig.3. TheNew- lapseproducedthesystemwiththesteepestinnerprofile tonian profile is very steep within R , while MOND and e (γ 1.7), the dMOND end-product has inner logarith- dMONDprofilesaresignificantlyflatterthere. Aswellas ∼ mic slope close to zero, while MOND collapses led to σ ,σ decreasesfordecreasingradiusintheinnerregion r los intermediate cases, with γ ranging from 1.2 (κ = 1) ofmodelD.ThekinematicalpropertiesofM2andM4are ∼ to 1.5 (κ = 4). We also note that the ratio r /r intermediate between those of M1 and of N: overall we c h ∼ (indicating the position of the knee in the density pro- find only weak dynamical non-homology among MOND file)increasessystematicallyfromdMONDtoNewtonian end-products. The empty symbols in Fig. 2 (right col- simulations. umn) refer to a test Newtonian simulation run with the The circularized projected density profiles of the end- FVFPS treecode (with 4 105 particles). The struc- × products are analyzed as described in Section 3.1. The tural and kinematical properties of the end-product of best-fit Sersic indices m , m and m (for projections this simulation are clearly in good agreementwith those a b c along the axes a, b, and c, respectively) are reported in of the end-product of simulation N (solid lines), which Table2,togetherwiththe1σuncertaintiescorresponding started from the same initial conditions. to∆χ2 =1;therelativeuncertaintiesonthebest-fitSer- 4.2. Phase-space properties of MOND collapses sicindicesareinallcasessmallerthan5percentandthe averageresiduals between the data and the fits are typi- To explore the phase-space evolution of the systems cally0.05< ∆SB <0.2,whereSB 2.5log[I(R)/Ie]. during the collapse and the following relaxation we con- The fitting∼rhadialira∼nge 0.1 < R/Re≡<−10 is comparable sider time snapshots of the particles radial velocity (vr) with or larger than the typ∼ical rang∼es spanned by ob- vs. radius as in Londrillo et al. (1991). In Fig. 4 we plot servations (e.g., see Bertin, Ciotti & Del Principe 2002). fiveofthesediagramsforsimulationsD,M1andN:each Inagreementwithpreviousinvestigations,wefoundthat plotshowsthephase-spacecoordinatesof32000particles the Newtoniancollapseproduceda systemwellfitted by randomly extracted from the corresponding simulation, the de Vaucouleurs (1948) law. MOND collapses led to and, as in Fig. 1, times are normalizedto the dynamical systems with Sersic index m < 4, down to m 2 in the time t (see Section3.1). At time t=0.5t allparti- ∼ dyn dyn case of the dMOND collapse. Figure 3 (bottom) shows cles are still collapsing in simulation N, while in MOND the circularized (major-axis) projected density profiles simulations a minority of particles have already crossed fortheend-productsofsimulationsD,M1andNtogether thecenterofmass,asrevealedbytheverticaldistribution with their best-fit Sersic laws (m =2.87, m = 3.20, and ofpointsatr 0inpanelsDandM1. Att=t (time dyn ∼ m=4.21,respectively),andthecorrespondingresiduals. ofthe peak of2K/W inthe three models), sharpshells | | Curiously,NLC06 foundthat low-msystems canbe also in phase space are present, indicating that particles are obtained in Newtonian dissipationless collapses in the movinginandoutcollectively andphasemixing has not presence of a pre-existing dark-matter halo, with Sersic taken place yet. At t = 4t is already apparent that dyn index value decreasing for increasing dark-to-luminous phase mixing is operating more efficiently in simulation mass ratio. N than in simulation M1, while there is very little phase mixing in the dMOND collapse. At significantly late 4.1.2. Kinematics times (t = 44tdyn), when the three systems are almost virialized(2K/W 1;seeFig.1),phasemixingiscom- We quantify the internal kinematics of the collapse | |∼ pleteinsimulationN,butphase-spaceshellsstillsurvive end-productsbymeasuringtheangle-averagedradialand in models M1 and D. Finally, the bottom panels show tangential components (σ and σ ) of their velocity- r t the phase-space diagrams at equilibrium (t = 150t ), dyn dispersion tensor, and the anisotropy parameter β(r) ≡ when phase mixing is completed also in the MOND and 1 0.5σ2/σ2. These quantities are shown in Fig. 2 for − t r dMOND galaxies: note that the populated region in the simulations D, M1, and N. We note that the σ pro- r (r,v ) space is significantly different in MOND and in r file decreases more steeply in the Newtonian than in the Newtoniangravity,consistently with the sharper decline MONDend-products,whileitpresentsaholeintheinner of radial velocity dispersion in the Newtonian system. regionsofthedMONDsystem. Inaddition,thedMOND Thus, our results indicate that phase mixing is more galaxy is radially anisotropic (β 0.4) even in the cen- ∼ effectiveinNewtoniangravitythaninMOND3.Itisthen tral regions, where models N and M1 are approximately isotropic (β 0.1). All systems are strongly radially 3 Ciotti, Nipoti & Londrillo(2007) found similarresults in“ad ∼ Dissipationless collapses in MOND 9 Fig. 4.—Phase-space(radial-velocityvs. radius)diagramsforsimulationsD,M1,andNatvarioustimes. vr isinunitsofv∗d(simulation D),andv∗n (simulationsM1andN;seeTable1). interestingtoestimateinphysicalunitsthephase-mixing tional systems is the energy distribution N(E) (i.e. the timescales of MOND systems. From Table 1 it follows number of particles with energy per unit mass between that t∗n 4.7(r∗/kpc)3/2(M∗/1010M⊙)−1/2Myr = E and E+dE; e.g., Binney & Tremaine 1987; Trenti & ≃ 29.8κ−3/4(M∗/1010M⊙)1/4Myr for a0 = 1.2 Bertin2005). Independently of the force law,the energy × 10−10ms−2. For example, in the case of model M1, perunitmassofaparticleorbitingatxwithspeedvina adoptingM∗ =1012M⊙ (andr∗ = GM∗/a0 34kpc), gravitationalpotentialφ(x)isE =v2/2+φ(x),andE is ≃ shells in phase space are still apparent after 8.3Gyr constant if φ is time-independent. In Newtonian gravity p ∼ ( 44t ). Simulation M1 might also be interpreted φisusuallysettozeroatinfinityforfinite-masssystems, dyn ≃ as representing a dwarf elliptical galaxy of, say, so E <0 for bound particles; in MOND all particles are M∗ = 109M⊙ (and r∗ = GM∗/a0 1.1kpc). In this bound, independently of their velocity,becauseφ is con- ≃ case44t 1.5Gyr. WeconcludethatinsomeMOND fining, and all energies are admissible. This difference dyn ∼ p systems substructures in phase space can survive for is reflected in Fig. 5, which plots the initial (top) and significantly long times. final (bottom) differential energy distributions for simu- In addition to the (r,v ) diagram, another useful di- lationsD,M1,andN.Giventhattheparticlesareatrest r agnostictoinvestigatephase-spacepropertiesofgravita- at t=0, the initial N(E) depends only on the structure of the gravitational potential, and is significantly differ- hoc”numericalsimulationsinwhichtheangularforcecomponents ent in the Newtonian and MOND cases. We also note werefrozentozero,sothattheevolutionwasdrivenbyradialforces only. In fact, while phase mixing is less effective both in MOND that N(E) is basically the same in models D and M1 at andinNewtoniansimulationswithrespecttothesimulationshere t =0, because model M1 is initially in dMOND regime. reported,thephasemixingtimescaleinMONDisstillconsiderably In accordance with previous studies, in the Newtonian longerthaninNewtoniangravity. case the final differential N(E) is well representedby an 10 Nipoti, Londrillo & Ciotti Fig. 5.—Initial(top)andfinal(bottom)differentialenergydistributions. TheenergyperunitmassE isinunitsofE∗d/M∗ (modelD), and E∗n/M∗ (models M1 and N). The energy zero points inmodels D and M1 are such that the most bound particles of the M1 and N end-products havethesameenergy,andthehighestenergyparticlesofmodelsDandM1att=0havethesameenergy. exponential function over most of the populated energy Sanders 2000). For example, when interpreting the FP range (Binney 1982; van Albada 1982;Ciotti 1991;Lon- tilt in Newtonian gravity one can invoke a systematic drillo et al. 1991; NLC06). In contrast, in model D the and fine-tuned increase of the galaxy dark-to-luminous final N(E) decreases for increasing energy, qualitatively mass ratio with luminosity (e.g., Bender, Burstein & preservingits initialshape. Inthe caseofsimulationM1 Faber 1992; Renzini & Ciotti 1993; Ciotti, Lanzoni & it is apparent a dichotomy between a Newtonian part Renzini1996),whileinMONDthe tiltshouldbe related at lower energies (more bound particles), where N(E) tothecharacteristicaccelerationa . Note,however,that 0 is exponential, and a dMOND part at higher energies, in MOND as well as in Newtonian gravity other impor- wherethefinalN(E)resemblestheinitialone. Weinter- tantphysicalproperties may help to explainthe FP tilt, pretthisresultasanothermanifestationofalesseffective suchasa systematicincreaseofradialorbitalanisotropy phase-spacereorganizationinMONDthaninNewtonian with mass or a systematic structural weak homology collapses. (Bertin et al. 2002). Due to the relevance of the subject, we attempt here to derive some preliminary hints. 4.3. Comparison with the observed scaling relations of In particular, for the first time, we can compare with the scaling relations of elliptical galaxies MOND models elliptical galaxies producedbyaformationmechanism,yetassimpleasthe It is not surprising that galaxy scaling relations rep- dissipationless collapse. resent an even stronger test for MOND than for New- In this Section we consider the end-products of sim- tonian gravity, due to the absence of dark matter and ulations M1, M2, and M4. As already discussed in the existence of the critical acceleration a with a uni- 0 Section 3, each of the three systems corresponds to a versalvalueintheformertheory(e.g.,seeMilgrom1984;

Dissipationless collapses in MOND

Similar Dissipationless collapses in MOND