EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 506 / No 2 (November I 2009) A&A, 506 2 (2009) 589-599 Full HTML
Free Access
Issue
A&A
Volume 506, Number 2, November I 2009
Page(s) 589 - 599
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/20079113
Published online 27 August 2009

A&A 506, 589-599 (2009)

Numerical computation of isotropic Compton scattering

R. Belmont1,2

1 - CESR (Centre d'Étude Spatiale des Rayonnements), Université de Toulouse [UPS], 9 Av. du Colonel Roche, 31028 Toulouse, France
2 - CNRS, UMR5187, 31028 Toulouse, France

Received 20 November 2007 / Accepted 7 July 2009

Abstract
Context. Compton scattering is involved in many astrophysical situations. It is well known and has been studied in detail for the past fifty years. Exact formulae for the different cross sections are often complex, and essentially asymptotic expressions have been used in the past. Numerical capabilities have now developed to a point where they enable the direct use of exact formulae in sophisticated codes that deal with all kinds of interactions in plasmas. Although the numerical computation of the Compton cross section is simple in principle, its practical evaluation is often prone to accuracy issues. These can be severe in some astrophysical situations but are often not addressed properly.
Aims. In this paper we investigate numerical issues related to the computation of the Compton scattering contribution to the time evolution of interacting photon and particle populations.
Methods. An exact form of the isotropic Compton cross section free of numerical cancellations is derived. Its accuracy is investigated and compared to other formulae. Then, several methods to solve the kinetic equations using this cross section are studied.
Results. The regimes where existing cross sections can be evaluated numerically are given. We find that the cross section derived here allows for accurate and fast numerical evaluation for any photon and electron energy. The most efficient way to solve the kinetic equations is a method combining a direct integration of the cross section over the photon and particle distributions and a Fokker-Planck approximation. Expressions describing this combination are given.

Key words: radiation mechanisms: general - plasmas - methods: numerical - galaxies: active - X-rays: binaries - X-rays: galaxies

Introduction

Compton scattering has become a major process in astrophysics, and for more than half a century, it has been applied to a large variety of astrophysical situations, particularly in high energy, low density media. It was first introduced in the context of electron cosmic rays interacting with thermal stellar photons (Follin 1947; Donahue 1951; Feenberg & Primakoff 1948, etc.). It was then used to model gases of high energy electrons in the strong radiation field of X-ray or relativistic sources (Levich & Syunyaev 1971). It is also responsible for the up-scattering of low energy photons from the microwave background radiation (Zeldovich et al. 1972). More recently, much effort has been spent on high energy media in relativistic sources such as $\gamma$-ray bursts, AGN, and X-ray binaries. Particular focus is now on the synchrotron self-Compton radiation that is thought to play a major role for example in AGN (Tavecchio et al. 2001; Saugé & Henri 2004).

The many varied applications result in different energy regimes for the photons and the particles. The energy of particles spans from the very cold electrons of pulsars winds to the hot gas in the corona of X-ray binaries ($E\approx$ 100 keV, i.e. $E/m_{\rm e}c^2\approx 1$), and to high energy cosmic rays (up to E > 100 TeV, i.e.  $E/m_{\rm e}c^2 >10^{8}$). The energy of photons also spans a wide range: in the case of AGN for instance, it goes from low energy synchrotron photons ( $\lambda = 1$ m, i.e. $h\nu/m_{\rm e}c^2 \approx10^{-12}$) to $\gamma$-ray emission (with E > 10 TeV, i.e.  $h\nu/m_{\rm e}c^2 > 10^7$ for some sources). General formulae must be able to deal simultaneously with very small and very large energies.

The exact cross section is quite complex and does not give any physical insight into the fundamental process. Moreover, often, only a fraction of the particle and photon spectra contribute significantly to the Compton scattering in one given astrophysical situation. For these reasons most of the theoretical work has been done on deriving limiting and/or averaged forms that can be easily included in the modeling of astrophysical objects. For example, much work is now based on the diffusion approach for particles, valid when large angle scattering can be neglected, that is, when the energy of the particles does not vary significantly in one single interaction (work initiated by Kompaneets 1957). Analytical expressions averaged over Maxwellian distributions and thus only valid in purely thermal plasmas have also been used. Many other approximations have been proposed in general papers, reviews and books (see Rybicki & Lighman 1979; Blumenthal & Gould 1970; Nagirner & Poutanen 1994; Zeldovich 1975; Katz 1987, for example). Nevertheless the relevant energy cannot always be confined to a narrow range where simple approximations can be made. Also more general expressions that can be used regardless of a specific problem are required for generic codes designed to address various situations. In these cases an exact description without limitation on the photon and electrons energies is required. A few exact formulae have been proposed for the isotropic Compton cross section (Nagirner & Poutanen 1994; Brinkmann 1984; Jones 1968, NP94 hereafter). Such formulae are analytically equivalent. When evaluated numerically however, they behave differently with respect to accuracy issues. For example, the original formula by Jones (1968) fails to describe the up-scattering of low energy photons which represents a typical astrophysical situation. Although these problems are well known, the Jones' cross section is still widely used in the literature.

Such general expressions are not very helpful for physical understanding. Nevertheless, they are complementary to asymptotic forms for they provides quantitative results that can be compared to observations. It has become particularly true recently, since more and more codes have been developed to deal with photon-photon, photon-particle and particle-particle interactions in high energy plasmas (e.g. Vurm & Poutanen 2009; Belmont et al. 2008; Coppi 1992; Poutanen & Svensson 1996; Nayakshin & Melia 1998). Even with an accurate, exact cross section, numerical computation of the Compton contribution to the evolution of the particle and photon distributions is problematic. When the distributions and cross section are discretized in energy bins, the grid resolution puts severe constraints on the code accuracy. Typically, high energy particles are scattered down by numerous and inefficient scattering events off soft photons. If the resolution is insufficient, these individual events are not captured by the numerical scheme and the global evolution of the particle distribution is not described accurately. Nayakshin & Melia (1998) proposed that a combination of this method with a Fokker-Planck approximation could enable good accuracy. This idea was applied in a recent numerical code by Belmont et al. (2008) (see also Vurm & Poutanen 2009). Here we investigate the method and quantify its accuracy.

In the first section, we derive a new form of the distribution of scattered photons. The accuracy of this formula is discussed and compared to other expressions. In the second section, numerical errors of several kinetic methods are estimated and compared.

1 Evaluating the isotropic distribution of Compton-scattered photons

1.1 Exact spectrum

As long as stimulated scattering is not considered, the evolution of a photon population interacting by Compton scattering with an electron distribution is described by the following equation (see Rybicki & Lighman 1979, for example):

                              $\displaystyle % \frac{\partial f_\omega(\vec{k}) }{\partial t}$ = $\displaystyle \int {\rm d}\vec{p}_0 f_{\rm e}(\vec{p}_0) \int {\rm d}\vec{k}_0 ... ...\rm d}\sigma}{{\rm d}\vec{k}} (\vec{p}_0,\vec{k}_0 \rightarrow \vec{k}) \right.$  
    $\displaystyle \left. - ~ f_\omega(\vec{k}) c\frac{{\rm d}\sigma}{{\rm d}\vec{k}_0}(\vec{p}_0, \vec{k} \rightarrow \vec{k}_0)\right].$ (1)

Here $f_{\omega}={\rm d}{\cal N}_{\omega}/{\rm d}^3\vec{r}/{\rm d}^3\vec{k}$ and $f_{\rm e}={\rm d}{\cal N}_{\rm e}/{\rm d}^3\vec{r}/{\rm d}^3\vec{p}$ are the photon and lepton distributions in their respective 6-dimension momentum space. The first term in the brackets describes the contribution of all photons of momentum $\vec{k}_0$ being scattered to momentum $\vec{k}$ after one interaction. The second one describes photons of momentum $\vec{k}$ being scattered to any other momentum. Both are proportional to the Compton differential cross section ${\rm d}\sigma/{\rm d}\vec{k}(\vec{p}_0,\vec{k}_0 \rightarrow \vec{k})$ giving the probability of one photon of momentum $\vec{k}_0$ being scattered to a momentum $\vec{k}$ by an electron of momentum $\vec{p}_0$. This cross section can be interpreted as the distribution of scattered photons resulting from the interaction of mono-energetic photons with one single electron and it will be often be referred to as such hereafter.

When the photon and particle distributions are isotropic, this equation can be integrated over all directions. By noting ${\rm d}^3\vec{p}/(m_{\rm e}c)^3=p^2 {\rm d}p {\rm d} \vec{\Omega}_p$ and ${\rm d}^3\vec{k}/(m_{\rm e}c)^3=\omega^2 {\rm d}\omega {\rm d} \vec{\Omega}_\omega$ (where $\omega=h\nu/m_{\rm e}c^2$), it yields:

                   $\displaystyle % \frac{\partial N_\omega(\omega)}{\partial t}$ = $\displaystyle \int \hspace{-.3cm} \int c \frac{{\rm d}\bar{\sigma}}{{\rm d}\omega}(p_0,\omega_0\rightarrow\omega) {\rm d}N_\omega(\omega_0){\rm d} N_{\rm e}(p_0)$  
    $\displaystyle - ~ N_\omega(\omega) \int c\bar{\sigma}_0(p_0,\omega_0) {\rm d} N_{\rm e}(p_0)$ (2)

