Indirect detection of long-lived particles in a rich dark sector with a dark vector portal

Simplified models of light new physics provide a convenient benchmark for experimental searches for new physics signatures, including dark matter (DM). However, additional detection modes can arise in less simplified and more realistic scenarios where new degrees of freedom are invoked. In this study, we introduce a non-minimal model based on a popular dark photon portal to DM where the mediator mass is obtained by interactions with the dark Higgs boson which acts as a long-lived particle. We further add to this scenario a new heavy DM species secluded from the Standard Model. In this model, which involves light and heavy particles in the dark sector, we find some new interesting phenomenological features that lead to complementary probes in intensity frontier searches for light long-lived particles, indirect detection searches for dark matter, and cosmic microwave background surveys. We also find possible non-local effects in the DM indirect detection searches that could significantly affect the usual detection strategies.


I. INTRODUCTION
The nature of the particle content of the dark matter (DM) sector of the Universe remains unknown.The dominant paradigm of the last several decades has been based on the assumption that DM is made up of some weakly interacting massive particles (WIMPs) in the approximate mass range of a few GeV to a few TeV.WIMP DM could naturally be produced in the early Universe via the popular thermal freeze-out mechanism, yielding a correct relic density for WIMP interactions set by a fraction of electroweak interactions of the Standard Model (SM); see, e.g., Ref. [1][2][3] for a review. 1 However, the lack of experimental signal in a wide array of worldwide searches has led one to explore alternative candidates for DM, its production mechanisms, and detection methods; see e.g., Refs [5,6].
For example, the standard freeze-out mechanism can be generalized to broad ranges of DM mass and coupling constants [7]; see also a discussion in Ref. [8].While too large values of DM couplings are constrained by unitarity, lowering them by even several orders of magnitude below the electroweak strength can still render the correct relic density, however, for DM mass much below the GeV scale.Some light mediator particle(s) beyond the Standard Model (BSM) spectrum are present in various alternative scenarios.They specify the interaction between the SM and some DM particle with the mass in the GeV range, or below [9,10].Experimental searches for such light DM (LDM) via direct detection (DD) in-voke new strategies [3,5], while further signatures of associated frameworks may be possible with searches for unstable light mediator species; see Refs [11][12][13][14] for recent reviews.
Such studies are usually performed within a simplified model framework where the DM particle is assumed to be part of some secluded (dark) sector that interacts with the SM sector only via one mediator particle, for instance, a dark photon or a dark Higgs boson; see Refs [15][16][17] for recent reviews.In more UV-complete models, further experimental probes could be possible and provide other complementary ways of their experimental tests, which triggers interest in considering non-minimal scenarios.
In models with a richer dark sector, there are typically more BSM particles with masses spanning a wide range of values, even up to several orders of magnitude.Suppose such a rich dark sector remains secluded from the SM sector, with only a weak portal communicating between them.In that case, the dark couplings of the BSM species are only mildly constrained and can lead to new phenomenological effects.This is especially relevant for cosmological tests and indirect detection (ID) searches of DM that can probe scatterings in the dark sector with the visible signals generated via subsequent decays of BSM species.Notable examples of such effects include, e.g., secluded [10,18] or Sommerfeld-enhanced [19] DM annihilations into the dark sector particles that modify ID bounds and can be constrained by Cosmic Microwave Background (CMB) radiation measurements. 2 In this study, we focus on such an interplay between light and heavy dark species in models that can simul- 2 Collider probes of light new physics can also be modified in the presence of large dark coupling constants via, e.g., possible secondary production of new BSM species in front of or inside the detectors [20].
taneously be probed by experiments targeting signatures related to heavy DM and signals associated with light and long-lived particles (LLPs).As we will show, already in a simple extension of the most minimal scenarios in which we add a heavy DM candidate and much lighter secluded dark species, one can find new complementary detection prospects in future CMB surveys and DM ID searches, in addition to standard LLP searches.
We will focus on a popular dark photon portal, which is often considered to be a mediator between the SM and even very complicated BSM sectors [10,15,18,21].
In particular, such a portal can naturally be related to the hidden valley scenarios predicting the existence of some LLPs [22].Motivated by this, we assume that the relevant new light species has a mass below 1 GeV up to a few tens of GeV and suppressed couplings to the SM sector.The remaining non-minimal BSM content of our model will further contain particles at mass scales up to 10 TeV, or so, including a heavy secluded DM candidate that could be targeted in future ID observations.As we will see, in our two-component DM scenario there are new effects that are not present in the minimal setup and that consequently lead to new experimental signatures.In particular, dark-bremsstrahlung annihilations of the dominant heavy DM particle might produce non-trivial spectra of SM species and be additionally affected by striking non-local DM ID signatures that both require modified search strategies.The presence of very long-lived light mediator species might also leave imprints on CMB observations.In addition, in our model, accelerator-based intensity-frontier experiments could unveil the existence of light dark sector portal species.The combination of the above signatures clearly differs from simplified models with light, GeV-scale DM.
This paper is organized as follows.In Sec.II, we introduce the BSM model of our interest.In Sec.III, we examine the relic abundance of both stable and long-lived dark species present in this scenario.In Sec.IV, we discuss past and future possible bounds on this model.Finally, we present our results in Sec.V and our conclusions in Sec.VI.The technical details of our analysis, including formulae for the relevant cross sections and decay widths, and a discussion of the photon spectrum in the DM ID, among others, have been moved to Appendices A and B. In Appendix C, we highlight possible unique non-local effects in DM ID that could appear in the presence of sufficiently long-lived light mediators.

