Tidal evolution of dark matter annihilation rates in subhalos

TL;DR

This paper develops a high-resolution simulation-based model to predict how tidal forces affect dark matter annihilation rates in subhalos, improving the understanding of signals from dark matter substructure.

Contribution

It introduces a new model that predicts the evolution of subhalo annihilation rates based on physically motivated parameters, enhancing accuracy over previous estimates.

Findings

Model accurately predicts the impact of tidal forces on annihilation rates.

Parameters include energy injection, tidal force ratio, and tidal heating distribution.

Model has broader implications for understanding subhalo evolution beyond annihilation signals.

Abstract

Dark matter halos grow by hierarchical clustering as they merge together to produce ever larger structures. During these merger processes, the smaller halo can potentially survive as a subhalo of the larger halo, so a galaxy-scale halo today likely possesses a rich abundance of substructure. This substructure can greatly boost the rate of dark matter annihilation within the host halo, but the precise magnitude of this boost is clouded by uncertainty about the survival prospects of these subhalos. In particular, tidal forces gradually strip material from the subhalos, reducing their annihilation signals and potentially destroying them. In this work, we use high-resolution idealized $N$-body simulations to develop and tune a model that can predict the impact of this tidal evolution on the annihilation rates within subhalos. This model predicts the time evolution of a subhalo's annihilation rate as a function of three physically motivated parameters of the host-subhalo system: the energy injected into subhalo particles per orbit about the host, the ratio of stretching to compressive tidal forces, and the radial distribution of tidal heating within the subhalo. Our model will improve the accuracy of predictions of the magnitude and morphology of annihilation signals from dark matter substructure. Additionally, our parametrization can describe the time evolution of other subhalo properties, so it has implications for understanding aspects of subhalo tidal evolution beyond the annihilation rate.

Figures

Equations

$$\rho (r)=\frac{{\rho }_s}{(r/r_s)(1+r/r_s{)}^2},$$

$$J\equiv \int {\rho }^2\mathrm{d}V,$$

$${\mathbf{F}}_{\mathrm{tidal}}(\mathbf{r})=-\frac{\mathrm{d}F}{\mathrm{d}R}(\mathbf{r}\cdot \hat{\mathbf{R}})\hat{\mathbf{R}}-F(R)\frac{\mathbf{r}-(\mathbf{r}\cdot \hat{\mathbf{R}})\hat{\mathbf{R}}}{R}$$

$$r_s\lesssim 0.1R_c.$$

$$\frac{M/m}{\ln (M/m)}\gtrsim 10(1+z{)}^{3/2}.$$

$$\tilde{x}\equiv \frac{R_c{\rho }_s}{R_s{P}_s}$$

$$\frac{\mathrm{d}\ln J}{\mathrm{d}\ln t}=-b{n}^{1-c}$$

$$\frac{1}{J}\frac{\mathrm{d}J}{\mathrm{d}n}=-b{n}^{-c}.$$

$$\frac{\mathrm{d}\ln J}{\mathrm{d}\ln t}\simeq -\min \left\{b{n}^{1-c},B\right\},$$

$$r_{n+1}^{\gamma }-r_n^{\gamma }+f{r}_{n+1}^{\alpha }=0,$$

