Chapter 5 Dark Matter and New Physics Beyond the Standard Model with LHAASO

Over the past century, a growing number of astrophysical and cosmological observations has led to overwhelming evidence for the so-called "dark matter" (DM), which is non-luminous and non-relativistic (hence dubbed cold). While it is not unusual in astrophysics to detect unseen objects gravitationally first or exclusively, the existence of DM cannot be explained in the current standard model (SM) of particle physics plus general relativity.

Given the universality of the gravitational force, in order to identify the nature of DM, it is required to detect other interactions of DM. Candidates for cold DM have been proposed in an extremely wide range of masses and interaction strengths. These parameters notably control its main observable: the DM average cosmological abundance, now measured to percent precision thanks to Cosmic Microwave Background (CMB) data, which must be correctly reproduced within a specific production mechanism.

Despite the huge parameter space, over the past few decades the dominant DM searches have overwhelmingly focused on the so-called Weakly interacting massive particles (WIMPs) paradigm, which features stable particles with masses in the electroweak scale region (give or take one order of magnitude) and couplings to the SM of electroweak strength, produced as thermal relics in the early universe. This has been promoted especially by

(1). The rather predictive thermal freeze-out production mechanism.

(2). A widespread theoretical expectation of new particle physics (such as Supersymmetry, Extra Dimensions etc.) at or slightly above the weak scale in order to fix the hierarchy problem of the Higgs mass. Such new physics produces (or can easily accommodate) also a WIMP DM particle as a byproduct.

(3). The numerous approaches that such a possibility would open for experimental searches, with technologically accessible and mature technologies: via nucleus recoils induced by DM scattering in underground experiments (direct detection), via direct searches at colliders, or indirectly detecting the SM byproducts of DM annihilation (the same process setting their relic abundance in the early universe) in space or ground-based observations of astrophysical fluxes (indirect detection).

Null results of all the searches at point (3). have often pushed the surviving WIMP candidates to heavier masses, notably above the kinematical production threshold of collider searches. There is thus a renewed interest of traditional DM searches above TeV scale, where both forthcoming Imaging Cherenkov Telescopes (notably CTA) and the Large High Altitude Air Shower Observatory (LHAASO) may contribute.

On the other hand, in recent years the trend has been to tackle the DM problem on its own, exploring alternatives to the baseline production mechanism (point (1).), to break more and more free from the "theoretical injunction" (point (2).) in building DM models. More and more studies explore alternative detection strategies, channels, and energy windows, compared to the limited framework of point (3). above. In that sense, the window roughly spanning the energies 0.1-10 PeV has become an interesting new target for DM modeling, for a variety of reasons. We can loosely classify them as:

1. Models explicitly built to account for a DM candidate, with mass and interaction scales emerging rather from the broader possibilities allowed for a viable production than from theoretical bias. In the standard freeze-out picture, an upper limit on the DM mass, of order 100 TeV, follows from the unitarity (i.e. the sum of probabilities equals one) upper limit on the DM annihilation cross section [1], see [2] for a recent appraisal. For the sake of logical organisation, let us classify the models of heavy DM in terms of the assumptions they drop, in the derivation of that result.

1.1 A straightforward option is that DM was never in chemical equilibrium with any bath. In this case DM can be produced in the right abundance by a multitude of mechanisms [3], just to name a few via freeze-in [4], via (p)reheating and/or inflaton decays (see the recent [5] and references therein), at phase transitions [6, 7], at inflation production and decay of heavy Gravitinos or Polonyi fields into the lightest supersymmetric particles (LSP) [8–10] etc. Such DM candidates can easily have the masses needed to give signals at LHAASO, e.g. via their decays.

1.2 An option that allows to keep DM at equilibrium with some (SM or 'dark') bath is to modify the standard radiation-dominated cosmology. A first possibility to achieve this is by an early matter-dominated phase, where the decay of massive particles dilute the density of preexisting relics, thus allowing for very large DM masses [11]. Models where the DM abundance is set by annihilations into lighter mediators ('Secluded DM' [12]) realise the above picture if mediators are long-lived enough [13, 14], and allow for DM masses as heavy as a few EeV [15] ('Homeopathic DM' [15]). Importantly for LHAASO, they have been shown to give indirect signals at high-energy telescopes [15]. Early vacuum-energy domination offers a second possibility to dilute the overabundance of heavy DM, via the resulting period of exponential expansion of the universe ('Supercool DM' [16]).

1.3 One can even not depart from radiation domination and keep DM at equilibrium with a bath, provided it is a dark one. If that bath is at a temperature much smaller than the SM one, then the unitarity limit on the cross section can allow for DM masses as heavy as tens of PeV if DM freezes-out while relativistic [17].

1.4 Models of heavy DM have also been proposed where both DM is at chemical equilibrium with the SM bath, and radiation domination is not modified. They require the DM abundance to be set either by long chains of nearest-neighbour interactions [18], allowing for DM as heavy as $ \sim 10^{14} $ GeV, or by processes like DM $ \zeta^+ \to \zeta \zeta $ with ζ any slightly lighter particle, allowing for DM masses as large as $ \sim 10^9 $ GeV ('Zombie DM') [19].

