Broadband Multiwavelength Study of LHAASO-detected Active Galactic Nuclei

Recently, the Large High Altitude Air Shower Observatory (LHAASO) collaboration presented the first catalog of γ-ray sources using 508 days of LHAASO data, from 2021 March to 2022 September. This catalog contains four blazars and a possible LINER-type active galactic nucleus (AGN) counterpart. In this work, we establish averaged multiwavelength spectral energy distributions (SEDs) by combining data from the Fermi-Large Area Telescope, Swift, Zwicky Transient Facility, and Wide-field Infrared Survey Explorer (WISE) covering the same period as the LHAASO detection. In general, these five AGNs are found in low states at all wavelengths. To study the multiwavelength properties of these AGNs, several jet emission models, including the one-zone leptonic model, the one-zone leptonic and hadronuclear (pp) model, the one-zone proton-synchrotron model, and the spine-layer model, are applied to reproduce their averaged SEDs. We find that the one-zone leptonic model can reproduce most of the SEDs, except for the high-energy tail of the LHAASO spectra of Mrk 421 and Mrk 501. To improve the fitting, emission from pp interactions is favored in the framework of a one-zone model. The spine-layer model, which can be treated as a multizone scenario, can also provide good spectral fits. The influence of different extragalactic background light models on fitting a LHAASO energy spectrum is also discussed.


Introduction
Very-high-energy (VHE, 0.1 ∼ 100 TeV) γ-rays are one of the most important messengers for the investigation of the most extreme phenomena in the Universe.More than 90 extragalactic sources have been detected in the VHE band, the majority of which are jetted active galactic nuclei (AGNs; de Naurois 2021).These VHE AGNs include blazars with powerful jets pointing toward the observer, and radio galaxies, which are considered as the misaligned counterparts of blazars (Blandford & Rees 1978;Urry & Padovani 1995).In addition, other subclasses with GeV detection (e.g., narrow-line Seyfert 1 galaxies; Luashvili et al. 2023), or jets are potential VHE emitters as well.
Recently, the Large High Altitude Air Shower Observatory (LHAASO) collaboration presented the first LHAASO catalog of VHE γ-ray sources, in which four blazars and a possible AGN counterpart are reported (Cao et al. 2023).The four blazars, i.e., Mrk 421, Mrk 501, 1ES J1727 + 502, and 1ES 2344 + 514, are classified as the high-synchrotron-peaked (HSP; Abdo et al. 2010) type and are known extragalactic VHE emitters that have been extensively studied and included in the TeVCat7 website.The possible AGN counterpart is the LINER-type AGN, namely NGC 4278.Although the Fermi telescope does not detect GeV photons from NGC 4278,8 the parsec-scale jet discovered by radio observation is still a possible site for the acceleration of relativistic particles (Ly et al. 2004;Giroletti et al. 2005).The new LHAASO observations shed light on the VHE radiation mechanism of AGNs.
The physical origin of the VHE emission of jetted AGNs is complex and under debate.Due to the lack of strong external photon fields for HSP blazars, the most commonly used interpretation is the synchrotron self-Compton process (SSC; Bloom & Marscher 1996;Marscher & Travis 1996;Abdo et al. 2011a).Since the Klein-Nishina (KN) effect softens the SSC spectrum, the SSC model predicts a soft VHE spectra naturally.However, the intrinsic hard VHE spectra of some AGNs imply a different physical interpretation.Several models are proposed, such as the spine-layer model (Ghisellini et al. 2005;Acciari et al. 2020a), the proton-synchrotron model (Aharonian 2000;Mücke & Protheroe 2001;Mücke et al. 2003;Cerruti et al. 2015;Xue et al. 2023), and the ultra-high-energy cosmicray propagation model (Essey et al. 2011;Prosekin et al. 2012;Das et al. 2020Das et al. , 2022)).The reported minute-scale variability at VHE band also implies a multizone origin of the multiwavelength emission of the jet (e.g., Begelman et al. 2008;Liu et al. 2023).On the other hand, many associations between Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
In this work, to understand the radiation mechanism of these LHAASO-detected AGNs comprehensively, we build averaged multiwavelength spectral energy distributions (SEDs) by combining observations of Wide-field Infrared Survey Explorer (WISE) in the infrared band, Zwicky Transient Facility (ZTF) in the optical band, Swift in the X-ray band and the ultraviolet band, and Fermi Large Area Telescope (LAT) data in the γ-ray band.Several models are applied to fit SEDs, especially the VHE spectra.This paper is organized as follows: In Section 2, we describe multiwavelength observations and data reduction.The model description and fitting results are presented in Section 3. Finally, we end with discussions and conclusions in Section 4. The cosmological parameters H 0 = 70 km s −1 Mpc −1 , Ω 0 = 0.3, and Ω Λ = 0.7 are adopted.

Observations and Data Reduction
In this section, we present the multiwavelength observations of Mrk 421, Mrk 501, 1ES 1727 + 502, 1ES 2344 + 514, and NGC 4278 from 2021 March 5 to 2022 September 30 during the operation of the Water Cherenkov Detector Array (WCDA) of LHAASO and describe the process of data reduction.Detailed information on these five LHAASO AGNs is given in Table 1.

LHAASO
LHAASO consists of two VHE emission detector arrays, WCDA sensitive to γ-rays with energies between 100 GeV and 30 TeV and a 1.3 km 2 array (KM2A) sensitive to γ-rays with energies above 10 TeV (Ma et al. 2022).The VHE data were collected from the first LHAASO catalog.All five AGN sources were detected only by WCDA and not by KM2A.Then, the parameters of the power-law SED of the VHE spectra obtained by WCDA were given in Table 1 of Cao et al. (2023), and all of them were evaluated without correcting the spectra for extragalactic background light (EBL) absorption.The 95% statistic upper limits of the KM2A component were also given in the same table.

Fermi-LAT
The LAT on board the Fermi mission is a pair-conversion instrument that is sensitive to GeV emission (Atwood et al. 2009).Data are analyzed with the fermitools version 2.2.0.A binned maximum likelihood analysis is performed on a region of interest (ROI) with a radius of 10°centered at the R.A. and decl. of each source.The recommended event selections for data analysis are "FRONT+BACK" (evtype = 3) and evclass = 128.We apply a maximum zenith angle cut of =  z 90 zmax to reduce the effect of the Earth albedo background.The standard gtmktime filter selection with an expression of (DATA_QUAL > 0 && LAT_CONFIG == 1) is set.A source model is generated containing the position and spectral definition for all the point sources and diffuse emission from the 4FGL-DR3 catalog (Abdollahi et al. 2022) within 15°of the ROI center.The analysis includes the standard Galactic diffuse emission model (gll_iem_v07.fits) and the isotropic component (iso_-P8R3_SOURCE_V3_v1.txt), respectively.We bin the data in count maps with a scale of 0°.1 pixel −1 and set ten logarithmically spaced bins per decade in energy.An energy dispersion correction is made when event energies extending down to 100 MeV are taken into consideration.The spectral parameters of weak sources located within 10°of the center of the ROI are fixed during the maximum likelihood fitting.In a few cases, we fix or delete some sources to obtain a convergent fit.We divide this SED into six equal logarithmic energy bins in the 0.1-100 GeV, and an additional bin in the 100-800 GeV for these LHAASO sources.We built GeV lightcurves using about 8 day intervals between 0.1 and 100 GeV photons, shown in Figure 1.For the data points with poorly measured fluxes (where the likelihood test statistic, hereafter TS, < 10 or the nominal uncertainty of the flux is larger than half the flux itself), upper limits at the 95% confidence level are given.The TS and spectral index can be found in Table 1.

