https://orcid.org/0000-0001-5501-6008
ABSTRACT
We describe a simple model that explains the qualitative and (approximate) quantitative features of the distribution of orbital velocities of merging pairs of dark matter haloes. Our model considers a primary dark matter halo as a perturber in a background of secondary haloes with velocities described by linear theory. By evaluating the ensemble of secondary haloes on orbits within the perturbing halo’s ‘loss cone’ we derive the distribution of orbital parameters of these captured haloes. This model is able provide qualitative explanations for the features of this distribution as measured from N-body simulations, and is in approximate quantitative agreement with those measurements. As the velocity dispersion of the background haloes is larger on smaller scales our model predicts an overall increase in the characteristic velocities of merging haloes, relative to the virial velocities of those haloes, in lower mass systems. Our model also provides a simple explanation for the measured independence of the orbital velocity distribution function on redshift when considered at fixed peak height. By connecting the orbital parameter distribution to the underlying power spectrum our model also allows for estimates to be made of the effect of modifying that power spectrum, for example by including a truncation at large wavenumber. For plausible warm dark matter models, we find that this truncation has only a small effect on the predicted distributions.
1 INTRODUCTION
Understanding the dark matter subhalo population – in particular its distributions in mass, host halo-centric radius, etc. – is of key importance for many studies aimed at constraining the nature of dark matter (Gilman et al. 2020; Bose et al. 2020; Nadler et al. 2021). While N-body and hydrodynamical simulations allow a direct calculation of these distributions, the computational cost of generating large numbers of high-resolution realizations (needed to quantify statistical properties) remains high. Semi-analytic models for the subhalo population (Taylor & Babul 2001; Benson et al. 2002; Zentner et al. 2005; Pullen, Benson & Moustakas 2014; Yang et al. 2020; Jiang et al. 2020), while being more approximate, achieve much lower computational costs. These models solve coupled sets of differential equations which describe the orbital position and velocity, the bound mass of the subhalo, and its density profile.
A key input to these models is the initial conditions for these differential equations. Typically, these are taken to be the orbital parameters of the subhalo at the point where it first crosses the virial radius of its host. Under the usual assumption of an isotropic distribution of merging haloes, these initial conditions are fully specified by the combination of the radial and tangential velocity of the subhalo (or any two equivalent quantities, such as the energy and angular momentum).
Distributions of orbital velocities for merging subhaloes have been measured by several authors (Benson 2005; Wetzel 2011; Jiang et al. 2015; Li et al. 2020). In these works, the distributions were obtained by identifying merging halo pairs in cosmological N-body simulations of cold dark matter, and constructing 2D histograms of their orbital velocities at the virial radius. These works show that radial and tangential velocities are strongly correlated – with the distribution peaking along the ridge line of constant total velocity (Jiang et al. 2015) – with mean values slightly lower than the virial velocity of the host, and with only weak dependence on the host mass, mass ratio, and redshift (when expressed in virial units).
The goal of this work is to gain some insight into the physics that sets the distribution of orbital parameters of merging dark matter haloes. Our motivation is to provide both some understanding of the physical origin of this distribution, and a simple model for how this distribution may be expected to respond to changes in the matter power spectrum (such as might arise in non-cold dark matter models). Given that, by definition, mergers between haloes happen far into the non-linear regime of gravitational structure formation, we do not expect our results to necessarily be quantitatively precise.
Our approach is to consider the more massive (‘primary’) halo embedded in a background of ‘secondary’ haloes which have a cosmological distribution of velocities, and to consider which of those haloes will be gravitationally captured by the primary, i.e. those whose orbits have pericentric distances that lie within the virial radius of the primary – similar to ‘loss cone’ calculations employed for supermassive black holes (Merritt 2013). Evaluating the velocities of these captured secondary haloes at the primary halo virial radius and integrating over all possible such haloes lead to a model of the orbital parameter distribution function.
The remainder of this paper is organized as follows. In Section 2, we describe our method for computing the orbital parameter distribution function. In Section 3, we show results from this method, and compare them to previous determinations of the distribution function. We also explore how the distribution function is expected to change in the low mass halo regime, examine any evolution with redshift, and explore the effects of a cut-off in the matter power spectrum. Finally, in Section 4, we discuss the outcomes and implications of our work.
Throughout this work will adopt a cosmological model described by |$(H_0,\Omega _\mathrm{m},\Omega _\Lambda ,\Omega _\mathrm{b})=(67.66\hbox{km/s},0.307,0.693,0.0482)$| and a matter power spectrum characterized by (σ8, ns) = (0.8228, 0.96) (Planck Collaboration XVI 2014), although our model applies equally well to other cases.
2 METHODS
In the following, we refer to the perturber halo as the ‘primary’ halo, with any halo from the background population referred to as the ‘secondary’ halo. Properties of these haloes are identified using subscripts ‘p’ and ‘s’ respectively. Without any loss of generality, we consider Mp, v ≥ Ms, v, where Mv is the virial mass of the halo, defined as the mass within a sphere enclosing a mean density contrast of Δv as computed under the standard spherical collapse model (e.g. Percival 2005).
Throughout this work we work in virial units, such that all lengths are measured in units of the primary halo virial radius, rp, v, and all velocities in units of the primary halo virial velocity, Vp, v = (GMp, v/rp, v)1/2.
Note that the initial factor of |Vr| in equation (1) arises from the fact that secondary haloes cross the virial radius of the primary halo (i.e. they merge with it) at a rate proportional to their radial velocity at the virial radius (Benson 2005). Note also that we include contributions from haloes which are initially ingoing (i.e. |$V^\prime _\mathrm{r} \lt 0$| and outgoing |$V^\prime _\mathrm{r} \gt 0$|).
The velocity |$(V^\prime _\mathrm{r},V^\prime _\theta)$| is computed from (Vr, Vθ) and r under the assumption that the secondary halo moves in a Keplerian potential2 corresponding to the mass of the primary halo. We exclude from the integral radii r < 1 (i.e. within the virial radius of the primary halo), and those for which the time-of-flight from r to the primary halo virial radius exceeds the age of the Universe (this latter condition has little effect on the results).
Evaluating the above distribution then only requires models for ξps(r) and |$p(V^\prime _\mathrm{r},V^\prime _\theta |r)$|. We assume that ξps(r) = b(Mp)b(Ms)ξlin(r), where b(M) is the bias of haloes of mass M, and ξlin(r) is the linear theory two-point correlation function, which is computed via the Fourier transform of the power spectrum. For the bias, b(M), we use the model of Tinker et al. (2010).
It is well-established that the cosmological pairwise velocity distribution function is, in fact, non-Gaussian, showing extended tails and skewness (Sheth & Diaferio 2001; Scoccimarro 2004) due to non-linear contributions to the pairwise velocities of haloes arising from scales much smaller than the halo separation. The core of the distribution is generally well represented by a Gaussian (e.g. Cuesta-Lazaro et al. 2020), with the non-Gaussian components affecting only the tails. Therefore, for simplicity, we do not consider these non-Gaussian features of the pairwise velocity distribution here.
Using the above, the orbital parameter distribution at the virial radius can be evaluated for any combination of primary and secondary haloes.
3 RESULTS
We use our model to compute the distribution of (Vr, Vθ) for various combinations of primary and secondary halo mass at z = 0. The upper-left panel of Fig. 1 shows the resulting distribution function for Mp = 1012M⊙ and Ms = 1011M⊙. The rainbow colour scale shows the results of our model, while contours show the model of Li et al. (2020) at levels of 0.01, 0.10, 0.30, and 0.60 for comparison. The upper right and lower left-hand panels show the marginal distributions of Vr and Vθ respectively, with results from our model shown by the blue lines, and those of Li et al. (2020) by the green lines. Finally, the lower right panel shows the distribution of Vr/V where |$V=(V_\mathrm{r}^2+V_\theta ^2)^{1/2}$|, again with results from our model shown by the blue lines, and those of Li et al. (2020) by the green lines.
Distributions of orbital parameters for merging haloes of masses Mp = 1012M⊙ and Ms = 1011M⊙ at z = 0. Upper left-hand panel: The joint distribution of radial and tangential velocities of merging haloes. The colour map shows results from our model, while the contours indicate the distribution of Li et al. (2020). Upper right-hand panel: The distribution of radial velocities of merging haloes. The blue curve shows results from our model, while the green line is the distribution of Li et al. (2020). Lower left-hand panel: As the upper right-hand panel, but for the tangential velocity. Lower right-hand panel: As the upper right-hand panel, but for the radial velocity expressed as a fraction of the total velocity.
We find that a value of α = 0.8 in equation (5) gives a reasonable match to the N-body results from Li et al. (2020). Had we instead taken α = 1 the predicted distributions of Vr and Vθ would be slightly broader, and the peak of the Vθ distribution shifted to slightly higher velocity, but the location of peak of the Vr distribution would be largely unchanged. Given the simplifying assumptions of our model, a value of α = 0.8 seems acceptable. There is good qualitative agreement in the locus of our distribution function with that of Li et al. (2020), although in detail they do not agree. There is also good agreement in the marginal distributions. For Vr our model predicts approximately the correct shape and width of the distribution, although it is shifted to slightly too large values. For Vθ there is excellent agreement between our model and the results of Li et al. (2020) in both the peak location, width of the distribution, and overall shape. Similarly, our model accurately matches the form of the distribution of Vr/V found by Li et al. (2020). In Appendix A, we show similar results for Mp = 1014M⊙ and Ms = 1011M⊙, where our model performs less well due to underestimating the dispersion of the secondary halo population.
Given this degree of success in matching measurements of the orbital parameter distribution in N-body simulations, we can utilize our model to provide some insights into the properties of the distribution. For example, the ridge of the distribution of (Vr, Vθ) approximately tracks a line of constant |$V^2=V_\mathrm{r}^2+V_\theta ^2 \equiv V_\mathrm{l} \approx 1.5$|, as was first noted by Jiang et al. (2015). This therefore corresponds to a locus of constant energy (since all orbits are measured at the same radius, the primary halo virial radius, and so at the same potential).
Orbits with infall velocities, V > Vl, in excess of the velocity of this locus are either unbound, or close to unbound. The radial velocity of such orbits remains non-zero out to large distances, unlike orbits with V < Vl which reach zero radial velocity (i.e. apocenter, rapo) at radii of 3–4 virial radii or less. As such, the contribution of these orbits is strongly down-weighted by the distribution function, |$p^\prime (V^\prime _\mathrm{r},V^\prime _\theta)$|, at the evaluation radius as |$V^\prime _\mathrm{r}$| is significantly offset from |$\bar{V}^\prime _\mathrm{r}$|. (For very large radii this effect is further enhanced by the Hubble expansion.) This is the cause in the decline in the distribution for V > Vl.
For V < Vl, the decline in the distribution function is driven by the available volume. These orbits have apocenters within at most a few virial radii, with the apocentric radius decreasing as V decreases. As we evaluate the contribution of each orbit only over the radial range 1(≡ rv) to rapo, there is little volume contributing to these regions of the distribution function.
The overall width of the distribution is controlled by the velocity dispersion of the background haloes (and, therefore, is fundamentally determined by the matter power spectrum). Without the smoothing effect of the velocity dispersion, the distribution would have a much sharper cut off beyond Vl, and similar (but somewhat less sharp) cut off for V < Vl. Artificially increasing/decreasing σ(r) results in a corresponding increase/decrease in the width of the distributions of Vr and Vθ, along with a positive/negative shift in the location of the peak of the Vθ distribution. The location of the peak in the Vr distribution is largely unaffected by such changes in σ(r) – it is instead mostly determined by the arguments given above for the location of the ridge line in the (Vr, Vθ) distribution.
In Fig. 2, we show summary statistics (mean, |$\bar{V}$|, and dispersion, σ, of the radial and tangential orbital velocities) as predicted by our model (solid lines), and as measured by Li et al. (2020; dashed lines). For the Li et al. (2020) results, we plot lines only over the range of halo masses considered in that work. Over the range of halo masses where Li et al. (2020) measured the orbital velocity distribution our model is in reasonable agreement with their results. Notably, over the range of primary halo masses considered by Li et al. (2020) our model predicts a weak dependence of these summary statistics on Mp, although not as weak as found by Li et al. (2020). Li et al. (2020) also find that these summary statistics are independent of the primary–secondary halo mass ratio for Ms/Mp < 10−2, with a weak dependence at higher mass ratios. Our model similarly predicts no dependence on mass ratio for Ms/Mp < 10−2. At higher mass ratios our model predicts some dependence, but much weaker than that found by Li et al. (2020).
Summary statistics of the orbital velocity distribution function as predicted by our model (solid lines), and as measured by dashed lines (Li et al. 2020). For the Li et al. (2020) results, we plot lines only over the range of halo masses considered in that work. Statistics shown are the mean (|$\bar{V}$|; left-hand panels), and dispersion (σ; right-hand panels), of the radial (upper panels) and tangential (lower panels) orbital velocities), as a function of the secondary-to-primary halo mass ratio, Ms/Mp. Lines colours indicate different primary halo masses, as indicated in each panel.
Most strikingly, our model predicts that, at primary halo masses below those considered by Li et al. (2020), the mean and dispersion in both radial and tangential velocities will increase significantly. This result can be understood by considering how the velocity dispersion of secondary haloes varies as a function of primary halo mass. Fig. 3 compares halo virial velocities (shown in dimensionful units in this figure) with the linear theory velocity dispersion for 10:1 primary–secondary mass ratio haloes separated by the Lagrangian virial radius of the primary halo (dashed lines) as a function of halo mass.3 At the highest masses shown the velocity dispersion of the secondary haloes is small compared to the primary halo virial velocity. However, as halo mass decreases the virial velocity decreases faster than the velocity dispersion of the secondary haloes. Consequently, for low mass primary haloes the secondary halo velocity dispersion becomes comparable to, and can exceed, the primary halo virial velocity. This increased secondary halo velocity dispersion (when expressed relative to the primary halo virial velocity) is the cause of the increase in the mean and dispersion of Vr and Vθ for low mass primary haloes. We find that the ratio σ/Vv is the primary driver of the shape of the orbital velocity distribution function.
![Virial velocities (solid lines) and halo-halo pairwise velocity dispersions for 10:1 mass ratio haloes separated by the Lagrangian virial radius of the larger halo (dashed lines) as a function of halo mass [equation (5) with the environmental factor of equation (12) applied]. Line colours correspond to different redshifts as indicated in the figure. In this figure, velocities are shown in dimensionful units.](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/mnras/505/2/10.1093_mnras_stab1413/1/m_stab1413fig3.jpeg?Expires=1750100690&Signature=o0sb98f9RQ9962Jqkv7umzzoE0AYVZAAxENDjHdBrTtEs7PjPlpOsO7ksA2W41v6WgBZScwDLiYZ2VDIN~sDgg6sDvpXbYgfYEzrtxHOhWWPRT26bSx6CXWsNaHH6FhqCzZLKvO8xZPalFsATJSlV7XAz26flNsrTExTcK-7ka~-ZYFVd5d8jp9f9aYpxcQrsDdHvx1fGBtPdlICLl2U70qMg-6wzrmLDc02icIt2X1W6clL5w10YldDuiiP5X5yzDIWaT9U7nzCTvw1kQpkruytTCiJrnu-3RjdCYZ7x0izmS-QvAumUa9hLa0zGOrCiemso3Gj2fMvnuGsKrMxEg__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
Virial velocities (solid lines) and halo-halo pairwise velocity dispersions for 10:1 mass ratio haloes separated by the Lagrangian virial radius of the larger halo (dashed lines) as a function of halo mass [equation (5) with the environmental factor of equation (12) applied]. Line colours correspond to different redshifts as indicated in the figure. In this figure, velocities are shown in dimensionful units.
3.1 Redshift dependence
In Fig. 3, we also show results for different redshifts. At z > 0 virial velocities are larger for fixed Mp, while the velocity dispersion of the secondary haloes is lower. Li et al. (2020) considered redshift dependence in the orbital velocity distribution, and showed that, when considered at fixed peak height, |$\nu = \delta _\mathrm{c}/\sqrt{S}$|, where δc is the critical linear theory overdensity for halo collapse, there was little redshift dependence. We therefore recast Fig. 3 in terms of peak height, and plot the velocity dispersion of the secondary haloes now in dimensionless units (i.e. in units of the primary halo virial velocity). The results are shown in Fig. 4.
Halo-halo pairwise velocity dispersions, in units of the primary halo virial velocity, for 10:1 mass ratio haloes separated by the Lagrangian virial radius of the larger halo (dashed lines) as a function of peak height, νp. Line colours correspond to different redshifts as indicated in the figure.
3.2 Truncated power spectra
As the orbital velocity distribution in our model is controlled primarily by the ratio σ/Vv, it will have some dependence on the matter power spectrum (which determines σ). In Fig. 5, we examine how σ is changed if we adopt a truncated power spectrum consistent with a 3 keV thermal warm dark matter particle. For such a particle, the half-mode mass of the power spectrum (Mhm; i.e. the mass corresponding to the wavenumber at which the transfer function is suppressed by a factor of 2 compared to the CDM case) is Mhm = 4.1 × 108M⊙ (Schneider et al. 2012). Comparing with Fig. 3 it is clear that, while there is some reduction in σ at low masses, it is quite small. In particular, at z = 0, we find that σ is reduced by around 30 per cent compared to the CDM case, even though this mass scale lies well below the half-mode mass. The reason for this is that, since the velocity field is sourced by the gradient of the density field, it gives a stronger weighting to large-scale modes – in equation (6) the power spectrum is weighted by k2 when evaluating the variance of the density field, but by k0 when evaluating the dispersion of the velocity field.
Virial velocities (solid lines) and halo–halo pairwise velocity dispersions for 10:1 mass ratio haloes separated by the virial radius of the larger halo (dashed lines) as a function of halo mass in a WDM model. Line colours correspond to different redshifts as indicated in the figure.
Fig. 6 shows the mean and dispersion of the radial velocity distribution as a function of primary halo mass and the secondary-to-primary halo mass ratio for this WDM power spectrum (solid lines), and for the equivalent CDM power spectrum (dashed lines). There is almost no change in these summary statistics due to the truncation of the power spectrum in the WDM case, with the exception of a small reduction in the mean radial velocity for Mp = 108M⊙. Since a truncation of the power spectrum on these scales is already strongly ruled out – for example, Nadler et al. (2021) find that mWDM > 9.7 keV at 95 per cent confidence (corresponding to Mhm < 107.4M⊙) – we conclude that the effects of plausible changes in the power spectrum are unlikely to have a strong effect on the distribution of orbital velocities on mass scales of current astrophysical interest.
Summary statistics of the orbital velocity distribution function as predicted by our model for a thermal WDM power spectrum (solid lines), and for the equivalent CDM power spectrum (dashed lines). Statistics shown are the mean (|$\bar{V}$|; left-hand panel), and dispersion (σ; right panel), of the radial orbital velocities, as a function of the secondary-to-primary halo mass ratio, Ms/Mp. Lines colours indicate different primary halo masses, as indicated in each panel.
4 DISCUSSION
We have described a simple model for the joint distribution function of the orbital velocities, (Vr, Vθ), of merging dark matter halo pairs. Our model considers a dark matter halo perturber, described by a Keplerian potential, embedded within a background of secondary haloes with velocities as predicted by linear perturbation theory, and assuming a Gaussian distribution of pairwise velocities for the background haloes.5 By considering secondary haloes on orbits in the ‘loss cone’ of the primary halo we derive the distribution of merging velocities. Despite its simplicity, and the fact that halo merging occurs far into the non-linear regime of structure formation, our model predictions are in excellent qualitative, and reasonable quantitative agreement with distributions measured from N-body simulations (Li et al. 2020). Our model, coupled with a description of the density profile of the inflowing stream of matter on to a primary halo, also produces good agreement with measurements of halo merger rates (see Appendix B).
Assuming a potential corresponding to an extended profile for the primary halo (rather than a Keplerian potential) would lead to the velocity distribution to be shifted to higher velocity. For example, if we consider an Einasto profile (with concentration parameter, c = 10, and shape parameter, α = 0.18, for specificity) then the potential difference from infinity to the virial radius is approximately |$1.46 V_\mathrm{V}^2$| (as opposed to just |$V_\mathrm{V}^2$| for the Keplerian potential) which would lead to velocities being around 20 per cent larger. However, the primary halo potential is growing via the accretion of secondary haloes (and, perhaps, some smooth accretion of mass). In a simple model in which there is no crossing of mass shells during the growth of the primary (e.g. Fillmore & Goldreich 1984) then the mass contained within the current orbital position of an infalling secondary would be constant (as the secondary comoves with the overall mass shell), and the Keplerian potential approximation would be valid. This assumption will break down for at least two reasons: (1) once the secondary reaches the backsplash radius of the primary there is shell crossing; and (2) since we are considering a distribution of secondary orbits, not all of them can comove with the mass shell.
We expect then that the Keplerian approximation is reasonable for the mean infall velocity (with some inaccuracy arising from shell-crossing close to the primary halo), but will be less accurate for haloes in the tails of the distribution (which experience significant shell-crossing during their infall). An improvement of our model could consider the assembly of the primary halo (perhaps using a secondary infall model; Fillmore & Goldreich 1984; Bertschinger 1985), and track the evolution of each secondary halo orbit in this evolving potential.
Our model connects the distribution of orbital velocities to the underlying matter power spectrum. This allows us to explore predictions for how the distribution of orbital parameters is expected to change for smaller primary halo masses, where we predict a significant increase in the width of these distributions, extending to higher velocities (when expressed in units of the characteristic virial velocity of the primary halo). Furthermore, this simple model shows that the quantity σ/Vv– the ratio of the background velocity dispersion to the primary halo virial velocity – is the primary driver of the shape and width of the orbital velocity distribution. We show that, when considered at constant peak height, ν, this quantity has a very weak dependence on redshift, in agreement with measurements from N-body simulations (Li et al. 2020).
The increased width of the merging halo velocity distributions for lower mass primary haloes may have some important consequences. As these distributions set the initial conditions for further evolution of the subhalo within the potential of its host, they may affect the radial distribution of subhaloes for example. However, primary halo mass scales on which measurements of the subhalo distribution are currently being made (of order 1012M⊙ for studies of dwarf galaxies around the Milky Way, and of order 1013M⊙ for studies of gravitational lensing around massive elliptical galaxies) are in the regime where σ/Vv is small, and, in any case, results from N-body simulations have been explored on these mass scales.
Other possible consequences of the increased width of the merging halo velocity distributions for lower mass primary haloes arise from the connection between halo mergers and halo structure. For example, both halo spin (Vitvitska et al. 2002; Benson, Behrens & Lu 2020) and concentration (Johnson, Benson & Grin 2021) appear to be connected to the orbital parameters of halo merger events. However, in cases where the merging halo velocity distribution width is increased many of the mergers will be unbound, and so may correspond to ‘fly-by’ encounters which have little lasting impact on the primary halo. We intend to explore the consequences of our model on halo structure in a future work.
Further exploiting this direct connection to the power spectrum, we explore how the distribution of orbital velocities is changed in the case of a truncated power spectrum associated with the thermal warm dark matter model. We find that, given current constraints on the mass of a thermal warm dark matter particle, there are no significant changes in the distribution of orbital parameters on scales of current astrophysical interest. As such, it is reasonable to utilize fitting functions for the distribution of orbital velocities of merging dark matter haloes that were derived from simulations of CDM even in the case of plausible WDM models.
In conclusion, our model provides physical insight into the origins of the distribution of orbital velocities of merging halo pairs, and a clear connection to the underlying matter power spectrum. While prior measurements of the orbital parameter distribution from N-body simulations have shown only a very weak dependence on halo mass, our model predicts a stronger dependence on mass scales smaller than those which have been explored in these prior studies. Our model therefore provides guidance on the extent of the parameter space over which fitting functions, often provided by those prior studies, can be considered to be valid. The code used to generate the distribution functions used in this work is publicly available as part of the Galacticus toolkit, which also allows those distribution functions to be utilized as part of other calculations made with the toolkit (e.g. evolving populations of subhaloes).
ACKNOWLEDGEMENTS
We thank Shaun Cole, Benedikt Diemer, Xiaolong Du, Fangzhou Jiang, and Annika Peter for helpful and insightful discussions. Computing resources used in this work were made available by a generous grant from the Ahmanson Foundation.
DATA AVAILABILITY
The data underlying this article will be shared on reasonable request to the corresponding author. The code used to generate the distribution functions used in this work is publicly available as part of the https://github.com/galacticusorg/galacticusGalacticus toolkit.
Footnotes
We could choose to integrate over any of the three variables, r, |$V^\prime _\mathrm{r}$|, or |$V^\prime _\theta$| in this equation, providing the other two variables were appropriately related to the variable integrated over, and the appropriate Jacobian used in place of the one used here. While integrating over r would lead to a more compact expression for the Jacobian, we find that integrating over |$V^\prime _\mathrm{r}$| is computationally more robust and efficient.
We discuss the validity of this assumption in Section 4.
We do not show comparisons to results from N-body simulations in this figure. The velocity dispersions shown are the linear theory expectation, computed at the Lagrangian virial radius of the primary halo (as this is the key quantity which goes into our model). What could be measured from N-body simulations is the velocity dispersion at the corresponding Eulerian separation but, by construction, that is far into the non-linear regime, making a comparison to what is plotted here not meaningful. The model presented in this work takes the linear theory expectation for the velocity distribution, and propagates it into the deeply nonlinear regime. Therefore, the confrontation with N-body results in Fig. 1 is the relevant comparison.
This behaviour can be understood as follows. At fixed density, velocities in any self-gravitating system will scale as V ∝ r. In cosmological linear theory the characteristic density is simply the mean density of the universe, independent of scale. For virialized dark matter haloes the characteristic density is the virial density, also independent of scale. As such, any scale dependence is expected to cancel out in the ratio σ/Vv. Additionally, |$\sigma _\mathrm{V}^2 \propto S$| since both are derived from an integral over the power spectrum. Since |$S \propto \nu _\mathrm{p}^{-2}$| this implied |$\sigma _\mathrm{V} \propto \nu _\mathrm{p}^{-1}$|.
This assumption of Gaussianity could be relaxed, for example using the distribution functions found by Cuesta-Lazaro et al. (2020), although the improved precision in the characterization of the background halo velocity distribution may not be warranted given the other approximations made in our model.
Specifically, at the separations of concern here the inflow stream is spatially coincident with the backsplash halo population of the primary halo. Any determination of the bias of haloes in the inflow stream would need to separate out these two populations.
REFERENCES
APPENDIX A: ORBITAL VELOCITY DISTRIBUTIONS
Fig. A1 shows distributions of orbital velocities as in Fig. 1 but for a much larger primary halo mass of Mp = 1014M⊙. The secondary halo mass is unchanged from Fig. 1 at Ms = 1011M⊙.
Distributions of orbital parameters for merging haloes of masses Mp = 1014M⊙ and Ms = 1011M⊙ at z = 0. Upper left-hand panel: The joint distribution of radial and tangential velocities of merging haloes. The colour map shows results from our model, while the contours indicate the distribution of Li et al. (2020). Upper right-hand panel: The distribution of radial velocities of merging haloes. The blue curve shows results from our model, while the green line is the distribution of Li et al. (2020). Lower left-hand panel: As the upper right-hand panel, but for the tangential velocity. Lower right-hand panel: As the upper right-hand panel, but for the radial velocity expressed as a fraction of the total velocity.
As expected from the summary statistics shown in Fig. 2 our model performs less well here compared to the results shown in Fig. 1. Specifically, while the location of the peak in the radial velocity is well-matched to the results of Li et al. (2020) the dispersion around this mean is too small in our model. This also leads to the distribution of tangential velocities being shifted toward low velocities compared to Li et al. (2020). In the top left-hand panel, the joint distribution of (vr, vθ) retains the ridge-line of approximately constant orbital energy, but the distribution is too narrow around this locus.
Overall, this suggests that our model underestimates the ratio σ/VV in this regime. Possible causes of this underestimate include the limitations of linear perturbation theory, and the approximate nature of our model for the environmental dependence of the secondary halo velocity dispersion.
APPENDIX B: HALO MERGER RATES
Using our model, which predicts |$\bar{v}_\mathrm{r}$|, we can predict the expected merger rates of haloes. For n(r, Mp, Ms) we assume that the density profile of the inflow stream around the primary halo follows the form found by Diemer & Kravtsov (2014), and assume no bias between haloes and mass in this stream such that |$n(r,M_\mathrm{p},M_\mathrm{s}) = n(M_\mathrm{s}) \rho _\mathrm{i}(r)/\bar{\rho }$| where n(Ms) is the usual halo mass function, ρi(r) is the density profile of the inflow stream from Diemer & Kravtsov (2014), and |$\bar{\rho }$| is the mean density of the universe. The assumption of no bias in the inflow is unlikely to be true in detail but has not been extensively studied.6
The results of this simple model are shown, as a function of secondary halo mass and for several different primary halo masses, in Fig. B1, where they are compared to the measured merger rates from Fakhouri et al. (2010). At secondary-primary mass ratios less than around 0.1 the agreement between the results of this simple model and the Fakhouri et al. (2010) measurements are very good – particularly for 1012M⊙ < Mp < 1014M⊙, the regime where our model predictions for |$\bar{v}_\mathrm{r}$| are in good agreement with N-body measurements. At lower values of Mp, our model predicts merger rates slightly larger than those of Fakhouri et al. (2010).
The merger rate, at z = 0, of secondary haloes of mass Ms is shown for several different values of primary halo mass, Mp (as indicated by line colour). Dashed lines show merger rates measured by Fakhouri, Ma & Boylan-Kolchin (2010), while solid lines show estimates based on the results of this work.
At higher mass ratios our model fails to capture the upturn in merger rates measured by Fakhouri et al. (2010). While a detailed analysis of this discrepancy is beyond the scope of this work there are several plausible contributing factors. First, as can be seen in Fig. 2, our model underestimates |$\bar{v}_\mathrm{r}$| for large Ms/Mp, which would lead to an underestimate of the merger rate. Perhaps more importantly the perturbative nature of our model (which assumes that the primary halo is a perturber in the background of secondary haloes) must break down when the secondary halo mass becomes comparable to the primary halo mass. Finally, as mentioned above, our assumption of no bias for the density profile of secondary haloes in the inflow stream is likely to be incorrect and may contribute to the discrepancy at large Ms/Mp.
Using these estimates of halo merger rates the model described in this work could be used to compute the mass growth rates of primary haloes. To do so would require assessing which merging secondaries are actually bound to the primary. Benson (2005) showed that around 18 per cent of merging secondaries are on unbound orbits, but that the majority of these become bound (either due to the effects of dynamical friction or because the primary halo continues to grow after the secondary becomes a subhalo within it). Therefore, to a reasonable approximation all merging secondary haloes contribute to the long-term mass growth of the primary. However, the results of Benson (2005) apply to haloes of masses of around 1011M⊙ and greater, and so may not be valid for lower mass haloes where this model described in this work predicts a much larger fraction of unbound orbits. Computing primary halo growth rates using this approach will therefore require a careful treatment of the orbital evolution of each subhalo and the growth of the primary halo potential.