2. Models at some high energy scale aimed at solving other issues of the Standard Model, that can accommodate DM at similar scales via one of the production mechanisms mentioned at point (1). "Traditional" constructions like supersymmetry (SUSY) and confinement constitute an only-partly explored playground in this respect, as we now sketch with a few examples. Models with a large SUSY-breaking scale can predict heavy DM in the form of a Gravitino, see e.g. [20] and references therein. Moving to lower SUSY-breaking scales, in gauge mediation new stable baryon-like composite particles, belonging to the strongly-interacting SUSY breaking sector, provide viable DM candidates [21–23]. SUSY also naturally accommodates dilution of relics via early matter domination, see e.g. [23]. Confining theories predict stable particles that are viable heavy DM candidates [24], and likewise can feature an early matter dominated phase associated to glueballs, which are naturally long-lived [25, 26]. If the confining phase transition is supercooled, then these models realise instead the picture of Supercool Dark Matter [27].

3. Models inspired by phenomenological puzzles in the 0.1-10 PeV energy range, notably the origin of the astrophysical neutrino flux discovered by IceCube [28–33]. If interpreted as a byproduct of DM, it is more straightforward to reconcile them with a decaying DM [34–51]. In these scenarios, DM is a non-thermal relic, whose abundance is set by particle decays, freeze-in, or inflationary-related mechanisms. Often, it is linked to the neutrino mass generation mechanism. It is important to remark that the sensitivity of LHAASO telescope, for large angular diffuse gamma-ray flux, will be an important discriminator to distinguish among astrophysical or heavy DM decay explanations for IceCube neutrinos [52, 53].

Similar to the WIMP paradigm, it seems to us that the PeV regime may establish itself as one of the new paradigms in the future, and therefore enjoy the synergy of theory and experiment that has made WIMP searches so prominent in the recent past. LHAASO offers a particularly important opportunity, in that respect. Note that many heavy DM models predict signals also at other experiments, from telescopes to colliders, from gravitational wavedetectors to cosmological surveys, thus widening the context of relevance for the LHAASO physics program.

The number of the particles of type i per unit energy E per unit volume per unit time injected at position x and time t via annihilation of the self-conjugated DM particles writes

For the decaying DM particles with lifetime $ \tau_{\rm{DM}} = \Gamma_{\rm{DM}}^{-1} $, the number of injected particles writes

In general, knowing the signal at the emission point is not enough, since we need the flux detectable at the Earth. For charged particles, one key difficulty is that they do not retain directionality due to deflections in interstellar magnetic fields: CR trajectories are similar to random-walks typically described via a diffusion-loss equation. Although we will not cover this channel in detail, note that it is of some interest for LHAASO as well, since at very least energetic $ e^\pm $ final states are responsible for secondary gamma-rays, notably via Inverse-Compton up-scattering of background photons, i.e. galactic radiation fields and the ubiquitous CMB.

Even limiting oneself to prompt photons, it is important to notice that absorption is significant at LHAASO energies. In particular, LHAASO will be the first experiment to be sensitive to gamma-ray absorption even for Galactic sources. If we define the optical depth to photon-photon pair production τ, the differential flux (number of particles per unit time, energy, surface, and solid angle) writes as an integral of the above term Q along the line of sight (below, b and l indicate latitude and longitude in Galactic coordinates, respectively)

For a spherically symmetric DM distribution (as typically assumed for the halo of our Galaxy) one has $ \rho_{\rm{DM}}(s,b,l) = \rho_{\rm{DM}}[r(s,b,l)] $ where

and $ R_{\odot} = 8.249\pm0.045_{\rm{syst}}\pm0.009_{\rm{stat}}\, $ kpc [54] is the distance of the Sun from the Galactic Center. A benchmark density profile is the Navarro-Frenk-White (NFW) density profile

with typical range of the scale radius $ 10\,{\rm{kpc}}\lesssim r_0\lesssim 20\, $ kpc and the $ \rho_0 $ corresponding to a local abundance of $ \rho_{\rm DM}(R_{\odot})\simeq0.4\, $ GeV/cm³ (for recent determinations, see e.g. [55, 56].) The profile becomes more uncertain towards the inner Galaxy. Yet, along most directions, the uncertainties on the flux for the case of decaying DM are within a factor two.

The best sensitivity for indirect searches of DM is naively expected by telescopes with the largest "grasp" $ G = A\Omega $ [53, 57–59], with A being the area and Ω the solid angle field of view. The current HAWC and, especially, the forthcoming LHAASO are the detectors with the largest grasp in the very-high-energy γ-band and are therefore well suited for the DM search.

Dwarf spheroidal satellites (dSphs), which are large Galactic DM substructures with high DM densities, are ideal targets for detecting the gamma-ray signals induced by DM annihilation or decay. This search is almost background-free, because of the lack of gamma-rays from astrophysical processes in dSphs. Since LHAASO has strong background rejection power (~1%) and large field of view (FOV) (~ 2 sr), it is promising for LHAASO to investigate the properties of massive DM particles via the gamma-ray observations of dSphs. In this subsection, we discuss the LHAASO sensitivities to gamma signals induced by DM annihilation from dSphs, which are investigated in Ref. [60].

For DM annihilation in a point-like source, the expected gamma-ray flux in an energy bin is given by

For a certain annihilation channel, the initial photon spectrum $ \dfrac{{\rm{d}}N_{\gamma}}{{\rm{d}}E} $ can be obtained using the $\tt{PPPC4DM}$ package [61, 62]. The J-factor is defined as the integral of the DM density squared along the line of sight within a solid angle of $ \Delta \Omega = 2\pi (1-\cos \alpha_{\rm{int}}) $. Using the results of the kinematic observation of dSphs, their J-factor can be derived from Jeans analysis [63–65].