Swift-XRT
We make use of the Swift X-Ray Telescope (XRT) data products generator9 (xrt_prods) to obtain 0.3-10 keV X-ray light curves and spectra.Version 1.10 of the xrt_prodsmodule is released as part of swift toolsv3.0.This facility allows the creation of publication-ready X-ray light curves and spectra.Processing is performed using HEASOFTv6.29.Note.Columns from left to right: (1) the source name, (2) R.A., (3) decl., (4) the redshift of the source, (5) the SMBH mass in units of the solar mass, M e , (6) test statistic (Fermi-LAT), ( 7) spectral index (Fermi-LAT), ( 8) the type of AGNs.Instrumental artifacts such as pile up and the bad columns on the CCD are corrected as suggested by Evans et al. (2007Evans et al. ( , 2009)). 10hese spectra and X-ray light curves are produced by specifying the same covering times as the optical band data.Other settings have adopted default values from xrt_prods.Those obtained spectra are not the single observed spectrum but the average spectra over the entire considered time window whose timescale for the time binning is from MJD = 59278, to MJD = 59852, which are observed by photon-counting (PC) mode and windowed timing (WT) mode.The WT mode spectra are taken for the Mrk 421 and Mrk 501 because there are only a few short exposure observations in PC mode and have longer exposure observations on the same day in WT mode.The PC mode spectra are chosen for 1ES 1727 + 502, 1ES 2344 + 514, and NGC 4278 due to the reasons similar to the above or no observations in WT mode.After downloading those average spectra, we chose XSPEC (version 12.9) to fit them, and fit the spectra of the two modes separately.The specific fitting process is as follows.In order to obtain smaller flux errors, we apply the grpphacommand to rebin channels, setting a minimum number of groups greater than 29.The group min of NGC 4278 is equal to 4 due to insufficient photon counts.The power-law (po) model is often considered for fitting of X-ray spectrum.It is good to reference the logarithmic parabolic (logpar) model (Massaro et al. 2004).Thus, the models of both TBabs  2 and Figure 2.For NGC 4278, the model of TBabs * TBabs * cflux * po with index of 1.157 ± 0.327 is chosen to fit the spectrum while the model of TBabs * TBabs * cflux * logpar cannot be completely excluded.Also, the model of TBabs * TBabs * cflux * (po +bbody) does not significantly improve the fitting compared to the results from po or logpar.After the fitting is completed, we use the eeufspec command to convert them into unfolded spectra whose flux value can be transformed into νF ν of SED.The spectra absorption is corrected by multiplying the ratio of nonabsorbed and absorbed model values.For Mrk 421 and Mrk 501, no absorption correction was made because the absorption is weak when E is greater than 2 keV.

Swift-UVOT
The ultraviolet and optical data can be obtained by the Swift Ultraviolet and Optical Telescope (UVOT), which is equipped with broadband ultraviolet (UVW1, UVM2, and UVW2) and optical (V, B, and U) filters (Roming et al. 2005).Based on the HEASARC Archive Search Web,12 Swift-UVOT images from the five sources between MJD = 59278 and MJD = 59852 were retrieved.NGC 4278 only has one UVOT observation (observation ID: 03109562002), but it is not located in the detection window.Therefore, no UVOT images are available for this source.There is only data in the ultraviolet band for Mrk 421.The latest version of HEASoft 6.32.1 and calibration files CALDB version 20211108 are used during data processing.According to UVOT analysis threads, 13 we check whether level 2 images < > ( [ ] ) sw obsid u filter sk.img are correctly aligned to the World Coordinate System.The small-scale sensitivity check is performed by default by the software.The uvotimsum command is then used to sum extensions within an image, and the uvotsource command is used to perform aperture photometry with a circular source region of 5″ radii and a circular (annular) background region of 15″-40″ (inner) radii.The output results include the magnitude of the AB system, corrected count rate, and the flux density in mJy, etc.The corrected count rate is converted into magnitude of the AB system using the new AB zero-points (Breeveld et al. 2011), and magnitude of the AB system into flux density by considering zero-point flux density.The flux density is corrected for Galactic extinction.The specific process is as follows: a reddening coefficient of E(B − V ) is obtained from Schlafly & Finkbeiner (2011) with R V = 3.1.Then, the extinction value A V is calculated using the dust extinction laws; those of Fitzpatrick (1999) are chosen.Based on this extinction curve, we obtain the extinction value for the Swift ultraviolet and optical bands, and then multiply our flux by ( ) 0.4 Swift band to perform the extinction correction.

ZTF
The optical magnitudes in g, r, and i bands are collected from the 17th ZTF public data release14 (Masci et al. 2019).If the parameter catflags for a ZTF image has a value less than 32768 (i.e., does not contain bit 15), the photometry at that epoch is probably usable (Masci et al. 2019).Thus, in order to obtain good observation data, we require catflags score = 0 for other sources apart from Mrk 421, which does not have data with catflags score = 0.It should be noted that these data with catflags score = 4096 are chosen for Mrk 421.We convert the g, r, i magnitudes into fluxes following Xiong et al. (2020).In addition, the Galactic extinctions in the g, r, and i bands are corrected, and the NASA/IPAC Extragalactic Database15 provides extinction values (also see Schlafly & Finkbeiner 2011).
We construct the SED using the average magnitudes and average errors during the selected period.

WISE
The WISE (Wright et al. 2010) telescope has been operating a repetitive all-sky survey since 2010, except for a gap between 2011 and 2013.The WISE telescope visits each location every half a year and takes >10 exposures during 1 day.Although initially four filters were used, most of the time, only two filters, named W1 and W2, are used at the moment.The central wavelengths of the two filters are 3.4 and 4.6 μm.We collected the magnitudes of the five sources by point-spread functions fitting from the NASA/IPAC InfRared Science Archive (IRSA). 16Following Jiang et al. (2021), we selected magnitudes with good image quality (qi_fact > 0) and unaffected by charged particle hits (saa_sep > 0), scattered moon light (moon_masked <1), or artifacts (cc_flags = 0); and then binned the magnitudes every half a year since we did not detect any intraday variabilities.The magnitudes are in the Vega system; then, we can convert WISE Vega magnitudes to flux density units with F ν = F ν0 × 10 −0.4 m , where the zero magnitude flux density (F ν0 ) for the W1 and W2 bands is 309.5 and 171.8 Jy, respectively, with m being the calibrated WISE magnitude.

