Martin Stref
email: martin.stref@umontpellier.fr

Julien Lavalle
email: lavalle@in2p3.fr

Modeling dark

Keywords: numbers: 12.60.-i, 95.35.+d, 96.50.S-, 98.35.Gi, 98.70.Sa

published or not. The documents may come

I. INTRODUCTION

While the long-standing issue of the origin of dark matter is still pending, many experiments involved in this quest have recently reached the sensitivity to probe the relevant parameter space for one of the most popular particle candidates, the WIMP, which finds specific realizations in many particle physics scenarios beyond the standard model (e.g. [1][2][3]). Among different search strategies, indirect DM searches (e.g. Refs. [4][5][6][7]) are becoming quite constraining for WIMPs annihilating through s-waves. This is particularly striking not only for indirect searches in gamma rays (e.g. [8][9][10]), but also in the antimatter cosmic-ray spectrum [11], both with positrons (e.g. [12]) and antiprotons (e.g. [13]). For indirect searches, the way the Galactic dark matter halo is modeled is a fundamental piece in deriving constraints or testing detectability. For direct DM searches, whether the local DM density is smooth or may contain inhomogeneities has also important consequences (see e.g. [14]).

A generic cosmological consequence of the WIMP scenario (among other CDM candidates) is the clustering of dark matter on very small, subgalactic scales, when the universe enters the matter domination era (e.g. [15][16][17][18][19][20][21][22][23][24], and [25] for a review). Both analytic calculations (see a review in e.g. [26]) and cosmological simulations (e.g. [27][28][29][30][31]) show that many of these subhalos survive in galaxies against tidal disruption, and further constrain their properties. Consequently, the DM halo embedding the MW, if made of WIMPs, is not a smooth distribution of DM, but instead exhibits inhomogeneities in the form of many subhalos or their debris. In the context of selfannihilating DM candidates, this leads to the interesting consequence of enhancing the average annihilation rate with respect to the smooth-halo assumption [16]. Generic methods to account for a subhalo population in the DM annihilation signal predictions were originally presented in [32,33] for gamma rays, and in [34,35] for antimatter cosmic rays.

While subhalos are now very often included when deriving constraints from the Galactic or extragalactic diffuse gamma-ray emissions (see e.g. [33,[36][37][38][39], and a review in [40]), this is still barely the case for the antimatter channels (e.g. [13,41]). In the latter case, although it was shown that subhalos could not enhance the predictions by orders of magnitude [35,36], the precision achieved by current experiments (see e.g. [42][43][44][45] for antiproton measurements) implies that even small changes in the predicted fluxes could still have strong impact on constraints on the WIMP mass. In this paper, our aim is to provide a dynamically self-consistent model of a Galactic subhalo component in order to improve the constraints derived on s-wave annihilating WIMPs.

The paper develops as follows. The main part of our study is described in Sect. II, where we introduce the dark halo setup including both a smooth and a subhalo component, and where we discuss the tidal effects induced by both baryons and dark matter. We then discuss the mass profiles, the dark matter annihilation profile, and the corresponding differential and integrated annihilation boost factors in Sect. III, which can be used in indirect detection studies. In that part, we also quantify the theoretical uncertainties coming from using different Galactic mass models, different tidal cut-off criteria, or other subhalo population properties. We conclude and present our perspectives in Sect. IV.

II. THE MILKY WAY DARK HALO AND ITS SUBHALO POPULATION

In this section, we propose a self-consistent method to constrain the subhalo population of the MW dark halo and to derive therein the DM annihilation rate including all components. This method subscribes to two main principles: (i) accounting for existing dynamical constraints in the MW; (ii) starting from general assumptions, then comparing to and calibrating on highresolution cosmological simulations only a posteriori. In the following, any halo mass m will, unless specified otherwise, express the mass contained within a sphere of radius r 200 such that

m = m(r 200 ) = m 200 = 4 π 3 (200 × ρ c ) r 3 200 , (1)
where ρ c is the critical density of the universe, which we compute from the best-fit Hubble parameter obtained by the Planck collaboration (combined analysis), H 0 = 67.74 km/s/Mpc.

A. Dark halo model

The most basic and obvious assumption one can make about the DM distribution in the Galaxy is that the DM density profile ρ tot can be splitted into two components, one smooth, ρ sm , and another made of subhalos, ρ sub , such that at any position x

ρ tot ( x) = ρ sm ( x) + ρ sub ( x) , (2)
such that the total dark mass is given by

M 200 = V200 dV ρ tot ( x) , (3)
where V 200 is the spherical volume delineated by the associated pseudo-virial radius R 200 . Furthermore, to get reliable predictions for DM annihilation signals, it is important to account for existing dynamical constraints on the DM profile -it is barely justified to import physical quantities (scale radius, scale density, density at 8 kpc, etc.) directly from cosmological simulations which are by no means reliable descriptions of our Galaxy, only using generic properties makes sense.

There has been increasing interest on this aspect in the recent years (e.g. [46][47][48][49][50][51][52][53][54][55][56]), such that modeling the dark halo in the context of DM searches can strongly benefit from the obtained results. We stress that global dynamical studies provide constraints on ρ tot , but not on ρ sm and ρ sub separately.

From cosmological structure formation (see e.g. [57][58][59]), we know that galactic halos form rather late (z ∼ 6) with respect to the smallest-scale halos expected in the WIMP scenario (z ∼ 80). It is therefore reasonable to assume that the smooth and subhalo components follow the same spatial distribution when the Galactic halo forms. Then, as the Galaxy evolves, several changes occur: (i) further subhalos are accreted, and (ii) subhalos may experience mergers, stellar encounters, and tidal disruptions. Since the former phenomenon also concerns the smooth component, it should not modify the overall picture (subhalos may be considered as test particles among others). However, the latter must be taken into account, since it will reduce the subhalo number density in regions close to the terrestrial observers. This approximate trend is actually what is found in very high-resolution cosmological simulations, where the subhalo number density is shown to depart from the overall DM distribution essentially in the central regions of galaxies [29,60,61]. In the same references, the global DM profile (including subhalos) is found to be consistent with the seminal earlier results obtained by Navarro, Frenk and White [62] (hereafter NFW) and subsequent refinements (e.g. [60,61,[63][64][65]). Inner cored profiles can also be found as a result of efficient feedback originating in star formation and supernova explosions [66,67].

All this suggests the following method to try to build a self-consistent dark halo with a substructure component: (i) assume a global DM halo profile ρ tot constrained by dynamical studies; (ii) start with a subhalo population tracking the smooth halo, such that both ρ sub ∝ ρ tot and ρ sm ∝ ρ tot ; (iii) plug in tidal disruption such that the mass contained in disrupted subhalos and in the pruned part of the survivors is transferred to the smooth halo component; (iii) compare/cross-calibrate the final result with/onto high-resolution cosmological simulations. Before we translate this method in terms of equations for DM searches, we need to figure out how to express the mass density profile ρ sub associated with subhalos. In practice, the smooth DM component will merely be determined from Eq. ( 2) as ρ sm = ρ tot -ρ sub , after having set ρ sub .

B. Accounting for dynamical constraints

As a template and dynamically constrained global dark halo, we will use the best-fit MW mass model obtained by McMillan [48] (M11 hereafter), which turns out to be fully consistent with more recent studies (e.g. [51,52,68]), while rather simple to implement. This model was derived from a Bayesian analysis run upon several observational data sets, photometric as well as kinematic, restricting to the terminal velocity curves measured for longitudes |l| > 45 • -this model does not address the complex structure of the very central regions of the MW, nor does it include any atomic or molecular gas component (we will use mass models including gas components in Sect. III C).

Global dark halo and baryons

M11 assumes a spherically symmetric NFW profile, given in terms of the general αβγ parameterization [69,70] as

ρ tot (r) = ρ αβγ (r) ≡ ρ s (r/r s ) -γ {1 + (r/r s ) α } - (β-γ) α , (4)
with (α, β, γ) = (1, 3, 1) for an NFW profile. The M11 best-fit values for the scale density ρ s and the scale radius r s are given in Tab. I. For the sake of comparison, we also introduce the Einasto dark matter profile [63,71]:

ρ ein (r) = ρ s exp - 2 α e r r s α e -1 . (5)
This profile halo was used in a dynamical study complementary to and consistent with M11, presented in Ref. [72] (CU10 hereafter). The associated parameters are also given in Tab. I. Irrespective of the MW mass model, these dark matter profiles will also be used to describe the subhalos. In the following, we will use M11 as our reference case. Since we also aim at considering the baryonic components when dealing with tidal effects (see Ref. [73] for a recent review), we provide the axisymmetric M11 bulge-disk density model below (with the convention r 2 = R 2 + z 2 ), where subscript b refers to the bulge and d to the disk:

ρ b (R, z) = ρ b (1 + r /r b ) α b exp - r r c b 2 , ρ d (R, z) = Σ d 2 z d exp - R R d - |z| z d , (6)
where r ≡ R 2 + (z/q) 2 ), q is the axial ratio, Σ d is the disk surface density, and the other parameters are scale parameters. All parameters are given in Tab. II, where a two-component disk is explicit (thin and thick disks) -note that the above disk parameterization can also be relevant to additional gas components (see Sect. III C). Since the model was not fitted against observational data featuring the central regions of the Galaxy, the bulge parameters but ρ b are actually fixed to those obtained in Ref. [74]. Note that such a disk profile can also be relevant to describe gaseous components, which have not been included in M11.

It will turn useful to have a spherical approximation of the disk density when dealing with global tides (see Sect. II D 1). We readily derive it by demanding that the disk mass inside a sphere of radius r equals the actual disk mass inside an infinite cylinder of radius R. It reads

ρ d,sph (r) = Σ d 2 r exp - r R d . (7)
One may find similar expressions with

R d ↔ R 2 d + z 2 d
(e.g. Ref. [START_REF] Binney | Galactic Dynamics: Second Edition[END_REF]), but using one or another has absolutely no impact in this study.

The overall subhalo component

The very presence of subhalos in the Galactic host halo leads to strong DM inhomogeneities, so defining a global regular mass density function for subhalos implicitly implies averaging over a certain volume. In the following, we will assume that subhalos are independent objects described over a phase space w n that includes their position x, mass m, and concentration c (we define these parameters in Sect. II C), such that their number density reads

d n N dw n = N sub K w dP V dV dP m dm dP c dc . (8)
Parameter K w is a normalization constant determined by the following closure relation:

dw n d n P w dw n = K w ⇔ dw n d n N dw n = N sub , (9)
where N sub is the total number of subhalos over the whole phase space embedded in the host dark halo. Each individual probability distribution function (pdf) dP wi /dw i , where w = V, m, c (V is the physical volume), is defined such that it is normalized over its phase-space subvolume δW i as wi,max∈δWi wi,min∈δWi

dw i dP wi dw i = 1 . (10)
We emphasize that as long as these individual pdfs are uncorrelated, K w = 1, but this is generally not the case.

In particular, when tidal effects are considered, then each subhalo is featured by a tidal radius r t which depends on its initial mass m, its position x in the Galactic halo, and its concentration c -we will detail the individual pdfs and discuss tidal disruption of subhalos in Sect. II D. Therefore, tidal effects will induce an explicit correlation between the pdfs, making the subhalo phase space intricate and non-trivial, and leading to K w = 1. [CU10], and [68] [M16]). ). † : The CU10 HI and HII gas disks are inferred from old data points, so we adopt the same parameterization as in M16 for simplicity -this has negligible impact on the final results.

