A stratified jet model for AGN emission in the two-flow paradigm

High-energy emission of extragalactic objects is known to take place in relativistic jets, but the nature, the location, and the emission processes of the emitting particles are still unknown. One of the models proposed to explain the formation of relativistic ejections and their associated non-thermal emission is the two-flow model, where the jets are supposed to be composed of two different flows, a mildly relativistic baryonic jet surrounding a fast, relativistically moving electron-positron plasma. Here we present the simulation of the emission of such a structure taking into account the main sources of photons that are present in active galactic nuclei (AGNs). We reproduce the broadband spectra of radio-loud AGNs with a detailed model of emission taking into account synchrotron and inverse-Compton emission by a relativistically moving beam of electron-positron, heated by a surrounding turbulent baryonic jet. We compute the density and energy distribution of a relativistic pair plasma all along a jet, taking into account the synchrotron and inverse-Compton process on the various photon sources present in the core of the AGN, as well as the pair creation and annihilation processes. We use semi-analytical approximations to quickly compute the inverse-Compton process on a thermal photon distribution with any anisotropic angular distribution. The anisotropy of the photon field is also responsible for the bulk acceleration of the pair plasma through the"Compton rocket"effect, thus imposing the plasma velocity along the jet. As an example, the simulated emerging spectrum is compared to the broadband emission of 3C273. In the case of 3C273, we obtain an excellent fit of the average broadband energy distribution by assuming physical parameters compatible with known estimates.


Introduction
Radio-loud active galactic nuclei (AGNs) are known to be powerful emitters of non-thermal radiation. In the simplest leptonic models, this emission is most commonly attributed to the presence of highly relativistic leptons accelerated in a relativistically moving magnetized jet, emitting synchrotron radiation and inverse-Compton photons on various sources of soft photons. However several features are still unclear, such as the composition of the jet (leptonic e + /e − or baryonic p + /e − ), the bulk acceleration mechanism, the heating mechanism of the relativistic non-thermal particles, and the precise size and location of the emitting zones. One-zone models are the most simple and widespread models for reproducing the AGN jet emission. They assume a spherical, homogeneous emission zone with a minimal number of free parameters: the radius of the zone, the magnetic field, the bulk Lorentz factor, and parameters describing the particle distribution. They present the advantage of being simple and give a good first approximation of the physical conditions in jets. However, they encounter several limitations. Among others, they have difficulties in reproducing the low-energy (radio) part of the spectral energy distribution (SED; most probably emitted by the farthest part of the jet). They assume the very stringent condition that the whole non-thermal emission must be produced in a single zone, which can be an issue for explaining the multi-wavelength variability of the sources. Furthermore, they do not offer any clue on the jet formation mechanism and its parameters outside this zone. Also the strong Doppler boosting associated with highly relativistic jets is incompatible with the detection of high-energy emission from unbeamed radio galaxies seen at large angles, since their emission should be strongly attenuated by a Doppler factor smaller than one.
Facing these weaknesses, models considering more complex structures have been proposed involving stratified inhomogeneous jets. For instance, the blob-in-jet models (Katarzynski et al. 2001;Hervet et al. 2015) propose a structure where blobs move at high relativistic speed in the jet. Those blobs can be responsible for some of the emission (especially at high energies) whereas the ambient jet can explain the rest of the spectrum, for example in the radio band. Spine/sheath models, implying the existence of two flows at different velocities, can provide a picture in agreement with both theoretical and observational constraints. Ghisellini et al. (2005) developed a model based on the same idea and showed that this kind of structure could help reduce the necessity for very high bulk Lorentz factor values as the two emitting zones interact radiatively, enhancing each other's emission.
But the original idea of a double jet structure stemmed from Sol et al. (1989), who coined the name "two-flow model" (more details in Sect. 2). In their original paper, the authors proposed for the first time a double jet structure for AGN jets. In this model, an outer jet or collimated wind is ejected from an accretion disc, with a mildly relativistic velocity (υ ∼ 0.5c). In the empty funnel of this jet, a fast inner electron-positron beam is formed and moves at much higher relativistic speeds (bulk Lorentz factor Γ b ≈ 10). The pairs are supposed to be continuously heated by the surrounding baryonic jet through MHD turbulence.
The model has several advantages compared to models assuming a single fluid. The first advantage is that the problem of the energy budget of relativistic jets is reduced, since only a minor component of leptons needs to be accelerated to high bulk Lorentz factors. The protons of the surrounding jets are not supposed to be highly relativistic. The second is that this model can provide a simple way to solve the discrepancy between the required high Lorentz factors to produce the observed gamma-ray emission and the slower observed motion in jets at large scales (e.g. Henri & Saugé 2006 and references therein). Furthermore, as the power is carried out mainly by the non-relativistic jet, it escapes the Compton drag issue. As explained below, the pair beam is only gradually accelerating thanks to the anisotropic inverse-Compton effect (or "Compton rocket" effect), and its velocity never exceeds the characteristic velocity above which the aberrated photon field starts to cause a drag (it actually remains at this characteristic velocity). Its density increases all along the jet due to the gamma-gamma pair production process, so the dominant emission region can be at large distances from the central black hole, avoiding the problem of gamma-gamma absorption by the accretion disk photon field. This model therefore offers a natural explanation of the main characteristics of the high-energy source deduced from observations. Noticeably, the dynamical effects of radiation on relativistic particles (which are intrinsically strongly dissipative) are very difficult to incorporate both in analytical and in numerical (M)HD simulations. The picture we present here is thus markedly different from most models available in the literature, since the pair component dynamics is mainly governed (at least at distances relevant for high-energy emission) by these radiative effects. However we argue that these effects are unavoidable since the cooling time of relativistic leptons is indeed very short at these distances, and these effects must be taken into account in any physically relevant model implying a relativistic pair plasma, which is in turn likely to exist given the high density of gamma-ray photons observed from radio-loud gamma-ray-emitting AGNs.
The first numerical model of non-thermal emission based on these ideas was proposed by Marcowith et al. (1995),who considered only the inverse-Compton process on accretion-disk photons. Assuming a power-law particle distribution and a stratified jet, the authors could derive the inverse-Compton emission from a plasma of relativistic leptons illuminated by a standard accretion disc as well as the opacity to pair production and the pair production rate . They showed that the spontaneous generation of a dense e + /e − pair beam continuously heated by the baryonic jet was indeed possible and was able to reproduce the gamma-ray emission of EGRET blazars.
Based on this work, Marcowith et al. (1996) further studied the possibility of accelerating the e + /e − pair beam in the framework of the two-flow through the Compton rocket mechanism introduced by O'dell (1981). In this mechanism, the motion of the relativistic pair plasma is entirely controlled by the local anisotropy of the soft photon field, which produces an anisotropic inverse-Compton emission transferring momentum to the plasma. The acceleration saturates when the photon field, aberrated by the relativistic motion, appears nearly isotropic (vanishing flux) in the comoving frame. Renaud & Henri (1998) continued this work on the acceleration via the Compton rocket effect. Considering a relativistic pair plasma following a powerlaw energy distribution coupled with the photon field from a standard accretion disc and extended sources (BLR or dusty torus); they computed the value of the terminal bulk Lorentz factor and showed that values up to Γ b = 20 are achievable for extragalactic jets, the most probable values being of the order of 10 in good agreement with VLBI motions studies (e.g. Lister et al. 2009). Subsequent works studied the possibility of explaining the spectra by a pile-up (relativistic Maxwellian) distribution (Saugé & Henri 2004), which better reproduces the high-energy spectra of BL Lacs that peak in the TeV range. A quasi-Maxwellian distribution is close to a monoenergetic one, however the spatial convolution of such a distribution whose parameters vary along the jet can mimic a power-law over a limited range of frequencies. These authors also proposed a time-dependent version of the model that was later shown to be able to successfully reproduce the rapid flares observed in some objects, such as PKS 2155-304 (Boutelier et al. 2008).
In a recent study, Vuillaume et al. (2015) (hereafter Vu15) further studied the acceleration through the Compton rocket effect. The complex photon field of an AGN was considered, carefully taking into account the spatial distribution of extended sources (standard accretion disc, dusty torus, and broad line region). The evolution of the resulting bulk Lorentz factor along the jet is then computed self-consistently and it appears that due to the complexity of the surrounding photon field, it can display a relatively tangled relation with the distance to the base of the jet. Moreover, variations of the bulk Lorentz factor (and thus of the Doppler factor) can induce complex variability of the observed emission in space and in time.
The goal of this paper is to present the complete calculation of the pair-beam non-thermal spectra assuming the above configuration, that is, 1) the pair plasma is assumed to be described by a pile-up distribution and is generated in situ by γ − γ interactions, 2) its motion is controlled at short and intermediate distances (up to ∼ 10 3 r g ) by the anisotropic Compton rocket effect in the complex photon field generated by an accretion disk, a broad line region, and a dusty torus, and 3) particles emit non-thermal radiation through synchrotron, synchrotron self-Compton (SSC), and external Compton (EC) in the various photon fields. Part of the high-energy γ−ray photons can also be absorbed to produce new pairs. We assume that these new pairs are continuously accelerated along the jet. In the two-flow paradigm (Sol et al. 1989), such acceleration is expected by the MHD turbulence generated by the surrounding baryonic flow. Subsequently, the density and internal energy of the pair plasma are computed self-consistently, and the emitted photon spectrum can be evaluated and compared to observations. The overall layout is as follows. In Section 2 we present the main theoretical interest of the two-flow paradigm. Then in Section 3 we detail our model and the numerical methods we use. In section 4 we apply this model to the bright quasar 3C273 Article number, page 2 of 17 and show the model can reproduce its jet emission from radio to gamma-rays.

