Indirect Detection of Sub-GeV Dark Matter Coupling to Quarks

We consider photon signals arising from the annihilation or decay of low-mass (sub-GeV) dark matter which couples dominantly to quarks. In this scenario, the branching fractions to the various kinematically accessible hadronic final states can largely be determined from chiral perturbation theory. Several of these final states yield striking spectral features in the sub-GeV photon spectrum. New experiments, such as e-ASTROGAM and AMEGO, are in development to improve sensitivity in this energy range, and we discuss their potential sensitivity to this class of models.

Subject to these selection rules, the only mesons which we will be interested in are π 0 (m π 0 ∼ 135 MeV), π ± (m π ± ∼ 140 MeV) and η (m η ∼ 548 MeV). Since kaons must appear in pairs to conserve strangeness, they can only be produced if √ s is very close to 1 GeV. ρ ± , ρ 0 can also be produced in conjunction with pions, but only if √ s 0.91 GeV. Although there are some regions of parameter space where these particles can be relevant, they are kinematically inaccessible in most of the parameter space we consider, so we will ignore them.
The decay of π ± generally produce muons and neutrinos, with only a small contribution to the photon spectrum. The mesons most relevant to the photon spectrum are π 0 , which decays to γγ with a ∼ 99% branching fraction, and η, which decays to γγ with a ∼ 39% branching fraction. Note, η also decays to 3π 0 with a ∼ 33% branching fraction, and to π + π − π 0 with a ∼ 23% branching fraction. As a result, a single annihilation or decay process can produce a relatively large number of neutral pions. Although this results in a larger number of secondary photons, their energy spectrum is much more complicated. Moreover, the astrophysical gamma ray backgrounds tend to grow rapidly at lower energy; since the photons arising from π 0 decay tend to be much less energetic than those arising from η decay, they will compete against a much larger background. As such, we will focus only on the photons arising from η decay. If √ s < 1 GeV, then the only kinematically allowed final states including at least one η are π 0 η, π 0 π 0 η and π + π − η. Since the relevant mesons are all pseudo-Nambu Goldstone bosons (pNGBs) of chiral symmetry breaking, their interactions with dark matter can be described using chiral perturbation theory, in which the dark matter is treated as a spurion whose interactions break flavor symmetry. This spurion is set equal to the dark sector operator which couples to quark bilinears in the fundamental Lagrangian. For the case of dark matter decay, the spurion will be the dark matter field, whereas for the case of dark matter annihilation, the spurion will be the appropriate dark matter bilinear, weighted by an energy scale associated with the energy scale of dark sector interactions. In the chiral Lagrangian, this spurion then couples to the octet of pNGBs, and the form of this interaction is determined by the Lorentz, parity, and flavor transformation properties of the quark bilinears to which the spurion couples in the fundamental Lagrangian. In the Standard Model, one already introduces scalar, pseudoscalar, vector, and axial-vector spurions in order to describe quark masses and electroweak interactions. As a result, one can use data to determine the coefficients of the operators coupling these spurions to the meson octet, order by order in the chiral Lagrangian. We will thus only consider the cases in which dark matter couples to scalar, pseudoscalar, vector, or axial-vector quark currents, and the corresponding spurions will be denoted as s, p, v µ and a µ , respectively.
We will work to lowest order in the chiral Lagrangian; this will be a good approximation for our purposes, especially if the momenta of the final state particles is small. The effective Lagrangian is thus given by where The field η 8 appearing above is linear combination of the physical η and η , with η 8 = η cos θ P +η sin θ P , and θ P ∼ 11.5 • (cos θ P ∼ 0.98). Henceforth, for simplicity, we will simply equate η 8 with the physical η. As expected, the constants F and B are determined from data; F ∼ F π ∼ 92 MeV is the pion decay constant, and B = m 2 π /(m u + m d ) + O(m 2 ). Henceforth, we will take m π 0 ∼ m π ± = m π , as this approximation will only have a non-negligible effect very near the threshold for π ± production. The chirality-violating real-valued spurions are given by where, for simplicity, we have considered the case of dark matter annihilation. If dark matter instead decays, the dark matter bilinear would be replaced by a single dark matter field. The dimensionless coefficients are equal to the coefficients of the operator coupling the dark matter bilinear to the appropriate quark current in the fundamental Lagrangian.
We can now expand this Lagrangian, in order to determine the coefficients of the operators which couple dark matter to the light mesons. Keeping only the terms which couple the spurions to at most three mesons, we find We immediately see that, at this order in the chiral Lagrangian, there are contact interactions which couple s to two-body L = 0 states and couple v i to two-body L = 1 states. But the contact interactions couple p and a 0 to three-body L = 0 states, and couple a i to three-body L = 1 states. But p and a µ can also couple directly to π 0 and η, allowing dark matter to annihilate through a mediator in the s-channel. In the case of dark matter decay, we would instead find mixing between dark matter and either π 0 or η, and this mixing is can be constrained by data. But this coupling of dark matter to a single meson will vanish if the coefficients are flavor-universal. A few other features are apparent from the Lagrangian in eq. 4. The contact interactions between dark matter and the final states we consider are actually independent of the dark matter coupling to strange quarks. We can see this by noting that all of the terms in eq. 4 involve a single trace; thus, if a trace involves more than one insertion of Φ, then α s S,P,A,V can only appear along with Φ 3i and Φ j3 . These terms are only relevant if the final state has more than one η or kaon, and these states are not kinematically accessible. To study contact interactions, we may thus truncate Φ to the upper-left 2 × 2 block, which we denote asΦ = (η/ √ 6)σ 0 + (π 0 / √ 2)σ 3 + π + σ + + π − σ − . Note, however, that a dark matter coupling to strange quarks does permit dark matter to annihilate to a three-body final state via an intermediate η in the s-channel.
As each spurion is a diagonal 3 × 3 matrix, it will be convenient to parameterize the flavor structure of spurions by expressing them as a linear combination of M 1 = diag(1, 1, 1), M 2 = diag(1, −1, 0), M 3 = diag(−1, −1, 2). In particular, spurions proportional to M 1 and M 3 will have the same contact interactions, but the spurion proportional to M 3 will also couple to mesons through an intermediate η, while the spurion proportional to M 1 will have only contact interactions. A spurion proportional to M 2 will have contact interactions with mesons, as well as a coupling through an intermediate π 0 .
Φ is a linear combination of an isospin singlet (η) and an isospin triplet (π ± , π 0 ), and the isospin of the final state is determined by the isospin of the spurion associated with the initial state. In particular, if the spurion describing the initial state is proportional to M 1 or M 3 , then the final state has I = 0, whereas if the spurion is proportional to M 2 , then the final state has I = 1, I 3 = 0 (I is the isospin quantum number). It is straightforward to determine which final states can be produced by various choices of spurions using the transformation properties of various particles under C, P and I. Note that a ππ pair with I 3 = 0 transforms as C : (−1) L , P : (−1) L , where L is the angular momentum of the two-body state (for example, see [15]). Symmetry of the wavefunction under particle interchange I: For each choice of spurion (s, p, vµ, aµ) and flavor structure (M1, M2, M3), we list the J P C quantum numbers of the initial state, the possible final states, and the average number of primary η and π 0 particles produced per annihilation/decay event. For final states with more than one pion, we indicate the total isospin of the pion state. We include only spurions which couple to final states including an η.
then requires I = L mod 2. We summarize the connection between various spurions and the relevant final states in Table I. Note that, at lowest order in the chiral Lagrangian, the v i spurion couples only to the π + π − state, whose decays yield few secondary photons. (The two-body final state must have J P C = 1 −− , and must thus be a ππ-state with I = L = 1; such a state has no π 0 π 0 contribution.) Final states involving π 0 or η arise only at higher order in the momentum expansion, and will have a small branching fraction. As such, we will henceforth ignore the case in which dark matter couples to vector quark currents. Similarly, there are several spurions which only couple to ππ or πππ final states, and we ignore them as well.
We are left with only the spurions s M2 , p M1,M3 and a M1,M3