II. MODEL
The BSM model that we put forward in this article contains a fairly rich dark sector that is coupled to the visible (SM) sector via the kinetic mixing between a dark photon A ′ and the SM hypercharge gauge boson.The dark photon is assumed to be a gauge boson of the secluded U (1) ′ group.Its mass is driven by the vacuum expectation value (vev) of a new scalar field, the dark Higgs boson σ, which is charged under U (1) ′ .Additional light DM particle η is often introduced in the model.It can be coupled to the SM via the dark photon or scalar portal.This is one of the prototype models predicting the existence of LLPs at the MeV-GeV scale, which is discussed in intensity frontier studies [10,18,21]; see, e.g., Refs [23][24][25] for recent analyses.
In our case, to illustrate the interplay between heavy and light new physics, we extend this popular scenario by adding one more heavy DM particle χ.The dark matter sector of the model then comprises of two scalar fields η and χ.The former mixes with the dark Higgs boson.The latter has negligible interaction strength with all the light dark species besides η, to which it couples via an additional heavy scalar field ϕ.This field only plays an auxiliary role in our analysis and could be replaced with a direct renormalizable contact operator, |χ| 2 |η| 2 .The dominant heavy χ DM component is highly secluded from the SM, although it's still thermally produced in the early Universe.The light dark sector species guarantee thermal contact with the SM in a cascade process.A scenario of this type is phenomenologically distinct from traditional WIMPs, as we discuss below.In Fig. 1, we schematically present the connection between the SM and the dark sectors in the model.The fields charged under the dark gauge group U (1) ′ are shown in green.
The Lagrangian of the model can be written as where L SM is the SM Lagrangian, L DS corresponds to the dark sector, and L portal describes interactions between the SM and the dark sector In the following, we assume ϵ ≪ 1 and typically m A ′ < m Z such that the mixing between the dark photon and the SM Z boson can be neglected [26].We also assume, for simplicity, that the mixing between the dark and SM Higgs bosons is negligible; see Ref. [27] for a recent review of the BSM Higgs portal.We parameterize both fields in the unitary gauge after the spontaneous symmetry breaking of both the electroweak and the dark gauge symmetries as follows where the SM Higgs vacuum expectation value is v h = 246 GeV, while the corresponding quantity for the dark Higgs field is denoted by v D .We denote resulting physical eigenvector states by H and h D , respectively.
The following Lagrangian describes the interactions between the dark sector species As we will discuss below, invoking two components of DM scalar fields allows one to relax otherwise stringent bounds on heavy WIMP-like scalars coupled to the SM sector via light mediator species; see Refs [28,29].In the following, we set q ′ h D = 1, while we take q ′ η = q ′ χ = 0. Hence, the η DM component does not couple to the dark photon directly but only indirectly via the dark Higgs boson.The couplings of η to the latter are described by the dark scalar potential terms As a result of the spontaneous breaking of the dark U (1) ′ gauge symmetry, the dark bosons A ′ and h D obtain their masses as The dark sector spans several orders of magnitude in mass, as shown on the left side of Fig. 2.Even for all the dark charges set to unity, the model is still characterized by 9 free parameters.In order to organize our discussion and better highlight interesting phenomenological prospects of this scenario, we assume below the following mass hierarchy in the BSM sector of the model With this assumption, we find below that the dominant DM component is the heavier scalar field χ which decouples earlier from thermal plasma than the lighter species η.In particular, the decoupling of χ proceeds via annihilations through an exchange of the intermediate heavy scalar, χχ → (ϕ * ) → ηη.In addition, 2 → 3 annihilation processes are possible with an outgoing dark scalar emitted from the η leg, χχ → (ϕ * ) → ηηh D .The latter process will play a crucial role in our discussion of ID observables.On the right side of Fig. 2, we show the respective Feynman diagrams for the key processes leading to indirect signals of DM in our model.As can be seen, for m h D ≲ m A ′ ,3 the dark Higgs boson produced in the 2 → 3 process subsequently decays into SM fermions via a triangle loop involving off-shell dark vectors; see Eq. ( 6). 4 As a result, the decay width is naturally suppressed, and h D may have a very large lifetime of astrophysical relevance, The approximate result in Eq. ( 7) is derived from a full expression given in Eq. (A2).As can be seen from Eq. ( 7), the dark scalar h D can travel galactic-scale distances before decaying and producing visible signals, h D → f f .The decay final-states are f f = e + e − , µ + µ − , τ + τ − , or hadrons formed due to hadronization of the light quark pairs q q.Instead, we require A ′ to decay before the Big Bang Nucleosynthesis (BBN) in the early Universe [30], i.e., τ A ′ ≲ 0.1 s.In fact, the lifetime of the dark photon is often within the reach of the intensity frontier searches for light new physics.We provide a complete list of decay width and annihilation cross sections relevant to our discussion in Appendix A.
III. RELIC DENSITY a. Boltzmann equations We begin the discussion of specific features of our model by examining the relic density of both DM components, the dominant one χ and a minor one η -for which we require Ω χ h 2 + Ω η h 2 ≃ 0.12 [31] -as well the abundance of the unstable longlived dark Higgs h D .To this end, we numerically solve a set of Boltzmann equations that extend the familiar assisted freeze-out mechanism discussed for two-component dark sectors [32] to compute the relic densities of three dark species:5 where Y i is the respective yield, i = χ, η, h D , and x = m η /T .We stress again that we assume the lifetime of the dark photon to be small such that A ′ decays before the era of BBN.Therefore, further analysis does not require a precise determination of its abundance.We assume, however, that A ′ remains in thermal contact with the SM radiation around the time of h D freeze-out.In fact, we will focus below on scenarios with a low mass splitting between these two dark species,6 The parameters λ i depend on the annihilation cross sections of the processes that play the dominant role in determining the abundances in the dark sector.Given mass ordering shown in Eq. ( 6), the abundance of the main DM component depends on while λ's for the other leading processes are where m Pl = 2.44×10 18 GeV is the reduced Planck mass, s ≡ s(T ) is the entropy density and H ≡ H(T ) is the Hubble rate.The effective number of degrees of freedom for the entropy and energy densities of the thermalized SM-DS plasma at temperature T is denoted by g * s (T ) and g * (T ), respectively.The equilibrium comoving yield of each particle is defined as Y eq i = n eq i /s, and it explicitly reads: Here, r i = m i /m η and the number of internal degrees of freedom of particle i is denoted by g i .The resulting relic density is obtained from , where s 0 is the present-day entropy density and Y i,0 is the final yield of the dark species i after its freeze-out.In all cases that we consider, χ is the dominant DM component, Ω χ h 2 ≃ 0.12, which undergoes ordinary freeze-out almost independently of the other species because it is much heavier than them; see Eq. (6).The values of the heavy (auxiliary) scalar ϕ mass and the coupling constants µ χ and µ η , see Eq. (A6), are chosen such that the total χχ annihilation cross section is equal to the thermal value ⟨σv⟩ χχ→ηη + ⟨σv⟩ χχ→ηηh D = 2.2 × 10 −26 cm 3 /s.
The abundances of lighter species -η and h D -depend on the couplings λ h D η and g D , respectively, see Eq. (A10) and Eq.(A11).These couplings must be large, λ h D η , g D ≳ 0.1.This is because we require the relic abundance of h D to be suppressed to avoid stringent cosmological constraints; see Sec.IV, and we are interested in regions of the parameter space predicting a large value of the cross section relevant to DM indirect detection, χχ → ηηh D .As a result, we typically find in our analysis Our discussion assumes that kinetic equilibrium is established until chemical decoupling.As discussed in [10], this can be obtained for sufficiently large values of the kinetic mixing parameter This also guarantees that A ′ decays early, before the era of BBN. b.Balanced and standard freeze-out regimes In a large part of the allowed parameter space of the model, the lighter species η and the long-lived h D , undergo the so-called balanced freeze-out [33], which results in phenomenological features different from the ordinary freezeout mechanism.In this regime, the relic density of the sub-dominant DM component η or a long-lived dark scalar h D is characterized by the following scaling Ω ∝ 1/ ⟨σv⟩, instead of the usual scaling Ω ∝ 1/ ⟨σv⟩ characteristic for standard freeze-out.This can be obtained when the annihilation cross sections of the sub-dominant DM components are orders of magnitude greater than the annihilation cross section of the dominant DM component.In the balanced freeze-out case, we predict the following simple scaling between the relic densities of the dominant and sub-dominant components which can be formally obtained by setting dY η /dx = 0 in Eq. ( 8).We note that this relation takes into account the possible production of ηs in annihilations of χ into both two-and three-body final states, χχ → ηη(h D ).A similar relation holds for the relic abundance of the From this expression, it is also clear that the balanced freezeout regime for h D points towards large values of g D ≳ 0.1, as λ h D ∼ g 4 D in Eq. (A11), and small dark Higgs boson mass m h D /m χ ≪ 1.This helps to suppress the relevant relic density.The same holds for η when λ h D η ≳ 0.5.8).We fixed the dark sector masses at the following values: m χ = 1.5 TeV, m η = 250 GeV, m h D = 3 GeV, and the mass splitting ∆ h D A ′ = 0.01; see Eq. ( 12).We also fixed the coupling constants g D = 0.4 and λ h D η = 0.5, while the value of µ χ , µ η , and m ϕ are chosen such that ⟨σv⟩ χχ→ηη is at thermal value and Ω χ h 2 ≃ 0.12, as indicated with the horizontal black dotted line.The brown horizontal line denotes the BBN limit on late-time decaying h D with τ h D ≲ 10 12 s.
By solving the relevant Boltzmann equations, we have numerically verified that the aforementioned simple expressions hold in appropriate regions of the model's parameter space and for benchmark scenarios presented below.We illustrate the evolution of the yield m i Y i for all the three dark species as a function of x in the right panel of Fig. 3.We assume the following mass benchmark: m χ = 1.5 TeV, m η = 250 GeV, m h D = 3 GeV, and a small mass splitting between the dark Higgs boson and dark photon ∆ h D A ′ = 0.01.As discussed in Sec.II, we keep the dark vector mass larger than the scalar mass to avoid direct decays of h D into the A ′ species.At the same time, for small mass splittings, the relic density of h D in the early Universe is still determined by the annihilation process h D h D → A ′ A ′ in the forbidden DM regime [34].
In the figure, we also fix the relevant coupling constants g D = 0.4 and λ h D η = 0.5.As one can see, for sufficiently large g D , the h D abundance may become strongly suppressed with respect to the total DM relic density.In particular, for the assumed values of the model parameters shown in the figure, we obtain Ω h D ≲ 10 −5 Ω χ .With such a small abundance, one can avoid stringent BBN bounds for even very long-lived dark Higgs bosons.We illustrate relevant BBN constraint in the plot assuming τ h D ≲ 10 12 s.

