A Unified Model for Multi-epoch Neutrino Events and Broadband Spectral Energy Distribution of TXS 0506+056

The blazar TXS 0506+056 has been proposed as a high-energy neutrino emitter. However, it has been shown that the standard one-zone model cannot produce sufficiently high neutrino flux due to constraints from the X-ray data, implying more complex properties of the radiation zones in the blazar than that described by the standard one-zone model. In this work we investigate multi-epoch high-energy muon neutrino events associated with the blazar TXS 0506+056 occured in 2014–2015, 2017–2018, 2021–2022 and 2022–2023, respectively. We applied the so-called “stochastic dissipation model” to account for the neutrino-blazar associations detected in the four epochs simultaenously. This model describes a scenario in which the emission of the blazar arise from the superimposition of two components: a persistent component related to the quasi-stable state of the blazar and a transient component responsible for the sudden enhancement of the blazar’s flux, either in electromagnetic radiation or in neutrino emission. The latter component could form at a random distance along the jet by a strong energy dissipation event. Under such assumption, the multi-epoch broadband spectral energy distribution (SED) can be well explained and the expected number of high-energy neutrino events is statistically realistic. The expected number of neutrino events in half-year is around 8.2, 0.07, 0.73 and 0.41, corresponding to the epoch in 2014–2015, 2017–2018, 2021–2022 and 2022– 2023, respectively. Hence, our model self-consistently explains the episodic neutrino emission from TXS 0506+056.