0
. The spurion s M2 can only produce the final state π 0 η, as this is the only two-body 0 ++ final state with I = 1, I 3 = 0. On the other hand the spurions p M1,M3 and a M1,M3 0 can each produce two three-body final states: π 0 π 0 η and π + π − η. But the branching fractions to these states (2/3 and 1/3, respectively) are determined by the fact that the two outgoing pions must be in an I = 0 state. But in any case, the branching fractions are not really relevant to our analysis, as we are only focussing on photons arising from the decay of the η.
Unlike the spurions (p, a 0 ) M1 , the the spurions (p, a 0 ) M3 can also couple to an intermediate η in the s-channel. But although this additional coupling will affect the total annihilation cross section (or decay width, in the case where the spurion is a single dark matter field), it will not affect the corresponding branching fractions, which are determined by the quantum numbers of the initial state. Thus, for our purposes, the spurions proportional to M 1 and M 3 can be treated identically. Note that, although we have obtained these results from considerations of symmetry, they can of course be obtained explicitly from the chiral Lagrangian.

III. PHOTON SPECTRA
Given a particular choice of Lorentz and flavor structure for the dark matter coupling to quarks, it is straightforward to determine the resulting secondary photon spectrum. Since the secondary photons arise from the process η → γγ, the resulting photon spectrum is simply related to energy spectrum of the parent η.

