The effect of stellar encounters on the dark matter annihilation signal from prompt cusps Free

ABSTRACT

Prompt cusps are the densest quasi-equilibrium dark matter objects; one forms at the instant of collapse within every isolated peak of the initial cosmological density field. They have power-law density profiles, ρ ∝ r^−1.5 with central phase-space density set by the primordial velocity dispersion of the dark matter. At late times, they account for |$\sim 1~{{\ \rm per\ cent}}$| of the dark matter mass but for |$\gt 90~{{\ \rm per\ cent}}$| of its annihilation luminosity in all but the densest regions, where they are tidally disrupted. Here we demonstrate that individual stellar encounters rather than the mean galactic tide are the dominant disruptors of prompt cusps within galaxies. Their cumulative effect is fully (though stochastically) characterized by an impulsive shock strength |$B_* = 2\pi G\int \rho _*({\bf x}(t))\, \mathrm{d}t$| where ρ_*, the total mass density in stars, is integrated over a cusp’s entire post-formation trajectory. Stellar encounters and mean tides have only a small effect on the halo annihilation luminosity seen by distant observers, but this is not true for the Galactic halo because of the Sun’s position. For a 100 GeV WIMP, Earth-mass prompt cusps are predicted, and stellar encounters suppress their mean annihilation luminosity by a factor of two already at 20 kpc, so that their annihilation emission is predicted to appear almost uniform over the sky. The Galactic centre γ-ray excess is thus unaffected by cusps. If it is indeed dark matter annihilation radiation, then prompt cusps in the outer Galactic halo and beyond must account for 20–80 per cent of the observed isotropic γ-ray background in the 1–10 GeV range.

1 INTRODUCTION

The nature of dark matter is still unknown. One of the most popular candidates for dark matter searches is a weakly interacting massive particle (WIMP), which could be produced thermally in the early universe. If dark matter is a WIMP, it may have a significant self-interaction cross-section that allows dark matter to self-annihilate in regions where the dark matter density is sufficiently high (Arcadi et al. 2018; Roszkowski, Sessolo & Trojanowski 2018).

Detection of the secondary products of such self-annihilation – also known as indirect dark matter detection – is one of the most promising ways of learning more about the nature of dark matter. Excitingly, the Fermi large area telescope (Fermi LAT) has detected a galactic centre excess (GCE) in gamma-ray radiation in the spectral range of |$1-10 \,\mathrm{G}\mathrm{e\mathrm{V}}$| (Hooper & Goodenough 2011). Although there has been a long-ongoing debate about the precise properties of the signal (see Slatyer 2021, chapter 6, for a review), it seems so far that the morphology and the spectrum of the signal could be consistent with a dark matter self-annihilation signal from the central regions of the Milky Way’s dark matter halo (e.g. Di Mauro 2021). Additionally, recent high-resolution hydrodynamical simulations show that the Milky Way likely exhibits a halo with the right spatial structure to explain the GCE (Grand & White 2022). However, there are also other astrophysical processes that might explain the GCE so that it cannot yet be clearly attributed to dark matter. A few recent studies have suggested that the shape of the GCE may be better explained by templates that trace the stellar mass distribution of the Galactic bulge rather than a spherical dark matter halo (Bartels et al. 2018; Macias et al. 2018, 2019; Abazajian et al. 2020). However, a consensus conclusion about the morphology of the GCE remains elusive; for example, a study by Di Mauro (2021) assuming different background models finds consistency with radiation from a close to spherically symmetric dark matter halo.

Recently, Delos & White (2023) have highlighted a different aspect of the interpretation of possible indirect detection signals from the haloes of the Milky Way and of other galaxies. Recent advances in the modelling of the formation of the smallest nonlinear objects in a WIMP cosmology have led to a clearer understanding of their origin, structure, and late-time abundance, leading to a re-evaluation of the expected dark matter self-annihilation signal at late times.

The primordial density field is smooth below the dark matter free-streaming scale (e.g. |$\sim 1 \,\rm {pc}$| comoving for a typical WIMP) and so exhibits a large number of density peaks at this scale (∼10⁴–10⁵ per solar mass of dark matter). The first nonlinear structures begin to form around redshift 30 through monolithic collapse of these peaks. Their collapse history differs radically from that of traditional, more massive haloes of the kind characterized by Navarro, Frenk & White (1996) (hereafter: NFW). While NFW haloes assemble through hierarchical accretion and merging over time-scales, which are comparable to their age, a prompt cusp forms almost instantaneously at the moment of first collapse of a density peak, and it contains dark particles with orbital periods much shorter than the collapse time. As a result, the density profiles of prompt cusps also differ from the NFW form, following a steep r^−1.5 density profile between an outer boundary set by the curvature of the initial density peak and an inner core determined by the physical nature of the dark matter. To distinguish these dense ‘first’ objects from traditional NFW haloes, we follow Delos & White (2023) in referring to them as ‘prompt cusps’.

Delos & White (2022) argue that prompt cusps should still be extremely abundant substructures today. In regions where their number is not significantly reduced by subsequent evolution, every solar mass of dark matter should contain tens of thousands of them, and we should expect ∼10¹⁶ cusps associated with the dark matter halo of the Milky Way. Since they are denser in their centres than traditional NFW haloes, prompt cusps survive tidal effects better and produce a substantially larger dark matter annihilation signal. In fact, the annihilation signal from r^−1.5-cusps is logarithmically divergent with radius, limited by the inner and outer boundaries of the power-law profile, which are set by the primordial dark matter phase-space density (e.g. Macciõ et al. 2012) and by the initial peak extent, respectively. This raises the signal for indirect detection compared to most previous studies. Thus, Delos & White (2022) predict that these cusps should enhance the total annihilation luminosity of NFW haloes by factors ranging between 100 and 2500, depending on halo concentration, and additionally that they should significantly alter the morphology of the signal which, except in the densest regions, is proportional to the first power of the mean dark matter density rather than to its square, as usually assumed.

As Delos & White (2022) note, this ‘cusp’-component impacts a possible annihilation interpretation of the GCE in two significant ways. (1) If disruption of prompt cusps is ignored, their emission dominates that of the smooth halo component beyond about 5 degrees from the Galactic centre, resulting in a profile in disagreement with observation. Accounting for truncation and disruption by Galactic tides reduces cusp emission from the inner Galaxy but leaves it still dominant beyond 10 degrees, reducing but not eliminating the contradiction with observation. (2) If the GCE is nevertheless due to annihilation, emission from prompt cusps in the outer Galactic halo and external to the Milky Way must constitute a major contribution to the isotropic γ-ray background (IGRB). However, the IGRB appears to be almost completely attributable to emission from star-forming galaxies and active galactic nuclei (AGN; Blanco & Hooper 2019), leaving little space to add an additional component. The resulting upper limit on the prompt cusp contribution to the IGRB allows Delos & White (2022) to strengthen constraints on the self-annihilation cross-section and the thermal relic mass of a hypothetical WIMP dark matter particle – effectively excluding thermal relic WIMPs with masses |$M \le 10 \,\mathrm{T}\mathrm{e\mathrm{V}}$|⁠.¹

Delos & White (2022) explicitly considered only the effect of the smooth tidal field of the Milky Way on prompt cusps. However, stellar encounters can also be important, inducing strong impulsive shocks in dark matter substructures and potentially even disrupt them. Such encounters should be very frequent in the central region of the galaxy and therefore can affect both the predictions of the last paragraph.

The main goal of this paper is to evaluate quantitatively the effect of stellar encounters on the structure, survival, and predicted annihilation signal from prompt cusps. We will find that prediction (1) is seriously affected. After accounting for the effect of stellar encounters, the tension between the GCE emission profile and that predicted including cusp emission vanishes. For prediction (2), we will find that reduced emission from cusps in the inner Galactic halo leads to slightly lower IGRB predictions for the annihilation cross-section required to produce the GCE. These predictions remain in tension with claims that the IGRB is almost entirely due to other sources.

We note that there is already a large literature on the effect of stellar encounters on NFW subhaloes (Angus & Zhao 2007; Goerdt et al. 2007; Green & Goodwin 2007; Schneider, Krauss & Moore 2010; Delos 2019a; Kavanagh et al. 2021; Facchinetti, Stref & Lavalle 2022; Shen et al. 2022). Much of this is based on the incorrect assumption that systems disrupt when the total injected energy exceeds their initial binding energy, and hence needs to be read with care (Aguilar & White 1985; van den Bosch et al. 2018). A particularly clear and general treatment has been presented by Delos (2019a), and we will often refer to this article as representative of stellar encounters with NFW haloes. Such results cannot be applied to prompt cusps, since their power-law structure dramatically enhances their resilience to tidal effects in comparison to the inner regions of NFW profiles (see Stücker et al. 2023, for a discussion of different power-law profiles). The only previous study to consider the effect of stellar encounters on prompt cusps is Ishiyama, Makino & Ebisuzaki (2010), but unfortunately, this presented a very limited treatment and used incorrect assumptions when scaling with encounter strength, leading to the erroneous conclusion that cusps would never be disrupted by stellar encounters. We will discuss this in more detail below.

The structure of this paper is as follows. In Section 2, we briefly introduce both the theoretical basis for predicting the distribution of prompt cusp structural properties, and the physical formalism needed to describe the effects of their encounters with stars. We also present a novel, simple, and very general scheme for calculating the full distribution of impulsive stellar shocks experienced by cusps (or normal subhaloes) as they pass through a galaxy (e.g. through the Milky Way’s disc or bulge). In Section 3, we numerically integrate cusp orbits in a realistic Milky Way model, and we infer the parameters describing their shock histories. In Section 4, we use idealized N-body simulations to estimate the disruptive effects of stellar encounters, developing simple formulae that predict the structure and annihilation signal of cusps that have gone through arbitrarily many encounters. Additionally, we show how this can be supplemented to include the effects of stripping by the smooth Galactic tidal field. In Section 5, we present our main results, predictions for the profile of the prompt cusp contribution to the annihilation signal of the Galactic halo both as seen by a distant observer and as seen from the Earth, together with an assessment of how prompt cusps affect the form and relative amplitude of the GCE and the IGRB. In Section 6, we will discuss the implications of these findings for indirect dark matter searches.

We make almost all codes used in this study available through an online python repository² so that our methods can easily be used in future studies and the results of this paper can be reproduced independently.