19 dSphs with large J-factors are considered in the analysis. These dSphs are located in the LHAASO FOV with favored declination angles. Compared with the HAWC investigation [66], four more dSphs, namely Draco II, Leo V, Pisces II, and Willman 1, are taken into account. Notice that Refs. [67, 68] provide two sets of J-factors with different choices for the integration angle $ \alpha_{\rm{int}} $. One set is derived with a constant $ \alpha_{\rm{int}} = 0.5^\circ $, while the other is derived with varying $ \alpha_{\rm{int}} = \theta_{\max} $. $ \theta_{\max} $ is the maximum angular radius of a certain dSph and can be determined by the observation of the outermost member star. To get a large signal-to-background ratio, the J-factor of each dSph is taken to be that derived with a smaller $ \alpha_{\rm{int}} $ as $ \min\{\theta_{\max},0.5\} $.

The main backgrounds in the analysis result from the mis-identified cosmic-ray particles. The expected background number in an energy bin can be calculated as

where $ \dfrac{{\rm{d}} \Phi_p}{{\rm{d}} E} $ is the flux of primary cosmic-ray proton and is assumed to be a single power-law, the observational time T is taken to be one year, and $ E_{\max}/E_{\min} $ is assumed to be 3. In order to include the contributions from heavier nuclei in the primary cosmic-rays, an additional factor $ \zeta_{cr} = 1.1 $ is introduced. The background number is calculated within a cone around the direction of dSph with a solid angle $ \Delta\Omega = 2\pi\times[1-\cos(\max\{\alpha_{\rm{int}}, \theta_c \})] $, where $ \theta_c $ is the angular resolution of LHAASO. With increasing photon energy, $ \theta_c $ varies from $ 2^\circ $ to $ 0.1^\circ $. The effective area of LHAASO $ A_{\rm{eff}}^{p} $ is taken from Ref. [69], which depends on energy and zenith angle $ \theta_{\rm{zen}}(t) $. In order to show the visibility of dSphs, the values of $ r_{\rm{eff}} $, which is defined as the ratio of the observation time with $ \theta_{\rm{zen}}< 60^\circ $ to total time, are listed in Table 1.

Source	RA. (deg)	DEC. (deg)	$ r_{\rm{eff}} \;{\rm{ (year)}}$	$ \theta_{\max}\;{\rm{ (deg)}}$	$ \log_{10}J_{\rm{obs}} $ ($ \text{GeV}^{2}{\rm{cm}}^{-5} $)
$ {\rm{Bo\ddot{o}tes \; I}}$	$ 210.02 $	$ 14.50 $	$ 0.352 $	$ 0.47 $	$ 18.2\pm0.4 $
Canes Venatici I	$ 202.02 $	$ 33.56 $	$ 0.398 $	$ 0.53 $	$ 17.4\pm0.3 $
Canes Venatici II	$ 194.29 $	$ 34.32 $	$ 0.399 $	$ 0.13 $	$ 17.6\pm0.4 $
Coma Berenices	$ 186.74 $	$ 23.90 $	$ 0.377 $	$ 0.31 $	$ 19.0\pm0.4 $
Draco	$ 260.05 $	$ 57.92 $	$ 0.442 $	$ 1.30 $	$ 18.8\pm0.1 $
Draco II*	$ 238.20 $	$ 64.56 $	$ 0.451 $	−	$ 18.1\pm2.8 $
Hercules	$ 247.76 $	$ 12.79 $	$ 0.348 $	$ 0.28 $	$ 16.9\pm0.7 $
Leo I	$ 152.12 $	$ 12.30 $	$ 0.346 $	$ 0.45 $	$ 17.8\pm0.2 $
Leo II	$ 168.37 $	$ 22.15 $	$ 0.372 $	$ 0.23 $	$ 18.0\pm0.2 $
Leo IV	$ 173.23 $	$ -0.54 $	$ 0.303 $	$ 0.16 $	$ 16.3\pm1.4 $
Leo V	$ 172.79 $	$ 2.22 $	$ 0.314 $	$ 0.07 $	$ 16.4\pm0.9 $
Pisces II$ ^{\star} $	$ 344.63 $	$ 5.95 $	$ 0.327 $	$ - $	$ 16.9\pm1.6 $
Segue 1	$ 151.77 $	$ 16.08 $	$ 0.357 $	$ 0.35 $	$ 19.4\pm0.3 $
Sextans	$ 153.26 $	$ -1.61 $	$ 0.299 $	$ 1.70 $	$ 17.5\pm0.2 $
Triangulum II*	$ 33.32 $	$ 36.18 $	$ 0.403 $	−	$ 20.9\pm1.3 $
Ursa Major I	$ 158.71 $	$ 51.92 $	$ 0.432 $	$ 0.43 $	$ 17.9\pm0.5 $
Ursa Major II	$ 132.87 $	$ 63.13 $	$ 0.449 $	$ 0.53 $	$ 19.4\pm0.4 $
Ursa Minor	$ 227.28 $	$ 67.23 $	$ 0.455 $	$ 1.37 $	$ 18.9\pm0.2 $
Willman 1*	$ 162.34 $	$ 51.05 $	$ 0.430 $	−	$ 19.5\pm0.9 $

Table 1.  The astrophysical properties of the 19 selected dSphs. The columns for each dSph denote the name, right ascension (RA.), declination (DEC.), effective time ratio ($ r_{\rm{eff}} $), maximum angular radius (${\rm{\theta_{max}}} $), and the median value and standard deviation of the J-factor. The J-factor and $ {\rm{\theta_{max}}} $ of most of the dSphs are taken from Ref. [67]. The J-factors of four dSphs marked with asterisks are not provided by Ref. [67]; they are taken from Ref. [68]. From [60].