A. η Injection Spectra
If dark matter annihilation or decay, with center-of-mass energy √ s, produces a π 0 η final state, then the energy spectrum of the outgoing η is given by If dark matter annihilation or decay produces a ππη final state, then the energy spectrum of the η can be determined by integrating the squared matrix element over the three-body phase space, which can be expressed in terms of the energies of any two of the three final state particles. In particular, we have where Note, if the final state is π 0 π 0 η, then the phase space integral has an additional 1/2 combinatoric factor. We will take i = 1 to denote the η, and i = 2, 3 to denote the two pions (either π + π − or π 0 π 0 ). Then we find Note that the combinatoric factor cancels in this expression. In general, the energy spectrum of a meson which is part of a three-body final state depends on the energy dependence of the matrix element. But one can verify from the chiral Lagrangian that, for the spurions we consider here, the matrix element is independent of x i . We thus find The injection spectrum for η is a single bump, which vanishes at the kinematic endpoints x min corresponds to the limit in which the η is produced with no boost, and the two pions are produced back-to-back, while x max corresponds to the limit in which both pions move with the same boost, in the opposite direction to the η.

B. Secondary Photon Spectrum
Given the injection spectrum of the parent η, the secondary photon spectrum is given by [16] A few points about the photon spectrum are immediately apparent [16]. This spectrum is log-symmetric about the energy scale E * = m η /2, and has a global maximum at E γ = E * . Moreover, the photon spectrum is non-increasing as E γ either increases or decreases away from E * . If dark matter decay or annihilation produces a π 0 η, than the η is monoenergetic with an energy determined by √ s; we will denote the boost of the η by γ η = E η /m η , with β η = 1 − (1/γ η ) 2 . From eq. 10, we see that the photon spectrum is uniform between the limits (m η /2)[γ η (1 − β η )] and (m η /2)[γ η (1 + β η )], and vanishes outside these limits.
On the other hand, if a three-body final state is produced, then the photon spectrum will be peaked at (m η /2) and will fall off monotonically in either direction until it vanishes at the limits (m η /2)[γ η (1 ± β η )], where γ η is the maximum boost factor for the η which is kinematically allowed. Following the results in [16], we see that the fact that dN η /dx 1 vanishes as x 1 → 2m η / √ s implies that the maximum in the photon spectrum at E γ = m η /2 is a smooth peak. For illustration, we plot in Figure 1 the photon spectrum arising from the decay process η → γγ for a ηππ final state with √ s = 850 MeV.