IV. CURRENT AND FUTURE CONSTRAINTS FROM ASTROPHYSICS, COSMOLOGY AND COLLIDERS
Although most of the dark sector species in the model considered here remain secluded from the SM, their indirect couplings via a light dark vector portal A ′ still induce experimental constraints.Below, we briefly summarize the current constraints on our scenario.We also discuss future experimental and observational prospects.

A. Accelerator-based searches
Light secluded dark photons are among the primary targets of intensity frontier searches for sub-GeV new physics as they correspond to one out of only a few available simple renormalizable portals to the dark sector.Based on available data, one can put upper bounds on the kinetic mixing parameter ϵ.Of particular relevance to our study are constraints obtained from CHARM [35,36], E141 [37,38], NuCal [39,40], Orsay [38,41] experiments, and the newest bounds obtained by the FASER [42] and NA62 [43,44] collaborations.
As discussed above, rare meson decays provide stringent bounds on light, sub-GeV dark vector bosons from various past accelerator-based experiments.Similar searches in the future are expected to constrain further the available parameter space of the model.For the kinetic mixing parameter range of our interest, important bounds are expected to come from the ongoing and proposed LLP searches in beam-dump and collider experiments; see Refs [11,12,14] for reviews.In the following, we will show selected such bounds from the DarkQuest [45], FASER2 [46,47], LHCb [48], NA62 [43], and SHiP [49] experiments.
We note that direct searches for heavier dark species at the LHC will be very challenging since we assume that they couple to the SM only via the A ′ portal and the corresponding kinetic mixing parameter is suppressed.Similarly, the dark Higgs boson in our setup remains only very weakly coupled to the SM, and it avoids acceleratorbased constraints.

B. Astrophysical and cosmological bounds
Additional constraints on light dark vectors can be obtained from their possible impact on astrophysical and cosmological observations.The dark photon of our interest, however, avoids bounds due to possible modifications of the supernovae cooling rate and the neutrino emission from SN1987A [50], and BBN constraints on late-time decaying unstable BSM species [30].
Important bounds arise due to the residual relic abundance of potentially very long-lived dark Higgs h D .Following Ref. [51], we implement relevant BBN constraints on late-time electromagnetic energy injection.
We also use the results of Ref. [52] for the CMB bounds based on the combined data from the Planck [31] and COBE/FIRAS [53] satellite missions.In the latter case, we follow Refs [54,55] and modify the CMB constraints from Ref. [52] by taking into account a finite fraction of the electromagnetic energy transferred to the intergalactic medium, f eff < 1, characterizing different SM final states present in the model.The final shape of the cosmological bounds derived this way depends on the interplay between the h D lifetime; see Eq. ( 7), and its relic abundance obtained by solving the Boltzmann equations discussed in Sec.III.
Future surveys are expected to improve the bounds on CMB spectral distortions significantly and to essentially exclude BSM scenarios predicting the mediator lifetime between τ h D ∼ 10 5 s and 10 12 s [52].The relevant bounds will constrain the relic abundance of unstable, long-lived species to be not larger than a tiny fraction of the total DM relic density, Ω h D ≲ 10 −6 Ω DM , but, depending on τ h D , even much more stringent limits can be derived of order 10 −12 Ω DM .Instead, for the CMB anisotropy data, the expected improvement over current bounds is less pronounced but will also result in more stringent constraints by about a factor of a few in Ω h D .This is relevant for the large lifetime regime, τ h D > 10 12 s.
In the following, we will present expected combined bounds from the Planck [31] and the proposed Primordial Inflation Explorer (PIXIE) [56] satellite mission.In addition, we will also employ more stringent constraints expected from the combination of the future data from the Polarized Radiation Imaging and Spectroscopy Mission (PRISM) [57], ground-based CMB Stage-4 (CMB S-4) searches [58] and the space mission LiteBIRD [59].We stress that, while the CMB data remain complementary to DM ID searches for mediator lifetimes τ h D ≳ 10 5 s, they will provide the best way of probing scenarios with extremely long-lived h D s that would predict muchsuppressed ID signal rates.
We note that both the current and future CMB bounds can be additionally affected by the impact of late-time DM annihilations, χχ → ηηh D , producing a subdominant abundance of very long-lived dark Higgs bosons in between the χ DM freeze-out and the CMB epoch.The typical energy of such h D s results in the relevant boost factors of the order γ h D ∼ 10 to 100 for the benchmark scenarios considered below.Hence, such subdominant 2 → 3 DM annihilation processes occurring even around the matter-radiation equality could produce light Higgs bosons that are not completely redshifted before the recombination period.The impact of the CMB bounds and future searches could then be modified in the regime of the large h D lifetime close to the current bounds, e.g., τ h D ≳ 10 11 s in the plots shown in Sec.V.While this effect is not expected to affect our discussion below significantly, we leave a relevant detailed analysis for future studies.

