© ESO 2022

1. Introduction

Active galactic nuclei (AGN) with jets, which are the brightest among the persistent extragalactic sources, present broadband emission covering the whole electromagnetic spectrum. Modelling this emission with the radiative processes of relativistic charged particles allows one to derive their non-thermal energy distributions, to infer their acceleration mechanism and site, and to investigate the jet composition (see Cerruti 2020 for a review centred on the radiative processes and Böttcher 2019 for a broader overview of the topic). The operation in the last two decades of space-borne and ground-based gamma-ray instruments has extended the spectral energy distributions (SEDs) of these sources with simultaneous data up to TeV energies (Madejski & Sikora 2016). Physical interpretation of their multi-wavelength (MWL) observations has traditionally been provided, even in the publications by the experiments and observatories gathering the data, by a few individuals with proprietary and private modelling software. The latter, though advancing the theoretical understanding of the field, is not accessible to the community at large and produces results difficult to reproduce or verify (Portegies Zwart 2018). Moreover, the operation of the next generation of gamma-ray observatories such as the Cherenkov Telescope Array (CTA) and the larger sample of simultaneous complete SEDs they will facilitate (CTA Consortium 2019, Chs. 8 and 12), calls for the modelling and interpretation effort to be opened to a wider number of astrophysicists. For all these reasons, it is desirable to develop a new generation of open-source AGN modelling software.

It has become increasingly common to build open-source astrophysical libraries in the ecosystem formed by the NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), and astropy (Astropy Collaboration 2013, 2018) libraries. The astropy project gathers ∼40 packages built on these fundamental dependencies in a community of ‘affiliated packages’ committed to create reusable and interfaceable astronomical python software¹. Notable among these are Gammapy (Deil et al. 2017), for the analysis of high-energy astrophysical data, adopted as science analysis tool for CTA, and naima (Zabalza 2015), the first python package modelling the non-thermal SED of astrophysical objects. Designed for sources with unbeamed radiation, naima computes the radiative processes assuming comoving densities of interacting particles. It is therefore appropriate for interpreting the emission of galactic sources like supernova remnants, but inadequate when considering jetted extragalactic sources. In this case, relativistic transformations have to be applied on the interacting particle (or photon) densities, that are usually moving at relativistic speed with respect to each other. Additionally, the emission from the jet has to be boosted, accounting for its beaming in a cone inclined at an angle θ_s with respect to the observer line of sight. jetset (Tramacere 2020) offers an open-source modelling framework for extragalactic and galactic jetted sources: blazars (Massaro et al. 2006; Tramacere et al. 2009, 2011) and microquasars (Rodi et al. 2021). The jetset core libraries compute the photon spectra due to leptonic (due to electrons and positrons) radiative processes and pp interactions and evolve the particle distributions accounting for their acceleration and cooling. Written in C, they are wrapped with a high-level python interface relying on NumPy and astropy, connecting the framework to this environment. The high-level interface includes routines for MWL SED data handling and fitting using other python data-analysis packages as backend.

We created agnpy with the objective of providing an open-source python package for jetted AGN modelling entirely built in the NumPy, SciPy, and astropy ecosystem, interfaceable with data-analysis tools such as Gammapy and completing the science cases covered by naima and jetset. agnpy gathers and implements the leptonic radiative processes most commonly considered in jetted AGN, its main objective being modelling the broadband SED, from radio to gamma rays, produced by these sources. Expanding the package’s scope beyond the non-thermal SED computation, we included additional classes describing the AGN components emitting thermal and line radiation. agnpy is able to compute the absorption due to γγ pair production on the soft photon fields of the different AGN components and provides the absorption due to the extragalactic background light (EBL). Though not designed for time evolution, the package offers the possibility of a self-consistent constraint of the model parameters. Rather than solving the differential equation regulating the evolution of the electron energy distribution (EED), agnpy assigns its maximum and break Lorentz factors using a simple parametrisation of the acceleration, escape, and radiation processes. Throughout all the physical processes described in the package, the angle of the jet axis to the observer line of sight, θ_s, is always considered among the model parameters, such that the use of agnpy is not restricted to the blazar case (θ_s → 0) but can also be extended to radio galaxies.

This article is structured as follows: in Sect. 2 we describe the package implementation and content; in Sect. 3 we expose some of its basic functionalities and in Sect. 4 we validate its results against the literature and other modelling software. In Sect. 5 we discuss reproducibility and performance of our software. We gather our conclusion and thoughts for future expansions in Sect. 6. The appendixes collect the mathematical derivations used in the code implementation, the benchmark model parameters used for the validation and some technical considerations.

2. Package description

The basic formulae and notation used in the package follow Dermer & Menon (2009). Assumptions and references employed for the specific radiative processes implementation are specified in what follows. The multidimensional integration needed to compute the emission of the electron distributions or of the photon fields (see Appendix A) is performed with NumPy. astropy provides all of the physical constants, handles unit transformations and distance calculations. The plotting utilities for visualisation of results rely on matplotlib (Hunter 2007). The code development is hosted on GitHub² and its documentation on readthedocs³. Each agnpy release tagged in GitHub is archived on the Zenodo platform that offers a digital object identifier (DOI) for each release of the code. Version 0.1.3 (Nigro et al. 2021) was used for this publication. Like Gammapy and naima, agnpy is one of the packages affiliated with the astropy project. The input quantities and the results of all the computations of agnpy are provided as astropy Quantities (or arrays of them, as in the case of a SED), allowing a seamless integration of other package using the same library upstream or downstream of the calculations performed with agnpy.

In what follows we detail the content of the individual modules constituting the package. Figure 1 offers an overview of their relations and the computations they perform.

Fig. 1. Scheme representing agnpy’s modules, their relations and the computations they perform. agnpy.emission_regions, defining the region responsible for the non-thermal radiation (containing the EED), is the initialising argument of the classes modelling radiative processes and computing their SEDs (agnpy.synchrotron and agnpy.compton). Thermal and line emitters, collected in agnpy.targets, are initialising arguments of the modules for the external Compton and the γγ absorption calculations (agnpy.absorption). The targets’ thermal SEDs and energy densities (u) can also be computed. The emission region is passed to the class for self-consistent modelling (agnpy.constraints) and to the absorption class to compute the γγ absorption on synchrotron photons.

2.1. Emitters

agnpy contains modules describing thermal and non-thermal emitters.

agnpy.emission_regions describes the emission regions responsible for the non-thermal radiative processes. The emission region contains the EED: the distribution of Lorentz factors, γ′, of the e^± accelerated in the source, ${\underline{n}}_{\text{e}}^{\prime }({\gamma }^{\prime })[{\mathrm{cm}}^{-3}]$ (primed quantities refer to the reference frame comoving with the blob, see Appendix A for more details on the notation). At the moment the only non-thermal emission region implemented is a spherical plasmoid streaming along the jet (blob). The electron distribution is considered homogeneous and isotropic in the blob’s rest frame.

agnpy.targets describes line and thermal emitters providing the photon fields target for inverse Compton scattering or γγ absorption. The line and thermal emitters modelled in the package are: the cosmic microwave background (CMB); a monochromatic point-like source behind the jet (to be used mostly for consistency checks, see Sect. 4); a Shakura & Sunyaev (1973) accretion disc, modelled as a geometrically thin disc whose emission is parameterised following Dermer & Schlickeiser (2002) and Dermer et al. (2009); a broad line region (BLR) simplified to an infinitesimally thin sphere reprocessing a fraction of the disc radiation in a monochromatic (line) emission (see Sect 3.5 in Finke 2016, and references therein)⁴; a dust torus (DT), simplified to a ring reprocessing a fraction of the disc radiation in a single-temperature black body (BB) emission (see Sect. 3.6 of Finke 2016). The thermal emission of some of the targets can also be computed: one can evaluate the multi-temperature BB SED produced in optical-UV by the accretion disc or the IR-peaking BB SED representing the DT emission. These BB SEDs approximate the actual thermal emission and are conventionally adopted in radio-to-gamma-ray MWL modelling to check the flux level of the thermal components against the dominant non-thermal ones (see e.g., Ghisellini & Tavecchio 2009). They are inadequate for precise modelling of the SED of those sources in which the thermal components can be directly observed (in particular non-jetted AGN). The energy density, u [erg cm⁻³], of every target can also be computed as a function of the distance from the central black hole (BH), r, both in the reference frame comoving with the blob and in the one comoving with the centre of the galaxy (see Fig. 1).

2.2. Radiative processes

The package contains a module for each leptonic radiative process (see Blumenthal & Gould 1970; Rybicki & Lightman 1986; Dermer & Menon 2009, for canonical references). All the formulae used in the radiative processes calculations are provided in Appendices A.3–A.5. Here we describe the fundamental physical assumptions considered in their implementation.

agnpy.synchrotron computes synchrotron spectra (with self-absorption), following Finke et al. (2008). It is assumed that the electron distribution is immersed in a large-scale random magnetic field (Crusius & Schlickeiser 1986).

agnpy.compton models the inverse Compton (IC) scattering either of the synchrotron photons produced by the same accelerated electrons, ‘synchrotron self-Compton’ (SSC), or of the soft photon fields produced by the AGN line and thermal emitters or by the CMB, ‘external Compton’ (EC). In both cases, the full Compton cross section, including the Klein-Nishina regime, is considered. SSC spectra are computed as in Finke et al. (2008), with the target synchrotron radiation assumed uniform in the blob (Jones 1968). EC spectra are computed according to Dermer et al. (2009) and Finke (2016). The electron distribution is transformed to a frame comoving with the target photon field (following the approach of Georganopoulos et al. 2001) and then convolved with the energy distribution of the latter and the Compton cross section. The head-on approximation of the cross section is used: it is assumed that the scattered photons have roughly the same direction of the scattering electrons (as in Dermer & Schlickeiser 1993, 2002).

agnpy.absorption allows one to compute the absorption due to γγ pair production (Gould & Schréder 1967) on the soft photon fields of the line and thermal emitters, as well as on the synchrotron photons produced in the emission region. Concerning γγ pair production with the EBL, absorption values computed according to the models of Franceschini et al. (2008), Finke et al. (2010), and Domínguez et al. (2011) are available in the package resources. They are interpolated to provide an energy- and redshift-dependent opacity function.

2.3. Parameter constraints

The package does not include any routine for the solution of the differential equation describing the EED time evolution (see e.g., Cerruti 2020, Eq. (7)). Nonetheless it contains an additional module, agnpy.constraints, constraining its spectral parameters such that the modelling can be considered self-consistent. Three time scales (computed in the reference frame of the emission region and therefore indicated with primed quantities, see notation in Appendix A) are considered in the code: acceleration, cooling and dynamical or ballistic. The acceleration time scale of electrons with energy E′ is defined by the acceleration parameter ξ (see e.g., Bednarek & Sitarek 2007) as

$$\begin{array}{c}\hfill t_{\mathrm{acc}}^{\prime }=E^{\prime }{\left(\frac{\mathrm{d}{E}^{\prime }}{\mathrm{d}{t}^{\prime }}\right)}_{\mathrm{acc}},\mathrm{where}{\left(\frac{\mathrm{d}{E}^{\prime }}{\mathrm{d}{t}^{\prime }}\right)}_{\mathrm{acc}}=\xi c{E}^{\prime }/R_{\mathrm{L}},\end{array}$$(1)

where c is the speed of light and R_L is the Larmor radius of the electron. The cooling time scale due to energy losses via synchrotron radiation or inverse Compton in the Thomson regime can be computed as (Blumenthal & Gould 1970)

$$\begin{array}{c}\hfill t_{\mathrm{cool}}^{\prime }=E^{\prime }{\left(\frac{\mathrm{d}{E}^{\prime }}{\mathrm{d}{t}^{\prime }}\right)}_{\mathrm{cool}},\mathrm{where}{\left(\frac{\mathrm{d}{E}^{\prime }}{\mathrm{d}{t}^{\prime }}\right)}_{\mathrm{cool}}=\frac{4}{3}{\sigma }_{T}c{u}^{\prime }{\gamma }^{\prime 2},\end{array}$$(2)