where $N_{\omega}(\omega) = 4\pi \omega^2 f_\omega(\vec{k})$ and $N_p(p) = 4\pi p^2 f_p(\vec{p})$ are the angle-integrated distributions. Here

\begin{displaymath}% \bar{\sigma}_0(p_0,\omega_0) = \int \frac{{\rm d} \bar{\sigma}}{{\rm d}\omega}(p_0,\omega_0\rightarrow\omega) {\rm d}\omega \end{displaymath} (3)

is the total scattering cross section, and

\begin{displaymath}% \frac{{\rm d}\bar{\sigma}}{{\rm d}\omega}(p_0,\omega_0\righ... ...}\omega {\rm d}\vec{\Omega}_\omega} {\rm d}\vec{\Omega}_\omega \end{displaymath} (4)

is the angle-averaged cross section for isotropic Compton scattering.

\begin{figure} \par\includegraphics[width=8.5cm,clip]{9113fg01.eps} \end{figure} Figure 1:

Angles and notations. The incoming photon direction $\vec{n}_0$ is defined by the incident angle $\theta _0$ and an azimuthal angle $\psi _0$ that orientates the global picture (not shown here). The scattered photon direction $\vec{n}$ is defined by the scattering angle and the azimuthal angle $\phi $. In coordinates $(\vec{e}_1,\vec{e}_2,\vec{e}_3)$, $\vec{\beta}_0/\beta_0=(\mu_0,\sin{\theta_0},0)$ and $\vec{n}=(c_\alpha,\sin{\alpha}\cos{\phi},\sin{\phi})$, so that $\mu = \mu_0 c_\alpha+ \sin{\alpha}\sin{\mu_0}\cos{\phi}$.
Open with DEXTER

Equation (4) can actually be simplified further. Because of symmetry, the differential cross section ${\rm d}\sigma/{\rm d}\omega/{\rm d}\vec{\Omega}_\omega$does not depend independently on the direction of the incoming photon and lepton, but only on the angle between their directions: $\mu_0 = \cos{\theta_0}$, $\theta_0=(\vec{\beta}_0,\vec{n}_0)$ where $\beta_0=p_0/\gamma_0$ is the lepton velocity in units of the speed of light and $\vec{n}_0$ is the direction of the incoming photon (see Fig. 1). Then, Eq. (4) reduces to:

\begin{displaymath}% \frac{{\rm d}\bar{\sigma}}{{\rm d}\omega}(\vec{p}_0,\omega_... ...\omega {\rm d}\vec{\Omega}_\omega} {\rm d}\vec{\Omega}_\omega. \end{displaymath} (5)

Also, although the result of a single scattering event is described by 6 variables (energy and direction of the scattered photon and particle), these variables are not independent. The conservation of the total energy-momentum 4-vector sets 4 constrains, so that only 2 independent variables are necessary to describe the scattering outcome. Although other choices can be made, these are often chosen to be the angles describing the scattered photon direction $\vec{n}$: the cosine of the scattering angle $c_\alpha=\cos{\alpha}$, $\alpha=(\vec{n}_0,\vec{n})$ and an azimuthal angle $\phi $ (see Fig. 1). Then, the differential cross section ${\rm d}\sigma/{\rm d}\omega/{\rm d}\vec{\Omega}_\omega$ is redundant and includes a delta-function that ensures that the total energy-momentum is conserved:

\begin{displaymath}% \frac{{\rm d}\bar{\sigma}}{{\rm d}\omega}(\vec{p}_0,\omega_... ...mega - \bar{\omega}) {\rm d}\mu_0 {\rm d}\phi {\rm d}c_\alpha. \end{displaymath} (6)

Here the latter condition is given by:

\begin{displaymath}% \bar{\omega} = \omega_0 \frac{1-\beta_0\mu_0}{1-\beta_0\mu +\frac{\omega_0}{\gamma_0}(1-c_\alpha)} \end{displaymath} (7)

where $ \mu=\mu_0c_\alpha + \sqrt{1-\mu_0^2}\sqrt{1-c_\alpha^2} \cos{\phi}$ is the cosine of the angle between the direction of the scattered photon and the incoming electron.

The angle differential cross section is a well known result of quantum mechanics. In the rest frame of the electron, the general Klein-Nishina cross section does not depend on the azimuthal angle and reads (Heitler 1954; Rybicki & Lighman 1979):

