© The Authors 2023

1. Introduction

Since 2014, neutrinos of astrophysical origin have been detected by IceCube in the energy range 10⁴ GeV ≲ E_ν ≲ 10⁶ GeV (IceCube Collaboration 2014). However, the sources corresponding most of these neutrinos could not be directly associated with any known sources, such as the typically proposed active galactic nuclei (AGN; e.g., Stecker et al. 1991), gamma-ray bursts (GRBs; e.g., Waxman & Bahcall 1997), tidal disruption events (TDEs; e.g., Dai & Fang 2017), and starburst galaxies (e.g., Loeb & Waxman 2006). The few exceptions are as follows: First, the neutrino event IC-170922A was associated with the blazar TXS 0506+056, which was almost simultaneously detected in gamma rays by Fermi, VERITAS, and MAGIC (IceCube Collaboration 2018). Second, another AGN, the Seyfert galaxy NGC 1028, was recently detected by IceCube without a gamma-ray counterpart (IceCube Collaboration 2022). And third, the TDEs AT2019dsg (Stein et al. 2021) and AT2019fdr (Reusch et al. 2022) were detected in the optical/UV bands by the Zwicky Transient Facility, and they are very plausibly the IceCube counterparts of the muon track events IC191001A and IC200530A, respectively.

Still, the origin of a great majority of the observed astrophysical neutrinos cannot be identified yet, and it can only be recognized that the sources should be extragalactic, given that the incoming directions in the sky are consistent with an isotropic emission. Among the candidate sources, choked gamma-ray bursts (CGRBs) arise as an interesting possibility. A CGRB is generated when a jet launched inside a collapsing massive star fails to emerge from the stellar mantle, yielding no gamma-ray counterpart. This phenomenon is expected to take place in some core-collapse supernovae (SNe), such as those of the types II (MacFadyen et al. 2001) and Ib/c. Although the fraction of such sources presenting jets is still uncertain (Piran et al. 2019), the link between SNe and CGRBs is conjectured in a similar manner to the case of the confirmed connection between type Ib/c SNe and long GRBs with flat spectrum after the observation of several SNe at the same location previously detected as long GRB (Hjorth & Bloom 2012).

Similar to the proposed scenarios of neutrino production in GRB jets (e.g., Waxman & Bahcall 1997; Murase & Nagataki 2006; Petropoulou et al. 2014), the study of CGRBs as hidden neutrino sources has been realized under different considerations (Mészáros & Waxman 2001; Razzaque et al. 2004; Murase & Ioka 2013; Senno et al. 2016; He et al. 2018; Fasano et al. 2021). Recently, Chang et al. (2022) performed a statistical analysis using data from observed SNe Ib/c, and although they found no significant correlation with IceCube data regarding muon track events, they concluded that choked jets in SNe can still be the main contributors to the diffuse neutrino flux.

We therefore have the relevant physical scenario well-established and we explore it to study the flavor composition of the emitted neutrinos from CGRBs in more detail and with the aim of obtaining specific signatures consistent with their assumed description. In general, it has been recognized that analyzing the flavor composition of the neutrino fluxes can be useful (Beacom et al. 2003; Bustamante et al. 2015), and particularly in highly magnetized sources, the synchrotron cooling of muons can gradually modify the emitted flavor ratio (f_e : f_μ : f_τ)_s from (1 : 2 : 0) to (0 : 1 : 0) at high energies (Lipari et al. 2007; Baerwald et al. 2012). Along this line, a detailed statistical analysis was performed in Bustamante & Tamborra (2020) using the available neutrino data along with a generic description of the neutrino sources, characterized by their comoving magnetic field B and Lorentz factor Γ. In Fiorillo et al. (2021), they also studied the effect of different distributions of target photons and the impact on the flavor ratios for different benchmark parameters describing typical sources.

In this work, we focus on the particular case of CGRBs and consider different aspects that we believe can play a role and should be taken into account in the neutrino production model. First, we solve a stationary transport equation for primary electrons as well as for protons, which are both supposed to be injected in the region with internal shocks of the outflow. The synchrotron emission of these electrons is an additional contribution to the soft photon background that serves as a target for the high energy particles. The other contribution is the commonly assumed one, i.e., the thermalized emission of another electron population, which is accelerated by the reverse shocks at the jet head. A fraction of these photons is considered to escape to the internal shock region and also serve as targets for the high energy particles there. Specifically, the possible interactions with the background photons include: pγ interactions (which produce pions and kaons) and inverse Compton (IC) interactions with electrons as well as with muons from the decay of pions. Taking these cooling processes into account, plus the synchrotron one, we also solve a transport equation to obtain the distributions of pions, muons, and kaons. This allows us to finally compute the neutrino emission consistently with the cooling processes affecting all the particle species involved in the production process. Hence, the generated neutrino flux is sensitive to important physical parameters, such as the jet power, magnetic field, and proton-to-electron ratio. Although the values of the latter, as well as that of the energy budget (Eichler et al. 2010), are uncertain, we aim to explore plausible combinations of them in order to determine their effect on the neutrino emission and flavor ratios.

In particular, for energies above ∼10⁵ − 10⁶ GeV, the electron flavor ratio to be observed on Earth gradually decreases, the muon flavor ratio increases, and the tau one increases only slightly. Experimental data from present and future detectors will therefore help in determining whether CGRBs, as they are normally conceived, can significantly contribute to the diffuse neutrino flux.

The rest of the manuscript is organized as follows. In Sect. 2, we present the basic scenario of the model adopted, and in Sect. 3 we describe the method of calculation used to characterize the particle distributions and neutrino emissivity. In Sect. 4, we present our results for the neutrino flux and flavor ratios as a function of the energy and also integrated over the energy. Finally, in Sect. 5, we conclude with a discussion.

2. Basic scenario

In this section, we briefly describe the model applied, which is based in the same set of assumptions adopted for CGRBs in several works (e.g., Mészáros & Waxman 2001; Murase & Ioka 2013; Senno et al. 2016; Fasano et al. 2021). In particular, our basic scenario is the same as in Fasano et al. (2021), since we assume similar values for the jet power, Lorentz factor, and the density of the stellar envelope. The main differences arise in the treatment applied to describe the high-energy particles involved in the neutrino production process, which in our case implies the use of a transport equation to account for the cooling processes, including πγ and muon IC interactions. Additionally, as mentioned above, we also include high energy electrons as a population of particles that is co-accelerated with protons by the internal shock.

Typically, the jet of a CGRB has a Lorentz factor Γ ∼ 100 − 300 and a half-opening angle θ_op ≈ 0.2 rad, and it carries a power L_j, which corresponds to an isotropic equivalent power

$$\begin{array}{c}\hfill L_0=2L_{\mathrm{j}}/(1-cos\theta )={10}^{50-51}\mathrm{erg}\mathrm{s}\mathrm{C}-1.\end{array}$$(1)

Fixing a variability timescale δt, the distance from the central black hole (BH) to the position where internal shocks develop is

$$\begin{array}{c}\hfill r_{\mathrm{is}}=2{\mathrm{\Gamma }}^{2}c\delta t\simeq 6\times {10}^{12}\mathrm{cm}\delta t_{-2}{\mathrm{\Gamma }}_2,\end{array}$$(2)

where Γ₂ = Γ/100 and δt₋₂ = δt/(0.01 s). There, the comoving number density of cold protons is

$$\begin{array}{c}\hfill n_{\mathrm{j}}^{\prime }=\frac{L_0}{4\pi {\mathrm{\Gamma }}^2{r}_{\mathrm{is}}^2{m}_{\mathrm{p}}{c}^3}\simeq 4.9\times {10}^{12}{\mathrm{cm}}^{-3}{L}_{0,51}{\mathrm{\Gamma }}_2^{-6}\delta t_{-2}^{-2}.\end{array}$$(3)

The jet is magnetized, and the magnetic energy density is usually taken to be a fraction ϵ_B of the kinetic energy density $m_{\text{p}}{c}^2{n}_{\text{j}}^{\prime }$, that is:

$$\begin{array}{c}\hfill B^{\prime }=\sqrt{\frac{2{ϵ}_{\mathrm{B}}{L}_0}{{\mathrm{\Gamma }}^2{r}_{\mathrm{is}}^{2}c}}=1.36\times {10}^5\mathrm{G}{ϵ}_{\mathrm{B},-1}^{1/2}{L}_{0,51}^{1/2}{\mathrm{\Gamma }}_2^{-2}\delta t_{-2}^{-1}.\end{array}$$(4)

We also consider the case of equipartition, ϵ_B = 1, to explore a situation with a higher magnetic field and the consequences for neutrino production.

As the jet propagates through the massive star, inside the extended hydrogen envelope up to a radius r_ext ≈ 10^13.5 cm, a forward shock and a reverse shock develop, and the jet can be stalled. This has been studied through simulations and analytically by Bromberg et al. (2011) and Mizuta & Ioka (2013). Specifically, the part of the jet that is affected by the reverse shock is usually called the “jet head”, and in the present context, it has been found that the Lorentz factor of this region is Γ_h ≃ 1 (e.g., He et al. 2018), and its position is given by (e.g., Mizuta & Ioka 2013; Murase & Ioka 2013)

$$\begin{array}{c}\hfill r_{\mathrm{h}}\simeq 1.3\times {10}^{13}\mathrm{cm}{t}_{\mathrm{j},3}^{3/5}{\left(\frac{L_0}{{10}^{50}{\mathrm{erg}\mathrm{s}}^{-1}}\right)}^{1/5}{\left(\frac{{\theta }_{\mathrm{op}}}{0.2}\right)}^{-4/5}{\rho }_{\mathrm{ext},-7}^{-1/5}.\end{array}$$(5)

Here, the lifetime of the jet is t_j, with t_j, 3 = t_j/(1000 s), and the density of the envelope ρ_ext at r_ext appears scaled with respect to the typical value, that is, ρ_ext, −7 = ρ_ext/(10⁻⁷g cm⁻³). Before reaching the jet head, the jet is expected to stop expanding and become collimated, keeping a cylindrical shape up to the position of the jet head r_h. In these conditions, a collimation shock is generated at a distance from the BH (Bromberg et al. 2011; Mizuta & Ioka 2013)

$$\begin{array}{c}\hfill r_{\mathrm{cs}}\simeq 8.2\times {10}^{12}\mathrm{cm}{t}_{\mathrm{j},3}^{2/5}{\left(\frac{L_0}{{10}^{50}{\mathrm{erg}\mathrm{s}}^{-1}}\right)}^{3/10}{\rho }_{\mathrm{ext},-7}^{-3/10}.\end{array}$$(6)

It is then assumed that a population of electrons in the shocked jet head emits synchrotron and IC radiation that thermalizes due to a high optical depth (e.g., He et al. 2018),