C. Dark matter indirect detection
The model of our interest can also be probed in DM ID experiments.The corresponding discovery prospects rely on annihilation rates that depend only on unsuppressed couplings in the secluded dark sector.On the other hand, future DD searches remain much less promising.The lighter DM species η couple to the SM only indirectly.Their suppressed couplings to light quarks and negligible abundance typically yield very low expected signal rates in DD experiments.This is also true for the heavier DM species χ for which the scattering rates off nuclei appear only via intermediate η and ϕ fields.
ID signatures in our model can arise from annihilations of both DM species, the heavier dominant scalar χ and the lighter one η.The latter process, however, is highly suppressed by the tiny abundance of η.This is not the case for annihilations of the dominant DM component.At the leading order, though, the main annihilation channel is purely into invisible final states χχ → ηη.Instead, the dominant contribution to ID signal rates appears at the next-to-leading order due to the 2 → 3 process χχ → ηηh D shown in Fig. 2.This is especially important for a growing value of the dark coupling constant λ h D η which increases the probability of emitting the final-state h D off the η leg; see Eq. (A12). 7Notably, as we have argued above, larger values of λ h D η are also preferred by cosmological and collider constraints.
We note that DM ID prospects in our model rely on the 2 → 3 annihilation process producing a continuous (not monochromatic) spectrum of the meta-stable h D mediator energies.This results in smearing the spectrum of final-state SM particles; see Appendix B for discussion about the photon spectrum.Hence, one effectively avoids bounds from the searches for the peaked spectral features in the positron data; see Refs [60,61].On top of this, interesting non-local effects can further affect DM ID event rates in this scenario [62][63][64][65][66].We discuss them in more detail in Appendix C.
A particularly promising way of probing these scenarios is via searches for a diffuse DM-induced γ-ray flux.The corresponding differential flux of photons coming to a detector from the angular region in the sky ∆Ω is given by where on the left-hand side, we have denoted that this result corresponds to the "standard" regime in which the long-lived mediator h D decays after traveling distances much shorter than a kpc.In our analysis below, we employ the Einasto DM profile, where ρ s = 0.079 GeV/cm 3 , r s = 20 kpc and α = 0.17 [67].In Eq. ( 18) the photon spectrum from 2 → 3 annihilations of χ and subsequent decays of h D is denoted by (dN γ /dE γ ) χ ; see Appendix B for further discussion.Importantly, the above chain of processes results in a much softer photon spectrum than one could expect based on the DM mass.In particular, a typical h D energy after the initial χ annihilation is of the order The peak energy of photons produced as a result of h D decays is then further shifted towards smaller energies, E γ ≲ few tens of GeV.However, it contains a high-energy tail extending towards larger E γ .In Sec.V, we will present the expected sensitivity reach of the forthcoming DM ID experiment targeting heavy DM species, namely the Cherenkov Telescope Array (CTA) [68]; see also Refs [69,70].Instead of performing full detector simulations for the scenario of our interest, which is beyond the scope of our study, we illustrate the impact of the CTA on the parameter space of the model with some approximate expected bounds.We obtain them by employing the CTA sensitivity plots for the secluded DM regime following Ref.[71], which we, however, modify by taking into account characteristic h D energies obtained after the 2 → 3 annihilation process.We expect this simplified procedure to encompass the most essential effects of our study correctly.Notably, the applicable CTA bounds for secluded DM only mildly depend on the DM mass in the regime between 100 GeV ≲ m DM ≲ few TeV.