\begin{displaymath}% \frac{{\rm d}^2 \sigma'}{{\rm d}\phi' {\rm d}c_\alpha'} = \... ...\rho'^2 \left( \rho' + \frac{1}{\rho'} -1 +c_\alpha'^2 \right) \end{displaymath} (8)

where $\sigma_{\rm T}$ is the Thomson cross section, $\rho'(c_\alpha) = \omega'/\omega'_0 = 1/(1+\omega'_0(1-c_{\alpha}'))$ is the ratio of the scattered to incoming photon energy. Here, all variables denoted with ' refer to quantities in the rest frame of the particle.

Equations (6) and (8) can be shown to be equivalent to Eq. (12) of Jones (1968) and Eqs. (6.2.1), (8.1.1) and (8.1.3) of NP94. Getting the angle-averaged cross section is now quite straightforward. However, including the Klein-Nishina cross section in Eq. (6) requires changes of frame and the triple integral quickly leads to cumbersome expressions. Different integration variables can be used and the triple integral can be performed in several orders. Here we follow the original method by Jones (1968). All quantities are expressed in the electron rest frame which gives:

\begin{displaymath}% \frac{{\rm d}\bar{\sigma}}{{\rm d}\omega}= \int \frac{{\rm ... ...{{\rm d}\mu_0'}{2} \frac{{\rm d}\phi'}{2\pi} {\rm d}c_\alpha'. \end{displaymath} (9)

Integration is then performed over $\mu_0'$, $c_\alpha'$ and then $\phi'$. The latter integration is made with the variable $1+\beta_0\mu'(\phi')$. The resulting formula by Jones (1968) contains several misprints[*]. Using slightly different notations, the correct formula for the scattered distribution reads (see also Pe'er & Waxman 2005):

\begin{displaymath}% \frac{{\rm d}\bar{\sigma}}{{\rm d}\omega}(p_0,\omega_0;\ome... ...m T}}{8\gamma_0p_0\omega_0^2}\left[ F(x_-,x_+) \right]_{a}^{b} \end{displaymath} (10)

with
                             F = $\displaystyle \frac{L_+}{\omega \omega_0} + \frac{1}{2L_+}\left( x_--\frac{1}{x_-} -4 -(\gamma/\omega+\gamma_0/\omega_0) \right)$  
    $\displaystyle + ~\frac{1}{2L_-}\left( x_--\frac{1}{x_-} - 4 + (\gamma_0\omega_0+\gamma\omega) \right)$  
    $\displaystyle + ~\frac{c_+}{\lambda_+}\left( \frac{1}{L_+} - {\rm acosh} \left(\sqrt{\lambda_+x_+/2}\right)/\sqrt{\lambda_+} \right)$  
    $\displaystyle + ~ \frac{c_-}{\lambda_-} \left( \frac{1}{L_-} - \left\{ \begin{a... ...)/\sqrt{\vert\lambda_-\vert} ~~\rm {if}~ \lambda_->0 \end{array}\right. \right)$ (11)

where $\gamma_0$ and $\gamma=\gamma_0+\omega_0-\omega$ are the Lorentz factors of the incoming and scattered particles respectively and $p=\sqrt{\gamma^2-1}$ is the momentum of the scattered particle. Here the following notations have been used:

$\displaystyle \lambda_+ = p^2_0 + 2\gamma_0\omega_0 + \omega_0^2$ (12)
$\displaystyle \lambda_- = p_0^2 -2\gamma_0\omega~ + \omega^2$ (13)
$\displaystyle L_\pm = \sqrt{\lambda_\pm \mp 2/x_\pm}$ (14)
$\displaystyle c_\pm = 1+\omega\omega_0+2\lambda_\pm.$ (15)

The cross section is proportional to the difference $\left[F\right]_a^b$ of function F evaluated at the upper and lower boundaries of the last integration (b and a respectively). At these boundaries, the variables $x_\pm$ have the following expressions:

  $\displaystyle x_+^a = {\rm min}\left[(\gamma_0+p_0 )/\omega_0; \; (\gamma+p)/\omega~ \right], \; x_-^a = 1/(\omega\omega_0x_+^a)$ (16)
  $\displaystyle x_-^b = {\rm min}\left[(\gamma_0+p_0)/\omega; \; (\gamma+p)/\omega_0 \right], \; x_+^b = 1/(\omega\omega_0x_-^b).$ (17)

The critical photon frequencies for which the arguments of the ${\rm min}$ and ${\rm max}$ functions are equal define distinct domains of the scattered distribution. For the lower integration boundary a, both arguments are equal at a unique frequency: $\omega=\omega_0$. For the higher boundary, both arguments are equal at two critical frequencies: $\omega=\omega_0$ and $\omega=\omega_{\rm c}$, where

\begin{displaymath}% \omega_{\rm c} = \omega_0 \frac{\gamma_0+p_0}{\gamma_0-p_0 + 2\omega_0}\cdot \end{displaymath} (18)

The scattered distribution is obviously continuous at these frequencies. However, higher derivatives are not, as it can be seen in Fig. 4 for example.

The scattered distribution is bounded in energy and the maximal and minimal energies guarantee that all radicals remain real (see Jones 1968, for more details):

  $\displaystyle \omega_{\rm min} = \omega_0 \frac{\gamma_0-p_0}{\gamma_0+p_0 + 2\omega_0}$ (19)
  $\displaystyle \omega_{\rm max} = \left\{ \begin{array}{ll} \omega_0+\gamma_0-1 ... ... \omega_{\rm c} & ~~~{\rm if }~~~\omega_0 < \omega_0^*(p_0), \end{array}\right.$ (20)

where

\begin{displaymath}% \omega_0^*=\frac{1}{2}(1+p_0-\gamma_0). \end{displaymath} (21)

When, $\omega_0>\omega_0^*$, the photon can gain up to the entire electron energy. However, this absolute limit is not necessarily reached. When $\omega_0<\omega_0^*$ (as for Compton up-scattering for instance), only a fraction of the electron energy is at most transfered to the photon.

1.2 Cancelation issues

As was first discussed by Jones (1968), evaluation of Eq. (11) suffers from accuracy issues due to machine round-off errors. Since most astrophysical problems involve an energy range spanning many orders of magnitude, Eq. (11) includes combinations of very small and very large terms. In this form, most large terms must cancel out, which obviously represents a numerical challenge. An example of accuracy issue is shown in Fig. 4.

Cancelation errors appear when the relative difference of two terms is smaller than the machine precision. Here we focus on computations in double precision, i.e. when reals are coded in 8 octets (64 bits) and their mantissa is coded in 52 bits. The relative error due to machine precision is then of the order of $\epsilon = 2^{-52} \approx 10^{-15}$.

We have checked the accuracy of the scattered distribution for a large range of photon and lepton energies and identify the domain of the $(\omega _0,p_0)$ plane where Eq. (11) gives accurate results. Results are shown in Fig. 2. We find that Compton scattering of mild photons ( $\omega_0\sim1$) by mid-relativistic particles ($p_0\sim 1$) is well described by this formula but problems arise for lower or higher energies. In particular, the scattered distribution cannot be evaluated in the most typical astrophysical case, namely the up-scattering of soft photons by high energy particles.

In the latter case, Jones (1968) derived asymptotic expansions in $\omega_0/\gamma_0$ ${\ll }$ 1. However, a large number of orders must be kept when $\omega_0/\gamma_0\la1$ and the numerical evaluation often becomes time consuming. Moreover, other situations are also affected by numerical issues, for which such an expansion is not relevant. Here, we rather concentrate on analytical manipulation of the original formula to get an exact expression free of cancelation issues.

1.3 Re-writing the exact cross section

Numerical accuracy issues result from two kinds of cancelations: when the relative difference of function F evaluated at the two integration boundaries is very small and when some terms constituting F should cancel out. Here we deal with both successively.

1.3.1 Cancelation issues when |F$^{\rm b}$ - F$^{\rm a}$| ${\ll }$ F$^{\rm a}$, F$^{\rm b}$

This happens typically when $x_\pm^a\approx x_\pm^b$, as for example, at the distribution upper and lower photon energies  $\omega_{\rm min}$ and  $\omega_{\rm max}$ where the distribution vanishes, but where Fa and Fb can remain large.

Differences of hyperbolic functions between the upper and lower boundaries can be computed analytically by using the following trigonometric relations:

$\displaystyle {\rm acosh} (\!\sqrt{a}) - {\rm acosh} (\!\sqrt{b}) = {\rm asinh}\left(\sqrt{a(b-1)}-\sqrt{b(a-1)} \right)$  
$\displaystyle {\rm asinh} (\!\sqrt{a}) - {\rm asinh} (\!\sqrt{b}) = {\rm asinh}\left(\sqrt{a(1+b)}-\sqrt{b(1+a)} \right)$  
$\displaystyle {\rm asin} (\!\sqrt{a}) - {\rm asin} (\!\sqrt{b}) = {\rm asin}\left(\sqrt{a(1-b)}-\sqrt{b(1-a)} \right).$  

The difference of the two hyperbolic functions between the two integration boundaries a and b simply reads: ${\rm asinh}(\Delta \sqrt{\lambda_\pm})$ if $\lambda_\pm >0$, or ${\rm asin}(\Delta\sqrt{-\lambda_\pm} )$ if $\lambda_\pm<0$, where
                          $\displaystyle % \Delta$ = L+(x+b)-L+(x+a) (22)
  = L-(x-b)-L-(x-a) (23)
  = $\displaystyle {\rm min}\left[p,p_0\right] + {\rm min}\left[ 0 ; \omega+\omega_0-p-p_0 \right]/2.$ (24)

Although the last term of Eq. (11) is not divergent it is numerically ill-behaved when $\lambda_- \ll 1$ and it is best to use the following auxiliary function:
S(x) = $\displaystyle \left\{ \begin{array}{c} \frac{1}{x} \left( \frac{{\rm asinh} (\s... ...})}{\sqrt{-x}} -\frac{1}{\sqrt{1+x}} \right) {~\rm if~} x<0 \end{array} \right.$ (25)

which is evaluated as a series expansion for small argument:

\begin{displaymath}% S(x) = \sum_{n=0}^{\infty} \frac{(2n+1)!(-x)^n}{4^n (n!)^2(2n+3)} ~~{\rm for} ~~\vert x\vert<1. \end{displaymath} (26)

Differences of all other terms in Eq. (11) between the two integration boundaries can also be computed analytically to factorize $\Delta$ and get the following quite simple expression for the scattered distribution:

\begin{displaymath}% \frac{{\rm d}\bar{\sigma}}{{\rm d}\omega} = \frac{3 \Delta}{8\omega_0^2\gamma_0p_0} G \end{displaymath} (27)

with
                       G = $\displaystyle 2 + (1+\omega\omega_0)\left[ z_+ + z_--\frac{2}{\omega\omega_0}\right]$  
    $\displaystyle + ~ 2 \left[A \right]_-^+ - \Delta^2(1+\omega\omega_0)\left[zA\right]_-^+$  
    $\displaystyle + ~ \Delta^2\left[ cz^{3/2}S\left(\lambda \Delta^2 z\right)\right]_-^+$ (28)

where $ A_\pm = (1/z_\pm+\lambda_\pm \Delta^2)^{-1/2} $ and $z_\pm=x_\pm^ax_\pm^b$. The scattered distribution is proportional to $\Delta$. This factor holds as the difference between the two integration boundaries a and b and vanishes at the minimal and maximal photon energies  $\omega_{\rm min}$ and  $\omega_{\rm max}$. For high particle energy, it is best computed as:
                       $\displaystyle % \Delta$ = $\displaystyle {\rm min} \left\{ {\rm min}(p;p_0);~ \frac{\gamma_0+\gamma+p_0+p}{2(p+p_0)} \right.$  
    $\displaystyle \times ~ \left. \left( {\rm min}(\omega;\omega_0) - \frac{{\rm max}(\omega;\omega_0)}{(\gamma+p)(\gamma_0+p_0)} \right) \right\}\cdot$ (29)

Here we note that an accurate evaluation of the scattered electron momentum p is crucial. When the latter is low (typically when p<10-3), the simple derivation $p=\sqrt{(\gamma_0+\omega_0-\omega)^2-1}$ is not accurate enough and an alternative derivation must used:

\begin{displaymath}% p = p_0 - \frac{(\omega-\omega_0)(p+p_0)(\gamma+\gamma_0)}{2p_0(p+p_0)-(\omega-\omega_0)(\gamma+\gamma_0)}\cdot \end{displaymath} (30)

Equations (27)-(28) are accurate over a larger domain than the original formula and can be used in many astrophysical situations.

\begin{figure} \par\includegraphics[width=8.5cm,clip]{9113fg02.eps} \end{figure} Figure 2:

Accuracy domains. Equation (11) is accurate for photons energy and electron momentum  $(\omega _0,p_0)$ in the region delimited by the solid line. The original formula derived by Jones (1968), and the cross section by Pe'er & Waxman (2005), have basically the same properties and are accurately evaluated in the same region. The published formula by NP94 is accurate in the domain determined by the dashed line. These boundaries are for double precision computing. The star shows the point for which the scattered distribution is plotted in Fig. 4.
Open with DEXTER

1.3.2 Cancelation when |z+ - z-| ${\ll }$ z+, z-

However, this new expression reveals other differences hidden in the original one and that result in cancelation issues, namely when $\vert z_+-z_-\vert \ll z_+,z_-$. In such cases, the differences must be computed analytically. Although alternate expressions for the differences can often be accurate in more cases, it is sufficient to use them when the relative differences are smaller than $\epsilon=10^{-10}$.

By using the definition of z+ and z-, the differences in the first three terms of Eq. (28) can be computed as:

$\displaystyle % z_+ + z_--\frac{2}{\omega\omega_0} = \left( \frac{1}{\omega}-\f... ...amma+p+\gamma_0+p_0)^2\omega\omega_0}{(\gamma_0+p_0)(\gamma+p)(p+p_0)^2}\right]$     (31)

  $\displaystyle \left[A \right]_-^+ = \frac{-A_-^2A_+^2}{A_++A_-} (\Delta^2 \left[\lambda\right]_-^+ - \omega^2\omega_0^2 \left[z\right]_-^+ )$ (32)
  $\displaystyle \left[zA\right]_-^+ = \frac{z_++z_-}{2}\left[A \right]_-^+ + \frac{A_++A_-}{2}\left[z\right]_-^+$ (33)

respectively, where
                            $\displaystyle % \left[ \lambda\right]_-^+$ = $\displaystyle (\omega+\omega_0)(\gamma+\gamma_0)$ (34)
$\displaystyle \left[z\right]_-^+$ = $\displaystyle -\left\vert\frac{1}{\omega}-\frac{1}{\omega_0} \right\vert$ (35)
    $\displaystyle \times ~ \mbox{min}\left[ \frac{(\gamma+p+\gamma_0+p_0)^2}{(p+p_0)(\gamma+p)(\gamma_0+p_0)},\frac{1}{\omega}+\frac{1}{\omega_0} \right]\cdot$ (36)

In the same manner, the last term can be computed as:

\begin{displaymath}% \left[ cz^{3/2}S\right]_-^+ = \frac{S_++S_-}{2}\left[cz^{3/... ...]_-^+ + \frac{c_+z_+^{3/2}+c_-z_-^{3/2}}{2} \left[S\right]_-^+ \end{displaymath} (37)

where
                    $\displaystyle % \left[cz^{3/2}\right]_-^+$ = $\displaystyle \frac{z_+^{3/2}+z_-^{3/2}}{2}\left[ \lambda\right]_-^+$ (38)
    $\displaystyle + ~\frac{c_++c_-}{2} \left[ z\right]_-^+\frac{1+z_+/z_-+z_-/z_+}{z_+^{1/2}/z_-+z_-^{1/2}/z_+}\cdot$ (39)

The numerical cancelation of $\left[S \right]_-^+$ fails either when $\lambda_+ z_+ \approx \lambda_- z_- $ or when $\Delta^2 \lambda_\pm z_\pm \ll 1$ and these two cases must be dealt with separately. In the former case, first order Taylor developments of S around the mid value $\Delta^2 (\lambda_+ z_++\lambda_- z_-)/2$ are sufficient whereas in the latter case (typically when $\Delta^2 \lambda_\pm z_\pm \ll 10^{-2} $), series expansions of S in zero can be used (Eq. (26)). In both cases, the difference is proportional to $\left[\lambda z \right]_-^+$ which must be computed as:

\begin{displaymath}% \left[\lambda z \right]_-^+ = \frac{z_++z_-}{2}\left[\lambda\right]_-^+ + \frac{\lambda_++\lambda_-}{2}\left[z\right]_-^+ \end{displaymath} (40)

for low photon energy ( $\omega_0 \ll 1$). With the additions of this subsection, Eq. (28) is accurately evaluated for any particle and photon energy.

1.4 Comparisons to other formulae

This expression was checked extensively against previous formulae. The exact differential cross section was first compared to the original Jones' formula (identical to 11 when the misprints are corrected) and to the formula by NP94[*] (Eqs. (8.1.7)-(8.1.22)). The latter was derived by integrating Eq. (9) in a different order and behaves numerically better than the Jones' cross section for Compton up-scattering of soft photons (see Fig. 2). Results show perfect agreement of the three formulae when the numerical evaluation of the three cross sections is accurate. The first moments of the distribution:

$\displaystyle \sigma_0(p_0,\omega_0) = \int \frac{{\rm d}\bar{\sigma}}{{\rm d}\omega} {\rm d}\omega$ (41)
$\displaystyle \sigma_1(p_0,\omega_0) = \int (\omega-\omega_0) \frac{{\rm d}\bar{\sigma}}{{\rm d}\omega} {\rm d}\omega$ (42)
$\displaystyle \sigma_2(p_0,\omega_0) = \int (\omega-\omega_0)^2 \frac{{\rm d}\bar{\sigma}}{{\rm d}\omega} {\rm d}\omega$ (43)

- namely the total cross section, the mean energy and the dispersion of the scattered photon energy - were also computed numerically by integrating over the scattered distribution and they were compared to analytical equations[*] (NP94, Eqs. (3.3.1)-(3.1.10)). Examples are shown in Fig. 3. Again, good agreement is obtained. Small differences are observed but they result from numerical errors in the integration of the scattered distribution.

\begin{figure} \par\hspace*{-1mm}\includegraphics[width=8.35cm,clip]{9113fg03.ep... ....eps}\vspace*{3mm} \includegraphics[width=8.3cm,clip]{9113fg05.eps} \end{figure} Figure 3:

First three moments of the scattered distribution as a function of the energy of the incoming photon: the total cross section, the mean energy of the scattered photon and the dispersion about this mean value. The solid and dashed lines show the numerical integration of the expression derived in this paper and analytical results given in NP94 respectively. The 3 curves are for p0=10-2,102 and 104. For each panel, relative errors are also shown.
Open with DEXTER

Compared to previous ones, the cross section derived here is found to be accurate over a wider energy range. An example of the scattered distribution is shown in Fig. 4. As summed-up in Fig. 2, the original cross section by Jones (1968) is inaccurate in typical astrophysical situations involving low energy photons up-scattered by high energy particles. The cross section derived by NP94 is accurate in most astrophysical cases. When computed directly from the published equations, it fails, however, for very high energy particles or very high energy photons. Although these cases may not be physically relevant since the Klein-Nishina cross section drops at these high energies, numerical issues can generate infinite values, propagate and affect numerical solutions. Nevertheless, it must be noted that a few changes in the way some quantities are computed enable good accuracy over a much larger domain (Poutanen, private communication). The cross section derived here is accurate for all photon and particle energies and overcomes these numerical issues.

\begin{figure} \par\includegraphics[width=8.2cm,clip]{9113fg06.eps} \end{figure} Figure 4:

The scattered distribution for soft photons ( $h\nu =10^{-5}m_{\rm e}c^2$) and mild relativistic particles (p=1). Solid and dotted lines show the distributions obtained with Eq. (28), and with Jones' formula respectively.
Open with DEXTER

2 Solving the kinetic equations

Once the Compton cross section is computed accurately, the Compton contribution to the evolution of interacting particles and photons is described by the set of equations:

                       $\displaystyle % \partial_t N_{\rm e}(p)$ = $\displaystyle \int\hspace{-.3cm} \int c \frac{{\rm d}\bar{\sigma}}{{\rm d}p}(p_0,\omega_0 ; \omega(p)) {\rm d}N_{\rm e}(p_0) {\rm d}N_\omega(\omega_0)$  
    $\displaystyle - ~ N_{\rm e}(p) \int c \bar{\sigma}_0(\omega_0,p) {\rm d} N_\omega(\omega_0),$ (44)


                   $\displaystyle % \partial_t N_\omega(\omega)$ = $\displaystyle \int \hspace{-.3cm} \int c \frac{{\rm d}\bar{\sigma}}{{\rm d}\omega}(p_0,\omega_0 ; \omega) {\rm d}N_{\rm e}(p_0) {\rm d}N_\omega(\omega_0)$  
    $\displaystyle - ~ N_\omega(\omega) \int c \bar{\sigma}_0(\omega,p_0) {\rm d}N_{\rm e}(p_0),$ (45)

where $\omega(p) = \omega_0 +\gamma_0-\gamma$ is the energy of the scattered photon and ${\rm d}\bar{\sigma}(p_0,\omega_0;\omega(p))={\rm d}\bar{\sigma}(p_0,\omega_0;\omega)$ is the differential cross section of the interaction $(p_0,\omega_0)\rightarrow(p,\omega)$. In this section we show that computing these integrals numerically can also lead to accuracy issues. We present an alternative approach that combines this integration with a Fokker-Planck treatment and overcomes most numerical issues. The basic idea of such a combination was proposed by Nayakshin & Melia (1998) for particles only but was not investigated in detail. A more precise analysis of this treatment was presented in Belmont et al. (2008) (see also Vurm & Poutanen 2009). Here we investigate the regimes where the various methods are valid and we focus on their numerical accuracy.

2.1 Integral approach

The simplest way to compute the time evolution of the lepton and photon populations is to discretize their distributions in energy bins. Integration over energy is then performed bin by bin by summing the contributions from all other bins. This process is time consuming since a double integral must be performed at each time step. Moreover it can be inaccurate in many situations: when the width of the scattered distribution is smaller than the bin size, the energy change of the scattered particles (or equivalently photons) is too small to be resolved by the numerical scheme. The particles (or photons) do not scatter far enough in the energy space to reach other bins and there is no numerical evolution of the distribution. It is well known that most numerical issues arise from what happens below the grid scale. In this peculiar case however, particles (or photons) can undergo many scattering events per unit time when the density of scattering centres is large. Errors associated with individual scattering events add and they can lead to accuracy issues. This is typically what happens in simulations of X-ray binaries: the soft disc photons are up-scattered by the high energy particles constituting the corona. This interaction is known to efficiently cool the corona down. Whereas the photons experience large energy changes, the relative variation of the particle energy is very small and their cooling may not be captured by a discrete integral with finite resolution. The energy is then not conserved since the photon field gains energy while the particle population does not lose any[*].

The accuracy of this integral method depends on the comparison of the energy bin size to the width of the scattered distribution. The latter depends on the momentum p0 of the incoming particle and the energy $\omega _0$ of the incoming photon. Most astrophysical cases involve large ranges of energy (typically many orders of magnitudes) which implies the need to use logarithmic grids, the resolution of which is also energy-dependent. The accuracy of the numerical integration thus strongly depends on the region in the lepton momentum-photon energy space  $(p_0,\omega_0)$the simulation has to deal with. In this subsection, the regimes where the integral approach remains accurate are investigated. Hereafter the subscript 0 is dropped for simplicity.

2.1.1 Photons

Numerical integration for the evolution of the photon population is accurate in the limit where the width of the scattered distribution $\Delta\omega(p,\omega)=\omega_{\rm {max}}(p,\omega)-\omega_{\rm {min}}(p,\omega)$ is much larger than the local photon bin size  $\delta \omega$. For linear grids, requiring that the scattered photon distribution spans at least n bins would simply read[*]: $\Delta\omega > (n-1) \delta\omega$. For logarithmic grids, the same condition reads: $\omega_{\rm {max}}/\omega_{\rm {min}} > r_\omega^{n-1} $ where $r_\omega=\omega^{i+1}/\omega^{i} $ is the logarithmic increment of the photon grid. By denoting $R_\omega$ the number of bins per decade energy, i.e. the resolution of the photon grid, and $\epsilon_{\rm {I,\omega}}=10^{(n-1)/R_\omega}-1$, the condition for the integral approach to be valid reads:

\begin{displaymath}% \Delta \omega(p,\omega)/\omega_{\rm {min}}(p,\omega) > \epsilon_{{\rm I},\omega}. \end{displaymath} (46)

When the resolution is large enough: $\epsilon_{\rm {I},\omega} \approx 2 (n-1)/R_\omega$. However, as n is typically of several, $(n-1)/R_\omega$ can be larger than unity for low resolution runs and the exact relation must be used. The condition 46 is satisfied in a well defined region of the $(p,\omega)$ space. Using Eqs. (19) and (20) to compute the scattered distribution width, one can write an equation in $\omega$ for the boundary of this region. It is the solution of a 2nd order polynomial and the direct numerical integration is found to be accurate for all photon and particle energies except in the region where:

\begin{displaymath}% 0 < p < p_{\rm I}^{\rm c} \quad {\rm and} \quad {\rm min}\l... ...,\omega_{\rm I}^-(p)\right] < \omega < \omega_{\rm I}^+(p) \\ \end{displaymath} (47)

with

\begin{displaymath}% p_{\rm I}^{\rm c} = \frac{(\epsilon_{{\rm I},\omega}+4) +(\... ...rt{1+\epsilon_{{\rm I},\omega}}}{4\epsilon_{{\rm I},\omega}}, \end{displaymath} (48)


                             $\displaystyle % \omega^\pm_{\rm I}$ = $\displaystyle \frac{1}{4} \left[ 2(1+p-2\gamma) + \epsilon_{{\rm I},\omega}(\gamma-p) \right.$  
    $\displaystyle \pm ~ \left. \sqrt{\epsilon_{{\rm I},\omega}}\sqrt{4(\gamma-p-1)+\epsilon_{{\rm I},\omega}(1-2p\gamma+2p^2)} \right],$ (49)


\begin{displaymath}% \omega_{\rm I}^{\rm c} = -\gamma + \frac{\epsilon/2}{(\gamma-p)\epsilon-2p}\cdot \end{displaymath} (50)

The inverse equation $p_{{\rm I}}(\omega )$ for this boundary would be more practical since it would be mono-valued and could be used directly in the combined approach (see hereafter). It is the solution of a high order polynomial and is most easily found numerically.

Boundaries for the integral approach are shown in Fig. 5 for $R_\omega \approx 9200,930,97,13,3.8,2.0,1.3$ and n=5. For comparison typical runs have a resolution $R_\omega=1{-}30$. In the case $R_\omega$ = 13, only the scattering of low energy photons ( $h\nu < 200$ keV) off sub-relativistic particles (p < 0.2) is inaccurately described by the integral approach. Although it is not general, the integral method for the photon equation can thus be used to address many of astrophysical situations, such as those involving high energy particles.

\begin{figure} \par\includegraphics[width=8.2cm,clip]{9113fg07.eps} \end{figure} Figure 5:

Photon equation: boundaries $p_{{\rm I}}(\omega )$ for the integral approach (solid lines) and $p_{{\rm FP}}(\omega )$ for the Fokker-Planck approximation (dashed lines) for $\epsilon _{{\rm I},\omega }=10^{-3},10^{-2},10^{-1},1,10,10^2,10^3$ and $\epsilon _{\rm FP}=10^{-3},10^{-2},10^{-1},1$ respectively. The integral approach and the FP approximation are valid respectively right and left of the related boundaries. The dotted line shows $\omega _0^*(p)$ (see Eq. (21)).
Open with DEXTER

2.1.2 Particles

Identically, the integral approach for particles is valid when the width of the scattered distribution is larger than the particle bin size. For logarithmic grids this condition is: $p_{\rm {max}}/p_{\rm {min}} > r_p^{n-1}$ where $p_{\rm {max/min}}=\sqrt{(\gamma+\omega-\omega_{\rm {min/max}})^2-1}$ are the largest and smallest momenta of the scattered particle distribution and rp=pi+1/pi is the logarithmic increment of the particle grid. By defining Rp and $\epsilon_{{\rm I},p}$ as for the photons, it yields:

\begin{displaymath}% \Delta p(p,\omega)/p_{\rm {min}}(p,\omega) > \epsilon_{\rm {I},p}. \end{displaymath} (51)

Contrary to the photon case, there is no simple equation for the boundary defining this region and the implicit Eq. (51) must be solved numerically. The boundary  $\omega _{{\rm I}}(p)$ is computed easily since there is only one positive solution to the equation on $\omega$ and it happens to be always smaller than  $\omega _0^*(p)$. A simple approximation for this boundary, accurate in all regimes to better than 30$\%$, is:

\begin{displaymath}% \omega_{\rm I}(p) \approx \frac{\epsilon_{{\rm I},p} p/2}{2... ...}-\sqrt{1+\epsilon_{{\rm I},p}})p+\epsilon_{{\rm I},p}+2}\cdot \end{displaymath} (52)

When the resolution is good enough ( $R_\omega \gg 1$, $\epsilon_{{\rm I},p} \ll 1$), condition 51 reduces to that used in Eq. (78) of Vurm & Poutanen (2009) for electron down-scattering: $\omega > \epsilon_{{\rm I},p}/(4\gamma)$ for relativistic particles.

Figure 6 shows this boundary for various resolutions ( $R_p \approx 9200,930,97,13,3.8,2.0,1.3$ and n=5). Contrary to the equation for photons, solving the evolution equation for particles with a simple numerical integration requires very high resolution to guarantee a correct accuracy. For example, soft photons of energy $\omega=10^{-5}$ scattering off high energy particles ( $\gamma =10^2$) are accurately described by the integral approach only if Rp > 900, which is far beyond current desktop computer capabilities and the constraint is even more severe when mildly relativistic particles are involved.

\begin{figure} \par\includegraphics[width=8.2cm,clip]{9113fg08.eps} \end{figure} Figure 6:

Particle equation: boundaries $\omega _{{\rm I}}(p)$ for the integral approach (solid lines) and $\omega _{{\rm FP}}(p)$ for the Fokker-Planck approximation (dashed lines) for $\epsilon _{{\rm I},p}=10^{-3},10^{-2},10^{-1},1,10,10^2,10^3$ and $\epsilon _{\rm FP}=10^{-3},10^{-2},10^{-1},1$ respectively. The integral approach and the FP approximation are valid respectively above and below the related boundaries. The dotted line shows  $\omega _0^*(p)$.
Open with DEXTER

2.1.3 Accuracy

If conditions (46) and (51) are not satisfied over the entire simulation domain, numerical computation of the Compton scattering may lead to inaccuracy. In particular, energy is not conserved when logarithmic grids are used (as explained before), and its measure provides a good way to estimate numerical errors[*]. Figure 7 shows the error on energy when soft mono-energetic photons are up-scattered by high energy particles. The numerical computation was performed with the code presented in Belmont et al. (2008), on a single time step, using an explicit scheme on logarithmic energy grids, and turning off all processes/injection/escape but Compton scattering. The contribution of Compton scattering was forced to be computed by the integral approach for the equation on particles. The contribution to the equation on photons was also computed with the integral approach but given the energy range and the grid resolution considered here, this method is accurate in this case. Both initial distributions were set to zero except in one energy bin, and they were normalised to unity to keep the number of scattering events constant as the photon energy is varied. The error was measured by computing the total energy lost by particles  $\partial_t E_p$ and comparing it to the energy gained by photons  $\partial_t E_\omega$.

\begin{figure} \par\includegraphics[width=8.2cm,clip]{9113fg09.eps} \end{figure} Figure 7:

Numerical error on energy conservation, for a single time step, as a function of the incident photon energy $\omega _0$, for particles of momentum p0=102,103,104. In the case p0=103, the boundaries between zones where the scattered distribution spans  2,3,4,5,6,7,8,9, and at least 10 energy bins are over-plotted (vertical dashed lines). The error is computed as the relative error between the total energy losses of the photon population and the total energy losses of the electron population: $\delta E /E = (\vert\partial _t E_p\vert - \vert\partial _t E_\omega \vert)/(\vert\partial _t E_p\vert + \vert\partial _t E_\omega \vert)$. For this run the resolutions are: Rp=64/6 and $R_\omega =256/12$.
Open with DEXTER

For low energy photons, the scattering takes place far below the boundary in the $(p,\omega)$ plane of Fig. 6. The scattered distribution typically spans less than 2 bins. The energy change of particles is not captured at all by the numerical scheme and the error is 100%. As the photon energy increases, the scattered distribution widens and the error decreases. From the cases presented here, we find that an error less than 1% corresponds to scattered distributions that span more than 5 bins. When n>10, the error is dominated by other weak numerical errors.

By the choice of particle and photon energies, we have only focused here on the error made in the equation for particles. A similar study can be done for photons. However, it would lead to very similar results and contrary to the case of particles, integration errors for photons are less relevant to astrophysical applications.

2.2 Fokker-Planck approach

Alternatively, Compton scattering can be described by the Fokker-Planck method. This approximation assumes that the cross section and the initial particle (or photon) distribution are only weakly dependent on the energy of the incoming particle (or photon) on the scale of the scattered distribution width. When the width of the scattered distribution is small, a second order Taylor expansion of the integral equation can be performed, which gives the well-known Fokker-Planck equations:

$\displaystyle \partial_t N_\omega = \partial_\omega \left(A_\omega N_\omega\right) + \frac{1}{2} \partial^2_{\omega^2} \left(D_\omega N_\omega\right)$ (53)
$\displaystyle \partial_t N_{\rm e} = \partial_p \left(\frac{\gamma}{p} A_\gamma... ...{\gamma}{p} \partial_p \left(\frac{\gamma}{p} D_{\rm e} N_{\rm e}\right)\right)$ (54)

with

$\displaystyle A_\omega ; D_\omega = \int_0^{\infty} c \bar{\sigma}_{1;2}(p;\omega) {\rm d}N_{\rm e}(p)$ (55)
$\displaystyle A_p ; D_p = \int_0^{\infty} c \bar{\sigma}_{1;2}(p;\omega) {\rm d} N_\omega(\omega).$ (56)

As for the integral approach, the validity of the FP approximation depends on the energy of the incoming particle and photon. Although the real criteria should in principle depend on the initial distributions and should be computed at each time step of a time dependant simulation, it is more convenient to define an approximate but general criteria.

2.2.1 Particles

The FP method for particles assumes that the quantity $N_{\rm e}(\gamma_0){\rm d}\bar{\sigma}(\gamma_0,\omega_0;\omega)$ does not vary much with $\gamma_0$ in the width $\Delta \gamma$ of the scattered distribution. There is no unique way to define such a condition. By assuming that the typical energy scale for the variation of this quantity is the initial kinetic energy: $\gamma_0-1$, then the FP approach is found to give a correct description of the particle evolution if:

\begin{displaymath}% \Delta\gamma(\omega,\gamma)/(\gamma-1) < \epsilon_{\rm FP} \end{displaymath} (57)

where $\epsilon_{\rm FP} \ll 1$ is a small parameter that describes the accuracy of the FP approximation: the smaller  $\epsilon_{\rm FP}$ the more valid the FP approach. The equation for the boundary of this region of the $(p,\omega)$ space reads:

\begin{displaymath}% \omega < \omega_{{\rm FP}}(p) \end{displaymath} (58)

where

\begin{displaymath}% \omega_{{\rm FP}} = \frac{\gamma}{2}\left( \sqrt{\frac{1-\e... ...gamma)\beta}{1-\epsilon_{\rm FP}(1-1/\gamma)/\beta}}-1\right). \end{displaymath} (59)

In the limit $\epsilon_{\rm FP} \ll 1$, a simple approximation for all particle energy is:
$\displaystyle % \omega_{{\rm FP}} \approx \frac{\epsilon_{\rm FP}\beta}{4(\gamma+1)},$     (60)

which reduces to $\omega_{{\rm FP},p} \sim \epsilon_{\rm FP}/(4\gamma)$ in the ultra-relativistic regime ($p \gg 1$). Boundaries  $\omega _{{\rm FP}}(p)$ are shown in Fig. 6 for $\epsilon_{\rm FP}=10^{-3},10^{-2},10^{-1}$ and 1. The Fokker-Planck approximation for particles is a good approximation when very low energy photons are up-scattered by mid-relativistic particles. It fails easily when the particle energy becomes large. For example, by setting $\epsilon _{\rm FP}=0.1$, the scattering of photons off particles with energy p=103 can only be described by the FP method when the photon energy is less than $\omega=10^{-5}$. However, when one does not need such a precision and sets $\epsilon_{\rm FP} \la 1$, the FP method is valid over a large range of energy.

2.2.2 Photons

Similarly, the FP method for photons assumes that the quantity $N_\omega(\omega_0){\rm d}\bar{\sigma}(\gamma_0,\omega_0;\omega)$ does not vary much with $\omega _0$ in the width  $\Delta\omega$ of the scattered distribution. By assuming that the typical energy scale for the variation of this quantity is of the order of the incoming photon energy $\omega _0$, the FP approach for photons is valid if:

\begin{displaymath}% \Delta\omega(\omega,\gamma) /\omega < \epsilon_{\rm FP}. \end{displaymath} (61)

When solved as a function of $\omega$, the equation for the boundary is again the solution of a second order polynomial and the condition 61 is satisfied in the region defined by:

\begin{displaymath}% p< p^{\rm c}_{\rm FP} \quad \rm {and} \quad \rm {min}\left[... ...),\omega_{\rm FP}^-(p)\right] < \omega < \omega_{\rm FP}^+(p) \end{displaymath} (62)