where σ_T is the Thomson cross section and u′ is the energy density of the magnetic field (for synchrotron cooling), or the radiation field (for inverse Compton cooling in the Thomson regime), computed in the reference frame of the emission region. We note that, due to relativistic beaming, the radiation density of the photon field of emitters external to the blob has to be transformed to its comoving reference frame (see Appendix A.2 for these transformations). Finally, the ballistic or dynamic time scale is determined from the blob size via

$$\begin{array}{c}\hfill t_{\mathrm{bal}}^{\prime }=R_{\mathrm{b}}^{\prime }/c.\end{array}$$(3)

This time scale determines the escape time of the radiation from the blob. As the electrons are confined in the blob by the magnetic field they are naturally not directly affected by such an escape condition. Nevertheless Eq. (3) also determines the time scale in which the blob crosses a distance corresponding to its size. In other words, in the reference frame comoving with the blob, the ambient matter traverses a distance corresponding to the radius of the emission region. Both the escape of the radiation and the change of the blob’s surrounding conditions can cause breaks in the spectrum. It should be noticed that the exemplary broken power law implicitly assumed will not accurately describe the evolving electron distribution in all cases. In the case of hard injection spectra, fast cooling can result in an additional pile-up in the EED (see e.g., Saugé & Henri 2004)

3. Applications

In the following section we illustrate some applications of the package functionalities, performing typical estimations necessary for modelling the emission of jetted AGN.

3.1. Computing the SED for a radiative process

As stated in the introduction, the main objective of the package is the computation of the broadband SED, from radio to gamma rays, due to the radiative process in jetted AGN. In the code snippet in Fig. 2a we illustrate how, with few lines of python, one can evaluate the SED due to synchrotron radiation of a power law of electrons (Fig. 2b). SED computations for further radiative processes are illustrated in Sects. 3.5, 4.1 and 4.2.

Fig. 2. SED due to synchrotron radiation computed with agnpy. Default parameters for the emission region are given in Table B.1 (SSC model column); fluxes and frequencies are given in the observer reference frame. (a) Example script using agnpy to compute and display a synchrotron SED. (b) Output of the example script.

3.2. Emission and densities of target photon fields

When considering the inverse Compton scenario in presence of different photon targets, one might want to determine which one produces the highest energy density u [erg cm⁻³] (see Appendix A.2 for a definition, Sect. 3 of Ghisellini & Tavecchio 2009 and Sect. 8 of Finke 2016 for a more complete discussion). In the Thomson regime, for example, the higher energy density will be the one dominating the electrons’ Compton cooling, Eq. (2), as long as the radiation can be considered quasi-isotropic. agnpy estimates the energy densities of line and thermal emitters as a function of the distance from the central black hole (BH) along the jet axis, u(r), in the reference frames comoving with the galaxy or the blob. Figure 3 illustrates, as an example, the energy densities of the photon fields generated by the CMB, a Shakura-Sunyaev disc, a spherical shell BLR emitting the Lyα line, and a ring DT, whose parameters are given in Table B.2. The energy densities are computed in the reference frame comoving with the emission region (EC model column in Table B.1). The energy density provided by the magnetic field of the blob, U_B = B²/(8π), and by the synchrotron radiation are plotted for comparison⁵. We observe that the energy densities of the AGN thermal and line emitters drop as we move at a distance larger than their typical sizes (for target parameters refer to Table B.2). For distances even larger than the dust torus radius r > 10²⁰ cm ≈ 32 pc the only relevant external photon field is the CMB, though, in the particular case realised in Fig. 3, the synchrotron radiation will be the one dominating the cooling.

Fig. 3. Energy densities of several photon fields: accretion disc (red solid line), BLR (blue solid line), DT (yellow solid line), CMB (black solid line), as a function of the distance from the central BH, computed with agnpy. The energy density of synchrotron radiation (green dotted line) and of the magnetic field (green dashed line) are plotted for comparison. The emitters parameters are specified in Table B.2. A BLR emitting the Lyα line is considered. The energy densities are expressed in the reference frame comoving with the second blob in Table B.1 (EC model).

3.3. Absorption

High-energy photons might be absorbed via γγ pair production by the very same soft photon fields that constitute the target for Compton scattering, or by the EBL on their path to Earth. agnpy estimates the γγ opacity, τ_γγ, due to the soft photon fields of synchrotron radiation, of the BLR and the DT, and of a monochromatic point source behind the jet (though the latter is used only for cross-checks, as illustrated in Sect. 4). As an illustrative example, Fig. 4 displays the total γγ opacity, computed with agnpy, due to photon fields internal and external to the blob. For the internal photon fields, the opacity due to synchrotron radiation of the second blob in Table B.1 (EC model) is represented, multiplied by a factor 10³ to be comparable with the opacity due to the external photon fields. For the latter, a DT and a BLR constituted of two concentrical shells, one emitting the Lyα and the other the Hα line, are considered. The luminosity and radius of the Hα shell are obtained from those of the Lyα, following the scaling of the stratified BLR model in Finke (2016). Parameters of line and thermal emitters are given in Table B.2. The blob is placed within the BLR, at a distance r = 1.1 × 10¹⁶ cm = 0.1 × R_Lyα. We note that internal and external absorption correspond to different attenuation factors of the final photon spectrum (see Appendix A.5). It should be noticed that in the case of high opacity for high-energy photons, the secondary e^± pairs can produce further IC emission in the ambient radiation field (see e.g., Bednarek & Kirk 1995; Blandford & Levinson 1995). This cascading effect is not calculated by agnpy.

Fig. 4. Total γγ opacity (black solid line and grey shaded area) produced by photon field external to the blob: BLR (red and blue dashed lines) and DT (yellow dashed line). The BLR is assumed to be composed of two concentric shells one emitting the Lyα (red) and the other the Hα line (blue). Line and thermal emitters parameters are specified in Table B.2. The green line and shaded area represent the opacity due to the synchrotron photons produced by the second blob in Table B.1 (EC model). The latter is multiplied by 103 to be comparable with the opacity due to the external photon fields. The blob is placed well within the BLR, at r = 1.1 × 1016 cm = 0.1 × RLyα. Frequencies are given in the observer reference frame.

The package resources additionally include tables of the absorption value, exp(−τ_γγ), computed according to the EBL models of Franceschini et al. (2008), Finke et al. (2010), and Domínguez et al. (2011). These are interpolated and made available to the user to be evaluated, for a given redshift, over an array of frequencies.

3.4. Constraining the Lorentz factor

The module for constraining the spectral parameters provides a number of functions that compute the Lorentz factors of the electrons for which the different time scales listed in Eqs. (1)–(3) match. Typically, a maximum Lorentz factor, ${\gamma }_{\max }^{\prime }$, is obtained from the comparison between a given cooling and acceleration time scales. A break Lorentz factor, ${\gamma }_{\text{b}}^{\prime }$, corresponding to a change in slope in the EED, is instead obtained equating the cooling and the ballistic time scale. A practical use case would be to compute the ${\gamma }_{\max }^{\prime }$ values for all the radiative processes considered in a particular scenario and then choose the lowest one as the most constraining. For example, considering the blob for the EC model in Table B.1, and the synchrotron, SSC, and EC on DT radiative processes, we obtain: ${\gamma }_{\text{max},\hspace{0.33em}\text{synch}}^{\prime }=1.6\times {10}^8$, ${\gamma }_{\text{max},\hspace{0.33em}\text{SSC}}^{\prime }=2.2\times {10}^{10}$, and ${\gamma }_{\text{max},\hspace{0.33em}\text{EC}}^{\prime }=2.4\times {10}^8$. We would then choose ${\gamma }_{\text{max},\hspace{0.33em}\text{synch}}^{\prime }$ as a constraint on the maximum energy. The process responsible for the limit on the maximum energy of the electrons will likely produce a break in the electron spectrum connected with the change of the conditions inside the blob if the corresponding ${\gamma }_{\text{b}}^{\prime }$ is within the $[{\gamma }_{\min }^{\prime },{\gamma }_{\max }^{\prime }]$ range, as it is in this case: ${\gamma }_{\text{b},\hspace{0.33em}\text{synch}}^{\prime }=7.4\times {10}^3$. Both spectral parameters are in fact in rough agreement with those chosen for the broken power-law EED in the EC model in the table.

3.5. Fitting the SED of blazars

agnpy is a package for numerical modelling; routines for data handling and fitting are, therefore, outside of its scope. Moreover, such functionalities are available in many other packages and re-implementing them will, in addition to creating duplication, break the modular approach we discussed in the introduction. To ease the interface between agnpy and other high-level data-analysis package we expose, for each radiative process, a function computing the SED that depends on the frequency and on all the model parameters simultaneously (i.e. νF_ν(ν; Θ), where Θ is a vector of the model parameters one wants to determine through some statistical procedure). Such functions can be easily wrapped by any fitting routine, leaving substantial freedom to the user to create an optimal fit model from any arbitrary combination of emission regions and radiative processes. Once the model is defined, users can choose which parameters of the model have to be fixed and which one left free, depending on the completeness of their data set and on information already available on the source. We provide, in the online documentation tutorials and in the code associated with this publication, two examples using both Gammapy (Deil et al. 2017) version 0.18.2 (Deil et al. 2021) and sherpa (Freeman et al. 2001; Doe et al. 2007) version v4.14.0 to wrap agnpy’s functions to fit the observed MWL SED of the blazars Mrk 421 (a BL Lac) and PKS 1510-089 (a flat spectrum radio quasar, FSRQ). While the low energy continuum of both source types is interpreted as due to synchrotron radiation of electrons, the high-energy one is commonly explained as due to SSC for BL Lacs and to EC (with a possible additional contribution from SSC process) for FSRQs. These two science cases offer us the possibility to test the fitting procedure with most of the radiative processes implemented in the package.

3.5.1. Mrk 421

We modelled the radio-to-gamma-ray broadband SED of Mrk 421 observed in 2009 by the MWL campaign described in Abdo et al. (2011) with the synchrotron and SSC radiation of a single blob accelerating a broken power-law EED. The model is described by nine parameters, six for the EED (k_e, p₁, p₂, ${\gamma }_{\text{b}}^{\prime }$, ${\gamma }_{\min }^{\prime }$, and ${\gamma }_{\max }^{\prime }$; see the parametrisation in Appendix A.1), and three for the blob ($R_{\text{b}}^{\prime }$, B, and δ_D). We fixed the minimum and maximum Lorentz factors of the EED, ${\gamma }_{\min }^{\prime }$ and ${\gamma }_{\max }^{\prime }$, and constrained the radius of the blob, $R_{\text{b}}^{\prime }$, with the observed variability time scale t_var = 1 d and the Doppler factor δ_D, via $R_{\text{b}}^{\prime }=c{t}_{\text{var}}{\delta }_{\text{D}}/(1+z)$. To limit the span of the parameters space, we fitted the log₁₀ of the parameters whose range is expected to cover several orders of magnitudes, such as the electron normalisation, k_e, the break Lorentz factors of the electron distribution, ${\gamma }_{\text{b}}^{\prime }$, and the magnetic field, B. Both for Gammapy and sherpa we minimised a χ² statistic. For sherpa we selected the Levenberg-Marquardt minimisation algorithm (Levenberg 1944; Marquardt 1963), while for Gammapy we used iminuit (Dembinski et al. 2020), a python wrapper for the minuit2⁶ library. The SED points below 10¹¹ Hz were not considered in the fit as they represent radio flux measurements with a larger integration region and are contaminated by the extended jet emission. We assumed 30% systematic uncertainties on very-high-energy gamma-ray flux measurements (E ≥ 100 GeV, Aleksić et al. 2016), 10% on the high-energy ones (0.1 GeV ≤ E < 100 GeV, Abdo et al. 2011), 10% in the X-ray band (Madsen et al. 2017), and a conservative 5% on the low-energy fluxes, from the UV to the radio band (Bohlin et al. 2014). The values of the model parameters obtained with the different software are reported in Table 1a. The best-fit model obtained with Gammapy, illustrated in Fig. 5 along with its individual components, returned a statistics χ² = 271.2 with 80 degrees of freedom. We observe a good agreement between the fitted model and the data over the entire energy range. We refrain from providing uncertainties on the best-fit values due to the simplistic treatment of the flux points uncertainties: the systematic uncertainties we considered represent an educated guess from the literature and we also neglected any correlation between flux points obtained by the same instrument. We remark that the objective of this exercise is to illustrate the seamless interface of analysis and modelling tools and not to perform a statistical analysis nor draw any scientific conclusion from the modelling.