2 THEORY

As described in the introduction, the effect of stellar encounters on NFW subhaloes has already been studied extensively, although in many cases using an incorrect criterion for subhalo disruption. However, there has been no realistic study of the effects for power-law cusps. Furthermore, even in the NFW case, most studies have not realistically modelled the full distribution of encounter parameters. Here we present the theoretical considerations necessary to treat full encounter histories in an accurate, general, but simple, way.

For this, we present in Section 2.1 the distribution of initial cusp properties as derived from the statistics of peaks in the initial Gaussian density field, in Section 2.2 how cusp structure can be modified by an inner core to account for the upper limit on phase-space density, in Section 2.3 the impulsive shocks induced by stellar encounters, in Section 2.4 a calculation of the total number of encounters expected on passing through a star distribution with arbitrary stellar mass and velocity distributions, and in Section 2.5 the statistical distribution of shock histories that follows from these considerations.

2.1 Cusps

Several numerical studies have found that peaks on the dark matter free-streaming scale collapse promptly to form dense cusps (Diemand, Moore & Stadel 2005; Ishiyama et al. 2010; Anderhalden & Diemand 2013; Ishiyama 2014; Polisensky & Ricotti 2015; Angulo et al. 2017; Ogiya & Hahn 2018; Delos et al. 2018a, b; Delos, Bruff & Erickcek 2019; Colombi 2021; White 2022; Delos & White 2023). These cusps have a density profile,

parametrized by a normalization A and an outer radius r_cusp, which limits the extent of the profile. Delos & White (2023) find that both parameters can be predicted well for a given cusp from the properties of the initial density peak from which it forms. Specifically,

where ρ₀ is the mean dark matter density of the universe today, a_col is the scalefactor when the peak first collapses, and R is the Lagrangian size of the initial density peak defined as |$R = \sqrt{|\delta /\nabla ^2 \delta |}$|⁠, where δ(x) is the linear overdensity field as a function of comoving position. The time of collapse can be estimated to sufficient accuracy by an ellipsoidal collapse model based on the triaxial structure of the initial peak. A more detailed description can be found in Delos & White (2022).

We follow the descriptions of Delos & White (2022) to sample a distribution of cusps. For this, we use a dark matter power spectrum generated by the Boltzmann code class (Blas, Lesgourgues & Tram 2011) up to a resolved wavenumber of |$k = 10^{4} \,\rm {\mathit {h}}\mathrm{M}^{-1} \rm {pc}$| at z = 30.6. Beyond that scale, we use the analytic prescriptions of Hu & Sugiyama (1996) but normalized so that it matches the class spectrum at |$k = 10^{4} \,\rm {\mathit {h}}\mathrm{M}^{-1} \rm {pc}$|⁠. We multiply this spectrum using the exponential power spectrum cutoff description of Bertschinger (2006) for a WIMP with mass |$m_{\rm {WIMP}} = 100 \,\mathrm{G}\mathrm{e\mathrm{V}}$| and decoupling temperature |$T_{\rm {d}} = 30 \,\mathrm{M}\mathrm{e\mathrm{V}}$| (corresponding to a decoupling scale-factor of a_d = 5.33 × 10⁻¹²). We use the free-streaming length |$k_{\rm {FS}} = 1.06 \,\rm {pc}^{-1}$| as calculated by Delos & White (2022). The distribution of initial density peaks can be sampled analytically given the initial power spectrum as described by Bardeen et al. (1986). We have created our own implementation of the sampling of peaks, which is published in the code repository of this paper. However, we tested it against the peak sampling implementation that was published by Delos et al. (2019). We map the distribution of peaks onto a distribution of cusps by using equations (2–3) with the ellipsoidal collapse correction for a_col as explained by Delos & White (2022) based on the approximation from Sheth, Mo & Tormen (2001).

We show the resulting distribution of cusps in Fig. 1 (dashed contours). However, more relevant than the number distribution of cusps is the annihilation weighted distribution. Under the assumption of a velocity-independent cross-section, the annihilation rate of a cusp should be proportional to

where we have neglected any contributions from r < r_core and r > r_cusp, where r_core is the core radius enforced by the phase-space density constraint (Delos & White 2023). We will explain in Section 2.2 how to calculate and treat the core radius more precisely.

The weighted distributions are shown as solid contours in Fig. 1. We can see that for this specific dark matter particle, the typical cusp relevant for the annihilation rate has |$A \sim 8 \times 10^{-4} \,\rm {{\rm M}_{\odot }}\rm {pc}^{-3/2}$| and |$r_{\rm {cusp}} \sim 5 \times 10^{-3} \,\rm {pc} \approx 10^{3} \,\rm {AU}$| (in physical units). It collapses at redshift z ∼ 30.

2.2 Phase-space cores

The fine-grained phase-space density of dark matter is constant as a function of time (Liouville’s theorem). The coarse-grained phase-space density – defined as the average over some finite phase-space volume – can therefore never exceed the fine-grained value established at dark matter freeze-out (Tremaine & Gunn 1979). Let us denote the maximum of this fine-grained phase-space density by f_max. The value of f_max is set in the early universe and depends strongly on the type of dark matter considered. For the WIMP model considered above, for example, it is |$f_{\rm {max}} = 9.98 \cdot 10^{13} \frac{\,\rm {{\rm M}_{\odot }}}{\,\mathrm{k}\mathrm{m}^{3} \mathrm{s}^{-3}\rm {pc}^{3}}$|⁠, whereas for thermal relic warm dark matter with |$m_X = 3.5 \rm {keV}$|⁠, it is |$f_{\rm {max}} = 6.75 \cdot 10^{-2} \frac{\,\rm {{\rm M}_{\odot }}}{\,\mathrm{k}\mathrm{m}^{3} \mathrm{s}^{-3}\rm {pc}^{3}}$| (see e.g. Delos & White 2023)³

The power-law profile of equation (1) corresponds, for an isotropic velocity distribution, to a phase-space distribution function,

(see e.g. Stücker et al. 2023). This description must break down at energies where f(E) > f_max. To estimate the radial scale where this happens we can consider at each radius the largest reachable phase-space density – given by f(ϕ(r)) where ϕ(r) is the potential. Inserting this into equation (5) and inverting for r, we find

Thus, in the above WIMP model, typical cusps relevant for the annihilation calculation have thermal core size, |$r_{\rm {core}} \sim 10^{-5} \,\rm {pc} \sim 2 \,\rm {au}$|⁠.

The detailed shape of the density profile near the core radius is unclear. It would be desirable to run simulations with the actual primordial velocity distribution of a WIMP, so that phase-space cores are created self-consistently, and then to measure their density profile. Such simulations would be computationally demanding and have not yet been performed, but they are not far beyond current capabilities. (Consider that typically r_cusp ∼ 500r_core so that simulations that resolve the cusp profile well are typically only an order of magnitude in linear scale from the core radius.)

To obtain a profile that has the desired maximum phase-space density and connects smoothly to the appropriate power-law at larger radii, we make the following Ansatz for the phase-space density,

where E_core is the core energy scale where the phase-space density f_max would be reached in the fiducial power-law profile:

so that together with equation (6) the profile is fully specified through A and f_max. Note that for f_max → ∞, the pure power-law profile is recovered.

We can integrate over the velocity components of the phase-space density to find the density,

as a function of the potential, ϕ, which is normalized so that ϕ(0) = 0 and ϕ(r → ∞) → ∞. Combined with Poisson’s equation, this forms a nonlinear second-order differential equation for the potential,

We do not know how to solve this differential equation analytically and therefore solve it through numerical integration starting at r = 0 using |$\phi (0) = \partial _r \phi (0) = 0$|⁠. As a result, we find ρ(r) and ϕ(r) and we display these and the phase-space density f(ϕ(r)) in Fig. 2.

As expected, the density profile reaches a well-defined maximum at r ≲ r_core, given by

A good approximation to the density profile is given by

which we also show as an orange line in Fig. 2. This deviates at most by 20 per cent from the actual density around r ∼ 3r_c where the latter is slightly enhanced with respect to the fiducial power-law. This enhancement is similar to the well-known behaviour of the isothermal sphere where the cored solution rises above the asymptotic r⁻² power-law solution at radii comparable to the core radius before aymptoting to constant density at smaller radii. (e.g. Binney & Tremaine 2008).

We note that while there are no simulations that have a phase-space density constraint consistent with the expected free-streaming scale, there have been simulations by Macciõ et al. (2012) with artificially large initial velocity dispersion and hence a lowered upper limit on phase-space density. Although it is difficult to compare these directly with our profile, it seems that the final cores do saturate the phase-space bound and that they transition relatively quickly to the asymptotic power-law behaviour – both consistent with the profile that we propose here.

We find that the annihilation radiation from the cored profile inside radius r can be approximated for r > 10r_core by

This is marginally larger than the annihilation radiation one obtains when assuming the power-law profile down to r_core [compare equation (4)] and a uniform density at smaller radii, in which case the 0.531 gets replaced by 0.333. This is due to the slight enhancement of the profile around 3r_core.

When needed, we will use the numerical cored profile throughout this study.

2.3 Impulsive encounters

We consider a star with mass M_* passing a cusp on a linear orbit at constant relative velocity v and minimal distance b. We can approximate the tidal forces acting on cusp particles through a multipole expansion of the potential up to second order,

where ϕ_s denotes the self-potential of the cusp, |$\boldsymbol{x}_0$| is the location of centre of the cusp, and |$\bf{T}_0(t)$| is the tidal field of the star evaluated at |$\boldsymbol{x}_0$|⁠. Here, we have neglected zeroth- and first-order terms, since they do not affect the internal dynamics of the cusp (e.g. Stücker, Angulo & Busch 2021). This ‘distant-tide’ approximation is valid for particles that are close enough to the centre, |$||\boldsymbol{x}-\boldsymbol{x}_0|| \ll b$|⁠. We will see that for encounters with stars with masses of order |$\,\rm {{\rm M}_{\odot }}$| the distant tide approximation is excellent for all particles that remain bound to the cusp (which typically have small |$\Delta \boldsymbol{x}$|⁠) and is only violated for particles that are kicked so strongly that they will leave the system. Thus, it is safe to adopt this approximation in all our calculations.