IV. SENSITIVITY
We will consider the sensitivity to dark matter annnihilation/decay which can be provided by experiments which are sensitive to O(10 − 1000) MeV gamma rays. Following [9], we consider constraints arising from diffuse emission in Galactic halo, and from emission from one particular dwarf (Draco). One can also consider limits from the Galactic Center, and these have been considered in a related context in [12]. However, there is great systematic uncertainty regarding astrophysical foregrounds/backgrounds from the direction of the Galactic Center. As a result, we do not consider this target in this work.
Since we are considering experiments which are in still in the design phase, we will make no serious attempt to optimize our analysis strategy. Instead, we will simply choose some reasonable cuts which will give us a good estimate for the sensitivity which can be obtained. We will consider, as a benchmark, an experiment with an exposure of 3000 cm 2 yr and a fractional 1 − σ energy resolution of = 0.3 throughout the entire range of interest. We will also assume an angular resolution of 1 • throughout the energy range of interest. The smallest target we will consider here is Draco, with an angular size of 1.3 • , and we assume the angular resolution is smaller than this size.
The differential photon flux arising from dark matter annihilation or decay can be written as where andJ ann.,dec. is the average J-factor of the target for either annihilation or decay (we assume that the dark matter particle is its own anti-particle). We consider diffuse emission in the region |b| > 20 • , where b is the latitude in Galactic coordinates. In this region, we will take the averaged J-factors to be given by [17] J ann. dif. = 3.5 × 10 21 GeV 2 cm −5 sr −1 , J dec. dif. = 1.5 × 10 22 GeV cm −2 sr −1 .
For Draco, we will take the averaged J-factors to be [18] J ann. Draco = 6.94 × 10 21 GeV 2 cm −5 sr −1 , J dec. Draco = 5.77 × 10 21 GeV cm −2 sr −1 , with uncertainties estimated at ∼ 60% [18]. The isotropic flux observed by COMPTEL (0.8 − 30 MeV) and EGRET (30 MeV-1 GeV) [19] can be well fit [9] to the function We will assume that the energy spectrum actually reported by the experiment is given by where is a smearing function which accounts for the fractional energy resolution of the instrument. If we denote by I exp. the exposure of the instrument, and by ∆Ω the solid angle viewed, then the number of events due to dark matter annihilation or decay expected to be observed within the energy window E − ≤ E obs. ≤ E + is given by while the number of expected events actually observed (based on the fit to COMPTEL and EGRET) is given by A. Bounds from diffuse emission Sensitivity to diffuse emission from dark matter annihilation or decay is largely controlled by systematic uncertainties in the astrophysical background. We will adopt the following criterion for estimating the sensitivity of a given instrument to diffuse emission: a model can be excluded if, within any energy bin of size set by the energy resolution (E + − E − = (E + + E − )), N S > αN O , where α is a constant set by the systematic uncertainty of the background. For example, a conservative bound would arise from setting α ∼ 1, and would be appropriate if one had little confidence in any background model; in this case, a model could only be excluded if there was an energy bin in which the estimated number of dark matter events exceeds the entire number of observed events (including statistical uncertainty). If one were confident that the background were smooth, then α would instead be determined by small fluctuations which could be accommodated by the uncertainty in the fit to the observed flux. For relatively large spectral features, this uncertainty would yield α ∼ 0.15 [9,12], but would decrease to ∼ 0.02 for narrow features [12].
This analysis then yields the following constraint: ann.,dec. dif.
where E 0 is the center of the energy bin. For a two-body final state, each meson produces a secondary photon spectrum which is constant over some energy range. In this case, E 0 should be chosen so that upper edge of the energy bin (E 0 (1 + )) lies at the highest energy such that the secondary photon spectrum is non-vanishing. For a three-body final state, one should choose E 0 = m η /2. Note that the sensitivity of an experiment to diffuse emission is determined by the signal-to-background ratio, and is thus independent of the exposure and angular resolution. Near threshold, when the spectral features are sharp, sensitivity will scale as −1 . But as √ s increases and the spectral features become large compared to the bin size, the dependence on disappears. In particular, at this level, future experiments would give no improvement over current bounds from the EGRET flux measurement.
But it is important to note that the exposure and angular resolution of future experiments can indirectly affect sensitivity to diffuse emission. We have assumed an isotropic flux equal to that observed by COMPTEL and EGRET. But if a future experiment is able to resolve a significant number of point sources which can then be subtracted from the isotropic background, the remaining background may be significantly smaller. Since sensitivity to diffuse emission is determined by the signal-to-background ratio, any reduction in the normalization of the observed isotropic flux (assuming the same spectral shape), will improve sensitivity by the same factor.
Note also that we have not included any uncertainty in the J-factor. But since N S scales linearly withJ, any deviation of the actual J-factor from our estimate will simply rescale our bound by the same factor.

B. Bounds from dSphs
Sensitivity to emission from a dwarf spheroidal galaxy, on the other hand, is largely controlled by statistical uncertainties. In this case, the observed isotropic flux can be treated as a estimate for the background in a search for emission from the dwarf (a more accurate estimate can be made for any particular dSph by measuring the flux from the direction of the dSph, but slightly off-axis [20][21][22][23][24]). This background includes all emission from dark matter annihilation or decay outside the dwarf but along the line of sight, as well as emission from astrophysical processes. Given an estimate for the expected background, and a measurement of the number of photons seen from the direction of a dSph, one can determine if any particular model is consistent with the data to any particular statistical confidence level. We will adopt the following criterion for estimating the sensitivity of a given instrument to emission from a dSph: if the number of events observed from the dSph in some energy range is equal to the expected number of background events (N O ), then a model can be ruled at confidence level n − σ if N S > n √ N O , where N S is the expected number of signal events in the same energy range.
We can express this constraint as Ξ ann.,dec. s −1J ann.,dec. Draco where E 0 is again the center of the energy bin. Here, we see that sensitivity scales with (I exp. ∆Ω) 1/2 .

