Neutrino Flavor Ratios Modified by Cosmic Ray Secondary-acceleration

Neutrino Flavor Ratios Modified by Cosmic Ray Secondary-acceleration

Norita Kawanaka norita@astron.s.u-tokyo.ac.jp Department of Astronomy, Graduate School of Science, University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo, 113-0033, Japan    Kunihito Ioka Theory Center, Institute of Particle and Nuclear Studies, KEK, Tsukuba 305-0801, Japan Department of Particle and Nuclear Physics, SOKENDAI (The Graduate University for Advanced Studies), Tsukuba 305-0801, Japan
July 14, 2019
Abstract

Acceleration of ’s and ’s modifies the flavor ratio at Earth (at astrophysical sources) of neutrinos produced by decay, , from () to () at high energy, because ’s decay more than ’s during secondary-acceleration. The neutrino spectrum accompanies a flat excess, differently from the case of energy losses. With the flavor spectra, we can probe timescales of cosmic-ray acceleration and shock dynamics. We obtain general solutions of convection-diffusion equations and apply to gamma-ray bursts, which may have the flavor modification at around PeV – EeV detectable by IceCube and next-generation experiments.

pacs:
13.85.Tp, 98.70.Sa, 98.38.Mz

I Introduction

The origin of high energy cosmic rays (CRs) has been a long-standing problem in astrophysics. Especially, CRs with energy above are considered to come from extra-Galactic sources such as active galactic nuclei (AGNs) and gamma-ray bursts (GRBs). In these sources we expect the production of high energy neutrinos () through interactions of accelerated protons with the ambient photons ( interactions) or gas ( or interactions) waxmanbahcall97 (); murase+06 (); muraseioka13 (). Detection of these neutrinos can provide us new information about high energy cosmic ray sources as well as the acceleration processes.

In cosmic ray accelerators, high energy neutrinos are mainly produced from the decay of charged pions: and . Therefore, the flavor ratios of these neutrinos are expected to be

(1)

at the sources, where denotes the flux of and (, or ). The observed flavor ratios become after the neutrino oscillations during the propagation to the Earth learnedpakvasa95 (). However, this argument may be too naive because we should take into account the finiteness of the decay timescale of pions and muons . For example, if the cooling timescale (kashtiwaxman05, ; lipari+07, ; tamborraando15, ) or acceleration timescale (murase+12, ; klein+13, ; winter+14, ; reynoso14, ) of a pion or a muon is shorter than the decay timescale, the spectral shape of neutrinos produced from the decay of those particles would be significantly modified. Especially, because the decay times are different between pions and muons, the energy dependence of neutrino fluxes would be different from flavor to flavor. The observed neutrino flavor ratio may be also modified by neutron decay anchordoqui+04 () and new physics such as neutrino decay beacom+03 (); baerwald+12 (), sterile neutrinos athar+00 (), pseudo-Dirac neutrinos beacom+04 (); esmaili10 (), Lorentz or CPT violation hooper+05 (), quantum-gravity anchordoqui+05 () and secret interactions of neutrinos iokamurase14 ().

Recently the high energy neutrino detector, IceCube, has discovered neutrinos aartsen+14 () which are confirmed to be non-atmospheric at the level of . The IceCube team has also analyzed the flavor composition of astrophysical neutrinos in the energy range of and demonstrate consistency with (although the best fit composition is ). In the near future, the next generation IceCube-Gen2 and KM3Net experiments will enable the precise study of the energy spectrum of high energy neutrinos and their flavor composition.

There is no study about how the flavor ratio of observed neutrinos as well as its energy dependence are modified by the acceleration of pions and muons, although several authors have investigated the neutrino energy spectrum under the secondary-acceleration murase+12 (); klein+13 (); winter+14 (); reynoso14 (). In addition, their approaches are based on one-zone klein+13 (); winter+14 () or two-zone reynoso14 () approximations, that is, they do not consider the spatial distribution of secondary particles (pions and muons) and their transport across the shock.