The two-flow paradigm: the hypothesis and the reasoning behind it
The two-flow paradigm is based on an original idea from Sol et al. (1989) (see Section 1).The model has evolved since then but the core hypothesis remains the same: an AGN "jet" is actually made of two interacting flows. The outer one is a MHD jet, or wind, fuelled by the accretion disc. It originates from and is self-collimated by the Blandford & Payne (BP) process (Blandford & Payne 1982) and is much like the jets found in other objects such as young stars or neutron stars. This baryonic jet is therefore mass loaded by the disc and mildly relativistic (β ≈ 0.5). On the rotation axis, where the angular momentum tends to zero, no MHD force counteracts the gravity from the central object and it is expected that the density strongly decreases near the axis, thus leaving an empty funnel. A lighter jet made of relativistic leptons (e − /e + ) can then be created there through γ − γ interaction between the surrounding soft and highenergy photons. This leptonic plasma is accelerated through the Compton rocket effect (as explained below) and travels at highly relativistic speeds. It is assumed to be confined, collimated, and continuously reheated by the surrounding MHD jet.

Interaction of highly relativistic flows
The first benefit of the two-flow hypothesis is to alleviate the problem of the confinement of a highly relativistic flow. Self-confinement of a jet can take place due to the magnetic field exerting a magnetic pressure (from the Lorentz force) balancing internal pressure and centrifugal force. However, the self-confinement of a highly relativistic jet through this process is quite inefficient. This has been shown first in numerical simulations by Bogovalov 2001 andBogovalov &Tsinganos 2001 and then demonstrated based on theoretical arguments by Pelletier 2004. Therefore the collimation of relativistic flows necessarily requires an external process like for example external pressure from the ambient (interstellar) medium.. In the two-flow paradigm, the outer self-confined massive and powerful MHD sheath confines the spine by ram pressure, providing an easy solution to the important problem of confinement of highly relativistic flows.
Moreover, the interface of the two flows can be subject to Kelvin-Helmholtz instabilities producing turbulence. This turbulence can then accelerate particles in the spine through secondorder Fermi processes. In that picture, because the MHD sheath carries most of the power, it can be seen as an energy reservoir continuously energizing the particles through turbulence. This continuous source of energy gives rise to two very interesting phenomena.
The first one is the most important feature of the pair creation process. As new pairs are created through γ − γ absorption, and immediately accelerated through turbulence to reach high energies, they can emit γ-rays that will create more pairs. The pair-creation process being very efficient and highly non-linear, a copious amount of new pairs can be created even from an initial, very low density, therefore constituting the spine. Moreover, above a certain threshold of energy, the process can runaway and give rise to episodes of rapid flares. This has been demonstrated by Renaud (1999), Saugé & Henri (2004), and Boutelier et al. (2008).
The second phenomenon is the possibility to accelerate the spine jet to relativistic motion through the anisotropic Compton rocket effect as discussed below.