C. Results
In Figure 2, we present lower bounds on the lifetime of decaying dark matter, for the case of either a π 0 η (blue) final state arising from a scalar interaction, or a ππη (red) final state, arising from a spurion which is either pseudoscalar, or the timelike component of an axial vector. Solid lines denote conservative bounds on (α = 1) diffuse photon emission, and dashed lines denote 2σ-bounds on photon emission from Draco. For the case in which dark matter annihilates, we present similar upper bounds on the annihilation cross section in Figure 3. In this figure, we also present recent bounds from Planck [25] (dotted black) on f ef f. σ A v /m X , where we have chosen f ef f. ∼ 0.4 as a rough approximation over the energy range and final states of interest [26]. Note, Planck bounds on decaying dark matter only provide a lower bound on the lifetime of ∼ O(10 24 s), which does not appear in Figure 2. For simplicity, in both figures we plot on the horizontal axis the center-of-mass energy ( √ s) of the process. This is equal to the dark matter mass in Figure 2, and equal to twice the dark matter mass in Figure 3. Since we have only included the photons arising directly from η decay, we will have missed a small number of photons arising from π 0 decay which are energetic enough to overlap the energy bins we consider; our bounds are thus conservative.
Although we have chosen specifications for an instrument which would roughly match the design of e-ASTROGAM, we have noted that the sensitivity to diffuse emission is largely independent of the exposure and angular resolution. As such, these conservative bounds on diffuse emission represent actual bounds which one can place on dark matter models using data from EGRET. Moreover, the sensitivity scales directly with the diffuse J-factor, and inversely with the magnitude of the background. If there are new estimates of the J-factor, or if point source measurements lead to a reduction of the diffuse background flux, then these limits can be rescaled appropriately.
The sensitivity from dSphs, on the other hand, scale directly with the J-factor, but only with the background to the −1/2 power. Moreover, the sensitivity from dSphs scales as the square root of the exposure. Unfortunately there is no dSph analysis available from EGRET, but for any future experiment with any exposure, the estimated sensitivity can be derived from these results by an appropriate rescaling. In fact, the estimated sensitivity for the benchmark case we consider is slightly worse than current bounds from Planck. However, if the J-factor for Draco is a little larger than the estimate we have used, then indirect detection could be competitive with current Planck bounds. Moreover, a stacked analysis of many dSphs would surely increase the sensitivity of indirect detection methods.

FIG. 2:
We plot lower bounds on the dark matter lifetime for the case in which the final state is π 0 η (blue) or ππη (red). We plot conservative limits on diffuse emission (solid) and 2σ-limits on emission from Draco (dashed).

FIG. 3:
We plot upper bounds on the dark matter annihilation cross section for the case in which the final state is π 0 η (blue) or ππη (red). We plot conservative limits on diffuse emission (solid) and 2σ-limits on emission from Draco (dashed). We also plot bounds arising from Planck data [25] (black dotted).
We note here that, in the case where dark matter annihilates through a scalar interaction to a π 0 η final state, the annihilation cross section is p-wave suppressed. In this case, the astrophysical factor associated with dark matter annihilation is not the standard J-factor associated with velocity-independent dark matter annihilation, but instead depends on the dark matter velocity-distribution [27][28][29][30]. A determination of the effective J-factor for Draco in the case of p-wave suppressed dark matter annihilation is beyond the scope of this work; instead, the associated velocity-suppression is absorbed into the quantity σ A v , which implicitly includes a convolution of the annihilation cross-section with dark matter velocity distributions in Draco. In this case, the bound from Planck cannot be directly applied, since the typically velocity of dark matter particles in the early Universe was far different that in the current epoch.
In the case of either dark matter or dark matter annihilation, the sensitivity arising from a search for emission from Draco exceeds the sensitivity of a search for diffuse emission. But this improvement is more pronounced for the case of dark matter annihilation, as in this case the photon signal scales as the square of the dark matter density in the dwarf. Similarly, bounds from Planck are the most constraining for the case of dark matter annihilation; this improvement is driven by the fact that the annihilation rate scales as the square of the density, which was much larger at the time of recombination than in the present epoch.
By comparison, we can consider the prospects for a search for dark matter decay or annihilation in the Galactic Center, which we define by the range |b| < 5 • , | | < 30 • . In this region of the sky, the observed flux is about a factor of 20 larger than the observed diffuse flux at O(100 − 1000 MeV) [12]. The average J-factor for decay for the Galactic Center is a factor of 3-4 larger than for the Galactic halo at high latitudes [17]. As a result, the sensitivity of an instrument to dark matter decay in the GC would be worse than from diffuse emission due to decay throughout the halo. On the other hand, the average J-factor for annihilation for the GC is about a factor of 30 larger than for the Galactic halo at high latitudes [17], assuming an NFW profile. Thus, there is little to be gained in searching the GC for dark matter annihilation or decay, in a conservative analysis. Note, however, that in a related context, the authors of [12] found much better prospects for studying emission from sub-GeV dark matter annihilation in the GC. We believe that this discrepancy is largely due to the fact that we have assumed a conservative analysis, in which a model can only be excluded if the expected number of signal events exceeds the observed number of events. More optimistic assumptions about one's ability to understand the backgrounds from the GC were made in [12].