SED Modeling
These five LHAASO AGNs do not show a significant flare17 during the same observation period of LHAASO in all bands, as shown in Figure 1.Therefore, we use the averaged flux of each band to construct SEDs.As mentioned above, the first LHAASO catalog reports five AGNs including four HSP blazars and one LINER-type AGN.In the case of blazars, the multiwavelength emission, apart from the obvious thermal peaks in the infrared and optical bands, is undoubtedly from the jet.For the LINER-type AGN NGC 4278, the jet/radiatively inefficient accretion flow and the thin accretion disk may alternately dominate as the origin of the radiation, depending on the strength of its X-ray emission (Younes et al. 2010).In our data analysis, the discovered hard X-ray power-law spectrum favors the jet origin.In order to better understand the radiation mechanisms of these LHAASO-detected AGN, we consider four popular jet models to reproduce the SEDs: the one-zone SSC model, the one-zone SSC+pp model, the onezone proton-synchrotron model, and the spine-layer model.
The synchrotron radiation, the inverse Compton (IC) radiation, and the pp interactions radiation are calculated using the naima Python package (Zabalza 2015).We also consider the absorption of γ-ray photons due to the soft photons in the radiation zone (Xue et al. 2022) and the EBL (Domínguez et al. 2011) during propagation in intergalactic space.In addition, the energy of the absorbed γ-ray photons in the radiation region is converted to lower energies through the cascade process.The calculation of the cascade spectrum is applied as proposed in Böttcher et al. (2013).
In the infrared and optical bands of SEDs, most sources exhibit a clear hump, which differs significantly from the trend in the other bands.This is normally suggested as the emission from the host galaxy.We assume that the host galaxy is a 13 Gyr old elliptical galaxy for all of the AGNs (Raiteri et al. 2014) and use the SWIRE template18 (Polletta et al. 2007) to generate the spectrum of host galaxy in the fitting.The host galaxy contribution is based on the results of Polletta et al. (2007) for Mrk 501 and 1ES 2344 + 514.The flux of the host galaxy in the R band of ∼1 mJy is used for 1ES 1727 + 502 (Nilsson et al. 2007).For NGC 4278, the host galaxy contribution can be clearly distinguished from the continuous radiation components, so it can be obtained by fitting.The spectrum of infrared, optical, and UV data from Mrk 421 in Figure 3 has a power-law shape, indicating a weak contribution from the host galaxy.

One-zone SSC Model
The one-zone SSC model is the simplest and most commonly used model in the study of jet emission.In this paper, we assume a broken power-law injection electron density distribution.By taking into account the radiative cooling and the escape of the electrons, the steady-state electron density distribution can be calculated with (Xue et al. 2019a) and g e,max are the minimum and maximum electron Lorentz factors of the distribution, γ e,b is the break electron Lorentz factor, s e,1 and s e,2 are the low-energy and the high-energy indexes of the broken power-law spectrum, L e inj is the electron injection luminosity, R is the radius of radiation zone, m e is the rest mass of the electron, c is the speed of light, t esc = R/c is the escape timescale, t cool = 3m e c/(4σ T γ e (u B +f KN u ph )) is the electron cooling timescale, σ T is the Thomson scattering cross section, f KN is the factor accounting for KN effects (Moderski et al. 2005) is the energy density of the magnetic field, B is the magnetic field strength, and u ph is the energy density of the soft photons. 19The observed emission will be Doppler boosted by a factor δ 4 , where 1 is the Doppler factor, Γ is the the bulk Lorentz factor, β Γ c is the velocity of the jet, and θ is the viewing angle of the jet.
Current observational data do not provide good constraints on all parameters in the modeling.Therefore, we fix some of the less sensitive parameters in fitting to reduce the number of free parameters.All the fitting parameters can be found in Table 3.The relativistic jet of the blazars is close to the line of sight of the observer, so we set the viewing angles for Mrk 421, Mrk 501, 1ES 1727 + 502, and 1ES 2344 + 514 to 1°.8 uniformly.Giroletti et al. (2005) suggest that the viewing angle of NGC 4278 is uncertain, and their study reports that it could have a small viewing angle 4 , alternatively a large viewing angle.We therefore divide it into two cases, one with the same viewing angle (θ = 1°.8,hence NGC 4278 a ) as the blazar sources and the other with a larger viewing angle (θ = 30°, hence NGC 4278 b ).
The fitting results of the one-zone SSC model are shown in Figure 3.In the case of four blazars, it can be found that the low-energy component of the SED can be reproduced by the superposition of host galaxy emission and electron synchrotron radiation, except that the model slightly underestimates the UV data.This may be due to the use of nonsimultaneous data.For example, there is an absence of Swift-UVOT observations from MJD 59500 to MJD 59550 for Mrk 421, in which period its flux for the optical band is at its lowest state during the entire observation period, and may overestimate its flux in the UV band.The high-energy component can be fitted very well by the SSC model.However, the high-energy tail of the VHE The light blue data points are infrared data from WISE, the dark blue data points are optical data from ZTF, the orange and green data points are X-ray data from Swift-XRT's PC mode and WT mode, respectively, and the purple data points are γ-ray data from Fermi-LAT, the red strap shows the observation of WCDA, and the red upper limit point is from KM2A.The gray data points are historical data from the SSDC SED Builder Tool (https://tools.ssdc.asi.it/) of the Italian Space Agency (Stratta et al. 2011).Notes.The minimum electron Lorentz factor g e,min is set to 1 × 10 2 because it is insensitive in fitting.The "L" sign indicates that the parameters do not exist in the one-zone SSC or SSC+pp scenario.
, where m is the number of the observational data points, n is the number of free parameters in the fitting model, ŷi is the expected value from the model, y i is the observed data, and σ i is the standard deviation for each data point.For the WCDA spectrum with only the power-law bow-tie, we divide it into an average of 30 bins in the logarithmic energy space to calculate the chi-square value.d The subscript "WCDA" indicates that we use only WCDA data to calculate the chi-square value.
e We consider two radiation zones here.The initial zone maintains the same leptonic parameters as the SSC+pp model, while the second zone primarily considers proton-synchrotron radiation.The parameters for the latter zone are shown here.For more detailed discussion, please refer to Section 3.3.f Although the proton injection luminosity is much smaller than the Eddington luminosity, the magnetic power is comparable to the Eddington luminosity in this case.g We present two strategies for modeling with the spine-layer model.Here, we consider using the EC process to fit the high-energy hump independently and show its parameters.A more detailed discussion can be found in Section 3.4. .This means that the IC radiation spectrum is steeper above ∼0.2TeV because of KN effects, as shown in Figure 3. Therefore, the high-energy tail of LHAASO spectra cannot be fitted with the one-zone SSC model, unless very extreme parameters are considered.
In the case of NGC 4278, due to the lack of GeV γ-ray data, the fitting parameters have a larger space to choose from.Nevertheless, in order to explain the spectra in both the X-ray and VHE bands simultaneously, extreme parameters are required.For example, the model requires a very hard lowenergy slope (s e,1 = 1) or a very large minimum electron Lorentz factor (approaching γ e,b ) to explain the very hard X-ray spectra.If we consider that its X-ray radiation is produced by the electron-synchrotron process, a large γ e,b (around 10 6 ) is required.In addition, the critical energy E KN obs is about 0.7 GeV for the case of NGC 4278 a and 7 MeV for the case of NGC 4278 b , which requires the low-energy slope s e,1 to be close to 1 and the high-energy slope s e,2 3 in order to counteract the impact of the KN limit and fit the VHE spectra.For NGC 4278 b , it is also necessary to set g e,max to 5 × 10 7 .