Bulk Lorentz factor of the spine
The question of the actual speed, or bulk Lorentz factor Γ b of jets, as well as their acceleration, is central to the understanding of their physics. This is a long-standing and debated issue in the community (Henri & Saugé 2006) with conflicts between theoretical arguments (Aharonian et al. 2007, Tavecchio et al. 2010, Begelman et al. 2008) and observations (Piner & Edwards 2004;Giroletti et al. 2004;Piner & Edwards 2014;Lister et al. 2013). Some of the attempts to solve this issue come in the form of structured jets (Chiaberge et al. 2000;Ghisellini et al. 2005;Georganopoulos & Kazanas 2003).
In the two-flow paradigm, the question of the acceleration of the spine is solved by the Compton rocket effect as proposed by O'dell (1981). In this process, inverse-Compton scattering in a strongly anisotropic photon field induces a motion of the emitting plasma in the opposite direction. O'dell (1981) showed already that a purely leptonic hot plasma must be dynamically driven by this process, up to relativistic speeds. Phinney (1982) opposed the fact that this process is actually quite inefficient at accelerating a pair blob as the pairs cool down very quickly through the inverse-Compton scattering, therefore killing it before it can be effective.
However, in the two-flow paradigm, as the pairs are continuously re-accelerated by the turbulence all along the jet, this argument does not hold and the Compton rocket process becomes an efficient source of plasma thrust (see Section 3.6). This effect is likely to be very efficient in type II AGNs (FSRQ and FR II galaxies) where an accretion disk with high luminosity is present. The ambient photon field in the first hundreds of r g will therefore be highly anisotropic and a pair plasma will be rapidly accelerated to a characteristic velocity for which the net aberrated flux, evaluated in the comoving frame, vanishes. As the cooling timescale for near-Eddington luminosities is much shorter than the other dynamical timescales, this effect will dominate over all other terms such as pressure gradients and magnetic effects. Hence the velocity of the plasma will be very close to this characteristic equilibrium velocity. If the photon field is mainly external (accretion disc and secondary reemission processes), and the scattering occurs in the Thomson regime, the characteristic velocity depends only on the angular distribution of photons and not on the plasma particle distribution, as explained in Section 3.6. We note however that the situation may be different for type I objects (FR I and BL lacs) where the radiation field is dominated by the jet itself; in this case , there is no simple calculation of the equilibrium velocity. In the following, we only consider FSRQ with an intense external photon field. This does not mean that the two-flow model is invalid for type I objects, but only that the velocity of the pair spine is not easily computed for these objects. In the same way, the application to other kinds of objects such as microquasars and GRBs is probably different given the very different physical conditions of the photon source.

Numerical modelling
The numerical model developed here is based on the one described in Boutelier et al. (2008). The jet is stratified along the axis coordinate z (see Fig. 1). The jet axis is viewed at an angle θ i with the line of sight, with µ i = cos θ i . Numerically, we define a slicing with an adaptive grid step as described below and physical conditions are computed at each slice, starting from initial conditions given at the basis of the jet. Each slice then acts as a one-zone and emits synchrotron radiation, synchrotron self-Compton (SSC), and external Compton (EC) radiation, computed as in a one-zone approximation with a spherical blob of radius R(z). The photon field from external sources (the accretion disc, the dusty torus, and the broad line region) is computed all along the jet. This determines the external inverse-Compton emission as well as the induced Compton rocket force on the plasma as described below. The opacity to high-energy photons inside and outside the jet is computed numerically, and is used to compute a pair creation term at each step.

Geometry of the jet
To compute the total jet emission, one needs to know the physical conditions all along the jet. Some of these conditions are derived from the computation of the spatial evolution equation along z (see Sect. 3.5) and only three physical parameters, that is, the inner radius of the jet R(z) (acting as an outer radius for the pair beam), the magnetic field B(z), and the heating term Q(z), are imposed on a global scale. Here we assume their evolution to be described by power-laws: This law describes a paraboloid with a radius close to R 0 at a distance Z 0 from the black-hole. The constant (R i /R 0 ) 1/ω allows to avoid divergence issues at z = 0 by setting a minimal jet radius R(z = 0) = R i corresponding to the disc inner radius. The starting and ending altitudes of the jet are free parameters. The index ω defines the jet opening. One must have ω < 1 to keep the jet collimated.
The magnetic field is supposed to be homogeneous and isotropic at every altitude z in the plasma rest frame. Its evolution is described by: The index λ gives the structure of the magnetic field and is confined by two extremes: λ = 1 corresponds to a pure toroidal magnetic field as the conservation of the magnetic field circulation gives B ∼ 1/R. -λ = 2 corresponds to a pure poloidal magnetic field as the conservation of the magnetic flux in this case gives B ∼ 1/R 2 .
The particle acceleration is a central part of the two-flow model as it is assumed that the particles are continuously heated by the outer MHD structure (see Sect. 2) which acts as an energy reservoir, compensating for radiation losses. Due to the lack of a precise expression for the acceleration rate per particle, Q acc , we use the following expression: The particle acceleration decreases as a power-law of index ζ up to an altitude z = Z c where it drops exponentially. This altitude Z c physically corresponds to the end of the particle acceleration (through turbulence) in the jet. Because of this exponential cut-off, whatever the index ζ is, the total amount of energy provided to accelerated particles remains finite. However, as Z c could be as large as desired (even as large as the jet), we consider ζ > 1 to be physically more satisfactory. This way, even an integration to infinity of Q acc would converge. Similarly to the jet radius expression (equation 1), the constant R i /R 0 avoids numerical issues for very small z.
The jet is then sliced along z. As we see below, the particle density can be subject to abrupt changes in case of short and intense events of pair creation which is a non-linear process. It is therefore essential to have an adaptive slicing as the physical conditions are computed in the jet. The condition have been chosen to ensure variation rates of the particle mean energy and of the particle flux of less than 1 between two computation steps.
Frequencies follow a logarithmic sampling between ν min and ν max in the observer frame. The sampling is therefore different at each altitude in the jet depending on the Lorentz factor. This ensures that local emissivities are computed for the same sampling (that of the observer) when transferred to the observer frame.

Geometry of the external sources of photons
There are several possible sources of soft photons in an AGN but three have an actual influence on the external Compton emission and are taken into account in our model: the accretion disc, the dusty torus, and the broad line region (BLR). In order to correctly compute the corresponding inverse-Compton radiation, the anisotropy of the sources is taken into account as detailed below.

The accretion disc
The geometry of the disc (see Fig. 2) is described by its internal radius R in and its external radius R out . It is then sliced along its radius r (with a logarithmic discretization) and its azimuthal angle ϕ (with a linear discretization). Therefore, each slice has a surface dS = dϕ rdr + dr 2 2 . From the jet axis at an altitude z, each slice is seen under a solid angle at an angle θ s = arccos(z/ √ r 2 + z 2 ) with the z axis. We consider here a standard accretion disc. Then the radial distribution of the temperature is given by T disc (r) (Shakura & Sunyaev 1976): with G the gravitational constant, σ the Stefan-Boltzmann constant, M the black-hole mass,Ṁ the accretion rate, and R isco = 6MG c 2 the innermost stable circular orbit. Its emissivity is equal to one if not specified otherwise.
The luminosity of one face of the disc is then given by which converges to L disc =Ṁ c 2 24 for R in = R isco and R out R in .