Further, we can assume the impulsive approximation that considers the limit that particles move very little within the cusp during the time of the encounter. In this case, the total change in the velocity of a particle due to the encounter can be approximated by

where we have defined the shock tensor K.

It is easy to see that the impulsive approximation is an excellent approximation here. The dynamical time-scale at radius r of our cusp is given by

where M(< r) is the enclosed mass profile. The internal dynamical time-scale is shortest at the core-radius, where it is of order |$2 \times 10^{4} \,\rm {yr}$| for the strongest cusps (which dominate the annihilation distribution). The impact parameter of the weakest relevant encounters is of order |$b \sim 10^{4} \,\rm {au}$| with v ∼ 200 kms⁻¹, which gives an encounter time-scale of order |$t_{\rm {enc}} = b/v \sim 200 \,\rm {yr}$| which is two orders of magnitude smaller than the dynamical time-scale of the quickest particles in the cusp. Stronger encounters (which are more relevant) will happen on even shorter time-scales and most particles orbit on longer time-scales so that the ratio should be even larger in practice, and we can safely assume the impulsive limit for all of our calculations.

If we assume that the star is a point mass with tidal field

and we assume its trajectory (without loss of generality) to be along the y-direction with the closest encounter at the coordinate (b, 0, 0), then the shock tensor is given by

(e.g. Aguilar & White 1985) where we have defined the shock parameter B. We note that B has dimensions of the inverse of time, but to simplify intuitive understanding, we will typically state it in units of |$\,\mathrm{k}\mathrm{m}\mathrm{s}^{-1} \rm {pc}^{-1}$|⁠.

It is clear that the only relevant parameter for describing the effect of a distant and impulsive stellar tidal shock on a prompt cusp is the tidal shock parameter B. The individual values of b, M_*, and v matter only insofar that they determine the value of B.

We can get a feeling for what range of B values will be relevant for typical cusps by considering the values,

B_cusp indicates the strength of a tidal shock that is needed to induce a velocity change at the radius r_cusp as large as the circular velocity and B_core is the analogue quantity at the core radius. We show the distributions of these two parameters for the fiducial |$100 \,\mathrm{G}\mathrm{e\mathrm{V}}$| WIMP in Fig. 3. We can expect that tidal shocks with B ≳ B_cusp will significantly alter the profile of the cusp. Tidal shocks with B ≳ B_core may possibly lead to complete disruption. Shocks with B ≪ B_cusp will leave the cusp largely unaffected. Typical cusps that are relevant for the annihilation rate have values of B_cusp and B_core of order |$0.3 \,\mathrm{k}\mathrm{m}\mathrm{s}^{-1} \rm {pc}^{-1}$| and |$30 \,\mathrm{k}\mathrm{m}\mathrm{s}^{-1} \rm {pc}^{-1}$|⁠, respectively. The impact parameters that are needed to reach such shock parameters for |$M_* = 1 \,\rm {{\rm M}_{\odot }}$| and v = 200 kms⁻¹ are |$1 \times 10^{-2} \,\rm {pc} \approx 2000 \,\rm {au}$| and |$1 \times 10^{-3} \,\rm {pc} \approx 200 \,\rm {au}$| respectively. We will see that tidal shocks of this order are not only possible but also quite likely. It is therefore important to make precise quantitative calculations to evaluate the number of such encounters and the effect they have on the aggregate annihilation luminosity from cusps.

Finally, it is useful to introduce the characteristic spatial scale associated with the action of shocks on a cusp,

which is the radius where the change in velocity induced by the shock is of order the circular velocity. We will see in Section 4 that the majority of particles beyond r_B will leave the system. The distant tide approximation is a good approximation if min(r_B, r_cusp) ≪ b so that it holds for all particles that remain bound. We find that typical encounter scenarios with stars for prompt cusps have min(r_B, r_cusp) ≪ 0.1b for any impact parameter b, so that the distant tide approximation is always valid. While we do not focus on NFW haloes in this study, we note that many of our subsequent derivations are of interest for studies of NFW substructures as well. For these, we estimate that |$r_{B, \rm {NFW}} \lt 0.1 b$| (defined through the circular velocity criterion for an NFW profile) holds for typical closest impact parameters with b ≲ 1000 au for haloes with virial masses |$m_{200\rm {c}} \lesssim 1 M_\odot$| (and concentration c ∼ 30). Many of our results below can therefore also be applied to NFW subhaloes with masses below 1.0 M_⊙, but additional care must be taken if larger masses are considered since the distant tide approximation may then fail.

The goal of the remaining parts of this section is to determine the distribution of shock parameters B for a cusp that is moving through an arbitrary stellar distribution (e.g. a component of the Milky Way). We note that this calculation is both simpler and more accurate than previous calculations in the literature, which have focused on estimating the distribution of the impact parameter b.

2.4 The expected number of encounters

We want to estimate the expected number of encounters with a tidal shock parameter that is greater than B, N(> B), for a cusp that is orbiting through an arbitrary stellar distribution. The stellar distribution can be described by a mass-dependent phase-space (number) density

which is normalized so that an integral over a phase-space volume and over a stellar mass interval gives the number of stars within that volume and mass range. An integral over phase space alone gives the stellar mass function, which is allowed to vary spatially.

If our cusp was passing through a uniform medium without velocity dispersion and with number density n_* then the expected number of encounters in a time-interval dt with impact parameters in the range [b, b + db] is given by

where v is the relative velocity with respect to the medium. Equation (25) arises by considering stars when they get closest to our cusp, which is when they pass through the plane orthogonal to the velocity vector. The relevant velocity here is the relative velocity between the stars and the cusp; higher velocities lead to more frequent encounters.

For arbitrary phase-space distributions, we have to consider each velocity and mass bin individually. The number of encounters within time interval dt, impact parameter range (b, b + db) and encounter velocity bin d³v with stars that have masses in the range |$(M, M + \, \mathrm{d} M)$| is given by

where |$\boldsymbol{v}$| is the encounter velocity and v_rel is the relative velocity between our cusp and the zero point of the stellar phase-space distribution. Here, we have assumed that the phase-space distribution function does not vary significantly over the distance b. This is an excellent approximation, since typical encounters of interest will have |$b \ll 1 \,\rm {pc}$|⁠, whereas the stellar distribution varies on much larger scales (⁠|$\gg 100 \,\rm {pc}$|⁠).

Now, it would be straightforward to estimate, for example, the total number of encounters with impact parameters smaller than some given b by integrating (26) over the corresponding variables. In general, this would turn out to depend both on the phase-space distribution of the stars and the stellar mass function.

However, as discussed in the previous subsection, we are actually interested in the distribution of shock parameters B. This can be evaluated by integrating (26) under the constraint that B = GM/vb², which leads us to

where δ_D is the Dirac delta function. We evaluate the integral, integrating in b first:

where we have used the substitution |$B^{\prime } = \frac{2MG}{v b^2}$| to evaluate the integral.

Sorting and evaluating individual terms gives us

where we have used the fact that the mass weighted integral over the phase-space number density gives the stellar mass density ρ_*. Further we have defined the time integral of the stellar density along the cusp trajectory χ_*. It is further convenient to define the characteristic shock parameter,

so that

Therefore, we expect on average one encounter with B > B_*.

It is worth noting that this result is surprisingly independent of the phase-space distribution function of stars and the stellar mass function but depends only on the stellar mass density. The mass function is irrelevant, because at a fixed mass density, a reduction in mass leads to an increase in number density, thus making close encounters more likely to the same degree that it decreases the strength of such encounters. A similar coincidence holds for the velocity dependence. The number of stars that are encountered increases with the velocity, while at the same time, the weight of each encounter decreases with the velocity to such a degree that the two effects cancel exactly. We call these effects the encounter conspiracy.

It is clear that these two simplifications could, in principle, break down at some scale. Mass function independence breaks down if the distant tide approximation fails – if we make our perturbers less massive, the encounters get closer for a given value of B. For very small masses e.g. |$10^{-6} \,\rm {{\rm M}_{\odot }}$| the closest approach needed for a significant perturbation would be of order 1 au, smaller than core radius of a cusp; the distant tide approximation would then certainly fail. The integral over the mass function in equation (30) should thus have a lower limit, in principle. The independence of velocity breaks down for very small encounter velocities. If an encounter takes longer than the orbital times within the cusp, then the cusp will react adiabatically, with no long-term changes in energy except for particles that leave the system in the adiabatic limit. Thus our approximations fail for stars that are moving almost at the same velocity as the cusp. However, neither of these problems has a significant effect in practice, since almost all stellar mass is in objects of mass within an order of magnitude or so of |$1.0~\,\rm {{\rm M}_{\odot }}$| and very few encounter velocities are smaller than, say, 10 kms⁻¹.

We note that the effects of encounters with other massive objects, such as planets or other prompt cusps would be overestimated if the calculation of this section were applied. However, the mass density in planets is so much lower than that in stars that planets are quite irrelevant in this context. The mass density in prompt cusps is non-negligible at large halocentric radii, but when their extended profile is taken into account, we find that the strongest possible shocks are far below |$B \ll 10^{-4} \,\mathrm{k}\mathrm{m}\mathrm{s}^{-1} \rm {pc}^{-1}$| so that cusps cannot significantly shock other cusps (cf. Fig. 3). Therefore, stars pose the only significant contributor to the distribution of encounter shocks.

2.5 Shock histories

We assume that all aspects of the problem follow Poisson statistics – for example that stars are drawn through a Poisson process from the continuous phase-space distribution and that stellar masses are drawn through a Poisson process from the stellar mass function. Then also the shock parameter distribution has to follow Poisson statistics. That means the probability of having exactly k encounters with shock strength bigger than B is given by

In particular, the probability of having at least one encounter with shock strength bigger than B is given by

The probability density function of the strongest encounter is therefore

It is straightforward to derive the corresponding functions for the second, third, etc. strongest shock. However, the distribution of e.g. the strongest and the second strongest shock are not independent and parametrizing the joint distribution is rather cumbersome. When considering a large number of encounters, it is more convenient to draw actual realizations. This can easily be done by mapping B onto another random variable that follows a uniform distribution x ≔ B_*/B,

