Modeling the Multiwavelength Spectral Energy Distributions of the Fermi-4LAC Bright Flat-spectrum Radio Quasars

In this paper, we present a long-term multiwavelength investigation focusing on 12 distinct samples of Fermi-4LAC bright flat-spectrum radio quasars (FSRQs). Detailed variability and spectral analyses of γ-ray, X-ray, and ultraviolet/optical data obtained by the Fermi Large Area Telescope, the Swift X-ray Telescope, and the Swift Ultraviolet and Optical Telescope were performed over a period of about 14 yr, spanning from 2008 October to 2022 October. These analyses provide insights into characterizing the variations within different activity states. To efficiently reproduce the multiwavelength simultaneous/quasi-simultaneous spectral energy distributions (SEDs) of the samples, we propose a novel approach for constraining the model parameters. By analyzing the parameters of the energy spectral curvature (β), the peak frequency (ν pk), the peak luminosity (L pk), the Compton dominance parameter (A C), and the variability timescale (t var) in different activity states, we can estimate the values of the jet radiation region parameters for the samples. Subsequently, we utilize the synchrotron-self-Compton and external Compton processes, employing a logarithmic-parabolic spectral shape to approximate the observed spectra of the sample sources, while considering the induced regime for the physical parameters. The model results show that: (1) by effectively reproducing SEDs in various active states of bright FSRQs, the parameters within the emission region were reasonably constrained; (2) compared to other active states, the emission region of the jet exhibits a reduced radius during the high state, while the magnetic field strength increases during the low state; and (3) for bright FSRQs in a high-activity state, there is an enhancement of the Doppler factor, often exhibiting a tendency toward energy equipartition.


Introduction
Blazars are an extreme subclass of radio-loud active galactic nuclei (AGNs).Their relativistic jets align with the observer's line of sight while producing Doppler-boosted broadband radiation ranging from radio to very-high-energy (VHE) γ-rays.(Urry & Padovani 1995;Abdo et al. 2010).Flat-spectrum radio quasars (FSRQs) and BL Lacertae (BL Lac) objects are two subclasses of blazars.FSRQs exhibit stronger and broader emission lines (equivalent width >5 Å), especially in the optical-ultraviolet (UV) spectrum attributed to the potential existence of an accretion disk, broadline regions (BLRs; Liu et al. 2006Liu et al. , 2008;;Bai et al. 2009), and a dusty molecular torus (DT; Sbarrato et al. 2012;Zheng et al. 2017;Paliya et al. 2018;Roy et al. 2021).According to the classification by Abdo et al. (2010), the majority of FSRQs are generally categorized as low-synchrotron-peaked (LSP; e.g., 10 peak syn 14 n < Hz) sources, and the inverse Compton (IC) peak often manifests at relatively low energies.The Large Area Telescope (LAT) on board the Fermi Gamma-Ray Space Telescope has been observing the entire sky in the γ-ray band (20 MeV to over 300 GeV) for more than 12 yr (Atwood et al. 2009), providing essential observations for the exploration of high-energy radiation from γ-ray sources.Currently, the fourth catalog of AGNs detected by an incremental version of the fourth Fermi-LAT γ-ray source catalog (4LAC-DR3; 4 Ballet et al. 2020; Abdollahi et al.  2022) has been released, furnishing a substantial number of γ-ray sources alongside their associated observational data.During the initial three months of the Fermi-LAT sky survey operation, 132 bright sources were detected at Galactic latitude b 10 | | > , exhibiting test statistics (TSs) exceeding 100 (∼10σ).These sources were subsequently designated as the LAT Bright AGN Sample (Abdo et al. 2009(Abdo et al. , 2010)).For specific luminous FSRQ populations characterized by relatively high-energy γ-ray emissions, investigating their broadband characteristics to establish reasonable physical constraints remains a valuable avenue for research (Abdo et al. 2009(Abdo et al. , 2010;;Sikora et al. 2009;Zhang et al. 2012;Kang et al. 2014Kang et al. , 2020;;Paliya et al. 2018;Anjum et al. 2020;Roy et al. 2021;Zhou et al. 2021;Zhu et al. 2021).
Numerous previous studies have focused on the multiwavelength variability and spectral analysis of blazars in recent years (Sikora et al. 2009;Böttcher et al. 2013;Dermer et al. 2014;Liao et al. 2014;Xiong et al. 2017Xiong et al. , 2020;;Paliya et al. 2018;Kang et al. 2021;Sahakyan 2021;Sahakyan et al. 2022).With the continuous advancement of computer algorithms, rigorous optimization methods-e.g., maximum likelihood estimation (Foreman-Mackey et al. 2013;Brien 2019), the least-squares method (Hogg et al. 2010;Mankuzhiyil et al. 2011;Fan et al. 2016;Kang et al. 2016), and/or Markov Chain Monte Carlo (MCMC) simulation (Foreman-Mackey et al. 2013;Yan et al. 2013;Qin et al. 2018;Kang et al. 2021)-have emerged for modeling the SEDs to achieve optimal results within the leptonic scenario.Modeling and analyzing the SEDs at different time periods enables the connection between observational properties and the underlying physical processes occurring within the jet.However, when dealing with multiple free parameters, the uncertainty in the range of each parameter leads to extensive computational time and reduced fitting efficiency.Furthermore, there is a possibility of encountering parameter combinations that are inconsistent with the observed conditions (Ghisellini et al. 2010;Tanaka et al. 2016).To address this issue, we propose the utilization of an analytical expression incorporating extended EC components, specifically the parameters of the energy spectral curvature (β), peak frequency (ν pk ), peak luminosity (L pk ), the Compton dominance parameter (A C ), and variability timescale (t var ), based on the logarithmic-parabolic (LP) model-fitting algorithm developed by Zhou et al. (2021).This approach aims to estimate the parameters of the jet radiation region during different active states of 12 Fermi-4LAC bright FSRQs and effectively reproduce their multiwavelength quasi-simultaneous SEDs.
The paper is structured as follows: in Section 2, we present the sample description.In Section 3, the data reduction, processing, and analysis are presented, including the presentation of the observed multiwavelength light curve and SEDs.In Section 4, the parameter constraints, model, and the results are described.A summary of our results is provided in Section 5. Throughout the paper, we assume the Hubble constant to be H 0 = 67.8km s −1 Mpc −1 , the matter energy density to be Ω M = 0.31, the radiation energy density to be Ω r = 0, and the dimensionless cosmological constant to be Ω Λ = 0.69 (Planck Collaboration et al. 2016).