Fig. 5. MWL SED of Mrk 421 (black points) from Abdo et al. (2011) fitted wrapping agnpy with Gammapy. The red solid line represents the best fit model, the coloured dashed lines its individual components: yellow for synchrotron, blue for SSC. Fluxes and frequencies are given in the observer reference frame.

3.5.2. PKS 1510-089

The MWL SED of PKS 1510-089, characteristic of a high optical and gamma-ray state observed in 2015 (Ahnen et al. 2017) was modelled with the synchrotron, SSC and EC radiation of a single blob accelerating a broken power-law EED. We included, in the total flux model, the disc and the DT thermal emissions since for FSRQs, while often outshined by synchrotron emission, they are still non-negligible. We considered the DT photon field as target for the EC and neglected any γγ absorption on photon fields external to the blob. Independent studies on the location of the emission region indeed finds it to be outside the BLR, without any evidence of absorption of the gamma-ray spectrum on its photon field (Costamante et al. 2018; Meyer et al. 2019). For this reason we set the emission region location at a distance r = 6 × 10¹⁷ cm, outside the BLR (following Ahnen et al. 2017). We accounted for self-absorption of the synchrotron radiation. For the χ² minimisation we used the same settings described for Mrk 421. In this case the number of model parameters is larger than those of the simple synchrotron and SSC model, as they include the parameters of the disc and DT thermal components. We fixed all the parameters related to the thermal emitters: the disc luminosity L_disc and accretion efficiency η, the mass of its BH M_BH, its inner and outer radii R_in(out). We considered a torus with radius R_DT reprocessing a fraction ξ_DT of the disc efficiency and radiating at temperature T_DT. The values adopted, listed in Table 1b, were obtained from previous modelling of the source (Aleksić et al. 2014; Ahnen et al. 2017) following the scaling relations in Ghisellini & Tavecchio (2009). As before, we fixed ${\gamma }_{\min }^{\prime }$ and ${\gamma }_{\max }^{\prime }$ of the EED, but due to the lower number of flux measurements compared to the Mrk 421 case, we were not able to fit simultaneously the EED parameters along with δ_D and B and therefore chose to fix the value of δ_D = 25. The blob size was again constrained by the time variability, fixed to t_var = 0.5 d. As for Mrk 421, spectral points below 10¹¹ Hz were not included in the fit and we assumed for the flux points the same systematic uncertainties previously specified. The value of the model parameters obtained with the different software are reported in Table 1a. The best fit model obtained with Gammapy, illustrated in Fig. 6 along with its thermal and non-thermal individual components, returns a statistics χ² = 230.5 with 36 degrees of freedom. Again we refrain from providing errors on the parameters.

Fig. 6. MWL SED of PKS 1510-089 (black points) from Ahnen et al. (2017) fitted wrapping agnpy with Gammapy. The red solid line represents the best fit model, coloured dashed lines represent its non-thermal components: yellow for synchrotron, blue for SSC, green for EC on DT; grey lines represent its thermal components: dot-dashed for the multi-temperature BB of the accretion disc, dotted for the BB of the DT. Fluxes and frequencies are given in the observer reference frame.

4. Validation

In order to validate the calculations performed by agnpy we compared them against bibliographic references and against jetset (version 1.2.0rc10). Moreover, we implemented internal consistency checks that examine the compatibility of different implementations of a similar scenario. One example is provided by the processes dependent on the distance between the emission region (the blob) and the central BH: EC and γγ absorption on AGN photon fields. The photon field produced by any emitter centred at the BH position can be approximated, for distances much larger than its size, r ≫ R_emitter, with the photon field produced by a point source (with the same luminosity) at the BH position. For very large distances, we can therefore verify the consistency of the EC spectrum (or the γγ opacity) generated by a given line or thermal emitter with the EC spectrum (or the opacity) generated by a point-source approximating it. All the comparisons illustrated in this section are embedded in the test system part of the package continuous integration (CI) cycle. Before a change is merged in the master branch of the code repository, these benchmark SEDs are automatically generated and it is numerically assessed that their deviations from the references are within a certain accuracy margin.

4.1. Synchrotron and synchrotron self-Compton

To cross-check the synchrotron and SSC radiation calculation we reproduced with agnpy one of the models in Fig. 7.4 of Dermer & Menon (2009) (whose parameters are given under the SSC model column of Table B.1). We used jetset, which for these processes relies on the same assumptions as agnpy (homogeneous magnetic field and homogeneous synchrotron radiation), to produce the SEDs for the same model parameters. Figure 7 illustrates the comparison between the reference SEDs and agnpy. In the bottom panels the deviation of agnpy from a given reference is represented with the quantity $\frac{\nu F_{\nu ,\mathrm{agnpy}}}{\nu F_{\nu ,\mathrm{ref}.}}-1$. We observe an absolute deviation for the synchrotron and SSC processes well within 20% both from Dermer & Menon (2009) and from jetset, for most of the frequency range considered.

Fig. 7. Comparison of the synchrotron (left panels) and SSC (right panels) radiation computed with agnpy (red solid line) and jetset (blue dashed line) and taken from Fig. 7.4 of Dermer & Menon (2009) (black dashed line). The bottom panels display the deviation of agnpy from the bibliographic reference or jetset. Fluxes and frequencies are given in the observer reference frame.

4.2. External Compton

For the EC, we performed comparisons against the literature reproducing the SEDs for the scattering of the accretion disc, BLR, and DT photons in Figs. 8, 10, and 11, respectively, of Finke (2016). We compared agnpy against jetset considering EC on the aforementioned photon fields and additionally on the CMB. For the Compton scattering calculation jetset relies on the same Georganopoulos et al. (2001) approach used in agnpy, but the two software packages rely on different parametrisations and assumptions for the target photon fields, therefore larger differences are to be expected in their comparison. As already mentioned, for the EC scenario it is possible to perform an additional internal consistency check approximating, for large distances from the BH, the photon field of the line or thermal emitter with the one produced by a point-like source.

Figures 8–10 present, for each AGN photon field, two comparisons. In the first, considering a distance between the blob and the BH smaller than the emitter size, r < R_emitter, we compare agnpy against the corresponding bibliographic reference and against jetset. In the second, considering a significantly larger distance, r ≫ R_emitter, instead of the literature, we compare against the SED computed with agnpy approximating the emitter with a point-like source. We notice that the point source we consider is monochromatic (Appendix A.2.2) hence the approximation is exact only in the case in which the emitter itself is also considered monochromatic in the scattering calculation; this is true for the BLR and DT, but not for the disc. The latter has an emission dependent on the radial coordinate ϵ₀(R) (see Eq. (A.17)). In considering EC scattering on its radiation (Fig. 8, right panel) we propose a comparison against two point sources: the first emitting with the same energy emitted at the inner disc radius ϵ₀ = ϵ₀(R_in) and the second emitting with the same energy emitted at the outer disc radius ϵ₀ = ϵ₀(R_out). Figure 11 shows instead the comparison against jetset for EC on the CMB.

Fig. 8. Left: comparison, for a distance smaller than the emitter size, of the EC scattering on the disc photon field computed with agnpy (red solid line) and jetset (dashed blue line) against the results in Finke (2016) (black dashed line). Right: comparison, for a distance much larger than the emitter size, of the SED obtained with agnpy considering the full disc energy density (solid red line) or approximating it with a monochromatic point-source (black line) with the same photon energy at the inner (dashed line) or outer (dotted line) disc radius. The SED computed by jetset at the same distance is displayed with the dashed blue line. The bottom panels display the deviation of agnpy from the reference (or the approximation). Fluxes and frequencies are given in the observer reference frame.

Fig. 9. Same as in Fig. 8, but considering EC scattering on the BLR photon field.

Fig. 10. Same as in Figs. 8 and 9, but considering EC scattering on the DT photon field.

Fig. 11. Comparison of the EC scattering on the CMB at z = 1 computed with agnpy (red solid line) and jetset (blue dashed line). Fluxes and frequencies are given in the observer reference frame.

Fig. 12. Left: comparison, for a distance smaller than the emitter size, of the γγ absorption on the BLR photon field computed by agnpy (red solid line) against the results in Finke (2016) (black dashed line). Right: comparison, for a distance much larger than the emitter size, of the opacity obtained with the full calculation using the BLR energy density (solid red line) against the calculation approximating the target with a monochromatic point source (black dashed line). To make the opacity for the point-source non null a viewing angle of θs = 20° is assumed. The bottom panels display the deviation of agnpy from the reference (or the approximation). Frequencies are given in the observer reference frame.

For the EC on disc we obtain a deviation within 30%, both from the literature and the point-source approximation, for most of the frequency range considered. The accretion disc radiation in jetset is modelled following Frank et al. (2002), considering a multi-temperature black body with a temperature profile similar to the one of Dermer & Schlickeiser (2002). Deviation from jetset are larger over the entire energy range, with the peaks of the spectra deviating more than 50%. For the EC on BLR again we observe a deviation within 30% from the literature and even smaller from the point-source approximation. The BLR radiation in jetset is modelled, following Donea & Protheroe (2003), with a thick spherical shell re-emitting the BB radiation of the disc. The different BLR models result in a factor 2 difference between agnpy and jetset, over most of the energy range considered. For large distances the deviation is reduced to within 50%. Similarly to the previous cases, for the EC on DT we observe a very small deviation, mostly within 20%, from the literature and point-source approximation and a larger one, within 50%, from jetset. The latter accounts in the EC scattering for the whole BB emission of the DT (approximated instead as monochromatic in the agnpy calculation). For the EC on CMB, no literature reference was available and the two software agree within 30%.

For all the EC validations, the parameters of the blob are given under the EC model column of Table B.1 and those of the targets in Table B.2. The distances r are specified in the plot titles.

4.3. γγ absorption

Similar tests performed for the EC were repeated to validate the γγ absorption produced by the BLR and DT photon fields, though in this case no comparison with other software was possible. Again, for the distances smaller than the size of the line or thermal emitter, r < R_emitter, we compared against the corresponding bibliographic reference, that is Fig. 14 of Finke (2016). For large distances, r ≫ R_emitter, we performed an internal comparison, approximating the target with a point source. We notice that, in the case of a source aligned with the observer, the absorption from a point-source behind the jet would be zero (see Eq. (A.52)). Therefore to perform a meaningful comparison for the absorption due to a point-source behind the jet approximating the emitter, we assumed a viewing angle of θ_s = 20° in the case of large distances. Again, with the exception of the extremes of the frequency range considered, we achieve a deviation smaller than 10% in all cases. We did not encounter any bibliographic or software reference to validate the absorption on synchrotron photons.

5. Reproducibility and performance

The issue of reproducibility in science is often addressed from the point of view of data availability, reduction, and analysis (e.g., Chen et al. 2018), but never from the point of view of physical interpretation of results. Nonetheless, since the latter is dependent on increasingly complex software, it is appropriate to apply the same standards we require for modern data-analysis software to modelling tools. For this reason, beside exposing the source code and opening the development of our software, we took the effort of thoroughly validating its results against bibliographic references, other software and internal consistency checks. Our package offers a modelling tool easily interfaceable with other python-based packages for astrophysical data analysis (Gammapy, sherpa), thus providing a modular approach to perform a complete and reproducible astrophysical analysis, from data reduction to physical interpretation, in the same computational environment. In order to also make the results in this publication reproducible, we made available all the scripts used to produce them in a GitHub repository⁷, that we also archived on Zenodo (Nigro et al. 2022). We provide an anaconda environment freezing the version of the packages (and dependencies) employed and a Docker container⁸ archiving the computational environment used to produce all the results in this paper.