In this work, we investigate the acceleration of pions and muons produced by protons by solving their convection-diffusion equations around the shock front taking into account their decay into other particles (i.e., pions into muons and muon neutrinos, and muons into muon neutrinos and electron neutrinos). The shock acceleration of secondary particles has been discussed in the context of the positron excess blasi09 (); blasiserpico09 (); mertschsarkar09 (); ahlers+09 () observed by PAMELA/Fermi LAT/AMS-02 (adriani+11, ; ackermann+12, ; aguilar+13, ) (see also (ioka10, ; kawanaka+10, ; kawanaka12, )). We develop their formalism by including the decay of secondary particles during their acceleration, and evaluate the energy spectrum of neutrinos produced from those particles.

This paper is organized as follows. In Section 2 we formulate a general model for the shock acceleration of pions and muons, which are produced from shock-accelerated protons via photomeson interactions, and give versatile expressions of the neutrino spectra for later application. In Section 3, we show that the acceleration of pions/muons would be possible in low-power GRBs occurring inside their progenitors, and apply our model to that case to compute the neutrino spectra and flavor ratios. Our results and discussions are summarized in Section 4. In Appendix A, we summarize general solutions of the convection-diffusion equations for pions (decaying secondary particles) and in Appendix B for muons (decaying tertiary particles).

Ii Model

In this section, we describe the shock acceleration of secondary pions and muons that are generated from the protons accelerated at the shock, and investigate the energy spectrum of neutrinos produced by those pions and muons without specifying a particular source. Hereafter we neglect the energy loss of particles due to synchrotron emission and inverse Compton scattering for simplicity (see discussions in Sec. 4).

In the shock rest frame, the transport and shock acceleration of particles decaying into other kinds of particles with a timescale can be described by the convection-diffusion equation,

(2)

where () is the equilibrium distribution function of accelerated particles per unit spatial volume and per unit volume in momentum space, is the diffusion coefficient, and is the the velocity of the background fluid. We assume that the shock is non-relativistic and that the distribution functions are stationary and isotropic except for the shock front for simplicity. When the shock is relativistic, we should solve the relativistic version of the convection-diffusion equation taking into account the anisotropy of the particle momentum distribution. The anisotropy is the order of , and therefore in the mildly-relativistic shock, it would be less than a factor of two.

The shock front is set at , and the upstream (downstream) region corresponds to (). If we ignore the third term of the right-hand side, which describes the loss of particles due to their decay, this is a well-known equation for the usual diffusive shock acceleration of particles blandfordeichler87 (); malkovdrury01 (). The decay timescales of a pion and a muon are the functions of their energy,

(3)
(4)

where is the energy of a particle at the shock rest frame and, and are the decay timescales of a pion and a muon at their rest frames, respectively.

The fourth term of the right-hand side of Eq.(2), , is the distribution function of particles injected per unit time. Here we consider that charged pions are produced in interactions, . In this case, should be given from the distribution function of primary protons. On the other hand, in the case of muons, since they are produced by the decay of pions, should be given from the distribution of pions (see below).

Hereafter, we assume the Bohm-type diffusion,

(5)

where is the charge of a particle, is the magnetic field strength and is the gyrofactor, which is equal to unity in the Bohm limit drury83 (). The fluid velocity is given by

(6)

where and are constants and the compression ratio is .

One should solve Eq.(2) taking the following boundary conditions into account:

(7)
(8)
(9)

where (iii) comes from the integration of Eq.(2) across the shock front. This condition yields the differential equation for the distribution function at the shock front with respect of (see Appendix A).

The detail of the general solution of Eq. (2) is presented in the Appendix A. In the following subsections we briefly describe the properties of the derived distribution functions of pions and muons.

ii.1 pion acceleration

In this subsection we show the pion distribution function evaluated from Eq.(2). Pions are produced through the interactions between the protons accelerated in the shock and the ambient photons. The distribution of protons is given by

(10)

where is the proton distribution function at the shock front, which is proportional to . This expression (10) is a well-known solution of the convection-diffusion equation (2) blandfordeichler87 (); malkovdrury01 (); drury83 ().

For simplicity we assume that the ambient photon field is uniform. Since pions are produced from shock-accelerated protons, the production spectrum at the source of pions is proportional to that of primary protons, which can be described as

