Hadronic Processes at Work in 5BZB J0630 − 2406

Recent observations are shedding light on the important role that active galactic nuclei play in the production of high-energy neutrinos. In this study, we focus on one object, 5BZB J0630 − 2406, which is among the blazars recently proposed as associated with neutrino emission during the ﬁ rst 7 yr of IceCube observations. Modeling the quasi-simultaneous, broadband spectral energy distribution, we explore various scenarios from purely leptonic to leptohadronic models, testing the inclusion of external photon ﬁ elds. This theoretical study provides a complementary testing ground for the proposed neutrino – blazar association. Despite being historically classi ﬁ ed as a BL Lac, our study shows that 5BZB J0630 − 2406 belongs to the relatively rare subclass of high-power ﬂ at-spectrum radio quasars. Our results indicate that interactions between protons and external radiation ﬁ elds can produce a neutrino ﬂ ux that is within the reach of the IceCube detector. Furthermore, the spectral shape of the X-ray emission suggests the imprint of hadronic processes related to very energetic protons.


INTRODUCTION
Active galactic nuclei (AGN) are among the most energetic and powerful objects in the universe.They are powered by a supermassive black hole (SMBH), in some of them a relativistic jet can be present and detected across the electromagnetic spectrum, from the radio to the very high energy (VHE) band.Efficient particle acceleration mechanisms such as magnetic reconnection (Blandford et al. 2017) and shock acceleration (Pelletier et al. 2019;Lemoine et al. 2019) can occur in distinct regions of the jet.Standing and moving features observed in AGN jets (Marshall et al. 2002;Jorstad et al. 2013;Lister et al. 2021), also reproduced numerically as shocks (Fromm et al. 2016;Fichet de Clairfontaine et al. 2021, 2022), support the idea that the observed gaetan.fichet-de-clairfontaine@physik.uni-wuerzburg.demulti-wavelength (MWL) electromagnetic emission can be explained trough radiative cooling of those accelerated particles along the jet.While to date the overwhelming majority of AGN can be explained invoking the leptonic framework, from the theoretical point of view relativistic jets harbored in AGN may be capable of accelerating hadrons.Within such scenarios, high-energy (HE; TeV/PeV energies) neutrinos are a natural by-product.
Located at the geographic South Pole, the IceCube observatory is the most sensitive high-energy (≳ TeV) neutrino detector currently operating.Since the beginning of its science operations in 2008, it allows us to study the putative neutrino counterparts of astrophysical sources.Although no firm associations between individual high-energy IceCube events and cosmic sources have been established to date, several claims of associations with AGN have been made at different statistical levels, e.g. the blazar TXS 0506+56, IceCube Collaboration et al. (2018a,b), which represent a subclass of AGN with the jet pointed directly at the observer.
Other approaches to studying such correlations exploit the time-integrated neutrino information over a given period of time, as in the case of observational evidence of neutrino emission in the direction of the Seyfert galaxy NGC 1068 (IceCube Collaboration et al. 2022), or through stacking analyses of populations of sources (e.g.Aartsen et al. 2017a;Padovani et al. 2016;Aartsen et al. 2017b;Abbasi et al. 2021;Plavin et al. 2021;Abbasi et al. 2023).
A recent work reports evidence for a statistically significant correlation between blazars listed in the fifth data release of the Roma-BZCat catalog (5BZCat, see Massaro et al. 2015), and a sample of IceCube hotspots, i.e. anisotropies in the distribution of IceCube events (Buson et al. 2022a,b).The study is based on the 7-yr IceCube southern sky map published by the IceCube collaboration and highlights 10 objects as candidate HE neutrino emitters, i.e.PeVatron blazars.Consistent findings are reported when expanding the investigation to the Northern Hemisphere with the latest sky map released by the IceCube collaboration (Buson et al. 2023).Investigations employing different analysis methodologies have provided mixed results.For instance, Bellenghi et al. (2023) performed an independent analysis of the public IceCube dataset, confirming the association established by Buson et al. (2022a,b) with the 7-yr Southern Hemisphere dataset.However, when analyzing the 10-year dataset, they did not observe a similar correlation.Such 10-yr sky map appears overall different from the one published by IceCube collaboration with the same dataset (Aartsen et al. 2020), and the one used by (Buson et al. 2023) for the 10-yr Northern Hemisphere analysis.The discrepancies can be attributed to the different likelihood formalism and the coarseness of the detector response matrices employed, that lead to an overall worst sensitivity, as acknowledged in Bellenghi et al. (2023).
Given the unascertained electromagnetic / neutrino relation, the statistical analysis was based solely on the positions of the objects.No a-priori selection was applied to the blazar sample, neither based on the objects' classification nor on their electromagnetic properties.
This paper aims to provide an initial characterization of the underlying physics of PeVatron blazars from the theoretical perspective (Fichet de Clairfontaine et al. 2023).
To this extent, we focus on one object for which broad, MWL simultaneous observations are available, namely 5BZB J0630−2406 (a.k.a.TXS 0628−240, WISE J063059.51−240646.2), while future works will address the physical properties of the sample (Azzollini et al. 2023).
Although the redshift of the source is unknown due to the lack of emission lines in the optical spectrum, a lower limit has been established z ≥ 1.239 (Shaw et al. 2013), and ab-sorption lines have been observed (Mg I, Fe II, Al II) that indicate the presence of the host galaxy (Rau et al. 2012a).
Theoretical modeling of the MWL emission of candidate neutrino emitters have shown the key role of external fields in the neutrino production (Dermer et al. 2014;Oikonomou et al. 2019) and markers of hadronic processes imprinted in the observed electromagnetic spectral energy distribution(SED, Reimer et al. 2019;Cerruti et al. 2018;Gao et al. 2019;Petropoulou et al. 2020;Rodrigues et al. 2021).Here, we model the SED of the object of interest, exploiting simultaneous and quasi-simultaneous MWL observations.
The paper is organized as follows.Sect. 2 introduces the neutrino observations that provided evidence of neutrino emission in the direction of 5BZB J0630−2406.The Sect. 3 displays the MWL observations and data reduction.The Sect. 4 discusses the temporal variability of this source.We present our approach to reproduce the observed SED with a one-zone leptonic and lepto-hadronic numerical models, with a complete description in Sect. 5.Then, in Sect.6, we present the results of our parameter space exploration in different scenarios, including scenarios where the X-ray band is dominated by either synchrotron radiation or by hadronic cascade processes related to secondaries particles.The characteristics of the source derived from the various scenarios will be detailed, as well as the expected neutrino event rates.In Sect.7, we discuss our findings and their implications in the context of recent neutrino -blazar associations.Throughout the paper, all primed quantities are evaluated in the rest frame of the relativistic jet.We also assume a flat ΛCDM cosmology with H 0 = 69.6 km•s −1 •Mpc −1 , Ω 0 = 0.29 and Ω Λ = 0.71.
2. NEUTRINO OBSERVATIONS 5BZB J0630−2406 has been proposed as candidate counterpart of the IceCube hotspot IC J0630−2353 (Buson et al. 2022a).The study presented in Buson et al. (2022a) is based on a Southern Hemisphere analysis (δ < −5) of the allsky map that spans the IceCube observations from 2008 to June 2015.This 7-yr sky map encodes the time-integrated information regarding long-term point-source neutrino emitters and is built using events of energy-proxy ≳ 100 TeV.In the 7-yr sky map, the direction of the IceCube hotspot IC J0630−2353 evidences an anisotropy in the spatial distribution of events (see Fig. 1) and, thus, may hint for the presence of astrophysical sources.The object is located at an angular separation of 0.28 • from the hotspot IC J0630−2353, within the association radius of 0.55 • derived for this dataset in Buson et al. (2022a).The evidence for a spatial correlation between IC J0630−2353 and the blazar 5BZB J0630−2406 suggests that this blazar may contribute to the observed anisotropy.