One-zone SSC+pp Model
The pp model has a potential to produce the observed TeV spectra of blazars without exceeding the Eddington luminosity, which is difficult to avoid in the pγ model (Xue et al. 2019b).Our recent study (Xue et al. 2022) shows that the pp interactions in the jet have the potential to generate VHE emission that can be deteced by LHAASO.Therefore, we incorporate the pp interactions into the one-zone SSC model to reproduce the SEDs.
In the pp modeling, we assume a power-law injection proton density distribution.By taking into account the radiative cooling and the escape of the protons, then, the steady-state proton density distribution can be calculated with p is the injection proton density distribution, γ p is the proton Lorentz factors in the range of g p,min to g p,max , s p is the slope of the power-law spectrum, L p inj is the proton injection luminosity, m p is the rest mass of the proton, and t cool (γ p ) is the cooling timescale of the proton.More specifically, t cool (γ p ) is dominated by the pp interactions in the SSC+pp scenario and can be approximated by H , where K pp ≈ 0.5 is the inelasticity coefficient, n H is the number density of cold protons in the jet, .
To maximize the efficiency of the pp interaction within a reasonable parameter range, analytical calculations suggest that the power of the cold protons in the jet should be set as half the Eddington luminosity (Li et al. 2022;Xue et al. 2022).Then, we can get the number density of cold protons , where L Edd = 1.26 × 10 38 M BH /M e is the Eddington luminosity, and M BH is the SMBH mass.We may estimate the maximum proton energy by equating the acceleration timescale with the escape timescale (Xue et al. 2019a).Then, the maximum proton energy can be calculated by where e is the elementary charge, and α is the factor representing the deviation from the highest acceleration rate.
We employ α = 1100, which corresponds to the shock speed measured in the upstream frame of ∼0.07c in the situation of shock acceleration (Rieger et al. 2007).We set the minimum proton energy g = 1 p,min , and the slope of the power-law spectrum s p = 1.5.These two parameters have no effect on the fitting results, only on the required proton injection luminosity.Finally, the only free parameter of the pp model remains the proton injection luminosity L p inj , which is shown in Table 3.The fitting results can be found in Figure 4.The radiation produced by the pp interactions fits perfectly the high-energy tail of the LHAASO spectrum of Mrk 421 and Mrk 501, which previously could not be explained by the one-zone SSC model.The dotted curve in Figure 4 shows the neutrino flux produced by the pp interactions, which should be comparable to the photon flux produced at the same time (as shown in two cases of NGC 4278).The sudden drop of photon flux in four cases of blazars is due to the absorption of photon-photon interactions occurring in the jet and in the intergalactic propagation.The fitting parameters in Table 3 show that Mrk 421, Mrk 501, 1ES 2344 + 514, NGC 4278 a , and NGC 4278 b can be fitted without exceeding the Eddington luminosity, and a quite low proton injection luminosity is required in the cases of 1ES 2344 + 514 and NGC 4278 a .The proton injection luminosity used in the 20 It should be noted that the VHE spectra of TeV AGNs are sometimes fitted as log-parabolas, whereas the first LHAASO catalog only reports powerlaw SEDs.
fitting of 1ES 1727 + 502 is 60 times that of the Eddington luminosity, because of its low SMBH mass (5.62 × 10 7 M e ).For comparison, the fiducial SMBH mass of the BL Lac objects is 10 8.5−9 M e (Shaw et al. 2013;Xiao et al. 2022).Moreover, based on the fitting results, it can be found that the cascade process makes a negligible contribution to the energy spectrum in the one-zone SSC+pp model.
From the chi-squared test results in Table 3, it can be seen that, for the cases of Mrk 421, Mrk 501, NGC 4278 a , and NGC 4278 b , the goodness of fitting has been significantly improved after the introduction of the pp model (especially considering only the chi-squared values of the WCDA data).While, for the two cases of 1ES 1727 + 502 and 1ES 2344 + 514, the introduction of the pp model in the fitting has no significant advantage.
It is necessary to evaluate the possible neutrino emission when high-energy protons are introduced into the model.Therefore, we calculate the neutrino flux that could be produced in the process of the pp interactions.The timeintegrated neutrino sources' searches with 10 yr of IceCube data collected between 2008 and 2018 report that the best-fit number of astrophysical neutrino events for Mrk 421 is 2.1, with a local pre-trial p-value of 0.42, and the number of astrophysical neutrino events for Mrk 501 is 10.3, with a pvalue of 0.25 (Aartsen et al. 2020).It should be noted that neutrino events from Mrk 421 and Mrk 501 have not yet been detected by IceCube.These best-fitting numbers of astrophysical neutrino events both correspond to very large p-values, which means that this result can only be treated as a rough upper limit (e.g., Abe et al. 2023).The neutrino event rate can be obtained by , where n E ,min and n E ,max are the lower and upper bounds of the neutrino energy, respectively, A eff is the IceCube point-source effective area for (anti)muon neutrinos (Carver 2019), δ decl is the decl., and f n ( ) E is the differential neutrino energy flux.Then, the expected neutrino event rates21 from the sources are 0.15 events yr -1 for 1ES 1727 + 502, 0.21 events yr -1 for 1ES 2344 + 514, 0.11 events yr -1 for NGC 4278 a , 0.10 events yr -1 for NGC 4278 b , 0.47 events yr -1 for Mrk 421, and 1.06 events yr -1 for Mrk 501, which are slightly exceed the upper limits for Mrk 421 and Mrk 501.Thus, although the SSC+pp model reproduces the SEDs best and has the smallest χ 2 /d.o. f for Mrk 421 and Mrk 501, it remains to be investigated whether pp interactions can reasonably explain the high-energy tail of the LHAASO spectra, after obtaining the simultaneous SEDs or the variability in the VHE band.