D. Constraints from the LHC and Direct Detection
If dark matter couples to light quarks, then one might hope to constrain these scenarios using LHC mono-anything searches for processes such as pp → XX + jet [31,32]. For example, if a dark matter bilinear coupling to quarks is represented by a p 1 spurion, then one would need α 1 P /Λ 2 ∼ (O(100) GeV) −2 in order for the annihilation cross section to be O(10 −2 ) pb, which is the approximate limit obtained from Planck data. If the contact approximation were valid at the energies of the LHC, then such a coupling would already be ruled out by LHC searches. But this result is highly model-dependent, and for sub-GeV dark matter, there is no real reason to expect the contact approximation to be valid at O(TeV) energies. By reducing the scale of the couplings α and the energy scale Λ, one can essentially evade current bounds from LHC mono-anything searches.
Direct detection experiments, such as CRESST [33], are now probing the mass range we consider here. If a dark matter bilienar interacts with quarks through a pseudoscalar interaction, then the dark matter-nucleon scattering cross section is v 4 -suppressed and spin-dependent, and CRESST would be unlikely to see a signal. If a dark matter blinear instead couples through a scalar interaction, then dark matter can also have velocity-independent spin-independent scattering with nuclei. But since dark matter annihilation through a scalar interaction is p-wave suppressed, one needs α 2 S /Λ 2 ∼ (O(1 − 10) GeV) −2 in order for the annihilation cross section to be O(10 −2 ) pb. For such models, the dark matter-proton scattering cross section would be much larger than a picobarn. However, for the scalar interaction generated by spurion s 2 , dark matter interactions are maximally isospin-violating [34][35][36], which would suppress the naive sensitivity of CRESST to these models. Such models could also be probed by the effect of scattering on CMB [37,38], but current bounds are not constraining. A detailed study of the direct detection prospects for these models for current and upcoming direct detection experiments would be very interesting, but is beyond the scope of this work.

V. CONCLUSIONS
We have considered the indirect detection of sub-GeV dark matter annihilation or decay. If dark matter couples to quarks, then the hadronic final states and branching fractions are largely determined by symmetry and kinematics, and can be derived in chiral perturbation theory. In particular, striking photon signals can be produced by the process η → γγ. Especially for the case of dark matter decay, the current lower bounds which can be obtained from EGRET data already exceed those obtained from Planck by orders of magnitude. Future data from an experiment such as e-ASSTROGAM, looking at dwarf spheroidal galaxies, can provide an even greater improvement.
In this work, we have utilized chiral perturbation theory at lowest order, and have focussed on the photons arising directly from the decay process η → γγ, since the signal is well understood and the background is small. But dark matter annihilation or decay in this energy range generally produces a larger number of pions, especially after including η decay, but the kinematics are more difficult. More generally, for slightly larger energies, a much wider range of final states are accessible and become relevant for indirect detection. But as increasing interest is shown in sub-GeV dark matter, it would be interesting to perform a more comprehensive study of the hadronic final states which can be produced.
In a similar vein, it is worth noting that if sub-GeV dark matter couples primarily to light quarks, then it can potentially be produced at proton beam fixed-target or beam-dump experiments, such as NA62 [39] or SeaQuest [40], or proposed experiments such as SHiP [41], or related proposed experiments such as FASER [42]. In particular, one might hope that dark matter could be produced in the rare decays of heavier mesons. But to determine the available signals at such experiments, a more detailed study beyond lowest order in chiral perturbation theory would be necessary.