MULTIWAVELENGTH OBSERVATIONS AND DATA ANALYSIS
In the following, we provide a description of the MWL data collection and reduction.In Figure 2, we present the MWL light curves, from near-infrared (NIR) to γ rays.The orange vertical line denotes the epoch of the simultaneous observations, used in the SED modeling presented in Sect. 5.

Radio observations
Archival radio data are available from the Academy of Sciences Radio Telescope (RATAN) 1 , the GaLactic and Extragalactic All-sky Murchison Widefield Array (GLEAM, Hurley-Walker et al. (2017)) and the Australia Telescope 20-GHz (AT20G, Mahony et al. (2011)) catalog.They are not used in the modeling and rather considered as upper limits for the fit.At such lower frequencies, the radio flux is likely originated from extended regions of the jet, as the electrons cool down via low-energy synchrotron over a long period of time while propagated through the jet.

Optical observations
Over the past decade 5BZB J0630−2406 has been monitored at optical wavelengths by several programs.The Katzman Automatic Imaging Telescope (KAIT, Filippenko et al. 2001;Li et al. 2003) performed optical photometric observations in the Clear band, which is close to the R band (6410 Å), within the Fermi-LAT AGN monitoring program (Cohen et al. 2014).The data are calibrated to Landolt R band, and corrected for the galactic reddening with A R = 0.138 mag, following Schlafly & Finkbeiner (2011).The magnitudes are converted to flux units using the zero points (Bessell et al. 1998).For earlier observations shown as black circles in Fig. 2, around MJD 55500, a manual comparison was made with calibrated data at corresponding points in time.
5BZB J0630−2406 has also been monitored by the All-Sky Automated Survey for Supernovae (ASAS-SN, Shappee et al. 2014;Kochanek et al. 2017) in the V-band and g-band.
The magnitude corrections are applied using the reddening coefficient E(B − V) = 0.054739 according to Schlafly & Finkbeiner (2011) with ratio of the extinction reddening A λ /E(B − V) for each filter taken from Fitzpatrick (1999).The conversion to flux units follows the same methodology as explained previously for KAIT.The Zwicky Transient Facility (ZTF, Masci et al. (2018)) monitored the source in the r-band and g-band, as well as the Small and Moderate Aperture Research Telescope System (SMARTS, Bonning et al. 2012) in the R-band.The collected magnitudes are corrected for reddening and converted to flux units consistent to KAIT and ASAS-SN.
Swift-UVOT observations are available in three optical filters (U, B, and V) between 2009-02-01 (MJD 54863) and 2014-11-10 (MJD 56971).The UVOT data above ν ≥ 10 15 Hz are affected by the Lyα forest absorption and hence not used in this study.We extracted source counts using apertures with a radius of 5", while the background counts were obtained from an annulus region centered on the source position where no other sources were visible.To compute the magnitudes, we used the UVOTSOURCE tool in HEASOFT.Extinction corrections were applied following Schlafly &Finkbeiner (2011) andFitzpatrick (1999) for each filter.Then, using the zero points derived from Breeveld et al. (2011), we subsequently converted the magnitudes to fluxes following the procedure described in Poole et al. (2008).

