TIME-DEPENDENT MODELING OF RADIATIVE PROCESSES IN HOT MAGNETIZED PLASMAS

The dimensionless four-momentum of a photon is $\underline{x}=\lbrace x, {\bm x}\rbrace = x \lbrace 1,{\bm \omega}\rbrace$, where ω is the unit vector in the photon propagation direction and x ≡ hν/m_ec². The photon distribution can be described by the occupation number $\tilde{n}_{\rm ph}$ or by the photon number density per linear and logarithmic interval of photon energy:

$$\begin{eqnarray} {N_{\rm ph}} &=& \int {N_{\rm ph}}(x) \; {d}x = \int {n_{\rm ph}}(x) \; {d}\ln {x} = \frac{2}{\lambda _{\rm C}^3} \int {d}^2 \omega \int \tilde{n}_{\rm ph}({\bm x})\ x^2\ {d}x, \end{eqnarray} \tag{ 1 }$$

where λ_C = h/m_ec is the Compton wavelength. Functions N_ph(x) and $\tilde{n}_{\rm ph}$ are used in general forms of kinetic equations and n_ph(x) is convenient for numerical work.

The dimensionless electron (positron) four-momentum is defined as $\underline{p}= \lbrace \gamma, {\bm p}\rbrace = \lbrace \gamma, p{\bOmega}\rbrace = \gamma \lbrace 1, \beta {\bOmega}\rbrace$, where ${\bOmega}$ is the unit vector in the electron propagation direction, γ, β, and $p=\beta \gamma =\sqrt{\gamma ^2-1}$ are the electron Lorentz factor, dimensionless velocity, and momentum, respectively. We can use subscripts + and − to distinguish between positrons and electrons. The electron/positron distributions can be defined in a number of alternative ways (normalized to their number density):

$$\begin{eqnarray} &&N_{\pm } = \int N_{\pm }(\gamma) {d}\gamma = \int {n_{\pm }}(p) \; {d}\ln {p} = \frac{2}{\lambda _{\rm C}^3} \int {d}^2 \Omega \int \tilde{n}_{\pm }({\bm p}) \; p^2 {d}p. \end{eqnarray} \tag{ 2 }$$

The occupation number $\tilde{n}_{\pm }(p)$ and the density per unit Lorentz factor are useful quantities used in general kinetic equations, while the electron density per logarithmic momentum interval n_±(p) is more appropriate for numerical work. For the processes, where the distinction between electrons and positrons is unnecessary, we use the sum of the distributions, for example, n_e = n_- + n₊.

The relativistic kinetic equation (RKE) describing the evolution of the occupation number $\tilde{n}_1 ({\bm p}_1)$ of species 1 (electron or photon) as a result of binary collisions can be written in the covariant form (de Groot et al. 1980)

$$\begin{eqnarray} &&\underline{p}_1 \cdot \underline{\nabla }\tilde{n}_1({\bm p}_1) = \int \frac{{d}^3 p_2}{\epsilon _2}\frac{{d}^3 p_3}{\epsilon _3} \frac{{d}^3 p_4}{\epsilon _4} \delta (\underline{p}_1+\underline{p}_2-\underline{p}_3-\underline{p}_4) W_{12\rightarrow 34}[ \tilde{n}_3({\bm p}_3) \tilde{n}_4({\bm p}_4)- \tilde{n}_1({\bm p}_1) \tilde{n}_2({\bm p}_2)],\qquad\quad \end{eqnarray} \tag{ 3 }$$

where $\underline{\nabla }=\lbrace \partial /c \partial t, {\bm \nabla}\rbrace$ is the four-gradient, _i is the zeroth component of the corresponding four-momentum, and W_12→34 = W_34→12 is a Lorentz scalar transition rate, which possesses the obvious symmetry. In this equation, the nonlinear terms related to fermion degeneracy and induced photon scattering are omitted. As it stands, the right-hand side of Equation (3) accounts for the rate of one particular process. To determine the full evolution of $\tilde{n}_1$ we should therefore sum up the collisional integrals accounting for all relevant processes.

In the frame where the particle distributions are isotropic (we call this frame E), the kinetic equation can be represented in the form (skipping subscript 1)

$$\begin{eqnarray} &&\begin{array}{c} \dsty \frac{\partial \tilde{n}(p)}{\partial t} + \frac{1}{p^2}\frac{\partial }{\partial p}\left\lbrace \dot{\epsilon } \epsilon p \tilde{n}(p) - \frac{1}{2} \frac{\partial }{\partial \epsilon } \left[ D(\epsilon) \epsilon p \tilde{n}(p) \right] \right\rbrace = \left. \frac{D \tilde{n}(p)}{Dt} \right|_{\rm coll},\end{array} \end{eqnarray} \tag{ 4 }$$

where the momentum derivative term accounts for continuous energy gain/loss processes, while the right-hand side contains all discontinuous processes such as scattering, emission, absorption, and escape. The quantities $\dot{\epsilon }$ and D() account for systematic particle heating/cooling and diffusion in energy space, respectively. Both are generally energy-dependent for the processes we are considering here. For the following discussion it is convenient to decompose the kinetic equations in terms of the contributions from different physical processes as

$$\begin{eqnarray} &&\frac{\partial {n_{\rm ph}}(x)}{\partial t} = \dot{n}_{\rm ph,syn} (x)+ \dot{n}_{\rm ph,cs}(x) + \dot{n}_{\rm ph,pp} (x) - \frac{{n_{\rm ph}}(x)}{t_{\rm esc}} + Q_{\rm ph}, \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} &&\frac{\partial {n_{\pm }}(p)}{\partial t} = \dot{n}_{\rm \pm,syn}(p) +\dot{n}_{\rm \pm,cs}(p) + \dot{n}_{\rm \pm,pp} (p) + \dot{n}_{\rm \pm,Coul} (p)- \frac{{n_{\pm }}(p)}{t_{\rm \pm, esc}} + Q_{\rm \pm }, \end{eqnarray} \tag{ 6 }$$

where syn, cs, pp, and Coul stand for synchrotron, Compton scattering, pair production (and annihilation), and Coulomb scattering, respectively. The terms describing physical processes can contain both differential and integral parts, depending on the nature of the process and the way we find most convenient to treat it. Thus the equation for photons has the form:

$$\begin{eqnarray} &&\frac{\partial {n_{\rm ph}}(x)}{\partial t} = -\frac{\partial }{\partial \ln x} \left[ A_{\rm ph}(x) {n_{\rm ph}}(x) - B_{\rm ph}(x) \frac{\partial {n_{\rm ph}}(x)}{\partial \ln x} \right]+ \int K_{\rm ph}(x,x_1) {n_{\rm ph}}(x_1) \: {d}\ln x_1 - \frac{{n_{\rm ph}}(x)}{t_{\rm ph}} + S_{\rm ph}. \end{eqnarray} \tag{ 7 }$$

Here, the differential term is responsible for Compton scattering in diffusion approximation, while the integral term with kernel K_ph describes scattering that can be resolved on the grid. The sink term ∝1/t_ph describes photon absorption (by synchrotron and pair production) and scattering as well as the escape, while S_ph gives the contribution from pair annihilation, synchrotron emission, and other (e.g., blackbody) photon injections.

Similarly for electrons and positrons, we write

$$\begin{eqnarray} &&\frac{\partial {n_{\pm }}(p)}{\partial t} = -\frac{\partial }{\partial \ln p} \left[ A_{\rm e}(p) {n_{\pm }}(p) - B_{\rm e}(p) \frac{\partial {n_{\pm }}(p)}{\partial \ln p} \right]+ \int K_{\rm e}(p,p_1) {n_{\pm }}(p_1) \: {d}\ln p_1 - \frac{{n_{\pm }}(p)}{t_{\rm \pm }} + S_{\rm \pm }, \end{eqnarray} \tag{ 8 }$$

where coefficients A_e and B_e describe electron cooling, heating, and diffusion as a result of synchrotron emission and absorption, Compton scattering in Thomson limit, Coulomb scattering as well as possible diffusive particle acceleration. The integral term with kernel K_e describes Compton scattering in Klein–Nishina limit into the bin and the sink term ∝1/t_± gives the scattering from the bin as well as the electron escape and pair annihilation. The source term S_± contains pair production as well as a possible electron injection term.

As we are studying radiative processes in a simple one-zone framework neglecting the radiative transport effects, we must include an escape term in Equation (5) to allow for the fact that photons can leave the emission region of finite size R and produce the radiation flux that is actually observed. The typical escape timescale is usually estimated from random walk arguments resulting in t_esc ∼ R(1 + τ_sc)/c, where τ_sc is the scattering opacity. Such form accounts for the fact that if multiple scatterings are important (τ_sc>1), photons have to "diffuse" out of the medium and the escape time is prolonged by a factor τ_sc. However, it does not account for the fact that if the medium is absorptive, a typical photon cannot diffuse further than the thermalization length l_⋆ = [α_a(α_a + α_sc)]^−1/2 before it is destroyed (α_a and α_sc are extinction coefficients due to absorption and scattering, respectively). To incorporate both effects, we employ the solution of a simple radiative diffusion problem in a sphere of radius R with constant emissivity, absorptivity, and monochromatic scattering. The escape timescale is estimated by comparing the emergent flux to the radiation density inside the source. While clearly an oversimplification, such estimation nevertheless has the desired properties mentioned above.

Defining the effective optical thickness of the medium as $\tau _*= \sqrt{3 \tau _{\rm a}(\tau _{\rm a}+ \tau _{\rm sc})}$, where τ_a = α_aR and τ_sc = α_scR are optical thicknesses due to absorption and scattering, respectively, we find

$$\begin{equation} t_{\rm esc} \,{=}\, \frac{2R}{3c} \left\lbrace 1 \,{+}\, \frac{\sqrt{3}}{2\sqrt{1-\lambda }} \left[\frac{\tau _*(1-{\rm e}^{-2\tau _*})}{\tau _* (1+{\rm e}^{-2\tau _*}) - (1-{\rm e}^{-2\tau _*})} \,{-}\, \frac{3}{\tau _*} \right] \right\rbrace\!, \end{equation} \tag{ 9 }$$

where λ = α_sc/(α_a + α_sc) is the single-scattering albedo. If the medium is translucent (τ_* ≪ 1), Equation (9) reduces to a more familiar form

$$\begin{equation} t_{\rm esc}= \frac{2R}{3c} \left(1 + \frac{3}{10} \tau _{\rm sc}\right). \end{equation} \tag{ 10 }$$

In our simulations, α_a includes cyclosynchrotron absorption and photon–photon pair production, and α_sc is the extinction coefficient for Compton scattering.

2.4.1. Compton Scattering of Photons

The explicitly covariant form of RKE for Compton scattering of photons ignoring nonlinear terms is (Pomraning 1973; Nagirner & Poutanen 1994, hereafter NP94)

$$\begin{eqnarray} &&\underline{x}\cdot \underline{\nabla }\tilde{n}_{\rm ph}({\bm x}) = \frac{r_{\rm e}^2}{2} \frac{2}{\lambda _{\rm C}^3} \int \frac{{d}^3 p}{\gamma } \frac{{d}^3 p_1}{\gamma _1} \frac{{d}^3 x_1}{x_1}\:\delta (\underline{p}_1 + \underline{x}_1 - \underline{p}- \underline{x}) F [ \tilde{n}_{\rm ph}({\bm x}_1) \tilde{n}_{\rm e}({\bm p}_1)- \tilde{n}_{\rm ph}({\bm x}) \tilde{n}_{\rm e}({\bm p})], \end{eqnarray} \tag{ 11 }$$

where r_e is the classical electron radius, F is the Klein–Nishina reaction rate (Berestetskii et al. 1982)