To provide the reader an estimate of agnpy’s performance in its current state, each script associated with this publication times the most relevant computations performed with the package functions. In Table 2 we report the execution times, t_exec, for the radiative processes computations in the validation scripts (those producing Figs. 7–13). They are obtained the on a macOS system (Monterey 12.0.1) with a 1.4 GHz processor and 8 GB RAM. For the calculation of a given processes, t_exec will depend on the number of dimensions of the integral and on the number of points chosen along each axis when the latter is discretised. We indicate in the table the number of points chosen over each dimension when computing the SED (or the absorption). The number along each dimension is chosen to obtain smooth spectra; we provide further discussion on the effects of the resolution of the discretisation on the final spectra in Appendix C. The fastest computations are given by one-dimensional integrals, t_exec ≈ 10 ms, such as those estimating the spectra due to synchrotron radiation (Eq. (A.30)) and EC on point-like source (Eq. (A.43)). Two-dimensional integrals like those computing the SSC (Eq. (A.33)) and EC on DT spectra (Eq. (A.46)) require up to 100 ms. Finally, the largest three-dimensional integrals such as the spectra of EC on disc (Eq. (A.45)), EC on BLR (Eq. (A.44)), EC on CMB (Eq. (A.42)) can take up to few s. We remark that the execution times are measured on a single run and might significantly differ from those reported for different machines. They should not differ significantly, mostly at the ms level, from execution to execution, on the same machine. The χ² minimisation for the Mrk 421 case (6 model parameters, 86 flux points) requires ≈15 s, when using sherpa and ≈120 s when using Gammapy. The χ² minimisation for the PKS 1510-089 case (5 model parameters, 41 flux points) requires ≈23 s when using sherpa and ≈85 s when using Gammapy.

Fig. 13. Same as in Fig. 12, but considering γγ absorption on the DT photon field.

6. Conclusions

In this article we presented agnpy: a python package modelling the radiative processes in jetted active galaxies. Built on NumPy, SciPy, and astropy, agnpy provides another element to the modular astrophysical software system envisioned in Portegies Zwart (2018). Not only it can be easily wrapped by other astrophysical libraries for data analysis living in the same ecosystem, such as sherpa and Gammapy, but it also complements the science cases covered by other modelling tools such as naima and jetset.

We present several examples of applications, the most useful being the calculation of the spectrum produced by a given radiative process and the computation of the energy densities and the γγ absorption on the photon fields of several AGN components. The features connected with line and thermal emitters and with absorption are not included in any other jetted AGN modelling software, nor is the possibility to describe any jetted radio galaxy (including non blazars). Even though not capable, in its current state, of properly describing the time-evolution of the EED, agnpy provides a class constraining its energies based on a simple parametrisation of the acceleration, escape, and radiation processes. In the spirit of the interfaceability previously discussed, we illustrate, in the code associated with this publication and in the online documentation, how to wrap agnpy with sherpa and Gammapy in two examples, fitting the MWL SED of the blazars Mrk 421 and PKS 1510-089.

We provide an extensively and continuously validated modelling tool. Extensively because its results have been validated against bibliographic references, other modelling software and through internal consistency checks. Continuously because the same checks illustrated in Sect. 4 constitute only a part of the tests executed during the CI cycle. For all the radiative processes considered, we obtain results that deviate between 10 and 30% from the bibliographic references and jetset, when considering the same physical assumptions. Larger differences, still within a factor 0.5 − 2, are observed when comparing with the calculations of jetset that assume a different parametrisations of the photon fields. To our knowledge, agnpy is the first open-source modelling package for radio-to-gamma-ray non-thermal emission openly providing software and bibliographic comparisons. It would be interesting, provided the availability of other authors, to expand them to results from other publications. Benchmark models for non-thermal radiative processes could be created, as already successfully done, for example, for radiative transfer of photons in dust structures (Ivezic et al. 1997; Gordon et al. 2017). Concerning performance, the most simple computations can be executed on current machines in few ms, while the most complex can take up to few s. Current performance results in a reasonable execution time, of the order of few seconds, for fitting routines wrapping the code.

Concerning future expansion of the package, we outline two possible paths. The first would be, given the recent multi-messenger observations of blazars (IceCube collaboration 2018), and the growing interest in hadronic radiative processes, their implementation in agnpy. The second might be the addition of a procedure to correctly evolve the EED properly accounting for electrons acceleration, radiation and escape (Kardashev 1962; Chiaberge & Ghisellini 1999).

With the advent of the next generation of high-energy observatories and the amount of data they will facilitate, open-source modelling tools will become essential to extending the interpretation effort to the largest possible number of astrophysicists. Preparing for widespread access to these data, high-energy astrophysicists have started to develop standardised data formats and open-source analysis tools. In the next era of computational astrophysics the same standardisation and reproducibility standards should be required from modelling software.

Acknowledgments

We are grateful to Charles Dermer and Justin Finke for patiently addressing our questions and for providing the spectral points from their publications. We are indebted with Justin Finke and Matteo Cerruti for reviewing this manuscript. We thank Jose Enrique Ruiz for providing comments on an earlier version of the paper and for helping to set up the reproducibility assets. We thank the journal referee for carefully examining the manuscript and the software and for helping us improve the paper. This work was supported by the European Commissions Horizon 2020 Program under grant agreement 824064 (ESCAPE – European Science Cluster of Astronomy & Particle Physics ESFRI Research Infrastructures), by the ERDF under the Spanish Ministerio de Ciencia e Innovación (MICINN, grant PID2019-107847RB-C41), and from the CERCA program of the Generalitat de Catalunya. J.S. and P.G. are supported by Polish Narodowe Centrum Nauki grant 2019/34/E/ST9/00224.

References

Appendix A: Notations and formulae

This appendix describes the notation and formulae used in the code implementation and in this paper. As some of these formulae are not explicit in the references given for the specific radiative process implementations (see Section 2.2), we reworked them ourselves. We collect them here hoping they will be useful for cross-checks and future work.

Following the notation in Dermer & Menon (2009) we express the energy of the electrons in terms of their Lorentz factor γ. For the photons we use a dimensionless energy, ϵ, expressed in units of electron rest-mass energy: ϵ = hν/(m_ec²), where h is the Planck constant, ν the frequency of the photon, and m_ec² the electron rest-mass energy. Among function arguments, a semicolon separates continuous from parametric dependency, for example the one-dimensional Heaviside function between the points x = a and x = b is written as H(x; a, b).

The notation we adopt makes frequent use of implicit differentials. We mark them by underlining the function symbol and moving the differential variable among the function arguments, for example the density of electrons with Lorentz factor between γ and γ + dγ will be indicated with ${\underline{n}}_{\mathrm{e}}(\gamma )\text{d}\gamma \equiv \frac{\text{d}{n}_{\mathrm{e}}}{\text{d}\gamma }\text{d}\gamma$. Similarly, the specific spectral energy density at a given distance r from a photon field, that is the density of photons with energy between ϵ and ϵ + dϵ and solid angle between Ω and Ω + dΩ, will be indicated with $\underline{u}(ϵ,\mathrm{\Omega };r)\text{d}ϵ\text{d}\mathrm{\Omega }\equiv \frac{\partial u}{\partial ϵ\partial \mathrm{\Omega }}(r)\text{d}ϵ\text{d}\mathrm{\Omega }$.

The emission region is assumed to be a spherical plasmoid (blob) with radius R_b, and tangled uniform magnetic field B. It is considered homogeneously exposed to the line and thermal emitters radiation (i.e. it is approximated as a point in the EC calculations). The blob moves along the jet axis with (constant) relativistic velocity ℬ and bulk Lorentz factor Γ, and the jet is oriented at an angle θ_s with respect to the observer line of sight. The Doppler factor, regulating the relativistic transformations to and from the blob’s reference frame, is

$$\begin{array}{c}\hfill {\delta }_{\mathrm{D}}=\frac{1}{\mathrm{\Gamma }(1-\mathcal{B}cos{\theta }_s)}.\end{array}$$(A.1)

In this appendix we express physical quantities in three reference frames. Unprimed quantities (ϵ, γ) are measured in a reference frame comoving with the galaxy and whose origin is the galactic BH. Primed quantities (ϵ′, γ′) refer instead to a reference frame comoving with the blob. Transformations of energies between these two reference frames read ϵ/ϵ′=γ/γ′=δ_D. The cosine of an angle μ, transforms instead as μ′=(μ − ℬ)/(1 − ℬμ). Hatted quantities refer to the inertial frame of the observer. An observed photon energy, $\widehat{ϵ}$, converts to the galactic reference frame via $ϵ=\widehat{ϵ}(1+z)$ and to the blob reference frame via ${ϵ}^{\prime }=\widehat{ϵ}(1+z)/{\delta }_{\mathrm{D}}$. We indicate the energy flux (or SED) with the short notation f_ϵ = νF_ν [erg cm⁻² s⁻¹] where ϵ is the dimensionless energy corresponding to ν. Observed energy fluxes (i.e. measured at Earth) relate to the intrinsic ones produced in the blob frame via ${\widehat{f}}_{\widehat{ϵ}}={\delta }_{\mathrm{D}}^4{f}_{{ϵ}^{\prime }}^{\prime }$. All of the plots in the main text of this paper present observed frequencies and fluxes; nonetheless we do not use the hatted notation not to overburden the notation appearing in them. We employ the proper notation in this appendix to mark the transformations between the different reference frames.

A.1. Non-thermal electron spectra

The electron energy density (EED) in the blob reference frame, ${\underline{n}}_{\text{e}}^{\prime }[{\mathrm{cm}}^{-3}]$, is parametrised as a function of the Lorentz factor only, γ′, thus assuming electrons with a uniform and isotropic distribution in the emission region. The total volume density of electrons is $n_{\mathrm{e},\mathrm{tot}}^{\prime }={\int }_1^{\infty }\mathrm{d}{\gamma }^{\prime }{\underline{n}}_{\mathrm{e}}^{\prime }({\gamma }^{\prime })$. The differential number of particles, ${\underline{N}}_{\text{e}}^{\prime }({\gamma }^{\prime })$, can be obtained multiplying the particle density, ${\underline{n}}_{\mathrm{e}}^{\prime }({\gamma }^{\prime })$, by the volume of the emission region, $V_{\text{b}}^{\prime }$. The total energy in electrons, $W_{\text{e}}^{\prime }$, is thus defined as

$$\begin{array}{c}\hfill W_{\mathrm{e}}^{\prime }=m_{\mathrm{e}}{c}^2{\displaystyle {\int }_1^{\infty }\mathrm{d}{\gamma }^{\prime }{\gamma }^{\prime }{\underline{N}}_{\mathrm{e}}^{\prime }({\gamma }^{\prime }).}\end{array}$$(A.2)

The simplest parametrisation available in agnpy for the electron energy density is a power law

$$\begin{array}{c}\hfill {\underline{n}}_{\mathrm{e}}^{\prime }({\gamma }^{\prime })=k_{\mathrm{e}}{\gamma }^{\prime -p}H({\gamma }^{\prime };{\gamma }_{\min }^{\prime },{\gamma }_{\max }^{\prime }),\end{array}$$(A.3)

where k_e is the normalisation constant, p the spectral index and the Heaviside function ensures null values outside the range ${\gamma }_{\text{min}}^{\prime }\le {{\gamma }^{\prime }}^{}\le {\gamma }_{\text{max}}^{\prime }$. Another parametrisation available is a power law with a spectral index changing from p₁ to p₂ at the Lorentz factor ${\gamma }_{\text{b}}^{\prime }$, commonly known as broken power law

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\underline{n}}_{\mathrm{e}}^{\prime }({\gamma }^{\prime })=k_{\mathrm{e}}& [{\left(\frac{{\gamma }^{\prime }}{{\gamma }_{\mathrm{b}}^{\prime }}\right)}^{-p_1}H({\gamma }^{\prime };{\gamma }_{\min }^{\prime },{\gamma }_{\mathrm{b}}^{\prime })\hfill \\ \hfill & +{\left(\frac{{\gamma }^{\prime }}{{\gamma }_{\mathrm{b}}^{\prime }}\right)}^{-p_2}H({\gamma }^{\prime };{\gamma }_{\mathrm{b}}^{\prime },{\gamma }_{\max }^{\prime })].\hfill \\ \hfill \end{array}\end{array}$$(A.4)

A power law with an exponential cutoff at the Lorentz factor ${\gamma }_{\text{c}}^{\prime }$ is also available