$$\begin{array}{cc}& {\tau }_{\mathrm{h}}\approx (4{\mathrm{\Gamma }}_{\mathrm{rel}}+3)n_{\mathrm{j}}(r_{\mathrm{h}}){\sigma }_{\mathrm{T}}\hfill \\ \hfill & \simeq 1170{\mathrm{\Gamma }}_2^{-2}{L}_{0,51}^{3/5}{\rho }_{\mathrm{ext},-7}^{2/5}{r}_{\mathrm{ext},13.5}^{5/2}{t}_{\mathrm{j},3}^{3/5},\hfill \end{array}$$(7)

where Γ_rel ≈ Γ/(2Γ_h) is the Lorentz factor of the jet head with respect to the jet. Specifically, a fraction ϵ_e ≈ 0.1 of the kinetic energy of the flow in the jet head is supposed to be carried by electrons so that the total radiation density of the thermalized emission is (e.g., Mészáros & Waxman 2001; He et al. 2018):

$$\begin{array}{c}\hfill a{T_{\mathrm{h}}}^4={ϵ}_{\mathrm{e}}4({\mathrm{\Gamma }}_{\mathrm{rel}}+3)({\mathrm{\Gamma }}_{\mathrm{rel}}-1)n_{\mathrm{j}}{m}_{\mathrm{p}}{c}^2,\end{array}$$(8)

which implies that

$$\begin{array}{c}\hfill T_{\mathrm{h}}={\left[{ϵ}_{\mathrm{e}}({\mathrm{\Gamma }}_{\mathrm{rel}}+3)({\mathrm{\Gamma }}_{\mathrm{rel}}-1)\left(\frac{L_{\mathrm{iso}}}{a\pi {\mathrm{\Gamma }}^2{r}_{\mathrm{h}}^{2}c}\right)\right]}^{\frac{1}{4}}.\end{array}$$(9)

Hence, the distribution of these photons in the jet head frame is

$$\begin{array}{c}\hfill N_{\mathrm{ph},\mathrm{h}}(E_{\mathrm{ph},\mathrm{h}})=\frac{2E_{\mathrm{ph},\mathrm{h}}^2}{{(hc)}^3[exp\left(\frac{E_{\mathrm{ph},\mathrm{h}}}{k{T}_{\gamma }}\right)-1]}[\mathrm{e}r{g}^{-1}c{m}^{-3}s{r}^{-1}],\end{array}$$(10)

and it is normally assumed that a fraction f_esc = 1/τ_h of them escape to the internal shock region. However, instead of making the approximation that all the escaping photons carry the same boosted energy ≈Γ_rel(2.8k_BT) in the internal shock frame, we transform the photon distribution to the internal shock frame and integrate on a solid angle element, as we describe in the next paragraph.

In the internal shock frame, the photon energy E_ph is such that E_ph, h = Γ_relE_ph(1 − β_relx), where x is the cosine of the angle between the relative velocity of the flow, β_relc, and the photon momentum in the internal shock frame. Considering the Lorentz invariant (Dermer & Schlickeiser 2002),

$$\begin{array}{c}\hfill \frac{1}{E\sqrt{E^2-m^{2}c}}\frac{\mathrm{d}\mathcal{N}}{\mathrm{d}V\mathrm{d}E\mathrm{d}\mathrm{\Omega }},\end{array}$$(11)

we can obtain the distribution of photons in the internal shock frame as

$$\begin{array}{c}\hfill N_{\mathrm{ph}}(E_{\mathrm{ph}},x)=\frac{N_{\mathrm{ph},\mathrm{h}}({\mathrm{\Gamma }}_{\mathrm{rel}}{E}_{\mathrm{ph}}(1-{\beta }_{\mathrm{rel}}x))}{{\mathrm{\Gamma }}_{\mathrm{rel}}^2{(1-{\beta }_{\mathrm{rel}}x)}^2}[\mathrm{e}r{g}^{-1}c{m}^{-3}s{r}^{-1}].\end{array}$$(12)

We then obtain the differential photon density as

$$\begin{array}{c}\hfill n_{\mathrm{ph}}(E_{\mathrm{ph}})=2\pi {\displaystyle {\int }_{-1}^{x_{\mathrm{is}}}\mathrm{d}x{N}_{\mathrm{ph}}(E_{\mathrm{ph}},x)[\mathrm{e}r{g}^{-1}c{m}^{-3}],}\end{array}$$(13)

where we perform the integration on the solid angle subtended by the internal shock region with respect to the jet head position (where the thermal photons are produced), ΔΩ_is = 2π(x_is + 1). Here, x_is = cos(π − ξ), with the angle ξ being such that (see Fig. 1)

$$\begin{array}{c}\hfill tan\xi =\frac{r_{\mathrm{is}}tan{\theta }_{\mathrm{op}}}{(r_{\mathrm{h}}-r_{\mathrm{is}})/\mathrm{\Gamma }}.\end{array}$$

Fig. 1. Schematic view (not to scale) of the jet propagating inside a star.

Therefore, the differential photon density of Eq. (13) characterizes the external photon background that could be relevant, in particular, for the neutrino producing pγ interactions of the protons accelerated by the internal shocks.

3. High energy particles and their interactions

Particle acceleration is normally expected to take place at the internal shocks mentioned above, provided they do not become radiation dominated. As discussed by Murase & Ioka (2013), this can be satisfied if the comoving size of the shocked region r_is/Γ is smaller than the photon interaction length $1/(n_{\text{j}}^{\prime }{\sigma }_{\text{T}})$, that is, if

$$\begin{array}{c}\hfill {\tau }_{\mathrm{is}}=\frac{r_{\mathrm{is}}}{\mathrm{\Gamma }}{n}_{\mathrm{j}}^{\prime }{\sigma }_{\mathrm{T}}\simeq 0.196L_{0,51}{\mathrm{\Gamma }}_2^{-5}\delta t_{-2}^{-1}<1,\end{array}$$(14)

and this condition is fulfilled in the cases considered in the present work, namely, for Γ = 100 and Γ = 300, as can be seen in the left panel of Fig. 2. In the figure, we plot the optical depth as a function of Γ in the cases of L₀ = 10⁵⁰ erg s⁻¹ and L₀ = 10⁵¹ erg s⁻¹ for both δt = 0.001 s and δt = 0.01 s. In particular, for the typically adopted value Γ = 100, the optical depth is less than one for both L₀ = 10⁵⁰ erg s⁻¹ and L₀ = 10⁵¹ erg s⁻¹, and the variability timescale is δt = 0.01 s. If this timescale is shorter, such as δt = 0.001 s (e.g., Senno et al. 2016), and setting L₀ = 10⁵¹ erg s⁻¹, we obtain low optical depths for higher values of Γ.

Fig. 2. Left panel: optical depth at the internal shock region as a function of the Lorentz factor of the jet for (L0, δt)=(1050 erg s−1, 0.001 s) and (L0, δt)=(1051 erg s−1, 0.01 s), both marked by the black curve, while the blue line corresponds to (L0, δt)=(1050 erg s−1, 0.01 s) and the orange one to (L0, δt)=(1051 erg s−1, 0.001 s). Right panel: positions of the jet head (rh) as a function of the Lorentz factor of the jet for L0 = 1050 erg s−1 and L0 = 1051 erg s−1 marked by long-dashed blue and orange lines, respectively. Solid blue and orange lines correspond to the positions of collimation shock (rcs) for L0 = 1050 erg s−1 and L0 = 1051 erg s−1, respectively. The positions of the internal shocks (ris) for δt = 0.001 s and δt = 0.01 s are indicated with black short-dashed and solid lines, respectively.

As an additional condition to avoid radiation dominated shocks, we required that the position where the particles are accelerated by internal shocks be placed inside the collimation shock, that is, r_is < r_cs, as the post-shock region would otherwise be dominated by radiation (Murase & Ioka 2013). As can be seen in the right panel of Fig. 2, this condition can be satisfied for a range of values of Γ that depend on jet power and other parameters, such as envelope density and jet duration. Adopting typical values for these parameters (e.g., Fasano et al. 2021), t_j = t_j, 3, ρ_ext = ρ_ext, −7, we observe that the internal shock region is inside the collimation shock, in particular for Γ = 100 and δt = 0.01 s as well as for Γ = 300 along with δt = 0.001 s. Hence, we retain these two combinations of basic parameters, which are consistent with an efficient particle acceleration, and we explore the impact of different values of other key physical parameters, such as jet power, proton-to-electron ratio, and jet magnetization, on neutrino production.

We therefore assume the particle acceleration rate to be

$$\begin{array}{c}\hfill t_{\mathrm{acc}}^{-1}(E_i^{\prime })=\eta \frac{eBc}{E_i^{\prime }},\end{array}$$(15)

where i = {e, p} stands for the primary electrons and protons, respectively, and we consider the acceleration efficiency coefficient to take a high value, namely, η = 0.1. In addition, we assume that the accelerated particles carried a power L_p + L_e = ϵ_relL_k, that is, a fraction ϵ_rel ≈ 0.1 of the total kinetic power carried by the jet. As mentioned above, we explore two different possibilities for the proton-to-ratio a_ep = L_p/L_e: an equal power injected for both protons and electrons (a_ep = 1) and a situation similar to what is inferred for the cosmic-ray flux arriving on Earth (a_ep = 100).

The key parameters of the model and the values assumed in the present work are shown in Table 1. We remark that given the adopted typical duration of a CGRB event, t_j = 1000 s, a higher energy budget is clearly required in the case of L₀ = 10⁵¹ erg s⁻¹. However, such conditions can still be attained, for instance, in blue supergiant (BSG) progenitors, where the corresponding jet breakout time can be estimated as (Murase & Ioka 2013)

$$\begin{array}{c}\hfill t_{\mathrm{bo}}\simeq 9800\mathrm{s}{L}_{0,51}^{-1/3}{\left(\frac{{\theta }_{\mathrm{op}}}{0.2}\right)}^{2/3}\left(\frac{R_{\ast }}{{10}^{13.5}\mathrm{cm}}\right).\end{array}$$(16)

Since this value is well above the assumed jet duration time, the jet would effectively be choked.

We proceed, then, by computing the injection of primary particles in the comoving frame as

$$\begin{array}{c}\hfill Q_i^{\prime }(E_i^{\prime })=\frac{\mathrm{d}{\mathcal{N}}_i}{\mathrm{d}{E}_i^{\prime }\mathrm{d}{\mathrm{\Omega }}^{\prime }\mathrm{d}{V}^{\prime }\mathrm{d}{t}^{\prime }}=K_i{E}_i^{\prime }exp(-\frac{E_i^{\prime }}{E_{i,\max }^{\prime }}),\end{array}$$(17)

where $E_{i,\max }^{\prime }$ is a maximum energy obtained by the balance of the acceleration rate with the energy loss rate, which we describe below. The normalizing constant K_i is therefore obtained from the following relation:

$$\begin{array}{c}\hfill \mathrm{\Delta }V{\displaystyle {\int }_{4\pi }\mathrm{d}\mathrm{\Omega }{\int }_{E_{i,\min }}^{\infty }\mathrm{d}{E}_i{E}_i{Q}_i(E_i)=L_i,}\end{array}$$(18)

where $\mathrm{\Delta }V=4\pi r_{\text{is}}^2\mathrm{\Delta }d$ is the volume of the region with internal shocks, and Δd ≃ cδt_v is the corresponding thickness in the central BH rest frame. The injection in this frame is obtained by taking into account that

$$\begin{array}{c}\hfill \frac{E_i}{p^2}\frac{\mathrm{d}{\mathcal{N}}_i}{\mathrm{d}V\mathrm{d}p\mathrm{d}\mathrm{\Omega }dt}\propto \frac{Q_i}{\sqrt{E_i^2-m_i^2{c}^4}}\end{array}$$

is a Lorentz invariant (Dermer & Schlickeiser 2002) and that the energy in the comoving frame is related to the one in the BH frame by the usual Lorentz transformation $E_i^{\prime }=\mathrm{\Gamma }(E_i-\beta \mu \sqrt{E_i^2-m_i^2})$, where μ = cos θ_i is the cosine of the angle between the particle momentum and the jet velocity in the BH frame. As for the minimum energy of the accelerated particles, our results are not very sensitive to its exact value, but we fixed it at $E_{i,\min }^{\prime }=2m_i{c}^2$ in the comoving frame, as this accounts for the fact that diffusive shock acceleration is effective only for suprathermal particles, that is, particles with energies well above those corresponding to the mean energy of the thermal distribution (Bosch-Ramon et al. 2006).

We obtain the particle distribution by solving the steady-state transport equation

$$\begin{array}{c}\hfill \frac{\mathrm{d}\left[b_{i,\mathrm{loss}}^{\prime }(E_i^{\prime })N_i^{\prime }(E_i^{\prime })\right]}{\mathrm{d}{E}_i^{\prime }}+\frac{N_i^{\prime }(E_i^{\prime })}{T_{i,\mathrm{esc}}}=Q_i^{\prime }(E_i^{\prime }),\end{array}$$(19)

where the $b{\prime }_{i,\mathrm{loss}}$ is the total energy loss for particles of the type i, which is related to the sum of the cooling rates $t_{i,j}^{\prime -1}$ corresponding to all the processes j as

$$\begin{array}{c}\hfill b_{i,\mathrm{loss}}^{\prime }={\frac{\mathrm{d}{E}_i^{\prime }}{\mathrm{d}{t}^{\prime }}|}_{\mathrm{loss}}\equiv E_i^{\prime }{\displaystyle \sum _j{t}_{i,j}^{\prime -1}(E_i^{\prime }).}\end{array}$$(20)

The escape timescale is conservatively considered to be the Bohm diffusion time corresponding to the size Δd′, such that

$$\begin{array}{cc}& T_{\mathrm{esc}}^{-1}=\frac{2E_i^{\prime }}{3{(\mathrm{\Delta }{d}^{\prime })}^2(e{B}^{\prime }c)}\hfill \\ \hfill & \simeq 1.7\times {10}^{-5}{\mathrm{s}}^{-1}\left(\frac{E_i^{\prime }}{{10}^4\mathrm{GeV}}\right){\left(\frac{{ϵ}_{\mathrm{B}}}{0.01}\right)}^{-\frac{1}{2}}{\left(\frac{L_0}{{10}^{51}\mathrm{erg}{\mathrm{s}}^{-1}}\right)}^{-\frac{1}{2}}\delta t_2^{-1}.\hfill \end{array}$$(21)

This escape rate is generally well below the cooling rates involved, which are presented in the next subsection.

The solution of Eq. (19) is formally given by

$$\begin{array}{cc}& N_i^{\prime }(E_i^{\prime })=\frac{1}{b_{i,\mathrm{loss}}(E_i^{\prime })}{\displaystyle {\int }_{E_i^{\prime }}^{\infty }\mathrm{d}{E}^{\prime }{Q}_i^{\prime }(E^{\prime })}\hfill \\ \hfill & \times exp{\displaystyle [-{\int }_{E_i^{\prime }}^{E^{\prime }}\frac{\mathrm{d}E″}{T_{\mathrm{esc}}(E″)b_{i,\mathrm{loss}}(E″)}]}\hfill \end{array}$$(22)

and is applied to the different particle species (e,  p,  π^±,  μ^±) while taking into account the corresponding cooling processes and injections in each case.

3.1. Cooling processes

In this section, we discuss the interactions and corresponding cooling rates for the high energy charged particles at the internal shock region. As discussed above, there is background of external soft photons (i.e., produced outside the internal shock region). These photons, characterized by the differential density of Eq. (13), enter the internal shock zone and can serve as targets for IC and pγ interactions. In this work, we also account for the possible interaction of high energy particles with synchrotron photons emitted by electrons in the internal shock zone itself.

In general, for any charged particle of mass m_i, the synchrotron cooling rate is given by

$$\begin{array}{c}\hfill t_{i,\mathrm{syn}}^{-1}(E_i^{\prime })=\frac{4}{3}{\left(\frac{m_{\mathrm{e}}}{m_i}\right)}^3\frac{{\sigma }_{\mathrm{T}}{B}^2}{m_{\mathrm{e}}c8\pi }\frac{E_i^{\prime }}{m_i{c}^2}.\end{array}$$(23)

We apply this equation to electrons, protons, pions, and muons in order to obtain their corresponding particle distribution in the internal shock zone.

The differential density of synchrotron photons is computed as

$$\begin{array}{c}\hfill n_{\mathrm{ph},\mathrm{syn}}=\frac{{ϵ}_{\mathrm{syn}}(E_{\mathrm{ph}})}{E_{\mathrm{ph}}}\frac{\mathrm{\Delta }{d}^{\prime }}{c},\end{array}$$(24)

where Δd′=Δd/Γ is the comoving thickness of the internal shock region, and ϵ_syn is the power per unit volume per unit energy of the photons,

$$\begin{array}{c}\hfill {ϵ}_{\mathrm{syn}}(E_{\mathrm{ph}})=\left(\frac{1-e^{-{\tau }_{\mathrm{SSA}}(E_{\mathrm{ph}}))}}{{\tau }_{\mathrm{SSA}}(E_{\mathrm{ph}})}\right){\displaystyle {\int }_{m_{\mathrm{e}}{c}^2}^{\infty }\mathrm{d}E4\pi P_{\mathrm{syn}}{N}_{\mathrm{e}}^{\prime }(E).}\end{array}$$(25)

Here, P_syn(E_ph) is defined as in Blumenthal & Gould (1970), and τ_SSA stands for the optical depth given by Rybicki & Lightman (2004).

The synchrotron emission of electrons in the internal shock region can provide an additional target for IC interactions of the electrons themselves, that is, the so-called synchrotron self-Compton (SSC) interactions. We also note that IC interactions of muons can be relevant, since they provide an additional cooling process (which may become even more significant than the synchrotron one) that decreases the neutrino output at high energies, as we show below in the next subsection. The IC cooling rate for a charged lepton of mass m_l is computed following the standard expressions of Blumenthal & Gould (1970), which can account for both the Thomson and the Klein-Nishina regimes:

$$\begin{array}{cc}& t_{\mathrm{IC}}^{-1}(E_l^{\prime })=\frac{3m_l^2{c}^4{\sigma }_{\mathrm{T}}}{4E^3}{\displaystyle {\int }_{E_{\mathrm{ph}}^{(\min )}}^E\frac{\mathrm{d}{E}_{\mathrm{ph}}}{E_{\mathrm{ph}}}{n}_{\mathrm{ph}}(E_{\mathrm{ph}})}\hfill \\ \hfill & \times {\displaystyle {\int }_{E_{\mathrm{ph}}}^{\frac{{\mathrm{\Gamma }}_l}{{\mathrm{\Gamma }}_l+1}E}\mathrm{d}{E}_{\gamma }{F}_{\mathrm{IC}}(q)[E_{\gamma }-E_{\mathrm{ph}}],}\hfill \end{array}$$(26)

where $E_{\mathrm{ph}}^{(\min )}$ is the lowest energy of the available background of target photons, q = E_γ(Γ_l(E_l − E_γ)), with ${\mathrm{\Gamma }}_l=4E_{\text{ph}}{E}_l/(m_l^2{c}^4)$, and the function F_IC(q) is given by

$$\begin{array}{c}\hfill F_{\mathrm{IC}}(q)=2qlnq+(1+2q)(1-q)+\frac{1}{2}(1-q)\frac{{(q{\mathrm{\Gamma }}_l)}^2}{1+{\mathrm{\Gamma }}_{l}q}.\end{array}$$(27)

In Fig. 3, we present the obtained acceleration and cooling rates for the injected primary electrons in the cases of L₀ = 10⁵⁰erg s⁻¹ in the left panels and of L₀ = 10⁵¹erg s⁻¹ in the right panels. The top panels correspond to Γ = 100 and the bottom ones to Γ = 300. Except for the IC cooling rate due to the external photons, the rest of the rates are sensitive to the magnetic field, for which we consider three different values determined by fixing the parameter ϵ_B at 0.01, 0.1, and 1.

Fig. 3. Cooling and acceleration rates for electrons at the internal shocks of CGRBs for L0 = 1050 erg s−1 and L0 = 1051 erg s−1 in the left and right panels, respectively. Top panels correspond to Γ = 100 and the bottom ones to Γ = 300. The different processes are indicated with the following curves: Dashed lines are for ϵB = 1, solid lines for ϵB = 0.1, and short-dashed lines for ϵB = 0.01. The blue curves refer to acceleration, orange ones to synchrotron, and green ones to external IC. The circles refer to SSC with (aep = 1, ϵB = 0.1) and squares to SSC with (aep = 100, ϵB = 0.1). With dashed gray lines, we also show the auxiliary results used to obtain the SSC rate with the iterative method mentioned in the text.

As for the SSC process, since the density of the target synchrotron photons depends on the electron distribution, the SSC cooling rate does as well. It is therefore a non-trivial issue to address this process self-consistently through a single expression for the cooling rate (e.g., Schlickeiser 2009, in the Thomson regime). To handle this situation, we obtain a first approximation of the electron distribution $N_{\text{e},\text{1}}^{\prime }$ while assuming a nil value for the SSC cooling rate. With $N_{\text{e},\text{1}}^{\prime }$ we compute a first approximation of the SSC cooling rate $t_{\mathrm{SSC},1}^{-1}$ using Eq. (26), and we use it to obtain a second approximation for the electron distribution, $N_{\text{e},\text{2}}^{\prime }$, and then $t_{\mathrm{SSC},2}^{-1}$. We repeatd this iterative method several times until the ratio between $t_{\mathrm{SSC},\mathrm{j}}^{-1}$ and $t_{\mathrm{SSC},\mathrm{j}-1}^{-1}$ was reduced to less than a few percent at low energies. For illustration of this, we include in Fig. 3 our successive approximations of the SSC cooling rate obtained in the case of ϵ_B = 0.1, using short-dashed gray lines for a_ep = 1 and long-dashed gray lines for a_ep = 100. As can be seen in the figure, for a higher injected power in electrons and Γ = 100 (i.e., in the top-right panel), the SSC becomes dominant as compared to the external IC rate, and more iterations are necessary, particularly for a_ep = 1. In the cases of Γ = 300, the external IC rate is higher than for Γ = 100. This is due to an increased boosting effect of the corresponding target density of photons.

