Time dependent photon and neutrino emission from Mkr 421 in the context of the one-zone leptohadronic model

We apply a recently developed time-dependent one-zone leptohadronic model to study the emission of the blazar Mrk 421. Both processes involving proton-photon interactions, i.e. photopair (Bethe-Heitler) and photopion, have been modeled in great detail using the results of Monte Carlo simulations, like the SOPHIA event generator, in a self-consistent scheme that couples energy losses and secondary injection. We find that TeV gamma-rays can be attributed to synchrotron radiation either from relativistic protons or, alternatively, from secondary leptons produced via photohadronic processes. We also study the variability patterns that each scenario predicts and we find that while the former is more energetically favored, it is the latter that produces, in a more natural way, the usual quadratic behavior between X-rays and TeV gamma-rays. We also use the obtained SEDs to calculate in detail the expected neutron and neutrino fluxes that each model predicts.


Introduction
Blazars, a subclass of Active Galactic Nuclei (AGN) emit non-thermal, highly variable radiation practically in all wavebands.The standard scenario is that an active region inside a blazar's jet energizes particles to high energies.As the jet points, within a small angle, towards our direction, the resulting photon emission is boosted due to relativistic beaming.In the case where the accelerated particles include both electrons and protons, then electron synchrotron radiation is responsible for the blazar emission from the radio up to UV/X-ray regime, while the high energy (GeV-TeV) emission is attributed to processes induced by hadrons which include proton synchrotron radiation [1,2] and pion-related cascades [3], producing thus the characteristic double hump blazar Spectral Energy Distribution (SED) in a νF ν vs ν diagram.
While hadronic models are used to fit the SED of AGN, they have not been used so far as to address issues like the SED time variability which is one of the cornerstones in understanding AGN emission.Furthermore, while they are also used to predict AGN neutrino fluxes, they do this mostly by assuming a generic proton and photon distribution, an assumption which inevitably leads to simplifications.
In the present paper we will revisit the neutrino emission and we will address the SED time variations expected from an AGN.For this we will use a recently developed a e-mail: amastich@phys.uoa.grb e-mail: maroulaaki@gmail.comc e-mail: st.seleukos@gmail.comnumerical code, which we will present in §2, that treats the radiative transfer problem in a region where both electrons and protons are accelerated.This code allows us to calculate the proton distribution at source simultaneously with the radiated photon and neutrino spectra.In §3 we will seek successful fits to the SED of the well monitored BL Lac object Mrk 421 which will give us also the production spectra of neutrons and neutrinos.In §4 we will address the X-ray -TeV correlations expected from our model and we will conclude with a discussion in §5.