$$\begin{equation} F = \left(\frac{1}{\xi } - \frac{1}{\xi _1}\right)^2 + 2 \; \left(\frac{1}{\xi } - \frac{1}{\xi _1}\right) + \frac{\xi }{\xi _1} + \frac{\xi _1}{\xi }, \end{equation} \tag{ 12 }$$

and $\xi = \underline{p}_1\cdot \underline{x}_1= \underline{p}\cdot \underline{x}$ and $\xi _1 = \underline{p}_1\cdot \underline{x}= \underline{p}\cdot \underline{x}_1$ are the scalar products of four-vectors.

We assume the existence of a reference frame where the particle and photon distributions are approximately homogeneous and isotropic. Under the spacial homogeneity assumption we can write Equation (11) as

$$\begin{eqnarray} &&\left. \frac{D\tilde{n}_{\rm ph}({\bm x})}{Dt} \right|_{\rm coll,cs} = {-} c \: \sigma _{\rm T}\: \overline{s}_0({\bm x}) \: {N_{\rm e}}\: \tilde{n}_{\rm ph}({\bm x}) + c\sigma _{\rm T}{N_{\rm e}}\frac{1}{x} \int \frac{{d}^3 x_1}{x_1} \: R_{\rm ph}({\bm x}_1 \rightarrow {\bm x}) \: \tilde{n}_{\rm ph}({\bm x}_1). \end{eqnarray} \tag{ 13 }$$

The scattering cross-section (in units of Thomson cross-section σ_T) is given by

$$\begin{eqnarray} &&\overline{s}_0({\bm x}) = \frac{3}{16\pi } \frac{2}{\lambda _{\rm C}^3{N_{\rm e}}} \frac{1}{x} \int \frac{{d}^3 p}{\gamma } \frac{{d}^3 p_1}{\gamma _1} \frac{{d}^3 x_1}{x_1} \tilde{n}_{\rm e}({\bm p}) F \: \delta (\underline{p}_1 + \underline{x}_1 - \underline{p}- \underline{x}) \end{eqnarray} \tag{ 14 }$$

and the redistribution function is

$$\begin{eqnarray} &&R_{\rm ph}({\bm x}_1 \rightarrow {\bm x}) = \frac{3}{16\pi } \frac{2}{\lambda _{\rm C}^3{N_{\rm e}}} \int \frac{{d}^3 p}{\gamma } \frac{{d}^3 p_1}{\gamma _1} \: \tilde{n}_{\rm e}({\bm p}_1) F \; \delta (\underline{p}_1 + \underline{x}_1 - \underline{p}- \underline{x}). \end{eqnarray} \tag{ 15 }$$

For isotropic particle distributions in frame E, Equation (13) can be written as

$$\begin{eqnarray} &&\left. \frac{D \tilde{n}_{\rm ph}(x)}{Dt} \right|_{\rm coll,cs} = {-} c \: \sigma _{\rm T}\: \overline{s}_0(x) \: {N_{\rm e}}\: \tilde{n}_{\rm ph}(x) + c\sigma _{\rm T}{N_{\rm e}}\frac{4\pi }{x} \int x_1 {d}x_1 \: \overline{R}_{\rm ph}(x,x_1) \: \tilde{n}_{\rm ph}(x_1), \end{eqnarray} \tag{ 16 }$$

where the redistribution function averaged over the cosine of the scattering angle μ = x  ·  x₁/(xx₁) = ω  ·  ω₁ is expressed via an integral over the electron distribution (NP94):

$$\begin{eqnarray} &&\overline{R}_{\rm ph}(x,x_1) = \frac{1}{2} \int _{-1}^{1} R_{\rm ph}({\bm x}_1 \rightarrow {\bm x}) \: {d}\mu = \frac{3}{16} \frac{2}{\lambda _{\rm C}^3{N_{\rm e}}} \int _{\gamma _{\star }(x,x_1)}^{\infty } \overline{R}_{\rm ph}(x,x_1,\gamma _1) \tilde{n}_{\rm e}(p_1) \ {d}\gamma _1. \end{eqnarray} \tag{ 17 }$$

Here,

$$\begin{eqnarray} &&\overline{R}_{\rm ph}(x,x_1,\gamma _1) = \frac{1}{4\pi ^2}\: p_1 \int \frac{{d}^3 p}{\gamma } {d}^2\Omega _1 \: {d}^2\omega _1 \: F \: \delta (\underline{p}_1 + \underline{x}_1 - \underline{p}- \underline{x}), \end{eqnarray} \tag{ 18 }$$

and the lower limit of the second integral in Equation (17) comes from the condition of energy and momentum conservation:

$$\begin{eqnarray} \gamma _{\star }(x,x_1)= \left\lbrace \begin{array}{@{}l@{\quad }l@{}} [ x - x_1 + (x+x_1)\sqrt{1+1/xx_1} ]/2 &\mbox{if \;$|x-x_1| \ge 2xx_1$,} \\ 1 + \left(x - x_1 + |x - x_1| \right)/2 & \mbox{if \;$|x-x_1| \le 2xx_1$.} \end{array}\right. \end{eqnarray} \tag{ 19 }$$

The integrals in Equation (18) can be calculated analytically (Brinkmann 1984; NP94) to obtain a fully general expression for $\overline{R}_{\rm ph}(x,x_1,\gamma _1)$ valid in all regimes (see Appendix B). This is an alternative form of the function derived by Jones (1968).

Since the total number of particles is conserved in Compton scattering, multiplying the right-hand side of Equation (16) by x² integrating over dx must give zero, implying a relation between the redistribution function and the extinction coefficient (NP94)

$$\begin{equation} \overline{s}_0(x) = \frac{4\pi }{x} \int _0^{\infty } \overline{R}_{\rm ph}(x_1,x) \; x_1 \; {d}x_1. \end{equation} \tag{ 20 }$$

This can also be inferred directly from the definitions (14) and (15).

In the kinetic equation (Equation (5)) for the photon density n_ph(x) the Compton term is obtained by multiplying Equation (16) by 8πλ⁻³_Cx³.

2.4.2. Compton Scattering of Electrons and Positrons

The description of Compton scattering for electrons and positrons is very similar to that for photons. In the linear approximation the RKE reads

$$\begin{eqnarray} &&\underline{p}\cdot \underline{\nabla }\tilde{n}_{\pm }({\bm p}) = \frac{r_{\rm e}^2}{2} \frac{2}{\lambda _{\rm C}^3} \int \frac{{d}^3 x}{x} \frac{{d}^3 x_1}{x_1} \frac{{d}^3 p_1}{\gamma _1} \: \delta (\underline{p}_1 + \underline{x}_1 - \underline{p}- \underline{x}) F [ \tilde{n}_{\rm ph}({\bm x}_1) \tilde{n}_{\pm }({\bm p}_1)- \tilde{n}_{\rm ph}({\bm x}) \tilde{n}_{\pm }({\bm p}) ]. \end{eqnarray} \tag{ 21 }$$

Neglecting spatial gradients, Equation (21) becomes

$$\begin{eqnarray} &&\left. \frac{D \tilde{n}_{\pm }({\bm p})}{Dt} \right|_{\rm coll,cs} = {-} c \: \sigma _{\rm T}\: \overline{s}_0({\bm p}) \: {N_{\rm ph}}\: \tilde{n}_{\pm }({\bm p}) + c \sigma _{\rm T}{N_{\rm ph}}\frac{1}{\gamma } \int \frac{{d}^3 p_1}{\gamma _1} \: R_{\rm e}({\bm p}_1 \rightarrow {\bm p}) \: \tilde{n}_{\pm }({\bm p}_1), \end{eqnarray} \tag{ 22 }$$

where the scattering cross-section for electrons is

$$\begin{eqnarray} &&\overline{s}_0({\bm p})= \frac{3}{16\pi } \frac{2}{\lambda _{\rm C}^3{N_{\rm ph}}} \frac{1}{\gamma } \int \frac{{d}^3 x}{x} \frac{{d}^3 x_1}{x_1} \frac{{d}^3 p_1}{\gamma _1} \tilde{n}_{\rm ph}({\bm x}) F \: \delta (\underline{p}_1 + \underline{x}_1 - \underline{p}- \underline{x}) \end{eqnarray} \tag{ 23 }$$

and the redistribution function

$$\begin{eqnarray} &&R_{\rm e}({\bm p}_1 \rightarrow {\bm p}) = \frac{3}{16\pi } \frac{2}{\lambda _{\rm C}^3{N_{\rm ph}}} \int \frac{{d}^3 x}{x} \frac{{d}^3 x_1}{x_1} \: \tilde{n}_{\rm ph}({\bm x}_1) F \; \delta (\underline{p}_1 + \underline{x}_1 - \underline{p}- \underline{x}). \end{eqnarray} \tag{ 24 }$$

Making use of the isotropy of the problem, we can rewrite the kinetic equation in frame E for isotropic distribution $\tilde{n}_{\pm}(p)$:

$$\begin{eqnarray} &&\left. \frac{D \tilde{n}_{\pm }(p)}{Dt} \right|_{\rm coll,cs} = {-} c \: \sigma _{\rm T}\: \overline{s}_0(p) \: {N_{\rm ph}}\: \tilde{n}_{\pm }(p) + c \sigma _{\rm T}{N_{\rm ph}}\frac{4\pi }{\gamma } \int p_1 {d}\gamma _1 \: \overline{R}_{\rm e}(p,p_1) \: \tilde{n}_{\pm }(p_1), \end{eqnarray} \tag{ 25 }$$

where the electron redistribution function averaged over cosine of the electron scattering angle μ_e is

$$\begin{eqnarray} &&\overline{R}_{\rm e}(p,p_1)= \frac{1}{2} \int_{-1}^{1} R_{\rm e}({\bm p}_1 \rightarrow {\bm p}) \: {d}\mu _e = \frac{3}{16} \: \frac{2}{\lambda _{\rm C}^3{N_{\rm ph}}} \int_{x_{\star }(\gamma,\gamma _1)}^{\infty } \overline{R}_{\rm e}(\gamma,\gamma _1,x_1) \tilde{n}_{\rm ph}(x_1) \ {d}x_1, \end{eqnarray} \tag{ 26 }$$

where

$$\begin{equation} \overline{R}_{\rm e}(\gamma,\gamma _1,x_1) = \frac{1}{4\pi ^2} \: x_1 \int \frac{{d}^3 x}{x} {d}^2\omega _1 \: {d}^2\Omega _1 \: F \: \delta (\underline{p}_1 + \underline{x}_1 - \underline{p}- \underline{x}) \end{equation} \tag{ 27 }$$

and

$$\begin{equation} x_{\star }(\gamma,\gamma _1) = [\gamma - \gamma _1 + |p - p_1|]/2. \end{equation} \tag{ 28 }$$

The relation between the redistribution function and the extinction coefficient is

$$\begin{equation} \overline{s}_0(p) = \frac{4\pi }{\gamma } \int \overline{R}_{\rm e}(p_1,p) p_1 \, {d}\gamma _1. \end{equation} \tag{ 29 }$$

Not surprisingly, there turns out to be a relation between the quantities $\overline{R}_{\rm ph}$ and $\overline{R}_{\rm e}$ (proved in Appendix A), namely,

$$\begin{equation} pp_1 \overline{R}_{\rm e}(\gamma,\gamma _1,x_1) = xx_1 \overline{R}_{\rm ph}(x,x_1,\gamma _1), \end{equation} \tag{ 30 }$$

together with the energy conservation condition x + γ = x₁ + γ₁. Through Equations (30) and (18) we have a generally valid expression also for $\overline{R}_{\rm e}(\gamma,\gamma _1,x_1)$.

In the kinetic equation (Equation (6)) for the electron and positron densities n_±(p) the Compton terms can be obtained by multiplying Equation (25) by 8πλ⁻³_Cp³.