The dusty torus
The dusty torus (see Fig. 3) is assumed to be in radiative equilibrium with the luminosity received from the accretion disc. The torus is sliced according to θ t ∈ θ t min , θ t max and ϕ ∈ [0, 2π] . Each slice is illuminated by the disk and is in radiative equilibrium. It emits as a grey-body with a temperature T t (θ t , ϕ) and an emissivity ε(θ t , ϕ). As most of the disc luminosity comes from its inner parts and considering R in (D t − R t ), the disc appears point-like from a torus surface point. Therefore, each slice of the torus of surface dS t and at an angle θ t from the X axis verifies the relation: with ω d the angle between the Z-axis and the emission direction from the disc: cos ω d = a sin θ t / 1 − 2a cos θ t + a 2 and a = R t /D t . The user can then either fix a constant temperature or a constant emissivity and compute the other parameters as a function of θ t to preserve the radiative equilibrium. If not specified otherwise, ε(θ t ) = 1 and the temperature is kept free.

The broad line region
The broad line region is modelled by an isotropic, optically and geometrically thin spherical shell of clouds situated at a distance R blr from the central black-hole. Like other sources, it is sliced angularly into different parts in order to perform the numerical integration. We choose a linear discretization along ω ∈ [ω max , ω min ] and along ϕ ∈ [0 : 2π]. Observed BLRs display a complex emission with a continuum and broad absorption lines but Tavecchio & Ghisellini (2008) showed that modelling the spectrum of the BLR with a grey-body spectrum at T = 10 5 K provides a good approximation to the resulting inverse-Compton spectrum. We followed this idea using a temperature of T blr = 10 5 K and an overall luminosity L blr , which is a fraction α blr of the disc luminosity: Article number, page 5 of 17 A&A proofs: manuscript no. vuillaume2018 Emissivity of the BLR is then given by:

Emission processes
As the spine is supposed to be filled by electrons/positrons, only leptonic processes need to be considered here: synchrotron and inverse-Compton radiation are computed all along the jet. The radiation is computed in the plasma rest frame and then boosted by the local Doppler factor δ(z) into the observer's rest frame -with β b = 1 − 1/Γ 2 b and i obs the observer angle relative to the jet axis as defined in Fig. 5.
In the two-flow paradigm, particles are supposed to be accelerated by a second-order Fermi process due to turbulence inside the jet. Thus a Fokker-Planck equation governs the evolution of the particle distribution which evolves by diffusive acceleration and radiative losses. (Schlickeiser 1984) showed that generic solutions of this equation were pile-up (or quasi-Maxwellian) distributions. Most of the particles are then concentrated around some characteristic energyγ where the acceleration and cooling time are similar. We adopt the following simplified form of such a distribution: where n 0 (z) = n e (γ)dγ.
As this acceleration mechanism does not produce power laws for the particles energy distributions (contrary to the first-order Fermi process at work in shocks), the reproduction of large-scale power laws in the spectra of the sources is not straight forward. Power-law shapes can however be reproduced by a summation of a number of pile-up (from different emission zones) with different parameters. However, the summation of the emission coming from each slice of the jet can be demanding in terms of computing time. A certain level of approximation is thus required to achieve computation of the model in a reasonable amount of time. The following subsections present the calculation of each type of radiation for a pile-up distribution.

Synchrotron radiation
Using the expression of the synchrotron-emitted power per unit of frequency and solid angle ( (Blumenthal & Gould 1970)), (Saugé 2004) showed that the synchrotron emissivity of a pile-up distribution can be written as: heBγ 2 m 2 e c 3 , y = / c and To fasten numerical computation, analytical approximations of the function Λ(y) can be done in different regimes. These approximations are presented in the appendix (equations A.1, A.2 and A.3).
Finally, one obtains the synchrotron spectrum in the optically thin regime ( > abs ): The transition between the optically thin and the optically thick regime happens at the frequency abs defined by For a pile-up distribution of particles, the distribution temperature k B T e = γ m e c 2 /3 =γm e c 2 does not depend on . In this case, the optically thick regime at low frequency ( < abs ) is described by the Rayleigh-Jeans law: Finally the synchrotron spectrum resulting from a pile-up distribution has three main parts: -< abs is the optically thick part of the spectrum described by a power-law of index 1. abs < < c is the optically thin part of the spectrum described by a power-law of index -2/3. c < is where the spectrum falls exponentially.

Synchrotron self-Compton radiation
We assume the synchrotron self-Compton (SSC) emission to be co-spatial with the synchrotron emission, and therefore neglect the synchrotron emission coming from other parts of the jet. We treat the Thomson and the Klein-Nishina (KN) regimes separately, following the analytical approximations proposed by (Saugé 2004). The results of these approximations are recalled here while the details are given in the appendix. We first consider the distinction between the Thomson and the KN regimes relative to the synchrotron peak c .
If c 1/ 1 , all synchrotron photons producing SSC photons of energy 1 are scattered in the Thomson regime. In this case, one can show that the SSC photon production rate is given by (see appendix A.2): with x dt e −t /t being the exponential integral function. Interestingly,G is a function of a single variable and can therefore be tabulated to speed up calculation in the Thomson regime.
For c > 1/ 1 , we must take into account KN corrections. However, synchrotron photons verifying < 1/ 1 are still in the Thomson regime and photons verifying > 1/ 1 are in the KN regime. In order to include KN corrections only when necessary, the emissivity in this regime is divided into two contributions, Thomson and Klein-Nishina, respectively given by J th ssc and J kn ssc : with Γ being the incomplete gamma function and with u min = 1 maxγ 2 and u max = 1 minγ 2 , The SSC photon production rate in the Klein-Nishina regime is then given by dn kn ssc ( 1 ) d 1 dt = n 1 n 0 cσ T h J th scc ( 1 ) + J kn ssc ( 1 ) .
The continuity between the two regimes given by Eqs. 15 and 19 is assured by an interpolation formula that gives the complete expression of the SSC radiation: We used some examples to verify that the value n = 6 gives a correct approximation of the full cross section. However the final results are relatively insensitive to the choice of n since the various contributions are smoothed by the spatial convolution of the different parts of the jet.