The model
In this section we will present the principles of the adopted numerical code.While details can be found in [4], for the sake of completeness we summarise here its basic points.We have assumed that high energy protons and electrons are injected uniformly with a constant rate inside a spherical volume of radius R containing a tangled magnetic field of strength B. Both species lose energy by various physical processes: Protons lose energy by synchrotron radiation, photopair and photopion production.All photopion interactions, including the decay of produced pions and muons, have been modelled using the results of the SOPHIA Monte Carlo code [5], while photopair interactions are modelled using the Monte Carlo results of [6].The stable secondary products of these processes include electron/positron pairs (which will be hereafter collectively referred to as electrons), photons, neutrons and neutrinos.Pions and muons, which are unstable secondary products, are allowed to lose energy through synchrotron radiation before their decay.Electrons (primary and secondary) lose energy by synchrotron radiation and inverse Compton scattering, while gamma-rays can be absorbed on the softer photons, produced mainly by electron synchrotron radiation, creating more electron-positron pairs.Neutrons, not being confined by the magnetic fields in the source, can either escape or interact with the photons before decaying back to protons.Finally neutrinos constitute the most uncomplicated component as they escape from the source essentially with their production spectrum.
In order to follow the evolution of this non-linear (see [7]) system we use the kinetic equation approach; thus we write five time-dependent equations for each stable species, namely for protons, electrons, photons, neutrons and neutrinos.The various rates are written in such a way as to ensure self-consistency, i.e. the amount of energy lost by one species in a particular process is equal to that emitted (or injected) by another.That way one can keep the logistics of the system in the sense that at each instant the amount of energy entering the source through the injection of protons and electrons should equal the amount of energy escaping from it in the form of photons, neutrons and neutrinos (in addition to the energy carried away due to electron and proton escape from the source).
The above approach has various advantages: since it uses a particle injection rate, one can calculate the efficiency of the hadronic model -note that the most usual approach of using a particle distribution function in order to calculate the resulting spectra cannot give such estimates.Furthermore, it can calculate in an exact way -under the assumption of a spherical geometry -the photon, neutron and neutrino transport and thus compare directly their respective fluxes.
The free parameters used are the injected luminosities of protons and electrons, which are expressed in terms of compactnesses: where i denotes protons or electrons and σ T is the Thomson cross section; the upper Lorentz factors of the injected protons and electrons, γ p,max and γ e,max respectively; the power law indices p p and p e of injected protons and electrons, respectively; the radius R of the region; its magnetic field, B; its Doppler factor, δ; and the escape time for both particles, which is assumed to be the same (t p,esc = t e,esc ).Finally, the lower Lorentz factors are assumed to be γ p,min = γ e,min = 1.Depending on the assumptions made about the time dependence of the parameters, the above scheme can be used to derive both steady-state and time-dependent solutions -see [8].Thus, if all parameters are constant in time the system will eventually reach a steady-state -note however that hadronic plasmas can become supercritical and in such cases they can exhibit limit cycle behaviour [9] even if the parameters are time independent.If, on the other hand, the system is in the subcritical regime and we allow for one or more parameters to have some explicit time dependence, then the system will not reach a steady state but it will show continuous temporal variations, which will reflect the corresponding ones imposed on the input parameters and will have an impact on the produced spectrum.In analogy to similar efforts in leptonic models ( [10], [11]), we will examine such a case in §4.

Photon and neutrino spectra from Mrk 421
Being one of the closest blazars to Earth, Mrk 421 has been the target of multiple multiwavelength campaigns over the years.Observations by [12] in 2001 produced excellent sets of time-dependent data both at the X-ray (RXTE) and TeV (Whipple and HEGRA) regimes during a six-day period.In the present section we focus on the pre-flare steady state emission of Mrk 421 in order to examine the resulting neutrino and neutron distributions in relation to the observed photon spectra.Interestingly enough we have found two sets of very different proton injection parameters which give very good fits to the data.In both cases, optical and X-rays are fitted by the primary injected leptonic component, while the origin of the GeV-TeV emission is different between the two models.We will discuss each one separately, since they imply very different physical conditions in the source.

The Leptohadronic-pion (LHπ) model
In this model the TeV data are fitted by the synchrotron radiation of electron/positron pairs that result both from charged pion decay and from the absorption of γ−rays from neutral pion decay.We will refer to this model as 'LHπ' (Leptohadronic-pion).Its initial parameters are shown in Table 1 while its fit to the SED is depicted in Fig. 1 along with its neutrino spectrum.The combination of a low magnetic field with a high proton injection compactness results in a compressed proton synchrotron component and prominent photopair and photopion components.Thus the SED does not have the usual double hump appearance as synchrotron photons from the photopair secondaries produce a distinctive broad hump at MeV energies.In general, we find that the inclusion of photopair -a process that is not usually taken into account -as an extra mechanism for injecting secondaries, alters the shape of the radiated photon spectra as for most 05005-p.2Fermi observations [13] (grey points) are not simultaneous with the rest of the data.The 40-String IceCube limit for muon neutrinos [14] is plotted with a dash-dotted line.

