Bjet_MCMC: A New Tool to Automatically Fit the Broadband Spectral Energy Distributions of Blazars

Multiwavelength observations are now the norm for studying blazars’ various states of activity, classifying them, and determining the possible underlying physical processes driving their emission. Broadband emission models became unavoidable tools for testing emission scenarios and setting the values of physical quantities such as the magnetic field strength, Doppler factor, or shape of the particle distribution of the emission zone(s). We announce here the first public release of a new tool, Bjet_MCMC, that can automatically fit the broadband spectral energy distributions (SEDs) of blazars. The complete code is available on GitHub and allows for testing leptonic synchrotron self-Compton models with or without external inverse-Compton processes from the thermal environment of supermassive black holes (accretion disk and broad-line region). The code is designed to be user-friendly and computationally efficient. It contains a core written in C++ and a fully parallelized SED fitting method. The original multi-SSC zone model of Bjet is also available on GitHub but is not included in the Markov Chain Monte Carlo fitting process at the moment. We present the features, performance, and results of Bjet_MCMC, as well as user advice.

1. One-zone "pure SSC" models, which were mostly relevant for the blazar subclass of high-frequency peaked BL Lacs (HBLs).This class of source was first defined in 1995 as "high-energy cutoff BL Lacs" to replace previous instrumentally based classifications such as X-ray-and radioselected BL Lacs (Padovani & Giommi 1995).We note here that the first SSC theoretical setup can be traced as far back as the late 1960s (Ginzburg & Syrovatskii 1965, 1969).2. Models with thermal external inverse-Compton (EIC) from the interaction of high-energy particles of the jet with the thermal ambient radiation field surrounding the nucleus due to the accretion disk emission reprocessed by the broad-line region (BLR).These SSC+EIC models were primarily used in the modeling of flat-spectrum radio quasars (FSRQs) such as 3C 279 (Sikora et al. 1994;Ghisellini & Madau 1996;Inoue & Takahara 1996).
In all cases, the high-energy emission zone is assumed to be characterized by a compact spherical zone, further referenced as a "blob," relatively close to the nucleus and moving along the jet at relativistic speed.This blob is isotropically filled with highenergy particles (usually simplified as an electron population) and a tangled magnetic field.The blob radiation in its reference frame is also considered isotropic.Most of the SSC models follow the spherical radiation transfer formula set by Gould (1979).The particle distribution within the blob, ρ(E), is characterized by a power law-like spectrum with many possible flavors presenting an average index α ∼ 2-3 (considering ρ(E) ∝ E − α ).Such a distribution is mostly justified by a process of diffuse shock acceleration (e.g., Section 21.4 in Longair 1994).
A known issue of one-zone models is that they usually poorly picture broadband SEDs below the infrared energy range.It is understood that most of the radio emission of jetted AGN is produced by large emission zones, which can be observed in radio very long baseline interferometry as parsec-to-kiloparsec radio core and radio knots.The very extended emission (>100 kpc) is mostly contributing below the centimeter-wavelength energy range.
The C++ code foundation that was later used to build Bjet was developed in the early 2000s for a study of the blazar Mrk 501 (Katarzyński et al. 2001).In order to explain the lowfrequency radiation of Mrk 501, the authors assumed an inhomogeneous model in a conical geometry with a constant bulk Lorentz factor and a power-law decrease of the magnetic 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.
field and particle density along the jet.The SSC of an inhomogeneous jet, discretized in homogeneous slices, was first developed by Marscher (1980) and Ghisellini et al. (1985).The approach of Katarzyński et al. (2001) consisted of two distinct models, one for the blob, "Sblob," and one for the conical jet, "Sjet."Bjet's primary goal was to merge these two models into a consistent multizone framework for the study of the low-frequency peaked BL Lac (LBL) AP Librae.AP Librae is a blazar that displays a multiwavelength SED with features inconsistent with one-zone models, such as a very broad and relatively flat inverse-Compton emission component from X-rays to very high energies (VHE; E > 100 GeV).This issue was tackled with the self-consistent multizone model Bjet that includes radiative interaction between the blob and the conical jet, such as synchrotron self-absorption, radiative absorption by pair creation, and the EIC emission (produced by the interaction of the blob's particles) onto the jet photons (Hervet et al. 2015); see Figure 1, right panel.
In addition to AP Librae, Bjet has been used for modeling several jetted AGN emitting in VHE, such as PKS 0625-354, HESS J1943+213, 1ES 1215+303, PKS 1222+216, and TON 599 (Archer et al. 2018;HESS Collaboration et al. 2018;Valverde et al. 2020;Adams et al. 2022).Since its first use, it has undergone multiple improvements in optimizing the computation time, as well as its scientific completeness, such as 1. better radiation transfer for large angles with the line of sight, 2. EIC from the blob's particles onto the direct disk radiation, and 3. a radial density profile of the BLR for the thermal EIC and its associated gamma-ray absorption by pair creation.The density profile is based on Nalewajko et al. (2014).
A geometrical scheme of Bjet (not to scale) is presented in the left panel of Figure 1.In this paper, we do not review the details of the radiative processes and formulae used in the code.They are described in Katarzyński et al. (2001), Hervet et al. (2015), and, for complementary details, the PhD thesis by Hervet (2015; in French).