External Compton radiation
The calculation of the inverse-Compton emission on a thermal distribution of photons (a.k.a. external Compton) is the most demanding in terms of computation time since it requires an integration over the energy and spatial distributions of the incoming photons which are not produced locally. We note that the anisotropy of the incoming radiation is important and has to be properly taken into account for the computation of the emissivity and the bulk Lorentz factor of the plasma through the Compton Rocket effect (see section 3.6). Since all our external emission sources (disk, BLR and torus) can be approximated by a blackbody energy distribution(or the sum of several blackbodies), some approximations are possible and have been proposed by Dubus et al. (2010), Khangulyan et al. (2013), and Zdziarski & Pjanka (2013) (hereafter ZP13). The method we developed and used in this paper is less precise than the one proposed by ZP13 but is at least twice faster; and up to ten times faster in some cases. It approximates the Planck's law for black-bodies emission by a Wien-like spectrum. Details of the calculation and comparison with ZP13 are done in Appendix B.
The number of inverse-Compton photons per energy unit and time unit produced by the scattering of a thermal photon field on a single particle of energy γm e c 2 is then given by: with x = 1 γ , θ 0 being the angle between the incoming photon and the particle direction of motion, H = x (1 − x) , and dt being the exponential integral. This relation takes into account the anisotropy of the emission through the angle θ 0 and the full Klein-Nishina regime.
Equation (21) needs to be integrated over the pile-up distribution to get the complete spectrum. This can be long in the general case. However, when applicable (we chose the conservative limit of < 0.01), which greatly simplifies in the Thomson regime (see Appendix B) and one obtains: with with χ(s) being a single variable function. As such, it can be computed once and then tabulated and interpolated over when required. As a result, when the emission occurs in the Thomson regime, the computation of the inverse-Compton spectra from the scattering of a thermal soft photon field on a pile-up distribution of electrons can be done much faster than usual with complete numerical integration.

Photon-photon absorption in the jet and induced pair creation
High-energy photons produced by inverse Compton will locally interact with low-energy photons produced by the synchrotron process (external radiations can be locally neglected). This photon-photon interaction induces an absorption that needs to be taken into account to compute the emitted flux and the pair creation process loading the jet with leptons.