Since the main backgrounds arise from cosmic-ray particles, it is crucial to improve the $ \gamma/p $ discrimination in the analysis. Ref. [70] provides the energy-dependent quality factor $ Q\equiv\varepsilon_{\gamma}/\sqrt{\varepsilon_{p}} $ for WCDA, where $ \varepsilon_{\gamma} $ and $ \varepsilon_p $ are survival ratios of gamma-rays and primary protons, respectively. For $ \varepsilon_{\gamma}$~50% and $ E>0.6\,\text{TeV} $, $ \varepsilon_{p} $ varies within a range of 0.04%-0.11%. In Here $ \varepsilon_{p} $ is taken to be 0.278% for $ \varepsilon_{\gamma} $~40.13% for a conservative calculation.

In order to calculate the LHAASO sensitivity, the mimic observations with only backgrounds are assumed. The event count N in each mimic observation is randomly generated around the expected background number B by Poisson sampling. The impact of the statistical uncertainty of the J-factor is also considered as Refs. [71, 72]. For one mimic observation, the likelihood is defined as

where i and j denote the i-th energy bin and j-th dSph, respectively, $ {\cal{L}}(S|B, N) $ is the Poisson distribution with expected signal number from DM annihilation S, ${\rm{log}}_{10}(J_{{\rm obs},{\it j}}) $ is the observed mean value of the J-factor, and $ \sigma_{j} $ is the corresponding standard deviation. The value of ${\rm{log}}_{10}(J_{\it j}) $ is adjusted to maximize $ {\cal{L}}_{j} $ for given $ \langle\sigma v\rangle $ and $ m_{\rm{DM}} $. In order to derive the upper-limit on $ \langle \sigma v\rangle $ at 95% C.L., it is required that the log-likelihood including the contribution of DM annihilation with increasing $ \langle \sigma v\rangle $ decreases by 2.71/2 from its maximum. Furthermore, a combined likelihood $ {\cal{L}}^{\rm{tot}} = \prod_{j}{\cal{L}}_{j} $ can be used to derive an improved upper-limit on $ \langle \sigma v\rangle $ in the joint analysis with many dSphs. As an example, the individual sensitivity for single dSph and the combined sensitivity for all the 19 dSphs for the $ b\bar{b} $ annihilation channel from one mimic observation are shown in Figure 1.

Figure 1.  (color online) The expected one-year sensitivities to the DM annihilation cross section $\left\langle {\sigma v} \right\rangle $ at 95% C.L. for the ${\rm{b\bar{b}}}$ annihilation channel in one mimic observation. The solid red line represents the combined sensitivity derived in the joint analysis for the 19 selected dSphs; while other lines represent the sensitivities for single dSph. From [60].

It is found that the sensitivities derived from the observations of three dSphs, namely Segue 1, Ursa Major II, and Triangulum II, are much better than other selected dSphs. Therefore, the combined sensitivity in the joint analysis is dominantly determined by these three dSphs, which have large J-factors and favorable locations in the LHAASO FOV. Note that Triangulum II has almost the largest J-factor among 19 dSphs, but the statistical uncertainty of its J-factor is also large. This is because that Triangulum is an ultra-faint dSph and the corresponding kinematic data is not sufficient. As a consequent, although Triangulum is located very close to the center of LHAASO FOV, it cannot utterly determine the combined sensitivity. This means that the uncertainty of the J-factor should be including in the analysis; otherwise, the sensitivity given by some dSphs may be overestimated.

Since the event number of very high energy photon is small, a series of mimic observations should be performed to reduce the impact of statistical fluctuation. The median values and the two-sided 68% and 95% containment bands of the combined sensitivities derived from 500 mimic observations for five typical DM annihilation channels, including $ b\bar{b} $, $ t\bar{t} $, $ \mu^{+}\mu^{-} $, $ \tau^+\tau^- $, and $ {\rm{W^+W^-}} $, are shown in Figure 2. In order to compare with the available limits from other gamma-ray observations of dSph, the results of the HAWC combined limit [66], Fermi-LAT combined limit [73], HESS combined limit [74], VERITAS Segue 1 limit [75], and MAGIC Segue 1 limit [76] are slso shown.

Figure 2.  (color online) The combined one-year LHAASO sensitivities at 95% C.L. for five DM annihilation channels, including $ {\rm b}\bar{\rm{b}} $, $ {\rm{t\bar{\rm{t}}}} $, $ \mu^{+}\mu^{-} $, $ \tau^+\tau^- $, and $ {\rm{W^+W^-}} $. The red solid lines represent the median values. The yellow and green bands represent the corresponding two-sided 68% and 95% containment bands, respectively. The HAWC combined limits [66], Fermi-LAT combined limit [73], VERITAS Segue 1 limit [75], HESS combined dSph limit [74] and MAGIC Segue 1 limit [76] are also shown for comparison. From [60].

Since the initial photon spectrum of the $ \tau^{+}\tau^{-} $ annihilation channel is hard, the sensitivity for this channel, which can reach $ \sim 10^{-24}\; {\rm{cm}}^{3}\; s^{-1} $ for $ m_{\rm{DM}}>1 $ TeV, is better than other channels. For the $ \tau^+\tau^- $, $ {\rm{W^+W^-}} $, and $ b\bar{b} $ channels, the LHAASO sensitivities are better than the current limits for $ m_{\rm{DM}} $ larger than $ \sim2 \;{\rm{TeV}} $, $ \sim3 \;{\rm{TeV}} $, and $ \sim8 \;{\rm{TeV}} $, respectively. For the $ \mu^{+}\mu^{-} $ and $ t\bar{t} $ channels, LHAASO has good sensitives to explore the DM signals for from $ m_{\rm{DM}} $ in the range of $ \sim $ 1 -100 TeV.

