# On the production of heavy axion-like particles in the accretion disks of gamma-ray bursts

###### Abstract

Heavy axion-like particles have been introduced in several scenarios beyond the Standard Model and their production should be possible in some astrophysical systems. In this study, we re-examine the possibility that this type of particle can be generated in the accretion disks of gamma-ray bursts (GRB), which are the most powerful events in the universe. If the produced axions decay into photons or pairs at the correct distances, a fireball is generated. We calculate the structure of transient accretion disks in GRBs (density, temperature and thickness profiles) considering the effect of heavy axion emission as well as the rest of the relevant standard cooling processes. This allows us to obtain the values of the coupling constant such that the axions do not become trapped, and we also compute the heavy axion luminosity emitted from the entire disk. We show that for the couplings within the ranges found, the mechanism for powering GRBs based on heavy axion production and decay is an alternative to the standard picture based on magnetohydrodynamic processes and neutrino-antineutrino annihilation. Alternatively, the mechanism fails if heavy axions are produced in the disk but their decay takes place further away. Still, the decay products (gamma rays or electrons and positrons) should leave observable signatures, which are not observed for different ranges of values of the coupling constants, depending on the mass of the heavy axion.

###### keywords:

new physics, axions, gamma-ray bursts^{†}

^{†}journal: Physics Letters B

## 1 Introduction

Gamma-ray bursts (GRBs) are among the most powerful events in the universe since the Big-Bang, emitting radiation at a rate , which is observed on Earth as flashes of gamma rays with energies and durations ranging from a fraction to hundreds of seconds (for a review, see (meszaros2006, )). GRBs that last less than two seconds are considered short GRBs and they are thought to be produced by the merging of compact objects in a binary system pacynski1986 (). Longer GRBs are associated with the collapse of a massive star into a black hole (woosley1993, ). It is generally accepted that an initial fireball of , gamma rays, and baryons is generated close to the black hole, which then expands to reach ultrarelativistic velocities (cavallorees1978, ). Hence, the observed radiation is considered to correspond to synchrotron and/or inverse Compton (IC) emission of electrons that have been accelerated to high energies in shocks created by relativistic plasma shells that collide (reesmeszaros1994, ; piran1999, ). GRBs can also be detected by observing their lower energy afterglow emission, which occurs from hours to days after the initial detection and it can last for months in some cases. This has allowed to measure the corresponding redshift , thereby confirming the extragalactic origin of GRBs (e.g. see afterglowdetection1997 (); afterglowdetection1998 ()).

As for the central engine, the black hole is assumed to be surrounded by a transient, hot, and dense accretion disk, which is considered to be cooled via advection and neutrino emission (pwf1999, ; kohri2002, ; dimatteo2002, ; kohri2005, ; chen2007, ). The high densities and temperatures in this disk as well as in the generated fireball have motivated the study of different particle physics beyond the Standard Model loeb1993 (); bertolami1998 (); demir1999 (); viaaxion2000 (); mirror2004 (); tu2015 (). In this study, we further explore the possibility that heavy axion-like particles can be emitted from such accretion disks (viaaxion2000, ; mirror2004, ; nucleonbrem2005, ).

The existence of the axion was proposed to solve the strong CP problem, which (for example) is reflected by the fact that the electric dipole moment of the neutron is unnaturally small lowenergyfrontier2010 (). As a possible solution to this problem, the standard axion arises as a pseudo-Nambu–Goldstone boson that spontaneously breaks the Peccei–Quinn symmetry () at a scale pq1977 (). This cancels out the CP-violating term of the QCD Lagrangian and the axion acquires a mass given by

(1) |

Heavy axion-like particles^{1}^{1}1Throughout the text we refer to “heavy axion-like particles” just as “heavy axions”. arise if the above condition for the mass is relaxed, and although the strong CP problem may not be solved, it is an interesting possibility that such heavy axions may indeed exist, as has been proposed in the context of several theoretical scenarios beyond the Standard Model (grandunification1997, ; extradim1999, ; strongcp2000, ). The Lagrangian terms describing the heavy axion interaction with matter are given by (viaaxion2000, )

(2) |