For protons, we compute the pγ cooling rate as (Atoyan & Dermer 2003):

$$\begin{array}{cc}& t_{\mathrm{p}\gamma }^{-1}(E_{\mathrm{p}}^{\prime })={\displaystyle {\int }_{E_{\mathrm{th}}/2{\gamma }_{\mathrm{p}}}^{\infty }\mathrm{d}{E}_{\mathrm{ph}}\frac{c{n}_{\mathrm{ph}}(E_{\mathrm{ph}})}{2{\gamma }_{\mathrm{p}}^2{E}_{\mathrm{ph}}^2}}\hfill \\ \hfill & \times {\displaystyle {\int }_{E_{\mathrm{th}}}^{2E_{\mathrm{ph}}{\gamma }_{\mathrm{p}}}\mathrm{d}{E}_{\mathrm{r}}{\sigma }_{\mathrm{p}\gamma }(E_{\mathrm{r}})K_{\mathrm{p}\gamma }(E_{\mathrm{r}})E_{\mathrm{r}},}\hfill \end{array}$$(28)

where ${\gamma }_{\text{p}}=E_{\text{p}}^{\prime }/(m_{\text{p}}{c}^2)$, E_th = 2m_ec² is the threshold energy corresponding to e⁺e-pairs (Bethe-Heitler process), or E_th ≃ 150 MeV in the case of pion production. The corresponding cross sections σ_pγ and inelasticity coefficients K_pγ are taken as in Begelman et al. (1990).

For completeness, we also consider the cooling due to pp interactions with a rate

$$\begin{array}{c}\hfill t_{\mathrm{pp}}^{-1}(E_{\mathrm{p}}^{\prime })\simeq \frac{1}{2}{n}_{\mathrm{j}}^{\prime }c{\sigma }_{\mathrm{pp}}(E_{\mathrm{p}}^{\prime }),\end{array}$$(29)

where the cross section σ_pp is taken from Kelner et al. (2006). In turn, the adiabatic cooling rate, which can be estimated as

$$\begin{array}{c}\hfill t_{\mathrm{ad}}^{-1}\approx \frac{c}{r_{\mathrm{is}}},\end{array}$$

gives a generally negligible effect for high energy protons in the present context.

In Fig. 4, we show the obtained acceleration and cooling rates for protons in the cases of L₀  = 10⁵⁰erg s⁻¹ and L₀  =  10⁵¹erg s⁻¹ in the left and right panels, respectively. Again, the top panels correspond to Γ = 100 and the bottom ones to Γ = 300. As can be seen in the figure, the pγ interactions are always the dominant process. In the plots, we also include the contributions to the cooling rate due to the external photons produced in the jet head and the contribution due to the synchrotron emission of electrons in the internal shock region for ϵ_B = 0.1. It can also be seen that the pγ interactions with external photons are dominant for Γ = 300, while for Γ = 100, the interactions with electron synchrotron photons can become relevant only for very high energy protons if a_ep = 1.

Fig. 4. Cooling and acceleration rates for protons at the internal shocks of CGRBs for L0 = 1050 erg s−1 and L0 = 1051 erg s−1 in the left and right panels, respectively. Top panels correspond to Γ = 100 and the bottom ones to Γ = 300. The different processes are indicated with the following curves: dashed lines are for ϵB = 1, solid lines for ϵB = 0.1, and short-dashed lines for ϵB = 0.01. The blue curves refer to acceleration, orange ones to synchrotron, magenta ones to pp interactions, and black ones to adiabatic cooling. With green, we mark the pγ cooling rates corresponding to a target of external photons (solid lines) and to internal pγ interactions with electron synchrotron photons as targets for aep = 1 (long-short dashed lines) and aep = 100 (long-short-short dashed lines).

We show in Figs. 5 and 6 the obtained electron and proton distributions for L₀ = 10⁵⁰erg s⁻¹ and L₀ = 10⁵¹erg s⁻¹. The top and bottom panels correspond to Γ = 100 and Γ = 300, respectively. In the figures, it can be seen that the electron distributions decrease with ϵ_B. This is because the synchrotron cooling becomes more relevant as the magnetic field increases. The maximum electron energy also decreases with the magnetic field, since the acceleration rate grows with B, but the synchrotron cooling rate grows with B². For Γ = 300 and L₀ = 10⁵⁰ erg s⁻¹, the cases with different a_ep only differ by a constant factor, reflecting only the fact that a different power is injected in the electrons since the SSC process does not play a significant role. In the other cases, SSC can be relevant for a_ep = 1 and high electron energies.

Fig. 5. Electron distributions for L0 = 1050 erg s−1 (left panel) and L0 = 1051 erg s−1 (right panel). The top panels correspond to Γ = 100 and the bottom ones to Γ = 300. The green, orange, and blue curves refer to the cases of ϵB = 0.01, ϵB = 0.1, and ϵB = 1, respectively. The dashed curves correspond to aep = 1 and the solid ones to aep = 100.

Fig. 6. Proton distributions for L0 = 1050 erg s−1 (left panel) and L0 = 1051 erg s−1 (right panel). The top panels correspond to Γ = 100 and the bottom ones to Γ = 300. The green, orange, and blue curves refer to the cases of ϵB = 0.01, ϵB = 0.1, and ϵB = 1, respectively. The dashed curves correspond to aep = 1 and the solid ones to aep = 100.

For protons, the synchrotron cooling is not so relevant, and the effect of increasing the magnetic field leads to a faster acceleration rate. Only in the cases with L₀ = 10⁵¹ erg s⁻¹ and Γ = 100 can this be compensated by a faster cooling via pγ interactions with an increased target of photons due significant electron synchrotron emission, particularly for a_ep = 1. In the rest of the cases, such additional interactions do not play a relevant role.

3.2. Pions, muons, and kaons

Pions are produced mainly through pγ interactions, and we compute the corresponding injection following Hümmer et al. (2010):

$$\begin{array}{cc}& Q_{\mathrm{p}\gamma \to {\pi }^{\pm }}^{\prime }(E_{\pi }^{\prime })={\displaystyle \sum _{\mathrm{it}=\{\mathrm{LR},\mathrm{NR}\}}{N}_{\mathrm{p}}^{\prime }\left(\frac{E_{\pi }^{\prime }}{{\chi }_{\mathrm{it}}}\right){\int }_{E_{\mathrm{th},\pi }/2}^{\infty }{n}_{\mathrm{ph}}\left(\frac{y{\chi }_{\mathrm{it}}{m}_{\mathrm{p}}{c}^2}{E_{\pi }^{\prime }}\right)\frac{M_{\pi }^{\mathrm{it}}}{2y^2}\mathrm{d}y}\hfill \\ \hfill & \times {\displaystyle {\int }_{E_{\mathrm{th},\pi }}^{2y}\mathrm{d}{E}_{\mathrm{r}}{E}_{\mathrm{r}}{\sigma }_{\mathrm{p}\gamma }(E_r),}\hfill \end{array}$$(30)

where the interaction channels “it” refers to are the low and high energy resonances (LR and HR). The corresponding parameters χ_LR, HR and $M_{\pi }^{\mathrm{LR},\mathrm{HR}}$ are listed in Hümmer et al. (2010). Additionally, charged kaons are produced in pγ collisions, although with a subdominant contribution at high energies. The injection of K⁺ can be computed as

$$\begin{array}{cc}& Q_{K^{+}}^{\prime }(E_K^{\prime })=N_{\mathrm{p}}^{\prime }\left(\frac{E_K^{\prime }}{{\chi }_K}\right){\displaystyle {\int }_{E_{\mathrm{th},K}/2}^{\infty }{n}_{\mathrm{ph}}\left(\frac{y{\chi }_K{m}_{\mathrm{p}}{c}^2}{E_K^{\prime }}\right)\frac{1}{2y^2}\mathrm{d}y}\hfill \\ \hfill & \times {\displaystyle {\int }_{E_{\mathrm{th},K}}^{2y}\mathrm{d}{E}_{\mathrm{r}}{E}_{\mathrm{r}}{\sigma }_{\mathrm{p}\gamma }(E_{\mathrm{r}}),}\hfill \end{array}$$(31)

where χ_K = 0.35 and E_th, K = 1 GeV. For the negative pions K⁻, we make the approximation Q_K⁻ ≈ 0.5Q_K⁺, which yields results consistent with the ones obtained by Lipari et al. (2007). The neutral kaons K_L are also produced at an intermediate level between the K⁻ and the K⁺ injections. However, we neglected the K_L contribution, since they only decay to neutrinos via three body channels, which implies that the neutrino energy is low and therefore the corresponding emissivity is expected to be negligible compared to those of the rest of the channels (see the discussion in Hümmer et al. 2010; Petropoulou et al. 2014).

Charged pions can also interact with matter and photons. In the former case, we obtain $t_{\pi \mathrm{p}}^{-1}$ making the simple approximation that ${\sigma }_{\pi \mathrm{p}}\approx \frac{2}{3}{\sigma }_{\mathrm{pp}}$, which is inspired by the fact that while protons are made up of three valence quarks, pions are formed by only two of them (Gaisser 1990). As for pion-photon interactions, we adopt an approach similar to Lipari et al. (2007) and consider the production of the resonances ρ[770], a₁[1260], a₂[1320], and b₁[1235], since they all have a possible decay channel to πγ according to the Particle Data Group (PDG) review (Particle Data Group 2022). Therefore, for low values of center-of-mass energies, $\sqrt{s}<2\mathrm{GeV}$, we describe the cross section as the sum of the corresponding contributions using the Breit-Wigner expression for the mentioned resonances (Wigner 1946). For higher values of s, we adopted the general parametrization derived for hadronic interactions (Block & Halzen 2004) along with the fitting constants derived in the PDG review (Particle Data Group 2016) to account for the total cross sections for the processes pp, pγ, and πp. These can be related to the πγ cross section if we apply the factorization relation that follows from the Regge theory (e.g., Gribov 2023)

$$\begin{array}{c}\hfill {\sigma }_{\pi \gamma }=\sqrt{{\sigma }_{\pi \pi }{\sigma }_{\gamma \gamma }}\simeq {\sigma }_{\pi \mathrm{p}}\sqrt{\frac{{\sigma }_{\gamma \gamma }}{{\sigma }_{\mathrm{pp}}}},\end{array}$$(32)