$$\begin{array}{c}\hfill {\underline{n}}_{\mathrm{e}}^{\prime }({\gamma }^{\prime })=k_{\mathrm{e}}{\gamma }^{\prime -p}exp(-\frac{{\gamma }^{\prime }}{{\gamma }_{\mathrm{c}}^{\prime }})H({\gamma }^{\prime };{\gamma }_{\min }^{\prime },{\gamma }_{\max }^{\prime }).\end{array}$$(A.5)

Finally, a curved power law, with spectral index changing logarithmically with the Lorentz factor, commonly known as log-parabola, can also be considered

$$\begin{array}{c}\hfill {\underline{n}}_{\mathrm{e}}^{\prime }({\gamma }^{\prime })=k_{\mathrm{e}}{\left(\frac{{\gamma }^{\prime }}{{\gamma }_0^{\prime }}\right)}^{-(p+q{log}_{10}({\gamma }^{\prime }/{\gamma }_0^{\prime }))}H({\gamma }^{\prime };{\gamma }_{\min }^{\prime },{\gamma }_{\max }^{\prime }).\end{array}$$(A.6)

A.2. Energy densities of line and thermal emitters

The (integral) energy density of a photon field at a given distance r above the central BH, u(r) [erg cm⁻³], (see Section 3.2) can be computed by integrating the specific spectral energy density over energy and solid angle

$$\begin{array}{c}\hfill u(r)={\displaystyle \int }\mathrm{d}ϵ{\displaystyle \int }\mathrm{d}\mathrm{\Omega }\underline{u}(ϵ,\mathrm{\Omega };r).\end{array}$$(A.7)

Let us consider the reference frame with the BH in its centre and the scattering geometry illustrated in Figure A.1. An illustrative disc-like emitter with variable radius R illuminates the blob, which is streaming along the jet axis and is located at a distance r above the BH. Photons from the emitter hit the blob with zenith angle θ and azimuth angle ϕ. The jet axis at the blob position forms a viewing angle θ_s with the line of sight of the observer. In the following section we calculate the energy densities for all the emitters (targets) implemented in agnpy, both in the galactic (u) and in the blob (u′) reference frames. The transformation laws for spectral energy densities between these two reference frames can be found in Dermer & Schlickeiser (2002, Section 4) and Dermer & Menon (2009, Chapter 5).

Fig. A.1. Geometry used for the energy density and external Compton scattering calculations. A disc-like emitter with radial coordinate R is represented as an example.

A.2.1. Isotropic monochromatic

An isotropic monochromatic radiation field with dimensionless energy ϵ₀ has specific spectral energy density

$$\begin{array}{c}\hfill \underline{u}(ϵ,\mathrm{\Omega })=\frac{u_0\delta (ϵ-{ϵ}_0)}{4\pi },\end{array}$$(A.8)

such that the energy density in the galactic frame is trivially u = u₀, while in the blob frame

$$\begin{array}{c}\hfill u^{\prime }=u_0{\mathrm{\Gamma }}^2(1+\frac{{\mathcal{B}}^2}{3}),\end{array}$$(A.9)

in both cases the integral energy density is independent from the distance, as visible in Figure 3 for the CMB case.

A.2.2. Monochromatic point source behind the blob

Let us consider a monochromatic (ϵ = ϵ₀) source of luminosity L₀ at the centre of our reference frame (the BH position). Its specific spectral energy density is

$$\begin{array}{c}\hfill \underline{u}(ϵ,\mathrm{\Omega };r)=\frac{L_0}{4\pi c{r}^2}\frac{\delta (\mu -1)}{2\pi }\delta (ϵ-{ϵ}_0),\end{array}$$(A.10)

where the condition μ = cos θ = 1 ensures that the photons are coming from behind the blob. The integral energy density in the galactic frame is

$$\begin{array}{c}\hfill u(r)=\frac{L_0}{4\pi c{r}^2},\end{array}$$(A.11)

and in the blob frame

$$\begin{array}{c}\hfill u^{\prime }(r)=\frac{1}{{\mathrm{\Gamma }}^2{(1+\mathcal{B})}^2}\frac{L_0}{4\pi c{r}^2}=\frac{u(r)}{{\mathrm{\Gamma }}^2{(1+\mathcal{B})}^2}.\end{array}$$(A.12)

Where Γ and ℬ are the Lorentz factor and velocity (in units of c) of the blob. We notice that Equations (A.11) and (A.12) can be used as a cross-check for the correct calculation of the photon energy densities produced by emitters with a more complex geometry. For a large distance from the BH, where any emitter is observed under a very small angle θ ≪ θ_s, its photon field tends to the one produced by a point-like source. Hence Equations (A.11) and (A.12) represent limits to the energy densities of any emitter, the limits being exact for monochromatic emitters.

A.2.3. Shakura-Sunyaev accretion disc

Let us consider an accretion disc with inner and outer radii R_in and R_out, respectively. From the sketch in Figure A.1 we see that the cosine of the photon incoming angle μ = cos θ and the coordinate along the disc radius R are related via

$$\begin{array}{c}\hfill \mu ={(1+\frac{R^2}{r^2})}^{-1/2},R=r\sqrt{{\mu }^{-2}-1}.\end{array}$$(A.13)

The minimum and maximum incoming angles of photons from the disc impinging on the emission region are therefore

$$\begin{array}{c}\hfill {\mu }_{\min }={(1+\frac{R_{\mathrm{out}}^2}{r^2})}^{-1/2},{\mu }_{\max }={(1+\frac{R_{\mathrm{in}}^2}{r^2})}^{-1/2}.\end{array}$$(A.14)

The specific spectral energy density of photons coming from a Shakura & Sunyaev (1973) accretion disc can be parametrised, following Dermer & Schlickeiser (2002) and Dermer et al. (2009), as:

$$\begin{array}{c}\hfill \underline{u}(ϵ,\mathrm{\Omega };r)=\frac{3GM\dot{m}}{{(4\pi )}^{2}c{R}^3}\phi (R)\delta (ϵ-{ϵ}_0(R))\end{array}$$(A.15)

where φ(R) represents the variation of energy flux along the radius

$$\begin{array}{c}\hfill \phi (R)=1-\sqrt{\frac{R_{\mathrm{in}}}{R}}\end{array}$$(A.16)

and ϵ₀(R) is the monochromatic approximation for the mean photon energy emitted from the disc at radius R

$$\begin{array}{c}\hfill {ϵ}_0(R)=2.7\times {10}^{-4}{\left(\frac{l_{\mathrm{Edd}}}{M_8\eta }\right)}^{1/4}{\left(\frac{R}{R_{\mathrm{g}}}\right)}^{-3/4}.\end{array}$$(A.17)

In the above equations, G is the gravitational constant, M = M₈ × 10⁸ M_⊙ is the mass of the black hole in solar mass units, ṁ [g s⁻¹] the black hole mass accretion rate, η the fraction of gravitational energy converted to radiant energy (L_disc = ηṁc²), l_Edd the ratio of the disc luminosity to the Eddington luminosity (L_Edd = 1.26 × 10⁴⁶ M₈ erg s⁻¹), and R_g = GM/c² is the gravitational radius of the BH. When using Equation (A.16) and Equation (A.17), we make explicit the fact that the photons emitted at a given disc radius R will have different incidence angles depending on the blob position r, by replacing φ(R)→φ(μ; r) and ϵ₀(R)→ϵ₀(μ; r).

The integral energy density in the galactic frame is

$$\begin{array}{c}\hfill u(r)=\frac{3}{8\pi c}\frac{GM\dot{m}}{r^3}{\displaystyle {\int }_{{\mu }_{\min }}^{{\mu }_{\max }}\mathrm{d}\mu \frac{\phi (\mu ;r)}{{({\mu }^{-2}-1)}^{3/2}}.}\end{array}$$(A.18)

In the blob frame the energy density is

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill u^{\prime }(r)=& \frac{3}{8\pi c}\frac{GM\dot{m}}{r^3}\hfill \\ \hfill & {\displaystyle {\int }_{{\mu }_{\min }}^{{\mu }_{\max }}\mathrm{d}\mu \frac{\phi (\mu ;r)}{{\mathrm{\Gamma }}^6{(1-\mathcal{B}\mu )}^2{(1+\mathcal{B}{\mu }^{\prime })}^4{({\mu }^{-2}-1)}^{-3/2}}.}\hfill \\ \hfill \end{array}\end{array}$$(A.19)

A.2.4. Spherical shell broad line region

Following Finke (2016), we model the BLR as an infinitesimally thin, monochromatic shell. A complex model of a BLR emitting several lines can be built as a superposition of these simpler objects. Figure A.2 presents the geometry of a spherically symmetric emitter reprocessing the accretion disc radiation. The cosine of the incidence angle of photons produced on the surface of the reprocessing material and crossing the jet at a distance r from the BH is μ_* = cos θ_*. The spherical coordinates (r_re, μ_re, ϕ_re) represent the reprocessing material (with μ_re = cos θ_re). From geometry

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mu }_{\ast }^2& =1-{\left(\frac{r_{\mathrm{re}}}{x}\right)}^2(1-{\mu }_{\mathrm{re}}^2),\hfill \\ \hfill x^2& =r_{\mathrm{re}}^2+r^2-2r{r}_{\mathrm{re}}{\mu }_{\mathrm{re}}.\hfill \end{array}\end{array}$$(A.20)

Fig. A.2. Geometry used for the energy density and external Compton scattering calculations in case of a spherically symmetric emitter reprocessing the accretion disc radiation.

The specific spectral energy density of a spherical shell BLR reads (Finke 2016)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \underline{u}(ϵ,\mathrm{\Omega };r)=\frac{{\xi }_{\mathrm{line}}{L}_{\mathrm{disc}}}{{(4\pi )}^{2}c}\delta (ϵ-{ϵ}_{\mathrm{line}}){\displaystyle {\int }_{-1}^1\frac{\mathrm{d}{\mu }_{\mathrm{re}}}{x^2}\delta (\mu -{\mu }_{\ast }),}\end{array}\end{array}$$(A.21)

where ξ_line is the fraction of the disc luminosity (L_disc) reprocessed by the BLR, ϵ_line is the dimensionless energy of the emitted spectral line, μ_* and μ_re are obtained from Equation (A.20) with r_re = R_line, the radius at which we assume the monochromatic line is emitted. The integral energy density in the galactic frame is

$$\begin{array}{c}\hfill u(r)=\frac{{\xi }_{\mathrm{line}}{L}_{\mathrm{disc}}}{8\pi c}{\displaystyle {\int }_{-1}^1\frac{\mathrm{d}{\mu }_{\mathrm{re}}}{x^2}.}\end{array}$$(A.22)

The monochromatic spherical shell BLR can be approximated, at a large distance from the BH, as a point source behind the blob, so Equation (A.22) should reduce to Equation (A.11) in the limit r ≫ R_line. As x² → r²

$$\begin{array}{c}\hfill u(r)=\frac{{\xi }_{\mathrm{line}}{L}_{\mathrm{disc}}}{8\pi c}{\displaystyle {\int }_{-1}^1\frac{\mathrm{d}{\mu }_{\mathrm{re}}}{r^2}=\frac{{\xi }_{\mathrm{line}}{L}_{\mathrm{disc}}}{4\pi c{r}^2},}\end{array}$$(A.23)

which is Equation (A.11), the energy density of a monochromatic point source behind the blob, with luminosity L₀ = ξ_lineL_disc. The integrated energy density in the blob frame reads

$$\begin{array}{c}\hfill u^{\prime }(r)=\frac{{\xi }_{\mathrm{line}}{L}_{\mathrm{disc}}}{8\pi c}{\mathrm{\Gamma }}^2{\displaystyle {\int }_{-1}^1{(1-\mathcal{B}{\mu }_{\ast })}^2\frac{\mathrm{d}{\mu }_{\mathrm{re}}}{x^2}.}\end{array}$$(A.24)

We can check again that in the limit of large distance the latter reduces to the equivalent expression for the point-source behind the blob, Equation (A.12). For r ≫ R_line, x² → r² and μ_* → 1, therefore

$$\begin{array}{c}\hfill u^{\prime }(r)={\mathrm{\Gamma }}^2{(1-\mathcal{B})}^2\frac{{\xi }_{\mathrm{line}}{L}_{\mathrm{disc}}}{8\pi c}\frac{2}{r^2}=\frac{1}{{\mathrm{\Gamma }}^2{(1+\mathcal{B})}^2}\frac{{\xi }_{\mathrm{line}}{L}_{\mathrm{disc}}}{4\pi c{r}^2}.\end{array}$$(A.25)