where are the coupling constants for interactions with fermions , with , and
parameterizes the coupling to photons. In particular, the most stringent bounds on the coupling constants and refer to the standard (light) axion, for which both constants are related and relation (1) holds lowenergyfrontier2010 (); khlopov1978 (); khlopov1991 (); khlopov1992 ().
For instance, bounds on were obtained based on the duration of SN 1987A masso1995 (), i.e. for , and a more recent analysis jaeckel2017 () concluded that for the same range of masses.
In addition, studying the production of heavy axions inside a supernova core leads to constraints on , as discussed by (newconstraints2011, ) in the case where the axions produced can escape freely outside^{2}^{2}2In particular, the results obtained by nucleonbrem2005 () imply that the mean free path for axions of mass MeV will be larger than the core radius km if for a core density and temperature .. Further bounds on the coupling to photons can be derived from cosmological considerations, i.e., heavy axions created in the early universe would decay, thereby affecting the cosmic microwave background (CMB) and the extragalactic background light (EBL), and such decays could also dilute the neutrino density or create a diffuse photon background cadamuro2012 (); millea2015 ().

In this study, as in viaaxion2000 (), we use a phenomenological approach and maintain all of the coupling constants as independent. We consider heavy axions of mass MeV which can couple to photons, electrons, and nucleons, adopting for the latter values in the range , and we make no distinction between neutrons and protons. We give two reasons why this range is suitable for consideration. First, there is an allowed window for hewett2012 (); pdg2016 (), and second, the existing bound excluding the range is derived from the observed duration of SN 1987A, although based on data with poor statistics, and without a complete understanding of the dynamics of the explosion (e.g., see cadamuro2012 ()). Hence, it may be interesting to explore further astrophysical phenomena in which such values of can also be tested. Under the assumptions stated above, the leading process for heavy axion production is nucleon–nucleon bremsstrahlung (raffeltlab1996, ), and considering non-zero couplings to leptons and/or photons allows the possibility of decays to these particles, which may have observable signatures under appropriate conditions. Studying these effects in the context of highly plausible scenarios that are considered appropriate for the generation of GRBs, our work provides complementary results to the existing constraints, which can be expressed as restrictions on the and planes for different values of .

In particular, we focus on the production of heavy axions in the accretion disks of GRBs. This possibility was studied by (viaaxion2000, ; mirror2004, ; nucleonbrem2005, ), and in the present study we compute the profiles for the density, temperature, and thickness of the accretion disk by taking into account the cooling rates of all the relevant mechanisms, including advection, neutrino emission, and heavy axion production as the new ingredient. We consider the effects of the non-zero masses of the axions and pions, following the approach described by nucleonbrem2005 (), and we then find the ranges of values for the coupling constant for which axions can escape from the disk. If heavy axions decay close to the central engine , then it is possible that a fireball is formed and this comprises a mechanism for powering the GRB, as discussed in the works mentioned. We shall also explore the possibility that the decays take place further away from the central engine, and leaving the bursts to be powered by a magnetohydrodynamical process (e.g. Blandford-Znajek bz1977 (); lee1999 ()) or by neutrino–antineutrino annihilation (pwf1999, ; zalamea2011, ). In this case, the heavy axions produced are still found to leave signatures which are not observed, because the energy dependence of the prompt spectrum would be different compared with that detected, e.g., if heavy axions decay preferably to photons. If they decay to electron–positron pairs, then the latter can produce a flux of gamma rays via IC interactions on the CMB, and this would be observed on Earth, although with a very different spectrum compared with the typical ones of GRBs.

The remainder of this work is organized as follows. In Section 2, we present the calculation of the accretion disk structure and we discuss on the values of the coupling constant for which heavy axions can escape from the disk for different accretions rates. In Section 3, we compute the contributions to the gamma-ray flux that would arise from the decaying axions and we compare them with a typical GRB photon flux. Finally, we conclude with a brief discussion in Section 4.

## 2 Structure of GRB accretion disks with heavy axion production

In both short and long GRBs, it is expected that matter falls to the newly formed Kerr black hole through a transient, hot and dense accretion disk pwf1999 (); kohri2002 (); dimatteo2002 (); kohri2005 (); chen2007 (). In these disks, energy can be efficiently liberated by advection and via neutrino emission, and the corresponding profiles for temperature, density and thickness (or scale height) can be obtained as a good approximation using steady state models dimatteo2002 (); chen2007 (); reynoso2006 (); janiuk2010 (); nonzero2016 (). In this study, we expand the model presented previously reynoso2006 () to include the new process of interest comprising the production of heavy axion-like particles. Similarly to previous studies (janiuk2010, ; nonzero2016, ; janiuk2014, ), we use the notation of Riffert and Herold riffert1999 () for the correcting factors to account for general relativistic effects due to the rotating black hole with mass and dimensionless spin parameter :