The electron RKE accounting for pair-production and annihilation processes can be written as (Nagirner & Loskutov 1999)

$$\begin{eqnarray} &&\underline{p}_{\,-}\cdot \underline{\nabla }\tilde{n}_{-}({\bm p}_{-})= \frac{r_{\rm e}^2}{4} \frac{2}{\lambda _{\rm C}^3} \int \frac{{d}^3 p_{+}}{\gamma _{+}} \frac{{d}^3 x_1}{x_1} \frac{{d}^3 x}{x}\:\delta (\underline{p}_{\,-} \,{+}\, \underline{p}_{\,+} \,{-}\, \underline{x}_1 \,{-}\, \underline{x}) F_{\gamma \gamma } [ \tilde{n}_{\rm ph}({\bm x}_1) \tilde{n}_{\rm ph}({\bm x}) - \tilde{n}_{-}({\bm p}_{-}) \tilde{n}_{+}({\bm p}_{+}) ], \end{eqnarray} \tag{ 31 }$$

where we used subscripts ∓ to explicitly show the momenta and the occupation number of electrons and positrons. Assuming homogeneity, we obtain

$$\begin{eqnarray} &&\left. \frac{D\tilde{n}_{-}({\bm p}_{-})}{Dt} \right|_{\rm coll, pp} = {-} c \: \sigma _{\rm T}\: \overline{s}_{\rm pa}({\bm p}_{-}) \: N_{\rm +}\: \tilde{n}_{-}({\bm p}_{-}) + c \: \sigma _{\rm T}\: {N_{\rm ph}^2}\: \frac{\lambda _{\rm C}^3}{2} \: j_{\rm pp}({\bm p}_{-}), \end{eqnarray} \tag{ 32 }$$

where the pair annihilation cross-section (in units of σ_T) is given by

$$\begin{eqnarray} &&\overline{s}_{\rm pa}({\bm p}_{-}) = \frac{3}{32\pi } \frac{2}{\lambda _{\rm C}^3N_{\rm +}} \frac{1}{\gamma _{-}} \int \frac{{d}^3 p_{+}}{\gamma _{+}} \frac{{d}^3 x_1}{x_1} \frac{{d}^3 x}{x} \tilde{n}_{+}({\bm p}_{+}) \: F_{\gamma \gamma } \delta (\underline{p}_{\,-}+ \underline{p}_{\,+}- \underline{x}_1 - \underline{x}) \end{eqnarray} \tag{ 33 }$$

and the pair production rate by

$$\begin{eqnarray} &&j_{\rm pp}({\bm p}_{-}) = \frac{3}{32\pi } \left(\frac{2}{\lambda _{\rm C}^3{N_{\rm ph}}}\right)^2 \frac{1}{\gamma _{-}} \int \frac{{d}^3 p_{+}}{\gamma _{+}} \frac{{d}^3 x_1}{x_1} \frac{{d}^3 x}{x} \tilde{n}_{\rm ph}({\bm x}) \: \tilde{n}_{\rm ph}({\bm x}_1) \: F_{\gamma \gamma }\: \delta (\underline{p}_{\,-}+ \underline{p}_{\,+}- \underline{x}_1 - \underline{x}). \end{eqnarray} \tag{ 34 }$$

The relativistically invariant reaction rate F_γγ is (Berestetskii et al. 1982)

$$\begin{equation} F_{\gamma \gamma }= \frac{\xi }{\xi _1} + \frac{\xi _1}{\xi } + 2 \; \left(\frac{1}{\xi } + \frac{1}{\xi _1}\right) - \left(\frac{1}{\xi } + \frac{1}{\xi _1}\right)^2, \end{equation} \tag{ 35 }$$

where $\xi = \underline{p}_{\,-}\cdot \underline{x}= \underline{p}_{\,+}\cdot \underline{x}_1$ and $\xi _1 = \underline{p}_{\,-}\cdot \underline{x}_1= \underline{p}_{\,+}\cdot \underline{x}$.

Assuming again isotropic particle distributions in frame E, we can write Equation (34) as

$$\begin{eqnarray} &&j_{\rm pp}(p_{-}) = 3\pi \left(\frac{2}{\lambda _{\rm C}^3{N_{\rm ph}}}\right)^2 \frac{1}{\gamma _{-}p_{-}} \int _{x^{(L)}}^{\infty } \tilde{n}_{\rm ph}(x) \: {d}x \int _{x_1^{(L)}}^{\infty } \tilde{n}_{\rm ph}(x_1) \: {d}x_1 \: R_{\rm \gamma \gamma }(\gamma _{-},x,x_1), \end{eqnarray} \tag{ 36 }$$

where we have defined

$$\begin{eqnarray} &&R_{\rm \gamma \gamma }(\gamma _{-},x,x_1) = \frac{1}{2} \frac{1}{(4\pi)^2} x x_1 p_{-}\int \frac{{d}^3 p_{+}}{\gamma _{+}} \: {d}^2\omega \: {d}^2\omega _1 F_{\gamma \gamma }\: \delta (\underline{p}_{\,-}+ \underline{p}_{\,+}- \underline{x}_1 - \underline{x}). \end{eqnarray} \tag{ 37 }$$

The cross-section becomes

$$\begin{equation} \overline{s}_{\rm pa}(p_{-}) = 4\pi \frac{2}{\lambda _{\rm C}^3N_{\rm +}} \int _0^{\infty } p_{+}^2 {d}p_{+}\: \sigma _{\rm pa}(\gamma _{+},\gamma _{-}) \: \tilde{n}_{+}(p_{+}), \end{equation} \tag{ 38 }$$

where

$$\begin{eqnarray} &&\sigma _{\rm pa}(\gamma _{+},\gamma _{-})= \frac{3}{8} \frac{1}{(4\pi)^2} \frac{1}{\gamma _{-}\gamma _{+}} \int \frac{{d}^3 x_1}{x_1} \frac{{d}^3 x}{x} \, {d}^2\Omega _{+}\: F_{\gamma \gamma } \delta (\underline{p}_{\,-}+ \underline{p}_{\,+}- \underline{x}_1 - \underline{x}). \end{eqnarray} \tag{ 39 }$$

The treatment of positrons is identical if we switch the subscripts − and + in Equations (31)–(39).

The photon kinetic equation accounting for pair production/annihilation processes is

$$\begin{eqnarray} &&\underline{x}\cdot \underline{\nabla }\tilde{n}_{\rm ph}({\bm x}) = \frac{r_{\rm e}^2}{2} \frac{2}{\lambda _{\rm C}^3} \int\! \frac{{d}^3 x_1}{x_1} \frac{{d}^3 p_{-}}{\gamma _{-}} \frac{{d}^3 p_{+}}{\gamma _{+}} \: \delta (\underline{p}_{\,-}+ \underline{p}_{\,+}- \underline{x}_1 - \underline{x}) F_{\gamma \gamma }[ \tilde{n}_{-}({\bm p}_{-}) \tilde{n}_{+}({\bm p}_{+}) - \tilde{n}_{\rm ph}({\bm x}_1) \tilde{n}_{\rm ph}({\bm x}) ]. \end{eqnarray} \tag{ 40 }$$

Neglecting the spatial derivatives in the left-hand side, this becomes

$$\begin{equation} \left.\frac{D \tilde{n}_{\rm ph}({\bm x})}{Dt} \right|_{\rm coll, pp} = - c \: \sigma _{\rm T}\: \overline{s}_{\rm pp}({\bm x}) \: {N_{\rm ph}}\: \tilde{n}_{\rm ph}({\bm x}) + c \: \sigma _{\rm T}\: {N_{\rm -}}\: N_{\rm +}\: \frac{\lambda _{\rm C}^3}{2}\ j_{\rm pa}({\bm x}), \end{equation} \tag{ 41 }$$

where the pair-production cross-section is

$$\begin{eqnarray} &&\overline{s}_{\rm pp}({\bm x}) = \frac{3}{16\pi } \frac{2}{\lambda _{\rm C}^3{N_{\rm ph}}} \frac{1}{x} \int \frac{{d}^3 x_1}{x_1} \frac{{d}^3 p_{-}}{\gamma _{-}} \frac{{d}^3 p_{+}}{\gamma _{+}} \tilde{n}_{\rm ph}({\bm x}_1) F_{\gamma \gamma }\: \delta (\underline{p}_{\,-}+ \underline{p}_{\,+}- \underline{x}_1 - \underline{x}) \end{eqnarray} \tag{ 42 }$$

and the emissivity due to pair annihilation

$$\begin{eqnarray} &&j_{\rm pa}({\bm x}) = \frac{3}{16\pi } \left(\frac{2}{\lambda _{\rm C}^3}\right)^2 \frac{1}{{N_{\rm -}}N_{\rm +}} \frac{1}{x} \int \frac{{d}^3 x_1}{x_1} \frac{{d}^3 p_{-}}{\gamma _{-}} \frac{{d}^3 p_{+}}{\gamma _{+}} \tilde{n}_{-}({\bm p}_{-}) \: \tilde{n}_{+}({\bm p}_{+}) \: F_{\gamma \gamma }\: \delta (\underline{p}_{\,-}+ \underline{p}_{\,+}- \underline{x}_1 - \underline{x}). \end{eqnarray} \tag{ 43 }$$

Note that unlike the electron equation, the photon equation is nonlinear owing to the fact that the cross-section (42) depends explicitly on the photon distribution.

Under the isotropy assumption Equations (42) and (43) in frame E become

$$\begin{eqnarray} &&j_{\rm pa}(x) = 6\pi \left(\frac{2}{\lambda _{\rm C}^3}\right)^2 \frac{1}{{N_{\rm -}}N_{\rm +}} \frac{1}{x^2} \int _{\gamma _{+}^{(L)}}^{\infty } \tilde{n}_{+}(p_{+}) \: {d}\gamma _{+} \int _{\gamma _{-}^{(L)}}^{\infty } \tilde{n}_{-}(p_{-}) \: {d}\gamma _{-}R_{\rm \gamma \gamma }(\gamma _{-},x,x_1), \end{eqnarray} \tag{ 44 }$$

where we have to substitute x₁ = γ₋ + γ₊ − x from the energy conservation condition, and

$$\begin{equation} \overline{s}_{\rm pp}(x) = 4\pi \frac{2}{\lambda _{\rm C}^3{N_{\rm ph}}} \int _{1/x}^{\infty } x_1^2 {d}x_1 \: \sigma _{\rm pp}(x,x_1) \: \tilde{n}_{\rm ph}(x_1), \end{equation} \tag{ 45 }$$

where

$$\begin{eqnarray} &&\sigma _{\rm pp}(x,x_1) = \frac{3}{4} \frac{1}{(4\pi)^2} \frac{1}{x x_1} \int \frac{{d}^3 p_{-}}{\gamma _{-}} \frac{{d}^3 p_{+}}{\gamma _{+}} \: {d}^2 \omega _1 \: F_{\gamma \gamma } \delta (\underline{p}_{\,-}+ \underline{p}_{\,+}- \underline{x}_1 - \underline{x}). \end{eqnarray} \tag{ 46 }$$

Explicit expressions for the rate R_γγ(γ₋, x, x₁) (derived by Svensson 1982, see also Boettcher & Schlickeiser 1997 and Nagirner & Loskutov 1999) and the cross-sections σ_pa(γ₊, γ₋), σ_pp(x, x₁) as well as the lower integration limits in Equations (36) and (44) are given in Appendix D.

The pair-production terms in Equations (5) and (6) take the form

$$\begin{eqnarray} \dot{n}_{\rm ph,pp} (x) & = & -\! c \alpha _{\rm pp}(x) {n_{\rm ph}}(x) + \epsilon _{\rm pa}(x), \end{eqnarray} \tag{ 47 }$$