Motivations for Bjet_MCMC
As mentioned above, Bjet has been used in multiple papers since its first development in 2015 and has shown its capability in modeling various types of blazars (HBLs, intermediate frequency-peaked BL Lacs (IBLs), LBLs, and FSRQs) and a radio galaxy candidate (PKS 0625-354).The main purpose of the project Bjet_MCMC is to provide this tool to the scientific community through an open GitHub project.5A static version on Zenodo is also available (Hervet et al. 2023).Users can have access to the full Bjet code and perform one-zone pure SSC, one-zone SSC with thermal nucleus interactions (EIC + pair absorption), and multi-SSC zones with thermal nucleus interactions (blob + jet + nucleus).
Another motivation was to make a tool that automatically fits the multiwavelength SEDs and is user-friendly and computationally efficient.SSC models are notorious for being challenging for standard χ 2 minimization methods.They have high dimensionalities, parameter degeneracies, local minima, and model-dependent parameter boundaries.To illustrate this last point, let us consider the particle distribution spectrum within the blob, which is set in our model as a broken power law, where g min , γ brk , and g max are the Lorentz factors of the radiating particles at the minimum, break, and maximum energy of their distribution.In this equation, ( ) , where ( ) N e 1 is the particle density factor set as . Keeping g min , γ brk , and g max free in a standard χ 2 minimization algorithm will certainly create issues, since the condition g g g < < min brk max has to be respected.These constraints let us consider Markov Chain Monte Carlo (MCMC) methods as the most suited to perform an SED fit and explore the parameter space.Indeed, MCMCs have the advantage of building the best solution from posterior probability distributions, which is by default less impacted by discontinuity or nonlinearity of the parameter space.We note here that this MCMC fitting approach of SSC models is relatively known in the community and has been implemented and used in multiple studies (e.g., Tramacere et al. 2011;Zabalza 2015;Qin et al. 2018;Jiménez-Fernández & van Eerten 2021).
In this paper, we do not intend to describe the general statistical concepts behind the MCMC methods.The literature on the subject is quite vast; we can recommend MacKay (2003), for example.

Implementing the MCMC Method
For our project, we used the MCMC emcee package, which is a handy Python tool allowing a relatively simple implementation (Foreman-Mackey et al. 2013).6 emcee requires a user-defined probability function to evaluate the goodness of a fit.It then automatically builds the posterior probability density function following a given number of steps and walkers and a defined burning sample.There are multiple flavors available of what type of move a walker can do on the parameter space; we use the default "StretchMove," developed by Goodman & Weare (2010).A few other moves-or combinations of moves -were tested but did not display significant improvements compared to the proposed default method.

