Constraints on metastable superheavy dark matter coupled to sterile neutrinos with the Pierre Auger Observatory

Dark matter particles could be superheavy, provided their lifetime is much longer than the age of the universe. Using the sensitivity of the Pierre Auger Observatory to ultra-high energy neutrinos and photons, we constrain a specific extension of the Standard Model of particle physics that meets the lifetime requirement for a superheavy particle by coupling it to a sector of ultra-light sterile neutrinos. Our results show that, for a typical dark coupling constant of 0.1, the mixing angle $\theta_m$ between active and sterile neutrinos must satisfy, roughly, $\theta_m \lesssim 1.5\times 10^{-6}(M_X/10^9~\mathrm{GeV})^{-2}$ for a mass $M_X$ of the dark-matter particle between $10^8$ and $10^{11}~$GeV.

Dark matter particles could be superheavy, provided their lifetime is much longer than the age of the universe.Using the sensitivity of the Pierre Auger Observatory to ultra-high energy neutrinos and photons, we constrain a specific extension of the Standard Model of particle physics that meets the lifetime requirement for a superheavy particle by coupling it to a sector of ultra-light sterile neutrinos.Our results show that, for a typical dark coupling constant of 0.1, the mixing angle   between active and sterile neutrinos must satisfy, roughly,   ≲ 1.5 × 10 −6 (  /10 9 GeV) −2 for a mass   of the dark-matter particle between 10 8 and 10 11 GeV.
Despite its countless successes, the Standard Model (SM) of particle physics is known to be incomplete.Missing elements include, among others, dark matter (DM) and neutrino masses.* spokespersons@auger.org; http://www.auger.orgA minimal extension to the SM is to add right-handed neutrinos inert to SM interactions ("sterile" neutrinos) that mix with left-handed ones active under weak interactions ("active" neutrinos) via the Brout-Englert-Higgs mechanism.Sub-eV mass eigenvalues result from mixing angles made sufficiently small by highly-enough massive sterile neutrinos (seesaw mechanism [1,2]).Other extensions may use additional sterile neu-trinos, as their number is in general largely under-constrained.This is the case of the extension presented in [3], considered in this study as a beyond-standard-model (BSM) benchmark.In this model, which includes DM and neutrino masses, a second sterile neutrino is assumed to render metastable a superheavy particle referred to as  with mass   and lifetime   , the decay products of which can leave unique signatures in the data of the Pierre Auger Observatory [4].
A superheavy particle that is metastable attracts particular attention if it is to be a viable DM candidate.The constraints imposed on   are indeed demanding:   ≳ 10 22 yr, see e.g.[5][6][7][8][9][10][11] for recent estimates in various channels.Compliance with these constraints indeed calls for adjusting a reduced coupling constant   down to a tiny level and, simultaneously, the multiplicity of the final state to a large value [12,13].This is in general challenging for theoretical constructions, unless the decay rules are based on non-perturbative effects [14,15].In the BSM benchmark [3], however, the lifetime does not resort to such a fine tuning.The DM particle  interacts only with sub-eV and superheavy (10 12−14 GeV) sterile neutrinos, of masses   and   respectively, via Yukawa couplings   and   .In the mass-eigenstate basis, neutrinos are then quasi-active or quasi-sterile, depending on the mixture of active and sterile neutrinos governed by a small mixing angle   ≃ √ 2  /  , with  the electroweak scale and   the mass of the known neutrinos.To leading order in   , quasi-active neutrinos are produced from quasi-sterile ones subsequent to the decay of .Consequently, the coupling   controls the dominant decay channels and allows for trading a factor (  / P ) 2 (with  P the Planck mass) for a (    /) 2 one in the decay width of .This trading enables the reduction of the width by a factor ∼ 10 −25  2  for a benchmark value   = 10 9 GeV, leading to the required lifetimes well beyond the age of the universe.
Ultra-high energy (UHE) neutrinos and photons are expected to emerge from the cascade of the  decay.In this Letter, we use the sensitivity of the Pierre Auger Observatory to such neutrinos and photons to impose constraints on the active-sterile neutrino mixing angle   for DM masses   ≳ 10 8 GeV.After recalling the main features of the Observatory that enable the detection of neutrinos and photons, we present the main decay channels of the  particle.We then calculate the number of neutrinos and photons expected to be observed at the Observatory as a function of   and   .Their non-observation allows us to constrain   for   ≳ 10 8 GeV.In a last section, we comment on the complementary bounds obtained on   from DM abundance and on   from cosmological observations.The Pierre Auger Observatory is the largest ground-based observatory that exploits extensive air showers to study UHE cosmic rays [4].Several detection techniques are combined.A surface of 3,000 km 2 is covered with a ground array of particle detectors separated by 1,500 m surrounded by fluorescence detectors spread over four sites.This baseline configuration is complemented with low-energy enhancements: a first nested array of 24 km 2 with particle detectors separated by 750 m and overlooked by three additional fluorescence telescopes with an elevated field of view, and a second one of 1.14 km 2 with detectors separated by 433 m, some of them on top of buried scintillation modules used to measure the number of highenergy muons.The fluorescence detectors provide, during moonless nights, direct observation of the longitudinal shower profile, which allows for the measurement of the energy and the primary-mass sensitive depth of the shower maximum,  max .On the other hand, the ground-level and underground detectors sample the shower particles with a permanent duty cycle.Although showers are observed at a fixed slice in depth with this technique, their longitudinal development is embedded in the signals detected.
Neutrinos of all flavors can be distinguished from nuclei through showers developing deeply in the atmosphere at large zenith angles (down-going detection mode) [16], while tau neutrinos can undergo charged-current interactions and produce a  lepton in the crust of the Earth that eventually decays in the atmosphere, inducing an upward-going shower (Earthskimming detection mode) [17].In both cases, neutrinos can be identified at the Observatory with signals in the groundlevel particle detectors that are spread over time, by contrast to narrow ones expected from the much more abundant nucleiinduced showers at large zenith angles [18][19][20].Photoninduced showers have also distinct salient features [21,22].Their first interactions are of electromagnetic nature, and the transfer of energy to the hadron/muon channel is reduced with respect to the bulk of hadron-induced showers.This results in a lower number of secondary muons.Additionally, as the development of photon showers is delayed by the typically small multiplicity of electromagnetic interactions, their  max is deeper in the atmosphere than for showers initiated by hadrons.
Based on these characteristics, a number of analyses have been designed to make the most of the sensitivity of the different detectors of the Observatory to neutrinos above ≃ 10 8 GeV [23] and to photons above 5 × 10 7 GeV [24][25][26][27].The non-observation of point sources and diffuse fluxes allowed the derivation of upper bounds that constrain various models very effectively.In the following, we draw constraints on the BSM benchmark [3] from these non-observations.
In the BSM benchmark [3], the DM candidate is a pseudoscalar particle  coupled to the sterile neutrinos sector alone.The interaction can be represented diagrammatically as with  the four-momentum of , and  the Dirac matrices.
Sterile neutrinos, which are right-handed, are denoted as .
Two types of Majorana sterile neutrinos are actually introduced.The first type,   , corresponds to the right-handed neutrino of the seesaw mechanism, with a superheavy mass   .The second type,   , with mass   , is necessary to allow the  particle to decay in an electroweak channel.In the basis of interaction eigenstates, the couplings between sterile neutrinos, active neutrinos and SM Higgs isospinor scalar fields through Yukawa parameters   and   give rise, after symmetry breakdown, to Dirac masses Assuming the hierarchy   <  D   ≪   , the three mass eigenstates associated to the eigenvalues ( 1 ,  2 ,  3 ) are well approximated by the following mixing between  L (active neutrinos),   and   (and their conjugate partners with superscripts c): with the mixing angle   , assumed to be ≪ 1, defined as For small mixings, the first (second) mass eigenstate  1 ( 2 ) is almost the light sterile (active) neutrino with a small admixture of the active (sterile) one, and that the third mass eigenstate  3 is almost the superheavy sterile neutrino.Note that to avoid complications that would be irrelevant to the question of DM, the flavor couplings in the active neutrino sector are considered diagonal in this study.The three mass eigenvalues read as The third eigenstate  3 appears to decouple from the two other states.While the superheavy particle   is essential to provide the mass to the active neutrino through a mixing angle   =   /2  identical to that of the standard seesaw mechanism [1,2], we shall ignore  3 hereafter.From the constraints on SM neutrino masses   ≤ 0.12 eV inferred from cosmological observations [28], we use hereafter  2 ≃ 0.04 eV as a benchmark.As for  1 , following [3], we use  1 = 10 −4 eV to fix the ideas; we shall see below that as long as  1 ≪  2 , results are not sensitive to the specific choice.That  1 ≪  2 is required for the lifetime of  to be much larger than the age of the universe.
To leading order in   , the interaction described by Eqn. 1 gives rise, in the basis of mass eigenstates, to the two-body decay  →  1  1 .However, in the relevant parameter space such that  1 ≪  2 ≪   ≪   ≪   , the total width Γ  is dominated by three-body channels stemming from the diagram depicted in Eqn. 1 and the interaction between active/sterile neutrinos and the Higgs isodoublet with Yukawa coupling gives the most important contribution to the width [3]: Due to the structure of the interaction depicted in Eqn. 1, a factor (  / P ) 2 expected from dimensional arguments [12] is traded for a factor ( 2 /) 2 .Although subdominant, there are two other three-body decay channels of interest, namely  →   1  2 and  →   1  with respective widths with  the weak-isospin gauge constant and  W the Weinberg angle.These two widths share the same structure as Γ  ℎ 1  2 modulo smaller numerical factors.
UHE neutrinos and photons are expected as byproducts of the decay of any particle with a mass much larger than the electroweak scale.Ultimately, they result from splitting-particle effects due to soft or collinear (real) radiative corrections enhanced by large logarithmic factors at high scale.The probability that a particle  at a scale   fragments to produce a particle  at a scale   carrying a fraction  of the initial energy is described by a fragmentation function (FF)    (;   ,   ).The FFs are evolved starting from measurements at the electroweak scale up to the energy scale   .In the QCD sector, the evolution is governed by DGLAP equations to account for the splitting function that describes the emission of parton  by parton .The resulting prompt spectra of photons and neutrinos have been derived in [29] and in several subsequent studies.Similarly, electroweak cascading can be described by evolution equations valid for a spontaneously broken theory [30].Seminal works of [29,31,32] in the QCD sector and of [33] in the electroweak one have provided the calculation of the FFs to derive the prompt flux of high-energy secondaries from the decay of a particle at high scale.We use hereafter the up-to-date HDMSpectra tool [34] to calculate the energy spectra of UHE neutrinos and photons from the decay of the  particle in the  → ℎ 1  2 ,  →   1  2 and  →   1  channels.
The number of neutrinos   () expected to be observed above an energy threshold  results from the integration over the sky of the directional exposure E  (, n) of the Observatory and of the flux of neutrinos emitted isotropically in proportion to the DM density accumulated in galaxy halos [35]: The first contribution is from the Milky Way halo, in which the DM energy density is parameterized by a profile function (x).x ⊙ is the position of the Solar system in the Galaxy,  is the distance from x ⊙ to the emission point, and n ≡ n(ℓ, ) is a unit vector on the sphere pointing to the longitude ℓ and latitude , in Galactic coordinates.There are uncertainties in the determination of the galactic-halo profile.We use here the traditional NFW profile as a reference [36], where  is the distance to the Galactic center,   = 24 kpc, and   is fixed by the DM density in the solar neighborhood, namely  ⊙ = 0.44 GeV cm −3 [37].To quantify the systematics stemming from the uncertainties on this profile, we repeat the analysis using other ones [38][39][40].The second contribution in Eqn. 12, which is obtained by integration over redshift  and which amounts to about 10% of the first one, is from all other galaxies. c is the critical energy density, Ω  is the DM abundance,   (, ) is the neutrino opacity of the universe as calculated in [41], and the Hubble rate  () depends on that observed today,  0 , and on the total matter abundance, In both contributions, the particle  decays into a (SM) daughter particle  whose FF leads to neutrinos.Thus, the differential decay width into neutrinos, dΓ  /d, results from the convolution of the differential decay width of  to  with the FF of  into neutrinos.For a single flavor, it reads as with  = 2  /  .In the ℎ 1  2 channel, the differential decay width reads as for neutrino final states, while it reads as for Higgs final states [3].Similar expressions hold in the cases of the   1  2 and   1  channels, with the corresponding decay widths.Finally, the differential widths entering into Eqn.14 account for all detectable flavors for the  2 species, where an explicit flavor index  is re-introduced.In the down-going detection mode, the three flavors contribute explicitly as , (17) with flavor indices  and .Additionally, the non-zero probability for neutrinos in the final state for "no-splitting" leads to an extra contribution with  = 0.13.The contribution to the Earth-skimming detection mode is obtained similarly, considering only the  flavor.Besides, the differential decay width into photons is also calculated following the same procedure, with proper FFs.Note that the expected photon fluxes to be compared to the flux limits are almost-entirely determined by the contribution of the Milky Way halo (due to their attenuation over inter-galactic scales).The resulting energy spectra are shown in Fig. 1 for   = 10 10 GeV.The high-energy enhancements are shaped by the non-zero probabilities for splitting a few times only at high scales.
Constraints in the planes (  ,   ) and (  ,   ) can be derived from the non-observation of neutrinos at the Observatory.First, 90% C.L. lower limits on the lifetime   are obtained by setting, for a specific value of   ,   () or   () to the 90% C.L. upper-limit numbers corresponding to the number of background-event candidates in the absence of signal [42,43].More details about the upper-limit numbers for each specific analysis searching for neutrinos or photons are provided in the supplemental material.Subsequently, a scan in   is carried out.It leads to a curve in the plane (  ,   ) that pertains to the energy threshold  considered.By repeating the procedure for several thresholds, a set of curves is obtained, reflecting the sensitivity of a specific energy threshold to some range of mass   .The union of the excluded regions finally provides the constraints in the (  ,   ) plane.Results are shown in Fig. 2 (top panel); lifetimes within the cross-hatched region are excluded.The region in full red pertains to a particular value of a Yukawa coupling    = 10 −5 , the meaning of which will be explained below.To illustrate the contribution from each secondary at our disposal, we show as the dotted the contribution to the constraints stemming from neutrinos alone; an analysis of the IceCube exposure dedicated to the benchmark-scenario decay channels would likely provide better sensitivity for exploring masses   ≲ 10 8.5 GeV.The lower limit on   is then transformed into an upper limit on   using the expressions of the total width of the particle .Results are shown in the bottom panel for separate values of   : color-coded regions pertain to    = 10 −5 while their extension (in cross-hatched) would require smaller values of    .Systematic uncertainties on   constraints amount overall to ≃ ±15%; they are dominated by those on the neutrino exposure [23].The restricted ranges of   for different   values come from the requirements not to overclose the universe with DM, while the exclusion hatched band comes from not altering the expansion history of the universe with the presence of ultra-light species such as sterile neutrinos   .We now briefly explain how these constraints are obtained.
In addition to its couplings to the DM sector and to the SM one through the Higgs isodoublet, the sterile neutrino   is also coupled to an inflationary sector in the BSM benchmark [3].This coupling, governed by a unique Yukawa parameter    for every  1 neutrinos, yields to a "radiative" production of  via a diagram similar to that depicted in Eqn. 8 (substituting  by the inflaton Φ in the initial state, and ℎ and  2 by  and  1 in the final states).Such a mechanism leads to a direct production of DM during the reheating period that can be sufficient, in general, to match the right amount of DM observed today [44].In the BSM benchmark [3], values for    are then required to range preferentially around 10 −5 .To infer the DM density   produced mainly during the reheating epoch, we also consider the minimal setup of gravitational production of  particles through the annihilation of SM (inflaton) particles as in [45] (as in [46]).In these conditions,  particles can be produced as long as the collision rate of particles is larger than the expansion rate  and/or as long as the inflaton field oscillates.By contrast,   is prohibitively low to allow any thermal equilibrium for DM.The collision term in the Boltzmann equation is then approximated as a source term only.Overall, the Boltzmann equation reads as Here, the sum on the right hand side represents the contributions from the SM and inflationary sectors.Using, on the one hand, the evolution of the SM-matter and inflaton densities derived in [47] and [46] respectively, and, on the other hand, the SM+SM →  +  and Φ+Φ →  +  reaction rates   derived in [48] and [46] respectively, the present-day relic abundance of DM can then be related, using Eqn.19 in the same way as in [15], to the mass   , the Hubble rate at the end of inflation  inf , and the reheating efficiency  quantifying the duration of the reheating period ( = 1 for an instantaneous reheating) [45].As a result, viable couples of values for ( inf ,   ) scale as  inf ∝  2  up to a maximum value for   , which depend on  and   .This scaling is a consequence of the domination of the radiative-production process over the gravitational one as long as the allowed values of  inf are too small for a given   value to generate significant particle production by gravitational interactions.For larger masses, the contribution from gravitational interactions added to the radiative production of  leads to an overproduction of DM that overcloses the universe, and there is thus no longer solution.This explains why the color-coded regions extend up to some maximum values of   in Fig. 2, for a benchmark value of    = 10 −5 .To the right of the regions shown in cross-hatched,    would need to be smaller.
Another constraint on   stems from the upper bound on the departure from 3 of the effective number of neutrino degrees of freedom  eff , Δ eff =  eff − 3, inferred from cosmological observations.From Friedmann equation, the Hubble parameter is governed by the energy-density content, which, during the radiation-dominated epoch, receives contributions from all relativistic species such as photons, active neutrinos and possible sterile ones.A significant value of Δ eff would thus change the dynamics of the expansion during the radiation era and thus impact on the cosmological microwave background measurements as well as on the big bang nucleosynthesis.Therefore, the 95%-confidence-level current limits on Δ eff < 0.3 [49] provide constraints on the temperature   1 of  1 neutrinos: This relationship is valid at any time before the radiation era started, in particular at the time each neutrino species decoupled from the thermal bath when the expansion rate became larger than the interaction rate    of the species   .Given that the ratio /   scales as (/ ★  ) 3 [50] (denoting as  ★  the decoupling temperature of   ), and given that the neutral-current-interaction rate of the  1 species is that of  2 but reduced by  2  , the relation  ★ 1 =  − 2 /3   ★ 2 holds.As soon as  dropped below  ★ 1 , the temperature of the freely expanding  1 species fell as   1 = ( ★1 /) ★ 1 (with  ★1 the scale factor at the time of the decoupling of the  1 species).At temperatures  <  ★ 1 , other particles that are still in thermal equilibrium, such as electrons and positrons, annihilate once  becomes smaller than their mass and thus heat the rest of the relativistic species (photons) relative to the neutrinos.Entropy conservation in this "freeze-out" then gives rise to the relationship with N ★ the number of degrees of freedom at  =  ★ .Combining Eqn.20 and Eqn.21 with N ★2 = 43/4 (accounting for two photon polarization states, three species of neutrinos and three of anti-neutrinos, and electrons and positrons each with two spin states) and  ★ 2 = 2 MeV [50], we obtain that N ★1 ≳ 35.95 and thus, using [51] to convert N into a temperature, the following bound on   , This bound is shown as the cross-hatched band in Fig. 2.1 In conclusion, we have shown that the data of the Pierre Auger Observatory provide stringent constraints on the angle   mixing sub-eV sterile and active neutrinos in the context of an extension to the SM that couples the sterile neutrinos to a superheavy DM candidate [3].For a typical dark coupling constant of 0.1, the mixing angle   must satisfy, roughly,   ≲ 1.5 × 10 −6 (  /10 9 GeV) −2 for a mass   of the darkmatter particle between 10 8 and 10 11 GeV.Future sensitivity to the effective number of neutrino degrees of freedom Δ eff from cosmological observations will be complementary, as they will probe   below 10 −3 independent of the mass   .In particular, a value of Δ eff departing significantly from 3 would call for sterile neutrinos and would fix the range of the mixing angle   .Our constraints would then be decisive in determining whether the new physics thus revealed in the neutrino sector is intimately related to that of superheavy DM.
culation is more involved.The efficiency  ES depends on the energy of the emerging  leptons   , the position of the signal pattern on the ground, the zenith and azimuth angles, and the altitude ℎ dec of the decay point of the  above ground. ES is averaged over the decay channels of the  with their corresponding branching ratios and is subsequently folded with the projected area of the detector, with the differential probability  exit (, ,   ) = d exit /d  of a  emerging from the Earth with energy   , as well as with the differential probability  dec (  , , ℎ dec ) = d dec /dℎ dec that the  decays at an altitude ℎ dec .The directional exposure is then The sky map shown in Fig. 3 (bottom panel) reflects that in local coordinates, the integrand of Eqn.24 is non-zero in a narrow zenithal range of a few degrees only just below the horizon.In the case of photons, the directional exposure is obtained through the time and area integration of the detection efficiency   and selection efficiency   projected onto the direction perpendicular to the arrival direction of the photon, Four different analyses, differing in the detector used, have been developed to cover the wide energy range probed at the Observatory.In particular, two of these analyses benefit from a direct estimate of the depth of shower maximum by the fluorescence detector as one of the discriminating variables.The use of the fluorescence-detector datasets introduces for these analyses an explicit dependence in time for   due to the limited duty cycle on moonless nights that propagates into a small dependence in right ascension  for E  .In addition to the detection efficiency, the   factor accounts for the dependencies of the selection process aimed at separating hadronicinduced showers from photon-induced ones.By contrast to the neutrino case, searches for photons are not backgroundfree and, for some of the energy thresholds explored, there is an irreducible small number of charged cosmic-ray events that resemble photon events despite the selection process.The values used for   () result from the upper limits on the number of candidates accounting for this irreducible background; they are listed in Table I.For illustration, we show in Fig. 4 an example of directional exposure to photons at 10 17 eV (top panel) and at 10 19 eV (bottom panel).The former is built on the 1.14 km 2 array with a separation of detectors of 433 m optimized to study the range of energies around 10 17 eV, while the latter on the 3,000 km 2 one, optimized for higher energies.This explains the large increase of exposure observed at high energies.

Figure 1 .
Figure1.Energy spectrum of neutrinos and photons from the decay of the pseudo-scalar particle  within the BSM benchmark[3] (  = 10 10 GeV).

Figure 2 .
Figure 2. Top: Constraints on   as a function of   for a value of the couplings of sterile neutrinos with the inflationary sector    = 10 −5 .The dotted line illustrates the constraints stemming from neutrino secondaries alone.Bottom: Constraints on   as a function of   for three different values of the coupling constant   , still for    = 10 −5 .The hatched-red region   ≥ 9×10 −4 is excluded from the constraint on Δ eff (see text).

Figure 4 .
Figure 4. Directional exposure of the Pierre Auger Observatory in Galactic coordinates to photons at 10 17 eV (top) and 10 19.5 eV (bottom).

Table I .
Upper limits on the number of photon candidates for each energy threshold.