where we also use ${\sigma }_{\pi \pi }={\sigma }_{\pi \text{p}}^2/{\sigma }_{\text{pp}}$. The resulting cross section is shown in Fig. 7, and it was adopted to estimate the πγ cooling rate with an expression analogous to that of Eq. (28) for pγ interactions, taking K_πγ ≈ 0.5 for the inelasticity.

Fig. 7. Pion-photon cross section adopted to obtain $t_{\pi \gamma }^{-1}$. The bottom x-axis refers to the square root of the center-of-mass energy s1/2, and the top one marks the corresponding photon energy in the pion rest frame, Er.

We show all the obtained cooling rates for pions in Fig. 8, with L₀ = 10⁵⁰ erg s⁻¹ and L₀ = 10⁵¹ erg s⁻¹ in the left and right panels, respectively, and Γ = 100 in the top panels and Γ = 300 in the bottom ones. In these plots, the decay is shown to be faster than the cooling for energies $E_{\pi }^{\prime }\lesssim {10}^5\text{GeV}$ for Γ = 100, and for $E_{\pi }^{\prime }\lesssim 2\times {10}^4\text{GeV}$ for Γ = 300. As for the pion cooling, we unified the πγ cooling rate in the figure and present the sum of the contribution due to collisions with internal photons from the synchrotron emission of electrons in the internal shock zone as well as the one corresponding to external photons produced in the jet head. These πγ interactions provide the dominant cooling mechanism for high energy pions regardless of the value of a_ep, except for the cases with Γ = 100 and L₀ = 10⁵¹ erg s⁻¹ at energies greater than 10⁵ GeV, which is where pion synchrotron can dominate. We note that πγ interactions actually refer to inelastic collisions characterized by the cross section obtained in Eq. (32). For high energies, where the pion decay rate becomes subdominant, such pγ interactions give the dominant cooling rate as compared to the corresponding pion-IC interactions. The latter take place at a rate similar to the muon IC cooling rate, and refer to the elastic case where there is only a pion and a photon in the final state.

Fig. 8. Pion decay and cooling rates for L0 = 1050 erg s−1 (left panel) and L0 = 1051 erg s−1 (right panel). The top panels correspond to Γ = 100 and the bottom ones to Γ = 300. The blue curves correspond to the pion decay rate, and the green curves represent the πγ interaction with external and internal photons with aep = 1 (dashed lines) and with aep = 100 (solid lines). The pion synchrotron cooling is marked with orange curves for ϵB = 0.01 (short-dashed lines), ϵB = 0.1 (solid lines), and ϵB = 1 (long-dashed lines). The solid black curves mark the adiabatic cooling rate.

The case of the charged kaons is quite similar to that of the pions. They decay about twice as fast as pions, and the K^±γ interactions are assumed to occur at the same rate as for pions. However, the synchrotron cooling rate for kaons is a factor (m_π/m_K)³ ≃ 2.3 × 10⁻² lower than the one for pions.

Muons are produced by the decay of pions, and we compute their injection $Q_{\mu }^{\prime }(E_{\mu }^{\prime })$ as in Reynoso & Romero (2009), making use of the formulas of Lipari et al. (2007), which depend on the distributions of pions $N_{\pi }^{\prime }$ and account for the kinematics of the decay. The relevant cooling rates for muons are shown in Fig. 9 for L₀ = 10⁵⁰ erg s⁻¹ and L₀ = 10⁵¹ erg s⁻¹, with Γ = 100 in the top panels and Γ = 300 in the bottom ones. In the figure, it can be seen that decay takes place faster than the cooling for $E_{\mu }^{\prime }\lesssim {10}^4\text{GeV}$, while for higher energies, the synchrotron and IC become significant. It can also be seen that the IC cooling of muons at high energies is sensitive to the value of a_ep, except for the cases with Γ = 300 and L₀ = 10⁵⁰ erg s⁻¹.

Fig. 9. Muon decay and cooling rates for L0 = 1050 erg s−1 (left panel) and L0 = 1051 erg s−1 (right panel). The top panels correspond to Γ = 100 and the bottom ones to Γ = 300. The blue curves correspond to the muon decay rate, and the green curves represent the IC interaction with external and internal photons with aep = 1 (dashed lines) and with aep = 100 (solid lines). The muon synchrotron cooling is marked with orange curves for ϵB = 0.01 (short-dashed lines), ϵB = 0.1 (solid lines), and ϵB = 1 (long-dashed lines). The solid black curves mark the adiabatic cooling rate.

The pion, muon, and kaon distributions are obtained using Eq. (22), making the replacement $T_{\mathrm{esc}}^{-1}\to T_{\mathrm{esc}}^{-1}+T_{i,\mathrm{d}}^{-1}$ where the decay timescales for each of these particle types are

$$\begin{array}{cc}& T_{\pi ,\mathrm{d}}(E_{\pi }^{\prime })\simeq 2.6\times {10}^{-8}\mathrm{s}\left(\frac{E_{\pi }^{\prime }}{m_{\pi }{c}^2}\right)\hfill \end{array}$$(33)

$$\begin{array}{cc}& T_{\mu ,\mathrm{d}}(E_{\mu }^{\prime })\simeq 2.2\times {10}^{-6}\mathrm{s}\left(\frac{E_{\mu }^{\prime }}{m_{\mu }{c}^2}\right)\hfill \end{array}$$(34)

$$\begin{array}{cc}& T_{K,\mathrm{d}}(E_K^{\prime })\simeq 1.2\times {10}^{-8}\mathrm{s}\left(\frac{E_K^{\prime }}{m_K{c}^2}\right).\hfill \end{array}$$(35)

We show the obtained $N_{\pi }^{\prime }(E_{\pi }^{\prime })$, $N_{\mu }^{\prime }(E_{\mu }^{\prime })$, and $N_K^{\prime }(E_K^{\prime })$ in the top, middle, and bottom panels of Fig. 10 for Γ = 100 and of Fig. 11 for Γ = 300. When comparing the resulting pion distributions, we observe that the different outcomes are sensitive to the proton distribution corresponding to each case as well as to the relevance of both πγ and synchrotron emission at high energies. The same goes for the muon distributions, which directly depend on the corresponding pion distribution, and since they have a slower decay rate, the IC and synchrotron losses affect their distribution at lower energies than for pions.

Fig. 10. Distributions of pions (top), muons (middle), and charged kaons (bottom) for Γ = 100. The left plots correspond to L0 = 1050 erg s−1 and the right plots to L0 = 1051 erg s−1. The green, orange, and blue curves refer to the cases of ϵB = 0.01, ϵB = 0.1, and ϵB = 1, respectively. The dashed curves correspond to aep = 1 and the solid curves to aep = 100.

Fig. 11. Distributions of pions (top), muons (middle), and charged kaons (bottom) for Γ = 300. The left plots correspond to L0 = 1050 erg s−1 and the right plots to L0 = 1051 erg s−1. The green, orange, and blue curves refer to the cases of ϵB = 0.01, ϵB = 0.1, and ϵB = 1, respectively. The dashed curves correspond to aep = 1 and the solid ones to aep = 100.

We note that for Γ = 100 and L₀ = 10⁵¹ erg s⁻¹, the synchrotron cooling can be important at high energies for muons and even pions. Specifically, in the top-right panel of Fig. 10, the obtained pion distributions for ϵ_B = 1 and ϵ_B = 0.1 extend up to similar energies despite the fact that in the former case, protons are accelerated to higher energies, but this is compensated by more efficient pion synchrotron losses than in the latter case. Similarly, in the middle-right panel of the same figure, the muon distributions extend to lower energies as ϵ_B increases, as a consequence of more severe synchrotron losses. These losses are also significant for Γ = 100 and ϵ_B = 1, as can be seen in the middle-left panel of Fig. 10.

In the other cases, since synchrotron cooling of pions or muons does not dominate, increasing ϵ_B just leads to a distribution that extends to higher energies due to a more efficient acceleration of parent protons. As for kaons, their distributions peak at higher energies than for pions, in agreement with Hümmer et al. (2010), and K^±γ interactions are the most relevant cooling mechanism. Therefore, even for Γ = 100, L₀ = 10⁵¹ erg s⁻¹, and ϵ_B = 1, synchrotron cooling does not significantly affect their distribution, as can be seen for a_ep = 100 in the bottom-right panel of Fig. 10. If a_ep = 1, instead, the cases with ϵ_B = 0.1 and ϵ_B = 1 become similar because the increase of the magnetic field is compensated by faster Kγ interactions with photons of the electron synchrotron emission. We also note that different values of the proton-to-electron ratio affect the muon distributions more significantly if Γ = 100, but only slightly for Γ = 300 and L₀ = 10⁵¹ erg s⁻¹.

3.3. Neutrinos

Once the pion, muon, and kaon distributions are obtained, we can compute the injection or emissivity of neutrinos plus antineutrinos, as we are interested in studying the overall flavor composition. The contribution to the emissivity of muon neutrinos and antineutrinos resulting from the pion decays π⁺ → μ⁺ν_μ and ${\pi }^{-}\to {\mu }^{-}{\overline{\nu }}_{\mu }$ can be computed as (e.g., Lipari et al. 2007):

$$\begin{array}{c}\hfill Q_{\pi \to {\nu }_{\mu }}^{\prime }(E_{\nu }^{\prime })={\displaystyle {\int }_{E_{\pi ,\min }}^{\infty }\mathrm{d}{E}_{\pi }\frac{N_{\pi }^{\prime }(E_{\pi })}{E_{\pi }(1-r_{\pi })T_{\pi ,\mathrm{d}}(E_{\pi })},}\end{array}$$(36)

where $E_{\pi ,\min }=\frac{E_{\nu }^{\prime }}{1-r_{\pi }}$, with r_π = (m_μ/m_π)². A similar expression corresponds to the kaon decays K⁺ → μ⁺ν_μ and $K^{-}\to {\mu }^{-}{\overline{\nu }}_{\mu }$ weighted by the branching ratio Br_{K → μνμ} ≃ 0.63, making the replacements $N_{\pi }^{\prime }\to N_K^{\prime }$, T_π,d → T_K,d, and r_π → r_K.

The contribution to the muon flavor at the production stage generated by the decays of muons ${\mu }^{+}\to e^{+}{\nu }_{\mathrm{e}}{\overline{\nu }}_{\mu }$ and ${\mu }^{-}\to e^{-}{\overline{\nu }}_{\mathrm{e}}{\nu }_{\mu }$ is accounted by the following emissivity (Lipari et al. 2007):