(11)

where is the momentum of a primary proton generating a secondary pion with a momentum , and these momenta are approximately related in a linear way:

(12)

where is the ratio of the energy of a pion to that of a primary proton waxmanbahcall97 ().

We can solve Eq. (2) for pions by using Eqs. (55), (56), and (57) with in Eq. (11), and obtain the pion distribution functions in the upstream and downstream as

(13)
(14)

The pion distribution functions at the shock front can be evaluated from Eqs. (59), (60) and (61) as

(15)

where and are numerical factors, being independent of (since both and are proportional to ):

(16)
(17)

where is the compression ratio. We can see that the distribution function of pions at the shock front, Eq. (15), becomes harder than their production spectrum by [], and this is similar to Eq.(6) of blasi09 (), where the acceleration of secondary positrons produced in the supernova remnant shock is discussed. The difference is that we take into account the decay of particles, the third term of the right-hand side of Eq. (2), in the convection-diffusion equation of secondary particles, which is reflected as numerical factors and .

To see the effects of the transport and acceleration of pions in the downstream region from Eq. (14), we divide into two components: and . The former component represents the pions that are reaccelerated at the shock, being proportional to . On the other hand, the latter component represents the pions that are produced from the protons and advected in the downstream region, being proportional to :

(18)
(19)

Hereafter, since the number of upstream pions, , is subdominant compared to that of downstream pions, , we only discuss the contribution of the downstream pions to the neutrino spectra.

In the limit of (i.e., the lifetime of a pion is much longer than the acceleration timescale, ), Eq.(15) has the same form as Eq.(6) of blasi09 (). In this limit, we have , . Therefore, the distribution function at the shock front in Eq.(15) can be approximated as

(20)

where is the power-law index of the production spectrum, , and in the strong shock limit with an adiabatic index . We can see that, since we assume , the resulting spectrum is proportional to , being harder than the production spectrum. This can be interpreted as the result of the secondary-acceleration: since pions produced from shock-accelerated protons can cross the shock front before their decay, they would gain the energy and their spectrum would become harder. We should also note that is proportional to . In this limit of , the pion distribution function in the downstream region (14) can be approximated as

(21)

where we use . Here the first term corresponds to the pions reaccelerated at the shock, and its damping length scale is identical to the distance over which a pion is advected with the fluid during its lifetime. The second term represents the pions that are produced from protons in the downstream region and simply advected further downward. On the other hand, in the upstream region, Eq.(13) can be approximated as

(22)

where we use .

ii.2 muon acceleration

Using the results of the last subsection, we can evaluate the distribution function of muons that are produced from the decay of pions. The production spectrum at the source of muons should be given based on the distribution function of pions as follows:

(23)

where is the ratio of the momentum of a muon to that of a primary pion . Using the method shown in the Appendix A, we can solve the muon convection-diffusion equation and derive . The detail of the solutions is shown in the Appendix B.

In the limit of (i.e., the lifetime of a pion is much longer than the acceleration timescale), the muon distribution function at the shock front can be approximated as

(24)

where we assume the power-law spectrum of pion production, . Since and , we can see that this spectrum Eq.(24) is proportional to , which is similar to the pion distribution function at the shock shown in Eq. (20). This can be interpreted as follows. As stated in Eq.(23), the production spectrum of muons is proportional to . Since the lifetime of a pion is proportional to , the production spectrum at the shock is proportional to . Injected muons are reaccelerated at the shock and, according to the similar discussion to that on pions, the muon spectrum at the shock would become harder than their injected spectrum by , which comes from the dependence of on . We should also note that is proportional to .

In a similar way to the pion distribution function, the muon distribution function in the downstream can be divided into two components, and as follows:

(25)
(26)

and, as in the case of pions, the number of upstream muons, is subdominant compared to that of downstream muons .

ii.3 neutrino spectra

From the pion and muon distribution functions calculated above, the neutrino spectrum can be obtained as follows:

(27)
(28)
(29)