INTRODUCTION
The origin of high-energy neutrinos is still unclear (Aartsen et al. 2013;IceCube Collaboration 2013).As the most powerful persistent emitters of electromagnetic (EM) radiation in the universe, blazars have long been suggested as potential accelerators of high-energy cosmic-ray protons (or nuclei) and sources of high-energy neutrinos (e.g., Mannheim 1993Mannheim , 1995;;Stecker & Salamon 1996;Halzen & Hooper 2002;Atoyan & Dermer 2001;Murase et al. 2014).The discovery of a possible correlation between high-energy neutrino events and blazars supports this speculation.(IceCube Collaboration et al. 2018a,b).Many researchers have studied these events in the framework of the so-called photohadronic interaction model (Ansoldi et al. 2018;Keivani et al. 2018;Murase et al. 2018;Padovani et al. 2018Padovani et al. , 2019;;Cerruti et al. 2019;Gao et al. 2019;Reimer et al. 2019;Rodrigues et al. 2019;Xue et al. 2019;Abbasi et al. 2021) and hadronuclear interaction model (Sahakyan 2018;Liu et al. 2019;Banik & Bhadra 2019).In these hadronic processes, the produced secondary EM particles can induce EM cascades, contributing to EM radiation from the optical band to the X-ray band or even to the MeV gamma-ray band (Xue et al. 2021;Reimer et al. 2019).On the other hand, the low-energy radiation of a blazar can be the target radiation field for the photohadronic interactions Corresponding author: Ruo-yu Liu; Junfeng Wang ryliu@nju.edu.cn;jfwang@xmu.edu.cn(Kelner & Aharonian 2008;Murase et al. 2014;Murase & Stecker 2022).In this sense, the neutrino emission from a blazar is related to its broadband radiation.
The first blazar associating with neutrino events at a significance above 3σ, TXS 0506+056, is a known gamma-ray blazar at redshift z = 0.3365 (Paiano et al. 2018).On 2017 September 22, a 290 TeV muon neutrino event (known as IC-170922A), was detected by IceCube (IceCube Collaboration et al. 2018a) from the direction of the blazar, when the blazar underwent a GeV gamma-ray flare.The chance coincidence of the association is rejected at a significance level of approximately 3σ.IceCube also found indication of a neutrino outburst with 13 ± 5 muon neutrino event in 32  PeV from this blazar over ∼ 160 days from 2014-2015 at a significance of 3.5σ level when examining over nine years of archival data (IceCube Collaboration et al. 2018b).Great efforts have been made to understand the theoretical correlation between the neutrino emission and the blazars (Mannheim 1995;Stecker & Salamon 1996;Atoyan & Dermer 2001;Murase et al. 2014;Murase & Stecker 2022).An consensus has been reached that the standard one-zone model cannot produce sufficiently high neutrino flux due to constraints from the X-ray data (Keivani et al. 2018;Gao et al. 2019;Rodrigues et al. 2019), and the detection of IC-170922A has to be interpreted as the Eddington bias (Strotjohann et al. 2019) or an upward fluctuation in the detection, given an anticipated event rate of 0.01/year.The problem for the 2014-2015 outburst is more serious.Rodrigues et al. (2019) found that the standard one-zone model can at most generate 4.9 events in half a year.This discrepancy cannot be ascribed to the Eddington bias because of the large difference in the event number between expectation and observation.More complex models invoking two or even more radiation zones are needed if we want to reproduce a more comfortable neutrino event rate to account for these association, especially for the 2014-2015 outburst (Xue et al. 2021;Rodrigues et al. 2019).Later, an independent search for point-like sources in the northern hemisphere using ten years of IceCube data revealed that TXS 0506+056 is coincident with the second hottest-spot of neutrino event excess (Aartsen et al. 2020).In addition, interestingly, a cascade neutrino event with the estimated energy of 224 ± 75 TeV, GVD-210418CA, from the direction of TXS 0506+056 was detected by the new Baikal-GVD neutrino telescope in April 2021 (Baikal-GVD Collaboration et al. 2022).On September 18, 2022, a ∼ 170 TeV neutrino, IC-220918A, was reported in spatial coincident with TXS 0506+056 again by IceCube (Becker Tjus et al. 2022a).Although the latter two neutrino events have not been studied in detail, the maximum expected neutrino event rate under the standard one-zone model would be similar to that for IC-170922A, i.e., ∼0.01/year, because the X-ray fluxes in the latter two epochs are comparable to (or slightly lower than) that in the epoch IC-170922A.If the latter two neutrino events are truly associated with TXS 0506+056, the standard one-zone model can be ruled out, because the joint possibility would be unreasonably low if all these associations are explained as lucky detections.
The motivation of this work is to explore whether the multi-messenger emission can be self-consistently accounted for in all four epochs.Also, as can be seen later, the SED in these four epochs are different.Hence, it is also interesting to investigate if there exists an inherent relation among the four epochs with neutrino detection.Wang et al. (2022) proposed the so-called "stochastic dissipation model" to explain the orphan and multiwavelength flares from blazars in a unified framework.They speculate that a blazar's emission is composed of contribution from numerous radiation zones (or blobs), formed by certain energy dissipation events, along the jet.Therefore, the blazar's emission consists of two components during a flare: (i) a quasi-stable component that arises from the superposition of numerous but comparatively weak dissipation zones along the jet, forming a persistent emission of the blazar which is statistically identical among different epochs (Meyer et al. 2019); (ii) a transient component, which is responsible for the sudden enhancement of the blazar flux, generated at a random distance along the jet by a strong energy dissipation event (Giannios & Spruit 2006;Nalewajko & Begelman 2012;Begelman 1998;McKinney & Blandford 2009).In this work, we explore the broadband SED and the neutrino emission of TXS 0506+056 in multi-epochs in the framework of this model, assuming the presence of high-energy protons in the dissipation events.
The rest of the paper is organized as follows.The observation and data reduction process are described in Section 2. The model and method for calculating the radiation is described in Section 3. In Section 4 we show the SED modeling and results.Finally we present discussion and summarize our findings in Section 5. Throughout the paper, H 0 =71 km s −1 Mpc −1 , Ω m = 0.27, and Ω Λ = 0.73 are adopted.

OBSERVATIONS AND DATA REDUCTION
For the 2014-2015 neutrino outburst and the IC-170922A flare, the multiwavelength data are from Rodrigues et al. (2019) and IceCube Collaboration et al. (2018a), respectively.For the periods covering the 2021-2022 neutrino event GVD-210418CA and the 2022-2023 neutrino event IC-220918A, the multiwavelngth data are processed by ourselves in this work.