We can then create a realization of a history of shocks with B > B_min, by sampling k uniformly distributed random variables x_i on the interval [0, B_*/B_min] and then transforming them as B = B_*x⁻¹. The number k must itself be drawn from a Poisson distribution with mean 〈k〉 = B_*/B_min. We note that there will be infinitely many encounters as B_min → 0 so that it is in general not possible to sample the full distribution, but only its truncated form. However, this is sufficient, since encounters with very small shock parameters become irrelevant. In particular, we will show in Section 4 that the joint effect of k encounters is more or less equivalent to a single encounter with an effective shock strength,

with p = 1.2.

We sample a large number of shock histories and show the strongest three encounters together with the distribution of B_eff in Fig. 4. The distribution of B_eff is already reasonably well approximated when only considering shocks with B > 0.1B_* i.e. approximately the 10 strongest shocks. It is almost fully converged when considering all shocks with B > 10⁻³B_* i.e. the 1000 strongest shocks. Typically the effective shock parameter is not too much larger than the strongest shock. Its median lies at 4.5B_*, whereas the median of the strongest shock lies at 1.46B_*. However, the low-end tail is shifted upwards quite a bit so that there are almost no cases with B_eff < B_*. The high-end tail is virtually identical to the distribution of the strongest shock.

3 SHOCK DISTRIBUTION FOR ORBITS IN THE MILKY WAY

In this section, we evaluate numerically the quantities that are relevant for describing the full distribution of shock parameters for a cusp that is orbiting in the Milky Way. We showed in the last section that, for a given orbit, the time integral of the stellar mass density along the cusp’s trajectory,

is sufficient to parametrize the full distribution of encounter shocks. Our aim is thus to infer realistic estimates of χ_* and of the corresponding characteristic shock parameter B_* = 2πGχ*.

We assume that the orbital distribution of cusps follows that of dark matter particles within the Milky Way’s halo.

3.1 Milky Way potential

For the baryonic components of the Milky Way, we assume the prescriptions used in the Phat ELVIS simulations (Kelley et al. 2019) at the specific time z = 0. The observational parameters underlying these simulations were taken from McMillan (2017) and Bland-Hawthorn & Gerhard (2016).

The stellar disc and the gas disc are each modelled through a superposition of three Miyamoto & Nagai (1975; hereafter MN) potentials as described by Smith et al. (2015). The parameters of the MN potentials have been tuned to recreate the mass distribution of an exponential disc up to 1 per cent accuracy. For the stellar disc, we use a total mass of |$M_{\mathrm{d}} = 4.1 \times 10^{10} \,\rm {{\rm M}_{\odot }}$|⁠, a scale radius |$R_{\mathrm{d}} = 2.5 \,\mathrm{k}\rm {pc}$|⁠, and a height parameter of |$z_{\mathrm{d}} = 0.35 \,\mathrm{k}\rm {pc}$|⁠. For the gas disc, we use |$M_{\mathrm{d}} = 1.9 \times 10^{10} \,\rm {{\rm M}_{\odot }}$|⁠, |$R_{\mathrm{d}} = 7.0 \,\mathrm{k}\rm {pc}$|⁠, and |$z_{\rm {d}} = 0.08 \,\mathrm{k}\rm {pc}$|⁠.

We approximate the stellar bulge through a Hernquist (1990) distribution with mass |$M_{\rm {B}} = 9 \times 10^{9} \,\rm {{\rm M}_{\odot }}$| and with scale length |$r_{\rm {a}} = 500 \,\rm {pc}$|⁠.

For the Milky Way’s dark matter halo, we assume that in the absence of baryonic components it would correspond to an NFW halo with concentration c = 8.7 and mass |$m_{200c} = 10^{12} \,\rm {{\rm M}_{\odot }}$|⁠. However, the baryons increase the dark matter density in the inner regions. We model this contraction using the semi-analytic approach presented in Cautun et al. (2020). Previously, we used the analytic procedure of Sellwood & McGaugh (2005), but this produces a slightly denser result in the central regions. We prefer the Cautun et al. (2020) approach, since it has been tested in detail both against state-of-the-art simulations and against observations of the Milky Way.

We show the result of the contraction of the halo together with profiles of all the different components in Fig. 5. We note that the contracted halo (in purple) is significantly denser in the centre than the uncontracted version (brown dashed line). This contraction is only moderately relevant for getting the correct potential structure of the Milky Way. However, it is very relevant for sampling self-consistent orbits of the dark matter distribution. To sample orbits, we use the (Eddington 1916) inversion method on the density profile of the contracted halo but inside the joint potential of all components. We have checked that if we sample particles from the contracted halo profile and integrate their orbits in the joint potential, the density profile is stable and evolves very little at later times.

3.2 Orbits

We create 10⁵ test particles from the contracted halo up to |$400 \,\mathrm{k}\rm {pc}$|⁠. We use an adaptive sampling method so that their initial number density is proportional to 1/r times that of the halo. This way we have a good sampling, both at small radii r ≲ 1 kpc and close to the virial radius r ≳ 200 kpc. Throughout this paper, when showing distributions, we always correct this non-uniform sampling using the appropriate weights. We integrate particle orbits in the full 3D Milky Way potential (including bulge, stellar disc, gas disc, and contracted halo) using 10⁵ constant time-steps over a period of |$10 \,\mathrm{G}\rm {yr}$|⁠. Along each orbit, we evaluate χ_* according to equation (39) using the stellar density inferred from the Laplacian of the potential of the stellar components (stellar disc and bulge) ρ_* = Δϕ_*/(4πG). Note that this can create negative densities in a (very) few locations, since the MN potential of the disc can create negative densities. We therefore clip χ_* after the integration at a minimum of |$10^{-8} {\rm M}_{\odot }\rm/ {pc}^{2} \mathrm{/}(\mathrm{k}\mathrm{m}\mathrm/{s}).$|⁠. This threshold has no impact on any of our results. The actual minimum should be set by the diffuse stellar halo (and partially by encounters with dark substructures). However, none of these aspects matter since, as we will see in Section 4.6, the effect of the smooth tidal field is much larger than this at radii where stellar encounters are infrequent.

We find the distributions of χ_* and B_* shown in Fig. 6. Here and in later plots, we assign each particle 1000 times with its final value of χ_* at different radii chosen uniformly in time over its full orbit history. To help with an intuitive understanding of the χ_* distribution, we have multiplied the χ_* axis by a value of 200 kms⁻¹ so that the left y-axis would correspond to the total stellar column density if the cusp always encountered stars at 200 kms⁻¹. Fig. 6 shows that a typical cusps orbiting around the solar radius (⁠|$8 \,\mathrm{k}\rm {pc}$|⁠) encounters a stellar column density of about |$10^{4} \,\rm {{\rm M}_{\odot }}\rm {pc}^{-2}$|⁠. These numbers can easily be understood, as orbits at these radii will have passed through the Galactic disc about 10² times and each disc passage adds a column density of order |$10^{2} \,\rm {{\rm M}_{\odot }}\rm {pc}^{-2}$|⁠.

Further, we show as an alternative y-axis in Fig. 6 the distribution of B_*. Recall that B_* indicates the value of B such that we expect on average one encounter with B > B_*. Therefore, we expect typically one encounter with |$B \gtrsim 1 \,\mathrm{k}\mathrm{m}\mathrm{s}^{-1} \rm {pc}^{-1}$| for cusps orbiting around the solar radius.

For each orbit, we sample the strongest encounters corresponding to the probability density in equation (36). In this way, we find the distribution of strongest shocks shown in Fig. 7. We also show the corresponding encounter distance for encounters with mass |$1 \,\rm {{\rm M}_{\odot }}$| and velocity 200 kms⁻¹. This figure is not used in our final results, but it is useful to gain intuition for the relevant quantities.

We conclude that inside the central |$10 \,\mathrm{k}\rm {pc}$| virtually all cusps will have had at least one stellar encounter that might significantly decrease the annihilation luminosity and possibly even disrupt them. For a reliable evaluation, we have to conduct numerical experiments that evaluate how strongly cusps react to shocks of a given amplitude.

4 THE EFFECT OF IMPULSIVE ENCOUNTERS

To evaluate the effects of impulsive encounters on prompt cusps, we ran several sets of N-body simulations addressing scenarios of increasing complexity. The first set considers pure power-law profiles experiencing a single encounter. The second set considers cored profiles also with a single encounter. The third set considers both cored and power-law profiles experiencing multiple stellar encounters. For each case, we develop simple descriptions for the annihilation luminosity expected from the remnants.

Finally, at the end of this section, we briefly discuss how the effect of the smooth tidal field can be incorporated simply in this formalism.

4.1 Simulation setup

We consider simulations starting from two different initial profiles. Power-law simulations start from a Poisson realization of a pure r^−3/2 density profile with the phase-space distribution from equation (5). In this case, there is no relevant length or mass scale so that the results can be rescaled to any desired tidal shock strength. For cored simulations, we use the phase-space distribution from equation (8) and the density profile that is self-consistently created from this distribution, as described in Section 2.2. For these, we set r_core = 1 so that length scales are given in units of the core radius.

To enhance the dynamic range that can be resolved in these simulations and to minimize the effects of the (numerically required) outer truncation radius, we use a similar non-uniform mass sampling strategy (Delos 2019b). Specifically, we arrange that the same number of particles N_split have pericentres in each of the radial ranges (0, 1), (1, 10), (10, 100), and (1000, 10000), which then have particle masses increasing by about a factor 30 between intervals. We explain this in more detail in Appendix A1, where we also present some numerical stability tests. An implementation of this method is available in the code repository of this article. We do not find any artefacts due to the non-uniform mass sampling, most likely because we use quite high resolution and because our pericentre criterion minimizes the intrusion of higher mass particles into higher-resolution regions.

For simulations with a single encounter (presented in Sections 4.2 and 4.3), we apply at time t = 0 a single kick according to equation (16) using the tensor from equation (19), with shock amplitude B. The shock introduces the characteristic spatial scale,

which is the radius where the tidal shock causes a (maximal) velocity change equal to the circular velocity in a pure power-law profile. By definition, r_B equals r_cusp and r_core for shocks with strength B_cusp and B_core, respectively, so that realistic shocks can easily produce r_B values that lie at any point of the cusp (compare Fig. 7). Further, the kick introduces a characteristic time-scale – corresponding to the dynamical time at this radius:

After the initial shock, we evolve the simulation for several dynamical times t_B until the profile no longer evolves at the radii of interest (e.g. ∼10t_B around r_B). Note that for typical shocks with |$B \sim 1 \,\mathrm{k}\mathrm{m}\mathrm{s}^{-1} \rm {pc}^{-1}$|⁠, we have |$t_B \sim 1 \,\mathrm{M}\rm {yr}$|⁠.

4.2 Truncation of pure power-law profiles

We ran two high-resolution simulations of shocked power-law profiles that each have N_split = 2²² particles per radial interval, so that they have in total about 2 × 10⁷ particles. These two simulations used different amplitudes for the shock parameter, varying by a factor of 4, leading to differing shock radii, r_B = 20 and r_B = 128. We remind the reader that, in principle, the result of a shocked power-law can be rescaled arbitrarily, so these two simulations differ only in how the physical scale r_B compares to resolution parameters.

We define the transfer function,

where ρ(r) and ρ₀(r) are the final and initial density profiles, respectively, and we present transfer functions for our simulations in Fig. 8 as functions of r/r_B. We show only the radial range that has reliably converged to a final, stable post-shock profile, and we use different output times to probe different parts of the transfer function. In Appendix A2, we provide convergence tests to determine these scales. At the small end, this is the largest radius where two-body relaxation is irrelevant, while at the large end, we limit to scales where at least 10 dynamical time-scales have passed. As additional evidence that the simulations are converged, we also show as grey lines in Fig. 8 reference simulations evolved to the same times but with no shock. Clearly, these are stable and show no discernable numerical evolution.

The simulations of different shock strengths line up very well if radii are scaled by r_B. In turn, since B ∝ b⁻², r_B scales with impact parameter as b^8/3. This is very different than the scaling of b^8/11 that Ishiyama et al. (2010) assumed (without further explanation) to extrapolate their results and to argue that stellar encounters should be irrelevant for the survival of cusps. Such a scaling is clearly incorrect and is also inconsistent with the profiles of the simulations that Ishiyama et al. (2010) presented themselves. These simulations, in fact, seem consistent with a b^8/3 scaling and agree qualitatively with what we find here.

Inspired by Delos (2019b), we fit the simulated transfer functions jointly using a function of the form,

where we find α = 1.256 and β = 0.639 as the best-fitting parameters. Interestingly, the value of β that we find here for our power-law profile is not too far from the one that Delos (2019b) found (β = 0.78) for an NFW initial structure. The small difference presumably reflects the different central slopes, −1.5 and −1.

We note that our measured profiles disagree slightly with our fit at larger radii r ≳ 2r_B. However, this is quite irrelevant for the calculation of the annihilation luminosity where those contributions are downweighted by a factor T²(r). For our fitted function, we find that the annihilation rate is given by

where |$\operatorname{Ei}$| is the exponential integral and r_min is a lower limit of integration that is necessary to obtain a finite result. The approximation that is used in the last line, r_cusp → ∞, is already very accurate for r_cusp ≳ r_B. Therefore, if r_B ≤ r_cusp, it is fine to neglect the initial boundary of the system, since the truncation through the shock sets the actual boundary. Further, it is insightful to consider the limit r_min ≪ r_B in which case

where γ ≈ 0.577 is the Euler–Mascheroni constant. This approximation holds (at 5 per cent accuracy or better) for radii r_min < 0.1r_B. One way to interpret this result is that the annihilation rate of the tidally truncated power-law corresponds to the annihilation rate of a power-law profile that is sharply truncated at 10 per cent of r_B [compare equation (4)]. This result is easily understood, since this is approximately the radius where T² reaches 50 per cent. We expect this approximation to be inaccurate if r_core ≳ 0.1r_B, since then shock effects are significant at the core radius where they may be amplified. We will consider such cases in the next subsection and derive a general formula for the annihilation rate of shocked cusps. However, we can make a rough estimate of the cored annihilation rate, by assuming the power-law + transfer profile up to the core radius and downweighting the core contribution by a factor of T²(r_core)

4.3 Truncation of cored profiles

We also ran several simulations with cored initial conditions and with different shock parameters. These are designed to probe the scales where the shock starts affecting the core. Each uses N_split = 2²⁰ and therefore in total about 5 × 10⁶ particles. We evaluate each simulation at a time where the final profile has reached a stable result. This requires a longer integration time for simulations where the central density is reduced significantly. We give the simulation parameters, durations, and final annihilation luminosities in Table 1.

We show the transfer functions of these simulations in Fig. 9. We can see that as r_core/r_B → 0, the pure power-law behaviour is recovered; at large radii, they have the same transfer function as the power-law case. At small radii, cored profiles are suppressed additionally and can even disrupt completely for sufficiently strong shocks.

The simulation with r_core/r_B = 0.666 is the first case that exhibits complete disruption – this means that after 50t_B, no stable remnant is left, and densities keep decreasing everywhere. We checked that this case still disrupts if a four times higher particle number and a smaller softening are used and that simulations with even larger shocks disrupt also. In addition, we a provide a video of one non-disrupting case and one disrupting case online.⁴ It seems clear that although it is impossible to disrupt centrally divergent power-law profiles, it is indeed possible to disrupt cored profiles. This agrees with theoretical (Amorisco 2021; Stücker et al. 2023) and numerical investigations (van den Bosch et al. 2018; Errani & Peñarrubia 2020; Errani et al. 2022) in the context of NFW haloes.

For the cored profiles, we do not attempt to fit a functional form to the transfer functions but instead directly calculate the annihilation radiation from the spherically averaged density profiles. This works very well numerically, because the part of the profile that is responsible for the bulk of the annihilation radiation (most significantly around 3r_core) is resolved very well.⁵ We describe the binning and integration procedures in Appendix A3 and show that all obtained annihilation luminosities are accurate at the 3  per cent level – except for the last non-disrupted case, r_core/r_B = 0.543, where the relative error is larger but still less than 40 per cent.

We list the annihilation rates we obtain in Table 1 and plot them in Fig. 10. In this figure, we also show the two power-law simulations with annihilation rates estimated by individually fitting their transfer functions and (arbitrarily) assuming a core radius, r_core = 0.05, for the annihilation estimate from equation (46). However, these two simulations could be arbitrarily rescaled along the black line.

We are able to fit the post-shock annihilation rates of our cored simulations using the function,

with free parameters a and b. By construction, this function approaches the behaviour of the power-law truncation fit of equation (46) for B ≪ B_core. However, it has two degrees of freedom a and b, which can modify the behaviour for large values of B. We find that a = 0.708 and b = 5.98 gives a reasonable fit to all annihilation rates at the |$10~{{\ \rm per\ cent}}$| level – both in the power-law regime B ≪ B_core and when the shock hits the core B ≳ 0.1B_core. The function reaches zero at B_dis ≈ 0.65B_core. We assume that all cusps that had an encounter with B > B_dis are disrupted and have zero contribution to the annihilation rate.

4.4 The effect of multiple encounters

As a final step, we need to understand what happens to the remnant if it is exposed to multiple shocks. To test this, we set up a large number of simulations (for both power-law and cored cases) with a variety of different shock histories. For each of these simulations, we apply a shock approximately every 10 dynamical time-scales t_B, and we use up to five shocks. We always compute the annihilation rate when 10t_B have passed since the last shock. We note that the value of t_B is not very clearly defined in the case with different shock amplitudes; we therefore choose a reference value B_ref to define a t_B that seems appropriate for each individual history. We have checked that this shock interval is large enough to ensure that the results would not change by slightly decreasing it or by increasing it further. We list all shock histories in Section 2. We have created several manually chosen shock histories (e.g. equal shock strengths, descending shocks, and ascending shocks), and we have also sampled a few shock histories from the full distribution of shock histories with B_* = B_ref, but only keeping the five strongest shocks – see Section 2.5. We sorted some of these histories accidentally, but we also created an additional set of four cored simulations with unsorted shock histories. However, neither the annihilation rate nor the transfer functions depend much on the order in which shocks happen.

We present transfer functions for power-law simulations with multiple encounters in Appendix A4. Most importantly, we note that cases with multiple encounters are still well approximated by the transfer function of equation (43), but with different cut-off radii, and we note that for the final transfer function, the order of shocks does not matter, only their amplitudes. These results are both consistent with findings in a previous study of NFW haloes (Delos 2019b).

We show the annihilation rates of cusps that have gone through multiple shocks in Fig. 11. The data points in this figure are obtained as follows. For power-law cases, we compute the annihilation rate by first fitting a transfer function to each encounter individually using equation (43) and then calculating the annihilation rate according to equation (46). Here we rescale the power-law results to two different core radii r_core = 0.05 and r_core = 0.2 so that each power-law simulation appears twice in Fig. 11. We note again that these simulations could be rescaled to be anywhere along the black line. For cored profiles, we calculate the annihilation rates as explained in the last section. For the B-axis, we calculate an effective shock parameter that summarizes the whole history of encounters of a cusp through a single number,

which is the p-norm of the shock history. Different values of p would give different importance to stronger versus weaker shocks. For p = 1, the value of B_eff would correspond to the sum of all shock parameters. For p < 1, multiple shocks would have an enhanced effect. For p → ∞, only the strongest shock would matter. For reference, we show how such cases would appear in Appendix A4. However, we have found that p = 1.2 gives excellent predictions, and this is the value of p that we use for calculating the B-values in Fig. 11. We note that some previous studies (of NFW subhaloes) have assumed that multiple shocks can be treated by adding the changes in binding energy (Shen et al. 2022). This would imply p = 2 and would give clearly wrong results for the case of prompt cusps and probably also for NFW haloes. The results of Delos (2019b) indicate that multiple shocks also have a significantly enhanced effect (p ≪ 2) for NFW haloes. However, it is not clear that our effective description will work equally well for NFW profiles, due to their more complicated form.

The blue line in 11 shows the prediction obtained by treating the shock histories as a single shock with effective shock parameter B_eff and inserting this into equation (47). This predicts the annihilation rates correctly to within |$20~{{\ \rm per\ cent}}$|⁠.⁶ We show in Appendix A4 that this is much more accurate than results that would be obtained by only considering the strongest shock p → ∞ (errors up to factor of a few) or by considering the sum of shocks as the effective shock parameter (errors up to 50 per cent).

4.5 Summary of the effect of stellar encounters