One-zone Proton-synchrotron Model
The proton-synchrotron emission in the framework of onezone model is often suggested as a possible interpretation of the high-energy component of HSP AGNs (e.g., Aharonian 2000;Böttcher et al. 2013;Acciari et al. 2020a), although extreme physical parameters, such as a super-Eddington jet power and a strong magnetic field, are usually introduced (see, Cerruti et al. 2015;Petropoulou & Dermer 2016;Xue et al. 2023).In this subsection, before fitting SEDs of these five LHASSO AGNs, we search the proton-synchrotron modeling parameter space with an analytical method proposed in our recent work (Xue et al. 2023).In this method, the parameter space is limited by three constraints, which are the total jet power (dominated by the injection power of relativistic protons and the power carried in magnetic field) does not exceed the Eddington luminosity, relativistic protons can be accelerated to the required maximum energy, and the emitting region is transparent to VHE photons, respectively.In addition, observation results suggest that the magnetic field in the inner jet of AGNs is typically lower than 10 G (O'Sullivan & Gabuzda 2009; Pushkarev et al. 2012;Hodgson et al. 2017;Kim et al. 2022), and the bulk Lorentz factor of the jet is lower than 30 (e.g., Hovatta et al. 2009).If a reasonable parameter space that satisfies observational constraints can be found (i.e., B  10 G and Γ  30), we fit their SEDs.
There are two strategies to fit the LHAASO spectra with proton-synchrotron emission.The first one is to use protonsynchrotron emission to account for the entire high-energy component.The second one is to use proton-synchrotron emission to fit the high-energy tail of the LHAASO spectra, while the rest of the high-energy component is still attributed to the leptonic SSC emission.For the first strategy, the index s p of injected proton energy distribution can be obtained by the photon index Γ index of Fermi-LAT spectrum, i.e., s p = 2Γ index − 1. for NGC 4278.Xue et al. 2023 find that the reasonable parameter space might be found if considering a relative large blob radius (R  10 16 cm), and the size of parameter space is inversely proportional to R. So, here, we only check if the parameter space can be found when R = 10 16 cm.For the second strategy, we default = E peak syn 14 TeV, and s p = 2 for all five AGNs, since there is no constraint on s p .Since the default peak energy is quite large, it is only necessary to check if the parameter space exists in a large emitting region.Here, we set R = 10 17 cm.
The results of the parameter space scans are shown in Figure 5.In the first strategy (the upper six panels), it can be seen that no valid parameter space is found for all five AGNs (SED fitting with strong magnetic fields is given in Appendix A).Among them, the super-Eddington jet power is needed for four blazars because soft proton indexes are suggested by the Fermi-LAT spectra.For NGC 4278, there is a large parameter space to get the sub-Eddington jet power; however, this space is in conflict with the Hillas condition (black curves with arrows).In the second strategy (the lower six panels), it can be seen that only 1ES 2344 + 514 can find a reasonable parameter space.However, with R = 10 17 cm, if we set B = 7 G, and Γ = 10, the energy density of low-energy component U syn ≈ 8.2 × 10 −7 erg cm −3 would be much lower than that of magnetic field U B ≈ 1.9 erg cm −3 .Therefore, in the framework of the one-zone model, it is impossible to fit the GeV data of 1ES 2344 + 514 with SSC emission when using the proton-synchrotron emission to explain the LHAASO spectrum.A second emitting region has to be introduced.In Figure 6, we show that the SED of 1ES 2344 + 514, including the LHAASO spectrum, can be explained as a superposition of leptonic emission from first emitting zone and protonsynchrotron emission from second emitting zone.The leptonic emission from the first emitting zone is the same as that obtained in Section 3.1.For the second emitting zone, we set R = 10 17 cm, B = 9 G, and Γ = 12 as indicated by the obtained parameter space.
Overall, the one-zone proton-synchrotron model seems to be difficult to interpret the currently observed LHAASO spectrum within a reasonable parameter space (i.e., B  10 G and Γ  30).If we introduce a second emitting zone, the proton-synchrotron emission is only applicable to 1ES 2344 + 514.However, the fitting results of the two-zone proton-synchrotron model do not show any advantages over the one-zone SSC model, in terms of either the fitting results or chi-squared test.As shown in Appendix A, the chi-square test shows that the one-zone proton-synchrotron model has no advantage in the Figure 5.The ratio of L jet /L Edd in the Γ − B diagram for the one-zone proton-synchrotron model.The upper six panels show the results that applying protonsynchrotron emission to explain the entire high-energy component, and the lower six panels show the results that applying proton-synchrotron emission to explain the high-energy tail of LHAASO spectra.The black curves with arrows represent the parameter space that satisfied the Hillas condition.The vertical blue curves with arrows show the lower limit of Γ that allows the escape of maximum energy γ-ray photons.The vertical and horizontal purple curves show the space that B  10 G and Γ  30.The white dashed contours denote specific values of log (L jet /L Edd ) associated with the color bar.
fitting of the SEDs over the one-zone SSC model after removing the parameter constraints.

Spine-layer Model
The observed limb-brightening at the parsec (e,g., Giroletti et al. 2004Giroletti et al. , 2006;;Piner et al. 2010) and the kiloparsec scales (e.g., Owen et al. 1989;Laing et al. 2011) suggests that the jet could be structured with a fast spine surrounded by a slower layer.Based on this observation, the spine-layer (or structured) jet model is proposed and applied to account for the rapidly variable VHE emission (e.g., Ghisellini et al. 2005) and to reproduce the quiescent state SED (Tavecchio & Ghisellini 2016).In this subsection, we apply the spine-layer model to fit the SEDs.
This model consists of two components, a relatively small cylinder that is the spine (denoted by the subscript "s") and another hollow cylinder wraps around the spine as the layer (denoted by the subscript "l").Similar to the conventional twozone model, this model also requires two sets of parameters.The difference is that the spine and the layer influence each other, and there is a relationship between these two sets of parameters.We basically use the same settings as the one-zone SSC model for each component.There are three differences from before: 1.The radiation zone changes from spherical to a cylinder (spine) or hollow cylinder (layer), so the radius of the radiation zone R is changed to the cross-sectional radius R c , and we add a parameter of the length of the cylinder L. In addition, all calculations in relation to the shape of the radiation zone must also be replaced.For example, the volume of a sphere in Equation (1) must be replaced by that for a cylinder or hollow cylinder.2. We set R l = 1.2R s to reduce the number of free parameters, which follows Ghisellini et al. (2005).The spectral indexes of the electron spectrum in the layer are set to be the same as that in the spine, because they cannot be constrained in fitting.

Photons produced in one component can enter another
component and, as soft photons, enhance the IC emission in both components.And this process of scattering the soft photons coming from another component is commonly known as the external Compton (EC) process.The energy density of the soft photons from another component is calculated as suggested in Ghisellini et al. (2005).
Finally, there are 15 free parameters, which can be found in Table 3.In this spine-layer model, we consider three strategies for modeling: 1. We consider using the EC process to fit the high-energy hump independently.As discussed in Section 3.1, it is difficult to reproduce the entire VHE data using the onezone SSC model due to the KN effect.In the spine-layer model, if soft photons come from another component, Equation ( 3) is rewritten as s l s l is the relative Lorentz factor between the spine and the layer, and β s and β l are velocities for the spine and the layer, respectively.In the EC process, the soft photons provided by another component would not be constrained by the fitting of the low-energy hump as in the SSC process.Therefore, the KN effect could be weakened as long as another component can provide low-energy soft photons, as shown in Equation ( 6).This indicates that the EC process might improve the fitting result of the SSC process on VHE observations.We then fit the SEDs where all the observed radiation is produced in one component, and the other component provides only soft photons.The fitting results are shown in Figure 7.It can be seen that the EC process with slight/no KN effect still cannot effectively reproduce the high-energy tail of LHAASO spectrum of blazars.This is because the radiation spectrum produced by the IC process is not in the standard power-law form.The morphology of the spectrum produced by the IC process approaches a smooth curve because it is influenced by both electron and soft photon spectra, the KN effect, and the absorption of photon-photon interactions.This is evident in the fitting result of Mrk 421, where the EC radiation spectrum is curved compared to the power-law spectrum observed by LHAASO.In the case of NGC 4278 a , the multiwavelength radiation is explained by emission produced in the spine, while for NGC 4278 b it is explained by an emission from the layer, because the relativistic beaming effect reduces the flux in the case with a larger viewing angle and a larger bulk Lorentz factor.The Doppler factor δ = 0.30 is small for the spine of NGC 4278 b .2. The superposition of radiation from two components seems to be a plausible strategy to explain the VHE spectra, as it has minimal parameter constraints.The fitting results are shown in Figure 8.In the cases of four blazars, the X-ray spectra are explained by synchrotron emission produced in the spine, and the γ-ray radiation is from the superposition of emission from the spine and the layer.As shown in Figure 8, the EC emission spectrum produced in the spine rapidly decreases near TeV or sub-TeV due to the KN effect.The EC radiation spectrum produced in the layer can be extended to higher energies, because a larger break electron Lorentz factor is set in the layer.Although the spectral index of the radiation spectrum for each component is different from the observed spectral index in the VHE band due to the influence of the KN effect, the VHE spectrum can still be reproduced by superimposing the radiation of the two regions.In the two cases of NGC 4278, a very hard electron spectrum is still required in fitting.The lowenergy and the high-energy humps are both explained by the radiation superposition from two components.3. We consider fitting the SED with a superposition of multiradiation processes from one emitting region.To be specific, the low-energy hump is fitted by synchrotron radiation, and the high-energy hump is fitted by the radiation superposition of the SSC and the EC processes.
To reproduce the high-energy hump by the radiation superposition of the SSC and the EC processes, we consider a scenario similar to that shown in Figure 6.In this strategy, the peak energy of SSC or EC radiation is required to reach ∼14 TeV.If the high-energy hump originates from the IC process, the threshold peak energy The characteristic photon energy in the observer frame The first fitting strategy shows no advantage over the onezone SSC model, either from the fitting results or from the chisquare results.The second fitting strategy is not recommended either.Despite its seemingly superior reproduction of SEDs to the naked eye, particularly for the LHAASO spectra of Mrk 421 and Mrk 501, it often results in larger χ 2 /dof.This is primarily due to the fact that it employs nearly twice as many fitting parameters (n) as the one-zone SSC model.When considering different models to fit the same SED, the number of observed data points (m) remains constant.According to the formula in footnote c of Table 3, an increase in the number of fitting parameters (n) leads to a decrease in the denominator (m − n), which ultimately results in an increase in χ 2 /dof.

Can Emission from pγ Interactions Interpret the LHAASO Spectra?
As shown in Section 3.1, the one-zone SSC model can fit most of the multiwavelength spectra, except for the highenergy tail of the LHAASO spectra.To obtain a better fit, we comprehensively test the contributions from pp interactions, proton-synchrotron emission, and the spine-layer model.On the other hand, since various soft photon fields exist in AGNs environment, many works dedicate themselves to studying the electromagnetic and neutrino emissions from pγ interactions.As suggested by many recent studies (e.g., Sahu et al. 2021a;Alfaro et al. 2022;Sahu & Valadez Polanco 2022), here, we analytically discuss if emission from π 0 decay in the pγ interactions can improve the fitting of LHAASO spectra.
Using the δ− approximation, the relation between the energy of π 0 decay VHE photons E VHE obs and the energy of target photons E tar obs both in the observer's frame can be obtained, when considering the peak cross section of photopion interactions due to the !+ (1232) resonance.By taking L Edd as the maximum proton injection luminosity, 10 −28 cm 2 as the photopion cross section weighted by inelasticity, the lower limit of flux of the target photons can be estimated by As shown in Figure 3, the model-predicted fluxes at ∼1 MeV for these LHAASO AGNs are ∼10 −12 erg s −1 cm −2 .If considering all emission processes occur in one single region, it can be seen that, even when the emission from pγ interactions is only used to account for the high-energy tail of the LHAASO spectra, the model-predicted flux is still 2 orders of magnitude lower than the flux required.If one attempts to use the emission from pγ interactions to account for the whole LHAASO spectrum, the required flux of target photons would be more than 3 orders of magnitude higher than the model-predicted flux.Even taking into account that the pγ interactions and the leptonic processes may occur in different regions, i.e., a multizone case, the required flux at MeV band (> ∼ 10 −10 erg s −1 cm −2 ) still far exceeds all the existing AGN observations.Therefore, using the emission from the pγ interactions to interpret the LHAASO spectra can be confidently ruled out.

The Influence of Different EBL Models
The observed γ-ray spectra of extragalactic sources, especially in the VHE band, are softened by the interactions of the γ-ray photons with the EBL.The energy spectrum of EBL is difficult to obtain by direct observation, so many researchers have used various methods to estimate it.To evaluate the influence of different EBL models on the fitting results, we apply five EBL models to Mrk 421, showing the optical depths (left panel in Figure 9) and the intrinsic VHE spectra (middle panel in Figure 9), respectively.
The energy of γ-ray opacity equal to unity varies in different EBL models.Within these models, the energy is focused between the 7 and 10 TeV range for Mrk 421.The EBL model (green line in Figure 3), used for the calculation in Section 3, shows a relatively moderate optical depth in the energy range of the WCDA.Furthermore, various EBL models are also applied to calculate the corresponding intrinsic VHE spectra, based on the observed spectrum given by WCDA.The middle panel of Figure 9 shows that all models indicate the presence of a new component beyond ∼10 TeV, which is consistent with our fitting results.From Figure 9, it can be seen that all EBL models have an equivalent effect on the VHE spectra below 4 TeV.Thus, the intrinsic VHE spectrum can be estimated by extrapolating the 1-4 TeV spectrum up to higher energies, if we assume that the VHE radiation comes from one single component.The right panel of Figure 9 displays the hypothetical extended intrinsic spectrum (represented by the black line) and the expected observed spectrum after absorption.It appears challenging to reproduce the entire VHE spectrum from the radiation of a single component, even without considering the KN effect in the one-zone SSC model.
The best-fit power-law bow-tie of the LHAASO data cannot be broken into energy bins to evaluate how the EBL impacts it over the entire energy range, so corrections other than a scaling factor would result in a break in the power law, which are not captured by the best fit.The right panel of Figure 9 shows that the power-law spectrum is unlikely to be broken before 7 TeV by the influence of absorption from the EBL model of Domínguez et al. (2011), which is used in the calculation of Section 3. From Figure 3, it can be found that the energy spectrum predicted by the one-zone SSC model already deviates from the LHAASO bow-tie at about 2 TeV, and the model-predicted flux is only about half the observed flux at 7 TeV.This suggests that there are factors other than EBL absorption that cause the power-law spectrum to bend, which is consistent with the conclusion in Section 3.1.In the Section 4.1, the theoretical analysis predicts that, if the pγ interactions are used to explain the observation of the LHAASO data at 14 TeV, it is necessary to increase the model-predicted flux of 1 MeV by more than 100 times.From the left panel of Figure 9, it can be seen that the optical depths of different EBL models at 14 TeV range from 1.2 to 1.7, and the correction factor for the flux is between 0.30 and 0.18.Thus, even if taking into account the effect of EBL absorption on the LHAASO best-fit power-law bow-tie, it is necessary to increase the 1 MeV flux by more than an order of magnitude.The above discussion indicates that applying different EBL models does not affect our conclusions.In addition, the observational data from the γ-ray telescope can help to constrain the absorption optical depth induced by EBL and to constrain the EBL model (e,g., Aharonian et al. 2007;Fermi-LAT Collaboration et al. 2018;Abeysekara et al. 2019;Acciari et al. 2019).However, this would require a more abundant or simultaneous set of observation data.As it is beyond the scope of this paper, we will not discuss it further.

Comparison with Previous TeV AGNs' Studies
The full broadband SEDs modeling has been the main tool for the blazar study.The VHE γ-ray observation from Imaging Air Cherenkov Telescopes provides a strong constraint to the jet models.The four LHAASO blazars are all known VHE emitters, of which Mrk 421, Mrk 501, and 1ES 2344 + 514 are the earliest detected extragalactic VHE γ-ray sources (Punch et al. 1992;Quinn et al. 1996;Chadwick et al. 1999).The VHE γ-ray of 1ES 1727 + 502 has also been observed for more than 10 yr (Aleksić et al. 2011(Aleksić et al. , 2014)).Therefore, these four blazars have been extensively studied with numerous multiwavelength SEDs from different periods.
Most previous studies on these four blazars have shown that the one-zone SSC model can reasonably reproduce the SEDs (e.g., Albert et al. 2007Albert et al. , 2022;;Anderhub et al. 2009;Tavecchio et al. 2010;Abdo et al. 2011a;Acciari et al. 2011aAcciari et al. , 2011bAcciari et al. , 2011c;;Bartoli et al. 2011;Aleksić et al. 2014;Furniss et al. 2015;Bartoli et al. 2016;Prince et al. 2022).In contrast, some studies suggest that the one-zone SSC model cannot fit the SEDs because it underestimates the TeV γ-ray flux, e.g., the high-state TeV flux of Mrk 421 observed by the Whipple Observatory between 2006 April 16 and 20 (Błażejowski et al. 2005) (Finke et al. 2008).These diverse observations and fitting results suggest that the realized radiation process from these sources may be very complex, and more observations from different time periods and energy bands are key to further research.
In addition to the one-zone SSC model, other models are also used to fit the SEDs of these sources.Aleksić et al. (2013) found that both one-zone and two-zone SSC models can well reproduce the SED observed by 1ES 2344 + 514 in late 2008.Abe et al. (2023) suggested that the observed different patterns of variability of Mrk 501 would naturally be expected from the two-zone model.Abdo et al. (2011b)   SED of Mrk 421 with a VHE flare observed on 2017 February 4 is better reproduced by a two-zone leptonic model than by a one-zone leptonic model.Sahu et al. (2020Sahu et al. ( , 2021b)), Sahu & Valadez Polanco (2022) found that the one-zone photohadronic model is inadequate to explain the multi-TeV flaring events from the transient extreme HSP-like sources of Mrk 501, Mrk 421, and 1ES 2344 + 514, and they proposed a two-zone photohadronic model as an effective methodology.Manzoor et al. (2023) suggested that an additional emission mechanism other than the SSC process is required to explain the TeV observations of Mrk 421 by MAGIC in 2013 February, because its VHE spectra are remarkably harder than the X-ray spectra.Hu et al. (2023) included that the SED observed from Mrk 421 in 2013 January, in particular the hard X-ray excess, could have been generated as a result of the two-injection scenario.
These sources exhibit a number of observational features, particularly in the flaring state, that cannot be explained by the one-zone SSC model.Higher sensitivity observations at higherenergy bands will hopefully verify the above assumptions and models, and LHAASO has great potential in this regard.In this work, we collect multiwavelength data from five LHAASO AGNs during the same observation period as LHAASO.Based on theoretical and fitting analysis, we suggest that the one-zone SSC model is capable of reproducing most of the SED, with the exception of the VHE tail in the cases of Mrk 421 and Mrk 501.The inability of the VHE tail is mainly due to the collective effect of the KN effect, the EBL absorption, and the parameter constraints for other bands observations.This is well demonstrated in the case of NGC 4278, which is very close to us and has almost no EBL absorption.In addition, its multiwavelength data have very weak parameter constraints.Therefore, when we consider more extreme parameters, the one-zone SSC model can reproduce its SED, especially the LHAASO spectrum.We suggest that the high-energy tail of the LHAASO data of Mrk 421 and Mrk 501 cannot be fitted with the one-zone SSC model, unless very extreme parameters are considered.This is similar to the conclusion of Katarzyński et al. (2005), which suggests that the Thomson scattering into VHE photon energies requires unacceptably large Doppler factors.
To reproduce the SEDs of LHAASO AGNs, we apply the pp model, the proton-synchrotron model, and the spine-layer model.The results of fitting and the chi-square test suggest that the one-zone model, upon incorporating pp interactions, effectively accounts for all observations in the SEDs, especially the tail of VHE observation.In addition, a multizone model is also feasible if we consider the superposition of radiation generated by different regions to explain VHE observations, as demonstrated in the spine-layer model presented in Section 3.4.Despite its seemingly superior reproduction of SEDs to the naked eye, particularly for the LHAASO spectra of Mrk 421 and Mrk 501, it often results in larger χ 2 /d.o. f.
Our analysis results indicate that the proton-synchrotron model and the pγ model are difficult to explain the SEDs without considering very extreme parameters.Of all the sources, only the SEDs of 1ES 2344 + 514 can be reproduced using the two-zone proton-synchrotron model.A very large magnetic field (>10 G) must be introduced to fit the SEDs of the other LHAASO AGNs, whether in one-zone or two-zone proton-synchrotron models.The low interaction efficiency of pγ model, brought about by the lack of suitable soft photon fields, prevents it from reproducing the SEDs within reasonable parameters.
NGC 4278 is the most possible association with 1LHAASO J1219 + 2915.Moreover, it is also found to be positionally consistent with the γ-ray transient source 1FLT J1219 + 2907 detected by Fermi-LAT (Baldini et al. 2021).Therefore, although Fermi-LAT did not detect high-energy radiation from NGC 4278 during the 500 day period of the LHAASO detection, it remains a candidate with great potential for the VHE source.Unfortunately, the data reduction and SED fitting in this paper do not allow us to determine further whether the VHE radiation comes from NGC 4278.

Outlook
Our results suggest that VHE observations are crucial to constrain the jet model.As mentioned in Section 3.1, the SSC process of HSP enters the KN regime in the VHE band.Detailed observations in the VHE band can verify or rule out the origin of the one-zone SSC model more precisely.Furthermore, Figure 9 shows that different EBL absorption models have a significant different influence on VHE observation beyond 7 TeV.Therefore, by conducting further observations of extragalactic sources exceeding 7 TeV, we can constrain the EBL model better.This method has already been extensively applied in other γ-ray telescopes (e.g., Fermi-LAT Collaboration et al. 2018;Abeysekara et al. 2019;Acciari et al. 2019).
Multiwavelength variability can provide a different perspective to study the emission origin.For example, long-term monitoring is carried out for Mrk 421, as it is one of the closest BL Lac objects.Its VHE variability displays a highly complex behavior.Most observations have found a strong correlation between flares in the VHE band and the X-ray band (Fossati et al. 2008;Acciari et al. 2021;Arbet-Engels et al. 2021;MAGIC Collaboration et al. 2021).Some observations have reported that variations in the VHE band correlated with X-rays, but not with the optical (Giebels et al. 2007) and the other bands (Aleksić et al. 2015).Some variability studies indicate that the correlation between the X-ray band and the VHE band shows different behavior (Acciari et al. 2020b), andAbeysekara et al. (2020) find that the flux relationship changes from linear to quadratic, to no correlation, and to anticorrelation over the decline epochs.Błażejowski et al. (2005) report the inconsistency of X-ray band and VHE-band flare times.Taken together, these phenomena are difficult to explain using the one-zone SSC model.Therefore, the observations of variability in VHE band and the corresponding simultaneously SEDs are very important to investigate the radiation mechanisms and the physical properties of blazars.2018,2022), naima (Zabalza 2015).