γ-Ray Observations by Fermi-LAT
We first obtained the Pass 8 data from the Fermi Science Support Center (FSSC).The data are selected within a 15 • region of interest (ROI) centered on the radio position of TXS 0506+056 (R.A. = 77.3582• , Decl.= 5.69315 • ; Beasley et al. (2002)).For the 2021-2022 neutrino events, the temporal coverage of the data is from 2021 January 01 to 2021 December 31 (MJD 59215-59580) of about 1.0 yr.For the 2022-2023 neutrino events, the temporal coverage of the data is from 2022 January 01 to 2022 December 31 (MJD 59580-59944) of about 1.0 yr.A binned likelihood analysis is performed for the γ-rays of 4FGL J0509.4+0542 using the publicly available software Fermitools v.2.0.8.In our analysis only the γ-ray photons in the energy range of 0.1-300 GeV and satisfying the standard data quality selection criteria "(DATA QUAL>0)&&(LAT CONFIG == 1)" are considered.A zenith angle cut of 90 • is set to avoid the γ-ray contamination causing by the Earth limb.We bin the data with a pixel size of 0.2 • in space and 25 logarithmical energy bins.The background models include all γ-ray sources listed in the 4FGL-DR3 Catalog and the Galactic diffuse component (gll iem v07.fits) as well as the isotropic emission (iso P8R3 SOURCE V3 v1.txt).The P8R3 SOURCE V3 set of instrument response functions (IRFs) is used.
The model for fitting the spectrum of 4FGL J0509.4+0542 in our analysis is a single power-law function, where Γ γ is the photon spectral index, and E b is the scale parameter of photon energy (Massaro et al. 2004).The spectral parameters of all sources lying within 8 • are left free, whereas the parameters of those sources lying beyond 8 • are fixed to their 4FGL-DR3 values.Also, the normalization parameters of the standard Galactic and isotropic background templates are set free in the likelihood fit.The significance of the γ-ray detection is quantified by adopting the maximum likelihood test statistic (TS), TS = 2(L 1 -L 0 ), where L 1 and L 0 are the maximum likelihood values for the models with and without the target source, respectively.The sources with TS < 16 are removed from further analysis.The results of the analysis and SED fitting are shown in Table 1 and Figure 7-8.The TS maps are shown in Figure 11-12.

X-Ray Observations with Swift-XRT
The X-ray telescope (XRT) on board the Swift satellite observed the source in the photon counting mode with exposure times of 21.1 ks and 2.5 ks, in 2021-2022 and 2022-2023 respectively.We perform standard filtering and data analysis using HEASOFT (v 6.29).To extract the source spectrum, we select a circular region of 30 ′′ , centered at the target.The background is estimated from an annular cell.The source spectrum is binned to have at least 20 counts per bin.The SED fitting process is performed in XSPEC and the χ 2 minimization technique is adopted for spectral analysis.The associated errors are calculated at the 90% confidence level.The spectrum is fitted by a single power law absorbed by two absorption components, one is absorption at z = 0 with the neutral hydrogen column density fixed at Galactic value N H gal = 1.55 × 10 21 cm −2 (Willingale et al. 2013), the other is an extragalactic foreground absorption N H int at a redshift of z = 0.3365 with column density set free.An extragalactic foreground absorption column density of N H int would be obtained by SED fitting.The complete process of building Swift-XRT products are described in https://www.swift.ac.uk/user objects/.The results of the analysis and SED fitting are shown in Table 2 and Figure 9-10.