$$\begin{eqnarray} \dot{n}_{\rm \pm,pp} (p_{\pm }) & = & -\!c \alpha _{\rm pa}(p_{\pm }) {n_{\pm }}(p_{\pm }) + \epsilon _{\rm pp}(p_{\pm }). \end{eqnarray} \tag{ 48 }$$

By comparing with Equations (32) and (41) we find the absorption coefficients and emissivities to be

$$\begin{eqnarray} \alpha _{\rm pp}(x) & =& \sigma _{\rm T}\: \overline{s}_{\rm pp}(x) \: {N_{\rm ph}},\ \epsilon _{\rm pa}(x) = 4\pi \: c \: \sigma _{\rm T}\: {N_{\rm -}}\: N_{\rm +}\: x^3 \: j_{\rm pa}(x), \end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray} \alpha _{\rm pa}(p_{\pm }) & =& \sigma _{\rm T}\: \overline{s}_{\rm pa}(p_{\pm }) \: N_{\mp },\ \epsilon _{\rm pp}(p_{\pm }) = 4\pi \: c \: \sigma _{\rm T}\: {N_{\rm ph}^2}\: p_{\pm }^3 \: j_{\rm pp}(p_{\pm }). \end{eqnarray} \tag{ 50 }$$

The kinetic equations describing synchrotron radiation need to be written in frame E, where we assume there is only tangled magnetic field (and no electric field). Using the Einstein coefficients and the cross-sections describing synchrotron emission and absorption (Ghisellini & Svensson 1991), we get the collision terms for these processes in the electron/positron and photon equations (see also Ochelkov et al. 1979):

$$\begin{eqnarray} \frac{D}{Dt} [\gamma p \tilde{n}_{\pm }(p)] \bigg|_{\rm coll, syn} &=& \int _0^{\infty } {d}x \int _{\gamma }^{\infty } {d}\gamma _1 \gamma _1 p_1 \frac{P(x,\gamma _1)}{x} \delta (\gamma _1 - \gamma - x) \{ \tilde{n}_{\pm }(p_1) [1 + \tilde{n}_{\rm ph}(x)] - \tilde{n}_{\pm }(p) \tilde{n}_{\rm ph}(x) \}\nonumber\\ && - \int _0^{\infty } {d}x \int _{1}^{\gamma } {d}\gamma _1 \gamma p \frac{P(x,\gamma)}{x} \delta (\gamma - \gamma _1 - x) \{ \tilde{n}_{\pm }(p) [1 + \tilde{n}_{\rm ph}(x)] - \tilde{n}_{\pm }(p_1) \tilde{n}_{\rm ph}(x) \}, \end{eqnarray} \tag{ 51 }$$

$$\begin{eqnarray} &&\frac{D}{Dt} [ x^2 \tilde{n}_{\rm ph}(x) ]\bigg|_{\rm coll, syn} = \int _{1}^{\infty } {d}\gamma \int _{1}^{\gamma } {d}\gamma _1 \: \gamma p \: \frac{P(x,\gamma) }{x} \delta (\gamma - \gamma _1 - x) \lbrace \tilde{n}_{\rm e}(p)[1+\tilde{n}_{\rm ph}(x)] - \tilde{n}_{\rm e}(p_1) \tilde{n}_{\rm ph}(x) \rbrace. \end{eqnarray} \tag{ 52 }$$

Here, P(x, γ) is the angle-integrated cyclosynchrotron spectrum of a single electron, normalized to the electron cooling rate:

$$\begin{equation} \int _{0}^{\infty } P(x, \gamma) \, {d}x = - \dot{\gamma }_{\rm s}= \frac{4 }{3} \frac{\sigma _{\rm T}U_{\rm B}}{m_{\rm e}c} p^2, \end{equation} \tag{ 53 }$$

where U_B = B²/(8π) is the magnetic energy density. One can readily verify that Equations (51) conserve the total number of electrons and positrons, and that the total energy is conserved by Equations (51) and (52).

Under the physical conditions that we are interested in, the average energy (or momentum) of an emitted or absorbed photon is much lower than the energy (momentum) of the electron taking part in the process. The standard way is therefore to treat synchrotron processes as continuous cooling or heating for electrons and as an emission or absorption process for photons.

We write the photon terms in the form

$$\begin{equation} \left.\frac{D \tilde{n}_{\rm ph}(x)}{Dt} \right| _{\rm coll, syn} = -c \alpha _{\rm s}(x) \tilde{n}_{\rm ph}(x) + \frac{\lambda _{\rm C}^3}{8\pi } \frac{\epsilon _{\rm s}(x)}{x^3}, \end{equation} \tag{ 54 }$$

where α_s and _s are cyclosynchrotron absorption and emission coefficients, respectively. In the kinetic Equation (5) for the photon density n_ph(x) the corresponding term can be obtained by multiplying Equation (54) by 8πλ⁻³_Cx³:

$$\begin{equation} \dot{n}_{\rm ph,syn} (x) = -c \alpha _{\rm s}(x) {n_{\rm ph}}(x) + \epsilon _{\rm s}(x). \end{equation} \tag{ 55 }$$

The emissivity _s gives the number of photons emitted per logarithmic dimensionless energy interval dln x, per unit volume and time and can be identified by comparing the corresponding terms in Equations (52) and (54):

$$\begin{equation} \epsilon _{\rm s}(x) = \frac{8\pi }{\lambda _{\rm C}^3} \int P(x,\gamma) \; p^2 \tilde{n}_{\rm e}(p) \; {d}p = \int P(x,\gamma) \; {n_{\rm e}}(p) \; {d}\ln {p}. \end{equation} \tag{ 56 }$$

Similarly, by comparing the terms proportional to $\tilde{n}_{\rm ph}$ we identify the absorption coefficient (e.g., Rybicki & Lightman 1979):

$$\begin{eqnarray} &&\alpha _{\rm s}(x) = \frac{1}{4\pi c x^3} \int {d}^3p [ \tilde{n}_{\rm e}(p_1) - \tilde{n}_{\rm e}(p)] P(x, \gamma) = {-} \frac{1}{c \ x^2} \int \frac{{d}\tilde{n}_{\rm e}(p)}{{d}p} P(x,\gamma) \: \gamma p \ {d}p, \end{eqnarray} \tag{ 57 }$$

where $p_1=\sqrt{(\gamma - x)^2-1}$ is the electron momentum corresponding to energy γ₁ = γ − x and the second expression is obtained by expansion to the first order in x ≪ γ.

In terms of the electron number density n_e(p) the absorption coefficient takes the form:

$$\begin{equation} \alpha _{\rm s}(x) = \frac{\lambda _{\rm C}^3}{8\pi c} \frac{1}{x^2} \int \frac{\gamma P(x,\gamma)}{p^2} \left[ 3 {n_{\rm e}}(p) - \frac{{d}{n_{\rm e}}(p)}{{d}\ln {p}} \right] \: {d}\ln {p}. \end{equation} \tag{ 58 }$$

The synchrotron processes for electrons can be treated as continuous using the Fokker-Planck equation. It can be obtained from Equation (51) employing the delta-function to take the integral over γ₁ and expanding γ₁p₁ P(x, γ₁) and n_±(p₁) near p to the second order in the small "parameter" x. Collecting the terms and finally integrating over the photon energy x we obtain

$$\begin{eqnarray} &&\frac{\partial }{\partial t} [\gamma p \tilde{n}_{\pm }(p)] = {-}\frac{\partial }{\partial \gamma } \left[ \dot{\gamma }_{\rm s}\ \gamma p\ \tilde{n}_{\pm }(p) - H(p)\ \gamma p\ \frac{\partial \tilde{n}_{\pm }(p)}{\partial \gamma } \right]+ \frac{1}{2} \frac{\partial ^2}{\partial \gamma ^2} [ H_0(p) \gamma p \tilde{n}_{\pm }(p) ], \end{eqnarray} \tag{ 59 }$$

where

$$\begin{eqnarray} &&H(p) = \int P(x,\gamma) \ \tilde{n}_{\rm ph}(x) \; x \ {d}x = \frac{\lambda _{\rm C}^3}{8\pi } \int \frac{P(x,\gamma) }{x} \ {n_{\rm ph}}(x) \ {d}\ln {x},\quad H_0(p) = \int P(x,\gamma) \; x \ {d}x. \end{eqnarray} \tag{ 60 }$$

To get the total electron energy gain/loss rate, one has to multiply Equation (59) by 8πλ⁻³_Cγ dγ and integrate. Multiplying Equation (54) by 8πλ⁻³_Cx³dx and integrating gives the corresponding rate for photons. Using Equations (53), (56), (57), and (60), we can verify that energy conservation is maintained when switching from Equation (51) to the continuous approximation (59).

Note that the last term on the right-hand side of Equation (59) is missing in similar equations derived previously (McCray 1969; Ghisellini et al. 1988). It corresponds to the diffusion due to spontaneous emission, but does not contribute to the electron cooling/heating. However, in most cases we expect its contribution to be negligible compared to the other terms. It is of the order x/γ smaller than the cooling term with $|\dot{\gamma }_{\rm s}|$ and, when electrons are mildly relativistic, self-absorption becomes important, $\tilde{n}_{\rm ph}\gg 1$ and the term containing H dominates. Therefore, we neglect the term with H₀ in our simulations. Thus, for the distributions n_±(p), Equation (59) takes the form

$$\begin{equation} \dot{n}_{\rm \pm,syn} (p) = -\frac{\partial }{\partial \ln {p}} \left[ A_{\rm e,syn}(p) {n_{\pm }}(p) - B_{\rm e,syn}(p) \frac{\partial {n_{\pm }}(p)}{\partial \ln {p}} \right], \end{equation} \tag{ 61 }$$

where

$$\begin{eqnarray} A_{\rm e,syn}(p) &=& \left(\dot{\gamma }_{\rm s}+ 3 \frac{\gamma }{p^2}H(p) \right) \;\frac{\gamma }{p^2},\quad B_{\rm e,syn}(p) = H(p)\frac{\gamma ^2}{p^4}. \end{eqnarray} \tag{ 62 }$$

It is worth mentioning here that other emission/absorption processes, e.g., bremsstrahlung, can be implemented analogously to the synchrotron radiation, once the emissivity function of a single electron P(x, γ) (which now may depend on the particle distribution) is specified.

The RKE accounting for electron (positron) evolution due to Coulomb scatterings is

$$\begin{eqnarray} &&\underline{p}\cdot \underline{\nabla }\tilde{n}_{\pm }({\bm p}) = r_{\rm e}^2 \frac{2}{\lambda _{\rm C}^3} \int \frac{{d}^3 p_1}{\gamma _1} \frac{{d}^3 p^{\prime }_1}{\gamma ^{\prime }_1} \frac{{d}^3 p^{\prime }}{\gamma ^{\prime }} \: \delta (\underline{p}_1 + \underline{p}- \underline{p}^{\prime }_{1}- \underline{p}^{\prime }) F_{\rm Coul}[ \tilde{n}_{\rm e}({\bm p}^{\prime }_1) \tilde{n}_{\pm }({\bm p}^{\prime }) - \tilde{n}_{\rm e}({\bm p}_1) \tilde{n}_{\pm }({\bm p}) ]. \end{eqnarray} \tag{ 63 }$$

The invariant reaction rate for Møller scattering (i.e., e⁻e⁻ and e⁺e⁺) is given by (Berestetskii et al. 1982)

$$\begin{equation} F_{\rm Coul}= \left(\frac{\xi _1}{\xi - 1} + \frac{\xi }{\xi _1 - 1} + 1 \right)^2 +\frac{1 - 4\xi \xi _1}{(\xi - 1)(\xi _1 - 1)} + 4 \end{equation} \tag{ 64 }$$