$$\begin{array}{cc}& Q_{\mu \to {\nu }_{\mu }}^{\prime }(E_{\nu }^{\prime })={\displaystyle \sum _{i=1}^4{\int }_{E_{\nu }^{\prime }}^{\infty }\frac{\mathrm{d}{E}_{\mu }}{E_{\mu }}\frac{N_{{\mu }_i}^{\prime }(E_{\mu })}{T_{\mu ,\mathrm{d}}(E_{\mu })}}\hfill \\ \hfill & \times [\frac{5}{3}-3x^2+\frac{4}{3}{x}^3+(3x^2-\frac{1}{3}-\frac{8x^3}{3})h_i],\hfill \end{array}$$(37)

where $x=E_{\nu }^{\prime }/E_{\mu }$, and the i index refers to the different muons depending on their charge and helecity (i.e., ${\mu }_{1,2}={\mu }_L^{-,+}$, ${\mu }_{3,4}={\mu }_R^{-,+}$), with h_1, 2 = −h_3, 4 = −1. Similarly, the injection of electron neutrinos and antineutrinos can be computed as:

$$\begin{array}{cc}& Q_{\mu \to {\nu }_{\mathrm{e}}}^{\prime }(E_{\nu }^{\prime })={\displaystyle \sum _{i=1}^4{\int }_{E_{\nu }^{\prime }}^{\infty }\frac{\mathrm{d}{E}_{\mu }}{E_{\mu }}\frac{N_{{\mu }_i}^{\prime }(E_{\mu })}{T_{\mu ,\mathrm{d}}(E_{\mu })}}\hfill \\ \hfill & \times [2-6x^2+4x^3+(2-12x+18x^2-8x^3)h_i].\hfill \end{array}$$(38)

In Figs. 12 and 13, we show the different contributions to neutrino emissivities for L₀ = 10⁵⁰ erg s⁻¹ and L₀ = 10⁵¹ erg s⁻¹, respectively, both for a_ep = 1. In these figures, the left panels correspond to Γ = 100 and the ones on the right to Γ = 300. The cases with ϵ_B = {0.01, 0.1, 1} are shown in the top, middle, and bottom panels, respectively, and the dashed lines refer to the outcomes that would be obtained if losses of pions, muons, and kaons were neglected. It can be seen that even in the cases of low magnetic field (i.e., Fig. 12 for L₀ = 10⁵⁰ erg s⁻¹, ϵ_B = 0.01, and Γ = 300), the neutrino emission is reduced at high energies due to interactions of the parent particles with the soft photon background, as mentioned above. In the cases of higher magnetic fields, synchrotron losses generate further depletion of the neutrino emissivity at high energies. In particular, the muon contribution is the most affected. Kaons are only moderately affected, and the neutrinos resulting from their decays give a small correction at the highest energies, similarly to what has been obtained in other previously studied cases (e.g., Hümmer et al. 2010, 2012).

Fig. 12. Emissivities of muon neutrinos and antineutrinos from the decay of pions (blue curves), muons (orange curves), and kaons (green curves) for L0 = 1050 erg s−1. The dashed curves represent the obtained output if all losses of the parent particles are neglected, while the solid curves mark the results if all the losses discussed are considered. The left plots correspond to aep = 1 and the right plots to aep = 100.

Fig. 13. Emissivities of muon neutrinos and antineutrinos from the decay of pions (blue curves), muons (orange curves), and kaons (green curves) for L0 = 1051 erg s−1. The dashed curves represent the obtained output if all losses of the parent particles are neglected, while the solid curves mark the results if all the losses discussed are considered. The left plots correspond to aep = 1 and the right plots to aep = 100.

4. Diffuse neutrino flux and flavor ratios

After considering in detail the cooling effect of the parent particles that generate the neutrinos, we are next interested in obtaining the diffuse neutrino flux produced by the cumulative CGRBs taking place at different redshifts in the universe. This would allow us to assess the contribution of these sources to the total neutrino flux that reaches the Earth, as well as their flavor composition. We assume that the generation rate of CGRBs R_CGRB(z) is proportional to the star formation rate ρ_⋆(z) (e.g., He et al. 2018; Fasano et al. 2021):

$$\begin{array}{c}\hfill R_{\mathrm{CGRB}}(z)=A_{\mathrm{CGRB}}{\rho }_{\star }(z),\end{array}$$(39)

where A_CGRB is a normalization constant. Following Madau & Dickinson (2014), we adopt

$$\begin{array}{c}\hfill {\rho }_{\star }(z)=0.015M_{\odot }{\mathrm{yr}}^{-1}{\mathrm{Mpc}}^{-3}\frac{{(1+z)}^{2.7}}{1+{[(1+z)/2.9]}^{5.6}}.\end{array}$$(40)

The source evolution function R_CGRB(z) is used to weight the neutrino spectrum produced by a single typical CGBRB,

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{\nu ,z}(E_{\nu ,z})}{\mathrm{d}{E}_{\nu ,z}}=\mathrm{\Delta }\mathrm{\Omega }{Q}_{\nu ,z}(E_{\nu ,z})\mathrm{\Delta }V{T}_{\mathrm{CGRB}}\left(\frac{\mathrm{\Delta }\mathrm{\Omega }}{4\pi }\right),\end{array}$$(41)

where E_ν, z = E_ν(1 + z) is the local neutrino energy at redshift z corresponding to an energy E_ν for an observer on Earth. We introduce the ratio ΔΩ/(4π) to account for the fact that these sources are beamed, and we use an isotropic equivalent luminosity to obtain the neutrino emissivity $Q_{\nu }^{\prime }(E_{\nu }^{\prime })$ at the comoving frame of source. The latter is then transformed to the local frame according to

$$\begin{array}{cc}& Q_{\nu ,z}(E_{\nu ,z})\simeq D{Q}_{\nu }^{\prime }(E_{\nu }^{\prime })\hfill \\ \hfill & \simeq 2\mathrm{\Gamma }{Q}_{\nu }^{\prime }(E_{\nu ,z}/(2\mathrm{\Gamma })).\hfill \end{array}$$(42)

Given that the comoving neutrino density can be expressed as (e.g., Murase 2007)

$$\begin{array}{c}\hfill \mathrm{d}{n}_{\nu }(E_{\nu })=R_{\mathrm{CGRB}}(z){(1+z)}^3\frac{\mathrm{d}t}{\mathrm{d}z}\mathrm{d}z\frac{\mathrm{d}{N}_{\nu }^{\prime }(E_{\nu }^{\prime })}{\mathrm{d}{E}_{\nu }^{\prime }}\mathrm{d}{E}_{\nu }^{\prime }{(1+z)}^{-3},\end{array}$$(43)

and considering that the differential neutrino flux at z = 0 is $\phi (E_{\nu })=\frac{c}{4\pi }\frac{\mathrm{d}{n}_{\nu }}{\mathrm{d}{E}_{\nu }}$, it follows that