The similar analysis procedures can be applied for decaying DM. The flux of gamma-ray signals from DM decays depend on the lifetime of DM $ \tau_{\rm {DM}} $ and the D-factor, which is the integral of the DM density along the line of sight. The sensitivities to $ \tau_{\rm {DM}} $ for five DM decay channels from the LHAASO gamma-ray observation of dSphs are investigated in Ref. [77]. 19 dSphs within the LHAASO FOV are studied in the individual and combined analyses. Two dSphs, namely Draco and Ursa Major II, would significantly affect the combined sensitivity, due to their large D-factors and suitable locations in the FOV of LHAASO. For $ m_{\rm{DM}}\sim 100 $ TeV, the LHAASO sensitivities to $ \tau_{\rm {DM}} $ can reach $ \sim 10^{-27} $ s.

In this subsection, we focus on the gamma-ray signals from DM decays in the Galactic halo. For particles with negligible absorption, such as neutrinos, a further, quasi-isotropic contribution to the flux is due to annihilations or decays in the whole universe. This term is in general negligible for annihilating DM (unless one assumes rather extreme DM halo clumpiness), but it is definitely comparable to the Galactic term for decaying DM, and should be taken into account. In particular, this component reduces the scale of variations of the signal across the sky in the neutrino channel. For gamma-rays in the range of LHAASO, however, the extragalactic sky is fully opaque (see e.g. [78]), hence the extragalactic DM contribution is degraded in energy below the pair-production threshold on the EBL, at $ E\lesssim 10^2-10^3, $ GeV. Still, the diffuse gamma ray background measured by Fermi-LAT provides an upper limit to the gamma-ray flux at TeV and PeV energies.

Different strategies are possible for the search of the DM decay signal. The analysis by the HAWC collaboration [79] has adopted an approach in which a signal from the direction around the Galactic Center (more precisely, the region of the Fermi Bubble) is searched for, and the rest of the sky is considered for the background estimate. An alternative possibility is to search for a somewhat weaker (by a factor of two, on average) signal, but extending across the entire sky. An advantage of the latter approach is the larger exposure available for the full-sky search, while a disadvantage is the stricter requirements on the charged-particle vs. gamma-ray separation, due to the modest if not negligible angular variation of the signal. A simple rescaling suggests that a full-sky exposure would provide an increase of the DM signal-to-noise ratio by a factor of $ \simeq 2 $ on one-year observation time span, compared to the Fermi Bubble region exposure, in spite of the lower average flux. The use of the full (or large) sky exposure, rather than of limited sky region around the Galactic Centre, is important also in the view of uncertainties on the Galactic diffuse γ-emission unrelated to the DM decay flux. This diffuse emission provides a background on top of which the DM decay signal is searched for. Even if it is possibly sub-dominant compared to the residual charged particle background in γ-telescopes, it might still be stronger than the DM decay signal.

Intriguingly, the birth of high-energy neutrino astronomy provides a benchmark region in parameter space to search for a possible DM decay signal. The IceCube experiment, completed in 2011, continues observing a flux of high energy ($ \gtrsim 10 $ TeV) neutrinos significantly in excess with respect to the expected background from atmospheric neutrinos and muons [28–33]. The source(s) of these neutrinos is yet unknown, although based on their almost uniform angular distribution an extragalactic origin or a galactic halo origin is favored. Directional analyses with various classes of astrophysical objects and catalogs are not showing any correlation leading to the conclusion that the contribution of well-known objects, such as blazars, to the observed diffuse neutrino flux is $ \lesssim{\cal{O}} $(10%) [80]. In this context, it has been quite natural to consider unconventional sources for these neutrinos. Also, since neutrinos provide for the first time a window on the 0.1-10 PeV Universe, it may not be so bizarre that new classes of sources can pop up. The potential to answer long-standing problems such as the nature of DM by investigating this energy regime has only recently been entertained.

A decaying DM scenario has gained some attention, mainly due to its minimal assumptions and its testability in future gamma-ray experiments. Interestingly, not only the first PeV-scale events discovered [34], but the whole observed flux of neutrinos by IceCube can be interpreted in this scenario [35], although a multi-component flux arising from both the conventional astrophysical sources and DM also has been investigated, starting from [38]. In a phenomenological approach, the properties of the required DM particle can be deduced from a fit to the neutrino data. In this case, the free parameters are the decay lifetime, the mass and the branching ratios of the DM decay to various standard model particles. The ballpark lifetime is $ \sim10^{27} -10^{28} $ s and the mass has to be $ \gtrsim $ few PeV in order to interpret the highest energy observed events in IceCube (obviously, assuming the multi-component neutrino flux, the DM mass can take any value in our range of interest $ \gtrsim10 $ TeV). The highest energy events are typically accounted for via 'hard' leptonic final states, while lower energy events are fitted via soft channels including e.g. gauge bosons and quarks. Part of the flux can also be accommodated via some astrophysical component.

For any decay channel of PeV-scale DM particles explaining or contributing to the IceCube neutrino flux, gamma rays are unavoidably associated decay products, and their Galactic fraction can be observed by LHAASO. The following processes contribute to the expected gamma ray flux at Earth: i) A prompt flux is at very least due to the electroweak corrections ii) A secondary flux is induced by the unavoidable prompt (as well as secondary) charged leptons, via the Inverse-Compton process onto the CMB and star-light which lead to a spectrum of high energy gamma rays where the Galactic part of it contribute to the total flux. Of course the exact spectral shape of the flux depends on the magnitude and profile of the magnetized halo in our Galaxy, which are yet not known very well.

