Issue 
A&A
Volume 620, December 2018



Article Number  A41  
Number of page(s)  16  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/201731899  
Published online  28 November 2018 
A stratified jet model for AGN emission in the twoflow paradigm
^{1}
Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, LAPP, 74000 Annecy, France
email: thomas.vuillaume@lapp.in2p3.fr
^{2}
Univ. Grenoble Alpes, CNRS, IPAG, 38000
Grenoble, France
Received:
5
September
2017
Accepted:
27
September
2018
Context. Highenergy 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 nonthermal emission is the twoflow 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).
Aims. We try to reproduce the broadband spectra of radioloud AGNs with a detailed model of emission taking into account synchrotron and inverseCompton emission by a relativistically moving beam of electron positron, heated by a surrounding turbulent baryonic jet.
Methods. We compute the density and energy distribution of a relativistic pair plasma all along a jet, taking into account the synchrotron and inverseCompton process on the various photon sources present in the core of the AGN, as well as the pair creation and annihilation processes. We use semianalytical approximations to quickly compute the inverseCompton 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 3C 273.
Results. In the case of 3C 273, we obtain an excellent fit of the average broadband energy distribution by assuming physical parameters compatible with known estimates. The asymptotic bulk Lorentz factor is lower than what is observed by superluminal motion, but the discrepancy could be solved by assuming different acceleration profiles along the jet.
Key words: galaxies: active / radiation mechanisms: nonthermal / gamma rays: general / quasars: individual: 3C 273 / galaxies: jets
© ESO 2018
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
Radioloud active galactic nuclei (AGNs) are known to be powerful emitters of nonthermal 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 inverseCompton 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 nonthermal particles, and the precise size and location of the emitting zones. Onezone 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 lowenergy (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 nonthermal emission must be produced in a single zone, which can be an issue for explaining the multiwavelength 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 highenergy 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 blobinjet 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 “twoflow 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 electronpositron 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 magnetohydrodynamic (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 gammaray 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 nonrelativistic jet, it escapes the Compton drag issue. As explained below, the pair beam is only gradually accelerating thanks to the anisotropic inverseCompton 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 gammagamma pair production process, so the dominant emission region can be at large distances from the central black hole, avoiding the problem of gammagamma absorption by the accretion disk photon field. This model therefore offers a natural explanation of the main characteristics of the highenergy 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 highenergy 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 gammaray photons observed from radioloud gammarayemitting AGNs.
The first numerical model of nonthermal emission based on these ideas was proposed by Marcowith et al. (1995), who considered only the inverseCompton process on accretiondisk photons. Assuming a powerlaw particle distribution and a stratified jet, the authors could derive the inverseCompton 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 gammaray 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 twoflow 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 inverseCompton 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 pileup (relativistic Maxwellian) distribution (Saugé & Henri 2004), which better reproduces the highenergy spectra of BL Lacs that peak in the TeV range. A quasiMaxwellian distribution is close to a monoenergetic one, however the spatial convolution of such a distribution whose parameters vary along the jet can mimic a powerlaw over a limited range of frequencies. These authors also proposed a timedependent version of the model that was later shown to be able to successfully reproduce the rapid flares observed in some objects, such as PKS 2155304 (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 selfconsistently 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 pairbeam nonthermal spectra assuming the above configuration, that is, (1) the pair plasma is assumed to be described by a pileup 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 (BLR), and a dusty torus, and (3) particles emit nonthermal radiation through synchrotron, synchrotron selfCompton (SSC), and external Compton (EC) in the various photon fields. Part of the highenergy γray photons can also be absorbed to produce new pairs. We assume that these new pairs are continuously accelerated along the jet. In the twoflow 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 selfconsistently, and the emitted photon spectrum can be evaluated and compared to observations.
The overall layout is as follows. In Sect. 2 we present the main theoretical interest of the twoflow paradigm. Then in Sect. 3 we detail our model and the numerical methods we use. In Sect. 4 we apply this model to the bright quasar 3C 273 and show the model can reproduce its jet emission from radio to gammarays.
2. The twoflow paradigm: the hypothesis and the reasoning behind it
The twoflow paradigm is based on an original idea from Sol et al. (1989; see Sect. 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 selfcollimated 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.
2.1. Interaction of highly relativistic flows
The first benefit of the twoflow hypothesis is to alleviate the problem of the confinement of a highly relativistic flow. Selfconfinement 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 selfconfinement of a highly relativistic jet through this process is quite inefficient. This has been shown first in numerical simulations by Bogovalov (2001) and Bogovalov & 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 twoflow paradigm, the outer selfconfined 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 paircreation process being very efficient and highly nonlinear, 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.
2.2. 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 longstanding 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, 2014; Giroletti et al. 2004; 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 twoflow 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, inverseCompton 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 inverseCompton scattering, therefore killing it before it can be effective.
However, in the twoflow paradigm, as the pairs are continuously reaccelerated 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 Sect. 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 nearEddington 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 Sect. 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 twoflow 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.
3. 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 onezone and emits synchrotron radiation, synchrotron selfCompton (SSC), and external Compton (EC) radiation, computed as in a onezone approximation with a spherical blob of radius R(z). The photon field from external sources (the accretion disc, the dusty torus, and the BLR) is computed all along the jet. This determines the external inverseCompton emission as well as the induced Compton rocket force on the plasma as described below. The opacity to highenergy photons inside and outside the jet is computed numerically, and is used to compute a pair creation term at each step.
Fig. 1. Schematic view of the model developed in the twoflow paradigm. 
3.1. 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 powerlaws:
This law describes a paraboloid with a radius close to R_{0} at a distance Z_{0} from the blackhole. 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 decreases as a powerlaw 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 cutoff, 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 (Eq. (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 nonlinear 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.
3.2. 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 BLR. In order to correctly compute the corresponding inverseCompton radiation, the anisotropy of the sources is taken into account as detailed below.
3.2.1. 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 . From the jet axis at an altitude z, each slice is seen under a solid angle
Fig. 2. Sketch of the disc radial and azimuthal splitting. A slice at (r, φ)∈([R_{in},R_{out}],[0, 2π]) is seen under a solid angle dΩ from an altitude z in the jet. 
at an angle 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 StefanBoltzmann constant, M the blackhole mass, Ṁ the accretion rate, and 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 for R_{in} = R_{isco} and R_{out} ≫ R_{in}.
3.2.2. 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} ∈ [θ_{tmin},θ_{tmax}] and φ ∈ [0, 2π]. Each slice is illuminated by the disk and is in radiative equilibrium. It emits as a greybody 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 pointlike 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:
Fig. 3. Dusty torus sliced as described in the text. Each slice is seen from an altitude z in the jet under a solid angle dΩ. 
Fig. 4. BLR, an optically and geometrically thin shell of clouds seen under a solid angle dΩ from an altitude z in the jet. The BLR is sliced according to ω ∈ [ω_{max},π/2] and φ ∈ [0, 2π]. The BLR absorbs and reemits part of the disk luminosity. 
with ω_{d} the angle between the Zaxis and the emission direction from the disc: 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.
3.2.3. The broad line region
The BLR is modelled by an isotropic, optically and geometrically thin spherical shell of clouds situated at a distance R_{blr} from the central blackhole. 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 greybody spectrum at T = 10^{5} K provides a good approximation to the resulting inverseCompton 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: .
Emissivity of the BLR is then given by:
3.3. Emission processes
As the spine is supposed to be filled by electrons/positrons, only leptonic processes need to be considered here: synchrotron and inverseCompton 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)=(Γ_{b}(z)(1 − β_{b}(z)cosi_{obs})^{−1} into the observer’s rest frame – with and i_{obs} the observer angle relative to the jet axis as defined in Fig. 5.
Fig. 5. Sketch of the computation of the absorption from external sources. Photons from all external sources (accretion disc, BLR and torus) interact with photons emitted by the jet. For every emitting zone at z_{0}, one needs to integrate the absorption over the path of the gamma photons to the observer and over incoming directions of thermal photons. 
In the twoflow paradigm, particles are supposed to be accelerated by a secondorder 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 pileup (or quasiMaxwellian) 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 firstorder Fermi process at work in shocks), the reproduction of largescale power laws in the spectra of the sources is not straight forward. Powerlaw shapes can however be reproduced by a summation of a number of pileup (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 pileup distribution.
3.3.1. Synchrotron radiation
Using the expression of the synchrotronemitted power per unit of frequency and solid angle (Blumenthal & Gould 1970; Saugé 2004) showed that the synchrotron emissivity of a pileup distribution can be written as:
with , 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 (Eqs. (A.1)–(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 pileup distribution of particles, the distribution temperature 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 pileup distribution has three main parts:

ϵ < ϵ_{abs} is the optically thick part of the spectrum described by a powerlaw of index 1.

ϵ_{abs} < ϵ < ϵ_{c} is the optically thin part of the spectrum described by a powerlaw of index −2/3.

ϵ_{c} < ϵ is where the spectrum falls exponentially.
3.3.2. Synchrotron selfCompton radiation
We assume the synchrotron selfCompton (SSC) emission to be cospatial 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
with being the exponential integral function. Interestingly, 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 and :
with Γ being the incomplete gamma function and with and
The SSC photon production rate in the Klein–Nishina regime is then given by
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.
3.3.3. External Compton radiation
The calculation of the inverseCompton 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 Sect. 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. (2014), 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 blackbodies emission by a Wienlike spectrum. Details of the calculation and comparison with ZP13 are done in Appendix B.
The number of inverseCompton 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 , θ_{0} being the angle between the incoming photon and the particle direction of motion, , and 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 pileup 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 inverseCompton spectra from the scattering of a thermal soft photon field on a pileup distribution of electrons can be done much faster than usual with complete numerical integration.
3.4. Photon–photon absorption in the jet and induced pair creation
Highenergy photons produced by inverse Compton will locally interact with lowenergy photons produced by the synchrotron process (external radiations can be locally neglected). This photonphoton 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.
3.4.1. 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 crosssection 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 angleaveraged 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 occurring at x = x_{max}, we make the approximation:
with
with σ_{Th} being the Thomson crosssection.
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 planeparallel approximation for photons absorbed insitu gives their escape probability :
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.
3.4.2. Pair creation
The pair production is a direct consequence of the γ − γ absorption in the jet. A photon of dimensionless highenergy ϵ > 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 rate ṅ_{prod} is given by
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 pileup distribution without considering the perturbation to the energy distribution introduced by the pair creation.
3.5. Evolution of the particle distribution
The pileup distribution has two parameters, the particle density n_{0}(z) and the particles characteristic energy (the particles mean energy being given by ). Both parameters are not imposed but computed consistently all along the jet as detailed below.
3.5.1. Evolution of the particles’ characteristic energy
The particles’ characteristic energy 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.
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 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: U = U_{B} + U_{syn} + U_{ext}.
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 cutoff frequency ν_{KN}). As the synchrotron emission is isotropic and cospatial with the SSC emission, one obtains
with 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 inverseCompton. 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 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.
3.5.2. 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.
In the absence of pair production, the particle flux Φ_{e}(z, t)=∫n_{e}(γ, z, t) πR^{2}(z) Γ_{b}(z) β_{b}(z) c dγ is conserved.
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 Eq. (32). In the stationary case, one simply solves .
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 nonlinear and can lead to explosive events. In the twoflow paradigm, this can lead to fast and powerful flares as the particle density will increase as long as the energy reservoir is not emptied.
3.6. 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 inverseCompton 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 inverseCompton 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 inverseCompton 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 inverseCompton 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 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 only take into account the geometry of the external photon fields in the Thomson regime. We make the assumption that most of the inverseCompton 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 being the relaxation length to equilibrium that can be written as
with 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 (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 . 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 twoflow interactions).
3.7. Absorption outside the jet
Between the jet and the observer, a photon encounters several radiation fields causing absorption through photonphoton 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).
3.7.1. Extragalactic background light
Extragalactic background light (EBL) is composed of several sources: the CMB, the diffuse UVoptical 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}.
3.7.2. Absorption in the jet vicinity
The immediate vicinity of the jet is composed of several sources of soft photons.
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 Eq. (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 Eq. (24) over the path followed by the photons from z_{0} to infinity (here M relates to the position along that path):
Terms in Eq. (44) have been detailed at Eq. (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 photonemitted altitudes, is given in Fig. 6. The figure represents isolines of photonemitted 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.
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. 
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
4. 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. (1999); 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 Sect. 3.6). We detail our modelling of this source in the following section.
4.1. 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 bestfit 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 Eq. (5)) and with the two free parameters: 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} (). The modelling of the disc is given in orange in Fig. 7.
Fig. 7. SED of 3C 273 – modelling compared to data from Turler et al. (1999). The synchrotron emission is shown in blue, the SSC in green, and the external Compton in purple. The torus emission is shown in red and the multicolor accretion disc in orange (filled orange curve), complete with a powerlaw describing the hot corona emission between 0.02 and 200 keV (dashed orange line). Different emission zones in the jet are represented with different dotted lines. The emission below 10^{9} Hz not reproduced by the model is interpreted as the emission from the jet hot spot, not modelled here. 
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).
Parameters corresponding to 3C 273 modelling.
Emission attributed to a hot corona is observed in the SED of 3C 273 in the Xray 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} Hz (0.02–200 keV) as an extension of the accretion disc represented in Fig. 7. These photons contribute as well to the inverseCompton 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 bestfit of the model on the average data from Turler et al. (1999) 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 ζ.
4.2. Results and discussion
With the seven free parameters, the fit gives a reduced . The model accurately reproduces the entire SED of 3C 273 from radio to gammarays 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 , more representative of the quality of the fit obtained.
The mass and accretion rate of the blackhole can be deduced from the radius R_{S} and the luminosity L_{disc} determined previously, depending on the BH spin. For a Schwarzschild (nonspinning) blackhole, M_{•} = R_{isco}c^{2}/6G = 1.8 × 10^{9} M_{⊙} and the reduced accretion rate . For a maximum spinning Kerr blackhole however, M_{•} = R_{isco}c^{2}/2G = 5.4 × 10^{9} M_{⊙} and . These values of mass are in agreement with those derived from observations; Paltani & Turler (2005) determined M_{•} = (5.69 − 8.27)×10^{9} M_{⊙} using the reverberation method on Lyα and C_{IV} lines and M_{•} = (1.58 − 3.45)×10^{9} M_{⊙} using Balmer lines whereas Peterson et al. (2004) determined a mass of M_{•} = (8.8 ± 1.8)×10^{9} M_{⊙} using a large reverberationmapping database.
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 powerlawlike 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 highenergy emission is produced by inverseCompton 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 Xrays and soft γrays are produced further, at z > 10^{3}R_{S} by SSC.
Figure 8 shows the evolution along the jet of the different model parameters. One can see first that the jet reaches ballistic motion (Γ_{b}(z)=Γ_{∞}) at z = 10^{4}R_{S}. From this point, there is no 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.
Fig. 8. Evolution of the physical parameters as a function of the distance in the jet for the modelling of 3C 273 
Concerning the geometrical parameters, the spine geometry (see Eq. (1)) follows a parabola with a radius . The jet opening can be evaluated by and goes from 2 × 10^{−3} rad below 1 parsec to about 10^{−5} rad at 1 kpc 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 twoflow 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 highenergy emission is emitted. This is in good agreement with the twoflow 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 nonlinear process, in turn creating flare episodes. A timedependent model such as the one described in Boutelier et al. (2008) is currently used to describe these flares.
The particles mean Lorentz factor () 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:
This plot shows that the jet is highly magnetised very close to the blackhole (below 100R_{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:

Our description of the acceleration with a powerlaw might be too simple. It is quite clear that jets are more complex objects and that some reacceleration sites are present far in the jet, as shown by the complex images of the jet (see also for example HST1 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.

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.
5. Conclusion
Emission modelling is essential to understanding the AGN jet physics, and yet despite many flaws the onezone 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 BLR are also modelled, including their spatial geometry. In the jet, we consider emission from synchrotron and inverseCompton (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 3C 273 from radio to gammarays. 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 highenergies while lowfrequency 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 Xray binaries. Further work is being planned for these applications.
Acknowledgments
We acknowledge funding support from the French CNES agency and CNRS PNHE. IPAG is part of Labex OSUG@2020 (ANR10 LABX56).
References
 Abramowitz, M. 1974, Handbook of Mathematical Functions, With Formulas, Graphs, and Mathematical Tables (New York: Dover Publications, Incorporated) [Google Scholar]
 Aharonian, F., Akhperjanian, A. G., BazerBachi, A. R., et al. 2007, ApJ, 664, L71 [Google Scholar]
 Begelman, M. C., Fabian, A. C., & Rees, M. J. 2008, MNRAS, 384, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 [NASA ADS] [CrossRef] [Google Scholar]
 Blumenthal, G. R., & Gould, R. J. 1970, Rev. Mod. Phys., 42, 237 [NASA ADS] [CrossRef] [Google Scholar]
 Bogovalov, S. V. 2001, A&A, 371, 1155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bogovalov, S., & Tsinganos, K. 2001, Astron. Astrophys. Trans., 20, 303 [NASA ADS] [CrossRef] [Google Scholar]
 Boutelier, T. 2009, tel.archivesouvertes.fr [Google Scholar]
 Boutelier, T., Henri, G., & Petrucci, P.O. 2008, MNRAS, 390, L73 [NASA ADS] [CrossRef] [Google Scholar]
 Celotti, A., Padovani, P., & Ghisellini, G. 1997, MNRAS, 286, 415 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chiaberge, M., Celotti, A., Capetti, A., & Ghisellini, G. 2000, A&A, 358, 104 [NASA ADS] [Google Scholar]
 Coppi, P. S., & Blandford, R. D. 1990, MNRAS, 245, 453 [NASA ADS] [Google Scholar]
 Dubus, G., Cerutti, B., & Henri, G. 2010, A&A, 516, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Franceschini, A., Rodighiero, G., & Vaccari, M. 2008, A&A, 487, 837 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Georganopoulos, M., & Kazanas, D. 2003, ApJ, 594, L27 [NASA ADS] [CrossRef] [Google Scholar]
 Ghisellini, G., Tavecchio, F., & Chiaberge, M. 2005, A&A, 432, 401 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Giroletti, M., Giovannini, G., Feretti, L., et al. 2004, ApJ, 600, 127 [NASA ADS] [CrossRef] [Google Scholar]
 Gould, R. J., & Schréder, G. P. 1967, Phys. Rev., 155, 1404 [NASA ADS] [CrossRef] [Google Scholar]
 Haardt, F., Fossati, G., Grandi, P., et al. 1998, A&A, 340, 35 [NASA ADS] [Google Scholar]
 Henri, G., & Pelletier, G. 1991, ApJ, 383, L7 [NASA ADS] [CrossRef] [Google Scholar]
 Henri, G., & Saugé, L. 2006, ApJ, 640, 185 [NASA ADS] [CrossRef] [Google Scholar]
 Hervet, O., Boisson, C., & Sol, H. 2015, A&A, 578, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jones, F. C. 1968, Phys. Rev., 167, 1159 [NASA ADS] [CrossRef] [Google Scholar]
 Katarzynski, K., Sol, H., & Kus, A. 2001, A&A, 367, 809 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Khangulyan, D., Aharonian, F. A., & Kelner, S. R. 2014, ApJ, 783, 100 [NASA ADS] [CrossRef] [Google Scholar]
 Lister, M. L., Cohen, M. H., Homan, D. C., et al. 2009, AJ, 138, 1874 [NASA ADS] [CrossRef] [Google Scholar]
 Lister, M. L., Aller, M. F., Aller, H. D., et al. 2013, AJ, 146, 120 [NASA ADS] [CrossRef] [Google Scholar]
 Mahadevan, R., Narayan, R., & Yi, I. 1996, ApJ, 465, 327 [NASA ADS] [CrossRef] [Google Scholar]
 Marcowith, A., Henri, G., & Pelletier, G. 1995, MNRAS, 277, 681 [NASA ADS] [Google Scholar]
 Marcowith, A., Pelletier, G., & Henri, G. 1996, A&AS, 120, 563 [NASA ADS] [Google Scholar]
 O’dell, S. L. 1981, ApJ, 243, L147 [NASA ADS] [CrossRef] [Google Scholar]
 Paltani, S., & Turler, M. 2005, A&A, 435, 811 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pearson, T. J., Unwin, S. C., Cohen, M. H., et al. 1981, Nature, 290, 365 [NASA ADS] [CrossRef] [Google Scholar]
 Pelletier, G. 2004, ArXiv eprints [arXiv:astroph/0405113] [Google Scholar]
 Peterson, B. M., Ferrarese, L., Gilbert, K. M., et al. 2004, ApJ, 613, 682 [NASA ADS] [CrossRef] [Google Scholar]
 Phinney, E. S. 1982, MNRAS, 198, 1109 [NASA ADS] [CrossRef] [Google Scholar]
 Piner, B. G., & Edwards, P. G. 2004, ApJ, 600, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Piner, B. G., & Edwards, P. G. 2014, ApJ, 797, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Reimer, A. 2007, ApJ, 665, 1023 [NASA ADS] [CrossRef] [Google Scholar]
 Renaud, N. 1999, tel.archivesouvertes.fr [Google Scholar]
 Renaud, N., & Henri, G. 1998, MNRAS, 300, 1047 [NASA ADS] [CrossRef] [Google Scholar]
 Rieke, G. H., & Weekes, T. C. 1969, ApJ, 155, 429 [NASA ADS] [CrossRef] [Google Scholar]
 Rybicki, G. B., & Lightman, A. P. 1979, Radiative Processes in Astrophysics (New York: WileyInterscience), 393 [Google Scholar]
 Rybicki, G. B., & Lightman, A. P. 1986, Radiative Processes in Astrophysics (New York: WileyVCH), 400 [Google Scholar]
 Saugé, L. 2004, tel.archivesouvertes.fr [Google Scholar]
 Saugé, L., & Henri, G. 2004, ApJ, 616, 136 [NASA ADS] [CrossRef] [Google Scholar]
 Schlickeiser, R. 1984, A&A, 136, 227 [NASA ADS] [Google Scholar]
 Schmidt, M. 1963, Nature, 197, 1040 [NASA ADS] [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1976, MNRAS, 175, 613 [NASA ADS] [CrossRef] [Google Scholar]
 Sol, H., Pelletier, G., & Asseo, E. 1989, MNRAS, 237, 411 [NASA ADS] [Google Scholar]
 Tavecchio, F., & Ghisellini, G. 2008, MNRAS, 386, 945 [NASA ADS] [CrossRef] [Google Scholar]
 Tavecchio, F., Ghisellini, G., Ghirlanda, G., Foschini, L., & Maraschi, L. 2010, MNRAS, 401, 1570 [NASA ADS] [CrossRef] [Google Scholar]
 Turler, M., Paltani, S., Courvoisier, T. J. L., et al. 1999, A&AS, 134, 89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vuillaume, T., Henri, G., & Petrucci, P. O. 2015, A&A, 581, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zdziarski, A. A., & Pjanka, P. 2013, MNRAS, 436, 2950 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Computation of the synchrotron and synchrotronself Compton radiation
A.1. Synchrotron radiation
The function Λ presented in Eq. (10) can be approximated analytically in the different regimes:

A development in Taylor series for small y gives:

Mahadevan et al. (1996) proposed an approximation for intermediates values of y:

For large values of y, with a saddlepoint method, (Saugé 2004) obtains:
A.2. Synchrotronself Compton radiation
In this section, we consider the dimensionless energy (as introduced before) ϵ = hν/m_{e}c^{2}. The photon production rate resulting from the scattering of a synchrotron photon field n_{ph}(ϵ) on a population of particles n(γ) cnan be written as:
with K_{jones} being the Compton kernel for an isotropic source of soft photons, considering the full Klein–Nishina cross section in the headon approximation as evaluated by Jones (1968). Numerical evaluation of this equation is time consuming but approximations may be used to speed up calculations. To do so, we consider the Thomson and the Klein–Nishina regimes separately:
# Emission in the Thomson regime
In this regime, Blumenthal & Gould (1970) and Rybicki & Lightman (1979) showed that the Jones kernel can be simplified:
with f(x)=2xlnx + x + 1 − 2x^{2}.
Injecting this simplified kernel into Eq. (A.4) and considering a pileup distribution (see Eq. (9)) for the particles, one obtains
with
Using the approximation for the function f proposed by Rybicki & Lightman (1979) f(x)≈2(1 − x)/3, one can show that the integral can be analytically solved and gives
with
where Γ(a, x) is the incomplete gamma function introduced by Abramowitz (1974).
Finally, the SSC photon production rate in the Thomson regime can be written as
with the function
and the exponential function
As one can see, in Eq. (A.10), the function depends only on the ratio . This function can therefore be evaluated and tabulated to decrease the computing time.
# Emission in the Klein–Nishina regime
In this regime, the complete Jones kernel must be considered. However, Rieke & Weekes (1969) showed that some approximations were possible, introducing the differential crosssection σ(ϵ_{1}, ϵ_{0}, γ), one can write
in the Thomson (1) and the Klein–Nishina (2) regimes, respectively,
Then the photon production rate can be written as
where g_{i} must be calculated in the two diffusion regimes.
In the Thomson regime, one can show easily that for a pileup distribution,
Injecting this expression into (A.15) and assuming a power law for the soft photon spectrum dn_{ph}/dϵ = n_{s}ϵ^{−s} over the range [ϵ_{min}, ϵ_{max}] after some manipulations one obtains
with Γ the incomplete gamma function and with and
In the Klein–Nishina regime, with a pileup distribution injected into (A.15), one obtains
With the assumption of a power law for the soft photon spectrum, the production rate becomes
with
As shown before, the synchrotron spectrum resulting from a pileup distribution is a broken power law, dn_{ph}/dϵdt = n_{s}ϵ^{−s}, with s = −1 in the optically thick part [ϵ_{min}, ϵ_{abs}] and s = 2/3 in the optically thin part [ϵ_{abs}, ϵ_{max}]. One must therefore consider two cases depending on the position of the absorption energy ϵ_{abs} relative to the Thomson/Klein–Nishina threshold 1/ϵ_{1}. The SSC spectrum in the Klein–Nishina regime is finally given by

If ϵ_{abs} < 1/ϵ_{1},
which gives

If ϵ_{abs} > 1/ϵ_{1},
which gives
Appendix B: Computation of the inverseCompton scattering over a thermal distribution of photons
The number of scattered photons of energy per element of time, per energy, per energy of the incoming photon ϵ′ and per elementary solid angle is equal to the number of photons at the energy ϵ′=hν′ verifying multiplied by their probability to scatter in an emission solid angle (see Jones 1968):
The distribution of photons, given by , should be given by Planck’s law in the case of a radiative black body. The approximation that we propose for the inverseCompton scattering is based on an approximation of Planck’s law that we call modified Wien’s law:
leaving the parameters A(T) and to be determined.
One can see the modified Wien’s law as a mix of Rayleigh–Jeans law and of Wien’s law. An (arbitrary) criteria to determine A(T) and is to conserve the total emitted power and the position of the emission peak in the Thomson regime. In this case, one can show that
with α ≈ 2.821 being a numerical constant.
When replacing with in Eq. B.1 and using the highly relativistic approximation (γ ≫ 1) and the headon approximation (Blumenthal & Gould 1970), one can show that the number of inverseCompton photons scattered in the particle rest frame per energy unit and per time unit (denoted by ′ subscript) for a thermal distribution of soft photons of energy γm_{e}c^{2} is given by
with , θ_{0} being the angle between the incoming photon and the particle direction of motion, and being the exponential integral that can be easily computed by any numerical package.
B.1. Approximations in the Thomson regime
In the Thomson regime, the expression (B.4) can be simplified, especially with and and one obtains
Therefore, the spectrum of the inverseCompton scattering of a thermal distribution of photons on a single particle in the Thomson regime is a function of a single variable.
Thomson regime and pileup particle energy distribution
The kernel computed above can then be integrated over a particle distribution. In the case of a pileup distribution this integration simplifies and one gets directly the inverseCompton spectrum:
with
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 inverseCompton spectra from the scattering of a thermal soft photon field on a pileup distribution of electrons can be done much faster than usually with complete numerical integration.
B.2. Validating the approximation
Here we compare our analytical approximation with a complete numerical integration of the inverseCompton spectrum. The result is displayed in Fig. B.1.
Fig. B.1. Comparison of the different inverseCompton spectra of a thermal soft photon field (T = 5780 K) scattering on a monoenergetic particle distribution of Lorentz factor γ = 1e2 with the relative errors in dotdotdashed thin lines. The numerically integrated spectra with the thermal distribution described by Planck’s law is shown in black. The analytically integrated spectra with the thermal distribution described by the modified Wien’s law (given by Eq. (21)) is shown in blue. As a comparison, we also provide the spectrum from Zdziarski & Pjanka (2013) in orange. 
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, another large reduction in computation time is also achieved with the one variable function χ(s) in the case of a pileup 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.
All Tables
All Figures
Fig. 1. Schematic view of the model developed in the twoflow paradigm. 

In the text 
Fig. 2. Sketch of the disc radial and azimuthal splitting. A slice at (r, φ)∈([R_{in},R_{out}],[0, 2π]) is seen under a solid angle dΩ from an altitude z in the jet. 

In the text 
Fig. 3. Dusty torus sliced as described in the text. Each slice is seen from an altitude z in the jet under a solid angle dΩ. 

In the text 
Fig. 4. BLR, an optically and geometrically thin shell of clouds seen under a solid angle dΩ from an altitude z in the jet. The BLR is sliced according to ω ∈ [ω_{max},π/2] and φ ∈ [0, 2π]. The BLR absorbs and reemits part of the disk luminosity. 

In the text 
Fig. 5. Sketch of the computation of the absorption from external sources. Photons from all external sources (accretion disc, BLR and torus) interact with photons emitted by the jet. For every emitting zone at z_{0}, one needs to integrate the absorption over the path of the gamma photons to the observer and over incoming directions of thermal photons. 

In the text 
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. 

In the text 
Fig. 7. SED of 3C 273 – modelling compared to data from Turler et al. (1999). The synchrotron emission is shown in blue, the SSC in green, and the external Compton in purple. The torus emission is shown in red and the multicolor accretion disc in orange (filled orange curve), complete with a powerlaw describing the hot corona emission between 0.02 and 200 keV (dashed orange line). Different emission zones in the jet are represented with different dotted lines. The emission below 10^{9} Hz not reproduced by the model is interpreted as the emission from the jet hot spot, not modelled here. 

In the text 
Fig. 8. Evolution of the physical parameters as a function of the distance in the jet for the modelling of 3C 273 

In the text 
Fig. B.1. Comparison of the different inverseCompton spectra of a thermal soft photon field (T = 5780 K) scattering on a monoenergetic particle distribution of Lorentz factor γ = 1e2 with the relative errors in dotdotdashed thin lines. The numerically integrated spectra with the thermal distribution described by Planck’s law is shown in black. The analytically integrated spectra with the thermal distribution described by the modified Wien’s law (given by Eq. (21)) is shown in blue. As a comparison, we also provide the spectrum from Zdziarski & Pjanka (2013) in orange. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.