Probability Function and Free Parameters
Our probability function is based on the χ 2 value of the model on all considered SED spectral points.Asymmetric error bars are fully implemented in our χ 2 calculation.We highlight here that flux upper limits are not considered in the fit but can still be included in the input SED data file for display purposes only.We advise to merge constraining upper limits together into larger energy bins until we have a statistically significant data point before fitting the SED.As the emcee package requires a log probability, our probability function is defined as c = -P ln 2 2 .The MCMC method is implemented for the single-zone SSC + thermal EIC model (blob + accretion disk + BLR).It includes up to 13 free parameters, as detailed in the sections below.The full Bjet model currently has 23 free parameters when including the SSC jet.We quickly realized that the computation time required to fit the full multizone model is not reasonable for user-friendly usage. 7In Bjet_MCMC, the user can decide to fix or free any of the 13 parameters.For a pure SSC model fit, the user can deactivate the EIC option to save computation time.

Defining the 1σ Parameter Space and Contour on the SED
The posterior distribution of probability allows us to define the parameter range corresponding to the 1σ confidence level (∼68%) around the best solution.We follow the general solution proposed by Avni (1976) and Lampton et al. (1976) based on the χ 2 cumulative distribution function c cdf 2 or, more precisely, the percent point function c ppf 2 , which returns the χ 2 value associated with a probability P and a number of degrees of freedom k of the c cdf 2 .In this approach, we consider all models that have as within 1σ of the best value, where c min 2 is our best solution, and , where k is the number of free parameters in our model that can range from 1 to 13.
Bjet_MCMC draws the 1σ contour in the SED associated with the parameter uncertainties.Solutions to get the exact contours were ruled out as too demanding of computing power and disk space.For example, one can save all SED points for all models tested during the MCMC process or run through thousands of randomly picked models within the 1σ parameter space.In order to save computation time, this contour is built by picking up models with extremum parameter values at the 1σ confidence level.We call it the "min-max method."This allows us to build relatively good contours by rerunning the code Bjet with only twice the number of free parameters, which typically takes up to a few minutes with all 13 parameters free.Hence, we must warn the user that the contours on the SED plots are an approximation and do not extend to the full theoretical space covered by the parameter uncertainties.

Parameter Boundaries
The pure SSC model has nine free parameters by default, in addition to the following three fixed parameters.
1.The redshift z is set by the user.2. The cosmology is set by default as a flat ΛCDM with H 0 = 69.6 km s −1 Mpc −1 , ΩM = 0.286, and ΩΛ = 0.714 (Bennett et al. 2014).3. The angle between the jet direction and the observer line of sight θ is fixed at 0°.57 to satisfy the Doppler-boosted regime d q < sin 1, with the Doppler factor δ 100.
In the MCMC implementation, we set minimum and maximum values for each of the parameters, as shown in Table 1.We intentionally use wide ranges for the parameters to make as few assumptions as possible.Many parameters are in log scale to ease the walkers' moves through multiple orders of magnitude.
Within our MCMC method, parameter boundaries mean that the posterior probability = P ln 0 for the fit with any parameter outside the given parameter range.The MCMC algorithm acts as an acceptance/rejection method only based on the change of the posterior probability value from one walker step to the other.If a walker moves outside the parameter range, the move will be automatically rejected, and the walker will try again.From the emcee package, a good acceptance rate is about 0.2.Given the additional parameter constraints developed below, Bjet_MCMC shows an acceptance rate on the order of ∼0.05-0.1 from previous tests.
As highlighted in Table 1, multiple parameters have fluctuating boundaries intertwined with other parameter values.The particle spectrum, for example, must follow the condition of g g g < < min brk max , and n 2 > n 1 .The fastest observed variability Dt obs,min is also used to constrain the Doppler factor and size of the emitting region.From the simple argument that in the jet frame, the fastest variation cannot happen faster than the time the light takes to cross the blob radius, we apply the condition