Figure 3, updated from Ref. [52], shows the spectrum of gamma ray yield from the decay of DM with mass 4 PeV and final state branching ratios given by: $ \ell^{\pm} W^{\mp}:\stackrel{(-)}{\nu_{\ell}} Z:\stackrel{(-)}{\nu_{\ell}} h = 1:2:2 $. The solid curves show the prompt flux accounting for gamma-ray absorption; different colors represent different directions in the sky. Even this Galactic flux suffers from absorption due to the pair production on CMB and star-light, with the suppression reaching ~70% for the Galactic center line of sight. Dashed curves show the flux due to inverse-Compton photons, for various assumptions for the constant halo magnetic field, $ B_{\rm{halo}} $, possibly pervading the thick diffusive halo of the Galaxy up to large distances. While uncertain, it is particularly important in the range above 100 TeV. Note how upper bound from two-decade old experiments CASA-MIA [81] and KASCADE [82] are within one order of magnitude of the expected flux (even less, if they had been sensitive to regions closer to the Galactic center), while LHAASO should provide a definite test of this scenario (the black line indicates its 1 yr nominal sensitivity).

Figure 3.  (color online) The gamma-ray flux from DM decay from various directions, with $m_{\rm{DM}} = 4$ PeV and $\tau_{\rm{DM}} = 10^{28}$ s, and branching ratios reported in the text. The solid colored curves show the prompt flux, including the absorption of gamma-rays; different colors represent different directions in the sky. The dashed curves show the IC flux, for various assumptions for the constant halo magnetic field, $B_{\rm{halo}}$, possibly pervading the thick diffusive halo of the Galaxy up to large distances. The green and brown bar lines show the upper bound on gamma-ray flux from CASA-MIA [81] and KASCADE [82], respectively. The black line is an indicative 1 yr LHAASO sensitivity.

Meaningful bounds can also be obtained thanks to diffuse gamma-ray data by Fermi-LAT in the GeV band, see e.g. [83], exploiting the cascading effect on the extragalactic part of the flux previously mentioned. However, such constraints are rather indirect, depending on the datasets used, the final state considered, and the different assumptions for the contributions to the astrophysical background. LHAASO would allow one to probe the scenario directly and unambiguously, achieving great sensitivity also to sub-leading DM contributions. This was explicitly illustrated in Ref. [84], where the authors estimated the LHAASO sensitivity reach for the decaying DM search, following the approach of Ref. [79]: In each energy bin they compared the DM decay flux levels for different values of $ m_{\rm DM},\tau_{\rm DM} $ with the residual charged particle background levels and calculate by how much is the $ \chi^2 $ of the fit of the signal+background data is inconsistent with the backgorund-only model in all energy bins. In this way they found the minimal detectable DM decay flux as a function of the DM mass for the model of Ref. [84] of DM decaying into quark-antiquark pair, in turn converted into a maximal measurable DM decay time. The results are shown in Fig. 4. It clearly illustrates how LHAASO will explore DM lifetimes up to $ \tau_{\rm DM}\sim 3\times 10^{29} $ s over a wide DM mass range $ m_{\rm DM}>10 $ PeV. In the mass range $ 10 $ TeV$ <m_{\rm DM}<10 $ PeV LHAASO will provide a factor of 3-to-10 improvement of sensitivity compared to HAWC. In any case, LHAASO will fully test models where a non-negligible fraction of the IceCube astrophysical neutrino flux is generated by DM decays.

Figure 4.  (color online) Sensitivity of LHAASO for the measurement of dark matter decay time (for DM decaying into quarks). Yellow band shows the range of decay times for which DM decays give sizeable contribution to the IceCube neutrino signal [84]. Blue and grey shaded regions show the existing bounds imposed by HAWC [79] and ultra-high-energy cosmic ray experiments [85]. and dashed cureves are from the HAWC search of the DM decay signal in the Fermi Bubble regions [79]. From [53].

Let us remind that the LHAASO sensitivity to new physics is not limited to DM candidates. As an example, light pseudo-Goldstone bosons (PGBs) can efficiently be converted by the interaction with high energy photons that belong to a magnetic field background, according to the Primakoff effect [86, 87]. PGBs arise ubiquitously in physics, in particular from several possible spontaneously broken global symmetries beyond the SM. Notable examples are provided by the QCD axion, the familon/archion [88, 89], the majoron [90, 91], the baryo-majoron [92] and and the exoticon [93]. All PGBs interact with the SM particles through terms of the Lagrangian that are suppressed by an energy scale related to the spontaneous symmetry breaking and the decay constant f, as $ f^{-1}J^{\mu}\partial_{\mu}a $, J standing for the Noether current related to the spontaneously broken global symmetry, and a denoting the light PGB field.

Similar light particles are also suggested from string-inspired phenomenology: moduli scalar fields from string compactifications, as well as anomalous Peccei-Quinn symmetries $ U_{P.Q.}(1) $ from intersecting D-branes models, provide a rich spectrum of different axion-like particles (ALPs) beyond QCD axion mass and couplings [94–96]. PQ symmetries, generated from string theory, are explicitly violated by worldsheet instantons, brane instantons, new gauge instantons and gravitational instantons, which are beyond QCD instantonic effects [94–96].

