---
author:
- |
    M. Giannotti$^a$, L.D. Duffy$^b$ and R. Nita$^a$\
    Physical Sciences, Barry University,\
    Miami Shores, FL 33161, U.S.A,\
    Los Alamos National Laboratory,\
    Los Alamos, NM 87545, U.S.A.\
    E-mail:
bibliography:
- 'heavyaxions.bib'
title: 'New constraints for heavy axion-like particles from supernovae'
---

Introduction
============

In this work, we analyze the gamma ray flux expected on Earth from the decay of heavy pseudoscalar particles produced in supernovae (SNe). The analysis of this flux and comparison with existing experimental observations provide new limits on these particles for the mass range from $\sim 100\,{\mbox{eV}}$ to $\sim 1\,$GeV.

Pseudoscalar particles occur in many extensions of the standard model. Arguably, the best known pseudoscalar is the axion [@Peccei:1977hh; @Peccei:1977ur; @Weinberg:1977ma; @Wilczek:1977pj], whose existence is a natural consequence of the Peccei-Quinn (PQ) mechanism for solving the strong CP problem. The axion emerges as a pseudo-Nambu-Goldstone boson from the axial Peccei-Quinn (PQ) symmetry, $U_{PQ}(1)$, spontaneously broken at the PQ (or axion) scale, ${f_{a}}$. This scale determines the strength of the axion interactions, as can be seen from the Lagrangian terms, $$\label{a-couplings}
\mathcal{L}_{\rm int} = -ig_{ak}a \bar\psi_k\gamma_5\psi_k +
{\frac{g_{a\gamma }}{4}} aF_{\mu\nu}\tilde F^{\mu\nu} + .... \, .$$ In Eq. (\[a-couplings\]), $g_{ak}=C_km_k/{f_{a}}$ is the axion-fermion coupling, with $m_k$ the fermion mass, and ${g_{a \gamma}}=C_\gamma \alpha/(2\pi{f_{a}})$ is the axion-photon coupling, where $\alpha$ is the fine structure constant. The constants $C_i$ are model-dependent parameters typically, but not necessarily, of order 1.

In the original axion models [@Weinberg:1977ma; @Wilczek:1977pj; @Kim:1979if; @Shifman:1979if; @Zhitnitsky:1980tq; @Dine:1981rt], the axion mass is related to the PQ constant via $$\label{mass}
    m_a\simeq 6.0\;\mu\mathrm{eV}\left(\frac{f_a}{10^{12}\;\mathrm{GeV}}\right)^{-1}\,.$$ These *standard axion* models have been studied extensively in the past 30 years (see, e.g., ref. [@Kim:2008hd] for a recent review). A pseudoscalar particle which possesses some, but not all, of the properties of the standard axion is known as an *axion-like* particle. These particles have recently been the subject of considerable interest (see, e.g., ref. [@Jaeckel:2010ni] for a recent review).

For the standard axion, experimental and astrophysical constraints exclude values of the PQ constant below $10^9\, {\mbox{GeV}}$ [@Raffelt:1996wa]. By Eq. (\[mass\]), the axion is then very light, with $m_a\lesssim 10^{-2}$ eV (see refs. [@Kim:2008hd; @Duffy:2009ig] for recent discussions). However, the existence of a heavy pseudoscalar particle is not excluded if relation (\[mass\]) is relaxed. Several models for these *heavy-axion-like* particles (HALPs) have been proposed in the past. Examples include models of large extra dimensions [@ArkaniHamed:1998rs; @ArkaniHamed:1998nn; @Dienes:1999gw; @DiLella:2000dn], where heavy HALPs emerge as Kaluza-Klein (KK) modes of the standard *light* axion, models with Plank-scale induced mass terms [@Berezhiani:1999qh] and models of axions interacting with a hidden sector [@Rubakov:1997vp; @Berezhiani:2000gh; @Gianfagna:2004je]. It is, therefore, phenomenologically interesting to explore the axion parameter space relaxing relation (\[mass\]).