We have obtained Equation (A.12), again with L₀ = ξ_lineL_disc.

A.2.5. Ring dust torus

Following Finke (2016), we model the DT as a geometrically thin ring of radius R_DT. The incidence angle and the distance of the blob from the reprocessing material in Equation (A.20) can be simplified to the ring case by imposing θ_re = π/2 and fixing the radial coordinate to the ring radius: r_re = R_DT

$$\begin{array}{cc}\hfill {\mu }_{\ast }& =\frac{r}{x},\hfill \\ \hfill x^2& =R_{\mathrm{DT}}^2+r^2.\hfill \end{array}$$(A.26)

The specific spectral energy density for a ring DT reads (Finke 2016)

$$\begin{array}{c}\hfill \underline{u}(ϵ,\mathrm{\Omega };r)=\frac{{\xi }_{\mathrm{DT}}{L}_{\mathrm{disc}}}{8{\pi }^{2}c{x}^2}\delta (\mu -\frac{r}{x})\delta (ϵ-2.7\mathrm{\Theta }),\end{array}$$(A.27)

where ξ_DT is the fraction of the disc luminosity (L_disc) reprocessed by the DT and we are approximating the torus black-body thermal emission as monochromatic at the peak energy ϵ_DT = 2.7Θ. The dimensionless temperature Θ = k_BT_DT/(m_ec²), with k_B Boltzmann constant and T_DT the peak temperature of the thermal emission. The integral energy density in the galactic frame is

$$\begin{array}{c}\hfill u(r)=\frac{{\xi }_{\mathrm{DT}}{L}_{\mathrm{disc}}}{4\pi c{x}^2}.\end{array}$$(A.28)

This can be straightforwardly reduced to Equation (A.11) (with L₀ = ξ_DTL_disc) for r ≫ R_DT (in which case x² → r²). In the frame comoving with the blob the energy density reads

$$\begin{array}{c}\hfill u^{\prime }(r)={\mathrm{\Gamma }}^2{(1-\mathcal{B}r/x)}^2\frac{{\xi }_{\mathrm{DT}}{L}_{\mathrm{disc}}}{4\pi c{x}^2}.\end{array}$$(A.29)

This expression can also be easily reduced to the case of the point-source behind the blob, Equation (A.12), considering that in the limit of large distances $\frac{x}{r}\to 1$ and Γ²(1 − ℬ)² = 1/[Γ²(1 + ℬ)²].

A.3. Synchrotron radiation

The assumptions used in the synchrotron radiation calculation are detailed in 2.2; here we report the formulae implemented in agnpy.

A.3.1. Spectral energy distribution

The synchrotron SED is computed, following Dermer & Menon (2009) and Finke et al. (2008), as

$$\begin{array}{c}\hfill {\widehat{f}}_{\widehat{ϵ}}^{\mathrm{synch}}=\frac{{\delta }_{\mathrm{D}}^4}{4\pi d_{\mathrm{L}}^2}{ϵ}^{\prime }\frac{\sqrt{3}{\mathrm{e}}^{3}B}{h}{\displaystyle {\int }_1^{\infty }\mathrm{d}{\gamma }^{\prime }{\underline{N}}_{\mathrm{e}}^{\prime }({\gamma }^{\prime })R(x),}\end{array}$$(A.30)

where d_L is the luminosity distance of the source, ϵ′ the energy of the emitted synchrotron photons in the blob reference frame, e the electron charge, and B the magnetic field. R(x) represents the pitch-angle-averaged synchrotron spectral power of a single electron, expressed as a function of the variable $x=(4\pi {ϵ}^{\prime }{m}_{\text{e}}^2{c}^3)/(3\text{e}Bh{\gamma }^{\prime 2})$. In agnpy the approximation in Eq. (D7) of Aharonian et al. (2010) is used for R(x).

A.3.2. Synchrotron self-absorption

The synchrotron self-absorption (SSA) effect (see Dermer & Menon 2009, Section 7.8) is taken into account by mutliplying the SED in A.30 by the absorption factor

$$\begin{array}{c}\hfill \frac{3}{{\tau }_{\mathrm{SSA}}}(\frac{1}{2}+\frac{exp(-{\tau }_{\mathrm{SSA}})}{{\tau }_{\mathrm{SSA}}}-\frac{1-exp(-{\tau }_{\mathrm{SSA}})}{{\tau }_{\mathrm{SSA}}^2})\end{array}$$(A.31)

where the energy-dependent SSA opacity is

$$\begin{array}{c}\hfill {\tau }_{\mathrm{SSA}}({ϵ}^{\prime })=-\frac{R_{\mathrm{b}}^{\prime }}{8\pi m_{\mathrm{e}}{ϵ}^{\prime 2}}{\left(\frac{{\lambda }_{\mathrm{C}}}{c}\right)}^3{\displaystyle {\int }_1^{\infty }\mathrm{d}{\gamma }^{\prime }R(x)\left[{\gamma }^{\prime 2}\frac{\partial }{\partial {\gamma }^{\prime }}\left(\frac{{\underline{n}}_{\mathrm{e}}^{\prime }}{{\gamma }^{\prime 2}}\right)\right]}\end{array}$$(A.32)

and where λ_C = h/(m_ec) is the Compton wavelength.

A.4. Inverse Compton radiation

As for the synchrotron case, we refer to 2.2 for the assumption used in the inverse Compton calculation and provide here the formulae implemented in agnpy. Some detail is given here to mathematically formalise the assumption going in the external Compton calculations.

A.4.1. Synchrotron self-Compton

The synchrotron self-Compton SED is computed, following Dermer & Menon (2009) and Finke et al. (2008), as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\widehat{f}}_{{\widehat{ϵ}}_{\mathrm{s}}}^{\mathrm{SSC}}=& \frac{{\delta }_{\mathrm{D}}^4}{4\pi d_{\mathrm{L}}^2}\frac{3}{4}c{\sigma }_{\mathrm{T}}{ϵ}_{\mathrm{s}}^{\prime 2}{\displaystyle {\int }_0^{\infty }\mathrm{d}{ϵ}^{\prime }\frac{{\underline{u}}_{\mathrm{synch}}^{\prime }({ϵ}^{\prime })}{{ϵ}^{\prime 2}}}\hfill \\ \hfill & {\displaystyle {\int }_0^{\infty }\mathrm{d}\gamma \frac{{\underline{N}}_{\mathrm{e}}^{\prime }({\gamma }^{\prime })}{{\gamma }^{\prime 2}}{F}_{\mathrm{C}}(q^{\prime },{\mathrm{\Gamma }}_{\mathrm{e}}^{\prime }),}\hfill \end{array}\end{array}$$(A.33)

here σ_T is the Thomson cross section, ϵ′ is the energy of the synchrotron photon, and ${ϵ}_{\text{s}}^{\prime }$ the energy of the scattered photon, both in the blob reference frame. The specific spectral energy density of synchrotron photons ${\underline{u}}_{\mathrm{synch}}^{\prime }$ can be obtained from Equation (A.30) as

$$\begin{array}{c}\hfill {\underline{u}}_{\mathrm{synch}}^{\prime }=\frac{3}{4}\frac{3d_{\mathrm{L}}^2{\widehat{f}}_{\widehat{ϵ}}^{\mathrm{synch}}}{c{R}_{\mathrm{b}}^{\prime 2}{\delta }_{\mathrm{D}}^4{ϵ}^{\prime }};\end{array}$$(A.34)

the factor 3/4 is due to averaging the radiation in a sphere. F_C is an integration kernel representing the Compton cross section for electrons and photons with a uniform spatial distribution (Jones 1968; Blumenthal & Gould 1970)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill F_{\mathrm{C}}(q^{\prime },{\mathrm{\Gamma }}_{\mathrm{e}}^{\prime })=& [2q^{\prime }ln{q}^{\prime }+(1+2q^{\prime })(1-q^{\prime })\hfill \\ \hfill & +\frac{1}{2}\frac{{({\mathrm{\Gamma }}_{\mathrm{e}}^{\prime }{q}^{\prime })}^2}{(1+{\mathrm{\Gamma }}_{\mathrm{e}}^{\prime }{q}^{\prime })}(1-q^{\prime })]H(q^{\prime };\frac{1}{4{\gamma }^{\prime 2}},1)\hfill \\ \hfill \end{array}\end{array}$$(A.35)

with

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{e}}^{\prime }=4{ϵ}_{\mathrm{s}}^{\prime }{\gamma }^{\prime }\mathrm{and}{q}^{\prime }=\frac{{ϵ}_{\mathrm{s}}^{\prime }/{\gamma }^{\prime }}{{\mathrm{\Gamma }}_{\mathrm{e}}^{\prime }(1-{ϵ}_{\mathrm{s}}^{\prime }/{\gamma }^{\prime })},\end{array}$$(A.36)

and the Heaviside functions on q′ expresses kinematic limits on ϵ′ and γ′.

A.4.2. External Compton

In order to simplify the calculation of the Compton spectra for photon targets external to the blob, we employ the ’head-on’ approximation (Reynolds 1982; Dermer & Schlickeiser 1993; Dermer et al. 2009): we assume that γ′≫1 and that, due to relativistic transformations of the angles, photons are scattered back approximately in the same direction as the incident electrons. Under this condition

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\sigma }_{\mathrm{C}}}{\mathrm{d}{ϵ}_{\mathrm{s}}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{s}}}\approx \frac{\mathrm{d}{\sigma }_{\mathrm{C}}}{\mathrm{d}{ϵ}_{\mathrm{s}}}\delta ({\mathrm{\Omega }}_{\mathrm{s}}-{\mathrm{\Omega }}_{\mathrm{e}}),\end{array}$$(A.37)

where ϵ_s is the scattered photon energy, Ω_s its direction that the delta ensures to be equal to the electron direction, Ω_e. To compute the Compton spectrum, the electron and photon densities have to be transformed to the same reference frame and then folded with the cross section. In agnpy, we follow the approach proposed by Georganopoulos et al. (2001) and convert the electron distribution to the galactic frame. For an isotropic electron distribution, ${\underline{N}}_{\mathrm{e}}(\gamma ,{\mathrm{\Omega }}_{\mathrm{e}})={\delta }_{\mathrm{D}}^3{\underline{N}}_{\text{e}}^{\prime }(\gamma /{\delta }_{\mathrm{D}})$. Given these assumptions, the energy flux produced by the blob electrons Compton-scattering an external radiation field reads (Dermer & Menon 2009)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\widehat{f}}_{{\widehat{ϵ}}_{\mathrm{s}}}^{\mathrm{EC}}=& \frac{c\pi r_{\mathrm{e}}^2}{4\pi d_{\mathrm{L}}^2}{ϵ}_{\mathrm{s}}^2{\delta }_{\mathrm{D}}^3{\displaystyle \int \mathrm{d}\mathrm{\Omega }{\int }_0^{{ϵ}_{\mathrm{high}}}\mathrm{d}ϵ\frac{\underline{u}(ϵ,\mathrm{\Omega };r)}{{ϵ}^2}}\hfill \\ \hfill & {\displaystyle {\int }_{{\gamma }_{\mathrm{low}}}^{\infty }\mathrm{d}\gamma \frac{{\underline{N}}_{\mathrm{e}}^{\prime }(\gamma /{\delta }_{\mathrm{D}})}{{\gamma }^2}{\mathrm{\Xi }}_{\mathrm{C}},}\hfill \end{array}\end{array}$$(A.38)

where r_e is the electron radius, ϵ_s the dimensionless energy of the scattered photon in the galactic frame and $\underline{u}$ the specific spectral energy density of the considered photon field. The Compton cross-section in the head-on approximation reduces to the integration kernel

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathrm{\Xi }}_{\mathrm{C}}=& y+y^{-1}+\frac{2{ϵ}_{\mathrm{s}}}{\gamma \overline{ϵ}y}+{\left(\frac{{ϵ}_{\mathrm{s}}}{\gamma \overline{ϵ}y}\right)}^2,\hfill \\ \hfill \mathrm{with}y=& 1-\frac{{ϵ}_{\mathrm{s}}}{\gamma },\hfill \\ \hfill \mathrm{and}\overline{ϵ}=& ϵ\gamma (1-cos\psi ),\hfill \\ \hfill \end{array}\end{array}$$(A.39)