We conclude that we can estimate the annihilation rate of cored power-law profiles that have gone through complicated shock histories simply by using equation (47) with an effective shock parameter calculated as the 1.2-norm of all the shocks. To account for the initial boundary of cusps when B ≪ B_cusp, we additionally use the boundary term from equation (44) so that we have

where we have defined

This is the main result of this section. Note that this function works accurately in all regimes – it recovers the correct suppression for B ≫ B_cusp, but it also recovers equation (10) for weak shocks with B ≪ B_cusp.

We note that the results we find here suggest that the impact of stellar encounters on the annihilation rate will be quite dramatic for cusps in the inner part of the Milky Way. We show in Fig. 12 the distribution of the reduction in annihilation luminosity relative to the initial luminosity as a function of B_*, the characteristic shock strength of a cusp’s trajectory. Here, we assume B_core = 100B_cusp, which is a typical ratio between these two parameters. We sample a large number of shock histories, considering all shocks with B > 10⁻³B_*, and evaluate the expected luminosities according to equation (50) using the effective shock parameter. For comparison, we show also the luminosities that would be obtained by considering only the strongest shock. Virtually, all cusps with B_* ≳ 0.3B_core get completely disrupted and cusps with B_* ≳ B_cusp are already dramatically affected. As expected (compare also Appendix A4), there is a significant difference between considering only the strongest shock and the full history of shocks. However, even when considering only the strongest shock, the suppression is quite strong.

4.6 The effect of smooth tides

In addition to stellar encounters, the smooth tidal field of the Milky Way can also induce mass-loss, and so a reduction in annihilation luminosity. We do not discuss the effect of smooth tides in great detail here, since this effect was already adequately incorporated by Delos & White (2022), based on the work by Stücker et al. (2023). However, since the effect of smooth tides should be included in addition to that of stellar encounters, we need to make a few modifications to their approach. We will discuss these modifications in Appendix  B. In summary, we include the effect of smooth tides by applying the adiabatic-tides model (Stücker et al. 2023) to the cusps that remain after stellar shock truncation. We find that the joint effect on annihilation radiation of a stellar shock with effective strength B_eff and the smooth tidal field is approximately equivalent to a pure shock with an effective shock parameter,

where λ is the largest eigenvalue of the tidal tensor of the spherically averaged Galactic mass distribution at the pericentre of the cusp’s orbit.

In this way, we are able to incorporate the effect of smooth tides simply through a redefinition of the shock parameter used in equation (50). In Fig. 13, we compare the effective shock parameter from the full history of shocks to the tidal contribution |$B_{\lambda } = \sqrt{42.2 \lambda }$| and to the total effect. At radii |$r \le 20 \,\mathrm{k}\rm {pc}$|⁠, the effect of encounters dominates, whereas at larger radii, the effect of the smooth tide dominates. Thus, viewed from the Earth’s position just 8 kpc from the Galactic centre, truncation by stellar encounters has a large effect on the angular distribution of the prompt cusp annihilation signal. This must be taken into account when evaluating whether these cusps affect interpretation of the GCE measured by Fermi LAT. On the other hand, stellar encounters have much less effect on the radiation seen by a distant observer since most of the mass of the Milky Way’s dark matter halo (and so most of its prompt cusps) are at larger radius.

5 RESULTS

We combine results from Section 3 on the distribution of stellar shock parameters and from Section 4 on the effect of shocks and smooth tides on prompt cusps to estimate the spatial distribution of the annihilation signal of cusps in the Milky Way.

For this, we assume the cusp population expected for a WIMP with mass |$100 \,\mathrm{G}\mathrm{e\mathrm{V}}$| and decoupling temperature |$30 \,\mathrm{G}\mathrm{e\mathrm{V}}$|⁠. This leads to values B_core and B_cusp through equations (22) and to initial annihilation rates through equation (14). We use the 10⁵ orbits and the final B_* values inferred in Section 3, but we create a different realization of a cusp and its shock history for each of 1000 different radii sampled uniformly in time along each orbit’s trajectory. For the shock history, we consider all shocks stronger than B > 10⁻⁴B_* (approximately the 10 000 strongest shocks) when evaluating the effective shock parameter B_eff. Additionally, we keep track of the smallest radius r_peri that each cusp has reached, and we evaluate the tidal field of the spherically averaged mass distribution at this point to obtain the value of λ. With this, we infer the effective shock parameter B_{eff, λ} as in equation (52) and evaluate the expected final annihilation luminosity according to equation (50). Thus, in total, we obtain 10⁸ pairs of radii and annihilation luminosities.

We show the distribution of the ratio between initial and final luminosities in the top panel of Fig. 14 as a function of radius. The percentiles of this distribution give an idea of how dramatic the effect of stellar encounters on cusps is. Typical cusps inside of the central |$10 \,\mathrm{k}\rm {pc}$| are disrupted by stellar encounters. Within the central |$3 \,\mathrm{k}\rm {pc}$|⁠, less than 2.5 per cent of cusps survive, and almost all cusps reduce their luminosities by dramatic factors (e.g. 99.7 per cent of cusps reduce their luminosity at least by a factor 5).

However, more relevant than the percentiles of the distribution is the annihilation weighted mean, since this gives the ratio between the initial annihilation profile and the final one. We show this as the red line in Fig. 14. The mean is clearly dominated by the least disrupted cusps. This is so because cusps that contribute more annihilation radiation are also more resilient to tides (compare Fig. 3) and further, since the bulk of the distribution gets completely disrupted, only the most resilient cusps with the least invasive shock histories contribute to the mean. Even so, the mean is dramatically suppressed with respect to the case where encounters are neglected. At 8 kpc, the average luminosity is already suppressed by a factor of 10. At smaller radii, the effect is even more dramatic; only about 10⁻³ of the luminosity expected in the absence of encounters remains near the galactic bulge at about 1 kpc. We show as a blue dashed line the mean annihilation reduction that would be obtained if only the effect of the smooth tide were considered (as in Delos & White 2022). This greatly overpredicts the luminosity in the central regions but is accurate at larger radii |$r \gtrsim 15 \,\mathrm{k}\rm {pc}$|⁠. The orange dashed line shows the mean if only the strongest shock is considered instead of the full history. This leads to significantly smaller, although still substantial, suppression.

In Fig. 15, we show the annihilation profile as a function of radius. Unlike the prediction when only the effect of smooth tides is considered, stellar encounters lead to a non-monotonic profile that decreases towards the centre. Its maximum is slightly outside the solar radius at around |$10 \,\mathrm{k}\rm {pc}$|⁠. At larger radii, the effect of stellar encounters becomes irrelevant so that the original profile is recovered at |$r \gtrsim 40 \,\mathrm{k}\rm {pc}$|⁠.

Let us briefly consider how these effects alter the annihilation luminosity of the Milky Way, as seen by a distant observer. For this, we simply have to sum up the luminosity of all cusps out to some truncation radius, which we assume to be |$r_{\rm {200b}} = 340 \,\mathrm{k}\rm {pc}$| (the radius where the enclosed mean density is 200 times the background density). Additionally, we add the contribution of the smooth halo, which is, however, very small (⁠|$\sim 2~{{\ \rm per\ cent}}$|⁠). We list the resulting total luminosities in Table 3. Clearly, the joint effect of tidal stripping and stellar encounters onto the total luminosity is relatively small, only reducing the estimated luminosity by 23 per cent. Stellar encounters do not much affect the total luminosity, since most of the cusps do not get close to the stars. As a result, it is not necessary to consider stellar and tidal disruption effects on cusps when inferring the total contribution of extragalactic sources to the IGRB. However, we observe the dark matter distribution of our own Galaxy from a highly biased view-point, which greatly increases the sensitivity to what happens in its inner regions.

We show in Fig. 16 the radiation profile of the Milky Way as it would be observed from the solar radius |$r = 8.2 \,\mathrm{k}\rm {pc}$| if dark matter has a significant self-annihilation cross-section. Here, we have inferred for each component the line-of-sight J-factor column-densities,

and then multiplied them by a normalization factor,

that makes the smooth halo profile agree with the GCE [this is the same factor that Delos & White (2022) have used]. Additionally, we have overplotted the data points of the observed GCE as given by Di Mauro (2021). Many error bars are too small to see, and some data points are only upper limits indicated by arrows.⁷ Stellar encounters affect the profile so strongly that there is only a very small variation with observation angle left. While the unperturbed profile and the profile including tides alone both seem in tension with the observed shape of the GCE, stellar encounters alleviate this tension and the shape of the GCE is again compatible with dark matter as an explanation. The signal from cusps in the Galactic halo makes an almost isotropic contribution – only deviating by about 30 per cent from its sky-average at its brightest angle (40°). When making conclusions from the shape of the signal profile, including the effects of stellar encounters is crucial.

However, if the GCE is due to dark matter, we still expect a significant contribution to the γ-ray background from prompt cusps. We list the sky-averaged flux that we predict for each component in Table 4. It is most interesting to compare these numbers to the observed IGRB

where the error indicates the systematic uncertainty due to foreground modelling, which we have assumed to accumulate linearly when integrating the spectrum, as measured by Ackermann et al. (2015). Here, we only consider the integrated flux of the IGRB, but in principle, the full spectral shape can be compared to the DM prediction. The observed IGRB only considers contributions from galactic latitudes b > 20°. To approximately mimic this, we have only included contributions from the Galactocentric angles larger than θ > 20° in the sky-averages listed in Table 4. Additionally, we have indicated what fraction of the IGRB each component would comprise. Finally, we have also listed the extra-galactic flux here that we estimate in the same manner as Delos & White (2022) – neglecting the effect from tidal fields, which might reduce the number by around 20 per cent as indicated by Table 3.

In comparison to the smooth tide prediction (equivalent to Delos & White 2022), the predicted background signal due to cusps in our own Milky Way goes down by a factor 2.0 and the signal is by a factor 3.5 smaller than the unperturbed one. After taking this correction into account, the total dark matter signal (smooth halo + Milky Way cusps + extra-galactic cusps) is dominated by the extra-galactic signal and goes down by roughly one third. Delos & White (2022) argue that the morphology of the signal that is expected from annihilation from prompt cusps and from dark matter decay are essentially the same. Therefore, they use constraints on the decay of dark matter from Blanco & Hooper (2019) and rescale them to obtain constraints on the dark matter annihilation cross-section. These constraints are proportional to the predicted background, thus the upper limit on the cross-section has to be increased by about one third. However, we note that these constraints depend quite strongly on how well the astrophysical contributions to the diffuse γ-ray background are understood (Blanco & Hooper 2019) and the assessment of uncertainties has not been very rigorous so far. For example, the constraints by Blanco & Hooper (2019) vary by a factor of about 4 simply by considering different assumptions about how correlated the error-bars are. Therefore, to infer reliable constraints, it will necessary to reanalyse the IGRB with a more sophisticated treatment of the systematic and statistical uncertainties.