Optical Observations with ALeRCE
The Automatic Learning for the Rapid Classification of Events (ALeRCE, Förster et al. 2021) broker light curve classifier (LCC, Sánchez-Sáez et al. 2021) is processing the alert stream from the Zwicky Transient Facility (ZTF) and which aims to become a Community Broker for the Vera C. Rubin Observatory and its Legacy Survey of Space and Time (LSST), as well as other large etendue survey telescopes.The goal of the LCC is to provide a fast classification of transient and variable objects by applying a balanced random forest algorithm.The LCC produced ariability features and colors obtained from AllWISE and ZTF photometry in g and r band.The complete set of features are described in http://alerce.science/features/.In this work, we extracted variability features from the g and r bands to follow the neutrino events of TXS 0506+056 during the periods of 2021-2022 and 2022-2023, respectively.6  2022-2023 2.20 ± 0.002 (7.19 ± 0.004) × 10 −11 868   3. MODEL SETUP AND METHOD Liu et al. (2023) and Wang et al. (2022) proposed the stochastic dissipation model to explain various observational features of blazars.The model suggests that there are numerous discrete radiation zones along the jet owing to certain dissipation events, such as magnetohydrodynamic instabilities (Giannios & Spruit 2006;Nalewajko & Begelman 2012), internal collisions (Deng et al. 2015;Joshi & Böttcher 2011;Böttcher & Dermer 2010), recollimation shocks (Nalewajko 2012;Fromm et al. 2016), magnetic reconnections in the jet (Begelman 1998;McKinney & Blandford 2009;Giannios et al. 2009;Petropoulou et al. 2016) and so on, which we do not specify in this work.The probability of a dissipation event occurring at a distance r from the supermassive black hole (SMBH) in unit time and unit length is phenomenologically characterised with a function p(r) ∝ r −α .The value of α and the injection electron spectral index into each blob can largely determine the spectral shape expected from a blazar.
The stochastic dissipation model posits that the emission of the blazar jet can be decomposed into two components (Wang et al. 2022;Tan et al. 2023).One component is the superposition of emission from numerous but comparatively weak dissipation zones along the entire jet, which is considered to be a quasi-stable, background emission.Another component arises from a strongly dissipating zone in the jet that dominates the emission in the flare state.The historical SED data taken from the SSDC SED builder1 representing a long-time average emission of the blazar is then considered as the low-state background emission.Since we here focus on the emission from a flaring zone in the jet to account for the the production of neutrinos, the modelling of the background radiation spectrum is simply characterized with a polynomial function (Wang et al. 2022).Liu et al. (2023) developed a time-dependent model for the low-state background emission, and the result is compatible with the polynomial function approximation.
To research the neutrino events of TXS 0506+056 in different epochs, we model the broadband SED in each epoch, based on which we can calculate corresponding neutrino detection rate in the epoch.In our model we consider the electron synchrotron radiation, the inverse Compton scattering off the synchrotron radiation field (i.e., synchrotron self-Compton, SSC, Harris & Krawczynski 2006) and off the external radiation (i.e., external Compton, EC) of electrons.We also take into account the proton synchrotron radiation, the photopion process, and the Bethe-Heitler process, as well as the emission of pairs generated in the electromagnetic cascade initiated by these processes.
We assumed that the flaring zone is a blob which has a spherical geometry with a radius R, filled with a uniformly entangled magnetic field B. The jet's bulk Lorentz factor Γ is assumed to remain constant up to 0.1 pc and we approximate the Doppler factor by δ = Γ, δ 0 is the Doppler factor at 0.1 pc.Considering a truncated conical jet structure and a dissipation zone of radius R comparable to the transverse radius of the jet at certain distance r from the SMBH (Wang et al. 2022), we have where R 0 is the transverse radius of the jet at 0.1 pc.Then we can assume that the magnetic luminosity is constant along the jet, which is consistent with results of the VLBA survey (O'Sullivan & Gabuzda 2009; Sokolovsky et al. 2010), the magnetic field strength can be approximated (Wang et al. 2022) as a function of r where B 0 is the magnetic field strength of the dissipation zone for r = 0.1 pc.We assumed that relativistic electrons or protons are injected in the blob with a broken power-law distribution (Ghisellini et al. 2010;Wang et al. 2020) where t e(p) is the minimum value among t cool and t e(p),dyn .
is the energy density of the magnetic field and U ph is the energy density of the soft photons.Here σ T is the Thomson scattering cross section and κ KN is a factor accounting for the Klein-Nishina effect (Moderski et al. 2005).For protons, we need to include the energy loss caused by the photopion production and the Bethe-Heitler (BH) pair production when calculating the cooling timescale.These two processes will be further described later.t e(p),dyn = R/c is the electron(proton) dynamical timescale, here we assumed the radius of the blob is R = r sin θ, r is the distance between the blob and the black hole.In this work all timescales are evaluated in the comoving frame of the blob.
Under the unified model, the synchrotron, SSC and EC radiations are calculated following Wang et al. (2022).If the blob (flaring zone) is close to the supermassive black hole (SMBH), the radiation of the broad line region (BLR) and the dust torus (DT) could enter the blob with a significant amplification of the energy density due to the Doppler effect.The radiation field of BLR is assumed to be a blackbody emission peaking at a frequency 2 × 10 15 Γ Hz (Tavecchio & Ghisellini 2008), and the radiation field of DT is assumed to be a blackbody emission peaking at a frequency 3 × 10 13 Γ Hz (Cleary et al. 2007) in the jet comoving frame, respectively.The energy density of external photons from the BLR and DT can be approximated (Hayashida et al. 2012) by and where η BLR = η DT = 0.1 are the fractions of the disk luminosity L d reprocessed into the BLR and DT radiation, respectively.We assumed the characteristic distance of BLR in the AGN frame is r BLR = 0.1 L d /10 46 erg s −1 1/2 pc and the characteristic distance of DT is r DT = 2.5 L d /10 46 erg s −1 1/2 pc.The disk luminosity L d is assumed to be 5 × 10 44 erg s −1 (Xue et al. 2021).Meanwhile, the X-ray photons emitted by the corona surrounding the accretion disk (Heckman & Best 2014) are also considered in our model, if the blob is near or at the jet base with a distance comparable to only a few times larger than the Schwarzschild radius (r Sch ∼ 10 14 (M SMBH /3 × 10 8 M ⊙ )cm) of the central SMBH, the X-ray photons would interact with gamma-ray photons emitted by the blob (Righi et al. 2019).
We assumed the emission of the corona has a power-law spectrum (Ricci et al. 2018;Xue et al. 2021) here we simply assume that its spectral index α = 1 and L 1 keV is set as 5 × 10 43 erg s −1 for TXS 0506+056 (Xue et al. 2021).The corresponding energy density in the comoving frame can be written as where r is the distance between the central black hole and the dissipation region (or the blob).Eq. ( 10) is valid for r ≫ r Sch , where as for r < 10r Sch we assume u corona =u corona (r = 10r Sch ).
Note that the relativistic protons can also interact with these radiation field via the photopion process and Bethe-Heitler process.The photopion process is the production of pions in proton-photon interactions: via which unstable pions are produced and pions will decay into Another important process is the Bethe-Heitler process, this process will lead to the production of electron/positron pairs, i.e., p Kelner & Aharonian ( 2008) developed a semi-analytical method to calculate spectra of secondary particles in these two processes and we followed the method in Kelner & Aharonian (2008).The cooling timescales of both electrons and protons via the aforementioned processes in the model are shown in Figure 1 for references.
The maximum Lorentz factor of protons in the emission region can be approximated by equating the acceleration and the cooling or dynamical timescales t acc = min {t pγ , t BH , t dyn , t p,syn } , ( The acceleration timescale can be evaluated by (Rieger et al. 2007) where α(≥ 1) is the parameter characterising the particle acceleration efficiency (α = 1 being the most efficient case), which depends on the spectrum of magnetic turbulence and on the velocity of the scattering center in the case of Fermi-type acceleration (Sikora 2011).Here we employ α = 10 in our model.The kinetic luminosity of the nonthermal particles and magnetic field can be written as (Celotti & Ghisellini 2008) where U i is the energy density of nonthermal particles or magnetic field, here i = {e, p, B}, represents electrons, protons and magnetic field respectively.

SED MODELING AND RESULTS
In this section, we use the stochastic dissipation model to reproduce the mutiwavelength emission of TXS 0506+056 in the four epochs with neutrino detection.For the photopion process and the BH process, radiation fields of the blazar serve as targets, including the emission of the BLR, the dusty torus and the corona, the radiation of primary electrons, as well as the radiation of secondary pairs developed in the electromagnetic cascade initiated by the two processes.The BH process and photopion process including the cascade initiated by these processes are calculated following the method shown in (Kelner & Aharonian 2008;Böttcher et al. 2013).
To reduce the number of free parameters, value of γ e,min in different neutrino events of TXS 0506+056 are fixed to 50.Our model has four parameters (R 0 , B 0 , n e,1 , and δ 0 ) that are common among different zone, and four parameters ( L inj e , n e,2 , γ break and r) that are unique to each dissipation zone.γ e,max and γ p,max are obtained by equating the acceleration and the cooling or dynamical timescales.Then we assumed protons are injected in the blob with a powerlaw distribution.The spectral index n p and the minimum Lorentz factor γ p,min of the proton distribution are fixed to 2 and 1 respectively, while the maximum injection luminosity of protons is obtained when the pair cascade emission saturates the X-ray or γ-ray data.We have a total number of 20 parameters in the model, with 4 common for all four epochs and 4 parameters that are separate for each of them.The best-fit and uncertainty of the model parameters are obtained via the Markov Chain Monte Carlo (MCMC) method with the PYTHON package EMCEE (Foreman-Mackey et al. 2013).

Epoch I: 2014-2015
The SED fitting result is shown in Figure 2 and the derived parameters are listed in Table 3.The best-fit blob distance is 3 × 10 −4 pc, or approximately 10 r rsc , from the SMBH, which is located within the corona.
At this distance the X-ray photons emitted by the hypothetical corona surrounding the accretion disk is very dense and absorb the γ-ray photons (Kamraj et al. 2018;Righi et al. 2019;Ricci et al. 2018) efficiently, so the observed γ-ray emission is mainly contributed by the background emission.This explains why the blazar is in the gamma-ray quiescent state during the neutrino outburst.Moreover, the X-ray corona of SMBH could be a promising neutrino production site provided that protons are accelerated in the region.The number of predicted muon neutrino events detected by IceCube during half a year is 8.2 based on the best-fit parameters, which is basically consistent with IceCube's measurement (IceCube Collaboration et al. 2018b).Xue et al. (2021) used a two-zone radiation model of blazars to explain the broadband SED of TXS 0506+056 during 2014-2015, the derived distance between the black hole and the inner blob is also approximately several times the Schwarzschild radius.In their model, the magnetic field of the inner blob is taken to be a free parameter and is only 2 G.The value is much smaller than the best-fit value 83.3 G in our work, which is based on a B ∝ r −1 relation and a simultaneous fitting to SEDs of other three epochs.The strong magnetic field makes synchrotron radiation comparably efficient as the EC process, or even the dominant radiation process for electrons above GeV energies considering the    KN effect (see Figure 6).The synchrotron radiation of pairs developed in the EM cascade creates a bump around MeV band in the SED.This is consistent with the results in some previous literature for the neutrino outburst (Xue et al. 2021) and could be hopefully detected by the next-generation MeV instrument as a test for the model.At GeV band, the high opacity of the γγ absorption due to the X-ray corona photon is very high for the blob.Therefore, the blazar's GeV emission is dominated by the background component, which is consistent with observation.

Epoch II: 2017-2018
The best-fit distance of the blob for this epoch is located at 0.09 pc from the SMBH and is the largest one among the four epochs.This distance is far beyond the size of the corona and comparable to the typical radius of the BLR.The external radiation field in this epoch is the BLR radiation and the DT radiation.The absorption of gamma-ray photons in the blob is not as important as that for Epoch I.The corresponding magnetic field in this blob is 0.28 G, which is similar to those found in many previous studies for the multi-messenger emission of IC-170922A (Xue et al. 2019;Liu et al. 2019;Das et al. 2022;Cerruti et al. 2019;Gao et al. 2019;Gasparyan et al. 2022).
The EM radiation of the blob is important in this epoch, where the optical flux is comparable to that of the background emission and the gamma-ray emission dominates.This is consistent with the observed flux enhancement in the optical band and the gamma-ray flare around the arrival time of IC-170922A (IceCube Collaboration et al. 2018a).The neutrino emission is produced by the photopion process between protons and BLR radiation.The number of predicted muon neutrino events detected by IceCube during half a year is approximate 0.07 based on the best-fit parameters, which is an acceptable value.
The SED fitting result is shown in Figure 3 and the derived parameters are listed in Table 3.The maximum Lorentz factor of electrons and protons are 7.2 × 10 7 and 6.2 × 10 8 respectively.The kinetic luminosity of the magnetic field, electrons and protons in the blob are 1.09 × 10 44 erg/s, 3.56 × 10 44 erg/s and 2.61 × 10 46 erg/s respectively.The total luminosity is close but still below the Eddington luminosity of the SMBH.

Epoch III: 2021-2022
The SED fitting result is shown in Figure 4 and the derived parameters are listed in Table 3.The magnetic field of the blob is approximately 20.8 G in this case, the distance between the black hole and the blob (flaring zone) is approximately 0.0012 pc.This value is 4 times the distance of the flaring blob in Epoch I.The predicted number of neutrino events, if assuming the muon neutrino effective area of IceCube, is 0.73 for half a year.Note that this neutrino event, GVD-210418CA, is a cascade event detected by Baikal-GVD.The effective area of Baikal-GVD for cascade event is several times smaller than that of IceCube for the track channel (Baikal-GVD Collaboration et al. 2022), so we expect the detection rate of neutrino in this epoch by Baikal-GVD is at the order of ∼ 0.1, which is not unreasonable for the detection.From the perspective of SED fitting, as shown in Figure 4 and Figure 2, we can see that the observed X-ray flux in Epoch III is lower than the X-ray flux upper limit in Epoch I.This leads to a restraint on neutrino flux in Epoch III than that in Epoch I, because otherwise the observed X-ray flux would be overshot by the cascade emission initiated by proton-photon interactions.As shown in Table 3, the total kinetic luminosity is 8.20 × 10 45 erg/s.This value is about three times smaller than that in Epoch I. Also, given the larger distance, the energy density of the X-ray corona has decreased by one order of magnitude.Therefore the photopion production efficiency decreases, as can be seen in Figure 1 by comparing the ratio of the cooling timescale of photopion to the dynamical timescale in this two Epochs.
On the other hand, where the X-ray photons emitted by the corona is still very intense at this distance and absorb the γ-ray photon above GeV emitted by the blob efficiently (Kamraj et al. 2018;Righi et al. 2019).As a result, the blob only has a comparable contribution to the gamma-ray flux around 0.1 GeV, leading to a slightly soft gamma-ray spectrum.The blazar thus appears in a low state in the GeV band during this epoch, which is consistent with the observation.The SED in Epoch IV is similar to that in Epoch III.The X-ray flux in this epoch is slightly lower than that in Epoch III while the gamma-ray flux around 0.1 GeV is a little bit higher.As a result, the best-fit distance of the blob is   12) the total luminosity of electrons and protons of the blob; Col.( 13) the numbers of neutrino events (during half a year).We assumed protons are injected in the blob with a power-law distribution.The spectral index and the minimum Lorentz factor of the proton distribution are fixed to 2 and 1 respectively, while the maximum Lorentz factor of the proton is obtained by equating the acceleration and the cooling or dynamical timescales.The maximum injection luminosity of protons is obtained when the pair cascade emission saturates the X-ray or γ-ray data. .a: calculated with IceCube's effective area for the track channel, see Section 4.3 for more discussion.
a few times farther away from the SMBH (0.003 pc from the SMBH), resulting in a slightly lower photopion efficiency and gamma-ray opacity.Comparing with the Epoch III, the blob's contribution at 0.1 GeV is a little bit higher but the emission in the Fermi-LAT band is still dominated by the background component.The lower X-ray flux leads to a lower neutrino flux, yielding 0.41 muon neutrino events in half a year.The kinetic luminosity of the magnetic field, electrons and protons in the blob are 1.09 × 10 44 erg/s, 1.82 × 10 43 erg/s and 1.12 × 10 46 erg/s respectively.The SED fitting result is shown in Figure 5 and the derived parameters are listed in Table 3.