Sample Description
The Fermi-LAT collaboration presented the fourth Fermi-LAT γ-ray Source Catalog (4FGL-DR1) in 2020, which covers 8 yr of observations from the Fermi-LAT Sky Survey, spanning 2008 to 2016 (Abdollahi et al. 2020).The fourth catalog of AGNs detected by Fermi-LAT (Ajello et al. 2020), based on 4FGL-DR1, comprises a total of 3207 AGNs.Fermi-LAT conducted its third data release recently, which is based on 12 yr of E > 50 MeV γ-ray data, including 4FGL-DR3 (Abdollahi et al. 2022) and 4LAC-DR3 (comprising 3814 AGNs; Ajello et al. 2022).The blazars in 4LAC-DR3 include 792 FSRQs, 1458 BL Lac objects, and 1493 blazar candidates of uncertain type.We filtered the FSRQs in 4LAC-DR3 based on the source significance (Signif_Avg in 4LAC-DR3 FITS Format) and obtained 39 bright FSRQ sources with TS > 100.To ensure the adequacy and completeness of the multiwavelength data for bright sources, we compiled the observation IDs (including UV, optical, and X-ray) from the High Energy Astrophysics Science Archive Research Center (HEASARC) website. 5ltimately, we selected 12 bright FSRQs with observation ID counts exceeding 50 as our samples, and the relevant parameters for these samples are displayed in Table 1.We obtained the radio and infrared flux data for the samples from both the ASI Space Science Data Center (SSDC),6 the NASA/ IPAC Extragalactic Database (NED),7 and literatures (Glauch et al. 2022).Additionally, SSDC provided historical flux data of various periods for the samples across multiple wavelength bands.Since we require broadband data for samples in various activity states, we utilize the threads of the Neil Gehrels Swift Observatory and Fermi-LAT for data processing to derive longperiod multiwavelength light curves and quasi-simultaneous SED data during different activity states.In the following section, we will provide a detailed description of the data processing procedures applied to the sample.