The Neil Gehrels Swift observatory
Overall, Swift-XRT (Gehrels et al. 2004)   period of the IceCube observations.To analyze the Swift-XRT data, the standard analysis tools provided by the HEA-SOFT2 V.6.31.1 software package were used.
The XRT observations were all in photon counting (PC) mode, and since the count rates were always below 0.5 counts s −1 , no pile-up correction was necessary.To extract the source events, a circular region with a radius of 45" was chosen within the 0.3-10 keV energy range.The background events were extracted from the annular region around the source position.Given the low photon statistics, we performed the spectral analysis with the Cash statistics method by rebinning the spectrum to ensure a minimum of one count per bin.

XMM-Newton and NuSTAR data reduction
5BZB J0630−2406 was targeted quasi-simultaneously by XMM-Newton (obsID: 0740820401; exposure: 9 ks) and NuSTAR (obsID: 60001140002; exposure: 66 ks) on October 17-18, 2014.We retrieved the XMM-Newton pn source and background spectra and associated matrices from the 4XMM catalog 3 (Webb et al. 2020).The source spectrum has been binned with 1 count per bin, to avoid having empty bins that can affect the spectral fit.To generate the NuSTAR spectra, we followed the standard data reduction procedure.Specifically, the data have been processed using the NuSTAR Data Analysis Software (NUSTARDAS) version 2.1.1.The raw event files are calibrated by the nupipeline script, using the response file from the Calibration Database (CALDB) version 20210202.The source and background spectra are extracted from a 45 ′′ (≈ 50% of the encircled energy fraction at 10 keV) circular region, centered at the optical position of the source and in a nearby (∼3 ′ separation) region which was visually inspected to avoid any possible contamination, respectively.Using nuproducts scripts, we then generated source and background spectra files, along with the corresponding ARF and RMF files.Finally, the NuSTAR spectra are grouped with 1 count per bin, using grppha.This procedure has been performed on both the NuSTAR focal plane modules, FPMA and FPMB.

XMM-Newton and NuSTAR spectral fitting results
We fitted the XMM-Newton and NuSTAR spectra of our target using the XSPEC (Arnaud 1996) software, version 12.12.1.We fixed the metal abundance to Solar metallicity using the abundances from Wilms et al. (2000), while the photoelectric cross-sections for all absorption components are those derived by Verner et al. (1996).The Galactic absorption column density is fixed to N H,gal = 7.5 × 10 20 cm −2 (Kalberla et al. 2005).To maximize the spectral statistic, we analyze the data with the Cash statistic (Cash 1979), which uses a Poisson likelihood function and is hence most suitable for low numbers of counts per bin.As mentioned in the previous section, we bin the spectra with 1 count per bin.The XMM-Newton pn spectrum is fitted in the 0.3-10 keV band, while the two NuSTAR FPMA and FPMB spectra are fit in the 3-70 keV band.
Following the standard procedure for fitting X-ray spectra of blazars, we first fit our data with a simple power-law model.We also add a cross-normalization constant to take into account the fact that the NuSTAR observation is significantly longer than the XMM-Newton one, and to model systematic cross-instrument offsets in flux.Finally, we include a column density at the redshift of the source, N H,z,ISM , 3 More details at http://xmm-catalog.irap.omp.eu/source/207408204010001.assuming a redshift of 1.239, to account for a possible contribution of the interstellar medium (ISM) to the absorption of X-ray photons.We measure a best-fit X-ray photon index Γ X = 3.03 ± 0.12, with an ISM column density N H,z,ISM = 3.0 +1.1 −1.0 × 10 21 cm −2 , and a XMM-Newton-NuSTAR crossnormalization C XMM−NuS = 1.90 +0.41  −0.34 .The best-fit statistic for this fit is Cstat/d.o.f.= 1377.8/1408.
We then fit the data with a log-parabola model, which has been shown to accurately describe the X-ray spectral shape of blazars (see, e.g., Bhatta et al. 2018).In particular, Middei et al. (2022) analyzed the NuSTAR spectra of 126 blazars and found that in some cases a log-parabola model provides a more statistically accurate description of a blazar X-ray spectrum.The log-parabola model is described by the following equation, where E 0 is the reference energy, which we fix as 5 keV (in the observer frame) following the results of previous works (e.g., Massaro et al. 2004;Baloković et al. 2016;Middei et al. 2022), while α and β are the photon index and the curvature parameter, respectively (e.g., Massaro et al. 2004), and K is the model normalization.
Finally, we performed a fit with a broken power-law model, which is also commonly used to describe the X-ray spectra of blazars.The parameterization of the broken power-law is, (2) Here, Γ 1 and Γ 2 are the low and high-energy photon indexes, K is the normalization parameter, and E b is the (restframe) energy of the break.The best-fit parameters are −0.37 The best-fit statistic for this fit is Cstat/d.o.f.= 1365.5/1406.The difference in Cstat between this model and the simple power-law one is ∆Cstat = 1377.8− 1365.5 = 12.3.Since the broken power-law model has 2 fewer degree of freedom than the simple power-law one, the broken power-law model is statistically preferred at the > 99.8 % confidence level (i.e., at a ≳ 3.1 σ significance).In summary, a two-component model is statistically favored with respect to a simple power-law one, while from a statistical point of view the log-parabola model and the broken power-law one are fully consistent.Differently from a previous study (Ackermann et al. 2016), our more sensitive analysis supports the presence of a break in the X-ray spectrum.We report in Fig. 3 the source light curves from both the NuSTAR cameras and, as it can be seen, no clear variability trend is observed within the NuSTAR observation.More in detail, we fit the two light curves with a constant and obtain best-fit count rates r FPMA = 0.45 cts•s −1 and r FPMB = 0.44 cts• s −1 , with reduced chi square (χ 2 /d.o.f.) FPMA = 9.49/20 and (χ 2 /d.o.f.) FPMB = 0.86/19.We also note that the XMM-Newton-NuSTAR cross-normalization values we measured in the log-parabola and broken power-law fits, while slightly higher than expected in quasi-simultaneous observations, have been observed in other XMM-Newton-NuSTAR quasisimultaneous observations (e.g., Marchesi et al. 2019) and may be explained by a simple cross-instrument calibration offset.