where , and are the ratios of the energy of a muon neutrino, an anti-muon neutrino, and an electron neutrino to that of their primary particles, respectively. Since each lepton produced from the decay of a pion (, , and ) carries approximately equal energy (i.e., 1/4 of that of the primary pion), we set , and . The volume integral should contain the surface area integral on the shocked matter plus the integral along the normal direction of the shock. Especially, defining the dynamical timescale as time for the shock to cross the system, the latter integral should be from to .

We divide the neutrino energy spectrum into two components according to the decomposition of the pion/muon distribution functions shown in Eqs. (18), (19), (25) and (26):

(30)
(31)

and and are defined in similar ways.

We should also consider neutrino oscillations during the propagation from the source to the Earth. When neutrinos propagate over the distances much longer than ( is the squared mass difference: , ), the observed fluxes of neutrinos () should be described as

(32)

where is the neutrino mixing matrix and the subscript represents the mass eigenstate of neutrinos. The matrix elements of can be described by the mixing angles , , and , and the Dirac phase . Based on fogli+12 (), we adopt , , , and .

Iii Applications to Low-Power GRBs

Now we consider long GRBs as neutrino sources. GRBs are thought to produce high-energy neutrinos waxmanbahcall97 (). In the standard model, the emission of long GRBs is believed to be produced by relativistic jets launched when a massive star collapses and a stellar-mass black hole is formed. In order for the jet to be observed as a GRB, it should penetrate the stellar envelope, otherwise the jet would stall inside the star and the gamma-ray emission would not be observed reesmeszaros94 (). Their prompt emission is often interpreted as synchrotron emission from non-thermal electrons accelerated at internal shocks. It is natural to consider the proton acceleration and the associated production of high energy neutrinos via / interactions waxmanbahcall97 ().

However, IceCube gave stringent upper limits on GRBs abbasi+11 (); abbasi+12 () and has ruled out the typical long GRBs as the main source of the observed diffuse neutrino events hummer+12 (); he+12 (); gao+13 ().

Instead of ordinary GRBs, we investigate high-energy neutrino production by low-power GRBs such as low-luminosity GRBs (LLGRBs) and ultralong GRBs (ULGRBs), which are still not strongly constrained by IceCube. In these low-power GRBs the high energy neutrinos may be produced inside the progenitor star meszaroswaxman01 (); razzaque+04 (); horiuchiando08 (). While a jet is penetrating in a stellar envelope, it becomes slow and cylindrical by passing through the collimation shock. The internal shocks would also occur when there is spatial inhomogeneity in a jet. Murase & Ioka muraseioka13 () recently investigated such high energy neutrino production expected from LLGRBs soderberg+06 () and ULGRBs gendre+13 (); levan+14 (), which have longer durations () and lower luminosity () compared to those of typical long GRBs. It has been suggested that ultra-long GRBs have bigger progenitors like blue supergiants (BSGs) with radii of suwaioka11 (); kashiyama+13 (); nakauchi+13 (). We apply our model of neutrino production in such GRB jets inside stars, taking into account the secondary-acceleration and decay of pions/muons that are produced by shock-accelerated protons via interactions. The internal shocks of GRBs are considered to be mildly-relativistic in the shock rest frame.

Let us evaluate the important timescales in our model by considering the internal shock scenario of GRBs. When two moving shells are ejected with comparable Lorentz factors of order of from the central engine during the time separation , these shells collide and make an internal shock at the radius . Here the magnetic field energy density can be estimated as , where is the magnetic luminosity. Then we can estimate the acceleration timescale , synchrotron cooling timescale as functions of the energy of a particle, and the dynamical timescale at the shock rest frame as follows:

(33)
(34)
(35)
(36)

where is the energy of a particle ( or ) at the shock rest frame, , , , and are the masses of a charged pion and a muon, respectively.

From Eqs. (3) and (4), a pion can be accelerated at the source before its decay when , i.e.,

(37)

while a muon can be accelerated before its decay when

(38)

Note that, since both and () are proportional to the energy of a pion (a muon), these conditions are independent of the energy of particles. Under these conditions, pions (muons) can be accelerated at the shock before they decay and therefore their spectra would become harder. On the other hand, we can see that the synchrotron cooling timescale would be shorter than the acceleration timescale when the energy in the shock rest frame is higher than , where