Observe that in some cases, for example in the models [@Berezhiani:2000gh; @Gianfagna:2004je], the HALP is a QCD axion, which effectively solves the strong CP-problem. Therefore, the term *axion-like* particle could be somewhat misleading. However, in general, these heavy pseudoscalar particles are not responsible for the solution of the strong CP problem, and cannot be considered QCD axions. Rather than focusing on the specific models, in this work we keep a phenomenological approach, and consider a general HALP, described by two independent parameters, the mass, $m_a$, and the PQ scale, ${f_{a}}$. We assume that this particle interacts with nucleons and photons as described by Eq. (\[a-couplings\]) with $C_{N}=1$ and $C_\gamma=0.7$. Our normalization choice for $C_\gamma$ is that of a DFSZ axion [@Zhitnitsky:1980tq; @Dine:1981rt] and is used in order to enable comparison with the existing literature. For the sake of simplicity, we neglect a possible HALP-lepton interaction. As we discuss later, the existence of such an interaction would not significantly change our results.

Presently, the strongest constraints on such HALPs are derived from the analysis of the effects of the photons from HALPs decay on some cosmological observations such as the extragalactic photon background, the cosmic microwave background, and the deuterium abundance [@Masso:1997ru]. This argument is based on the hypothesis that the HALPs are only allowed to decay into photons. The analysis applies to the mass region between 1 keV and 1 GeV, and ${f_{a}}$ is restricted to being greater than $10^{10}{\mbox{GeV}}$ (see Fig. 1 in ref. [@Masso:1997ru]) for an axion-photon coupling with $C_\gamma=0.7$.

In this paper, we extend the bound on ${f_{a}}$ to greater values for axion masses in a similar range, using Fermi-LAT observations of supernovae photons.

Weakly interacting pseudoscalars may be produced via nucleon bremsstrahlung in the core of supernovae. Due to their weak interactions, they would freely stream out and eventually decay into photons. The decay probability depends on both the HALP’s mass and coupling with photons. In general, a fraction of the HALPs produced would decay before reaching Earth, generating an observable gamma-ray flux. The signal, peaked at energies of the order of the SN core temperature, would be stretched on a much longer timescale than the corresponding neutrino signal, and so the decay photons would still be observable today. We have analyzed the gamma-ray signal expected from the decay of massive HALPs produced by SN 1987A and Cassiopeia A (Cas A). We find that, for a large region of the axion parameter space, the gamma ray flux from the directions of SN 1987A and Cas A should be easily observable by the Fermi Large Area Telescope (Fermi-LAT). Observations of this signal would provide information on the properties of these particles. The lack of such observations allows us to place limits on the HALP-photon coupling (or equivalently, $f_a$) as a function of the HALP mass. We find that for ${f_{a}}=10^{12}\,{\mbox{GeV}}$, the HALP mass is constrained in the range from $100\,{\mbox{eV}}$ to $1 \,{\mbox{GeV}}$. The most stringent constraint on the Peccei-Quinn scale applies to HALP masses of approximately 100 MeV and restricts ${f_{a}}$ to values greater than about $10^{15}\,{\mbox{GeV}}$.

This paper is organized as follows. In the next section, we analyze the HALP production from the SN core, calculate the corresponding photon flux from their decay, and compare our results to the Fermi-LAT observations, thereby obtaining limits for HALP parameters. In section \[sec:Conclusions\], we review and discuss our results. Finally, in the appendix we provide a quantitative analysis of the effect of a lepton decay channel, neglected in the rest of the paper.