. Finally, a last condition is applied, specifying that the blob diameter cannot be larger than the jet cross section.From a radio study of a large sample of blazars, it can be noted that the intrinsic jet half-opening angle α jet/2 is no more than 5°(e.g., Hervet et al. 2016).Using this conical jet approximation, we set the condition a = jet 2 4. Pure SSC Validation on the HBL 1RXS J101015.9-311909 The blazar 1RXS J101015.9-311900 is an HBL with a redshift of z = 0.143 that has been discovered emitting up to a few TeV by the H.E.S.S. Collaboration after an observing campaign between 2006 and 2010 (H.E.S.S. Collaboration et al. 2012).The multiwavelength SED of this source has been successfully fitted with a one-zone SSC model by Cerruti et al. (2013).In their study, they developed a fitting algorithm that relies on a strong parameterization of the SED features such as slopes and peaks and extended the approach made by Tavecchio et al. (1998).They also probed the parameter space of their SSC model by discretizing it in a grid, where each free parameter is divided into 10 points, with the exception of the particle index n 1 , which is only divided into 3 points.This kind of probing method quickly becomes very computationally heavy for large dimensionality and has only six parameters that can be considered free to mitigate the computation time.On these parameters, they used a different definition of the blob density, which does not allow a straightforward comparison with our results.
Nevertheless, this multiwavelength SED is ideal for a validation test of Bjet_MCM and allows a partial comparison of the results with the work of Cerruti et al. (2013).For this model, we considered 100 walkers, 5000 steps, and a burning phase of 200 steps.We ran it over 15 parallelized threads over a full time of 5 hr 45 minutes.First, we can visually compare the two models in Figure 2. We see notable differences but a relatively equally good fit at first glance.However, we achieved a better c red 2 considering the full data set, or only the X-ray to gamma-ray data set as used by Cerruti et al. (2013).These values are reported in Table 2.We observe a good agreement within the errors on our free parameters.A notable difference is that the best value found for n 1 with Bjet_MCMC is outside the probed parameter range of Cerruti et al. (2013).Overall, this comparison fully validates our approach for single-zone pure SSC models by providing the best SED fit and parameter characterization to date on the blazar 1RXS J101015.9-311909.
Figure 2. Results of Bjet_MCMC on the SED fit of the blazar 1RXS J101015.9-311909.The model applied was a one-zone pure SSC, including a gamma-ray absorption on the extragalactic background light (EBL) following the model of Franceschini & Rodighiero (2017).The SED data points have been shared by Cerruti et al. (2013), while their model line itself has been manually digitized from the original paper.

Fitting the FSRQ PKS 1222+216 with SSC +
Thermal EIC In order to further display the capabilities of Bjet_MCMC, we performed a test on the FSRQ PKS 1222+216 (z = 0.432), which is known to display bright gamma-ray outbursts (e.g., Tavecchio et al. 2011;Adams et al. 2022).It has been noticed that this blazar was a good candidate for SSC + thermal EIC models, as EIC was proposed as the main source of VHE gamma rays.However, without a proper fitting method, it is challenging to fully discard the one-zone pure SSC.This point may actually be the most critical in the relevance of Bjet_MCMC, as it is likely the first public code to date that can provide a fit with a full SSC+EIC model (13 free parameters).As pure SSC and SSC+EIC models are nested, Bjet_MCMC can provide a statistical test that allows one to reject the pure SSC hypothesis.
It has to be noted that MCMC methods can still struggle in complicated parameter spaces and be stuck in local minima.Eventually, a proper MCMC algorithm should always end up in the real best solution.The computation time needed can, however, be unachievable with standard computers.
In this section, we compare the results of Bjet_MCMC on the 2014 flare SED of PKS 1222+216, published by Adams et al. (2022).In their paper, the SED model was manually crafted through a "fit-by-eye" model using Bjet.Being fitted with the same core code (with only minor updates since then), we can have a proper parameter-by-parameter check on how the results of Bjet_MCMC differ from the previous model.
For this fit, we used the same MCMC setup as for 1RXS J101015.9-311900(100 walkers, 5000 steps, 200 steps of burning phase, 15 computing threads).Activating the interaction with the thermal nucleus emission significantly increases the computation time for each step.The full MCMC chain took a total of 12 hr to run.After noticing multiple χ 2 local minima, we fixed two parameters to the Adams et al. (2022) value: the BLR covering factor ò BLR = 1 × 10 −2 and the accretion disk luminosity L disk = 2.8 × 10 46 erg s −1 .As seen in Figure 3, the best fit is visually convincing.We also observe in Appendix Figure 6 that the MCMC walkers display a good general χ 2 convergence.However, we still have hints of local minima, such as in Appendix Figure 6 (top left), where a small fraction of walkers get stuck away from the best fit.They also appear as "islands" in the parameter space corner plot (see Appendix Figure 7).
As shown in Table 3, the results of Bjet_MCMC show a better fit χ 2 compared to the model of Adams et al. (2022).It is interesting to notice that the best solution of Bjet_MCMC does not favor any significant contribution for EIC emission in gamma rays.However, it does not rule out a strong EIC emission either.In the study of Adams et al. (2022), it was estimated that the distance to the supermassive black hole (SMBH) should be at least about 1 pc to avoid too much gamma-gamma absorption from the BLR.This is relatively consistent with the estimated distance of D BH from our fit, which is found above 0.96 pc from the SMBH.