V. RESULTS
The non-minimal content of our model results in 9 free parameters that contain the masses of the dark species and various dark sector couplings; see Sec.II.Our goal in this paper has been to point out some new signatures that arise in the model, rather than to perform a thorough investigation of its rich phenomenology.We, therefore, limit our discussion to only slices of this multidimensional parameter space.Below, we first justify our choice for the fixed parameters, and later we present the results of our analysis in a convenient two-dimensional (m A ′ , g D ) plane.
As discussed above, we work in the regime of a small mass splitting between the dark Higgs boson and the dark photon, ∆ h D A ′ = 0.01 in Eq. (12).We employ three different values of the kinetic mixing parameter ϵ = 10 −4 , 10 −5 , and 10 −6 .We set m χ = 1.5 TeV for the heavy DM particle and m η = 250 GeV for the lighter one.In our plots, the η DM component is kept heavier than the dark Higgs boson, m h D < m η .This allows for efficient annihilations ηη → h D h D and suppresses the η abundance to negligible levels.Otherwise, it will significantly contribute to the total DM relic density and could be overproduced.We fit the auxiliary parameters in the dark sector, m ϕ , µ χ , and µ η , such that at each point in the parameter space shown in the figures below, the heavy scalar DM obtains the correct value of the thermal relic density, Ω χ h 2 ≃ 0.12.In contrast, Ω η h 2 is negligible.
The values of the dark coupling constant g D are only limited by the perturbativity bound.We take it to be α D /4π < 1, which originates from the fact that such combination enters into 1-loop corrections.This condition is equivalent to g D < 4π.The astrophysical and cosmological bounds discussed in Sec.IV constrain this coupling constant from below.In particular, lower values of g D would lead to a too large abundance of h D , and its late-time decays would violate BBN and CMB limits.Instead, collider bounds affect the region of the parameter space with too light dark photons below a few hundred MeV.All the bounds are shown as gray-shaded regions in the figures.
We present expected sensitivity reach lines of future searches in Fig. 4. The top (central, bottom) panels correspond to ϵ = 10 −4 (10 −5 , 10 −6 ).On the left, we fix the coupling strength between the dark Higgs boson and η DM component to λ h D η = 0.5.In the right panels, we assume it saturates the perturbativity bound λ h D η = 4π.We note that in our analysis, we assumed for simplicity that the bare value of µ h D η vanishes, and instead µ h D η = λ h D η v D is obtained only after the U (1) dark spontaneous symmetry breaking.A non-zero bare value of µ h D η could lower the value of λ h D η corresponding to the regime of long-lived h D significantly below ∼ O(1).As can be seen, in each case, the allowed region of the parameter space corresponds to the A ′ mass in between tens or hundreds MeV and m A ′ ≃ m η = 250 GeV.The dark coupling constant varies in this region between g D ∼ (0.01 − 0.1) and the perturbativity bound.
Scenarios with the sub-GeV dark photon can be constrained in the future in collider-based searches, as we show with vertical, colorful lines in the plots.These projected constraints and past bounds do not depend on the coupling strengths g D and λ h D η , while they vary for different values of the kinetic mixing parameter ϵ.Future LHCb searches are expected to provide the most severe such bounds for ϵ = 10 −4 and 10 −5 .For ϵ = 10 −6 , the dominant bounds will be found by DarkQuest, FASER2, and SHiP experiments.
Complementary bounds from future CMB surveys are indicated in the plots with light blue and red dashed lines.They will exclude too low values of g D , for which the h D lifetime grows; see Eq. (7).Instead, these means do not constrain the region in the model's parameter space with large values of g D , close to the perturbativity limit.The relic abundance of h D remains, in this case, highly suppressed beyond the expected sensitivity of the future CMB and BBN analyses.The projected constraints and past bounds become more severe for decreasing ϵ, which is due to its impact on the h D lifetime in a loop-induced decay process; see Eq. (7).In particular, for ϵ = 10 −6 , future CMB-S4 measurements will be able to probe almost the entire remaining available parameter space of this model.
We also show in the plots regions of the parameter space that the future DM ID searches in CTA can probe.These also constrain the available parameter space from below, i.e., they will exclude too low values of g D .This is because the 2 → 3 annihilation cross section determining the DM ID signal rates decreases for increasing value of this dark coupling constant; see Eq. (A12) in which σ χχ→ηηh D ∝ λ h D η /g D and we fix λ h D η in the plots.We note, however, that for natural values of the dark coupling constant g D ≲ 1, some regions of the available parameter space of the model might be probed in future DM ID searches.This especially concerns scenarios with increasing m A ′ ≃ m h D .Instead, the 2 → 3 annihilation cross section is suppressed for decreasing m A ′ .The projected CTA bounds will only mildly constrain the available parameter space of the model for λ h D η = 0.5 or ϵ ≲ 10 −5 .Instead, they become more severe with a growing value of these coupling strengths and become the dominant indirect bounds for ϵ = 10 −4 and λ h D η = 4π.The dependence on λ h D η is driven by the aforementioned 2 → 3 annihilation cross section.On the other hand, the DM ID signal rates do not depend on ϵ.However, with decreasing value of this parameter, these bounds become less competitive with CMB and BBN constraints.
For sufficiently low values of g D , the lifetime of the dark Higgs boson is very large, and non-local effects in DM ID become important.These can suppress DM ID signal rates, especially for small regions of interest around the Galactic Center.We show it in the plots with dashed, dotted, and dash-dotted purple lines; see Appendix C for detailed discussion.

VI. CONCLUSIONS
Light sub-GeV portal models of dark matter have gained much interest in recent years and have sparked both experimental and theoretical activity in the field.Most studies so far have focused on simplified frameworks with only a limited number of new species added to extend the SM.These typically correspond to the most popular interaction operators and are supposed to encompass essential phenomenological aspects of such simplified scenarios.While this approach allows one to do a manageable comparison among many experimental proposals, it is also understood that new effects might be observed in more elaborate, and more realistic, models incorporating a larger number of BSM species with a mass hierarchy that can be quite complex and in particular span several orders of magnitude.In this paper we have presented an example of such a scenario.
The presence of the light portal particle that thermally   connects the SM sector and a heavy DM in the early Universe leaves important observational imprints that can be distinct from those in both more popular heavy WIMP models and from scenarios predicting the existence of only light DM.In particular, new detection prospects can appear, in this case, due to the interplay between intensity frontier searches, cosmological observations, and DM indirect detection.
In this study, we have found such effects and exposed them in a model which extends a popular scenario with a dark photon portal to DM and employs the dark Higgs boson vev to drive the light dark sector masses.We have discussed the consequences of adding a new heavy DM particle to this setup coupled predominantly to the lighter DM particle, ensuring its secluded nature differs from traditional WIMPs.This scenario remains beyond the reach of current and near-future direct detection experiments and collider searches for heavy DM.However, we have shown here that the best way of probing this model type, besides intensity frontier searches for LLPs, is to employ DM-induced signatures in future ID and CMB experiments.
The presence of unstable subdominant dark species produced in annihilations of much heavier dominant DM particles can lead to further striking signatures.We have illustrated this in the case of the dark Higgs boson h D that can feature a very large loop-suppressed lifetime and an astrophysically interesting decay length γβcτ h D ∼ 1 kpc.This can lead to interesting non-local effects in DM ID searches, for instance, enhanced DM signal rates from the GC and simultaneously suppressed corresponding rates expected from the dwarf galaxies.A thorough experimental testing of such scenarios will require going beyond the traditional approach to DM ID.
In searching for new physics, it remains essential to encompass a broad range of possible BSM scenarios and to study all potential complementarities among different types of experiments.An attractive framework based on the popular thermal DM paradigm has led to numerous studies in the past years.Models of this class that predict the existence of a light sub-GeV portal to heavy secluded DM deserve special attention in such efforts.As shown here, their characteristic signatures are absent in more popular simplified scenarios but may be successfully probed in extensive experimental programs in the coming years.
a. Acknowledgements We would like to thank Brian Batell for useful remarks and comments on the manuscript.We would like to thank Luc Darmé, Iftah Galon, Andrzej Hryczuk, Arvind Rajaraman for useful discussions.In our analysis, we employed python module vegas [72], which uses an algorithm developed in Refs [73,74]  Below, we provide expressions for all the decay widths and scattering cross sections relevant for our analysis.
1. Decay widths a. Decay width of A ′ It is given by where where m e − is the electron mass and B e = B(A ′ → e + e − ) is the branching fraction of a decay into an electronpositron pair, which can be found, e.g., in Ref. [75].b.Decay width of h D The dark Higgs plays the role of the LLP in our setup, which can be naturally achieved due to loop suppression of the decay width into SM fermions [18]: where the I loop function is given by ) and an analogous expression holds for h D decays into pairs of other SM charged states.
c. Decay width of ϕ While it is not essential for our study, we also provide for completeness the decay width of the heavy scalar ϕ, which can decay into both the χ χ and η η pairs: (A4) In the assumed mass scheme, Eq. ( 6), and for µ χ = µ η , µ χ /m ϕ ∼ 0.1, we obtain a very small lifetime of ϕ, which then does not have any implications on our discussion of the cosmological bounds

