The origin of the astrophysical neutrino flux discovered by IceCube remains largely unknown. Several individual neutrino source candidates were observed. Among them is the gammaray flaring blazar TXS 0506056. A similar coincidence of a highenergy neutrino and a gammaray flare was found in blazar PKS 0735178. By modeling the spectral energy distributions of PKS 0735178 we expect to investigate the physical conditions for neutrino production during different stages of the source activity. We analyze the multiwavelength data during the selected periods of time. Using numerical simulations of radiation processes in the source we study the parameter space of onezone leptonic and leptohadronic models and find the bestfit solutions that explain the observed photon fluxes. We show the impact of model parameter degeneracy on the prediction of the neutrino spectra. We show that the available mutliwavelength data are not sufficient to predict the neutrino spectrum unambiguously. Still under the condition of maximal neutrino flux we propose a scenario in which 0.2 neutrino events are produced during the 50 days flare.
['Astrophysics - High Energy Astrophysical Phenomena']
Multiepoch leptohadronic modeling of neutrino source candidate blazar PKS 0735178
{'Multi-epoch leptohadronic modeling of neutrino source candidate blazar PKS 0735+178': '1, ⋆ 2 , 1 3 , 4 , 1 1 5 3 , 6 3  \n- and J. Tammi  \nA. Omeliukh , S. Garrappa , V. Fallah Ramazani , A. Franckowiak , W. Winter , E. Lindfors , K. Nilsson , J. Jormanainen 3 , 6 , 4 , F. Wierda 6 , A. V. Filippenko 7 , W. Zheng 7 , M. Tornikoski 4 , A. Lähteenmäki 4 , 8 , S. Kankkunen 4 , 8 , 4  \n- 1 Ruhr University Bochum, Faculty of Physics and Astronomy, Astronomical Institute (AIRUB), Universitätsstraße 150, 44801 Bochum, Germany\n- 2 Weizmann Institute of Science, Department of Particle Physics and Astrophysics, Herzl Street 234, 76100 Rehovot, Israel\n- 3 Finnish Centre of Astronomy with ESO (FINCA), University of Turku, FI-20014 Turku, Finland\n- 4 Aalto University Metsähovi Radio Observatory, Metsähovintie 114, 02540 Kylmälä, Finland\n- 5 Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, D-15738 Zeuthen, Germany\n- 6 Department of Physics and Astronomy, University of Turku, FI-20014, Turku, Finland\n- 7 Department of Astronomy, University of California, Berkeley, CA 94720-3411, USA\n- 8 Aalto University Department of Electronics and Nanoengineering, P.O. BOX 15500, FI-00076 AALTO, Finland  \nReceived XX; accepted XX', 'ABSTRACT': 'Context. The origin of the astrophysical neutrino flux discovered by IceCube remains largely unknown. Several individual neutrino source candidates were observed. Among them is the gamma-ray flaring blazar TXS 0506 + 056. A similar coincidence of a highenergy neutrino and a gamma-ray flare was found in blazar PKS 0735 + 178.  \nAims. By modeling the spectral energy distributions of PKS 0735 + 178, we expect to investigate the physical conditions for neutrino production during di ff erent stages of the source activity.  \nMethods. We analyze the multi-wavelength data during the selected periods of time. Using numerical simulations of radiation processes in the source, we study the parameter space of one-zone leptonic and leptohadronic models and find the best-fit solutions that explain the observed photon fluxes.  \nResults. We show the impact of model parameter degeneracy on the prediction of the neutrino spectra. We show that the available mutli-wavelength data are not su ffi cient to predict the neutrino spectrum unambiguously. Still, under the condition of maximal neutrino flux, we propose a scenario in which 0.2 neutrino events are produced during the 50 days flare.  \nKey words. Galaxies: BL Lacertae objects: individual - Neutrinos - Methods: numerical - Radiation mechanisms: nonthermal', '1. Introduction': "The detection of a di ff use high-energy astrophysical neutrino flux (IceCube Collaboration 2013) started a new era of neutrino astronomy. However, the nature of the extragalactic neutrino flux remains unclear.  \nBlazars are a subclass of active galactic nuclei (AGNs) with a relativistic jet pointing close to the observer's line of sight. The relativistic boosting of the jet radiation and the overall high power of the jet emission make them prominent neutrino source candidates (see, e.g., review by Giommi & Padovani 2021). The multi-wavelength spectral energy distribution (SED) of blazars typically shows a two-bump structure. While the low-energy bump likely originates from synchrotron radiation of relativistic electrons in the blazar jets, the origin of the high-energy bump is still debated, as contributions from both leptonic and hadronic processes are possible.  \nBlazars can be probed as high-energy neutrino emitters in searches for temporal and spatial associations of blazar flares with high-energy neutrino events. The most promising neutrino blazar candidate is the blazar TXS 0506 + 056, which was in a gamma-ray flaring state during the arrival of a 300 TeV neu-  \ntrino detected by IceCube, resulting in a significance at the 3 σ level (IceCube Collaboration et al. 2018a). An additional flare of TeV neutrinos at the 3 . 5 σ level was found from the direction of TXS 0506 + 056 during a gamma-ray quiet state (IceCube Collaboration et al. 2018b).  \nIn addition to TXS 0506 + 056, IceCube has also detected high-energy events in spatial coincidence with other individual blazars of di ff erent classes, among which are PKS 1424-418 (Kadler et al. 2016; Gao et al. 2017), GB6 J1040 + 0617 (Garrappa et al. 2019), PKS 1502 + 106 (Franckowiak et al. 2020; Rodrigues et al. 2021), and PKS 0735 + 178 (Sahakyan et al. 2022). No statistical correlation between gamma-ray blazars and neutrinos was found in the IceCube analysis (Aartsen et al. 2017; Abbasi et al. 2023). The correlations of neutrino events with blazar positions were claimed by Buson et al. (2022, 2023) but disfavored by Bellenghi et al. (2023). Additionally, a statistical correlation between neutrinos and radio blazars was claimed by Plavin et al. (2021, 2023) but was not confirmed by Abbasi et al. (2023).  \nIn early December 2021, multiple neutrino events detected by IceCube (IceCube Collaboration 2021), Baikal-GVD (Dzhilkibaev et al. 2021), the Baksan Underground Scintillation Telescope (Petkov et al. 2021), and KM3NeT (Filippini et al.  \n2022) were in temporal and spatial coincidence with the largest ever observed flare of the blazar PKS 0735 + 178 in the gammaray band. The source was also flaring in the optical, infrared, ultraviolet, and X-ray bands. Among neutrino source candidates, PKS 0735 + 178 is the only source for which multiple neutrino events from di ff erent detectors were observed.  \nBlazar PKS 0735 + 17 was among the first sources to be designated as 'classical BL Lac' (Carswell et al. 1974), owing to the fact that its spectrum does not display strong, broad emission lines. Therefore, the determination of the distance of the host galaxy is a challenge. Falomo et al. (2021), assuming that the host galaxy belongs to a small group, estimated the redshift of PKS0735 + 178 as z ≈ 0 . 65. Nilsson et al. (2012) performed deep optical imaging and derived z = 0 . 45 ± 0 . 06 in agreement with the absorption redshift of 0.424 (Carswell et al. 1974) within the uncertainties. We adopt z = 0 . 45 as the redshift of PKS 0735 + 178 throughout this work.  \nThe possible neutrino emission from PKS 0735 + 178 was already discussed by Sahakyan et al. (2022), Acharyya et al. (2023), and Prince et al. (2023). The source emission in these works was modeled either during the flare in December 2021 that lasted over one month or based on all available data nonsimultaneous for this source. Multi-epoch leptonic modeling was performed by Bharathan et al. (2024) but does not characterize the source behavior in a multi-messenger context. In this work, we describe the multi-wavelength and multi-messenger emission from PKS 0735 + 178 during both quiescent and active states, which is motivated by the fact that the blazar activity is highly variable. This allows us to consider a scenario in which physical conditions for the neutrino production vary, not requiring sustaining the same extreme conditions (e.g., high proton injection power) over months of the enhanced source activity.  \nThe previous works that utilize highly advanced numerical modeling methods (e.g., Böttcher et al. 2013; Cerruti et al. 2015; Petropoulou et al. 2020a; Cerruti et al. 2019; Oikonomou et al. 2019, in addition to works already mentioned above) usually propose only a single set of radiation model parameters to explain the blazar emission. Radiation models typically have a high number of free parameters, with the possible degeneracy of some of them. In the current study we, for the first time, aim to address this issue and explore the parameter space of one-zone leptonic and leptohadronic models in a more complete way than current state-of-the-art models. We apply a novel approach to search for the best-fit solutions.  \nThe paper is structured as follows. Section 2 presents the available multi-wavelength and neutrino data used for this study. In Section 3 we describe the radiative models that were used to explain the multi-messenger emission as well as our method for searching the best-fit values. Section 4 provides the results of the modeling, which are further discussed and summarized in Section 5.  \nFor the calculations in the paper, we adopt a flat Λ CDMcosmological model with parameters H0 = 70 km s -1 and Ω m = 0 . 3.", '2. Data collection and data analysis': 'The detection of multiple neutrinos from the direction of PKS 0735 + 178 triggered a campaign of multi-wavelength observations of this source. We analyze the available data and build multi-wavelength light curves as well as SEDs for this source.', '2.1. Neutrino emission': 'On 8 December 2021 at 20:02:51.1 UTC, the IceCube real-time alert system (Aartsen et al. 2017) detected a 172 TeV track-like event, IceCube-211208A, with a 30% probability to be an astrophysical neutrino (IceCube Collaboration 2021). Four hours after the IceCube event, Baikal-GVD reported the detection of a 43 TeV cascade event with a probability of ∼ 50%of being an astrophysical neutrino (Dzhilkibaev et al. 2021). Additionally, four days prior to the IceCube event, the Baksan Underground Scintillation Telescope detected a muon neutrino with energy larger than 1 GeV (Petkov et al. 2021) consistent with the location of PKS 0735 + 178 and the directions of the neutrinos detected by IceCube and Baikal-GVD.  \n·  \nFig. 1: Localization of neutrino arrival directions and the position of PKS 0735 + 178. The gray markers indicate the position of Fermi -LAT sources in the 4FGL.  \n<!-- image -->  \nFigure 1 shows the spatial uncertainty contours of the reconstructed neutrino detections in all three neutrino telescopes. Blazar PKS 0735 + 178 lies just outside of the IceCube-211208A 90% uncertainty region. The neutrino still can be associated with the blazar owing to the possible existence of IceCube systematic uncertainties and the fact that 10% of the real counterparts are expected to be outside of it as argued by Sahakyan et al. (2022).  \nThe KM3NeT neutrino detector found no neutrinos from the direction of PKS 0735 + 178 within a time window of one day before and after the IceCube event, but an additional search over the entire month of December 2021 revealed an 18 TeV neutrino event in ARCA, a component designed for the detection of highenergy neutrinos, with a weak association ( p -value 0.14) and no neutrinos in ORCA, a component targeted to low-energy neutrinos (Filippini et al. 2022).', '2.2. Broad-band light curves': 'The multi-wavelength light curves of PKS 0735 + 178 from 2006 to 2022 are illustrated in Fig. 2. In December 2021, the source showed the highest activity in gamma rays over the time of its monitoring with Fermi -LAT. Simultaneously with the gammaray flare, the source was flaring in optical and X-ray bands.', '2.3. Spectral energy distributions': 'Based on the source activity in the gamma-ray band (see Section 2.3.1), we define the following states that will be used for further studies: quiescent state, beginning of the flare (neutrino arrival),  \nTable 1: Selected time periods used for the SED generation.  \nflare peak, and end of the flare. The dates of the selected time periods are shown in Table 1. The corresponding SEDs were composed based on the observations made within those periods.', '2.3.1. Gamma rays': 'We use the gamma-ray observations of PKS 0735 + 178 from the Large Area Telescope (LAT) on board of the Fermi satellite collected between August 2008 and May 2022 presented by Garrappa et al. (2024). The bottom panel of Figure 2 shows the adaptively binned light curve of the gamma-ray flux integrated between 100 MeV and 1 TeV. The continuous all-sky coverage of Fermi -LAT shows a detailed picture of the temporal activity of the source in the 14 years, with the flaring activity temporally coincident with the arrival of IC211208A being by far the brightest activity observed in gamma rays. During this period, from the Bayesian Blocks (Scargle et al. 2013) representation of the light curve in Garrappa et al. (2024), the shortest significant ( ≥ 2 σ ) variability timescale can be derived as the minimum of the quantity (e.g., Meyer et al. 2019)  \nt var , i j = Fi + Fj 2 GLYPH<13> GLYPH<13> GLYPH<13> GLYPH<13> GLYPH<13> GLYPH<13> ti -t j Fi -Fj GLYPH<13> GLYPH<13> GLYPH<13> GLYPH<13> GLYPH<13> GLYPH<13> , (1)  \nwhere Fi and Fj are the fluxes in the bins at the times before ( ti ) and after ( t j ) each block edge that denotes a significant variation. We find the shortest variability timescale at ∼ 14 . 4 days. This timescale will be considered to constrain the radius of the emission region in Sec. 3.  \nIn addition, we use the gamma-ray SEDs from Garrappa et al. 2024 calculated in di ff erent time windows selected from the Fermi -LAT light curve. Three periods have been chosen with respect to the bright gamma-ray flaring activity of the source: the first time window is simultaneous with the IceCube neutrino arrival (MJD 59556-59559), the second one includes the peak of the gamma-ray flare (MJD 59565-59567), and the third one is during part of the decaying phase of the flare (MJD 5957359600). As comparison to the period of flaring activity, we consider also the quiescent state SED from MJD 55219-55233 (see also Table 1).  \nVERITAS and H.E.S.S. observations shortly after the IceCube neutrino detection put upper limits on the very-high-energy gamma-ray emission from PKS 0735 + 178 (Acharyya et al. 2023).', '2.3.2. X-rays': 'The X-ray Telescope (XRT; Burrows et al. 2004) onboard the Neil Gehrels Swift Observatory (Swift) observed the source 24 times between 20 December 2009 and 6 January 2022 in photoncounting (PC) mode. The multi-epoch event lists were down-  \nloaded from the publicly available SWIFTXRLOG ( Swift -XRT Instrument Log) 1 . Following the standard Swift -XRT analysis procedure described by Evans et al. (2009), the data were processed using the configuration described by Fallah Ramazani et al. (2017) for blazars and assuming a Galactic column density of hydrogen of 4 . 42 × 10 20 cm -2 reported by Willingale et al. (2013). The spectra of each observation were fitted using the maximum-likelihood-based statistic for Poisson data (Cash statistics; Cash 1979). The X-ray integral fluxes between 0.3 and 10 keV of these observations are shown in Fig. 2.  \nNo spectral variability was observed within the Swift -XRT data between 31 January 2010 and 17 February 2010 (observation ID 00090099001 and 00090099002). These two observations were combined to produce an average X-ray spectrum of the source during the quiescent state. Similarly, no spectral variability was observed between 25 December 2021 and 6 January 2022 (observation ID 00036372023, 00036372024, and 00036372025). These three observations were combined to produce an average X-ray spectrum of the source after the flaring state. In order to produce the X-ray spectra for the periods of the neutrino arrival and the gamma-ray peak (as defined above), the Swift -XRT data on 10 December 2021 (observation ID 00036372014) and 17 December 2021 (observation ID 00036372021) were used.  \nThe X-ray spectrum above 2 keV from NuSTAR observation on December 11, 2021 was taken from Acharyya et al. (2023).', '2.3.3. Ultraviolet and Optical': 'Simultaneously with Swift -XRT observations, PKS 0735 + 178 was observed with the UVOT instrument (Roming et al. 2005) onboard Swift . The source was observed in the W 1 band on 31 January 2010, and in the U , B , W 1, M 2, and W 2 bands on 10, 17, 25, 30 December 2021 and 6 January 2022. In addition, it was observed with the V filter on 6 January 2022. The data were reduced using the analysis pipeline from the Space Science Data Center (SSDC 2 ).  \nPKS 0735 + 178 has been monitored in the optical R band as part of the Tuorla blazar monitoring program 3 since 2005. The monitoring observations were performed using a 35 cm Celestron telescope coupled to the KVA (Kunglinga Vetenskapsakademi) telescope located on the Canary Island of La Palma, Spain. Data analysis was performed using standard procedures with the semi-automatic pipeline developed in Tuorla (Nilsson et al. 2018).  \nThe optical light curve is complemented by data obtained with the 76 cm Katzman Automatic Imaging Telescope (KAIT) as part of the Lick Observatory Supernova Search (LOSS; Filippenko et al. 2001) in the Clear (open) band (close to the R band; see Li et al. 2003). All images were reduced using a custom pipeline 4 detailed by Stahl et al. (2019). Several nearby stars were chosen from the Pan-STARRS1 5 catalog for calibration; their magnitudes were transformed into Landolt magnitudes using the empirical prescription presented by Equation 6 of Tonry et al. (2012). The light curve is shown in Figure 2.  \nFor the quiescent state SED, we used the average R -band flux on MJD 55231, 55232, and 55233 where no significant variabil-  \nFig. 2: The multi-wavelength light curves of PKS 0735 + 178 from 2006 to 2022. The gray vertical line corresponds to the time of the IceCube neutrino detection. The light curves show that the source was in the exceptional high state in all bands during the neutrino arrival.  \n<!-- image -->  \nity was detected in the data obtained by the Tuorla blazer monitoring. During the neutrino arrival time, the optical measurement by KAIT on MJD 59559 was used for the SED modeling. For the SEDmodeling of the gamma-ray flare peak, the optical and nearinfrared (NIR) data are obtained from Lindfors et al. (2021). Finally, for the post-flare period, the NIR data from Carrasco et al. (2021) and optical data from KAIT on MJD 59585 and 59600 were used.  \nThe observed UV, optical, and NIR fluxes are corrected for Galactic extinction using the values from Schlafly & Finkbeiner (2011). The contribution of the host-galaxy flux in UV, optical, and NIR bands is negligible when compared with the statistical uncertainty of the flux measurements (Nilsson et al. 2012).', '2.3.4. Radio': 'The Metsähovi radio telescope (located at Kirkkonummi, Finland), operating at 37 GHz, has been observing the source for three decades. We selected radio data obtained after mid-2006 for our study. The instrument and data-reduction procedures are described by Teräsranta et al. (1998). The long-term radio light curve is presented in Figure 2. Additionally, we use MOJAVE (Lister et al. 2018) measurements of the 15 GHz radio flux from two days after the arrival of IceCube-211208A.', '3.1. Simulation framework of one-zone models': 'Within the one-zone framework, we assume that the observed photon and neutrino fluxes originate from the emission of highly energetic particles inside a compact zone (blob) in the jet. Radio fluxes cannot be explained within the one-zone model framework since the compact region is optically thick to lowfrequency radio emission owing to synchrotron self-absorption. In addition, the typical size of the emission region obtained from the gamma-ray variability timescales is two to four orders of magnitude smaller than the size of the radio core which is the dominant emitter of the radio flux. For our modeling, radio fluxes will be treated as upper limits.  \nWe utilize the open-source code AM 3 (Klinger et al. 2023) for simulating the radiation processes and interactions undergone by relativistic electrons and protons. AM 3 numerically solves the system of coupled di ff erential equations that describes the evolution of the particle spectra, in a fully timedependent and self-consistent manner. The absorption of highenergy gamma rays due to extragalactic background light (EBL) is accounted for in all models, based on Domínguez et al. (2011).', '3.2. Search for the best solutions': 'Self-consistent modeling of the radiation processes in one-zone models utilizes a high number of free parameters (usually between seven and fourteen). Some of the parameters can be ten-  \ntatively constrained from the observables, such as an upper limit on the size of the emission region (from the gamma-ray variability timescale) or the blob Lorentz factor (from radio observations). The radiation power of the jet is fundamentally limited by the Eddington luminosity. Still, in most cases, those constraints leave large regions of possible parameter values and the problem of finding the best-fit values arises.  \nThe radiation processes are highly sensitive to the values of the model parameters. A small change in a parameter value can result in an SED whose fluxes di ff er by orders of magnitude or has a drastically di ff erent shape due to a di ff erent radiation process starting to contribute significantly. Thus, the goodnessof-fit is not a smooth function of the model parameters, which introduces challenges to numerical di ff erentiation. Additionally, common algorithms based on gradient descent methods encounter performance problems in high-dimensional parameter spaces (Han & Neumann 2006).  \nWe propose here a simple and e ff ective two-step algorithm for searching for the best-fit solutions. As a first step, we generate a large number of di ff erent SEDs with values of the model parameters that are equally spaced between their boundary values. When choosing an equal number of points for each parameter, we get N m di ff erent models, where m is the number of model parameters (dependent on the model) and N is the number of points per dimension (dependent on the available computational resources). This gives a rough localization of the best solutions. As a second step, we locally minimize the goodness-of-fit value within a narrow region in the parameter space. Since we know that all neighbor points of the current best-fit solution have a lower goodness-of-fit value, we adopt their parameter values as new boundaries for searching of the new best fit. We apply the migrad and symplex algorithms from Minuit (James & Roos 1975; Dembinski et al. 2020) for this local minimization because the parameter space is now better constrained, ensuring the convergence of the algorithms.  \nA search of the parameters with algorithms like Markov Chain Monte Carlo (MCMC) or genetic algorithm typically requires simulating 10 6 -10 7 SEDs for finding one best-fit solution. The proposed grid scan requires only a single generation of a set of SEDs (in the source frame) which can be reused for different sources, thus using the computational resources only to recalculate the goodness-of-fit. The local minimization is computationally e ff ective as well, requiring only ∼ 200-300 SED generations. The disadvantage of the method is an exponential increase in the required computational resources when the number of model parameters increases.', '3.3. Leptonic models': "We start with the simplest models, which assume that all radiation from the NIR to gamma rays originates solely from leptonic processes.  \nWe assume that electrons are accelerated to a simple powerlaw spectrum 6 dN / d γ ' e ∝ γ '-α e e with spectral index α e, spanning a range of Lorentz factors from γ ' min e to γ ' max e . The energy spectrum of the electrons is normalized to the total electron luminosity parameter, L ' e . These particles are then injected into a single spherical blob of size R ' b (in the comoving frame of the jet) moving along the jet with Lorentz factor Γ , where they encounter a homogeneous and isotropic magnetic field of strength B ' . We  \nassume the jet is observed at an angle θ obs = 1 / Γ b relative to its axis, resulting in a Doppler factor of δ D = Γ b.  \nTo explore the seven-dimensional parameter space of the leptonic models, we do a simple grid scan probing 10 points per dimension resulting in 10 million simulated models. The parameters of the leptonic models and the boundaries of the parameter space are given in Table A.1.  \nWe find the best-fit parameters in two steps. First, for each of the four SEDs (each corresponding to a di ff erent period of blazar activity) we evaluate the reduced χ 2 using the 10 million simulated models. It is calculated as  \nχ 2 = 1 N -N par + 1 X i ( F data i -F model i ) 2 σ 2 i , (2)  \nwhere N is the number of data points, N par the number of free parameters in the model, F data are the observed fluxes, F model are the predicted fluxes by the model, and i is the summation index which corresponds to the observed frequency values. The models that predict higher flux values than Fermi -LAT sensitivity-based upper limits are rejected because the model fluxes are expected to significantly decrease in the corresponding energy range.  \nFor the blazar states during the 2021 flare, an additional constraint on the radius of the emission region applies:  \nR ' b ≤ δ D c t var 1 + z . (3)  \nWe use the value of t var = 14.4 days for the 2021 flare from Sect. 2.3.1, and the values of R ' b and δ D are free model parameters.  \nWe select the best-fitting model for a particular SED based on the value of its reduced χ 2 . Afterward, we locally minimize the reduced χ 2 with the parameter boundaries corresponding to the neighbor parameter values on the grid to obtain a final solution. If the values of R ' b and δ D do not satisfy Eq. 3, the model is rejected.", '3.4. Leptohadronic models': "In an alternative scenario, we assume that both electrons and protons are accelerated in the source to power-law spectra dN / d γ ' e , p ∝ γ ' -α e , p e , p with spectral indices α e , p, spanning a range of Lorentz factors from γ ' min e , p to γ ' max e , p . The energy spectra of the electrons and protons are normalized to the corresponding total electron and proton luminosities, L ' e and L ' p . Similarly to the leptonic case, electrons and protons undergo interactions and radiate inside of the spherical blob of size R ' with a homogeneous and isotropic magnetic field of strength B ' moving along the jet with Lorentz factor Γ .  \nThe application of the same two-step algorithm (grid scan followed by a local minimization) described in Section 3.3 would now require significantly more computational resources, since the number of model parameters increased from seven to ten.  \nTo circumvent this issue, we follow the approach adopted by Rodrigues et al. (2024b,a) in which IR, optical, UV, and gammaray fluxes are first fit with a purely leptonic model followed by the selection of hadronic parameters to fit the X-rays. Most of the contribution to optical and GeV gamma-ray emission typically originates from leptons, even in the models where protons are added to the emission zone, as shown by leptohadronic models for several neutrino candidate blazars (e.g., Cerruti et al.  \n2019; Rodrigues et al. 2021; Sahakyan et al. 2022). We fit first IR, optical, UV, and gamma rays via the grid scan method (using the same boundaries as in Table A.1) and locally minimize them using Minuit . X-ray fluxes at this step are treated as upper limits. Afterward, we add protons to the AM 3 simulation. All leptonic parameters are fixed to the previously found values, while the four hadronic parameters are varied according to the grid probing 10 4 combinations. The boundaries of parameters in the hadronic parameter space used for the grid scan are shown in Table A.2. Adding Lp = 0 to the parameter space ensures a possible convergance to the purely leptonic solutions. As a final step, all ten parameters of the leptohadronic models are locally minimized with Minuit allowing for slight variations (with the boundaries being parameter value ± 20%) in each of them.", '4. Results': 'We apply the above described method of the analysis of the parameter space of a one-zone model to the SEDs of PKS 0735 + 178 during the selected periods of observations.', '4.1. Purely leptonic fits': "When analyzing the reduced χ 2 values resulting from the grid scan, we conclude that many di ff erent combinations of parameters yield similar goodness-of-fit values for a particular dataset. We select two physically di ff erent solutions for each SED and explore the parameter space around them to illustrate model parameter degeneracy. The number of solutions that provide a similarly good fit to the observed data can be higher, especially if the boundaries of the parameters are expanded beyond what was selected for this study.  \nThe results of the leptonic modeling are shown in Fig.3. The corresponding parameter values are presented in Table A.3. We show two best-fit results that correspond to distant regions in the parameter space also implying physically di ff erent solutions. For each epoch, we refer to the two di ff erent solutions as 'slow' and 'fast' based on their Lorentz factors. The slow solutions have Lorentz factor values between two and six, while for the fast solution those values span between twelve and eighteen. The different values of the magnetic field, electron luminosities, and Lorentz factors lead to the di ff erent shapes of the SED peaks. Still, both slow and fast solutions explain the observed data well, with the values of reduced χ 2 being 1.4 and 1.1, respectively. In the period of neutrino arrival, the fast solution is characterized by a substantially higher magnetic field (4 G compared to 0.7 G for the slow model), a softer electron energy spectrum, a smaller emission-region size, and a lower electron luminosity. Despite the fact that model parameters lie in the di ff erent regions of the parameter space, both models explain well the observed photon fluxes (the major contribution to the large value of χ 2 / n.d.f. comes from the five UV data points with ∼ 1% error bars). Similar di ff erences between slow and fast solutions are observed for the gamma flare peak SED. For the post-flare SED, on the contrary, our method failed to explain the gamma-ray fluxes. Both slow and fast solutions fit optical, UV, and X-rays. While the slow solution produces a significantly lower gamma-ray flux than observed, the fast solution matches the flux level but fails to reproduce the correct shape.  \nSmall deviations in values of certain model parameters can still keep the model valid. To explore how the change of each pair of parameters influences the goodness of fit (the rest of the parameters are fixed to their best-fit values), we build contour maps of the projections of the parameter space. Fig. 4 shows the  \ncontour maps constraining the parameters of the neutrino arrival SED fast solution. Similar plots for the remaining leptonic models can be found in Appendix B. Since our original grid resolution (ten points per dimension) is insu ffi cient to draw the contours, we perform additional scans in which all parameters but two are fixed to their best-fit values, probing 100 points per dimension. The color bar corresponds to the values of reduced χ 2 with values over 100 being saturated.  \nThe best-fit fast solution for the period of neutrino arrival suggest that the sychrotron emission is produced by a population of electrons with a soft spectrum (with a power-law index of 2.7). Therefore, the spectrum is dominated by low-energy electrons and the high-energy end of the spectrum is not well constrained, which can be seen in the third row of plots in Fig. 4. Since the peak of the synchronton emission cannot be estimated from the SED (which was also noticed by Acharyya et al. 2023), the minimum electron energy also has a wide range of possible values. The Lorentz factor, magnetic field, and electron luminosity, on the contrary, produce good fits only in narrow regions of the parameter space.", '4.2. Leptohadronic fits': "Since the leptonic grid scan revealed the existence of multiple solutions with close χ 2 / n.d.f value, we modeled the leptohadronic SEDs according to the grid-scan-based procedure described in Section 3.4 also for the slow and fast solutions.  \nBefore the local minimization, we selected the models with the lowest reduced χ 2 value for each SED and each solution type (slow or fast). In the leptonic case, the slow and fast solutions had similar values of reduced χ 2 . When adding protons to the radiation zone, the values of reduced χ 2 for the slow leptohadronic solutions became more than twice worse than that for the fast leptohadronic solutions. We also notice that among the best-fit slow leptohadronic solutions, the hadronic contribution was either zero or heavily suppressed owing to favoring the lowest proton energies or the lowest proton luminosities. When searching for the leptohadronic solutions, we started with the assumption that hadrons mostly contribute to the observed X-ray fluxes as described in Sec. 3. Generally, higher values of Lorentz factors produce a deeper gap in the X-rays, which can be filled with the emission from hadronic processes. Therefore, the slow solutions naturally suggest a leptonic origin of the X-rays, with any additional contribution from hadronic processes leading to overshooting the observed X-ray fluxes. The slow leptohadronic models converge to purely leptonic models with the lowest proton energies or luminosities producing a better fit than any other combination of the hadronic parameters. Since we are interested in the leptohadronic models, all the following steps were carried out only with the fast leptohadronic solutions for all selected periods of the blazar activity.  \nAs for the fast leptohadronic solutions, the analysis revealed the degeneracy of the hadronic contributions in all fast leptohadronic models. We find that along with the best-fit solution, there exist many other solutions with the same leptonic parameters and di ff erent hadronic contributions; all these solutions have very similar goodness-of-fit values (less than 2% di ff erence in the reduced χ 2 values).  \nAn example of the degeneracy of the hadronic contributions is shown in Figure 5. Although all shown leptohadronic models produce similar levels of photon fluxes that agree with the data, the hadronic contributions in each model are subdominant and di ff er both in spectral shapes and in energies. This results in a possible range of neutrino spectra located in the shadowed  \nFig. 3: Purely leptonic best fits for the SEDs of PKS 0735 + 178 during the selected periods. Slow and fast solutions are named according to their Lorentz factors and correspond to two models with di ff erent parameters. The parameters are given in Table A.1  \n<!-- image -->  \n.  \nFig. 4: Values of the reduced χ 2 around the best-fit fast solution for the neutrino arrival SED. Plots for other models and SEDs are shown in Appendix B.  \n<!-- image -->  \ngray area in Fig. 5. The parameters of the leptohadronic models shown in Fig. 5 are presented in Table A.4. The degeneracy comes from the fact that the total energy of the protons can be deposited into the emission zone in di ff erent ways. For example, a high number of low-energy protons can produce a similar amount of radiation as a small number of high-energy protons. The neutrino spectra for models 1-9 in Fig. 5 span from lower energies to higher energies as the maximum proton energy for each model increases. Models 2, 3, and 4 all have the same minimal and maximal proton energies and di ff er only by the value of spectral index, which results into the di ff erence in the required proton luminosity being two orders of magnitude between models 2 and 4. Models 1-3 with lower proton energies have high ratio L ' p / L Edd, while models 7 and 8 have L ' p / L Edd ≈ 10 -4 . The hadronic interactions of the proton population in model 5 produce a photon spectrum that peaks around the X-ray gap in the SED, which limits the possible proton luminosity and leads to lower neutrino fluxes than other models.  \nFig. 5: Multiple leptohadronic models that explain the observed photon fluxes during the neutrino arrival. The solid curves correspond to the electromagnetic radiation. The dotted curves represent the contribution of the hadronic process in total photon fluxes, and the dashed curves represent the predicted all-flavor neutrino spectra. The values of the hadronic parameters corresponding to each model number from the legend are shown in Table A.4.  \n<!-- image -->  \nThe collected data for the neutrino arrival SED are located in the typical frequency ranges of available instruments (NIR, optical, UV, X-rays, and gamma rays). They also constrain well the trough in the X-rays. The fact that the observed data are not su ffi cient to estimate the expected neutrino emission from the source is critical for modeling neutrino source candidates. The same photon energy ranges were used for obtaining the data for other sources associated with neutrino emission such as TXS 0506 + 056 or PKS 1502 + 106. The results of our modeling show that the SED alone is not su ffi cient to unambiguously predict the neutrino spectrum.  \nAmong the solutions with various hadronic contributions, we search for the one that predicts the highest number of detected neutrinos in IceCube. We do this by calculating the number of detected neutrino events as  \nN ν = 1 3 T Z Φ ν ( E ) A e ff ( E , θ ) dE , (4)  \nArticle number, page 8 of 15  \nwhere A e ff ( E , θ ) is the e ff ective area of the neutrino detector which depends on the neutrino energy and source declination, Φ ν ( E ) is the all-flavor neutrino flux, and T is the exposure time. A coe ffi cient 1 / 3 is introduced to account for neutrino mixing during the propagation, while we only consider IceCube's muon track channel. We calculate the muon neutrino event rates in IceCube using the experiment's e ff ective area (IceCube Collaboration 2017). For the exposure time, we use the di ff erent flare periods as follows: beginning of the flare (neutrino arrival SED) - 15 days, gamma flare peak - 8 days, and post flare - 27 days.  \nThe expected number of neutrino events is calculated for all models with δ ( χ 2 / n.d.f.) = χ 2 model / n.d.f. -χ 2 best / n.d.f. < 1 corresponding to the various hadronic contributions. One model resulting in the maximum number of neutrino events is selected. Those maximal neutrino models for all analyzed periods of activity of PKS 0735 + 178 are shown in Figure 6. The corresponding parameters of the leptohadronic models are shown in Table A.5.  \nUnder the condition of maximization of the neutrino fluxes in IceCube, we predict at most 0.04 muon neutrino event per year during the quiescent state and 0.2 muon neutrinos during the 50 days 2020-2021 flare.", '5. Discussion': "We now discuss the obtained results and place them in the context of previous studies of this source. The model presented for the neutrino arrival SED in Fig. 6 predicts a slightly higher number of neutrino events than the hybrid model of Sahakyan et al. (2022) which has a similar physical setup. Meanwhile, the proton luminosity in our model is four orders of magnitude lower. This is explained by the fact that the maximum proton energy in our models is not fixed. Protons of higher energy produce more neutrinos through photo-pion production, and the resulting cascades contribute mostly to the gamma-ray emission and thus are not limited by the X-ray emission. This compensates for the fact that the neutrino peak is shifted toward high energies in comparison with the IceCube maximum sensitivity which is around 100 TeV (Aartsen et al. 2019).  \nOur best fit for the gamma-ray flare peak SED, on the contrary, produces almost three times more neutrino events in IceCube than the external field model of Sahakyan et al. (2022). In this case, lower proton energies and higher proton luminosities were preferable, leading to a scenario in which more neutrinos are produced (proportional to the proton luminosity) in the energy range of IceCube's peak sensitivity.  \nThe model for the post-flare SED cannot explain the observed gamma-ray fluxes. The fact that our algorithm failed to find any leptonic fit (and subsequently leptohadronic fit as well) that could simultaneously explain optical and gamma-ray data may be an indication of the presence of external fields. The hard shape of the gamma-ray spectrum, also observed in other sources, may be explained by synchrotron self-Compton (SSC) and external Compton interactions with photons from the broadline region (Fichet de Clairfontaine et al. 2023). In our case of multi-epoch modeling, this result can indicate an evolution of the emission region. The emission zone could be located between the supermassive black hole (SMBH) and the broad-line region (BLR) at the moment of the neutrino detection and during the peak of the gamma-ray flare. With the blob moving along the jet in the post-flare period, the emission zone can enter the BLR. In our models, introducing external photon fields created by the BLR adds an additional parameter (the disk temperature) and increase the complexity of numerical modeling with more  \nFig. 6: Best-fit models for the selected periods of PKS 0735 + 178 activity under the condition of maximization of neutrino events in IceCube.  \n<!-- image -->  \ninteractions. This creates a computational challenge and is left for future work.  \nFor the quiescent-state SED, even under the condition of maximization of the neutrino events in IceCube, the best fit can produce at most 0.04 neutrino events per year, consistent within Poisson fluctuations with 0.0 neutrino events found in the IceCube point-source analysis using 10 yr of IceCube data between 2008 and 2018 (Aartsen et al. 2020).  \nAcharyya et al. (2023) argue that the SSC model cannot explain the SED of PKS 0735 + 178, and instead propose a model with external photon fields which induces higher neutrino rates owing to more targets for p γ interactions. As shown in the results, our SSC model for the neutrino arrival period explains the observed photon fluxes. It produces a factor of two fewer neutrino events than the external field model of Acharyya et al. (2023) but lowers the proton power requirement by two orders of magnitude. The neutrino emission during the gamma flare peak, on the contrary, produces a factor of two more neutrino events than the external field model by Acharyya et al. (2023).  \nThe leptohadronic one-zone models have a long-standing problem of too high proton power requirement to explain the observed neutrino emission (see, e.g., Gasparyan et al. 2022; Petropoulou et al. 2020b). The fundamental jet luminosity is constrained by the Eddington limit L Edd = 1 . 3 × 10 38 ( M / M ⊙ ) erg s -1 (Ghisellini et al. 2010). The mass of the SMBH in PKS 0735 + 178 was estimated based on the optical intraday variability timescales (Gupta et al. 2012) as 1 . 89 × 10 8 M ⊙ (corresponding to L Edd = 2 . 5 × 10 46 erg s -1 ) assuming a jet Doppler factor of δ ≈ 3 . 5. If the Doppler factor of the emission zone is higher, which is a case for our leptohadronic models, the corresponding value of the SMBH mass increases by a factor of five, leading to L Edd ≈ 10 47 erg s -1 . Even though the quiescent and neutrino arrival SED are explained well with sub-Eddington proton luminosities, the gamma flare peak and post-flare SEDs with maximum neutrino rates still require L ' p = 10-100 L Edd.  \nWe analyze the multi-messenger behavior of the source during the di ff erent activity states. For the purely leptonic models, the existance of mutliple solutions does not allow us to trace the exact evolution of the radiation zone since the causal relation between di ff erent solutions at any epoch is not clear. However, we still see that for any solution the flare is characterized by higher values of the magnetic field strength and electron luminosity, simultaneous increase of which suggests energy equipartition.  \nUnder the condition of neutrino rates maximization, the different periods during the flare have enhanced magnetic field strength and more compact emission-zone regions. The neutrino flare could be linked to the increase of minimum proton energy and hardening of the proton spectrum leading to increased p γ interaction rates.", '6. Summary and conclusions': "In this work, we modeled the mutli-epoch emission of PKS 0735 + 178 including di ff erent stages of flare evolution. The grid-scan-based approach showed the degeneracy of leptonic and leptohadronic model parameters. Under the condition of maximization of the neutrino detection rates, we find leptohadronic solutions for di ff erent stages of the flare that produce 0.2 neutrino events during the 50 day flare, a higher number of neutrino events than any of the previously proposed models. We also show that in the quiescent state the maximum neutrino rate, which in our model is limited by the X-ray emission, is two orders of magnitude lower than during the flare. Still, the post-flare SEDcannot be reproduced with a simple one-zone leptohadronic model and likely requires external field photons for the explanation of the gamma rays, which may be a sign of spatial evolution of the emission zone.  \nThe demonstrated degeneracy of the hadronic parameters in the one-zone leptohadronic models creates a significant challenge for understanding potential neutrino sources. The prop-  \nerties of the proton spectrum cannot be constrained by the observed photon fluxes, leading to a great uncertainty in the subsequent neutrino emission. To break this degeneracy, nextgeneration neutrino telescopes such as IceCube-Gen2 (Aartsen et al. 2021), KM3NeT (Adrián-Martínez et al. 2016), or P-ONE (Malecki 2024) are needed. Additionally, multi-wavelegth polarization can potentially constrain the hadronic component (Zhang et al. 2024), thus highlighting the importance of future X-ray and MeV polarimeters such as eXTP (Zhang et al. 2016), COSI (Tomsick et al. 2019), and AMEGO-X (Caputo et al. 2022).  \nAcknowledgements. We thank Xavier Rodrigues for valuable discussions. This research has made use of data from the MOJAVE database that is maintained by the MOJAVE team (Lister et al. 2018). This publication makes use of data obtained at the Metsähovi Radio Observatory, operated by the Aalto University. A.O. was supported by DAAD funding program 57552340. A.F. acknowledges support from the DFG via the Collaborative Research Center SFB1491 Cosmic Interacting Matters - From Source to Signal . A.V.F.'s research group at UCBerkeley acknowledges financial assistance from the Christopher R. Redlich Fund, Gary and Cynthia Bengier, Clark and Sharon Winslow, Alan Eustace (W.Z. is a Bengier-Winslow-Eustace Specialist in Astronomy), and numerous other donors. KAIT and its ongoing operation were made possible by donations from Sun Microsystems, Inc., the Hewlett-Packard Company, AutoScope Corporation, Lick Observatory, the U.S. National Science Foundation, the University of California, the Sylvia & Jim Katzman Foundation, and the TABASGO Foundation. Research at Lick Observatory is partially supported by a generous gift from Google.", 'References': "- Aartsen, M., Ackermann, M., Adams, J., et al. 2017, Astroparticle Physics, 92, 30\n- Aartsen, M. G., Abbasi, R., Ackermann, M., et al. 2021, Journal of Physics G Nuclear Physics, 48, 060501\n- Aartsen, M. G., Abraham, K., Ackermann, M., et al. 2017, ApJ, 835, 45\n- Aartsen, M. G., Ackermann, M., Adams, J., et al. 2020, Phys. Rev. Lett., 124, 051103\n- Aartsen, M. G., Ackermann, M., Adams, J., et al. 2019, European Physical Journal C, 79, 234\n- Abbasi, R., Ackermann, M., Adams, J., et al. 2023, ApJ, 954, 75\n- Acharyya, A., Adams, C. B., Archer, A., et al. 2023, ApJ, 954, 70\n- Adrián-Martínez, S., Ageron, M., Aharonian, F., et al. 2016, Journal of Physics G Nuclear Physics, 43, 084001  \nBellenghi, C., Padovani, P., Resconi, E., & Giommi, P. 2023, ApJ, 955, L32  \nBharathan, A. M., Stalin, C. S., Sahayanathan, S., Bhattacharyya, S., & Mathew,  \nB. 2024, MNRAS, 529, 3503  \n- Böttcher, M., Reimer, A., Sweeney, K., & Prakash, A. 2013, ApJ, 768, 54\n- Burrows, D. N., Hill, J. E., Nousek, J. A., et al. 2004, in Proc. SPIE, Vol. 5165, XRay and Gamma-Ray Instrumentation for Astronomy XIII, ed. K. A. Flanagan & O. H. W. Siegmund, 201-216\n- Buson, S., Tramacere, A., Oswald, L., et al. 2023, arXiv e-prints, arXiv:2305.11263\n- Buson, S., Tramacere, A., Pfei ff er, L., et al. 2022, ApJ, 933, L43\n- Caputo, R., Ajello, M., Kierans, C. A., et al. 2022, Journal of Astronomical Telescopes, Instruments, and Systems, 8, 044003\n- Carrasco, L., Recillas, E., Escobedo, G., et al. 2021, The Astronomer's Telegram, 15148, 1\n- Carswell, R. F., Strittmatter, P. A., Williams, R. E., Kinman, T. D., & Serkowski, K. 1974, ApJ, 190, L101\n- Cash, W. 1979, ApJ, 228, 939\n- Cerruti, M., Zech, A., Boisson, C., et al. 2019, MNRAS, 483, L12\n- Cerruti, M., Zech, A., Boisson, C., & Inoue, S. 2015, Monthly Notices of the Royal Astronomical Society, 448, 910\n- Dembinski, H., Ongmongkolkul, P., Deil, C., et al. 2020\n- Domínguez, A., Primack, J. R., Rosario, D. J., et al. 2011, MNRAS, 410, 2556\n- Dzhilkibaev, Z. A., Suvorova, O., & Baikal-GVD Collaboration. 2021, The As-\n- tronomer's Telegram, 15112, 1  \nEvans, P. A., Beardmore, A. P., Page, K. L., et al. 2009, MNRAS, 397, 1177  \nFallah Ramazani, V., Lindfors, E., & Nilsson, K. 2017, A&A, 608, A68  \nFalomo, R., Treves, A., & Paiano, S. 2021, The Astronomer's Telegram, 15132,  \n1  \n- Fichet de Clairfontaine, G., Buson, S., Pfei ff er, L., et al. 2023, ApJ, 958, L2 Filippenko, A. V., Li, W. D., Tre ff ers, R. R., & Modjaz, M. 2001, in Astronomical Society of the Pacific Conference Series, Vol. 246, IAU Colloq. 183: Small Telescope Astronomy on Global Scales, ed. B. Paczynski, W.-P. Chen, &C. Lemme, 121\n- Filippini, F., Illuminati, G., Heijboer, A., et al. 2022, The Astronomer's Telegram, 15290, 1\n- Franckowiak, A., Garrappa, S., Paliya, V., et al. 2020, The Astrophysical Journal, 893, 162\n- Gao, S., Pohl, M., & Winter, W. 2017, ApJ, 843, 109\n- Garrappa, S., Buson, S., Franckowiak, A., et al. 2019, ApJ, 880, 103\n- Garrappa, S., Buson, S., Sinapius, J., et al. 2024, arXiv e-prints, arXiv:2401.06666\n- Gasparyan, S., Bégué, D., & Sahakyan, N. 2022, Monthly Notices of the Royal Astronomical Society, 509, 2102\n- Ghisellini, G., Tavecchio, F., Foschini, L., et al. 2010, Monthly Notices of the Royal Astronomical Society, 402, 497\n- Giommi, P. & Padovani, P. 2021, Universe, 7, 492\n- Gupta, S., Pandey, U., Singh, K., et al. 2012, New Astronomy, 17, 8\n- Han, L. & Neumann, M. 2006, Optimization Methods and Software, 21, 1 IceCube Collaboration. 2013, Science, 342, 1242856\n- IceCube Collaboration. 2017, The Astrophysical Journal, 835, 151\n- IceCube Collaboration. 2021, GRB Coordinates Network, 31191, 1\n- IceCube Collaboration, Aartsen, M. G., Ackermann, M., et al. 2018a, Science, 361, eaat1378\n- IceCube Collaboration, Aartsen, M. G., Ackermann, M., et al. 2018b, Science, 361, 147\n- James, F. & Roos, M. 1975, Computer Physics Communications, 10, 343 Kadler, M., Krauß, F., Mannheim, K., et al. 2016, Nature Physics, 12, 807\n- Klinger, M., Rudolph, A., Rodrigues, X., et al. 2023, AM 3 : An OpenSource Tool for Time-Dependent Lepto-Hadronic Modeling of Astrophysical Sources\n- Li, W., Filippenko, A. V., Chornock, R., & Jha, S. 2003, Publications of the Astronomical Society of the Pacific, 115, 844\n- Lindfors, E., Hovatta, T., Pursimo, T., et al. 2021, The Astronomer's Telegram, 15136, 1\n- Lister, M. L., Aller, M. F., Aller, H. D., et al. 2018, ApJS, 234, 12 Malecki, P. 2024, Universe, 10\n- Meyer, M., Scargle, J. D., & Blandford, R. D. 2019, ApJ, 877, 39\n- Nilsson, K., Lindfors, E., Takalo, L. O., et al. 2018, A&A, 620, A185\n- Nilsson, K., Pursimo, T., Villforth, C., et al. 2012, A&A, 547, A1\n- Nilsson, K., Pursimo, T., Villforth, C., et al. 2012, Astronomy & Astrophysics, 547, A1\n- Oikonomou, F., Murase, K., Padovani, P., Resconi, E., & Mészáros, P. 2019, MNRAS, 489, 4347\n- Petkov, V. B., Novoseltsev, Y. F., Novoseltseva, R. V., & Baksan Underground Scintillation Telescope Group. 2021, The Astronomer's Telegram, 15143, 1 Petropoulou, M., Murase, K., Santander, M., et al. 2020a, ApJ, 891, 115\n- Petropoulou, M., Oikonomou, F., Mastichiadis, A., et al. 2020b, ApJ, 899, 113\n- Plavin, A. V., Kovalev, Y. Y., Kovalev, Y. A., & Troitsky, S. V. 2021, ApJ, 908, 157\n- Plavin, A. V., Kovalev, Y. Y., Kovalev, Y. A., & Troitsky, S. V. 2023, MNRAS, 523, 1799\n- Prince, R., Das, S., Gupta, N., Majumdar, P., & Czerny, B. 2023, Monthly Notices of the Royal Astronomical Society, 527, 8746\n- Rodrigues, X., Garrappa, S., Gao, S., et al. 2021, The Astrophysical Journal, 912, 54\n- Rodrigues, X., Karl, M., Padovani, P., et al. 2024a, arXiv e-prints, arXiv:2406.06667\n- Rodrigues, X., Paliya, V. S., Garrappa, S., et al. 2024b, A&A, 681, A119 Roming, P. W. A., Kennedy, T. E., Mason, K. O., et al. 2005, Space Sci. Rev., 120, 95\n- Sahakyan, N., Giommi, P., Padovani, P., et al. 2022, Monthly Notices of the Royal Astronomical Society, 519, 1396\n- Scargle, J. D., Norris, J. P., Jackson, B., & Chiang, J. 2013, ApJ, 764, 167\n- Schlafly, E. F. & Finkbeiner, D. P. 2011, ApJ, 737, 103\n- Stahl, B. E., Zheng, W., de Jaeger, T., et al. 2019, MNRAS, 490, 3882  \nArticle number, page 10 of 15  \nTeräsranta, H., Tornikoski, M., Mujunen, A., et al. 1998, A&AS, 132, 305  \n- Tomsick, J., Zoglauer, A., Sleator, C., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 98\n- Tonry, J. L., Stubbs, C. W., Lykke, K. R., et al. 2012, ApJ, 750, 99\n- Willingale, R., Starling, R. L. C., Beardmore, A. P., Tanvir, N. R., & O'Brien, P. T. 2013, MNRAS, 431, 394\n- Zhang, H., Böttcher, M., & Liodakis, I. 2024, ApJ, 967, 93\n- Zhang, S. N., Feroci, M., Santangelo, A., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9905, Space Telescopes and Instrumentation 2016: Ultraviolet to Gamma Ray, ed. J.-W. A. den Herder, T. Takahashi, & M. Bautz, 99051Q", 'Appendix A: Tables of parameter space boundaries and model parameters': "For the search of the best-fit parameters of leptonic one-zone models, a grid-scan method was used. The region of parameter space that was studied is defined by the boundaries presented in Table A.1. The intervals where equally divided into 10 points (including the boundaries) in linear space for B ' , Γ b, and α e , and in decimal logarithmic space for all other parameters.  \nTable A.1: List of leptonic model parameters and boundaries of the parameter space.  \nDuring the search for best-fit parameters of one-zone leptohadronic models (see Section 3.4), the grid scan of four hadronic parameters was performed. The subregion of the parameter space that was studied is defined by the boundaries presented in Table A.2. The intervals of γ ' min p and γ ' max p where equally divided into 30 points (including the boundaries) in decimal logarithmic space. The L ' p interval was divided into 29 points (30 together with zero) in decimal logarithmic space. Five equally distant in linear space points from the range of α p were selected including boundaries. The unphysical combinations where γ ' min p > γ ' max p were excluded.  \nTable A.2: List of hadronic model parameters and boundaries of the parameter space.  \nTable A.3 shows the best-fit parameter values of leptonic models. The solutions are divided into two cathegories, slow and fast, based on their bulk Lorentz factor. The parameters are given for each of the four activity states of PKS 0735 + 178.  \nTable A.4 shows values of hadronic parameters of the models presented in Fig.5. The model number corresponds to the legend of Fig. 5.  \nTable A.4: Parameters of the leptohadronic models shown in Fig. 5. The leptonic parameters in all models are R ' b = 10 15 . 57 cm, B ' = 2.2 G, Γ b = 15.1, γ ' min e = 10 3 . 46 , γ ' max e = 10 4 . 34 , α ' e = 2.42, L ' e = 10 42 . 65 erg s -1 .  \nTable A.5 presents the best-fit parameter values of the leptohadronic models under the condition of neutrino rate maximization. They correspond to the models shown in Fig. 6.", 'Appendix B: Paramater space of one-zone leptonic models': 'Similarly to Fig. 4, Figures B.1-B.6 show the value of reduced χ 2 for variations of every pair of leptonic model parameters when the remaining parameters are fixed to the best fit. The bestfit value is marked with a pink circle in all plots.  \nTable A.3: Best-fit leptonic model parameters.Table A.5: Best-fit leptohadronic model parameters with neutrino event rate maximization condition.  \nFig. B.1: Values of the reduced χ 2 around best-fit slow solution for the neutrino arrival SED.  \n<!-- image -->  \nFig. B.2: Values of the reduced χ 2 around best-fit slow solution for the quiescent state SED.  \n<!-- image -->  \nFig. B.3: Values of the reduced χ 2 around best-fit fast solution for the quiescent state SED.  \n<!-- image -->  \nFig. B.4: Values of the reduced χ 2 around best-fit slow solution for the γ ray flare peak SED.  \n<!-- image -->  \nFig. B.5: Values of the reduced χ 2 around best-fit fast solution for the γ ray flare peak SED.  \n<!-- image -->  \nFig. B.6: Values of the reduced χ 2 around best-fit slow solution for the after-flare SED.  \n<!-- image -->  \nFig. B.7: Values of the reduced χ 2 around best-fit fast solution for the after-flare SED.  \n<!-- image -->'}