Fermi-LAT Observations
The γ-ray data used in this paper were collected by Fermi-LAT from 2008 October to 2022 October (MJD 54740 to MJD 59854).We have downloaded and analyzed the newest Pass 8 SOURCE data in the energy range from 100 MeV to 1 TeV using the standard reduction methodology recommended by the Fermi Science Support Center. 8The events classified as evclass = 128 and evtype = 3 within a circular region of interest (ROI) of 10°centered on the γ-ray positions of the targets are subject to analysis using the Fermi Science Tools (version 2.0.8) with the P8R3_SOURCE_V3 instrument response function.As suggested by the Fermi-LAT team, we apply the quality cuts, (DATA_QUAL == 1)&& (LAT_CONFIG == 1), and also implement a zenith angle cut at 90°to eliminate events originating from the Earth's limb.The good time intervals were selected using the gtmktime task.Then a count map of the ROI was generated using the gtbin task, while the exposure map was created with the help of both the gtltcube and gtexpmap tasks.The Galactic diffuse emission component, as well as the isotropic background templates,9 were utilized with P8R3_SOURCE_V3 and iso_P8R3_SOUR-CE_V3_v1.txt, respectively.The XML files were generated using the user-contributed tool make4FGLxml.py.The diffuse source responses were generated using the gtdiffrsp tool, and the flux was obtained by executing the gtlike task.The gtlike tool was employed to estimate the γ-ray flux and the photon index of the sources.The significance of the γ-ray detection was quantified using the maximum likelihood TS = L 2 log( ) D , where L represents the ratio of likelihood values for models with and without a γ-ray point source (Mattox et al. 1996;Atwood et al. 2009).The γ-ray light curves were calculated using the standard unbinned likelihood data reduction procedure as described in the documentation,10 along with the corresponding quality cuts provided therein.The γ-ray photon indices and normalized parameters are allowed to vary freely during the model fitting for sources located within a 10°r adius from the center of the ROI.For sources located outside a 10°radius from the center of the ROI, their parameters are held constant at the values provided in 4FGL-DR3.The normalization of the diffuse background components remains free throughout the model-fitting process.Furthermore, we computed the γ-ray photon spectra of Fermi-LAT sources using Fermitools (v11r5p3) and the open-source Python package Fermipy (Wood et al. 2017).We utilized the appropriate instrument response function P8R3_SOURCE_V3, the galactic interstellar emission model gll_iem_v07 (i.e., "gll_iem_v07.fits"), and a newly introduced isotropic spectral template "iso_P8R3_SOURCE_V3_v1. txt" (Abdollahi et al. 2020).During the data processing, we also utilized the easyFermi community tool, 11 an open-source graphical interface built upon Fermitools and Fermipy, to efficiently generate γ-ray photon spectra and light curves for the samples (de Menezes 2022).For easyFermi data processing of energy spectra and the light curve, we selected LAT 0.1-300 GeV Pass 8 events recommended by the Fermi-LAT collaboration from a circular ROI with an angular radius of less than 15°.During the XML model-fitting process, which encompassed power-law, log-parabola, and power-law with an exponential cutoff spectral models, for sources located within a 5°radius from the center of the ROI, we allowed for variations in the γ-ray photon index and normalization parameters.When generating γ-ray photon spectra using easyFermi tools, we utilized the default settings: a minimum significance of 4.0 and a minimum separation of 0°. 5 within the ROI for finding extra sources.Additionally, we divided the energy range from 0.1 to 300 GeV into 10 logarithmically equal energy bins.Data points with TS values less than 10 were considered as 95% confidence upper limits.

