Neutrino production in Population III microquasars
NNeutrino production in Population III microquasars
Agust´ın M. Carulli a , Mat´ıas M. Reynoso a , Gustavo E. Romero b,c a IFIMAR (CONICET-UNMdP) and Departamento de F´ısica, Facultad de Ciencias Exactas y Naturales, Universidad Nacional de Mar del Plata, Funes 3350,(7600) Mar del Plata, Argentina b Instituto Argentino de Radioastronom´ıa (CCT-La Plata, CONICET; CICPBA)C.C. No. 5, (1894) Villa Elisa, Argentina c Facultad de Ciencias Astron´omicas y Geof´ısicas, Universidad Nacional de La PlataPaseo del Bosque s / n, (B1900FWA) La Plata, Argentina. Abstract
Microquasars (MQs) are binary systems composed by a star feeding mass to a compact object through an accretion disk. Thecompact object, usually a black hole, launches oppositely directed jets which are typically observed in our galaxy through theirbroadband electromagnetic emission. These jets are considered potential galactic neutrino sources. MQs can also have been formedby the first generations of stars in the universe, i.e., Population III (Pop III) stars, which are considered essential contributors tothe ionization processes that took place during the period of “cosmic reionization”. In the present work, we develop a model thataccounts for the main particle processes occurring within Pop III MQ jets, with the aim to obtain the di ff use neutrino flux at theEarth. We define di ff erent zones within the jets of Pop III MQs where particle interactions occur, and primary particles (i.e protonsand electrons) are injected. We solve a transport equation for each zone, including the relevant cooling and escape processes, whichinclude p γ and pp interactions. Once we obtain the primary particle distributions, we compute the pion and muon distributions, aswell as the neutrino output produced by their decays. Finally, we obtain the di ff use neutrino flux by integration over the redshift,the line-of-sight angle, and the MQs lifetime. We find that, for a range of parameters suitable for Pop III MQ jets, the most relevantsite for neutrino production in the jets is the base of the inner conical jet. Additionally, if protons accelerated at the forward shockformed at terminal jet region can escape from the outer shell, they would produce further neutrinos via p γ interactions with thecosmic microwave background (CMB). The latter contribution to the di ff use neutrino flux turns out to be dominant in the range10 GeV (cid:46) E ν (cid:46) GeV, while the neutrinos produced in the inner jet could only account for a small fraction of the IceCube fluxfor E ν ∼ GeV. The co-produced multiwavelength photon background is also computed and it is checked to be in agreementwith observations.
Keywords:
Radiation mechanisms: non-thermal – Neutrinos; X-rays: binaries
1. Introduction
The first stars formed in the early universe ( z ∼ − ∼ M (cid:12) − M (cid:12) according to star formation simulations [1].Their low metallicity and mass imply that nearly ∼
90% of themcollapsed to black holes. Moreover, theoretical results suggestthat these stars often formed binary systems ( ∼
50% of themaccording to Ref. [2] ). In this context, it is expected that Pop-ulation III Microquasars (Pop III MQs) can indeed have beengenerated [3, 4]. Composed by a black hole ( ∼ M (cid:12) in thiswork) and a Pop III companion star, Pop III MQs are expectedto be super-acreeting systems, i.e., much more powerful thantypical galactic MQs [4]. Two of the latter have been detectedat high energy gamma rays (Cyg X-1 and Cyg X-3). In addition,the MQ SS433, which is the only super-accreting binary in theGalaxy, was also found to emit very high energy gamma rays[5]. Given their nature, MQs have long been considered as po-tential high energy neutrino sources [6, 7, 8, 9, 10, 11, 12, 13].Pop III MQs are also considered as one of the possiblemain contributors to the process of cosmic reionization [14, 4]. Around 0 .
37 Myr after the Big Bang, a period known as re-combination took place: the plasma of electrons and protonswent through a phase transition that coupled them together forthe first time to form neutral atomic hydrogen. Therefore, anera named as the “Dark Ages” arose, during which no ob-jects capable of producing radiation had yet been formed. Af-terwards, the universe passed through a period known as the“Epoch of Reionization”, which began thanks to the ultravi-olet (UV) radiation produced by the first formed stars. Thisradiation might have been capable of ionizing the intergalac-tic medium (IGM) within the boundaries of the haloes wherethe stars were born. However, the radiation emanated from thistype of sources could not ionize farther away, meaning that an-other ionization mechanism should have taken place in orderto explain the ionization at longer distances. As proposed byRef. [14], X-ray radiation produced by jets arising from accret-ing BH is essential for the ionization process, since the meanfree path of X-rays is much longer than the corresponding toUV radiation. Under these considerations, the reionization ca-pabilities of Pop III MQs jets have also been studied taking intoaccount more complete models [4].
Preprint submitted to Astroparticle Physics January 11, 2021 a r X i v : . [ a s t r o - ph . H E ] J a n ince Pop III MQs would have been formed at high red-shifts, gamma rays of very high energy ( E (cid:38)
100 GeV) thatcould have been produced in their jets would have been ab-sorbed during their propagation to Earth. On the other hand,neutrinos would be completely unabsorbed, and hence, it isinteresting to compute the neutrino output from these sourcesin order to assess their detection with neutrino telescopes suchas IceCube. Previous works have addressed the possible neu-trino emission arising from other phenomena related to PopIII stars, such as supernova remnants [15, 16] and gamma-raybursts [17, 18]. Here, we concentrate on a neutrino contributionwhich has not been previously considered, i.e., the correspond-ing to the MQ evolutionary phase of binary systems composedby Pop III stars.Specifically, we develop a model that accounts for the mainparticle processes occurring within Pop III MQs jets, with thepurpose of obtaining a contribution to the di ff use neutrino fluxthat would arrive on the Earth. To do so, we take into ac-count contributions of neutrino production via the decay of pi-ons and muons resulting from pp collisions between high en-ergy protons and cold proton targets, p γ interactions betweenhigh energy protons and soft photons produced by electron syn-chrotron, plus protons interacting with photons from the CMB.We compare the resulting neutrino fluxes with the best fits avail-able for the di ff use flux of astrophysical neutrinos obtained ex-perimentally by IceCube [19, 20]. We also compare the resultswith the upper limit from the Pierre Auger Observatory [21]and the expected sensitivity for GRAND [22], which wouldbe sensible to neutrinos of higher energies (10 GeV (cid:46) E ν (cid:46) GeV).This work is organized as follows: in the next section, wedescribe the jet model and the calculation procedure appliedobtain to particle distributions. In the following section, wepresent the di ff use neutrino flux for di ff erent combinations ofparameters, and we also obtain the accompanying flux of mul-tiwavelength photons. Finally, in the last section, we discussthe results and give our concluding remarks.
2. The model
The present model is based on the one presented in Ref.[4], where the companion star provides the mass that is trans-ferred to the central back hole through an accretion disk inan extremely super-critical regime. The critical accretion rateof the black hole is associated to the Eddington luminosity L edd ( ˙ M crit = L edd / c ) and thus can be written as ˙ M crit = π GM BH m p / ( σ T c ), where M BH is the mass of the black hole, m p is the proton mass, and σ T is the Thomson scattering cross-section. The Eddington luminosity is such that for sphericallysymmetric accretion, the radiation pressure would exactly bal-ance the e ff ect of gravity, implying that higher accretion ratesare not attainable. However, if accretion proceeds through adisk as, then it is possible that the emission is directed mostlyperpendicularly to the disk plane, thus not cancelling to theaccretion process [23]. This means that it is perfectly pos-sible to have super-critical systems with accretion rates muchhigher than the critical one, as it occurs in the mentioned case of SS433. Is is also considered that this regime is appropriatefor Pop III MQs due to the large amount of mass accreted bythe black hole through overflow of the Roche lobe [4]. The MQphase lasts τ MQ ∼ × yr according to Ref .[24], and we willtake it into account in order to compute the total neutrino emis-sion along the whole life of the MQs. It is also expected thatunder this regime of accretion, a large fraction of the accretedmaterial has to be ejected in powerful winds and jets that canreach a kinetic power as high as L k ∼ erg s − . Therefore,Pop III MQs are sources expected to be significantly more pow-erful than their galactic counterparts. Since their formation andexistence is well-motivated, we consider that it is worth exam-ining their possible neutrino production potential. For details ofthe accretion disk, the reader is referred to Ref. [4], while herewe concentrate on the mechanisms of neutrino generation, aswell as the associated multiwavelength photons resulting fromparticle acceleration in the jets of Pop III MQs. Figure 1: Schematic view of a Population III microquasar. The black hole ac-cretes material from the Population III companion star. The jets arise oppositelydirected from the plane containing the accretion disk, and form an angle i j withthe line of sight. Relevant zones, such as the cocoon and the bow shock, arealso shown. The inner microquasar jets are modeled as two oppositely di-rected outflows arising from the vicinity of the central blackhole (BH) of mass M BH (cid:39) M (cid:12) . In Fig. 1, we show anschematic view of the main components of the system. Wesuppose that each jet is accelerated through the conversion ofmagnetic to bulk kinetic energy, according to the basic mecha-nism discussed by Ref. [25]. In this scenario, the flow is initialymagneticaly dominated in the vicinity of the BH and subjectto di ff erential collimation decreasing along the jet. It can beshown that this leads to an acceleration of the flow, increasingits Lorentz factor (see Ref. [25] for details). Energy equiparti-tion is expected to hold at an inner position z ∼ (10 − R g ,where R g = GM BH / c is the gravitational radius of the accretor.For greater distances along the jet ( z j > z ), the magnetic energydrops as the jet gradually accelerates, and we assume that the jetreaches its final Lorentz factor Γ at a distance z acc (cid:29) z . At this2oint, we consider that the jet becomes conical, and thus sub-ject to uniform collimation, which prevents further magneticacceleration of the flow. The magnetic energy at z acc is thensupposed to be a small fraction q m of the kinetic energy, i.e.,the jet becomes matter-dominated, which is a necessary con-dition for the development of shocks. Hence, the first regionof particle acceleration that we consider is placed at z acc witha size ∆ z b ∼ (1 − R j , where R j = z acc tan ψ is the radius ofthe jet and ψ its half-opening angle. The jet continues its prop-agation with a constant velocity, and this is consistent with amagnetic field dependence on the distance along the jet ( z j ) asRefs. [25, 26, 4] B ( z j ) = B acc (cid:32) z acc z j (cid:33) , (1)where B acc is the magnetic field at z acc . Following Ref. [27], weconsider further possible zones throughout the jets where somemechanism of particle acceleration can take place [37, 38, 39].The inner jets propagate expanding laterally up to the recon-finement point z rec , where the pressure in the jets equals that ofthe external medium. Beyond the reconfinement point, the jetscontinue their propagation keeping a constant radius until theyreach the terminal regions where they are finally stopped by theexternal medium.Since our goal is to obtain the total contribution of neutrinosfrom Pop III MQs to the di ff use neutrino flux, this is to be foundby adding up the individual contributions from the followingdi ff erent zones where emission can take place (see Fig.2 for anschematic view): • Base zone • Conical jet • Reconfinement zone • Cocoon • Shell or bow shock • External regionThe first zone we consider, where accelerated particles areinjected, is close to the base of the jet (base zone, for short),and it is placed in the inner jet at the distance z acc from the BH.There, the kinetic energy density is in sub-partition with themagnetic energy density, i.e. ρ m = q m ρ k . This means that thekinetic energy is smaller by a fraction q m than it would be underequipartition. The jet kinetic density is: ρ k = L k [( Γ − Γ π R v j ] , (2)where L k is the jet kinetic power and v j is its bulk velocity.The particles that escape from the base zone are injected intoa larger conical region at the inner jet, which extends up tothe reconfinement point. In order to account for the propaga-tion e ff ect, a convection term is included in the correspondingtransport equation, as we discuss below. At the terminal re-gion, the jets interact with the external medium and both a for-ward shock (bow shock) and a reverse shock are produced. The former propagates into the IGM forming a shell, and a reverseshock that goes inward the jet creating a region called cocoon,which exerts a pressure directed towards the jet. This pressure isable to stop the lateral expansion of the jet and forces the cone-shaped part of the jet to become a cylinder for larger distances.The third emission zone we consider is the reconfinement zone,placed at a distance to the BH that can be computed as Refs.[28, 27] z rec ∼ (cid:118)(cid:116) L k v j ( γ + (cid:16) Γ j − (cid:17) π c P coc (3)where γ = / v j is the jet velocity, and P coc is the pressure by the cocoon: P coc = m p n IGM ( z ) v . (4)In the latter expression, n IGM is the number density of IGM mat-ter, n IGM ( z ) = H π G m p Ω M (1 + z ) (5)where Ω M is the is the matter density and H is the Hubbleconstant, both corresponding to z =
0. The velocity of the bowshock is v bs = l bs t MQ . (6)In turn, the distance from the BH to the bow shock is given by l bs = (cid:32) L k m p n IGM (cid:33) t MQ (7)where t MQ is the microquasar age and L k is the jet power.The remaining two emission zones considered in each MQ jetare placed at the cocoon and the bow shock (shell). However,we still consider that the protons escaping from the shell areinjected outside the system in an external zone, where they canfurther interact with the CMB. Here we discuss the main cooling processes that a ff ect therelativistic particles at the zones mentioned above. We considera density of cold protons in each zone of the model given by: n p = ρ k m p c , (8)where ρ k is the kinetic energy density given by Eq.2. These coldprotons are considered to be targets for the relativistic protons,and the corresponding rate of pp interactions is [29]: t − pp ( E p ) = n p c σ (inel) pp ( E p ) K pp , (9)where the inelasticity coe ffi cient is K pp ≈ / σ (inel) pp is taken to be as in Ref. [30].For each particle type considered (electrons, protons, pions,and muons), the synchrotron cooling rate is, t − ( E i ) = (cid:32) m e m i (cid:33) σ T B m e c π E i m i c , (10)3 igure 2: Emission zones considered in the jets model. m e is the mass of an electron, and m i and E i are the massand the energy of particles of type “ i ”, respectively. The parti-cle types considered are electrons ( i = e ), protons ( i = p ), pions( i = π ), and muons ( i = µ ). B is the magnetic field correspond-ing to each of the zones.The synchrotron emission of electrons with an energy distri-bution N e produces a background density of photons which canbe approximated by the following expression in the comovingreference frame: n ph ( E ph ) = (cid:15) syn ( E ph ) E ph R j c , (11)where E ph is the photon energy, and (cid:15) syn is the power per unitvolume per unit energy of the photons, (cid:15) syn ( E ph ) = (cid:32) − e − τ SSA ( E ph )) τ SSA ( E ph ) (cid:33) (cid:90) ∞ m e c dE π P syn N e ( E ) . (12)Here, P syn ( E ph , E ) is the power per unit energy of synchrotronphotons with energy E ph emitted by an electron of energy E . Itis defined as [31, 33]: P syn ( E ph , E ) = √ e Bm e c h E ph E cr (cid:90) ∞ E ph E cr d ζ K / ( ζ ) (13)where K / ( ζ ) is the modified Bessel function of order 5 / E cr = √ heB π m e c (cid:32) Em e c (cid:33) . (14)The e ff ect of synchrotron self-absorption (SSA) is taken intoaccount with the factor between parentheses in Eq. (12). Thiscorrects the synchrotron emissivity by accounting for the possi-bility that low energy synchrotron photons may be reabsorbedby electrons [32]. The corresponding optical depth τ SSA isgiven by τ SSA ( E ph ) = (cid:90) ∞ m e c dE α SSA ( E ph ) , (15)where the SSA coe ffi cient is [33]: α SSA ( E ph ) = c h π E (cid:90) ∞ m e c dEP syn ( E ph , E ) E × (cid:34) N e ( E − E ph )( E − E ph ) − N e ( E ) E (cid:35) . (16)Inverse Compton (IC) interactions of relativistic electronswith soft photons are considered in the model, in particular,the CMB photons are the most relevant target for the termi-nal jet zones, while for the electrons in the base zone, thesynchrotron emission of the electrons themselves become thedominant target, giving rise to the so-called synchrotron self-Compton (SSC) process. In order to obtain t − , we apply asuccessive approximation method as we discuss below, since itis necessary to know the particle distribution of the synchrotron emitting electrons: t − ( E e ) = m e c σ T E e (cid:90) E e E (min)ph dE ph n ph ( E ph ) E ph × (cid:90) Γ e Γ e + E e E ph dE γ F ( q ) (cid:104) E γ − E ph (cid:105) . (17)Here, E (min)ph is the lowest energy of the available background ofphotons produced by synchrotron of electrons, q = E γ ( Γ e ( E e − E γ )), with Γ e = E ph E e / ( m e c ), and the function F ( q ) is ob-tained following Ref. [31]: F ( q ) = q ln q + (1 + q )(1 − q ) +
12 (1 − q ) ( q Γ e ) + Γ e q . (18)The cooling of relativistic protons by p γ interactions is givenby Ref. [34]: t − p γ ( γ p ) = (cid:90) ∞ (cid:15) th mpc Ep d (cid:15) cn ph ( (cid:15) ) m p c E p (cid:15) (cid:90) (cid:15) Epmpc (cid:15) th d (cid:15) r σ p γ ( (cid:15) r ) K p γ ( (cid:15) r ) (cid:15) r , (19)where (cid:15) th ≈
150 MeV, σ p γ is the inelastic cross-section for pho-topion and photopair creation, K p γ is the inelasticity coe ffi cient(taken as in Ref. [34]), and n ph ( E ) represents the density oftarget photons.Electrons can also be cooled down by Bremsstrahlung,though it is small compared with the other processes for theparameters considered: t − ( E e ) = α FS r e c n p ln (cid:32) E e m e c − (cid:33) . (20)The adiabatic cooling rate for a gas of relativistic particles inan expanding volume at a rate dV / dt is [35]: t − ( E ) = V dVdt . (21)For a conical jet at a distance z j from the central source, consid-ering an element of volume V = π R dz j and a lateral expansionvelocity as dR j / dt = v j tan ψ , yields [36]: t − = β cz j . (22)where β is the bulk velocity of the jet at that position in units of c . In the case of spherical expansion (as considered in the bowshock) is: t − = β c ∆ z j , (23)where ∆ z j is the size of the zone considered.
3. Relativistic particles at the di ff erent zones ymbol Description Base Reconfinement Cocoon Shell Units L p , j + L e , j Power injected 10 erg s − Γ j Lorentz factor 1 . −
10 1 . −
10 1 . − < .
003 1 q m Magnetic parameter 5 × − . . . α Injection index 1 . − . . − . . − . . − . R j Radius of emitter 4 . × − . × . × . × . × cm z j Injection point 4 . × − . × . × . × . × cm B Magnetic field 2 × − . × . × − × − × − GTable 1: Main parameters of the model for the four zones where primary particles are injected as a power-law in the energy. The conical jet and the external zones arenot listed because injection there is determined by the escaping particles from the base and from the shell, respectively. The values for the reconfinement, coccoon,and shell correspond to z = t MQ = . × yr. As mentioned above, populations of relativistic primary par-ticles (electrons and protons) can be accelerated to very highenergies by some mechanism such as shock acceleration. Ateach zone ” j ”, the power injected in the form of relativisticparticles ( L e , j + L p , j ) is taken to be fraction q rel of the total jetkinetic power ( L k ). We refer to such particles as primary, be-cause these are ones that can initiate the radiation and emissionprocesses that give rise to the production of other, secondaryparticles which include photons, pions, muons, and neutrinos.We consider a steady-state one-zone treatment where theemission region is spatially homogeneous, and the injection andcooling rates are independent of time. This approach is appliedto the base zone, reconfinement region, cocoon, shell, and alsoto the external zone considered. Instead, for the extended con-ical part of the inner jet, we apply an inhomogeneous transportequation which also accounts for the convection e ff ect, as dis-cussed below.In the terminal regions (reconfinement, cocoon, and shell),the jet becomes a ff ected by the IGM, which first makes thelateral expansion to cease, and ultimately stops the jet propa-gation. At the reconfinement point z rec given by Eq. (3), rec-ollimation shocks can give rise to further particle acceleration(see e.g. Ref. [27]). The accelerated particles are injected into acylinder-shaped zone which extends up to the position of the re-verse shock, where the cocoon forms ( z j (cid:39) l bs ). As mentioned,this zone corresponds to the reconfined jet, which has a size ∆ z cyl = l bs − z rec and a radius R rec = z rec tan ψ , and where addi-tional pp , p γ interactions can take place. Both the bow shockand the reverse shock can generate particle acceleration, and thecorresponding emission zones are the shell and the cocoon, re-spectively. Their radius are l bs / ∆ z bs (cid:39) R bs , that of the cocoon is ∆ z coc = R rec .In Table 3, we present typical values of the main parameters ofour model adopted throughout the work. In the emission zones placed at the jet base, reconfinementregion, cocoon, and shell, primary particles are injected as apower law in the energy at the comoving frame, for energiesgreater than E i , min = m i c : Q i , j ( E i ) = K i , j E − α i e ( − E i / E max , i , j ) . (24) Here, “ i ” refers to electrons ( i = e ) and protons ( i = p ), α is theinjection index and K i , j is a constant fixed by normalization onthe total power injected in electrons and protons, L i , j = π ∆ V j , com (cid:90) ∞ E i , min dE i E Q i , j ( E i ) . (25)This expression is applied in the comoving reference, E i , min isthe minimum energy of injection, ∆ V j , com = Γ∆ V j is the comov-ing volume of the corresponding zone, and ∆ V j is the Lorentzcontracted volume as seen from the BH frame. The maximumenergies E max , i , j , which appear in the exponential cut-o ff , cor-respond, in principle, to the balance energies for which thetotal rate of cooling plus escape is equal to the accelerationrate. The latter is obtained by applying the general requirementthat the timescale for energy gain is greater than r gyr / c , where r gyr = E i / ( e B ) is the gyroradius for a particle with energy E i and charge e . Therefore, the acceleration rate is expressed as t − ( E i ) = η ecBE i , (26)where η < ffi ciency coe ffi cient that depends on the de-tails of the acceleration mechanism [29]. One further require-ment is that the particles can only remain confined inside thezone if their gyroradius does not exceed the size of the acceler-ation region. Thus, it must be fulfilled that E i / eB ( z j ) < R j ( z j ),which is known as Hillas criterion, E H = eB ( z j ) R j ( z j ). There-fore, if the balance energy E max , i mentioned above happens tobe higher than E H , then we simply set E max , i = E H . We show inFig. 3 the proton and electron cooling rates at the base zone andat the shell obtained for redshift z = t MQ ∼ . × yr .We choose as representative values ∆ z j = R j , α = Γ = . q m = × − , which are included in the ranges indi-cated in Table 3. We also show in the figure two di ff erent casesfor the escape rates. One is determined by a constant escapetimescale in the comoving frame, T esc (cid:39) Γ j ∆ z / v j , where Γ j isthe Lorentz factor of the zone considered and v j is its veloc-ity. The other case considered corresponds to a Bohm di ff usiontimescale, T B ( E ) = ( ∆ z ) / [2 D B ( E i )], where the di ff usion coef-ficient is D B ( E i ) = r gyr c /
3, so that T B ( E i ) = e B ( ∆ z ) E i c . t - [ s - ] −6 −4 −2 E e [GeV] −2 −1 electrons [base zone] a cc e l e r a t i o n s y n c h r o t r o n s y n c h r o t r o n SSCconst. escapeadiabatic b r e m ss t r a h l un g B o h m e s c a p e t - [ s - ] −18 −16 −14 −12 −10 −8 −6 −4 E p [GeV] protons [shell] a cc e l e r a t i o n const. escapeadiabaticpp p γ B o h m e s c a p e t - [ s - ] −18 −16 −14 −12 −10 −8 −6 −4 −2 E e [GeV] −2 −1 electrons [shell] a cc e l e r a t i o n s y n c h r o t r o n s y n c h r o t r o n SSC const. escapeadiabaticIC B o h m e s c a p e b r e m ss t r a h l un g t - [ s - ] −6 −5 −4 −3 −2 −1 E p [GeV] protons [base zone] a cc e l e r a t i o n const. escapeadiabaticppp γ B o h m e s c a p e Figure 3: Proton (electron) cooling rates for the base zone and bow shock are shown on the left (right) panels. The bow shock rates correspond to a redshift z = t MQ ∼ . × yr . We adopt ∆ z j = R j , α = Γ = . η = .
1, and q m = × − .
7e then calculate the distributions of primary particles N i , j in the zone “ j ” by solving the steady-state transport equation: d (cid:104) b i N i , j ( E i ) (cid:105) dE + N i , j ( E i ) T esc = Q i , j ( E i ) , (27)where i = e stands for electrons and i = p for protons. b i ≡ dE i / dt = − E i t − embody the continuous energy lossesof the particles due to the cooling processes that occur in thezone, i.e., synchrotron, IC (or SSC), adiabatic expansion, pp and p γ interactions. In the case of the base zone, we performsuccessive approximations in order to obtain the electron dis-tribution N e , b , since the SSC cooling rate cannot be neglected.First, we obtain N (0) e , b as a solution of the transport equation with-out considering SSC interactions. Then, we calculate the SSCcooling rate with the obtained N (0) e , b , and then we include it inthe transport equation to obtain a new approximation N (1) e , b . Weiterate this process until it converges to the correct N e , b .The solution of Eq. (27) is given by: N i , j ( E i ) = (cid:90) ∞ E i dE (cid:48) Q i , j ( E (cid:48) ) | b i ( E (cid:48) ) | exp (cid:34) − (cid:90) E (cid:48) E i dE (cid:48)(cid:48) T esc | b i ( E (cid:48)(cid:48) ) | (cid:35) , (28)and we show in Fig. 4 the distributions N e , b and N p , b obtainedfor the base zone using with the same parameter values as forFig. 3. It can be seen in the left panel of Fig.4 that the protondistribution is higher in the case of Bohm escape as comparedto the faster constant escape case. The dependence with the en-ergy is still the same, since the dominant cooling process in theBohm escape case is adiabatic cooling, which is constant. Inthe case of electrons (right panel of Fig. 2.1), there is no dif-ference between the Bohm and the constant escape cases sincesynchrotron emission largely dominates.Once we obtain the solution of Eq. (27), it is possible tocompute the injection of charged pions produced by p γ and pp interactions. We compute the pion injection using accurateapproximations to the SOPHIA code for p γ interactions [40]given in Ref. [41]: Q p γ → π ± ( E ) = (cid:90) ∞ E p dE p E p N p , b ( E p ) (cid:90) ∞ (cid:15) th n ph ( E ph ) R π ( E p , E ph ) , (29)where R π ( E , E ph ) is a function that depends on the cross section σ p γ and includes the di ff erent channels for pion production, asdiscussed in Ref. [41].As for the injection due to pp interactions, we compute it as Q pp → π ± ( E ) = n p c (cid:90) ∞ E N p , b ( E p ) F π ( E p , E ) σ pp ( E p ) , (30)where σ pp ( E p ) is the pp cross section and F π is a fitting func-tion given in Ref.[30] to reproduce the outputs of the SIBYLLsimulation code [42]. In order to obtain the pion distribution,we include the corresponding decay term in the trasport equa-tion: d (cid:104) b i N i , j ( E i ) (cid:105) dE i + N i , j ( E i ) T esc + N i , j ( E i ) T i , d ( E i ) = Q i , j ( E i ) , (31) where T i , d is the particle lifetime and “ i ” refers to pions ( i = π )and muons ( i = µ ). Once N π, j is obtained, we can compute theinjection of the muons generated by pion decays Q µ, j applyingthe formulae given in Ref. [43], which account for the kine-matics of the decay process. We then plug the muon injectioninto Eq. (31) and compute the corresponding muon distribution N µ, j using an analogous expression to Eq. (28). The result isshown in Fig. 5, along with the pion distribution N π, j , both cor-responding to the base zone with the same parameter values asin the previous figures. The electrons injected at the base zone su ff er severe syn-chrotron losses for the values of jet power and magnetic fieldconsidered, i.e., the radiative cooling rate dominates over theescape rate. Therefore, the power injected by the electrons atthe base zone is completely radiated there before the electronsescape. Conversely, as it can be seen in Fig. 3, the escaperate of protons dominates, along with the adiabatic cooling rate,which means that protons do not participate very e ffi ciently inthe cooling processes at the jet base. Thus, a significant fractionof them escape from the base zone to continue their propagationalong the rest of the cone-shaped inner jet. The position in thejet where the conical jet zone begins is z eoi = z acc + ∆ z b , whichis the end of the base zone. In turn, the end of the extendedconical region is determined by the reconfinement point z rec ,where, as mentioned above, the pressure exerted by the cocoonchanges the geometry of the jet into a cylinder. Supposing thatthe total rate of protons escaping from the base zone is equal tothe total rate of protons injected in the conical jet zone, N p , b ( E p ) t − Γ∆ V b = (cid:90) Q p , c ( E p ) dV c , it follows that the injection term in the second zone can be ex-pressed by: Q p , c ( z j , E p ) = ∆ V b Γ π R ( z eoi ) N p , b ( E p ) t − p , esc δ ( z j − z eoi ) . (32)In order to obtain the distribution of protons along the ex-tended conical jet region, we consider a more general transportequation with a convection term. It is convenient to expresses itusing spherical coordinates, so that the transport equation reads[44]: v j Γ r ∂ ( r N i , c ) ∂ r − ∂ (cid:0) b i N i , c (cid:1) ∂ E + N i , c T i , d = Q i , c , (33)where the convection term is the first one on the left memberand r is the radius in spherical coordinates with the origin in theBH. The term of decay is omitted for protons and it is kept in thecase of pions and muons in the form T i , d ( E i ) = T i , d (cid:16) E i m i c (cid:17) , where T i , d is the lifetime of the particle at rest. We solve Eq. (33)applying the method of the characteristic curve as described inthe appendix Appendix A.8 E p N p , b [ G e V c m - s r - ] E p [GeV] −2 −1 proton distribution [base zone] constant escapeBohm escape E e N e , b [ G e V c m - s r - ] −4 −2 E e [GeV] −2 −1 electron distribution [base zone] Figure 4: Primary particle distributions as a function of the energy at the base zone in the cases of a constant escape rate (solid lines) and Bohm escape rate (dashedlines). E π N π , b [ G e V c m - s r - ] −1 E π [GeV] pion distributions [base zone] pions from pγ, Bohm esc.pions from pγ, const. esc.pions from pp, Bohm esc.pions from pp, const. esc. E μ N μ , b [ G e V c m - s r - ] −1 E μ [GeV] muon distributions [base zone] muons from pγ, Bohm esc.muons from pγ, const. esc.muons from pp, Bohm esc.muons from pp, const. esc. Figure 5: Distributions of secondary pions and muons as a function of energy at the base zone in the cases of a constant escape rate (solid lines) and Bohm escaperate (dashed lines).
9n the case of protons, T p , d → ∞ , and after integrating andsimplifying Eq.(A.8), we obtain N p , c ( r , E p ) = ∆ V b N p , b ( E (cid:48) ( r eoi )) t − π R v j H ( r eoi − r min ) × (cid:18) r eoi r (cid:19) + C a (1 + C a ) r + C a ( r + C a − r + C a eoi ) C b E p − (1 + C a ) r r + C a eoi (34)This result is shown in Fig. 6 as a function of the position andthe energy in the case of a MQ of age t MQ = . × yr atredshift z =
8. As it can be seen, the distribution in the case ofa Bohm di ff usion escape is lower and with a flatter dependenceon the energy. This is because the Bohm escape is proportionalto the energy and hence this dependence a ff ects the injection atthe conical region. -20 -15 -10 -5 proton distribution [conical jet] l o g [ z j / c m ] E p [ G e V ] E p N p , c [ G e V c m - s r - ] const. escapeBohm escape Figure 6: Proton distributions at the conical part of the inner jet for a con-stant escape rate from the base zone (magenta) and for a Bohm escaperate (green) . The results correspond to a MQ at redshift z = t MQ ∼ . × yr. The dominant neutrino production process for such protons isthrough pp interactions. This is because the electrons injectedin the base zone radiate practically all their power there beforethey can be injected in the conical zone. In other words, theescape rate for electrons at the base zone is orders of magnitudebelow the synchrotron cooling rate (see Fig. 3, top-right panel),meaning that the electrons that can escape to be injected in theconical region carry only a negligible power. Therefore, thereis no significant electron synchrotron radiation to act as targetfor p γ interactions in the conical region.The corresponding distributions of the produced pions and ofthe muons generated by pion decays are found using Eq. (A.8)along with the injection given by Q i , c ( r , E ) = Q i ( z j , E ) + ∆ V b Γ π R ( z eoi ) N i , b ( E ) t − p , esc δ ( z j − z eoi ) , (35) where “ i ” refers to pions ( i = π ) and muons ( i = µ ). The firstterm accounts for the injection produced along the jet and thesecond term corresponds to the injection due to the escape fromthe base zone. The obtained distributions of pions and muonsare to be used in the calculation of the neutrino emission, asdiscussed below. Since we are interested in capturing all the relevant neutrinoproducing processes that can be triggered by high energy pro-tons accelerated in Pop III MQs, and, in particular, taking intoaccount that escape is dominant in the outermost zone of thesystem, we simply consider the injection of the protons escap-ing from the shell into an external zone at the IGM. In this way,we account for the possibility that protons that are accelerated atthe bow shock and escape from the shell could, in turn, generateneutrinos by p γ interactions on the CMB if they are energeticenough to produce pions.The corresponding energy density of CMB photons is [45]: n ph , CMB ( z , E ph ) = π E ( hc ) (cid:104) exp (cid:16) E ph k B T (1 + z ) (cid:17) − (cid:105) , (36)where T = .
725 K. Taking the size of this external zone to be ∆ z ext =
10 Mpc is adequate to consider the CMB as constantand homogeneous within. And this also leads to an escape ratewhich is lower than the p γ cooling rate for energies above ∼ × GeV, where pion production is activated. This can beseen in the left panel of Fig. 7 for MQs of age t MQ = × yrat a redshift z = t − = t − p γ, CMB , t − = c / ∆ z ext , and: Q p , ext ( E ) = ∆ V bs ∆ V ext N p , bs ( E ) t − , bs , (37)where ∆ V ext is the volume of the external zone. Eq. (37) en-sures that the total power injected matches the total power thatescapes from the shell carried by protons. The resulting dis-tributions of protons are shown in the right panel of Fig.7 forvarious redshifts, t MQ (cid:39) . × yr. We also show the casesfor a constant escape from the shell and for a escape term as-suming Bohm di ff usion. In the latter case, it can be seen thatthe proton distributions are below the ones corresponding to aconstant escape except at the highest energies.The distributions of secondary pions and muons are obtainedfollowing the procedure described above and are used to com-pute the expected neutrino output, as is discussed in the nextsection.
4. Neutrino and electromagnetic emission
In each zone considered in the model, pp and p γ interactionslead to the production of charged pions, which decay to neu-trinos and muons, and the latter also decay yielding neutrinos.The accompanying broadband photon emission co-produced by10 t - [ s - ] −17 −16 −15 −14 −13 −12 E p [GeV] p γ and escape rates [external zone] pγ, z=10pγ, z=8pγ, z=6pγ, z=3escape E p N p , e x t [ G e V c m - s r - ] −27 −26 −25 −24 −23 −22 E p [GeV] proton distributions [external zone] z=3, Bohm esc.z=3, const. esc.z=5, Bohm esc.z=5, const. esc.z=8, Bohm esc.z=8, const. esc.z=10, Bohm esc.z=10, const. esc. Figure 7: Proton cooling rates at the external zone for Pop III MQs of age t MQ (cid:39) . × yr at di ff erent redshifts (left panel), and the corresponding protondistributions (right panel) for a constant escape rate from the shell (solid lines) and for a Bohm escape rate (dashed lines). the high energy particles in our model is also computed consis-tently to check that it is not in conflict with any existing bound,as we show below.The emissivity ν µ + ¯ ν µ from direct pion decays can be ob-tained following Ref. [43]: Q π → ν µ ( E ν ) = (cid:90) ∞ E dE π T − π, d ( E π ) N π ( E π ) Θ (1 − r π − x ) E π (1 − r π ) , (38)where Θ ( x ) is the step function, x = E ν / E π , and the pion life-time is T π, d = . E π m π c × − s. The contribution from muondecays to ν µ + ¯ ν µ is [43] Q µ → ν µ ( E ν ) = (cid:88) i = (cid:90) ∞ E dE µ E µ T − µ, d ( E µ ) N µ i ( E µ ) × (cid:34) − x + x + (cid:32) x − − x (cid:33) h i (cid:35) , (39)where x = E ν / E µ , the muon lifetime is T µ, d = . E µ m µ c × − s, µ , = µ − , + L and µ , = µ − , + R . Here, L and R indicate the helicityof the muons, that is h i = h i = − ν e + ¯ ν e , the emissivity from thedecay of muons is given by [43]: Q µ → ν e ( E ν ) = (cid:88) i = (cid:90) ∞ E ν dE µ E µ T − µ, d ( E µ ) N µ i ( E µ , t ) × (cid:104) − x + x + (cid:16) − x + x − x (cid:17) h i (cid:105) . (40)The higher the Lorentz factors of the plasma in the di ff erentemission zones, the more boosted in the direction of the jet theobserved flux would be. Certainly, counter-jets are de-boosted.Nevertheless, we still account for their contributions, since theycan be significant, particularly in the cases of the shell and ex-ternal zone, where bulk velocities are lower. The comoving emissivities above can then be transformed tothe local frame at rest with the central BH to give Q (cid:48) ν ( E (cid:48) ν ) = D i j Q ν (cid:32) E (cid:48) ν D i j (cid:33) + D π − i j Q ν (cid:32) E (cid:48) ν D π − i j (cid:33) , (41)where E (cid:48) ν = D i j E com ν is the energy in the BH frame and theDoppler factor corresponding to a viewing angle i j is D i j = (cid:104) Γ j (1 − β j cos i j ) (cid:105) − , (42)and β j is the bulk velocity of the zone j in units of c . We remarkthat the first term on the right member of Eq. (41) includes thecontribution of the jet, whereas the second term accounts forthe counter-jet contribution.The neutrino spectrum, corresponding to one MQ at redshift z for which the angle of the jet with the line of sight is i j can becomputed as: dN (cid:48) ν dE (cid:48) ν d Ω (cid:48) = Q (cid:48) ν ( E (cid:48) ν ) dV dt MQ (43)In order to compute the di ff use neutrino flux due to all possi-ble Pop III MQs that existed along the history of the universe,we consider the rate of their formation per unit mass to be afraction of the corresponding rate of Pop III star formation: dR MQ ( z ) dM = f BH f bin dR PopIII ( z ) dM , (44)which represents the number of MQs generated per unit time,per unit volume and per unit mass of the Pop III stars produced( M ). We consider that a fraction f BH (cid:39) . M min , MQ = M (cid:12) produced BHs of about half ofits mass, according to Ref.[46], and a fraction f bin (cid:39) . dR PopIII ( z ) dM ∝ M − b , with b = (0 − b =
0) [48], which leads toa greater number of high mass systems and hence to a higherneutrino emissivity overall.We proceed to normalize, at each redshift z , the distribution dR PopIII ( z ) dz using the total mass generated in Pop III stars ˙ M PopIII ( z )according to Ref. [49], (cid:90) M max M min M dR
PopIII ( z ) dM dM = ˙ M PopIII ( z ) , (45)where we suppose that the possible range of masses for the starsis M min (cid:39) . M (cid:12) and M max = M (cid:12) . The total generation rateof Pop III MQs created by the evolution of stars with massesabove M min , MQ is a fraction of the total generation rate of PopIII stars ˙ M PopIII ( z ). The latter is shown for illustration in Fig.8and we also show it weighted by H | dtdz | , where (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) dtdz (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = H (1 + z ) (cid:112) (1 + z ) Ω m + Ω Λ , (46)with H =
70 km s − Mpc − , Ω m = . Ω Λ = . [ M ☉ M p c - y r - ] −7 −6 −5 −4 −3 −2 z dtdzH0-RMQ-s-1cm-3RMQ-cm-3-s-1H0dtdzRMQinMsun2yr2Mpc3RinMsun2yr2Mpc3 M PopIII H |dt/dz| M PopIII
Pop III star formation rate vs redshift . .
Figure 8: Formation rate of Pop III stars, ˙ M PopIII ( z ), in blue, adopted fromRef.[49]. The product H | dtdz | ˙ M PopIII ( z ) is shown in red, as this is useful for thecalculation of the di ff use neutrino flux. The di ff erential density of the produced neutrinos using asimilar expression to the given by Refs. [50, 51], i.e., dn ν ( E ν ) dM = (cid:90) d Ω (cid:48) dR MQ ( z ) dM (1 + z ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) dtdz (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) dz dN (cid:48) ν dE (cid:48) ν d Ω (cid:48) dE (cid:48) ν (1 + z ) − . (47)Eq. (47) accounts for the contributions of the jet and counter-jet, assuming that the orientation is distributed isotropically,that is, dR MQ dM is independent of i j . According to Fig. 8, it canbe concluded that the dominant neutrino contribution arises forreshifts z ≈ −
8, where H | dtdz | ˙ M PopIII ( z ) peaks. Integration over the total MQ life T MQ , the solid angle, andvolume of the emitting zone j , yields, in units of [energy − ], thespectrum of the muonic flavor of neutrinos and antineutrinos, ν µ + ¯ ν µ : dN (cid:48) ν dE (cid:48) ν = π (cid:90) ∆ V j dV j (cid:90) T MQ dt MQ (cid:90) π di j sin( i j ) × (cid:104) Q (cid:48) ν µ ( E (cid:48) ν ) P ν µ → ν µ + Q (cid:48) ν e ( E (cid:48) ν ) P ν e → ν µ (cid:105) . (48)Here, P ν µ → ν µ (cid:39) .
453 is the probability that the generated ν µ or¯ ν µ keep the same flavor, and P ν e → ν µ (cid:39) .
171 is the probabilitythat ν e or ¯ ν e oscillate into ν µ or ¯ ν µ . These probabities are derivedfrom the unitary mixing matrix U α j , which is determined bythree mixing angles, θ (cid:39) . ◦ , θ (cid:39) . ◦ , and θ (cid:39) ◦ ,and the CP-violating phase δ CP ≈ ◦ [52]. The values usedfor the probabilities correspond to a normal mass ordering ofthe massive neutrinos ( ν , ν , ν ), i.e., m < m < m .Considering that the emission from any redshift z is the samein all directions, the di ff erential neutrino density can be relatedto the di ff erential neutrino flux with as d Φ ν dE ν = c π n ν . Hence, thefinal expression for the di ff use neutrino flux originated in PopIII MQs is given by d Φ ν ( E ν ) dE ν = c π (cid:90) z max z min dz (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) dtdz (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:90) M max , MQ M min , MQ dM dR MQ ( z ) dM dN (cid:48) ν dE (cid:48) ν , (49)where M min , MQ = M (cid:12) and M max , MQ = M (cid:12) limit the rangeof masses of Pop III stars supposed to lead to the formation aMQs, and z min = z max =
25 indicate the limiting values ofredshift along which Pop III MQs were distributed according to˙ M PopIII ( z ).Taking into account the dominating cooling processes de-scribed above, we can make some simple order-of-magnitudeestimates of the neutrino flux that could be expected at theEarth. As it can be seen from Fig. 3, p γ interactions at the basezone are most e ff ective at high energies, E p ∼ GeV, andescape dominates otherwise for the adopted parameters in thecase of a constant escape rate. The escaping protons cool dom-inantly by adiabatic expansion, but they can still undergo pp interactions along the rest of the conical part of the jet. Consid-ering that, on average, ∼
20% of the energy of the parent protongoes to the produced neutrinos [34], we can estimate that the av-erage power carried by the final neutrinos generated in the innerjet to be L ν, jet ≈ . L p , b for a neutrino energy range between E (cid:48) ν, ≈ . × GeV and E (cid:48) ν, ≈ . × GeV. This estimationaccounts for the fact that p γ interactions do not dominate overthe whole mentioned energy range, but only for the most ener-getic protons. Considering a typical lifetime T MQ ≈ × yrfor Pop III MQs, a simplistic ∼ E (cid:48)− ν single-source spectrum ofneutrinos of all flavors can be obtained as dN (cid:48) ν dE (cid:48) ν (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) jet ≈ L ν, jet T MQ log E (cid:48) ν, E (cid:48) ν, E (cid:48)− ν (50) (cid:39) . × GeV − (cid:32) L p , b L i (cid:33) (cid:32) T MQ T . (cid:33) (cid:32) E (cid:48) ν GeV (cid:33) − , where L i = × erg s − and T . = × yr. Similarly,since the protons accelerated at the shell escape to the external12 ν d Φ ν / d E ν [ G e V c m - s r - s - ] − − − − − − E ν [GeV] totalbaseconical jetshellexternal zone [muon neutrinos], T esc = const., Γ j = acc = 2× 10 R g , Δz j = 5R j E ν d Φ ν / d E ν [ G e V c m - s r - s - ] − − − − − − E ν [GeV] totalbaseconical jetshellexternal zone [muon neutrinos], T = T B , Γ j =1.67, z(E) acc = 2× 10 R g , Δz j = 5R jesc IceCube IceCube Auger G R A N D G R A N D IceCubeIceCube Auger
Figure 9: Di ff use neutrino flux including the major contributions in the cases of a constant escape rate (left panel) and a Bohm escape rate (right panel). The shadedregion corresponds to varying the index b characterizing the mass distribution of MQs between b = b = zone, and those with energies E p (cid:38) × GeV photo-producepions e ffi ciently by interactions with the CMB, the power car-ried by the neutrinos generated is roughly L ν, ext ≈ . × L p , bs ( E p > × GeV) (cid:39) . × × erg s − (cid:32) L p , bs . L i (cid:33) . Therefore, an estimate for a typical spectrum of the neutrinoproduced in the external zone is dN (cid:48) ν dE (cid:48) ν (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ext ≈ L ν, ext T MQ log E (cid:48) ν, E (cid:48) ν, E (cid:48)− ν (51) (cid:39) . × GeV − (cid:32) L p , ext . L i (cid:33) (cid:32) T MQ T . (cid:33) (cid:32) E (cid:48) ν GeV (cid:33) − , with E ν, = . × GeV and E ν, = . × GeV.The integrals on M and z of Eq.(49) can be estimated makingthe rough approximation that the rate of generated MQs is suchthat H (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) dtdz (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) dRdM ∆ M ∼ f BH f bin ˙ M PopIII M (cid:12) (cid:39) . × − Mpc − yr − for redshifts between z = z =
10 (see Fig. 8). Assum-ing that flavor mixing leads to an approximate equal ratio forthe three neutrino flavors, the contributions to the di ff use fluxfrom the inner jet and the external zone are roughly given by E ν d Φ ν dE ν (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) jet ≈ × − GeV cm − sr − s − (cid:32) L p , b L i (cid:33) (cid:32) T MQ T . (cid:33) for 250GeV (cid:46) E ν (cid:46) . × GeV E ν d Φ ν dE ν (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ext ≈ − GeV cm − sr − s − (cid:32) L p , bs . L i (cid:33) (cid:32) T MQ T . (cid:33) for 4 × GeV (cid:46) E ν (cid:46) × GeV . While useful as order of magnitude estimations, these ex-pressions clearly do not account for the exact dependence ofthe particle distributions, injections, and intervening coolingrates, so that, for instance, the e ff ect of synchrotron losses bypions and muons at the inner jet were not included at that point.Another important issue is accounting for MQs with di ff erentBH masses which arise, as explained above, by the gravita-tional collapse of a Pop III stars with masses between 50 M (cid:12) and100 M (cid:12) in binary systems. We address this point by consideringthat the power injected in relativistic particles is proportionalto the BH mass, and hence to M . Therefore neutrino emissionis also proportional to M for all the processes except for p γ atthe jet base, since there the target photons correspond to syn-chrotron emision by the electrons, which is also proportionalto M . Hence, neutrino production the jet base is considered toscale as ∝ M . We apply these scalings to perform the integra-tion over M in Eq.(49) making use of our central result obtainedfor M BH = M (cid:12) , i.e. for a Pop III star mass M (cid:39) M (cid:12) .This can be performed for di ff erent cases of mass distributions dR MQ dM ∝ M − b , with b = (0 −
2) as discussed above. In Fig. 9we show the results obtained with the full numeric code for thedi ff use neutrino flux of ν µ + ¯ ν µ in the case of b =
1, and the grayshaded region indicates the possible range of the flux values be-tween the lowest flux corresponding to b = b =
0. We also show individually the most significantcontributions among the di ff erent emission zones consideredfor the cases of constant escape rates (left panel) and for Bohmescape rates (right panel). We also include the fit obtained forIceCube data [19, 20], as well as the upper limits of higher en-ergy neutrinos given by Auger [21] and IceCube [53]. For ref-erence, we also show the expected sensitivity for GRAND [22],but other planed detectors will be sensible to UHE neutrinosas well, such as IceCube-Gen2 [54], PUEO [55], RNO-G [56],Trinity [57], and BEACON [58].13 ν d Φ ν / d E ν [ G e V c m - s r - s - ] − − − − − E ν [GeV] Δ z j = 1 R j Δ z j = 5 R j Δ z j = 10 R j Δ z j = 20 R j [muon neutrinos] - α =2, Γ j =1.67, z acc =2×10 R g , Δ z j =(1– 20)R j E ν d Φ ν / d E ν [ G e V c m - s r - s - ] − − − − − E ν [GeV] z acc = 10 R g z acc = 2×10 R g z acc =4×10 R g [muon neutrinos] - α =2, Γ j =1.67, z acc =(1– 4)×10 R g , Δ z j =5R j E ν d Φ ν / d E ν [ G e V c m - s r - s - ] − − − − − E ν [GeV] Γ j = 1.25 Γ j = 1.67 Γ j = 3 Γ j = 5 Γ j = 10 [muon neutrinos] - α =2, Γ j =1.25–10, z acc =2×10 R g , Δ z j =5R j E ν d Φ ν / d E ν [ G e V c m - s r - s - ] − − − − − E ν [GeV] α =1.8 α =2 α =2.2 [muon neutrinos] - α =1.8–2.2, Γ j =1.67, z acc =2×10 R g , Δ z j =5R j IceCube IceCubeIceCube IceCubeIceCubeIceCubeIceCube IceCube Auger AugerAugerAuger G R A N D G R A N D G R A N D G R A N D Figure 10: Di ff use neutrino flux resulting with di ff erent combinations of parameters in the case of a constant escape rate.
14n Fig. 10, we plot the results corresponding to the di ff useneutrino flux if four key parameters of the model are varied,adopting for illustration a constant escape rate. In the top leftpanel, we show the fluxes obtained for di ff erent values of theposition of the emitter in the jet base ( z acc ), while in the topright panel, the size of the base zone ( ∆ z b ) is varied. Likewise,the resulting flux is shown for di ff erent values of the jet Lorentzfactor ( Γ ) in the bottom left panel, and with di ff erent values ofthe index of injection of primary particles ( α ). Although forsimplicity we have kept the ratio of magnetic to kinetic energyat the base as constant ( q m = × − ), varying z acc leads to dif-ferent values of the magnetic field at the base zone, since its sizeis set in reference to the expanding jet radius. Then, in the topleft panel of Fig.10, ∆ z b = R j (cid:39) . z acc is assumed in the threecases, and we obtain B acc (cid:39) . z acc R g , with R g (cid:39) . × cm.In the top right panel, we fix z acc = × R g and changingthe size of the base zone basically modifies linearly the escaperate and the p γ cooling rate. Therefore, for instance, for thehighest value considered ( ∆ z b = R g ), both rates are low lead-ing to a less e ff ective neutrino production in comparison withthe other cases for smaller sizes. In particular, for the smallestvalue adopted ( ∆ z b = R g ), it can also be seen that the maximumneutrino energy is lower, and this is because the accelerationrate is the same for all the cases of that panel and the maximumproton energy is correspondingly lower for the high rates of es-cape and p γ collisions. We note that, given the values of themagnetic field considered at the base in general, the electroncooling is so fast that no significant synchrotron emission takesplace outside the injection region at the jet base. Therefore, ifthe volume of this zone is increased, the density of synchrotrondecreases and p γ become less e ff ective.In the bottom left panel of Fig.10, it can be seen that the con-tribution from the inner jets decreases as the bulk Lorentz factorof the jet increases. This can be understood as a consequenceof the fact that under the assumptions made, the neutrino emis-sivity in the comoving frame is Q ν ∝ Γ − , as is shown in Ap-pendix Appendix B. Therefore, when transformed to the BHframe, Eq.(41) implies that Q (cid:48) ν ∝ Γ − and this is reflected inthe final possible fluxes to arrive at the Earth. Varying Γ alsoimplies varying the magnetic field at the base, since the mag-netic energy is proportional to the kinetic one, and the latter is ∝ Γ ( Γ −
1) (Eq.2). Therefore, for
Γ = . B acc ≈ . × G,and for
Γ =
10 we have B acc ≈ . × G. In the bottom rightpanel, we show the di ff use neutrino flux obtained for other val-ues of the spectral index of primary particle injection: α = . α = .
2. As it can be seen, the prospects for detection fallfor steeper injections.The flavor ratios of neutrinos has become an interesting ob-servable which can bring information on the nature of the pro-duction mechanism operating at the sources [62, 59, 61, 60, 63].In Fig. 11, we show the neutrino flavor ratios that are obtainedwithin our model for Pop III MQs, in the case of a Bohm escaperate, but the result is very similar for a the constant escape case.The e ff ect caused by the magnetic field is manifest for the en-ergy window ∼ (3 × − )GeV, where high energy muonsat the inner jet are a ff ected by synchrotron losses and a deficit of electron neutrinos are produced. For still higher energies, theflux is dominated by the contribution from the escaping protonsinteracting with the CMB, and no magnetic field e ff ects are ex-pected. This also happens for neutrino energies E ν (cid:46) GeV,for which synchrotron losses of pions and muons at the innerjet are not significant. f ν E ν [GeV] f ν f ν f ν e μτ neutrino flavor ratios, α =2, Γ =1.67, z acc =2×10 R g , Δ z b =5R j Figure 11: Neutrino flavors ratios as a function of energy in the case of a Bohmescape rate.
For completeness, we compute the di ff use background fluxof multiwavelength photons that are co-produced along withthe neutrinos, using a expression analogous to Eq. (49). Themain contributing processes are synchrotron emission, IC inter-actions, pp , and p γ collisions, and we compute the correspond-ing emissivities following, e.g., Refs. [26, 10, 64]. We show inFig. 12 the photon spectra obtained for MQs at redshift z = ff use flux obtained isshown in Fig. 13, where it can be seen that the flux level is wellbelow that of the extragalactic background of multiwavelengthphotons [65]. We include the flux corrected by γγ absorption onthe CMB and EBL through an exponential factor e − τ γγ , wherethe optical depth τ γγ is integrated following Ref. [66].
5. Discussion
In this work, we have applied a model that allows to obtain adi ff use neutrino flux produced by a distribution of Pop III MQsduring their lifetime at a wide range of redshifts ( z = − E γ d N γ / d E γ [ e r g ] log(E γ /eV) −4 −2 0 2 4 6 8 10 12 14 16 18 20 [photons] base zone - T esc =const.e-syn SSC p - s y n μ - s y n π - s y n pp p γ E γ d N γ / d E γ [ e r g ] log(E γ /eV) −4 −2 0 2 4 6 8 10 12 14 16 18 20 [photons] conical, shell, and external zones - T esc =const.e-syn (sh) I C ( s h ) p - s y n ( c o n ) μ - s y n ( c o n ) π - s y n ( c o n ) pp(con) p γ (ext)p γ (sh)pp(sh)pp(ext) E γ d N γ / d E γ [ e r g ] log(E γ /eV) −4 −2 0 2 4 6 8 10 12 14 16 18 20 [photons] base zone - T esc =T B (E)e-syn SSC p - s y n μ - s y n π - s y n pp p γ E γ d N γ / d E γ [ e r g ] log(E γ /eV) −4 −2 0 2 4 6 8 10 12 14 16 18 20 [photons] conical, shell, and external zones - T esc =T B (E)e-syn (sh) IC (sh) p - s y n ( c o n ) μ - s y n ( c o n ) π - s y n ( c o n ) pp(con)p γ (ext)p γ (sh) pp ( s h ) pp ( e x t ) Figure 12: Total photon spectra emitted by MQs at redshift z = γ d Φ γ / d E γ [ e r g c m - s r - s - ] − − − − − log(E γ /eV) − − − − γ rays diffuse flux [multiwavelenth photons] Figure 13: Di ff use flux of multiwavelength photons from Pop III MQs in thecase of a constant escape rate as compared to existing data of the extragalacticphoton background, adapted from Ref. [65] . In particular, we assumed a high acceleration e ffi ciency andalso that the injected power in electrons ( L e , j ) is the same asthat in protons ( L p , j ) at each emission zone “ j ” considered. Weexplored di ff erent combinations of parameters, varying the po-sition and size of the emitter at base zone of the inner jet, theLorentz factor of the jet, and the index of injected particles. Outof the possibilities mentioned, which involve plausible valuesfor the parameters, we select the case shown in right panel ofFig. 9 as a representative one for e ffi cient neutrino production.The main contributions to the di ff use neutrino flux arise at theinner jets for energies ∼ [10 − ]GeV, and at the externalzone for higher energies ( ∼ [10 − ]GeV).The inner jets are relevant sites for neutrino production forthe following reasons: first, at the base of the jet the magneticfield derived is strong ( B ∼ G), which enhances the acceler-ation e ffi ciency and also the production of the low energy pho-tons generated by electron synchrotron, thus favoring p γ inter-actions. And second, the target of cold protons at this regionis highly dense, n p ∼ cm − at z j = z acc , which favors the pp interactions at the jet base and at the conical part of the jet.The contribution from the inner jet still do not reach the level ofthe detected neutrino flux according to a global fit of IceCube[19], and the best neutrino fit of astrophysical ν µ + ¯ ν µ [67], atenergies ∼ × GeV. In order to have a higher neutrinoflux, the typical jet power could be increased, and / or the typicalMQ lifetime, but this would be inconsistent with the results ob-tained by simulations in Ref. [24]. The mentioned parameterswere taken from this reference an the fraction of power that isinjected in relativistic particles was considered with the typicalvalue q rel = . ∼ GeV (see Fig.3), thus generating neutrinos peaking inthe energy range ∼ (10 − )GeV. However, since protonescape from the shell actually dominates over p γ interactionswithin it, the great majority of the protons indeed escape andare injected into the IGM. The external zone considered allowsto account for the possibility that further p γ interactions withthe CMB take place outside the MQs, and we found that thesegive the major contribution at the highest energy part of theobtained neutrino output. This contribution does not violatethe upper limits given by Pierre Auger Observatory and Ice-Cube, but could still be at the reach of future detectors suchas GRAND, as shown in Figs. 9 and 10. Furthermore, forenergies from ∼ to 10 GeV, this contribution to the dif-fuse neutrino flux overlaps the energy range that the flux ex-pected from cosmogenic neutrinos produced by the interactionof ultra-high energy cosmic rays (UHECRs) with photon tar-gets from the CMB. Cosmogenic neutrinos are sensitive to thechemical composition of UHECRs, namely, their expected fluxis higher for higher proton content in UHECR with respect toheavier nuclei. On the other hand, since no significant heavynuclei contribution is expected from Pop III MQs because theseelements are released by supernova explosions and Pop III starsare the first generation of stars in the universe and have zerometallicity. This means that these stars basically burn hydrogento helium so that heavier nuclei are not present in the accreet-ing matter, and hence can not be accelerated in jets of Pop IIIMQs. The simple approach applied to obtain this contribution isstill adequate as long as over the interaction length, the photonbackground can be considered as constant, and this conditionis satisfied. Since, as mentioned no significant contribution ofheavy nuclei is present, hence it is in principle not necessaryto account for a cascade of nuclear reactions. We also do notcompute any electromagnetic cascade that would develop byinteractions with the CMB. However, the emission by the dom-inant processes allows to conclude that there is no conflict withdata (see Fig. 12), and this is enough for the purposes of thepresent work.Detailed studies for cosmogenic neutrino production accountfor in-source nuclear cascades [68, 69] to characterize the cor-rect level of neutrino flux consistent with di ff erent chemicalcompositions. Therefore, if cosmic ray data finally establisheda chemical composition consistent with a very weak flux of ac-companying cosmogenic neutrinos, and if future neutrino ob-servations yield a di ff use signal above the predicted level, thenthe posibility that the sources of these neutrinos could be PopIII MQs should not be ruled out . Conversely, in case of a futurenon-detection of the high energy part of the neutrino flux pre-dicted, this would require that either Pop III MQs themselvesdid not generate at the rate here assumed, and / or that their e ffi -ciency for accelerating protons at their shells should be boundto a lower value than the assumed in this work. For instance,lower values of the e ffi ciency of acceleration η would shift thebumps of the main contributions towards lower energies, andin the case of the external zone neutrino production by protonsaccelerated at the shell could even be supressed if the pion pro-duction threshold is not reached.Future neutrino observations with new generation instru-17ents such as IceCube-gen2, GRAND, PUEO, RNO-G, Trin-ity, and BEAC will be useful to probe the flux neutrinos fromPop III MQs at the highest energies. This will also help to ob-tain more accurate measurments of the flavor composition alongan extense energy range, which would constrain neutrino pro-ducing models such as the presented in this work and yield morelight on the origin of astrophysical neutrinos. Appendix A. Solution of the inhomogeneous transportequation with convection and decay
Here we describe the steps followed to solve Eq.(33) usingthe method of the characteristics. We assume the boundary con-dition N c ( r , E i ) | r −→ =
0, i.e., the escaping particles have a van-ishing distribution at r (cid:28) z eoi . Since the dominant cooling pro-cesses are adiabatic expansion and synchrotron emission, werewrite the transport equation as:1 r ∂ ( r N i , c ) ∂ r − (cid:34) C a Er + C b E r (cid:35) ∂ N i , c ∂ E − (cid:34) C a r + C b Er − C c E (cid:35) N i , c = Q i , c Γ v j , (A.1)with the constants are given by: C a = Γ (A.2) C b = (cid:32) m e m i (cid:33) σ T c B π z m e c Γ v j m p c (A.3) C c = m i c T Γ v j (A.4)The solution to the characteristic equation dE i dr = − C a E i r − C b E i r , (A.5)gives the characteristic curve E (cid:48) ( r (cid:48) ; r , E i ) = (1 + C a ) E i r (cid:48) r + C a (1 + C a ) r (cid:48) + C a r + C b E i ( r (cid:48) + C a − r + C a ) . (A.6)Using the curve corresponding to each pair of values ( r , E i ), wesolve the following ordinary di ff erential equation, dN i , c dr (cid:48) = Q Γ v j + (cid:32) C a r (cid:48) + C b E (cid:48) ( r (cid:48) ) r (cid:48) − C c E (cid:48) ( r (cid:48) ) − r (cid:48) (cid:33) N i , c . (A.7)to obtain the particle distribution along the inner jet as: N i , c ( r , E i ) = (cid:90) rr ini dr (cid:48) Q ( r (cid:48) , E (cid:48) ( r (cid:48) )) Γ v j × exp (cid:40)(cid:90) rr (cid:48) dr (cid:48)(cid:48) E (cid:48) ( r (cid:48)(cid:48) ) r (cid:48)(cid:48) (cid:2) C a r (cid:48)(cid:48) E (cid:48) ( r (cid:48)(cid:48) ) + C b (cid:0) E (cid:48) ( r (cid:48)(cid:48) ) (cid:1) − E (cid:48) ( r (cid:48)(cid:48) ) r (cid:48)(cid:48) − r (cid:48)(cid:48) / T i , d (cid:105)(cid:111) . (A.8)Here, r ini = max ( r acc , r min ), where r min is the value for whichthe characteristic curve goes to infinity: r min ( r , E i ) = r (cid:34) + r (1 + C a ) C b E i (cid:35) − + Ca . (A.9) Appendix B. Analytical estimate of the p γ cooling ratewith photons from electron synchrotron astargets In this appendix we estimate the cooling rate t − p γ in the casethat the target photon density the synchrotron emission of elec-trons given by Eq.(11). If synchrotron cooling dominates forelectrons, as is the case for the magnetic field values adopted,we can approximate the photon density by supposing that thesame power injected in electrons is radiated. Since the cor-responding electron distribution is N e ∝ E − e for a simplifiedinjection of electrons Q e ∼ K e E − e , with K e = L e π Γ∆ V b log E e , max E e , min and E e , max = m e c (cid:115) π e ησ T B acc , considering that the synchrotron emission is concentrated in theenergy range given by E min(max)ph = √ heB acc π m e c (cid:32) E min(max) m e c (cid:33) leads to an emissivity Q e , syn ≈ K e E − . The density of suchphotons is, then: n ph = π Q e , syn R j c ≈ K e π R j cE , (B.1)in units of [energy − length − ]. In order to estimate the t − p γ , weapply the approximation for the cross section given by Atoyan& Dermer (2003) [34], i.e., σ p γ ( E r ) = E r < . µ barn for 0 . < E r < . µ barn for E r ≥ . , (B.2)where the low energy range corresponds to the single pion ( p + γ → π + n ) with an inelasticity K = .
2, and for higher energiesthe multipion channel dominates ( p + γ → p + π + + π − + π )with K = .
6. In the case of the single-pion channel, t − p γ, ( γ p ) = K (cid:90) ∞ e th2 γ p dE ph cn ph ( E ph )2 γ p E (cid:90) E ph γ p E th dE r σ E r H ( E − E r ) t − p γ, ( γ p ) ≈ K e K σ R j π γ p (cid:90) E γ pE γ p dE ph E − (cid:16) E γ p − E (cid:17) + K e R j π ( E − E )2 γ p (cid:90) E ph , max E γ p dE ph E − t − p γ, ( γ p ) ≈ K e K σ R j π ( E − E ) γ p E E = . × − s − (cid:32) E p GeV (cid:33) (cid:32) L e L i (cid:33) (cid:32) Γ . (cid:33) − . (B.3)18or the multipion channel, we find, t − γ, ( γ p ) = (cid:90) ∞ E γ p dE ph cn ph ( E ph )2 γ p E (cid:90) E ph γ p E dE r σ E r = (cid:90) E ph , max E γ p dE ph cn ph ( E ph )2 γ p E σ (cid:16) E γ p − E (cid:17) ≈ K e σ R j πγ p E ≈ × − s − (cid:32) E p GeV (cid:33) (cid:32) L e × erg s − (cid:33) (cid:32) Γ . (cid:33) − . Hence, the total interaction rate can be approximated by t − p γ ( E p ) ≈ × − s − (cid:32) E p GeV (cid:33) (cid:32) L e L i (cid:33) (cid:32) Γ . (cid:33) − , (B.4)and this matches the result shown in Fig.3 for E p (cid:46) , whilethe exact result flattens at higher energies because the targetphoton density used in that case is corrected by synchrotronself-absorption. Appendix B.1. Estimation of neutrino emissivity
In the case of a constant escape rate t − = c / ( Γ∆ z ), the dis-tribution of protons at the base is roughly given by N p ( E p ) ≈ L p ∆ z b π c ∆ V b log (cid:18) E p , max m p c (cid:19) E − p , with E p , max ≈ × GeV. The emissivity of pions producedby p γ interactions can be approximated using the collision fre-quency ω p γ = t − p γ / K p γ . In turn, following Ref. [34], the neu-trino emissivity can be obtained approximately by supposingthat the pion energy is equally distributed among the four finaldecay products (after the muon decay). For the π + decay of thesingle pion channel, considering that the pion production takesplace half of the times as compared to the neutron production,this leads to: Q (1) ν e ( E ν ) = Q (1)¯ ν µ ( E ν ) = Q (1) ν µ ( E ν ) ≈ N p (20 E ν ) ω (20 E ν ) . For the decays the π + decay in the multipion channel, we have: Q (1) ν e ( E ν ) = Q (1)¯ ν µ ( E ν ) = Q (1) ν µ ( E ν ) ≈ N p (20 E ν ) ω (20 E ν ) , and for the π − decay of the multipion channel: Q (2) ν µ ( E ν ) = Q (2)¯ ν µ ( E ν ) = Q (2)¯ ν e ( E ν ) ≈ N p (20 E ν ) ω (20 E ν ) . The total emissivities for ν µ + ¯ ν µ and ν e + ¯ ν e are then: Q ν µ + ¯ ν µ ( E ν ) = Q ν e + ¯ ν e ( E ν ) ≈ ω (20 E ν ) N p (20 E ν ) + ω (20 E ν ) N p ( E ν ) ≈
150 GeV − cm − sr − s − (cid:18) E ν GeV (cid:19) − × (cid:32) L p L i (cid:33) (cid:32) L e L i (cid:33) (cid:32) Γ . (cid:33) − . (B.5)In this case, these simplified expressions yield results within afactor ∼ . Acknowledgments
AMC and MMR are supported by grants PIP 0046 (CON-ICET) and 15 / E870EXA912 /
18 (Universidad Nacional de Mardel Plata). GER is supported by grant PIP 0338 (CON-ICET), PICT 2017-2865 (ANPCyT), and by the Ministerio deEconom´ıa y Competitividad (MINECO) under grant AYA2016-76012-C3-1-P and PID 2019-105510GB-C31.
References [1] M. Fraser, A. R. Casey, G. Gilmore, A. Heger, and C. Chan, The massdistribution of Population III stars, MNRAS (2017) 418-425 (2017)[2] A. Stacy and V. Bromm, Constraining the statistics of Population III bina-ries, MNRAS (2013) 1094-1107[3] G. E. Romero and P. Sotomayor Checa, Population III microquasars, Int.J. Mod. Phys. D (2018) 1844019, 1-7[4] P. Sotomayor Checa and G. E. Romero, Model for population III micro-quasars, A& A (2019) A76 (2019).[5] A. U. Abeysekara, A. Albert, R. Alfaro, C. Alvarez, J. D. lvarez, R. Arceo,J. C. Arteaga-Velzquez, D. Avila Rojas, H. A. Ayala Solares, E. Belmont-Moreno, and et al., Very-high-energy particle acceleration powered by thejets of the microquasar SS 433, Nature (2018) 82-85.[6] A. Levinson and E. Waxman, Probing microquasars with TeV neutrinos,Phys. Rev. Lett. (2001) 171101.[7] C. Distefano, D. Guetta, E. Waxman, and A. Levinson, Neutrino flux pre-dictions for known galactic microquasars, ApJ (2002) 378.[8] F. Aharonian, L. Anchordoqui, D. Khangulyan, and T. Montaruli, Micro-quasar LS 5039: a TeV gamma-ray emitter and a potential TeV neutrinosource, Journal of Physics: Conference Series (2006) 408-415.[9] M. M. Reynoso, G. E. Romero, and H. R. Christiansen, Production ofgamma rays and neutrinos in the dark jets of the microquasar SS433, MN-RAS (2008) 1745-1754.[10] M. M. Reynoso and G. E. Romero, Magnetic field e ff ects on neutrinoproduction in microquasars, A& A (2009) 1.[11] J. F. Zhang, Y. G. Feng, M. C. Lei, Y. Y. Tang, and Y. P. Tian, High-energyneutrino emission from low-mass microquasars, MNRAS (2019) 2468-2474.[12] L. A. Anchordoqui, H. Goldberg, T. C. Paul, L. H. M. da Silva, and B. J.Vlcek, Estimating the contribution of galactic sources to the di ff use neutrinoflux, Phys. Rev. D (2014) 123010.[13] M. M. Reynoso and A. M. Carulli, On the possibilities of high-energyneutrino production in the jets of microquasar SS433 in light of new obser-vational data, Astropart. Phys. (2019) 25-32.[14] I.F. Mirabel, M. Dijkstra, P. Laurent, A. Loeb, and J.R Pritchard, Stellarblack holes at the dawn of the universe, A&A (2011) A149.[15] V. Berezinsky and P. Blasi, UHE neutrinos from Pop III stars: conceptand constraints, Phys. Rev. D (2012) 3003.[16] D. Xiao, P. Mszros, K. Murase and Z. Dai, Revisiting the Contributionsof Supernova and Hypernova Remnants to the Di ff use High-Energy Back-grounds: Constraints on Very-High-Redshift Injections, Astrophys.J. (2016) 133.[17] R. Schneider, D. Guetta and A. Ferrara, Gamma Ray Bursts from the FirstStars: Neutrino Signals, MNRAS (2002) 173[18] S. Gao, K. Toma and P. Meszaros, High Energy Neutrino Emission fromthe Earliest Gamma-Ray Bursts, Phys. Rev. D (2011) 103004.[19] M. G. Aartsen, K. Abraham, M. Ackermann, J. Adams, J. A. Aguilar,M. Ahlers, M. Ahrens, D. Altmann, T. Anderson, M. Archinger, and et al.,A combined maximum-likelihood analysis of the high-energy astrophysicalneutrino flux measured with IceCube, ApJ (2015) 98.[20] M. G. Aartsen, M. Ackermann, J. Adams, J. A. Aguilar, M. Ahlers, M.Ahrens, I. A. Samarai, D. Altmann, K. Andeen, T. Anderson, and et al.,Extending the search for muon neutrinos coincident with gamma-ray burstsin IceCube data, ApJ (2017) 112.[21] D. G´ora, The Pierre Auger Observatory: Review of Latest Results andPerspectives, Universe 4 (2018) 128[22] J. ´Alvarez-Mu˜niz,et al. [GRAND Collaboration], The Giant Radio Arrayfor Neutrino Detection (GRAND): Science and Design, Sci. China Phys.Mech. Astron. (2020) no.1, 219501
23] K. Ohsuga, S. Mineshige, Why is supercritical disk accretion feasible?,ApJ (2017) 1283. doi:10.1086 / (2017) 10.1103 / phys-revd.96.063014.[25] S. S. Komissarov, Magnetic acceleration of relativistic jets, Memorie dellaSocieta Astronomica Italiana (2011) 95-103.[26] G. E. Romero and G. S. Vila, The proton low-mass microquasar: high-energy emission, A&A (2008) 623-631.[27] P. Bordas, V. Bosch-Ramon, J. M. Paredes, and M. Perucho, Non-thermalemission from microquasar / ism interaction, A&A (2009) 325-334.[28] C. R. Kaiser and P. Alexander, A self-similar model for extragalactic radiosources, MNRAS (1997) 215-222.[29] M. C. Begelman, B. Rudak, and M. Sikora, Consequences of relativisticproton injection in active galactic nuclei, ApJ (1990) 38-51.[30] S. R. Kelner, F. A. Aharonian, and V. V. Bugayov, Energy spectra ofgamma-rays, electrons and neutrinos produced at proton-proton interactionsin the very high energy regime, Phys. Rev. D (2006) 034018, [Erratum:Phys. Rev.D (2009) 039901].[31] G. R. Blumenthal and R. J. Gould, Bremsstrahlung, synchrotron radiation,and compton scattering of high-energy electrons traversing dilute gases,Rev. Mod. Phys. (1970) 237-270.[32] G. E. Romero, M. Boettcher, S. Marko ff , et al., Relativistic Jets in ActiveGalactic Nuclei and Microquasars. Space Sci. Rev. (2017) 5-61.[33] G. B. Rybicki and A. P. Lightman, Radiative processes in astrophysicsWiley-VCH (2004).[34] A. M. Atoyan and C. D. Dermer, Neutral beams from blazar jets, ApJ (2003) 79-96.[35] M. S. Longair, High Energy Astrophysics, 3rd ed. Cambridge UniversityPress (2011).[36] V. Bosch-Ramon, J. M. Paredes, G. E. Romero, and D. F. Torres, A micro-quasar model applied to unidentified gamma-ray sources, A& A (2006)1081-1087.[37] R. D. Blandford and J. P. Ostriker, Particle acceleration by astrophysicalshocks, ApJL (1978) L29–L32[38] F. M. Rieger, V. Bosch-Ramon and P. Du ff y, Fermi acceleration in astro-physical jets,Astrophys. Space Sci. (2007), 119-125[39] L. Sironi, U. Keshet and M. Lemoine, Relativistic Shocks: Particle Ac-celeration and Magnetization, Space Sci. Rev. (2015) 519-544[40] A. M¨ucke, R. Engel, J. Rachen, R. Protheroe, and T. Stanev, Montecarlosimulations of photohadronic processes in astrophysics, Computer PhysicsCommunications (2000) 290-314.[41] S. H¨ummer, M. Ruger, F. Spanier, and W. Winter, Simplified models forphotohadronic interactions in cosmic accelerators, ApJ (2010) 630-652.[42] R. S. Fletcher, T. K. Gaisser, P. Lipari, and T. Stanev, Sibyll: An eventgenerator for simulation of high energy cosmic ray cascades, Phys. Rev. D (1994) 5710.[43] P. Lipari, M. Lusignoli, and D. Meloni, Flavor Composition and En-ergy Spectrum of Astrophysical Neutrinos, Phys. Rev. D (2007) 123005(2007).[44] A. A. Zdziarski, P. Pjanka, M. Sikora, and M. Stawarz, Jet models forblack-hole binaries in the hard spectral state, MNRAS (2014) 2238-2254.[45] R. Ru ffi ni, G. V. Vereshchagin, and S. S. Xue, Cosmic absorption of ultrahigh energy particles, Astrophys. Space Sc. 361 (2016) 82.[46] A. Heger A., S. E. Woosley, The Nucleosynthetic Signature of PopulationIII, ApJ (2002) 532.[47] B. Liu and V. Bromm, When did Population III star formation end?, Mon.Not. Roy. Astron. Soc. (2020) 2839-2854[48] S. Hirano and V. Bromm, Formation and survival of Population III stellarsystems, Mon. Not. Roy. Astron. Soc. (2017) 898-914[49] R. S. de Souza, N. Yoshida, and K. Ioka, Populations iii.1 and iii.2gamma-ray bursts: constraints on the event rate for future radio and x-raysurveys, A&A (2011) A32.[50] S. Ando and K. Sato, Relic neutrino background from cosmological su-pernovae, New J. Phys. 6 (2004) 170.[51] K. Murase, High energy neutrino early afterglows from gamma-ray burstsrevisited, Phys. Rev. D (2007) 123001.[52] I. Esteban, M. Gonzalez-Garcia, M. Maltoni, et al. The fate of hints: up-dated global analysis of three-flavor neutrino oscillations. J. High Energ. Phys. 2020 (2020) 178.[53] M. G. Aartsen, et al. [IceCube Collaboration], Di ff erential limit on theextremely-high-energy cosmic neutrino flux in the presence of astrophysi-cal background from nine years of IceCube data, Phys. Rev. D (2018)062003[54] M. G. Aartsen, et al. [IceCube-Gen2 Collaboration], IceCube-Gen2: TheWindow to the Extreme Universe (2020) arXiv:2008.04323 [astro-ph.HE][55] P. Allison, et al., The Payload for Ultrahigh Energy Observations (PUEO):A White Paper (2020) arXiv:2010.02892 [astro-ph.IM][56] J. A. Aguilar, et al., The Next-Generation Radio Neutrino Observa-tory – Multi-Messenger Neutrino Astrophysics at Extreme Energies (2019)arXiv:1907.12526 [astro-ph.HE].[57] A. Nepomuk Otte, et al., Trinity: An Air-Shower Imaging System for theDetection of Ultrahigh Energy Neutrinos, Bulletin of the AAS (2019) 7,arXiv:1907.08732 [astro-ph.IM].[58] Wissel, Stephanie, et al., Expanding the Reach of Tau Neutrino Tele-scopes with the Beamforming Elevated Array for COsmic Neutrinos (BEA-CON), Astro2020: Decadal Survey on Astronomy and Astrophysics, APCwhite papers, no. 191; Bulletin of the American Astronomical Society ,Issue 7 (2019) id. 191.[59] O. Mena, S. Palomares-Ruiz, and A. C. Vincent, Flavor composition ofthe high-energy neutrino events in IceCube, Phys. Rev. Lett. (2014)091103.[60] M. G. Aartsen, et al. [IceCube Collaboration], Flavor ratio of astrophysi-cal neutrinos above 35 TeV in IceCube, Phys. Rev. Lett. , (2015) 171102.[61] M. Bustamante, J. F. Beacom and W. Winter, Theoretically palatable fla-vor combinations of astrophysical neutrinos , Phys. Rev. Lett. (2015)161302.[62] J. F. Beacom, N. F. Bell, D. Hooper, S. Pakvasa and J. Weiler, MeasuringFlavor Ratios of High-Energy Astrophysical Neutrinos, Phys. Rev. Lett. D (2005) 093005.[63] M. Bustamante and M. Ahlers, Inferring the flavor of high-energy astro-physical neutrinos at their sources, Phys. Rev. Lett. (2019) 241101.[64] M. M. Reynoso, M. C. Medina, and G. E. Romero, A lepto-hadronicmodel for high-energy emission from fr i radiogalaxies, A& A (2011)A30.[65] A. Cooray, Extragalactic background light measurements and applica-tions, Royal Society Open Science (2016) 150555.[66] A. Domnguez, J. R. Primack, D. J. Rosario, F. Prada, R. C. Gilmore, S. M.Faber, D. C. Koo, R. S. Somerville, M. A. Prez-Torres, P. Prez-Gonzlez, andet al., Extragalactic background light inferred from aegis galaxy-sed-typefractions, MNRAS (2010) 2556-2578.[67] M. G. Aartsen, et al. [IceCube Collaboration], An All-sky Search forThree Flavors of Neutrinos from Gamma-ray Bursts with the IceCube Neu-trino Observatory, ApJ (2016) 115.[68] D. Biehl, D. Boncioli, A. Fedynitch and W. Winter, Cosmic-Ray and Neu-trino Emission from Gamma-Ray Bursts with a Nuclear Cascade, Astron.Astrophys. (2018), A101[69] L. Morejon, A. Fedynitch, D. Boncioli, D. Biehl and W. Winter, Improvedphotomeson model for interactions of cosmic ray nuclei, JCAP (2019),007(2019),007