Fermi-Large Area Telescope observations
The study utilized gamma-ray data taken from Ackermann et al. (2016), and obtained from Fermi-LAT observations spanning the period between August 4, 2008, and January 31, 2015 within the energy range of 100 MeV to 500 GeV.This time period is contemporaneous to the simultaneous Xray data introduced in the previous section.

TEMPORAL VARIABILITY
The long-term optical light curves presented in Fig. 2 highlight a similar flux variability pattern over the monitored bands.According to the 4FGL-DR4 catalog (Ballet et al. 2023), the source is variable at γ-rays on ≳ year-timescales.To characterize statistically-significant variations in the γray light curve, we adopt the Bayesian algorithm available in Astropy 5 (Astropy Collaboration et al. 2013).To determine the optimal value of the prior for the number of blocks, we use the empirical relation evaluated in Scargle et al. (2013) for the probability to falsely report a detection of a change point, setting it to 0.05.Applying this approach to the yearly binned light curve, the first ∼ 7-yr of LAT observations, up to 5 More details at http://docs.astropy.org/en/stable/api/astropy.stats.bayesianblocks.html∼ MJD 57250, are overall consistent with a steady state, indicating that if variability is present at these frequencies it may be below the sensitivity of the LAT data.During the more recent LAT observations, since ∼ MJD 57250, the source is undergoing a long-term enhanced state.
Optical evidence of year-long modulations of the flux, on timescales of ∼ 3/4 yr.This is consistent with the behavior traced by the > 15-yr γ-ray light curve, including a major flux enhancement observed around MJD 57800.The sparse X-ray/UV data mimic the overall variability pattern at other frequencies.The high-cadence KAIT monitoring, with sampling as short as a ∼ 3 days cadence, evidences statistically significant (≳ 6σ) changes between consecutive observations (Ackermann et al. 2016).This suggests that at least some of the jet's emission arises in compact R ≲ 10 16 cm regions, in the observer's frame.vations, i.e. 2008-2015.Within this time range, MWL contemporaneous observations of 5BZB J0630−2406 were collected with the gamma-ray burst Optical/Near-Infrared Detector (GROND) instrument at the 2.2 m MPG telescope at the ESO La Silla Observatory (Greiner et al. 2008), the Swift Niel Gehrels Observatory, XMM-Newton and NuSTAR (Harrison et al. 2013) satellites.They were performed around October 17, 2014 , i.e. ∼ MJD 56948, close to the time period of the 7 yr neutrino observations.These simultaneous (≲ 1 day) multiwavelength data constrain the X-ray part of the synchrotron component.To constrain the second hump of the SED, we employ contemporaneous observations carried out by the Fermi-LAT in the MeV-GeV range.As discussed in Sec. 4, the source of interest is characterized by long-term variability in the Fermi-LAT band, and no significant variations in the flux are observed during the first 7-yr of monitoring.Therefore, for the SED modeling we employ the LAT spectrum integrated over the first 7 years of observations (from August 04, 2008 to August 15, 2015, i.e.MJD 54682 − 57250) available from the literature (Ackermann et al. 2016).We exclude the UVOT data above ν ≥ 10 15 Hz as they are affected by the Lyα forest absorption.Fig. 5 and Fig. 6 show the contemporaneous broadband SED built with these data (black points).Further archival radio, optical, and near-IR observations are displayed for comparison (gray points), and not included in the SED modeling.The red downward triangle represents a limit on accretion disk luminosity from Ghisellini et al. (2012).