Swift Observations
The Neil Gehrels Swift Observatory, a space-based observatory, is equipped with three main instruments: the Ultraviolet and Optical Telescope (UVOT; Roming et al. 2005), the X-ray Telescope (XRT; Burrows et al. 2005), which operates in the 0.3-10 keV band, and the Burst Alert Telescope (Barthelmy et al. 2005), sensitive to the 15-150 keV band.Swift is an ideal satellite for simultaneous/quasi-simultaneous observations of blazars in the X-ray, optical, and UV bands.The XRT data were acquired in both photon counting (PC) mode and windowed timing (WT) mode, with individual exposures ranging from 0.22 to 14.35 ks, resulting in a total exposure time of approximately 0.76 Ms.We retrieved data from 2008 October to 2022 October from the HEASARC website. 12For data analysis at level I, we followed the standard threads outlined in the Swift-XRT analysis documentation. 13We processed all Swift data using HEASoft 6.26.1.The xrtpipeline task was executed, selecting a circular source region with a radius of approximately 19 pixels (45″) and an annular background region with an inner radius of approximately 43 pixels (100″) and an outer radius of approximately 85 pixels (200″).The light curves and spectra were generated using xselect with level II data from the the WT and PC modes, respectively.To enhance the signal-to-noise ratio (S/N) of the merged X-ray spectrum, the spectra from the samples were binned using grppha, with a minimum of 20 photons per bin for WT-mode spectra and a minimum of 10 photons per bin for PC-mode spectra.In the reduction of the X-ray spectra, the response matrix files swxwt0to2s6_20131212v015.rmf for the WT mode and swxpc0to12s6_20130101v014.rmf for the PC mode were employed, and the standard ancillary response files were generated using xrtmkarf.The individual X-ray spectra were fitted using XSPEC 12.12.0, 14an X-ray spectral fitting package, which employed an absorbed and redshifted powerlaw (zPL) model 15 (e.g., ) with neutral hydrogen column density (N H ), while ignoring channels with energies below 0.3 keV and above 10 keV.Galactic absorption (i.e., within the Milky Way) is included in all of these spectral models (i.e., tbabs×zPL, where tbabs represents the absorption due to the interstellar medium).

Archival Data
We have been collecting multifrequency flux observations from each sample since 2008 October, using data from the SSDC.The SSDC combines flux data across the radio and γray (including TeV) bands, from various missions, external services (e.g., NED, the Sloan Digital Sky Survey, and the United States Naval Observatory), experiments, as well as catalogs and archival data.In multiwavelength SED modeling, we compare historical flux data with quasi-simultaneous data to refine the fitting.