Annihilation cross sections
a. χχ → ηη The relic density of χ is driven by annihilations into two lighter species η, as shown on the left in Fig. 5.The relevant cross section is given by The relic density of η is obtained via its annihilations into the light mediators, A ′ and h D , as shown in Figs. 6 and 7, and are described by the following formula: Figure 7: Feynman diagrams for the η annihilation processes into a pair of dark Higgs boson.which, in the limit of m η ≫ m A ′ , m h D , simplifies to c. ηη → h D h D The cross section for η annihilations into light dark Higgs bosons reads which, in the limit of m η ≫ m h D , simplifies to The metastable relic density of the dark Higgs is obtained thanks to its annihilations into A ′ A ′ pairs, as shown in Fig. 8.As we are working in the forbidden DM regime, m h D < m A ′ , we give the formula for the cross section for the opposite reaction, We connect the two processes by using the detailed balance condition, e. χχ → ηηh D For small values of the dark gauge coupling g D , the dark bremsstrahlung process becomes important, because the coupling between a pair of ηs and a h D is proportional to the dark Higgs vev, v D = m A ′ /g D .The amplitude of the process depends on, i.a., the momenta of all the particles and the total energy in the CM frame s = (p 1 + p 2 ) 2 : where the incoming χ particles have momenta p 1 , p 2 , the outgoing η particles have momenta p 4 , p 5 and the outgoing h D has momentum p 3 .The total cross section is given by8 (A13) To obtain this expression, we have already performed a trivial integral over the azimuth angle ϕ, which gives a factor of 2π.The angles written without the asterisk are evaluated in the CM frame of a pair of particles with momenta p 1 and p 2 .In turn, quantities with the asterisk ( * ) are evaluated in the CM frame of pair of particles with momenta p 4 and p 5 .For example, θ is the scattering angle between the particle with momentum p 3 and the collision axis of the particles with momenta p 1 and p 2 in their CM frame, while θ * and ϕ * denote the angles of the particle with momentum p 3 in the frame where ⃗ p 4 + ⃗ p 5 = 0.In our case, we set: In order to obtain the quantities with the asterisk, one needs to use the following boost factor The invariant mass m 45 and the energy of the pair E 45 are given by We also write down the expressions for all components of p 4 and p 5 : where Finally, we can express p 4 and p 5 in the CM frame as a function of (E 3 , θ, θ * , ϕ * ) by applying the following boost and rotation transformations where the usual boost and rotation matrices are given by In this appendix, we provide further details about our modeling of the DM-induced photon spectrum in the model under study.For the primary spectra of γ rays, we rely on the PPPC data [77]; see also Ref. [78] for useful discussion.We denote the relevant spectrum obtained in the rest frame of the dark Higgs boson by (dN γ /dE ′ γ ) h D .We then boost the photon spectrum into the χ rest frame.For each fixed value of the h D energy, E h D , we obtain the corresponding photon spectrum from the χχ → ηηh D process as where E γ denotes the photon energy in the χ rest frame. 9n order to obtain the final γ-ray spectrum (dN γ /dE γ ) χ , we then convolute this with the actual continuous h D spectrum from 2 → 3 annihilations, χχ → ηηh D ; see Eq. (A13) for the relevant differential cross section.In order to simplify our analysis, we evaluate this as a discretized weighted average where ⟨σv⟩ E h D corresponds to the 2 → 3 annihilation cross section integrated over a limited range of the outgoing h D energies (within the energy bin centered around In this case, the dark Higgs boson decays promptly afbeing produced in the 2 → 3 annihilation process of χ DM.Instead, the black dash-dotted, red dotted, and blue dashed lines correspond to the long-lived mediator case with the relevant typical decay length of boosted h D equal to 0.1, 3, and 30 kpc, respectively.E h D ), while ⟨σv⟩ is the total such cross section.We used 20 equally spaced bins in x h D = E h D /E χ in the logarithmic scale, which has been numerically verified to be sufficient for the results presented in Sec.V.
An example of such a photon spectrum is shown with the black solid line in Fig. 9.We have obtained it for the fixed masses of the dark sector species m χ = 1.5 TeV, m η = 250 GeV, m h D = 5 GeV, and ∆ h D A ′ = 0.02.The coupling constants were fixed as follows: ϵ = 10 −6 , α D = 0.1, and λ h D η = 0.5.As expected, even though the assumed χ mass is above TeV, the resulting photon spectrum is much softer and peaked at E γ ∼ 100 GeV.
In Fig. 9, we further present several such photon spectra corresponding to very long-lived mediators with the decay length of order dh D ∼ 0.1, 3 and 30 kpc.We present the photon fluxes for the |b|, |l| < 12 • region around the GC.As discussed in Appendix C, for dh D ≪ d RoI ≃ 2.3 kpc the impact of non-local effects on the observed spectrum is very small and the photon spectrum resembles the one obtained in the short lifetime regime.This is indicated by the black dash-dotted line in the figure.Instead, for larger values of the decay length, dh D ∼ d RoI , we see a relative increase in the photon flux observed from this RoI due to anisotropy effects, as shown with the red dotted line.Finally, for very large decay lengths, dh D ≫ d RoI , the flux becomes suppressed, which is illustrated with a blue dash-dotted line.The suppression is more pronounced for energetic photons as they originate from more boosted dark vectors.These can more easily escape the RoI before decaying.As a result, the observed photon spectrum is even more shifted towards lower energies.As discussed in the main text, the phenomenology of the type of model of our interest differs from both traditional WIMPs and light DM species.Especially interesting such effects can be expected in indirect DM searches where secluded heavy DM can be detected thanks to non-negligible annihilation rates into lighter dark sector species.Importantly, the light BSM species can even be very long-lived and constrained by cosmology.This opens up new detection prospects for this kind of LLPs in the lifetime regime that remain highly complementary to intensity frontier searches.
In DM ID searches, this can lead to striking non-local effects [62][63][64][65][66].In this case, the interplay between different methods of performing DM ID observations could be challenging to interpret in standard WIMP scenarios, while it could indicate the existence of a dark sector comprising both heavy WIMP-like DM particles and much lighter and long-lived species.Additional non-local effects that can modify signals include enhanced DMinduced signal rates from extensive regions outside the Galactic Center (GC) and decreased signals from individual small regions, e.g., around dwarf spheroidal galaxies (dSphs).Below, we first discuss these effects for general WIMP DM and then present the possible impact on the parameter space of the model under study.
In the non-local case, the DM-induced γ-ray flux is partially driven outside the dense region around the GC by a very long-lived mediator h D .As a result, both the energy spectrum and morphology of the signal might be affected.The respective photon flux reads where dh D = cτ h D γ h D β h D is the decay length of h D in the Galactic frame and vectors ⃗ r χ , ⃗ r h D , and ⃗ r GC correspond, respectively, to the position of the χ annihilation, the h D decay, and the GC with respect to the detector on Earth.As can be seen, compared with the standard case, Eq. ( 18), in the non-local DM ID an additional integration appears over the position of the initial χ annihilation.This takes into account the fact that the longlived mediator h D produced in 2 → 3 processes at ⃗ r χ can travel long-distances before decaying at position ⃗ r h D .In particular, the initial position ⃗ r χ can lie outside the region of interest (RoI) in a given DM ID analysis.The mediator decay probability decreases exponentially with the growing distance |⃗ r h D − ⃗ r χ |.Hence, typically only a limited region in the Galaxy around the RoI contributes to the observed DM-induced photon flux, although this depends on the value of dh D .
In Eq. (C1), we also employ anisotropy factors that depend on the angle θ defined as the angle between the h D boost direction and the detector, In the case of a two-body decay of h D → XX, the function f (θ) reads in which θ is the relevant angle in the h D rest frame and we obtain cos θ = cos θ+ for θ ≤ π/2 and cos θ = cos θ− otherwise, where In the simplest case (not relevant to our model), in which the traveling mediator directly decays into a pair of photons, X = γ, one would reproduce a known expression for radiative beaming, f (θ) = γ 2 h D (β h D cos θ + 1) 2 .The anisotropy factors appear due to the boost of the decaying h D in the Galactic frame.They affect the final observed photon flux at Earth since decaying mediators preferentially come from the direction of the GC in each given region in the Galaxy.However, in the non-relativistic and local limit, we obtain f (θ) → 1.
We illustrate the impact of non-local effects on the total integrated flux of γ-rays in a toy DM model with m DM = 100 GeV and a long-lived mediator mass m med = 10 GeV in the left panel of Fig. 10.Here, we assume for simplicity that the mediator decays directly into a pair of photons.We show with the red lines the ratio between the flux obtained in the long and short mediator lifetime regimes for several different RoIs.In particular, the solid red line corresponds to a large region around the GC characterized by |b|, |l| < 12 • .10While this is a larger region than for typical CTA analyses, we employ it to highlight the difference between the searches better focused on the closed vicinity of the GC and the possible extended Galactic center survey present in the non-local DM ID regime; see, e.g., Refs [68,70] for further discussion about CTA analyses.The large RoI employed in our study extends to roughly d RoI ∼ R 0 sin b ≃ 2.3 kpc distance away from the GC, where R 0 = 9 kpc is the distance between the Earth and the GC.
As can be seen in the figure, for dmed ≲ d RoI the impact of non-local effects on the observed spectrum is very small.The photon spectrum resembles the one obtained in the short lifetime regime denoted by the horizontal black solid line, Φ nonlocal /Φ stand.≃ 1.In contrast, for very large decay lengths of A ′ the expected DM-induced photon flux coming from the RoI drops down much below the standard expectations.This is due to the efficient escape of mediators away from the GC before they decay.The decrease of the flux is roughly linear in growing dmed , as can be seen from Eq. (C1) in the limit of dmed ≫ |⃗ r med − ⃗ r χ | ∼ d RoI .
The relative increase of the DM-induced photon flux for intermediate values of dmed ∼ d RoI can be understood as follows [64].In this case, most of the mediators produced close to the GC decay within the RoI.Only a small fraction of dark photons produced at the GC would generate photons traveling toward the Earth.In the standard case, the signal from other mediators will be lost.However, dark mediators can partially overcome this in the non-local regime by traveling away from the GC before decaying.At these remote positions, they can produce photons in the direction of the Earth, which would not be seen had they been produced close to the GC.As a result, the DM ID signal rates from distant positions within the RoI receive contributions from both local DM annihilations and annihilations occurring close to the GC.This increase is even more pronounced for a modified RoI around the GC in which we exclude the innermost region with the 2 • size aperture.We present this with the dashed red line in the figure.In this case, the photon flux in the non-local regime gains even more from the dark mediators traveling inside the RoI from the innermost region close to the GC.
In addition, we also show in the figure the expected photon fluxes for much smaller RoIs.Here, with a red dotted line, we present the results for a small RoI characteristic for the CTA Galactic center survey, in which we have additionally excluded part of the region very close to the GC, i.e., we assume 0.3 • < |b| < 1 • and |l| < 1 • .Instead, with a red dash-dotted line, we present the flux for an even smaller region around the GC with |b|, |l| < 0.5 • .The size of this region encompasses a typical DM halo size for dwarf galaxies in the Fermi-LAT analyses [79].
As can be seen, for both small RoIs, the relative growth of the flux for smaller dmed is hard to reconstruct, while for the decay length of order several kpcs, the flux is already suppressed.We note that if dSphs are modeled as point-like sources in the analysis with the 0.1 • × 0.1 • bin size [80], the impact of non-local effects is even stronger.Last but not least, we stress that in the non-local regime, the DM-induced photons escaping from small RoIs could also affect the analysis based on local background expectations around each of the dSphs.This could further ameliorate the relevant bounds.
As can be seen, for dmed ∼ a few kpc, the difference between the impact of non-local effects for DM ID focusing on extensive regions around the GC and for searches targeting small RoIs can reach up to a factor of a few in the predicted photon flux.One can then expect a weakening of the relevant DM constraints derived based on stacked dwarf analyses, while DM-induced signals could be stronger for searches using larger RoIs around the GC.This might open new possibilities in explaining persisting anomalies in DM searches; see, e.g., Ref. [66] for the relevant discussion about the Galactic Center Excess (GCE) and non-local DM ID effects.However, as mentioned above, such effects would also modify the morphology of the DM-induced signals, which could then no longer follow the original DM profile but appear less cuspy.In particular, as shown in Ref. [66], when compared to the standard short-lived regime, DM solutions to the GCE employing non-local effects struggle to improve the global χ 2 fit for this anomaly.
In the right panel of Fig. 10, we illustrate the impact of non-local effects on DM ID searches for the aforementioned toy model.To this end, we compare the expected sensitivity reach of the CTA in a secluded WIMP DM scenario presented in Ref. [71] with the relevant reach obtained for a very long-lived mediator with τ med = 10 9 s and fixed m med = 10 GeV.Here, we assume that the mediator decays into light quarks.In the figure, larger values of the mass of annihilating DM m DM imply larger values of the boost factor of the mediator and the corresponding decay length dmed ≃ (m DM / 1 TeV) (10 GeV/m med ) × 1 kpc.This results in an effective suppression of the DM-induced signal from a small RoI around the GC for m DM ≳ a few hundred GeV, as indicated with the dotted red line in the figure.We also show there with the red dashed line the expected sensitivity for a larger RoI |b|, |l| < 12 • , which, for illustrative purposes, has been re-scaled to match the sensitivity of the small RoI in the limit of low m DM .As can be seen, for the extended RoI, the weakening of the future bounds corresponds to larger DM masses than for small RoI.In this case, we also observe a relative improvement of the bound in the intermediate region of m DM ∼ a few TeV.This is due to the excess DM-induced photon flux for dmed ∼ d RoI , as discussed above.
In Fig. 11, we show the results of a similar analysis of the BSM model described in the main text.In the left panel, we fix the heavy DM mass to m χ = 1.5 TeV, while it is increased to 15 TeV in the right panel.The solid purple lines in the plots correspond to projected CTA bounds obtained when neglecting non-local effects.We also show in the plots dashed and dotted lines obtained for two different RoI, as discussed above, for which the impact of increasing h D decay length is considered.As can be seen, for small values of m A ′ ≃ m h D in the plot, which assumes ∆ h D A ′ = 0.01, the projected CTA bounds become much weaker.This is expected based on the growing escape probability of h D s outside the RoI.The same effect is also shown in Fig. 4 in the main text.One also finds the aforementioned enhanced signal rate effect in the plots for selected values of the dark sector masses.
Comparing left and right panels, we see that the nonlocal effects become important at even larger m A ′ ≃ m h D for increasing m χ .This is because more boosted light dark Higgs bosons are then produced in the 2 → 3 annihilation process, χχ → ηηh D .In addition, all the DM ID projected bounds are shifted to higher values of g D in this case.This is primarily driven by an increase in the typical energy carried away by h D after χ annihilation, which corresponds to the region with an improved CTA sensitivity.The negative impact of growing g D on the dark Higgs boson lifetime in Eq. ( 7) partially compensates the effect of increasing the boost factor of h D in estimating the non-local effects.As a result, the shift in m A ′ ≃ m h D is lower than a naively expected order of magnitude difference driven by the change in m χ .

