Hazma: A Python Toolkit for Studying Indirect Detection of Sub-GeV Dark Matter
Adam Coogan Logan Morrison Stefano Profumo
hep-ph astro-ph.HE physics.comp-ph

G-flocking: Flocking Model Optimization based on Genetic Framework
    Li Ma, Weidong Bao, Xiaomin Zhu, Member, IEEE, Meng Wu, Yuan Wang, Yunxiang Ling, and Wen Zhou
This work was supported in part by the National Natural Science Foundation of China under Grants 61872378,  91648204, 71702186, in part by Postgraduate Research Innovation Project in Hunan Province under grant CX2018B021, in part by the Scientific Research Project of National University of Defense Technology under Grants  ZK17-03-48
    Received: date / Accepted: date
=================================================================================================================================================================================================================================================================================================================================================================================================================================================





§ INTRODUCTION



The search for particle debris from dark matter (DM) annihilation or decay has thus far largely centered on DM masses in the GeV-TeV scale, for a variety of reasons. First, if the DM shares electroweak interactions with the Standard Model, as in the weakly-interacting massive particle (WIMP) scenario, then, by the Lee-Weinberg limit, its mass is expected to be more than a few GeV<cit.>[Note that exceptions exist to the Lee-Weinberg limit, see e.g. <cit.>]. Second, the general expectation for the scale of new physics based on the “small hierarchy problem” is that new physics, and thus new massive particles possibly including the particle making up the cosmological DM, should appear around the electroweak scale. Finally, the GeV scale is testable with an array of currently-operating gamma-ray and cosmic-ray observatories, including, but not limited to, the Fermi Large Area Telescope (LAT) <cit.>, Cherenkov Telescope Arrays such as MAGIC <cit.>, HESS <cit.>, VERITAS <cit.>, and the Alpha-Magnetic Spectrometer (AMS-02) <cit.>.

On the theory front, the calculation of the detailed expected particle spectrum of the debris resulting from DM annihilation or decay has thus focused on the GeV-TeV regime. State-of-the-art codes utilize simulations describing the results of high-energy collisions of elementary particles yielding jets and leptons, which in turn decay and produce stable final-state particles. Many such codes, such as DarkSUSY <cit.>, micrOMEGAs <cit.> and PPPC4DM <cit.> utilize tabulated results from PYTHIA <cit.>, one of the most widely-used programs for performing these simulations. Such results are reliable at center-of-mass energy scales at or above roughly 5 GeV <cit.>, but not at lower energies, where, for instance, strongly-interacting particles form hadronic bound states and are no longer described by parton showers, fragmentation and decay. It is well known that the resulting spectra of gamma rays, electrons and positrons and antiprotons, are dramatically different in that case.

For a variety of reasons it is now quite timely to offer the community a reliable computational package that provides the spectra of particles resulting from lighter, sub-GeV DM annihilation or decay. First and foremost, MeV astronomy will soon be revolutionized with a new generation of telescopes such e-ASTROGAM <cit.>, AMEGO <cit.> and others, including concept telescopes such as the Advanced Energetic Pair Telescope (AdEPT) <cit.>, the PAir-productioN Gamma-ray Unit (PANGU) <cit.>, and the Gamma-Ray Imaging, Polarimetry and Spectroscopy (“GRIPS”) <cit.>. Second, the persistent absence of any conclusive astrophysical signal from DM in the GeV-TeV range has furthered theoretical and phenomenological interest in the mass range below the GeV, providing additional motivation to investigate the details of DM decay or annihilation processes. Lastly, at present no code exists that allows users to readily study gamma-ray and cosmic-ray production from DM particles annihilating dominantly into hadronic bound states. 

With these motivations, we here introduce a Python toolkit, hazma[Hazma is a small rounded Pokemon, with light green spikes running down its back and tail, making it appear somewhat dinosaurian. Its body resembles a yellow hazardous materials suit, with a face resembling a respirator or gas mask, and a zipper-like marking running down its stomach. It has two stubby legs, the feet of which are green. Hazma is one of the few stable Nuclear types, the others being Nucleon and Urayne. Its leaded skin makes it immune to nuclear radiation. In the aftermath of a nuclear accident, groups of Hazma will appear and feed on the radioactive gas, eventually cleaning the air of the area over time. <cit.>], that computes spectra of gamma rays and cosmic-ray electrons and positrons from the decay of muons and pions, calculates constraints from gamma-ray observations and the cosmic microwave background, and allows users to compute composite spectra for selected built-in models of DM-parton interactions.

From a field theoretic standpoint, the description of the interactions of fundamental fields with hadrons is performed in the context of chiral perturbation theory (ChPT, see e.g. Ref. <cit.> for a review). A full account of mapping a fundamental, parton-level Lagrangian onto its ChPT counterpart will be given elsewhere <cit.>; here, however, we do provide selected examples of how models where the DM interacts with a mediator of specified spin and parity produces ChPT vertices.

As for any effective theory, ChPT possesses a certain range of validity which depends upon the size of some dimensionless parameter, here the ratio of the meson momentum to a scale Λ_ ChPT≡ 4π f_π∼ 1 GeV. Below we will describe the range of dark matter and mediator masses for which our EFT framework can be reliably used to compute annihilation cross sections and mediator decay rates. The mass ranges dictate which combination of mesons we include in the computational package we hereby present.

In the light DM mass limit, we also found that the standard approach for studying radiative emission from leptonic final states is problematic. In short, utilizing the Altarelli-Parisi splitting function to calculate the final state radiation spectrum assumes that radiating particle's center-of-mass frame energy is much larger than its mass. For 𝒪(100MeV) dark matter annihilating into muons, this is not the case. As a result we compute the exact spectrum for a few model cases. We also provide spectra for the final state radiation off of charged pions, and account for radiative decays of all relevant particles (e.g. π^+→μ^+νγ, π^+→ e^+νγ etc).

A general issue with light, sub-GeV DM models is that constraints from perturbations to the cosmic microwave background (CMB) are generically very strong if the DM freezes out as a thermal relic. Of course there exist a broad variety of workarounds and caveats (see e.g. <cit.>), but any light DM model is prone to CMB constraints. hazma implements such constraints by including functionality for computing electron and positron spectra from dark matter annihilation.