Multiwavelength Light Curves and SEDs
Based on the data processing of Swift-UVOT, Swift-XRT, and Fermi-LAT mentioned above, 12 sample sources can be processed uniformly within a specific period to obtain their quasi-simultaneous observations in the optical-UV, X-ray, and γ-ray bands.Given that the long-term multiwavelength light curve of the source is instrumental in enabling us to assess and analyze flares (or outbursts, also referred to as flux surges) as well as the quiescent state of the source during various time intervals, we processed the 14 yr γ-ray light curve of each source within the sample, corresponding to their respective Swift observation IDs.We categorized each source into high states (i.e., outbursts), low states (i.e., quiescence), and average states using the following methods.
delineating the high-state (H) period.If there is no corresponding outburst in the γ-ray light curve during the X-ray outburst interval, we consider the X-ray outburst to be associated with a period of the γ-ray outburst (with the fastest rising flux).This ensures the concomitant occurrence of outbursts in both X-rays and γ-rays within the same H period. 2. We initially identify the period characterized by the lowest flux in the X-ray light curve as the low-state (L) interval.If this period coincides with an outburst in the γray light curve, we designate the period with the lowest flux in the γ-ray light curve to represent the L for X-rays.This ensures that both X-ray and γ-ray fluxes remain relatively stable and low during the L period.During the L period, due to the diminished flux of γ-ray photons, the number of photons collected within a few days is not optimal.Consequently, we consider a longer period, spanning from one month before to one month after, for generating γ-ray spectra in such scenarios.3. We use Fermitools and Fermipy, respectively, to calculate the average γ-ray spectrum from 14 yr of Fermi data for each source in the energy range of 100 MeV to 1 TeV, with each spectrum divided into 10 energy bins.Moreover, xselect and ximage, components of the HEASOFT software, are employed to combine the observational events and images from all Swift observation IDs for each source, followed by the use of xrtmkarf to generate the response file.Finally, we utilize grppha to derive the average X-ray spectrum.The Cash statistical (Cash 1979) method was employed for the analysis of ungrouped data, as the count numbers were relatively low for some observations.Nevertheless, for observations with a high count rate, a cross-check was conducted by rebinning the data to have at least 10 counts per bin, followed by fitting using the χ 2 minimization technique.Additionally, in accordance with the UVOT data analysis threads, 18 we employ the uvotimsum tool to merge different images collected for the same filter and the same Swift observation ID, before running uvotsource.Subsequently, the UVOT images summed for each of the six filters for each source are processed and subjected to reddening correction to yield averaged UV-optical data.
We combine these averaged γ-ray, X-ray, UV, and optical data sets as our average-state (A) period.
In the process of handling multiwavelength quasi-simultaneous data, due to the absence of radio and infrared data during various activity phases, we incorporate historical data provided by the SSDC to constrain the low-energy portion of the SED.For each source in the SSDC data set, we process it by taking the arithmetic mean of multiple fluxes at the same frequency to obtain low-energy A data; we use the maximum flux among multiple fluxes at the same frequency for the low-energy H data; and we select the minimum flux among multiple fluxes at the same frequency for the low-energy L data.In addition, for the H and L data, we employed easyFermi with the default energy range of 100 MeV to 300 GeV for photon events to generate their γ-ray photon spectra.We conducted multiwavelength quasi-simultaneous data for each source within the sample across various activity states.The corresponding results of this data processing are presented in Tables 1-2 and Figure 1 (See Figures A1-A12 in the Appendix for the results of all samples).
When processing the 14 yr γ-ray light curves from Fermi-LAT, we adopted a 30 days binning interval for γ-ray flux.Furthermore, during periods of high activity, we can constrain the parameters of the radiation region by analyzing the outburst profiles to determine the variability timescale of the sources.Consequently, we analyzed the light curves of sample sources within a 6 hr interval, and the processing results are shown in the top panel of Figure 1(a) and Table 1.As indicated in Table 1, we used the easyFermi tool to generate γ-ray light curves with 6 hr time bins before and after extracting the H states for each source.To quantitatively determine the rise and decay timescales of flux variations, we employ the following double exponential functional form to fit the time profiles of γray flares during high-activity states (Zheng et al. 2013(Zheng et al. , 2014;;Gasparyan et al. 2018): where F c is the quiescent flux, t pk denotes the time of the flare peak (F 0 ), and t r and t d are the rise and decay times, respectively.Each light curve within the sample was fitted using the nonlinear optimization Python package lmfit. 19The results of the fitting are shown in Table 3.

Results
Figure 1(a) (see also Appendix Figures A1-A12) displays the multiwavelength light curves of the sample spanning from 2008 to 2022, covering the γ-ray, X-ray, UV, and optical bands.In the first panel, the data were extracted from Fermi-LAT (30 days binning).An enlarged view in the top panel denotes the γ-ray light curve at a 6 hr binning at the H state.The second panel displays the γ-ray photon index.The third panel displays data obtained from Swift-XRT, while the fourth panel presents the X-ray photon index in the 0.3-10.0keV range.The fifth and sixth panels exhibit data from Swift-UVOT, including flux measurements in the V, B, and U filters, as well as the UVW1, UVM2, and UVW2 filters.The purple and gray shaded areas represent H and L activity states, respectively.
Figure 1(b) (see also Appendix Figures A1-A12) displays the fitting results of the LP model algorithm under different activity states.Building upon Zhou et al. (2021), we introduced MCMC techniques (based on the emcee20 software package, as described in Foreman-Mackey et al. 2013) to fit multiwavelength observations across various activity states, obtaining phenomenology-related spectral parameters: the energy spectral curvature (β), the peak frequency (ν pk ), the peak luminosity (L pk ), and the Compton dominance parameter (A C ).
To prevent anomalous fitting of low-energy peaks resulting from the lack of simultaneous radio data, we employed flux values in the 1.4-5 GHz range from SSDC to constrain the fitting when using the LP model.The corresponding results are shown in Table 4.  Here, "C-stat.,""Γ x ," and "K p " denote the values of the Cash statistics, the photon index of the X-ray, and the normalization factor in the zPL fits, respectively.The fitting parameters and goodness evaluations were obtained from xspec 12.10.1,(https://heasarc.gsfc.nasa.gov/xanadu/xspec/)an X-ray spectrum fitting program."dof" are the degrees of freedom." F log x " denotes the flux in 0.3-10 keV in units of erg cm −2 s −1 ." F log g ," "K γ ," "P1," "P2," and "TS" represent, respectively, the γ-ray flux in the 0.1-300 GeV energy range in units of erg cm −2 s −1 , the normalization factor of the XML spectral model, the characteristic parameter 1 of the XML spectral model (α in logParabola(LP) or index1 in PLSuperExpCutoff2(PLEC)), the characteristic parameter 2 of the XML spectral model (β in LP or index2 in PLEC), and the TS value from the likelihood analysis (https://fermi.gsfc.nasa.gov/ssc/data/p6v11/analysis/scitools/source_models.html)."Type" denotes the spectral models of the Fermi-LAT sources available for use in gtlike.

Constraining the Parameters of the Jet Radiation Region
We aim to establish a self-consistent connection between observational parameters related to phenomenology, derived from data analysis, and radiation region parameters associated with physical models.Similar to the derivation approach in Zhou et al. (2021), within the framework of the lepton model, we introduce external radiation field components and derive from the perspective of the electron energy distribution (EED) in the radiation region.We consider a homogeneous spherical radiating region with a radius of R b ¢ , containing ultrarelativistic electrons, where a uniform and randomly oriented magnetic field strength B¢ exists (Ghisellini & Maraschi 1989;Ghisellini et al. 2010;Zhang et al. 2012;Kang et al. 2014Kang et al. , 2016Kang et al. , 2021;;Chen 2017;Zheng et al. 2019).In the context of stochastic acceleration with mono-energetic and instantaneous injection, the EED is characterized by the LP distribution, expressed in the form (Massaro et al. 2006;Paggi et al. 2009Paggi et al. , 2011;;Zheng et al. 2018b;Tan et al. 2020;Zhou et al. 2021) where s denotes the electron spectrum index, r denotes the electron spectrum curvature, N 0 denotes the normalization constant in units of cm −3 , 0 g ¢ represents the initial electron where t var denotes the observed variability timescale; β syn represents the energy spectral curvature of the Syn bump; L pk syn and L pk SSC denote the peak luminosities of Syn and the SSC energy spectrum, respectively, measured in erg s −1 ; A C is the dominant parameter of Compton in the EC radiation field; and u ext denotes the energy density of the external radiation field, expressed in erg cm −3 .u ext can be calculated as

=
´(see, e.g., Paliya 2015; Gasparyan et al. 2018).During the L and A states, we simplify the given variability timescales using t var = 1 day in the source coordinate system (Ulrich et al. 1997;Nalewajko 2013;Kang et al. 2014).Additionally, based on the γγ attenuation argument, we can numerically estimate the lower limit of the Doppler factor during the H state as follows (Dondi & Ghisellini 1995;Finke et al. 2008;Paliya 2015;Ding et al. 2019;Geng et al. 2022): where d L represents the luminosity distance to the source with a redshift z and σ T is the Thomson cross section.During the flare, the energy of the highest-energy photon is E ,max g , and the contemporaneous X-ray flux is f x  .The estimation results of the jet radiation region parameters under various activity states are presented in Table 3  , respectively.These values closely match the estimated range for FSRQs provided by Ghisellini et al. (1998Ghisellini et al. ( , 2010) ) and Ghisellini & Tavecchio (2015) As long as we provide the EED within the blob in the comoving coordinate system of the jet, we can calculate the various radiation components of the jet, encompassing Syn, SSC, and EC (Yan et al. 2012;Kang et al. 2014;Zheng et al. 2019).
Assuming an isotropic distribution of electrons described by Equation (2), the theoretical Syn spectrum can be calculated as follows:  where e represents the electron charge, B¢ is the magnetic field strength is the energy density of the Syn photons in units of erg cm −3 (Yan et al. 2012;Zheng et al. 2019Zheng et al. , 2020)) represents the energy density of the external photon field, measured in units of erg cm −3 .In the EC case, we assume that the seed photons predominantly originate from the BLR and the DT.For the dissipation region within the BLR, the seed photon energy density is given by u , where τ BLR ∼ 0.1 represents the fraction of the disk luminosity, L disk , that is re-emitted by the broad lines.Here, f p is the Planck function for blackbody radiation and r BLR is the typical size of the BLR (Inoue & Takahara 1996;Dermer & Schlickeiser 2002;Ghisellini & Tavecchio 2008, 2009;Sikora et al. 2009;Potter & Cotter 2013;Kang et al. 2014).In the case where the dissipation region is outside the BLR, the seed photons from the BLR decrease quickly and the main contribution comes from the DT.Similar to u BLR , we have , where τ DT ∼ 0.5 and r DT represents the typical size of the DT.The external radiation field is characterized by an isotropic blackbody with a temperature of T BLR/DT = hν p /(3.93k B ), where ν p represents the maximum frequency of the seed photons in the ν − νF ν space.In the case of the BLR, this is expressed as 2 10 Hz  et al. 2007;Ghisellini & Tavecchio 2008;Kang et al. 2014Kang et al. , 2016Kang et al. , 2021)).
Hence, the observed spectrum for the total emissions can be expressed as where τ γγ (ν, z) is the absorption optical depth due to the interaction of VHE γ-ray photons with the photons from the extragalactic background light (Franceschini et al. 2008;Franceschini & Rodighiero 2017;Zheng et al. 2019).
We compared the Swift-UVOT data with archival data and found a thermal radiation bump in the UV to optical bands.For this portion of the data, we consider adopting the approach proposed by Zheng & Yang (2016), which involves calculating the radiation spectrum originating from the accretion disk for interpretation.A simple Shakura-Sunyaev disk spectrum can be calculated as follows: where ò is the photon energy that is emitted by the accretion disk.While the value of max  depends on the spin of the black hole and relative Eddington luminosity, we set a typical characteristic temperature of the UV bump in Seyfert galaxies with 10 eV max  . The values for the accretion disk luminosity (L d ) that are listed in Table 4 were obtained from the following literature: Ghisellini et al. 2010, Zhang et al. 2015, and Paliya et al. 2021.
The parameters of the jet radiation region, as estimated in Section 4.1, serve as the initial values for modeling the theoretical photon spectrum.This enables us to reproduce the multiwavelength SED under different active states, utilizing Equations (8) to (12).There are 12 parameters in our SED modeling: R b ¢ , B¢, δ, pk g ¢ , min g ¢ , max g ¢ , s, r, N 0 ¢ , r BLR , r DT , and L d .
To reduce the degrees of freedom for the parameters and effectively achieve optimized results, we imposed the  = ´ (Kaspi et al. 2007;Ghisellini & Tavecchio 2008;Bentz et al. 2009 g ¢ = ´(3) data points with energies exceeding 10 GeV were excluded from our analysis due to their VHE nature; and (4) we optimized the fitting results by adjusting the initial parameter values of the jet radiation region, guided by the minimum χ 2 process (Kang et al. 2014;Kang 2017).Therefore, after imposing the above restrictions, the number of free parameters in the actual fit is reduced to seven (R b ¢ , B¢, δ, pk g ¢ , s, r, and N 0 ¢ ), and all of them vary within a local range.The SED modeling results corresponding to these parameters are presented in Table 5 and Figure 1(c).Figure 1(c) (see also Appendix Figures A1-A12) displays the multiwavelength SEDs reproduced by our one-zone leptonic model, which includes Syn radiation, SSC, and EC components, under the H, L, and A period.The distribution of parameters in the jet radiation region under various active states is presented in Figure 3.The analysis of Figure 3 yields the following results: (1) during the H state, the radius of the jet radiation region in the samples is smaller compared to other active states, and during the L state, the magnetic field strength is greater compared to other active states; (2) based on the 1σ confidence ellipses, we provide reasonable ranges for the parameters of the emission region in different active-state samples as follows: 0.88 log , where N is the number of observed data points and dof is the degrees of freedom, i.e., the number of free parameters used for the model; y i denotes the expected values from the model and y i denotes the observed data.It is known that both the particle and field energy significantly contribute to the overall energy content of jet plasma.The interplay between these factors is intricately linked to the mechanisms governing jet formation, particle acceleration, and radiation (Dermer et al. 2014;Zheng et al. 2017Zheng et al. , 2018a)).As a means to elucidate this relationship, we introduce the equipartition fraction, denoted as u u  The parameter distributions and ranges of the jet radiation region are obtained from our modeling of the SEDs in various activity states.The orange dots and orange curves represent the parameters during the "L"; the blue dots and blue curves represent the parameters during the "H"; and the green dots and the gray curves represent the parameters during the "A."The dotted line represents the confidence ellipse of 1σ (Kang et al. 2021).The red cross marks the intersection of the 1σ candidate region obtained from Gaussian fits for the parameter (Foreman-Mackey 2016).