where

\begin{displaymath}% p^{\rm c}_{\rm FP} = \frac{(4-\epsilon_{\rm FP}) - (4+\epsilon_{\rm FP})\sqrt{1-\epsilon_{\rm FP}}}{4\epsilon_{\rm FP}}\cdot \end{displaymath} (63)


                           $\displaystyle % \omega^\pm_{\rm FP}$ = $\displaystyle \frac{1}{4} \left[ 2(1-p-\gamma) + \epsilon_{\rm FP}(\gamma+p) \right.$  
    $\displaystyle \pm ~ \left. \sqrt{\epsilon_{\rm FP}}\sqrt{4(1-\gamma-p)+\epsilon_{\rm FP}(1+2p\gamma+2p^2)} \right],$ (64)


$\displaystyle % \omega^{\rm c}_{\rm FP}= \frac{p}{\epsilon_{\rm FP}} ~\left(1-\... ...eta+\sqrt{(1+\epsilon_{\rm FP}/\beta)^2-\epsilon_{\rm FP}^2/p^2} } \right)\cdot$     (65)

As for the integral approach, the inverse equation for this boundary  $p_{{\rm FP}}(\omega )$ is more practical to use and is best evaluated numerically.

Boundaries $p_{\rm {FP}}$ are plotted in Fig. 5 for $\epsilon_{\rm FP}=10^{-3},10^{-2},10^{-1}$, and 1. Contrary to the equation for particles, the photon evolution is only poorly described by the FP approximation and only low energy photons ( $\omega<0.1$) scattering off low energy particles (p<0.1) undergo small angle scattering that can be accurately described by the FP approximation.