One zone model
We explore different scenarios to describe the MWL emission of 5BZB J0630−2406 using the time-dependent code AM 3 (Gao et al. 2017).This code is able to solve the system of coupled differential equations that describe the transport of particles through the relativistic jet.Relativistic electrons and protons are assumed to be accelerated initially and injected into a single zone where they can radiate their energy and interact with external photon fields.For more details on the implementation of the synchrotron, inverse Compton scattering, photo-photon pair annihilation, hadronic processes (i.e., Bethe-Heitler pair production and photo-pion production), we refer the reader to Gao et al. (2017).We also note that the synchrotron self-absorption opacity of the blob is treated only for the electrons (we expect a negligible contribution from the proton population in the radio band).
As represented in the Fig. 4, the emission region is modeled as a sphere of radius R ′ b that moves at relativistic speed with a Lorentz factor of Γ b .From Earth, the blob is observed with an observation angle θ obs = 1/Γ b , with quantities Doppler-shifted according to the Doppler factor δ D = Γ b .A magnetic field with a strength of B ′ is assumed to be present in the blob and isotropic, i.e. with no preferential direction.The blob is located at a distance of R diss from the supermassive black hole, also known as the dissipation radius, considered static here.Since the energy density of external photon fields varies with R diss , this parameter is a crucial aspect of our model.In our model, relativistic electrons and/or protons are injected into the blob.For the electrons, we use a broken power-law distribution given by, where γ ′ e,min and γ ′ e,max are the cut-off values, and γ ′ e,brk is the break energy.The choice of this distribution is supported by a pre-fit analysis on the observed SED that it used for initial guesses only (as it can be seen in Ghisellini et al. (2012)).For protons, we assume a simple power-law distribution, The parameters N e and N p are the normalization factors.We ensure that the constraints on γ ′ e,brk and γ ′ e/p,max are compatible by respectively equating the synchrotron cooling timescale τ syn with the adiabatic timescale τ ad = 2R ′ b /c and the acceleration timescale, τ acc 6 with the shortest cooling timescale.For γ ′ e/p,max , we also check that the Hillas criterion is satisfied.
Given the expected large redshift, we incorporate the effect of absorption by the extragalactic background light (EBL) in our model.To achieve this, we use the Python library ebltable 7 , which is based on the model presented by Domínguez et al. (2011).
The presence of a luminous accretion disk and radiation fields has been suggested by previous literature studies ( Ghisellini et al. 2012;Padovani et al. 2012).External radiation field radiation can either interact directly with the source in the jet or be re-processed by a broad line region (BLR).A dust torus is also present in this model, as its presence has been considered in numerous previous studies (Finke 2016;Murase et al. 2014;Oikonomou et al. 2021).An upper limit on the disk luminosity has been derived in Ghisellini et al. (2012), (5) It was obtained assuming an empirical relation between the BLR and gamma-ray luminosities (L BLR ∼ 4 × L 0.93 γ ).Therefore, the model presented accounts for both external radiation fields originated either from the accretion disk and the dust torus, modeling them in the observer frame as a single temperature black body emission (Dermer & Menon 2009) for simplicity, as visible in Rodrigues et al. (2019).This choice is also motivated by the lack of optical lines, as it implies a higher uncertainty on the black hole mass and the disk type.
Following Ghisellini et al. (2017), we assume in this study that the luminosity and the size of the dust torus can be derived from the disk luminosity and the dust torus temperature, where L disk , T disk , and T torus are free parameters in the model.Our model assumes that the emission from the accretion disk is generated by a single-temperature blackbody that interacts with the source jet through a Doppler de-boosting effect, which is proportional to Γ 2 b .In fact, we consider a scenario in which a single-temperature accretion 6 Here, we define the acceleration timescale as, where we assume η = 0.1 (Cerruti et al. 2015). 7Available at https://github.com/me-manu/ebltabledisk emits an isotropic radiation "behind" the blob, as the jet is pointing at the observer.The emission of the dust torus is boosted in the blob frame and also consider isotropic.Additionally, the emission can be scattered by the BLR, which we assume is modeled as a thin shell located at a distance R BLR from the central engine and radiating the luminosity L BLR .Both L BLR and R BLR are derived following Ghisellini & Tavecchio (2008), We assume that the BLR reprocesses 10% of the disk emission isotropically in the rest frame of the black hole (Sbarrato et al. 2012).In fact, most of the reprocessed flux is emitted in emission lines (e.g.Lyman α) and not as a thermal continuum (which accounts for 1% (Blandford & Levinson 1995;Murase et al. 2014)).In fact, we consider that most of them will lie on the disk emission range (Lyman α is situated at ∼ 10 15 Hz) and are outshined by the non-thermal continuum from the jet (Rodrigues et al. 2021).In the blob frame, this radiation field is enhanced by a factor δ 2 .This Doppler factor depends on Γ b but also on the dissipation radius R diss , so that the energy density perceived by the blob will be lower outside the BLR radius, according to scaling factors calculated in Ghisellini & Tavecchio (2009).
To explore the parameter space, we develop a parameter search algorithm which is parallelized, with initial and mutation guesses stored on the first and successive central processing units (CPU), respectively.Gaussian noises are applied to each parameter as mutations.We use the least squares function from the scipy.optimize8library to optimize the parameters.For each parameter set and CPU, the algorithm calls the AM 3 code, which returns the simulated SED.The optimization is based on minimizing the residuals.Finally, we select the parameter sets associated with the best χ 2 /d.o.f.value.Here, the d.o.f.term refers to the number of degrees of freedom.If the χ 2 /d.o.f.value is not acceptable, we use this latest set of parameters found as input for the next generation, with the other solutions being mutations of this one.The process stops when the χ 2 /d.o.f.value is sufficiently low, i.e., χ 2 /d.o.f.≤ 2.More details on the parameter space research, and on the evaluation of initial guesses, are show in Appendix B.