Perhaps even more intriguing than the implications for constraining dark matter, the predicted contribution to the IGRB offers an independent test for the dark matter interpretation of the GCE. If the GCE is due to annihilation, then we predict an additional approximately isotropic γ-ray signal that would comprise about 40 per cent of the observed 1–10 GeV background for a |$100 \,\mathrm{G}\mathrm{e\mathrm{V}}$| WIMP. If we additionally consider that different WIMP models can lead to a factor ∼2 variation in predicted cusp J-factors (consider Delos & White 2022, Fig. 7 and equation 2), the total dark matter contribution should range between 20 per cent and 80 per cent of the observed IGRB. Given current uncertainties in the (apparently dominant) contribution of star-forming galaxies and AGN to the observed signal, it is unclear whether there is room for such a component. A careful reevaluation of these other contributions to the γ-ray background is clearly well motivated. The combination of total flux, the spatial morphology, and the spectrum of the IGRB can be used simultaneously to constrain the existence of such an additional component. Our work here can be used to infer templates for such an analysis. The detection or exclusion of this additional component would confirm or contradict the dark matter interpretation of the GCE and so substantially advance our understanding of dark matter itself. Firm exclusion could rule out an annihilation interpretation of the GCE, while robust detection would strongly support this interpretation, since it would be a remarkable coincidence to find the GCE and the additional IGRB contribution at the right relative level yet due to different astrophysics.

6 DISCUSSION

In this article we have modelled the effect of stellar encounters on the expected dark matter annihilation signal from prompt cusps orbiting within the Milky Way’s halo, presenting several advances in the treatment of such encounters.

First, we have developed a new method for inferring the full history of the impulsive shocks experienced by a dark substructure as it moves through the Galaxy. This method is both simpler and more general than previous approaches. While we have focused on prompt cusps in this article, our results on the distribution of shock histories (Sections 2 and 3) could be applied to conventional NFW subhaloes also, as long as the dominant shocks are in the distant-encounter regime (masses ≲1 M_⊙).

Secondly, we have performed idealized N-body simulations to infer how stellar encounters affect prompt cusps. Here, we have for the first time considered the phase-space core that any (otherwise) centrally divergent profile must exhibit. Only when this core is considered can a prompt cusp disrupt. Our simulations allow us to derive accurate formula to describe the structural effect of encounters with arbitrary combinations of shock history, core radius, truncation radius, characteristic density, and smooth tidal field. As a result, we are able to account for the joint effect of smooth tides and of any number of stellar encounters along the entire trajectory of any cusp.

While our model is relatively complete and comprehensive, a few of its assumptions could still be improved. We have assumed a static Milky Way potential, and simply integrated cusp orbits for |$10 \,\mathrm{G}\rm {yr}$| within it. A more accurate treatment would consider evolution of the host potential and of the stellar population it contains and would follow cusps from their initial formation and growth through their accretion onto precursor objects, and finally onto the Milky Way itself. While the early stages of this process remain quite uncertain, we believe that our procedure should be relatively accurate and should be conservative for effects at later times since the great majority of Milky Way stars formed (approximately) in situ in the disc and the bulge, and most of them are less than |$10 \,\mathrm{G}\rm {yr}$| old.

Another uncertain point is the precise profile of the phase-space core of prompt cusps. Here, we used a simple heuristic approach to obtain a stable profile that is consistent with the phase-space density constraint, but a rigorous investigation of the central profiles with N-body simulations would be desirable. However, we expect the results of our study to be quite robust to this uncertainty; annihilation rates depend relatively weakly on the precise radial scale and shape of the core because most cusps in the central region of the Galaxy have had encounters that put them far beyond the threshold for disruption. To significantly alter our predictions for the Galactic centre would require increasing the resilience of cusps to shocks (B_core) by 1–2 orders of magnitude, which is not plausible given the strict phase-space density constraint.

When we apply our modelling to the prompt cusp population in the Milky Way’s halo, we find that stellar encounters have a dramatic effect on any cusps that enter the star-dominated regions – the vast majority gets disrupted within the central |$10 \,\mathrm{k}\rm {pc}$|⁠, leading to a radial annihilation radiation profile that peaks near |$10 \,\mathrm{k}\rm {pc}$|⁠, dropping strongly at smaller radii. While this has little effect on the luminosity of the Milky Way as seen by a distant observer (⁠|$\lesssim 20~{{\ \rm per\ cent}}$|⁠), it strongly affects the annihilation flux observable from Earth. Stellar encounters destroy cusps so efficiently in the inner Galactic halo that the surface brightness of cusp annihilation radiation is predicted to vary only slightly between the centre and anticentre directions. As a result, it has no noticeable effect on the GCE, which therefore remains compatible with production by annihilation of the smooth dark matter distribution in the inner few kpc. This removes the issue raised by Delos & White (2022) who included cusp disruption due to the smooth Galactic tide, but not that due to stellar encounters, and hence found a total annihilation profile apparently incompatible with the GCE profile of Di Mauro (2021, see Fig. 14).

This opens up an intriguing new possibility. If the GCE is indeed due to dark matter annihilation, this implies an approximately isotropic annihilation signal from prompt cusps that would have an amplitude in the range |$20~{{\ \rm per\ cent}} - 80~{{\ \rm per\ cent}}$| of the observed IGRB, depending on the mass and decoupling temperature of the WIMP. The results of Blanco & Hooper (2019) suggest that a signal of this amplitude is inconsistent with the observed IGRB, since the latter appears to be explained entirely by emission from star-forming galaxies and AGN. We think, however, that our current results warrant a careful reevaluation of those of Blanco & Hooper (2019) since it seems conceivable that some of the fitted templates might absorb a near-isotropic dark matter annihilation signal, or that the statistical treatment of systematic errors is not yet robust enough to exclude such a signal. If such a signal is firmly ruled out, it becomes unlikely that the GCE can be ascribed to annihilation radiation. On the other hand, if it were robustly detected, this would confirm the annihilation interpretation of the GCE.

ACKNOWLEDGEMENTS

JS thanks Sten Delos for answering quickly and clearly all questions regarding the distribution and properties of cusps. JS thanks Mattia Di Mauro for providing data related to the GCE. JS and RA thank all members of the cosmology group at Donostia International Physics Center for daily discussions and for the motivating research environment. JS and RA acknowledge the support of the European Research Council through grant number ERC-StG/716151. GO was supported by the National Key Research and Development Program of China (No. 2022YFA1602903) and the Fundamental Research Fund for Chinese Central Universities (No. 226-2022-00216).

DATA AVAILABILITY

The code used to generate all results of this study other than the simulation-based analysis of Section 4 is available in an online repository under https://github.com/jstuecker/cusp-encounters. Some results are based on the adiabatic-tides code, which is also publicly available under https://github.com/jstuecker/adiabatic-tides. The simulation data from Section 4 will be shared on reasonable request to the corresponding author.

Footnotes

References

APPENDIX A: NUMERICAL CONVERGENCE

A1 Non-uniform mass initial conditions

To simulate the behaviour of our cusps, we have to truncate them at some radius r_max. We first tried to run numerical experiments with uniform mass sampling and a sharp truncation radius. With particle numbers of order 10⁶, one can typically resolve scales down to 10⁻³–10⁻² r_max before two-body relaxation effects become substantial. However, at the same time, a sharply truncated r^−3/2 power-law is far from equilibrium around the truncation radius. We found that truncation still has significant numerical effects below 10⁻¹r_max (depending on how many dynamical times are simulated). This behaviour is much worse than for profiles that are steeper in their outer regions, such as Hernquist or NFW profiles. As a result, only a small range of radii is reliably resolved and a lot of care has to be taken to get numerically converged results over a substantial radial range. Unfortunately, the transfer function we wish to measure spans several orders of magnitude in radius.

We have therefore decided to abolish the uniform mass sampling and instead divide the profile into several populations of particles with differing mass. This allows us to get reliably converged results over a dynamical range of 10⁴ in radius.

We define five populations of particles where the ith population is defined so that all of its particles have their orbital pericentres in the range,

and we choose |$\boldsymbol{r}_{\rm {split}} = (0, 1, 10, 100, 1000, 10000)^T$|⁠.

In each interval, we choose to use the same number of particles N_split. In principle, one can create such a realization by setting up several full profiles at different resolutions and then discarding all particles that do not fulfill the criterion of each population. However, this would be very inefficient for generating the innermost populations.

Instead, we can directly sample particles from the cumulative energy, pericentre distribution function. The density contribution at radius r from particles that have pericentres in a radial range r₁ < r_peri < r₂ and energy smaller than E is given by

where the boundaries L₁ and L₂ are chosen so that the integral goes only over particles which fulfill the pericentre criterion:

We find that

with

where the terms containing r₁ have to be set to zero if r < r₁ and the terms containing r₂ have to be set to zero for r < r₂.

Knowledge of this density function is enough to sample radii and energies of particles in each group. Radii can be sampled through inverse-distribution function sampling with the cumulative mass profile that can be obtained by integrating ρ and using the limit E → ∞, but maintaining a finite pericentre range. Energies can be sampled by inverse-distribution function sampling among energies using the normalized version of ρ(r, r₁ < r_peri < r₂, <E). Finally, angular momenta can be sampled by considering the cumulative angular momentum distribution function,

An implementation of this scheme for creating initial conditions for cored cusps is also publicly available in the code repository of this paper.

We show the radial density profiles of the initial conditions that have been created in this manner for a profile with r_core = 1 and N_split = 2¹⁸ as solid lines in Fig. A1. The individual components are shown as coloured lines, and the sum is shown as the black lines. Clearly, this gives an excellent initial representation of the density profile. The dashed lines in the same figure show the density profile after the simulation has been evolved for 100 times the dynamical time-scale at the core radius. The final profile still agrees excellently with the initial one, and we can see that particles of different masses have mixed very little. Note that the particle number here is still a factor 4–16 lower than for the simulations that we actually use in the main text so that we can be fairly confident that all simulations are stable and well converged.