EPJ Web of Conferences
relevant fitting parameters it has approximately the same significance as the widely used photopion process.The neutrino spectrum is sharply peaked at an energy about 30 times lower than the maximum proton energy in the νF ν vs ν diagram.At lower energies it follows an approximate power law with index p ν , which is harder than the power law of the initial protons by a factor of ∼ 1, in accordance with the approximate relation p ν ≃ (p p − 0.5)/2.5 from [4].Neutrinos from neutron decay peak at an energy two orders of magnitude lower than those from meson decay and their luminosity is similarly smaller.In this case their contribution is noticeable as a small peak to the left of the main neutrino spectrum peak in Fig. 1.
The LHπ model produces a substantial neutrino flux that is of the same order as the TeV γ−rays and peaks at energies around E ν,peak = 3 PeV.Its spectral parameters, compared to those of other particle types, are summarized in Table 2.Note that the luminosities are given in the comoving frame.The expected neutrino flux is just under the sensitivity of the IC-40 detector and should be producing observable neutrinos when its activity is high.Recent observations of PeV-energy neutrinos by the IceCube collaboration [15] are in good agreement with this prediction, although the calculated flux is still too low to offer any hint at a spectral shape.

The Leptohadronic-synchrotron (LHs) model
In this case, TeV γ−rays are produced by proton synchrotron radiation.We will refer to this model as 'LHs' (Leptohadronic-synchrotron).Its initial parameters are also shown in Table 1 while its fit to the SED and its neutrino spectrum are depicted in Fig. 2  Although neutrinos from this model are well below current instrumental sensitivities, the ultra-high energy protons produced by the escaping neutron decay could be observable.Neutrons resulting from photopion interactions are an effective means of facilitating proton escape from the system, as they are unaffected by its magnetic field and their decay time is high enough to allow them to escape freely before reverting to protons [16][17][18][19].A further advantage is that they are unaffected by adiabatic energy losses that the protons may sustain in the system as it expands [20].Those effects make them excellent originators of ultra-high energy cosmic rays (UHECR).In Fig. 3 we show the neutron spectrum as it would appear at Earth, having first decayed into protons.Their propagation in a uniform intergalactic pG magnetic field and The Innermost Regions of Relativistic Jets and Their Magnetic Fields their energy losses from interactions with the cosmic microwave and infrared background are modelled using CR-Propa [21].Since at this energy range those losses are the dominant factor in reshaping the proton spectrum and diffusion effects are minimal, we have restricted ourselves to one-dimensional simulations.The spectra (black dots) are compared to observations from Auger [22] (open triangles), HiRes-I [23] (open squares) and Telescope Array [24] (x's).The resulting UHECR spectra, neglecting any energy losses and propagation effects, are plotted for comparison reasons with a dashed line.As Mrk421 is a northern source, its contribution to UHECR could explain some of the discrepancy between the Auger and HiRes-I/Telescope Array data (if this is, indeed, an intrinsic difference) since its position makes it invisible to Auger.Interestingly enough, the same applies for every other highfrequency-peaked BL Lac within a distance of z = 0.05 [25].
We should note here that this is the minimum contribution of UHECR that this model predicts.It is quite possible that the main high energy population of protons inside the source contributes as well once it escapes.However, this will depend on whether the protons suffer strong adiabatic losses before escaping the blazar jet.

SED time variability
In the previous section we showed that two different hadronic models can give good fits to the SED of Mrk 421 as this was obtained by [12].Thus, it is possible that each fitting combination will have distinct temporal signatures which are worth investigating.Since the observing sessions mentioned above have given good correlations between the X-rays and TeV γ−rays, it would have been interesting to examine the expected correlations that  each model produces and deduce which one comes to better agreement with the observations.Motivated by the results of long-term variability studies of Mrk 421 and other prototype blazars (e.g.[26]), we have introduced, for the temporal variations, a randomwalk type of change in the injection luminosity of particles.Thus, we use where ξ is a uniformly distributed random number in the range (0,10).The integer parameter a k is then scaled in such a way as to produce the following change in the injection parameters where ℓ inj k is the value of the injection parameter at time t k , ℓ inj 0 is the value of the injection parameter that gave the low-state fit and f a multiplication factor -in all the examples shown in the present paper we have chosen f = 0.05.For the values of the injection parameter between two successive crossing times, we have chosen a linear interpolation scheme.
Previous analyses of Mrk 421 have noted the presence of lags in its spectral evolution (e.g.[12,27,28]).In our case, where two primary particle populations are being injected into the source, these lags can be simulated by introducing a shift N in the temporal variation of ℓ inj e compared to ℓ inj p , i.e. (ℓ inj e ) i = (ℓ inj p ) i+N , where the subscript denotes that quantities are calculated at time t i (in units of t cr ).
Figure 4 depicts the TeV vs. X-ray flux obtained in the LHπ case.In both panels, the black line corresponds to N = 0, i.e. ℓ this way we imposed less correlated variations on particle injection, having a Pearson's correlation coefficient r = 0.63.We have used a variety of positive correlated ℓ inj p , ℓ inj e by introducing various shifts.In all cases, which we do not present here, we derived large correlation coefficients (r 0.5) which show that the TeV and X-ray fluxes retain their strong correlation.The respective diagram for the LHs model and N = 80 are shown in Fig. 5.
In general, the introduction of a shift loosens the F TeV −F X correlation.The effect, however, is more evident in the proton synchrotron case where the TeV/X-ray flux correlation is almost destroyed for N = 80 (r = −0.15).On the other hand, in the LHπ model, we find a larger absolute value for the correlation coefficient.In other words, if the γ-rays are modelled by the emission of secondaries produced by photohadronic processes and therefore any variations of F TeV reflect only indirectly the changes in the proton injection rate, the correlation between the fluxes in the X-rays and TeV γ-rays is partially retained.On the other hand, in the LHs model the variability observed in both X-rays and γ-rays reflects directly the variability pattern of the particle injection rate.Thus, if there is a degree of decorrelation in the injection it will be seen also in the flux-flux diagram.For the LHπ case (top panel) the introduction of a shift decreases the maximum/minimum flux values in both energy bands, since the γ-ray luminosity depends also on the number density of soft target photons, e.g.synchrotron photons in the X-rays.In the LHs case on the other hand, where the γ-rays are the product of proton synchrotron radiation, the γ-ray flux does not depend on the X-ray photons, as long as the X-ray photon number density is low enough as not to cause significant γγ absorption.Thus, the range of flux variations remains approximately the same.
Comparing the above synthetically derived correlations to the observed ones that showed, in most cases, a quadratic trend between TeV γ−rays and X-rays, we deduce that the LHπ model can fit in a more natural way this trend.On the other hand, the LHs requires a rather fine tuning, i.e. it requires that the injection of protons should vary quadratically with respect to the one of the electrons -for more details see [8].

Summary/Discussion
In the present paper we have used an one-zone leptohadronic model to obtain fits to the SED of the BL Lac object Mrk 421. as obtained in 2001 using contemporaneous optical, X-ray and TeV γ−ray data [12].Even if the assumed geometry of the source within an one-zone model is simple enough, we have applied a sophisticated method to calculate in a self-consistent and time-dependent way the photopair and photopion interactions -see [4].Consequently, we have used the so-obtained fitting parameters in a two-fold way: On the one hand we have calculated the expected neutron and neutrino fluxes and we have followed the propagation of protons, which result from the escaping neutron decay to derive the expected cosmic ray contribution of Mrk 421 at Earth.On the other, we have introduced perturbations in the injection particle lumimosities in order to study the temporal variations of the system and compare them with observed ones from Mrk 421.
We found that two sets of very different parameters produce good fits to the SED of the source.In the first case, γ−rays are produced from synchrotron radiation of secondaries resulting from photopion interactions, therefore they are pion-induced (the LHπ model).In the other case, the γ−rays originate from proton synchrotron radiation (the LHs model).
The LHs model requires proton acceleration at E p,max = 1.4 × 10 20 eV which agrees well with the requirement for UHECR accelerators.It is a very economic model in the sense that it requires a relatively low jet luminosity and it is magnetically dominated since u B ≫ u p -see Table 2. Propagation of the protons produced from the escaping neutron decay results in a UHECR flux at Earth that is very close to the measurements of HiRes, Auger and Telescope Array at energies around 30 EeV.However, due to the fact that our Cosmic Ray (CR) spectrum is peaked at high energies, it fails to reproduce the observed CR spectrum at energies below 10 EeV -see Fig. 3.This situation cannot be improved if we assume that the other Northern Hemisphere nearby BL Lac objects produce the same spectral shape of CRs as Mrk 421 and normalize their output to the photon luminosity of this source.The combination of the sources' low luminosities and CR propagation through larger distances minimizes their contribution to the CR flux at Earth.Furthermore, the LHs model produces a low neutrino flux, since photohadronic processes are suppressed to a low level, that is many orders of magnitude below the IceCube sensitivity threshold -see Fig. 2. Finally, it requires a fine-tuning between the injection of electrons and protons to obtain the usual quadratic behaviour between X-rays and TeV γ-rays -see Section 4.
The LHπ model, on the other hand, requires a large, but not unacceptable, jet power which is heavily particle 05005-p.5 The Innermost Regions of Relativistic Jets and Their Magnetic Fields dominated -see Table 2. Good fits to the SED of the source are obtained assuming that the protons are accelerated up to E p ≃ 30 PeV, therefore they cannot contribute to the UHECR flux.However, the produced neutrino flux is tantalizingly close to the IC-40 sensitivity limit for Mrk 421 [14] and peaks within the energy range where the first PeV neutrino detection from the IceCube collaboration was reported [15].A shortcoming of the LHπ model is its low efficiency.If we define the total efficiency ξ as the ratio of the sum of escaping luminosities in photons, neutrons and neutrinos over the sum of the injected luminosities, then we find ξ ≈ 2 × 10 −5 .However it can expalin in a more natural way the observed X-ray -TeV correlations as the X-rays become targets for the photopion induced gamma-rays.
The novelty of the present work is that instead of modelling a generic blazar, we have focused on a specific source.Thus, the distribution of ultra-high energy protons that we use in order to derive the variability and the neutrino/CR spectra is specifically determined by fitting the SED of the high-peaked blazar Mrk 421.As such it is limited to these sources.Efforts to fit the SED of the more luminous, low-peaked blazars would probably have resulted in different parameters as the role of the Broad Line Region, which in the case of Mrk 421 is marginal, could have been central.We plan to extent our analysis to these sources in the near future.

Figure 1 .
Figure 1.Spectra of photons (thin line) fitting the March 22nd/23rd 2001 observation of Mrk 421 (black points), neutrinos of all flavours (grey line) and muon neutrinos (thick line) escaping from the source, according to the photopion model LHπ.Fermi observations[13] (grey points) are not simultaneous with the rest of the data.The 40-String IceCube limit for muon neutrinos[14] is plotted with a dash-dotted line.

Figure 2 .
Figure 2. As in Fig. 1 but for the LHs model

Figure 3 .
Figure 3. Cosmic ray spectra obtained within the LHs model after taking into account the effects of propagation using the numerical code CRPropa[21]; the uncorrected high-energy neutron spectrum is plotted with a dashed line.The cosmic ray energy spectrum, as observed by Auger[22], HiRes-I[23] and Telescope Array[24] is overplotted with grey open triangles, black open squares and grey x's respectively.

Figure 4 .
Figure 4. Plot of the TeV vs. X-ray fluxes obtained within the LHπ model, after varying ℓ inj p and ℓ inj e with N = 0 (black line) and N = 80 (grey line).Inset legends show the Pearson's correlation coefficients for each case.

Figure 5 .
Figure 5.As in Fig. 4 for the case of the LHs model.

Table 1 .
Initial parameters for the two fits.

Table 2 .
. The high magnetic field coupled with a low proton injection compactness results in a compressed photopair and photopion component Neutrino parameters compared to the parameters of the protons from which they originated.Energy densities are at steady state; those of the other particles and the magnetic field are included.