Orphan galaxies in semi-analytic models
Facundo M. Delfino, Claudia G. Scoccola, Sofia A. Cora, Cristian A. Vega-Martinez, Ignacio D. Gargiulo
MMNRAS , 1–14 (2020) Preprint 4 February 2021 Compiled using MNRAS L A TEX style file v3.0
Orphan galaxies in semi-analytic models
Facundo M. Delfino, , ★ Claudia G. Scóccola, , Sofía A. Cora , , ,Cristian A. Vega-Martinez , and Ignacio D. Gargiulo , Facultad de Ciencias Astronómicas y Geofísicas, Universidad Nacional de La Plata, Observatorio Astronómico, Paseo del Bosque,B1900FWA La Plata, Argentina Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET), Rivadavia 1917, Buenos Aires, Argentina, Instituto de Astrofísica de La Plata (CCT La Plata, CONICET, UNLP), Observatorio Astronómico, Paseo del Bosque B1900FWA, La Plata, Argentina, Instituto de Investigación Multidisciplinar en Ciencia y Tecnología, Universidad de La Serena, Raúl Bitrán 1305, La Serena, Chile, Departamento de Astronomía, Universidad de La Serena, Av. Juan Cisternas 1200 Norte, La Serena, Chile
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We present an updated model for the evolution of the orbits of "orphan galaxies" to be used in theSAG semi-analytical model of galaxy formation and evolution. In cosmological simulations,orphan galaxies are those satellite galaxies for which, due to limited mass resolution, halofinders lose track of their dark matter subhalos and can no longer be distinguished as self-bound overdensities within the larger host system. Since the evolution of orphans dependsstrongly on the orbit they describe within their host halo, a proper treatment of their evolutionis crucial in predicting the distribution of subhalos and satellite galaxies. The model proposedtakes into account the dynamical friction drag, mass loss by tidal stripping and a proximitymerger criterion, also it is simple enough to be inexpensive from a computational point ofview. To calibrate this model, we apply it onto a dark matter only simulation and compare theresults with a high resolution simulation, considering the halo mass function and the two-pointcorrelation function as constraints. We show that while the halo mass function fails to put tightconstraints on the dynamical friction, the addition of clustering information helps to betterdefine the parameters of the model related to the spatial distribution of subhalos. Using themodel with the best fit parameters allows us to reproduce the halo mass function to a precisionbetter than 5 per cent, and the two point correlation function at a precision better than 10 percent.
Key words: galaxies:formation – galaxies:evolution – galaxies: haloes – methods: numerical
The observed Universe is successfully described by the Λ CDMmodel. According to this concordance model, at early times, theUniverse underwent a period of exponential expansion, called Infla-tion, in which the primordial perturbations in the metric were set-tled. These metric fluctuations produced perturbations in the matterdensity field that are characterised by the matter power spectrum.The large scale structure (LSS) we see today is the result of thegravitational growth of these tiny matter perturbations. It is cur-rently accepted that structure formation proceeds in a hierarchicalway, with small structures being the first ones to collapse and reach astate close to virial equilibrium. Larger structures, like massive darkmatter (DM) haloes, form later by mergers of pre-existing virialisedhaloes, and by accretion of diffuse dark matter (Frenk & White2012). Within the dark matter haloes, the first galaxies are born andevolve (White & Rees 1978).Galaxies are highly non-linear objects that are the result of a ★ E-mail: fdelfi[email protected] complex formation mechanism, which involves several astrophysi-cal processes and spans a wide range of spatial scales (for a reviewof the theory of galaxy formation and evolution see, e.g. Benson2010; Silk & Mamon 2012; Wechsler & Tinker 2018). Due to thenon-linear nature of these processes, it is not feasible to treat themusing analytic methods and hence the use of N-body simulationsis required. The most common approaches to generate simulatedgalaxy populations are
Hydrodynamic Simulations , Halo Occupa-tion Distributions and
Semi-Analytic Models , each of them havingadvantages and disadvantages (see e.g. Vogelsberger et al. 2020).On large scales, hydrodynamic simulations can take into ac-count the effects of both dark matter content and baryons. However,for scales below the resolution limit of the grid, the astrophysi-cal processes are simulated using semi-analytic recipes. Since thenumber of physical processes to be solved is very large, these sim-ulations are usually extremely demanding on computing power. Atpresent, the state-of-the-art of this class of simulations correspondsto the Virgo Consortium’s EAGLE project (Schaye et al. 2015;Crain et al. 2015) and The Next Generation Illustris project (Illus-trisTNG, Springel et al. 2018; Nelson et al. 2018; Pillepich et al. © a r X i v : . [ a s t r o - ph . C O ] F e b Delfino et al. tidal heating ), making it more susceptible to tidal stripping (Zentner &Bullock 2003; Gan et al. 2010; Pullen et al. 2014). Ram pressurestripping (RPS, Gunn & Gott 1972) is another mechanism thatproduces material loss, affecting only the gas content of the satellitegalaxy that moves within a high density environment. Indeed, agalaxy that moves in an eccentric orbit experiences periods of strongRPS as it approaches the pericentre combined with periods of lowRPS as it passes through the apocentre. On the other hand, a slowlydecaying galaxy will experience a continuous increase in RPS asit falls into denser regions (Brüggen & De Lucia 2007). All theseprocesses depend on the position and velocities of satellite galaxies,so a proper treatment of orphan satellite orbits is crucial.In this work, we present an updated treatment for the orbitsof orphan satellite galaxies used in the semi-analytic model SAG(Semi-Analytic Galaxies, Cora 2006; Lagos et al. 2008; Tecce et al.2010; Orsi et al. 2014; Muñoz Arancibia et al. 2015; Gargiulo et al.2015; Cora et al. 2018). In the previous version of SAG, whena subhalo is no longer identified due to a merger with a largerstructure, its corresponding galaxy becomes an orphan. Assuminga circular orbit, the orphan is initially located at a distance of its hostgiven by the virial radius and with a velocity determined by the virialvelocity of the host. The evolution of the radial distance is estimatedusing dynamical friction, with randomly generated positions andvelocities. Finally, the orphan merges with the central galaxy in atime according to the dynamical friction time-scale (Pujol et al.2017). Here we present an improved model for the evolution oforphan galaxies which includes analytic prescriptions for the tidalstripping (TS) and dynamical friction (DF) mechanisms.Using different SAM and HOD models, Pujol et al. (2017)studied the impact that different orphan satellite treatments haveon the clustering signal. In general, it is found that models that donot include orphans present a lower clustering at low scales. Thisresult is consistent with the work of Kitzbichler & White (2008).These authors use a SAM model to show that, without the additionof orphan satellites, the clustering at low scales is much lower thanthat observed in actual galaxy surveys. In addition, Guo & White(2014) show that the inclusion of orphan satellites improves theabundance of subhaloes at low masses and removes the resolutiondependence observed between SHAM catalogues built from simu-lations of different resolutions. Based on these results, we propose toapply the orphan evolution model to a dark matter only simulation,and compare the results to a higher resolution simulation with thesame settings, considering the halo mass function and the two-pointcorrelation function as constraining relations. In the higher reso-lution simulation, there are subhaloes that would not be identifiedby the halo finder algorithm in the lower resolution simulation. Weassume that the difference in those constraining relations will begiven by the orphan galaxies, and use these quantities to tune theparameters of the orbit model.Finally, it is important to remark that, in this study, we only con-sider DM subhalos and in no case do we apply the semi-analyticalmodel. Therefore, when we refer to “orphan galaxies”, “orphansatellites” or simply “orphans”, we always refer to the subhaloscorresponding to orphan galaxies and not to the orphan galaxiesthemselves.This paper is organised as follows: in Section 2, we introducethe orbit model for the orphan galaxies and its many ingredients.Section 3 describes the statistical tools used to analyse galaxy clus-tering and gives details about the simulations used in this work,mdpl2 and smdpl. In Section 4, we describe the method used tocalibrate the free parameters of the model and discuss our results inSection 5. We present our conclusions in Section 6.
MNRAS000
MNRAS000 , 1–14 (2020) rphan galaxies in SAMs As many other SAMs, the SAG model of galaxy formation andevolution takes as input the properties of dark matter haloes andtheir merger trees extracted from cosmological N-body simulations(see e.g. Roukema et al. 1997; Springel et al. 2001b; Croton et al.2006). Based on this information, the model assigns a central galaxyto each new halo that appears in the simulation and then simulates itsevolution. The central galaxies of main halos are referred to as type0 galaxies, while the central galaxies of satellite halos are labelledas type 1 . Finally, orphan galaxies (i. e. galaxies that have lost theirhost halo) are referred to as type 2 galaxies. For type 0 and type1 galaxies, where the halo finder is still able to identify their darkmatter haloes, effects such as dynamical friction or tidal stripping aregiven in a self-consistent way by the base N-body simulation. On theother hand, for type 2 (orphans) galaxies, since they lose their DMhalo, we cannot follow their evolution. To avoid this problem, whena subhalo is no longer detected by the halo finder algorithm, themodel continues integrating its orbit numerically, taking as initialconditions the position, velocity, mass and radius of the halo at theinstant of last identification, as given by the N-body simulation.The simplest model of a subhalo moving within its host halocan be approximated by a point mass without internal structure or-biting in a static potential. This model does not take into account im-portant aspects of the composition and evolution of satellite haloessuch as the internal structure of the subhalo or the interactions withthe material that forms the host halo. For example, the interactionbetween the subhalo and the matter of its host gives rise to a dynam-ical friction force that causes the evolution of the orbits to deviatefrom that of the simplest model. In addition, during their evolu-tion, satellite haloes may experience mass loss due to tidal strippingmechanism or gravitational shocks. Therefore, to describe the sub-halo evolution accurately, it is necessary to take into account all theprocesses mentioned above.Within our semi-analytic model SAG, we consider each orphangalaxy as a particle whose properties may change over time. The or-bits of this type of satellite galaxies are estimated in a pre-processingstep, that is, before applying SAG to the underlying cosmologicalDM simulation. To determine the orbit of an orphan galaxy, weconsider each subhalo as a particle of the same mass moving ina smooth spherical potential generated by its host halo. In orderto take into account the effect of the host halo over the satellite,at each instant, we compute the effect of dynamical friction usingChandrasekhar’s formula and we also take into account the loss ofmaterial by considering a tidal stripping model. If the mass of a sub-halo falls below a certain resolution limit we consider it disrupted;if, at any instant, the satellite-host distance is less than a fractionof the viral radius of the host we consider the satellite galaxy to bemerged with its host. Below, we describe these processes in moredetail.
When a subhalo of total mass 𝑀 sat moves through a large col-lisionless system composed of particles of mass 𝑚 << 𝑀 sat , itperturbs the particle field creating an over-dense region behind it.This “wake” pulls the subhalo in the opposite direction causing adrag force called dynamical friction . Therefore, we can separate theforce experienced by an orphan galaxy orbiting within a massivehalo into two contributions: one due to the potential of the centralhalo, and a higher-order correction due to the background particles(the dynamical friction term). The first part is given by the simple expression 𝑓 𝑖 = − 𝜕 Φ / 𝜕𝑥 𝑖 , that relates the force acting on a particleat a given position with the potential of the main system Φ at that po-sition. The dynamical friction force is given by the Chandrasekharformula (Chandrasekhar 1943; Binney & Tremaine 2008), i. e. F df ( 𝑟 ) = − 𝜋𝐺 𝑀 𝜌 host ( 𝑟 ) ln Λ 𝑉 (cid:20) erf ( 𝑋 ) − 𝑋 √ 𝜋 exp (− 𝑋 ) (cid:21) V 𝑉 , (1)where 𝑟 is the position of the satellite relative to its host halo, V isthe subhalo velocity, 𝑉 = | V | , 𝑋 = 𝑉 /(√ 𝜎 ) with 𝜎 the velocitydispersion of dark matter particles, 𝜌 host represents the densitydistribution of the host halo, ln Λ is the Coulomb logarithm and erfis the Gauss error function.Assuming a NFW profile (Navarro et al. 1997) for the densitydistribution 𝜌 host and an isotropic velocity distribution, the velocitydispersion 𝜎 is given by Łokas & Mamon (2001) 𝜎 ( 𝑠 ) 𝑉 = 𝑔 ( 𝑐 ) 𝑠 ( + 𝑐𝑠 ) ∫ ∞ 𝑠 (cid:20) 𝑔 ( 𝑐𝑠 ) 𝑠 ( + 𝑐𝑠 ) (cid:21) 𝑑𝑠 , (2)where 𝑐 is the concentration parameter, 1 / 𝑔 ( 𝑥 ) = log ( + 𝑥 ) − 𝑥 /( + 𝑥 ) , and 𝑠 is the distance normalised by the virial radius. Furtherdetails on the NFW density profile can be found in Appendix A andin Łokas & Mamon (2001). For simplicity, in this work we use thefollowing approximation for 𝜎 , which is accurate to 1 per cent for 𝑥 in the range 0 . −
100 (Zentner & Bullock 2003) 𝜎 ( 𝑥 ) (cid:39) 𝑉 max . 𝑥 − . + . 𝑥 − . , (3)where 𝑥 = 𝑐𝑠 , and 𝑉 max is the maximum circular velocity related tothe virial velocity via 𝑉 max (cid:39) 𝑉 vir √︁ . 𝑐 𝑔 ( 𝑐 ) .The argument of the Coulomb logarithm can be expressedas Λ = 𝑏 max / 𝑏 min where 𝑏 max and 𝑏 min are the maximum andthe minimum impact parameters for gravitational encounters be-tween the satellite and the background objects (Binney & Tremaine2008). Typically, 𝑏 min corresponds to a close encounter, then 𝑏 min (cid:39) 𝐺 𝑀 sat / 𝑉 where 𝑉 is a velocity typical of the encounter,such as the rms velocity of the background particles. The choiceof the value for 𝑏 max is more ambiguous, and for a finite system istaken to be the characteristic scale of the system.It should be noted that the derivation of Chandrasekhar’s for-mula assumes a massive particle moving in a homogeneous mediumcomposed by an infinite number of low-mass particles with aMaxwellian velocity distribution. However, in the literature, sev-eral works show that this equation is applicable to more generalcontexts, where these hypotheses are not satisfied, if the Coulomblogarithm is chosen appropriately (Weinberg 1986; Cora et al. 1997;Velazquez & White 1999).There has been much debate in the literature about the appro-priate choice of Coulomb logarithm. For example, Springel et al.(2001b) uses an approach given by ln Λ (cid:39) ln ( + 𝑀 cen / 𝑀 sat ) where 𝑀 cen and 𝑀 sat are the masses of the central halo and the satellitesubhalo, respectively. On the other hand, some authors use otherdefinitions that allow them to reproduce results from N-body simu-lations (Hashimoto et al. 2003; Zentner & Bullock 2003; Petts et al.2015, 2016; Ogiya & Burkert 2016). In particular, Hashimoto et al.(2003), hereafter H03, propose a variable Coulomb logarithm. Thischoice avoids the strong circularization effect that is observed when MNRAS , 1–14 (2020)
Delfino et al. comparing these models with the results obtained from N-bodysimulations. Following this, we use the expressionln
Λ = (cid:40) ln ( 𝑟 / 𝑏𝑅 sat ) 𝑟 > 𝑏𝑅 sat 𝑟 ≤ 𝑏𝑅 sat , (4)where 𝑟 is the distance from the satellite subhalo to the centre ofthe host halo, 𝑅 sat is the virial radius of the satellite and 𝑏 is afree parameter. Note that in H03, the Coulomb logarithm is givenby ln Λ = ln ( 𝑟 / . 𝜖 sat ) for 𝑟 > . 𝜖 sat , where 𝜖 sat is the softeninglength corresponding to a Plummer sphere. Here we assume a NFWprofile for the satellite, thus we introduce the virial radius of thesubhalo and leave 𝑏 as a free parameter to be adjusted. As mentioned above, a subhalo orbiting within its host system issubjected to tidal forces. When tidal forces are greater than thegravitational force of the satellite itself, the material becomes un-bound and the satellite loses mass. The dynamical friction force isproportional to 𝑀 (see equation 1), and hence the magnitude ofthe deceleration experienced by a satellite halo is proportional to 𝑀 sat . As a result, mass loss can have a major impact on the orbitalevolution of the satellite halo. For this reason, it is necessary toestimate the amount of mass lost by tidal stripping.We estimate the tidal radius as the distance at which the self-gravity force and the tidal forces cancel out; material outside thisdistance become unbound and could be stripped out from the satel-lite. The tidal radius is given by 𝑟 t = (cid:18) 𝐺 𝑀 sat 𝜔 − 𝑑 Φ / 𝑑𝑟 (cid:19) / , (5)where 𝑀 sat is the mass of the satellite, 𝜔 is its angular velocityand Φ characterise the potential of the host system (King 1962;Taylor & Babul 2001; Zentner & Bullock 2003). This equation isderived under the assumption that the satellite moves in a circularorbit and the potential of the main system is spherically symmetric.But, even under these restricted assumptions, the tidal limit cannotbe represented as a spherical surface, because some particles of thesatellite within 𝑟 t are unbound while other particles outside 𝑟 t mayremain bound to the suhalo (Binney & Tremaine 2008).In general, satellites do not move in circular orbits and thepotential of the host system is not spherically symmetric. As anapproximation, we can still apply equation 5 to eccentric orbits, inwhich case we estimate an instantaneous tidal radius by using thecorresponding instantaneous values, i.e. 𝜔 = | V × r |/ 𝑟 , where V and r are the instantaneous position and velocity of the subhalo.Another aspect that remains unclear is the rate at which the ma-terial located outside the instantaneous tidal radius 𝑟 t is going tobe removed. Following Zentner et al. (2005), we absorb all thesecomplicated details in a free parameter to be adjusted by externalconstraints. Then, the rate of mass loss of the satellite subhalo byTS is given by 𝑑𝑀 sat 𝑑𝑡 = − 𝛼 𝑀 sat ( > 𝑟 t ) 𝑇 orb . (6)Here 𝑇 orb = 𝜋 / 𝜔 , with 𝜔 the instantaneous angular velocity ofthe satellite, and 𝛼 is treated as a free parameter. The value ofthe parameter 𝛼 differs from author to author. For example, Taylor& Babul (2001) and Zentner & Bullock (2003) choose a value 𝛼 =
1; on the other hand, Peñarrubia & Benson (2005) assume aninstantaneous stripping , which effectively implies 𝛼 → ∞ . Finally,some authors (e.g. Zentner et al. 2005, Pullen et al. 2014) varythe value of 𝛼 in order to reproduce the halo mass function ofsimulations. In this paper, we will follow the latter approach. According to hierarchical structure formation models, mergers playa critical role in the formation and evolution of galaxies. Thus, acriterion to determine whether an orphan satellite is merged withits host is another important aspect to take into account. In the caseof type 2 galaxies, since the halo finder cannot follow its evolution,if we want to estimate the time that orphan satellites remain in thesimulation before merging we need to rely on alternative methods.One possibility, is to use an analytical expression to estimatethis using a dynamical friction timescale 𝑡 df (Binney & Tremaine2008). Then, an orphan satellite galaxy remains in the simulation aslong as the time it has orbited around its central galaxy is less than 𝑡 df . Instead of using an analytical estimate, in this paper, following(Boylan-Kolchin et al. 2008), we assume that a subhalo that is notlonger detected (and its corresponding galaxy) merges with its hosthalo when it loses 99 percent of its initial angular momentum.When that condition is satisfied, the subhalo is considered merged.In addition to this condition on angular momentum, we will alsodefine a proximity criterion for mergers: we consider a subhalo tobe merged when the subhalo-host distance is smaller than a fraction 𝑓 of the virial radius of the main system; i.e. if 𝑟 sat < 𝑓 𝑅 host , (7)where we treat 𝑓 as a free parameter of the model. Therefore, weconsider a subhalo to be merged with its host halo when any of theabove criteria is met; either the subhalo loses 99 percent of its initialangular momentum, or it gets very close to the centre of the hosthalo, as determined by equation 7. Once the subhalo of a galaxy is no longer detected by the halo finderalgorithm, we cannot follow its position evolution and we flag itas a type 2 galaxy (orphan). From that moment on, we integrateits orbit numerically taking as initial conditions the last knownvalues of position, velocity, mass and radius. The subhalo orbitbetween two snapshots is divided into time intervals of length 𝛿𝑡 .Here, the time spacing between snapshots is given by the base DMsimulation. At each time interval, we compute the forces that actover the subhalo according to equations 1, 3 and 4, and perform theevolution of positions and velocities using a kick-drift-kick (KDK)leapfrog scheme.Between snapshots, we took outputs of the orbit evolution ofthe subhalo at time intervals given by Δ 𝑡 ( Δ 𝑡 > 𝛿𝑡 ). For each outputwe compute the tidal radius 𝑟 t (equation 5). If 𝑟 t is smaller than thecurrent radius of the subhalo, we assume that a certain amount ofmaterial that is outside 𝑟 t is unbound and can be removed by tidalstripping. At this stage, we remove an amount of material equivalentto Δ 𝑡 / 𝑇 strip = 𝛼 Δ 𝑡 / 𝑇 orb from the unbound mass, that is Δ 𝑀 sat = 𝛼𝑀 sat ( 𝑟 > 𝑟 t ) Δ 𝑡𝑇 orb (8)where 𝑀 sat ( 𝑟 > 𝑟 t ) is the mass of the satellite outside a radius 𝑟 t MNRAS000
1; on the other hand, Peñarrubia & Benson (2005) assume aninstantaneous stripping , which effectively implies 𝛼 → ∞ . Finally,some authors (e.g. Zentner et al. 2005, Pullen et al. 2014) varythe value of 𝛼 in order to reproduce the halo mass function ofsimulations. In this paper, we will follow the latter approach. According to hierarchical structure formation models, mergers playa critical role in the formation and evolution of galaxies. Thus, acriterion to determine whether an orphan satellite is merged withits host is another important aspect to take into account. In the caseof type 2 galaxies, since the halo finder cannot follow its evolution,if we want to estimate the time that orphan satellites remain in thesimulation before merging we need to rely on alternative methods.One possibility, is to use an analytical expression to estimatethis using a dynamical friction timescale 𝑡 df (Binney & Tremaine2008). Then, an orphan satellite galaxy remains in the simulation aslong as the time it has orbited around its central galaxy is less than 𝑡 df . Instead of using an analytical estimate, in this paper, following(Boylan-Kolchin et al. 2008), we assume that a subhalo that is notlonger detected (and its corresponding galaxy) merges with its hosthalo when it loses 99 percent of its initial angular momentum.When that condition is satisfied, the subhalo is considered merged.In addition to this condition on angular momentum, we will alsodefine a proximity criterion for mergers: we consider a subhalo tobe merged when the subhalo-host distance is smaller than a fraction 𝑓 of the virial radius of the main system; i.e. if 𝑟 sat < 𝑓 𝑅 host , (7)where we treat 𝑓 as a free parameter of the model. Therefore, weconsider a subhalo to be merged with its host halo when any of theabove criteria is met; either the subhalo loses 99 percent of its initialangular momentum, or it gets very close to the centre of the hosthalo, as determined by equation 7. Once the subhalo of a galaxy is no longer detected by the halo finderalgorithm, we cannot follow its position evolution and we flag itas a type 2 galaxy (orphan). From that moment on, we integrateits orbit numerically taking as initial conditions the last knownvalues of position, velocity, mass and radius. The subhalo orbitbetween two snapshots is divided into time intervals of length 𝛿𝑡 .Here, the time spacing between snapshots is given by the base DMsimulation. At each time interval, we compute the forces that actover the subhalo according to equations 1, 3 and 4, and perform theevolution of positions and velocities using a kick-drift-kick (KDK)leapfrog scheme.Between snapshots, we took outputs of the orbit evolution ofthe subhalo at time intervals given by Δ 𝑡 ( Δ 𝑡 > 𝛿𝑡 ). For each outputwe compute the tidal radius 𝑟 t (equation 5). If 𝑟 t is smaller than thecurrent radius of the subhalo, we assume that a certain amount ofmaterial that is outside 𝑟 t is unbound and can be removed by tidalstripping. At this stage, we remove an amount of material equivalentto Δ 𝑡 / 𝑇 strip = 𝛼 Δ 𝑡 / 𝑇 orb from the unbound mass, that is Δ 𝑀 sat = 𝛼𝑀 sat ( 𝑟 > 𝑟 t ) Δ 𝑡𝑇 orb (8)where 𝑀 sat ( 𝑟 > 𝑟 t ) is the mass of the satellite outside a radius 𝑟 t MNRAS000 , 1–14 (2020) rphan galaxies in SAMs and 𝑇 orb is the same as in equation 6. Then, we update the massand radius to the new values. In this procedure, we assume that theoriginal density profile is not altered by the stripping process. Nostripping takes place if 𝑟 t is greater than the original radius of thesatellite halo. We continue evolving the subhalo of an orphan galaxyuntil a merger with its host eventually occurs, following the criteriaexplained in the previous subsection (see equation 7). If the subhalomass is reduced below a threshold limit, the satellite is considereddisrupted. In the previous section, we introduced a model for the evolution ofthe orbits of orphan galaxies in SAG. This model depends on freeparameters that have to be determined by external constraints. Inthis work, we propose to use the halo mass function (HMF) andthe two-point correlation function (2PCF) as constraints for the freeparameters of the orbit model, ( 𝑏, 𝑓 , 𝛼 ) , introduced in equations 4,6 and 7, respectively.We apply our model to two dark matter only N-body simu-lations, with the same cosmological parameters but different massand force resolutions. It is worth noting that in the higher resolu-tion simulation, the majority of halos that disappear in the lowerresolution simulation, leaving an orphan galaxy, will be detectedby the halo finder. The results of the orbit model applied to thelow resolution simulation can be compared with those obtainedfrom the high resolution one, which is considered as the referencesimulation. Then, we vary the free parameters of the model, un-til we find a combination for which the low resolution simulationconverges to the high resolution one, both for the halo mass func-tion ( 𝜙 = 𝑑𝑛 / 𝑑 log 𝑀 , with 𝑛 the numerical density) and for thetwo-point correlation function ( 𝜉 ).Below, we describe the set of dark matter only simulationsused in this work and the computation of the 2PCF. In order to findoptimal parameter values, we need to run the model several times.Since running the code over the full boxes is a numerically expensivetask, we make a parameter exploration over a set of calibrationboxes, which are selected subvolumes of the full simulations with2PCF and HMF similar to those of the full boxes. The details of thisprocedure is covered in the rest of this section. mdpl2 and smdpl simulations In this subsection, we describe the two N-body simulations withdifferent mass and force resolutions that we use for our study ofthe evolution of the orbits of orphan satellites. smdpl and mdpl2are dark matter only N-body simulations that are part of the Multi-Dark cosmological simulation suite . These simulations follow theevolution of 3840 particles within a box, and are characterisedby Planck cosmological parameters: Ω m = . Ω Λ = . Ω b = . 𝑛 s = .
96 and 𝐻 = ℎ km s − Mpc − , where ℎ = .
678 (Planck Collaboration et al. 2016). Both simulationshave been carried out with l-gadget-2 code, a version of the pub-licly available code gadget-2 (Springel et al. 2001a; Springel 2005)whose performance has been optimised for simulating large num-bers of particles. smdpl simulation has a box size of 400 ℎ − Mpcwhich implies a particle mass of 9 . × ℎ − M (cid:12) , while mdpl2simulation has a box size of 1 ℎ − Gpc which implies dark matter particles of 1 . × ℎ − M (cid:12) . Table 1 shows the numerical andcosmological parameters for the simulations. For more details aboutthis set of cosmological simulations see Klypin et al. (2016).These simulations were analysed with the rockstar halo finder(Behroozi et al. 2013a), and merger trees were constructed usingconsistent-trees (Behroozi et al. 2013b). The virial mass of thesestructures is defined as the mass enclosed by a sphere of radius 𝑅 vir , so that the mean density is equal to Δ =
200 times the criticaldensity of the universe 𝜌 c , i.e. 𝑀 vir = / 𝜋𝑅 Δ 𝜌 c . Dark matterhaloes can exist over the background density or lie within anotherdark matter halo. To differentiate them, the former are referred to asmain host haloes, whereas the latter are called subhaloes or satellitehaloes.Figure 1 shows the HMF of the full sample of haloes forthe smdpl ( 𝜙 SM in solid line) and mdpl2 ( 𝜙 MD in dashed line)simulations at redshift 𝑧 =
0. The halo mass function for the mdpl2simulation presents a break at a mass of ∼ . ℎ − M (cid:12) , whichestablishes the minimum mass from which we can guarantee thatwe have completeness in the number of haloes for both simulations(vertical dashed line in Figure 1). In the lower panel, the fractionaldifference between the HMF of the mdpl2 simulation with respectto the one of the smdpl simulation is shown. In the mass range10 . − . ℎ − M (cid:12) , the fractional difference is of the orderof 0 .
2, hence there are 20 percent more low mass haloes in smdplas compared to mdpl2. On the other hand, for halo masses greaterthan 10 . ℎ − M (cid:12) , the fractional difference is always below 0 . ∼ . ℎ − M (cid:12) , while this calculation for smdplgives a minimum halo mass of ∼ ℎ − M (cid:12) . Since completenessis important when comparing the two point correlation functionsof different simulations, in this paper we will consider only haloeswith masses greater than 10 . ℎ − M (cid:12) for both simulations (seevertical dashed line in Figure 1). Given a set of points, the probability of finding an object in aninfinitesimal volume 𝑑𝑉 is, 𝑑𝑃 = 𝑛𝑑𝑉 , where 𝑛 is the mean numberdensity and 𝑁 = 𝑛𝑉 is the number of objects in a finite volume 𝑉 . Then, the 2PCF is defined as the excess probability of findingone of them inside a small volume 𝑑𝑉 and the other in a smallvolume 𝑑𝑉 , separated by a distance 𝑟 (Peebles 1980; Martínez &Saar 2001) 𝑑𝑃 = 𝑛 ( + 𝜉 ( 𝑟 )) 𝑑𝑉 𝑑𝑉 . (9)In practice, for a catalogue of 𝑁 particles and volume 𝑉 with periodicboundary conditions, the 2PCF can be estimated by counting thenumber of pairs of objects, 𝑁 s ( 𝑟 ) , in a shell of volume, 𝑉 s ( 𝑟 ) , atdistance 𝑟 from each other using 𝜉 ( 𝑟 ) = 𝑛 𝑉 𝑁 s ( 𝑟 ) 𝑉 s ( 𝑟 ) − . (10) MNRAS , 1–14 (2020)
Delfino et al. simulation box 𝑁 p 𝑚 p 𝜖 Ω m Ω b Ω Λ 𝜎 𝑛 s 𝐻 ℎ − Gpc ℎ − M (cid:12) ℎ − kpc km s − Mpc − MDPL2 1.0 3830 . × . .
307 0 .
048 0 .
693 0 .
829 0 .
96 67 . . × . .
307 0 .
048 0 .
693 0 .
829 0 .
96 67 . Table 1.
Numerical and cosmological parameters for the N-body simulations. − − − φ [ h M p c − ] SMDPLMDPL2 M h [ h − M (cid:12) ] − . − . . ( φ M D − φ S M ) / φ S M Figure 1.
Upper panel:
Halo mass function of mdpl2 (red dashed line) andsmdpl (black solid line) simulations for redshift 𝑧 =
0. The vertical linedenotes halo masses of 10 . ℎ − M (cid:12) . Note that for masses below thislimit, mdpl2 simulation presents a break which establishes the minimummass from which we can guarantee that we have completeness in the numberof haloes for both simulations. Lower panel:
Fractional difference betweenmdpl2 and smdpl halo mass functions. The dashed horizontal line indicatesa fractional difference of 0 .
05. Note that for masses below than 10 ℎ − M (cid:12) we have approximately 20 percent more haloes in smdpl than in mdpl2. However, if the boundaries are not periodic, this equation cannot beapplied. This is the case for a real survey, or when we restrict theanalysis to a small region within a large periodic simulation. In thesecases, the common method consists in comparing the observed data(D), with catalogues composed of random points (R) that reproducethe same geometry and artefacts of the original catalogue. In general,random catalogues contain at least 10 times more objects than thedata catalogue, in order to reduce the noise level. In such cases, the2PCF can be computed using the Landy-Szalay estimator (Landy& Szalay 1993) 𝜉 ( 𝑟 ) = 𝐷𝐷 ( 𝑟 ) − 𝐷𝑅 ( 𝑟 ) + 𝑅𝑅 ( 𝑟 ) 𝑅𝑅 ( 𝑟 ) ; (11) here 𝐷𝐷 ( 𝑟 ) , 𝐷𝑅 ( 𝑟 ) and 𝑅𝑅 ( 𝑟 ) are the normalized data-data, data-random and random-random pairs, respectively. If 𝑁 d and 𝑁 r arethe number of objects in D and R then 𝐷𝐷 ( 𝑟 ) = 𝑑𝑑 ( 𝑟 ) 𝑁 d ( 𝑁 d − )/ , (12) 𝑅𝑅 ( 𝑟 ) = 𝑟𝑟 ( 𝑟 ) 𝑁 r ( 𝑁 r − )/ , (13) 𝐷𝑅 ( 𝑟 ) = 𝑑𝑟 ( 𝑟 ) 𝑁 d 𝑁 r , (14)where 𝑑𝑑 ( 𝑟 ) is the number of objects pairs separated by a distance 𝑟 in D, 𝑟𝑟 ( 𝑟 ) is the number of pairs separated by a distance 𝑟 inR and 𝐷𝑅 ( 𝑟 ) is the cross-correlation statistic, the number of pairsseparated by a distance 𝑟 with one point taken from D and theother from R (for more details see e.g. Vargas-Magaña et al. 2013).To compute correlation functions, throughout this work we use thepublicly available python package corrfunc (Sinha & Garrison2020).Figure 2 shows the 2PCF for halo masses greater than10 . ℎ − M (cid:12) for the smdpl (solid line) and mdpl2 (dashed line)simulations at redshift 𝑧 =
0. From this figure, we see that the clus-tering of smdpl is greater than that of mdpl2 for all scales, this effectis more significant at lower scales (between 0 . − . ℎ − Mpc). Thesuppression observed in the amplitude of the correlation functionof mdpl2 could be due to a low force resolution. In cosmologicalsimulations, the resolution of the gravitational force is defined bythe value of the softening length 𝜖 . However, there are studies (seee.g. Jenkins et al. 1998) that show that the softening length intro-duces considerable suppression on the clustering signal only forseparations lower than 2 𝜖 . In our case, mdpl2 is characterised by 𝜖 = . ℎ − Mpc, then for scales greater than 0 . ℎ − Mpc theeffect of the softening length should be very small. Thus, the soften-ing length can not account for the large clustering suppression seenin mdpl2 simulation for separations between 10 − − − ℎ − Mpc.According to the halo model, the two-point correlation functioncan be decomposed into two contributions: a 1-halo term and a 2-halo term. The 1-halo term involves correlations between haloes be-longing to the same system, i.e. correlations between central haloesand their corresponding satellites and correlations between all satel-lites that belong to a system. On the other hand, the 2-halo terminvolves correlations between haloes belonging to different systems(Cooray & Sheth 2002). When the contributions from these twoterms are added together, the resulting correlation function shouldroughly follow a power law (see eg. Coil 2013). In general, the 2-halo term dominates at large scales (greater than 1 ℎ − Mpc), whilethe 1-halo term dominates at scales lower than 1 ℎ − Mpc. Notethat the characteristic scale ∼ ℎ − Mpc, which is of the order ofthe size of the main systems, indicates the transition between thetwo regimes. In van den Bosch et al. (2013), it is shown that in-creasing the number of satellites increases mainly the 1-halo term,
MNRAS000
MNRAS000 , 1–14 (2020) rphan galaxies in SAMs ξ ( r ) SMDPL M h ≥ . h − M (cid:12) MDPL2 M h ≥ . h − M (cid:12) − − r [ h − Mpc] − . . ( ξ M D − ξ S M ) / ξ S M Figure 2.
Upper panel:
Two-point correlation function for the smdpl (blacksolid line) and mdpl2 (red dashed line) simulations at redshift 𝑧 = Lowerpanel:
The fractional difference between mdpl2 and smdpl two-point corre-lation functions. For separations greater than 0 . ℎ − Mpc, both simulationspresent a similar clustering. However, for smaller scales, the clustering func-tions differ considerably (fractional difference is greater than 0 . . because subhaloes are located within host systems of typical sizes (cid:46) ℎ Mpc. Therefore, the addition of orphan galaxies, which aresatellites, will enhance the clustering at small scales (Kitzbichler &White 2008).The fraction of satellite haloes with respect to total (satellites+ centrals) for the smdpl simulation is 0 . . . − ℎ − M (cid:12) insmdpl than in mdpl2. Therefore, we conclude that the discrepancyobserved between the two-point correlation functions of these sim-ulations at small scales (see Figure 2) is due to the greater fractionof satellite haloes in smdpl compared to mdpl2. To compensate forthis lack of low-mass subhaloes we introduce the orphan galaxies,integrating their orbits until they merge with their host halo (seeSection 2).These additional satellite galaxies will modify both the HMFand the 2PCF at small scales. Comparing these functions in smdpland mdpl2 simulations, and assuming that their differences couldbe solved by following in a semi-analytic way the dynamics of thosesubhaloes that will host orphan galaxies, we find a way to constrainthe parameters of the orbit model for such satellites. The goal is to make a fast exploration of the parameters of theorbit model, to find regions of the parameter space where there isconvergence between the simulations (i.e., similar HMF and 2PCF)after applying the orbit model. Since running the model over thefull simulations is computationally very expensive, we are inter-ested in finding a set of boxes that are relatively small in volume butrepresentative of the characteristics of the full simulations. We con-sider sub-volumes with the following box sizes: smdpl 50 ℎ − Mpc,mdpl2 50 ℎ − Mpc and mdpl2 100 ℎ − Mpc. In addition, when se-lecting these sub-volumes we need to estimate their “goodness” (inthe sense that they reproduce the characteristics of the full simu-lations). Here, we exemplify how to obtain the best box for mdpl250 ℎ − Mpc; the other cases are analogous.We begin partitioning mdpl2 into 8000 (20 ) disjoint subsam-ples with a box size of 50 ℎ − Mpc. For each of these boxes, wecalculate both 𝜙 and 𝜉 . In the case of the correlation function 𝜉 ,we consider only haloes with masses greater than 10 . ℎ − M (cid:12) ,i.e. where the mdpl2 simulations is complete. Using these results,we compute the mean values of the HMFs and the 2PCFs, 𝜙 and 𝜉 , respectively. Analogously, we obtain the mean values 𝜙 and 𝜉 for both mdpl2 100 ℎ − Mpc (with 10 subvolumes) and smdpl 50(with 8 subvolumes).Figures 3 and 4 show (in circle-dashed line) the mean valuesfor the HMF and 2PCF, respectively. We also plot (in continuosline) the HMF and 2PCF corresponding to the full simulations. Forboth figures, we have the following cases: smdpl 50 ℎ − Mpc (left),mdpl2 50 ℎ − Mpc (centre) and mdpl2 100 ℎ − Mpc (right). Errorbars indicate the standard deviation of the subvolume sets. In Figure3, we note that, as it happens in the mdpl2 full simulation, the boxesalso present a lack of low mass haloes (centre and right). The verticaldashed line indicates the minimum mass we consider in our study(10 . ℎ − M (cid:12) ). Figure 4 shows the clustering signal computedconsidering only haloes with masses greater than 10 . ℎ − M (cid:12) .We note that for very low separations ( (cid:46) . ℎ − Mpc), the scatterin 𝜉 is high, which is reflected in the magnitude of the error bars.This effect is mainly due to the fact that at these scales, the numberof pairs is scarce and therefore Poissonian error dominates.In order to find the best mdpl2 50 ℎ − Mpc box, we need tocharacterise how representative each of our subsamples is. To esti-mate how much each box deviates from the full mdpl2 simulation(regarding to the 2PCF and the HMF), and pre-select some boxesthat can be good candidates, we compute MAPE (mean absolutepercentage error) estimates, i.e., for a given sub-box 𝑗 , we computethe sum over all binsMAPE (cid:104) ℎ ( 𝑗 ) (cid:105) = 𝑛 𝑛 ∑︁ 𝑖 = | ℎ ( 𝑗 ) 𝑖 − ℎ 𝑖 || ℎ 𝑖 | , (15)where ℎ represents either of the relevant constraining functions (i.e.,the 2PCF 𝜉 or the HMF 𝜙 ), 𝑗 indicates the given sub-box, and 𝑛 is thenumber of bins where the function is computed. In this expression, ℎ 𝑖 indicates the value of function ℎ in a given bin 𝑖 correspondingto the full simulations.As mentioned above, the full mdpl2 simulation is not completefor low masses; furthermore, the correlation function computed inboxes becomes too noisy for scales below 0 . ℎ − Mpc. Takingthis into account, for each box we only consider MAPE [ 𝜙 ] errorsfor masses higher than 10 . ℎ − 𝑀 (cid:12) and MAPE [ 𝜉 ] errors forseparations in the range 0 . − ℎ − Mpc. We choose candidatesfor the best mdpl2 50 ℎ − Mpc box as the ones that simultaneously
MNRAS , 1–14 (2020)
Delfino et al. − − φ [ h M p c − ] SMDPLSMDPL 50 h − MpcSMDPL mean + std M h [ h − M (cid:12) ] − . . . ( φ − φ S M ) / φ S M MDPL2MDPL2 50 h − MpcMDPL2 mean + std M h [ h − M (cid:12) ] ( φ − φ M D ) / φ M D MDPL2MDPL2 100 h − MpcMDPL2 mean + std M h [ h − M (cid:12) ] ( φ − φ M D ) / φ M D Figure 3.
HMF for the selected boxes compared to the full simulations and the mean of the boxes. Results are shown for the smdpl 50 ℎ − Mpc (left), mdpl250 ℎ − Mpc (centre) and mdpl2 100 ℎ − Mpc (right).
Upper panels:
The continuos black line correspond to the full simulations. The circle-dashed black lineshows the mean value for the halo mass function computed from the small boxes. The error bars correspond to the standard deviation at each halo mass bin. Redtriangle-dashed lines correspond to the best boxes found.
Lower panels:
Fractional differences in HMF taking the full simulations as a reference. The verticalline denotes the mass cut at ∼ . ℎ − M (cid:12) . In the case of mdpl2 50 ℎ − Mpc, the HMF of the selected box is slightly higher than the mean for all masses. ξ ( r ) SMDPLSMDPL 50 h − MpcSMDPL mean + std − − r [ h − Mpc] − . . . ( ξ − ξ S M ) / ξ S M MDPL2MDPL2 50 h − MpcMDPL2 mean + std − − r [ h − Mpc] ( ξ − ξ M D ) / ξ M D MDPL2MDPL2 100 h − MpcMDPL2 mean + std − − r [ h − Mpc] ( ξ − ξ M D ) / ξ M D Figure 4. ℎ − Mpc (left), mdpl250 ℎ − Mpc (centre) and mdpl2 100 ℎ − Mpc (right).
Upper panels:
The continuous black line shows the mean value for the 2PCF computed from the smallboxes. The error bars correspond to the standard deviation. Red dashed lines correspond to the best boxes found.
Lower panels:
Fractional differences in 2PCFtaking the full simulations as a reference. The vertical dashed line indicates the scale 0 . ℎ − Mpc. This limit corresponds to the smaller scale we use to makethe comparison, since below this limit, boxes have very few pairs and the correlation function becomes too noisy. MNRAS000
Fractional differences in 2PCFtaking the full simulations as a reference. The vertical dashed line indicates the scale 0 . ℎ − Mpc. This limit corresponds to the smaller scale we use to makethe comparison, since below this limit, boxes have very few pairs and the correlation function becomes too noisy. MNRAS000 , 1–14 (2020) rphan galaxies in SAMs minimises both MAPE [ 𝜙 ] and MAPE [ 𝜉 ] . From these candidateswe choose our final best box by visual inspection, prioritizing the2PCF.Figure 3 shows in triangle-dashed red line the HMF of the bestsubvolumes for smdpl 50 ℎ − Mpc (left), mdpl2 50 ℎ − Mpc (cen-ter) and mdpl2 100 ℎ − Mpc (right). The method applied to find thebest subvolumes for smdpl 50 ℎ − Mpc and mdpl2 100 ℎ − Mpc, isanalogous to what we have done for mdpl2 50 ℎ − Mpc. Likewise,Figure 4 shows in triangle-dashed red lines the 2PCF of the best sub-volumes for the same box sizes. These figures show that we are ableto find smaller boxes in which to perform the parameter exploration,that are representative of the complete simulation boxes.
In this section, we make an exploration of the parameter space of theorbit model for orphan galaxies and see how the HMF and the 2PCFchange when these parameters are varied. We recall the meaning ofthe parameters under consideration: 𝑏 characterises the dynamicalfriction through the Coulomb logarithm (equation 4), 𝛼 controls thetidal stripping (equation 6) and 𝑓 is related to the merger criterion(equation 7).Figure 5 shows the results obtained after applying the orbitmodel for orphan satellites to the mdpl2 50 ℎ − Mpc sub-volumeintroduced in Section 3. Here we plot fractional differences takingthe full smdpl simulation as a reference, and selecting all haloeswith masses greater than 10 . ℎ − M (cid:12) at redshift 𝑧 =
0. Thetop panels show fractional differences in HMF while the bottompanels show fractional differences in 2PCF. The solid (black) lineindicates the fractional difference between mdpl2 and smdpl fullsimulations, while the dotted (black) line indicates the fractionaldifference for the mdpl2 50 ℎ − Mpc box. The dashed lines withsymbols correspond to fractional differences considering differentparameter combinations of the model. In all panels, the dashed-circle (black) line correspond to the same parameters ( 𝑏 = . , 𝑓 = . , 𝛼 = . ) . From left to right we vary 𝑏 (left panels), 𝑓 (centrepanels) and 𝛼 (right panels), in each case we leave the remainingparameters fixed. 𝑏 The parameter 𝑏 enters the dynamical friction force term throughthe Coulomb logarithm (see equations 1 and 4). Decreasing 𝑏 isequivalent to increasing the value of the Coulomb logarithm ln Λ andthis leads to a greater deceleration of subhaloes due to dynamicalfriction drag. Since decreasing 𝑏 leads to a greater decelerationforce, then we have a greater number of subhalo mergers and alower number of total subhaloes. On the other hand, increasing 𝑏 decreases the value of ln Λ , we have smaller deceleration and thesubhaloes decay at a slower rate. Thus, if we increase 𝑏 we have alower number of mergers and a greater number of total subhaloes.As we have seen in the previous paragraph, decreasing 𝑏 resultsin a smaller number of total subhaloes, and 𝜙 decreases over theentire mass range. On the other hand, if we increase 𝑏 we have agreater number of total subhaloes and 𝜙 increases. The HMF fordifferent values of 𝑏 , it is shown in the upper-left panel of Figure 5.The 2PCF for different values of 𝑏 is shown in dashed lines in thelower-left panel of Figure 5. Since increasing 𝑏 results in a lowernumber of total satellite haloes, we have a smaller fraction of satellitehaloes and the clustering signal decreases mainly at low scales. Ifwe decrease 𝑏 , we obtain a greater fraction of satellite haloes and the clustering increases at low scales. Note that since decreasing 𝑏 implies that the satellite haloes decay faster, one might think thereshould be an increase in clustering at lower scales. However, thisdoes not occur because most of these haloes end up merging withits main system and therefore do not contribute to the 2PCF. 𝑓 The parameter 𝑓 is related to the merger criterion, increasing (de-creasing) the value of 𝑓 reduces (increases) the distance at whichsatellites merge with their hosts. Therefore, a higher value of 𝑓 im-plies shorter merging times, a greater number of subhalo mergersand a lower number of total satellite haloes. A lower value of 𝑓 implies larger merging times, a smaller number of mergers and ahigher number of total satellite haloes.The central panels of Figure 5 shows, in dashed lines, theresults of the satellite model for different values of parameter 𝑓 . Aswe have seen, a higher value of 𝑓 implies a lower number of satelliteshaloes and the HMF decreases, while a lower value of 𝑓 impliesa higher number of total satellite haloes and the HMF increases.This is shown in the upper-central panel of Figure 5. For the 2PCF,we note that increasing (decreasing) 𝑓 decreases (increases) thevalue of the correlation function 𝜉 . This is shown in the lower-central panel of Figure 5. This effect is more relevant for separationsbelow log ( 𝑟 / ℎ − Mpc ) ≤ − ( 𝑟 / ℎ − Mpc ) ≥ −
1. Since increasing 𝑓 decreases the numberof satellite haloes, then we have a stronger clustering signal at lowscales than for intermediate scales, because the 1-halo term of 𝜉 depends strongly on the satellite fraction. 𝛼 The parameter 𝛼 controls the rate at which the tidal stripping mech-anism removes mass from a subhalo. A greater value of 𝛼 impliesa more efficient tidal stripping, and more material is removed viaTS. A lower value of 𝛼 slows down the TS effect and we have lessmaterial removed from a satellite.The top-right panel of Figure 5 shows (in dashed lines) theHMF for different values of 𝛼 . As we increase 𝛼 , the HMF decreasesover the entire range of masses. Since 𝛼 controls the rate at whichthe tidal stripping mechanism removes mass from a subhalo, thena greater value of 𝛼 implies a more efficient TS, and 𝜙 decreases.Finally, the bottom-right panel of Figure 5 shows (in dashed lines)the 2PCF for different values of 𝛼 . Since increasing 𝛼 increases theefficiency of TS, this reduces the fraction of satellite haloes andthen the clustering signal decreases. To find the best fit parameters, we perform an exploration of theparameter space running the model over the calibration box mdpl250 ℎ − Mpc varying the parameters 𝑏 , 𝑓 and 𝛼 . For parameter 𝑏 ,we consider values in the range 0 . − .
50, we allow 𝑓 to varybetween 0 . − .
08 and 𝛼 takes values between 1 . − .
0. For eachrun output, we compute the corresponding HMF and the 2PCF atredshift 𝑧 = ( 𝑏 = . , 𝑓 = . , 𝛼 = . ) . MNRAS , 1–14 (2020) Delfino et al. . . . − . − . . . . ( φ − φ S M ) / φ S M b = 0 . f = 0 . α = 1 . b = 0 . f = 0 . α = 1 . b = 0 . f = 0 . α = 1 . . . . h /h − M (cid:12) ) b = 0 . f = 0 . α = 1 . b = 0 . f = 0 . α = 1 . b = 0 . f = 0 . α = 1 . . . . b = 0 . f = 0 . α = 2 . b = 0 . f = 0 . α = 1 . b = 0 . f = 0 . α = 5 . − . − . − . . − . − . . . . ( ξ − ξ S M ) / ξ S M − . − . − . . r/h − Mpc) − . − . − . . Figure 5.
Fractional differences in HMF (top panels) and in 2PCF (bottom panels) with respect to the smdpl full simulation for a sample obtained fromapplying the orbit model of orphan satellites to the mdpl2 50 ℎ − Mpc calibration considering different parameter combinations of the model (dashed lined withsymbols). The dashed-circle (black) line correspond to the same parameters ( 𝑏 = . , 𝑓 = . , 𝛼 = . ) in all panels. Solid and dotted lines correspond tothe fractional difference for mdpl2 and the mdpl2 50 ℎ − Mpc calibration box, respectively.
Left panels: effect of varying the parameter 𝑏 leaving 𝑓 and 𝛼 fixed. Parameter 𝑏 enters the dynamical friction force term through the Coulomb logarithm. Decreasing this parameter increases ln Λ and the DF effect, thusmore satellite merge with its host and 𝜙 decreases. Since we have a smaller satellite fraction, the clustering signal 𝜉 also decreases. Central panels: correspondto different values of the parameter 𝑓 leaving 𝑏 and 𝛼 fixed. The value of 𝑓 determines the minimum distance a satellite halo can be from its host before itis considered as merged with it. Increasing this parameter, implies a higher number of mergers and a lower number of satellite haloes. This also impacts theclustering signal 𝜉 , which decreases. Right panels: show the effect of varying 𝛼 leaving 𝑏 and 𝑓 fixed. Increasing the value of 𝛼 implies greater efficiencyof the TS process and greater mass loss and, consequently, a lower number of satellite haloes. As we have fewer satellites, then 𝜉 and 𝜙 decreases. For moredetails see text. We present the treatment for the evolution of the orbits of orphansatellites and the method used for calibrating the free parametersof the model. Here we discuss in more detail some aspects of ourresults.From the analysis presented above we have that, in general,the inclusion of the subhaloes of orphan galaxies helps to enhancethe HMF over the entire mass range. However, as we can see fromthe panels upper-panels of Figure 5, HMF depends more stronglyon parameter 𝛼 , while varying 𝑏 and 𝑓 has little impact on theoverall shape of the HMF. Since 𝑏 and 𝑓 are related mainly to theevolution of the position of the satellites, thus HMF fails to put atight constraint on those parameters. In particular, note that sincethe dynamical friction decelaration is proportional to the mass ofthe satellite (equation 1), then massive haloes are more sensitive tovariations of the DF model. However, the lack of massive subhaloes makes it difficult to put a strong constraint on parameter 𝑏 usingonly information from the halo mass function.On the other hand, we see that the 2PCF is sensitive to varia-tions of the three parameters (see lower panels of Figure 5), althoughthe behavior is different depending on the scales considered. Forvery small scales ( (cid:46) . ℎ − Mpc, see van den Bosch et al. 2013),the correlation function is dominated by the 1-halo term whichdepends strongly on the satellite fraction, then this region is verysensitive to variations on the dynamical friction model ( 𝑏 ) or themerging criterion ( 𝑓 ). At greater separations ( (cid:38) . ℎ − Mpc), the2-halo term begins to compete with the 1-halo term and the con-straining power of the correlation function is reduced. It is worthnoticing that if the strength of TS is high ( 𝛼 ∼ . 𝑏 and 𝑓 , the correlation function has greatconstraining power for all parameters of the model.As a result of the exploration of the model parameters, wehave found the following best fit: 𝑏 = . , 𝑓 = . , 𝛼 = .
43. For
MNRAS000
MNRAS000 , 1–14 (2020) rphan galaxies in SAMs 𝑏 = .
2, we can estimate the maximum and minimum values ofthe Coulomb logarithm. Since the dynamical friction term can onlyproduce deceleration, we impose that ln
Λ = 𝑟 ≤ 𝑏𝑅 sat , as inequation 4. On the other hand, in our simulations, the virial radius ofhaloes (and subhaloes) are in the range 10 − ℎ − kpc. Assuminga maximum satellite-host distance of the order of the virial radiusof main systems, then we have that, at most, 𝑟 / 𝑅 sat ∼ 𝑅 vir / 𝑅 sat ∼ ( ln Λ ) max ∼ .
2. These estimated values are inagreement with our simulations, where the maximum and minimumvalues for Coulomb logarithm are in the range 0 − . Λ =
2. Other authors choose the value ofthis parameter in order to reproduce the results of numerical simu-lations. For example, Velazquez & White (1999), by studying orbitsin numerical simulations, find values of this parameter between1 −
2. On the other hand, Yang et al. (2020), using a model similarto the one presented here, perform a parameter space explorationfinding values for ln Λ of the order of 1 .
5. Clearly, the maximumfor Coulomb logarithm we have found does not agree with the re-sults of these previous works. However, we remark that those resultswere obtained assuming a constant value for ln Λ , whereas in oursimulations we assume a Coulomb logarithm that varies with thesatellite-host distance.The parameter 𝑓 determines the minimum distance a satellitecan approach before merging with its host, therefore this value givesus an estimate of the size of the central galaxy that will inhabit theDM halo as provided by a SAM. For 𝑓 = .
04, the central galaxysize would be of the order 𝑅 g ∼ . 𝑅 vir . This value is compatiblewith the size-virial radius relation found by Kravtsov (2013). Thisrelation links the half mass radius of a galaxy with the virial radiusof its host according to 𝑟 / ∼ . − . 𝑅 .As we have seen, the parameter 𝛼 controls the efficiency ofTS mechanism. Zentner & Bullock (2003) assumes a fixed value 𝛼 =
1, while other authors adjust 𝛼 in order to match the results ofsimulations finding 𝛼 = . 𝛼 = . 𝛼 = ∞ (Peñarrubia& Benson 2005). In our parameter exploration, we have found that 𝛼 = .
43, which is much smaller than the results obtained by Z05and P14. Finally, our results suggest that a value 𝛼 > 𝛼 , as comparedto other works, could be the fact that in our model, the ln Λ isvariable, whereas in previous works, this value is usually taken asconstant. In our model, we initially have a larger value of ln Λ , andhence, the satellites decay more rapidly and come closer to the host,where the TS is stronger. On the other hand, in the case of fixedln Λ , the initial deceleration is smaller and therefore it takes moretime for the satellite to reach inner regions, where the density andTS are higher. To compensate for this, a higher value of 𝛼 would berequired, in order to have a higher mass loss.We also perform a convergence test using the other calibra-tion boxes described in Section 3.3. We run the model using thebest fit parameters over mdpl2 050 ℎ − Mpc (MD050), mdpl2100 ℎ − Mpc (MD100) and smdpl 50 ℎ − Mpc boxes (SM050).The results of this exercise are shown (in dashed lines with sym-bols) in Figure 6 and Figure 7 for HMF and 2PCF, respectively. Indotted lines with symbols we plot the calibration boxes without theorphan satellites. As a reference, we plot the mdpl2 full simulation(in dashed line) and the smdpl full simulation (in continuous line). In the lower panel of these figures, we show fractional differencestaking the smdpl full simulation as a reference.Figure 6 shows that, after applying the model, we have a goodagreement on the different calibration boxes after including theorphan satellites (dashed lines with symbols), compared to the boxeswithout including the orphans (dotted lines with symbols). From thisfigure, we also notice that smdpl 50 ℎ − Mpc and mdpl2 50 ℎ − Mpcpresent some “spikes” at masses of the order of 10 ℎ − M (cid:12) ; thisis a particular characteristic of the selected boxes. For lower halomasses, the inclusion of orphans in smdpl 50 ℎ − Mpc (in squares)has little impact on the halo mass function. For mdpl2 50 ℎ − Mpc(in circles) we note that the addition of orphan satellites helps toenhance the HMF at low masses, in this case the agreement withthe HMF of smdpl full is within 5 percent over the entire massrange. Finally, for mdpl2 100 ℎ − Mpc (in triangles), we obtain agood agreement (within 5 percent) with the other cases except formasses below 10 ℎ − M (cid:12) where the difference is of the order of10 percent.Figure 7 shows the effect of running the model over the cali-bration boxes using the best fit parameters for the 2PCF. In general,the addition of orphan satellites in mdpl2 50 ℎ − Mpc and mdpl2100 ℎ − Mpc, considerably improves the two point correlation func-tion over the entire range of scales (dasehd lines with symbols),compared to the boxes without the addition of orphans (dotted lineswith symbols). However, this enhancement is more important atsmall scales (0 . − . ℎ − Mpc) where the 1-halo term domi-nates, which mainly depends on the fraction of satellite haloes. Ingeneral, fractional errors for 2PCF are within 10 percent for thethree cases considered. Note that for both mdpl2 subvolumes, theclustering is strongly suppressed with a discrepancy greater than50 percent for scales close to 0 . ℎ − Mpc. From these results, wecan conclude that the inclusion of orphans improves the HMF and2PCF of the simulation with lower resolution (mdpl2) as comparedto smdpl.With this analysis, we have shown that, on the one hand, thesmdpl simulation has a complete population of satellite galaxiesabove 10 . ℎ − M (cid:12) , and can be used as a reference for the ex-ploration of the parameters of the orbit model. Indeed, if we ap-ply the orbit model to a smaller fraction of smdpl, namely smdpl50 ℎ − Mpc, we would obtain the same HMF and 2PCF (see Figure6 and 7). On the other hand, the mdpl2 50 ℎ − Mpc box is suf-ficiently good to test the parameters, and we do not need a largerbox. This can be verified by taking a larger box, namely mdpl2100 ℎ − Mpc, in which the results, when the best fit parameters areused, are consistent with those of the mdpl2 50 ℎ − Mpc box (theobserved differences are due to the selected box).
In this paper, we present an updated model for the orbital evolutionof dark matter subhaloes. The model used includes tidal strippingeffects and dynamical friction. It also takes into account the possiblemerger of the haloes using a proximity criterion. We have charac-terised our model by three free parameters ( 𝑏, 𝑓 , 𝛼 ) . The proposedmodel describes the main processes that affect the orbital evolutionof a satellite halo and it is simple enough to be relatively cheap froma computational point of view (at least two order of magnitude fasterthan a full N-body simulation of the same resolution and size).In order to calibrate the free parameters of the model, satellitehalo evolution studies in the literature used the HMF as a constraint.Figure 5 shows that 𝜙 is sensitive to variations in the efficiency of the MNRAS , 1–14 (2020) Delfino et al. − − φ [ h M p c − ] SM050MD050MD100SM050 + modelMD050 + modelMD100 + modelMDPL2SMDPL M h [ h − M (cid:12) ] − . . . ( φ − φ S M ) / φ S M Figure 6.
Upper panel:
Halo mass functions. The dotted lines with symbolscorrespond to calibration boxes (see Section 3.3), i.e. without the orphansatellites. The dashed lines with symbols correspond to the calibration boxesafter applying the orbit evolution model using the best parameters (see Sec-tion 4.4). Here SM050 correspond to smdpl 50 ℎ − Mpc best box, MD050to mdpl2 50 ℎ − Mpc and MD100 indicates 100 ℎ − Mpc simulations. Thesimple dashed line denotes mdpl2 full and continuous line correspond tothe smdpl full simulation.
Lower panel:
Fractional differences taking thesmdpl full simulation as a reference. We see that after applying the model,for smdpl 50 ℎ − Mpc and mdpl2 50 ℎ − Mpc the difference is below 5percent for the entire mass range. tidal stripping (parameter 𝛼 ). On the other hand, the halo mass func-tion fails to put tight constraints on the dynamical friction model,since this effect is most important for massive satellites, which areless abundant in the simulation (see right panels Figure 5). In thiswork, we have introduced the 2PCF as another constraint to the non-linear evolution of the orbits of orphan satellite haloes. In contrastwith the HMF, the correlation function includes information aboutthe distribution of subhaloes around their host. In this sense, the in-clusion of clustering information helps to constrain the dynamicalfriction model and merging criterion (see left and centre panels ofFigure 5).We have performed an exploration of the parameter space ofthe model, this gives us the following best fit set ( 𝑓 = . , 𝑏 = . , 𝛼 = . ) . The value 𝑓 = .
04 is compatible with the valueof half mass radius relation presented by Kravtsov (2013). Theparameter 𝑏 gives us a variable Coulomb logarithm which takesminimum - maximum values in the range ln Λ ∼ − .
2. Otherworks in the literature usually use ln
Λ = −
2; since these valuesare taken to be constant, we cannot compare them directly with ourresults. On the other hand, 𝛼 = .
43, the parameter that controls the ξ ( r ) SM050MD050MD100SM050 + modelMD050 + modelMD100 + modelSMDPLMDPL2 − r [ h − Mpc] − . . ( ξ − ξ S M ) / ξ S M Figure 7.
Upper panel:
Two-point correlation functions. The dotted lineswith symbols correspond to the calibration boxes, i.e. without the orphansatellites. The dashed lines with symbols correspond to the samples obtainedfrom the best boxes after applying the orbit model using the best parame-ters. Here, SM050 correspond to SMDPL 50 ℎ − Mpc best box, MD050 toMDPL2 50 ℎ − Mpc and MD100 indicates MDPL2 100 ℎ − Mpc satellites.The simple dashed line denotes mdpl2 full and continuous line correspondto the smdpl full simulation.
Lower panel: fractional differences in theHMFs taking smdpl full as a reference. After applying the model, we ob-tain an enhancement of the clustering for MDPL2 50 ℎ − Mpc and MDPL2100 ℎ − Mpc. In general, the agreement between the different calibrationboxes is within 10 percent except for scales below 0 . ℎ − Mpc, where theagreement is within 20 percent. efficiency of the TS mechanism, is lower compared to the valuesobtained in previous works. This might be related to the fact that weare assuming a variable Coulomb logarithm instead of a constantone, but this point requires further study.Finally, we remark that we have not included a tidal heatingterm in our orbital evolution model. Tidal heating is a mechanismproduced by rapid changes in tidal forces when a halo passes throughthe pericentre. These tidal shocks transfer energy to the satellite,then the subhalo expands and therefore a greater amount of matteris susceptible to be removed via tidal stripping. In a future work,we plan to include tidal heating effects in the orbital evolution ofsubhaloes and study individual orbits, taking special care on theevolution of masses and radius.
ACKNOWLEDGEMENTS
FMD and CGS are supported by the National Agency for thePromotion of Science and Technology (ANPCYT) of Argentina
MNRAS000
MNRAS000 , 1–14 (2020) rphan galaxies in SAMs grant PICT-2016-0081; and grants G140 and G157 from UNLP.CVM acknowledges support from ANID/FONDECYT throughgrant 3200918. SAC acknowledges funding from Consejo Nacionalde Investigaciones Científicas y Técnicas (CONICET, PIP-0387),
Agencia Nacional de Promoción de la Investigación, el DesarrolloTecnológico y la Innovación (Agencia I+D+i, PICT-2018-03743),and
Universidad Nacional de La Plata (G11-150), Argentina.The CosmoSim database used in this paper is a service by theLeibniz-Institute for Astrophysics Potsdam (AIP). The MultiDarkdatabase was developed in cooperation with the Spanish MultiDarkConsolider Project CSD2009-00064.
DATA AVAILABILITY
The data underlying this article will be shared on reasonable requestto the corresponding author. The simulations used in this paper arepublicly available at . REFERENCES
Amendola L., et al., 2013, Living Reviews in Relativity, 16, 6Baugh C. M., 2006, Reports on Progress in Physics, 69, 3101Behroozi P. S., Conroy C., Wechsler R. H., 2010, ApJ, 717, 379Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013a, ApJ, 762, 109Behroozi P. S., Wechsler R. H., Wu H.-Y., Busha M. T., Klypin A. A.,Primack J. R., 2013b, ApJ, 763, 18Benson A. J., 2010, Physics Reports, 495, 33Benson A. J., 2012, New Astron., 17, 175Berlind A. A., Weinberg D. H., 2002, ApJ, 575, 587Berlind A. A., et al., 2003, ApJ, 593, 1Binney J., Tremaine S., 2008, Galactic Dynamics: Second EditionBoylan-Kolchin M., Ma C.-P., Quataert E., 2008, MNRAS, 383, 93Brüggen M., De Lucia G., 2007, Monthly Notices of the Royal AstronomicalSociety, 383, 1336Chandrasekhar S., 1943, ApJ, 97, 255Coil A. L., 2013, The Large-Scale Structure of the Universe. p. 387,doi:10.1007/978-94-007-5609-0_8Conroy C., Wechsler R. H., Kravtsov A. V., 2006, ApJ, 647, 201Cooray A., Sheth R., 2002, Phys. Rep., 372, 1Cora S. A., 2006, MNRAS, 368, 1540Cora S. A., Muzzio J. C., Vergne M. M., 1997, MNRAS, 289, 253Cora S. A., et al., 2018, ArXiv e-prints: 1801.03883,Crain R. A., et al., 2015, MNRAS, 450, 1937Croton D. J., et al., 2006, MNRAS, 365, 11DESI Collaboration et al., 2016, arXiv e-prints, p. arXiv:1611.00036Dawson K. S., et al., 2015, preprint, ( arXiv:1508.04473 )Dawson K. S., et al., 2016, AJ, 151, 44Euclid Collaboration et al., 2019, arXiv e-prints, p. arXiv:1910.09273Frenk C. S., White S. D. M., 2012, Annalen der Physik, 524, 507Gan J., Kang X., van den Bosch F. C., Hou J., 2010, Monthly Notices of theRoyal Astronomical Society, 408, 2201Gargiulo I. D., et al., 2015, MNRAS, 446, 3820Gonzalez-Perez V., Lacey C. G., Baugh C. M., Lagos C. D. P., Helly J.,Campbell D. J. R., Mitchell P. D., 2014, MNRAS, 439, 264Gunn J. E., Gott J. Richard I., 1972, ApJ, 176, 1Guo Q., White S., 2014, MNRAS, 437, 3228Hashimoto Y., Funato Y., Makino J., 2003, ApJ, 582, 196Henriques B. M. B., White S. D. M., Thomas P. A., Angulo R. E., Guo Q.,Lemson G., Springel V., 2013, MNRAS, 431, 3373Jenkins A., et al., 1998, ApJ, 499, 20King I., 1962, AJ, 67, 471Kitzbichler M. G., White S. D. M., 2008, Monthly Notices of the RoyalAstronomical Society, 391, 1489 Klypin A., Yepes G., Gottlöber S., Prada F., Heß S., 2016, MNRAS, 457,4340Kravtsov A. V., 2013, ApJ, 764, L31Kravtsov A. V., Berlind A. A., Wechsler R. H., Klypin A. A., Gottlöber S.,Allgood B. o., Primack J. R., 2004, ApJ, 609, 35LSST Dark Energy Science Collaboration 2012, arXiv e-prints, p.arXiv:1211.0310LSST Science Collaboration et al., 2009, arXiv e-prints, p. arXiv:0912.0201Lagos C. D. P., Cora S. A., Padilla N. D., 2008, MNRAS, 388, 587Landy S. D., Szalay A. S., 1993, ApJ, 412, 64Laureijs R., et al., 2011, preprint, ( arXiv:1110.3193 )Levi M., et al., 2013, preprint, ( arXiv:1308.0847 )Łokas E. L., Mamon G. A., 2001, MNRAS, 321, 155Marinacci F., et al., 2018, MNRAS, 480, 5113Martínez V. J., Saar E., 2001, Statistics of the Galaxy Distribution. Chapmanand Hall/CRC, Boca Raton, FL, USAMoster B. P., Somerville R. S., Maulbetsch C., van den Bosch F. C., MacciòA. V., Naab T., Oser L., 2010, ApJ, 710, 903Muñoz Arancibia A. M., Navarrete F. P., Padilla N. D., Cora S. A., GawiserE., Kurczynski P., Ruiz A. N., 2015, MNRAS, 446, 2291Naiman J. P., et al., 2018, MNRAS, 477, 1206Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493Nelson D., et al., 2018, MNRAS, 475, 624Ogiya G., Burkert A., 2016, MNRAS, 457, 2164Orsi Á., Padilla N., Groves B., Cora S., Tecce T., Gargiulo I., Ruiz A., 2014,MNRAS, 443, 799Peñarrubia J., Benson A. J., 2005, MNRAS, 364, 977Peacock J. A., Smith R. E., 2000, MNRAS, 318, 1144Peebles P. J. E., 1980, The large-scale structure of the universePetts J. A., Gualandris A., Read J. I., 2015, MNRAS, 454, 3778Petts J. A., Read J. I., Gualandris A., 2016, MNRAS, 463, 858Pillepich A., et al., 2018, MNRAS, 475, 648Planck Collaboration et al., 2016, A&A, 594, A13Pujol A., et al., 2017, MNRAS, 469, 749Pullen A. R., Benson A. J., Moustakas L. A., 2014, ApJ, 792, 24Reddick R. M., Wechsler R. H., Tinker J. L., Behroozi P. S., 2013, ApJ, 771,30Roukema B. F., Quinn P. J., Peterson B. A., Rocca-Volmerange B., 1997,MNRAS, 292, 835Schaye J., et al., 2015, MNRAS, 446, 521Silk J., Mamon G. A., 2012, RAA, 12, 917Sinha M., Garrison L. H., 2020, MNRAS, 491, 3022Springel V., 2005, MNRAS, 364, 1105Springel V., Yoshida N., White S. D. M., 2001a, New Astron., 6, 79Springel V., White S. D. M., Tormen G., Kauffmann G., 2001b, MNRAS,328, 726Springel V., et al., 2018, MNRAS, 475, 676Taylor J. E., Babul A., 2001, ApJ, 559, 716Tecce T. E., Cora S. A., Tissera P. B., Abadi M. G., Lagos C. D. P., 2010,MNRAS, 408, 2008The LSST Dark Energy Science Collaboration et al., 2018, arXiv e-prints,p. arXiv:1809.01669Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G.,Gottlöber S., Holz D. E., 2008, ApJ, 688, 709Trujillo-Gomez S., Klypin A., Primack J., Romanowsky A. J., 2011, ApJ,742, 16Vale A., Ostriker J. P., 2004, MNRAS, 353, 189Vargas-Magaña M., et al., 2013, A&A, 554, A131Velazquez H., White S. D. M., 1999, MNRAS, 304, 254Vogelsberger M., Marinacci F., Torrey P., Puchwein E., 2020, Nature Re-views Physics, 2, 42Vollmer B., Cayatte V., Balkowski C., Duschl W. J., 2001, ApJ, 561, 708Wechsler R. H., Tinker J. L., 2018, ARA&A, 56, 435Weinberg M. D., 1986, ApJ, 300, 93White S. D. M., Rees M. J., 1978, MNRAS, 183, 341Yang S., Du X., Benson A. J., Pullen A. R., Peter A. H. G., 2020, MNRAS,498, 3902Zentner A. R., Bullock J. S., 2003, ApJ, 598, 49MNRAS , 1–14 (2020) Delfino et al.
Zentner A. R., Berlind A. A., Bullock J. S., Kravtsov A. V., Wechsler R. H.,2005, ApJ, 624, 505eBOSS Collaboration et al., 2020, arXiv e-prints, p. arXiv:2007.08991van den Bosch F. C., More S., Cacciato M., Mo H., Yang X., 2013, MNRAS,430, 725
APPENDIX A: INGREDIENTS OF THE ORBIT MODELFOR A NFW DENSITY PROFILE
In numerical studies, the detected dark matter haloes are well de-scribed by a double power-law function of the radius, first introducedby (Navarro et al. 1997). Here we describe the Navarro-Frenk-White(NFW) density profile, and we give expressions for the dynamicalfriction formula and tidal radius equation corresponding to this par-ticular profile. The Navarro-Frenk-White is given by the followingexpression 𝜌 ( 𝑟 ) = 𝜌 𝛿 char 𝑟𝑟 s (cid:16) + 𝑟𝑟 s (cid:17) . (A1)The parameters of this expression are the scale radius 𝑟 s and thecharacteristic density 𝛿 char . The concentration parameter is definedas 𝑐 = 𝑟 vir 𝑟 s , (A2)where 𝑟 vir is the virial radius of the halo, defined as the distancefrom the center of the halo within which the mean density is Δ times the critical density 𝜌 . The value of the virial overdensity isoften assumed to be 178, predicted by the spherical collapse model.However, numerical simulations typically use Δ = 𝑠 and the function 𝑔 ( 𝑐 ) , which often appears in calculations involving the NFW profile, 𝑠 = 𝑟𝑟 vir , (A3) 𝑔 ( 𝑐 ) = ( + 𝑐 ) − 𝑐 /( + 𝑐 ) . (A4)With the above definitions, the density equation becomes 𝜌 ( 𝑠 ) 𝜌 = Δ 𝑐 𝑔 ( 𝑐 ) 𝑠 ( + 𝑐𝑠 ) . (A5)The mass of the halo is usually defined as the mass within the virialradius 𝑀 vir = 𝜋𝑟 Δ 𝜌 𝑐 (A6)Then the mass profile in units of the virial mass is 𝑀 ( 𝑠 ) 𝑀 vir = 𝑔 ( 𝑐 ) (cid:104) log ( + 𝑐𝑠 ) − 𝑐𝑠 + 𝑐𝑠 (cid:105) (A7)Using the same definitions, we can express the gravitational poten-tial as Φ ( 𝑠 ) 𝑉 = − 𝑔 ( 𝑐 ) log ( + 𝑐𝑠 ) 𝑠 , (A8) where 𝑉 is the circular velocity at 𝑟 = 𝑟 vir .For a point particle moving within a halo with a NFW densityprofile, the 𝑖 -component of the acceleration is a = − 𝑔 ( 𝑐 ) (cid:20) log ( + 𝑐𝑠 ) 𝑠 − 𝑐 ( + 𝑐𝑠 ) 𝑠 (cid:21) s 𝑠 (A9)In this case, the Chandrasekhar formula for the dynamical frictionforce is given by a df = − 𝑀 sat ln Δ c 𝑣 𝑐 𝑔 ( 𝑐 ) 𝑠 ( + 𝑐𝑠 ) (cid:20) erf ( 𝑋 ) − 𝑋 √ 𝜋 exp − 𝑋 (cid:21) v 𝑣 . (A10)Finally, we need to compute the tidal radius, i.e. 𝑟 t = (cid:18) 𝐺 𝑀 sat 𝜔 − 𝑑 Φ / 𝑑𝑟 (cid:19) / . (A11)For 𝑑 Φ / 𝑑𝑟 we have 𝑑 Φ 𝑑𝑟 ( 𝑟 = 𝑟 sat ) = − 𝐺 𝑀 host ( < 𝑟 sat ) 𝑟 + 𝜋𝐺 𝜌 ( 𝑟 sat ) . (A12)Then using the above expressions for 𝑀 host (A7) and 𝜌 host (A5), weobtain the tidal radius. This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000