A2 Time evolution and convergence of shocked power-law profiles

We briefly discuss the convergence aspects of the shocked power-law profiles here. We show the profile of the r_B = 20 simulation with N_split = 2²² at different output times as solid lines in Fig. A2. The profile is very stable with time over large radial ranges, but there are differences at large radii where, at early times, the profile has not yet had enough time to relax to its final form. We find that a conservative criterion for the largest radius r_max, where the profile has relaxed to its final form, is given by

where t_sim is the run-time of the simulation, and we use the mass profile at the output time (not the initial time) for M(< r). We invert this equation numerically to find r_max at each time and mark this point in Fig. A2. Clearly, this is a very conservative criterion and cuts out all parts of the profile that are yet to reach their final form.

The innermost radius r_min, where we can trust the simulations, is given by the two-body relaxation radius that we can approximate through

(e.g. Binney & Tremaine 2008), where r₀ is an estimate of the largest two-body encounter radius and ϵ is the softening. Since the dependence on r₀ is rather weak, we just use r₀ = r_min. This formula only holds strictly for systems that are made of uniform mass particles. Therefore, to be extra conservative, we use here not the actual particle number here, but the particle number as if all mass inside of r_min was made up from particles of the second-highest resolution level (from the pericentre interval [1, 10]). Since most particles in the region are actually much lower mass (and none are higher mass), the profiles are actually converged at substantially smaller radii than the thus determined r_min, which we indicate in Fig. A2. Further, we show two lower resolution simulations with 4 and 16 times fewer particles to demonstrate that the simulations are actually very well converged in the indicated regimes. In Fig. 8, we show only the reliable regions inside the range [r_min, r_max].

A3 Annihilation rates

We describe here how we compute the J-factors of the cored power-law profiles that we presented in Section 4.3. In many cases, it is numerically problematic to directly infer annihilation rates from N-body simulations of centrally divergent profiles. This is so since annihilation rates scale as the density squared and are therefore particularly sensitive to the central regions of the profile, which are, however, in general, also the regions that suffer the most from numerical inaccuracies caused by two-body relaxation, sampling noise and the effects of force softening.

However, for the simulated cored profiles, none of these are problematic, since most of the annihilation radiation comes from radii r ≳ r_core, which are resolved extremely well in our simulations. Therefore, we estimate the annihilation radiation in the following manner: We infer the radial density profile by assigning particles to 50 logarithmically spaced radial bins in the range [0.1 r_core, 10⁴r_core] (i.e. 10 bins per dex). We then numerically integrate the J-factor by combining third-order spline interpolation with a large number of integration points that we place between r_min = 0.2r_core and r_max = 6r_B. Finally, we add the contribution from r ≲ r_min by assuming that the density is uniform in this range with a value given by the average density within that radius. We show the resulting annihilation rates as the blue line in Fig. A3. We then consider different variations of the numerical procedure: using a two times larger value of r_min = 0.4r_core; using a smaller value for r_max = 3r_B; using two times the number of bins (20 per dex); and by evaluating the profile at a slightly different simulation time (5t_B earlier). Note that the last variation tests that the final profile is both stable and robust to shot noise. We show these cases in Fig. A3, together with the residuals with respect to the reference case in the bottom panel. Clearly, none of these deviations has a significant impact on the annihilation radiation. In all cases, the associated relative error is less than 3 per cent, except for the case of the strongest shock considered where the error is significantly larger, but still less than 40 per cent. These results are more than accurate enough for the purposes of this paper.

A4 Transfer functions and multiple encounters

Here, we present the transfer functions obtained for power-law profiles that have gone through up to five shocks in different sequences. Each shock is applied from a random direction. We use the power-law simulations that are listed in Table 2 and show their transfer functions in Fig. A4. Additionally, we show the transfer function that is predicted by treating the full history as a single encounter with effective shock parameter B_eff, given by the p = 1.2 norm of the shock history. Clearly, this effective shock parameter predicts the transfer function of multiple encounters reasonably well.

Finally, we present alternative versions of Fig. 11 that use different values of p for calculating the norm of the shock history. The top panel of Fig. A5 assumes p = 1, so that B_eff corresponds to the linear sum of B. In this case, the prediction (blue line) overestimates the reduction in annihilation radiation by up to 50 per cent. The bottom panel adopts the value p = 100 – effectively selecting only the strongest shock to define B_eff. Clearly, considering only the strongest shock can dramatically underpredicts the reduction in annihilation luminosity. For example, at B_eff/B_core ∼ 0.1, there are cases where the annihilation luminosity based on the full shock history is 10 times smaller than predicted from just the strongest encounter.

APPENDIX B: THE EFFECT OF THE SMOOTH TIDAL FIELD

While the main text of this paper discusses the effect of stellar encounters in great detail, the smooth tidal field also affects the annihilation luminosity of prompt cusps orbiting in the Milky Way (see e.g. Delos & White 2022). As discussed in Stücker et al. (2023), the most important parameter determining this effect is the largest eigenvalue λ of the tidal tensor at the orbital pericentre of a prompt cusp’s orbit. After a sufficient amount of time (≳10 orbits), the orbiting cusp approaches an asymptotic structure, which changes little in subsequent evolution (Errani & Navarro 2021). Stücker et al. (2023) show that the asymptotic remnant can be reasonably approximated by the adiabatic-tides model. This is an analytic description of an object’s reaction when a tidal field is applied very slowly (in the adiabatic limit) and in a spherical approximation. The corresponding code can be found online.⁸ While actual N-body simulations in the Galactic context would of course be more accurate, the adiabatic-tides model allows easy exploration of a variety of different scenarios. We note that most studies of tidal stripping have focused on NFW subhaloes, which are much less resilient to tides than prompt cusps. As a result, most previously published results cannot be applied to the case of interest here.

For a pure power-law profile with slope −1.5, the truncation radius due to smooth tides scales with a characteristic length r_λ, which may be defined by equating the cusp’s attractive force and the disruptive force due to the tidal field:

For a pure power-law initial profile, the final profile reaches zero density at the tidal radius r_t = 0.24r_λ (Stücker et al. 2023). If we set up a pure power-law profile and evaluate the structure of the remnant using the adiabatic-tides model, we find that its annihilation luminosity is given by

so that, in this case, the effect of smooth tide on the annihilation luminosity is equivalent to a sharp truncation of the initial profile at 0.78 per cent of r_λ.

To model accurately the joint effect of stellar encounters and tidal stripping, the whole encounter and tidal history should, in principle, be considered. Here, we will simply assume that we can model the joint effect approximately by first considering the effect of all stellar encounters – by truncating the prompt cusp as described in Section 4 – and thereafter applying the pericentre tidal field λ. We choose this order because it is technically simpler for us, but we expect that consideration of the effects in reverse order would be just as valid and would give similar results.

We set up a profile following the shocked power-law form,

and we apply various tidal fields λ using the adiabatic-tides model. We note that this set-up has one relevant parameter – the ratio between the two truncation scales:

We note that the annihilation luminosity in equation (B3) is equivalent to that in equation (45) for a single stellar encounter with shock parameter

Thus, in the limit B_λ ≫ B (where the tidal field sets the truncation scale), we expect the annihilation luminosity to be equivalent to that after a single shock with strength B_λ, while in the limit B_λ ≪ B, it should be equivalent to a single shock of strength B. A reasonable guess for intermediate cases is for the joint effect to be given by a weighted average of the two scales,

In principle, any p-norm of B and B_λ would be a reasonable guess: we find that p = 2 works very well.⁹ In Fig. B1, we compare the transfer functions obtained from the adiabatic-tides calculations to an effective description, assuming that the joint effect of encounters and tides is equivalent to a single shock with B_{eff, λ}. This works reasonably well in the regime where ρ/ρ_power-law > 0.5. In the regime, where this ratio is smaller, the approximation is worse; the tidal truncation is complete at the tidal radius, whereas the encounter truncation has a much longer tail. However, the tail of the profile is almost irrelevant for the annihilation luminosity and is unlikely to be correct in the adiabatic-tides description (Stücker et al. 2023).

The bottom panel of Fig. B1 shows the J-factors obtained by integrating ρ² for the adiabatic remnants over the range [10⁻⁶r_B, ∞]. Note that we chose the lower bound so that it is well in the power-law regime of the profile for all the cases considered. Apart from that, it is, of course, arbitrary, and so the absolute offset on the J-axis is also arbitrary. We show additionally the two limiting cases of a pure encounter truncation and a pure tidal truncation, together with the effective description by a single shock with strength B_{eff, λ}. Clearly, this recovers the asymptotic limiting cases as well as the intermediate regime very well.

To verify that this gives a reasonable description of the joint effect of encounters and of the smooth tidal field in all relevant cases, we must also check that it applies to cored initial profiles. Here, we only consider the case where B_λ ≫ B so that we do not have to deal with two scales at the same time, but only with a single scale B_λ/B_core. We set up cored profiles as described in Section 2.2 and apply tidal fields of varying amplitudes through the adiabatic-tides model. We show the corresponding J-factors in Fig. B2.

Clearly, the effective description works very well for cored profiles also, as long as B_{eff, λ} ≲ 0.2B_core. Beyond that scale, the cored profiles reach an earlier disruption threshold of B_{dis, λ} = 0.33B_core, which is a factor two smaller than the encounter disruption threshold. While it would certainly be possible to improve our effective model further to incorporate this reduced disruption threshold, this is unnecessary, since such large tidal fields are not reached in the Milky Way – especially not at the large radii where smooth tides dominate over stellar encounters and |$B_\lambda \ll 1 \,\mathrm{k}\mathrm{m}\mathrm{s}^{-1} \rm {pc}^{-1}$| typically. We conclude that treating the joint effect of stellar encounters and smooth tides as that due to a single encounter with a shock parameter B_{eff, λ} is an excellent approximation within the adiabatic-tides framework. We note that this may slightly overestimate the effects of smooth tides, since the adiabatic-tides predictions are often a slight overprediction in practice (Aguirre-Santaella et al. 2023; Stücker et al. 2023). However, since the effects on cusp annihilation luminosity are in any case rather weak, this seems acceptable.