where cos ψ is the angle between the direction of the incident photon and the scattering electron (see Figure A.1)

$$\begin{array}{c}\hfill cos\psi =\mu {\mu }_{\mathrm{s}}+\sqrt{1-{\mu }^2}\sqrt{1-{\mu }_{\mathrm{s}}^2}cos\varphi .\end{array}$$(A.40)

The extremes of integration on the target energy ϵ_high and Lorentz factor of the electron γ_low are imposed by kinematic limits

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {ϵ}_{\mathrm{high}}& =\frac{2{ϵ}_{\mathrm{s}}}{1-cos\psi },\hfill \\ \hfill {\gamma }_{\mathrm{low}}& =\frac{{ϵ}_{\mathrm{s}}}{2}[1+\sqrt{1+\frac{2}{ϵ{ϵ}_{\mathrm{s}}(1-cos\psi )}}].\hfill \end{array}\end{array}$$(A.41)

Considering the specific spectral energy densities $\underline{u}$ given in Section A.2, we compute the Compton-scattered spectra corresponding to each of the targets implemented in agnpy.

If the target for external Compton scattering is an isotropic monochromatic radiation field, then plugging Equation (A.8) in (A.38), we obtain

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\widehat{f}}_{{\widehat{ϵ}}_{\mathrm{s}}}^{\mathrm{EC}}=& \frac{3}{2^7{\pi }^2}\frac{c{\sigma }_{\mathrm{T}}{u}_0}{d_{\mathrm{L}}^2}{\left(\frac{{ϵ}_{\mathrm{s}}}{{ϵ}_0}\right)}^2{\delta }_{\mathrm{D}}^3{\displaystyle {\int }_0^{2\pi }\mathrm{d}\varphi {\int }_{-1}^1\mathrm{d}\mu }\hfill \\ \hfill & {\displaystyle {\int }_{{\gamma }_{\mathrm{low}}}^{\infty }\mathrm{d}\gamma \frac{{\underline{N}}_{\mathrm{e}}^{\prime }(\gamma /{\delta }_{\mathrm{D}})}{{\gamma }^2}{\mathrm{\Xi }}_{\mathrm{C}}.}\hfill \\ \hfill \end{array}\end{array}$$(A.42)

For a point source behind the blob, $\underline{u}$ is given by Equation (A.10), and the Compton-scattered SED reads

$$\begin{array}{c}\hfill {\widehat{f}}_{{\widehat{ϵ}}_{\mathrm{s}}}^{\mathrm{EC}}(r)=\frac{3}{2^7{\pi }^2}\frac{{\sigma }_{\mathrm{T}}{L}_0}{d_{\mathrm{L}}^2{r}^2}{\left(\frac{{ϵ}_{\mathrm{s}}}{{ϵ}_0}\right)}^2{\delta }_{\mathrm{D}}^3{\displaystyle {\int }_{{\gamma }_{\mathrm{low}}}^{\infty }\mathrm{d}\gamma \frac{{\underline{N}}_{\mathrm{e}}^{\prime }(\gamma /{\delta }_{\mathrm{D}})}{{\gamma }^2}{\mathrm{\Xi }}_{\mathrm{C}},}\end{array}$$(A.43)

where r indicates the distance of the blob from the central BH (Figure A.1) and we remark the dependence of the SED on this quantity in the left side of the equation. We note that in this case μ = 1 ⇒ cos ψ = μ_s. Considering a Shakura-Sunyaev accretion disc as a target, $\underline{u}$ is Equation (A.15), one derives

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\widehat{f}}_{{\widehat{ϵ}}_{\mathrm{s}}}^{\mathrm{EC}}(r)=& \frac{3^2}{2^9{\pi }^3}\frac{{\sigma }_{\mathrm{T}}GM\dot{m}}{d_{\mathrm{L}}^2{r}^3}{ϵ}_{\mathrm{s}}^2{\delta }_{\mathrm{D}}^3{\displaystyle {\int }_0^{2\pi }\mathrm{d}\varphi }\hfill \\ \hfill & {\displaystyle {\int }_{{\mu }_{\min }}^{{\mu }_{\max }}\mathrm{d}\mu \frac{\phi (\mu ;r)}{{ϵ}_0^2(\mu ;r){({\mu }^{-2}-1)}^{3/2}}}\hfill \\ \hfill & {\displaystyle {\int }_{{\gamma }_{\mathrm{low}}}^{\infty }\mathrm{d}\gamma \frac{{\underline{N}}_{\mathrm{e}}^{\prime }(\gamma /{\delta }_{\mathrm{D}})}{{\gamma }^2}{\mathrm{\Xi }}_{\mathrm{C}},}\hfill \\ \hfill \end{array}\end{array}$$(A.44)

where the extremes of integration on the cosine of the zenith are given in Equation (A.14) and once more we have made explicit that the photons emitted at a given disc radius R will have different incidence angles, depending on the blob position r, by replacing φ(R)→φ(μ; r) and ϵ₀(R)→ϵ₀(μ; r). The Compton kernel Ξ_C is dependent on μ through cos ψ and through the target energy ϵ = ϵ₀(μ; r). The external Compton spectrum for a spherical shell BLR target photon field, $\underline{u}$ in Equation (A.21), is

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\widehat{f}}_{{\widehat{ϵ}}_{\mathrm{s}}}^{\mathrm{EC}}(r)=& \frac{3}{2^9{\pi }^3}\frac{{\sigma }_{\mathrm{T}}{\xi }_{\mathrm{line}}{L}_{\mathrm{disc}}}{d_{\mathrm{L}}^2}{\left(\frac{{ϵ}_{\mathrm{s}}}{{ϵ}_{\mathrm{line}}}\right)}^2{\delta }_{\mathrm{D}}^3{\displaystyle {\int }_0^{2\pi }\mathrm{d}\varphi {\int }_{-1}^1\frac{\mathrm{d}{\mu }_{\mathrm{re}}}{x^2}}\hfill \\ \hfill & {\displaystyle {\int }_{{\gamma }_{\mathrm{low}}}^{\infty }\mathrm{d}\gamma \frac{{\underline{N}}_{\mathrm{e}}^{\prime }(\gamma /{\delta }_{\mathrm{D}})}{{\gamma }^2}{\mathrm{\Xi }}_{\mathrm{C}},}\hfill \end{array}\end{array}$$(A.45)

where now the angle between the scattering electron and the target photon is $cos\psi ={\mu }_{*}{\mu }_{\mathrm{s}}+\sqrt{1-{\mu }_{*}^2}\sqrt{1-{\mu }_{\mathrm{s}}^2}cos\varphi$ with μ_* depending on the zenith angle of the reprocessing material μ_re according to Equation (A.20). Last, considering the photon field of a ring dust torus, $\underline{u}$ in Equation (A.27), we obtain the Compton-scattered SED

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\widehat{f}}_{{\widehat{ϵ}}_{\mathrm{s}}}^{\mathrm{EC}}(r)=& \frac{3}{2^8{\pi }^3}\frac{{\sigma }_{\mathrm{T}}{\xi }_{\mathrm{DT}}{L}_{\mathrm{disc}}}{d_{\mathrm{L}}^2{x}^2}{\left(\frac{{ϵ}_{\mathrm{s}}}{{ϵ}_{\mathrm{DT}}}\right)}^2{\delta }_{\mathrm{D}}^3{\displaystyle {\int }_0^{2\pi }\mathrm{d}\varphi }\hfill \\ \hfill & {\displaystyle {\int }_{{\gamma }_{\mathrm{low}}}^{\infty }\mathrm{d}\gamma \frac{{\underline{N}}_{\mathrm{e}}^{\prime }(\gamma /{\delta }_{\mathrm{D}})}{{\gamma }^2}{\mathrm{\Xi }}_{\mathrm{C}},}\hfill \\ \hfill \end{array}\end{array}$$(A.46)

where now $cos\psi =\frac{r}{x}{\mu }_{\mathrm{s}}+\sqrt{1-{\left(\frac{r}{x}\right)}^2}\sqrt{1-{\mu }_{\mathrm{s}}^2}cos\varphi$.

A.5. γγ absorption

The optical depth for γγ absorption for a photon with observed frequency ${\widehat{\nu }}_1$ on a target photon field with specific spectral energy density $\underline{u}(ϵ,\mathrm{\Omega };r)$ is

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\tau }_{\gamma \gamma }({\widehat{\nu }}_1)=& {\displaystyle {\int }_0^{\infty }\mathrm{d}l{\int }_0^{2\pi }\mathrm{d}\varphi {\int }_{-1}^1\mathrm{d}\mu (1-cos\psi )}\hfill \\ \hfill & {\displaystyle {\int }_0^{\infty }\mathrm{d}ϵ\frac{\underline{u}(ϵ,\mathrm{\Omega };l)}{ϵm_{\mathrm{e}}{c}^2}{\sigma }_{\gamma \gamma }(s),}\hfill \\ \hfill \end{array}\end{array}$$(A.47)

where cos ψ is the cosine of the angle between the emitted and the absorbing photon (same as in Equation (A.40)), ϕ and μ the polar angles over which we integrate, and l represents the distance from the emission region to the observer along which the absorption takes place. σ_γγ(s) represents the γγ absorption cross section

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\sigma }_{\gamma \gamma }(s)& =\frac{1}{2}\pi r_{\mathrm{e}}^2(1-{\beta }_{\mathrm{cm}}^2)[(3-{\beta }_{\mathrm{cm}}^4)ln\left(\frac{1+{\beta }_{\mathrm{cm}}}{1-{\beta }_{\mathrm{cm}}}\right)\hfill \\ \hfill & -2{\beta }_{\mathrm{cm}}(2-{\beta }_{\mathrm{cm}}^2)]\hfill \\ \hfill \end{array}\end{array}$$(A.48)

where r_e is the electron radius, ${\beta }_{\mathrm{cm}}=\sqrt{1-s^{-1}}$ and s = ϵ₁ϵ (1 − cos ψ) / 2 are kinematic variables representing the momentum of the produced e^± pair in the centre-of-momentum frame. ϵ₁ is the dimensionless energy of the photon emitted in the galactic frame (whose absorption we want to estimate), ϵ the dimensionless energy of the target photon. Both in Dermer et al. (2009) and Finke (2016) the integral in Equation (A.47) is simplified assuming that the emitted photons travel in the direction parallel to the jet axis (μ_s → 1 ⇒ cos ψ = μ), decoupling the cross section and the (1 − cos ψ) term from the integral on ϕ. The optical depths thus calculated are valid only for blazars (characterised by a small viewing angle to the observer line of sight). agnpy carries out the full integration, computing optical depths valid for any jetted AGN. For speed reasons, agnpy falls back to the simplified calculations in Dermer et al. (2009) and Finke (2016) in the case in which μ_s = 1 is specified.

A.5.1. Absorption on targets

Figure A.3 illustrates the geometry used for the opacity calculation in case of sources misaligned with the observer, those whose jet axis forms an angle θ_s to the line of sight. In Equation (A.47) the distance over which we are integrating runs along the line of sight and it is marked by the vector l in the figure. We can use simple vector algebra to relate it to the distance from the reprocessing material, x, and the distance of the blob from the BH centre, r. Assuming that the line of sight lies in the x − z plane, we can write

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathbf{x}& =\mathbf{r}+\mathbf{l}-{\mathbf{r}}_{\mathrm{re}}\hfill \\ \hfill \mathbf{r}& =(0,0,r)\hfill \\ \hfill \mathbf{l}& =(lsin{\theta }_s,0,lcos{\theta }_s)\hfill \\ \hfill {\mathbf{r}}_{\mathrm{re}}& =(r_{\mathrm{re}}sin{\theta }_{\mathrm{re}}cos{\varphi }_{\mathrm{re}},r_{\mathrm{re}}sin{\theta }_{\mathrm{re}}sin{\varphi }_{\mathrm{re}},r_{\mathrm{re}}cos{\theta }_{\mathrm{re}}),\hfill \\ \hfill \end{array}\end{array}$$(A.49)