and the scalar products of particles' four-momenta are defined as $\xi = \underline{p}\cdot \underline{p}^{\prime }$ and $\xi _1 = \underline{p}_1\cdot \underline{p}^{\prime }$. As discussed by Baring (1987) and Coppi & Blandford (1990), the corresponding rates for Bhabha e^±e^∓ scattering are nearly the same in the small-angle scattering approximation, we therefore do not distinguish between electrons and positrons in these equations.

Although the Coulomb process is collisional in nature, it is customary to treat it in the Fokker–Planck framework, i.e., as a continuous diffusive energy exchange mechanism. This is due to the well-known divergent nature of the Coulomb cross-section for small-angle scatterings with negligible energy exchange per event, while in the parameter regimes we are interested in, a large number of such scatterings dominate the energy gain or loss rate of a particle over a much smaller number of large-angle scatterings. In frame E, where the particle distributions are approximately homogeneous and isotropic, we can therefore write the Coulomb terms in the form of the Fokker–Planck equation in (scalar) momentum space

$$\begin{equation} \dot{n}_{\rm \pm,Coul} (p) = -\frac{\partial }{\partial \ln {p}} \left[ A_{\rm e,Coul}(p) {n_{\pm }}(p) - B_{\rm e,Coul}(p) \frac{\partial {n_{\pm }}(p)}{\partial \ln {p}} \right] \end{equation} \tag{ 65 }$$

with coefficients given by

$$\begin{eqnarray} &&A_{\rm e,Coul}(p) = \frac{\dot{\gamma }_{\rm Coul}\gamma }{p^2} - \frac{\partial }{\partial \gamma } \left(\frac{1}{2} \frac{\gamma D_{\rm Coul}}{p^2} \right),\quad B_{\rm e,Coul}(p) = \frac{1}{2} \frac{\gamma ^2 D_{\rm Coul}}{p^4}. \end{eqnarray} \tag{ 66 }$$

The energy exchange rate and the diffusion coefficient can be obtained by calculating the first and second moments of Equation (63) keeping only small-angle scatterings and are expressed as integrals over the particle distributions:

$$\begin{eqnarray} &&\dot{\gamma }_{\rm Coul}= \int a(\gamma,\gamma _1) {n_{\rm e}}(p_1) \, {d}\ln {p_1},\quad D_{\rm Coul}(p) = \int d(\gamma,\gamma _1) {n_{\rm e}}(p_1) \, {d}\ln {p_1}. \end{eqnarray} \tag{ 67 }$$

The rates a(γ, γ₁) and d(γ, γ₁) have been calculated by Nayakshin & Melia (1998) and are given in Appendix F.

The matrix $M_{i,i^{\prime }}$ of the linear system can be decomposed into two parts arising from the differential and integral terms in Equations (7) and (8). The differential part contributes a tridiagonal matrix, the form of the equation (e.g., for electrons), giving rise to it, is

$$\begin{equation} \frac{n_{i}^{k+1} - n_{i}^{k}}{\Delta t_k} = -\frac{1}{\Delta _p} \left[ F^{k+1/2}_{i+1/2} - F^{k+1/2}_{i-1/2} \right], \end{equation} \tag{ 71 }$$

where the momentum space flux is given by

$$\begin{equation} F^{k+1/2}_{i+1/2} = A^{k+1/2}_{i+1/2} \; n_{i+1/2}^{k+1/2} - B^{k+1/2}_{i+1/2} \; \frac{n_{i+1}^{k+1/2} - n_{i}^{k+1/2}}{\Delta _p}. \end{equation} \tag{ 72 }$$

The distribution function between time grid points is defined according to the Crank-Nicolson scheme as (omitting the momentum index)

$$\begin{equation} n^{k+1/2} = \frac{1}{2} (n^{k+1} + n^{k}). \end{equation} \tag{ 73 }$$

We also have to somehow define the distribution function between momentum grid points. Following Chang & Cooper (1970) we introduce a parameter δ_i so that (now omitting the time index)

$$\begin{equation} n_{i+1/2} = (1-\delta _i) n_{i+1} + \delta _i n_{i}, \quad \delta _i \in [0, 1]. \end{equation} \tag{ 74 }$$

The basic idea of the Chang and Cooper scheme is to employ this parameter to ensure that the differencing scheme converges to the correct equilibrium solution independently of the size of the grid step Δ_p. Assuming that the momentum space flux through the boundaries vanishes, the equilibrium solution tells us that it must vanish everywhere, i.e., F = 0. From Equations (72) and (74) we then have

$$\begin{equation} \frac{n_{i+1}}{n_{i}} = \frac{\delta _i A_{i+1/2} \Delta _p+ B_{i+1/2}}{B_{i+1/2} - (1-\delta _i) A_{i+1/2} \Delta _p}, \end{equation} \tag{ 75 }$$

while the exact solution gives (Chang & Cooper 1970)

$$\begin{equation} \frac{n_{i+1}}{n_{i}} = \exp \left[\frac{A_{i+1/2}}{B_{i+1/2}} \Delta _p\right]. \end{equation} \tag{ 76 }$$

We can see that using either centered differencing (δ = 1/2) or forward differencing (δ = 0), Equations (75) and (76) agree only to the first order in AΔ_p/B. To make the correspondence exact, one has to equate the two equations and solve for δ_i, to obtain

$$\begin{equation} \delta _i = \frac{1}{w_i} - \frac{1}{\exp (w_i) - 1}, \qquad w_i = -\frac{A_{i+1/2}}{B_{i+1/2}} \Delta _p. \end{equation} \tag{ 77 }$$

Aside from converging to the correct equilibrium solution, such choice of δ_i also guarantees positive spectra, as shown by Chang & Cooper (1970). Although this method applies to purely differential equations, we can still use it in our integro-differential equations to ensure that the differential part tends toward its own correct equilibrium solution, which would also be the correct solution for the full equation in the region where the differential terms happen to dominate.

Accurate numerical treatment of Compton scattering over a wide range of energies is not straightforward. This is caused by the well-known fact that at different energies the process takes place in different regimes. If the energy of a photon in electron rest frame is much smaller than the electron rest energy, the process takes place in the Thomson regime and the electron loses a very small amount of its energy in one scattering. Correspondingly, there is a sharp peak in the electron redistribution function R_e near p = p₁. We cannot therefore numerically resolve R_e on our finite grid and have to treat the energy loss process as continuous. On the other hand, for scattering in the Klein–Nishina regime the electron can lose a significant amount of its energy in one scattering. Wishing to include both regimes, we need a way to switch from the continuous approximation (implying a differential equation) to direct calculation of scattering through the integral terms. Similar treatment is required for photons, although the continuous approximation is only needed in the regime where the photon energy is much lower than the electron rest energy and the electron is nonrelativistic.

3.2.1. Scattering of Electrons: Separation of Regimes

Let us first look at the electron redistribution function (Equation (26)). We wish to know what is the lowest incoming photon energy x^±_⋆(p₁) that can cause a shift in electron momentum p₁ by |Δln p|. This energy is related to the lower limit (Equation (28)) of the integral in Equation (26). If the shift is small enough, we can write

$$\begin{eqnarray} &&x_{\star }^{\pm }(p_1) \approx x_{\star }(\gamma,\gamma _1)= \frac{1}{2} \left(\pm |\Delta \gamma | + |\Delta p| \right) \approx \frac{1}{2} p_1 |\Delta \ln {p_1}| \left(1 \pm \frac{p_1}{\gamma _1}\right), \end{eqnarray} \tag{ 78 }$$

where we have used p dp = γ dγ. The plus sign applies when the electron gains energy and the minus when it loses it. We see that for high-energy electrons, the minimum energy of photons for which we can resolve up- or downscattering is vastly different. However, since the upscattering (energy increase) of relativistic electrons is extremely inefficient, we concern ourselves only with being able to resolve their downscattering (i.e., cooling) and so use the minus sign in Equation (78). Choosing |Δln p₁| comparable to our grid step (we use somewhat arbitrarily 3Δ_p) in the electron equation, we then state that the scattering of electrons on photons with x₁ < x⁻_⋆(p₁) cannot be resolved.

We now split the redistribution function into two parts according to whether we can or cannot resolve it on our grid

$$\begin{equation} \overline{R}_{\rm e}(p,p_1) = \overline{R}_{\rm e}^{<}(p,p_1) + \overline{R}_{\rm e}^{>}(p,p_1), \end{equation} \tag{ 79 }$$

where for the first term the integral in Equation (26) is taken over x₁ < x⁻_⋆(p₁), and the second term is defined by integrating over the remaining x₁. To totally isolate scatterings that undergo on photons with energies below and above x⁻_⋆, we have to write the extinction coefficient as an analogous sum, $\overline{s}_0(p) = \overline{s}_0^{<}(p) + \overline{s}_0^{>}(p)$, where

$$\begin{equation} \overline{s}^{\lessgtr }_0(p) = \frac{4\pi }{\gamma } \int \overline{R}_{\rm e}^{\lessgtr }(p_1,p)\ p_1 \, {d}\gamma _1, \end{equation} \tag{ 80 }$$

in accordance with Equation (29). For the terms containing $\overline{R}_{\rm e}^{>}$ and $\overline{s}_0^{>}$ in the electron equation, we compute the integrals through the discrete sums, but the terms containing $\overline{R}_{\rm e}^{<}$ and $\overline{s}_0^{<}$ have to be accounted for by continuous energy exchange terms in the equation. Since we also want to treat thermalization by Compton scattering, these terms have to contain a second-order derivative of the electron distribution (a diffusive term). Therefore, we take the standard form of the Fokker–Planck equation

$$\begin{equation} \dot{N}_{\rm \pm,diff,cs}(\gamma) = -\frac{\partial }{\partial \gamma } \left\lbrace \dot{\gamma }_{\rm c} N_{\pm }(\gamma) - \frac{1}{2} \frac{\partial }{\partial \gamma } \left[ D_{\rm e}(\gamma) N_{\pm }(\gamma) \right] \right\rbrace, \end{equation} \tag{ 81 }$$

while the exact equation for the (<) terms comes from Equation (25), written here for N_±(γ),

$$\begin{eqnarray} &&\dot{N}_{\rm \pm,coll,cs}^{<} (\gamma) = {-}c \sigma _{\rm T} \overline{s}_0^{<}(p)\; {N_{\rm ph}} N_{\pm }(\gamma) + 4\pi c \sigma _{\rm T}{N_{\rm ph}} p \int \frac{{d}\gamma _1}{\gamma _1} \: \overline{R}_{\rm e}^{<}(p,p_1) \: N_{\pm }(\gamma _1). \end{eqnarray} \tag{ 82 }$$

In order to make a physically sensible correspondence between these two representations, we demand that the first three moments of Equations (81) and (82) were identical. Substituting Equation (80) into (82), we find

$$\begin{eqnarray} \int \dot{N}_{\rm \pm,coll,cs}^{<}(\gamma) \gamma ^i {d}\gamma &=& 4\pi c \sigma _{\rm T}{N_{\rm ph}} \int {d}\gamma \int {d}\gamma _1 \gamma ^i\! \left\lbrace -\frac{p_1}{\gamma } \overline{R}_{\rm e}^{<}(p_1,p) N_{\pm }(\gamma) + \frac{p}{\gamma _1} \overline{R}_{\rm e}^{<}(p,p_1) N_{\pm }(\gamma _1) \right\rbrace \nonumber \\ & =& 4\pi c \sigma _{\rm T}{N_{\rm ph}} \int {d}\gamma \int {d}\gamma _1 \frac{p_1}{\gamma } \left(\gamma _1^i - \gamma ^i \right) \overline{R}_{\rm e}^{<}(p_1,p) N_{\pm }(\gamma) = c \sigma _{\rm T}{N_{\rm ph}} \int {d}\gamma \; \overline{(\gamma _1^i - \gamma ^i)} \overline{s}^{<}_0(p) N_{\pm }(\gamma), \end{eqnarray} \tag{ 83 }$$