Summary
We conducted an investigation using quasi-simultaneous multiwavelength observations for 12 Fermi-4LAC bright FSRQs obtained through long-term light-curve and spectral analysis.By establishing a correlation between the observed parameters and the physical model parameters of the jet emission region, we effectively constrained the parameter ranges applicable to these bright FSRQs across different activity states.Subsequently, we effectively reproduced their multiwavelength SEDs in various activity states.Additionally, we quantified the energy equipartition within distinct activity states.The variations in the parameters of the jet emission region in different activity states may indicate changes in the conditions of the emission region.In the H state, we found that the jet's emission region has a smaller radius compared to other activity states, suggesting that it may be closer to the central engine during high activity.In the L state, the slightly larger magnetic field could reflect potential magnetic energy conversion mechanisms.
Due to limitations in the number of observed samples and the complexity of the data processing, our conclusions are more applicable to bright FSRQs.The task of estimating radiation region parameters in various activity states with a more extensive sample is left for future work.

Figure 1 .
Figure 1.Example (J1224.9+2122) of the multiwavelength light curves and SEDs.(a) Multiwavelength light curves covering γ-ray, X-ray, UV, and optical bands from 2008 to 2022 with corresponding epochs.(b) The outcomes of fitting multiwavelength SEDs using the LP model in various activity states.(c) Modeling multiwavelength SEDs in various active states using the leptonic model of SSC+EC with an LP spectral shape.
g¢ is the electron Lorentz factor.Throughout the article, quantities represented with a prime are calculated in the comoving coordinate system.Based onZhou et al. (2021) andChen & Bai (2011), in the Thomson scenario, we have obtained the following analytical of the emission blob.Additionally, h denotes the Planck constant.In this context, photon energy in the comoving frame, where m e is the rest mass of the electron and c is the speed of light.The coefficient for synchrotron self-absorption is expressed as(Rybicki & Lightman 1979) the angle between the magnetic field and the velocity of high-energy electrons, and K 5/3 (t) is a modified Bessel function of the second kind of order 5/3.The theoretical SSC/EC spectrum is given by