Heavy Axion-like Particles and Supernovae {#sec:HeavyAxionsAndSupernovae}
=========================================

In this section, we will first calculate the HALP luminosity from SNe, disregarding relation (\[mass\]), and then the flux of photons that we could expect from HALP decay on Earth. Combined, these calculations allow us to derive new limits on HALP properties, based on the possible detection of the emitted photons from HALP decay.

Generalization of the energy loss mechanism {#sec:GeneralizationOfTheEnergyLossMechanism}
-------------------------------------------

Massive pseudoscalars may be produced efficiently in a SN core, due to the high temperature and density. With temperatures of $T\sim 30\,{\mbox{MeV}}$ and densities of $\rho \sim 3\times 10^{14}~{\rm g~ cm}^{-3}$, the most efficient axion production mechanism in a SN core is the nucleon-nucleon bremsstrahlung process [@Raffelt:1996wa] $$\label{brem}
    N~N\to N~N~a~,$$ where $N$ is a nucleon. The bremsstrahlung production of massive axions has been studied in ref. [@Giannotti:2005tn], using the one pion exchange approximation. We use the same formalism to obtain our HALP emission rates. Setting $C_N\simeq 1$ to normalize the HALP-nuclear coupling, and a SN core of radius of $10$ km, we obtain the production rate of massive pseudoscalar in SN, $$\label{eq:emission_number}
    \frac{dN}{d\omega}\simeq \frac{1.25\times 10^{49}\,\left({\mbox{MeV}}\,{{\rm s}}\right)^{-1}}{{f_{12}}^2} 
     \left(\frac{T}{30\;{\mbox{MeV}}}\right)^{3/2}\left(\frac{\rho}{3\times 10^{14}\,{{\rm g}\,{\rm cm}^{-3} }}\right)^{2}
     \beta^3 \phi \left(\frac{\omega}{2T}\right)\,,$$ where $dN$ is the number of axions with energy between $\omega$ and $\omega+d\omega$, $\beta$ is the axion velocity as a fraction of the speed of light, and ${f_{12}}={f_{a}}/10^{12}{\mbox{GeV}}$. The function $\phi(x)$ is given by $$\begin{aligned}
\label{phi}
    \phi(x)= \frac52\,x^2\,e^{-x}\, K_1(x)\,,\end{aligned}$$ where $K_1$ is the modified Bessel function of the second kind. The normalization of $\phi$ is such that $\int_0^\infty \phi\,dx=1$. In the relativistic case ($\beta\simeq1$), all the energy dependence is contained in $\phi(x)$ which then represents the axion spectrum. The maximum of $\phi(x)$ occurs at $x\simeq 0.62$, approximately corresponding to an energy of $1.2 \,T$. At high energies, $$\phi \left(\frac{\omega}{2T}\right) \simeq \frac52\sqrt{\frac{\pi}{2}}\,
    \left(\frac{\omega}{2T}\right)^{3/2}e^{-\omega/T}\,,\qquad (\omega\gg T)
    \label{eq:phihit}$$ and thus, for $m_a\gg T$, the production rate is suppressed by an exponential factor. This is the origin of the Boltzmann suppression factor expected for the production of heavy particles in SNe. Additionally, in the nonrelativistic regime, the production rate of Eq. (\[eq:emission\_number\]) is suppressed by the velocity factor, $\beta^3$.

As we are interested in large values of the PQ constant, corresponding to weak coupling, we assume the HALPs free stream from the SN core and decay into photons only when they are free of the SN matter. In the free-streaming regime, the total energy loss rate via HALPs, or HALP luminosity, is $$\label{luminosity}
    L_a=\int_{m_a}^\infty \omega \,\left(\frac{d N}{d\omega}\right)d\omega\,,$$ with $dN/d\omega$ given by Eq. (\[eq:emission\_number\]). Now we can easily generalize the novel particle cooling limit of SN 1987A, originally derived for a massless axion, to a general HALP. Requiring that the novel particle luminosity of Eq. (\[luminosity\]), evaluated at typical core condition 1 s after the collapse ($T\sim 30\,{\mbox{MeV}}$, $\rho \sim 3\times 10^{14}~{\rm g~ cm}^{-3})$ does not exceed $3\times 10^{52}$ erg/s, we find the mass-dependent constraint on the PQ constant shown in Fig. \[fig:Contour1\]. Our result is in agreement with the standard bound, constraining $f_a$ to be greater than $10^9$ GeV for masses less than 30 MeV. As expected, the bound is strongly reduced for larger masses.

Constraints from decay into photons {#sec:ConstraintsFromDecayIntoPhotons}
-----------------------------------

In the rest of this section, we discuss the stronger phenomenological constraint on axion-like particles which arises from their subsequent decay into photons. Weakly coupled particles will freely stream out of the SN core, and, if massive enough, decay into photons before reaching Earth, generating an observable photon flux. This flux would be observable for a long time after the SN explosion, as massive particles propagate more slowly than neutrinos and photons produced directly in the SN, and their decay time is statistically distributed. Therefore, a flux of gamma photons would still be observable today.

We compare the predicted photon flux from decaying HALPs with the observations of gamma-rays from the Fermi Large Area Telescope, a detector designed to detect high-energy gamma-rays up to $300{\mbox{GeV}}$, with an angular resolution of $1^\circ$ or so at energies of a few $100{\mbox{MeV}}$ [@Atwood:2009ez].

We study both SN 1987A and Cas A. We are interested in the photon signal from the decaying-HALP observable today, after a time of some $10^{9}{{\rm s}}$ from the direct observation of the SN explosion. As the SNe distances are of the order $10^{12}{{\rm s}}$, decay-produced photons coming from an angle larger than $\sim 1^\circ$ will not have had time to reach Earth. Therefore, both sources would appear as point sources, if observed by Fermi-LAT.

Cas A, which seems to correspond to the Flamsteed’s SN of 1680, is the closest known core-collapse SN [@Pavlov:1999et] from which a flux of gamma rays in the range we are interested in has not been observed [@Hannestad:2001xi]. Recently observed gamma rays from Cas A [@Funk:2010cg] are too energetic to be caused by the mechanism that we are describing here. No gamma-ray counterpart of SN 1987A has yet been detected. As they are not observed at the energies we are interested in, we compare our predicted fluxes with the point source sensitivity of the Fermi-LAT, $< 6 \times 10^{-9} {\rm photons\,cm}^{-2} {\rm s}^{-1}$ for photons of energies $\omega_\gamma \geq 100$ MeV, at the latitudes of the SNe we are considering. As we shall see below, the analysis from Cas A will provide the most stringent bounds on the axion-photon coupling, because of its closer distance to Earth.

We now compute the gamma-ray flux on Earth from the decay of massive pseudoscalar particles produced in SNe. In the following, we assume that the HALP couples only to photons and nucleons. The definitions of these coupling strengths can be found after Eq. (\[a-couplings\]).

The HALP lifetime is then given by $$\label{tau}
    \tau\equiv  \tau_{a\to \gamma \gamma}
    \simeq 1.7\times 10^{17}\, {f_{12}}^2\, \left(\frac{m_a}{{\mbox{MeV}}}\right)^{-3}\, {\rm s}\,.$$ SN 1987A occurred at a distance of ${{d_{\rm SN}}}\simeq 51.4$ kpc from Earth or approximately $5.3\times 10^{12}$ s in natural units ($c=1$). Cas A is located ${{d_{\rm SN}}}\simeq 3.4$ kpc or $3.5\times 10^{11}$ s from Earth. Comparing these numbers with the lifetime of Eq. (\[tau\]), we see that a sizable fraction of HALPs produced can decay into photons before reaching Earth if their mass is large enough. The actual number of HALPs that decay into photons after a time $t_a$ from production is $$\label{eq:axion_decay}
    dN_{a\to\gamma\gamma}=\left[1-\exp{\left(-\frac{t_a}{\gamma\tau}\right)}\right]dN\,,$$ where $\gamma=\omega/m_a$ is the HALP Lorentz factor and $\omega$ is the HALP energy.

In order to calculate the expected photon flux today, we need to find the HALP decay time $t_a$ in Eq. (\[eq:axion\_decay\]). In natural units, the distance traversed by the HALP before decay, $\beta t_a$, and that by the photons produced by HALP decay which travel to Earth, $t_\gamma$, are geometrically related to the SN distance, ${{d_{\rm SN}}}$, via $$(\beta t_a)^2+t_\gamma^2+2(\beta t_a) t_\gamma \cos{\theta}=d_{\rm SN}^2 \,,$$ where $\theta$ is the photon emission angle. This can be solved for $t_a$ in terms of known parameters,

$$\label{eq:ta}
t_a=\frac{\left(d_{\rm SN}+\Delta t\right)\left(1-x\beta\right)
-\sqrt{
d_{\rm SN}^2\left(1-x\beta\right)^2
+\beta^2\Delta t\left(2d_{\rm SN}+\Delta t\right)\left(x^2-1\right)}}
{1-2x\beta+\beta^2} \,,$$

where $x=\cos{\theta}$ and $\Delta t=t_a+t_\gamma-{{d_{\rm SN}}}$ is the HALP *delay time*. This is the time between the signal seen on Earth from a photon directly produced in the SN and that from the photons produced via HALP decays. In this analysis, the delay time is that between direct observation of each SN and present day. For SN 1987A, $\Delta t\simeq 22$ years or $6.9\times10^{8}$ s and for Cas A, $\Delta t\simeq 320$ years or $1.0\times10^{10}$ s.

Notice that for both SN 1987A and Cas A, $\Delta t\ll {{d_{\rm SN}}}$. For the purpose of our calculation, we can therefore ignore the terms of order $(\Delta t/{{d_{\rm SN}}})^2$ or higher, which contribute less than $1\%$ to the time $t_a$. With this simplification $t_a$ becomes: $$\label{eq:ta2}
t_a\simeq\frac{\Delta t}{1-x\beta} +O\left(\frac{\Delta t}{{{d_{\rm SN}}}}\right)^2 \,.$$ We now proceed to calculate the resulting photon spectrum from HALPs decay. Here, we will make the simplifying assumption that the two photons are emitted with the same energy, $\omega_\gamma=\omega/2$. This assumption is reasonable in both the nonrelativistic and the ultrarelativistic limits, and numerical tests show that the effect of this approximation on the calculation of the photon flux is small. The resulting photon spectrum is given by $$\label{g_spectrum}
    \frac{dN_{\gamma}}{d\omega_\gamma}=
    4 \int_{-1}^{1} \left[1-\exp{\left(-\frac{t_a}{\gamma\tau}\right)}\right]\frac{dN}{d\omega} f(x) \,dx \,,$$ where $$\label{eq:angle}
    f(x)=\frac12 \frac{1-\beta^2}{(1-\beta x)^2}$$ is the normalized probability distribution for the photon emission angle [@Landau].

The expected HALP-produced photon flux on Earth can now be calculated from Eqs. (\[eq:ta\]), (\[g\_spectrum\]), and (\[eq:angle\]) in terms of the axion mass, the delay time $\Delta t$, and the SN distance ${{d_{\rm SN}}}$. Assuming isotropic emission, massive axion-like particles would produce a flux of high energy photons, $$\label{g_flux}
    \frac{dF_{\gamma}}{d\omega_\gamma}=\frac1{4\pi d_{\rm SN}^2}\frac{dN_{\gamma}}{d\omega_\gamma}\,.$$ In order to compare our prediction with the Fermi-Lat data, we need to calculate the expected number of photons with energy greater that $100~{\mbox{MeV}}$ generated by the HALPs decay, $$\label{Eq_F100}
    F_{100}=\int_{100{\mbox{MeV}}}^{\infty}\left(\frac{dF_{\gamma}}{d\omega_\gamma}\right)\,d\omega_\gamma \,,$$ as a function of the HALP mass, $m_a$, and the PQ constant, ${f_{a}}$.

Notice that since we are only interested in photons with energy greater than $100$ MeV, we can use the approximation (\[eq:phihit\]) in Eq. (\[Eq\_F100\]). Setting $T=30$MeV and $\rho=3\times 10^{14}\,{{\rm g}\,{\rm cm}^{-3} }$, and defining $\omega_1=\omega_\gamma/1{\mbox{MeV}}$ and $d_{\rm 87A}=51.4$ kpc, we find $$\begin{aligned}
\label{Eq_F100_2}
      F_{100}\simeq \frac{3\,{f_{12}}^{-2}}{{\rm cm}^2\,{{\rm s}}}  \left(\frac{d_{\rm 87A}}{{{d_{\rm SN}}}}\right)^2        \int_{100}^{\infty}   \beta^3  \omega_1^{3/2}e^{-\omega_1/15} \,d\omega_1
         \int_{-1}^{1} \left[1-e^{-t_a/\gamma\tau}\right]   f(x)dx   \,.\end{aligned}$$ Using Eq. (\[eq:ta2\]) for the axion decay time, the angular integral can be easily estimated in both the relativistic and nonrelativistic regimes.

In the relativistic regime (which applies to axions with mass $m_a\ll 200$ MeV, since we are interested only in photons of energy above 100 MeV) the calculation of $F_{100}$ can be analytically performed to the end, with the use of very reasonable approximations. In this regime $t_a/\gamma\tau \sim \Delta t/\gamma\tau\ll 1$, and so the term in the square brackets in Eq. (\[Eq\_F100\_2\]) can be substituted with $t_a/\gamma\tau$. In addition, we can neglect the $\beta^3\simeq1$ factor in Eq. (\[Eq\_F100\_2\]). With these approximations, the integral in $dx$ becomes straightforward, and Eq. (\[Eq\_F100\_2\]) reduces to $$\begin{aligned}
\label{Eq_F100_3}
         F_{100}\simeq \frac{6.8 \times 10^{-5}}{{\rm cm}^2\,{{\rm s}}}  
         \left(\frac{d_{\rm 87A}}{{{d_{\rm SN}}}}\right)^2\left(\frac{\Delta t}{\Delta t_{\rm 87A}}\right)\,\left(\frac{m_a}{1{\mbox{MeV}}}\right)^2\,{f_{12}}^{-4}\,,
           \quad (m_a\ll 200\,{\mbox{MeV}}),\end{aligned}$$ where $\Delta t_{\rm 87A}=22$ yr.

For axions heavier than $200$MeV, the lower integration limit in the photon energy in Eq. (\[Eq\_F100\_2\]) should be substituted with $2m_a$. Since this is much larger than the SN core temperature, the exponential in the integral dominates, and so the integral is peaked in the region of low energy. In this region, the exponential term is small, $\exp[-t_a/\gamma\tau]\ll 1$, and so $$\begin{aligned}
\label{eq_heavy}
      F_{100}\simeq \frac{3\,{f_{12}}^{-2}}{{\rm cm}^2\,{{\rm s}}}  \left(\frac{d_{\rm 87A}}{{{d_{\rm SN}}}}\right)^2        \int_{\frac{m_a}{2{\tiny {\mbox{MeV}}}}}^{\infty}   \beta^3  \omega_1^{3/2}e^{-\omega_1/15} \,d\omega_1
     \,,  \qquad (m_a\gg 200{\mbox{MeV}})\,.\end{aligned}$$ The above expression cannot be explicitly calculated because of the $\beta^3$ factor. However, since the integral has to be evaluated in the region where the exponential dominates, $\omega_\gamma\geq m_a\gg 15{\mbox{MeV}}$, the expression for $F_{100}$ is going to behave as $$F_{100}\sim {f_{12}}^{-2}\, e^{-m_a/\zeta} \,,  \qquad (m_a\gg 200\,{\mbox{MeV}}),$$ with $\zeta\sim 30{\mbox{MeV}}$, predicting an exponential suppression of the photon flux for large axion mass.

We show the results of the numerical integration of Eq. (\[Eq\_F100\]) in Fig. \[fig:graph\_1\] and \[fig:graph\_2\]. In Fig. \[fig:graph\_1\], we compare the photon flux from HALP decay, $F_{100}$, with the Fermi-LAT point source sensitivity as a function of the HALP mass, for the case of ${f_{a}}=10^{12}{\mbox{GeV}}$. For comparison, we also show the detection limit from the Energetic Gamma Ray Experiment Telescope (EGRET), an older gamma-ray telescope with point source sensitivity of approximately $10^{-7} {\rm photons\,cm}^{-2} {\rm s}^{-1}$, for photon energies of $\omega_\gamma \geq 100$ MeV, at high latitudes [@Sreekumar:1997un]. Where the predictions are above the Fermi-LAT point source sensitivity, we are able to provide an exclusion limit. In this case, the region of HALP mass between $10\,{\mbox{keV}}$ and $1\,{\mbox{GeV}}$ is excluded from the lack of observation of a photon flux from SN 1987A. The data from Cas A exclude a larger mass region, from approximately $100\,{\mbox{eV}}$ to $1{\mbox{GeV}}$.

Fig. \[fig:graph\_2\] shows the complete exclusion region obtained numerically, in the $m_a$-${f_{a}}$ plane. The excluded values are those for which Fermi would be able to observe our predicted spectrum of gamma-rays. The constraint from Cas A is more stringent, due to both a smaller value of ${{d_{\rm SN}}}$ and a larger value of $\Delta t$, with respect to SN 1987A (cfr. Eq. (\[eq\_rel\_approx\]) below).

The agreement with the analytical approximation in the relativistic regime is within a few percent, and the exponential suppression in the very massive regime is also evident.

The contour line in the low mass region ($m_a\lesssim 100{\mbox{MeV}}$) shown in Fig. \[fig:graph\_2\] is, according to Eq. (\[Eq\_F100\_3\]), a straight line with equation $$\label{eq_rel_approx}
    \log \left(\frac{f_a}{1{\mbox{GeV}}}\right)\simeq 13 + 
    \frac14\log \left[\left(\frac{d_{87A}}{{{d_{\rm SN}}}}\right)^2 \left(\frac{\Delta t}{\Delta t_{\rm 87A}}\right)\right]
    +\frac12\log \left(\frac{m_a}{1{\mbox{MeV}}}\right)\,.  $$ The above relation gives a quantitative prediction of the effect of the SN distance and the delay time on the expected photon flux.

From the figures, we see that the bounds on ${f_{a}}$ are peaked around the value $m_a\sim 100$ MeV, with ${f_{a}}$ restricted to being greater than approximately $10^{15}$ GeV. The bounds are strongly suppressed at lower and higher masses. This is not surprising, as the the HALP flux is maximized if the HALP can be easily produced, but is also heavy enough to efficiently decay. This differs from the standard energy loss (or novel particle cooling) arguments, which are particularly efficient for low mass particles.

Discussion and Summary {#sec:Conclusions}
======================

Here we briefly summarize and comment on the key points of our analysis and discuss possible future directions. The main motivation for this work is the analysis of the heavy-axion-like particles (HALP) parameter space, particularly in the case when the mass and PQ constant are not directly related. We have calculated the production rate of massive HALPs in a SN core and analyzed the photon flux expected from their decay into photons. The flux has then been compared to the results of the Fermi Large Area Telescope. Our main results are shown in Figs. \[fig:graph\_1\] and \[fig:graph\_2\]. For a HALP mass of $\sim 100\,$MeV, we find that ${f_{a}}$ must be greater than roughly $10^{15}\,{\mbox{GeV}}$. For the specific case of ${f_{a}}=10^{12}\,{\mbox{GeV}}$, we find that the values of the HALP mass between approximately $100\,{\mbox{eV}}$ and $1 \,{\mbox{GeV}}$ are excluded.

The idea of constraining radiatively decaying axions from gamma-ray observations of SN 1987 A was considered already shortly after the SN 1987 A event (see, e.g., [@Kolb:1988pe]). The original analysis was applied, however, only to photons from light axions decay, which reach the earth within a few seconds from the observation of the explosion. In general, the axion bounds on ${f_{a}}$ from stellar cooling are constant for masses up to a few $T$, where $T$ is the relevant stellar temperature, and are exponentially suppressed for higher masses. Our bound, instead, is maximal for HALP masses around $100$ MeV (a few times the SN core temperature of 30 MeV), but is reduced in the high and the low mass region. In fact, in the mechanism we have studied, two different effects play a role: HALP production in the SN core, and HALP decay into photons. The production rate of heavy HALPs is exponentially suppressed as the energy required to produce HALPs increases, but the decay rate into photons increases with the third power of the mass.

We have also investigated the SN 1987A novel cooling limit that applies when the HALP mass is independent of the PQ constant. The standard energy loss argument from the analysis of the observed SN 1987A neutrino signal becomes inefficient for axion masses greater than the SN core temperature, approximately $30$ MeV. This is quantitatively shown in Fig. \[fig:Contour1\]. This limit is in agreement with previous bounds derived from SN 1987A for the standard axion.

In deriving these new bounds, we have made a number of hypotheses. First, and most importantly, we have calculated the axion production rate in the one-pion-exchange (OPE) approximation, which is not completely justifiable. As shown in [@Hanhart:2000ae], this approximation might induce an overestimate of the production rate. Although this conclusion cannot be directly applied to the results in our paper, since the analysis in [@Hanhart:2000ae] does not take into account the effects of finite axion mass, and although our bound on ${f_{a}}$ depends only on the fourth-root of the production rate, it is evident that our results should be considered within some not easily quantifiable uncertainty, as many details of the SN physics and particles production in its core are yet not completely known (see also the discussion in [@Raffelt:1996wa], page 120).

We have also assumed that our axion-like particle interacts with nucleons and photons only. We have not included the possibility of axion decay into leptons, even when this would be kinematically possible ($m_a>2m_\ell$). Even if one assumes that the pseudoscalar-lepton coupling is absent in the Lagrangian (e.g. following hadronic axion models), the coupling with photons would necessarily induce a (suppressed) coupling with electrons, and so massive axions could have another decay channel. This effect would reduce the expected photon flux and somewhat relax our bounds. Our numerical tests discussed, in the appendix, show, however, that this would not be a very significant effect for large PQ constant. Finally, we have not considered the HALP couplings to nucleons and photons as independent. These couplings are still related in the same manner as standard axion models, via the PQ constant. This may not be the case for a more generic pseudoscalar. In principle, it is possible to tune the axion-photon coupling to very small values without changing the axion coupling with nuclei [@Kaplan:1985dv]. This would, of course, modify our bounds.

Other possible arguments could be used to constrain the region of the parameter space which we have analyzed. For example, for large PQ constant and small mass, axions would have a lifetime longer than the age of the universe. Therefore, another possible constraint could be derived from the analysis of the axion thermal production, and the calculation of the axion relic abundance today. A full calculation of the HALP relic abundance, based on references [@Turner:1986tb; @Masso:2002np], shows that this is the case for $f_a\lesssim 5\times 10^{12}{\mbox{GeV}}$, and $m_a>$ a few $10^{-4}{\mbox{MeV}}$, providing that the HALP lifetime exceeds the age of the universe. Therefore, from our analysis, it results that the SN argument which we used in this work is more efficient in constraining the region of high ${f_{a}}$. A full analysis of other possible constraints in that region would, however, be very interesting and should be considered in the future.

Our bounds provide strong implications for the models which predict heavy-axion-like particles, for example, the axion models in extra dimensions. A detailed analysis of these models is beyond the purpose of this paper, and will be the topic of future work.

We thank Brenda Dingus for guidance on Fermi observations.

Appendix: Inclusion of the lepton decay channel {#sec:CommentsOnTheInclusionOfLeptonDecayChannel .unnumbered}
===============================================

In our calculations we have not included the possibility of HALP decay into leptons, even when this would be kinematically possible ($m_a>2m_\ell$). Here we want to briefly comment on the impact of the axion-lepton decay on our results.

The effect of a new decay channel would be dual: *first*, the axion lifetime would be reduced, and this would affect the exponential factor in Eq. (\[g\_spectrum\]); *second*, the decay products would not just be photons. This last corresponds to a reduction of the expected photon flux by the factor $\tau/\tau_\gamma$, where $\tau_\gamma$ is the HALP lifetime into photons given by Eq. (\[tau\]), and $\tau=(\tau_\gamma^{-1}+\tau_\ell^{-1})^{-1}$ is the lifetime including the lepton decay channel. For a DFSZ axion, $$\tau_\ell=\frac{6.7\times 10^{10}\,{f_{12}}^2}{(m_a/{\mbox{MeV}}) \sqrt{1 - (2 m_\ell/m_a)^2}}\,{{\rm s}}\,. \nonumber$$ Both effects would be negligible in the case of a very massive particle. On the other hand, the reduction factor $\tau/\tau_\gamma$ is small for small axion mass, and could potentially cause a sizable reduction of the photons flux coming from relativistic pseudoscalars. However, in general, whenever the condition $t_a/\gamma\tau\ll 1$ is satisfied, the two effects described above cancel out, as can be easily understood by expanding the exponential factor in Eq. (\[g\_spectrum\]). This happens for large values of the PQ constant, and for relatively small axion mass, as confirmed by our numerical analysis in Fig. \[fig:leptons\]. Therefore, we expect the reduction induced by the lepton decay to be sizable only in the region of intermediate axion mass, and small PQ constant.

Our numerical results are shown in Fig. \[fig:leptons\]. The figure represents the percentage of the flux of photons with energy above 100 MeV expected on Earth if the lepton decay channel is open, with respect to the flux (\[Eq\_F100\]), expected if the decay into leptons is neglected.

In the analysis we are including the decay channel into electrons and muons, and assuming an axion-lepton coupling as for the DFSZ axion. In general, it is possible to consider models where the axion-lepton coupling is much smaller (hadronic axion models), reducing even more the effects of the lepton decay channel.

Notice that, since $t_a\simeq\Delta t$, the reduction of the flux (at large ${f_{a}}$) induced by the lepton decay channel is more pronounced for Cas A, which has a larger delay time.