The remainder of this paper is structured as follows. Sec. (<ref>) offers a high-level overview of the hazma code and introduces the Theory class, the main user-facing component of hazma. Section (<ref>) describes the effective field theory framework used to study sub-GeV dark matter, provides details of the scalar mediator (sec. (<ref>) and vector mediator (sec. (<ref>)) models, and describes the particle physics outputs of hazma, including cross section, decay widths, and branching fractions. Section (<ref>) explains how to calculate gamma-ray spectra from individual annihilation final states, while Sec. (<ref>) combines the latter with the scalar- and vector-mediator models to obtain the overall gamma-ray spectra from DM annihilation. Section (<ref>) describes the calculation of the positron spectra, and, finally, Sec. (<ref>) and Sec. (<ref>) describe, respectively, gamma-ray and CMB limits. Section (<ref>) concludes. Examples of how to use hazma are woven throughout the text. Appendix (<ref>) describes the installation process for hazma, App. (<ref>) review the basics of using hazma, and App. (<ref>) gives examples of more advanced applications of the code, such as incorporating new models. Appendices (<ref>-<ref>) describe modules provided with hazma for convenience which can perform arbitrary phase-space integrations and compute photon spectrum for ≥ 3-body final states.

Our code is available on GitHub at https://github.com/LoganAMorrison/Hazmahttps://github.com/LoganAMorrison/Hazma, and the manual is located at https://hazma.readthedocs.io/https://hazma.readthedocs.io/. The icons  and  provide the jupyter notebook and python script used to make each figure, which are paired with https://github.com/mwouts/jupytextjupytext. The code snippets appearing throughout this work are collected in a jupyter notebook: https://github.com/LoganAMorrison/Hazma/blob/master/notebooks/hazma_paper/snippets.ipynb. An archived version of hazma is available on Zenodo . Finally, we have created a Mathematica package, HazmaTools, to reproduce all calculations of cross sections, mediator decay widths, and final state radiation spectra used in Hazma. The package includes FeynRules <cit.> model files for the scalar and vector mediator models, and utilizes FeynArts <cit.> and FeynCalc <cit.>. It is available for download at https://github.com/LoganAMorrison/HazmaTools/https://github.com/LoganAMorrison/HazmaTools/.



  
Conventions: throughout this paper, unless otherwise noted, the units used are MeV for energies, masses and decay widths and c̆m̆^3 / s̆ for the thermally-averaged DM self-annihilation cross section . The numpy package <cit.> is referred to in some of the code snippets as np. Lines in code blocks beginning with >>> indicate the python command prompt. We have sometimes rounded or formatted the output from hazma to make the code snippets more readable.



§ STRUCTURE OF HAZMA





In this section, we describe the structure of the hazma codebase and the intended workflow. The general workflow for a user of hazma is shown in Fig. (<ref>). The light green boxes at the top indicate the types of physics models a user can study, the dark green boxes show possible analysis outputs, and the light green boxes at the bottom denoting inputs required for these outputs.

The user has several options for tapping into the resources provided by hazma. The easiest is to use one of the built-in models, where a user only needs to specify the parameters of the model (see below and Sec. (<ref>) and Sec. (<ref>) for details on the built-in models and their corresponding parameters). Alternatively, if the user is working with a model which is a specific version one of the included models (e.g., where the mediator's couplings to Standard Model particles are interrelated), they can define their own subclass of that model. By using inheritance, one retains all of the functionality of the built-in model (such as functions for computing final state radiation spectra, cross sections, mediator decay widths, etc.) while supplying the user with a simpler, more specialized interface to the underlying models. For a detailed explanations of how to set model parameters and make subclasses of built-in models, see App. (<ref>). Another option is for the user to define their own model. To do this, they need to define a class which contains functions for the gamma-ray and positron spectra, as well as the annihilation cross sections and branching fractions. In App. (<ref>) we provide a detailed example of how to do this for a toy model, utilizing various helper modules provided by hazma.

After choosing a model, the user represents it in hazma by creating an instance of the Theory class (purple box in Fig. (<ref>)). Various particle physics quantities (the DM self-annihilation cross section, gamma-ray spectra per DM annihilation, etc) can be computed using the values of the masses and couplings in the Theory. Computing constraints on the DM self-annihilation cross section requires additional inputs (green boxes at bottom of Fig. <ref>). To make it straightforward for the user to constrain models using existing data, hazma comes with flux measurements from several gamma-ray telescopes, along with Planck's upper bound on how much energy DM annihilations could inject into the CMB. For projecting the reach of future gamma-ray telescopes, the user must provide information characterizing the detector.



Fig. (<ref>) displays the modules contained in hazma as well as the general structure of the code, which are explained in detail below.



  
*theory
The primary goal of hazma is to provide users with a simple interface for setting indirect detection constraints on ⟨σ v ⟩ as a function of DM mass for sub-GeV dark matter models. This is done using the Theory abstract base class, which every dark matter model in hazma should inherit from. At a high level, it contains three methods for determining these constraints: one that uses existing gamma-ray data (Theory.binned_limit), one for computing the discovery reach for planned gamma-ray detectors (Theory.unbinned_limit), and one for deriving CMB constraints (Theory.cmb_limit). The class also provides methods for computing particle physics quantities such as annihilation cross sections (with annihilation_cross_sections), mediator decay widths (partial_widths), the continuum and monochromatic gamma-ray spectra from DM annihilation (total_spectrum, gamma_ray_lines), and the same for positron spectra (total_positron_spectrum, positron_lines), and more.

Custom models must implement functions for computing gamma-ray spectra, positron spectra and branching fractions to gain use of the methods in Theory. As mentioned above, Appendix (<ref>) shows how to implement a simple custom model. Two other modules may be useful for complex custom models: rambo, which performs N-body phase space integrations using a Monte Carlo algorithm, and gamma_ray, which computes gamma ray spectra for N-body final states. These are not used by the simplified models that ship with hazma, and are described in detail in Appendices (<ref>-<ref>).



  
*decay
This contains high-performance functions for computing the decay spectra from π^±, π^0 and μ^±. Details of how these are computed are found in Sec. (<ref>). The functions in this module allow the user to compute the decay spectra for particles decaying with arbitrary lab-frame energies. This requires computing the decay spectra in the rest frame of the parent-particle and performing a Lorentz boost, amounting to performing a change-of-variables along with a convolution integral. To achieve high computational performance, we perform all integrations in c using cython <cit.> and build extension modules to interface with python.



  
*positron_spectra
This module computes the electron/positron spectra from decays of π^± and μ^±, which are critical inputs for constraining dark matter models using CMB observations. See Sec. (<ref>) for details on how these spectra are computed. As in the decay module, the positron_spectra module allows users to compute the electron/positron spectra for arbitrary energies of the parent-particle. The procedure for computing the spectra for arbitrary parent-particle energies is identical to the procedure used for decay.


hazma also ships with various particle physics models of sub-GeV DM. These models are located in the modules scalar_mediator and vector_mediator, and contain all the relevant annihilation cross sections, branching fractions, decay spectra, FSR spectra and positron spectra. They can be used by instantiating the appropriate class (for example, the HiggsPortal class from scalar_mediator for the Higgs-portal model.) The user only needs to specify the parameters of the model. The particle physics frameworks used to construct these models, the specialized subclasses of these models included with hazma, and accessing functions to compute their annihilation cross sections and other particle physics quantities is the topic of the following section.



§ PARTICLE PHYSICS FRAMEWORK



Each of the models distributed with hazma contain two BSM particles: a dark matter particle χ and a mediator M that interacts with the DM as well as Standard Model particles. The Lagrangian can be expressed as

    Ł    = Ł_S̆M̆ + Ł_D̆M̆ + Ł_M̆ + Ł_Ĭn̆t̆(M),

which consists of the SM Lagrangian, the free Lagrangians for the Dark Matter (DM) and mediator (M), and the mediator's interactions with the DM and SM fields, collectively included in the term Ł_Ĭn̆t̆(M). For the models currently available in hazma the DM is taken to be a Dirac fermion, so that

    Ł_D̆M̆    = χ̅ (i -̣ m_χ) χ.

Both the dark matter and the mediator are taken to be uncharged under the Standard Model gauge group. The Lagrangian is defined in terms of the microscopic degrees of freedom of the Standard Model (quarks, leptons and gauge bosons). However, at the energy scale of interest for self-annihilations of non-relativistic MeV dark matter (√(s)∼ 1-100 MeV), quarks and gluons are not the correct strongly-interacting degrees of freedom. At these energies, QCD confines, and quarks and gluons group into bound states called mesons and baryons. Given that the DM interacts with quarks and gluons, at these energies interactions with mesons and baryons are induced. Since QCD is non-perturbative for energies √(s)≲GeV, we must use an effective field theory to describe the interactions of DM with mesons and baryons. We thus match our microscopic Lagrangian onto the effective Lagrangian for pions and other mesons using the techniques of chiral perturbation theory (ChPT) <cit.>. The models currently implemented in hazma utilize leading-order ChPT.

Before describing these models, we review the range of validity of ChPT, which is limited since it is an effective theory. ChPT is organized as an expansion in a small parameter, the meson momentum squared p^2 divided by the squared mass scale associated with loop diagrams, Λ_C̆h̆P̆T̆∼ 4 π f_π≈ 1.2GeV, where f_π = 92.2MeV is the pion decay constant. Derivatives of the meson fields and factors of the meson masses are counted as 𝒪(p). The expansion parameter p^2 is used in two ways. First, the chiral Lagrangian is organized as a sum of terms of decreasing importance

    ℒ_C̆h̆P̆T̆ = ℒ^(2) + ℒ^(4) + ℒ^(6) + …,

where the superscript indicates the number of derivatives or powers of meson masses. Second, the expansion parameter can be used to assess the relative size of an individual Feynman diagram's contribution in the calculation of a given observable. This is quantified using the chiral dimension, which is lower for more important diagrams. Generally the tree-level diagrams from ℒ^(2) have the lowest chiral dimension and provide the largest contributions to observables, followed by tree-level diagrams from ℒ^(4) and one-loop diagrams from ℒ^(2), and so on. As p^2 →Λ_C̆h̆P̆T̆^2, this scheme for organizing calculations in terms of the chiral dimension breaks down and ChPT becomes unreliable.

Assessing the precise range of validity of leading-order ChPT is complicated. A simple estimate suggests that for processes with meson energies of e.g. 500MeV the leading-order observables receive ∼ (500MeV / Λ_C̆h̆P̆T̆)^2 ∼ 20% corrections from next-to-leading-order (NLO) contributions <cit.>. In reality, the magnitude of NLO corrections is larger for annihilation final states with the same quantum numbers as resonances such as the ρ <cit.> and f_0(500) <cit.>, since the final state particles can rescatter. While resonances can be accounted for with unitarization techniques <cit.> or as explicit resonance fields <cit.>, we instead use leading-order ChPT, restricting to DM masses below 250 MeV in the case of annihilation into Standard-Model particles and mediator masses below 500 MeV for analysis of their decays to avoid resonances. An important consequence is that annihilation into kaons and heavier mesons is not considered herein since that would require DM masses far above this range. User-defined models (see App. (<ref>)) need not use leading-order ChPT, and this is an interesting topic for future work.

While the preceding discussion is quite general, we now specialize to the two models that come with hazma: scalar_mediator, which contains a real scalar mediator S, and vector_mediator, where the mediator is a vector V [In future work, we plan to add pseudoscalar and axial-vector mediators.]. These are subclasses of the abstract class Theory, and thus implement a variety of functions for computing physical quantities. In the following two subsections, we present Ł_Ĭn̆t̆(S) and Ł_Ĭn̆t̆(V) at the level of quarks and gluons as well as the Lagrangians obtained by performing the ChPT matching.[The interaction Lagrangians, matching procedure and a review of the chiral Lagrangian are explained in detail in a forthcoming companion paper <cit.>.] Snippets are provided to demonstrate how to construct each model and change its parameters. The third subsection shows how to access various particle physics quantities in hazma.



 §.§ Scalar Mediator



The free Lagrangian for a real scalar is

    Ł_S     = 1/2 (_̣μ S) (^̣μ S) - 1/2 m_S^2 S^2,

where m_S is the scalar's mass. The interactions with the light fundamental SM degrees of freedom read

    Ł_Ĭn̆t̆(S)    = -S ( g_Sχ + g_Sf∑_f y_f/√(2)f̅ f )                                                                                     
        + S/Λ( g_SGα_ĔM̆/4π F_μν F^μν + g_SFα_s/4π G_μν^a G^a μν).

The sum runs over fermions with mass below the GeV scale (f = e, μ, u, d, s). Note that the coupling g_Sf is outside the sum. The Yukawa couplings are defined to be y_f = √(2) m_f / v_h, with the Higgs vacuum expectation value (vev) defined as v_h = 246GeV. The parameter Λ is the mass scale at which S acquires (non-renormalizable) interactions with the photon and gluon.

After performing the matching onto the chiral Lagrangian and expanding to leading order in the pion fields, the resulting interaction Lagrangian is

    Ł_Int(S)    = 2 g_SG/9 Λ S [ (_̣μπ^0) (^̣μπ^0) + 2 (_̣μπ^+) (^̣μπ^-) ]                                                                        
        + 4 i e g_SG/9 Λ S A^μ[ π^- (_̣μπ^+) - π^+ (_̣μπ^-) ]
        - B (m_u + m_d)/6( 3 g_Sf/v_h + 2 g_SG/3 Λ) S [ (π^0)^2 + 2 π^+ π^- ]
        + B (m_u + m_d) g_SG/81 Λ( 2 g_SG/Λ - 9 g_Sf/v_h) S^2 [ (π^0)^2 + 2 π^+ π^- ]
        + 4 e^2 g_SF/9Λ S π^+ π^- A_μ A^μ
        - g_S χ S χ̅χ - g_Sf S ∑_ℓ=e,μy_ℓ/√(2)ℓ̅ℓ.

where B = m_π^±^2/(m_u + m_d) ≈ 2800 MeV <cit.>.

The parameters for the scalar model are attributes of the scalar_mediator class. Their names in hazma are

    (m_χ, m_S, g_Sχ, g_Sf, g_SG, g_SF, Λ) ↔ (𝚖𝚡, 𝚖𝚜, 𝚐𝚜𝚡𝚡, 𝚐𝚜𝚏𝚏, 𝚐𝚜𝙶𝙶, 𝚐𝚜𝙵𝙵, 𝚕𝚊𝚖).

The following snippet shows how to instantiate scalar_mediator, change the value of a parameter, and print its new value:
python
>>> from hazma.scalar_mediator import ScalarMediator
>>> sm = ScalarMediator(mx=150., ms=1e3, gsxx=1., gsff=0.1,
...                     gsGG=0.1, gsFF=0.1, lam=2e5)
>>> sm.gsff
0.1
>>> sm.gsff = 0.5
>>> sm.gsff
0.5

In addition to the general scalar_mediator model, hazma also comes with two subclasses of scalar_mediator, which represent UV completions: a Higgs-portal model <cit.> (HiggsPortal) and a model containing a heavy quark (HeavyQuark.) The Higgs-portal model assumes the scalar mediator interacts with the Higgs through gauge-invariant interactions resulting in the scalar mediator mixing with the Higgs. After diagonalizing the scalar-mediator/Higgs mass matrix, the scalar mediator and Higgs are replaced with:

    h    → hcosθ - Ssinθ, 
            S    → hsinθ + Scosθ

where θ is the scalar-mediator/Higgs mixing angle. This replacement results in the following cutoff scale and couplings of the scalar mediator with the Standard Model fermions, gluons and photon, the later of which require integrating out the τ, c, b, t, W and Z <cit.>:

    g_Sf   = sinθ,     g_SG   = 3sinθ,     g_SF    = -5/6sinθ,    Λ   = v_h.

where v_h = 246  GeV is the Higgs vacuum expectation value. The parameters for the HiggsPortal class are:

    (m_χ, m_S, g_Sχ, sinθ) ↔ (𝚖𝚡, 𝚖𝚜, 𝚐𝚜𝚡𝚡, 𝚜𝚝𝚑𝚎𝚝𝚊).


The heavy-quark model assumes the existence a new heavy, colored and charged fermion which enters the Lagrangian as:

    ℒ⊃ iQ̅DQ - m_QQ̅Q - m_Q/v_hhQ̅Q - g_SQSQ̅Q

with D_μ = ∂_μ -ieQ_QA_μ -ig_sG_μ^aλ^a/2 where Q_Q is the U_ĔM̆(1) charge of the heavy-quark. Integrating out the heavy quark induces a cutoff scale and effective couplings of the scalar mediator with gluons and photons:

    g_SG   = g_SQ,     g_SF   = 2Q_Q^2 g_SQ,    Λ   = m_Q.

Note that integrating out the heavy-quark also induced effective couplings to the SM fermions, but at two-loop order. We do not include these interactions.  The parameters for the HeavyQuark class are:

    (m_χ, m_S, g_Sχ, g_SQ, m_Q, Q_Q) ↔ (𝚖𝚡, 𝚖𝚜, 𝚐𝚜𝚡𝚡, 𝚐𝚜𝚀, 𝚖𝚀, 𝚀𝚀).

Both of these models can be imported and used in a similar fashion to the generic ScalarMediator class. The following snippet displays how to initialize these subclasses:
python
>>> from hazma.scalar_mediator import HiggsPortal, HeavyQuark
>>> hp = HiggsPortal(mx=150., ms=1e3, gsxx=1., stheta=1e-3)
>>> hq = HeavyQuark(mx=150., ms=1e3, gsxx=1., gsQ=1.0, mQ=1e6, QQ=1.)

Trying to set the underlying ScalarMediator model's attributes will result in an error when these are fully determined by the attributes of the subclass:
python
>>> hp.gsff        # can only be accessed
0.001
>>> hp.gsff = 0.1  # fully determined by hp.stheta
AttributeError: Cannot set gsff




 §.§ Vector Mediator


For the vector mediator the free part of the Lagrangian is

    Ł_V     = -1/4 V_μν V^μν + 1/2 m_V^2 V_μ V^μ,

where m_V is the mass of the vector. The interactions considered are

    Ł_Ĭn̆t̆(V)    = V_μ( g_Vχχ̅γ^μχ + ∑_f g_Vff̅γ^μ f ).

The sum again runs over the light fermions (f = e, μ, u, d, s), and V may have different couplings to each of these.

A kinetic mixing term between the mediator and photon is also possible:

    Ł_Ĭn̆t̆(V)   ⊃ - ϵ/2 V^μν F_μν.

However, this can be eliminated by redefining the photon field using A_μ→ A_μ - ϵ V_μ. Upon this field redefinition V acquires an ϵ-suppressed interaction with the SM fermions, which is captured by changing the fermion couplings

    g_Vf→ g_Vf - ϵ e Q_f,

where Q_f is the electric charge of the fermion f and e > 0 is the electron's charge. Because of this, we do not include ϵ as a coupling in the vector model, and instead provide a subclass KineticMixing to handle this case.

Matching onto the chiral Lagrangian and isolating the terms contributing at leading order to the quantities computed in hazma gives

    Ł_Int(V)    = -i (g_Vu - g_Vd) V^μ( π^+ _̣μπ^- - π^- _̣μπ^+ ) 
        + ( g_Vu - g_Vd )^2 V_μ V^μπ^+ π^- 
        + 2 e (Q_u - Q_d) (g_Vu - g_Vd) A_μ V^μπ^+ π^-
        + V_μ(g_Vee̅γ^μe + g_Vμμ̅γ^μμ)
        + 1/8π^2 f_πϵ^μνρσ (_̣μπ^0)
       ×{ e (2 g_Vu + g_Vd) [ (_̣ν A_ρ) V_σ + (_̣ν V_ρ) A_σ] .
       . + 3 (g_Vu^2 - g_Vd^2) (_̣ν V_ρ) V_σ}.

The contributions in the last three lines come from matching onto the anomalous terms in the chiral Lagrangian <cit.>, which is explained in detail in our companion paper <cit.>. While these terms come from the NLO chiral Lagrangian, the resulting final states contribute significantly to the gamma-ray spectra, so we include them here. In particular, in the case where V couples only to quarks, the only final state accessible below the two-pion threshold is π^0γ. In contrast, other terms from the NLO chiral Lagrangian only lead to corrections for scattering rates into final states accessible at tree-level, and are thus not included. 

The correspondence between the microscopic parameters for the vector model and attributes of the base vector_mediator class is[Note that gvss is not currently used since vector_mediator is based on leading-order ChPT.]

    (m_χ, m_V, g_Vχ, g_Vu, g_Vd, g_Vs, g_Ve, g_Vμ)                                                           
       ↔ (mx, mv, gvxx, gvuu, gvdd, gvss, gvee, gvmumu).

The following snippet shows how to instantiate vector_mediator:
python
>>> from hazma.vector_mediator import VectorMediator
>>> vm = VectorMediator(mx=150., mv=1e3, gvxx=1., gvuu=0.1,
...                     gvdd=0.2, gvss=0.0, gvee=0.4,
...                     gvmumu=0.5)


For the subclass KineticMixing which handles the important case of a kinetically-mixed mediator, the parameters are

    (m_χ, m_V, g_Vχ, ϵ) ↔ (mx, mv, gvxx, eps).

While the underlying parameters gvuu, …, gvmumu can be accessed by instances of KineticMixing, they cannot be set directly since they are fully determined by eps:
python
>>> from hazma.vector_mediator import KineticMixing
>>> km = KineticMixing(mx=150., mv=1e3, gvxx=1., eps=0.1)
>>> km.gvuu
-0.020187846690459792  # = -0.1 * 2/3 * sqrt(4 pi / 137)
>>> km.gvuu = 0.1
AttributeError: Cannot set gvuu


Finally, the subclass QuarksOnly is provided for analyzing the case where the mediator is hadrophilic and has only couplings to quarks. The parameters are those for the full vector-mediator model, with the leptonic couplings set to zero.



 §.§ Computing Particle physics quantities







Every model in hazma is required to implement functions listing the available annihilation final states and providing corresponding functions for computing annihilation cross sections. For example, for the models described above we have:
python
>>> ScalarMediator.list_annihilation_final_states()
['mu mu', 'e e', 'g g', 'pi0 pi0', 'pi pi', 's s']
>>> VectorMediator.list_annihilation_final_states()
['mu mu', 'e e', 'pi pi', 'pi0 g', 'pi0 v', 'v v']

These lists exclude final states where the annihilation is extremely suppressed by couplings or phase-space factors (e.g., π^+ π^- γγ, S π^0 π^0, S π^+ π^- for the scalar model). They also exclude final states arising from radiative processes (such as e^+ e^- γ). Depending on the couplings' values and the center of mass energy of the DM, some of these final state may be inaccessible.

The corresponding annihilation cross sections and branching fractions can be accessed using the annihilation_cross_sections function:
python
>>> from hazma.scalar_mediator import ScalarMediator
>>> sm = ScalarMediator(mx=180., ms=190., gsxx=1., gsff=0.1,
...                     gsGG=0.1, gsFF=0.1, lam=2e5)
>>> e_cm = 400.
>>> sm.annihilation_cross_sections(e_cm)  # MeV^-2
'mu mu': 1.018678054222354e-16,
 'e e': 3.892649866478948e-21,
 'g g': 2.4381469306302895e-21,
 'pi0 pi0': 7.975042237472552e-17,
 'pi pi': 1.5169347389281449e-16,
 's s': 1.78924656025887e-07,
 'total': 1.7892465635920504e-07

Note that in this example m_χ < m_S, but the center of mass energy is large enough to allow the process χ̅χ→ S S. The function Theory.annihilation_branching_fractions() is invoked the same way and returns a dict of branching fractions for annihilation into each final state X:

    BF(χ̅χ→ X) = σ_χ̅χ→ X/∑_Yσ_χ̅χ→ Y


Fig. (<ref>) shows the annihilation branching fractions for the scalar model. The different rows exhibit results for the Higgs portal and heavy quark UV completions discussed above, which will be referred to throughout this work. In the first case, the branching fraction is largest into whichever final state is closest to threshold since S couples to Standard Model states roughly proportionally to their Yukawas. The branching fraction into γγ is small for non-relativistic DM annihilations due to the derivative coupling between S and the photon field. For the heavy quark coupling pattern only hadronic final states are accessible since the mediator couplings exclusively to photons and gluons.

Fig. (<ref>) collects the branching fractions for the vector model, assuming it kinetically mixes with the photon (top row) or couples only to quarks (middle and bottom rows). The coupling-dependence of the vector's cross sections can be understood without performing detailed calculations. Since the vector model can be recast as a U(1) gauge theory (by using the Stuckelberg mechanism to generate the mass term and decoupling the Stuckelberg field), the vector's couplings to pions are the sums of its couplings to their constituents: g_Vπ^±π^± = g_Vuu - g_Vdd, g_Vπ^0π^0 = 0. The cross section σ_χ̅χ→π^+ π^- is thus proportional to (g_Vuu - g_Vdd)^2, explaining why that final state has a branching fraction of zero in the bottom row of the figure. The vector's coupling to π^0γ and π^0 V come from the anomalous part of the chiral Lagrangian. Since the anomaly is exact at one loop, the relevant diagrams in the former case can be evaluated with a proton running in the loop, making it clear that σ_χ̅χ→π^0γ∝ (2 g_Vuu + g_Vdd)^2. This cross section is suppressed due to the loop factors, but provides an important contribution to the gamma-ray spectrum. For the coupling choices used for the plots in this work, the spectrum contribution from the π^0 V final state is either negligible or identically zero.

Finally, the mediator decay partial widths are easily obtained with the partial_widths function:
python
>>> from hazma.scalar_mediator import ScalarMediator
>>> sm = ScalarMediator(mx=120., ms=280., gsxx=1., gsff=0.1,
...                     gsGG=0.1, gsFF=0.1, lam=2e5)
>>> sigmas = sm.partial_widths()
>>> sigmas  # MeV
'g g': 1.472617003459079e-13,
 'pi0 pi0': 3.261686997076484e-09,
 'pi pi': 1.8664869506864194e-09,
 'x x': 1.522436123428156,
 'e e': 4.798160511838726e-13,
 'mu mu': 5.792897882600423e-09,
 'total': 1.522436134349855




§ BUILDING BLOCKS OF MEV GAMMA-RAY SPECTRA



The total gamma-ray spectrum for dark matter annihilation consists of three parts:

    . dN/dE_γ|_χ̅χ   = . dN/dE_γ|_χ̅χ,line + . dN/dE_γ|_χ̅χ,dec. + . dN/dE_γ|_χ̅χ,FSR.

The first term accounts for photons from final states containing one or more monochromatic gamma rays, which result in spectral lines. The second term is the spectrum of photons produced through the radiative decays of final state particles. Decays can produce dramatic spectral features such as the neutral pion “box” as well as E_γ^-1 spectra. The third term accounts for photons produced through final state radiation (FSR) from annihilation into electromagnetically charged final states, and possibly internal bremmstrahlung (IB). FSR generically produces spectra that scale as dN/dE_γ∝ E_γ^-1 at low energies, in accordance with Low's theorem <cit.>, which could be distinguished from the softer dN/dE_γ∝ E_γ^-Γ, with Γ∼2 for typical astrophysical background.

The calculations required to account for these contributions are described in detail in this section, beginning with the trivial case of spectral lines and the associated hazma functions. The model-dependent FSR spectra are explained in detail since we compute them exactly for leptonic and hadronic final states. With regard to radiative decays, we present a comprehensive overview of the contributions to the π^± decay spectrum, which is currently absent from the literature and thus has not been considered in prior work on sub-GeV DM. Code snippets for computing different decay spectra are included throughout.

The section concludes by showing how to compute FSR and radiative decay spectra in hazma for the built-in models. A key feature of hazma is that the spectrum functions can be employed to analyze user-defined models, as demonstrated in a detailed example in App. (<ref>).



 §.§ Monochromatic Gamma Rays



The only final states containing monochromatic photons that are relevant for this work are π^0 γ and γγ:

    . dN/dE_γ|_χ̅χ,line (E_γ)     = Br(χ̅χ→π^0 γ)  δ(E_π^0 γ - E_γ) + 2 Br(χ̅χ→γγ)  δ(E_CM / 2 - E_γ),

where the energy of the line from the π^0 γ final state is E_π^0 γ≡ (E_CM^2 - m_π^0^2) / (2 E_CM) and E_CM is the center of mass energy. This contribution to the DM annihilation spectrum is thus obtained by computing the branching fractions appearing in this expression.

The Theory.gamma_ray_lines() method returns a dict with information about monochromatic gamma-ray lines produced in dark matter annihilations. For example, the scalar mediator model has a 2γ line and the vector model has a line from the π^0γ final state. Assuming heavy quark-type couplings for the former and quark-only couplings for the later, the line energies and branching fractions are found to be
python
>>> from hazma.scalar_mediator import HeavyQuark
>>> sm = HeavyQuark(mx=140., ms=1e3, gsxx=1., gsQ=0.1, mQ=1e3, QQ=0.1)
>>> sm.gamma_ray_lines(e_cm=300.)
'g g': 'energy': 150.0, 'bf': 1.285653242415087e-08
>>> from hazma.vector_mediator import QuarksOnly
>>> vm = QuarksOnly(mx=140., mv=1e3, gvxx=1., gvuu=0.1, gvdd=0.1, gvss=0.)
>>> vm.gamma_ray_lines(e_cm=300.)
'pi0 g': 'energy': 119.63552908740002, 'bf': 1.0




 §.§ Final State Radiation


Dark matter annihilating into charged SM particles generically produces photons via final state radiation (FSR), which is accounted for by computing

    . dN/dE_γ|_χ̅χ,FSR/IB   = 1/σ_χ̅χ∑_A dσ_χ̅χ→ A γ/dE_γ(E_γ)
       ≈∑_A [ σ_χ̅χ→ A/σ_χ̅χ] [ 1/σ_χ̅χ→ Adσ_χ̅χ→ A γ/dE_γ(E_γ) ]
       ≡∑_A Br(χ̅χ→ A) . dN/dE_γ|_χ̅χ→ A (E_γ),

where in the second line we assumed σ_χ̅χ→ Aγ≪σ_χ̅χ→ A and the sum is over final states containing electromagnetically-charged particles. The required branching fractions were computed in Sec. (<ref>).

It is well-known that the FSR physics factorizes from the hard process in the limit where the mass of the radiating particle is zero and the radiation is collinear with the radiating particle<cit.>. This leads to a universal spectrum for FSR for DM annihilating into charged particle-antiparticle pairs X̅ X, which, at leading order in α is: [Note that the approximation takes several forms in the literature, often omitting the -1 term, which are equivalent in the small-μ limit.]

    . dN/dE_γ|_χ̅χ→X̅ X   ≈2 α/π Q· P_γ X(x) ·[ log( 1 - x/μ_X^2) - 1 ].

with μ_X = m_X/Q is the ratio of the mass of the radiating particle divided by the center of mass energy. The form factor P_γ X(x) is called the splitting function, and depends only on the spin of the radiating particle <cit.> and the fraction of energy it takes from the splitting particle, x: [Note that the scalar splitting function in Ref. <cit.> is a factor of two smaller than in Ref. <cit.>. We have rederived the splitting function and agree with the result in the later reference.]

    P_γ X(x) = 2(1 - x)/x   for  X=s̆c̆ăl̆ăr̆
    1 + ( 1 - x )^2/x   for  X=f̆ĕr̆m̆ĭŏn̆,

where x ≡ 2 E_γ / Q.

This spectrum, which we call the Altarelli-Parisi (AP) approximation, has been employed in prior studies of indirect detection of sub-GeV DM <cit.>. However, the AP approximation insufficient for our analysis since it is only valid near the mass singularity, i.e. for μ_X≪1. This does not hold for MeV DM in the Milky Way halo annihilating non-relativistically into μ^+ μ^- (or π^+ π^-), since in this case Q ∼ m_μ∼𝒪(100 MeV). As a result of this assumption, the AP spectrum cuts off at x = 1- exp(1)μ^2, which is different from the exact kinematic threshold x = 1-4μ^2.

In this section we instead exactly compute the model-dependent FSR spectra for the e^+ e^- and μ^+ μ^- final states, as well as the π^+ π^- final state. The resulting spectra for leptonic final states are (see chapter 20 of <cit.> for details of how these calculations can be performed)

    . dN/dE_γ|_χ̅χ→ S^* →ℓ̅ℓ   = α/E_γπ (1 - 4 μ^2)^3/2[ ( 2 (1 - x - 6 μ^2) + (x + 4 μ^2)^2 ) log1 + √(1 - 4 μ^2/1-x)/1 - √(1 - 4 μ^2/1-x). 
       . - 2 (1-4μ^2)(1-x) √(1 - 4 μ^2/1-x)]
    . dN/dE_γ|_χ̅χ→ V^* →ℓ̅ℓ   = α/E_γπ√(1 - 4μ^2) (1 + 2 μ^2)[ (1 + (1-x)^2 - 4μ^2 (x + 2μ^2)) log1 + √(1 - 4 μ^2/1-x)/1 - √(1 - 4 μ^2/1-x). 
       . - (1 + (1 - x)^2 + 4μ^2 (1-x)) √(1 - 4 μ^2/1-x)].

The FSR spectrum for the π^+ π^- final state can be computed using the interaction Lagrangians derived in our companion paper <cit.> and reproduced in Eqn. (<ref>) and (<ref>). The final expressions for the spectra are

    .dN/dE_γ|_χ̅χ→ S^* →π^+π^-   = 2 α/E_γπ√(1 - 4 μ^2)[ (1 - x - 2μ^2) log1 + √(1 - 4 μ^2/1-x)/1 - √(1 - 4 μ^2/1-x) - (1 - x) √(1 - 4 μ^2/1-x)]
    . dN/dE_γ|_χ̅χ→ V^* →π^+π^-   = 2 α/E_γπ (1 - 4μ^2)^3/2[ (1 - x - 2μ^2)(1 - 4μ^2) log1 + √(1 - 4 μ^2/1-x)/1 - √(1 - 4 μ^2/1-x). 
       . - ( (1 - x)(1 - 4 μ^2) - x^2 ) √(1 - 4 μ^2/1-x)].

In the limit μ≪ 1, these spectra reduce to the AP approximation, up to the additional, non-collinear terms: α x / π and 2 α x / π for the processes χ̅χ→ S^* →ℓ̅ℓγ and χ̅χ→ V^* →π^+π^-γ respectively. The spectra all scale as E_γ^-1 in accordance with Low's theorem at low photon energies.

Computing the FSR spectra for annihilation into π^+ π^- requires including internal bremmstrahlung diagrams with the 4-point vertices S,V - π^+ - π^- - γ. In the scalar mediator case, this diagram is only present when g_SG≠ 0. The FSR spectrum is in accord with the AP approximation in the small μ limit. In the vector case, this diagram cannot be eliminated since it is required by gauge invariance.





In Fig. (<ref>) we plot the spectra for FSR in the scalar and vector models as well as the AP spectrum (blue, red and yellow curves respectively). As can be seen from the above equations, the specific values of the couplings, mediator masses and DM mass do not impact the shape of the spectrum since they cancel when dividing by σ_χ̅χ→ℓ̅ℓ. As expected, the AP approximation performs poorly near the dimuon threshold (upper left panel). The photon energy at which the spectrum cuts off is an order of magnitude too large, and while the normalization is within a factor of two for the scalar model's spectrum, it is an order of magnitude larger than the vector model's one. The situation improves at larger center of mass energies, and just above the ππ threshold the spectra all agree to within a factor of two. The pion FSR spectra are shown in Fig. (<ref>).

We demonstrate how to compute these spectra in hazma at the end of this section.



 §.§ Radiative Decay Spectra



Photons are produced when dark matter annihilates into a state A that subsequently undergoes radiative decay A → B γ:

    . dN/dE_γ|_χ̅χ,dec. (E_γ)    = ∑_A Br(χ̅χ→ A) . dN/dE_γ|_A (E_γ).

In addition to collecting well-known expressions for the neutral pion and muon radiative decay spectra, we carefully account for the three dominant contributions to the charge pion's radiative decay spectrum. Functions for computing radiative decay spectra are contained in the decay module, and described below. These can be utilized to compute the decay contribution to the annihilation spectrum for arbitrary final states.



  §.§.§ Decay Spectra in Different Frames


Before computing the radiative decay spectra for the μ^±, π^0 and π^±, we review how to boost a spectrum dN/dE_R from the parent particle's rest frame to obtain the “lab frame” spectrum dN/dE_L, where the particle has boost γ = E / m along the z-axis. Let the photon phase space coordinates be (E_R, cosθ_R) in the rest frame and (E_L, cosθ_L) in the lab frame, where θ_R and θ_L are the angles with respect to the z axis and we have assumed azimuthal symmetry. Since the number of photons in a small patch of phase space does not depend on which coordinates we choose, we have

    dN    = f(E_R, cosθ_R) dE_R dcosθ_R
       = dE_R/dE_L   dE_R/dcosθ_L
    dcosθ_R/dE_L   dcosθ_R/dcosθ_L f(E_R(E_L, cosθ_L), cosθ_R(E_L, cosθ_L)) dE_L dcosθ_L,

where we recognize the partial phase space density as f(E, cosθ) ≡ dN/(dE dcosθ) and the term with the two vertical lines is the Jacobian factor for the change of variables. If f(E_R, cosθ_R) is independent of cθ_R (which is the case for all the spectra we will consider), we can write:

    f(E_R, cθ_R) = 12dNdE_R.

The variables E_1, E_2, cθ_1 and cθ_2 are related via a Lorentz boost. The relationships are

    E_R    = γ E_L (1 - β cθ_L),    cosθ_R    = β - cosθ_L1 - βcosθ_L,

where β = √(1 - 1 / γ^2) is the particle's velocity in natural units. The spectrum in the lab frame is then obtained by integrating Eqn. (<ref>) over cosθ_L:

    dNdE_L    = ∫cosθ_L12γ(βcosθ_L - 1)dNdE_R.




  §.§.§ Neutral Pions

The dominant decay mode for neutral pions is π^0→γγ with a branching fraction of about 99%. Due to this decay modes' large branching fraction and the fact that it has two photons in the final state, we ignore the π^0's other decay modes. In the pion's rest frame, the gamma-ray spectrum is trivial:

    . dNdE_γ|_π^0(E_π^0 = m_π^0)     = 2×δ(E_R - m_π2)

where the factor of 2 comes from the fact that there are two photons in the final state.

Applying Eqn. (<ref>) gives that the gamma ray spectrum in the laboratory frame is

    . dNdE_γ|_π^0(E_π^0)    = ∫cosθ_L12γ(βcosθ_L - 1) 2×δ(E_R - m_π^02)
       = 2γβ m_π^0[θ(E_γ - E_-)- θ(E_γ -E_+)]

with E_± = m_π^0 / 2γ(1∓β), which is the characteristic box spectrum centered at m_π^0 / 2.

The boosted spectrum can be computed in hazma as follows:
python
>>> from hazma.decay import neutral_pion as dnde_pi0
>>> e_gams = np.array([100., 125., 150])  # photon energies
>>> e_pi0 = 180.                          # pion energy
>>> dnde_pi0(e_gams, e_pi0)
array([0.0165965, 0.0165965, 0.       ])




  §.§.§ Muons


The primary muon decay channel is μ^- → e^- ν̅_e ν_μ. The photon in the corresponding radiative decay μ^- → e^- ν̅_e ν_μγ has the following spectrum in the muon's rest frame <cit.>:

    dNdE_γ|_μ^±(E_μ = m_μ)    = α (1-x)/36 π E_γ[ 12 (3 - 2 x (1-x)^2) log(1-x/r) . 
       . + x (1-x) (46 - 55x) - 102 ],

where r ≡(m_e/m_μ)^2, x ≡ 2 E_γ / m_μ and the kinematic bounds on E_γ translate into 0 ≤ x ≤ (1 - r). To obtain the spectrum in the lab frame we substitute the above expression into Eqn. (<ref>) and evaluate the integral numerically. 

The lab-frame decay spectrum can be compute for arbitrary muon energies using:
python
>>> from hazma.decay import muon as dnde_mu
>>> e_gams = np.array([1., 10., 100.])  # photon energies
>>> e_mu = 130.                         # muon energy
>>> dnde_mu(e_gams, e_mu)
array([1.76076858e-02, 1.34063877e-03, 4.64775301e-08])




  §.§.§ Charged Pions


The dominant pion decay is π^+ →ℓ^+ ν_ℓ, where ℓ = e, μ and ℓ is produced approximately on shell, thanks to the narrow width approximation when ℓ = μ. Contributions to the pion's radiative decay spectrum come from initial state radiation from the photon or FSR from the lepton, as well as emission from virtual hadronic states (the “structure-dependent” component). An expression for d^2 Γ_π^+ →ℓ^+ νγ / dx dy is given in Eqn. (72) of <cit.> and Eq. (68.3) in the Review of Particle Physics <cit.> [Note that the Particle Data Group uses the convention f_π = 130 MeV, which is a factor of √(2) larger than the value used in this work.], where x ≡ 2 E_γ / m_π^+ and y ≡ 2 E_ℓ / m_π^+. After changing to the Mandelstam variables

    s    ≡ (p_π - p_γ)^2 = m_π^+^2 (1-x),
    
        t    ≡ (p_π - p_ℓ)^2 = m_π^+^2 ( 1 - y + m_ℓ^2/m_π^+^2),

integrating over 0 ≤ t ≤ (m_π^+^2 - s) (s - m_ℓ^2) / s, and dividing by the total π^+ decay width, we obtain an analytic expression for this channel's contribution to the π^+ decay spectrum:

    Br(π^+ →ℓ^+ ν_ℓ) ·. dN/dE_γ|_π^+ →ℓ^+ ν_ℓ (E_π^+ = m_π^+)    = Γ_π^+ →ℓ^+ ν_ℓ/Γ_π^+·α (f(x) + g(x))/24 π m_π^+ f_π^2 (r-1)^2 (x-1)^2 r x,

where r ≡ m_ℓ^2 / m_π^+^2 and

    f(x)     = (r+x-1) [m_π^+^2 x^4 (F_A^2+F_V^2 ) (r^2-r x+r-2 (x-1)^2) . 
       . - 12 √(2) f_π m_π^+ r (x-1) x^2 (F_A (r-2 x+1)+x F_V) . 
       . - 24 f_π^2 r (x-1) (4 r (x-1)+(x-2)^2)], 
    
        g(x)     = 12 √(2) f_π r (x-1)^2 log(r/1-x) [ m_π^+ x^2 (F_A (x-2 r) - x F_V) . 
       . + √(2) f_π(2 r^2-2 r x-x^2+2x-2) ].

The vector form factor is F_V(q^2) = F_V(0) (1 + a q^2) with F_V(0) = 0.0254, slope parameter a = 0.10 and q^2 = 1 - x; the axial form factor is F_A = 0.0119 <cit.>.

When ℓ = μ, the muon's subsequent radiative decay also contributes significantly to the pion's decay spectrum.  Since we can take the muon to be on shell by the narrow width approximation, this decay path's contribution is simply the product of the branching fraction for π^+ →μ^+ ν_μ and the muon decay spectrum Eqn. (<ref>) evaluated at the muon's energy. Our final expression for the charged pion decay spectrum is thus

    dNdE_γ|_π^+(E_π^+ = m_π^+)    = ∑_ℓ = e, μBr(π^+ →ℓ^+ ν_ℓ) ·. dN/dE_γ|_π^+ →ℓ^+ ν_ℓ (E_π^+ = m_π^+)
        + Br(π^+ →μ^+ ν_μ) ·dNdE_γ|_μ^±( E_μ = m_π^+^2 + m_μ^2/2 m_π^+).

This rest frame spectrum can be substituted into Eqn. (<ref>) and numerically integrated to obtain the total charged pion radiative decay spectrum. 



Fig. (<ref>) shows the contributions of each term and the total spectrum in the pion's rest frame. The muon's radiative decay is by far the most important component due to its large branching fraction. There is also more phase space available for the photon in comparison with π^+ →μ^+ ν_μγ since the other final state particles are massless: the former spectrum cuts off at ∼ 69MeV and the latter at (m_π^+^2 - m_μ^2) / (2 m_π^+) ≈ 27MeV. The contribution from π^+ → e^+ νγ is smaller than the others due to helicity suppression.

The charged pion decay spectrum is included in hazma.decay:
python
>>> from hazma.decay import charged_pion as dnde_pi
>>> e_gams = np.array([1., 10., 100.])  # photon energies
>>> e_pi = 150.                         # pion energy
>>> dnde_pi(e_gams, e_pi)
array([1.76949944e-02, 1.32675207e-03, 1.16607174e-09])

The individual contributions can also be computed by setting the mode argument. The orange, blue and green curves are obtained using mode="munu", "munug" and "enug" respectively.



 §.§ Continuum Gamma-Ray Spectra in hazma



Models in hazma provide a few methods for computing gamma-ray spectra at different levels of detail: accessing the total gamma-ray spectrum including all available final states, computing spectrum from specific final states or even decomposing the spectra into FSR or decay for individual channels and processes. The total_spectrum method gives the total continuum gamma-ray spectrum at specified photon energies and fixed center-of-mass energy:
python
>>> from hazma.scalar_mediator import HiggsPortal
>>> e_cm = 305.                         # DM center of mass energy
>>> e_gams = np.array([1., 10., 100.])  # photon energies
>>> hp = HiggsPortal(mx=150., ms=1e3, gsxx=0.7, stheta=0.1)
>>> hp.total_spectrum(e_gams, e_cm)
array([0.02484114, 0.00186874, 0.01827828])

Underlying this are methods for computing the gamma-ray spectra for individual final states. These are accessible through the spectra method, which also accepts a list of photon energies and a center of mass energies. It returns the total continuum spectrum and the contribution from each final state (ie, dN/dE_γ multiplied by the final state's branching fraction):
python
>>> hp.spectra(e_gams, e_cm)
'mu mu': array([4.87977497e-03, 3.83154301e-04, 1.22977086e-06]),
 'e e': array([4.19129970e-07, 3.93151015e-08, 2.22413275e-09]),
 'pi0 pi0': array([0.        , 0.        , 0.01827702]),
 'pi pi': array([1.99609455e-02, 1.48554303e-03, 2.28115016e-08]),
 's s': array([0., 0., 0.]),
 'total': array([0.02484114, 0.00186874, 0.01827828])
 
 The spectra and total_spectrum methods are provided by the Theory class, and are driven by the method spectrum_funcs which classes inheriting from Theory are required to implement. This returns a dict whose keys are strings corresponding to annihilation final states and whose values are methods returning the continuum gamma-ray spectrum for annihilations into that final state (ie, omitting the branching fraction factor).

The models built into hazma provide more fine-grained methods for studying different channels' spectra. Each of these has a name corresponding to the final state and takes an additional string argument specifying the spectrum type ("decay", "fsr" or "all", the default). The return value depends on the spectrum type argument. In the first case, the returned spectrum only accounts for the final state particles' radiative decays, which are model-independent. In the second case, the method returns the radiative decay spectrum, which is model-dependent. In the third, default case, the method returns the total continuum spectrum.

For example, in the scalar mediator model, the annihilation final states e^+ e^-, μ^+ μ^-, π^0π^0, π^+ π^- and S S contribute to the continuum annihilation spectrum. The following snippet shows how to call the spectrum methods for the instance of HiggsPortal created above. The methods all follow the same naming conventions and return dN/dE_γ in M̆ĕV̆^-1:
python
>>> from hazma.scalar_mediator import HiggsPortal
>>> e_cm = 305.                         # DM center of mass energy
>>> e_gams = np.array([1., 10., 100.])  # photon energies
>>> hp = HiggsPortal(mx=150., ms=1e3, gsxx=0.7, stheta=0.1)
>>> hp.dnde_ee(e_gams, e_cm, spectrum_type="fsr")
>>> array([0.05435176, 0.00509829, 0.00028842])
>>> hp.dnde_ee(e_gams, e_cm, spectrum_type="decay")
>>> array([0., 0., 0.])  # electrons don't decay
>>> hp.dnde_pipi(e_gams, e_cm, spectrum_type="fsr")
>>> array([1.01492577e-03, 5.00245201e-05, 0.00000000e+00])
>>> hp.dnde_pipi(e_gams, e_cm, spectrum_type="decay")
>>> array([3.53972084e-02, 2.65985674e-03, 4.16120296e-08])
>>> hp.dnde_pipi(e_gams, e_cm)
>>> array([3.64121342e-02, 2.70988126e-03, 4.16120296e-08])




§ GAMMA RAY SPECTRA FROM DM ANNIHILATION











Using the ingredients described in the previous section, we now compute the photon spectra from DM annihilation for the scalar and vector models with the same couplings as in Fig. (<ref>) and Fig. (<ref>). In what follows, the center of mass energy is set to the non-relativistic approximation E_CM = 2 m_χ( 1 + 1/2 v_MW^2 ), where v_MW = 10^-3 c is the fiducial velocity dispersion for the Milky Way.

Fig. (<ref>) displays the spectrum for the scalar model when annihilation into the mediator is kinematically prohibited (m_S ≳ m_χ). The dark matter mass is fixed for each column and the model parameters in each row. Individual final states' spectra are indicated by the colored curves. To visualize the di-photon final state's contribution, the black curve shows the total spectrum convolved with a 5% energy resolution function (the case for COMPTEL <cit.>), with a vertical dashed line highlighting the photons' energies (E_CM / 2).

For low DM masses, the electron and di-photon line are responsible for the spectrum. Above the muon threshold the Yukawa-suppressed electron contribution is negligible. Since the branching fractions into different pion species are identical (up to isospin-breaking corrections proportional to m_π^± - m_π^0), the charged pion spectrum is always sub-dominant to the neutral pion box for E_CM > 2 m_π. Depending on the ratio of the microscopic coupling to photons (g_SF) with the couplings to gluons and fermions (g_SG and g_Sf), the neutral pion box or di-photon line are generally the most distinctive spectral features for large DM masses, which can also be seen in the figure. Note that if g_Sff and g_SGG are considered as independent parameters, they can take values such that the cross section for annihilation into pions is zero, though this scenario is very fine-tuned.

Spectra for annihilation into mediators are shown in Fig. (<ref>). The panels are labeled with the DM mass and the different colored curves are the spectra computed with the scalar's mass set to the indicated values. For different mediator masses, the following processes dictate the shape of the spectrum:

    
  * m_S < 2 m_e: the mediator can only decay to two photons, so the annihilation spectrum is a box centered at m_S/2, just like for the neutral pion.
    
  * 2 m_e < m_S < 2 m_μ: the spectrum is a superposition of the box and the electron FSR spectra, and when the scalar can decay into two muons this dominates the spectrum.
    
  * m_S > 2 m_π: the spectrum consists of a peak from boosting the π^0 decay spectrum and a softer component from FSR and decays from the other final states into which S decays. The peak is more rounded if the scalar is highly boosted (m_S ≪ E_CM / 2) and boxy for m_S ∼ E_CM / 2. The difference in the case of heavy quark microscopic couplings is that the muon and electron final states are not present.


The components of the spectrum for DM annihilation into SM final states in the vector model are illustrated in Fig. (<ref>). In the kinetic mixing case the spectrum consists purely of a continuum, since the branching fraction for annihilation into π^0 γ is highly suppressed. In contrast, if the vector couples only to quarks the spectrum is composed of a monochromatic gamma ray line and box from the π^0 decay at low energies. If not forbidden by the vector's couplings (as is the case for the bottom row), as the center-of-mass energy grows above the two-pion threshold the continuum of photons from the π^+ π^- final state overwhelm those from π^0 γ.

Fig. (<ref>) exhibits spectra from DM annihilating into V V followed by the vectors' decays for a kinetically mixed vector (top row) and couplings to quarks only. As in the case of annihilations directly into SM particles, a power law spectrum is obtained in the kinetically mixed case. For the spectrum for quark-only couplings (middle and bottom rows) is sourced exclusively by the vector's decay to π^0 γ. The resulting spectrum is comprised of a box (from the boosting the monochromatic photon) and a rounded peak (from boosting the neutral pion's decay spectrum). The box component is offset towards higher energies than the rounded peak component, leading to the sawtooth pattern in the total spectrum. This feature would only by observable at a detector with very good energy resolution.

The colored curves in the Fig. (<ref>)-(<ref>) can be reproduced using the Theory.spectra() function described in Sec. (<ref>). The total spectrum convolved with an energy resolution function can be obtained using Theory.total_conv_spectrum_fn(). The arguments specify the range of photon energies over which to perform the convolved spectrum, the center of mass energy of the DM annihilation, the energy resolution function, and the number of points to use when computing the convolution:
python
>>> from hazma.scalar_mediator import HiggsPortal
>>> e_cm = 305.                  # DM center of mass energy
>>> e_min, e_max = 1., 100.      # define energy range
>>> energy_res = lambda e: 0.05  # 5
>>> hp = HiggsPortal(mx=150., ms=1e3, gsxx=0.7, stheta=0.1)
>>> dnde_conv = hp.total_conv_spectrum_fn(
...     e_min, e_max, e_cm, energy_res, n_pts=1000)

Decreasing n_pts from its default value of 1000 increases speed at the cost of accuracy. The return value (dnde_conv above) is an interpolator of type InterpolatedUnivariateSpline from the scipy.interpolate module. This can be used to compute the spectrum at a particular photon energy:
python
>>> dnde_conv(25.)
array(0.00048767)

This spectrum is can easily and efficiently be integrated over a range of photon energies between e_min and e_max:
python
>>> dnde_conv.integral(25., 85.)
0.815998464406668




§ ELECTRON AND POSITRON SPECTRA FROM DM ANNIHILATION



Dark matter annihilation in the mass range under consideration produces relativistic electrons and positrons. Such particles could be in principle directly detected as part of the cosmic radiation, albeit at energies well below the GeV solar modulation strongly affects cosmic-ray electron spectra. Additionally, relativistic electrons and positrons can be indirectly detected through (1) effects on photons of the cosmic microwave background and (2) the secondary emission of radiation via e.g. up-scattering of background photons (inverse Compton processes), synchrotron radiation, and bremsstrahlung. The calculation of the electron-positron spectra for MeV dark matter is therefore needed to reliably compute CMB limits and the spectrum of secondary radiation.

In this section we compute the positron spectrum for MeV dark matter, which follows along the lines of the gamma-ray ones.[The electron spectrum looks identical to the positron one since our dark matter particles do not carry lepton number.] The spectrum consists of a line-like spectrum from χ̅χ→ e^+ e^- and a continuum piece from the decays of unstable particles:

    dN/dE_e^+(E_e^+)     = Br(χ̅χ→ e^+ e^-)  δ(E_CM/2 - E_e^+) + . dN/dE_e^+|_χ̅χ,dec(E_e^+)

For the scalar and vector model, the decay piece receives contributions from the μ^+μ^- and π^+π^- final states.

The muon decays nearly 100% of the time through μ^+→ e^+ν_eν̅_μ. In the muon's rest frame, the positron's energy spectrum is <cit.> 


    dN/dE_e^+    = -4 √(x^2-4 r^2)(r^2 (4-3 x)+x (2 x-3))/m_μ

where r = m_e / m_μ and x = 2E_e^+ / m_μ. The spectrum can be boosted to other frames using Eqn. (<ref>). It is accessible in hazma using:
python
>>> from hazma.positron_spectra import muon as dnde_p_mu
>>> e_mu = 150.                      # muon energy
>>> e_p = np.array([1., 10., 100.])  # positron energies
>>> dnde_p_mu(e_p, e_mu)
array([4.86031362e-05, 4.56232320e-03, 4.45753994e-03])


The charged pion primarily decays through π^+ →μ^+ ν_μ, which yields a positron when the muon decays. The energy spectrum is thus obtained by boosting the muon spectrum from Eqn. (<ref>). An additional contribution comes from the monochromatic positron produced through π^+ → e^+ ν_e, which has a helicity suppressed branching fraction Br(π^+ → e^+ ν_e) = 1.23-4. The total spectrum can be computed in hazma with
python
>>> from hazma.positron_spectra import charged_pion as dnde_p_pi
>>> e_pi = 150.                      # charged pion energy
>>> e_p = np.array([1., 10., 100.])  # positron energies
>>> dnde_p_pi(e_p, e_pi)
array([3.84163631e-05, 3.85242442e-03, 2.55578895e-05])






Fig. (<ref>) exhibits representative positron annihilation spectra from the scalar model with Higgs portal couplings and vector model with kinetic mixing couplings (first and second rows), and several DM masses (different columns). Changing the couplings to prohibit annihilation into leptons has a straightforward impact on the total spectrum. Here annihilation into mediators is kinematically forbidden. Positron spectra for annihilation into mediators are shown for these models in Fig. (<ref>). When the mediator is light the spectrum is box-shaped since the mediator can only produce electrons by decaying into them directly. At higher energies the spectrum is a superposition of the rounded spectra from mediator decays into muons and pions and this box spectrum. For the Higgs portal model, the box is subdominant to the rounded part of the spectrum since the mediator decay width into e^+ e^- is Yukawa-suppressed. The top edge of the box is barely visible for the orange curve near 200 MeV in the top middle panel. For the vector model, the box component of the spectrum is visible as a sawtooth edge in the lower middle and right panels.

The interface for computing positron spectra mirrors the gamma-ray spectrum interface. The total_positron_spectrum and positron_spectra methods in Theory give the total positron spectrum and individual channels' contributions (scaled by branching fraction). The positron_lines method gives information about final states producing monochromatic positrons:
python
>>> from hazma.scalar_mediator import HiggsPortal
>>> e_cm = 305.                       # DM center of mass energy
>>> e_ps = np.array([1., 10., 100.])  # positron energies
>>> hp = HiggsPortal(mx=150., ms=1e3, gsxx=0.7, stheta=0.1)
>>> hp.total_positron_spectrum(e_ps, e_cm)
array([2.60677406e-05, 2.58192085e-03, 4.09383294e-04])
>>> hp.positron_spectra(e_ps, e_cm)
'mu mu': array([6.10588051e-06, 5.70145391e-04, 5.78482356e-04]),
 'pi pi': array([2.14034924e-05, 2.13725359e-03, 9.48371777e-05]),
 's s': array([0., 0., 0.]),
 'total': array([2.75093729e-05, 2.70739898e-03, 6.73319534e-04])
>>> hp.positron_lines(e_cm)
'e e': 'energy': 152.5, 'bf': 7.711433862697906e-06

For the built-in models, functions computing the positron spectra for annihilations into particular final states follow a similar naming scheme to the gamma-ray ones, and can be called directly.

A method (more precisely, an InterpolatedUnivariateSpline) to compute the total positron spectrum convolved with an energy resolution function is obtained using:
python
>>> e_cm = 305.                  # DM center of mass energy
>>> e_p_min, e_p_max = 1., 100.  # define energy range
>>> energy_res = lambda e: 0.05
>>> dnde_p_conv = hp.total_conv_positron_spectrum_fn(
...     e_p_min, e_p_max, e_cm, energy_res)
>>> dnde_p_conv(20.)
0.00864851382906508
>>> dnde_p_conv.integral(10, 100)  # integrate spectrum
0.6538810882108401




§ GAMMA RAY LIMITS



In hazma we include two routines to set constraints on the DM self-annihilation cross section  using current gamma-ray flux measurements and to make projections for planned detectors.

As is well-known, the primary gamma-ray flux from DM annihilating in a region of the sky subtending a solid angle ΔΩ is

    dΦ/dE_γ(E)    = ΔΩ/4π·/2 f_χ m_χ^2· J ·dN/dE_γ(E).

The factor f_χ is 1 if the DM is self-conjugate and 2 otherwise (2 for the models currently included with hazma). The J factor counts the number of pairs of DM particles in the observing region:

    J ≡1/ΔΩ∫_ΔΩ dΩ ∫_L̆ŎS̆ dl ρ(r(l, ψ))^2,

where ρ(r) is the DM density profile. The annihilation spectrum dN/dE_γ|_χ̅χ was computed in detail in the previous section. A detector with finite energy resolution observes the smoothed flux

    dΦ_ϵ/dE_γ(E)    ≡∫ dE'  R_ϵ(E | E') dΦ/dE_γ(E'),

where the spectral resolution function is the probability that a gamma ray with true energy E' is reconstructed with energy E, and can be approximated as a Gaussian <cit.>:

    R_ϵ(E|E')    ≡1/√(2π) ϵ(E')  E'exp[ -(E - E')^2/2 (ϵ(E')  E')^2].


In the limit-setting procedure for existing experiments, the integral of the smoothed spectrum in each energy bin is required not to exceed the observed value plus twice the upper error bar: Φ_ϵ^(i) |_χ̅χ < Φ^(i)|_ŏb̆s̆ + 2 σ^(i). While better limits can be obtained by assuming a background model, this introduces significant systematic uncertainties since the model is derived from precisely the observations under consideration and would thus only marginally improve the limits on  <cit.>.

For assessing the discovery reach by upcoming experiments, hazma uses an unbinned procedure. A detector can be characterized by an effective area A_ĕf̆f̆(E) and an (energy-independent) observing time T_ŏb̆s̆. Over the energy range [E_m̆ĭn̆, E_m̆ăx̆], the number of photons observed is then

    N_γ   = ∫_E_m̆ĭn̆^E_m̆ăx̆ dE  T_ŏb̆s̆  A_ĕf̆f̆(E)  dΦ_ϵ/dE_γ(E).

Assuming Poisson statistics, the minimum-detectable  is obtained by requiring the signal-to-noise ratio to be significant at the n_σ = 5 level: N_γ|_χ̅χ / √(N_γ|_b̆ğ) = n_σ. Holding the DM mass fixed, the lowest-possible value of  is found by optimizing over E_m̆ĭn̆ and E_m̆ăx̆ to maximize N_γ|_χ̅χ / √(N_γ|_b̆ğ). This energy range does not completely align with spectra features. For example, in the case of the π^0 spectrum, setting E_m̆ĭn̆ equal to the left boundary of the “box” would include too much of the dΦ/dE_γ|_b̆ğ∝ E_γ^-2 spectrum.

The unbinned analysis requires a background model, two of which are included with hazma. The simplest was obtained by fitting a power law to COMPTEL and EGRET data from high galactic latitudes over the energy range 0.8MeV - 10GeV <cit.>. This is suitable for projecting limits on  from dwarf galaxy observations. A model derived using GALPROP to fit existing MeV gamma-ray data in the region |b| < 10^∘, |l| < 30^∘ is also included to assess the sensitivity to emission from the Galactic Center.





These two limit-setting procedures can currently be applied in hazma to data from three existing experiments (COMPTEL, EGRET and Fermi-LAT) and one proposed detector (e-ASTROGAM). The energy ranges, approximate effective areas and energy resolutions for these are supplied in Tab. (<ref>). Each gamma-ray data set is paired with a target region, defined by the galactic coordinate ranges, solid angles and J-factors (assuming an NFW distribution for the galactic DM halo) in Tab. (<ref>). These are packaged into FluxMeasurements objects in hazma.

The limit-setting procedures are methods in the Theory class. For example, the following snippet imports a FluxMeasurement object containing information about EGRET observations and computes the  limits for the Higgs portal model at three DM masses in c̆m̆^3/s̆:
python
>>> from hazma.scalar_mediator import ScalarMediator
>>> from hazma.gamma_ray_parameters import egret_diffuse, TargetParams
>>> egret_diffuse.target  # target region information
TargetParams(J=3.79e+27, dOmega=6.584844306798711)
>>> sm = ScalarMediator(mx=150., ms=1e4, gsxx=1., gsff=0.1, gsGG=0.5,
...                     gsFF=0, lam=1e5)
>>> sm.binned_limit(egret_diffuse)
7.890021123433229e-27  # cm^3 / s

It is important to note that hazma only uses the model parameters (m_S, g_Sχχ, g_Sff, g_SGG, g_SFF and Λ) to compute the shape of the spectrum and treats  as an independent parameter (c.f. <ref>). To relate the model parameters to , the user can employ the method Theory.annihilation_cross_sections(). The FluxMeasurement objects for COMPTEL and Fermi-LAT observations are named comptel_diffuse and fermi_diffuse and can be imported similarly. App. (<ref>) describes how to create new flux measurement objects.

Computing unbinned limits requires specifying several quantities characterizing the behavior of a gamma-ray detector:

    
  * An effective area function returning A_ĕf̆f̆(E) in c̆m̆^2;
    
  * An energy resolution function returning ϵ(E') in %;
    
  * An observing time in seconds;
    
  * An object of type TargetParams containing the J-factor and ΔΩ for the target region;
    
  * A BackgroundModel object, which has a method dPhi_dEdOmega() returning the differential gamma-ray flux flux per solid angle for the background model, and an attribute e_range specifying the range of E_γ values for which the model is valid.

Using the same instance of scalar_mediator and DM masses as above, we can compute the projected limits for e-ASTROGAM using the 10^∘×10^∘ window around the Galactic Center from Tab. (<ref>) as the target:
python
>>> from hazma.gamma_ray_parameters import A_eff_e_astrogam
>>> from hazma.gamma_ray_parameters import energy_res_e_astrogam
>>> from hazma.gamma_ray_parameters import T_obs_e_astrogam
>>> from hazma.gamma_ray_parameters import gc_target, gc_bg_model
>>> gc_target.J                               # target J-factor
1.795e+29
>>> gc_target.dOmega                          # target extent
0.12122929761504078
>>> sm.unbinned_limit(A_eff_e_astrogam,       # effective area
...                   energy_res_e_astrogam,  # energy resolution
...                   T_obs_e_astrogam,       # observing time
...                   gc_target,              # target region
...                   gc_bg_model)            # background model
5.6701062876845636e-30  # cm^3 / s

Modifying these ingredients to study other planned gamma-ray detectors is described in App. (<ref>.

We do not currently use secondary radiation to set limits on our models. Such limits are highly dependent on (for example) the choice of observational target, the magnetic field intensity and spatial structure <cit.>; Ref. (<cit.>) calculates the secondary emission for select MeV dark matter annihilation final states. Since hazma can already compute the relevant input spectra for determining secondary spectra, we plan to incorporate this into future versions of the code.



§ COSMIC MICROWAVE BACKGROUND LIMITS



Dark matter annihilations at recombination inject ionizing particles into the photon-baryon plasma. The resulting changes in the residual ionization fraction and baryon temperature modify the CMB temperature and polarization power spectra, particularly at small scales, which puts a constraint on , the DM self-annihilation cross section at recombination <cit.>. The size of the effect of the changes depend on the energy per unit volume per unit time, which is related to the effective parameter

    p_ăn̆n̆≡ /.

The current 95% upper limit on p_ăn̆n̆ from Planck is 3.5-31cm^3s^-1MeV^-1 <cit.>.  is the fraction of energy per DM annihilation imparted to the plasma. This depends on how many electrons/positrons and photons are produced per DM annihilation (ie, the spectra computed earlier in this work) and how efficiently these particles re-ionize the CMB (quantified by the factors f_ĕf̆f̆^γ and f_ĕf̆f̆^e^+ e^-) <cit.>:

    f_ĕf̆f̆^χ = 1/E_C̆M̆∫_0^E_C̆M̆ / 2 dE ( 2 f_ĕf̆f̆^e^+e^- E  dN/dE_e^+ + f_ĕf̆f̆^γ E  dN/dE_γ).


If the DM self-annihilation cross section is velocity dependent (as is the case for the scalar mediator model, where the annihilation is p-wave), the DM velocity v_χ,C̆M̆B̆ at recombination is required to relate  to the present-day value of  in indirect detection targets. This velocity is also required for computing the center of mass energy in DM self-annihilation events at recombination. The velocity depends on the kinetic decoupling temperature, T_K̆D̆. Before kinetic decoupling DM initially is coupled to the plasma so that v_χ,C̆M̆B̆ = √(3 T_γ / ), but afterwards it cools more quickly (T_χ∝ z^2). The DM velocity is thus <cit.>

    v_χ,C̆M̆B̆ = √(3 T_χ/)≈ 2-4( T_γ/1eV) ( 1MeV/) ( 10^-4/x_k̆d̆)^1/2,

where x_k̆d̆≡ T_K̆D̆ / generally lies in the range 10^-6 - 10^-4. Since this velocity is much smaller than the characteristic velocity in the Milky Way, the CMB constraints on  for models with p-wave suppression are in general much weaker than the gamma-ray ones. The CMB constraints in the s-wave case are generally much stronger, but can be evaded in models where the DM is partially produced through the decay of a heavier dark-sector degree of freedom after recombination <cit.>.



The quantity  as well as the separate e^± and γ terms in Eqn. (<ref>) can be computed for a given kinetic decoupling temperature:
python
>>> from hazma.scalar_mediator import ScalarMediator
>>> sm = ScalarMediator(mx=150., ms=1e4, gsxx=1., gsff=0.1, gsGG=0.5,
...                     gsFF=0, lam=1e5)
>>> x_kd = 1e-4
>>> sm.f_eff(x_kd), sm.f_eff_ep(x_kd), sm.f_eff_g(x_kd)
(0.4364732027954761, 0.16674887821288106, 0.26972432458259504)

Fig. (<ref>) shows this factor as a function of DM and mediator mass for the scalar and vector models with the couplings we have been using as benchmarks. There are sharp variations in  arising from SM or two-mediator final states becoming accessible as the DM mass changes. Deriving the kinetic decoupling temperature for particular models rather than leaving it as a free parameter is left for future work (or the user).

The Theory.cmb_limit() function computes the constraint on  given a value of x_k̆d̆ and upper limit on p_ăn̆n̆:
python
>>> from hazma.cmb import p_ann_planck_temp_pol as p_ann
>>> p_ann
3.5e-31  # cm^3 / MeV / s
>>> x_kd = 1e-4
>>> sm = HiggsPortal(mx=150., ms=1e4, gsxx=1., stheta=0.01)
>>> sm.cmb_limit(x_kd, p_ann)
1.2398330863024796e-28  # cm^3 / s

As with the gamma-ray limits described in the previous section, the parameters of the theory (here m_S, g_Sχχ and s_θ) determine the shapes of the photon and e^± spectra, and Theory.cmb_limit() treats  as a free parameter. The user can find which model parameters reproduce this cross section using Theory.annihilation_cross_sections().



§ CONCLUSION



We have introduced hazma, a code for the calculation of gamma-ray and electron-positron spectra as well as constraints and projected reach for future telescopes for sub-GeV dark matter models. hazma is timely both because of forthcoming new observational capabilities at sub-GeV gamma-ray energies, and because of renewed interest in dark matter models beyond the WIMP paradigm and at lower-than-usually-considered masses. In future work we will use hazma to study the complementarity between indirect and other probes of sub-GeV dark matter.

hazma allows users to compute both spectra of gamma rays, electrons and positrons for individual annihilation final states containing arbitrarily many particles, as well as for user-implemented or a set of built-in dark matter models. In the latter case, the code computes the relevant branching fractions and the ensuing decay chains leading to photons, electrons and positrons, and calculates model-dependent emission processes such as e.g. internal bremsstrahlung. The interactions between the mediator particles in these models and the final state mesons is implemented using chiral perturbation theory techniques.

hazma contains several updates and refinements over how gamma-ray spectra have been calculated in the recent literature. These include several radiative decay spectra, such as from charged pion decay or from muon decay, final state radiation calculations going beyond the Altarelli-Parisi splitting function approximation for the built-in models, and the ability to include arbitrarily many mesons and leptons in the final states.

We presented numerous snippets of code, and all code used to produce the figures shown herein manuscript is available through clickable links. We also showed a few advanced usage features, such as adding new gamma-ray experiments for deriving constraints from existing data or assessing projected reach for planned detectors, and implementing user-defined models.

LM and SP are partly supported by the U.S. Department of Energy grant number de-sc0010107. We acknowledge with gratitude early collaboration on this project with Francesco D'Eramo. Francesco D'Eramo was the first to put forward the idea of using chiral perturbation theory in the context of accurately calculating the annihilation and decay products of MeV dark matter. He was also responsible for several of the early calculations, part of the writing of the text, and for overall educating us all on numerous topics in chiral perturbation theory. We thank Michael Peskin for important input and feedback, and Regina Caputo for discussions about MeV gamma-ray observatories. We thank Graham White as well as Nicholas Rodd for pointing out typos in an earlier version of the manuscript. Finally, we thank Eulogio Oset and Jose Antonio Oller for clarifying some points about unitarized chiral perturbation theory.





§ INSTALLATION



hazma was developed for python3. Before installing hazma, the user needs to install several well-established python packages: cython, scipy, numpy, and scikit-image. Theses are easily installed by using PyPI. If the user has PyPI installed on their system, then these packages can be installed using
bash
pip install cython, scipy, numpy, scikit-image, matplotlib

hazma can be installed in the same way, using:
bash
pip install hazma

This will download a tarball from the PyPI repository, compile all the c-code and install hazma on the system. Alternatively, the user can install hazma by downloading the package from <https://github.com/LoganAMorrison/Hazma.git>. Once downloaded, navigate to the package directory using the command line and run either
bash
pip install .

or
bash
python setup.py install


Note that since hazma makes extensive usage of the package cython, the user will need to have a c and c++ compiler installed on their system (for example gcc and g++ on unix-like systems or Microsoft Visual Studios 2015 or later on Windows). For more information, see the cython installation guide: <https://cython.readthedocs.io/en/latest/src/quickstart/install.html>



§ BASIC USAGE





 §.§ Using the Built-In Models



Throughout the text, we gave usage examples of how to use the models built into hazma. Here we give a compact review of how to use these models in hazma. All the models built into hazma have identical interfaces. The only difference in the interfaces is the parameters which need to be specified for the particular model being used. Thus, we will only show the usage of one of the models. The others can be used in an identical fashion. The example we will use is the KineticMixing model. To create a KineticMixing model, we use:
python
# Import the model
>>> from hazma.vector_mediator import KineticMixing
# Specify the parameters
>>> params = 'mx': 250.0, 'mv': 1e6, 'gvxx': 1.0, 'eps': 1e-3
# Create KineticMixing object
>>> km = KineticMixing(**params)

Here we have created a model with a dark matter fermion of mass 250 MeV, a vector mediator which mixes with the standard model with a mass of 1 TeV. We set the coupling of the vector mediator to the dark matter, g_Vχ = 1 and set the kinetic mixing parameter ϵ = 10^-3. To list all of the available final states for which the dark matter can annihilate into, we use:
python
>>> km.list_annihilation_final_states()
['mu mu', 'e e', 'pi pi', 'pi0 g', 'pi0 v', 'v v']

This tells us that we can potentially annihilate through: χχ̅→μ^+μ^-,e^+e^-,π^+π^-,π^0γ, π^0 V or V^μV^μ. However, which of these final states is actually available depends on the center of mass energy. We can see this fact by looking at the annihilation cross sections or branching fractions, which can be computed using:
python
>>> cme = 2.0 * km.mx * (1.0 + 0.5 * 1e-6)
>>> km.annihilation_cross_sections(cme)
'mu mu': 8.998735387276086e-25,
 'e e': 9.115024836416874e-25,
 'pi pi': 1.3013263970304e-25,
 'pi0 g': 1.7451483984993156e-29,
 'pi0 v': 0.0,
 'v v': 0.0,
 'total': 1.941526113556321e-24
>>> km.annihilation_branching_fractions(cme)
'mu mu': 0.4634877339245762,
 'e e': 0.46947732367713324,
 'pi pi': 0.06702595385888176,
 'pi0 g': 8.988539408840104e-06,
 'pi0 v': 0.0,
 'v v': 0.0

Here we have chosen a realistic center of mass energy for dark matter in our galaxy, which as a velocity dispersion of δ v∼ 10^-3. We can that the V^μV^μ final state is unavailable, as it should be since the vector mediator mass is too heavy. In this theory, the vector mediator can decay. If we would like to know the decay width and the partial widths, we can use:
python
>>> km.partial_widths()
'pi pi': 0.0006080948890354345,
 'pi0 g': 0.23374917012731816,
 'x x': 26525.823848648604,
 'e e': 0.0024323798404358825,
 'mu mu': 0.0024323798404358803,
 'total': 26526.0630706733

If we would like to know the gamma-ray spectrum from dark matter annihilations, we can use 
python
# Specify photon energies to compute spectrum at
>>> photon_energies = np.array([cme / 4])
>>> km.spectra(photon_energies, cme)
'mu mu': array([2.8227189e-05]),
 'e e': array([0.00013172]),
 'pi pi': array([2.1464018e-06]),
 'pi0 g': array([7.66452627e-08]),
 'pi0 v': array([0.]),
 'v v': array([0.]),
 'total': array([0.00016217])

Note that we only used a single photon energy because of display purposes, but in general the user can specify any number of photon energies. If the user would like access to the underlying spectrum functions, they can use:
python
>>> spec_funs = km.spectrum_functions()
>>> spec_funs['mu mu'](photon_energies, cme)
array([6.09016958e-05])
>>> mumu_bf = km.annihilation_branching_fractions(cme)['mu mu']
>>> mumu_bf * spec_funs['mu mu'](photon_energies, cme)
array([2.8227189e-05])

Notice that the direct call to the spectrum function for χχ̅→μ^+μ^- doesn't given the same result as km.spectra(photon_energies, cme)['mu mu']. This is because the branching fractions are not applied for the spec_funs = km.spectrum_functions(). If the user doesn't care about the underlying components of the gamma-ray spectra, the can simply call
python
>>> km.total_spectrum(photon_energies, cme)
array([0.00016217])

to get the total gamma-ray spectrum. The reader may have caught the fact that there is a gamma-ray line in the spectrum for χχ̅→π^0γ. To get the location of this line, the user can use:
python
>>> km.gamma_ray_lines(cme)
'pi0 g': 'energy': 231.78145156177675, 'bf': 8.988539408840104e-06

This tells us the process which produces the line, the location of the line and the branching fraction for the process. We do not include the line in the total spectrum since the line produces a Dirac-delta function. In order to get a realistic spectrum including the line, we need to convolve the gamma-ray spectrum with an energy resolution. This can be achieved using:
python
>>> min_photon_energy = 1e-3
>>> max_photon_energy = cme
>>> energy_resolution = lambda photon_energy : 1.0
>>> number_points = 1000
>>> spec = km.total_conv_spectrum_fn(min_photon_energy, max_photon_energy, 
...                                  cme, energy_resolution, number_points)
>>> spec(cme / 4)  # compute the spectrum at a photon energy of `cme/4`
array(0.00167605)

The km.total_conv_spectrum_fn computes and returns an interpolating function of the convolved function. An important thing to note here is that the km.total_conv_spectrum_fn takes in a function for the energy resolution. This allows the user to define the energy resolution to depend on the specific photon energy, which is often the case for gamma-ray telescopes.

Next we present the positron spectra. These have an identical interface to the gamma-ray spectra, so we only show how to call the functions and we suppress the output
python
>>> from hazma.parameters import electron_mass as me
>>> positron_energies = np.logspace(np.log10(me), np.log10(cme), num=100)
>>> km.positron_spectra(positron_energies, cme)
>>> km.positron_lines(cme)
>>> km.total_positron_spectrum(positron_energies, cme)
>>> dnde_pos = km.total_conv_positron_spectrum_fn(min(positron_energies),
...                                           max(positron_energies),
...                                           cme, 
...                                           energy_resolution,
...                                           number_points)

The last thing that we would like to demonstrate is how to compute limits. In order to compute the limits on the annihilation cross section of a model from a gamma-ray telescope, say EGRET, we can use:
python
>>> from hazma.gamma_ray_parameters import egret_diffuse
# Choose DM masses from half the electron mass to 250 MeV
>>> mxs = np.linspace(me/2., 250., num=100)
# Compute limits from e-ASTROGAM
>>> limits = np.zeros(len(mxs), dtype=float)
>>> for i, mx in enumerate(mxs):
...     km.mx = mx
...     limits[i] = km.binned_limit(egret_diffuse)

Similarly, if we would like to set constraints using e-ASTROGAM, one can use:
python
>>> # Import target and background model for the e-ASTROGAM telescope
>>> from hazma.gamma_ray_parameters import gc_target, gc_bg_model
>>> # Import telescope parameters
>>> from hazma.gamma_ray_parameters import A_eff_e_astrogam
>>> from hazma.gamma_ray_parameters import energy_res_e_astrogam
>>> from hazma.gamma_ray_parameters import T_obs_e_astrogam
>>> # Choose DM masses from half the electron mass to 250 MeV
>>> mxs = np.linspace(me/2., 250., num=100)
>>> # Compute limits
>>> limits = np.zeros(len(mxs), dtype=float)
>>> for i, mx in enumerate(mxs):
...     km.mx = mx
...     limits[i] = km.unbinned_limit(A_eff_e_astrogam,
...                                   energy_res_e_astrogam,
...                                   T_obs_e_astrogam,
...                                   gc_target,
...                                   gc_bg_model)




 §.§ Subclassing the Built-In Models



As mentioned in Sec. (<ref>), the user may not be interested in the generic models built into hazma, but instead a more specialized model. In this case, it makes sense for the user to subclass one of the models (i.e. create a class which inherits from one of the models.) As and example, we illustrate how to do this with the Higgs-portal model (of course this model is already built into hazma, but it works nicely as an example.) Recall that the full set of parameters for the scalar mediator model are:

    
  * m_χ: dark matter mass,
    
  * m_S: scalar mediator mass,
    
  * g_Sχ: coupling of scalar mediator to dark matter,
    
  * g_Sf: coupling of scalar mediator to standard model fermions,
    
  * g_SG: effective coupling of scalar mediator to gluons,
    
  * g_SF: effective coupling of scalar mediator to photons and
    
  * Λ: cut-off scale for the effective interactions.

In the case of the Higgs-portal model, the scalar mediator talks to the standard model only through the Higgs boson, i.e. it mixes with the Higgs. Therefore, the scalar mediator inherits its interactions with the standard model fermions, gluons and photon through the Higgs. In Sec. (<ref>), we gave the relationships between the couplings g_Sf, g_SG, g_SF, Λ and sinθ, v_h, where sinθ is the mixing angle for the scalar and Higgs and v_h is the Higgs vacuum expectation value. These relationships are:

    g_Sf   = sinθ,     g_SG   = 3sinθ,     g_SF    = -5/6sinθ,    Λ   = v_h.

Therefore, the Higgs-portal model doesn't need to know about all of the couplings of the generic scalar mediator model. It only needs to know about m_χ, m_S, g_Sχ and sinθ. Below, we construct a class which subclasses the scalar mediator class to implement the Higgs-portal model.
python
from hazma.scalar_mediator import ScalarMediator
from hazma.parameters import vh

class HiggsPortal(ScalarMediator):
    def __init__(self, mx, ms, gsxx, stheta):
        self._lam = vh
        self._stheta = stheta
        super(HiggsPortal, self).__init__(mx, ms, gsxx, stheta, 3.*stheta,
                                          -5.*stheta/6., vh)

    @property
    def stheta(self):
        return self._stheta

    @stheta.setter
    def stheta(self, stheta):
        self._stheta = stheta
        self.gsff = stheta
        self.gsGG = 3. * stheta
        self.gsFF = - 5. * stheta / 6.

    # Hide underlying properties' setters
    @ScalarMediator.gsff.setter
    def gsff(self, gsff):
        raise AttributeError("Cannot set gsff")

    @ScalarMediator.gsGG.setter
    def gsGG(self, gsGG):
        raise AttributeError("Cannot set gsGG")

    @ScalarMediator.gsFF.setter
    def gsFF(self, gsFF):
        raise AttributeError("Cannot set gsFF")

There are a couple things to note about our above implementation. First, our model only takes in m_χ, m_S, g_Sχ and sinθ, as desired. But the underlying model, i.e. the ScalarMediator model only knows about m_χ, m_S, g_Sχ, g_Sf, g_SG, g_SF and Λ. So if we update sinθ, we additionally need to update the underlying parameters, g_Sf, g_SG, g_SF and Λ. The easiest way to do this is using getters and setters by defining sinθ to be a property through the @property decorator. Then every time we update sinθ, we can also update the underlying parameters. The second thing to note is that we want to make sure we don't accidentally change the underlying parameters directly, since in this model, they are only defined through sinθ. We an ensure that we cannot change the underlying parameters directly by overriding the getters and setters for gsff, gsGG and gsGG and raising an error if we try to change them. This isn't strictly necessary (as long as the user is careful), but can help avoid confusing behavior. 



§ ADVANCED USAGE





 §.§ Adding New Gamma-Ray Experiments



Currently hazma only includes information for producing projected unbinned limits with e-ASTROGAM, using the dwarf Draco or inner 10^∘×10^∘ region of the Milky Way as a target. Adding new detectors and target regions is straightforward. As described in Sec. (<ref>), a detector is characterized by the effective area A_ĕf̆f̆(E), the energy resolution ϵ(E) and observation time T_ŏb̆s̆. In hazma, the first two can be any callables (functions) and the third must be a float. The region of interest is defined by a TargetParams object, which can be instantiated with:
python
>>> from hazma.gamma_ray_parameters import TargetParams
>>> tp = TargetParams(J=1e29, dOmega=0.1)

The background model should be packaged in an object of type BackgroundModel. This light-weight class has a function dPhi_dEdOmega() for computing the differential photon flux per solid angle (in M̆ĕV̆^-1sr) and an attribute e_range specifying the energy range over which the model is valid (in MeV). New background models are defined by passing these two the initializer:
python
>>> from hazma.background_model import BackgroundModel
>>> bg = BackgroundModel(e_range=[0.5, 1e4],
...                      dPhi_dEdOmega=lambda e: 2.7e-3 / e**2)


As explained in Sec. (<ref>), gamma-ray observation information from Fermi-LAT, EGRET and COMPTEL is included with hazma, and other observations can be added using the container class FluxMeasurement. The initializer requires:

    
  * The name of a CSV file containing gamma-ray observations. The file's columns must contain:
    
        
  * Lower bin edge (MeV)
        
  * Upper bin edge (MeV)
        
  * E^n d^2Φ / dE  dΩ (in M̆ĕV̆^n-1cm^-2s^-1sr^-1)
        
  * Upper error bar (in M̆ĕV̆^n-1cm^-2s^-1sr^-1)
        
  * Lower error bar (in M̆ĕV̆^n-1cm^-2s^-1sr^-1)
    
    Note that the error bar values are their y-coordinates, not their relative distances from the central flux.
    
  * The detector's energy resolution function.
    
  * A TargetParams object for the target region.

For example, a CSV file obs.csv containing observations
text
150.0,275.0,0.0040,0.0043,0.0038
650.0,900.0,0.0035,0.0043,0.003

with n=2 for an instrument with energy resolution ϵ(E) = 0.05 observing the target region tp defined above can be loaded using[If the CSV containing the observations uses a different power of E than n=2, this can be specified using the power keyword argument to the initializer for FluxMeasurement.]
python
>>> from hazma.flux_measurement import FluxMeasurement
>>> obs = FluxMeasurement("obs.dat", lambda e: 0.05, tp)

The attributes of the FluxMeasurement store all of the provide information, with the E^n prefactor removed from the flux and error bars, and the errors converted from the positions of the error bars to their sizes. These are used internally by the Theory.binned_limit() method, and can be accessed as follows:
python
>>> obs.e_lows, obs.e_highs
(array([150., 650.]), array([275., 900.]))
>>> obs.target
<hazma.gamma_ray_parameters.TargetParams at 0x1c1bbbafd0>
>>> obs.fluxes
array([8.85813149e-08, 5.82726327e-09])
>>> obs.upper_errors
array([6.64359862e-09, 1.33194589e-09])
>>> obs.lower_errors
array([4.42906574e-09, 8.32466181e-10])
>>> obs.energy_res(10.)
0.05




 §.§ User-Defined Models



In this subsection, we demonstrate how to implement new models in Hazma. A notebook containing all the code in this appendix can be downloaded from GitHub https://github.com/LoganAMorrison/Hazma/blob/b525ac3482ce355f4109c9ad6622470c14d0931f/notebooks/hazma_paper/hazma_example.ipynb. The model we will consider is an effective field theory with a Dirac fermion DM particle which talks to neutral and charged pions through gauge-invariant dimension-5 operators. The Lagrangian for this model is:

    ℒ   ⊃c_1Λχχπ^+π^- + c_2Λχχπ^0π^0

where c_1, c_2 are dimensionless Wilson coefficients and Λ is the cut-off scale of the theory. In order to implement this model in Hazma, we need to compute the annihilation cross sections and the FSR spectra. The annihilation channels for this model are simply χ̅χ→π^0π^0 and χ̅χ→π^+π^-. The computations for the cross sections are straight forward and yield:

    σ(χ̅χ→π^+π^-)    = c_1^2 √(1-4 μ _π^2)√(1-4 μ _χ^2)/32 πΛ^2
    σ(χ̅χ→π^0π^0)    = c_2^2 √(1-4 μ_π^0^2)√(1-4 μ_χ^2)/8 πΛ^2

where Q is the center of mass energy, μ_χ = m_χ/Q, μ_π = m_π^±/Q and μ_π^0 = m_π^0/Q. In addition to the cross sections, we need the FSR spectrum for χχ→π^+π^-γ. This is:

    dN(χ̅χ→π^+π^-γ)dE_γ   = α(2 f(x)-2(1-x-2 μ_π ^2)
       log(1-x-f(x)/1-x+f(x)))/π√(1-4 μ_π ^2) x

where 

    f(x)    = √(1-x)√(1-x-4 μ_π ^2)

We are now ready to set up the Hazma model. For hazma to work properly, we will need to define the following functions in our model:


    
  * annihilation_cross_section_funcs(): A function returning a dictionary of the annihilation cross sections functions, each of which take a center of mass energy.
    
  * spectrum_funcs(): A function returning a dictionary of functions which take photon energies and a center of mass energy and return the gamma-ray spectrum contribution from each final state.
    
  * gamma_ray_lines(e_cm): A function returning a dictionary of the gamma-ray lines for a given center of mass energy.
    
  * positron_spectrum_funcs(): Like spectrum_funcs(), but for positron spectra.
    
  * positron_lines(e_cm): A function returning a dictionary of the electron/positron lines for a center of mass energy.

In the interests of compartmentalization, we find it easiest place each of these components in its own class (a mixin in python terminology) which are then combined into a master class representing our model. Before we begin writing the classes, we need to import a few helper functions and constants from hazma:
python
import numpy as np # NumPy is heavily used
import matplotlib.pyplot as plt # Plotting utilities
# neutral and charged pion masses
from hazma.parameters import neutral_pion_mass as mpi0
from hazma.parameters import charged_pion_mass as mpi
from hazma.parameters import qe # Electric charge
# Positron spectra for neutral and charged pions
from hazma.positron_spectra import charged_pion as pspec_charged_pion
# Deay spectra for neutral and charged pions
from hazma.decay import neutral_pion, charged_pion
# The `Theory` class which we will ultimately inherit from
from hazma.theory import Theory

Now, we implement a cross section class:
python
class HazmaExampleCrossSection:
    def sigma_xx_to_pipi(self, Q):
        mupi = mpi / Q
        mux = self.mx / Q

        if Q > 2 * self.mx and Q > 2 * mpi:
            sigma = (self.c1**2 * np.sqrt(1 - 4 * mupi**2) * 
                     np.sqrt(1 - 4 * mux**2)**2 / 
                     (32.0 * self.lam**2 * np.pi))
        else:
            sigma = 0.0

        return sigma

    def sigma_xx_to_pi0pi0(self, Q):
        mupi0 = mpi0 / Q
        mux = self.mx / Q

        if Q > 2 * self.mx and Q > 2 * mpi0:
            sigma = (self.c2**2 * np.sqrt(1 - 4 * mux**2) *
                     np.sqrt(1 - 4 * mupi0**2) /
                     (8.0 * self.lam**2 * np.pi))
        else:
            sigma = 0.0

        return sigma

    def annihilation_cross_section_funcs(self):
        return 'pi0 pi0': self.sigma_xx_to_pi0pi0,
                'pi pi': self.sigma_xx_to_pipi

The key function is annihilation_cross_section_funcs, which is required to be implemented. Next, we implement the spectrum functions which will produce the FSR and decay spectra:
python
class HazmaExampleSpectra:
    def dnde_pi0pi0(self, e_gams, e_cm):
        # Decay spectra for the neutral pions
        return 2.0 * neutral_pion(e_gams, e_cm / 2.0)
    
    def __dnde_xx_to_pipig(self, e_gam, Q):
        # Unvectorized function for computing FSR spectrum
        mupi = mpi / Q
        mux = self.mx / Q
        x = 2.0 * e_gam / Q
        if 0.0 < x and x < 1. - 4. * mupi**2:
            dnde = (
                qe ** 2
                * (
                    2 * np.sqrt(1 - x) * np.sqrt(1 - 4 * mupi ** 2 - x)
                    + (-1 + 2 * mupi ** 2 + x)
                    * np.log(
                        (-1 + np.sqrt(1 - x) *
                         np.sqrt(1 - 4 * mupi ** 2 - x) + x) ** 2
                        / (1 + np.sqrt(1 - x) *
                           np.sqrt(1 - 4 * mupi ** 2 - x) - x) ** 2
                    )
                )
            ) / (Q * 2.0 * np.sqrt(1 - 4 * mupi ** 2) * np.pi ** 2 * x)
        else:
            dnde = 0
        
        return dnde

    def dnde_pipi(self, e_gams, e_cm):
        # Sum the FSR and decay spectra for the charged pions
        return (np.vectorize(self.__dnde_xx_to_pipig)(e_gams, e_cm) +
                2. * charged_pion(e_gams, e_cm / 2.0))
    
    def spectrum_funcs(self):
        return 'pi0 pi0':  self.dnde_pi0pi0,
                'pi pi':  self.dnde_pipi

    def gamma_ray_lines(self, e_cm):
        return 

The required functions here are spetrum_funcs and gamma_ray_lines. Note the the second __dnde_xx_to_pipig is an unvectorized helper function, which is not to be used directly. Next we implement the positron spectra:
python
class HazmaExamplePositronSpectra:
    def dnde_pos_pipi(self, e_ps, e_cm):
        return pspec_charged_pion(e_ps, e_cm / 2.)
    
    def positron_spectrum_funcs(self):
        return "pi pi": self.dnde_pos_pipi

    def positron_lines(self, e_cm):
        return 

Lastly, we group all of these classes into a master class and we're done:
python
class HazmaExample(HazmaExampleCrossSection,
                   HazmaExamplePositronSpectra,
                   HazmaExampleSpectra,
                   Theory):
    # Model parameters are DM mass: mx, 
    # Wilson coefficients: c1, c2 and
    # cutoff scale: lam
    def __init__(self, mx, c1, c2, lam):
        self.mx = mx
        self.c1 = c1
        self.c2 = c2
        self.lam = lam
    
    @staticmethod
    def list_annihilation_final_states():
        return ['pi pi', 'pi0 pi0']

Now we can easily compute gamma-ray spectra, positron spectra and limit on our new model from gamma-ray telescopes. To implement our new model with m_χ = 200 MeV, c_1 = c_2 = 1 and Λ = 100 GeV, we can use:
python
>>> model = HazmaExample(200.0, 1.0, 1.0, 100e3)

To compute a gamma-ray spectrum:
python
# Photon energies from 1 keV to 1 GeV
>>> e_gams = np.logspace(-3.0, 3.0, num=150)
# Assume the DM is moving with a velocity of 10^-3
>>> vdm = 1e-3
# Compute CM energy assuming the above velocity
>>> Q = 2.0 * model.mx * (1 + 0.5 * vdm**2)
# Compute spectra
>>> spectra = model.spectra(e_gams, Q)

Then we can plot the spectra using:
python
>>> plt.figure(dpi=100)
>>> for key, val in spectra.items():
...     plt.plot(e_gams, val, label=key)
>>> plt.xlabel(r'E_γ (MeV)', fontsize=16)
>>> plt.ylabel(r'dN/dE_γ (MeV^-1)', fontsize=16)
>>> plt.xscale('log')
>>> plt.yscale('log')
>>> plt.legend()

Additionally, we can compute limits on the thermally-averaged annihilation cross section of our model for various DM masses using 
python
# Import target and background model for the E-Astrogam telescope
>>> from hazma.gamma_ray_parameters import gc_target, gc_bg_model
# Choose DM masses from half the pion mass to 250 MeV
>>> mxs = np.linspace(mpi/2., 250., num=100)
# Compute limits from E-Astrogam
>>> limits = np.zeros(len(mxs), dtype=float)
>>> for i, mx in enumerate(mxs):
...     model.mx = mx
...     limits[i] = model.unbinned_limit(target_params=gc_target,
...                                      bg_model=gc_bg_model)




§ MULTI-PARTICLE PHASE SPACE INTEGRATION: RAMBO



The rambo module implements the RAMBO algorithm, developed by Kleiss, Stirling and Ellis <cit.>, for generating phase-space points using a Monte Carlo procedure. We advise the interested reader to consult the original paper for details on how this algorithm works. We will only described how the algorithm works at a high-level.

The RAMBO algorithm works by first democratically generations random, independent, isotropic, massless four-momenta, q_i^μ, with i running from 1 to N, N being the number of final state particles. The momenta q_i^μ are generated with energies distributed according to q_i^0exp(-q_i^0)dq_i^0. These four momenta are then transformed into four-momenta p_i^μ by performing a particular Lorentz boost such that ∑_i p_i^μ = P^μ where P^μ is the total four-momenta of the given process in the center-of-mass frame (i.e. P^0 = E_CM and P⃗ = 0⃗.) At this stage, a corresponding weight W_0(p_i) is computed to give the probability of the event with the particular four-momenta p_i^μ. Lastly, the four-momenta p_i^μ are transformed into four-momenta k_i^μ such that the k_i^μ have the correct masses (i.e. k_i^μk_i,μ = m_i^2) and a new weight W(k_i) is computed. This completes the algorithm. Below we give pseudo-code for the algorithm:

    
  * Generate the q_i^μ's with energies distributed according to q^0_iexp(-q^0_i)dq^0_i using:
    
    q_i^0   = -log(ρ_i^(3)ρ_i^(4)),     q_i^x   = q_i^0√(1-c_i^2)cosϕ_i,     q_i^y   = q_i^0√(1-c_i^2)sinϕ_i,     q_i^z   = q_i^0c_i.

    where c_i = 2ρ_i^(1) - 1, ϕ_i=2πρ_i^(2) and ρ_i^(1,2,3,4) are uniform random numbers on (0,1). See <cit.> for a proof that these are indeed distributed according to q^0_iexp(-q^0_i)dq^0_i.
    
  * Transform the q_i^μ's into p_i^μ's such that the ∑_ip_i^μ = P^μ with P^0 = E_CM and P⃗ = 0⃗. This is done using:
    
    p_i^0   = x(γ q^0_i+b⃗·q⃗_i),    p⃗_i   = x(q⃗_i+ b⃗q_i^0 + a (b⃗·q⃗_i)b⃗).

    where 
    
    b⃗   = -Q⃗/M,     x    = E_CM/M,    γ   = Q^0/M,     a    = 1/√(1+γ),
     
            Q^μ   = ∑_iq_i^μ,     M    = √(Q_μQ^μ)

    Assign the following weight to the event:
    
    W_0   = (π2)^N-1E_CM^2N-4(2π)^4-3NΓ(N)Γ(N-1)

    
  * Next we transform the p_i^μ's to k_i^μ's such that k_i^μk_i,μ = m_i. This is done via:
    
    k⃗_i   = ξp⃗_i,     k_i^0   = √(m_i +(ξ p_i^0)^2)

    where ξ is the solution to 
    
    E_CM   = ∑_i=1^N√(m_i^2 + (ξ p_i^0)^2)

    This equation can be solved efficiently using Newton's method with the initial guess ξ_0 = √(1 - (∑_i=1^Nm_i/E_CM)^2). Lastly, we re-weight the probability of the event using:
    
    W = W_0 E_CM(∑_i=1^N|k⃗_i|E_CM)^2N-3(∑_i=1^N|k⃗_i|^2k_i^0)^-1(∏_i=1^N|k⃗_i|k_i^0)

    We note that this fixes a typo from the original paper, where the overall factor of E_CM was missing. Additionally, we note that the weights clearly don't take into account the matrix element of the particle physics process. The matrix elements are included by simply multiplying the weights by the corresponding matrix element. 

This completes the pseudo-code for the RAMBO algorithm. 

In hazma, the rambo module can be used to perform various tasks. At the highest level, the function compute_annihilation_cross_section (compute_decay_width) is provided for computing cross-sections (decay widths) for 2→ N (1→ N) processes. For example, consider computing the partial decay width of μ^- → e^- ν_e ν_μ. One first declares a function that takes a list of the final state particles' four-momenta and muon's four-momentum, and returns the matrix element squared for the reaction:
python
# Import the fermi constant
>>> from hazma.parameters import GF
# Import a helper function for scalar products of four-vectors
>>> from hazma.field_theory_helper_functions.common_functions import ...    minkowski_dot as MDot
# Declare the matrix element
>>> def msqrd_mu_to_enunu(momenta):
...    pe = momenta[0] # electron four-momentum
...    pve = momenta[1] # electron-neutrino four-momentum
...    pvmu = momenta[2] # muon-neutrino four-momentum
...    pmu = sum(momenta) # muon four-momentum
       # squared matrix element
...    return 64. * GF**2 * MDot(pe, pvmu) * MDot(pmu, pve)

Then, the partial decay width can be computed using:
python
# Import function to compute decay width
>>> from hazma.rambo import compute_decay_width
# import masses of muon and electron
>>> from hazma.parameters import muon_mass as mmu
>>> from hazma.parameters import electron_mass as me
# Specify the masses of the electron and neutrinos
>>> fsp_masses = np.array([me, 0., 0.])
# compute the partial width
>>> partial_width = compute_decay_width(fsp_masses, mmu, num_ps_pts=50000,
...                                     mat_elem_sqrd=msqrd_mu_to_enunu)

Using 50000 phase-space points gives a result within 5% of the analytical value.

In addition, rambo includes a function for performing integrations over all variables except the energy of one of the final-state particles called generate_energy_histogram. This is useful for computing energy spectra for FSR from final states with large numbers of particles. The following demonstrates how it can be used:
python
>>> from hazma.rambo import generate_energy_histogram
>>> import numpy as np
>>> num_ps_pts = 100000 #number of phase-space points to use
# masses of final-state particles
>>> masses = np.array([100., 100., 0.0, 0.0])
>>> cme = 1000. # center-of-mass energy
>>> num_bins = 100 # number of energy bins to use
# computing energy histograms
>>> eng_hist = generate_energy_histogram(num_ps_pts, masses, cme,
...                                  num_bins=num_bins)
# plot the results
>>> import matplotlib as plt
>>> for i in range(len(masses)):
...     # pts[i, 0] are the energies of particle i
...     # pts[i, 1] are the probabilities
...     plt.loglog(pts[i, 0], pts[i, 1])                                   


At the lowest level, rambo includes a function for generating phase-space points called generate_phase_space. We also include a function integrate_over_phase_space which will perform the integral

    ∫(∏_i=1^Nd^3p⃗_i(2π)^312E_i)(2π)^4δ^4(P-∑_i=1^Np_i)|ℳ|^2

where P^μ is the total four-momentum, p_i^μ are the individual four-momenta for each of the N final-state particles and ℳ is the matrix element. 



§ GAMMA-RAY SPECTRA FROM MANY FINAL STATE PARTICLES: GAMMA_RAY



Since computing gamma-ray spectra is a model-dependent process, we include in hazma tools for computing gamma-ray spectra from both FSR and the decay of final state-particles. The gamma_ray module contains two functions called gamma_ray_decay and gamma_ray_fsr. The gamma_ray_decay function accepts a list of the final-state particles, the center-of-mass energy, the gamma-ray energies to compute the spectrum at and optionally the matrix element squared. Currently, the final-state particles can be π^0, π^±, μ^±, e^±, K^±, K_L and K_S where K stands for kaon. We caution that when including many final-state mesons, one needs to take care to supply a matrix element squared that is valid at the reaction's center-of-mass energy (see Sec. (<ref>) for a discussion on the validity of ChPT for large energies).

The gamma_ray_decay function works by first computing the energies distributions of all the final-state particles and convolving these with the decay spectra of the final-state particles. It can be used as follows:
python
>>> from hazma.gamma_ray_decay import gamma_ray_decay
>>> import numpy as np
# specify the final-state particles
>>> particles = np.array(['muon', 'charged_kaon', 'long_kaon'])
# choose the center of mass energy
>>> e_cm = 5000.
# choose list of the gamma-ray energies to compute spectra at
>>> e_gams = np.logspace(0., np.log10(e_cm), num=200, dtype=np.float64)
# compute the gamma-ray spectra assuming a constant matrix element
>>> spec = gamma_ray_decay(particles, e_cm, e_gams)



The gamma_ray_fsr function computes the gamma-ray spectrum from X→ Yγ, i.e.:

    dN(X→ Yγ)dE_γ = 1σ(X→ Y)dσ(X→ Yγ)dE_γ

where X and Y are any particles excluding the photon. As input this function takes a list of the initial state particle masses (either 1 or 2 particles), the final state particle masses, the center-of-mass energy, a function for the tree-level matrix element squared (for X→ Y) and a function for the radiative matrix element squared (X→ Yγ). The functions for the squared matrix elements must take is a single argument which is a list of the four-momenta for the final state particles. As an example, we consider the process of two dark-matter particles annihilating into charged pions, χχ̅→π^+π^-(γ) using the model from App. (<ref>). In App (<ref>), we gave the analytic expressions for the gamma-ray spectra. The tree-level and radiative matrix elements squared for this process are:

    |ℳ(χχ̅→π^+π^-)|^2    = c_1^2 (s-4 m_χ^2)/2 Λ^2
    
        |ℳ(χχ̅→π^+π^-γ)|^2    = 2 c_1^2 (4 μ_χ^2-1) Q^2 e^2Λ^2 (t-μ_π^2 Q^2)^2
       (u-μ_π^2 Q^2)^2
           ×((μ_π Q (t+u)-2 μ_π^3 Q^3)^2+s
       (t-μ_π^2 Q^2) (μ_π^2
       Q^2-u))

where Q is the center-of-mass energy, e is the electromagnetic coupling, μ_π,χ = m_π,χ/Q and 

    s    = (p_π,1 + p_π,2)^2,    t    = (p_π,1 + k)^2,     u    = (p_π,2 + k)^2

with p_π,1,2 are the four-momenta of the two final-state pions and k is the four-momenta of the final-state photon. Below, we create a class to implement functions for the tree and radiative squared matrix elements. Note that these functions take in an array of four-momenta. 
python
from hazma.field_theory_helper_functions.common_functions import     minkowski_dot as MDot
    
class Msqrd(object):
    def __init__(self, mx, c1, lam):
        self.mx = mx # DM mass
        self.c1 = c1 # effective coupling of DM to charged pions
        self.lam = lam # cut off scale for effective theory

    def tree(self, momenta):
        ppi1 = momenta[0] # first charged pion four-momentum
        ppi2 = momenta[1] # second charged pion four-momentum
        Q = ppi1[0] + ppi2[0] # center-of-mass energy
        return -((self.c1**2 * (4 * self.mx**2-Q**2)) / (2 * self.lam**2))

    def radiative(self, momenta):
        ppi1 = momenta[0] # first charged pion four-momentum
        ppi2 = momenta[1] # second charged pion four-momentum
        k = momenta[2] # photon four-momentum
        Q = ppi1[0] + ppi2[0] + k[0] # center-of-mass energy

        mux = self.mx / Q
        mupi = mpi / Q

        s = MDot(ppi1 + ppi2, ppi1 + ppi2)
        t = MDot(ppi1 + k, ppi1 + k)
        u = MDot(ppi2 + k, ppi2 + k)

        return ((2*self.c1**2*(-1 + 4*mux**2)*Q**2*qe**2 *
                 (s*(-(mupi**2*Q**2) + t)*(mupi**2*Q**2 - u) +
                  (-2*mupi**3*Q**3 + mupi*Q*(t + u))**2)) /
                (self.lam**2*(-(mupi**2*Q**2) + t)**2*
                (-(mupi**2*Q**2) + u)**2))

Next, we can compute the gamma-ray spectrum for χχ̅→π^+π^-γ using:
python
>>> from hazma.gamma_ray import gamma_ray_fsr
# specify the parameters of the model
>>> params = 'mx': 200.0, 'c1':1.0, 'lam':1e4
# create instance of our Msqrd class
>>> msqrds = Msqrd(**params)
# specify the initial and final state masses
>>> isp_masses = np.array([msqrds.mx, msqrds.mx])
>>> fsp_masses = np.array([mpi, mpi, 0.0])
# choose the center-of-mass energy
>>> e_cm = 4.0 * msqrds.mx
# compute the gamma-ray spectrum
>>> spec = gamma_ray_fsr(isp_masses, fsp_masses, e_cm, msqrds.tree,
                         msqrds.radiative, num_ps_pts=500000, num_bins=50)
# plot the spectrum
>>> import matplotlib.pyplot as plt 
>>> plt.figure(dpi=100)
>>> plt.plot(spec[0], spec[1])
>>> plt.yscale('log')
>>> plt.xscale('log')
>>> plt.ylabel(r'dN/dE_γ (MeV^-1)', fontsize=16)
>>> plt.xlabel(r'E_γ (MeV)', fontsize=16)


apsrev