Appendix Proton-synchrotron Modeling with Strong Magnetic Fields
Some studies show that the high-energy hump of SEDs can be fitted by the proton-synchrotron process with a strong magnetic field (Cerruti et al. 2015;Wang et al. 2023), although such a large magnetic field strength contradicts current observations.In the following, we will present the fitting results of the proton-synchrotron model, with a strong magnetic field.In this scenario, the leptonic modeling follows that given in Section 3.1, and the hadronic modeling basically follows that given in Section 3.2.There are five differences from before: 1.The power-law proton spectrum cannot fit the observation of LHAASO, and therefore, the injection proton density distribution is changed to a broken power-law spectrum, i.e.,    2. In the proton-synchrotron modeling, a higher maximum proton Lorentz factor is required to produce TeV emission than that in the pp interactions.Here, we set α = 1, which implies an extreme acceleration efficiency.3.In the proton-synchrotron model, a large magnetic field is needed to accelerate protons to higher energies and produce higher-energy emissions.We boldly fix the magnetic field B to 35 G for all of five AGNs.4. To maximize the efficiency of the proton-synchrotron process within a reasonable parameter space, we assume that the power of the magnetic field equals to half the Eddington luminosity.Then, the radius of radiation zone can be written as 5. During fitting, we find a significant degeneracy between g p,max and γ p,b .In order to reduce the number of free parameters, we set g g = 10 p,b p,max .
Finally, there are nine free parameters left, which can be found in Table 4.The fitting results are shown in Figure 10.It can be seen that the LHAASO observations are well reproduced for Mrk 421, Mrk 501, 1ES 2344 + 514, NGC 4278 a , and NGC 4278 b .In the case of 1ES 1727 + 502, however, it deviates significantly from the observations, which may be caused by the maximum energy that protons can reach.The characteristic photon energy in the observer's frame produced by the proton-synchrotron process can be calculated by  It is clear that α, B, and Γ are the three parameters that will affect the value of E p,max syn .To increase E p,max syn , α must be lowered, but α = 1 is already the theoretical minimum value.So to increase E p,max syn , alternative acceleration mechanisms with higher efficiency than Fermi first-order are needed.Similar to α, it is also needed to reduce Γ to get a larger E p,max syn .However, reducing Γ will also lead to the observed flux decrease due to the weakening of the beaming effect, unless we increase L p inj at the same time, but that would cause the jet power to exceed the Eddington luminosity.The proton injection luminosity in fitting of four blazars (shown in Table 3) is close to half of the Eddington luminosity, and the power of the magnetic field is assumed previously to be half of the Eddington luminosity, so the sum of the two is very close to the Eddington luminosity.