The QCD axion is proposed as a solution of the strong CP problem, dynamically shifting out the θ-parameter in $ {\cal{L}}_{\theta} = -\theta(\alpha_{s}/8\pi)G^{\mu\nu a}\tilde{G}^{a}_{\mu\nu} $, where $ G^{a}_{\mu\nu} $ is the gluon field strength and $ \tilde{G}^{\mu\nu a}\equiv \epsilon^{\mu\nu\lambda\rho}G_{\lambda\rho}^{a}/2 $. The QCD axion is related to a $ U_{P.Q.}(1) $, which is broken through anomalous triangle couplings to gluons: $ {\cal{L}} = (a/f-\theta)(\alpha_{s}/4\pi)G^{\mu\nu a}\tilde{G}^{a}_{\mu\nu} $. Non-perturbative topological fluctuations of the gluon field induce a periodic potential for a-field, the minimum of which is set at $ a = f\theta $, screening the CP violating θ-term [97–100]. It is worth to note that axions may be also thought as composite states of neutrinos, through new fundamental fifth force interactions [101] or chiral symmetry breaking topological terms sourced by non-perturbative gravitational effects [102, 103].

The ALP-photons interaction can be generated by anomalous triangle diagrams as $ {\cal{L}}_{a\gamma\gamma} = g_{a\gamma\gamma}aF_{\mu\nu}\tilde{F}^{\mu\nu} $ /4, where $ F^{\mu\nu} $ is the electromagnetic field strength. In term of electric and magnetic fields, $ {\cal{L}}_{a\gamma\gamma} $ can be rewritten as $ g_{a\gamma\gamma} a\vec{E}\cdot \vec{B} $. This interaction implies that an external magnetic field background may source photon-ALPs mixings, through the aforementioned Primakoff effect [86, 87]. This motivates searches for ALPs in laboratory experiments with high magnetic fields, including the CAST experiment in CERN facility [104, 105].

On the other hand, possible hints of $ a-\gamma $ oscillations may be observed from VHE gamma-ray sources, triggered by cosmic magnetic fields along the sight line [106–110]. An amplifier for this effect would be provided by the distance of sources, such as blazars. These are active galactic nuclei (AGN) with gamma beam pointing toward our line of sight, and representing the most distant and energetic gamma-ray sources observed. This certainly individuates blazars as natural candidates for the searches of ALP's effects. The mean free path of VHE photons above $ 100\;{\rm{GeV}} $ energy is related to the interactions with radiation backgrounds and pair production $ \gamma\gamma\rightarrow e^{+}e^{-} $. It is well known that VHE gamma rays have a non-negligible probability of disappearing into electron-positron pairs, because of the annihilation with extragalactic background light (EBL) [111–114].

For gamma-ray energies below $ 100\, {\rm{GeV}} $, the pair production effect is completely negligible, since the mean free path length $ \lambda_{\gamma}(E) $, with E the VHE photon energy, is comparable with the Hubble radius. On the other hand, for $ E>100\, {\rm{GeV}} $, the EBL effects quickly increase and $ \lambda_{\gamma}(E) $ rapidly decreases with the energy. The $ \gamma\rightarrow a \rightarrow \gamma $ oscillations can occur in the intergalactic space from the source to the Earth, turning the average $ \lambda_{\gamma} $ into an higher value than expected, competitively with electron-positron pair production [106–108].

Another possibility is that $ \gamma-a $ conversion happens inside or near blazar sources, and axions are converted back into photons inside the Milky Way [109, 110, 115–120]. This mechanism induces VHE gamma path length to be become larger than expected. Thus the Universe would appear more transparent in the very high energetic spectrum. A hardening of the observed blazar spectra, parametrically controlled by the $ a\gamma\gamma $ coupling, would provide a hint of a photon-ALP oscillation in $ 1-100\,{\rm{TeV}} $ range. The future LHAASO measurements of blazar gamma-ray spectra in this energy range will impose stronger limits on the ALP parameters than current results. For the blazars with variations in the spectra, such as Markarian 421, the combined analysis of observations in different periods can future improve the constraints [120].

The ALP interpretations of such observations depend on some astrophysical features. If the blazar is located at a rich galaxy cluster, the $ \gamma-a $ conversion could happen in the turbulent inter cluster magnetic field of $ {\cal{O}}(1)\mu {\rm{G}} $. Since there is no concrete model for this magnetic field, a series of random configurations is introduced in the analysis [115]. For blazars that are not resides in such environment, the $ \gamma-a $ conversion dominantly happen in the blazar jet magnetic field, whose strength can be determined by the multi-wave observations [115, 118]. The ALP interpretations also depend on the EBL model. The predicted blazar spectra with the ordinary EBL models such as given by Ref. [78], can well fitted most of current measurements. Thus the measurements of blazar gamma-ray spectra impose constraints on the ALP parameters. On the other hand, if there exist some exotic EBL contents at $ \sim {\cal{O}}(1)\mu {\rm{m}} $ as reported by some recent measurements such as CIBER [121] the additional EBL absorption effect would reduce the observed VHE photons. In this case, the $ \gamma-a $ oscillation compensates the lost of VHE photons and would provide better fitting to the measurements [116, 119].

The similar analysis can be applied for the VHE gamma-ray sources in the Galaxy [122, 123]. An advantage of such analysis is that the model of the Galactic magnetic field is more clear than those in the extragalactic space. The ASγ observation of sub-PeV gamma-rays from the Crab Nebula has been used to set limits on the ALP parameter space [123]. The future LHAASO detection for the photons with higher energies and larger statistic will improve the current limits.