DISCUSSION AND SUMMARY
In this work, we applied the stochastic dissipation model to simultaneously explain neutrino events and broadband SED of TXS 0506+056 detected in four different epochs, corresponding to the 2014-2015 neutrino outburst and the arrival time of three single muon neutrino event IC-170922A, GVD-210418CA, and IC-220918A.In this model, the emission of the jet is composed of two emission component (Wang et al. 2022).One emission component is the superposition of emission from numerous but comparatively weak dissipation zones along the entire jet, which is considered to be the same in all four epochs.Another component arises from a compact flaring zone (i.e., a blob with intense particle injection) that could form at any position in the jet randomly.In addition to radiation of electrons considered by Wang et al. (2022), we also assume protons are accelerated in the flaring zone and undergo synchrotron radiation and proton-photon interactions (including the Bethe-Heitler process and the photopion process), as well as the emission of pairs generated in the EM cascade initiated by these processes.Some of properties of the flaring zones in different epochs such as the magnetic field and the Doppler factor are related to each other, while the injection luminosity of relativistic electrons and protons and the distance from the SMBH are independent.In particular, the distance of the flaring zone determines the EM environment it is located and leads to different radiation properties.
The best-fit distance of the flaring zones in all four epochs are small, where the one during the Epoch I (2014-2015), Epoch III (2021-2022), Epoch IV (2022-2023) are close to the X-ray corona and the one during Epoch II (2017-2018) is beyond the X-ray corona but still within the BLR.These external radiation fields provide promising target photons for neutrino production.Compared to conventional one-zone models, this model can yield a higher neutrino flux with a more reasonable detection probability of these neutrino events, i.e., 8.2, 0.07, 0.73 (∼ 0.1 for Baikal-GVD), and 0.41 in half a year for Epoch I-IV respectively, while the required kinetic luminosity of the jet is below the Eddington luminosity.Then, we may expect an overall probability of detecting the three single-neutrino event in Epoch II, III, and IV, to be 0.07 × 0.1 × 0.41 ≈ 0.003, which is marginally within 3σ confidence level and not too small to be unreasonable.
We note that the neutrino flux is positively correlated to the X-ray flux of this blazar, which is similar to that predicted in conventional one-zone models.This is because the co-produced EM particles in the photopion production can initiate EM cascades and deposit energy in lower energy band via synchrotron radiation of electron/positron pairs generated in the cascade, while the X-ray flux radiated by primary electrons accelerated in the jet is weak for this blazar.In fact, the main observational constraints on the neutrino flux is from the measured X-ray flux.This also explains why the predicted neutrino detection rate in this model is much higher with respect to that in conventional one-zone models: the low-energy synchrotron bump in the SED is not fully contributed by the neutrino emission zone, so a weaker magnetic field can be employed in the neutrino emission zone which also reduces the cascade emission at the X-ray band.This is similar to reason of a higher neutrino flux in the two-zone model as discussed by Xue et al. (2019).
To conclude, the blazar TXS 0506+056 is a candidate of high-energy neutrino emitter.In addition to the 2014-2015 neutrino outburst and IC-170922A as reported by IceCube several years ago, two additional high-energy muon neutrino events, GVD-210418CA and IC-220918A, are reported recently in spatial association with the blazar (Baikal-GVD Collaboration et al. 2022;Becker Tjus et al. 2022b).If the neutrino events in all four epochs truly originate from TXS 0506+056, the conventional one-zone lepto-hadronic model, with which IC-170922A was explained as a lucky detection (Keivani et al. 2018;Cerruti et al. 2019;Gao et al. 2019), can be ruled out, because the probability of detecting all these neutrino events from the blazar under the one-zone model would be unreasonably low.We attempt to extend the stochastic dissipation model by including hadronic emissions and apply the model to this blazar.We found that SEDs of TXS 0506+056 in all four epochs with neutrino detection can be simultaneously explained and the predicted neutrino event are reasonable.The model supports blazars to be efficient high-energy neutrino emitters and sheds some light on the radiation models of blazars.
.0 Note-The column information are as follows: (1) the year of the observation of neutrino events; Col.(2) the best-fit photon index for a power-law model; Col.(3) the integrated γ-ray flux in the energy range of 0.1-300 GeV; Col.(4) the maximum likelihood test statistic.
column information are as follows: (1) the year of the observation of neutrino events; Col.(2) the exposure time of X-ray; Col.(3) the neutral hydrogen column density of the extragalactic foreground absorption; Col.(4) the best-fit photon index for a power-law model; Col.(5) the integrated X-ray flux in the energy range of 0.3-10 KeV; Col.(6) the reduced χ 2 .
or a power-law distribution, i.e., Q e(p) (γ e(p) ) = Q e(p),0 γ −n e(p),1 e(p) 1 + γ e(p) γ e(p),b (ne(p),2−ne(p),1) Q e(p) (γ e(p) ) = Q e(p),0 γ −n e(p) e(p) , (5) Here γ e(p)min and γ e(p)max are the minimum and maximum Lorentz factor of electrons or protons respectively, γ e(p),b is the break Lorentz factor of electrons or protons, Q e(p),0 is the normalization, n e(p),1 and n e(p),2 represent the lowenergy slope and the high-energy slope of the broken power-law distribution respectively.Q e(p),0 can be obtained from Q e(p) γ e(p) m e(p) c 2 dγ e(p) = L e(p),inj / 4/3πR 3 where c is the speed of light and m e(p) is the rest mass of electron or proton, L e(p),inj is the injection luminosity of electrons or protons.The steady-state density distribution of electrons(protons) can be written as N e(p) (γ e(p) ) = Q e(p) (γ e(p) )t e(p) ,