Figure 1 .
Figure 1.Multiwavelength light curves (LCs) of LHAASO sources.Panels from top to bottom in these six figures: LCs of WISE, ZTF, Swift-UVOT, Swift-XRT, Fermi-LAT, and TS of GeV detection.There is no Swift-UVOT figure for NGC 4278.The meaning of symbols is given in the legend of Mrk 421.Note that the y-axis of some panels does not start from zero.

Figure 2 .
Figure 2. The fitting results of the Swift-XRT spectra.

Figure 3 .
Figure3.One-zone SSC modeling.The meanings of line styles are given in the legend of Mrk 421.The light blue data points are infrared data from WISE, the dark blue data points are optical data from ZTF, the orange and green data points are X-ray data from Swift-XRT's PC mode and WT mode, respectively, and the purple data points are γ-ray data from Fermi-LAT, the red strap shows the observation of WCDA, and the red upper limit point is from KM2A.The gray data points are historical data from the SSDC SED Builder Tool (https://tools.ssdc.asi.it/) of the Italian Space Agency(Stratta et al. 2011).

c
a The SSC+pp model is fitted based on the SSC model for cases of Mrk 421, Mrk 501, and 1ES 1727 + 502, and its leptonic parameters are identical to those of the SSC model, barring an extra proton injection luminosity.In the other three cases, the SSC+pp model has different leptonic parameters than the SSC model.For more detailed discussion, please refer to Section 3.2.b In the SSC+pp model, we assume that the power of the cold protons in the jet is 0.5 L Edd .So the power of the jet should be at least + The corresponding chi-square value per degrees of freedom for each object is calculated by /