(3) | |||||

(4) | |||||

(5) | |||||

(6) |

where is the radius in cylindrical coordinates. These factors were derived from the conservation equation of the energy–momentum tensor corresponding to a viscous flow in a spacetime described by the Kerr metric. They are useful for correcting the standard expressions for the viscous shear, disk thickness, and heating rate of a Keplerian accretion disk, in order to make them valid for a disk around a rotating black hole. The Keplerian disk case is then recovered if these factors are set to one (see (riffert1999, ) for details).

The accretion rate in the disk is supposed to be constant () and mass conservation at a radius from the black hole implies that

(7) |

where is the mass surface density, is the disk mass density, is the radial velocity and is the half thickness of the disk. The latter can be written as (riffert1999, )

(8) |

and it is related to the viscous shear as

(9) |

where is the Shakura–Sunyaev viscosity coefficient shakurasunyaev1973 (). The total pressure is given by

(10) |

where the fraction of free nuclei is approximated by (kohri2005, )

with and . Electron neutrinos and antineutrinos are the ones which are produced more efficiently and they can become trapped, thereby contributing to the pressure. To describe their energy density we follow Ref. (dimatteo2002, ) and adopt the prescription:

where the neutrino optical depth is the sum of the scattering plus the absorptive contributions, , which are given in the following.

Accretion proceeds as the energy generated by friction is either advected toward the black hole or emitted by the disk. The heating rate due to viscosity can be written as

(11) |

The steady-state solutions are obtained requiring that the heating rate is equal to the total cooling rate at each radius, , including all the relevant cooling processes:

(12) |

The rate of photo-disintegration for heavy nuclei is given by pwf1999 ()

(13) |

and the cooling by advection can be approximated as narayanyi1994 ()

(14) |

Neutrino cooling occurs mainly through the electron–positron pair capture process, and , at a rate

(15) |

and also through electron–positron pair annihilation at rates

(16) | |||||

(17) |

Considering the corresponding inverse processes, the absorption optical depth can be approximated as

(18) |

and the scattering process as . Then, by using a simplified treatment for neutrino emission and transport based on a two-stream approximation (e.g., see (dimatteo2002, ; hubeny1990, )), we obtain the escaping energy rate in neutrinos as follows:

(19) |

Finally, we include the energy loss rate due to the emission of heavy axions by nucleon–nucleon bremsstrahlung according to (nucleonbrem2005, ), without neglecting the finite mass of the axions and pions,

(20) |

where is the intensity of the emitted axions,

(21) |

with

(22) |

and , where

(23) | |||||

(24) | |||||

(25) | |||||

(26) |

It is also useful to compute the mean free path for the inverse process, i.e., axion capture by nucleons, which can be obtained as (nucleonbrem2005, )

(27) |

Figure 1 shows the obtained cooling rate as a function of the density for and K for and as compared to the corresponding neutrino cooling rates. It can be seen from this plot that axion emission becomes more important than neutrino cooling for densities at the temperatures shown, which are typical values for the central part of the accretion disks in GRBs.

In order to calculate the density, thickness, and temperature of the disk as functions of the radius , we proceed as described by reynoso2006 () and first solve numerically Eqs. (8) and (9) to obtain and as functions of the density . Then, by using of the equation of state (Eq. 10), we can also obtain as a function of , and we employ this relation to evaluate the total cooling rate (Eq. 12) and equate it to the heating rate (Eq. 11 ), thereby obtaining the correct pair that satisfies the energy balance at each radius .

Parameter | Description | Values |
---|---|---|

black hole mass | ||

black hole spin | ||

viscosity parameter | ||

accretion rate | ||

axion–nucleon coupling | ||

heavy axion mass |

Figure 2 shows the results obtained for the profiles , , and in the left, middle, and right panels, respectively. The values of the parameters used are summarized in Table 1. The accretion rate is in the top panels, in the center panels, and in the bottom panels. In this plot, we assume that the mass of the heavy axions is . These results show that the density, temperature, and thickness do not change significantly with respect to the case with no axion production, at least for the coupling strength ranges considered. In the left panels, we also show the mean free path for the axions () evaluated at their mean energy (which turns out to be nucleonbrem2005 ()), compared with the thickness . Hence, we can conclude that for couplings as high as , heavy axions with mass will escape freely from the disk if , whereas for disks with higher accretion rates , axions will escape without interacting for couplings . For higher values of , axions will become trapped in the disk but we do not address such cases in the present study.