where similarly to the moments of the photon redistribution function (NP94), we defined the moments of the electron redistribution function

$$\begin{equation} \overline{\gamma _1^i} \overline{s}^{<}_0(p) \equiv \frac{4\pi }{\gamma } \int p_1 \gamma _1^i \, {d}\gamma _1 \; \overline{R}_{\rm e}^{<}(p_1,p). \end{equation} \tag{ 84 }$$

The zeroth moment (giving zero in the right-hand side of Equation (83)) is just a statement of particle number conservation, while the first moment gives the total rate at which the electrons gain (or lose) energy. The moments defined by Equation (84) can be calculated analytically using the exact Klein–Nishina scattering cross-section. For photons this was shown by NP94, while the extension of these calculations to the electrons is given in Appendix C.

The moments of the continuous approximation (Equation (81)) are

$$\begin{eqnarray} \int \dot{N}_{\rm \pm,diff,cs} (\gamma) \, {d}\gamma & =& 0, \end{eqnarray} \tag{ 85 }$$

$$\begin{eqnarray} \int \dot{N}_{\rm \pm,diff,cs} (\gamma) \gamma \, {d}\gamma & =& \int \dot{\gamma }_{\rm c} N_{\pm }(\gamma) \, {d}\gamma, \end{eqnarray} \tag{ 86 }$$

$$\begin{eqnarray} \int \dot{N}_{\rm \pm,diff,cs} (\gamma) \gamma ^2 \, {d}\gamma & =& \int \left[ 2\gamma \dot{\gamma }_{\rm c}+ D_{\rm e}(\gamma) \right] N_{\pm }(\gamma) \, {d}\gamma.\quad \end{eqnarray} \tag{ 87 }$$

Here, we have assumed that the distribution function N_±(γ) vanishes at the boundaries of integration. Exact correspondence with Equation (83) can be made if we identify

$$\begin{eqnarray} &&\dot{\gamma }_{\rm c}= c \sigma _{\rm T}{N_{\rm ph}} \overline{(\gamma _1 - \gamma)} \overline{s}^{<}_0(p),\quad D_{\rm e}(\gamma) = c \sigma _{\rm T}{N_{\rm ph}} \overline{(\gamma _1 - \gamma)^2} \overline{s}^{<}_0(p), \end{eqnarray} \tag{ 88 }$$

while for the zeroth moment the correspondence is automatic. These moments can be computed using Equations (C11) and (C12). Finally, we write Equation (81) through n_±(p) and in the form that can be included in the Chang and Cooper differencing scheme together with other terms

$$\begin{equation} \dot{n}_{\rm \pm,diff,cs} (p) = -\frac{\partial }{\partial \ln {p}} \left[ A_{\rm e,cs}(p) {n_{\pm }}(p) - B_{\rm e,cs}(p) \frac{\partial {n_{\pm }}(p)}{\partial \ln {p}} \right], \end{equation} \tag{ 89 }$$

where

$$\begin{eqnarray} &&A_{\rm e,cs}(p) = \frac{\dot{\gamma }_{\rm c}\gamma }{p^2} - \frac{\partial }{\partial \gamma } \left(\frac{1}{2} \frac{\gamma D_{\rm e}(\gamma)}{p^2} \right),\quad B_{\rm e,cs}(p) = \frac{1}{2} \frac{\gamma ^2 D_{\rm e}(\gamma)}{p^4}. \end{eqnarray} \tag{ 90 }$$

3.2.2. Scattering of Photons and Three-Bin Approximation

Insufficient resolution of numerical calculations can become an issue also for the scattering of photons if the electron energies are low enough. A photon will then exchange very little energy with an electron upon scattering and the redistribution function is strongly peaked near x = x₁. To overcome this we propose the following approach. We separate scatterings that take place within some narrow interval around the energy of the incoming photon from those invoking photon energy outside this interval. We then approximate the scatterings taking place within the central interval by a continuous process and account for this by differential terms calculated through the exact moments of the redistribution function.

To keep the correspondence to the electron equation, we rewrite the photon evolution Equation (16) in terms of N_ph(x):

$$\begin{eqnarray} \dot{N}_{\rm ph,coll,cs} (x) &=& 4\pi c \sigma _{\rm T}{N_{\rm e}}\bigg\lbrace \!\!\int _{\notin } {d}{x_1} \bigg[ {N_{\rm ph}}(x_1) \frac{x}{x_1} \overline{R}_{\rm ph}(x,x_1) - {N_{\rm ph}}(x) \frac{x_1}{x} \overline{R}_{\rm ph}(x_1,x) \bigg]\nonumber\\ && +\, \int _\in \: {d}{x_1} \bigg[ {N_{\rm ph}}(x_1) \frac{x}{x_1} \overline{R}_{\rm ph}(x,x_1) - {N_{\rm ph}}(x) \frac{x_1}{x} \overline{R}_{\rm ph}(x_1,x) \bigg] \bigg\rbrace, \end{eqnarray} \tag{ 91 }$$

where the extinction coefficient is expressed explicitly through $\overline{R}_{\rm ph}$ using Equation (20). Here, ∈ stands for the interval $[x {\rm e}^{- \delta _x},x {\rm e}^{\delta _x}]$ and ∉ means integration from 0 to ∞ excluding that interval. The width of the central region (2δ_x in log units) is somewhat arbitrary, but should include at least a couple of bins, with our choice being three, i.e., $\delta _x=\frac{3}{2} \Delta _x$.

For the second integral in Equation (91) we wish to write a continuous approximation similar to Equation (81)

$$\begin{equation} \dot{N}_{\rm ph,diff,cs} (x) = -\frac{\partial }{\partial x} \left\lbrace \dot{x}_{\rm c} {N_{\rm ph}}(x) - \frac{1}{2} \frac{\partial }{\partial x} [ D_{\rm ph}(x) {N_{\rm ph}}(x) ] \right\rbrace. \end{equation} \tag{ 92 }$$

Similarly to what was done for electrons, the coefficients in Equation (92) are determined from the requirement that the first three moments of the differential and integral equations coincide. The moments of the "central" part of Equation (91) (denoted by ∈) are

$$\begin{eqnarray} &&\int _0^{\infty } \dot{N}_{\rm ph,coll,cs}^{\in } (x) x^i\, {d}x = 4\pi c \sigma _{\rm T}{N_{\rm e}}\int _{0}^{\infty } {d}{x}\int _{\in } {d}{x_1} \: \left(x_1^i - x^i \right) \frac{x_1}{x} \overline{R}_{\rm ph}(x_1,x) {N_{\rm ph}}(x), \end{eqnarray} \tag{ 93 }$$

where the integration limits for x and x₁ in the first term were switched, because for constant δ_x the area on the (x, x₁) plane is the same. The moments of the differential equation are similar to what were obtained for electrons

$$\begin{eqnarray} \int _0^{\infty } \dot{N}_{\rm ph,diff,cs}^{\in } (x) \, {d}x & =& 0, \end{eqnarray} \tag{ 94 }$$

$$\begin{eqnarray} \int _0^{\infty } \dot{N}_{\rm ph,diff,cs}^{\in } (x) x \, {d}x & =& \int _0^{\infty } \dot{x}_{\rm c} {N_{\rm ph}}(x) \, {d}x, \end{eqnarray} \tag{ 95 }$$

$$\begin{eqnarray} \int _0^{\infty } \dot{N}_{\rm ph,diff,cs}^{\in } (x) x^2 \, {d}x & = & \int _0^{\infty } \left[ 2x \dot{x}_{\rm c}+ D_{\rm ph}(x) \right] {N_{\rm ph}}(x) \, {d}x. \end{eqnarray} \tag{ 96 }$$

Equations (93)–(96) give identical expressions for the first three moments of the "central" part of the equation if we identify

$$\begin{eqnarray} &&\dot{x}_{\rm c}= 4\pi c \sigma _{\rm T}{N_{\rm e}}\int _{\in } {d}{x_1} \: \left(x_1 - x \right) \frac{x_1}{x} \overline{R}_{\rm ph}(x_1,x),\quad D_{\rm ph}(x) = 4\pi c \sigma _{\rm T}{N_{\rm e}}\int _{\in } {d}{x_1} \: \left(x_1 - x \right)^2 \frac{x_1}{x} \overline{R}_{\rm ph}(x_1,x). \end{eqnarray} \tag{ 97 }$$

The 0th moment is identically zero for both Equations (93) and (94), implying particle conservation.

The integrals in Equations (97) are computed numerically at a finer grid. At low photon energies, the redistribution function can be narrower than the whole integration interval, and integration can present a problem. In this case, however, we can extend the integration limits in Equations (97) from 0 to ∞ and to calculate the moments of the redistribution function analytically (NP94). Using the limits on γ_⋆, given by Equation (19), one can show that scattering takes place entirely within the central interval ∈ for incident photons and electrons satisfying the following relations:

$$\begin{equation} x < \frac{\delta _x}{2}, \qquad p < p_{\star }^{-}(x) = \frac{\delta _x}{2} - x. \end{equation} \tag{ 98 }$$

We can write the moments of the redistribution function in a way similar to Equation (84):

$$\begin{equation} \overline{x_1^i} \overline{s}_0^{<}(x) \equiv \frac{4\pi }{x} \int x_1^{i+1} \, {d}x_1 \; \overline{R}_{\rm ph}^{<}(x_1,x), \end{equation} \tag{ 99 }$$

where the < superscript signifies that only electrons with p < p⁻_⋆(x) are taken into account. Equations (97) then (for x < δ_x/2) become

$$\begin{eqnarray} &&\dot{x}_{\rm c}= c \sigma _{\rm T}{N_{\rm e}}\ \overline{(x_1-x)} \overline{s}_0^{<}(x),\quad D_{\rm ph}(x) = c \sigma _{\rm T}{N_{\rm e}}\ \overline{(x_1-x)^2} \overline{s}_0^{<}(x), \end{eqnarray} \tag{ 100 }$$

and can be computed using Equations (C5) and (C6).

For numerical differencing Equation (92) has to be written in the form

$$\begin{equation} \dot{n}_{\rm ph,diff,cs}(x) = -\frac{\partial }{\partial \ln {x}} \left[ A_{\rm ph,cs}(x) {n_{\rm ph}}(x) - B_{\rm ph,cs}(x) \frac{\partial {n_{\rm ph}}(x)}{\partial \ln {x}} \right], \end{equation} \tag{ 101 }$$

where

$$\begin{eqnarray} &&A_{\rm ph,cs}(x) = \frac{\dot{x}_{\rm c}}{x} - \frac{\partial }{\partial x} \left(\frac{1}{2} \frac{D_{\rm ph}(x)}{x} \right),\quad B_{\rm ph,cs}(x) = \frac{1}{2} \frac{D_{\rm ph}(x)}{x^2}. \end{eqnarray} \tag{ 102 }$$

The numerical treatment of pair-production and annihilation processes in our code is fairly straightforward. The only potential difficulty can arise from the nonlinearity of the absorption term in the photon equation. To deal with this we have chosen the simplest possible approach: for calculating the pair-production opacity at each step we simply use the photon distribution from the previous step. The error caused by doing so is not expected to be significant in most cases. It is well known that a photon of energy x will most efficiently interact with photons of energy x₁ ≈ 3/x, thus if its energy is not very close to the electron/positron rest energy, the photon will most likely annihilate on another photon of a vastly different energy than its own. Therefore, we can visualize two separate populations of photons that pair-produce on each other, with the dividing energy at m_ec². The photon distribution from the previous step is then taken to be the "target" population on which the photons that are being evolved pair produce.