Figure 2 .
Figure 2. Distribution of estimated values for the parameters of the jet radiation region log est d , R log b,est ¢ , B log est ¢ , and log pk,est g¢ .
states of activity, we have computed the equipartition fraction and determined their ranges through 1σ Gaussian fits as follows:1.150.80 4.91 -+ for the H; 0.11 0.09 1.49 -+ for the L; and 0.42 0.31 6.79 -+for the A. Thus, at the H, energy equipartition is more closely achieved due to the larger value of A C in the source.

Figure 3 .
Figure3.The parameter distributions and ranges of the jet radiation region are obtained from our modeling of the SEDs in various activity states.The orange dots and orange curves represent the parameters during the "L"; the blue dots and blue curves represent the parameters during the "H"; and the green dots and the gray curves represent the parameters during the "A."The dotted line represents the confidence ellipse of 1σ(Kang et al. 2021).The red cross marks the intersection of the 1σ candidate region obtained from Gaussian fits for the parameter (Foreman-Mackey 2016).

Figure A1 .
Figure A1.The multiwavelength light curves and SEDs for J0210.7-5101(PKS 0208-512).(a) Multiwavelength light curves covering γ-ray, X-ray, UV, and optical bands from 2008 to 2022 with corresponding epochs.(b) The outcomes of fitting multiwavelength SEDs using the LP model in various activity states.(c) Modeling multiwavelength SEDs in various active states using the leptonic model of SSC+EC with an LP spectral shape.

Figure A3 .
Figure A3.The multiwavelength light curves and SEDs for J1159.5+2914(Ton 599).The explanations for each subfigure are the same as in Figure A1.

Figure A11 .
Figure A11.The multiwavelength light curves and SEDs for J2232.6+1143(CTA 102).The explanations for each subfigure are the same as in Figure A1.

Table 1
Summary of Observations for Swift-XRT, Swift-UVOT, and Fermi-LAT in Our Sample

Table 2
Summary of Fitted Results of the X-Ray Spectra by zPL (tbabs×zpowerlw) and the γ-Ray Spectra in Our Sample and Figure 2. The histogram distributions of the parameters log est 11 g