For different values of the heavy axion mass, their emission from the disk is also different, and in particular less intense for higher masses. This can be seen in Figure 3, where we plot the relative cooling parameter for axions and neutrinos as functions of the disk radius in the case of an accretion rate for a coupling , and for heavy axion masses . It can be seen from this plot that at the innermost regions of the disk, neutrino cooling becomes less efficient than at its maximum values reached further away from the center. This is due to increased neutrino trapping, which implies that the advection process becomes more significant, thereby leading to a decrease in the density and the temperature in agreement with the results by (chen2007, ; nonzero2016, ).

## 3 Heavy axion emission from GRB accretion disks: Implications

In this section, we study what is to be expected as observational consequences in the case that heavy axions are produced in GRB disks as described above. Again, as mentioned by (viaaxion2000, ; mirror2004, ; nucleonbrem2005, ), if the axions generated in the disk are to decay close enough (), the fireball can be formed more efficiently than by neutrino-antineutrino annihilation. The decay rates of these heavy axions to photons and to electron–positron pairs are

(28) | |||||

(29) |

so we find that, for instance, for axions of mass , the coupling constants must be and/or for the decays to occur at distances less than from the central black hole. Under these conditions, the standard phenomenology of the fireball model can then be employed to describe the burst, i.e., the created pairs or photons would form an optically thick plasma, which will expand due to radiation pressure and generate gamma-rays via internal shocks at larger distances where the flow becomes optically thin meszaros2006 (); piran1999 ().

We also consider the possibility that the decays can take place far from the central engine, such as in the interstellar medium or even outside the host galaxy. In this case, the GRBs will not be powered by axion decay but by other mechanism instead, such as neutrino–antineutrino annihilation or a magnetohydrodynamic mechanism (e.g. Blandford–Znajek). The high luminosity of the emitted axions, implies a high luminosity of the decay products, which would be directly observable in the case of photons, or via the IC radiation generated by the produced scattering on the CMB.

As mentioned in previous studies (viaaxion2000, ; mirror2004, ), axions would escape freely from short GRB accretion disks generated in compact merger events. Furthermore, it can be seen that they would also be capable of escaping from the collapsing star in long GRBs. To estimate the corresponding optical depth, we consider that the central temperature of a GRB progenitor is and the density is (e.g. woosley2002 ()). If we assume that the density drops on the envelope as (e.g. meszarosmurase2016 ()), then by taking , we find that the optical depth for axions of energies is much less than one in all of the cases studied:

(30) |

and thus in the present context, we can consider that also axions escape freely from the stellar envelope in long GRBs.

The luminosity of the emitted axions can be calculated as

(31) |

and an analogous expression is used for the neutrino luminosity . In the expression above, is the outer radius of the disk and the inner radius is taken as that corresponding to the last stable circular orbit,

where , , and .

We can also estimate the power that would be generated by the Blandford–Znajek process using the expression given by xue2013 ():

(32) |

where is the poloidal magnetic field near the horizon. The latter can be related to the pressure via , and the pressure is approximated by (see liu2015 ()) . For comparison, Table 2 shows the obtained values for , , and using different values of the accretion rate, mass of the heavy axion, and their coupling constant to nucleons. It is important to note that the efficiency of neutrino–antineutrino annihilations is always less than in the cases studied of black hole spin and accretion rates, as pointed out by liu2015 (). In general this process is less efficient than the Blandford–Znajek process, which varies here only with the accretion rate because we keep in all cases. The luminosity in heavy axions clearly increases with and with the accretion rate, since the latter implies a higher density and temperature.

Without considering the specific details of how the prompt emission is generated in GRBs, we can still rely on observations of the energy dependence of the detected gamma-ray flux and the observed luminosity. A standard fit to data on many bursts is the so-called Band flux band1993 ()

(33) |

with , and for the low- and high-energy indexes, we take , and based on the analysis of the Fermi collaboration using four years of data taking on many bursts with the Gamma-ray Burst Monitor fermigbm2014 (). The constant can be fixed by normalization on the total flux in a given energy band, and we consider the band and a flux which according to Ref. fermigbm2014 (), is close to the mean flux according to their analysis on the samples.