Figure 1 .
Figure1.The cooling timescale of protons and electrons in different epochs, as well as the opacity of the intrinsic γγ absorption of the radiation zone.The horizontal axis eV represents ernergy of protons for synchrotron, photopion and Bethe-Heitler process in the comoving frames, and also represents ernergy of electrons for synchrotron, SSC and EC process in the comoving frames.In addition, we also plot the optical depth of photons(the red dot-dashed line) in the figure.The vertical grey dashed lines show the maximum electron and proton energy allowed by the tacc =min t cool , t e(p),dyn .Parameters of the model are described in Table3, respectively.

Figure 2 .
Figure 2. Fitting to the broadband SED of TXS 0506+056 in 2014-2015.The red data are from Rodrigues et al. (2019), the gray data points are historical archival data.The pink bow tie shows the muon-neutrino flux as measured during the 14-15 flare (IceCube Collaboration et al. 2018b).

Figure 5 .
Photopion BH total column information are as follows: (1) the year of the observation of neutrino events; Col.(2) the transverse radius of the jet at 0.1 pc; Col.(3) the Doppler factor at 0.1 pc; Col.(4) the magnetic field strength of the dissipation zone for r = 0.1 pc; Col.(5) low energy spectral index of electron; Col.(6) the distance between black hole and blob (flaring zone); Col.(7) high energy spectral index of electrons; Col.(8) the break Lorentz factor of electrons; Col.(9) the injection luminosity of electrons; Col.(10) the magnetic field strength of the blob; Col.(11) the luminosity of magnetic field of the blob; Col.(

Table 1 .
Results of the Analysis of Fermi-LAT Observations

Table 2 .
Results of the Analysis of Swift-XRT Observations of TXS 0506+056

Table 3 .
The physical parameters derived from the unified models fits