Line emission from ALPs decaying everywhere across the Universe also contributes to the EBL [133]. Superposition of the line emission accumulated at different redshifts leads to a broad "bump-like" feature in the EBL spectrum with the width $ \Delta E\sim E $. Such feature could be directly measured by dedicated observations probing the EBL. Such measurements [121, 134–137] are, however, challenging in the eV photon energy range because of the presence of strong Zodiacal light background. The direct measurements typically find an excess EBL flux above the estimate from direct galaxy counts in the near-infrared [138–141]. Recent measurements by CIBER [121] estimate the excess flux at the level 29 to 42 $ {\rm{nW/m}}^2/sr $, depending on assumptions about Zodiacal flux model. Recent results from AKARI [136] between $ \sim 2\;\mu {\rm{m}} $ and $ \sim 5\;\mu {\rm{m}} $ also demonstrate an increase in EBL intensity at shorter wavelengths.

Attenuation of gamma-ray signal from distant sources by the effect of pair production on the EBL photons provides an alternative way of probing the EBL [142–146]. This method has been used to set limits on the normalisation of the EBL, under certain assumptions about its spectral shape. However, such an analysis does not account for possible unforeseen spectral features in the EBL, such as e.g. the ALP dark matter decay bump. On the other hand, presence of the bump leads to a characteristic dip in the observed gamma-ray spectra of distant sources, thus making this feature potentially observable.

In a recent study [147] combined Fermi/LAT and Cherenkov telescopes gamma-ray spectra of a large number of blazars were used to derive tight constraints on possible EBL bump features by fitting the gamma-ray spectra with models which include additional attenuation of the gamma-ray flux by absorption on the EBL bump. These constraints were then used to explore the parameter space of ALPs. Fig 5 shows constraints on the mass and two photon coupling constant obtained with current data and predictions for the sensitivity of LHAASO.

Figure 5.  (color online) Constraints of ALP parameters from laboratory and astrophysical probes. Yellow band and black solid line corresponds to the QCD axion models [124–127]. Blue region is the limit impost by CAST non-observation of solar axions [128]. Grey vertical band comes from high-resolution spectroscopy of intracluster medium in by optical telescope [129]. Black dashed line is the EBL limit estimate of [130]. Green horizontal line is the limit imposed by non-observation of excess energy loss in Horizontal Branch stars [131]. Light grey shaded region shows 95% confidence level excluded range of parameters from non-observation of absorption feature in the Fermi/LAT+VERITAS spectrum of 1ES 1218+304 [132]. Pink hatched and region show 68% confidence level preferred range of the parameters obtained by the spectral analysis of 1ES 1218+304. Red region is the same as pink but for the simulated CTA spectrum. Brown dashed and solid lines show the precision of the ALP mass and two-photon coupling which can be reached with 1- and 5-years of observations with LHAASO (95% confidence level).

The existing data from Imaging Atmospheric Cherenkov Telescopes and the direct measurements of the EBL flux indicate that a small bump in the EBL spectrum, which might be due to the ALP decays, exists in the wavelength range around 1 micron [132]. The range of ALP parameters which corresponds to this hint is shown by the pink hatched region in Fig. 5. The "dip" in the gamma-ray spectra imprinted by the ALP-induced EBL bump is only marginally detectable by the current generation telescopes. Next generation telescope CTA will be able to test the hypothesis of the EBL excess produced by the ALP decays and pinpoint the parameters of ALP particles. Red colored region in Fig. 5 shows the precision of measurement of ALP paerameters which can be achieved with CTA [132].

Non-observation of axion flux from the Sun by CAST experiment constrains $ g_{a\gamma\gamma} $ from above in the mass range below 1 eV, as shown in Fig. 5 [128]. The CAST limit will be improved by the next-generation facility IAXO [148], also shown in Fig. 5. Another strong constraint on the ALP parameters is imposed by non-observation of the effect of energy loss through emission of ALPs on stellar evolution [149, 150]. Recent constraints on two-photon coupling [131, 151] are shown by the horizontal line in Fig. 5, stem from non-detection of excessive energy loss through ALP emission which would affect evolution of the Horizontal Branch (HB) stars in globular clusters. Similarly to most of the astrophysical constraints, the stellar evolution argument relies on the assumption of validity of complex stellar evolution models controlled by large number of parameters which are not directly measured [152, 153]. Direct searches for ALP decays into two photons were performed using highresolution spectroscopy of galaxy clusters in the range of masses between 4.5 and 7.7 eV [129]. Non-observation of non-identified spectral emission lines from ALP decays imposes a constraint on $ g_{a\gamma\gamma} $ shown by grey shading in Fig. 5.

To estimate the sensitivity reach of LHAASO for the ALPs parameters we performed the same analysis as proposed in [132]. Six sources, which satisfy the criteria from [132] and located in LHAASO field of view were selected: Markarian 501, Markarian 421, 1ES 1218+304, 1ES 1215+303, 1ES 1959+650, 1ES 1011+496. Using gamma-ray sensitivity, effective area and gamma/hadron separation efficiency we estimated background rate and simulated spectra for the blazars in the sample for 1- and 5-years of observations with LHAASO in the energy range between 200 GeV and up to 20 TeV. One can see that 5-years LHAASO sensitivity is sufficient to explore potentially interesting area of ALPs masses around 1 eV and $ g_{a\gamma\gamma} \sim 10^{-10}\, \rm{GeV}^{-1} $ indicated from previous observations. In the case of absence of the signal LHAASO will improve current limitation on $ g_{a\gamma\gamma} \sim 10^{-10}\, \rm{GeV}^{-1} $ by a factor of 3-5 in the mass range from 1 eV to 5 eV. From Fig. 5 one can see that LHAASO will have enough sensitivity to cross-check the CTA measurement of the ALP parameters (if the existence of a bump in the EBL spectrum will be confirmed).