$$\frac{\Delta \ln r}{\Delta \ln n}\simeq \{\begin{array}{ll}-fn/\gamma ,& fn\ll 1,\\ -1/(\alpha -\gamma ),& fn\gg 1,\end{array}$$

$$\frac{\Delta \ln J}{\Delta \ln n}\simeq -\min \left\{bn,B\right\}$$

$$\frac{1}{J}\frac{\mathrm{d}J}{\mathrm{d}n}=-bq,\frac{1}{q}\frac{\mathrm{d}q}{\mathrm{d}n}=-c{q}^{1/c}.$$

$$\ln \frac{J}{J_{\mathrm{init}}}=b\left[a-\frac{1}{1-c}\left(n^{1-c}-1\right)\right]$$

$$E_b=-4\pi (\ln 2)G{\rho }_s{r}_s^2$$

$$x\equiv |E_b|/\Delta E_{\mathrm{imp}}$$

$$x\sim {\rho }_s/\bar{P}(<R_c)$$

$$b=b_0{x}^{-b_1},\text{if}{R}_c\ll R_s,$$

$$b=b_0{x}^{-b_1}\left[1+b_{2}f(y)\right],$$

$$y\equiv \bar{R}/R_s.$$

$$f(y)=\frac{2\ln (1+y)-y(2+3y)/(1+y{)}^2}{\ln (1+y)-y/(1+y)}.$$

$$a=a_0-a_{1}f(y)$$

$${\rho }_a\equiv \frac{V_p^2}{G{R}_p^2}={\eta }^2\frac{M(R_c)R_c}{R_p^4},$$

$${\rho }_t\equiv M(R_a)/R_a^3,$$

$$z\equiv r_a/r_t,$$

$$c=c_0{z}^{c_1}$$

$$\ln \frac{J}{J_{\mathrm{init}}}=b\left[a-\frac{1}{1-c}\left(n^{1-c}-1\right)\right].$$

$$a$$

$$b$$

$$c$$

$$n_1=b^{1/(c-1)}$$

Full text

I Introduction

Despite overwhelming evidence for the existence of dark matter (e.g., Refs. Zwicky (1933); Clowe et al. (2004); Ade et al. (2016)), its microphysical details remain unknown. Numerous models have been considered, but none have been experimentally confirmed (see Refs. Bertone and Tait (2018); Feng (2010); Bertone (2010); Bertone et al. (2005); Bergström (2000) for reviews). However, to explain the present abundance of dark matter, a large class of models, including the popular weakly interacting massive particle Jungman et al. (1996), propose that dark matter was pair produced from the thermal plasma in the hot early universe. In this scenario, the dark matter can annihilate back into standard-model particles today, leading to prospects for the detection of high-energy gamma rays or other annihilation products (e.g., Ref. Strigari et al. (2007)).

The rate of dark matter annihilation scales as the square of the dark matter density, so it is strongly sensitive to the spatial distribution of the dark matter. At galactic scales and above, this spatial distribution is well understood. Initially overdense patches in the Universe collapse into gravitationally bound dark matter halos, which thereafter merge to produce successively larger structures. Numerical simulations demonstrate that the spherically averaged mass distributions of the resulting dark matter halos are well described by the NFW density profile Navarro et al. (1996, 1997),

[TABLE]

which has scale parameters $r_{s}$ and $\rho_{s}$. Baryonic effects may subsequently alter this density profile (e.g., Ref. Brooks (2014)).

However, at subgalactic scales, the spatial distribution of the dark matter is less clear. As a halo is built up through hierarchical merging, it accretes smaller halos that survive as subhalos within the larger host. These subhalos gradually lose mass due to the influence of tidal forces from the host (e.g., Refs. Binney and Tremaine (1987); Mo et al. (2010)), surviving until either they are completely stripped or dynamical friction causes them to sink into the host’s center Chandrasekhar (1943). However, for sufficiently small subhalos dynamical friction is inefficient Mo et al. (2010). Moreover, numerous analyses have found that if the subhalos possess divergent central density, as in the NFW profile, then tidal forces may never fully strip them Goerdt et al. (2007); Penarrubia et al. (2010); Errani and Peñarrubia ; Berezinsky et al. (2008).

Thus, a galactic-scale dark matter halo is likely to possess a multitude of subhalos. This substructure can significantly boost the rate of dark matter annihilation, depending on the scale of the smallest halos and when they form Berezinsky et al. (2003); Diemand et al. (2005); Pieri et al. (2008); Berezinsky et al. (2008); Ishiyama et al. (2010); Anderhalden and Diemand (2013a); *anderhalden2013erratum; Ishiyama (2014); Sánchez-Conde and Prada (2014); Anderson et al. (2016); Gao et al. (2012); Springel et al. (2008); Stref and Lavalle (2017); Stref et al. (2019); Bartels and Ando (2015); Hiroshima et al. (2018) (see Ref. Ando et al. (2019) for a recent review). When the smallest halos are microhalos of roughly earth mass, annihilation rates may be boosted by a factor of about $10$ Stref and Lavalle (2017); Hiroshima et al. (2018) relative to those expected in the absence of substructure, assuming that these halos arise from primordial density fluctuations comparable to the large-scale fluctuations inferred from the cosmic microwave background (e.g., Ref. Ade et al. (2016)). Moreover, this boost can be raised by orders of magnitude by cosmological scenarios that amplify small-scale density fluctuations and thereby lead to earlier and more abundant microhalo formation. Such scenarios include a period of domination by a heavy species Erickcek and Sigurdson (2011); Barenboim and Rasero (2014); Fan et al. (2014); Erickcek (2015) or a fast-rolling scalar field Redmond et al. (2018) prior to nucleosynthesis, along with a variety of inflationary models Silk and Turner (1987); Salopek et al. (1989); Starobinskij (1992); Ivanov et al. (1994); Randall et al. (1996); Stewart (1997); Adams et al. (1997); Starobinsky (1998); Covi and Lyth (1999); Martin et al. (2000); Chung et al. (2000); Martin and Brandenberger (2001); Joy et al. (2008); Barnaby and Huang (2009); Barnaby (2010); Ben-Dayan and Brustein (2010); Gong and Sasaki (2011); Lyth (2011); Bugaev and Klimai (2011); Barnaby and Peloso (2011); Achúcarro et al. (2011); Cespedes et al. (2012); Barnaby et al. (2012). In these cases, the high density within these microhalos causes them to completely dominate any dark matter annihilation signal (e.g., Refs. Bringmann et al. (2012); Delos et al. (2018a); Blanco et al. ).

Unfortunately, all estimates of the substructure’s boost to annihilation rates are subject to uncertainties about the impact of tidal effects on subhalos. Cosmological simulations cannot resolve subhalos that are much smaller than the host, and those that are resolved are prone to artificial destruction van den Bosch (2017); van den Bosch et al. (2018); van den Bosch and Ogiya (2018). Numerous semianalytic models have been developed to describe the dynamical evolution of subhalos; Refs. Taylor and Babul (2001); Peñarrubia and Benson (2005); van den Bosch et al. (2005); Zentner et al. (2005); Kampakoglou and Benson (2007); Gan et al. (2010); Pullen et al. (2014); Jiang and van den Bosch (2016) model a subhalo’s loss of mass due to tidal stripping, and Refs. Hayashi et al. (2003); Penarrubia et al. (2010) predict the impact of this mass loss on the halo’s density profile. However, these models are typically tuned to the results of cosmological simulations, so their predictions are affected by the artificial subhalo disruption occurring therein. They cannot fully reproduce the results of idealized simulations van den Bosch et al. (2018); Ogiya et al. (2019).

Meanwhile, calculations of the dark matter annihilation rate in the substructure have employed a number of different treatments of tidal evolution. Some, such as Refs. Bartels and Ando (2015); Hiroshima et al. (2018), employ a combination of the models above to predict the time evolution of subhalo density profiles. Others, such as Refs. Springel et al. (2008); Erickcek (2015); Stref and Lavalle (2017); Stref et al. (2019), employ simpler models, often either truncating subhalos at a characteristic tidal radius or formulating a destruction condition for subhalos and assuming the survivors are unaltered. Still others, such as Refs. Bringmann et al. (2012); Anderhalden and Diemand (2013a); *anderhalden2013erratum; Ishiyama (2014); Sánchez-Conde and Prada (2014); Anderson et al. (2016); Delos et al. (2018a); Blanco et al. neglect the tidal disruption of the substructure altogether.

Our work is motivated by this context. Since cosmological simulations cannot resolve the smallest substructures, we follow Refs. Hayashi et al. (2003); Kazantzidis et al. (2004); Read et al. (2006); Penarrubia et al. (2010); van den Bosch and Ogiya (2018); Ogiya et al. (2019) in using idealized simulations of an $N$-body subhalo inside an analytic galactic potential. However, unlike these works, we focus on understanding the impact of tides on the subhalo’s annihilation rate, a goal that requires significantly better resolution than has been attained in previous studies. Moreover, previous works have focused on understanding the evolution of subhalos of scales resolvable in cosmological simulations, such as halos associated with dwarf galaxies within a galactic halo. Accordingly, they probe only the subhalo properties and orbits that are found in such simulations. For instance, Ref. Ogiya et al. (2019) only studies subhalos orbiting above the host’s scale radius. In contrast, we seek to probe the full range of subhalos down to the smallest microhalos, which span a far broader range of properties and orbits.

Using the results of 52 high-resolution $N$-body simulations, we develop a physically motivated model that can predict the time evolution of a subhalo’s annihilation rate due to tidal effects. In the process, we isolate three physical variables that determine this evolution:

1. (1)

The energy injected by tidal forces into subhalo particles over the course of each orbit about the host, in units of the particle’s binding energy to the subhalo; 2. (2)

The ratio of stretching (radial) tidal forces to compressive (tangential) tidal forces; 3. (3)

The range of radii in the subhalo across which material is heated by tidal forces, which is set by the shape of the subhalo’s orbit.

This model predicts the suppression of a subhalo’s annihilation rate, characterized by its $J$ factor111We assume the dark matter annihilation cross section is velocity independent in the nonrelativistic limit.

[TABLE]

as a function of its orbit about the host. To assist the application of our model, we supply convenient fitting functions.

This article is organized as follows. In Sec. II, we detail how we carry out our $N$-body simulations. Section III qualitatively discusses the trajectory of a subhalo’s $J$ factor, interpreting simulation trends physically and motivating our model. In Sec. IV, we develop our predictive model for the evolution of a subhalo’s $J$ factor, and Sec. V summarizes the model and discusses limitations and extensions. In Sec. VI, we compare the model’s predictions to those of previous semianalytic models. Section VII concludes, after which we supply a variety of appendixes. Appendix A supplies further details about our simulations, while Appendix B quantifies the range of subhalo sizes over which the results of these simulations are applicable. Appendix C presents fitting formulas and other computational details that aid in applying our model. Appendix D tests our model against a publicly available simulation library Ogiya et al. (2019). Finally, in Appendix E, we observe that our model can be adapted to describe the evolution of subhalo properties beyond the $J$ factor.

II Simulations

Owing to the difference in scales between a host and its smallest subhalos, the computational challenge in simulating subhalo evolution in a cosmological context is formidable. A number of previous works have addressed this problem by simulating an $N$-body subhalo inside an analytic host potential Hayashi et al. (2003); Kazantzidis et al. (2004); Read et al. (2006); Penarrubia et al. (2010); van den Bosch and Ogiya (2018); Ogiya et al. (2019); our approach is similar but differs in one key step. Instead of placing a subhalo in orbit about the host potential, we subject the subhalo directly to the time-dependent tidal force field experienced by an analytic orbit about the host. This procedure minimizes the impact of numerical precision errors that can result from differences in scale between the subhalo’s orbital and internal dynamics. In this section, we detail that procedure and present qualitative results.

We assume that both the host and the subhalo possess the NFW density profile given by Eq. (1). While there is evidence that many galactic halos possess constant-density cores instead of the NFW profile’s cusp Moore (1994), at least some galactic halos appear to be cuspy Read et al. (2018). Additionally, while microhalos are expected to form with $\rho\propto r^{-3/2}$ inner profiles Ishiyama et al. (2010); Anderhalden and Diemand (2013a); *anderhalden2013erratum; Ishiyama (2014); Polisensky and Ricotti (2015); Ogiya and Hahn (2018); Delos et al. (2018b, a); Angulo et al. (2017); Delos et al. (2019), it is likely that mergers will drive their inner cusps toward the $\rho\propto r^{-1}$ of the NFW profile Ogiya et al. (2016); Angulo et al. (2017); Delos et al. (2019).

To model the host’s tidal field we begin with an analytically computed orbit, described by the time-dependent vector $\bm{R}(t)$ pointing from host center to subhalo center. The tidal acceleration at position $\bm{r}$ relative to the subhalo center is222We experimented with using the full tidal force $\bm{F}_{\mathrm{tidal}}(\bm{r})=\bm{F}(\bm{R}+\bm{r})-\bm{F}(\bm{R})$, but because $r\ll R$ in our simulations, it offers no advantage; moreover, it is less numerically stable due to the subtraction of two close numbers.

[TABLE]

at linear order in $r/R$, where $F(R)$ is the force profile of the host, $R=|\bm{R}|$, and $\bm{\hat{R}}=\bm{R}/R$. We modified the Gadget-2 $N$-body simulation code Springel et al. (2001); Springel (2005) to include this tidal acceleration.

We prepare the initial $N$-body subhalo with an NFW profile by drawing particles from an isotropic distribution function computed using the fitting form in Ref. Widrow (2000). Additionally, we sample the subhalo’s central region at increased resolution; particles whose orbital pericenters are below $r_{s}/3$, where $r_{s}$ is the subhalo scale radius, have $1/64$ the mass and $64$ times the number density of the other particles. Appendix A demonstrates that there is no significant relaxation associated with the use of particles of different masses. We cut off the density profile at $r=500r_{s}$; subhalo particles this far out are stripped immediately, so as long as the cutoff radius is much larger than $r_{s}$, the precise choice makes no difference.333The natural place to cut off the density profile would be where the density reaches that of the subhalo’s background: the host. However, tidal forces automatically truncate a subhalo’s density profile at roughly the radius where its average density equals that of the host [see, e.g., Eq. (24)], so it is not necessary to tune a cutoff radius by hand. We represent the subhalo using a total of $8\times 10^{6}$ particles, and roughly 70% of them, carrying roughly 4% of the total mass, are high-resolution particles. All of our subhalos have $r_{s}\simeq 10^{-6}R_{s}$, where $R_{s}$ is the scale radius of the host, but as we will soon discuss, the precise choice of $r_{s}$ has no impact on dynamics.

For our simulations, we consider a variety of orbits about the host. An orbit in a spherically symmetric potential is characterized by two parameters: energy $E$ and angular momentum $L$ or, equivalently, a scale parameter and a shape parameter. For convenience, we use the circular orbit radius444Note that $R_{c}$ is roughly the time-averaged radius; for a power-law potential $\phi(R)\propto R^{n}$, $R_{c}=\langle R^{n}\rangle^{1/n}$. See Appendix C for a more precise relationship for NFW profiles. $R_{c}$, defined as the radius of the circular orbit with energy $E$, and the “circularity” $\eta=L/L_{c}$, where $L_{c}$ is the angular momentum of the circular orbit with the same energy. In each simulation the subhalo begins at its orbital apocenter.

Figure 1 illustrates a simulation executed through this arrangement. The host has scale radius $R_{s}=0.8$ kpc and scale density $P_{s}=$5\text{\times}{10}^{7} $M_{\odot}/\text{kpc}^{3}$, while the $N$-body subhalo is initially a microhalo with scale radius $r_{s}=$6\text{\times}{10}^{-7}\text{\ }\mathrm{k}\mathrm{p}\mathrm{c} that has $\rho_{s}/P_{s}=1285$ times the scale density of the host. The subhalo orbit has $R_{c}=0.15$ kpc and $\eta=0.5$. The simulation runs through 18 orbits about the host, and Fig. 2 plots the density profile of the subhalo at each apocenter. Consistently with the results of other works, such as Refs. Goerdt et al. (2007); Penarrubia et al. (2010); Errani and Peñarrubia ; Berezinsky et al. (2008), we find that this subhalo’s central cusp is highly resistant to disruption by the host’s tidal forces. 91% of the subhalo’s mass is stripped by simulation termination, but its central density profile is largely unscathed.

Our goal in this work is to understand how the annihilation signal decays due to tidal effects. For this purpose, we consider the $J$ factor, Eq. (2), integrated over the subhalo mass distribution. This $J$ factor is the factor in the annihilation rate that depends on mass distribution, and Appendix A discusses our procedure to extract it from the simulations. We also show in Appendix A that the resulting $J$-factor trajectories are converged with respect to simulation parameters.

Finally, we conclude this section by discussing the applicability of our simulation results. The linearized tidal force in Eq. (3) is valid for $r\ll R$, and in Appendix B we show that it yields accurate results in simulations as long as

[TABLE]

Our simulation results are only applicable if this condition is satisfied. Additionally, since we do not simulate the host halo’s dynamics, we cannot account for dynamical friction. Reference Ogiya et al. (2019) found that dynamical friction555Specifically, Ref. Ogiya et al. (2019) studied dynamical self-friction, or the dynamical friction that results from the subhalo’s own tidal tail. This friction can be considerably more efficient than that resulting from the host’s material alone Fujii et al. (2006); Fellhauer and Lin (2007). has minimal impact on a subhalo’s mass evolution for host-to-subhalo mass ratios $M/m\gtrsim 100$. However, since subhalos relevant to dark matter annihilation may orbit over significantly longer timescales than considered in Ref. Ogiya et al. (2019), it is also useful to have an analytic estimate for when dynamical friction can be neglected. It follows from the analysis in Ref. Mo et al. (2010) that for a subhalo that accretes onto a host at redshift $z$, dynamical friction can be neglected as long as the host-to-subhalo mass ratio $M/m$ satisfies

[TABLE]

Our results may be considered applicable as long as Eqs. (4) and (5) are satisfied, but we remark that if one is satisfied, then the other likely is too.

III Trends in the tidal evolution

In this section, we explore trends in the evolution of $J$ as a function of system and time and attempt to explain them physically. Our goal is to find the $J$ factor as a function of time $t$, orbital parameters $R_{c}$ and $\eta$, subhalo parameters $r_{s}$ and $\rho_{s}$, and host parameters $R_{s}$ and $P_{s}$. The dimensionality of this space is large, but some immediate simplifications are evident:

1. (1)

As long as the subhalo is much smaller than its orbit, or $r_{s}\ll R_{c}$, the value of $r_{s}$ has no impact on dynamics.666For fixed density, the internal velocities of particles in a subhalo are proportional to its radius. Since the tidal acceleration in Eq. (3) is also proportional to the radius, fractional velocity changes induced by tidal forces are independent of the subhalo’s radius. All of our simulated subhalos have $r_{s}\simeq 10^{-6}R_{s}$, which leads to $r_{s}\ll R_{c}$ for all orbits we consider. 2. (2)

If instead of time $t$ we use the orbit count $n=t/T$, where $T$ is the orbital period, then the overall density scale has no impact on dynamics, and only the ratio $\rho_{s}/P_{s}$ enters. 3. (3)

The overall size of the host-subhalo system is irrelevant, so only the ratio $R_{c}/R_{s}$ affects dynamics.

We have verified that all of these simplifications are borne out in our simulations. Hence, if $J_{\mathrm{init}}$ is the initial $J$ factor, then $J/J_{\mathrm{init}}$ is now a function of time $n=t/T$ and just three system parameters: $\rho_{s}/P_{s}$, $R_{c}/R_{s}$, and $\eta$. Notably, the tidal evolution is independent of the subhalo’s mass, a property also noted in prior works (e.g., Ref. Ogiya et al. (2019)).

III.1 Trends in the simulations

We first inspect the results of selected simulations in order to find trends in the behavior of $J$. As a further simplification, we focus on the $R_{c}\ll R_{s}$ regime. The host potential is self-similar in this regime, reducing the tidal evolution problem in two additional ways:

1. (1)

The orbital radius $R_{c}$ is degenerate with properties of the host and subhalo. For instance, reducing the orbital radius is equivalent to making the host denser. 2. (2)

Orbits with the same circularity $\eta$ have the same shape; they are rescaled versions of one another.

The first simplification further reduces the parameter space so that in this self-similar regime, there are only two parameters,777We reserve $x$ (without the tilde) for later use as a modified version of $\tilde{x}$.

[TABLE]

and $\eta$. Figure 3 shows the success of this parameter reduction; different systems with the same $\tilde{x}$ and $\eta$ follow precisely the same $J(n)$ trajectories. Meanwhile, the second simplification allows us to isolate the impact of these two parameters; we can vary the “reduced orbital radius” $\tilde{x}$ without altering the orbit’s shape.

We first investigate the impact of orbit shape. Figure 4 shows the trajectory of the $J$ factor for several values of $\eta$. We see immediately that the $J$ factor oscillates with the orbital period with a larger amplitude for more eccentric orbits. This trend is explained by noting that all tidal forces are compressive in the self-similar regime,888Tangential tidal forces are always compressive, while radial tidal forces are negligible when $R_{c}\ll R_{s}$; see Sec. IV.1 for further discussion. so the subhalo becomes most compact near the orbital pericenter. The subhalo’s $J$ factor, being proportional to its mass-weighted average density [e.g., Eq. (32)], is maximized at this point. The precise appearance of these oscillations can be complicated because a subhalo’s response to these tidal forces is delayed; for instance, double peaks in Fig. 3 arise because the subhalo and its unbound tidal stream are maximally compressed at different times. However, these oscillations are relatively unimportant. If subhalos are at random points in their orbits, then the $J$ factor averaged over an orbital period suffices to predict the aggregate signal from a population of subhalos. We discuss this point further in Sec. V.

More interesting trends arise in the broader time evolution. Figure 5 plots the running power-law index $\mathrm{d}\ln J/\mathrm{d}\ln t$ of the $J$ factor with time, and the equation

[TABLE]

describes the evolution of this index reasonably well as long as $|\mathrm{d}\ln J/\mathrm{d}\ln t|<\mathcal{O}(1)$. Here, $b>0$ and $c>0$ are constant parameters, and $c$ is smaller for more eccentric orbits. Figures 4 and 5 also show fits to the $J$-factor trajectories using this form, which determines $J(n)$ up to a constant multiple. The exponent in Eq. (7) is so defined because it leads to the more evocative expression

[TABLE]

If $c=0$, this equation tells us that the $J$ factor would decay by the same factor $\mathrm{e}^{-b}$ over each orbit. The parameter $c$, when $c>0$, accommodates some physical process by which tidal effects lose efficiency over time.

We show the impact of the orbital radius in Fig. 6, which plots the trajectories of $J$ and $\mathrm{d}\ln J/\mathrm{d}\ln t$ for a variety of reduced orbital radii $\tilde{x}$ ranging from $0.7$ to $230$. Evidently, $\tilde{x}$ affects the initial decay rate of the $J$ factor, described by the parameter $b$ in Eq. (7), without altering the rate at which the decay slows over time. This figure also shows more clearly that there is a steepness limit to the decay of the $J$ factor:

[TABLE]

where $B\sim\mathcal{O}(1)$.

III.2 Physical interpretation

Behavior similar to that of Eq. (9) can be reproduced in a toy model. Suppose the subhalo has potential ${\phi(r)\propto r^{\gamma}}$ up to an additive constant; for instance, an NFW profile would have $\gamma=1$ for ${r\ll r_{s}}$. Now discretize time, perhaps as a count of orbits, letting $r_{n}$ be the radius of the subhalo’s outer boundary at time $t_{n}$. Any material outside $r_{n}$ at time $t_{n}$ is free, fixing the additive constant in the potential such that $\phi_{n}(r)\propto r^{\gamma}-r_{n}^{\gamma}$ (where $\phi\geq 0$ implies freedom). Now suppose that at each time step, particles in the subhalo experience an injection of energy ${\Delta E\propto r^{\alpha}}$ due to tidal forces, and any radius with $\Delta E+\phi>0$ no longer belongs to the halo. For instance, if the time step is constant, then $\alpha=2$ since the energy injection is proportional to the square of the tidal force, which is in turn proportional to the radius. This tidal heating rule leads to the evolution equation $\phi_{n}(r_{n+1})+\Delta E(r_{n+1})=0$, or

[TABLE]

where $f$ incorporates all of the proportionality constants.

For simplicity, we may assume $r_{0}=1$, absorbing its dimensionful value into $f$. For $\alpha>\gamma$, Eq. (10) obeys

[TABLE]

where $\Delta$ denotes the discrete difference across time steps. If $J\propto r^{\beta}$, then

[TABLE]

with $b=\beta f/\gamma$ and $B=\beta/(\alpha-\gamma)$. The quantity $fn$ can be understood as (up to factors of order unity) the ratio of the tidal energy injection to the subhalo’s internal energy. The two separate regimes arise physically because when $fn\ll 1$, the total radius change $|r_{n}-r_{0}|\ll r_{0}$. Since the radius does not change appreciably, the efficiency of tidal heating does not change, so $r$ and $J$ drop by the same fraction in each orbit. However, when $fn\gg 1$, $|r_{n}-r_{0}|\sim r_{0}$. In this case, the radius is decreasing significantly, which implies that the density at the halo’s shrinking outer boundary is increasing and hence that the halo is becoming more difficult to strip.

This toy model has reproduced Eq. (7) with ${c=0}$, successfully explaining the apparent upper limit in $|\mathrm{d}\ln J/\mathrm{d}\ln t|$. We remark, however, that there is another, completely different, physical reason to expect an upper limit in $|\mathrm{d}\ln J/\mathrm{d}\ln t|$: an unbound tidal stream grows in length $L$ as $L\propto t$. Hence, its volume grows as ${V\propto t}$, so its $J$ factor drops as $J\propto M^{2}/V\propto t^{-1}$. Once a subhalo has been stripped to the point that its own $J$ factor is dwarfed by that of its tidal stream, the $J$ factor of the subhalo remnant decays as ${J\propto t^{-1}}$. The combination of these two processes—the increasing density of the subhalo as its radius drops and the $J$ factor of its tidal stream—may explain the behavior in Figs. 5 and 6 wherein $|\mathrm{d}\ln J/\mathrm{d}\ln t|$ initially shallows toward some value larger than $1$ before subsequently returning back to $1$. Note, however, that the precise evolution of the $J$ factor in the $|\mathrm{d}\ln J/\mathrm{d}\ln t|\sim 1$ regime is of little consequence. By this point, the subhalo has already lost most of its $J$ factor and contributes little to annihilation signals.

The physical explanation for the $c>0$ behavior observed in the simulations remains unclear. However, it is necessarily connected to how the shape of the subhalo’s density profile changes in response to tidal effects (see Fig. 2), which the toy model does not account for. In a more complete picture, the rate $\mathrm{d}J/\mathrm{d}n$ of tidal evolution should be sensitive only to the instantaneous host-subhalo system with no explicit dependence on the time $n$. Hence, it should be possible to replace the factor $n^{-c}$ in Eq. (8) with a function of the subhalo’s density profile (and other properties of the system). However, in the $|\mathrm{d}\ln J/\mathrm{d}\ln t|<1$ regime, the total change in $J$ is much smaller than $J$ itself, and if we neglect changes in the shape of the density profile, then any parameter of the density profile (e.g., $\rho_{s}$ or $r_{s}$) must experience a similarly small change. Since the factor $n^{-c}$ can change by an order of magnitude in the same regime, it is not possible, except in a very contrived way, to replace this factor with a function of the density profile.

Thus, the $c>0$ behavior must follow from changes in the density profile’s shape. As a result of these changes, the density profile picks up new parameters that can potentially vary wildly without significantly altering $J$, and the factor $n^{-c}$ can be replaced with a function of those parameters. For instance, by introducing a new parameter $q$, we can write

[TABLE]

This system no longer has explicit time dependence, but if $q=1$ when $n=1$, then it is equivalent to Eq. (8).

IV Modeling the tidal evolution

Motivated by the results of the previous section, we seek a model of the form

[TABLE]

for the case where $|\mathrm{d}\ln J/\mathrm{d}\ln t|<1$. The parameters $b$ and $c$ follow immediately from Eq. (7), and we have inserted another parameter $a$ to fix the overall normalization. Our goal is now to relate $a$, $b$, and $c$ to the parameters of the host-subhalo system. For this purpose, we use the results of 52 idealized $N$-body simulations that we carried out as described in Sec. II. The parameter space covered by these simulations is depicted in Fig. 7.

For each simulation we obtain the trajectory of the subhalo’s $J$ factor, stopping if ${|\mathrm{d}\ln J/\mathrm{d}\ln t|\geq 1}$ or otherwise at an arbitrarily chosen simulation termination time. As we discussed in the previous section, the evolutionary behavior changes markedly when ${|\mathrm{d}\ln J/\mathrm{d}\ln t|\geq 1}$, but precise predictions in this regime are unnecessary. Next, we convolve the $J$-factor trajectory in log space with a top-hat filter of width equal to the radial (apocenter-to-apocenter) orbit period.999For circular orbits the radial orbit period is ill defined, and we substitute its limit as the orbit approaches circular as obtained using the fitting form in Appendix C. This step suppresses the influence of the periodic oscillatory behavior observed in Sec. III; we are effectively finding the moving logarithmic average of $J$ over this period. Finally, we fit Eq. (14) to this smoothed trajectory of $J$, but we only employ times after the end of the first radial period (so the first point is at $n=1.5$, whose corresponding $J$ factor averages from $n=1$ to $n=2$). This restriction is intended to remove the influence of any transient effects associated with suddenly turning on the tidal field. Additionally, in case the smoothing procedure fails to fully suppress periodic effects, we minimize any resulting bias by ending the fit at an integer number of orbits (so for instance, we might end at $n=15.5$, corresponding to the average $J$ from $n=15$ to $n=16$). The number of radial orbits fit through this procedure is represented in Fig. 7 as the marker size; this number is a proxy for how much information that simulation provides.101010When performing fits, we weight a simulation spanning $n$ orbits by $\sqrt{n}$. Figure 8 illustrates the procedure; the smoothing filter suppresses oscillations quite effectively.

IV.1 Parameter $b$: The initial $J$-factor decay rate

From Fig. 6 we anticipate that $b$ should depend strongly on the orbital radius. For simplicity we first study the self-similar regime, $R_{c}\ll R_{s}$. The upper panel of Fig. 9 plots $b$ against the reduced orbital radius $\tilde{x}$ for the 36 of our simulations that satisfy $R_{c}/R_{s}<0.3$. While $b$ is strongly sensitive to $\tilde{x}$, there is also significant sensitivity to the orbital shape, parametrized by $\eta$. However, it turns out that we can eliminate the shape dependence of $b$ by defining the reduced orbital radius more carefully.

We first remark that $\tilde{x}\sim|E_{b}|/\Delta E$, where $E_{b}$ is the binding energy of a particle at the subhalo’s initial scale radius $r_{s}$ and $\Delta E$ is the energy injected into that particle by tidal forces over the subhalo’s orbital period. To see this, observe that the particle’s binding energy is

[TABLE]

(per mass). Meanwhile, the tidal acceleration on this particle is roughly $(r_{s}/R)F$, where $R$ is the subhalo’s orbital radius and $F$ is the host force (per mass) at radius $R$. In the self-similar regime, $F\sim GP_{s}R_{s}$. The total velocity injected into the particle is $\Delta v\sim FT$, where $T$ is the subhalo orbital period. Since $T\sim\sqrt{R/F}$, the energy injection (per mass) is $\Delta E\sim(\Delta v)^{2}\sim GP_{s}R_{s}r_{s}^{2}/R$, and since $R\sim R_{c}$, this leads to $|E_{b}|/\Delta E\sim\tilde{x}$.

With this motivation, we define

[TABLE]

as a more exact version of $\tilde{x}$. Here, $\Delta E_{\mathrm{imp}}$ is the energy injection per orbit on a particle at $r_{s}$ computed using the impulse approximation as in Ref. Gnedin et al. (1999). In this approximation, the subhalo particle is treated as stationary while Eq. (3) is integrated to find the velocity (and hence energy) injection. We supply a fitting formula for $\Delta E_{\mathrm{imp}}$ in Appendix C for convenience. As intuition, $x$ is of order the ratio

[TABLE]

between the subhalo’s density and the average host density within the subhalo’s orbital radius, a connection that follows from the observation that ${\Delta E_{\mathrm{imp}}/r_{s}^{2}\sim F(R_{c})/R_{c}}$ (see Appendix C). In the bottom panel of Fig. 9 we plot $b$ against $x$ for the self-similar regime. Evidently, our definition of $x$ captures most or all of the sensitivity of the parameter $b$ to the orbit shape, and

[TABLE]

with $b_{0}=0.58$ and $b_{1}=0.58$. This success is remarkable; the approximation that subhalo particles are stationary during the application of tidal forces can only be valid for highly eccentric orbits, and yet the impulsive energy calculation accurately predicts the tidal evolution for more circular orbits as well.

To complete our understanding of the parameter $b$ we must move beyond the self-similar regime. In the upper panel of Fig. 10, we plot $b$ against $x$ for all of our simulated subhalos. The color scale indicates the time-averaged orbital radius $\bar{R}$ in units of $R_{s}$. It appears that the effect of leaving the self-similar regime is to alter the normalization of $b$ while keeping the power-law sensitivity to $x$ unchanged. In particular, we may write

[TABLE]

for some function $f(y)$ and parameter $b_{2}$, where we define

[TABLE]

For convenience, we supply a fitting formula for $\bar{R}$ in Appendix C. While we could use the circular orbit radius $R_{c}$ instead, we favor $\bar{R}$ because its physical significance is clearer.

To define the function $f(y)$, we consider the physical impact of leaving the self-similar regime. The magnitudes of the tidal forces are altered, but this effect should be accounted for by the definition of $x$. However, the directions of the tidal forces also change. In particular, Eq. (3) implies that there are stretching tidal forces proportional to $\mathrm{d}F/\mathrm{d}R$ along the radial axis from the host and compressive tidal forces proportional to $F/R$ along the perpendicular directions. When the host force profile $F(R)$ is self-similar, the ratio between the stretching and compressive forces is fixed.111111In fact, for an NFW profile, $\mathrm{d}F/\mathrm{d}R=0$ when $R\ll R_{s}$. Beyond the self-similar regime, however, the ratio between these forces can change. With this motivation, we define $f(R/R_{s})\equiv(\mathrm{d}F/\mathrm{d}R)/(F/R)$ as this ratio. For the NFW profile this definition implies that

[TABLE]

In the bottom panel of Fig. 10, we plot $b/(b_{0}x^{-b_{1}})$ against $y$ for the purpose of tuning the parameter $b_{2}$ in Eq. (19). We find that this equation121212The force-ratio argument motivates any expression of the form $\left[1+b_{2}f(y)^{\alpha}\right]^{\beta}$, but we assume for simplicity that $\alpha=\beta=1$. works reasonably well, and we obtain $b_{2}=1.29$. The introduction of stretching tidal forces increases the efficiency of tidal effects, which is reflected as an increase in the decay rate $b$ of the $J$ factor.

IV.2 Parameter $a$: The $J$-factor normalization

We next handle the overall normalization of $J/J_{\mathrm{init}}$. According to Eq. (14), the $J$ factor changes by the factor $\mathrm{e}^{ab}$ after the first orbit, which is sensitive to a second parameter: $a$. Because of the way we defined this parameter, it turns that $a$ is almost wholly sensitive to $y=\bar{R}/R_{s}$ alone. Figure 11 plots $a$ against $y$ for all of our simulations, and we find that with only moderate scatter,

[TABLE]

with $a_{0}=0.44$ and $a_{1}=1.32$.

In some sense, the parameter $a$ describes the initial behavior of the subhalo as it equilibrates—to the extent that this is possible—into the tidal field generated by the host. When $y\ll 1$ all nonzero tidal forces are compressive, so the $J$ factor is initially slightly boosted ($a>0$). However, when $y\gtrsim 1$, the stretching tidal forces cause the $J$ factor to be initially suppressed ($a<0$). Note that the trajectory given by Eq. (14) is only valid after this equilibration takes place, so it is not valid for $n<1$.

IV.3 Parameter $c$: The loss of tidal efficiency

Finally, we address the parameter $c$ that characterizes the drop in the efficiency of tidal effects over time. As we found in Sec. III, $c$ is sensitive to the orbit shape; more eccentric orbits yield smaller values of $c$ while more circular orbits yield larger values. In the self-similar regime we could write $c$ as a function of $\eta$, since $\eta$ completely describes the orbit shape. However, beyond this regime, orbits with the same $\eta$ could have different shapes. Thus, to accurately describe the sensitivity of the parameter $c$ to the host-subhalo system, it is necessary to find the correct orbit-shape parametrization.

We argued in Sec. III that the loss of tidal efficiency encoded in $c$ is related to changes in the shape of the subhalo density profile. The connection to the orbital shape is that circular orbits tidally heat material more predominantly in the outskirts of the subhalo, while eccentric orbits can alter the density profile further inward. This tendency is illustrated in Fig. 12, which depicts the tidally altered density profiles of two subhalos with different orbit shapes. The subhalo on the circular orbit loses more material from its outskirts, while the subhalo on the eccentric orbit loses more material from its interior.

Differences in the radii at which material is heated can be understood in terms of adiabatic shielding (e.g., Refs. Spitzer (1987); Weinberg (1994a, b); Gnedin and Ostriker (1999)). Deep within the subhalo, the internal dynamical timescale is much shorter than the timescale over which the external tidal field changes. In this case, the conservation of adiabatic invariants prevents any energy injection by tidal forces; these radii are adiabatically shielded. Meanwhile, adiabatic shielding is connected to the shape of the subhalo’s orbit. The timescale over which tidal forces change is related to the timescale of the pericenter passage, which can be very short for highly eccentric orbits.

Up to factors of order unity, the subhalo’s internal dynamical timescale is $t_{\mathrm{dyn}}\sim(Gm(r)/r^{3})^{-1/2}$ at radius $r$, where $m(r)$ is the subhalo mass profile Binney and Tremaine (1987). Meanwhile, the pericenter passage timescale is $t_{p}\sim R_{p}/V_{p}$, where $R_{p}$ and $V_{p}$ are the radius and velocity at the pericenter, respectively. To make precise the connection between the orbit shape and the radii at which tidal heating is efficient, we define the adiabatic shielding radius $r_{a}$ as the radius at which $t_{\mathrm{dyn}}=t_{p}$. This definition motivates a characteristic density scale

[TABLE]

so that $r_{a}$ is the radius at which $m(r_{a})/r_{a}^{3}=\rho_{a}$. Note that we used the definitions of the circular orbit radius $R_{c}$ and orbit circularity $\eta$ to eliminate $V_{p}$ from Eq. (23); $M(R)$ is the host mass profile at radius $R$.

To quantify changes in the shape of the subhalo density profile, we can compare the radius $r_{a}$ below which material is shielded to the tidal radius $r_{t}$ above which all material is stripped. The tidal radius is the radius above which the tidal force from the host exceeds the gravitational force from the subhalo. There are several definitions of the tidal radius in the literature, but they are all related to the expression Klypin et al. (1999, 2015) $r_{t}=R[m(r_{t})/M(R)]^{1/3}$ by (possibly nonconstant) factors of order unity (see, e.g., Ref. van den Bosch et al. (2018)). The tidal radius is only well defined for circular orbits, but it is common to apply the concept to eccentric orbits as well van den Bosch et al. (2018). In particular, if we seek the radius above which all material is stripped, we can define the tidal radius $r_{t}$ using the orbital apocenter radius $R_{a}$. In this case, there is a characteristic density scale

[TABLE]

and $r_{t}$ is the solution to $m(r_{t})/r_{t}^{3}=\rho_{t}$.

Anticipating that the drop in the efficiency of tidal effects encoded in the parameter $c$ is a consequence of changes to the shape of the subhalo density profile, we may hypothesize that $c$ is sensitive to the ratio

[TABLE]

which ranges from 0 for radial orbits to 1 for circular orbits. For simplicity, in defining $z$ we employ the subhalo’s initial NFW mass profile [see Eq. (35)]. Figure 13 shows the relationship between $c$ and $z$; there is some scatter, but the trend is that131313We argued in Sec. III.2 that $c>0$ is connected to changes in the density profile’s shape. In this light, Eq. (26) implies that in the limit $z=0$ where the tidal energy injection is completely impulsive, the shape of the density profile does not change over successive orbits. This notion is consistent with the results of Ref. Delos (2019), which found that impulsive point-object encounters yield a universal density profile.

[TABLE]

with $c_{0}=0.73$ and $c_{1}=0.21$. Note that $z$ is not solely a function of the subhalo’s orbit. Because it depends on the subhalo mass profile $m(r)$, it is also sensitive to the density ratio $\rho_{s}/P_{s}$. We explored using $\rho_{t}/\rho_{a}$, a purely orbital parameter, instead of $r_{a}/r_{t}$. This parameter exhibited a similar power-law relationship with $c$, but it left significant residual sensitivity to the parameter $x$, which is related to $\rho_{s}/P_{s}$. Using $z=r_{a}/r_{t}$ mostly eliminates that sensitivity.

V Model summary and discussion

In the last section, we developed a model for the evolution of a subhalo’s $J$ factor due to tidal effects as a function of parameters of the host-subhalo system. As long as the $J$ factor decays slower than $|\mathrm{d}\ln J/\mathrm{d}\ln t|=1$, its trajectory is well fit by the expression

[TABLE]

Here, $n=t/T$ is the number of subhalo orbits, and $a$, $b$, and $c$ are parameters that depend on the host-subhalo reduced system parameters $x$, $y$, and $z$ through

[TABLE]

with ${a_{0}=0.44}$, ${a_{1}=1.32}$, ${b_{0}=0.58}$, ${b_{1}=0.58}$, ${b_{2}=1.29}$, ${c_{0}=0.73}$, and ${c_{1}=0.21}$. In Appendix C, we detail how to compute $x$, $y$, $z$, and $T$ from the subhalo parameters $r_{s}$ and $\rho_{s}$, host parameters $R_{s}$ and $P_{s}$, and orbital parameters $R_{c}$ and $\eta$.

Equation (14) applies only when ${|\mathrm{d}\ln J/\mathrm{d}\ln t|<1}$. The $J$ factor’s precise behavior when $|\mathrm{d}\ln J/\mathrm{d}\ln t|\gtrsim 1$ is of little consequence, as the subhalos in this regime contribute only minimally to aggregate annihilation signals. Nevertheless, it is useful to have an approximate treatment in this regime. As we discussed in Sec. III, when $|\mathrm{d}\ln J/\mathrm{d}\ln t|\geq 1$ it is a reasonable approximation to enforce $-\mathrm{d}\ln J/\mathrm{d}\ln t=1$, i.e., $J\propto n^{-1}$. We define

[TABLE]

as the orbit count at which $-\mathrm{d}\ln J/\mathrm{d}\ln t=1$. Additionally, when $b>1$ (so $n_{1}<1$), we cannot expect our treatment of the normalization of $J$ (Sec. IV.2) to be accurate. To handle these issues, we can write

[TABLE]

where the last case follows from continuity considerations. We further note that this equation is valid only when $n\geq 1$ (see Sec. IV.2).

One can now use our model to understand the emission from a host halo due to dark matter annihilation in subhalos. In particular, one can sample subhalos from an orbital distribution in $R_{c}$ and $\eta$ (e.g., Refs. Tormen (1997); Zentner et al. (2005); Khochfar and Burkert (2006); Wetzel (2011); Jiang et al. (2015); van den Bosch (2017)). Accounting for tidal evolution, each subhalo’s contribution to the dark matter annihilation signal is then scaled by the orbit-dependent function given by Eq. (28). In Ref. Delos et al. (shed) (in preparation), we will use this model to study the annihilation signature arising from the extreme-density microhalos that result from certain early universe scenarios. In this case, the orbital distribution of subhalos is the same as that of particles, and one may employ the host halo’s distribution function (e.g., Ref. Widrow (2000)) to sample subhalo orbits.

Our model does not include the periodic oscillations in the $J$ factor observed in Sec. III. These oscillations do not affect the overall annihilation rate in subhalos, but they still introduce a systematic biasing effect where subhalos at smaller radii have larger $J$, and this effect can alter the morphology of an annihilation signal. However, we remark that these oscillations only have a significant amplitude in the $R_{c}\ll R_{s}$ regime, when all tidal forces are compressive, and at small $x\lesssim 10$. Because of these restrictions, we anticipate that their impact is minor. However, in forthcoming work Delos et al. (shed) we will quantify the impact of these oscillations.

We also address another potential limitation to our model. The differential equation driving it, Eq. (8), has explicit time dependence in the factor $n^{-c}$, so the resulting tidal evolution is not completely determined by the system’s instantaneous state. Physically, we view $n^{-c}$ as a proxy for unknown physical variables [e.g., Eq. (13)], and as long as the host halo’s density profile and the subhalo’s orbit are static, this formulation poses no difficulty. Since halos grow from the inside outward, subsequent accretion is not expected to significantly alter the density profile of a host halo at the radii of already-present subhalos, so the host halo is generally expected to remain static. Moreover, if dynamical friction is negligible [see Eq. (5)], the subhalo’s orbit is also static. However, there is a scenario where a subhalo’s host is expected to change dramatically. If the host is itself a satellite of a larger host halo, then the subhalo may be tidally stripped from its host, becoming itself a satellite of the superhost. In this scenario, it is not obvious how to continue the subhalo’s tidal evolution.

If the initial host-subhalo system yields trajectory parameters $a^{\prime}$, $b^{\prime}$, and $c^{\prime}$ and the new host-subhalo system yields parameters $a$, $b$, and $c$, then a self-consistent way to treat this problem is to substitute the factor $n^{-c}$ in Eq. (8) with $(n+n^{\prime c^{\prime}/c})^{-c}$ and integrate the resulting expression. This treatment follows from the assumption that the parameter $q=n^{-c}$ in Eq. (13) is a function of the subhalo alone. Additionally, the $J$ factor should be rescaled by $\mathrm{e}^{ba-b^{\prime}a^{\prime}}$, a consideration motivated by the discussion in Sec. IV.2. However, it turns out that while this treatment works reasonably well for a portion of the $a^{\prime}$, $b^{\prime}$, $c^{\prime}$, $a$, $b$, $c$ parameter space, it does not accurately predict every scenario; the parameter $q$ in Eq. (13) is not a function of the subhalo alone. We leave a detailed investigation of this problem to future work.

As another caveat, the long-term accuracy of the trajectory in Eq. (14) relies on the assumption that the efficiency of tidal effects follows precisely the power law $n^{-c}$, as described by Eq. (8). While such a power law is a natural assumption [e.g, Eq. (13)] and is borne out in our simulations, it does not have a direct physical motivation; tidal heating models considered in Sec. III.2 and elsewhere Pullen et al. (2014) can only reproduce $c=0$. Without such motivation it is unclear that this power-law behavior should extend beyond the $n=20$ orbits of our longest simulations. Also, the system parameter $z=r_{a}/r_{t}$ that sets the power-law index $c$ is defined based on two concepts that are not themselves entirely well defined: the adiabatic shielding radius $r_{a}$ and the tidal radius $r_{t}$. Moreover, since our model does not predict the larger evolution of the subhalo density profile, we use the subhalo’s initial density profile to define $r_{a}$ and $r_{t}$ even though the density profile quickly begins to change. For these reasons, we anticipate that it is possible to find a better-motivated parameter to replace $z=r_{a}/r_{t}$.

Nevertheless, this model describes the results of our simulations with remarkable success. As further validation, we consider the library of idealized subhalo simulations, called DASH (for dynamical aspects of subhaloes), published by Ref. Ogiya et al. (2019). These simulations have a lower resolution than ours, but because of the extraordinary volume of this library, it still supplies a valuable test for our model. In Appendix D we verify that modulo substantial scatter and certain systematic effects associated with their lower resolution, the DASH simulations are consistent with our model.

VI Comparison to previous work

Numerous prior works have endeavored to model the impact of tidal effects on a subhalo’s dynamical evolution Taylor and Babul (2001); Hayashi et al. (2003); Peñarrubia and Benson (2005); van den Bosch et al. (2005); Zentner et al. (2005); Kampakoglou and Benson (2007); Gan et al. (2010); Penarrubia et al. (2010); Pullen et al. (2014); Jiang and van den Bosch (2016). In this section, we explore how our results compare to those of previous studies. Motivated primarily by simulations, our model is based on the notion that a subhalo’s $J$-factor evolution is determined by

[TABLE]

where $b$ and $c$ are functions of the host-subhalo system and $n$ counts the number of orbits. In contrast to our focus on the $J$ factor, previous works have largely focused on the evolution of a subhalo’s total bound mass $m_{\mathrm{bound}}$ and of its maximum circular velocity $v_{\mathrm{max}}$ and corresponding radius $r_{\mathrm{max}}$. However, the general form of our model is not specific to the $J$ factor, and we show in Appendix E that it can also describe the evolution of $v_{\mathrm{max}}$, $r_{\mathrm{max}}$, and $m_{\mathrm{bound}}$.

Despite the broad applicability of our model suggested by Appendix E, no prior work (to our knowledge) has proposed tidal evolution of the form given in Eq. (8). Broadly, prior models of tidal evolution fall into two main categories, although a given work may employ more than one:

1. (1)

Tidal stripping models, where material outside the characteristic tidal radius [e.g., Eq. (24)] is assumed to be stripped over some time period; 2. (2)

Tidal heating models, where energy injected by tidal forces heats subhalo material, causing it to rise and possibly become freed from the subhalo.

We found in Sec. IV.1 that the parameter $b$ in Eq. (8), which characterizes the rate of tidal evolution, is tightly sensitive to the energy injected by tidal forces (see Fig. 9). Additionally, we observed in Sec. III that the tidal evolution in our simulations closely resembles that predicted by a toy model of tidal heating. Consequently, we anticipate that of the two classes of models, tidal heating models should yield results most similar to those of our model. We will first compare the results of our model to those of the tidal heating model developed by Ref. Pullen et al. (2014), hereafter P14.

However, prior treatments of dark matter annihilation within subhalos predominantly treat the impact of the host halo’s tidal forces using models based on tidal stripping Bartels and Ando (2015); Hiroshima et al. (2018); Ando et al. (2019); Springel et al. (2008); Erickcek (2015); Stref and Lavalle (2017); Stref et al. (2019). Tidal stripping models cannot prescribe how to change a subhalo’s density profile below the tidal radius, but it is possible to apply a simulation-tuned prescription for how the density profile responds to mass loss Hayashi et al. (2003); Penarrubia et al. (2010). We will subsequently compare the results of our model to those of a tidal stripping model developed by Ref. Jiang and van den Bosch (2016) (hereafter J16), using the prescription of Ref. Penarrubia et al. (2010) (hereafter P10) to predict the subhalo’s density profile. This pair of models has been employed by Refs. Bartels and Ando (2015); Hiroshima et al. (2018) to predict dark matter annihilation rates in subhalos.

VI.1 Comparison to a tidal heating model

We first compare our model’s predictions to those of the analytic tidal heating model given in P14. In the tidal heating picture, energy injected by tidal forces causes subhalo material to move to higher radii, and P14 employed the assumption of virial equilibrium to predict this change in radius and consequently the subhalo’s new density profile. We follow the prescription in P14141414For simplicity, we compute the energy injection $\Delta E$ directly using the impulse approximation (Appendix C), neglecting additional corrections suggested in P14; these corrections will not qualitatively alter the results. to compute the evolution of a subhalo’s density profile, subsequently integrating it to obtain the $J$ factor. Figure 14 shows a sample of the resulting $J$-factor trajectories, and we compare those trajectories to our model’s predictions and to the results of our simulations. Generally, we find that for a model constructed from first principles, the P14 model is remarkably accurate. However, it does not fully capture the sensitivity of tidal evolution to system parameters, a matter we explore next.

The quantity $Q=\Delta E/r^{2}$ employed by P14 is related to host-subhalo system parameters by $Q/(G\rho_{s})=4\pi(\ln 2)N/x$, where $N$ is the number of orbits over which the energy injection is taken and $x$ is the system parameter (see Sec. IV). For our comparison, we take $N=1$ and iterate the calculation, assuming the halo revirializes during each orbit. Since the density profile evolution in P14 is only sensitive to the ratio $Q/(G\rho_{s})$, we see immediately that this model’s predictions are sensitive only to the system parameter $x$ and are insensitive to $y$ and $z$. Additionally, the $J$-factor evolution predicted by P14 turns out to be only sensitive to $x$ in the combination $n/x$, where $n=t/T$ is the number of orbits, so every system follows the same trajectory rescaled in time. In this respect, the P14 model is similar to the toy model we explored in Sec. III.2, which was only sensitive to the combination $fn$ of system parameters $f$ and orbit count $n$. In fact, the P14 model approximately obeys the toy model solution Eq. (12) with $b\simeq 3.2/x$ and $B=1$, but it can potentially transition between the $n/x\ll 1$ and $n/x\gg 1$ regimes extremely slowly, and all behavior seen in Fig. 14 is in the intermediate regime.

The combination of its single time-rescaled trajectory and its insensitivity to $y$ and $z$ leaves the P14 model unable to accurately predict tidal evolution in the full host-subhalo parameter space. We see evidence of this deficiency in Fig. 14, but we further note that we did not plot any subhalos in the $y\ll 1$ regime. In this regime, the P14 model dramatically overestimates the impact of tidal stripping since it does not account for the directions of tidal forces, which are encapsulated in the parameter $y$. While the P14 model yields reasonably accurate predictions over a small range of host-subhalo system parameters, our model can accurately predict the evolution of a much broader variety of systems.

VI.2 Comparison to a tidal stripping model

Finally, we compare our model to a semianalytic model of tidal stripping that has been employed in previous calculations of annihilation rates in the substructure Bartels and Ando (2015); Hiroshima et al. (2018). This semianalytic model uses the tidal stripping model in J16 to characterize a subhalo’s mass loss, subsequently using the results of P10 to connect this mass loss to the subhalo’s density profile and hence annihilation signal. In J16, the rate of mass loss for subhalos of mass $m$ inside a host halo of mass $M$, averaged over subhalo orbits, is modeled using

[TABLE]

Here, $\mathcal{A}$ and $\zeta$ are simulation-tuned parameters and $t_{\mathrm{dyn}}$ is the host’s dynamical timescale at its virial radius (e.g., Ref. Binney and Tremaine (1987)). For this comparison we adopt J16’s central values $\mathcal{A}=0.86$ and $\zeta=0.07$.

In P10, it is shown that the subhalo’s maximum circular velocity $v_{\mathrm{max}}$ and the radius $r_{\mathrm{max}}$ at which it is attained are related to the fraction $m/m_{\mathrm{acc}}$ of the subhalo’s mass that remains gravitationally bound, where $m_{\mathrm{acc}}$ is the subhalo’s virial mass at accretion. We confirm in Appendix E that these relations are reasonably accurate if $y>1$ and the subhalos have concentration $r_{\mathrm{vir}}/r_{s}\simeq 20$ at accretion. If we assume that subhalos possess NFW profiles, then the mass fraction $m/m_{\mathrm{acc}}$ predicted by J16 thereby determines each subhalo’s $J$ factor.

To compare our model, we employ the same subhalo orbital distribution considered in J16, which is drawn from Ref. Zentner et al. (2005). The circular orbit radius $R_{c}$ is taken to be uniformly distributed between $0.6R_{\mathrm{vir}}$ and $R_{\mathrm{vir}}$, where $R_{\mathrm{vir}}$ is the host’s virial radius. Meanwhile, the circularity $\eta$ is distributed proportionally to $\sin\pi\eta$, and we assume that the distributions of $R_{c}$ and $\eta$ are independent. By drawing subhalo orbits from this distribution, we are able to compute the orbit-averaged151515Specifically, we take the median $J/J_{\mathrm{init}}$ at each time, but using the mean or the logarithmic mean instead does not significantly alter the results. value of $J/J_{\mathrm{init}}$ using our model given in Eq. (28). In Fig. 15, we plot the resulting orbit-averaged $J$-factor trajectories along with those predicted by the semianalytic model of J16 and P10. We consider two different host-subhalo systems, listed on the figure, and since the semianalytic model is sensitive to the total virial masses of the host and subhalo, we employ the concentration parameter $c\equiv r_{\mathrm{vir}}/r_{s}$ to describe these systems; $c_{\mathrm{host}}$ is the host halo’s concentration, while $c_{\mathrm{sub}}$ is the subhalo’s concentration when it is accreted.

Compared to our model, Fig. 15 shows that the semianalytic model underestimates the impact of tidal forces early on while overestimating their impact at late times. These discrepancies arise from several sources. As we show in Appendix E, a subhalo’s $J$ factor after tidal evolution is about 30% smaller than what would be predicted from its parameters $r_{\mathrm{max}}$ and $v_{\mathrm{max}}$ assuming an NFW profile. However, this source of error is relatively minor. The main differences arise from the model in J16 given by Eq. (29). Since this model does not account for the subhalo’s density profile, it takes too long to strip the subhalo’s weakly bound outskirts (beyond $r_{s}$) that contribute little to annihilation rates.161616For the two cases shown in Fig. 15, it takes, respectively, 6 and 3 dynamical times for the J16 model to bring the subhalo’s bound mass below its initial $m_{\mathrm{max}}$, the mass enclosed within the radius $r_{\mathrm{max}}$ at which the maximum circular velocity is attained. As we find in Appendix E (see Fig. 21), the drop in the $J$ factor is minimal above this mass threshold. This behavior partially explains why the semianalytic model underestimates the early impact of tidal effects. Meanwhile, for small $\zeta\ll 1$, Eq. (29) describes nearly exponential decay, analogous to our model, Eq. (8), with $c=0$. Without the braking behavior contributed by $c>0$ (and attributed to changes in the shape of the subhalo’s density profile; see Sec. III) along with the limiting $|\mathrm{d}\ln J/\mathrm{d}\ln t|\sim 1$ behavior, the semianalytic model overestimates the impact of tidal effects at late times. For these reasons, our simulation-tuned model supplies significantly more accurate predictions of subhalo annihilation rates.

VII Conclusion

In this work, we used 52 idealized $N$-body simulations to develop a model that can predict the impact of a host halo’s tidal forces on the rates of dark matter annihilation within its subhalos. Our model is given by Eq. (28) and summarized in Sec. V, and it predicts the evolution of the subhalo’s $J$ factor, the factor in the annihilation rate that depends on mass distribution, as a function of the subhalo’s orbit and other properties of the host-subhalo system. These properties are distilled into three physically motivated variables $x$, $y$, and $z$ that characterize the energy injected by tidal forces, the ratio of stretching to compressive tidal forces, and the radial distribution of tidally heated material, respectively. Appendix C details how to compute these variables from standard properties of the host-subhalo system.

Our model is based on the notion that for sufficiently small changes in $J$, the $J$ factor evolves according to

[TABLE]

where $n=t/T$ is the time in units of the subhalo’s orbital period and $b$ and $c$ are parameters that depend on the system. If $c=0$, Eq. (8) states that the subhalo loses a fixed fraction $\mathrm{e}^{-b}$ of its $J$ factor in each orbit. The parameter $c\geq 0$ is motivated by simulation results and adds a braking mechanism to the $J$ factor’s decay. To our knowledge, a model of this form has not previously been put forward, even though we find that it can also describe other structural properties of the subhalo. We also find that our model predicts significantly different $J$-factor trajectories than prior semianalytic models. We further validate our model by testing it against the publicly available DASH library of subhalo simulations Ogiya et al. (2019), finding reasonable agreement.

Our model has limitations. As presented, it is restricted to host-subhalo systems in which both halos possess NFW density profiles. The NFW profile (possibly with minor corrections; e.g., Ref. Navarro et al. (2010)) arises generically in dark matter simulations of halos built by hierarchical clustering Navarro et al. (1996, 1997). However, the smallest subhalos, forming by direct collapse, exhibit steeper density profiles Ishiyama et al. (2010); Anderhalden and Diemand (2013a); *anderhalden2013erratum; Ishiyama (2014); Polisensky and Ricotti (2015); Angulo et al. (2017); Ogiya and Hahn (2018); Delos et al. (2018b, a, 2019). Additionally, the density profiles of many galactic halos (but not all Read et al. (2018)) are inferred to be shallower than the NFW profile, an observation that may be explained by baryonic effects or unknown dark matter properties (see Refs. Brooks (2014); Bullock and Boylan-Kolchin (2017) for reviews). Despite being developed using NFW profiles, we anticipate that the physical manner in which we defined the model parameters $x$, $y$, and $z$ implies that our model can be adapted to accommodate different host or subhalo density profiles.

Also, our model only accounts for tidal forces from the host halo. Subhalos can also be disrupted by encounters with other subhalos, but the results of Ref. van den Bosch et al. (2018) suggest that this effect is subdominant. More importantly, subhalos can be affected by baryonic content residing within the host, such as stars (e.g., Refs. Berezinsky et al. (2003, 2006); Goerdt et al. (2007); Zhao et al. (2007); Green and Goodwin (2007); Schneider et al. (2010); Ishiyama et al. (2010); Angus and Zhao (2007); Berezinsky et al. (2014); Delos (2019)) or a disk (e.g., Refs. Taylor and Babul (2001); Angus and Zhao (2007); Schneider et al. (2010); D’Onghia et al. (2010); Berezinsky et al. (2014); Zhu et al. (2016); Errani et al. (2016); Garrison-Kimmel et al. (2017); Stref and Lavalle (2017); Kelley et al. (2019); Hütten et al. (2019)). These effects are not included in our model and must be accounted for separately. However, we remark that many of the dwarf spheroidal galaxies, already some of the most promising targets for dark matter annihilation searches Strigari et al. (2007), have such little baryonic content (e.g., Ref. McConnachie (2012)) that it may be possible to neglect the influence of this content on their subhalos.

Despite these limitations, we anticipate that our model will prove useful in understanding the annihilation signals of dark matter substructure. In a subsequent paper Delos et al. (shed), we will explore the consequences of our model by using it to study microhalo-dominated annihilation signals in nearby dwarf galaxies. Such signals are expected to arise from certain cosmological scenarios, such as an early matter-dominated era prior to nucleosynthesis, and our model enables precise characterization of the magnitude and morphology of these signals.

Acknowledgements.

The simulations for this work were carried out on the Killdevil and Dogwood computing clusters at the University of North Carolina at Chapel Hill. This work was funded by NASA through the Fermi Guest Investigator Cycle 10 Award No. 80NSSC17K0751 (PI A. Erickcek). The author thanks Adrienne Erickcek and Tim Linden for helpful discussions. Key figures in this work employ the cube-helix color scheme developed by Ref. Green (2011).

Appendix A Simulation details

A.1 High- and low-resolution particles

As Sec. II notes, we sample the subhalo’s central region at increased resolution such that particles whose orbital pericenters are below $r_{s}/3$ have $1/64$ the mass and $64$ times the number density of the other particles. When simulation particles have different masses, it is possible for two-body interactions to artificially transfer energy from the heavy to the light particles. To verify that this effect is not significant in our simulations, we show in Fig. 16 the density profiles of light and heavy particles in a subhalo not exposed to tidal forces. Even after duration $t=318(G\rho_{s})^{-1/2}$, where $\rho_{s}$ is the subhalo’s scale density, there is no visible tendency for the heavy particles to sink to smaller radii.

A.2 Density profiles and $J$ factors

We obtain each subhalo’s density profile by binning it in factors of $1.1$ in the radius. At small radii, there is a resolution limit driven by three effects: force softening, Poisson noise, and artificial relaxation. Each effect is associated with a minimum resolved radius below which the density profile artificially flattens. For force softening, that radius is the distance $r_{\mathrm{soft}}=2.8\epsilon$, where $\epsilon$ is Gadget-2’s force-softening parameter, at which forces become non-Newtonian. For Poisson noise, we take it to be the radius $r_{100}$ enclosing 100 particles. To estimate the radius $r_{\mathrm{rel}}$ at which artificial relaxation becomes significant, we compute the relaxation time Binney and Tremaine (1987)

[TABLE]

at each radius $r$, where $M$ and $N$ are the mass and particle count interior to $r$, and $\Lambda=\max\{N,r/\epsilon\}$. If $\alpha t_{\mathrm{relax}}$ at radius $r$, with an appropriate proportionality constant $\alpha$, is shorter than the system age, then $r<r_{\mathrm{rel}}$. The proportionality constant $\alpha$ is tuned to predict the correct $r_{\mathrm{rel}}$ in a simulation of the subhalo without a host; in that case, any change to the density profile is artificial since the halo was built from an equilibrium distribution. From this calibration we use $\alpha=5$. Figure 16 shows how $r_{\mathrm{rel}}$ marks where the density profile begins to shallow due to relaxation effects.

The minimum resolved radius of the density profile is $r_{\mathrm{min}}=\max\{r_{\mathrm{soft}},r_{100},r_{\mathrm{rel}}\}$. For the purpose of accurately computing $J$ factors we extrapolate the density profile below $r_{\mathrm{min}}$ as $\rho=Ar^{-1}$, where $A$ is the average of $\rho r$ in the three smallest radial bins above $r_{\mathrm{min}}$, so that

[TABLE]

Effectively, this procedure produces a lower bound on $J$ under the assumption that larger radii are always stripped more than smaller radii. Below $r_{\mathrm{min}}$, we simply assume all radii are stripped equally. We can also compute an upper bound on $J$ by assuming that radii below $r_{\mathrm{min}}$ are completely unaffected (so $A=\rho_{s}r_{s}$), and this allows us to estimate the uncertainty in our $J$ factors. We find that by the termination time of each subhalo’s $J$-factor trajectory (as defined in Sec. IV), the uncertainty in the $J$ factor, taken as $J_{\mathrm{upper}}/J_{\mathrm{lower}}-1$, is 31% for one simulation (parameters $x=31$, $y=11$, $z=0.07$; see Sec. IV), smaller than 17% for the remaining 51 simulations, and smaller than 10% for 44 of them.

To understand the $J$ factors in the cases where $|\mathrm{d}\ln J/\mathrm{d}\ln t|\gtrsim 1$, another step is necessary. In this regime, the elongated tidal stream can contribute significantly to the $J$ factor, making the spherical integral Eq. (31) inaccurate. Thus, we also compute the $J$ factor as the sum over simulation particles

[TABLE]

where $m_{i}$ is the mass of particle $i$ and $\rho_{i}$ is its local density. The density $\rho_{i}$ is estimated as

[TABLE]

over the $N=50$ nearest particles $j$, where $r_{ij}$ is the distance to particle $j$, $h_{i}$ is the distance to the $N$th particle, and $W(r,h)$ is the cubic spline kernel defined as in Ref. Springel (2005).

Equation (32) underestimates the $J$-factor contribution at $r<r_{\mathrm{min}}$ due to artificial flattening of the density profile. To accommodate the extrapolation procedure in Eq. (31) that addresses this problem, an additional step is required. We find the bound remnant of the subhalo using a procedure similar to that in Ref. van den Bosch et al. (2018). Beginning with the assumption that all particles are bound, we iteratively compute the gravitational potential of each particle due to all other bound particles using a Barnes-Hut octree Barnes and Hut (1986) with $\theta=0.7$ and the same softening length as the simulation. Subsequently, we mark each particle as unbound if its total energy is positive and bound if its total energy is negative. At each step, we find the center-of-mass position and velocity of the 100 most bound particles and recenter the full system to be relative to this center of mass. All particles are initially marked as bound, and the procedure terminates when the count of bound particles converges.171717This halting condition is stricter than the one in Ref. van den Bosch et al. (2018).

By assuming that the bound remnant is spherically symmetric, we can estimate the $J$ factor both including spherical asymmetry and compensating for the flattening of the density profile below $r_{\mathrm{min}}$. If $J_{\mathrm{full}}$ is the $J$ factor of the full system computed using Eq. (32) and $J_{\mathrm{bd,rad}}$ and $J_{\mathrm{bd}}$ are the $J$ factors of the bound remnant computed using Eqs. (31) and (32), respectively, then

[TABLE]

A.3 Numerical convergence

In our simulations, we set Gadget-2’s force-softening length to be $\epsilon=0.003r_{s}$. This small value is intended to evade the artificial subhalo disruption observed by Ref. van den Bosch and Ogiya (2018). Meanwhile, the subhalo’s high-resolution particles (see Sec. II) have mass $$4.3\text{\times}{10}^{-7}$\rho_{s}r_{s}^{3}$. To check that numerical artifacts in our simulations are under control, we test the impact of changing the softening length and the particle resolution. Additionally, we test the impact of altering the (adaptive) integration time steps in order to ensure there are no artifacts arising from the application of the host’s tidal field over these discrete intervals. In Fig. 17, we plot the $J$-factor trajectory in a reference simulation (with system parameters $x=34$, $y=0.018$, and $z=0.15$; see Sec. IV) along with three simulations of the same system with different particle resolution, force softening, and integration time steps. We plot the upper and lower limits of the $J$-factor trajectory as discussed above. These limits overlap for all simulation parameters, suggesting that the simulation is converged.

Appendix B Subhalo size

In our simulations we applied the host halo’s tidal forces using the linearized expression given by Eq. (3), which is valid in the limit that the subhalo is much smaller than its orbital radius. Thus, our results are applicable in the $r_{s}\ll R_{c}$ limit. In this appendix, we explore precisely how far the applicability of our results can be taken. For this purpose we executed several simulations using the exact tidal force $\bm{F}_{\mathrm{tidal}}(\bm{r})=\bm{F}(\bm{R}+\bm{r})-\bm{F}(\bm{R})$ instead of the linearized version, and we compare the results of these simulations to those of a simulation that employed the linearized force. The subhalos in all of these simulations are cut off at radius $5r_{s}$ to avoid excessive overlap with the host’s center; this change does not affect the comparison since it applies equally to every simulation.

We subjected subhalos of different sizes to the same tidal evolution scenario with parameters $x=21$, $y=0.056$, and $z=0.18$ in the parametrization given in Sec. IV. Figure 18 shows the tidal evolution of the subhalo’s $J$ factor for different values of the ratio $r_{s}/R_{c}$ between the subhalo’s scale radius and the radius of its orbit about the host. We find that the tidal evolution begins to diverge markedly from that induced by the linearized tidal force when $r_{s}\gtrsim 0.1R_{c}$. Note that this analysis still neglects dynamical friction (including self-friction, due to the limited duration of these simulations); the influence of this effect is also sensitive to the subhalo’s size. Thus, our results are applicable if both $r_{s}\lesssim 0.1R_{c}$ and dynamical friction can be neglected (see Sec. II).

Appendix C Computational details

In this appendix, we present practical ways to compute the reduced variables $x$, $y$, and $z$ along with the orbital period $T$. For convenience, we include fitting formulas to approximate the necessary integrals. In what follows, the host is assumed to possess an NFW profile with scale radius $R_{s}$ and scale density $P_{s}$; its mass profile is

[TABLE]

its force profile is $F(R)=GM(R)/R^{2}$, and its potential profile is

[TABLE]

Meanwhile, the subhalo’s orbit about the host is parametrized by the circular orbit radius $R_{c}$ and circularity $\eta$, and as shorthand, we define $y_{c}\equiv R_{c}/R_{s}$.

C.1 Computing $x=|E_{b}|/\Delta E_{\mathrm{imp}}$

The binding energy $E_{b}$ of a particle at the subhalo’s scale radius $r_{s}$ is given by Eq. (15). Meanwhile, the total energy $\Delta E_{\mathrm{imp}}$ injected into a particle at radius $r$ by tidal forces over the course of a subhalo orbit is computed using the impulse approximation, as described in Ref. Gnedin et al. (1999). This energy depends on the particle’s full three-dimensional position within the subhalo, but we simplify the picture by averaging this energy over the sphere at radius $r$. Dimensionally, ${\Delta E_{\mathrm{imp}}/r^{2}\sim F(R_{c})/R_{c}}$, and we can approximate

[TABLE]

where $P_{1}(y_{c})$ is defined

[TABLE]

$P_{2}(y_{c})$ is defined

[TABLE]

and $P_{3}(y_{c})$ is defined

[TABLE]

For $\eta>0.04$, this expression is accurate to within 3% for $y_{c}<10$ and within 14% for $y_{c}<10^{3}$.

C.2 Computing $y=\bar{R}/R_{s}$

The time-averaged radius $\bar{R}$ of the orbit is approximately $R_{c}$, and, in fact, $\bar{R}/R_{c}\to 1$ as $R_{c}/R_{s}\to 0$. More broadly, the expression

[TABLE]

is accurate to within 0.3% for $\eta>0.04$ and $y_{c}<10^{3}$.

C.3 Computing $z=r_{a}/r_{t}$

We define the adiabatic shielding radius $r_{a}$ and the tidal radius $r_{t}$ as the solutions to $m(r_{a})/r_{a}^{3}=\rho_{a}$ and $m(r_{t})/r_{t}^{3}=\rho_{t}$, respectively, where $m(r)$ is the subhalo’s initial NFW mass profile [see Eq. (35)]. Here, $\rho_{a}$ and $\rho_{t}$ are functions of the subhalo’s orbit; in particular,

[TABLE]

where $R_{p}$ and $R_{a}$ are the orbital pericenter and apocenter radii, which may be obtained as the two solutions $R$ to

[TABLE]

C.4 Computing $T$

Dimensionally, the radial orbit period $T\sim t_{0}$, where ${t_{0}\equiv\sqrt{R_{c}/F(R_{c})}}$. More precisely, the expression

[TABLE]

is accurate to within 0.2% for $\eta>0.04$ and $y_{c}<10^{3}$.

Appendix D Comparison to the DASH library

Reference Ogiya et al. (2019) published a library called Dynamical Aspects of SubHaloes (DASH) of idealized subhalo simulations. This library includes the results of 2177 simulations, with different system parameters, of an $N$-body subhalo orbiting an analytic host potential. These simulations resolve significantly less of the subhalo density profile than do ours; as shown in Fig. 19, they can leave large fractions of the $J$ factor unresolved. Also, the DASH library covers a smaller parameter range in $x$, $y$, and $z$. Nevertheless, because of the extraordinary volume of this library, it can serve as a test for our model.

We use the procedure in Appendix A to find the $J$-factor trajectory of each DASH simulation,181818We use a larger $\alpha=20$ to find $r_{\mathrm{rel}}$ for the DASH simulations, obtained by recalibrating for these simulations. Note that larger $\alpha$ implies more optimism about simulation resolution. imposing an additional constraint that the trajectory halt when the maximum uncertainty in the $J$ factor is larger than a factor of 3. Next, we fit the parameters $a$, $b$, and $c$ to this trajectory as in Sec. IV. For the DASH simulations, Fig. 20 plots (in the same way as Figs. 10, 11, and 13) the trajectory parameters $a$, $b$, and $c$ against the system parameters $x$, $y$, and $z$. Superposed are our model predictions, as solid lines, using the parameters obtained in Sec. IV.

We first remark that all DASH simulations have $y>2$, so we cannot directly test Eq. (18) describing the behavior of $b$ in the $y\ll 1$ self-similar regime. Nevertheless, the upper-left panel of Fig. 20 shows that the DASH simulations exhibit roughly the same power-law behavior $b\propto x^{-0.58}$ predicted by Eq. (19) (the offset between our curve and the simulations here is not a discrepancy). The lower panels show the sensitivity of $a$ and $b$ to $y$. Because the DASH simulations only cover a small range of $y$, we cannot verify the functional form of each parameter in $y$. Also, there is substantial scatter, especially at large $x$. Nevertheless, our model predicts roughly the correct values of $a$ and $b$ for these simulations, although there is a tendency for the simulations to have smaller values of $a$ and larger values of $b$. Finally, although the scatter in $c$ is quite large, the relationship between $c$ and $z$ is approximately borne out in the DASH simulations.

The tendency for the DASH simulations to yield small $a$ and large $b$ can be understood as a resolution artifact. Below the resolution limit, we extrapolate the density profile in a way that always underestimates the $J$ factor (see Appendix A). This underestimation both increases the immediate loss of the $J$ factor, reducing $a$, and increases the rate at which $J$ decays (since artificial relaxation worsens the resolution over time), raising $b$.

Also, there is a tendency for systems at the large-$x$ end to exhibit large scatter in $a$, $b$, and $c$ as well as a precipitous drop in $b$ (sometimes even to $b<0$). This trend is also an unphysical artifact. In our simulations, we observed the same trend when $x\gtrsim 200$, which is why our simulation sample in Sec. IV only includes $x<200$. For the lower-resolution DASH simulations, the trend begins at $x\gtrsim 50$. The numerical difficulty with large $x$ is unclear, but it is likely connected to the fact that large $x$ implies the subhalo’s internal forces are much stronger than the external tides. The vast difference in the scales of these forces could lead to issues in numerical precision when the tiny tidal forces are added to the large internal forces.

Appendix E The broader density profile

References Hayashi et al. (2003) (hereafter H03) and Penarrubia et al. (2010) (hereafter P10) studied the tidal evolution of a subhalo’s density profile, focusing on the structural parameters $v_{\mathrm{max}}$, the maximum circular velocity within the subhalo, and $r_{\mathrm{max}}$, the radius at which this velocity is attained. Prior treatments of the annihilation rate in subhalos (e.g., Refs. Bartels and Ando (2015); Hiroshima et al. (2018)) have employed these works’ predictions of $r_{\mathrm{max}}$ and $v_{\mathrm{max}}$, along with the assumption that subhalos retain NFW profiles, to predict subhalo $J$ factors. To understand the connection between our work and these prior works, this appendix investigates the evolution of $r_{\mathrm{max}}$ and $v_{\mathrm{max}}$ in our simulations.

At each snapshot of our simulations, we find $r_{\mathrm{max}}$ as the radius $r<r_{t}$ that maximizes $v_{\mathrm{circ}}=\sqrt{Gm(r)/r}$, and $v_{\mathrm{max}}$ is the corresponding maximum. We only consider snapshots up to the point where $|\mathrm{d}\ln J/\mathrm{d}\ln t|=1$, as discussed in Sec. IV. Additionally, we halt the $r_{\mathrm{max}}$ and $v_{\mathrm{max}}$ trajectory when $r_{\mathrm{max}}$ becomes smaller than the resolution limit (see Appendix A), and we only include simulations whose trajectories cover at least five orbits about the host. This restriction reduces our simulation count to $41$, $33$ of which are in the self-similar regime ($R_{c}/R_{s}<0.3$).

E.1 Relations between structural parameters

Following H03 and P10, we first explore the relationship between a subhalo’s structural parameters and its total mass loss. These prior works parametrize the mass loss using the ratio $m_{\mathrm{bound}}/m_{\mathrm{acc}}$, where $m_{\mathrm{bound}}$ is the mass that remains bound to the subhalo and $m_{\mathrm{acc}}$ is its virial mass at accretion. However, this parametrization implies that the impact of tides is strongly sensitive to the subhalo’s initial concentration, and we propose that this sensitivity is unphysical since the outer layers may be stripped almost immediately upon accretion onto the host. To evade this problem, we instead parametrize the mass loss using the ratio $\tilde{m}\equiv m_{\mathrm{bound}}/m_{\mathrm{max,init}}$ of $m_{\mathrm{bound}}$ to the mass initially enclosed within $r_{\mathrm{max}}$; this ratio is initially larger than unity. We compute $m_{\mathrm{bound}}$ using the procedure in Appendix A, and Fig. 21 shows these relationships. For comparison we also plot the predictions of H03 and P10 assuming $m_{\mathrm{acc}}=4.5m_{\mathrm{max,init}}$, which corresponds to subhalo concentration $c_{\mathrm{sub}}\simeq 20$ at accretion.

We find that $r_{\mathrm{max}}$ is cleanly related to $m_{\mathrm{bound}}$ by a power law. Additionally, for both $r_{\mathrm{max}}$ and $v_{\mathrm{max}}$, the predictions of P10 (with $c_{\mathrm{sub}}\simeq 20$) work reasonably well as long as $\bar{R}>R_{s}$. However, the bottom panel of Fig. 21 shows that there is substantial scatter in the relationship between $J$ and $m_{\mathrm{bound}}$, and P10 does not accurately predict the $J$ factor if an NFW profile is assumed. The scatter partially results from the modest scatter in $v_{\mathrm{max}}$, since $J\propto v_{\mathrm{max}}^{4}$, but it also reflects that tidally altered density profiles differ significantly from NFW. Figure 22 investigates this effect further and shows that a subhalo’s $J$ factor is roughly 30% smaller than what would be predicted from $r_{\mathrm{max}}$ and $v_{\mathrm{max}}$ assuming an NFW profile.

E.2 Time evolution of structural parameters

We can also predict the evolution of the subhalo’s structural parameters more explicitly. Subjected to tidal forces, $r_{\mathrm{max}}$ and $v_{\mathrm{max}}$ follow qualitatively similar trajectories to the $J$ factor:

[TABLE]

where $n=t/T$ is the number of orbits [compare Eq. (14)]. Note that the total bound mass, $m_{\mathrm{bd}}$, behaves similarly; its trajectory follows from Eq. (45) by inverting the equation in the top panel of Fig. 21. As shown in Fig. 23, $b_{r}$, $c_{r}$, $b_{v}$, and $c_{v}$ appear to depend on the system parameters $x$, $y$, and $z$ in the same way that $b$ and $c$ did:

[TABLE]

However, the middle panels of Fig. 23 show that unlike $a$, the parameters $a_{r}$ and $a_{v}$ depend not only on $y$ but also on $x$. Moreover, they are only sensitive to $x$ in the self-similar regime ($R_{c}/R_{s}<0.3$). We fit the equation $a_{r}=a_{r0}\ln(x/a_{r1})$, and likewise for $a_{v}$, in the self-similar regime. Next, we fit $a_{r}-a_{r0}\ln(x/a_{r1})[1-f(y)/2]=-a_{r2}f(y)$, and likewise for $a_{v}$, using all simulations. The function $f(y)$ asymptotes at 2 for large $y$, so the combination $[1-f(y)/2]$ suppresses the $x$-dependent part of $a_{r}$ and $a_{v}$ at large $r$. Hence, we obtain

[TABLE]

as depicted in Fig. 23.

Broadly, there is more scatter in the trajectory parameters $a_{r}$, $b_{r}$, $c_{r}$, $a_{v}$, $b_{v}$, and $c_{v}$ of $r_{\mathrm{max}}$ and $v_{\mathrm{max}}$ than in the parameters $a$, $b$, and $c$ of the $J$ factor, when plotted against the system parameters $x$, $y$, and $z$ (cf. Figs. 10, 11, and 13). The source of this scatter is unclear, but it is likely that $r_{\mathrm{max}}$ and $v_{\mathrm{max}}$ are more sensitive than $J$ to additional effects beyond those accounted for by the parameters $x$, $y$, and $z$. Such heightened sensitivity is plausible for two reasons. First, $r_{\mathrm{max}}$ and $v_{\mathrm{max}}$ are more sensitive than $J$ to the density profile at large radii, which could depend on details of the tidal forces to which the inner profile is insensitive. Second, $r_{\mathrm{max}}$ and $v_{\mathrm{max}}$, being defined using the condition $\mathrm{d}(m(r)/r)/\mathrm{d}r=4\pi r\rho(r)-m(r)/r^{2}=0$, can be sensitive to fine details in the density profile $\rho(r)$. As an integrated quantity, the $J$ factor does not exhibit this sensitivity.