We then take a Band flux at the level mentioned above and compare it to the fluxes that would arise from the decay of heavy axions if they decay far from the central engine and before arriving on Earth. In particular, we consider the possibility that they decay primarily to photons at distances , such that for GRBs at redshift , i.e. at a luminosity distance . We consider this particular distance for illustration because the distribution of GRB with redshift exhibits its higher values for fermigbm2014 (). Given the above expressions for the decay rates, we can see that will hold, e.g., if for , and if for . Hence, the ranges of the couplings and that correspond to the channel as dominant are marked by the shaded regions in the upper plots in Fig. 4, for and in the left and right panels, respectively.

For simplicity, the flux of gamma rays produced in this case can be estimated assuming an isotropic emission as

(34) |

where and .

On the other hand, if the dominant decay channel is , then we consider decay lengths such that , i.e., up to a maximum of from the burst site, in order to consider that the interactions with the CMB are initiated at but not closer than to avoid remaining trapped in the magnetic field of the host galaxy. The shaded regions in the bottom plots in Fig. 4 indicate the ranges of the couplings and at which these decays are dominant, for and in the left and right panels, respectively.

In order to estimate the flux of scattered photons that would arrive on Earth, we follow the analytical treatment of electromagnetic cascades described by (venyakala2016, ), such that the spectrum of cascade photons initiated by an electron of energy can be generally expressed as

(35) |

Here, , , and the characteristic energies of the CMB and the EBL are and , respectively . We adopt this dichromatic approximation for the background photons, although the produced will interact mainly with the CMB radiation and the scattered photons are cascade sterile, because they are below the threshold for production. The parameter can be computed by considering energy conservation between the initial energy and the total energy of the scattered photons. We note that in the present context and for redshifts , the break energy is , so that most have energies , and hence the scattered photons have a spectrum . In these cases, it is found that , whereas if , then . According to venyakala2016 (), if we assume that the interactions occur in less than a timescale , then we can write the flux of IC scattered photons to be observed on Earth as

(36) |

Here, the intensity of decaying electrons and positrons is given by

(37) |

where is the Lorentz factor of the axion and is its velocity in units of . The minimum energy for the decay to of energy is given by

which follows when we consider that the maximum and minimum energies in the laboratory frame add up to , similarly to the discussion by stecker1971 ().

Figure 5 shows the flux contributions obtained in gamma rays from the decay of heavy axions in the cases of direct decay to photons and also if the decay to electron–positrons is dominant. In the former case, the flux to be observed would be very different from the typical Band flux, which appears in dashed grey lines, whereas the contribution from direct decays to photons is in blue. We note that even if we considered a Band flux twice as high as the level shown in the plots, the contributions of decaying axions would still introduce visible signatures into the spectra. And in this case, it appears reasonable to expect that most of the GRBs presents such a high flux or less (), according to the analysis of Fermi-GBM (see Fig. 11 in (fermigbm2014, )). Therefore, the features corresponding to the red and blue curves in Fig. 5 would have been observed in most of the GRBs, and thus for , we can exclude the ranges of and marked by the shaded regions in Fig. 4.

Now, in the cases of axions that decay dominantly to , the flux contributions due to IC on the CMB are found to be significant and more spread widely to lower energies (red curves in Fig. 5). Actually, this component of the flux is not expected to arrive from the same direction as the original burst emission, which is supposed to be beamed. This is because the electrons and positrons are deviated in the intergalactic magnetic field (e.g.(venyakala2016, )). Hence, for coupling values and within the shaded regions in the lower panels of Fig. 4, the mentioned flux component should have been clearly observed, i.e., not superimposed directionally and temporally to a Band-like flux. The lack of observations of any flux contribution such as the described allows us to exclude the aforementioned combination of values for , , and .

## 4 Discussion

In this study, we investigated the structure of accretion disks in GRB by considering the cooling term that would arise due to heavy axion production via the nucleon–nucleon bremsstrahlung process. We found values of the coupling constant to nucleons for which the axions produced can escape from the disk by comparing the mean free path with the disk thickness. For instance, for heavy axions with a mass , their coupling to nucleons can be and they could leave the disk without interacting for accretion rates , whereas for higher accretion rates (), the disk is denser and it is necessary that in order to have free streaming. In these cases, the structure of the disk does not depart significantly from the result corresponding to no axion production, although the slight changes are more noticeable for the higher values of considered.