MW mass model q α b r b r c b ρ b R d z d Σ d [kpc] [kpc] [M /pc
However, we can still self-consistently define the global subhalo mass density profile as

ρ sub ( x) = N sub m t ( x) dP V ( x) dV , with m t ( x) ≡ 1 K w mmax mmin dm dP m dm cmax cmin dc dP c dc m t (r t (c, m, x), m, c) , (11)
where m t is the subhalo mass contained within the tidal radius r t , to be contrasted with m which is the mass contained inside an approximate virial radius assuming a homogeneous background matter, usually r 200 (see Sect. II C for details). Symbol is not the average over the mass and concentration subpart of the phase space because of the normalization K w , which is calculated over the full phase space. The real mean mass (or any other quantity depending on mass and concentration) is actually given by

m t ( x) =
The dependence of the tidal radius r t on position, mass, and concentration will be discussed in Sect. II D. Notice that there is also spatial dependence hidden in the denominator above, as the minimal concentration will be shown to be spatial dependent in Sect. II D 4.

The total mass M tot sub in the form of subhalos is thereby given by

M tot sub = N sub host halo dV m t ( x) dP V ( x) dV . ( 13
)
It will also turn useful to define the total subhalo mass contained in a specific subhalo mass subrange

∆ m 12 = [m 1 , m 2 ] ⊂ [m min , m max ], M ∆ m 12 sub = N sub host halo dV m t ∆ m 12 ( x) dP V ( x) dV , (14)
with

m t ∆ m 12 ( x) ≡ 1 K w m2 m1 dm cmax cmin dc m t (r t (c, m, x), m, c) dP m dm dP c dc . (15)
From Eqs. ( 3) and ( 14), we can then define the total dark mass fraction in the form of subhalos within the

mass range ∆ m 12 , f ∆ m 12 sub ≡ M ∆ m 12 sub M 200 . ( 16
)
We will actually use this fraction to normalize our subhalo population and to calculate N sub , which we discuss in the next section.

Calibration of the subhalo component

The overall subhalo distribution being defined, we need to calibrate the subhalo mass content. To proceed, we will first rely on cosmological simulation results, which provide pictures of MW-like halos at redshift z = 0, with subhalo populations that have already experienced all relevant dark matter-only non-linear disruption or pruning processes (see e.g. [29,60,61,75,76]). Calibration from first principles is also possible, while more involved and subject to large theoretical uncertainties; this gives similar constraints though, as reviewed in Ref. [26]. Besides, it is well known that cosmological parameters have significant impact on the global and structural properties of subhalos, especially the matter abundance Ω m , the normalization of the power spectrum σ 8 , and the inflation spectral index n s (see e.g. [77][78][79][80] and [75,81,82]) -larger values of the former lead to more concentrated halos on all scales, while larger values of the latter increases the power on small scales. Therefore, we should favor references with input cosmological parameters not too far from the most recent estimates. In particular, the Planck mission [83] has provided combined constraints, Ω m 0.31, σ 8 0.82, and n s 0.97, directly relevant to the structuring of DM subhalos -note, though, that there are still mild tensions between different cosmological probes (see e.g. Ref. [84] for a recent illustration). This makes the Via Lactea II ultra-high resolution simulation [60] (VL2 hereafter) a rather conservative reference, since it was run with WMAP-3 best-fit parameters, Ω m 0.24, σ 8 = 0.74, and n s = 0.951 [85]. For comparison, the Aquarius simulation series [29] were run with Ω m = 0.25, σ 8 = 0.9, n s = 1, and with a spatial resolution similar to VL2.

We will use the VL2 results to calibrate the subhalo mass fraction defined in Eq. ( 16), but the method presented below can be used with any calibration source. In particular, the authors of VL2 provide the cumulative number of subhalos N VL2 (> v max ) as a function of the maximal velocity v max . Note that their census is made up to the host halo radius R 50 and not R 200 , as defined in Eq. (1). A very good fit to this measurement is obtained in the range v max ∈ [3 km/s, 20 km/s] = [v max,1 , v max,2 ] with the following parameterization [60]:

N VL2 (> v max ) = 0.036 v max v max,host 3 , (17)
with v max,host = 201 km/s the maximal velocity of the host halo. The maximal velocity is directly measured in simulations, and is related to the (sub)halo profile through the relation

v max = max G M (r) r = G M (r vmax ) r vmax , (18)
which defines the radius r vmax , and which can easily be computed for any choice of subhalo profile once its parameters are fixed (total mass m, concentration c, and position x in a host halo, if relevant). We can therefore calculate the effective mass fraction contained in subhalos within the mass range [m 1 = m(v max,1 ),

m 2 = m(v max,2 )] as f ∆ m 12 sub,VL2 = M ∆ m 12 sub,VL2 M VL2 200 ( 19
)
with

M ∆ m 12 sub,VL2 ≡ vmax,2 vmax,1 dv max m(v max,1 ) dN VL2 dv max .
This is an effective mass fraction since it is computed from the pseudo-virial mass m = m 200 instead of the tidal mass m t used in Eqs. ( 14) and ( 16), which is unknown here. The subtlety is that the subhalo population under scrutiny here has actually already experienced tidal effects, which we will ultimately have to account for.

Assuming that VL2 subhalos are well fitted by NFW profiles and taking a concentration function matching the VL2 results (we actually take the VL2 concentration fit proposed in Ref. [36]), we find that the total effective subhalo mass in the mass range

[m(v max,1 ) = 3.14×10 6 M , m(v max,2 ) = 1.25×10 9 M ] is M ∆ m 12
sub,VL2 = 2.24 × 10 11 M . Taking the global VL2 halo mass M 50 = 1.93 × 10 12 M , we obtain an effective mass fraction of f50 = 11.6%. If we further assume that the subhalo number density profile spatially tracks the global halo density profile in the outer halo regions, then the extrapolation to R 200 < R 50 is trivial, f200 f50 . For completeness, we can express this result in terms of a relative mass range, as different halo models come with different global masses. From the Einasto profile fitted on the VL2 host halo (see the caption of Fig. 2

∀ m 200 M 200 ∈ m 1 M 200 = 2.2 × 10 -6 , m 2 M 200 = 8.8 × 10 -4 . (20)
We note that this is fully consistent with the semianalytic result obtained in Ref. [75], which sets this fraction to ∼10% for subhalos in the mass range 10 -5 < m/M < 10 -2 , assuming that dN/dm ∼ ∝ m -2 (see also Refs. [36,[86][87][88]).

In practice, we will match the fraction f

∆ m 12
sub defined in Eq. (16) to the above f ∆ m 12 sub by replacing the tidal mass m t by m 200 in Eq. (15). An important subtlety is that we will only integrate over the subhalo population which has not been disrupted by tidal interactions with the host dark halo (so-called global tides in Sect. II D 1). Indeed, these interactions are at play in VL2, so this normalization procedure must take them into account. Note that the calculation of the phase space normalization factor K defined in Eq. ( 9) must also include these tidal cuts, which are position-mass-concentration-dependent. In practice, this is done by integrating the concentration function from a minimal concentration, c min (m, R) ≥ 1, which is set by the tidal disruption model and depends on the subhalo mass and its position in the Galaxy (see Sect. II D 4).

We emphasize that since VL2 is a DM-only simulation, the above normalization can only be used to calculate the total number of subhalos before plugging in tidal stripping from the baryonic components. This is actually very fortunate because this really allows us to predict the baryonic effects (at least those related to tidal stripping), instead of trying to reproduce them. Indeed, we stress that the way baryons are implemented in simulations is still highly debated in numerical cosmology (see e.g. Ref. [89]). We will deal with baryonic tides only in a second, independent step.

To summarize the normalization procedure, we first fix the total number of subhalos before baryonic tides by matching f ). In a second step, we plug in baryonic tides, which turn out to be dominated by disk shocking effects. It is easy to show that the final number of subhalos N will merely be given by N = (K /K) × N ≤ N , where K (K ≤ K) and N are the phase-space normalization and the total number of objects, respectively, before (after) including baryonic tides. This relies on matching the global subhalo mass density in the outskirts of the Galaxy, where baryonic effects can be neglected. Tidal effects will be discuss in detail in Sect. II D.

C. Global and internal subhalo properties

In this section, we specify the global and internal properties. The latter are mostly featured by the inner density profile ρ of a subhalo and its specific concentration c. For the density profile, we assume spherical symmetry and adopt the NFW shape given by Eq. ( 4), with (α, β, γ) = (1, 3, 1) as our default configuration, unless specified otherwise. We define the concentration parameter c as

c = c 200 = r 200 r -2 , ( 21
)
where r -2 is the radius at which the logarithmic slope d ln(ρ)/d ln(r) = -2. In the αβγ case, we have

κ ≡ r -2 r s = (β -2) (2 -γ) -1 α ∀ 0 < γ < 2 , β -γ α > 2 , (22)
such that r -2 = r s for an NFW profile. The same is readily obtained for an Einasto profile. The concentration parameter will play a significant role not only in ruling the subhalo annihilation rate, but also in characterizing the resistance of subhalos to tidal stripping. We now formulate the overall mass and the tidal mass,

m = m 200 = m(r 200 ) = 4 π r 3 s c 0 dx x 2 ρ(x r s ) m t = m(r t ) = 4 π r 3 s xt 0 dx x 2 ρ(x r s ) ζ(x t ) , (23)
where the dimensionless parameter x ≡ r/r s , and x t ≡ r t /r s (r t is the tidal radius). Function ζ(x t ) takes values 0 or 1 to account for the potential tidal disruption of the subhalo. We will specify this function as well as our definition of the tidal radius in Sect. II D. On the same vein, we also introduce the subhalo effective annihilation volume ξ,

ξ ≡ ξ 200 = ξ(r 200 ) = 4 π r 3 s c 0 dx x 2 ρ(x r s ) ρ 0 2 ξ t = ξ(r t ) = 4 π r 3 s xt 0 dx x 2 ρ(x r s ) ρ 0 2 ζ(x t ) , (24)
which provides a measure of the WIMP annihilation rate in a subhalo. It actually quantifies the volume a subhalo would have to supply its annihilation rate if it were an homogeneous sphere of reference density ρ 0 . In practice, we will set ρ 0 = ρ , unless specified otherwise. This is particularly convenient a choice in the context of indirect DM searches with antimatter cosmic rays [34][35][36]. It is similar to the definition of the J(ψ) luminosity factor in the context of gamma-ray searches [90].

We now introduce key physical quantities to describe bounded systems, which we will use when addressing the tidal effects in Sect. II D. We first define the gravitational binding energy, i.e. the minimum energy to unbound the system, as

E b (r t ) = 4 π G N rt 0 dr r ρ(r) m(r) , (25)
where r t is the subhalo tidal radius, ρ(r) its mass density at radius r, and m(r) its mass inside r; the binding energy is defined as positive. Alternatively, we also introduce the potential energy of a bounded system:

U g (r t ) = 2 π G N rt 0 dr r 2 ρ(r) φ(r) , (26)
where we have used the gravitational potential

φ(r) = φ(r, r t ) ≡ φ(r) -φ(r t ) (27) with φ(r) = -G N ∞ r dr m(r) r 2 ,
taking into account that subhalos have finite extensions set by their tidal radii r t . This potential takes an analytic expression for an NFW profile, easy to derive and available in any relevant textbook. Both the binding energy and the (absolute value of the) potential energy scale similarly with r t for NFW profiles, very roughly ∝ r 2 t when r t r s , and ∝ ln(r t ) when r t r s . In the following, we will provide more details on the overall global phase space characterizing our subhalo population model. We will discuss the concentration function in Sect. II C 1, the mass function in Sect. II C 2, the spatial distribution in Sect. II C 3, and tidal effects and induced correlations in Sect. II D.

An important point, we will assume that subhalos are independent of each other, which means that each physical quantity (mass, annihilation volume, etc.) can be dealt with as a random variable over the global phase space. This will allow us to compute different moments of any observable and thereby estimate the associated statistical uncertainty.

Concentration function

The concentration of DM (sub)halos has long been studied in the literature (see e.g. [29,61,81,[91][92][93][94][95][96][97][98]). In Ref. [97] (SCP14 hereafter), the authors compared the concentration model of Ref. [96] to various sets of cosmological simulation data, spanning a large range of subhalo masses, notably from ∼ 10 -6 M (from Refs. [27,30,99]), and also including the VL2 data. It turns out that in spite of the slightly different input cosmological parameters, these data can be relatively well described by the model within statistical errors -note that the rather large values of Ω m and σ 8 inferred from the recent Planck data would even favor a more optimistic modeling [81,82]. The authors of SCP14 also provide a fitting function of the central concentration value, inspired by Ref. [35], which is quite convenient for our purposes:

c(m, z = 0) = 5 i=0 c i ln m h -1 M i , (28)
with c i = [37.5153, 1.5093, 1.636 • 10 -2 , 3.66 • 10 -4 , 2.89237 • 10 -5 , 5.32 • 10 -7 ], which gives values from ∼ 65 at the lower subhalo mass edge ∼ 10 -10 M to ∼ 10 at the bigger mass edge ∼ 10 10 M . This is reminiscent from the fact that smaller objects have formed earlier, in a denser universe, and this further induces a larger luminosity-to-mass ratio ξ/m ∼ ∝ c 3 for lighter objects.

Furthermore, there is a scatter about this central value related to the fact that structure formation is a statistical theory of initial density perturbations. The associated pdf can be very well described by a log-normal distribution (see e.g. Refs. [93-95, 100, 101]):

dP c dc (c, m) = 1 K c exp -(c-c(m)) 2 2 σ 2 c c 2 π σ 2 c , (29)
where we will fix the variance in log space to σ c = σ dec c × ln (10), with σ dec c = 0.14, a mass-independent and rather generic value consistent with several detailed studies (e.g. [93,94,101]). Parameter K c = K c (m) allows a normalization to unity over the range considered in this work, that we set in practice to c ∈ [1, exp(ln(c(m)) + 8 σ c )]. The lower value c min = 1 is constrained by the definition of r -2 , which is no longer consistent when the halo extent is found smaller in the case of both NFW and Einasto profiles. This does not mean that subhalos for which one cannot specify r -2 are nonphysical, this is just a limit of our definition of the concentration itself [28,102]. However, this has no impact on the observables we will be dealing with along this article, for which only large values of the concentration will be relevant. Note that, according to Eqs. ( 28) and ( 29), the central concentration c and the averaged concentration c do not coincide:

c(m) = dc c dP c dc (c, m) (30)
c(m) e σ 2 c 2 1.05 c(m) = c(m) .
To summarize, once the density profile is fixed, the inner structure of a subhalo is fully determined from its mass m and its concentration c. The former gives r 200 from Eq. ( 1), and the latter provides the scale radius r s and the scale density ρ s from Eqs. ( 21) and (23).

Finally, we emphasize that concentration will play a crucial role in characterizing the resistance of subhalos to tidal effects, as we will discuss in more detail in Sect. II D.

Mass function

An important part of the subhalo phase space consists in the mass function. The Press and Schechter (PS) formalism and its extensions (see Refs. [57,58,75,77,79,[103][104][105]), in the frame of hierarchical structure formation and standard cosmology, provide the basic theoretical paradigm to understand why cosmological simulations exhibit power-law (sub)halo mass functions down to very small subhalo masses (see e.g. Refs. [29,60,61,86,106]). The mass index is actually related to the index of the power spectrum of primordial perturbations, and remains weakly constrained on the very small scales relevant to DM subhalos (for recent studies, see e.g. [107,108]). However, we will still assume that the mass function is regular over the whole subhalo mass range, as expected in standard cosmology, such that the initial mass pdf may be written as a simple power law,

dP m dm = K m m M -αm and mmax mmin dm dP m dm = 1 , (31)
where The mass index α m is typically expected 2 as a prediction of the PS theory with standard cosmological parameters, which is actually recovered in cosmological simulations [29,60,61,86,106]. In the following, we will assume 1.9 ≤ α m ≤ 2, unless specified otherwise. We emphasize that the actual subhalo mass distribution, which should incorporate tidal stripping and disruption, and depends on the tidal subhalo mass m t rather than on m 200 , is not directly described by Eq. ( 31). Indeed, tidal effects will make m t become position dependent, and thereby the subhalo mass range too. Nevertheless, the procedure presented in Sect. II A (see Sect. II B 2 and Sect. II B 3) includes all this self-consistently while still based on Eq. ( 31) as the initial mass function.

K m = K m (m min , m max )

Spatial distribution

The spatial distribution of subhalos in the Galaxy is a very important input in this work because it will allow to compute the local number density of subhalos, which will itself set the local annihilation boost factor, relevant, for instance, to indirect DM searches with antimatter cosmic rays. As for the mass function introduced above, we will also define the initial spatial distribution, which will further be distorted by tidal effects from the procedure defined in Sect. II A. As argued above, since small-scale subhalos have already virialized when the Galaxy forms, it is reasonable to match the initial spatial pdf to the global dark halo profile, such that

dP V ( x) dV = ρ tot ( x) M 200 , (32)
where M 200 is the global dark halo mass within R 200 , and ρ tot is the global DM density profile discussed in Sect. II B 1. This pdf is normalized to unity within a sphere or radius R 200 by construction. Of course, tidal effects will strongly distort this initial distribution because of tidal disruption, such that the effective and real spatial distribution of subhalos will eventually not look like Eq. (32). Actually, tidal effects will make this spatial distribution become mass-dependent, exactly as the actual mass function becomes spatial dependent, such that the mass and spatial distributions are fully intricate in practice (tidal effects are discussed in Sect. II D). Therefore, even though we do use Eq. (32) to formally describe the initial spatial distribution, the effective spatial distribution is still self-consistently determined through the procedure described in Sect. II B 2 and II B 3.

D. Tidal effects

Tidal effects play a fundamental role in shaping the phase space relevant to Galactic DM subhalos as defined in Eq. ( 8). As discussed above, they affect their mass, concentration, and spatial distributions, and will thereby distort and mix the pdfs defined in Eqs. ( 31) and ( 29) by pruning and disrupting subhalos. In the following, we describe in detail the way we implement these effects, which are critical to our final results.

Many studies have been, and are still being, carried out on this topic (e.g. Refs. [61,76,[109][110][111][112][113][114][115][116][117][118][119][120][121]). In this study, we will mostly consider two distinct effects: tidal stripping from the overall Galactic potential, and tidal shocking by the Galactic disk, which are known to be the most significant processes (see e.g. Ref. [26]).

Implicit in what follows, we will assume that any derived tidal radius cannot exceed r 200 , such that formally, throughout all this paper, we will always impose

r t = Min {r t , r 200 } . ( 33
)
1. Global tides from the host halo Tidal effects generated by the host Galactic halo induce a pruning of subhalos that can be accounted for by setting the actual spatial extent of a subhalo to its tidal radius. In the simplest approximation where both the host halo and its subhalo are considered as point-like objects, and taking into account the centrifugal force, the tidal radius can be defined as [76,[122][START_REF] Binney | Galactic Dynamics: Second Edition[END_REF][START_REF] Read | [END_REF]

r t• = r t• (R, m, M ) = m t 3 M 1/3 R , (34)
where M and m t are the point masses of the whole host galaxy and the subhalo, respectively, and R is the radial position of the subhalo in the host galaxy. Note that m t = m(r t• ) features the above equation, not m 200 . This formula can be generalized to the case of objects orbiting galaxies with continuous mass density profiles, more relevant to our case, as (see Ref. [START_REF] Binney | Galactic Dynamics: Second Edition[END_REF])

r t = r t (R, m t , ρ tot (R)) (35) =    m t 3 M (R) 1 -1 3 d ln M (R) d ln(R)    1/3 R ,
where M (R) is the host galaxy mass within a radius R, which depends on the global mass density profile ρ tot . This equation may be solved iteratively as it implies the tidal subhalo mass m t = m(r t ) defined in Eq. ( 23), and is shown to provide a rather good description of a subhalo radial extent in DM-only cosmological simulations (see e.g. Ref. [29]).

For completeness, we may also introduce an empirical tidal radius definition where we just delineate the subhalo by the radius at which its density equals the overall density locally, i.e.

r t such that ρ(r t ) = ρ tot (R) . ( 36
)
We finally stress that when baryons are included, they also contribute to ρ tot and thereby to M (R) in the equations above [for the baryonic disk, we will use the spherical approximation of the density, given in Eq. ( 7)]. We will discuss the impact of using one or the other definition in Sect. II D 4. Besides, note that although global tides from the host halo are indeed important in the outskirts of the Galaxy, other processes become more and more efficient in the inner regions, as the ratio of baryons to dark matter increases, as we will see below.

Baryonic disk shocking

An important source of destructive gravitational interaction arises during disk crossing, where subhalos can acquire a substantial amount of kinetic energy which can unbind them (see Refs. [110,114,116,125,126]). Termed disk shocking, this effect dominates over more local destructive effects like encounters with stars, and is actually found the most efficient subhalo disruption mechanism in the luminous part of spiral galaxies [26]. These effects are much more tricky to include than those discussed in Sect. II D 1.

Below, we discuss the physical steps that allow to account for disk shocking in a subhalo population model. We first review the seminal results obtained in Ref. [125] by Ostriker, Spitzer, and Chevalier, and further extended in e.g. Ref. [126], which were related to the study of Galactic stellar clusters.

We wish to evaluate the kinetic energy gained by a WIMP orbiting a subhalo only subject to the gravitational field of the Galactic disk during one crossing. Assuming the disk is an infinite slab (radial boundaries are sent to infinity), then the disk gravitational force field is directed along the axis perpendicular to the disk and of unitary vector e z , so the z coordinate is the only relevant here. This is a fair approximation when a subhalo is about to cross the disk. Setting x the full 3D WIMP position and x 0 the subhalo center position, the change in the WIMP velocity along the z axis and in the subhalo frame reads:

dv z dt = d( ˙ x -˙ x 0 ) dt • e z = g z,disk (Z) -g z,disk (Z 0 ) = g z,disk (Z 0 + δZ) -g z,disk (Z 0 ) ≈ δZ dg z,disk (z) dz , (37)
where we have defined δZ ≡ Z -Z 0 = ( x -x 0 ) • e z , and where the latest line is merely obtained from a Taylor expansion to first order. We have used the disk gravitational force field g z,disk , which can be inferred from the baryonic disk profile introduced in Eq. ( 6),

|g z,disk (R, z)| = 4 π G N z d ρ d (R, z) . (38)
Eq. ( 37) can further be integrated over the disk cross-ing time δt = t > -t < to get the net velocity change ∆v z ,

∆v z = t> t< dt δZ dg z,disk (z) dz (39) ≈ δt δZ (g z,disk (z(t > )) -g z,disk (z(t < )) z(t > ) -z(t < ) = δZ δt z(t > ) -z(t < ) = 1 V z × 2 g z,disk (z = 0) ,
where V z is the component of the subhalo velocity perpendicular to the disk. This approximation is licit as long as δZ does not vary much over the crossing time (i.e. the WIMP orbital time in the subhalo is much longer than the disk crossing time) and as long as the modulus of the gravitational force field remains close to its maximal value (but for the flip of sign when crossing z = 0). This is known as the impulsive approximation.

We can therefore derive the net average gain in kinetic energy per unit WIMP mass for a single disk crossing,

0 k (z) ≡ ∆E 0 k m χ = 1 2 (∆v z ) 2 = 2 g 2 z,disk (z = 0) z 2 V 2 z , (40)
which depends on the squared vertical coordinate z relative to the subhalo center.

A key assumption in deriving the previous results is that δZ does not vary significantly as the subhalo crosses the disk. This is very likely not verified for the innermost orbits, nor for the smallest objects, for which the impulsive approximation readily breaks down. Indeed, had subhalo particles enough time to circulate several times about the center as the object crosses the disk, conservation of angular momentum would prevent them to leave the system, and disk shocking would become inefficient. This is an example of the manifestation of adiabatic invariance, which was extensively studied in the context of stellar clusters in Refs. [126][127][128][129][130], from both analytic and numerical calculations. Following Ref. [126], capturing the results derived in Ref. [128] from the linear theory approximation, we introduce an adiabatic correction,

A(η) = (1 + η 2 ) -3/2 , ( 41
)
where η is the so-called adiabatic parameter, with η 1 for orbits close to the object's center, and η 1 close to the tidal radius. This gives A(η 1) → 0, and A(η 1) → 1, the latter case corresponding to the the parameter space for which the impulsive approximation holds. The adiabatic parameter is formally defined as

η(r, R) ≡ ω(r) τ (R) , ( 42
)
where ω is the orbital frequency that can be estimated from the inner dispersion velocity, ω = v 2 (r)/r, r being the distance to the subhalo center, and τ is the effective crossing time. The latter is given in terms of the half-height H of the disk, and of the vertical component of the subhalo velocity V z (R) at radius R in the Galactic frame. In the following, we will make use of the isothermal approximation, such that each Cartesian component of the velocity dispersion, for any system of mass m(r) inside a radius r, is related to the circular velocity according to

σ 2 v,i (r) = 1 2 v 2 c (r) = 1 2 G N m(r) r . (43)
Consequently, we get

ω(r) ≡ 3 G N m(r) 2 r 3 (44)
≈ 3.5 × 10 -2 Myr -1 m/6 × 10 -8 M (r/3.5 pc) 3 and

τ (R) ≡ H V z (R) = H 2 R G N M (R) (45)
≈ 0.45 Myr (H/100 pc) (V z /200 km/s) ,

where m(r) stands for the subhalo mass inside a radius r, while M (R) is the total Galactic mass inside a radius R. The orbital frequency is indicated for the typical mass a subhalo of m 200 = 10 -6 M has inside its scale radius r s ≈ 3.5 × 10 -8 M , taking a median concentration. This shows that except in the very central parts of subhalos where A(η) → 0, we will essentially have A(η) ∼ 1, corresponding to a maximal efficiency for disk shocking. Nevertheless, since m(r s )/r 3 s ∝ c 3 , we see that this efficiency will decrease as the concentration increases, protecting the most concentrated objects from disk-shocking effects. Actually, for a flat Galactic velocity curve of ∼ 200 km/s, we find assuming an NFW profile that to get η > 1, condition for the disk-shocking efficiency to start to be damped out, one needs x = r/r s 10 -3 c, regardless of the subhalo mass.

The adiabatic correction A(η) allows to modify the kinetic energy transfer defined in Eq. ( 46) in such a way that it is now valid over the full extent of any considered subhalo. This reads

k (z) ≡ 2 g 2 z,disk (z = 0) z 2 V 2 z A(η) , (46)
where the vertical subhalo velocity component V z (R) has been implicitly defined in Eq. (45). Finally, assuming circular orbits for WIMPs in a subhalo, one can easily express the average kinetic energy gain as a function of the radius r only, as

z 2 = (1/2) d cos θ r 2 cos 2 θ = r 2 /3. We get k (r) = 2 g 2 z,disk (z = 0) r 2 3 V 2 z (R) A(η) . (47)
The scaling with r is explicit, but for the quasiexponential suppression when r → 0 due to adiabatic 47) (here, we use the M11 mass model parameters). Also shown are their gravitational potentials. We plot the results in units of GN ρs r 2 s , which comes out as a natural scaling, to make the comparison more striking.

invariants: the gain in kinetic energy increases like the squared radius, and is maximal close to the tidal boundary of the subhalo. This scaling is shown as red curves in Fig. 1 for two different subhalo masses, 10 -6 (solid curve) and 10 6 M (dashed curve), and further compared to the moduli of their gravitational potentials, defined in Eq. (27).

The calculations presented above are at the basis of the methods we propose to follow to account for disk shocking, and thereby to further prune or destroy subhalos. Below, we discuss two different strategies, that we will call differential and integrated disk shocking to make the distinction clear. Common to both methods is the number of disk crossings, N cross , which is computed from the circular velocity of a subhalo in the Galactic frame (we implicitly assume circular orbits) and the age of the Galaxy T MW :

N cross (R) = G N M (R) R T MW π R . ( 48
)
Throughout this paper, we will set T MW = 10 Gyr. Beside the disk-shocking methods presented below, which are aimed to determine subhalo tidal radii in Galactic regions encompassing the baryonic disk, our tidal disruption criteria will be discussed in Sect. II D 4 where the disk-shocking methods will be further compared.

a. Tidal radius from differential disk shocking

The so-called differential disk shocking method will be our primary one, and relies on a comparison between the kick in velocity induced by disk shocking, as effectively described in Eq. ( 47), and the escape velocity,

v esc (r) = -2 φ(r) . ( 49
)
If the kick induced by disk shocking is such that the particle reaches the escape velocity, then it gets unbound to the system. Therefore, for each disk crossing, we will accordingly define the tidal radius as the radius at which the kick in velocity equals the escape velocity. In terms of energies, this reads:

r t such that k (r t ) = -φ(r t ) . (50)
This procedure must be applied at each crossing, such that it may somehow capture the dynamics of disk shocking. Indeed, hidden in φ [see Eq. ( 27)] is the radial boundary of the subhalo, which means that the above equation must be applied iteratively up to the number N cross of disk crossings given in Eq. ( 48). More explicitly, we have for the i th crossing

r t,i such that k (r t,i ) = -φ(r t,i , r t,i-1 ) . (51)
In practice, we start with the tidal radius inferred from the global tidal effects induced by the host halo and discussed in Sect. II D 1. This method can easily be applied to any subhalo model, irrespective of the inner density profile. It also provides a dynamical description of disk shocking, while only approximately. Indeed, this iterative procedure assumes that the internal structure of the shocked subhalo is not altered between two crossings, while part of the energy could actually be redistributed. Anyway, this picture is still consistent with adiabatic invariance, which partly protects the inner parts of subhalos against tidal pruning. An illustration of this differential disk shocking method is shown in Fig. 1, where we have plotted the diskshocking energy k (r) (red curves) and the gravitational potential modulus | φ(r)| as a function of the scaling variable r/r s (r s is the subhalo scale radius). We have considered two different NFW subhalos, 10 -6 (solid curves) and 10 6 M (dashed curves, respectively), both located at R . The corresponding gravitational potential moduli evaluated using two different radial boundaries for subhalos, one set to r 200 (blue curves), the other set to r 200 /10 (green curves), about which they are exponentially suppressed -the 1/r scaling expected beyond r s is poorly seen as the potential goes from ∝ cst to ∝ (1/r -1/r t ) ≈ e -r/rt /r very fast. These radial boundaries can be thought of as initial tidal radii before disk crossing. By virtue of Eq. ( 50), the tidal radius after one disk crossing will be set to the radius at which the kinetic energy and the potential curves intersect. Therefore, Fig. 1 nicely illustrates why the tidal stripping efficiency is much larger (i) in the outer regions of the system, and (ii) for more massive subhalos. This is due to 52) for the theory-based calculation, and Eq. ( 54) for the simulation-based estimate. Also shown are the binding and potential energies. We plot the results in units of GN ρ 2 s r 5 s , which leads to the same binding and potential energies for both subhalo masses.

the different scaling in r/r s , which is much sharper for the disk-shocking energy than for the gravitational potential. Not only do these results follow expectations, but they also allow to make specific calculations related to the phenomenology of subhalos.

b. Tidal radius from integrated disk shocking In contrast to what was presented above as a differential disk-shocking method, we can now try to integrate the kinetic energy gain over the whole subhalo -so the denomination integrated disk shocking method. Such a method was partly followed in Ref. [114], where the authors used the Eddington equation in the isothermal limit to convert the energy gain in phase space into a mass loss. Here, instead, we will use spherical symmetry, and simply assume that WIMPs take only circular orbits, such that the integrated kinetic energy gain can be expressed as

E k (r t , R) = 2 π rt 0 dr r 2 1 -1 d cos(θ) k (z, R) ρ(r) m χ , (52)
where k (z, R) is given by Eq. ( 46), and ρ is the inner subhalo mass density profile. Spherical symmetry merely implies that z 2 = r 2 cos 2 (θ), which makes the computation easy. This integrated energy gain can then be com-pared to the binding energy or to the potential energy, for each subhalo. For an NFW profile, the scaling goes from roughly ∝ r 4 t for r t r s , to ∝ r 2 t for r t r s . This is illustrated for a single disk crossing in Fig. 2, for two subhalos of 10 -6 (solid curves) and 10 6 M (dashdotted curves) located about the solar position -for completeness, we use two different concentration values for each subhalo: the median value (thin curves), and twice it (thick curves). The red curves show the integrated disk-shocking energy given in Eq. ( 52), as compared to the binding (green curve) and potential (blue curve) energies [see Eqs. ( 25) and ( 26), and the associated comments about the radial scaling]. Using units of G N ρ 2 s r 5 s allows to get a single curve for each of the latter energies, for both subhalo masses. The scaling of the disk-shocking energy with r t /r s is again sharper than that of the potential or binding energy, as expected. This implies that the central subhalo regions will be less prone to tidal stripping. We also see that as in the case of differential disk shocking, more massive subhalos will be more efficiently affected by tidal stripping.

Since we are dealing with integrated energies, we define the subhalo tidal radius after N cross disk crossings as

r t such that N cross E k (r t , R) = E b (r t ) , (53)
where we use the binding energy E b defined in Eq. ( 25) as a reference. c. Tidal radius from integrated disk shocking (fits on cosmological simulations) For the sake of comparison, we now introduce a result fitted on dark matter-only zoomed-in cosmological simulations, given in Ref. [116], wherein a baryonic disk potential was grown adiabatically to study the induced tidal disruption of subhalos. The qualitative features of this result were recently recovered in cosmological simulations including baryons, and discussed in Ref. [61]. The authors of Ref. [116] have tried to capture disk shocking effects by a simple and physically motivated ansatz, which, as they found, match rather well with their simulation results (see e.g. Ref. [START_REF] Binney | Galactic Dynamics: Second Edition[END_REF] for the dynamical grounds). They introduced an integrated-like disk-shocking energy Ẽk (r t , R) given by Ẽk (r t , R)

E b (r t ) = (1.84 r 1/2 ) 2 g 2 z,disk 3 σ2 v V 2 z , (54)
where r 1/2 is the radius containing half the subhalo mass, g z,disk is the disk gravitational force field given in Eq. ( 38), σ is an estimate of the internal dispersion velocity given by σ2 = 0.4 G N m t /r 1/2 , and V z is the velocity component perpendicular to the disk, that will be inferred from the approximation given in Eq. ( 43). This disk-shocking energy is shown as purple curves in Fig. 2, for the two subhalo prototypes introduced above. It still scales more sharply with r t (readily inferred as ∝ r 3 t from the equation just above) than the potential or binding energy, though less sharply than the integrated disk-shocking energy discussed in the previous paragraph while with similar amplitude around r t /r s 10. This means that this way to implement disk shocking will likely disrupt subhalos more efficiently, as gravitational stripping toward the central regions becomes more efficient. Still, we note that Eq. ( 54) relies on fits on simulation results, and could therefore be more specific to the subhalo mass range probed by cosmological simulations, which is still strongly limited by resolution issues. Anyway, the resulting subhalo tidal radius after N cross disk crossings can then be calculated by means of Eq. ( 53), merely by replacing Ẽk (r t , R) by Ẽk (r t , R). d. Disk-shocking summary We have introduced the so-called differential and integrated disk-shocking energies. For the latter, we have derived two expressions, one consistent with the differential one, another inspired by cosmological simulation and fully independent. These physical quantities allow us to derive the subhalo tidal radius r t after N cross disk crossings for any method. These calculations lead to different results, but common to all is the fact that r t does depend either on the subhalo mass m, its concentration c, its position in the Galaxy R, and its internal density profile. Our primary method will be the one based on the differential disk-shocking energy, as it relies on fewer assumptions. We will compare all these results in Sect. II D 4.

Subhalo mass-independence of xt = rt/rs

A striking property of all the tidal radius calculation methods discussed above, either those involving global tides and those involving disk shocking, is that the ratio x t = r t /r s turns out to be independent of the subhalo mass. Actually, x t depends only on the subhalo concentration c and on its radial position R in the Galaxy. If the latter dependence is rather easy to understand (tidal stripping depends on the position), the former is much less trivial.

For the global tides discussed in Sect. II D 1, it is easy to show that the methods based on the Jacobi limit can be formulated along

x t = ∆ 200 f (x t ) ∆ t (R) f (c) 1/3 κ c ⇔ x t [f (x t )] -1/3 = ∆ 200 ∆ t (R) f (c) 1/3 κ c ,
which makes it clear that x t is only a function of R and c. Here, κ = r -2 /rs is set by the choice of the inner profile, and function f can be defined on general grounds by means of the subhalo mass, m(x) = 4 π ρ s r 3 s f (x), where x = r/r s -for an NFW profile, it is simply f (x) = ln(1+x)-x/(1+x). We have also defined ∆ x = ρ rx /ρ c , i.e. the ratio of the average subhalo density within a radius r x to the critical density (∆ 200 = 200). In the case of the point-like Jacobi approximation corresponding to the tidal definition of Eq. ( 34), for instance, we have

∆ t (R) = 9 M/(4 π R 3 ρ c ) ,
where M is the whole host galaxy mass.

The demonstration for the method setting the tidal radius from equating the inner density to the outer density, given in Eq. ( 36), is trivial, and relies on the fact that the subhalo scale density ρ s , regardless of its profile and its mass, is only set by the concentration parameterfor an NFW profile, it reads

ρ s = ∆ 200 ρ c 3 c 3 f (c) .
If we write the density profile as ρ(r) = ρ(x = r/r s ) = ρ s u(x), then Eq. ( 36) translates into u(x) = ρ tot (R)/ρ s , which makes it clear again that x t depends only on c and R.

Finally, the cases of disk-shocking tidal effects are more subtle. In the differential method, x t can readily be shown to be a function of c and R only from Eq. ( 50). This is simply because the potential φ(r t ) ∝ r 2 s g(x t , c), where it is not necessary to specify function g, while the kinetic energy k (r t ) ∝ r 2 t ḡ(x t , c, R), function ḡ being unspecified too, such that equating them leads to an equation that involves only variables x t , c and R. This proves that x t only depends only on c and R. The reasoning is similar for the so-called integrated disk-shocking methods, and also leads to the dependence only on R and c of the associated x t .

Tidal disruption criterion and minimal concentration

Equipped with several tidal radius definitions, we can now define a tidal disruption criterion by specifying the function ζ(r t ) introduced in Eq. (23), where r t is the subhalo tidal radius. We remind that the latter depends on all the specific subhalo properties, and on its position in the host halo. In light of results obtained in Ref. [76], we may define the following very simple disruption function,

ζ x t ≡ r t r s ≡ θ (x t -ε t ) , ( 55
)
where θ is the usual dimensionless step function, r s = r s (m, c) is the subhalo scale radius, and parameter ε t sets the minimal value allowed for x t . This parameter very likely depends on the inner subhalo density profile, and could also depend on the specific process responsible for tidal stripping. Typical values found using dark matter-only simulations are ε t ≈ 2 (see Ref. [76]), but we may wonder whether simulations can efficiently capture the continuous limit due to their limited spatial/mass resolution. For definiteness, we will set ε t = 1 in the following, unless specified otherwise. This translates into a minimal bound on the subhalo concentration, c min (R), as the surviving subhalos are only those with scale radii such that r t /r s ≥ ε t . This concentration cut-off reads a transcendental equation that can be solved iteratively.

c min (R) = ε t κ r 200 (m) r t (c min (R), m, R) , (56)
Here, κ = r -2 /r s is fixed by the choice of density profile (κ = 1 for an NFW or an Einasto profile). In practice, we will further impose that

c min (R) = Max {c min (R); 1} . ( 57
)
We emphasize that c min does actually not depend on the subhalo mass, only on its location R in the Galaxy. This is because x t is only a function of the concentration c and R, as explained in Sect. II D 3. This concentration lower bound, c min (R), is the very variable that differentiates the tidal stripping methods discussed in Sect. II D 1 and in Sect. II D 2. We report our calculations of c min in Fig. 3, as a function of the dimensionless Galactic radius R/R 200 (R 200 = 237 kpc in the M11 model). The curves related to the global tides are shown as solid colored lines, while those associated with disk shocking are the non-solid ones. Note that we have also included the baryons in the calculation of the global-tide effects (see Sect. II D 1). We also stress that we did the calculation assuming two different inner subhalo profiles: an NFW profile (thick curves), and an Einasto profile (thin curves) -we took an index of α e = 0.17 for the latter. From the plot, there is no significant qualitative difference between these profiles, but that Einasto subhalos are very slightly more resistant to gravitational tides.

We note that the most approximate method for the global tides, the point-like Jacobi limit given in Eq. (34), is also the one that destroys subhalos most efficiently, even more efficiently than disk shocking in the central parts of the Galaxy. It can therefore be used for fast and conservative calculations, while it is highly sensitive to the estimate of the total mass of the Galaxy, which is often ambiguous as it depends on the choice for the virial radius. To make the discussion more quantitative, we recall that a 10 -6 M subhalo has a peak concentration of ∼ 60, which will serve as a reference value here. We see from the plot that the point-like tide method already affects such tiny objects already from 20 kpc and pushes them to exponentially high concentration inward. This means that at the solar position these objects have already been almost fully disrupted. The two other globaltides methods [given in Eqs. ( 35) and ( 36)], much more realistic, give similar results and lead to much less tidal stripping than the point-like approximation. Subhalos of 10 -6 M start to be strongly affected around 2-4 kpc from the MW center in these scenarios.

Disk-shocking effects start to play a role only from 20 kpc inward, as expected from the typical gravitational size of the Galactic disk. All disk-shocking methods lead to more stripping than global tidal effects, except for the point-like approximation discussed above. Here again, we see that the most approximate method, the integrated disk-shocking method fitted on cosmological simulations and given in Eq. ( 54), is the most efficient to destroy or prune subhalos. Besides being based on very crude approximations, we stress that it is also likely biased by the resolution limit inherent to cosmological simulations, where only subhalos with masses 10 4-7 can be tracked. These massive objects are much less concentrated than their lighter brothers and sisters, and more prone to stripping and disruption. In contrast, the less efficient method is the one based on integrated disk shocking and given in Eq. (53). Intermediate is the method the most motivated on theoretical grounds. Interestingly, the latter starts to deplete subhalos of 10 -6 M around the position of the Sun.

In summary, global tides tend to dominate the stripping beyond the disk, while disk shocking dominates inward. This was obviously expected, but we quantified and illustrated these effects rather exhaustively. Moreover, we showed that the point-like Jacobi approximation makes it irrelevant to include disk shocking, as it supersedes all over effects over the whole Galactic range. Nevertheless, as we discussed above, this point-like approximation is by far the worst to make, while being conservative. Obviously, in a consistent and complete model, one has to include all tides, those coming from global gravitational effects, and those coming from disk shocking. This is what we will do when discussing our final results in Sect. III.

Tidal selection of the most concentrated objects: shift of the average concentration

By depleting the lower tail of the concentration distribution, tidal effects modify the average concentration of subhalos as a function of their mass. This can be explicitly calculated by means of the first moment of the concentration function, given in Eq. (30). The increase in the average concentration merely comes from that tidal effects reduce the concentration range from below by c in (R) ≥ 1. In reality, the concentration function should not be truncated that sharply, but this truncation still captures the main physical effects at play.

We illustrate this in Fig. 4, where we report our calculations of c(m) as a function of the dimensionless Galactic radius X 200 = R/R 200 for three different subhalo masses, 10 -6 , 1, and 10 6 M , and for all the tidalstripping methods introduced above. The asymptotic values of c(m) at X 200 → 1 correspond to the average concentration computed in the range c min = 1, c max = ∞. As we go inward, tidal effects come into play and c min increases, leading to the increase in c(m) . Recalling that the concentration function is Gaussianly suppressed beyond the median value c 0 (m) ≈10-100, in the considered subhalo mass range, we can therefore read off from the plot that most of subhalos with masses larger than that of a given curve are tidally depleted as the curve exceeds ∼ 100. This trend is consistent with previ-ous studies performed from dark matter-only cosmological simulation results (see e.g. [28,76,88,119,121,131]), or from simple analytic approximations (see e.g. [118]), but these works did not include baryonic effects. Here we provide quantitative estimates for both baryonic and dark matter tidal effects, and comparisons between different approaches.

Impact of tidal effects on the calibration and normalization procedure

It may turn useful to summarize the way tidal effects are integrated in the full procedure in practice. As discussed in Sect. II B 3, we calibrate the subhalo population by considering only the so-called global tidal effects presented in Sect. II D 1. These global tidal effects translate into a function c 0 min (R) that cuts the concentration pdf from below, and allows us to determine both N 0 sub and the associated normalization of the whole subhalo phase space K 0 . This must be done without baryons at all, consistently with the fact that the calibration is based upon dark matter-only simulation results. Then, we compute the final phase-space normalization K that accounts for the baryonic tides (both the global-tide and the disk-shocking calculations), which are characterized by a new cut-off function c min (R). We obtain the final number N sub of subhalos by demanding that the overall subhalo mass density is unaffected at very large radii, far from the disk, where baryonic effects can be neglected. This can be rephrased as setting N sub = (K/K 0 ) N 0 sub .

III. CONCRETE RESULTS: MASS PROFILES, LUMINOSITY PROFILES, AND BOOST FACTORS

In the previous section, we have introduced the whole scheme to derive a Galactic subhalo population consistently with current dynamical constraints, assuming only a smooth halo component and spherical symmetry. This scheme was integrated in a C++ numerical code, and in this section, we will present our main results. In Sect. III A, we will first inspect the overall subhalo mass profile that comes out as a result of our tidal-stripping procedure. Then, in Sect. III B, we will show how the obtained subhalo distribution translates into annihilation profile, and subsequently quantify the associated annihilation boost factor. Throughout this section, we will also comment on the specific impacts of the initial subhalo mass distribution index α and the minimal subhalo mass (see Sect. II C 2).

A. Mass profiles

In Sect. II A, we have introduced the smooth and subhalo dark matter components from the overall density profile ρ tot in Eq. ( 2), that we rewrite with an explicit spherical symmetry to make the present discussion clear:

ρ tot (R) = ρ sm (R) + ρ sub (R) . (58)
We remind that ρ tot is subject to dynamical constraints, and we have adopted the M11 model as a template Galactic mass model (see Sect. II B). The smooth dark matter component ρ sm featuring above can only be derived a posteriori, after having determined the subhalo component ρ sub given in Eq. (11). Making the tidal cut-off c min (R) explicit, the latter reads:

ρ sub (R) = N sub K w dP V (R) dV mmax mmin dm cmax cmin(R) dc m t (r t (c, m, R), m, c) dP m dm dP c dc with K w = 4 π R200 0 dR R 2 dP V (R) dV mmax mmin dm cmax cmin(R) dc dP m dm dP c dc , (59)
where N sub is the total number of subhalos, m t is the subhalo mass contained in the tidal radius r t , and the P's define the global subhalo phase space, normalized to unity thanks to K w , and were introduced in Sect. II B 2.

Taking the M11 Galactic mass model as a reference, the associated prediction of the overall subhalo mass density profile is shown in Fig. 5, where we also represent the impact of the mass index (α m = 1.9/2 in the top/bottom row panels) and of the minimal subhalo mass (10 -6 /10 -10 M in the left/right column panels). In each panel, we give predictions for all the methods introduced in Sect. II D to compute the tidal stripping and associated subhalo disruption -the upper part of each plot shows the mass density profile, and the lower part the subhalo mass fraction, as functions of the dimensionless Galactic radius R/R 200 . Subhalo mass profiles relying on global tides only are reported as solid colored curves -see Sect. II D 1 -while those incorporating also diskshocking effects are shown as dashed (differential), dashdotted (integrated), and dotted (simulation-fit inspired method, respectively) -see Sect. II D 2. We recall that our reference model is based on the global tides evaluated for a smooth halo in the Jacobi limit (dubbed smooth host in the plots), given in Eq. ( 35), and on the differential disk-shocking tides given Eq. ( 51). It is represented as dashed red curves in each panel. Some illustrative FIG. 5: . numbers can also be found in Tab. III.

Overall, these results show that tidal effects strongly deplete the subhalo population in the central parts of the Galaxy, and underline the effects of disk shocking, which plays an important role. This leads to a coredlike spatial distribution inward, before the full disruption of subhalos in the very center ( 4 kpc). These are generic features observed in cosmological simulations, but our analytic procedure allows us to make predictions down to much lower spatial and mass scales, in a dynamically constrained and consistent frame. Going to more specific global tidal stripping configurations, we see that the global point-like Jacobi method, which is clearly too approximate as it does not account for the host halo and subhalo profiles details, disrupts almost all subhalos within R/R 200 0.1 ⇔ R 20 kpc, making disk-shocking effects even irrelevant. This leads to a negligible local (R/R 200 ∼ 0.3-0.4) subhalo mass fraction, typically 1%. The two other more physically motivated global methods provide slightly more optimistic Results for our subhalo model as embedded in the M11 Galactic configuration, when all tidal effects are included (global smooth Jacobi limit, and differential disk-shocking method). Here, we take a cut-off subhalo mass of 10 -10 M , and a tidal disruption efficiency of ε t = 1.

predictions, with a subhalo mass fraction ∼ 10%. When disk-shocking effects are further included, however, it decreases down to 1%. Still, we will see in Sect. III B that this low fraction is somewhat compensated for, in terms of annihilation rate, by a tidal selection of more concentrated objects. The impact of the minimal subhalo mass is only noticeable for a mass index of 2, as most of the mass fraction is then carried by the smallest objects, which are much more resilient to tidal effects. Since the minimal subhalo mass can in principle be determined by the interaction properties of WIMPs (more or less straightforwardly related to its mass), most of the theoretical uncertainties is then featured by α m . For or reference model, we see that while the mass fraction can vary by a factor of ∼ 2 between α m = 1.9 and α m = 2 in the outskirts of the Galaxy, its differential value is much more sensitive because of more efficient tidal selection of lighter subhalos. In the latter case, the variation can reach an order of magnitude. This should have impact on predictions for direct subhalo searches (see e.g. Refs. [132,133]).

Finally, we note that the amount of subhalo mass lost during disk crossings could in principle be quantified from our method, which may be used to size the impact of the smallest pruned subhalos on the high tail of the WIMP velocity distribution (see e.g. Refs. [134,135]). Indeed, disk shocking induces a net kinetic energy gain for the pruned WIMPs. This, however, goes beyond the scope of this work.

B. Annihilation rate profiles and boost factors

In this paragraph, we discuss the potential enhancement the presence of dark matter subhalos may induce in the WIMP annihilation rate in the Galaxy, usually dubbed as boost factor. Here, we will only determine the differential and integrated boost factors on the annihilation rate, not on the observable cosmic-ray or gammaray fluxes. We remind the reader that in terms of these fluxes, the annihilation boost factor is angular-dependent for gamma rays [32], while energy-dependent for antimatter cosmic rays [34,35].

In this work, we will assume that subhalos do not superimpose, such that we will not account for the potential existence of sub-subhalos, which might be relevant in the most massive subhalos as they have formed at later epochs than the lightest ones. We still stress that any inclusion of sub-subhalos should be consistent with the normalization and calibration procedures one subscribes to (in particular, the overall mass function should be recovered after all layers of subhalos have been accounted for).

With this assumption, the total annihilation rate can be derived starting from a discrete distribution of N sub subhalos, where the total dark matter density would be given by

ρ tot,N sub ( x) = ρ sm ( x) + N sub i M t,i δ( x -x i ) , (60)
where each subhalo i is point-like and allocated a tidal mass M t,i , and ρ sm is the smooth dark matter density. Squaring this equation, we get

ρ 2 tot,N sub ( x) = ρ 2 sm ( x) + N sub i≥j M t,i δ( x -x i ) M t,j δ( x -x j ) + 2 ρ sm ( x) N sub i M t,i δ( x -x i ) . (61)
Neglecting sub-subhalos formally implies that x i = x j ∀ i = j. If we now take the continuous limit and make spherical symmetry explicit, we obtain where m t and ξ t are the subhalo tidal mass and annihilation volume, defined in Eqs. ( 23) and ( 24) respectively. Symbol denotes the averaging over the concentration and mass parts of the subhalo phase space, made explicit in Eq. (11), which is position-dependent. Notice the crossing term above induced by the interaction between subhalos and the host halo. Usually assumed subleading and thereby neglected, it may actually dominate over the smooth contribution at large Galactic radii, as will be shown below.

ρ 2 tot (R) = ρ 2 sm (R) + N sub K w dP V (R) dV mmax mmin dm dP m (m) dm ∞ cmin(R) dc dP c (c, m) dc rt(R,c,m) 0 dr 4 π r 2 ρ 2 (r) + 2 ρ(r) ρ sm (R) = ρ 2 sm (R) + N sub dP V (R) dV ρ 2 ξ t (R) + 2 ρ sm (R) m t (R) (62
We now define the dimensionless WIMP luminosity L, which measures the spatial dependence of the annihilation rate, as

L(R) ≡ ρ 2 tot (R) ρ 2 , ( 63
)
where the normalization ρ is made at the solar position. We further introduce the differential annihilation boost factor B(R), and the integrated annihilation boost factor B(R), as

B(R) ≡ L(R) L no sub. (R) , (64)
B(R) = R 0 dr r 2 L(r) R 0 dr r 2 L no sub. (r)
.

Defined so, these boost factors are merely the multiplicative corrections to apply to the differential or integrated annihilation rate as computed by neglecting the subhalo component. Note that in principle, any WIMP signal prediction involving subhalos should be affected by a statistical variance, reflecting the possible fluctuations of the number of contributing objects [34,35]. We keep this aspect for further dedicated studies. We report our calculation results for the annihilation profiles in Fig. 6, where we again adopt M11 as the reference Galactic mass model -this is the translation of Fig. 5 in terms of annihilation profiles. Global tides are calculated from the smooth Jacobi method [see Eq. ( 35)], while disk-shocking effects are described from the differential method [see Eq. ( 51)]. Top (bottom) row panels correspond to a subhalo mass index of 1.9 (2, respectively). Left (right) column panels correspond to a minimal subhalo mass of 10 -6 M (10 -10 M , respectively). In each panel, the upper part shows the different components of the annihilation profile, and the lower part shows the differential and integrated boost factors as defined above. We display the impact of neglecting disk shocking as dashed curves, which demonstrate the importance of this effect in the central parts of the Galaxy.

A generic result is that the subhalo contribution dominates at large radii, typically from the edge of the Galactic disk for a subhalo mass index of α m = 1.9, even from much more inner regions in the case of α m = 2. For an overall NFW profile, this leads to a characteristic scaling of 1/r 2 toward the Galactic center, where the smooth halo dominates, progressively changing to 1/r 3 outward, when the luminosity profile tracks the subhalo spatial distribution. Interestingly, in the case of α m = 2, where the global subhalo luminosity is enhanced, a plateau arises in the overall luminosity profile at the transition between smooth-halo domination and subhalo domination. This is actually an imprint of disk-shocking effects, which delay the rising of the subhalo contribution. We will see later that this plateau does not depend on the tidal disruption efficiency, and is also preserved in the case of an overall Einasto profile. This striking feature might be used in gamma-ray searches.

Finally, we comment on our results for the differential and integrated boost factors. The so-called differential boost factor (reported as "local" in the plots) is mostly relevant to indirect dark matter searches with antimatter cosmic rays because of the limited horizon of the latter induced by propagation effects. It also represents the correction to apply to the integrand of the line-of-sight integral used in gamma-ray searches. On the other hand, the integrated boost is related to the absolute Galactic luminosity, and thereby to extragalactic gamma-ray searches -then the Galaxy appears as a template case characterizing other galaxies close in mass. We see that at the solar position (R/R 200 ∼ 0.03-0.04), the boost is locally < 2 for α m = 1.9, while it reaches ∼ 5 for α m = 2. Though moderate, these values may have some impact on the existing limits on WIMPs as the precision in the cosmic-ray data has strongly increased in the recent years. We also remark that the differential boost increases up to 10 3 -10 4 toward the edge of the Galaxy, which strongly affects, for instance, the diffuse gamma-ray signal on high Galactic latitudes, as known from long ago (see e.g. Refs. [32,36]). Regarding the integrated boost, the values obtained at the edge of the dark halo can represent useful calibration values for calculations of the dark matter contribution to the extra-galactic diffuse gamma-ray. These go from ∼ 3 for α m = 1.9, to ∼ 20 for α m = 2. This is fully consistent with the recent study in Ref. [121] (see Fig. 6 in this article), which is based on fits of cosmological simulations, and does not include baryonic effects -while global tides are merely the outcomes of the simulations themselves. That baryons play no role for the integrated boost at the whole Galactic scale should not come as a surprise, as they only affect the dynamics in the very central parts of the halo.

At this stage, we illustrated our results assuming a tidal efficiency of ε t = 1 [see Eq. ( 55)]. It is important to check their stability against changes in this parameter. In Fig. 7, we investigate the impact of ε t by computing the annihilation profiles for ε t = 0.5 (subhalos can be pruned down to r s /2 before getting disrupted), and for ε t = 0.5 (subhalos can be pruned only down to 2 × r s before getting disrupted). We adopt the configuration for which the plateau discussed above is visible, namely α m = 2. We see that the plateau is slightly smeared when ε t = 0.5, as disk-shocking effects are then less disruptive and smear the previously abrupt rising of the subhalo contribution. On the contrary, the plateau is much more salient when ε t = 2, as expected. Values of 0.5 are rather small compared to what is found in cosmological simulations (see e.g. Ref. [76]), but could be relevant, for instance, to the very concentrated cores of ultra-compact mini-halos, close to the minimal cut-off mass (see e.g. Ref. [26] for further discussion). Nevertheless, we see that except close to the smooth-halo/subhalo luminosity domination transition, changes in ε t have no significant impact on our predictions. In particular, the plateau feature does not seem to be spoiled.

C. Comparison between different Galactic mass models

Throughout this study, we have adopted the M11 Galactic mass model as a reference. It is interesting to check how our predictions are affected by changes in the Galactic model itself, while still trying to use dynamically constrained scenarios. To this aim, we will use two other Galactic models: (i) the one of Ref. [46] (CU10 hereafter), which is particularly interesting as it relies on an Einasto dark halo instead of an NFW one, and (ii) the upgraded and updated version of M11, recently released in Ref. [68] by the same author (M16 henceforth), still based on an NFW profile for the dark halo, but with additional constraints from novel kinematic data. Beside the details of the dark matter halo profile and the data sets used as dynamical constraints, changes with respect to M11 also come from the additional inclusion of a twocomponent gaseous disk (HI and HII) in both models. Including gas has actually no significant impact on the dark matter profile on large scales, as this mostly tunes the distribution of the inner overall gravitational potential among all, baryonic and dark, components. Rather, this mostly influences the understanding of the very local stellar dynamics. Note that since the gas components of CU10 are directly inferred from data points, we will instead, for convenience, use that of M16 in both M16 and CU10, which may be described from Eq. ( 6) (we further trade the original vertical sech 2 of M16 for an exponential function, which has no impact on the final result). The sets of parameters of these models are given in Tab. I and Tab. II.

Concrete comparisons are displayed in Fig. 6 in terms of annihilation and boost-factor profiles, where solid (dashed and dotted) curves correspond to M11 (CU10 and M16, respectively) predictions. We adopt the "bestcase" configuration where α m = 2 and the cut-off subhalo mass is 10 -10 M . The color code is the same as in Fig. 6. Luminosity profiles are measured in units of the squared dark matter density ρ 2 at the solar position R , which (barely) change from one model to another. The upper horizontal axes feature the Galactic radius in units of R 200 , which also varies between configurations -see Tab. I. We add another lower panel that provides the real annihilation profile ratio with respect to M11, where the luminosity is then evaluated at the corresponding M11 radius for each model -consequently, the lower horizontal axis features R/R 200 as inferred from M11 only.

We first notice that the difference between M11 and M16 is hardly visible, and amounts to 10% over the full Galaxy, M16 being slightly less luminous in terms of dark matter annihilation. This is not surprising as the only changes between M11 and M16 are the addition of a gaseous disk and a new set of constraining data. On the other hand, differences are more pronounced between CU10 and M11-M16: CU10 is brighter than M11-M16 in the central parts of the Galaxy, typically for R/R 200 < 0.1. This is obviously a consequence of the different halo shape, as Einasto profiles are known to be more luminous than NFW profiles within the scale radius (but for the divergence of NFW profiles at the very centers -see e.g. Ref. [36]), while having a faster luminosity decrease outward because of the exponential cut-off in the halo shape. This amounts to an increase of ∼ 20% in integrated luminosity at the solar position, and a decrease of the same order at the edge of the Galaxy. In terms of boost factors, which measure the impact of subhalos relative to the host halo and therefore can be directly compared between different mass models, we see that the difference is very moderate for the differential boost, leading to an integrated difference < 2 at the edge of the dark halo (CU10 leads to a slightly smaller integrated boost factor).

Finally, we remark that the plateau feature emphasized in Sect. III B as a signature of a sharp subhalo mass function also shows up in CU10, despite the different overall halo shape. This prediction is therefore robust against systematic uncertainties in the overall dark halo modeling, provided the smooth halo density continues increasing inward, where the subhalo population has been fully depleted -this plateau could convert into a bump for a cored smooth halo profile.

IV. CONCLUSIONS

We have proposed a method to model a galactic dark matter subhalo population consistently with dynamical constraints, focusing on the Milky Way -subhalos are unavoidable galactic components if dark matter is made of WIMPs or any other dark matter candidates with suppressed self-interactions and devoid of additional pressure. Dynamical consistency is important to make sense of contraints or discovery potentials of both direct and indirect dark matter searches (see e.g. Ref. [136] for an illustration in direct searches). We have assumed that subhalos initially track the host halo profile when the galaxy forms and then we have explicitly calculated the effects of tidal stripping and subsequent potential disruption induced not only by the overall gravitational poten-tial but also by baryonic disk crossing. We have developed and compared different theoretical approaches to deal with the latter, and retained the so-called differential disk-shocking method as our reference case, since built upon more accurate physical grounds. This method was inspired by previous works dedicated to the understanding of stellar clusters, in particular by Refs. [125-127, 129, 130]. These works' results were already used in other analytic studies (see e.g. Ref. [114]), but were dealt with in a significantly different way, leading to a different formulation of tidal mass losses without explicit links to the definition of the tidal radius. Nevertheless, even if it is difficult to make quantitative comparisons between the mentioned study and ours, it seems that both approaches are in agreement at least at the qualitative level. Our study was more aimed to quantify the impact of a subhalo population on the dark matter annihilation rate and to size the related theoretical uncertainties in a realistic Galactic mass configuration, while Ref. [114] was more concerned with the survival probability of subhalos against different types of tidal effects. On the whole, our method to include disk-shocking effects is likely simpler to implement in numerical calculations.

The main inputs of our model are (i) the Galactic mass model, (ii) the subhalo mass function, and (iii) the subhalo concentration function, for which we adopted consensual prescriptions. Further assumptions regard the choice of the inner subhalo profile. We considered a spherically symmetric host halo, hence a spherically symmetric subhalo distribution. Our model can in principle easily be extended to axisymmetric host halos, while its numerical implementation will then likely become much trickier. We stress that we calibrated the subhalo mass fraction using constraints from cosmological simulations without baryons. This is important to use dark matteronly simulations because baryonic components in hydrodynamic simulation are likely to differ significantly from those of the real Milky Way, while strongly affecting the dynamics of subhalos: it would then become impossible to disentangle the global from other tidal effects, and it would make it spurious to calibrate a subhalo population model a posteriori. This underlines the need to still continue running dark matter-only simulations with increased resolution and up-to-date cosmological parameters, even in a context where issues related to the impact of baryons on cold dark matter halos are certainly the most pressing ones.

Using the recent and constrained Galactic mass models from Refs. [48,68,72] (dubbed M11, CU10, and M16 -M11 being used as the template case), characterized by different assumptions on the dark halo profile while providing results consistent with each other, we computed the overall subhalo mass profiles and further made predictions for the induced annihilation profiles. We stress that these results incorporate a self-consistent calculation of the tidal radius for each subhalo, depending on its mass, concentration, and position in the Galaxy; individual mass and luminosity are calculated up to the tidal radius for each subhalo. We used different assumptions for the mass index α m and the cut-off subhalo mass. Since the latter could in principle be determined from WIMP interaction properties in specific scenarios, the main the-oretical uncertainty is actually carried by α m . We showed that the global or integrated boost factor could vary between ∼ 2 (for α m = 1.9) and ∼ 20 (for α m = 2, respectively) for all choices of Galactic mass models. This may provide interesting and complementary calibration points for estimates of the dark matter contribution to the extragalactic diffuse gamma-ray background (see e.g. Ref. [39]). We also derived differential boost factors (i.e. the boost factor profile) that could be used to revisit estimates of the dark matter contribution to the Galactic diffuse gamma-ray emission (e.g. Refs. [137][138][139][140]), or to the local antimatter cosmic-ray flux (e.g. Refs. [13,[141][142][143]). Interestingly, our model predicts a plateau in the overall annihilation rate in the case of a sharp mass function (α m = 2) that could lead to specific observable effects. This feature seems to persist within the considered theoretical and systematic uncertainties of the model.

The local subhalo population and induced boost factor, relevant to direct searches and antimatter searches, respectively, are very sensitive to α m . For antimatter searches, though, the precision achieved in the most recent measurements is such that even a moderate effect could have significant impact on the existing limits or discovery prospects. In the most optimistic case, when α m = 2, the enhancement can reach a factor of 10.

It would be interesting to test this model against cosmological simulations with baryons in the relevant subhalo mass range in the near future, but this is clearly beyond the scope of the present study. In any case, the model is easily tunable in terms of initial distribution functions, provided internal consistency with the dynamical constraints, which was the main purpose of this work. Finally, self-made predictions for direct and indirect dark matter searches are left to further studies.

  dc m t (r t (c, m, x), m, c)

  in Ref.[60]), we getM VL2 200 = M VL2 (R VL2 200 = 225.44 kpc) = 1.42 × 10 12M . This allows us to propose the following ansatz to normalize the subhalo population:

∆ m 12 sub∆ m 12 sub

 1212 defined in Eq. (16) to the constraint f given in Eq. (20) (replacing the tidal mass m t by m 200 in the definition of f ∆ m 12 sub

  allows the normalization of the pdf to unity over the mass range delineated by [m min , m max ]. Note that we implicitly assume m = m 200 .

FIG. 1 :

 1 FIG. 1:Average kinetic energy gain induced by a single encounter with the Galactic disk as a function of r/rs, for two NFW subhalos (10 -6 and 10 6 M ) located at R -see the definition in Eq. (47) (here, we use the M11 mass model parameters). Also shown are their gravitational potentials. We plot the results in units of GN ρs r 2 s , which comes out as a natural scaling, to make the comparison more striking.

5 s 2 s r 5 s)FIG. 2 :

 5252 FIG. 2: Integrated kinetic energy gain induced by theGalactic disk as a function of rt/rs, for two NFW subhalos (10 -6 and 10 6 M ) located at 8 kpc, assuming two different concentration values (the median one, and twice it) -see the definitions in Eq. (52) for the theory-based calculation, and Eq. (54) for the simulation-based estimate. Also shown are the binding and potential energies. We plot the results in units of GN ρ 2 s r 5 s , which leads to the same binding and potential energies for both subhalo masses.

FIG. 3 :

 3 FIG. 3: Minimal concentration as a function of the dimensionless Galactocentric radius R/R200, induced by different tidal effects. The solid curves show the global-tides effects discussed in Sect. II D 1, and the non-solid curves show the disk-shocking effects discussed in Sect. II D 2. See more comments in Sect. II D 4.

FIG. 4 :

 4 FIG. 4: Minimal concentration as a function of the dimensionless Galactocentric radius R/R200, induced by different tidal effects. The solid curves show the global-tides effects discussed in Sect. II D 1, and the non-solid curves show the disk-shocking effects discussed in Sect. II D 2. See more comments in Sect. II D 4.

3 ] 3 ] 2 M

 332 )halo mass profile with tides Whole DM profile Sub (point-like host) Sub (smooth host) Sub (ρsub(rt) = ρtot(R)) )halo mass profile with tides Whole DM profile Sub (point-like host) Sub (smooth host) Sub (ρsub(rt) = ρtot(R)) Stref & Lavalle 2016 DM (sub)halo mass profile with tides Whole DM profile Sub (point-like host) Sub (smooth host) Sub (ρsub(rt) = ρtot(R)) Stref & Lavalle 2016 DM (sub)halo mass profile with tides Whole DM profile Sub (point-like host) Sub (smooth host) Sub (ρsub(rt) = ρtot(R)) min = 10 -10 M t = 1

  FIG.8:.

TABLE I :

 I Dark matter halo parameters for different Galactic mass models (best-fit models of Refs.[48] [M11],[72]

TABLE II :

 II Baryonic component parameters for different Galactic mass models (best-fit models of Refs.[48] [M11],[72] [CU10], and[68] [M16]

TABLE III :

 III

  Annihilation/luminosity profiles for different assumptions. Upper part of each panel: Black curves show the profile when neglecting subhalos; blue solid curves show the overall profile when subhalos are included and all tidal effects considered; green curves show then the separate contribution of the smooth halo, red curves the contribution of subhalos, magenta curves the contribution of the subhalo-smooth halo cross product; dashed lines show the impact of neglecting disk-shocking effects. Lower part of each panel: differential (blue curves) and integrated (red curves) boost factors. Upper/lower row: αm = 1.9/2. Left/right column: minimal subhalo mass of 10 -6 /10 -10 M .

		10 2			Stref & Lavalle 2016		10 2		Stref & Lavalle 2016
	10 -1 10 0 10 1 L/ρ 2 [dimensionless]			Annihilation profile Without subhalos With subhalos: smooth contrib. With subhalos: sub. contrib. With subhalos: cross product With subhalos: total contrib.	L/ρ 2 10 -1 10 0 10 1 [dimensionless]		Annihilation profile Without subhalos With subhalos: smooth contrib. With subhalos: sub. contrib. With subhalos: cross product With subhalos: total contrib.
	Annihilation profile		Tidal effects Without disk shocking With disk shocking Tidal effects Without disk shocking With disk shocking		Annihilation profile	10 -4 10 -3 10 -2	Tidal effects Without disk shocking With disk shocking Tidal effects Without disk shocking With disk shocking
		Annihilation Boost factor	Boost factor Local boost Integrated boost Boost factor Local boost Integrated boost	m cut =10 -6 M α M =1.9 t =1	10 0 10 1 10 2 10 3 10 4		10 -5 Annihilation Boost factor	Boost factor Local boost Integrated boost Boost factor Local boost Integrated boost	m cut =10 -10 M α M =1.9 t =1	10 4 10 0 10 1 10 2 10 3
		10 -2		10 -1 R/R 200	10 0 10 -1		10 -2		10 -1 R/R 200	10 0 10 -1
		10 2			Stref & Lavalle 2016		10 2		Stref & Lavalle 2016
	10 -1 10 0 10 1 L/ρ 2 [dimensionless]			Annihilation profile Without subhalos With subhalos: smooth contrib. With subhalos: sub. contrib. With subhalos: cross product With subhalos: total contrib.	L/ρ 2 10 -1 10 0 10 1 [dimensionless]		Annihilation profile Without subhalos With subhalos: smooth contrib. With subhalos: sub. contrib. With subhalos: cross product With subhalos: total contrib.
	Annihilation profile	10 -4 10 -3 10 -2	Tidal effects Without disk shocking With disk shocking Tidal effects Without disk shocking With disk shocking		Annihilation profile	10 -4 10 -3 10 -2	Tidal effects Without disk shocking With disk shocking Tidal effects Without disk shocking With disk shocking
		10 -5 Annihilation Boost factor	Boost factor Local boost Integrated boost Boost factor Local boost Integrated boost	m cut =10 -6 M α M =2 t =1	10 0 10 1 10 2 10 3 10 4		10 -5 Annihilation Boost factor	Boost factor Local boost Integrated boost Boost factor Local boost Integrated boost	m cut =10 -10 M α M =2 t =1	10 4 10 0 10 1 10 2 10 3
		10 -2		10 -1 R/R 200	10 0 10 -1		10 -2		10 -1 R/R 200	10 0 10 -1
	FIG. 6:

ACKNOWLEDGMENTS

The PhD grant of MS is funded by the OCEVU Labex (ANR-11-LABX-0060), which also provided financial support to this project. We also acknowledge support from the CNRS program Défi InPhyNiTi, and from the European Union's Horizon 2020 research and innovation program under the Marie Sk lodowska-Curie grant agreements No 690575 and No 674896.