Purely leptonic solution
In our initial attempt to model the SED of 5BZB J0630−2406, we considered only a population of relativistic electrons.The simulated SED is shown in Fig. 5, and the fitted parameters can be found in Table 1 labeled L. In this figure, the labels SY and IC stand respectively for synchrotron and inverse Compton, from either the leptonic population (e ± ) or from γ − γ pair productions.In this scenario, the data points until the X-ray are explained by synchrotron emission only.The GeV gamma-ray data are partly explained by synchrotron self-Compton (SSC).The best solution found in that scenario is close to equipartition with u ′ e /u ′ b ∼ 3.9, where u ′ e is the electron energy density and u ′ b the magnetic energy density.It should be noted that we did not take into account the presence of cold protons in this model.Here, we derived a relatively high value for γ e,min , this can be due to the prior acceleration (injection) of a truncated power-law from a region closer to the black hole.The cooling process might also be inefficient for lower energies, leading to the development of the low energy tail outside the zone model here (Katarzyński et al. 2006), explaining the radio counterpart.In fact, the extended jet is expected to contribute to the integrated radio flux where those low energy electrons will cool down (Plavin et al. 2022).Such a model is beyond the scope of this paper.Specific interactions between the electron and proton populations can also explain high γ e,min values (Zech & Lemoine 2021).
To account for the limit proposed by Ghisellini et al. (2012), we included the black body emission of the accretion disk.The final disk luminosity obtained is L disk = 4.8×10 45 erg•s −1 , which explains the peaky feature observed at the highest γ-ray energies with external Compton interactions.In this model, the blob is located at a dissipation radius of R diss /R BLR = 1.7, indicating that the influence of the BLR radiation field is still significant.The dust torus black body is also present but subdominant.
From the accretion disk parameters obtained, we can derive a black hole mass of ∼ 10 10 M ⊙ and an associated Eddington luminosity of ∼ 3 × 10 48 erg • s −1 .As suggested in Sbarrato et al. (2012), we use the accretion regime η = L BLR /L Edd as physical criteria to distinguish BL Lacertae (BL Lacs) from FSRQs.We find an accretion regime of η ∼ 2 × 10 −4 , close to the boundary that separates FSRQs from BL Lacs (η = 5 × 10 −4 , Sbarrato et al. 2012).Finally, the shortest timescale variability τ var derived from this model is 3.7 days in the observer's frame, close to the derived value from optical variability.
Given the X-ray spectrum, we further tested the possibility to reproduce the break in the X-ray band with the purely leptonic model, but we were not able to find an acceptable solution (χ 2 d.o.f.≤ 2, by lowering γ e,min for example).This suggests that, if the break is present, it may be interpreted as the presence of additional processes.As we discuss in the following sections, a hadronic component is capable to successfully account for the broken spectral shape observed in the X-ray band as well as the putative neutrino emission.

Lepto-hadronic solution
We use the parameters found in the purely leptonic solution as starting point for the parameter space research of the lepto-hadronic scenario.The final parameters are displayed in the Tab. 1, labeled as LH, and the SED is displayed on the Fig. 6.Here BH stands from Bethe-Heitler reaction, while p − γ stands for photo-pion production.Protons are injected with a simple power-law index of 2.0 and the hadronic processes remain globally sub-dominant, except in the X-ray and MeV bands.Within the lepto-hadronic solution, we find values that are in agreement with those found for the purely leptonic model with an accretion regime of η ∼ 2 × 10 −4 , and an Eddington luminosity of ∼ 3 × 10 48 erg • s −1 .Similarly, the shortest timescale variability derived here, τ var ∼ 3.3 days, is consistent with the value derived from optical variability.The cascade component accounts for the X-ray flux, particularly at the highest energy.
The higher energy MWL peak is mainly explained by SSC and EC interactions with the BLR, consistently with the pure leptonic solution.Similarly, we find a dissipation radius close to 1.6 R BLR .The cascade component accounts for the hard X-ray data, leading to a steeper index for the electron broken power-law.Indeed, significant contributions from γ − γ pair production in the GeV energy band can also be observed, which are more pronounced than the ones observed in the leptonic model due to the additional photon fields from mesons disintegration.Although the model finds u ′ p /u ′ b ∼ 10 3 (where u ′ p is the proton energy density), far from the equipartition, the proton luminosity remains below the Eddington limit.
The model presented here may be considered a viable solution for the efficient neutrino production case found by our modeling.As the X-ray data provide upper-limits for the cascade component, a broad range of solutions can be obtained with less energetic protons and, hence, lower predicted neutrino fluxes.