Since we wish the numerical scheme to treat electrons and positrons identically (particularly when we are dealing with pure pair plasma), while at each step one of them has to be evolved first when the outcome of the other is yet unknown, we use a fully implicit scheme for the pair-annihilation terms.

The quantities R_γγ(γ₋, x, x₁), σ_pa(γ₊, γ₋), and σ_pp(x, x₁) defined by Equations (37), (39), and (46) are precalculated on a fine grid and thereafter averaged within the electron/positron and photon bins used by the code. The integrals in the Equations (36), (38), (44), and (45) for emissivities and absorption coefficients are calculated through discrete sums.

One of the main difficulties in numerically treating synchrotron processes in compact sources is that the optical thickness of the medium due to self-absorption might become extremely large at low energies compared to, say, Thomson optical thickness. Almost all photons that are produced are immediately absorbed, so very few escape. But the energy which those few carry away comes from the small net energy exchange rate between electrons and photons, which we need to keep track of to maintain the energy balance. Near the equilibrium, in the photon equation we have two large terms describing emission and absorption, which nearly exactly cancel out. A small error in either of them produces a significant error in the total energy transfer rate. In the electron equation this transfer rate is given by the difference in the synchrotron cooling and heating rates. To maintain the energy balance between the two equations, we need to ensure that in our numerical scheme these rates are seen identically by both equations.

In discretized form, the synchrotron processes for electrons/positrons are described by Equations (71) and (72), with n = n_±, A = A_e,syn, and B = B_e,syn. To obtain the total energy gain we have to multiply Equation (71) by γ_iΔ_p, sum over i, and sum the corresponding terms in the electron and positron equations. Assuming vanishing boundary currents, we have

$$\begin{equation} \frac{\Delta E_{\rm e}}{\Delta t_{k}\: \Delta V} = \sum _{i=1}^{i_{m} - 1} \Delta \gamma _{i+1/2} \left[ A_{i+1/2} \; n_{{\rm e},i+1/2} - B_{i+1/2} \frac{n_{{\rm e},i+1} - n_{{\rm e},i} }{\Delta _p} \right], \end{equation} \tag{ 103 }$$

where Δγ_i+1/2 ≡ γ_i+1 − γ_i and we have omitted the time index k + 1/2 for brevity. The exchange rate as seen by the photon equation can be evaluated by writing the integrals in emissivity and absorptivity Equations (56) and (58) as sums over the grid, multiplying Equation (55) by x_lΔ_x and summing over l:

$$\begin{eqnarray} \frac{\Delta E_{\rm ph}}{\Delta t_{k}\: \Delta V} &=& \sum _{l=1}^{l_m} \left[- c \: x_l \: \alpha _l \: n_{{\rm ph},\: l} + x_l \: \epsilon _{\, l} \right] \Delta _x = \frac{\lambda _{\rm C}^3}{8\pi } \sum _{l=1}^{l_m} \left[ \frac{n_{{\rm ph},\: l} }{x_l} \: \Delta _x\sum _{i=1}^{i_m} \frac{\gamma _i P_{l,i}}{p^2_i} \Delta _p\left(-3 n_{{\rm e},i} \: + \frac{n_{{\rm e},i+1} - n_{{\rm e},i} }{\Delta _p} \right) \right] \nonumber\\ && + \sum _{l=1}^{l_m} \left[ x_l \: \Delta _x\sum _{i=1}^{i_m} n_{{\rm e},i} \: P_{l,i} \: \Delta _p\right], \end{eqnarray} \tag{ 104 }$$

where P_l,i = P(x_l, γ_i). Changing the order of summation, identifying the sum over the photon distribution as the discretized version of the definition H(p), and noting that ∑_lP_l,ix_lΔ_x gives the electron cooling rate $-\dot{\gamma }_{{\rm s},i}$, we obtain

$$\begin{eqnarray} &&\frac{\Delta E_{\rm ph}}{\Delta t_{k}\Delta V} = \sum _{i=1}^{i_m} \Delta _p\bigg[{-} \bigg(\dot{\gamma }_{{\rm s},i} + \frac{3\gamma _i H_{i} }{p^2_i} \bigg) n_{{\rm e}, i} + \frac{\gamma _i H_{i} }{p^2_i} \frac{n_{{\rm e}, i+1} - n_{{\rm e}, i} }{\Delta _p} \bigg]. \end{eqnarray} \tag{ 105 }$$

To make Equations (103) and (105) identical (except for the sign) we have to make subtle changes in the definition of coefficients and the way integrals are numerically calculated. In Equation (105), we have to define the coefficients in between the electron momentum grid points, at i + 1/2, substitute the electron distribution n_e, i by n_e, i+1/2 (except in the derivative term), where the latter is calculated using the same Chang and Cooper coefficients δ_i as in the electron equation, and sum up to i = i_m − 1 instead of i_m. This amounts to defining the emission and absorption coefficients as

$$\begin{eqnarray} &&\epsilon _l = \sum _{i=1}^{i_m - 1} P_{l,i+1/2} \: n_{{\rm e},\: i+1/2} \: \Delta _p,\quad \alpha _l = \frac{\lambda _{\rm C}^3}{8\pi } \frac{1}{x_l^2} \sum _{i=1}^{i_m-1} \frac{\gamma _{i+1/2} \: P_{l,i+1/2}}{p_{i+1/2}^2} \left[ 3 n_{{\rm e},i+1/2} - \frac{n_{{\rm e},i+1} - n_{{\rm e},i} }{\Delta _p} \right] \: \Delta _p. \end{eqnarray} \tag{ 106 }$$

Also, the coefficients A and B entering the momentum space flux (Equation (72)) and thus also the electron energy exchange rate (Equation (103)) should be written as

$$\begin{eqnarray} &&A_{i+1/2} = \frac{\Delta _p}{\Delta \gamma _{i+1/2}} \left(\dot{\gamma }_{\rm s}+ 3 \frac{\gamma }{p^2} H \right)_{i+1/2} \quad \mbox{and}\quad B_{i+1/2} = \frac{\Delta _p}{\Delta \gamma _{i+1/2}} \left(\frac{\gamma }{p^2} H \right)_{i+1/2}, \end{eqnarray} \tag{ 107 }$$

which become identical to Equation (62) in the limit Δ_p → 0 and ensure that the energy exchange rates as seen by the electron and photon equations are the same.

The only discrepancy left is that we cannot use the same n^k+1/2_ph and n^k+1/2_e in both equations. This is because each of them contains a function n^k+1, which, in the equation that we evolve before, is not known for the other type of particle. The solution to this, at least in the average sense, is to regard the time grids for each equation as shifted by a half time step. Then n^k+1 obtained from one equation can be used as n^k+1/2 in the other and vice versa.

Coulomb scattering only redistributes the energy between different parts of the lepton population. It is easy to see that the total energy is conserved in the sum of two Equations (65) for electrons and positrons, provided that a(γ, γ₁) is antisymmetric, the latter simply reflects the energy conservation in two-body interactions. Similarly to synchrotron, our numerical treatment has to ensure that the conservation is exact, otherwise unphysical runaways can occur near the equilibrium.

The flux in momentum space in Equation (71) for Coulomb scattering is given by Equation (72) with coefficients expressed as (see Equation (66))

$$\begin{eqnarray} &&A_{i+1/2} = \left(\frac{\dot{\gamma } \gamma }{p^2} \right)_{i+1/2} - \frac{1}{2\Delta \gamma _{i+1/2}} \left[ \left(\frac{\gamma D}{p^2} \right)_{i+1} - \left(\frac{\gamma D}{p^2} \right)_{i} \right], \quad B_{i+1/2} = \frac{1}{2} \left(\frac{\gamma ^2 D}{p^4}\right)_{i+1/2}. \end{eqnarray} \tag{ 108 }$$

The total energy exchange rate is identical to Equation (103) for synchrotron.

Let us now look separately at terms containing $\dot{\gamma }$ and D. For $\dot{\gamma }$ we have

$$\begin{eqnarray} &&\left. \frac{\Delta E_{\rm e}}{\Delta t_{k}\: \Delta V} \right|_{\dot{\gamma }} = \sum _{i=1}^{i_{m} - 1} \Delta \gamma _{i+1/2} \left(\frac{\dot{\gamma } \gamma }{p^2} \right)_{i+1/2} \; n_{{\rm e},i+1/2}. \end{eqnarray} \tag{ 109 }$$

It is now easy to see that this quantity can be made to vanish if we write

$$\begin{eqnarray} &&\dot{\gamma }_{i+1/2} = \sum _{l=1}^{i_{m} - 1} a(\gamma _{i+1/2},\gamma _{l+1/2}) \; n_{{\rm e},l+1/2} \Delta _p\quad \mbox{and}\quad \left(\frac{\dot{\gamma } \gamma }{p^2} \right)_{i+1/2} \rightarrow \dot{\gamma }_{i+1/2} \frac{\Delta _p}{\Delta \gamma _{i+1/2}}, \end{eqnarray} \tag{ 110 }$$

provided that a is antisymmetric. The terms containing D(γ) in the energy exchange rate are

$$\begin{eqnarray} &&\left.\! \frac{\Delta E_{\rm e}}{\Delta t_{k}\: \Delta V} \right|_{D} = {-}\frac{1}{2} \sum _{i=1}^{i_{m} - 1} \left\lbrace \left[ \left(\frac{\gamma D}{p^2} \right)_{i+1} - \left(\frac{\gamma D}{p^2} \right)_{i} \right] n_{{\rm e},i+1/2} + \left(\frac{\gamma D}{p^2}\right)_{i+1/2} \left(n_{{\rm e},i+1} - n_{{\rm e},i} \right) \right\rbrace, \end{eqnarray} \tag{ 111 }$$

where we have redefined the coefficient B as

$$\begin{equation} B_{i+1/2} \rightarrow \frac{1}{2} \frac{\Delta _p}{\Delta \gamma _{i+1/2}} \left(\frac{\gamma D}{p^2}\right)_{i+1/2}. \end{equation} \tag{ 112 }$$

One can see that Equation (111) has the form of an integral over a full differential and, as such, should vanish provided that D = 0 at the boundaries. To ensure this numerically for any electron distribution we write explicitly n_e,i+1/2 = (1 − δ_i)n_e,i+1 + δ_in_e,i and demand that the coefficient in front of n_e,i in Equation (111) is equal to zero for every i. Rearranging terms, we get

$$\begin{eqnarray} \left.\! \frac{\Delta E_{\rm e}}{\Delta t_{k}\: \Delta V} \right|_{D} &=& -\!\!\frac{1}{2} \sum _{i=2}^{i_{m} - 1} n_{{\rm e},i} \left\lbrace \vphantom{\frac{\gamma D}{p^2}_{i-1/2}} \delta _i \left[ \left(\frac{\gamma D}{p^2} \right)_{i+1} - \left(\frac{\gamma D}{p^2} \right)_{i} \right]+ (1-\delta _{i-1}) \left[ \left(\frac{\gamma D}{p^2} \right)_{i} - \left(\frac{\gamma D}{p^2} \right)_{i-1} \right] \right.\nonumber \\ &&- \left. \left(\frac{\gamma D}{p^2}\right)_{i+1/2} + \left(\frac{\gamma D}{p^2}\right)_{i-1/2} \right\rbrace + n_{{\rm e},1} S^{-} + n_{{\rm e},i_{m}} S^{+}. \end{eqnarray} \tag{ 113 }$$

The expression in the curly brackets is identically zero if we set

$$\begin{equation} \left(\frac{\gamma D}{p^2}\right)_{i+1/2} = \delta _i \left(\frac{\gamma D}{p^2}\right)_{i+1} + (1-\delta _i) \left(\frac{\gamma D}{p^2}\right)_{i}, \end{equation} \tag{ 114 }$$