2.3 Combined approach

As can be seen in Figs. 5 and 6, the two methods happen luckily to give a correct description of Compton scattering in regions of the $(p,\omega)$ space that can be complementary. This suggests that combining the two methods is a good way to overcome numerical accuracy issues in computing the effects of Compton scattering. However, this can only be done if the validity regions for the two methods cover the entire simulation domain, that is, for given $\epsilon_{{\rm I},p}$, $\epsilon_{{\rm I},\omega}$ and $\epsilon_{\rm FP}$, if:

$\displaystyle p_{{\rm FP}}(\omega) \ge p_{{\rm I}}(\omega) \;\; \forall \omega,$ (66)
$\displaystyle \omega_{{\rm FP}}(p) \ge \omega_{{\rm I}}(p) \;\; \forall p$ (67)

for the equation on photons and particles respectively. Typically, this implies that $\epsilon_{{\rm I},\omega} \le \epsilon_{\rm FP}$ and $\epsilon_{{\rm I},p} \le \epsilon_{\rm FP}$, which sets a constraint on the grid resolution:

\begin{displaymath}% R_{\omega;p} \ge \frac{n-1}{\epsilon_{\rm {FP}}}\cdot \end{displaymath} (68)

As good accuracy for each methods requires n>5 and $\epsilon_{\rm FP} \ll 1$, this puts strong constraints on the resolution. For example, by setting $\epsilon _{\rm FP}=0.1$ and n=5 this leads to R>40. For particle grids that typically span 5 orders of magnitude, it implies the need to use 200-bin grids accessible to most desktop computers. However, for photon grids that typically span 15 orders of magnitude, it corresponds to 600-bin grids, which requires high computing power.