Escape probability in the jet
For a photon of dimensionless energy, 1 = hν 1 /m e c 2 in a photon field of photon density per solid angle, and per dimensionless energy, n ph ( 2 , Ω), the probability to interact with a soft photon of energy 2 on a length dl is given by: with σ(β) being the interaction cross-section given by (Gould & Schréder 1967): with µ being the cosine of the incident angle between the two interacting photons in their reference frame.
If we make the approximation that the synchrotron emission is locally isotropic in the plasma rest frame, the optical depth per unit length dl simplifies and is then given by (Coppi & Blandford 1990 with R pp the angle-averaged pair production rate (cm 3 .s −1 ). Analytical approximations of R pp have been proposed by Coppi & Blandford (1990) and Saugé & Henri (2004) but they still require a timeconsuming integration over all energies. In order to simplify numerical calculations, as the function R pp is peaked around its maximum R max pp occurring at x = x max , we make the approximation: with R max pp = 0.283 with σ T h being the Thomson cross-section.
The optical depth in the jet finally simplifies to Assuming that the absorption coefficient is constant at a given altitude z in the jet of radius R(z), the opacity can be calculated as Solving the transfer equation in the plane-parallel approximation for photons absorbed in-situ gives their escape probability P jet esc : This escape probability is useful for computing the pair production rate inside the jet. However, another absorption factor must be considered in the vicinity of the jet as discussed in Sect. 3.7, since the photon density does not vanish abruptly outside the jet.

Pair creation
The pair production is a direct consequence of the γ − γ absorption in the jet. A photon of dimensionless high-energy > 1 interacts preferentially with a soft photon of energy ≈ 1/ to form a leptonic pair e + /e − . Both created particles have the same energy γm e c 2 (with m e the electron mass) and one can write the energy conservation as + 1/ ≈ = 2γ. Therefore, at a given altitude of z, the pair production ratė n prod is given bẏ This pair creation is then taken into account in the evolution of the particle density. The created particles are not supposed to cool freely but they are rather constantly reaccelerated by the turbulence. We assume that the acceleration process is fast enough to maintain a pile-up distribution without considering the perturbation to the energy distribution introduced by the pair creation.

Evolution of the particle distribution
The pile-up distribution has two parameters, the particle density n 0 (z) and the particles characteristic energyγ(z)m e c 2 (the particles mean energy being given by 3γm e c 2 ). Both parameters are not imposed but computed consistently all along the jet as detailed below.

Evolution of the particles' characteristic energy
The particles' characteristic energyγ(z)m e c 2 results from the balance between heating from the turbulence and radiative cooling. The energy equilibrium for each particle in the comoving frame can be written as follows. ∂ ∂t γm e c 2 = 1 Γ b Q acc m e c 2 − P cool , with the acceleration parameter Q acc (s −1 ) (see Eq. 3) and P cool being the emitted power derived from synchrotron emission, SSC, and EC emission. We consider here that the cooling is efficient and occurs mainly in the Thomson regime, neglecting the cooling in the Klein-Nishina regime. In the case of isotropic distribution of photons, Rybicki & Lightman (1986) showed that the emitted power at the characteristic particles energy writes P = 4 3 cσ T U γ 2 − 1 with U being the total energy density of the photon field. One must compute the contribution of each energy density corresponding to each emission process: For the synchrotron emission, we consider an isotropic magnetic field. Rybicki & Lightman (1986) gives The power emitted through SSC is computed considering the effective energy density of the synchrotron photon field corresponding to the Thomson regime (using a cut-off frequency ν kn ). As the synchrotron emission is isotropic and co-spatial with the SSC emission, one obtains with I syn being the local synchrotron specific intensity.
The external photon field is not isotropic and one should integrate over all directions to derive the power emitted through external inverse-Compton. To ease numerical computation, we assume that the total emitted power at each altitude z can be approximated by the power emitted in the direction perpendicular to the jet axis and integrated over 4π: with dN dtd (t, ) being the emitted external Compton spectrum in the comoving frame. Finally, we compute the characteristic particle Lorentz factor by numerically solving the following energy equation in the comoving frame: This equation is solved in the stationary regime in each slice of the jet using a Runge-Kutta method of order 4.

Evolution of the particles density
The particle density evolves with the pair production and annihilation. As the particles move forward in the jet, one can apply flux conservation to compute the evolution of the density along the jet.
By generalizing the standard continuity equation, (Boutelier 2009) showed that for a stationary jet structure, the flux conservation can be written as follows.
with S (z) = πR 2 (z) being the section of the jet andṅ prod the pair production rate given by equation 32. In the stationary case, one simply solves ∂ ∂z Φ e = 1 c Sṅ prod .
As the new particles are created in the jet filled with turbulences, they are accelerated and can emit more radiation, creating more particles. This process is highly non-linear and can lead to explosive events. In the two-flow paradigm, this can lead to fast and powerful flares as the particle density will increase as long as the energy reservoir is not emptied.

Evolution of the bulk Lorentz factor of the jet
As shown by O'dell 1981, the radiative forces act on a hot plasma dynamic through the Compton rocket effect. In the present case, the timescale of the Compton rocket effect is equal to the inverse-Compton scattering time of a photon field of energy density U ph on a particle of energy γm e c 2 and thus can be evaluated as: In the inner regions of powerful AGNs, the photon field energy density at a distance z in the jet can be evaluated as with L edd = 4πGMm p c/σ T being the Eddington luminosity and R g = GM/c 2 the gravitational radius. This allows us to compare the inverse-Compton scattering time with the dynamical time of the system t dyn = z/c as The bulk Lorentz factor of the plasma is the result of a balance of the radiative and dynamical forces. However, as shown here, in the inner parts of an AGN, the inverse-Compton scattering time is several orders lower than the dynamical time. In this case, a purely leptonic plasma is strongly tide to the photon field and its dynamic must be imposed by the inverse-Compton scattering as others forces (such as MHD forces or the interaction with the external flow) acting on timescales of the order of the dynamical time are too slow to counteract this force.
Therefore, the bulk Lorentz factor of the jet Γ b (z) is not a free parameter in the model but is imposed by the Compton rocket force and thus by the external photon fields.
As shown in Vuillaume et al. (2015) (hereafter Vu15), as long as most of the emission happens in the Thomson regime, the computation of the resulting Γ b (z) depends exclusively on the distribution of the external photon field in the inner parts of the jet whereas its final value Γ ∞ reached at parsec scales when the bulk acceleration stops or becomes ineffective depends on the jet energetics. Indeed, acceleration through the Compton rocket effect ceases when the scattering time in the rest frame of the plasma becomes greater than the dynamical time.
We use the method described in Vu15 to determine the equilibrium bulk Lorentz factor Γ eq (z) imposed by the Compton rocket and and only take into account the geometry of the external photon fields in the Thomson regime. We make the assumption that most of the inverse-Compton scattering always happens in the Thomson regime and we can verify this assumption afterwards for each object modelling.
From there, the evolution of Γ b (z) is determined by solving the differential equation (Vu15) with l(z,γ) being the relaxation length to equilibrium that can be written as with H = 1 4π I ν s (Ω s )µ s dΩ s dν s being the Eddington parameter proportional to the net energy flux of the external photon fields on the jet axis. In order to ease numerical computation, we do not solve the complete differential equation. The relaxation length is determined at each step and as long as l(z,γ) z (equivalent to t ic t dyn ), one can consider that Γ b (z) = Γ eq (z). At high z in the jet, or whenγ decreases enough, one has l(z,γ) z. In this case, jet moves in a ballistic motion and its speed reaches a constant Γ ∞ . Numerical tests showed that Γ eq (z) reaches Γ ∞ when l(z)/z ≈ 0.6. At this point we fix Γ b (z) = Γ ∞ .
The bulk Lorentz factor then reaches an asymptotic value corresponding to ballistic motion of the plasma. In the present model, we assume that the jet follows such a ballistic motion and is not affected by external factors (such as interstellar medium or two-flow interactions).

Absorption outside the jet
Between the jet and the observer, a photon encounters several radiation fields causing absorption through photon-photon interaction. There are two main sources of external radiation that should be considered here: the immediate vicinity of the jet and the extragalactic background light (EBL) (the absorption from the host galaxy being negligible).

Extragalactic background light
Extragalactic background light (EBL) is composed of several sources: the CMB, the diffuse UV-optical background produced by the integrated emission from luminous AGNs and stellar nucleosynthesis, and the diffuse infrared background composed of the emission from stars' photospheres and thermal emission from dust. Here we use absorption tables from (Franceschini et al. 2008). This introduces another term, exp(τ ebl ), to the escape probability that depends on the source redshift z s .

Absorption in the jet vicinity
The immediate vicinity of the jet is composed of several sources of soft photons. Accretion disc BLR Dusty torus log(Z 0 /R S ) = 1 2 3 4 5 6 7 8 Fig. 6. Integrated τ ext created by the three external sources (disc, torus, BLR) and experienced by a photon of frequency ν leaving the jet from an altitude z 0 and travelling along the jet axis to infinity. Parameters of the sources: Disc: L d = 0.2L edd , R in = 3R S and R out = 5e4R S . Torus: R t = 5e4R S and D t = 1e5R S . BLR: R blr = 8e3R S , ω ∈ [0 : π/2], α blr = 0.01.
First, we consider the synchrotron photons from the jet itself. Marcowith et al. (1995) showed that these photons interact on a typical distance corresponding to the jet radius R, which introduces an absorption term exp(−τ jet ) with τ jet given by equation 30.
Once a photon has escaped the jet photon field, it enters the external photon field generated by external sources described in Sect. 3.2. The induced opacity is generally greater than one in FSRQ and therefore cannot be neglected.
As the external photon field resulting from the sources is highly anisotropic, one cannot use the approximation used in Sect. 3.4.1. To compute the opacity, τ ext , experienced by a photon emitted from the jet at an altitude z 0 , one needs to integrate the complete γγ absorption given by equation 24 over the path followed by the photons from z 0 to infinity (here M relates to the position along that path): Terms in equation 44 have been detailed at equation 24. The factor θ(M) is the angle between the direction of the photon from the jet and the photons from thermal sources as illustrated in Fig.  5.
As an example, the result of this numerical integration, carried out for a range of energies and photon-emitted altitudes, is given in Fig. 6. The figure represents isolines of photon-emitted altitude z 0 . In order to validate our approach, we chose an observing angle i obs = 0 and the same source parameters as in Reimer (2007) and find consistent results.
Considering the different sources of absorption described above, the total escape probability of photons produced in the jet is then given by and the specific intensity reaching the observer is thus given by I obs ( ) = I jet ( )P tot esc ( , z s ).

Example of application: the quasar 3C 273
3C 273 is the first quasar to have been identified thanks to its relatively close distance (redshift z = 0.158, Schmidt 1963) and has been extensively observed and studied. The broadband SED data used here come from Turler et al. (1998); they averaged over 30 years of observations and are more likely to represent the average state of the AGN. This quasar represents a good test for the model as this average state can be associated to a quiescent state that we can model. Moreover, it is a FSRQ with relatively wellknown external sources, allowing us to test the complete model with good constraints on these sources while imposing values of Γ b (z) through the Compton rocket effect (see section 3.6). We detail our modelling of this source in the following section.

Modelling
The accretion disc is easily visible in the SED in Fig. 7 as the big 'blue bump' in the optical band. For our modelling, me make a best-fit of this bump (between 2.4 × 10 14 Hz and 2.4 × 10 15 Hz) considering the temperature given by a standard accretion disc (see equation 5) and with the two free parameters: R S = 1 3 R isco and L disc , the total luminosity of the disc. We obtain R S = 5.1 × 10 14 cm and L disc = 1.7 × 10 46 erg.s −1 (χ 2 = 0.2). The modelling of the disc is given in orange in Fig. 7.
The signature of the dusty torus is perceptible as a bump in the infrared band. Parameters of the torus are chosen to be in agreement with the position and luminosity of this bump. The resulting size of the torus is R T = 10 4 R S and D T = 1.5 × 10 4 R S .
Strong emission lines are observed, indicating the presence of a BLR. To model the BLR, we chose its radius from the value of Paltani & Turler (2005) that derived a size of R blr /c = 986 days by studying the lag of the Balmer lines. The BLR opening is set to an ad hoc value ω max = 35°(there is no strong constraints on ω max , however, it must be large enough for the BLR to be outside the MHD jet). Its luminosity was set with α blr = 0.1 derived from observations (Celotti et al. 1997).
Emission attributed to a hot corona is observed in the SED of 3C 273 in the X-ray band. Following Haardt et al. (1998) who used a model of thermal comptonization to reproduce the hot corona emission, we add an emission following a power law of photon index Γ = 1.65 between 5 × 10 15 Hz and 5 × 10 19  as an extension of the accretion disc represented in Fig. 7. These photons contribute as well to the inverse-Compton emission and Compton rocket process.
Superluminal motion has been observed in 3C 273 by Pearson et al. (1981), proving the existence of relativistic motion in this jet. The deduced apparent velocity is v app ≈ 7.5c 1 . This imposes constraints on the observation angle and one can determine that i obs < 15°(or cos i obs > 0.96). Here we fixed i obs = 13°in agreement with these constraints.
A best-fit of the model on the average data from Turler et al. (1998) is performed. Such a fit is a difficult task in this case as the model is computationally expensive and requires many parameters. We developed a method using a combination of genetic algorithms and gradient methods. As the source parameters have been fixed by observations, the free parameters of the model used for the fit are R 0 , n 0 , B 0 , Q 0 , λ, ω and ζ. R S (cm) 5.3 × 10 14 Z c /R S 10 9 L disc (erg/s) 1.7 × 10 46 Table 1. Parameters corresponding to 3C 273 modelling. R S = 2GM/c 2 is the Schwarzschild radius deduced from a best-fit of the optical emission. R in and R out are the inner and outer radii of the disc. D T and R T are geometrical parameters of the torus. R blr is the BLR radius. Z 0 , R 0 , and Z c are geometrical parameters of the jet (see section 3.1). ω max corresponds to the BLR opening and α blr to its luminosity fraction with respect to that of the disc, L disc . The value z is the redshift of 3C273. The factors n 0 , B 0 , Q 0 , λ, ω and ζ are free parameters of the jet model as described in section 3.1

Results and discussion
With the seven free parameters, the fit gives a reducedχ 2 = 4.92. The model accurately reproduces the entire SED of 3C273 from radio to gamma-rays except for the far radio below 10 8 Hz. These points show a break in the spectrum and are interpreted as the jet hot spot on a very large scale. This hot spot is not supposed to be modelled here and it is actually not surprising that these points are not well fitted by the best fit model; if we do not take them into account, then the reduced χ 2 is reduced toχ 2 = 1.3, more representative of the quality of the fit obtained.
The mass and accretion rate of the black-hole can be deduced from the radius R S and the luminosity L disc determined previously, depending on the BH spin. For a Schwarzschild (non-spinning) black-hole, M • = R isco c 2 /6G = 1.8 × 10 9 M and the reduced accretion rateṁ =Ṁ/Ṁ edd = 0.08. For a maximum spinning Kerr black-hole however, M • = R isco c 2 /2G = 5.4 × 10 9 M andṁ = 0.027. These values of mass are in agreement with those derived from observations; 1 The authors compute an apparent velocity of ∼ 9.5c assuming a value for the Hubble constant H 0 = 55km.s −1 .M pc −1 which is now known to be closer to 70km.s −1 .  In the jet, the low energy (radio to optical) is produced by the synchrotron process. The synchrotron part of the inner jet (below 10 3 R S ) emits in the optical and is hidden by the emission from the accretion disc and the dusty torus. When moving further in the jet, the synchrotron peak shifts to lower frequencies and the further we go, the more the peak shifts. Finally, the whole jet from 10 3 R S to 10 9 R S is necessary to reproduce the power-law-like radio spectrum. Its slope is determined by different factors: the increase of the jet radius, the decrease of the magnetic field and of the particle heating, and the bulk Lorentz factor. The spectrum of 3C 273 shows a break at ∼ 10 9 Hz which is poorly fitted by our model and is thought to be the result of synchrotron emission produced by the extended radio structure (lobe+hot spot). It could also show the limitation of our assumption on ballistic motion of the jet at very large distances. A deceleration of the jet could result in the observed break in the spectrum.
The high-energy emission is produced by inverse-Compton processes, either on synchrotron (SSC) or on thermal photons (EC). Similarly to the synchrotron, the highest energies are produced close to the central engine and further regions emit at lower energies. In particular, the spectrum at ν > 10 21 Hz is produced by regions at z < 10 3 R S with a combination of SSC and EC. However, the X-rays and soft γ-rays are produced further, at z > 10 3 R S by SSC.    further pair creation, the density decreasing only by dilution (due to the increase of the jet radius), and all parameters follow a very smooth evolution, corresponding to the almost featureless spectrum below 10 13 Hz.
Concerning the geometrical parameters, the spine geometry (see equation 1) follows a parabola with a radius R(z) ∝ √ z.
The jet opening can be evaluated by tan θ jet = R z and goes from 2 × 10 −3 rad below 1 parsec to about 10 −5 rad at 1kpc which makes it a very narrow jet (these values concern only the spine jet here which is collimated by the outer MHD jet in the two-flow paradigm). However, the jet radius is important only to compute the SSC radiation, the other process depending only linearly on the number of particles: if the jet were to widen after the SSC emitting region, this would not affect the overall spectrum.
The particle density (middle curve of the top plot in the Fig. 8) evolves along the jet and goes through important pair creation under 10 3 R S . This part is also the part where most of the high-energy emission is emitted. This is in good agreement with the two-flow paradigm as flares (which can be extreme at high energies) are explained by intense phases of pair creation. Slight changes in the particle acceleration could induce large changes in the pair creation which is a highly non-linear process, in turn creating flare episodes. A time-dependent model such as the one described in Boutelier et al. (2008) is currently used to describe these flares.
The particles mean Lorentz factor (3γ) varies between 600 and 1800; it increases and decreases rapidly below 10 3 R S , following changes in the bulk Lorentz factor Γ b . This is due to changes in the cooling inferred by changes in Doppler aberration in the plasma rest frame and in the particle density along the jet. Further in the jet,γ slowly increases but due to an important decrease in the particle density, the particle energy density also decreases (dotted line in the last plot in Fig. 8).
The last plot in Fig. 8 represents the energy density in the particles and in the magnetic field as well as the equipartition ratio defined as the ratio of both energy densities: ξ = n e < γ > m e c 2 B 2 /8π .
This plot shows that the jet is highly magnetised very close to the black-hole (below 100 R S ) as much of the energy is carried by the magnetic field. Particles then carry much more overall energy and bring this energy very far along the jet, allowing for the far synchrotron emission.
The bulk Lorentz factor Γ b varies as a function of the distance in the jet and is equal to Γ eq (z) under 10 4 R S . Changes in Γ b (z) = Γ eq (z) are due to effects discussed in Vuillaume et al. (2015) and are the result of the Compton rocket effect. The jet reaches ballistic motion at z = 10 4 R S as the acceleration becomes ineffective and then Γ b (z) reaches its final value Γ ∞ .
The value of Γ ∞ ≈ 2.7 is very questionable in this modelling, as superluminal motion inferring v app > 7c has been observed, implying at least the same value for Γ b . Superluminal motion as seen by very long base interferometry (VLBI) corresponds to regions very far from the central engine and should thus be produced by components travelling at Γ ∞ . Therefore, Γ ∞ should be able to reach higher values. In the Compton rocket model, this implies that particles should be coupled to the photon field at larger distances than what is found here (well outside the BLR region). There are two simple ways to explain this lack of energy in particles: 1. our description of the acceleration with a power-law might be too simple. It is quite clear that jets are more complex ob-jects and that some re-acceleration sites are present far in the jet, as shown by the complex images of the jet (see also for example HST-1 in M87). Such far acceleration would give an energy boost to particles that would turn into jet acceleration on greater scales, pushing Γ ∞ to the observed values. 2. Our modelling of 3C 273 corresponds to a quiescent state.
During flaring states, it is expected that more energy is transferred to particles, here again providing a possibility to reach higher values of Γ ∞ . In this view, the fastest features observed in jets would correspond to past flaring episodes that accelerated the plasma to high values of Γ b -the latter now following ballistic motion.
On the other hand, it is very interesting that the model is able to reproduce the broadband spectrum of 3C 273 with Γ b (Z) < 3. As explained in section 2.2, there is substantial evidence pointing to low Lorentz factors in AGNs and high Lorentz factors are very difficult to explain if applied to the entire jet. Here we show that with a stratified jet model, high Lorentz factors are not necessary to reproduce the averaged broadband emission and that the quiescent state of an object could correspond to very moderate Γ b in opposition to episodes of flares that could induce higher Γ ∞ , as high as the ones observed.

Conclusion
Emission modelling is essential to understanding the AGN jet physics, and yet despite many flaws the one-zone model is still predominant. However, more complex models are being developed, including structured jets models. These models introduce more physics and improve our understanding of AGNs as they are often better at explaining the observations. In this work, we present a numerical model of a stratified jet based on the twoflow idea first introduced by Sol et al. (1989) and further developed by Henri & Pelletier (1991).
The numerical model is based on a prescription of the geometry of the jet providing the evolution of the radius, magnetic field, and particle acceleration. External thermal sources such as the accretion disc, the dusty torus, and the broad-line region are also modelled, including their spatial geometry. In the jet, we consider emission from synchrotron and inverse-Compton (over synchrotron and thermal photons) processes. The physical conditions such as the particles energy and density are consistently computed along the jet based on initial conditions at a fixed point. In particular, the bulk Lorentz factor is not a free parameter in our model but its evolution is entirely constrained by the external sources (accretion disc, dusty torus, BLR) through the Compton rocket process. Photon-photon absorption is taken into account in the jet (essential for pair production) as well as outside the jet, in its vicinity, and on the photon path to the observer.
Despite being very constrained, the model is able to accurately reproduce the broadband emission of 3C273 from radio to gamma-rays. As the physical conditions evolve along the jet, different zones contribute to different parts of the observed spectrum. Inner parts of the jet contribute more to high-energies while low-frequency radio is emitted further in the jet, at parsec scales.
In the future, the model could be extended to incorporate variability (as done in Boutelier et al. 2008) and can be applied to other types of compact object such as galactic X-ray binaries. Further work is being planned for these applications.
Thomson regime and pile-up particle energy distribution The kernel computed above can then be integrated over a particle distribution. In the case of a pile-up distribution n e (γ) = N e 2γ 3 γ 2 exp (−γ/γ) this integration simplifies and one gets directly the inverse-Compton spectrum: As χ(s) is a single variable function, it can be computed once, and then tabulated and interpolated over when required. This way, the computation of the inverse-Compton spectra from the scattering of a thermal soft photon field on a pile-up distribution of electrons can be done much faster than usually with complete numerical integration.

Appendix B.2: Validating the approximation
Here we compare our analytical approximation with a complete numerical integration of the inverse-Compton spectrum. The result is displayed in figure B.1. We see that the solution proposed by Zdziarski & Pjanka (2013) (hereafter ZP13) is, overall, better than ours with a relative error below 0.1% in the medium range. The emission peak is also closer to the complete calculation.
Nevertheless, our calculation presents the advantage of being faster than the one from ZP13; in our simulations it was shown to be more than twice as fast in the case studied here. Moreover, Article number, page 16 of 17 T. Vuillaume et al.: A stratified jet model for AGN emission in the two-flow paradigm. another large reduction in computation time is also achieved with the one variable function χ(s) in the case of a pile-up distribution in the Thomson regime. We note that the error made compared to the exact numerical integration will be smoothed in the complete model as the spectra are convolved on several parameters and on the integration along the jet.