Figure 1 :
Figure 1: Schematic illustration of the model.The dark sector (DS) is connected to the Standard Model (SM) through a light dark photon portal (A ′ ), which kinetically mixes with the SM hypercharge gauge boson.We denote the mixing parameter by ϵ.The mass of A ′ is determined by the vev of the additional dark Higgs boson h D .The vector field and the dark Higgs field charged under the secluded U (1) ′ group are shown in green.The dark matter sector indicated with the black color consists of the additional scalar fields η and χ.The former mixes with h D , while the latter is heavy and is the dominant DM component.It couples to η via the auxiliary scalar field (ϕ).

Figure 2 :
Figure 2: On the left side, a schematic illustration of the mass hierarchy across the energy scales for dark sector particles in the model is shown.The unstable mediators are denoted in dark-red, while the two stable DM species are in black.On the right side, Feynman diagrams for processes leading to DM indirect detection signatures are shown.From top to bottom, these include the 2 → 3 annihilation process of χ DM producing the dark Higgs boson h D and the loop-induced h D decay into the SM species.

Figure 3 :
Figure3: Comoving energy densities, (mY ) i , for twocomponent DM and the unstable dark Higgs boson undergoing thermal freeze-out obtained by numerically solving the set of Boltzmann equations given by Eq. (8).We fixed the dark sector masses at the following values: m χ = 1.5 TeV, m η = 250 GeV, m h D = 3 GeV, and the mass splitting ∆ h D A ′ = 0.01; see Eq.(12).We also fixed the coupling constants g D = 0.4 and λ h D η = 0.5, while the value of µ χ , µ η , and m ϕ are chosen such that ⟨σv⟩ χχ→ηη is at thermal value and Ω χ h 2 ≃ 0.12, as indicated with the horizontal black dotted line.The brown horizontal line denotes the BBN limit on late-time decaying h D with τ h D ≲ 10 12 s.