Computational Performance and General Usage Advice
Bjet_MCMC makes full use of the parallelized capability of the emcee package.By running several tests, we observed a roughly linear improvement in the computation time following the number of parallel threads used for the fitting process.We have not performed extensive testing to check if this linear relation was holding true at more than about a dozen threads.It is expected that I/O processes will diminish the relevance of large parallelization at some point.We recommend using a large number of computing threads, if available, likely at least four for the pure SSC and at least 15 for SSC+EIC, if a user wants to get results overnight.Bjet_MCMC will be the most relevant if used in a computing center with several tens of available computing threads.

Quality Checks and Length of MCMC Chains
We propose a few ways to estimate if a user gets enough walkers and steps to be confident in the output of Bjet_MCMC, with some warnings and advice.The favored test to check if the fit is optimal is to get a look at the "median χ 2 per step" plot.For example, one can see in Appendix Figure 4 for 1RXS J101015.9-311909 that the median χ 2 plateaued at about 2000 steps.We can confidently deduce that only 3000 steps would have been enough for this fit, as no further improvements are observed afterward.Now looking at the same plot from PKS 1222+216 in Appendix Figure 6, we observe a plateauing of the median χ 2 at the very end of the chain after an uneven convergence .This should lead to some caution in the estimation of the 1σ parameter space.A weak convergence may not be a big issue, but one should avoid drawing too firm a conclusion from the exact number of the error associated with the parameters.A good practice would be to add an extra 20% on the parameter errors to get a more conservative parameter range when the median χ 2 curve does not fully flatten out.If the median χ 2 curve does not show sign of asymptotic behavior, then the number of steps and/or walkers needs to be increased.
Note that the best χ 2 convergence gives a confidence estimation on the best model, while the median χ 2 convergence gives a confidence estimation of the associated parameter errors.The plot "χ 2 per step" gives a view of the entirety of the walkers.It is a relatively efficient way of checking for local χ 2 minima in the parameter space.As soon as most walkers converge toward the best solution, this should not have significant consequences in the results.If most of the walkers appear to be stuck in local minima, then the SED data set is not constraining enough for the complexity of the model.One can run a longer chain to hope for the MCMC method to eventually converge or reduce the model complexity by freezing parameters.
emcee considers unsafe a number of steps fewer than 50 times the integrated autocorrelation time τ corr (in steps).We find this boundary challenging to achieve in most of our tests, but it also does not systematically lead to better results when met.The value τ corr is given as an output of Bjet_MCMC.We generally advise a safe boundary of n steps 10τ corr .However, a general check on the χ 2 convergence plots and a visual check on the SED itself should always be given to assess the fit quality.