while the boundary terms S⁻ and S⁺ vanish if

$$\begin{equation} \left(\frac{\gamma D}{p^2}\right)_{1} = 0 \quad \mbox{and} \quad \left(\frac{\gamma D}{p^2}\right)_{i_m} = 0. \end{equation} \tag{ 115 }$$

Using Equations (110) in the first term in coefficient A and Equations (114) and (115) in the definition (112), we ensure precise energy conservation in the numerical scheme.

As a first test we compare our code to the well-known pair plasma code eqpair by Coppi (1992, 1999). eqpair also considers a uniform emission region into which high-energy electrons/pairs are injected, mimicking an unspecified acceleration mechanism. Some low-energy photons are also injected, emulating a source of external soft radiation (e.g., accretion disk). The high-energy pairs cool by Compton scattering and Coulomb energy exchange with colder thermal pairs. The Compton upscattered photons can produce electron–positron pairs which then upscatter more photons, etc., initiating a pair cascade. Once the pairs cool down to low enough energies, the timescale of the systematic energy losses becomes longer than that of diffusive processes, leading to relaxation into a low-energy thermal distribution. In eqpair, Coulomb collisions between particles are assumed to be the thermalizing mechanism. However, the thermalization process is not treated entirely consistently in this code in a sense that there exists only one thermal bin into which particles are put once they have cooled below a certain threshold energy, chosen to be γ = 1.3. The electron temperature associated with this thermal bin is nevertheless calculated self-consistently from energetic considerations. Furthermore, the code does not consider thermalization by synchrotron self-absorption, which can be an efficient mechanism if the medium is magnetized (Ghisellini et al. 1988, GHS98).

The setup of this test run is similar to what was used for Figure 1 in Coppi (1999). We switched off synchrotron processes in our code and left other processes. We inject a Gaussian distribution of pairs centered at γ_inj = 10³ and a low-energy blackbody distribution of photons. In addition, there is a background electron plasma present with optical depth τ_p = 0.1. There is no escape term for pairs, meaning that all injected pairs eventually annihilate transferring their energy to the radiation field. The power injected as nonthermal pairs is parameterized by compactness

$$\begin{equation} l_{\rm nth} = \frac{\sigma _{\rm T}}{m_{\rm e}c^3} \frac{L_{\rm nth}}{R}, \end{equation} \tag{ 116 }$$

where L_nth is the injected luminosity (including rest mass) and R is the linear dimension of the emission region. Similarly, we define the compactness of the injected soft radiation as

$$\begin{equation} l_{\rm s} = \frac{\sigma _{\rm T}}{m_{\rm e}c^3} \frac{L_{\rm s}}{R}, \end{equation} \tag{ 117 }$$

where L_s is the relevant luminosity. To mimic acceleration with less than 100% efficiency, additional power is supplied to low-energy electrons in the form of continuous heating, parameterized by l_th. In Coppi (1999), this energy was just given to the thermal bin, but since we do not have such bin in our code, we need to explicitly specify the form of this heating. This is done by stochastic acceleration prescription of the form

$$\begin{equation} \left.\frac{D \tilde{n}_{\pm }(p)}{Dt}\right|_{\rm {stoch.}} = \frac{1}{p^2} \frac{\partial }{\partial p} \left[ p^2 D_{\rm acc}(p) \frac{\partial \tilde{n}_{\pm }(p)}{\partial p} \right]. \end{equation} \tag{ 118 }$$

The momentum diffusion coefficient is assumed to take the form characteristic of stochastic acceleration by resonant interactions with plasma waves (Dermer et al. 1996), D_acc(p) ∝ p^q. We have chosen q = 2 in our calculations. The mean energy gain rate of a particle resulting from Equation (118) is

$$\begin{equation} \left. \left\langle \frac{{d}\gamma }{{d}t} \right\rangle \right|_{\rm stoch.} = \frac{1}{p^2} \frac{\partial }{\partial p} [ \beta p^2 D_{\rm acc}(p) ], \end{equation} \tag{ 119 }$$

where β = p/γ is the particle speed. We can see that for a power-law diffusion coefficient the gain rate is proportional to p^q−1 in the relativistic regime, while in the nonrelativistic regime it is proportional to p^q. Our choice q = 2 means that at high energies Compton losses always overcome gains by stochastic acceleration, the main effect of the latter process is therefore the heating of low-energy pairs.

The differential term given by Equation (118) is included in the Chang and Cooper scheme on the same grounds with other continuous terms. Therefore, before discretization it has to be written in the form compatible with Equations (71) and (72):

$$\begin{equation} \left.\frac{D {n_{\pm }}(p)}{Dt}\right|_{\rm {stoch.}} = -\frac{\partial }{\partial \ln {p}} \left\lbrace D_{\rm acc}(p) \frac{1}{p^2} \left[ 3 {n_{\pm }}(p) - \frac{\partial {n_{\pm }}(p)}{\partial \ln {p}} \right] \right\rbrace. \end{equation} \tag{ 120 }$$

The results of the test are shown in Figure 1. Varying the amount of stochastic heating (l_th) keeping all other parameters constant, we see that we can well reproduce the behavior of the spectrum in Figure 1 in Coppi (1999). Just as expected by Coppi (1999), the equilibrium electron distribution is hybrid: Maxwellian at low energies with a nonthermal high-energy tail. Note that we get such shape even if we switch off Coulomb scattering. The thermal-looking distribution is produced by the stochastic heating itself, which gives a Maxwellian slope at low energies irrespective of the shape of D_acc(p), while the location of the peak of the distribution is determined by the balance between heating and Compton cooling. The behavior of the spectrum in response to varying the power of stochastic heating seen in Figure 1(a) was analyzed in detail by Coppi (1999), we are not going to repeat it here.

Zoom In Zoom Out Reset image size

For the second test, we compared our results with these of GHS98. They studied electron thermalization by synchrotron self-absorption in the presence of Compton cooling. The electron cooling, heating, and diffusion due to the synchrotron were described by Equation (59; without the last term), while Compton scattering was assumed to take place in the Thomson regime and contribute only to systematic cooling. Furthermore, the treatment was not fully self-consistent since only the electron equation was actually solved. While the equilibrium synchrotron spectrum was self-consistently calculated at each time step from the formal solution of the radiative transfer equation, the Comptonized spectrum was not. Thus, only the synchrotron spectrum entered the electron heating rate by self-absorption, while the radiation energy density needed to account for Compton cooling was estimated from energetic considerations.

We ran our code with the same parameters used to obtain the results in Figures 1 and 2 in GHS98. The pair production/annihilation and Coulomb scattering have been switched off for this test. High-energy electrons are injected into the emission region, with the total power (including rest mass) parameterized by the injection compactness l_nth. The magnetic compactness is defined by

$$\begin{equation} l_{\rm B} = \frac{\sigma _{\rm T}}{m_{\rm e}c^2} R U_{\rm B}, \end{equation} \tag{ 121 }$$

where U_B is the magnetic energy density. In addition there is an external source of soft blackbody photons assumed to arise from reprocessing half of the hard radiation by cold matter in the vicinity of the emission region. The electron escape timescale is fixed at t_esc = R/c.

Zoom In Zoom Out Reset image size

In the first case, the injected electrons have a Gaussian distribution peaking at γ = 10. The evolution of this distribution is followed in time as it cools and thermalizes by Compton and synchrotron processes. We can see that our results shown in Figure 2 are almost identical to those presented in Figure 1 in GHS98. However, we would like to stress that we also compute self-consistently the photon spectrum. We see the partially self-absorbed synchrotron bump at small energies, then the blackbody photons and two Compton scattering orders at higher energies.

In the second case, we calculated the steady-state particle distributions for different injection compactnesses. The injected electron distribution (per unit ln p) is

$$\begin{equation} Q_{\rm e} = Q_0 \frac{p^3}{\gamma ^2} \exp {\left(-\frac{\gamma }{\gamma _{\rm c}} \right)}, \end{equation} \tag{ 122 }$$

where γ_c = 3.33. The resulting equilibrium electron distributions plotted in Figure 3(b) are again very similar those obtained by GHS98 in their Figure 2. The corresponding radiation spectra shown in Figure 3(a) are computed self-consistently and simultaneously with the electron distribution (while the spectra in Figure 4 of GHS98 are calculated a posteriori, i.e., after the equilibrium electron distribution has been determined). As discussed in GHS98, if the source is strongly magnetically dominated, the equilibrium distribution is almost purely Maxwellian. When the injection compactness increases, Compton losses become non-negligible and the electrons cool down to lower energies before they have time to thermalize. Note that at the highest compactness (l_nth = 100) the temperature of the Maxwellian part of the distribution inferred from Figure 3(b) deviates appreciably from the one obtained by GHS98. This is caused by the fact that at high compactness a significant fraction of the soft radiation is Compton upscattered to energies comparable to the energies of the Maxwellian electrons. These photons are therefore not effective in cooling the electrons any further. However, in GHS98 Compton cooling is accounted for through a term proportional to the radiation energy density, which includes all photons, and therefore overestimates the cooling rate. Overall, the simple prescription for Compton cooling without actually solving the photon equation appears to work well in the parameter regimes considered here.

Zoom In Zoom Out Reset image size

Finally, we compare our code to the large particle Monte Carlo code by Stern et al. (1995), with all the processes operating now. The setup is similar to that used in Stern & Poutanen (2004) for simulating the spectral evolution of gamma-ray bursts. They consider an initially optically thin distribution of electrons in a cylinder-shaped emission region. Arguing that impulsive first-order Fermi acceleration would result in cooling spectra that are too soft to be consistent with observations, energy is instead supplied to the electrons continuously, mimicking dissipation by plasma instabilities behind the shock front. As electrons are heated to relativistic energies in the prescribed background magnetic field, they emit synchrotron radiation, providing seed photons for Compton upscattering. The high-energy upscattered photons then initiate pair production.

In our simulation, we consider a spherical region permeated by magnetic field and start by heating a cold electron distribution (with initial Thomson optical depth τ₀ = 6 × 10⁻⁴) according to the stochastic acceleration prescription (120). No pair escape is allowed. The results of simulations are shown in Figure 4 and can be compared to a similar Figure 2 in Stern & Poutanen (2004). In both cases, the electrons are rapidly heated to about γ ∼ 100, as determined by the balance between stochastic heating and synchrotron cooling. As the photon field builds up, additional cooling by Compton scattering causes the electron "temperature" to start dropping. After about 1/3 of the light crossing time, the number of photons upscattered to the MeV range becomes large enough to start significant pair production. With the increasing pair density (at t = 1, opacity has grown by a factor of 20) the available energy per particle decreases, causing a further drop in the temperature of the now almost pure pair plasma. After about 10 light-crossing times the Thomson opacity is τ_T = 1.3 and the pair density reaches the value where the pair annihilation and creation rates are balanced and a steady state is attained.

Zoom In Zoom Out Reset image size

The spectral behavior seen in Figure 4(a) is similar to what was obtained by Stern & Poutanen (2004). The synchrotron peak rises first, being initially in the optically thin regime and thus following the evolution of the peak of the electron distribution according to x ∝ γ². The first Compton scattering order lags slightly behind synchrotron, while the second scattering order is initially in the Klein–Nishina regime and thus hardly visible at all. As the electron temperature drops and the peak of the first scattering order evolves to lower energies, the second order shifts to the Thomson regime and becomes comparable to and eventually dominant over the first order. At the same time the decreasing temperature and increasing pair opacity causes the synchrotron emission to switch to optically thick regime and the synchrotron luminosity to drop dramatically. The plasma becomes photon starved and the Comptonized spectrum hardens.