The luminosity in heavy axions can still be important and it would lead to a more efficient production of photons and/or than via neutrino–antineutrino annihilation. This can be seen in Table 2, especially for the highest values considered for . Hence, as proposed by viaaxion2000 () but having performed a more realistic calculation of the axion luminosity, if the axions produced decay to photons and/or at distances , then the initial fireball could be generated and give rise to the observed GRB. For example, this would be the case if for dominant decays to photons and for dominant decays to , with . However, the former possibility has been excluded by beam dump and collider experiments beamdump1988 (); bauer2017 (), constraints based on the primordial D/H ratio (millea2015, ), and the duration of SN 1987A. Figure 6 shows these and the remaining current bounds in the plane, which are valid in the case of heavy axions decaying to photons as the dominant channel. In this figure, a black dotted line denotes the required values of so that the typical decay distance of the heavy axions produced is , which shows that higher coupling values which would yield , but they have already been excluded for the relevant range of heavy axion mass.

On the other hand, the possibility of dominant decays to to power GRBs has not yet been excluded completely. This can be seen in Fig. 7, where the Boxerino bound extends up to , and the black dotted line denotes the minimum values of required for decays at distances shorter than , thereby allowing the possibility that heavier axions decay to to create the initial GRB fireball provided that , which has been excluded by beam dump experiments. We leave for future work a detailed study of any possible effects regarding axion production within the expanding fireball for these particular cases of and .

We also studied the cases in which heavy axions decay at longer distances. In the case that the decay to photons is dominant, the flux generated would lead to a clearly different energy dependence compared with that typically observed, which is usually fitted by the so-called Band model (Fig. 5). This would occur for taking the values mentioned above and for the couplings and within the shaded regions in Fig. 4 (top panels). In addition, we considered the possibility that the decay channel to electron–positron pairs is dominant. In these cases, a flux of gamma rays would be generated by IC scattering on the CMB, and these photons would arrive at the Earth from a different direction compared with that of the usual GRB prompt emission, which is explained by the deflection of electrons and positrons in the intergalactic magnetic field.

In Fig. 5, the GRB reference flux corresponds to a luminosity of in the energy band for a GRB at a redshift . In fact, similar plots would be obtained for higher redshifts, and it is reasonable to expect that accretion disks with higher values of become more feasible, so the conditions for heavy axion production would remain at a significant level even for .

The lack of observations of the signatures described along with the assumption that GRBs actually involve accretion disks with the physical conditions discussed pwf1999 (); kohri2002 (); kohri2005 (); chen2007 (); zalamea2011 (); janiuk2010 (); liu2015 (), imply that either and there are no restrictions on and , or if , then the values of and in the shaded regions of the plots in Fig. 4 should be excluded for and , as well as similar regions for intermediate values of . We note that because the IC emission generated by the would not be superimposed directionally and temporally to a normal GRB one, then it can be expected that even for lower levels of heavy axion production (i.e., ), these IC fluxes would have been observed for the aforementioned and values.

In particular, applying these arguments to models with negligible or absent decays to pairs and MeV, we can exclude the region of the space marked in cyan color and labeled as “GRB - ” in Fig. 6. This region overlaps with existing bounds derived from different astrophysical arguments, i.e., excessive energy loss in red giant stars would affect the observed counts in the horizontal branch of color–magnitude diagrams of globular clusters (“ HB”), the duration of the neutrino signal of SN 1987A (“SN”), and an unobserved delayed photon burst due to axion decay if they were produced in SN 1987A (“SN decay”). Our region also overlaps with part of excluded regions by cosmological considerations, i.e., heavy axion decays would have caused distortions in the CMB spectrum when the universe was opaque to photons (“CMB”), and the observed EBL flux cannot be exceeded (“EBL”) by axion decay to photons when the universe became transparent. In addition, the decays to photons in the early universe cannot cause dilution of the neutrino density or affect the abundances of primordial nuclei that are consistent with Big Bang nucleosynthesis (“BBN”). These constraints derived from cosmology, as discussed by (jaeckel2017, ), are model dependent in the sense that they are valid provided that the reheating temperature is relatively high ( GeV) in order to allow significant thermal production of heavy axions. Since there are cosmological models which involve lower reheating temperatures, the bounds based on the cosmology arguments mentioned above (marked in light gray in Fig. 6) would not apply, and then the results of the present work become more useful.

In the case of dominant decays to , the existing bounds in the space are shown in Fig. 7. It can be seen from this plot that the region studied here, marked with cyan color and labeled “GRB-