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 e + e − 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 g a N 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.


Introduction
Gamma-ray bursts (GRBs) are among the most powerful events in the universe since the Big-Bang, emitting radiation at a rate L γ ∼ 10 51−53 erg s −1 , which is observed on Earth as flashes of gamma rays with energies ∼ (10 keV − 50MeV) and durations ranging from a fraction to hundreds of seconds (for a review, see [1]). 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 [2]. Longer GRBs are associated with the collapse of a massive star into a black hole [3]. It is generally accepted that an initial fireball of e + e − , gamma rays, and baryons is generated close to the black hole, which then expands to reach ultrarelativistic velocities [4]. 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 [5,6]. 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 z, thereby confirming the extragalactic origin of GRBs (e.g. see [7,8]).
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 [9,10,11,12,13]. The high densities and temperatures in Email address: mreynoso@mdp.edu.ar (Matías M. Reynoso) this disk as well as in the generated fireball have motivated the study of different particle physics beyond the Standard Model [14,15,16,17,18,19]. In this study, we further explore the possibility that heavy axion-like particles can be emitted from such accretion disks [17,18,20].
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 [23]. As a possible solution to this problem, the standard axion arises as a pseudo-Nambu-Goldstone boson that spontaneously breaks the Peccei-Quinn symmetry (U PQ ) at a scale f a [24]. This cancels out the CP-violating term of the QCD Lagrangian and the axion acquires a mass given by m a 6 × 10 −4 eV 10 10 GeV Heavy axion-like particles 1 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 [26,25,27]. The Lagrangian terms describing the heavy axion interaction with matter are given by [17] L a = − 1 2 m 2 a a 2 − ig ai aψ i γ 5 ψ i − g aγ 4 aF µνF µν , where g ai are the coupling constants for interactions with fermions ψ i , (i = e, N) with N = (p, n), and g aγ parameterizes the coupling to photons. In particular, the most stringent bounds on the coupling constants g aN and g aγ refer to the standard (light) axion, for which both constants are related and relation (1) holds [23,28,29,30]. For instance, bounds on g aγ were obtained based on the duration of SN 1987A [31], i.e. g aγ 10 −6 GeV −1 for m a ∼ 1 − 10 MeV, and a more recent analysis [32] concluded that g aγ 10 −14 GeV −1 for the same range of masses. In addition, studying the production of heavy axions inside a supernova core leads to constraints on g aN , as discussed by [33] in the case where the axions produced can escape freely outside 2 . Further bounds on the coupling to photons g aγ 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 [34,35].
In this study, as in [17], we use a phenomenological approach and maintain all of the coupling constants as independent. We consider heavy axions of mass m a = 0.1 − 10 MeV which can couple to photons, electrons, and nucleons, adopting for the latter values in the range g aN ∼ 10 −7 − 10 −5 , 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 g aN ∼ (1 − 3) × 10 −6 [36,37], and second, the existing bound excluding the range (3 × 10 −9 g aN 10 −6 ) 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 [34]). Hence, it may be interesting to explore further astrophysical phenomena in which such values of g aN can also be tested. Under the assumptions stated above, the leading process for heavy axion production is nucleonnucleon bremsstrahlung [38], 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 g aγ −g ae and g aγ −m a planes for different values of g aN .
In particular, we focus on the production of heavy axions in the accretion disks of GRBs. This possibility was studied by [17,18,20], 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 2 In particular, the results obtained by [20] imply that the mean free path for axions of mass ∼ 1 MeV will be larger than the core radius R SN = 10 km if g aN < 10 −8 for a core density ρ SN ≈ 10 14 g cm −3 and temperature T SN ≈ 30 MeV. approach described by [20], and we then find the ranges of values for the coupling constant g aN for which axions can escape from the disk. If heavy axions decay close to the central engine d a < 10 9 cm, 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 [39,40]) or by neutrinoantineutrino annihilation [9,41]. 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 electronpositron 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 g aN 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.

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 [9,10,11,12,13]. 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 [11,13,42,43,44]. In this study, we expand the model presented previously [42] to include the new process of interest comprising the production of heavy axion-like particles. Similarly to previous studies [43,44,45], we use the notation of Riffert and Herold [46] for the correcting factors to account for general relativistic effects due to the rotating black hole with mass M bh and dimensionless spin parameter a * : where r 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 [46] for details). The accretion rate in the disk is supposed to be constant (∼ 0.1 − 10 M s −1 ) and mass conservation at a radius r from the black hole implies thatṀ where Σ = 2ρH is the mass surface density, ρ is the disk mass density, v r is the radial velocity and H is the half thickness of the disk. The latter can be written as [46] H Pr 3 ρGM bh and it is related to the viscous shear as where α is the Shakura-Sunyaev viscosity coefficient [47]. The total pressure is given by where the fraction of free nuclei is approximated by [12] X nuc ≈ 295.5ρ with ρ 10 = ρ 10 10 g cm −3 and T 11 = T 10 11 K . 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. [11] and adopt the prescription: where the neutrino optical depth is the sum of the scattering plus the absorptive contributions, τ ν l = τ s,ν l + τ a,ν l , 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 The steady-state solutions are obtained requiring that the heating rate is equal to the total cooling rate at each radius, Q − = Q + vis , including all the relevant cooling processes: The rate of photo-disintegration for heavy nuclei is given by [9] Q phot = 10 29 ρ 10 v r dX nuc dr H, (13) and the cooling by advection can be approximated as [48] Q − adv v r r Neutrino cooling occurs mainly through the electron-positron pair capture process, p + e − → n + ν e and n + e + → p +ν e , at a rate and also through electron-positron pair annihilation at rates Considering the corresponding inverse processes, the absorption optical depth can be approximated as and the scattering process as τ s = 2.7 × 10 −7 T 2 11 H. Then, by using a simplified treatment for neutrino emission and transport based on a two-stream approximation (e.g., see [11,49]), we obtain the escaping energy rate in neutrinos as follows: Finally, we include the energy loss rate due to the emission of heavy axions by nucleon-nucleon bremsstrahlung according to [20], without neglecting the finite mass of the axions and pions, where I a (E a ) = dN a dE a is the intensity of the emitted axions, 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 [20] Figure 1 shows the obtained cooling rate Q a as a function of the density for T = 5×10 10 K and T = 10 11 K for g aN = 2×10 −6 and m a = 2m e c 2 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 ρ 3 × 10 11 g cm −3 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 r, we proceed as described by [42] and first solve numerically Eqs. (8) and (9) to obtain H and P as functions of the density ρ. Then, by using of the equation of state (Eq. 10), we can also obtain T 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 (ρ, T ) that satisfies the energy balance at each radius r . Figure 2 shows the results obtained for the profiles ρ(r), T (r), and H(r) in the left, middle, and right panels, respectively. The values of the parameters used are summarized in Table 1. The accretion rate isṀ = 0.1 M s −1 in the top panels,Ṁ = 1 M s −1 in the center panels, andṀ = 3 M s −1 in the bottom panels. In this plot, we assume that the mass of the heavy axions is m a = 2m e . 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 (λ a ) evaluated at their mean energy (which turns out to be E a 2kT [20]), compared with the thickness H(r). Hence, we can conclude that for couplings as high as g aN 10 −5 , heavy axions with mass m a = 2m e will escape freely from the disk ifṀ 0.1M s −1 , whereas for disks with higher accretion ratesṀ = 1−3M s −1 , axions will escape without interacting for couplings g aN 10 −6 . For higher values of g aN , 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 q a = Q − a /Q − tot and neutrinos q ν = Q − ν /Q − tot as functions of the disk radius in the case of an accretion rateṀ = 1M s −1 for a coupling g aN = 2 × 10 −6 , and for heavy axion masses m a = {0.1MeV, 10MeV}. 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 [13,44].

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 [17,18,20], if the axions generated in the disk are to decay close enough (d a < 10 9 cm), the e + e − fireball can be formed more efficiently than by neutrino-antineutrino annihilation. The decay rates of these heavy axions to photons and to electronpositron pairs are so we find that, for instance, for axions of mass m a ∼ 1MeV, the coupling constants must be g aγ 10 −5.1 GeV −1 and/or g ae 10 −8.8 for the decays to occur at distances less than 10 9 cm 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 [1,6]. 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 e + e − scattering on the CMB.
As mentioned in previous studies [17,18], 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 T 0 10 10 K and the density is ρ 0 10 10 g cm −3 (e.g. [21]). If we assume that the density drops on the envelope as ρ e = ρ 0 r r 0 −1.5 (e.g. [22]), then by taking r 0 = 10 4 cm, we find that the optical depth for axions of energies ∼ 1 MeV is much less than one in all of the cases studied: 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 and an analogous expression is used for the neutrino luminosity L ν . In the expression above, R out = 200 r g is the outer radius of the disk and the inner radius is taken as that corresponding to the last stable circular orbit, and z 2 = 3a 2 + z 2 1 1 2 . We can also estimate the power that would be generated by the Blandford-Znajek process using the expression given by [50]: where B in is the poloidal magnetic field near the horizon. The latter can be related to the pressure via B 2 in /(8π) ≈ P in , and the pressure is approximated by (see [51]) P in ≈ 10 [30+1.22 a +log(Ṁ/(M s −1 ))] erg cm −3 . For comparison, Table 2 shows the obtained values for L a , L ν , and L BZ 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 10% in the cases studied of black hole spin and accretion rates, as pointed out by [51]. In general this process is less efficient than the Blandford-Znajek process, which varies here only with the accretion rate because we keep a = 0.9 in all cases. The luminosity in heavy axions clearly increases with g aN 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 [52] φ Band with E 0 = 100keV, and for the low-and high-energy indexes, we take α 0 = −1, and β 0 = 2.1 based on the analysis of the Fermi collaboration using four years of data taking on many bursts with the Gamma-ray Burst Monitor [53]. The constant A can be fixed by normalization on the total flux in a given energy band, and we consider the band (10keV − 40MeV) and a flux  Φ γ ≈ 2.5 ph cm −2 s −1 which according to Ref. [53], 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 d a = E a m a c 2 v a Γ a→γγ , such that 10 17 cm < d a < 10 28 cm for GRBs at redshift z = 1, i.e. at a luminosity distance d L = 3.6 × 10 28 cm. We consider this particular distance for illustration because the distribution of GRB with redshift exhibits its higher values for z ∼ 1 [53]. Given the above expressions for the decay rates, we can see that Γ a→γγ = Γ a→e + e − will hold, e.g., if g ae = (5.8 × 10 −4 GeV)g aγ for m a = 2m e , and if g ae = (3.5 × 10 −3 GeV)g aγ for m a = 10 MeV. Hence, the ranges of the couplings g aγ and g ae that correspond to the channel a → γγ as dominant are marked by the shaded regions in the upper plots in Fig. 4, for m a = 2m e and m a = 10 MeV 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 On the other hand, if the dominant decay channel is a → e + e − , then we consider decay lengths such that 10 20 cm < d a < 3 × 10 26 cm, i.e., up to a maximum of 400Mpc from the burst site, in order to consider that the interactions with the CMB are initiated at z 1 but not closer than ∼ 30 kpc to avoid e ± 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 g ae and g aγ at which these decays are dominant, for m a = 2m e and m a = 10 MeV 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 [54], such that the spectrum of cascade photons initiated by an electron of energy E e can be generally expressed as Here, E γ = m e c 2 /[ ebl (1 + z)], E X = 1 3 E γ /m e c 2 2 cmb (1 + z), and the characteristic energies of the CMB and the EBL are cmb = 6.3 × 10 −4 eV and ebl = 0.68 eV, respectively . We adopt this dichromatic approximation for the background photons, although the e ± produced will interact mainly with the CMB radiation and the scattered photons are cascade sterile, because they are below the threshold for e + e − production. The parameter K can be computed by considering energy conservation between the initial e ± energy and the total energy of the scattered photons. We note that in the present context and for redshifts z 2, the break energy is E X 15 MeV, so that most e ± have energies E e < E X , and hence the scattered photons have a spectrum ∝ E − 3 2 γ . In these cases, it is found that K = 1 2 √ E e /E X , whereas if E e < E e < E γ , then K = E e / [E X (2 + ln(E e /E X ))]. According to [54], if we assume that the interactions occur in less than a timescale ∼ H −1 hubble (z), then we can write the flux of IC scattered photons to be observed on Earth as Here, the intensity of decaying electrons and positrons is given by where γ a = E a m a c 2 is the Lorentz factor of the axion and β a is its velocity in units of c. The minimum energy for the decay to e ± of energy E e is given by which follows when we consider that the maximum and minimum e ± energies in the laboratory frame add up to E a , similarly to the discussion by [55]. Figure 4: Regions in the g aγ − g ae plane that correspond to dominant a → γγ decays in the top panels, and to dominant a → e + e − decays in the lower panels, and for m a = 2m e and m a = 10 MeV on the right and left plots, respectively. The shaded regions can be excluded if 10 −7 g aN 10 −5 , for which axions can be efficiently emitted from GRB accretion disks. The arrows denote that the regions can be extended arbitrarily in the directions indicated. 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 (Φ γ 5 ph cm −2 s −1 ), according to the analysis of Fermi-GBM (see Fig. 11 in [53]). 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 10 −7 g aN 10 −5 , we can exclude the ranges of g aγ and g ae marked by the shaded regions in Fig. 4. Now, in the cases of axions that decay dominantly to e + e − , 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. [54]). Hence, for coupling values 10 −7 g aN 10 −5 and (g ae , g aγ ) 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 g aN , g aγ , and g ae .

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 m a ∼ 1 MeV, their coupling to nucleons can be g aN < 10 −5 and they could leave the disk without interacting for accretion ratesṀ 0.1 M s −1 , whereas for higher accretion rates (Ṁ 1 − 3 M s −1 ), the disk is denser and it is necessary that g aN 10 −6 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 g aN considered.
The luminosity in heavy axions can still be important and it would lead to a more efficient production of photons and/or e + e − than via neutrino-antineutrino annihilation. This can be seen in Table 2, especially for the highest values considered for g aN . Hence, as proposed by [17] but having performed a more realistic calculation of the axion luminosity, if the axions produced decay to photons and/or e + e − at distances d a < 10 9 cm, then the initial fireball could be generated and give rise to the observed GRB. For example, this would be the case if g aγ 10 −5.1 GeV −1 for dominant decays to photons and g ae 10 −8.8 for dominant decays to e + e − , with m a ∼ 2m e . However, the former possibility has been excluded by beam dump and collider experiments [57,63], constraints based on the primordial D/H ratio [35], and the duration of SN 1987A. Figure 6 shows these and the remaining current bounds in the g aγ − m a 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 g aγ so that the typical decay distance of the heavy axions produced is d a→γγ = 10 9 cm, which shows that higher coupling values which would yield d a→γγ < 10 9 cm, but they have already been excluded for the relevant range of heavy axion mass.
On the other hand, the possibility of dominant decays to e + e − to power GRBs has not yet been excluded completely. This can be seen in Fig. 7, where the Boxerino bound extends up to m a 5.6 MeV, and the black dotted line denotes the minimum values of g ae required for decays at distances shorter than 10 9 cm, thereby allowing the possibility that heavier axions decay to e + e to create the initial GRB fireball provided that g ae 10 −7.5 , 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 g ae and m a .
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 g aN taking the values mentioned above and for the couplings g aγ and g ae 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 e ± 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 L γ = 10 52 erg s −1 in the energy band (1 keV − 40 MeV) for a GRB at a redshift z = 1. 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 g aN ∼ 10 −7 .
The lack of observations of the signatures described along with the assumption that GRBs actually involve accretion disks with the physical conditions discussed [9,10,12,13,41,43,51], imply that either g aN 10 −7 and there are no restrictions on g aγ and g ae , or if 10 −7 g aN 10 −5 , then the values of g aγ and g ae in the shaded regions of the plots in Fig. 4 should be excluded for m a = 2m e and m a = 10 MeV, as well as similar regions for intermediate values of m a . We note that because the IC emission generated by the e ± would not be superimposed directionally and temporally to a normal GRB one, then it can Figure 5: Contributions to the flux of gamma rays due to the decay of heavy axions produced in the accretion disk of a GRB at z = 1, and for accretion rateṡ M = 1M s −1 in the left panels andṀ = 3M s −1 in the right panels. The mass of the heavy axions m a is 0.1 MeV, 2m e , and 10 MeV in the top, middle and bottom panels, respectively. Blue lines correspond to φ a→γγ (E γ ) and red lines to φ IC e ± (E γ ), and a typical GRB flux is shown by dashed gray curves. Figure 6: Excluded regions in the m a − g aγ space adapted and combined from previous studies [36,37,34,32], including the region explored in this work based on heavy axion decays to photons outside the GRB. The latter region, marked with in a transparent cyan color and labeled as "GRBg aN ∈ (10 −7 , 10 −5 )" can be ruled out for g aN within the range of values indicated. The brown and light brown regions have been excluded by collider and beam dump experiments [57,63]. The orange region is excluded based on searches of axion-photon conversions of solar axions with the "CAST + Sumico" helioscope. The magenta region ("x ion ") is ruled out since heavy axion decays would lead to a too early reionization of the universe [56], while the dark gray region ("DM") is excluded because heavy axion dark matter would be excessively produced, as well as the light gray region labeled "X-rays" because heavy axion decays inside galaxies would imply unobserved X-ray features [34]. The region marked by black vertical lines corresponds to values of the standard "axion models" and they are shown for illustration. The dotted black line indicates the minimum values of g aγ required for heavy axions to power GRBs. The remaining labels are described in Section 4. Figure 7: Excluded regions in the m a − g ae space adapted and combined from previous studies [60,61,62,63], including the region explored in this work based on heavy axion decays to electron-positrion pairs outside the GRB. The latter region, marked in cyan color and labeled as "GRBg aN ∈ (10 −7 , 10 −5 )", can be ruled out for g aN within the range of values indicated. The orange region has been excluded by searches using the Boxerino detector [64], while the light brown region is excluded by beam dump experiments [57]. If galactic dark matter is composed enterely of heavy axions, then the gray region is excluded with the upper bounds given by the experiments indicated. We also include the region in transparent yellow, which is to be probed at LHC Run-2 searching for possible signals of Z → a γ → e + e − decays [63]. The dotted black line indicates the minimum values of g ae required for heavy axions to power GRBs.
be expected that even for lower levels of heavy axion production (i.e., g aN < 10 −7 ), these IC fluxes would have been observed for the aforementioned g aγ and g ae values.
In particular, applying these arguments to models with negligible or absent decays to e + e − pairs and m a ∼ 0.01 − 10 MeV, we can exclude the region of the g aγ − m a space marked in cyan color and labeled as "GRBg aN ∈ (10 −7 , 10 −5 )" 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 [32], are model dependent in the sense that they are valid provided that the reheating temperature is relatively high (T RH 120 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 e + e − , the existing bounds in the m a − g ae 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-g aN ∈ (10 −7 , 10 −5 )", has not been ruled out previously. We find that for g aN ∈ (10 −7 , 10 −5 ), this region should be excluded because it implies unobserved gamma ray flashes with an energy dependence that is clearly different from that corresponding to GRBs (i.e., the red lines in Fig. 5).
To summarize our conclusions, we have seen that the production of heavy axions in GRB accretion disks can lead to the formation of the GRB fireball via the decay channel a → e + e − , because an allowed region in the m a − g ae space is still compatible with this situation (Fig. 7). By constrast, heavy axions decaying to photons cannot be the origin of the fireball, because the required values for g aγ have been excluded (Fig. 6). Finally, we have found other combinations of values for the couplings g aN , g aγ , and g ae , which as discussed above, are inconsistent with the current understanding and observations of GRBs. This can be used as a complementary tool to constrain models that involve heavy axion-like particles.