For each equation (for photons and for particles), one can define an average boundary when the validity regions of both methods overlap:

\begin{displaymath}% \omega_{\rm c}(p) = \left( \omega_{\rm {I}} \omega_{{\rm FP... ...\rm c}(\omega) = \left(p_{\rm {I}} p_{{\rm FP}} \right)^{1/2}. \end{displaymath} (69)

Finally, a combination of the two methods is achieved by solving the following equations:
                       $\displaystyle % \partial_t N_{\rm e}(p)$ = $\displaystyle \partial_p \left( \frac{\gamma}{p} A^{\omega<\omega_{\rm c}}_{\rm e} N_{\rm e} \right)$  
    $\displaystyle + ~ \frac{1}{2}\partial_p \left( \frac{\gamma}{p} \partial_p \left(\frac{\gamma}{p}D_{\rm e}^{\omega<\omega_{\rm c}} N_{\rm e} \right) \right)$  
    $\displaystyle + ~ \int {\rm d}N_{\rm e}(p_0) \int_{\omega_{\rm c}(p_0)}^{\infty... ...rm d}\bar{\sigma}}{{\rm d}p}(p_0,\omega_0; \omega(p)) {\rm d}N_\omega(\omega_0)$  
    $\displaystyle - ~ N_{\rm e}(p) \int_{\omega_{\rm c}(p)}^{\infty} c \bar{\sigma}_0(\omega_0,p){\rm d}N_\omega(\omega_0),$ (70)


                  $\displaystyle % \partial_t N_\omega(\omega)$ = $\displaystyle \partial_\omega\left( A^{p<p_{\rm c}}_\omega N_\omega \right) +\frac{1}{2} \partial^2_{\omega^2} \left(D_\omega^{p<p_{\rm c}} N_\omega \right)$ (71)
    $\displaystyle + ~ \int {\rm d}N_\omega(\omega_0) \int_{p_{\rm c}(\omega_0)}^{\i... ...\rm d}\bar{\sigma}}{{\rm d}\omega}(p_0,\omega_0 ; \omega) {\rm d}N_{\rm e}(p_0)$  
    $\displaystyle - ~ N_\omega(\omega) \int_{p_{\rm c}(\omega)}^{\infty} c \bar{\sigma}_0(\omega,p_0) {\rm d}N_{\rm e}(p_0)$ (72)