Figure 4 .
Figure 4. One-zone SSC+pp modeling.The meanings of symbols and line styles are given in the legend of Mrk 421.

Figure 6 .
Figure 6.Two-zone proton-synchrotron Modeling.The meanings of symbols and line styles are given in the legend.

Figure 7 .
Figure 7.The first fitting strategy of the spine-layer model.The multiwavelength data are explained by the emission from one component.The meanings of symbols and line styles are given in the legend of Mrk 421.

Figure 8 .
Figure 8.The second fitting strategy (two zone) of the spine-layer model, the SEDs are reproduced by the superposition of emissions from two components.The meanings of symbols and line styles are given in the legend of Mrk 421.
, the flare-state VHE band (above 6 TeV) flux of Mrk 501 observed by ARGO-YBJ in 2011 October (Bartoli et al. 2012), the narrow spectral feature at ∼3 TeV of Mrk 501 observed by MAGIC in 2014 July 19 (MAGIC Collaboration et al. 2020a).Some other studies have found that high Doppler factors are necessary for fitting with the one-zone SSC model, e.g., Doppler factors 60 obtained by fitting the H.E.S.S. and Swift data of a TeV flare observed from PKS 2155-304 between 2006 July 28 and 30, Doppler factors 30 for fitting a TeV flare observed in 2001 from Mrk 421 presented the average SED of Mrk 421 in the low state between 2019 January 19 and June 1, and suggested that both the one-zone SSC model and the hadronic proton-synchrotron Blazar model (Mücke & Protheroe 2001) are able to describe the SED well.MAGIC Collaboration et al. (2020b) also found that a flaring state SED of 1ES 2344 + 514, observed in 2016 August, can be successfully described by both the one-zone SSC model and the proton-synchrotron model.Note that a larger magnetic field (B = 50 G for Mrk 421, and B ∼ 50 G for 1ES 2344 + 514) is required in their hadronic scenario.Mastichiadis et al. (2013) found that the observations of Mrk 421 on 2001 March can be naturally reproduced with the leptohadronic model.Sahu & Miranda (2016) concluded that the TeV flaring of Mrk 421 can be well explained by the photohadronic model.MAGIC Collaboration et al. (2020c) suggested that a colocated twozone model is a more reasonable explanation for the overall SEDs of five TeV blazers, including 1ES 2344 + 514 and 1ES 1727 + 502.MAGIC Collaboration et al. (2021) found that the
To reproduce the VHE spectra, the protons with maximum energy should emit at least 20 TeV photons (the energy range of WCDA data is 1-25 TeV).Substituting Equation (5) and

Figure 10 .
Figure 10.One-zone proton-synchrotron modeling.The meanings of symbols and line styles are given in the legend of Mrk 421.
Finally, B is the only parameter that can be adjusted to get a larger E p,max syn .1, Γ, and θ used in fitting into Equation (A4), we can obtain the minimum required magnetic field strength = .The jet is unlikely to have such a strong magnetic field.

Table 3
The Fitting Parameters SSC Model and SSC+pp Model a h Here are the parameters for another fitting strategy of the spine-layer model.The SEDs are reproduced by the superposition of emissions from two components.spectrum of LHAASO is poorly interpreted for Mrk 421 and Mrk 501.20This is caused by the KN effects, which steepens the spectrum naturally.When g , the IC scattering occurs from the Thomson regime into the KN regime, where E 0 is the energy of the soft photon in the comoving frame.Then, we can obtain the critical electron Lorentz factor Then, we derive values of s p of Mrk 421, Mrk 501, 1ES 1727 + 502, and 1ES 2344 + 514, which are 2.66, 2.56, 2.44, and 2.64, respectively.For NGC 4278, since only upper limits are given by Fermi-LAT, we default s p = 2.
under grant No. 12203043.D.R.X.acknowledges the science research grants from the China Manned Space Project with No. CMS-CSST-2021-A06, Yunnan Province Youth Top Talent Project (YNWR-QNBJ-2020-116), and the CAS Light of West China Program.F.K.P. acknowledges support by the National Natural Science Foundation of China (grant No. 12003002), the University Annual Scientific Research Plan of Anhui Province 2023 (2023AH050146), the Excellent Teacher Training Program of Anhui Province (2023), and the Doctoral Starting up Foundation of Anhui Normal University 2020 (903/752022).L. M.S. acknowledges the support from NSFC grant No. 12103002 and Anhui Provincial Natural Science Fondation grant No. 2108085QA43.This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester.Based on observations obtained with the Samuel Oschin Telescope 48 inch and the 60 inch Telescope at the Palomar Observatory as part of the Zwicky Transient Facility project.ZTF is supported by the National Science Foundation under grants Nos.AST-1440341 and AST-2034437 and a collaboration including current partners Caltech, IPAC, the Weizmann Institute for Science, the Oskar Klein Center at Stockholm University, the University of Maryland, Deutsches Elektronen-Synchrotron and Humboldt University, the TANGO Consortium of Taiwan, the University of Wisconsin at Milwaukee, Trinity College Dublin, Lawrence Livermore National Laboratories, IN2P3, University of Warwick, Ruhr University Bochum, Northwestern University and former partners the University of Washington, Los Alamos National Laboratories, and Lawrence Berkeley National Laboratories.Operations are conducted by COO, IPAC, and UW.This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

Table 4
The Fitting Parameters of Proton-synchrotron Model with Strong Magnetic Fields