DISCUSSION
Although still widely used, the historical classification of blazars based purely on observational characteristics, e.g. the optical spectrum and the location of the low-energy peak, has been put into question for a long time.Since the availability of large samples of blazars, thanks to Fermi-LAT observations, Ghisellini et al. (2011) has debated in favor of a more physical distinction for blazars based on the luminosity of the  (Ghisellini et al. 2012).
BLR measured in the Eddington units.Indeed, the BLR and Eddington luminosities are respectively related to the optical spectrum and the jet power (energy injected in the primary relativistic electrons).This can set a divide approximately where the disc transitions from a radiatively efficient to an inefficient regime.
Prior to being included in the sample of PeVatron blazars and proposed as a high-energy neutrino emitter, 5BZB J0630−2406 stood out in the literature due to its peculiarities.Historically, it has been classified as a BL Lac object due to the featureless optical spectrum and its high synchrotron peak, with ν sy pk ∼ 10 15 Hz.5BZB J0630−2406 was pinpointed as an exemplary blazar, it displays properties typical of "blue flat spectrum radio quasar" (Ghisellini et al. 2011;Padovani et al. 2012;Ghisellini et al. 2012, a.k.a."high-power high-synchrotron-peak blazars"), i.e. high-emitting power sources that are intrinsically FSRQs where their broad emission lines are swamped by the jet syn-chrotron emission.For reference, this typology of blazars has also been called "masquerading BL Lacs" (Giommi et al. 2013;Padovani et al. 2019a).In contrast to "true" highfrequency-peaked BL Lacs which have intrinsically poor radiation fields, these objects host powerful jets and radiatively efficient accretion.
The study presented here allows us to provide conclusive evidence to the earlier speculations regarding the peculiar nature of 5BZB J0630−2406.The SED modeling reveals a bright accretion disk with L disk ∼ 4 × 10 45 erg • s −1 .The presence of external fields that are partly re-processed by the BLR naturally explains the peaky feature in the γ-ray band.Its accretion regime, i.e. the energy injected into these external fields, is of the order of L BLR /L Edd ∼ 2 × 10 −4 , close to the values physically suggested for FSRQs (Ghisellini et al. 2011).The relatively high accretion regime of 5BZB J0630−2406 is supported also by the ratio of the γray luminosity L γ (0.1 − 100 GeV) and the Eddington lumi- The expected flux of µ-neutrinos as observed on Earth is shown in Fig. 7 and is close to the IceCube sensitivity (Aartsen et al. 2017b, estimated for a ∝ E −2 spectrum).Our leptohadronic model predicts N events = 4.82 +5.18  −3.82 over a livetime period of 7 years, as shown on Fig. 8. Fig. 8 shows the temporal evolution of the number of events N events assuming a constant flux over the livetime of the IceCube detector, accounting for the different strings configurations, with 3σ uncertainties.Assuming low counting statistics, e.g.Poisson statistics, we can evaluate the probability of N events detection against a null hypothesis of no detection, under the form of a p-value.Considering the first 7-year livetime Aartsen et al. (2017c), we find a p-value of 0.03, indicating a small tension with the null hypothesis of no detection.It should be noted that, in the most conservative scenario, our model still predicts a minimum of N min = 4.82 − 3.82 = 1 event, with a pvalue higher than 0.05.We find that 5BZB J0630−2406 contributes for ≤ 1% to the IceCube diffuse muon neutrino flux (IceCube Collaboration et al. 2022).Similar values were obtained for TXS 0506+056 (Aartsen et al. 2016;IceCube Collaboration et al. 2018c).
We can derive the baryonic loading in the blob frame, which represents the ratio between the proton luminosity and  the γ-ray luminosity, denoted by ξ = L ′ p /L ′ γ .Based on our lepto-hadronic model, we find a loading factor of ξ ∼ 10 3 .This value appears notably high when compared to the typical order of magnitude, which is around 100 as derived in Murase et al. (2014).Nevertheless, the value obtained here is consistent with the range derived in Palladino et al. (2019), as the intermediate BL Lac -FSRQ objects show a higher ξ, but still respecting the IceCube stacking limit (Aartsen et al. 2015).Appendix A discusses our findings in the context of literature studies that tackled 5BZB J0630−2406 as an ultrahigh-energy cosmic ray (UHECR) candidate.in the context of both a purely leptonic and mixed leptohadronic scenario.We summarize the main findings in the following.
• Despite being formally classified as a highsynchroton-peaked BL Lac object, based on its featureless optical spectrum, the intrinsic nature of 5BZB J0630−2406 is that of a high-power FSRQ.It hosts a standard accretion disk and BLR; the optical emission lines elude direct observation, being the optical band swamped by the non-thermal continuum.The combination of efficient external radiation fields and enhanced particle acceleration efficiency offers ideal conditions for the production of neutrinos.
• The presence of a bright accretion disk is confirmed by a peaky feature in the γ-ray spectrum at the highest observable LAT energies.In FSRQ sources, a similar spectral shape may be observed and is naturally explained by external Compton reprocessing of the disk/BLR radiation by the jet.
• The SED can be adequately modeled via both purely leptonic and mixed lepto-hadronic scenarios, suggesting that the hadronic component is sub-dominant except in the X-ray and in the MeV bands.Our results predict that future missions in the MeV band, such as AMEGO-x and ASTROGAM (Caputo et al. 2022; De • The analysis of the simultaneous XMM-Newton and NuSTAR spectra during a comparatively low state, provides evidence (≳ 3σ) of a break in the X-ray band.If the break was to be intrinsic to the object, a pure leptonic model faces challenges in reproducing it.On the other hand, the proposed lepto-hadronic model shows a turnover of the spectrum in the X-ray band, that marks the kick-in of the hadronic component contribution.Based on our theoretical modeling, the SED is overall leptonic-dominated.Therefore, observations in lower activity states, as the one studied here, may offer better chances to pinpoint the hadronic fingerprint.
• The relatively high accretion regime of 5BZB J0630-2406 is supported also by the ratio of the γ-ray luminosity L γ (0.1 − 100 GeV) and the Eddington luminosity, L γ /L Edd ≃ 0.15.As further reprove of the intrinsic FSRQ nature, the model reveals an efficient accretion regime of 2 × 10 −4 .Besides, a relatively large fraction of the γ-ray luminosity (∼ 15% of L Edd ) is observed.Similarities can be found with the object TXS 0506+056 (Padovani et al. 2019b), where the localization of the γ emission region is thought to be, as in our case, on the edge of the BLR avoiding a significant γ − γ absorption and an efficient neutrino production.
• The neutrino emission predicted within the leptohadronic framework is at reach of the IceCube detector, close to the flux sensitivity.During the 7-yr integration span of the data used in Buson et al. (2022a), we expect to observe N events = 4.82 +5.18 −3.82 neutrinos.Assuming that the object maintains a constant neutrino rate over time, more high-energy neutrinos may be expected with increased instrument exposure.However, the uncertainties on the predicted numbers are large and long-term variability in the MWL light curve is clearly present at almost all observable frequencies.
• The contribution of 5BZB J0630−2406 to the astrophysical neutrino diffuse flux is expected to be of the order of ∼ 1%.
Based on the theoretical predictions presented here, the PeVatron blazar 5BZB J0630−2406 is capable of producing neutrinos in the IceCube energy range and can plausibly contribute to the anisotropy observed in the distribution of IceCube events of the hotspot IC J0630−2353 (see Fig. 1, Buson et al. 2022a).5BZB J0630−2406 is a high-power, radiatively efficient blazar.Other objects in the PeVatron blazar sample display similar characteristics, i.e.TXS 0506+056, PKS 1424+240 and 5BZB J0035+1515 (Buson et al. 2022a,b).At the current status it remains unclear whether this peculiar, relatively rare characteristic describes the persistent behavior of their engine and/or is linked to different environment properties, or may be tracing temporary physical changes, such as changes in the state of the accretion mode, as suggested for "changing-look blazars" (Peña-Herazo et al. 2021), or changes in the location of the dissipation region (Ghisellini et al. 2013).Future investigation of the PeVatron blazar sample will enable us with a broader understanding of the neutrino/blazar physical relation.

Figure 1 .
Figure 1.Cut-out region of the 7-yr IceCube L-value map (Aartsen et al. 2017c) centered at the position of the hotspot IC J0630−2353, displayed in celestial coordinates.The position of the associated blazar 5BZB J0630−2406 is highlighted by a green cross.The black lines indicate the 1 degree × 1 degree coordinate grid.

Figure 2 .
Figure2.Multi-wavelength light curve of 5BZB J0630−2406 from γ-rays, X-rays to the optical band.The top panel shows the γ-ray light curve with one-yr binning, available from the 4FGL-DR4 catalog.The second panel displays measurements from Swift-XRT, XMM Newton and Nustar used for the modeling.The third panel shows values taken by Swift-UVOT in the three optical filters V, B, and U.The fourth panel shows the V-band and g-band data collected by ASAS-SN and the bottom panel shows the KAIT Clear band measurement, the SMARTS R-band and the ZTF r-band and g-band.The X-ray and optical light curves are corrected for galactic extinction.The orange line highlights the time of quasi-simultaneous data used in the SED analysis.GROND data also used for modeling but not included in this panel are discussed in section 3.2.

Figure 3 .
Figure 3. 5BZB J0630−2406 NuSTAR FPMA (blue circles) and FPMB (red squares) light curves.Both datasets are fit with a constant, shown with a blue (red) dashed line.No trend is observed in the light curves, suggesting no significant intra-observation variability in the NuSTAR data.
Building the quasi-simultaneous broadband SEDIn this work, we are mostly interested in the modeling of the MWL SED during the time span of the IceCube obser-

Figure 4 .
Figure4.Schematic view of the model.1: Supermassive black hole, 2: Accretion disk defined by a luminosity L disk and a temperature T disk , 3: Dust torus defined by a luminosity L torus and a temperature T torus , 4: Broad line region (BLR) that scattered direct emission, 5: Relativistic jet (not modeled here) and 6: Moving emission region ("blob") situated at a distance R diss from the black hole.Direct (indirect) external photon fields are represented by continuous (dashed) arrows.Scheme not to scale.

Figure 5 .
Figure5.Spectral energy distribution from the purely leptonic model of 5BZB J0630−2406 in the observer frame.The dotted-dashed curves represent synchrotron and inverse Compton from pair production due to gamma-gamma absorption.The black solid squares represent the contemporaneous MWL data used for the modeling.Archival data are shown by gray solid dots.The red downward triangle represents the accretion disk limit from(Ghisellini et al. 2012).

Figure 7 .
Figure 7. Neutrino flux scaled according to the livetime showed in Aartsen et al. (2017c) derived for both mixed lepto-hadronic modeling of 5BZB J0630−2406.The uncertainties are computed assuming Poisson statistics and 3σ levels.The dashed gray line represents the 7-years IceCube sensitivity, assuming E −2 and δ = −24 • .06.

Figure 8 .
Figure 8. Evolution of the number of observed events by the Ice-Cube detector in time according to various string configurations, and assuming a constant flux.The uncertainties are evaluated assuming Poisson statistics (3σ levels).The blue dashed line represents the integrated livetime used in Buson et al. (2022a).

Table 1 .
Parameters used for the leptonic and the mixed leptohadronic models.