where integrations in the integral part are performed only above the boundary energy  $\omega_{\rm c}$ and momentum $p_{\rm c}$ and where the FP coefficients are computed as integrals of the cross section only below these boundaries:

$\displaystyle A_\omega ; D_\omega = \int_0^{p_{\rm c}(\omega)}c \sigma_{1;2}(p;\omega) d N_{\rm e}(p)$ (73)
$\displaystyle A_p ; D_p = \int_0^{\omega_{\rm c}(p)} c \sigma_{1;2}(p;\omega) d N_\omega(\omega).$ (74)

Figure 8 shows the steady photon and particle distributions for a more realistic simulation. Particles are injected with a mono-energetic distribution at high energy ( $\gamma =10^2$) at a constant rate (the compactness of which is set to $l_{\rm in,e}=\sigma_T P /(m_{\rm e}c^3R)=10$, where P is the injected power) and escape at a constant rate with the speed of light to allow for a steady state. Soft photons are also injected with a black body spectrum of temperature $k_{\rm B}T=3$ $\times $ $10^{-4} m_{\rm e}c^2$ and a flux of compactness $l_{{\rm in},\omega}=\sigma_T L_{\rm soft} /(m_{\rm e}c^3R)=1$ (see e.g. Belmont et al. 2008, for detailed definitions of the compactness parameters). We start with an empty system and the code simultaneously evolves both distributions with time until they reach a steady state. The only process turned on is Compton scattering. Figure 8 shows the results of three kinds of runs where the contribution of Compton scattering to the evolution of the particle distribution was forced to be treated with the integral (several resolutions), the Fokker-Planck and the combined approaches. The distributions look very different.

\begin{figure} \par\hspace*{-1mm}\includegraphics[width=8.55cm,clip]{9113fg10.eps}\vspace*{3mm} \includegraphics[width=8.5cm,clip]{9113fg11.eps} \end{figure} Figure 8:

Particle distribution and spectra in steady state when Compton scattering for particles is treated with the integral method for different resolutions: Rp=512/6, 256/6, 128/6, and 64/6 (solid lines), when it is treated with the combined approach with $\epsilon _{\rm FP}=0.1$ and n=10 (dashed line), and with the Fokker-Planck approximation (dotted line). For these runs, the power injected in particles and photons is: $l_{\rm in, e}=10$, $l_{{\rm in},\omega }=1$ and the source size is set to R=107 cm. The temperature of the seed photons is $k_{\rm B}T/m_{\rm e}c^2=3$ $\times $ 10-4 and the energy of the injected particles is $\gamma =10^2$. For the combined and FP approaches, the particle resolution is: Rp=64/6. For all runs, the resolution of the photon grid is $R_\omega =512/9$.
Open with DEXTER

As expected, the integral approach fails to capture efficiently the particle cooling for low resolution runs: the lower the resolution, the steeper the slope of the high energy tail of the particle distribution. These excess high energy particles more efficiently up-scatter the soft photons and the resulting spectra are harder. In steady state, the total energy loss must balance the total injected power. Power is supplied through particles $l_{\rm in, e}=10$ and through photons: $l_{\rm in}=l_{\rm in,e}+l_{{\rm in},\omega}=11$. Energy is also lost both through particles and photons. The steady particle distribution has low density $\tau=\sigma_TR\int {\rm d}N_{\rm e} =2.3$ $\times $ 10-2 so that the energy lost through particle escapes is rather small ( $l_{\rm out,e}\approx 1$) and most of the energy escape as radiation. As numerical errors prevent particles from cooling as they should, the overall system actually gains energy artificially and the total energy loses are measured to be larger than the injected power: $l_{\rm out} = l_{\rm out,\omega}+l_{\rm out,e} =19.6, 14.2, 12.2$, and 11.4 from lower to higher grid resolution. The corresponding effective electron temperatures are $k_{\rm B}T_{\rm e}=6.23, 4.68, 3.95$, and $3.60~m_{\rm e}c^2$ (see Eq. (2.8), Coppi 1992).

The Fokker-Planck equation provides good energy conservation and the plasma luminosity balances the injected power: $l_{\rm out}=11.0$. However, energy conservation does not guarantee accurate results and significant deviation is observed in the low energy part of the particle distribution, where the Fokker-Planck approximation is expected to fail. In the particular case presented here, this has no effect on the emitted spectrum since the spectrum is dominated by the Comptonisation by high energy particles. Nevertheless it could strongly affect the predicted spectrum when other processes such as Coulomb collisions or pair annihilation are involved. The electron temperature is measured to be: $k_{\rm B}T_{\rm e}=3.40~m_{\rm e}c^2$.

The combined approach produces the best results. With low grid resolution (Rp=64/6), it gives the correct photon spectrum and the correct particle distribution, it conserves energy to better than $1\%$ accuracy ( $l_{\rm out}=11.0$) and it was checked that the shape of the steady distributions does not depend on the grid resolution (not shown here). The electron temperature is $k_{\rm B}T_{\rm e}=3.39~m_{\rm e}c^2$. In comparison, the integral approach requires more than 10 times higher resolution to produce similar results, which also means at least 10 times more computational time. Extensive tests show that the combined approach is very efficient in most cases. It takes advantage of the two approaches. The way the two methods are associated is numerically simple to implement and it can easily be shown that the total number of particles and photons as well as the total energy are conserved. Moreover even when conditions (66) and (67) are not exactly satisfied this method still provides good accuracy in many cases since, as long as these conditions are not strongly violated, only a small faction of particles and/or photons is treated inaccurately. Also, even if the evolution of a large fraction of particles and/or photons is not described correctly by the numerical scheme they might have only a minor contribution to the overall distribution evolution in some specific cases. For these reasons some simulations might reach a correct final accuracy, despite violating the accuracy condition. However, this is very problem-dependent and the accuracy should be checked very carefully when these conditions are not satisfied.

Conclusion

In this paper we have investigated numerical issues related to the computation of isotropic Compton scattering in numerical codes. We have derived a form of the distribution of Compton scattered photons that is free of numerical cancellations. This form is exact with no limitation on the photon and electron energies. The numerical accuracy resulting from the direct use of the cross section has been studied for various grid resolutions and equations have been proposed for the boundaries where this method must be replaced by an approximate Fokker-Planck treatment to guarantee sufficient accuracy. All results presented here have been included in and extensively tested with the new code developed by Belmont et al. (2008).

Acknowledgements
The author thanks J. Poutanen for providing the numerical routines corresponding to the formulae of NP94.