$$\begin{array}{cc}& {\phi }_{{\nu }_{\alpha }}(E_{\nu })=\frac{c}{4\pi H_0}{\displaystyle {\int }_0^{z_{\max }}\frac{\mathrm{d}z{R}_{\mathrm{CGRB}}(z)}{\sqrt{{\mathrm{\Omega }}_{\mathrm{\Lambda }}+{\mathrm{\Omega }}_{\mathrm{m}}{(1+z)}^3}}}\hfill \\ \hfill & [\frac{\mathrm{d}{N}_{{\nu }_{\mathrm{e}},z}(E_{\nu }(1+z)}{\mathrm{d}{E}_{\nu ,z}}{P}_{\mathrm{e}\alpha }+\frac{\mathrm{d}{N}_{{\nu }_{\mu },z}(E_{\nu }(1+z)}{\mathrm{d}{E}_{\nu ,z}}{P}_{\mu \alpha }].\hfill \end{array}$$(44)

Here, we adopt H₀ = 70 km s⁻¹ Mpc⁻¹, and we also include the effect of neutrino oscillation during their propagation by multiplying the emitted neutrino spectrum of the flavors e and μ by P_eα and P_μα, which represent the probabilities of oscillation to the flavor of interest to be observed on Earth, α = {e,μ,τ}. These probabilities depend on three mixing angles (θ₁₂, θ₂₃, θ₁₃) and a phase (δ_cp) that characterize that the mixing matrix U_αβ. We adopt the following values obtained by Esteban et al. (2020) through global fits to data from neutrino baseline experiments:

$$\begin{array}{cc}& {\theta }_{12}=32.4^{°}\hfill \end{array}$$(45)

$$\begin{array}{cc}& {\theta }_{23}={49}^{°}\hfill \end{array}$$(46)

$$\begin{array}{cc}& {\theta }_{13}=8.6^{°}\hfill \end{array}$$(47)

$$\begin{array}{cc}& {\delta }_{\mathrm{cp}}={197}^{°},\hfill \end{array}$$(48)

and they yield the following values for the oscillation probabilities: P_ee ≃ 0.552, P_eμ ≃ 0.171, P_eτ ≃ 0.276, and P_μτ ≃ 0.375.

In Fig. 14, we show the resulting muon neutrino and antineutrino fluxes for L₀ = 10⁵⁰ erg s⁻¹ and L₀ = 10⁵¹ erg s⁻¹, with Γ = 100 in the left panels and Γ = 300 in the right ones. We include the cases of different values of the magnetic fields and of the proton-to-electron ratio a_ep in the figure. The shape and relative importance of the resulting fluxes for the different cases studied can be understood by looking at the distributions of the pions, muons and kaons corresponding to each case.

Fig. 14. Diffuse flux of muon neutrinos and antimuons for L0 = 1050 erg s−1 (top panels) and L0 = 1051 erg s−1 (bottom panels). The results obtained with Γ = 100 and Γ = 300 are shown in the left and right panels, respectively. The green, orange, and blue curves refer to the cases of ϵB = 0.01, ϵB = 0.1, and ϵB = 1, respectively. The dashed curves correspond to aep = 1 and the solid ones to aep = 100.

The normalizing constant mentioned above is fixed at $A_{\text{CGRB}}=8\times {10}^{-8}{M}_{\odot }^{-1}$ for L₀ = 10⁵¹ erg s⁻¹, $A_{\text{CGRB}}=8\times {10}^{-7}{M}_{\odot }^{-1}$ for L₀ = 10⁵⁰ erg s⁻¹ and Γ = 100, and $A_{\text{CGRB}}=4\times {10}^{-7}{M}_{\odot }^{-1}$ for L₀ = 10⁵⁰ erg s⁻¹ and Γ = 300. These values were obtained in order to have our resulting illustrative fluxes for a_ep = 100 at the level of the fit obtained with IceCube 10-year data for muon tracks (Stettner 2019).

We note that even for $A_{\text{CGRB}}=8\times {10}^{-7}{M}_{\odot }^{-1}$, the corresponding local rate of choked jet events is R_CGRB(z = 0)≈12 yr⁻¹Gpc⁻³, which is within the range of plausible values when taking into account that a certain fraction f_jet of the total core-collapse SNe harbor jets (see, e.g., Senno et al. 2016; Denton & Tamborra 2018; He et al. 2018; Fasano et al. 2021). The precise fraction is actually unknown, but it plays a crucial role in establishing a link between CGRBs and such SNe and therefore in associating the location of neutrino events to a possible location of an SN. Recent studies in this area (e.g., Guetta et al. 2020; Chang et al. 2022) suggest that the mentioned fraction should be low, f_jet ≲ 0.1, and although no correlation was confirmed for SNe and neutrino events, CGRBs can still be the dominant sources of the detected diffuse neutrino flux.

To make a comparison of the neutrino fluxes of the different flavors, it is useful to obtain the energy dependent flavor ratios,

$$\begin{array}{c}\hfill f_{\alpha }(E_{\nu })=\frac{{\phi }_{{\nu }_{\alpha }}(E_{\nu })}{{\phi }_{{\nu }_{\mathrm{e}}}(E_{\nu })+{\phi }_{{\nu }_{\mu }}(E_{\nu })+{\phi }_{{\nu }_{\tau }}(E_{\nu })}.\end{array}$$

This includes all the effects of energy losses suffered by the parent particles. We plot the corresponding ratios in Figs. 15 and 16 for Γ = 100 and Γ = 300, respectively. We also consider, for reference, the standard estimate derived from the fact that for pion decays, a ratio of 2 : 1 of muon neutrinos to electron neutrinos is produced at the sources, and neglecting all losses and without considering the decay spectra of the different particles, it follows that

$$\begin{array}{cc}& f_{\mathrm{e},0}=2P_{\mathrm{e}\mu }+P_{\mathrm{ee}}\simeq 0.30\hfill \end{array}$$(49)

$$\begin{array}{cc}& f_{\mu ,0}=2P_{\mu \mu }+P_{\mathrm{e}\mu }\simeq 0.36\hfill \end{array}$$(50)

$$\begin{array}{cc}& f_{\tau ,0}=2P_{\mu \tau }+P_{\mathrm{e}\tau }\simeq 0.34.\hfill \end{array}$$(51)

Fig. 15. Energy dependent neutrino flavor ratios for L0 = 1050 erg s−1 (left panels) and L0 = 1051 erg s−1 (right panels) for Γ = 100. The top plots correspond to aep = 1 and the bottom ones to aep = 100. The green, orange, and blue curves correspond to the electron, muon, and tau flavor ratios, respectively. The long-dashed, solid, and short-dashed lines refer to the results for ϵB = 0.01, ϵB = 0.1, and ϵB = 1, respectively. The flavor ratios expected without considering losses are indicated with gray dotted lines for the electron flavor, gray long-dashed lines for the muon flavor, and gray long-short dashed lines for the tau flavor.

Fig. 16. Energy dependent neutrino flavor ratios for L0 = 1050 erg s−1 (left panels) and L0 = 1051 erg s−1 (right panels). Both plots are for Γ = 300. The green, orange, and blue curves correspond to the electron, muon, and tau flavor ratios, respectively. The long-dashed, solid, and short-dashed lines refer to the results for ϵB = 0.01, ϵB = 0.1, and ϵB = 1, respectively. The flavor ratios expected without considering losses are indicated with gray dotted lines for the electron flavor, gray long-dashed lines for the muon flavor, and gray long-short dashed lines for the tau flavor.

This estimate gives roughly (f_e, 0 : f_μ, 0 : f_τ, 0)∝(1 : 1 : 1) (e.g., Bustamante & Ahlers 2019; Bustamante & Tamborra 2020). Therefore, it can be seen in the mentioned figures that as the neutrino energy increases, the flavor ratios obtained as a result of our full calculation differ more from the expected outcome if losses are neglected. Specifically, the losses discussed (interactions with soft photons, and synchrotron emission) cause the electron flavor ratio to decrease to ≈0.17, the muon ratio to increase up to ≈0.45, and the tau flavor to increase slightly to 0.36. The smooth transition to these values takes place at different energies for the different cases considered.

For L₀ = 10⁵¹ erg s⁻¹ in particular, the flavor ratios are modified at energies lower than for L₀ = 10⁵⁰ erg s⁻¹, given that the interaction rates are higher in the former cases than in the latter ones. When comparing the top and lower panels of Fig. 15 for Γ = 100, the flavor ratios are shown to be affected at slightly lower energies for a_ep = 1 than for a_ep = 100, and it is the interactions with more electron synchrotron radiation that makes the difference.

As discussed above, for Γ = 100 and L₀ = 10⁵¹ erg s⁻¹, synchrotron cooling of pions and muons become more relevant for higher values of ϵ_B, and this affects the flavor ratios at lower energies. However, for Γ = 100 and L₀ = 10⁵⁰ erg s⁻¹, synchrotron cooling of muons is only relevant if ϵ_B = 1, and therefore the cases with lower magnetic fields have similar flavor ratios (long dashed and solid lines in left panels of Fig. 15) and are affected mostly by pion and muon interactions with soft photons.

As for the cases with Γ = 300, the flavor ratios we present in Fig. 16 refer to a_ep = 1. A different situation arises in the cases with L₀ = 10⁵⁰ erg s⁻¹ shown in the left panel as compared those with L₀ = 10⁵¹ erg s⁻¹ shown in the right panel. In the former cases, unlike the rest of the cases studied, the corresponding muon distributions extend only to energies below ∼10⁵ GeV (see Fig. 11, middle-left panel), and cooling is mostly dominated by the muon IC process, even for the case of ϵ_B = 1. In fact, at high energies, the IC cooling rate is in the Klein-Nishina regime (i.e., decreasing with the energy). This implies that increasing ϵ_B only increases the acceleration rate and therefore the maximum energy of the parent protons and of the produced pions and muons. However, since muons cool more slowly at higher energies in these particular cases (for Γ = 300 and L₀ = 10⁵⁰ erg s⁻¹), the resulting curves for the flavor ratios for low values of ϵ_B appear slightly more affected by interactions at high energies as compared those for higher magnetic fields. This situation is different for Γ = 300 and L₀ = 10⁵¹ erg s⁻¹ as wells as in the rest of the cases studied, as the IC cooling of muons does not always dominate if the magnetic field is increased because synchrotron cooling can become important at high energies, and this results in similar flavor ratios to the obtained, for instance, using Γ = 100 and L = 10⁵¹ erg s⁻¹.

Given the relatively low number of detected neutrino events of astrophysical origin, a more useful observable to be extracted from experiments may be the flavor ratios of the integrated fluxes. Depending on the energy range of integration, the effects discussed for the energy dependent ratios may also be appreciated in these flavor ratios. If we integrate the fluxes on the neutrino energy above a minimum energy E_ν,min

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{{\nu }_{\alpha }}(E_{\nu ,\min })={\displaystyle {\int }_{E_{\nu ,\min }}^{\infty }\mathrm{d}{E}_{\nu }{\phi }_{{\nu }_{\alpha }}(E_{\nu }),}\end{array}$$

the corresponding flavor ratios can be obtained as

$$\begin{array}{c}\hfill F_{\alpha }(E_{\nu ,\min })=\frac{{\mathrm{\Phi }}_{{\nu }_{\alpha }}(E_{\nu ,\min })}{{\mathrm{\Phi }}_{{\nu }_{\mathrm{e}}}(E_{\nu ,\min })+{\mathrm{\Phi }}_{{\nu }_{\mu }}(E_{\nu ,\min })+{\mathrm{\Phi }}_{{\nu }_{\tau }}(E_{\nu ,\min })}.\end{array}$$(52)

We plot the results in Figs. 17 and 18. The plots show that the effects of cooling can be better observed by concentrating on the neutrino fluxes of the different flavors at energies higher than ∼5 × 10⁵ GeV. Otherwise, the more abundant events of lower energies tend to mask the departure from the standard values of the flavor ratios corresponding to no cooling of the secondary particles.

Fig. 17. Flavor ratios of the integrated fluxes of neutrinos with energy above Eν,min for L0 = 1050 erg s−1 (left panels) and L0 = 1051 erg s−1 (right panels). All plots are for Γ = 100. The top plots correspond to aep = 1 and the bottom ones to aep = 100. The green, orange, and blue curves correspond to the electron, muon, and tau flavor ratios, respectively. The long-dashed, solid, and short-dashed lines refer to the results for ϵB) = 0.01, ϵB = 0.1, and ϵB = 1, respectively. The flavor ratios expected without considering losses are indicated with gray dotted lines for the electron flavor, gray long-dashed lines for the muon flavor, and gray long-short dashed lines for the tau flavor.

5. Discussion

We have studied in detail the processes leading to the production of high energy neutrinos in CGRBs by making use of steady-state kinetic equations to account for the cooling mechanisms of all the intervening particle populations, namely, electrons, protons, pions, kaons, and muons. Typically, the target photons considered for protons accelerated at the internal shocks are those that escape from the shocked jet head. In the present work, we also considered as targets the synchrotron photons emitted by electrons accelerated at the internal shocks as well as the protons. This implies a rise in the pγ cooling rate only in the cases with a powerful jet L₀ = 10⁵¹ erg s⁻¹ and Γ = 100, and therefore means a decrease in the maximum proton energy.

We have accounted for IC cooling to obtain the particle distributions of electrons as well as of muons. Pions and kaons produced in pγ interactions are also affected by interactions with the mentioned target of soft photons. We found that these interactions can indeed cause important losses of the mentioned particles and therefore a decrease in the neutrino flux at high energies, even in the cases where the magnetic field is relatively low and no significant losses by synchrotron emission are expected. As a consequence of these losses, the flavor ratios become modified: The electron flavor ratio decreases, and the muon flavor ratio increases. This effect gradually manifests in the energy dependent flavor ratios at high energies, as can be seen in Figs. 15–18, where the different curves represent possible signatures of key physical conditions in the jets of CGRBs, such as the magnetic field, proton-to-electron ratio, Lorentz factor, and the jet power. For these parameters, we adopted typically plausible values, according to several works describing the same general context. In particular, by requiring the flow to be optically thin at the internal shock region and that this region be placed inside the collimation shock, we made sure that the sets of parameters we used would imply shocks that are not mediated by radiation, since that allows for efficient particle acceleration.

Fig. 18. Flavor ratios of the integrated fluxes of neutrinos with energy above Eν,min for L0 = 1050erg s−1 (left panels) and L0 = 1051 erg s−1 (right panels). Both plots are for Γ = 100. The green, orange, and blue curves correspond to the electron, muon, and tau flavor ratios, respectively. The long-dashed, solid, and short-dashed lines refer to the results for ϵB) = 0.01, ϵB = 0.1, and ϵB = 1, respectively. The flavor ratios expected without considering losses are indicated with gray dotted lines for the electron flavor, gray long-dashed lines for the muon flavor, and gray long-short dashed lines for the tau flavor.