∆Figure 4 :
Figure4: The impact on the parameter space of the considered model of various past bounds (gray-shaded regions) and future searches is shown in the (m A ′ , g D ) plane; see the text for details.We fix other parameters of the model to m χ = 1.5 TeV, m η = 250 GeV, and m h D = 0.99 × m A ′ .We also fix the kinetic mixing parameter to ϵ = 10 −4 (10 −5 , 10 −6 ) in the top (central, bottom) panels, and we assume λ h D η = 0.5 (4π) in the left (right) panels.In each point in the figure, the parameters m ϕ , µ χ , and µ η are chosen such that the correct value of the heavy DM relic density is obtained, Ω χ h 2 ≃ 0.12.

Figure 8 :
Figure 8: Feynman diagrams for the dark vector annihilations into light dark Higgs bosons.

Figure 9 :
Figure9: The differential flux of photons produced in the DM-induced cascade annihilation process described in the text for the model under study is shown as a function of the photon energy.In the figure, we show the result for the standard scenario with the black solid line.In this case, the dark Higgs boson decays promptly afbeing produced in the 2 → 3 annihilation process of χ DM.Instead, the black dash-dotted, red dotted, and blue dashed lines correspond to the long-lived mediator case with the relevant typical decay length of boosted h D equal to 0.1, 3, and 30 kpc, respectively.
Appendix C: Non-local effects in γ-ray DM ID for χ

Figure 10 :
Figure 10: Left: The ratio of the integrated photon fluxes obtained for increasing decay length d of the mediator and in the standard regime of prompt decays shown with the horizontal black solid line.The figure has been preparedfor the toy model with the DM and mediator masses equal to m DM = 100 GeV and m med = 10 GeV, respectively, and assuming that the mediator decays directly into a pair of photons.We integrate the fluxes over the energy range between 0.1 GeV and 100 GeV.The solid (dashed, dotted, dash-dotted) red line corresponds to the DM ID observation region around the GC defined by different longitude and latitude limits, as indicated in the figure.Right: The CTA sensitivity for the secluded WIMP DM scenario in the (m DM , ⟨σv⟩) plane is shown with the black solid line following Ref.[71].For comparison, we also present with the red dotted line the expected such reach for the toy model with the long-lived mediator with τ med = 10 9 s and m med = 10 GeV which decays into light quarks.The red dashed line corresponds to a re-scaled sensitivity for a larger RoI, as indicated in the figure (see text for details).

Figure 11 :
Figure 11: Similar to Fig. 4 but here, zoomed-in plots are shown focusing on DM indirect detection bounds obtained based on γ-ray search toward the Galactic Center.In the plot, we fix λ h D η = 4π and assume m χ = 1.5 TeV (15 TeV) in the left (right) panel.The solid line corresponds to the projected CTA bound obtained without considering nonlocal effects.Dashed and dotted lines illustrate the impact of non-local effects on analyses for two different regions of interest around the GC, as indicated in the plots.
. KJ and LR are supported by the National Science Centre, Poland, research grant No. 2015/18/A/ST2/00748.LR and ST are supported by the grant "AstroCeNT: Particle Astrophysics Science and Technology Centre" carried out within the International Research Agendas programme of the Foundation for Polish Science financed by the European Union under the European Regional Development Fund. ST is supported by the National Science Centre, Poland, research grant No. 2021/42/E/ST2/00031.KJ is supported by the IBS under the project code, IBS-R018-D1.