References

  • Barbosa, D. D. 1982, ApJ, 254, 301 [NASA ADS] [CrossRef]
  • Belmont, R., Malzac, J., & Marcowith, A. 2008, A&A, 491, 617 [NASA ADS] [CrossRef] [EDP Sciences]
  • Blumenthal, G. R., & Gould, R. J. 1970, Rev. Mod. Phys., 42, 237 [NASA ADS] [CrossRef]
  • Brinkmann, W. 1984, J. Quant. Spectrosc. Radiat. Trans., 31, 417 [NASA ADS] [CrossRef]
  • Coppi, P. S. 1992, MNRAS, 258, 657 [NASA ADS]
  • Coppi, P. S., & Blandford, R. D. 1990, MNRAS, 245, 453 [NASA ADS]
  • Donahue, T. M. 1951, Phys. Rev., 84, 972 [NASA ADS] [CrossRef]
  • Feenberg, E., & Primakoff, H. 1948, Phys. Rev., 73, 449 [NASA ADS] [CrossRef]
  • Follin, J. W. 1947, Phys. Rev., 72, 743 [NASA ADS]
  • Heitler, W. 1954, International Series of Monographs on Physics (Oxford: Clarendon), 3rd edn.
  • Jones, F. C. 1968, Phys. Rev., 167, 1159 [NASA ADS] [CrossRef]
  • Katz, J. I. 1987, Frontiers in Physics (Reading: Menlo Park, Addison-Wesley)
  • Kompaneets, A. S. 1957; Sov. Phys., JETP 4, 730
  • Levich, E. V., & Syunyaev, R. A. 1971, Sov. Astron., 15, 363 [NASA ADS]
  • Nagirner, D. I., & Poutanen, J. 1994, Single Compton scattering, ed. D. I. Nagirner, & J. Poutanen, Astrophys. Space Phys. Rev. (Amsterdam: Harwood Academic Publishers), part 1, 9, 83 (NP94)
  • Nayakshin, S., & Melia, F. 1998, ApJS, 114, 269 [NASA ADS] [CrossRef]
  • Pe'er, A., & Waxman, E. 2005, ApJ, 628, 857 [NASA ADS] [CrossRef]
  • Poutanen, J., & Svensson, R. 1996, ApJ, 470, 249 [NASA ADS] [CrossRef]
  • Rybicki, G. B., & Lightman, A. P. 1986, Radiative Processes in Astrophysics, ed. G. B. Rybicki, & A. P. Lightman (Wiley-VCH), 400
  • Saugé, L., & Henri, G. 2004, ApJ, 616, 136 [NASA ADS] [CrossRef]
  • Tavecchio, F., Maraschi, L., Pian, E., et al. 2001, ApJ, 554, 725 [NASA ADS] [CrossRef]
  • Vurm, I., & Poutanen, J. 2009, ApJ, 698, 293 [NASA ADS] [CrossRef]
  • Zeldovich, I. B. 1975, Sov. Phys. Uspekhi, 115, 161 [NASA ADS]
  • Zeldovich, Y. B., Illarionov, A. F., & Syunyaev, R. A. 1972, Zhurnal Eksperimental noi i Teoreticheskoi Fiziki, 62, 1217 [NASA ADS]
 

Footnotes

... misprints[*]
Contrary to what is claimed in Coppi & Blandford (1990), all intermediate expressions are correct and only the final Eqs. (23)-(27) contains misprints.
... NP94[*]
The equations in NP94 contain some misprints too. In particular: Eq. (8.1.10) should read $\mu_+ = 1- (x-x_1)^2/D_m/x/x_1$ and Eq. (8.1.22) should read $H_1 = \left[ H_0-\Delta h A_1(h_-) + H/A(h_-)A(h_+)\right]/h_+ $
... equations[*]
A few misprints were corrected. In Eq. (3.3.5), $\Psi_{14}$ should read $(2\psi_{-11}-\psi_{-10}-\psi_{-12})/\xi$, the last expression for Eq. (3.3.9) should be $\pi^2/6+\ln^2{(2\xi)}/2 - g_*(1/2\xi)$, and in Eq. (3.3.10), $\psi_{00}$ should read: $3/(8\xi)\left[ g(\xi)-2/\xi+ (1/2+2/\xi+1/\xi^2)l_\xi-R_\xi/2-3/2\right]$.
... lose any[*]
We note that the energy is not conserved only with non-linear grids. For linear grids and identical resolution in the photon and electron energy grids ($\omega _0$ and $\gamma_0$ respectively), the energy is balanced to machine precision. In such case indeed, the bin size is uniform. As the energy width of the scattered distribution is by definition the same for both species, if the scattering of particles to neighbour energy bins is not captured by the numerical integration, the scattering of photon is not either, so that the integral approach fails simultaneously for both species and there is no net energy exchange between the two species. Although the integral approach is clearly not accurate in this case, the energy balance is computed to machine precision.
... read[*]
The -1 comes from the fact that when the width of the scattered distribution is exactly $n\delta\omega$, it is discretized in n+1 bins in all situations except when it is exactly centred at the centre of a bin.
... errors[*]
For linear grids, other accuracy indicators should be used.
Copyright ESO 2009

All Figures

  \begin{figure} \par\includegraphics[width=8.5cm,clip]{9113fg01.eps} \end{figure} Figure 1:

Angles and notations. The incoming photon direction $\vec{n}_0$ is defined by the incident angle $\theta _0$ and an azimuthal angle $\psi _0$ that orientates the global picture (not shown here). The scattered photon direction $\vec{n}$ is defined by the scattering angle and the azimuthal angle $\phi $. In coordinates $(\vec{e}_1,\vec{e}_2,\vec{e}_3)$, $\vec{\beta}_0/\beta_0=(\mu_0,\sin{\theta_0},0)$ and $\vec{n}=(c_\alpha,\sin{\alpha}\cos{\phi},\sin{\phi})$, so that $\mu = \mu_0 c_\alpha+ \sin{\alpha}\sin{\mu_0}\cos{\phi}$.
Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=8.5cm,clip]{9113fg02.eps} \end{figure} Figure 2:

Accuracy domains. Equation (11) is accurate for photons energy and electron momentum  $(\omega _0,p_0)$ in the region delimited by the solid line. The original formula derived by Jones (1968), and the cross section by Pe'er & Waxman (2005), have basically the same properties and are accurately evaluated in the same region. The published formula by NP94 is accurate in the domain determined by the dashed line. These boundaries are for double precision computing. The star shows the point for which the scattered distribution is plotted in Fig. 4.
Open with DEXTER
In the text

  \begin{figure} \par\hspace*{-1mm}\includegraphics[width=8.35cm,clip]{9113fg03.ep... ....eps}\vspace*{3mm} \includegraphics[width=8.3cm,clip]{9113fg05.eps} \end{figure} Figure 3:

First three moments of the scattered distribution as a function of the energy of the incoming photon: the total cross section, the mean energy of the scattered photon and the dispersion about this mean value. The solid and dashed lines show the numerical integration of the expression derived in this paper and analytical results given in NP94 respectively. The 3 curves are for p0=10-2,102 and 104. For each panel, relative errors are also shown.
Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=8.2cm,clip]{9113fg06.eps} \end{figure} Figure 4:

The scattered distribution for soft photons ( $h\nu =10^{-5}m_{\rm e}c^2$) and mild relativistic particles (p=1). Solid and dotted lines show the distributions obtained with Eq. (28), and with Jones' formula respectively.
Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=8.2cm,clip]{9113fg07.eps} \end{figure} Figure 5:

Photon equation: boundaries $p_{{\rm I}}(\omega )$ for the integral approach (solid lines) and $p_{{\rm FP}}(\omega )$ for the Fokker-Planck approximation (dashed lines) for $\epsilon _{{\rm I},\omega }=10^{-3},10^{-2},10^{-1},1,10,10^2,10^3$ and $\epsilon _{\rm FP}=10^{-3},10^{-2},10^{-1},1$ respectively. The integral approach and the FP approximation are valid respectively right and left of the related boundaries. The dotted line shows $\omega _0^*(p)$ (see Eq. (21)).
Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=8.2cm,clip]{9113fg08.eps} \end{figure} Figure 6:

Particle equation: boundaries $\omega _{{\rm I}}(p)$ for the integral approach (solid lines) and $\omega _{{\rm FP}}(p)$ for the Fokker-Planck approximation (dashed lines) for $\epsilon _{{\rm I},p}=10^{-3},10^{-2},10^{-1},1,10,10^2,10^3$ and $\epsilon _{\rm FP}=10^{-3},10^{-2},10^{-1},1$ respectively. The integral approach and the FP approximation are valid respectively above and below the related boundaries. The dotted line shows  $\omega _0^*(p)$.
Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=8.2cm,clip]{9113fg09.eps} \end{figure} Figure 7:

Numerical error on energy conservation, for a single time step, as a function of the incident photon energy $\omega _0$, for particles of momentum p0=102,103,104. In the case p0=103, the boundaries between zones where the scattered distribution spans  2,3,4,5,6,7,8,9, and at least 10 energy bins are over-plotted (vertical dashed lines). The error is computed as the relative error between the total energy losses of the photon population and the total energy losses of the electron population: $\delta E /E = (\vert\partial _t E_p\vert - \vert\partial _t E_\omega \vert)/(\vert\partial _t E_p\vert + \vert\partial _t E_\omega \vert)$. For this run the resolutions are: Rp=64/6 and $R_\omega =256/12$.
Open with DEXTER
In the text

  \begin{figure} \par\hspace*{-1mm}\includegraphics[width=8.55cm,clip]{9113fg10.eps}\vspace*{3mm} \includegraphics[width=8.5cm,clip]{9113fg11.eps} \end{figure} Figure 8:

Particle distribution and spectra in steady state when Compton scattering for particles is treated with the integral method for different resolutions: Rp=512/6, 256/6, 128/6, and 64/6 (solid lines), when it is treated with the combined approach with $\epsilon _{\rm FP}=0.1$ and n=10 (dashed line), and with the Fokker-Planck approximation (dotted line). For these runs, the power injected in particles and photons is: $l_{\rm in, e}=10$, $l_{{\rm in},\omega }=1$ and the source size is set to R=107 cm. The temperature of the seed photons is $k_{\rm B}T/m_{\rm e}c^2=3$ $\times $ 10-4 and the energy of the injected particles is $\gamma =10^2$. For the combined and FP approaches, the particle resolution is: Rp=64/6. For all runs, the resolution of the photon grid is $R_\omega =512/9$.
Open with DEXTER
In the text

Copyright ESO 2009

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.