Fig. A.3. Geometry used for the absorption calculation in case of arbitrary viewing angle θs.

where we have used the same spherical coordinates of Figure A.2 (r_re, θ_re, ϕ_re) to characterise the reprocessing material. The distance between a point on the reprocessing surface and a point on the line of sight l is

$$\begin{array}{c}\hfill \mathbf{x}=\left(\begin{array}{c}lsin{\theta }_s-r_{\mathrm{re}}sin{\theta }_{\mathrm{re}}cos{\varphi }_{\mathrm{re}}\\ -r_{\mathrm{re}}sin{\theta }_{\mathrm{re}}sin{\varphi }_{\mathrm{re}}\\ r+lcos{\theta }_s-r_{\mathrm{re}}cos{\theta }_{\mathrm{re}}\end{array}\right)\end{array}$$(A.50)

and its magnitude

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill x^2& =l^2+r_{\mathrm{re}}^2+r^2-2l{r}_{\mathrm{re}}sin{\theta }_{\mathrm{re}}sin{\theta }_{\mathrm{s}}cos{\varphi }_{\mathrm{re}}\hfill \\ \hfill & -2r_{\mathrm{re}}cos{\theta }_{\mathrm{re}}(r+l{\mu }_{\mathrm{s}})+2rl{\mu }_{\mathrm{s}}.\hfill \end{array}\end{array}$$(A.51)

Assuming μ_s = 1, at the blob position, l = 0, Equation (A.51) reduces to the distance calculation employed for the EC in case of reprocessing material, Equation (A.20).

Considering a point-like source behind the jet and inserting its specific spectral energy density, Equation (A.10), in Equation (A.47), we obtain

$$\begin{array}{c}\hfill {\tau }_{\gamma \gamma }=\frac{L_0}{4\pi {ϵ}_0{m}_{\mathrm{e}}{c}^2}{\displaystyle {\int }_0^{\infty }\mathrm{d}l\frac{(1-cos\psi ){\sigma }_{\gamma \gamma }(s)}{x^2},}\end{array}$$(A.52)

where the cosine of the angle between the photons is μ = (r + lμ_s)/x.

Considering a spherical shell BLR and replacing its specific spectral energy density, Equation (A.21), in the general opacity formula, Equation (A.47), we obtain

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\tau }_{\gamma \gamma }({\widehat{\nu }}_1)& =\frac{1}{4{\pi }^2}\frac{{\xi }_{\mathrm{line}}{L}_{\mathrm{disc}}}{{ϵ}_{\mathrm{line}}{m}_{\mathrm{e}}{c}^3}{\displaystyle {\int }_0^{\infty }\mathrm{d}l{\int }_0^{2\pi }\mathrm{d}{\varphi }_{\mathrm{re}}}\hfill \\ \hfill & {\displaystyle {\int }_{-1}^1\frac{\mathrm{d}{\mu }_{\mathrm{re}}}{x^2}(1-cos\psi ){\sigma }_{\gamma \gamma }(s).}\hfill \end{array}\end{array}$$(A.53)

In this case, we use for the distance x Equation (A.51) with r_re = R_line. The angles indicating the photon direction are μ = x_z/x and $\varphi =arctan\left(\frac{x_y}{x_x}\right)$. In the case μ_s = 1, one can see from Figure A.3 that l lies on the z axis and the integration along the line of sight can be carried out along the z axis, from r to ∞ (Finke 2016, Sec. 7). The integration over ϕ becomes trivial.

For the case of the DT, replacing its specific spectral energy density, Equation (A.27), in the general opacity formula, Equation (A.47), we obtain

$$\begin{array}{c}\hfill {\tau }_{\gamma \gamma }({\widehat{\nu }}_1)=\frac{1}{8{\pi }^2}\frac{{\xi }_{\mathrm{DT}}{L}_{\mathrm{disc}}}{{ϵ}_{\mathrm{DT}}{m}_{\mathrm{e}}{c}^3}{\displaystyle {\int }_0^{\infty }\mathrm{d}l{\int }_0^{2\pi }\mathrm{d}\varphi \frac{(1-cos\psi )}{x^2}{\sigma }_{\gamma \gamma }(s).}\end{array}$$(A.54)

The distance x is obtained from Equation (A.51) with r_re = R_DT and sin θ_re = 0. The angles indicating the photon direction can be obtained as remarked for the BLR. Considering the case μ_s = 1, the same simplifications described in the BLR case can be applied. Photoabsorption by fields external to the emission region results in an attenuation of the photon flux by a factor exp(−τ_γγ).

A.5.2. Absorption on synchrotron radiation

agnpy includes also a simplified calculations of optical depth of high energy photons on the synchrotron radiation field using an approach similar to the one in Finke et al. (2008), that is without taking into account the inhomogeneity and anisotropy of the radiation field. Contrary to Finke et al. (2008), however, we assume that the total distance over which the absorption can occur is equal to the diameter (rather than radius) of the blob, and that the two photons collide at the angle of π/2. The optical depth is then computed (in the blob frame) as:

$$\begin{array}{c}\hfill {\tau }_{\mathrm{synch}}({ϵ}_1^{\prime })=2R_b^{\prime }{\displaystyle {\int }_{{ϵ}_{min}^{\prime }}^{{ϵ}_{max}^{\prime }}\mathrm{d}ϵ\frac{{\underline{u}}_{\mathrm{synch}}^{\prime }}{{ϵ}^{\prime }{m}_{\mathrm{e}}{c}^2}{\sigma }_{\gamma \gamma }(s),}\end{array}$$(A.55)

where $s={ϵ}_1^{\prime }{ϵ}^{\prime }/2$ and ${\underline{u}}_{\mathrm{synch}}^{\prime }$ is the differential energy density of the synchrotron photons averaged over the blob volume, following Equation (A.34). It should be noticed that both the gamma-ray emission in the blob and the soft radiation on which it is absorbed fill homogeneously the blob. Therefore the attenuation factor, instead of exp(−τ), which is the case for the attenuation of the blob emission in large-scale external radiation fields, is computed as (see e.g. Finke et al. 2008)

$$\begin{array}{c}\hfill \frac{1-exp(-{\tau }_{\mathrm{synch}})}{{\tau }_{\mathrm{synch}}}.\end{array}$$(A.56)

Appendix B: Physical parameters for the examples and validation

The physical parameters used in all the examples and validations provided in the paper are listed in Tables B.1 and B.2. Table B.1 provides the parameters of the emission regions and their EEDs, Table B.2 provides instead those of the line and thermal emitters.

Appendix C: Effect of the resolution of the discretisation on the smoothness of the spectra

In this section we briefly illustrate the effect of the discretisation of the integral on the smoothness of the final spectrum. We choose as a characteristic case the SED for EC on DT, as it contains a double integral with an energetic and a spatial dimensions represented by the EED energy and the azimuth angle, respectively (see Equation (A.46)). Figure C.1 shows the results of computing the SED with a fixed dimension of frequency values (n_ν = 100) and changing in one case (left panel) the number of points along the EED energy dimension, n_γ logarithmically spaced values, and in the other (right panel) the number of points along the azimuth dimension, n_ϕ linearly spaced values. The bottom panels represent the deviations of the individual SEDs from the one with the denser integration grid (drawn in black). For this specific case it is evident that one can obtain a smoother spectrum by increasing the resolution along the EED energy dimension, while increasing the resolution along the azimuth dimension has no significant effect. We remark that the variation of the number of points along the spatial dimensions can have a different and more relevant impact if different parameters or a completely different physical process are considered. The ultimate resolution one can achieve depends on the machine used and its local memory, as larger integration grids will generate larger NumPy arrays in memory. Of course results might improve with the integration technique; we find the simple trapezoidal rule applied to logarithmically or linearly spaced integration grids acceptable for the processes described in this paper.

Fig. C.1. Effect of the resolution of the integration grid on the final spectrum. Left: variation of the number of points along the EED energy dimension, nγ. Right: variation of the number of points along the azimuth angle dimension, nϕ. In both bottom panels the deviation of the individual SEDs from the one with the denser integration grid (drawn with a dashed black line in the upper panel) is plotted.

All Tables

All Figures

Fig. 1. Scheme representing agnpy’s modules, their relations and the computations they perform. agnpy.emission_regions, defining the region responsible for the non-thermal radiation (containing the EED), is the initialising argument of the classes modelling radiative processes and computing their SEDs (agnpy.synchrotron and agnpy.compton). Thermal and line emitters, collected in agnpy.targets, are initialising arguments of the modules for the external Compton and the γγ absorption calculations (agnpy.absorption). The targets’ thermal SEDs and energy densities (u) can also be computed. The emission region is passed to the class for self-consistent modelling (agnpy.constraints) and to the absorption class to compute the γγ absorption on synchrotron photons.

In the text

	Fig. 2. SED due to synchrotron radiation computed with agnpy. Default parameters for the emission region are given in Table B.1 (SSC model column); fluxes and frequencies are given in the observer reference frame. (a) Example script using agnpy to compute and display a synchrotron SED. (b) Output of the example script.
In the text

Fig. 3. Energy densities of several photon fields: accretion disc (red solid line), BLR (blue solid line), DT (yellow solid line), CMB (black solid line), as a function of the distance from the central BH, computed with agnpy. The energy density of synchrotron radiation (green dotted line) and of the magnetic field (green dashed line) are plotted for comparison. The emitters parameters are specified in Table B.2. A BLR emitting the Lyα line is considered. The energy densities are expressed in the reference frame comoving with the second blob in Table B.1 (EC model).

In the text

Fig. 4. Total γγ opacity (black solid line and grey shaded area) produced by photon field external to the blob: BLR (red and blue dashed lines) and DT (yellow dashed line). The BLR is assumed to be composed of two concentric shells one emitting the Lyα (red) and the other the Hα line (blue). Line and thermal emitters parameters are specified in Table B.2. The green line and shaded area represent the opacity due to the synchrotron photons produced by the second blob in Table B.1 (EC model). The latter is multiplied by 103 to be comparable with the opacity due to the external photon fields. The blob is placed well within the BLR, at r = 1.1 × 1016 cm = 0.1 × RLyα. Frequencies are given in the observer reference frame.

In the text

	Fig. 5. MWL SED of Mrk 421 (black points) from Abdo et al. (2011) fitted wrapping agnpy with Gammapy. The red solid line represents the best fit model, the coloured dashed lines its individual components: yellow for synchrotron, blue for SSC. Fluxes and frequencies are given in the observer reference frame.
In the text

Fig. 6. MWL SED of PKS 1510-089 (black points) from Ahnen et al. (2017) fitted wrapping agnpy with Gammapy. The red solid line represents the best fit model, coloured dashed lines represent its non-thermal components: yellow for synchrotron, blue for SSC, green for EC on DT; grey lines represent its thermal components: dot-dashed for the multi-temperature BB of the accretion disc, dotted for the BB of the DT. Fluxes and frequencies are given in the observer reference frame.

In the text

	Fig. 7. Comparison of the synchrotron (left panels) and SSC (right panels) radiation computed with agnpy (red solid line) and jetset (blue dashed line) and taken from Fig. 7.4 of Dermer & Menon (2009) (black dashed line). The bottom panels display the deviation of agnpy from the bibliographic reference or jetset. Fluxes and frequencies are given in the observer reference frame.
In the text

Fig. 8. Left: comparison, for a distance smaller than the emitter size, of the EC scattering on the disc photon field computed with agnpy (red solid line) and jetset (dashed blue line) against the results in Finke (2016) (black dashed line). Right: comparison, for a distance much larger than the emitter size, of the SED obtained with agnpy considering the full disc energy density (solid red line) or approximating it with a monochromatic point-source (black line) with the same photon energy at the inner (dashed line) or outer (dotted line) disc radius. The SED computed by jetset at the same distance is displayed with the dashed blue line. The bottom panels display the deviation of agnpy from the reference (or the approximation). Fluxes and frequencies are given in the observer reference frame.

In the text

	Fig. 9. Same as in Fig. 8, but considering EC scattering on the BLR photon field.
In the text

	Fig. 10. Same as in Figs. 8 and 9, but considering EC scattering on the DT photon field.
In the text