Evolution of subhalo orbits in a smoothly-growing host halo potential
MMNRAS , 1–14 (2020) Preprint 5 February 2021 Compiled using MNRAS L A TEX style file v3.0
Evolution of subhalo orbits in a smoothly-growing host halo potential
Go Ogiya , ★ , James E. Taylor , , and Michael J. Hudson , Waterloo Centre for Astrophysics, University of Waterloo, Waterloo, ON N2L 3G1, Canada Department of Physics and Astronomy, University of Waterloo, 200 University Avenue West, Waterloo, Ontario N2L 3G1, Canada
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
The orbital parameters of dark matter (DM) subhaloes play an essential role in determiningtheir mass-loss rates and overall spatial distribution within a host halo. Haloes in cosmologicalsimulations grow by a combination of relatively smooth accretion and more violent mergers,and both processes will modify subhalo orbits. To isolate the impact of the smooth growth of thehost halo from other relevant mechanisms, we study subhalo orbital evolution using numericalcalculations in which subhaloes are modelled as massless particles orbiting in a time-varyingspherical potential. We find that the radial action of the subhalo orbit decreases over the firstfew orbits, indicating that the response to the growth of the host halo is not adiabatic during thisphase. The subhalo orbits can shrink by a factor of ∼ 𝑧 < ∼
3, indicating that any subhaloesaccreted earlier are unresolved in the simulations. We also discuss tidal stripping as a formationscenario for NGC1052-DF2, an ultra diffuse galaxy significantly lacking DM, and find thatits expected DM mass could be consistent with observational constraints if its progenitor wasaccreted early enough, 𝑧 > ∼ .
5, although it should still be a relatively rare object.
Key words: galaxies: haloes – galaxies: kinematics and dynamics – cosmology: dark matter– methods: numerical
In the standard paradigm for structure formation in the Universe, the Λ cold dark matter ( Λ CDM) cosmological model, small dark matter(DM) haloes form early on through gravitational collapse, and thenmerge to form larger structures subsequently. At the same time,baryonic gas cools within haloes over some mass scale, ignitinggalaxy formation. As a result, the hierarchy of observed galaxiesis formed (e.g., White & Rees 1978; Frenk & White 2012), i.e.,smaller DM subhaloes and associated satellite galaxies, orbitingwithin larger host systems consisting of a host DM halo and acentral galaxy. In this paper, we consider the orbital evolution ofsubhaloes in their host halo, given the mass assembly history of thelatter.The orbital parameters of subhaloes are essential in determin-ing their spatial distribution and mass-loss rates, and thus they arean important input to observational tests of the nature of DM. Forinstance, in a structure formation scenario based on an alternativeDM model, warm dark matter (WDM), the primordial density fluc-tuations on small scales are smoothed out by the free-streamingmotion of WDM particles, and fewer subhaloes are formed, while ★ E-mail: [email protected] (GO) the minimum halo mass is greater than in the Λ CDM cosmology(Bode et al. 2001; Angulo et al. 2013; Lovell et al. 2014; Bose et al.2016, and references therein). This should affect both the abundanceand the spatial distribution of subhaloes, since the formation of theearliest, densest structures would be suppressed. Tests of DM basedon observations of gravitational lensing (e.g., Dalal & Kochanek2002; Vegetti et al. 2012; Shu et al. 2015; Hezaveh et al. 2016),gaps in stellar streams (e.g., Carlberg 2012; Ngan & Carlberg 2014;Erkal et al. 2016; Ibata et al. 2020), and annihilation or decay signalsof DM particles (e.g., Strigari et al. 2007; Pieri et al. 2008; Hayashiet al. 2016; Hiroshima et al. 2018; Okoli et al. 2018) also dependstrongly on the expected abundance, spatial distribution and massfunction of DM subhaloes.Subhaloes evolve dynamically within the host halo. Since sub-haloes are extended objects, the gravity of the host halo works as atidal force and gradually strips a subhalo’s mass (e.g., King 1962;Spitzer 1987; Taylor & Babul 2004; Binney & Tremaine 2008). Thesubhalo orbit plays an essential role in determining the dynamicalevolution of subhaloes in the tidal force field. When a subhalo ison a radial orbit, its mass is significantly reduced by tidal strippingonce per orbit, at the pericentre where it feels the strongest tidalforce. On the other hand, a subhalo’s mass is reduced more gradu-ally when it is on a circular orbit. As a response to tidal stripping, theinternal structure of the subhalo is altered through re-virialization © a r X i v : . [ a s t r o - ph . GA ] F e b G. Ogiya et al. (e.g., Hayashi et al. 2003; Peñarrubia et al. 2010; Drakos et al.2017; Ogiya et al. 2019; Delos 2019; Green & van den Bosch 2019;Drakos et al. 2020). An accurate prediction of this tidal evolution isrequired to interpret observational tests of the nature of DM.When the host halo is spherical and static, the orbit of a sub-halo is specified by its energy 𝐸 , and angular momentum vector,or equivalently by the energy, orbital plane, and scalar angular mo-mentum 𝐿 . In this paper, we use two dimensionless parametersto characterise the orbit. The first one, characterising the orbitalenergy, is 𝑥 c ≡ 𝑟 c ( 𝐸 )/ 𝑟 , (1)where 𝑟 c ( 𝐸 ) and 𝑟 are the radius of a circular orbit of the or-bital energy, 𝐸 , and the virial radius of the host halo (see Eq. (3)below). The second parameter, the orbital circularity, characterisesthe angular momentum of subhalo’s orbit, 𝜂 ≡ 𝐿 / 𝐿 c ( 𝐸 ) , (2)where 𝐿 and 𝐿 c ( 𝐸 ) are the angular momentum of the subhalo orbitand the angular momentum of a circular orbit with the same energy.Orbits of 𝜂 = 𝑁 -body simulations (e.g., Tormen 1997; Zentner et al.2005; Khochfar & Burkert 2006; Wetzel 2011; Jiang et al. 2015;van den Bosch et al. 2018) and have shown that the distributions ofthese parameters peak at 𝑥 c ∼ . 𝜂 ∼ .
6, respectively.These studies generally measure the orbital properties of sub-haloes at the time of accretion, but their orbits will evolve sub-sequently, due to several different processes, including dynamicalfriction and self-friction. Dynamical friction is the drag force ex-erted by the wake that forms behind the subhalo in the density fieldof the host halo (Chandrasekhar 1943), due to the subhalo’s owngravity (e.g., Ogiya & Burkert 2016). Self-friction is another dragforce caused by tidally stripped material from the subhalo (Fujiiet al. 2006; Fellhauer & Lin 2007; van den Bosch & Ogiya 2018).Similar to dynamical friction, self-friction works more efficientlywhen the ratio of the subhalo mass to the host mass is larger (Ogiyaet al. 2019; Miller et al. 2020). Note that self-friction operates in allsimulations associated with tidal mass-loss, even if the host halo ismodelled with an analytical potential . These mechanisms couplewith each other non-linearly and drive orbital evolution.Due to the difficulty of developing fully analytical treatmentsfor the underlying physics, semi-analytical modelling, in whichfudge parameters in analytical formulae are calibrated to reproducethe simulation results, is a useful way to study the orbital evolution ofsubhaloes (e.g., Lacey & Cole 1993; Taylor & Babul 2001; Taffoniet al. 2003; Peñarrubia & Benson 2005; Zentner et al. 2005; Boylan-Kolchin et al. 2008; Jiang et al. 2008; Gan et al. 2010; Pullen et al.2014). While the fudge parameters are mainly calibrated for the for-mulation of dynamical friction (Chandrasekhar 1943), self-frictionis also taken into account since the simulations are associated withtidal mass-loss of subhaloes. Using idealised simulations, Milleret al. (2020) found that the impact of self-friction is sub-dominant( ∼ /
10 of dynamical friction) in mergers between a host halo anda subhalo.Other mechanisms, including violent relaxation driven by ma-jor mergers (Lynden-Bell 1967), interactions between subhaloes(e.g., Moore et al. 1996) and the growth of the host halo, make the In the simulations using of an analytical potential for the host, dynamicalfriction is absent because of the absence of the density wakes. potential field time-varying and can alter the orbits of subhaloes.While all of these, together with the drag forces of dynamical frictionand self-friction, contribute to the orbital evolution of subhaloes,their individual impacts have been less well studied and are thusuncertain. The fact that these mechanisms are also coupled makesunderstanding the orbital evolution of subhaloes in cosmologicallyrealistic situations particularly challenging.In the hierarchical structure formation scenario, DM haloesgrow through mergers with other haloes and smooth accretion fromadjoining filaments or the surrounding field. Both processes changethe potential of DM haloes with time, and alter the orbits of theirsubhaloes. To avoid the complexities raised by mergers, which in-troduce several new parameters and internal degrees of freedom tothe problem, in this paper we focus on the impact of the smoothgrowth on subhalo orbits. We use numerical calculations in which aspherical host halo is modelled with an analytical potential, and theorbits of subhaloes are traced with massless particles. This treat-ment allows us to isolate the impacts of the smooth host halo growthfrom those of the other mechanisms.The rest of the paper is organised as follows: § 2 describesthe setup and assumptions of our numerical models. In § 3, westudy the orbital evolution of subhaloes driven by the smooth massgrowth of the host halo, and derive an empirical model describingthe evolution of subhalo orbits. We discuss the spatial distribution ofsubhaloes in § 4 and the mass evolution of the possible progenitor ofthe DM-deficient galaxy NGC1052-DF2 in § 5, before summarisingour results in § 6. Throughout this paper, we use the cosmologicalparameters obtained by Planck Collaboration et al. (2016).
We use numerical calculations to study the orbital evolution ofsubhaloes driven by the smooth growth of the spherical host halopotential. In these calculations, subhaloes are treated as masslessparticles, i.e. interactions between subhaloes, dynamical frictionand self-friction are all neglected, such that the (time-varying) ana-lytic potential representing the host halo fully governs their motionafter accretion. Therefore we solve 𝑁 individual one-body prob-lems. The following part of this section describes the justificationfor, and details of, the numerical model.Here, we define basic quantities in our model. The virial massof the DM halo is given as 𝑀 vir ≡ 𝜋 Δ vir ( 𝑧 ) 𝜌 crit ( 𝑧 ) 𝑟 , (3)where Δ vir ( 𝑧 ) and 𝜌 crit ( 𝑧 ) are the virial overdensity and criticaldensity of the universe at given redshift 𝑧 . As seen in Eq. (3),the virial radius, 𝑟 vir , is the radius in which the mean density is Δ vir ( 𝑧 ) times of the critical density at 𝑧 . The virial overdensity of Δ vir ( 𝑧 ) =
200 is employed, and hereafter we denote the virial massand radius of the host halo as 𝑀 and 𝑟 . Throughout this paper,the DM host halo is taken to be spherical, with a Navarro-Frenk-White (NFW) density profile (Navarro et al. 1997), 𝜌 ( 𝑟 ) = 𝜌 s ( 𝑟 / 𝑟 s )( + 𝑟 / 𝑟 s ) . (4)Here 𝑟 , 𝜌 s and 𝑟 s are the distance from the host halo’s centre, thescale density, and the scale length of the host halo, respectively. Thehalo concentration, 𝑐 , is defined as 𝑐 ≡ 𝑟 / 𝑟 s . (5) MNRAS , 1–14 (2020) ubhalo orbits in smoothly-growing host haloes A key assumption in our numerical calculations is treatment ofsubhaloes as massless particles, which corresponds to neglectingdynamical friction, self-friction and subhalo-subhalo interactions.Here we justify this approximation.
The strength of dynamical friction and self-friction relative to grav-ity is roughly proportional to the mass of the subhalo (e.g., Chan-drasekhar 1943; Miller et al. 2020). For the low mass systems thatdominate the subhalo population (Giocoli et al. 2008; Springel et al.2008; Jiang & van den Bosch 2016, and references therein), theseeffects will be negligible, and thus treating subhalos as massless par-ticles is justified. On the other hand, the orbits of larger subhaloeswill shrink more than our numerical calculations predict, due to theneglected impacts of dynamical friction and self-friction.What is the condition on subhalo mass, to be able to neglectthese drag forces? Because self-friction only works subdominantly(Miller et al. 2020), we will only consider dynamical friction inwhat follows. According to Chandrasekhar’s formula for dynamicalfriction, the deceleration is given as a DF = − 𝜋𝐺 ln Λ 𝜌 h 𝑀 s 𝑣 v 𝑣 , (6)where 𝐺 , ln Λ and 𝜌 h are the gravitational constant, the Coulomblogarithm, and the mass density of the host halo. A subhalo witha mass of 𝑀 s moves in the host halo potential with velocity 𝑣 . Wedefine the orbital decay timescale of subhaloes accreted at 𝑧 acc as 𝜏 ( 𝑧 acc ) = 𝑣 /| 𝑎 DF | = 𝑓 (cid:2) 𝜋𝐺 ( ln Λ ) 𝜌 crit ( 𝑧 acc ) (cid:3) − / , (7)where 𝑓 ≡ 𝑀 h / 𝑀 s . Here the host halo mass is given as 𝑀 h ≡ 𝑀 ( 𝑧 acc ) = ( 𝜋 / ) 𝜌 crit ( 𝑧 acc ) 𝑟 ( 𝑧 acc ) (Eq. (3)) and we as-sume that the subhalo is on a circular orbit at the virial radius of thehost halo at 𝑧 acc , 𝑟 ( 𝑧 acc ) .Equating Eq. (7) to the lookback time corresponding to theaccretion redshift, 𝑧 acc , we derive the critical ratio of the host halomass to the subhalo mass, 𝑓 crit ( 𝑧 acc ) . Dynamical friction and self-friction will be negligible for subhaloes of mass 𝑀 s < 𝑀 h / 𝑓 crit ,while more massive systems will sink to the centre of the host haloand merge with it by the present time. In Fig. 1, we show 𝑓 crit as afunction of 𝑧 acc assuming ln Λ = 𝑓 crit is well fitted bylog 𝑓 crit ( 𝑧 acc ) = . 𝑧 . + log ( ln Λ / ) . (8)Because the lookback time is shorter and the critical density de-creases at low redshift, 𝑓 crit decreases as 𝑧 acc goes to zero. UsingEq. (8), we find that dynamical friction and self-friction do not alterthe orbit of subhaloes accreted at 𝑧 acc = ∼ . 𝑧 acc .The impact of dynamical friction estimated in Eq. (6) is in factan upper limit, given tidal mass-loss, which reduces the efficiency ofdynamical friction. Thus the estimated timescale 𝜏 will be a lowerlimit, and 𝑓 crit will be smaller than estimated. While the detailedmass-loss history for a subhalo will depend on its structure and orbit, 𝜏 gets three times longer in typical cases (Mo et al. 2010) and thus 𝑓 crit is reduced by a factor of three. In addition, only particles havingvelocity less than 𝑣 are actually expected to cause dynamical frictionwhile Eq. (6) supposes all particles in the host halo contribute. l o g ( f c r i t ) z acc Figure 1.
The ratio of the host halo mass to the subhalo mass, 𝑓 crit , abovewhich dynamical friction is negligible. The black curve is derived by equat-ing the orbital decay time Eq. (7) to the lookback time for that accretionredshift, 𝑧 acc . A Coulomb logarithm of ln Λ = 𝑓 crit decreases as 𝑧 acc goes to zero, dynamical friction becomes less important at recent times. Interactions between subhaloes are another possible mechanism todrive orbital evolution, whose impact is neglected in our numericalcalculations. To estimate how subhalo-subhalo interactions affectorbital evolution, we will use a toy model based on the argument in§ 1.2.1 of Binney & Tremaine (2008) The details of this model aredescribed in Appendix A.We consider two channels by which subhalo-subhalo interac-tion can alter the orbit of a given ‘subject’ subhalo. The first oneis the deflection of the subject subhalo’s orbit by an angle of morethan 90 degrees (a ‘strong’ deflection) due to a single close en-counter with another subhalo (the ‘perturber’). The second channelis through the cumulative effect of many weak encounters. Fig. 2studies how many strong deflections are expected to occur, and howstrong the effect of cumulative weak encounters should be, as afunction of accretion redshift. We find that while strong deflectionsare generally unlikely (solid curve), the cumulative impact of manyweak encounters may scatter the orbital parameters away from theirinitial values for subhaloes accreted at 𝑧 acc > ∼ < ∼ 𝑧 acc . We expect that these perturbers wouldhave little effect on a subject subhalo’s orbit, however; given theirlarger orbits and longer orbital periods, the subject subhalo couldrespond adiabatically to any changes that they generate in the mainpotential. More detailed studies based on cosmological 𝑁 -bodysimulations are needed to fully explore the effect of encounters onsubhalo orbits, but overall we expect it to be minor in most cases. MNRAS000
The ratio of the host halo mass to the subhalo mass, 𝑓 crit , abovewhich dynamical friction is negligible. The black curve is derived by equat-ing the orbital decay time Eq. (7) to the lookback time for that accretionredshift, 𝑧 acc . A Coulomb logarithm of ln Λ = 𝑓 crit decreases as 𝑧 acc goes to zero, dynamical friction becomes less important at recent times. Interactions between subhaloes are another possible mechanism todrive orbital evolution, whose impact is neglected in our numericalcalculations. To estimate how subhalo-subhalo interactions affectorbital evolution, we will use a toy model based on the argument in§ 1.2.1 of Binney & Tremaine (2008) The details of this model aredescribed in Appendix A.We consider two channels by which subhalo-subhalo interac-tion can alter the orbit of a given ‘subject’ subhalo. The first oneis the deflection of the subject subhalo’s orbit by an angle of morethan 90 degrees (a ‘strong’ deflection) due to a single close en-counter with another subhalo (the ‘perturber’). The second channelis through the cumulative effect of many weak encounters. Fig. 2studies how many strong deflections are expected to occur, and howstrong the effect of cumulative weak encounters should be, as afunction of accretion redshift. We find that while strong deflectionsare generally unlikely (solid curve), the cumulative impact of manyweak encounters may scatter the orbital parameters away from theirinitial values for subhaloes accreted at 𝑧 acc > ∼ < ∼ 𝑧 acc . We expect that these perturbers wouldhave little effect on a subject subhalo’s orbit, however; given theirlarger orbits and longer orbital periods, the subject subhalo couldrespond adiabatically to any changes that they generate in the mainpotential. More detailed studies based on cosmological 𝑁 -bodysimulations are needed to fully explore the effect of encounters onsubhalo orbits, but overall we expect it to be minor in most cases. MNRAS000 , 1–14 (2020)
G. Ogiya et al. N Δv /v i n d i c a t o r z acc Figure 2.
The expected number of close encounters (solid) and the impactof cumulative weak encounters (dotted) as a function of the redshift atwhich the subhalo was accreted into the host halo. A final virial mass of 𝑀 = 𝑀 (cid:12) is assumed, but the results are insensitive to 𝑀 . Whilestrong deflections by subhalo-subhalo interactions are not expected, anyinformation about the orbits of subhaloes accreted at 𝑧 acc > ∼ A number of previous studies have proposed models describing themean mass assembly history (MAH) of DM haloes, based on thebehaviour seen in cosmological 𝑁 -body simulations (e.g., Wech-sler et al. 2002; McBride et al. 2009; Wong & Taylor 2012; vanden Bosch et al. 2014), or on extended Press-Schechter theory (e.g.,Lacey & Cole 1993; van den Bosch 2002; Correa et al. 2015a). Inthis paper, given a virial mass for the host halo at the present time, 𝑀 ( 𝑧 = ) ≡ 𝑀 , we use the model of Correa et al. (2015b) toevaluate its prior evolution. The model requires a final halo concen-tration at 𝑧 = 𝑐 ( 𝑀, 𝑧 ) of Ludlow et al. (2016) to set this. Fig. 3 shows the evolutionof host halo mass and virial radius over time. DM haloes with lowerfinal mass (redder lines) grow earlier than those with higher finalmass (bluer lines), consistent with the basic pattern of hierarchicalstructure formation seen in cosmological 𝑁 -body simulations (e.g.,Fakhouri et al. 2010).Based on the model outlined above, we show the evolution ofthe radial profiles of the enclosed mass, 𝑀 ( < 𝑟 ) , and the gravi-tational potential, Φ ( 𝑟 ) , of the host halo with 𝑀 = 𝑀 (cid:12) inFig. 4. Note that the radial bins are given in (fixed) physical kpc. Inthe upper panel, we see that at the outskirts of the halo, close to thevirial radius (located roughly at the break point in the profiles), theenclosed mass 𝑀 ( < 𝑟 ) increases steadily with time. On the otherhand, close to the centre of the halo, the enclosed mass increasesonly until 𝑧 =
6, and then decreases with time. This behaviourseems inconsistent with some previous analyses of individual halosin cosmological 𝑁 -body simulations, notably (e.g., Diemand et al.2007). The DM halo analysed in Diemand et al. (2007) was isolated,however (with a last major merger at 𝑧 = . M [ M ⊙ ] M =10 M ⊙ M =10 M ⊙ M =10 M ⊙ M =10 M ⊙ M =10 M ⊙ M =10 M ⊙ r [ k p c ] Figure 3.
The growth of the virial mass ( upper ) and radius ( lower )of the host halo. Solid curves show models with a final virial mass of 𝑀 = , , , , and 10 𝑀 (cid:12) . Lower (Higher) masshaloes grow earlier (later). tion of the MAH model and 𝑐 ( 𝑀, 𝑧 ) relation that our calculationassumes. The presence or absence of major mergers may explain thedifference between the two pictures. We have tested several variantsof the MAH model and/or 𝑐 ( 𝑀, 𝑧 ) relation (Appendix B), but theevolution of the central mass structure is always qualitatively consis-tent with that shown in Fig. 4. The mass growth deepens the overallhalo potential, with more growth in the centre than in the outskirts(lower panel). The potential of NFW haloes is relatively flat in thecentre of the halo. Because of this, the potential at 𝑟 < ∼ 𝑟 s deepensalmost as much as that at 𝑟 =
0. Similar results hold for the modelswith other values of 𝑀 . The results shown in the lower panel arerobust even if the other MAH models and/or 𝑐 ( 𝑀, 𝑧 ) relations areemployed (Fig. B1). In our numerical calculations, the host halo centre is fixed at theorigin and the host-centric coordinates are used in what follows.Massless particles (subhaloes) fall into the host halo potential at theaccretion redshift, 𝑧 acc . The host halo potential at 𝑧 acc and the twoorbital parameters, 𝑥 c and 𝜂 (Eq. (1) and Eq. (2)) specify the orbital MNRAS , 1–14 (2020) ubhalo orbits in smoothly-growing host haloes z=0.0z=1.5z=3.0z=4.5z=6.0z=7.5 M ( < r ) / M z = ( < r ) | Φ ( r ) |/ v , Figure 4.
Evolving mass and potential profiles of the model of 𝑀 = 𝑀 (cid:12) . ( Upper ) The enclosed mass profile at a given redshift 𝑧 , 𝑀 ( < 𝑟 ) ,relative to the profile at 𝑧 = 𝑀 z = ( < 𝑟 ) . ( Lower ) The potential profileat a given 𝑧 , Φ ( < 𝑟 ) , relative to the virial velocity squared at 𝑧 = 𝑣 , ≡ 𝐺𝑀 / 𝑟 ( 𝑧 = ) . The radial bins are given in (fixed) physicalkpc. The profiles at 𝑧 are computed in the range of 𝑟 = [ . , 𝑟 ( 𝑧 ) ] where 𝑟 ( 𝑧 ) is the virial radius of the halo at 𝑧 . We note that the centraldensity of the DM halo was higher at earlier times, peaking around 𝑧 ∼ energy and angular momentum of massless particles at accretion.In the numerical calculations presented in § 3, we linearly samplethem over the range of 𝑥 c = [ . ] and 𝜂 = [ .
05 : 0 . ] in 10steps of Δ 𝑥 c = .
15 and Δ 𝜂 = .
09 (i.e. 121 combinations at given 𝑧 acc ). The two-dimensional parameter space covers the most of theinitial subhalo orbits found in cosmological 𝑁 -body simulations(e.g., Tormen 1997; Zentner et al. 2005; Khochfar & Burkert 2006;Wetzel 2011; Jiang et al. 2015; van den Bosch et al. 2018). Whennot specified, at 𝑧 = 𝑧 acc , massless particles are located at theapocentre, 𝑟 a , with zero radial velocity; we determine the amplitudeof the tangential velocity using the pair of orbital parameters andthe host halo potential at 𝑧 acc . The polar and azimuthal anglesat accretion are randomly drawn. The amplitudes of the azimuthand polar components of the particle’s velocity ( 𝑣 𝜃 and 𝑣 𝜙 ) are randomly determined and satisfy 𝜂𝐿 c ( 𝐸 )/ 𝑟 i = ( 𝑣 𝜃 + 𝑣 𝜙 ) / , where 𝑟 i is the initial distance of the subhalo to the host halo centre. We perform the numerical calculations from 𝑧 = . 𝑧 = 𝑀 . 𝑀 is sam-pled over the range from log [ 𝑀 / 𝑀 (cid:12)] =
10 to 15 in 15 steps of Δ log [ 𝑀 / 𝑀 (cid:12)] = / Δ 𝑧 = .
1. Note that the number ofmassless particles in the numerical calculations does not reflect theactual number of subhaloes in a typical halo; instead, the numericalcalculations sample subhalo orbits in the four dimensional spaceof i) accretion redshift, 𝑧 acc ; ii) final mass of the host halo, 𝑀 ; iii)orbital energy parameter at accretion, 𝑥 c , i ; and iv) orbital circularityat accretion, 𝜂 i . In the end, we have 76 redshift bins and the totalnumber of data points is ∼ . × . Neglecting interactions be-tween subhaloes, computing the gravity of the spherical host halopotential is straightforward. We update the host halo mass and struc-ture based on the models of Correa et al. (2015b) and Ludlow et al.(2016) every Δ 𝑡 ≈ . × yr. The particle orbits are integratednumerically using a second-order Leapfrog scheme with the samefixed timestep, Δ 𝑡 . Numerical calculations varying Δ 𝑡 , which con-trol the smoothness of the host halo growth and the accuracy oforbit integration, confirm the numerical convergence of the results.We have also checked that the results are insensitive to the choiceof Δ 𝑧 , Δ log [ 𝑀 / 𝑀 (cid:12)] , Δ 𝑥 c and Δ 𝜂 . If the mass growth of the spherical host halo is adiabatic, i.e. thetimescale for this growth is long compared to the typical orbital pe-riod for subhaloes, the evolution of the subhalo’s orbit is specifiedanalytically, based on the conserved actions, which are adiabaticinvariants. Models based on the adiabatic invariants provide a rea-sonable description of the properties of observed galaxies (e.g.,Blumenthal et al. 1986; Dutton et al. 2007), as well as those ofgalaxies in cosmological hydrodynamical simulations (e.g., Gnedinet al. 2004). One of the adiabatic invariants is the angular momentumof the subhalo orbit 𝐿 ; this is perfectly conserved in our numericalcalculations because of the assumption of a spherical host halo, andthe use of massless particles. The other one is the radial action, 𝐽 r = ∮ orbit 𝑣 r ( 𝑟 (cid:48) ) 𝑑𝑟 (cid:48) (9)where 𝑣 r ( 𝑟 ) is the radial velocity of the subhalo at 𝑟 .While 𝐿 and 𝐽 r are the actual adiabatic invariants for sphericalsystems, their proxies are usually used to discuss the adiabatic evo-lution of DM and stellar orbits. The standard model by Blumenthalet al. (1986) assumes that a test particle moves on a circular orbitin an evolving spherical system; the proxy, [ 𝐺 𝑀 ( < 𝑟 ) 𝑟 ] / , corre-sponds to the specific angular momentum of the test particle, whichis conserved in the evolution. However, the assumption of circularorbits is clearly unrealistic for most cases. Gnedin et al. (2004) pro-posed instead a modified proxy for 𝐽 r , 𝐾 ≡ [ 𝐺 𝑀 ( < ¯ 𝑟 ) ¯ 𝑟 ] / , where¯ 𝑟 is the orbit-averaged radius,¯ 𝑟 = ∮ orbit 𝑟 (cid:48) 𝑑𝑟 (cid:48) 𝑣 r ( 𝑟 (cid:48) ) (cid:30)∮ orbit 𝑑𝑟 (cid:48) 𝑣 r ( 𝑟 (cid:48) ) . (10) MNRAS000
1. Note that the number ofmassless particles in the numerical calculations does not reflect theactual number of subhaloes in a typical halo; instead, the numericalcalculations sample subhalo orbits in the four dimensional spaceof i) accretion redshift, 𝑧 acc ; ii) final mass of the host halo, 𝑀 ; iii)orbital energy parameter at accretion, 𝑥 c , i ; and iv) orbital circularityat accretion, 𝜂 i . In the end, we have 76 redshift bins and the totalnumber of data points is ∼ . × . Neglecting interactions be-tween subhaloes, computing the gravity of the spherical host halopotential is straightforward. We update the host halo mass and struc-ture based on the models of Correa et al. (2015b) and Ludlow et al.(2016) every Δ 𝑡 ≈ . × yr. The particle orbits are integratednumerically using a second-order Leapfrog scheme with the samefixed timestep, Δ 𝑡 . Numerical calculations varying Δ 𝑡 , which con-trol the smoothness of the host halo growth and the accuracy oforbit integration, confirm the numerical convergence of the results.We have also checked that the results are insensitive to the choiceof Δ 𝑧 , Δ log [ 𝑀 / 𝑀 (cid:12)] , Δ 𝑥 c and Δ 𝜂 . If the mass growth of the spherical host halo is adiabatic, i.e. thetimescale for this growth is long compared to the typical orbital pe-riod for subhaloes, the evolution of the subhalo’s orbit is specifiedanalytically, based on the conserved actions, which are adiabaticinvariants. Models based on the adiabatic invariants provide a rea-sonable description of the properties of observed galaxies (e.g.,Blumenthal et al. 1986; Dutton et al. 2007), as well as those ofgalaxies in cosmological hydrodynamical simulations (e.g., Gnedinet al. 2004). One of the adiabatic invariants is the angular momentumof the subhalo orbit 𝐿 ; this is perfectly conserved in our numericalcalculations because of the assumption of a spherical host halo, andthe use of massless particles. The other one is the radial action, 𝐽 r = ∮ orbit 𝑣 r ( 𝑟 (cid:48) ) 𝑑𝑟 (cid:48) (9)where 𝑣 r ( 𝑟 ) is the radial velocity of the subhalo at 𝑟 .While 𝐿 and 𝐽 r are the actual adiabatic invariants for sphericalsystems, their proxies are usually used to discuss the adiabatic evo-lution of DM and stellar orbits. The standard model by Blumenthalet al. (1986) assumes that a test particle moves on a circular orbitin an evolving spherical system; the proxy, [ 𝐺 𝑀 ( < 𝑟 ) 𝑟 ] / , corre-sponds to the specific angular momentum of the test particle, whichis conserved in the evolution. However, the assumption of circularorbits is clearly unrealistic for most cases. Gnedin et al. (2004) pro-posed instead a modified proxy for 𝐽 r , 𝐾 ≡ [ 𝐺 𝑀 ( < ¯ 𝑟 ) ¯ 𝑟 ] / , where¯ 𝑟 is the orbit-averaged radius,¯ 𝑟 = ∮ orbit 𝑟 (cid:48) 𝑑𝑟 (cid:48) 𝑣 r ( 𝑟 (cid:48) ) (cid:30)∮ orbit 𝑑𝑟 (cid:48) 𝑣 r ( 𝑟 (cid:48) ) . (10) MNRAS000 , 1–14 (2020)
G. Ogiya et al.
While they found that the modified form, 𝐾 (cid:48) = [ 𝐺 𝑀 ( < ¯ 𝑟 ) 𝑟 ] / ,was a better proxy in their simulations, in our calculations we findthat 𝐾 (cid:48) can fluctuate by a factor of > ∼
10 over a single orbital period,because of the variation in the instantaneous value of 𝑟 . Thus, in thediscussion below, 𝐾 is used.We will also consider another proxy of the adiabatic invariants, 𝐿 c = [ 𝐺 𝑀 ( < 𝑟 c ) 𝑟 c ] / . The advantage of using 𝐿 c is that it doesnot require orbit integration, unlike 𝐾 , and we can evaluate it at anyarbitrary phase of the orbit. Since 𝑟 c roughly corresponds to theapocentre, 𝑟 a , 𝐿 c is an approximation of another proxy consideredby Blumenthal et al. (1986), 𝐺 𝑀 ( < 𝑟 a ) 𝑟 a . This corresponds to 𝐽 when the orbit is purely radial.A way to estimate the adiabaticity of the orbital evolution ofsubhaloes is to compare the timescale of the change in the host halopotential, 𝑡 Φ , to the orbital period for subhaloes, 𝑡 orb . The former isgiven as 𝑡 Φ ( 𝑧, 𝑟 ) ≡ Φ ( 𝑧, 𝑟 )/ (cid:164) Φ ( 𝑧, 𝑟 ) , (11)where Φ ( 𝑧, 𝑟 ) is the potential profile of the host halo at redshift, 𝑧 ,and (cid:164) Φ ( 𝑧, 𝑟 ) is its time derivative. The latter is proportional to thedynamical time of the halo, 𝑡 dyn ( 𝑧, 𝑟 ) = √︄ 𝜋 𝐺 ¯ 𝜌 ( 𝑧, 𝑟 ) = 𝜋 √︄ 𝑟 𝐺 𝑀 ( 𝑧, 𝑟 ) , (12)where ¯ 𝜌 ( 𝑧, 𝑟 ) = 𝑀 ( 𝑧, 𝑟 )/( 𝜋𝑟 ) is the mean density within 𝑟 . Weparametrise the ratio of 𝑡 orb to 𝑡 dyn as 𝛼 . According to Ogiya et al.(2019), 𝑡 orb ∼ . 𝑥 . Gyr at 𝑧 =
0. Putting a typical value ataccretion, 𝑥 c = .
2, while comparing 𝑡 orb with 𝑡 dyn [ , 𝑟 ( )] , wederive 𝛼 = . 𝑡 Φ to 𝑡 orb . Each line colour representsthe results for a model with given value of 𝑀 . Measuring thetimescales at the halo outskirt ( 𝑟 ( 𝑧 ) ; solid), we find that 𝑡 Φ < 𝑡 orb at 𝑧 > ∼
3. While the ratio increases with decreasing of 𝑧 , it is still oforder unity in most cases. This indicates that the change in the hosthalo potential is not adiabatic for subhaloes orbiting around 𝑟 ,i.e. ones recently accreted into the host halo. When the timescales aremeasured at the scale radius ( 𝑟 s ( 𝑧 )/
2; dashed), we find 𝑡 Φ (cid:29) 𝑡 orb at all redshifts 𝑧 <
7. Thus the growth of the host halo shouldproduce adiabatic orbital changes for subhaloes in the halo centre.Such subhalos would generally have been accreted earlier, however,and would have experienced some non-adiabatic orbital evolutionsoon after accretion.We test the predicted adiabaticity of subhalo orbits with nu-merical calculations. To compute 𝐽 r ( 𝐾 ) in the 𝑁 p -th orbit, thenumerical integration of Eq. (9) (Eq. (10)) starts when the masslessparticle reaches the 𝑁 p -th apocentre, and continues until the mass-less particle reaches the ( 𝑁 p + ) -th apocentre. When evaluating theproxies 𝐾 and 𝐿 c , the mass profile of the host halo potential at thetime of the ( 𝑁 p + ) -th apocentre approach is used.Fig. 6 presents the evolution of the adiabatic invariant 𝐽 r (black), and its proxies 𝐾 (orange) and 𝐿 c (blue), as a functionof the number of orbital periods, 𝑁 p . We find that 𝐽 r decreases byup to ∼
10 percent over the first few orbital periods. This confirmsthat during this initial phase of evolution, the mass growth of thehost halo potential is not adiabatic for subhalo orbits, so numericalcalculations are needed to accurately model the orbital evolutionof subhaloes in the growing host halo potential. In the later phase, 𝐽 r remains constant, i.e. the host halo growth appears adiabatic forsubhalo orbits. These results are consistent with the expectationfrom the comparison of timescales in Fig. 5. We also find in Fig. 6that the change in the proxies in the non-adiabatic phase is greater dashed: r s (z)/2solid: r (z) M =10 M ⊙ M =10 M ⊙ M =10 M ⊙ M =10 M ⊙ M =10 M ⊙ M =10 M ⊙ t Φ / t o r b Figure 5.
Comparison between the timescale of the change in the potential, 𝑡 Φ , and the orbital period, 𝑡 orb . Lines show the models with a final virial massof 𝑀 = , , , , and 10 𝑀 (cid:12) . The black horizontalline indicates equality between the two timescales, for guidance; above thisline, the response to changes in the host potential should be adiabatic, whilebelow the line it will not be. For subhaloes orbiting close to the virial radiusof the halo at redshift 𝑧 , ( 𝑟 ( 𝑧 ) ; solid), 𝑡 Φ is shorter than or comparableto 𝑡 orb . Thus, the potential change is not adiabatic for recently accretedsubhaloes with apocentres ∼ 𝑟 ( 𝑧 ) . On the other hand, for subhaloesorbiting in the centre of the host halo ( 𝑟 s ( 𝑧 )/
2; dashed lines) the responseto changes in the potential should be adiabatic at all redshifts. M =10 M ⊙ , z acc =7.5 Γ ( N p ) / Γ ( ) p =10 M ⊙ , z acc =7.5 Γ ( N p ) / Γ ( ) p =10 M ⊙ , z acc =7.5 Γ ( N p ) / Γ ( ) p =10 M ⊙ , z acc =3.0Γ = J r Γ = KΓ = L c Γ ( N p ) / Γ ( ) p Figure 6.
Evolution of 𝐽 r (black), 𝐾 (orange), and 𝐿 c (blue), as a function ofthe number of orbital periods, 𝑁 p . Solid lines show the mean of the changein each quantity. Upper and lower dashed (dotted) lines show the 75th and25th (90th and 10th) percentiles of the distribution. The host halo mass at 𝑧 = 𝑀 , and the accretion redshift, 𝑧 acc , are indicated in each panel. 𝐽 r and its proxies change significantly over the first few orbits, but are wellconserved subsequently. MNRAS , 1–14 (2020) ubhalo orbits in smoothly-growing host haloes M =10 M ⊙ , z acc =7.5 r c ( N p ) / r c ( ) p M =10 M ⊙ , z acc =3.0 r c ( N p ) / r c ( ) p Figure 7.
Evolution of 𝑟 c as a function of the number of orbital periods, 𝑁 p .The upper (lower) panel shows results for 𝑀 = 𝑀 (cid:12) and 𝑧 acc = . 𝑟 c . Upper and lower dashed(dotted) lines show the 75th and 25th (90th and 10th) percentiles of thedistribution. In the first few orbits, 𝑟 r decreases significantly, but the decreaseis stalled thereafter. than that in 𝐽 r , while they remain almost constant during the lateradiabatic phase. Thus, the orbital evolution of subhaloes will bemiscalculated if one relies on either of 𝐾 or 𝐿 c evaluated at accre-tion. Also note that the net changes in 𝐽 r , 𝐾 and 𝐿 c are larger forlarger halo masses 𝑀 or earlier accretion times 𝑧 acc . In such cases,the ratio 𝑡 Φ / 𝑡 orb is smaller (Fig. 5), i.e. the change in the host halopotential is less adiabatic.In Fig. 7, we show the evolution of 𝑟 c as a function of orbitalperiod 𝑁 p , for the model with 𝑀 = 𝑀 (cid:12) . Since 𝑟 c correspondsroughly to the apocentric radius, Fig. 7 indicates that the smoothgrowth of the host halo shrinks subhalo orbits over the first feworbital periods. The contraction in orbital radius is similar to theevolution of 𝐽 r , with significant change over the first few orbitalperiods, and little change thereafter. In the upper panel ( 𝑧 acc = . 𝑁 p > ∼ 𝑧 acc = . ∼ 𝑟 < ∼ 𝑧 acc = . 𝑟 c decreases monotonically with time. At 𝑧 acc = .
0, the virial radius of the host halo is ∼
50 kpc (Fig. 4) andrecently accreted subhaloes have an apocentre of a similar value. In this radial range, the enclosed mass increases monotonically withtime, and thus the orbit only ever contracts. We confirm that theevolution of the pericentre and apocentre resembles the 𝑟 c -evolutionclosely, i.e. they are reduced by almost the same factor over the firstfew periods, while the orbital contraction then stops during the lateradiabatic phase.As indicated in Fig. 5, the growth of the host halo potential isless adiabatic when 𝑧 acc and/or 𝑀 is larger. This leads to the moresignificant reduction of 𝑟 c (up to a factor of ∼ 𝑟 c is reduced by a larger factorin orbits with higher 𝑥 c , i or 𝜂 i . This is because subhaloes on theseorbits spend more time in the outskirts of the host halo, where theadiabatic condition is strongly broken. In predicting the orbital evolution of subhaloes in the smoothly-growing host halo potential, 𝐿 c has the advantage relative to 𝐽 r and 𝐾 that it can be evaluated at any phase of the orbit. Since themost significant part of the orbital evolution occurs rapidly (Fig. 7),this advantage will be important for making accurate predictions.However, as shown in Fig. 6, the change in 𝐿 c during the evolutionis greater than those in 𝐽 r and 𝐾 . If one applies a model basedon 𝐿 c , without any further correction, the orbital evolution of sub-haloes will be mispredicted. Thus, to increase the accuracy of ourpredictions, we introduce a correction factor.As shown in the numerical calculations in § 3.1, most of thechange in 𝐿 c occurs during the first few orbital periods. In practice,we find that the amplitude of the change in 𝐿 c depends on all fourpossible variables, i) the accretion redshift, 𝑧 acc ; ii) the final massof the host halo, 𝑀 ; iii) the orbital energy parameter at accretion, 𝑥 c , i ; and iv) the orbital circularity at accretion, 𝜂 i . Motivated by thisobservation, the correction factor is defined as 𝐿 c ( 𝑁 𝜏 , 𝑧 acc , 𝑀 , 𝑥 c , i , 𝜂 i ) 𝐿 c ( , 𝑧 acc , 𝑀 , 𝑥 c , i , 𝜂 i ) = + 𝐴 exp ( 𝐵𝑧 acc + 𝐶𝑀 + 𝐷𝑥 c , i + 𝐸𝜂 i ) tanh ( 𝐹𝑁 𝜏 ) . (13)The number of orbital periods is estimated using the smooth func-tional form (Jiang & van den Bosch 2016), 𝑁 𝜏 = ∫ 𝑧𝑧 acc 𝑑𝑧 (cid:48) ( 𝑑𝑡 / 𝑑𝑧 (cid:48) ) 𝑡 dyn [ 𝑧 (cid:48) , 𝑟 ( 𝑧 (cid:48) )] . (14)We fit the results from numerical calculations for various 𝑀 byusing the curve_fit procedure in the scipy.optimize module. The fitting parameters are obtained as the averages in 𝑘 -fold cross-validation, adopting 𝑘 = 𝐴 = . × − , 𝐵 = . 𝐶 = . 𝐷 = . 𝐸 = − .
832 and 𝐹 = . 𝐿 c at a given 𝑧 . Giventhe mass profile of the host halo at 𝑧 , the predicted 𝐿 c is thenconverted into the two orbital parameters of interest, 𝑥 c and 𝜂 . InFig. 8, we show the distribution of the error in the estimate of theevolved orbital parameters. We find that they are reproduced at thefive percent level. Note that data points of 𝑧 = 𝑧 acc are excluded fromthe analysis. The prediction based on Eq. (13) is applicable at anyarbitrary time, and does not require the orbital integration, unlike 𝐽 r or 𝐾 . In addition, it is more accurate than an estimate of evolvedorbital parameters based on 𝐽 r ( 𝐾 ), which can change by ∼ MNRAS000
832 and 𝐹 = . 𝐿 c at a given 𝑧 . Giventhe mass profile of the host halo at 𝑧 , the predicted 𝐿 c is thenconverted into the two orbital parameters of interest, 𝑥 c and 𝜂 . InFig. 8, we show the distribution of the error in the estimate of theevolved orbital parameters. We find that they are reproduced at thefive percent level. Note that data points of 𝑧 = 𝑧 acc are excluded fromthe analysis. The prediction based on Eq. (13) is applicable at anyarbitrary time, and does not require the orbital integration, unlike 𝐽 r or 𝐾 . In addition, it is more accurate than an estimate of evolvedorbital parameters based on 𝐽 r ( 𝐾 ), which can change by ∼ MNRAS000 , 1–14 (2020)
G. Ogiya et al. ζ=x c C15 + L16 (fiducial)C15 + D19B14 + L16B14 + D19 ( d N / d ε ) d ε ζ=η ( d N / d ε ) d ε pred -ζ sim )/ζ sim −0.2 −0.1 0 0.1 0.2 Figure 8.
Error distributions when estimating the evolved values of theorbital parameters, 𝑥 c (upper) and 𝜂 (lower). Numerical calculations as-suming the MAH model by Correa et al. (2015b, C15) and the 𝑐 ( 𝑀 , 𝑧 ) relation by Ludlow et al. (2016, L16) (black line) are used to derive the fit-ting parameters in Eq. (13). Additional calculations using the MAH modelby van den Bosch et al. (2014, B14) and/or the 𝑐 ( 𝑀 , 𝑧 ) relation by Diemer& Joyce (2019, D19) are also shown for reference. The general evolution ofboth orbital parameters is well reproduced, regardless of the detailed modelsused in the numerical calculation. parameters in Eq. (13) are derived from the numerical calculationsemploying the MAH model by Correa et al. (2015b, C15) andthe 𝑐 ( 𝑀, 𝑧 ) relation by Ludlow et al. (2016, L16). The details ofthe structural evolution of the host halo and the orbital evolutionof subhaloes could in principle depend on this choice of models.To study the dependence, we have performed additional numericalcalculations employing the MAH model by van den Bosch et al.(2014, B14) and/or the 𝑐 ( 𝑀, 𝑧 ) relation by Diemer & Joyce (2019,D19). While the exact values of the fitting parameters in Eq. (13) dodepend on the models chosen, we find that the parameter set listedabove provides excellent accuracy in describing the evolution of 𝑥 c and 𝜂 , even in numerical calculations with the other models. z acc =[0, 7.5]z acc =[0, 5]z acc =[0, 4]z acc =[0, 3]z acc =[0, 2]z acc =[0, 1]z acc =[0, 0.5] n ( r ) [ a r b . un i t ] ∝ r n ( r ) / ρ ( r ) [ a r b . un i t ] −3 n ( r ) / n . ( r ) r/r Figure 9.
Spatial distribution of massless particles representing subhaloesin the numerical calculation of 𝑀 = 𝑀 (cid:12) . The snapshot at 𝑧 = 𝑧 acc . ( Top ) Radial profile ofthe number density of the particles, 𝑛 ( 𝑟 ) . ( Middle ) The ratio of 𝑛 ( 𝑟 ) to thehost halo density profile, 𝜌 ( 𝑟 ) . The dashed line is the scaling found in theAquarius simulations (Springel et al. 2008; Han et al. 2016). ( Bottom ) Theratio of 𝑛 ( 𝑟 ) to that of 𝑧 acc = [ , . ] , 𝑛 . ( 𝑟 ) , showing the accumulationhistory of subhaloes. The radial bins are normalised by the virial radiusof the host halo, 𝑟 , and 𝑛 ( 𝑟 ) is given in arbitrary units. The subhalodistribution found in the cosmological simulation matches the distributionof subhaloes accreted at 𝑧 acc < ∼ The spatial distribution of DM subhaloes has been an importantsubject of study with cosmological 𝑁 -body simulations (Ghignaet al. 2000; Diemand et al. 2004; Nagai & Kravtsov 2005; Ludlowet al. 2009; Gao et al. 2012; Hellwing et al. 2016, and referencestherein). These authors have shown that the radial number densityprofile of subhaloes within their host halo, 𝑛 ( 𝑟 ) , has a core ofconstant number density at the centre of the host halo and decreases MNRAS , 1–14 (2020) ubhalo orbits in smoothly-growing host haloes with distance from the centre of the host halo, 𝑟 , but that its slope isshallower than that of the host halo density profile, 𝜌 ( 𝑟 ) . Motivatedby this observation, Han et al. (2016) advocated a modification ofthe density profile of the host halo that represents the number densityof subhaloes as 𝑛 ( 𝑟 ) ∝ 𝑟 𝛾 𝜌 ( 𝑟 ) , (15)where 𝛾 is a parameter controlling the significance of the modifica-tion. They found that 𝛾 = . 𝛾 analytically, based on the number ofmerging subhaloes, and the rate at which subsequent tidal mass-losscan lead to complete disruption of subhaloes.Recent papers have cast doubt on the accuracy of tidal evolutionfor subhaloes in cosmological simulations, and thus on the actualrate of tidal disruption. For instance, van den Bosch et al. (2018)used analytical formulae from the literature to estimate the tidalmass-loss rate due to tidal shocking and stripping, and found thatneither mechanism can explain the rate of subhalo disruption seen incosmological simulations. van den Bosch & Ogiya (2018) showedthat artificial subhalo disruption occurs when the force softeningis inadequate, or the number of particles used to model a subhalois too small. They concluded that many of the subhalo disruptionsseen in cosmological simulations are artificial (see also Errani &Peñarrubia 2020).Motivated by this situation, we ran one more numerical cal-culation for a Milky Way-like host halo (final host halo mass of 𝑀 = 𝑀 (cid:12) ) like those in the Aquarius simulations, to studythe subhalo spatial distribution free from artificial disruptions. Theresult from the cosmological 𝑁 -body simulations that 𝑛 ( 𝑟 ) is in-dependent of the subhalo mass (Hellwing et al. 2016) justifies theuse of massless particles, since the impacts of dynamical frictionand self-friction would be negligible for subhaloes with low enoughmasses, as shown in § 2.1. Based on Fakhouri et al. (2010), themerger rate at 𝑧 , 𝑑𝑁 / 𝑑𝑧 scales as 𝑑𝑁𝑑𝑧 ∝ 𝑀 ( 𝑧 ) . ( + 𝑧 ) . . (16)While the halo mass definition in Fakhouri et al. (2010) is differentfrom 𝑀 (1-3 times greater than 𝑀 ; Jiang et al. 2014), weuse the scaling relation they found with 𝑀 . The simple scalingworks for subhaloes with low enough mass, the subject of this pa-per. In our additional numerical calculation, the number of subhalomassless particles accreted at 𝑧 is determined with Eq. (16) and 𝑁 tot = particles accrete in total from 𝑧 = . 𝑁 tot does not correspond to the actual number of subhaloes accretedinto a single host halo, we use a large number of massless parti-cles to improve the statistics. At 𝑧 = 𝑧 acc , massless particles areintroduced at 𝑟 i = 𝑟 ( 𝑧 acc ) , with an inward velocity based onresults from cosmological 𝑁 -body simulations. We draw 𝑥 c and 𝜂 by the rejection sampling from the probability distribution function(PDF) of the orbital parameters measured by Jiang et al. (2015).Since the PDF does not strongly depend on host-to-subhalo massratio nor on redshift, we use the fitting result for host haloes of thevirial mass of 10 𝑀 (cid:12) at 𝑧 = 𝑛 ( 𝑟 ) obtained from the particledata at 𝑧 =
0. When we constrain the accretion redshift to lower 𝑧 acc (bluer lines), the central number density gets lower and the sizeof the central core becomes larger. The orbital energy of subhaloesaccreted earlier is lower, and such subhaloes live in the centre of the host halo. In the middle panel, we show the ratio of 𝑛 ( 𝑟 ) to 𝜌 ( 𝑟 ) .We find that when taking subhaloes accreted at 𝑧 acc < ∼
3, the ratio 𝑛 / 𝜌 shows a power-law behaviour, and that the slope is consistentwith the value found by Han et al. (2016) ( 𝛾 = .
3; dashed line).This result implies that the cosmological simulations can resolveonly subhaloes accreted at 𝑧 acc < ∼
3, and that the profile found incosmological simulations may be biased. We show in the bottompanel the ratio of 𝑛 ( 𝑟 ) to that of 𝑧 acc = [ , . ] , 𝑛 . ( 𝑟 ) , indicatingthat subhaloes accreted earlier are located close to the centre, whileones accreted later dominate the outskirts of the host halo. Thispanel implies that inferring the accretion epoch of subhaloes fromtheir position is possible (see also Oman et al. 2013). For instance,looking at the pink line ( 𝑧 acc = [ , ] ), the ratio exceeds 0.5 at 𝑟 / 𝑟 ∼ .
1. Subhaloes at 𝑟 / 𝑟 < . 𝑟 / 𝑟 > .
1) are inferredto have been accreted at 𝑧 acc > 𝑧 acc < 𝑛 ( 𝑟 ) is almost cuspy at the centre (toppanel). Comparing 𝑛 ( 𝑟 ) to 𝜌 ( 𝑟 ) , the subhalo distribution has almostthe same radial dependence as the density profile of the host halo(middle panel). This was in fact implied by Han et al. (2016). Theytraced the position of unresolved (disrupted) subhaloes virtually, bytracking the most bound particle from each disrupted subhalo, andfound that their spatial distribution is similar to the density profile ofthe host halo in a broad radial range (0 . < ∼ 𝑟 / 𝑟 < ∼ Subhaloes accreted at higher 𝑧 acc have lower orbital energies(i.e. smaller 𝑥 c ) at a given 𝑧 . They have smaller pericentres andfeel stronger tidal forces from their host halo than those accretedlater. In addition, their orbital period is shorter, and thus the numberof pericentric passages in a fixed interval of time is larger. Thussubhaloes accreted earlier will experience more significant tidalstripping in their evolution.We here consider the tidal evolution of an ultra diffuse galaxy,NGC1052-DF2 (hereafter, DF2), in the group centred on the largeelliptical galaxy NGC1052. Recently, observations (van Dokkumet al. 2018; Wasserman et al. 2018; Danieli et al. 2019) have inferredthat the DM mass contained in DF2 is several hundred times smallerthan expected from the empirical models of galaxy formation andevolution (Moster et al. 2018; Behroozi et al. 2019). While thereis considerable debate on this interpretation, e.g., discussion of theoverall statistical confidence due to the small number of kinematictracers (Martin et al. 2018; Laporte et al. 2019), the details of thedata processing (Hayashi & Inoue 2018; Trujillo et al. 2019), andthe consistency with the orbital decay timescale due to dynamicalfriction (Nusser 2018; Dutta Chowdhury et al. 2019), consideringthe possible formation processes of such extreme galaxies is aninteresting and important step in advancing our understanding ofgalaxy formation and evolution (Leigh & Fragione 2019; Saleset al. 2019). For instance, Ogiya (2018) showed that DF2 couldbe a result of a violent tidal stripping event, but that a cored DMdensity profile, rather than a cuspy profile like the NFW density For instance, Han et al. (2016) selected subhaloes having at least 1,000particles in their analysis.MNRAS000
1) are inferredto have been accreted at 𝑧 acc > 𝑧 acc < 𝑛 ( 𝑟 ) is almost cuspy at the centre (toppanel). Comparing 𝑛 ( 𝑟 ) to 𝜌 ( 𝑟 ) , the subhalo distribution has almostthe same radial dependence as the density profile of the host halo(middle panel). This was in fact implied by Han et al. (2016). Theytraced the position of unresolved (disrupted) subhaloes virtually, bytracking the most bound particle from each disrupted subhalo, andfound that their spatial distribution is similar to the density profile ofthe host halo in a broad radial range (0 . < ∼ 𝑟 / 𝑟 < ∼ Subhaloes accreted at higher 𝑧 acc have lower orbital energies(i.e. smaller 𝑥 c ) at a given 𝑧 . They have smaller pericentres andfeel stronger tidal forces from their host halo than those accretedlater. In addition, their orbital period is shorter, and thus the numberof pericentric passages in a fixed interval of time is larger. Thussubhaloes accreted earlier will experience more significant tidalstripping in their evolution.We here consider the tidal evolution of an ultra diffuse galaxy,NGC1052-DF2 (hereafter, DF2), in the group centred on the largeelliptical galaxy NGC1052. Recently, observations (van Dokkumet al. 2018; Wasserman et al. 2018; Danieli et al. 2019) have inferredthat the DM mass contained in DF2 is several hundred times smallerthan expected from the empirical models of galaxy formation andevolution (Moster et al. 2018; Behroozi et al. 2019). While thereis considerable debate on this interpretation, e.g., discussion of theoverall statistical confidence due to the small number of kinematictracers (Martin et al. 2018; Laporte et al. 2019), the details of thedata processing (Hayashi & Inoue 2018; Trujillo et al. 2019), andthe consistency with the orbital decay timescale due to dynamicalfriction (Nusser 2018; Dutta Chowdhury et al. 2019), consideringthe possible formation processes of such extreme galaxies is aninteresting and important step in advancing our understanding ofgalaxy formation and evolution (Leigh & Fragione 2019; Saleset al. 2019). For instance, Ogiya (2018) showed that DF2 couldbe a result of a violent tidal stripping event, but that a cored DMdensity profile, rather than a cuspy profile like the NFW density For instance, Han et al. (2016) selected subhaloes having at least 1,000particles in their analysis.MNRAS000 , 1–14 (2020) G. Ogiya et al.
Table 1.
Summary of model parameters describing the mass evolution of the DF2 progenitor. Description of each column: (1) accretion redshift. (2) lookbacktime. (3) host halo mass. (4) concentration of the host halo. (5) subhalo mass. (6) concentration of the subhalo. (7) probability of orbital parameter setsreproducing the DF2 mass criterion, when neglecting the orbit contraction due to the smooth growth of the host halo. (8) probability of orbital parameter setsreproducing the DF2 mass criterion, when taking the orbit contraction into account. Since the host halo mass and structure at 𝑧 = 𝑧 acc 𝑡 lb [Gyr] 𝑀 h [ 𝑀 (cid:12) ] 𝑐 h 𝑀 s [ 𝑀 (cid:12) ] 𝑐 s 𝑃 no − cont 𝑃 with − cont . × . × . × . × . × . × − . × − . × . × . × − . × − x c =0.6, η=0.1z acc =1.0z acc =1.5z acc =2.0Ogiya (2018) b o un d m a ss [ M ⊙ ] lookback time [Gyr]−12 −10 −8 −6 −4 −2 0 Figure 10.
An example of the bound mass evolution of the DF2 progenitormodels. Red, blue and magenta lines show the cases in which the progenitoraccreted into the host system (main progenitor of NGC1052) at 𝑧 acc = . 𝑥 c = . 𝜂 = .
1. Horizontal dashed line is the upper DM mass limitfor DF2, inferred by van Dokkum et al. (2018). profile (Eq. (4)), is necessary to reproduce the observations (seealso Peñarrubia et al. 2010; Yang et al. 2020).A caveat to the proposal of Ogiya (2018) is that the structureof galaxies and the orbital parameters prior to tidal interactionswere based on observations and empirical relations at 𝑧 =
0. SinceDF2 is a member of the galaxy group, the progenitor of DF2 musthave accreted at 𝑧 acc > 𝑧 acc , 𝑀 s , from the current stellar mass of DF2, ∼ × 𝑀 (cid:12) (van Dokkum et al. 2018). The empirical modelalso predicts the mass growth history of the halo surrounding thehost galaxy, NGC1052, given its current stellar mass, ∼ 𝑀 (cid:12) (Forbes et al. 2017). Once the virial masses of the host- and sub-haloes at 𝑧 acc are obtained, we use the 𝑐 ( 𝑀, 𝑧 ) relation by Ludlowet al. (2016) to determine the structure of the two merging systems,assuming NFW density profiles. Table 1 lists the parameter setsobtained for 𝑧 acc = .
0, 1,5 and 2.0.We re-examine the tidal stripping model for a cuspy NFWdensity profile, using the DASH library of idealised 𝑁 -body simu-lations of minor halo mergers. Ogiya et al. (2019) ran more than2,000 high-resolution simulations covering a broad range in the pa-rameters determining the structure of the haloes and the orbit ofthe subhalo. They then trained a non-parametric machine learningmodel to reproduce the mass evolution of subhaloes in the tidalfield of the host halo potential, based on the random forest algo-rithm (Breiman 2001) as implemented in scikit-learn (Pedregosaet al. 2011) . Relative to the simulations, the final predictions ofthe model are accurate at the 0.1 dex level. Although the DASHsimulations consider mergers between pure dark matter structures,they should give a reasonable indication of the overall behaviour.The DASH simulations assume that both the host- and sub-haloeshave NFW density profiles prior to the merger, and that the ratio ofthe host halo to subhalo mass, 𝑓 ≡ 𝑀 h / 𝑀 s , is large. In such merg-ers, we can safely neglect the impacts of dynamical friction andself-friction. While the DASH library is only formally applicable tomergers with 𝑓 > ∼ 𝑓 ∼ −
50 (see Table 1). Note that dynamical fric-tion would not be completely negligible in those cases (Fig. 1). Weneglect the growth of the host halo potential, i.e. the structure ofthe host halo is fixed from 𝑧 = 𝑧 acc to 0, and orbital decay due todynamical friction and self-friction, as assumed in the simulationsperformed by Ogiya et al. (2019).To study the importance of orbital contraction driven by thesmooth growth of the host halo, two cases, with and without orbitalcontraction, are considered for each 𝑧 acc . In the models neglectingorbital contraction, the orbital parameters are unchanged since ac-cretion. In the other models, we virtually take orbital contractioninto account by using Eq. (13). While the orbital parameters are pre-dicted for 𝑧 =
0, they are applied at 𝑧 = 𝑧 acc . The orbital circularity, 𝜂 , is used as it is predicted, and 𝑥 c is scaled by multiplying the factorof 𝑟 ( 𝑧 = )/ 𝑟 ( 𝑧 acc ) to remove the effect of the growth of thevirial radius from the 𝑥 c -evolution. The orbit is assumed to shrinkinstantaneously at 𝑧 acc , while the actually orbit will in fact shrink https://cosmo.oca.eu/dash/ https://scikit-learn.org Because of the potential growth, the same values of 𝑟 c at two differentredshifts do not correspond to the same orbital energy. This is neglected inthe model for simplicity. MNRAS , 1–14 (2020) ubhalo orbits in smoothly-growing host haloes more gradually (Fig. 7). Thus, this treatment provides an upper limiton the estimated strength of the effect of the orbit contraction.van Dokkum et al. (2018) inferred the enclosed DM masswithin 3.1 and 7.6 kpc from the centre of DF2. On the other hand,the machine learning model predicts the mass gravitationally boundto the subhalo at given time measured from the beginning of accre-tion. The comparison is nonetheless justified, for two reasons. First,the process of tidal stripping removes the mass preferentially fromthe outskirts of the subhalo, and hardly changes its central densitystructure (e.g., Hayashi et al. 2003; Peñarrubia et al. 2010; Ogiyaet al. 2019). Second, the radius where the enclosed mass goes belowthe upper mass limit in the subhalo mass profile prior to accretionis ∼ . − . In Fig. 10, we show an example of the expected bound massevolution of the halo surrounding the DF2 progenitor accreted at 𝑧 acc = . 𝑥 c = .
6) radial ( 𝜂 = .
1) orbit is shown. These orbitalparameters are the same as employed in Ogiya (2018). As expected,the number of pericentric passages is larger (8-10) than that seenin Ogiya (2018), 3, and the bound mass goes below the upper masslimit inferred by van Dokkum et al. (2018, horizontal dashed line)in the models with 𝑧 acc = . 𝑥 c , i = [ . , . ] and 𝜂 i = [ .
01 : 0 . ] withthe interval of Δ 𝑥 c , i = Δ 𝜂 i = .
01 (i.e. 14,949 pairs of the orbitalparameters) for 𝑧 acc = .
0, 1.5 and 2.0, and weight each model basedon the PDF of Jiang et al. (2015). We use the fitting parametersfor host haloes of virial mass of 10 𝑀 (cid:12) at 𝑧 = 𝑥 c , i and 𝜂 i satisfying the masslimit when neglecting and when considering orbital contractionare shown in the seventh and eighth columns in Table 1. In themodel taking orbital contraction into account, we set 𝑥 c = . 𝑥 c < . 𝑧 acc > ∼ . ∼
10 when orbital contraction is considered, it remains fairly low,however, and thus DF-2 should be considered very uncommon atbest. A follow-up observational paper reported another DM deficitgalaxy in the same galaxy group, NGC1052-DF4 (van Dokkum An additional justification is provided by the model for the density profileof tidally stripped subhaloes by Green & van den Bosch (2019) that requiresthe subhalo concentration and the ratio of the bound mass of the subhaloat given time to the subhalo mass prior to the merger. We use the machinelearning model to compute the latter and derive the enclosed mass within3.1 kpc from the centre of the DF2 progenitor by the numerical integra-tion. While the predicted bound mass is larger than the enclosed mass, thedifference between them gets smaller when considering more significantmass-loss events and is reduced to ∼
10 percent in the events in which thebound mass goes below the upper mass limit of DF2. Using the enclosedmass instead of the bound mass in the analysis, the probabilities of orbitalparameter sets reproducing the DF2 mass criterion are the same as shownin Table 1. et al. 2019, but see also Monelli & Trujillo 2019). The probabili-ties we have obtained suggest that tidal stripping is very unlikely toexplain the existence of two or more DF-deficit galaxies in a singlegalaxy group. They would increase, however, when taking into ac-count other mechanisms that could shrink the orbit of DF2, such asdynamical friction and self-friction. The probability would increasefurther if the DM halo of DF2 has a cored density profile. A moredetailed study is needed to reach a firm conclusion regarding theformation of these mysterious objects.
Subhaloes are a promising probe of the nature of DM. Astrophysi-cal observations can constrain their masses and spatial distribution,but accurate theoretical predictions of their properties are needed tosupport these efforts. The mass evolution and spatial distribution ofsubhaloes within the host halo depend strongly on their orbital pa-rameters, which in turn are influenced by several different processes.Some of these (e.g. dynamical friction and self-friction) are negligi-ble for all but the most massive subhaloes, while others (e.g. violentrelaxation driven by major mergers) are hard to model analytically.In this paper we have considered two other processes, interactionsbetween subhaloes, and the smooth component of mass growth ofthe host halo, that occurs even in realistic systems formed by hi-erarchical merging. Based on a simple model of subhalo-subhalointeractions, we find that these can erase the memory of initial sub-halo orbital parameters by the present day, if systems are accreted at 𝑧 acc > ∼
5, while large changes due to strong encounters are rare andthus unimportant. For systems accreted at later times, the smoothchange in the background potential is the main effect driving orbitalevolution.To isolate this effect from the other mechanisms, we use nu-merical calculations with massless particles that represent subhaloesorbiting in the smoothly evolving potential of a host halo. Since thesubhalo particles do not have mass, they do not feel drag forces, andthe interactions between them are neglected. We find that the radialaction of subhalo orbits, 𝐽 r , which is a conserved quantity in spher-ical systems evolving adiabatically, decreases by ∼
10 percent soonafter accretion into the host halo, and stays almost constant there-after. The non-conservation indicates that the smooth mass growthof the host halo is not adiabatic for subhalo orbits, and as a resultthey shrink by a factor of ∼ 𝐽 r , 𝐾 and 𝐿 c , is larger than that in 𝐽 r , but theoverall form of the evolution is similar to that of 𝐽 r . We introducean analytic model for the evolution of the orbital parameters basedon 𝐿 c , since it can be evaluated at an arbitrary phase of the orbit.Our model based on the corrected 𝐿 c should accurately de-scribe the orbital evolution of subhaloes whose mass is small enoughcompared to that of the host, that dynamical friction and self-friction are negligible. For more massive subhaloes, these mech-anisms should be properly taken into account. While we have as-sumed the spherical host halo potential, DM haloes in cosmologicalsimulations are in fact triaxial (Jing & Suto 2002; Kuhlen et al.2007; Vera-Ciro et al. 2014, and references therein). The modellingfor spherical systems would be applicable for triaxial systems withsome modification of the proxy for 𝐽 r (e.g., Lithwick & Dalal 2011).The proper modification would depend on the triaxiality of the hosthaloes and further investigations are needed for accurate predictionsof the subhalo orbital evolution. Relaxing the simplifying assump-tions made in this paper would lead to a more realistic model oforbital evolution, although it would make analytical investigations MNRAS000
10 percent soonafter accretion into the host halo, and stays almost constant there-after. The non-conservation indicates that the smooth mass growthof the host halo is not adiabatic for subhalo orbits, and as a resultthey shrink by a factor of ∼ 𝐽 r , 𝐾 and 𝐿 c , is larger than that in 𝐽 r , but theoverall form of the evolution is similar to that of 𝐽 r . We introducean analytic model for the evolution of the orbital parameters basedon 𝐿 c , since it can be evaluated at an arbitrary phase of the orbit.Our model based on the corrected 𝐿 c should accurately de-scribe the orbital evolution of subhaloes whose mass is small enoughcompared to that of the host, that dynamical friction and self-friction are negligible. For more massive subhaloes, these mech-anisms should be properly taken into account. While we have as-sumed the spherical host halo potential, DM haloes in cosmologicalsimulations are in fact triaxial (Jing & Suto 2002; Kuhlen et al.2007; Vera-Ciro et al. 2014, and references therein). The modellingfor spherical systems would be applicable for triaxial systems withsome modification of the proxy for 𝐽 r (e.g., Lithwick & Dalal 2011).The proper modification would depend on the triaxiality of the hosthaloes and further investigations are needed for accurate predictionsof the subhalo orbital evolution. Relaxing the simplifying assump-tions made in this paper would lead to a more realistic model oforbital evolution, although it would make analytical investigations MNRAS000 , 1–14 (2020) G. Ogiya et al. far more complicated. To understand the statistical evolution ofsubhalo orbits, analysing the data from fully realistic cosmologicalsimulations would be a promising avenue.We study the spatial distribution of subhaloes in a Milky Way-sized host halo and find that in current high-resolution cosmologicalsimulations, the dominant fraction of subhaloes surviving at 𝑧 = 𝑧 acc < ∼
3, while those accreted earlierare unresolved. We also consider the implications of our numericalcalculations for the mass evolution of a DM-deficient galaxy, DF2.If the progenitor of DF2 was accreted into the host galaxy at highenough redshift, 𝑧 acc > ∼ .
5, with the orbital parameters of 𝑥 c ∼ . 𝜂 ∼ .
1, tidal stripping by the host galaxy potential canreproduce the upper mass limit inferred from observations. Therequired orbital parameters are in the tail of the PDF, however,making this scenario somewhat unlikely, even if the orbit contractiondriven by the smooth growth of the host halo is taken into account.Additional factors increasing the overall mass-loss efficiency, suchas a further reduction of the orbital size by dynamical friction andself-friction, or increase susceptibility to mass-loss due to a coredDM density profile, would be needed to explain the existence of twoor more DM-deficient galaxies in the galaxy group.
ACKNOWLEDGEMENTS
We thank the anonymous referee and Neal Dalal for providingthe insightful comments that greatly improved the paper. We aregrateful to the developers of SciPy and scikit-learn for makingtheir code publicly available. A part of numerical calculations wasperformed on the Graham cluster operated by Compute Canada( ). JET and MJH acknowledge financialsupport from NSERC Canada, through Discovery Grants.
DATA AVAILABILITY
The data underlying this article will be shared on reasonable requestto the corresponding author.
REFERENCES
Angulo R. E., Hahn O., Abel T., 2013, MNRAS, 434, 3337Behroozi P., Wechsler R. H., Hearin A. P., Conroy C., 2019, MNRAS, 488,3143Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. PrincetonUniversity PressBlumenthal G. R., Faber S. M., Flores R., Primack J. R., 1986, ApJ, 301, 27Bode P., Ostriker J. P., Turok N., 2001, ApJ, 556, 93Bose S., Hellwing W. A., Frenk C. S., Jenkins A., Lovell M. R., Helly J. C.,Li B., 2016, MNRAS, 455, 318Boylan-Kolchin M., Ma C.-P., Quataert E., 2008, MNRAS, 383, 93Breiman L., 2001, Machine Learning, 45, 5Browne M. W., 2000, Journal of Mathematical Psychology, 44, 108Bryan G. L., Norman M. L., 1998, ApJ, 495, 80Carlberg R. G., 2012, ApJ, 748, 20Chandrasekhar S., 1943, ApJ, 97, 255Correa C. A., Wyithe J. S. B., Schaye J., Duffy A. R., 2015a, MNRAS, 450,1514Correa C. A., Wyithe J. S. B., Schaye J., Duffy A. R., 2015b, MNRAS, 450,1521Dalal N., Kochanek C. S., 2002, ApJ, 572, 25Danieli S., van Dokkum P., Conroy C., Abraham R., Romanowsky A. J.,2019, ApJ, 874, L12Delos M. S., 2019, Phys. Rev. D, 100, 063505 Diemand J., Moore B., Stadel J., 2004, MNRAS, 352, 535Diemand J., Kuhlen M., Madau P., 2007, ApJ, 667, 859Diemer B., Joyce M., 2019, ApJ, 871, 168Drakos N. E., Taylor J. E., Benson A. J., 2017, MNRAS, 468, 2345Drakos N. E., Taylor J. E., Benson A. J., 2020, MNRAS, 494, 378Dutta Chowdhury D., van den Bosch F. C., van Dokkum P., 2019, ApJ, 877,133Dutton A. A., van den Bosch F. C., Dekel A., Courteau S., 2007, ApJ, 654,27Erkal D., Belokurov V., Bovy J., Sanders J. L., 2016, MNRAS, 463, 102Errani R., Peñarrubia J., 2020, MNRAS, 491, 4591Fakhouri O., Ma C.-P., Boylan-Kolchin M., 2010, MNRAS, 406, 2267Fellhauer M., Lin D. N. C., 2007, MNRAS, 375, 604Forbes D. A., Sinpetru L., Savorgnan G., Romanowsky A. J., Usher C.,Brodie J., 2017, MNRAS, 464, 4611Frenk C. S., White S. D. M., 2012, Annalen der Physik, 524, 507Fujii M., Funato Y., Makino J., 2006, PASJ, 58, 743Gan J., Kang X., van den Bosch F. C., Hou J., 2010, MNRAS, 408, 2201Gao L., Navarro J. F., Frenk C. S., Jenkins A., Springel V., White S. D. M.,2012, MNRAS, 425, 2169Ghigna S., Moore B., Governato F., Lake G., Quinn T., Stadel J., 2000, ApJ,544, 616Giocoli C., Tormen G., van den Bosch F. C., 2008, MNRAS, 386, 2135Gnedin O. Y., Kravtsov A. V., Klypin A. A., Nagai D., 2004, ApJ, 616, 16Green S. B., van den Bosch F. C., 2019, MNRAS, 490, 2091Han J., Cole S., Frenk C. S., Jing Y., 2016, MNRAS, 457, 1208Hayashi K., Inoue S., 2018, MNRAS, 481, L59Hayashi E., Navarro J. F., Taylor J. E., Stadel J., Quinn T., 2003, ApJ, 584,541Hayashi K., Ichikawa K., Matsumoto S., Ibe M., Ishigaki M. N., Sugai H.,2016, MNRAS, 461, 2914Hellwing W. A., Frenk C. S., Cautun M., Bose S., Helly J., Jenkins A.,Sawala T., Cytowski M., 2016, MNRAS, 457, 3492Hezaveh Y. D., et al., 2016, ApJ, 823, 37Hiroshima N., Ando S., Ishiyama T., 2018, Phys. Rev. D, 97, 123002Ibata R., Thomas G., Famaey B., Malhan K., Martin N., Monari G., 2020,arXiv e-prints, p. arXiv:2002.01488Jiang F., van den Bosch F. C., 2016, MNRAS, 458, 2848Jiang C. Y., Jing Y. P., Faltenbacher A., Lin W. P., Li C., 2008, ApJ, 675,1095Jiang L., Helly J. C., Cole S., Frenk C. S., 2014, MNRAS, 440, 2115Jiang L., Cole S., Sawala T., Frenk C. S., 2015, MNRAS, 448, 1674Jing Y. P., Suto Y., 2002, ApJ, 574, 538Khochfar S., Burkert A., 2006, A&A, 445, 403King I., 1962, AJ, 67, 471Kuhlen M., Diemand J., Madau P., 2007, ApJ, 671, 1135Lacey C., Cole S., 1993, MNRAS, 262, 627Laporte C. F. P., Agnello A., Navarro J. F., 2019, MNRAS, 484, 245Leigh N., Fragione G., 2019, arXiv e-prints, p. arXiv:1903.06717Lithwick Y., Dalal N., 2011, ApJ, 734, 100Lovell M. R., Frenk C. S., Eke V. R., Jenkins A., Gao L., Theuns T., 2014,MNRAS, 439, 300Ludlow A. D., Navarro J. F., Springel V., Jenkins A., Frenk C. S., Helmi A.,2009, ApJ, 692, 931Ludlow A. D., Bose S., Angulo R. E., Wang L., Hellwing W. A., NavarroJ. F., Cole S., Frenk C. S., 2016, MNRAS, 460, 1214Lynden-Bell D., 1967, MNRAS, 136, 101Martin N. F., Collins M. L. M., Longeard N., Tollerud E., 2018, ApJ, 859,L5McBride J., Fakhouri O., Ma C.-P., 2009, MNRAS, 398, 1858Miller T. B., van den Bosch F. C., Green S. B., Ogiya G., 2020, arXive-prints, p. arXiv:2001.06489Mo H., van den Bosch F. C., White S., 2010, Galaxy Formation and EvolutionMonelli M., Trujillo I., 2019, ApJ, 880, L11Moore B., Katz N., Lake G., 1996, ApJ, 457, 455Moster B. P., Naab T., White S. D. M., 2018, MNRAS, 477, 1822Nagai D., Kravtsov A. V., 2005, ApJ, 618, 557Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493MNRAS , 1–14 (2020) ubhalo orbits in smoothly-growing host haloes Ngan W. H. W., Carlberg R. G., 2014, ApJ, 788, 181Nusser A., 2018, ApJ, 863, L17Ogiya G., 2018, MNRAS, 480, L106Ogiya G., Burkert A., 2016, MNRAS, 457, 2164Ogiya G., van den Bosch F. C., Hahn O., Green S. B., Miller T. B., BurkertA., 2019, MNRAS, 485, 189Okoli C., Taylor J. E., Afshordi N., 2018, J. Cosmology Astropart. Phys.,2018, 019Oman K. A., Hudson M. J., Behroozi P. S., 2013, MNRAS, 431, 2307Peñarrubia J., Benson A. J., 2005, MNRAS, 364, 977Peñarrubia J., Benson A. J., Walker M. G., Gilmore G., McConnachie A. W.,Mayer L., 2010, MNRAS, 406, 1290Pedregosa F., et al., 2011, Journal of Machine Learning Research, 12, 2825Pieri L., Bertone G., Branchini E., 2008, MNRAS, 384, 1627Planck Collaboration et al., 2016, A&A, 594, A13Pullen A. R., Benson A. J., Moustakas L. A., 2014, ApJ, 792, 24Sales L. V., Navarro J. F., Penafiel L., Peng E. W., Lim S., Hernquist L.,2019, arXiv e-prints, p. arXiv:1909.01347Shu Y., et al., 2015, ApJ, 803, 71Spitzer L., 1987, Dynamical evolution of globular clustersSpringel V., et al., 2008, MNRAS, 391, 1685Strigari L. E., Koushiappas S. M., Bullock J. S., Kaplinghat M., 2007, Phys.Rev. D, 75, 083526Taffoni G., Mayer L., Colpi M., Governato F., 2003, MNRAS, 341, 434Taylor J. E., Babul A., 2001, ApJ, 559, 716Taylor J. E., Babul A., 2004, MNRAS, 348, 811Tormen G., 1997, MNRAS, 290, 411Trujillo I., et al., 2019, MNRAS, 486, 1192Vegetti S., Lagattuta D. J., McKean J. P., Auger M. W., Fassnacht C. D.,Koopmans L. V. E., 2012, Nature, 481, 341Vera-Ciro C. A., Sales L. V., Helmi A., Navarro J. F., 2014, MNRAS, 439,2863Wasserman A., Romanowsky A. J., Brodie J., van Dokkum P., Conroy C.,Abraham R., Cohen Y., Danieli S., 2018, ApJ, 863, L15Wechsler R. H., Bullock J. S., Primack J. R., Kravtsov A. V., Dekel A., 2002,ApJ, 568, 52Wetzel A. R., 2011, MNRAS, 412, 49White S. D. M., Rees M. J., 1978, MNRAS, 183, 341Wong A. W. C., Taylor J. E., 2012, ApJ, 757, 102Yang D., Yu H.-B., An H., 2020, arXiv e-prints, p. arXiv:2002.02102Zentner A. R., Berlind A. A., Bullock J. S., Kravtsov A. V., Wechsler R. H.,2005, ApJ, 624, 505van Dokkum P., et al., 2018, Nature, 555, 629van Dokkum P., Danieli S., Abraham R., Conroy C., Romanowsky A. J.,2019, ApJ, 874, L5van den Bosch F. C., 2002, MNRAS, 331, 98van den Bosch F. C., Ogiya G., 2018, MNRAS, 475, 4066van den Bosch F. C., Jiang F., Hearin A., Campbell D., Watson D., Padman-abhan N., 2014, MNRAS, 445, 1713van den Bosch F. C., Ogiya G., Hahn O., Burkert A., 2018, MNRAS, 474,3043
APPENDIX A: TOY MODEL FOR SUBHALO-SUBHALOINTERACTIONS
This appendix provides a simple estimate of the relative importanceof subhalo-subhalo interactions in the orbital evolution of subhaloes(the final results of the analysis are presented in § 2.1.2). The modelis a modified version of the argument presented in § 1.2.1 of Binney& Tremaine (2008).In this model for interactions, subhaloes are assumed to bepoint masses. Suppose that a subhalo moves on a straight path pass-ing through the centre of the host halo and passes by another subhalowith a mass of 𝑀 s (the ‘perturber’). For simplicity, during the inter-action, the perturber is fixed and the relative velocity between them is a constant, 𝑣 . We denote the impact parameter (the perpendicu-lar distance between the path and the perturber) as 𝑏 . The velocityperturbation in the perpendicular direction in an interaction is 𝑑𝑣 = 𝐺 𝑀 s / 𝑏𝑣. (A1)When 𝑏 is smaller than 𝑏 ≡ 𝐺 𝑀 s / 𝑣 , (A2) 𝑑𝑣 can be greater than 𝑣 , in which case and the orbit of the ‘subject’subhalo will be deflected by more than 90 degrees. We refer tosuch interactions as close or strong encounters. Even if 𝑏 > 𝑏 ,cumulative impacts of weak interactions with multiple perturberscan alter the orbit of the subject subhalo. Supposing that perturbersare isotropically distributed in the host halo, Δ 𝑣 ≡ (cid:205) 𝑑𝑣 =
0, while Δ 𝑣 ≡ (cid:205) 𝑑𝑣 = 𝑑𝑁 × 𝑑𝑣 >
0. Here 𝑑𝑁 represents the numberof perturbers in [ 𝑏 : 𝑏 + 𝑑𝑏 ] and 𝑑𝑁 = 𝜋𝑏 Σ ( 𝑏 ) 𝑑𝑏 . The columnnumber density of the subhalo is derived by the integration, Σ ( 𝑏 ) = ∫ √︃ 𝑟 − 𝑏 𝑛 (cid:0)√︁ 𝑧 + 𝑏 (cid:1) 𝑑𝑧, (A3)where 𝑛 ( 𝑟 ) is the number density profile of subhaloes in the hosthalo.Host haloes contain multiple subhalo populations in terms oftheir masses. The efficiency of perturbers in altering the subjectsubhalo’s orbit will depend on their mass, 𝑀 s . When 𝑀 s is larger,the impact of a single interaction is larger. For instance, 𝑑𝑣 is pro-portional to 𝑀 . On the other hand, subhaloes with smaller massesare more abundant than those with larger masses (e.g., Giocoli et al.2008; Springel et al. 2008; Jiang & van den Bosch 2016). Accord-ing to cosmological 𝑁 -body simulations, the subhalo mass functionroughly scales as 𝑑𝑁 / 𝑑 ln ( 𝑀 s / 𝑀 h ) ∝ ( 𝑀 s / 𝑀 h ) − in the limit of 𝑀 s / 𝑀 h (cid:28)
1, where 𝑀 h is the host halo mass. Therefore massiveperturbers are the main contributor in altering the subject subhalo’sorbit in a cumulative fashion. They also more efficiently alter thesubject subhalo’s orbit in the close encounter channel, as indicatedby Eq. (A2).Given this argument, we focus on encounters with massivesubhaloes. As discussed in detail in § 4, Eq. (15) models the spatialdistribution of recently accreted subhaloes within the host halo. Suchsubhaloes are more massive than those accreted earlier because ofthe nature of the hierarchical structure formation (larger structuresare formed, and merge, later) and because they have experiencedless tidal mass-loss. We assume that all subhalo populations followthe distribution of Eq. (15) with 𝛾 = .
3, and compute the columnnumber density of subhaloes of all populations (Eq. (A3)). Then thecolumn number density of subhaloes with 𝑀 s is given by Σ ( 𝑏, 𝑀 s ) = Σ ( 𝑏 ) 𝑁 ( 𝑀 s )/ 𝑁 sub , tot , (A4)where 𝑁 ( 𝑀 s ) and 𝑁 sub , tot are the number of subhaloes with 𝑀 s and the total number of subhaloes in the host halo, respectively.We estimate 𝑁 ( 𝑀 s ) and 𝑁 sub , tot for a host halo with given 𝑀 atgiven redshift assuming the model for the MAH of DM haloes byCorrea et al. (2015b, see § 2.2 for details) and the formulation forthe subhalo mass function by Jiang & van den Bosch (2016). The square of the relative velocity between the subject sub-halo and the perturber, 𝑣 , appears in computing 𝑑𝑣 and 𝑏 . Forsimplicity, the subject subhalo is assumed to be on a radial orbitpassing through the centre of the host halo and drag forces (dy-namical friction and self-friction) and the mass growth of the host They referred to it as the ‘evolved subhalo mass function’.MNRAS000
3, and compute the columnnumber density of subhaloes of all populations (Eq. (A3)). Then thecolumn number density of subhaloes with 𝑀 s is given by Σ ( 𝑏, 𝑀 s ) = Σ ( 𝑏 ) 𝑁 ( 𝑀 s )/ 𝑁 sub , tot , (A4)where 𝑁 ( 𝑀 s ) and 𝑁 sub , tot are the number of subhaloes with 𝑀 s and the total number of subhaloes in the host halo, respectively.We estimate 𝑁 ( 𝑀 s ) and 𝑁 sub , tot for a host halo with given 𝑀 atgiven redshift assuming the model for the MAH of DM haloes byCorrea et al. (2015b, see § 2.2 for details) and the formulation forthe subhalo mass function by Jiang & van den Bosch (2016). The square of the relative velocity between the subject sub-halo and the perturber, 𝑣 , appears in computing 𝑑𝑣 and 𝑏 . Forsimplicity, the subject subhalo is assumed to be on a radial orbitpassing through the centre of the host halo and drag forces (dy-namical friction and self-friction) and the mass growth of the host They referred to it as the ‘evolved subhalo mass function’.MNRAS000 , 1–14 (2020) G. Ogiya et al. halo are neglected. Given that the specific orbital energy of thesubject subhalo is 𝐸 = Φ ( 𝑟 ) where Φ ( 𝑟 ) is the gravitationalpotential of the host halo, the subject subhalo oscillates in the rangeof 𝑋 = [− 𝑟 ( 𝑧 acc ) : 𝑟 ( 𝑧 acc )] while 𝑌 = 𝑍 =
0. We define 𝑣 as 𝑣 = ∫ 𝑟 − 𝑟 [ 𝐸 − Φ ( 𝑋 )]N ( 𝑋 ) 𝑑𝑋 (cid:30)∫ 𝑟 − 𝑟 N ( 𝑋 ) 𝑑𝑋. (A5)Here, N ( 𝑋 ) is computed by N ( 𝑋 ) = ∫ √︃ 𝑟 − 𝑋 𝜋𝑅𝑛 (cid:0)√︁ 𝑅 + 𝑋 (cid:1) 𝑑𝑅. (A6)Note that in Eq. (A5) the local velocity squared of the subjectsubhalo is weighted with the number of perturbers at [ 𝑋 : 𝑋 + 𝑑𝑋 ] , N ( 𝑋 ) 𝑑𝑋 , so that Eq. (A5) corresponds to the averaged velocitysquared of the subject subhalo at the the closest approach to theperturbers. This treatment is justified by the fact that the velocityperturbation of single encounters (Eq. (A1)) corresponds to theproduct of the acceleration at the closest approach, 𝐺 𝑀 s / 𝑏 , andthe time duration of the interaction, 2 𝑏 / 𝑣 (Binney & Tremaine2008).Perturbers are counted as close or weak encounters based onEq. (A2), i.e. interactions with 𝑏 ≤ 𝑏 are close encounters whilethose with 𝑏 > 𝑏 are weak encounters. Perturbers located at √ 𝑌 + 𝑍 / 𝑟 = [ − : 1 ] with a mass of 𝑀 s / 𝑀 h = [ − : 1 ] are considered. While Jiang & van den Bosch (2016) adopted thedefinition of virial overdensity by Bryan & Norman (1998), wedefine 𝑀 h ≡ 𝑀 ( 𝑧 acc ) , i.e., Δ vir ( 𝑧 acc ) = 𝑟 ( 𝑧 acc ) are unevolved from 𝑧 = 𝑧 acc to 0. The number of close encountersand Δ 𝑣 are counted in a single crossing of the subject subhalo.We multiply them by the ratio of the lookback time to 𝑧 acc to thecrossing time of the subject subhalo, defined as twice the dynamicaltime of the host halo at 𝑧 acc , 𝑡 cross ( 𝑧 acc ) ≡ 𝑡 dyn [ 𝑧 acc , 𝑟 ( 𝑧 acc )] = √︄ 𝜋 𝐺 𝜌 crit ( 𝑧 acc ) . (A7) APPENDIX B: EVOLUTION OF HOST HALO PROFILES
In this appendix, we see how the enclosed mass and potential pro-files, 𝑀 ( < 𝑟 ) and Φ ( 𝑟 ) , depend on the MAH model and the 𝑐 ( 𝑀, 𝑧 ) relation. The fiducial case employs the MAH model by Correa et al.(2015b) and the 𝑐 ( 𝑀, 𝑧 ) relation by Ludlow et al. (2016) to specifythe internal structure of the host halo at a given redshift 𝑧 , and theradial profiles are shown in Fig. 4. Fig. B1 shows the evolution of theenclosed mass and potential profiles obtained by varying the MAHmodel (Correa et al. 2015b or van den Bosch et al. 2014) and/or the 𝑐 ( 𝑀, 𝑧 ) relation (Ludlow et al. 2016 or Diemer & Joyce 2019), asindicated in the upper panels. We find that the derived profiles areconsistent with the fiducial case; i) the mass at the halo outskirts andthe virial radius of the host halo increases, while the central densitydecreases with time at 𝑧 < ∼
6. ii) the potential in the centre of the halodeepens relative to the value in the outskirts. Thus, the radial pro-files are insensitive to the choice of the MAH model or the 𝑐 ( 𝑀, 𝑧 ) relation. We note that the predicted evolution of the mass profile atthe halo centre contrasts with the results found for individual halosin some cosmological 𝑁 -body simulations (Diemand et al. 2007).More detailed studies are needed to achieve a firm conclusion onthe evolution of the central mass structure of DM haloes. The main results of this paper do not depend on strongly on this behaviour,however, since the most significant orbital evolution of subhaloesoccurs while subhaloes move in outskirts of the host halo, and thecentral mass is much smaller than the virial mass. This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS , 1–14 (2020) ubhalo orbits in smoothly-growing host haloes c(M,z): Diemer & Joyce (2019)MAH: Correa et al. (2015) M ( < r ) / M z = ( < r ) van den Bosch et al. (2014)Ludlow et al. (2016) van den Bosch et al. (2014)Diemer & Joyce (2019)z=0.0z=1.5 | Φ ( r ) |/ v , z=3.0z=4.5 r [kpc]0.1 1 10 100 z=6.0z=7.5 r [kpc]0.1 1 10 100 Figure B1.
Same as Fig. 4 but varying the MAH model and/or the 𝑐 ( 𝑀 , 𝑧 ) relation, as indicated in the upper panels. The fiducial case shown in Fig. 4 employsthe MAH model of Correa et al. (2015b) and the 𝑐 ( 𝑀 , 𝑧 ) relation of Ludlow et al. (2016). ( Upper ) The enclosed mass profile at given 𝑧 , 𝑀 ( < 𝑟 ) , scaledwith that derived in the fiducial case at 𝑧 = 𝑀 z = ( < 𝑟 ) . ( Lower ) The potential profile at given 𝑧 , Φ ( < 𝑟 ) , scaled with the virial velocity squared at 𝑧 = 𝑣 , ≡ 𝐺𝑀 / 𝑟 ( 𝑧 = ) . The radial bins are given in (fixed) physical kpc. The profiles at 𝑧 are computed in the range of 𝑟 = [ . , 𝑟 ( 𝑧 ) ] where 𝑟 ( 𝑧 ) is the virial radius of the halo at 𝑧 . The potential profiles are insensitive to the choice of the MAH model or the 𝑐 ( 𝑀 , 𝑧 ) relation. While the details ofthe evolution of the central mass structure depends on the choices, the same trend, i.e., that the central density decreases with time at 𝑧 < ∼
6, is obtained in allcases. The growth of the halo outskirts is similar in all cases.MNRAS000