Advice on Building a Multiwavelength SED
Blazars showing great variability should be treated with caution when building their multiwavelength (MWL) SED.
Bjet_MCMC is a stationary model, which means that within the period considered for the SED, we assume an equilibrium between particle injection and cooling mechanisms.The code cannot model any flux variation or describe any time lags between flares observed at different energies.It is possible to model the SED of a flare with Bjet_MCMC, such as PKS 1222 +216 in Figure 3, but we need to ensure that the time period  Note.
a The value ò BLR was considered as fixed in the study of Adams et al. (2022).
used for modeling can be considered as steady.It is done, for example, by using a time period much smaller than the duration of the flare.Longer integration times are often used for lowersensitivity instruments, and nonsimultaneous data sets are often used to build a more detailed SED.There is no single way of building an MWL SED, since most of the time, it consists of finding the best trade-off between data quantity and quality.In any case, data sets, time periods, and assumptions used to build SEDs should be documented by users of Bjet_MCMC.The modeling results are as relevant as the prior assumptions used to build an SED.Some spectral points with very small errors may be overconstraining the fit to the detriment of another energy range.For example, a 10% model flux variation in the optical may be as constraining as a factor of 2 for gamma rays.From a statistical point of view, this is exactly how we expect the fit to behave following the error bars of spectral points.However, the flux "real" errors are often widely underestimated in low energies.The error in each spectral point should reflect the flux variations during the integration period used to build the SED.We advise users to opt for a conservative estimation of the spectral point errors.If the rms errors of flux variation are larger than the statistical errors, the rms errors should be favored when building the SED.It appears that both optical SEDs used in this paper likely have overconstraining error bars.

Other Considerations
A full version of the multizone Bjet code is available through the package Bjet_MCMC.As it contains total of 23 free parameters, it is not included in the MCMC method.The multizone model is recommended to be used only by scientists who have a deep knowledge of blazar emission models, as there are significant risks of having an inconsistent or unphysical set of parameters.This risk has been mitigated for the pure SSC and SSC+EIC models through parameter constraints, but the final assessment of interpreting the quality and relevance of the Bjet_MCMC results is the responsibility of the user.

Conclusion
Bjet_MCMC is a new tool in the growing family of SSC models of blazars.Its full version includes two SSC zones + EIC, based on Hervet et al. (2015).However, only the one-zone pure SSC and SSC+EIC are fully implemented in the automatic MCMC fitting method.Bjet_MCMC aims to be a user-friendly tool that only requires minimal input from the user, namely, a configuration file and an SED data file.It is fully parallelized and can take advantage of computing clusters.The code works as intended and produces consistent results.There are other publicly available SSC models at this time that contain SED fitting algorithms.Among the most known are AGNpy (Nigro et al. 2022), JetSet (Tramacere et al. 2011), andNaima (Zabalza 2015).We do not provide any comparison between Bjet_MCMC and these models in terms of consistency of results, performance, and capabilities.This will be addressed in further studies.However, it has to be noted that Bjet_MCMC appears at the current time to be the only public tool with an automatic fitting that can handle a full SSC+EIC model with up to 13 free parameters.It needs to be noted that excellent broadband energy coverage of the SED has to be built in order to obtain a good convergence of all parameters for the most extensive model version.This tool needs the user to have some knowledge of blazar emission processes to understand the limits of SSC models and infer a scientific interpretation from their outputs.Finally, we highlight that the model is still having frequent updates, and some information mentioned in this paper may be quickly outdated.One of our priorities is to explore routes to further reduce the computation time, such as reducing I/O and further optimizing the energy discretization of the radiative components.The most updated version with relevant information is publicly available on Github. 8

Figure 1 .
Figure 1.Left: scheme of the Bjet model, where dashed lines show the considered radiative transfers.Right: example of an application of the Bjet model on the SED of the blazar AP Librae (Hervet et al. 2015).

Figure 3 .
Figure3.Results of Bjet_MCMC on the SED fit of the blazar PKS 1222+216.The model applied was a one-zone SSC + thermal EIC from the disk and BLR radiative interaction, including a gamma-ray EBL absorption by the model ofFranceschini & Rodighiero (2017).The SED data points have been shared byAdams et al. (2022), while their model line itself has been manually digitized from the original paper.

Table 3
Best Fit and Parameter Comparison for PKS 1222+216