A treatment like the present one could be applied to other compact sources, such as tidal disruption events, where the magnetic field can be high and therefore significant cooling of pions and/or muons can take place (e.g., Senno et al. 2017). The cores of AGN could also be studied by adapting the present model, such as in the case of NGC 1068. As described by Murase (2022), the possible neutrino emission site could be relatively close to the central BH (e.g., in a corona). Another possibility is that it is located further away where protons would be accelerated through shocks generated by a wind from the accretion disk (Lamastra et al. 2016). In this latter case, a sub-gauss magnetic field would yield the expected no-loss flavor ratios, ≈(1 : 1 : 1), while in the case of a corona, the magnetic field would be expected to be higher (possibly B ≳ 1000 G), leading to flavor ratios that would depart from the no-loss result for high neutrino energies.

Still, in order to be able to measure the energy dependent flavor ratios, as well as the flavor ratios of the integrated fluxes for different minimum energies, more observation time and larger neutrino telescopes are necessary. In general, such resources would allow important information on the physical conditions at the neutrino sources to be obtained. In the specific case of CGRBs, these data would help probe the consequences of the generation of such explosive, gamma-ray hidden events, with the physical conditions normally expected for them, providing new clues as to whether these types of sources can indeed contribute significantly to the diffuse neutrino flux.

Acknowledgments

We thank ANCyT and Universidad Nacional de Mar del Plata for their financial support through grants PICT 2021-GRF-T1-00725 and EXA1102/22, respectively.

References

All Tables

All Figures

	Fig. 1. Schematic view (not to scale) of the jet propagating inside a star.
In the text

Fig. 2. Left panel: optical depth at the internal shock region as a function of the Lorentz factor of the jet for (L0, δt)=(1050 erg s−1, 0.001 s) and (L0, δt)=(1051 erg s−1, 0.01 s), both marked by the black curve, while the blue line corresponds to (L0, δt)=(1050 erg s−1, 0.01 s) and the orange one to (L0, δt)=(1051 erg s−1, 0.001 s). Right panel: positions of the jet head (rh) as a function of the Lorentz factor of the jet for L0 = 1050 erg s−1 and L0 = 1051 erg s−1 marked by long-dashed blue and orange lines, respectively. Solid blue and orange lines correspond to the positions of collimation shock (rcs) for L0 = 1050 erg s−1 and L0 = 1051 erg s−1, respectively. The positions of the internal shocks (ris) for δt = 0.001 s and δt = 0.01 s are indicated with black short-dashed and solid lines, respectively.

In the text

Fig. 3. Cooling and acceleration rates for electrons at the internal shocks of CGRBs for L0 = 1050 erg s−1 and L0 = 1051 erg s−1 in the left and right panels, respectively. Top panels correspond to Γ = 100 and the bottom ones to Γ = 300. The different processes are indicated with the following curves: Dashed lines are for ϵB = 1, solid lines for ϵB = 0.1, and short-dashed lines for ϵB = 0.01. The blue curves refer to acceleration, orange ones to synchrotron, and green ones to external IC. The circles refer to SSC with (aep = 1, ϵB = 0.1) and squares to SSC with (aep = 100, ϵB = 0.1). With dashed gray lines, we also show the auxiliary results used to obtain the SSC rate with the iterative method mentioned in the text.

In the text

Fig. 4. Cooling and acceleration rates for protons at the internal shocks of CGRBs for L0 = 1050 erg s−1 and L0 = 1051 erg s−1 in the left and right panels, respectively. Top panels correspond to Γ = 100 and the bottom ones to Γ = 300. The different processes are indicated with the following curves: dashed lines are for ϵB = 1, solid lines for ϵB = 0.1, and short-dashed lines for ϵB = 0.01. The blue curves refer to acceleration, orange ones to synchrotron, magenta ones to pp interactions, and black ones to adiabatic cooling. With green, we mark the pγ cooling rates corresponding to a target of external photons (solid lines) and to internal pγ interactions with electron synchrotron photons as targets for aep = 1 (long-short dashed lines) and aep = 100 (long-short-short dashed lines).

In the text

	Fig. 5. Electron distributions for L0 = 1050 erg s−1 (left panel) and L0 = 1051 erg s−1 (right panel). The top panels correspond to Γ = 100 and the bottom ones to Γ = 300. The green, orange, and blue curves refer to the cases of ϵB = 0.01, ϵB = 0.1, and ϵB = 1, respectively. The dashed curves correspond to aep = 1 and the solid ones to aep = 100.
In the text

	Fig. 6. Proton distributions for L0 = 1050 erg s−1 (left panel) and L0 = 1051 erg s−1 (right panel). The top panels correspond to Γ = 100 and the bottom ones to Γ = 300. The green, orange, and blue curves refer to the cases of ϵB = 0.01, ϵB = 0.1, and ϵB = 1, respectively. The dashed curves correspond to aep = 1 and the solid ones to aep = 100.
In the text

	Fig. 7. Pion-photon cross section adopted to obtain $t_{\pi \gamma }^{-1}$. The bottom x-axis refers to the square root of the center-of-mass energy s1/2, and the top one marks the corresponding photon energy in the pion rest frame, Er.
In the text

Fig. 8. Pion decay and cooling rates for L0 = 1050 erg s−1 (left panel) and L0 = 1051 erg s−1 (right panel). The top panels correspond to Γ = 100 and the bottom ones to Γ = 300. The blue curves correspond to the pion decay rate, and the green curves represent the πγ interaction with external and internal photons with aep = 1 (dashed lines) and with aep = 100 (solid lines). The pion synchrotron cooling is marked with orange curves for ϵB = 0.01 (short-dashed lines), ϵB = 0.1 (solid lines), and ϵB = 1 (long-dashed lines). The solid black curves mark the adiabatic cooling rate.

In the text

Fig. 9. Muon decay and cooling rates for L0 = 1050 erg s−1 (left panel) and L0 = 1051 erg s−1 (right panel). The top panels correspond to Γ = 100 and the bottom ones to Γ = 300. The blue curves correspond to the muon decay rate, and the green curves represent the IC interaction with external and internal photons with aep = 1 (dashed lines) and with aep = 100 (solid lines). The muon synchrotron cooling is marked with orange curves for ϵB = 0.01 (short-dashed lines), ϵB = 0.1 (solid lines), and ϵB = 1 (long-dashed lines). The solid black curves mark the adiabatic cooling rate.

In the text

	Fig. 10. Distributions of pions (top), muons (middle), and charged kaons (bottom) for Γ = 100. The left plots correspond to L0 = 1050 erg s−1 and the right plots to L0 = 1051 erg s−1. The green, orange, and blue curves refer to the cases of ϵB = 0.01, ϵB = 0.1, and ϵB = 1, respectively. The dashed curves correspond to aep = 1 and the solid curves to aep = 100.
In the text

	Fig. 11. Distributions of pions (top), muons (middle), and charged kaons (bottom) for Γ = 300. The left plots correspond to L0 = 1050 erg s−1 and the right plots to L0 = 1051 erg s−1. The green, orange, and blue curves refer to the cases of ϵB = 0.01, ϵB = 0.1, and ϵB = 1, respectively. The dashed curves correspond to aep = 1 and the solid ones to aep = 100.
In the text

Fig. 12. Emissivities of muon neutrinos and antineutrinos from the decay of pions (blue curves), muons (orange curves), and kaons (green curves) for L0 = 1050 erg s−1. The dashed curves represent the obtained output if all losses of the parent particles are neglected, while the solid curves mark the results if all the losses discussed are considered. The left plots correspond to aep = 1 and the right plots to aep = 100.

In the text

Fig. 13. Emissivities of muon neutrinos and antineutrinos from the decay of pions (blue curves), muons (orange curves), and kaons (green curves) for L0 = 1051 erg s−1. The dashed curves represent the obtained output if all losses of the parent particles are neglected, while the solid curves mark the results if all the losses discussed are considered. The left plots correspond to aep = 1 and the right plots to aep = 100.

In the text

Fig. 14. Diffuse flux of muon neutrinos and antimuons for L0 = 1050 erg s−1 (top panels) and L0 = 1051 erg s−1 (bottom panels). The results obtained with Γ = 100 and Γ = 300 are shown in the left and right panels, respectively. The green, orange, and blue curves refer to the cases of ϵB = 0.01, ϵB = 0.1, and ϵB = 1, respectively. The dashed curves correspond to aep = 1 and the solid ones to aep = 100.

In the text

Fig. 15. Energy dependent neutrino flavor ratios for L0 = 1050 erg s−1 (left panels) and L0 = 1051 erg s−1 (right panels) for Γ = 100. The top plots correspond to aep = 1 and the bottom ones to aep = 100. The green, orange, and blue curves correspond to the electron, muon, and tau flavor ratios, respectively. The long-dashed, solid, and short-dashed lines refer to the results for ϵB = 0.01, ϵB = 0.1, and ϵB = 1, respectively. The flavor ratios expected without considering losses are indicated with gray dotted lines for the electron flavor, gray long-dashed lines for the muon flavor, and gray long-short dashed lines for the tau flavor.

In the text

Fig. 16. Energy dependent neutrino flavor ratios for L0 = 1050 erg s−1 (left panels) and L0 = 1051 erg s−1 (right panels). Both plots are for Γ = 300. The green, orange, and blue curves correspond to the electron, muon, and tau flavor ratios, respectively. The long-dashed, solid, and short-dashed lines refer to the results for ϵB = 0.01, ϵB = 0.1, and ϵB = 1, respectively. The flavor ratios expected without considering losses are indicated with gray dotted lines for the electron flavor, gray long-dashed lines for the muon flavor, and gray long-short dashed lines for the tau flavor.

In the text

Fig. 17. Flavor ratios of the integrated fluxes of neutrinos with energy above Eν,min for L0 = 1050 erg s−1 (left panels) and L0 = 1051 erg s−1 (right panels). All plots are for Γ = 100. The top plots correspond to aep = 1 and the bottom ones to aep = 100. The green, orange, and blue curves correspond to the electron, muon, and tau flavor ratios, respectively. The long-dashed, solid, and short-dashed lines refer to the results for ϵB) = 0.01, ϵB = 0.1, and ϵB = 1, respectively. The flavor ratios expected without considering losses are indicated with gray dotted lines for the electron flavor, gray long-dashed lines for the muon flavor, and gray long-short dashed lines for the tau flavor.

In the text

Fig. 18. Flavor ratios of the integrated fluxes of neutrinos with energy above Eν,min for L0 = 1050erg s−1 (left panels) and L0 = 1051 erg s−1 (right panels). Both plots are for Γ = 100. The green, orange, and blue curves correspond to the electron, muon, and tau flavor ratios, respectively. The long-dashed, solid, and short-dashed lines refer to the results for ϵB) = 0.01, ϵB = 0.1, and ϵB = 1, respectively. The flavor ratios expected without considering losses are indicated with gray dotted lines for the electron flavor, gray long-dashed lines for the muon flavor, and gray long-short dashed lines for the tau flavor.

In the text