(39)
(40)

In order to evaluate the timescales of inverse Compton cooling and interaction, we should give the target photon spectrum at the local rest frame. In the case of the internal shock occurring inside a star, the accelerated particles mainly interact with photons that are produced in the jet head and escape back from there. Here we estimate the spectrum of the target photon field according to the procedure adopted in muraseioka13 (). At the head of the collimated jet, the photon temperature is given as

(41)

where is the total jet luminosity ( is the fraction of the magnetic energy), is the Boltzmann constant, is the Stefan-Boltzmann constant, is the radius where a jet becomes cylindrical through the collimation shock, and is the Lorentz factor of the collimated jet (note that this is different from the Lorentz factor of the precollimated jet ). The fraction of photons escaping the collimated jet is where is the comoving proton number density in the collimated jet, and is the Thomson cross section. Therefore, the number density of the target photons is given as

(42)

where is the comoving photon number density in the collimated jet, and is the Riemann zeta function. We assume that the escaping photon field has a thermal spectrum,

(43)

with the effective temperature of .

The photomeson production ( interaction) timescale can be evaluated as

(44)

where , is the cross section of pion production as a function of photon energy in the proton rest frame, is the average fraction of energy lost from a proton to a pion, and is the threshold energy (waxmanbahcall97, ). In the following discussion, we use the resonance approximation: is approximated to be a function with a peak at , where with the width of , and .

Figures 1 and 2 depict the acceleration timescales, cooling timescales via synchrotron emission and inverse Compton scattering, decay timescales of a pion and a muon, the timescale of interactions, and the dynamical timescale () in the internal shock occurring inside a star expected for an ultralong GRB (, , , , ). We can see that, with the current choice of parameters, the acceleration timescales of a pion and a muon are shorter than their lifetimes for arbitrary energy range, and that the decay timescale becomes longer than the dynamical timescale above the energy of for pions and for muons (at the shock rest frame). Note that, since synchrotron cooling timescale for pions and muons becomes shorter than acceleration timescale when the energy of particles is larger than , our formalism is not applicable in the energy range above . Note also that the efficiencies of pion/muon production would be suppressed in the energy range where the timescale of interactions is comparable () or shorter than the acceleration timescale. In the current work, this effect is not taken into account.

Figure 3 depicts the energy spectra of muon neutrinos and electron neutrinos expected from the internal shock inside a progenitor of an ultralong GRB. The flux from reaccelerated pions and muons, Eqs. (18) and (25), and that from advected pions and muons, Eqs. (19) and (26), are also shown (with dotted lines and dashed lines, respectively). We can see that the electron neutrino flux from advected muons drops above the energy of in the observer frame ( at the shock rest frame). This corresponds to the energy at which the muon decay timescale is equal to the dynamical timescale. Above this energy, only part of muons can decay into within the dynamical timescale muonadcool (). As for the muon neutrino flux from advected particles, it slightly drops around the energy where the electron neutrino flux drops because anti-muon neutrinos are generated from the decay of muons , and drops again at the energy where the pion decay timescale is equal to the dynamical timescale ( for the current parameter set) because muon neutrinos are generated from the decay of pions . We can interpret this behavior as follows. From Eqs. (19), (26), (27), (28) and (29), in the limit of and , the neutrino fluxes from advected particles can be approximated as

(45)
(46)

while in the limit of and they can be approximated as

(47)
(48)

where is the volume of the merged shell making the internal shock. Here we neglect the contribution from pions/muons in the upstream region () because it is subdominant compared to that from the downstream pions/muons. We can easily see that in the latter limit , the energy spectra of neutrino fluxes are softer than by because the decay timescale is proportional to .

On the other hand, the neutrino fluxes from reaccelerated pions/muons increase more as the acceleration timescales become longer. Under the condition , from Eqs. (18), (25), (27), (28) and (29), we can approximate the neutrino fluxes from reaccelerated pions/muons as

(49)
(50)

In the high energy limit, where the decay timescales of pions/muons are much longer than the dynamical timescale, each of these neutrino fluxes behaves asymptotically as

(51)