A&A 506, 589599 (2009)
Numerical computation of isotropic Compton scattering
R. Belmont^{1,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 FokkerPlanck
approximation. Expressions describing this combination are given.
Key words: radiation mechanisms: general  plasmas  methods: numerical  galaxies: active  Xrays: binaries  Xrays: 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 Xray or relativistic sources (Levich & Syunyaev 1971). It is also responsible for the upscattering 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 ray bursts, AGN, and Xray binaries. Particular focus is now on the synchrotron selfCompton 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 Xray binaries ( 100 keV, i.e. ), and to high energy cosmic rays (up to E > 100 TeV, i.e. ). The energy of photons also spans a wide range: in the case of AGN for instance, it goes from low energy synchrotron photons ( m, i.e. ) to ray emission (with E > 10 TeV, i.e. 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 upscattering 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 photonphoton, photonparticle and particleparticle 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 FokkerPlanck 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 Comptonscattered 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):
=  
(1) 
Here and are the photon and lepton distributions in their respective 6dimension momentum space. The first term in the brackets describes the contribution of all photons of momentum being scattered to momentum after one interaction. The second one describes photons of momentum being scattered to any other momentum. Both are proportional to the Compton differential cross section giving the probability of one photon of momentum being scattered to a momentum by an electron of momentum . This cross section can be interpreted as the distribution of scattered photons resulting from the interaction of monoenergetic 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
and
(where
), it yields:
where and are the angleintegrated distributions. Here
(3) 
is the total scattering cross section, and
is the angleaveraged cross section for isotropic Compton scattering.
Figure 1: Angles and notations. The incoming photon direction is defined by the incident angle and an azimuthal angle that orientates the global picture (not shown here). The scattered photon direction is defined by the scattering angle and the azimuthal angle . In coordinates , and , so that . 

Open with DEXTER 
Equation (4) can actually be simplified further. Because of symmetry, the differential cross section
does
not depend independently on the direction of the incoming photon and
lepton, but only on the angle between their directions:
,
where
is the lepton velocity in units of the speed of light and
is the direction of the incoming photon (see Fig. 1). Then, Eq. (4) reduces to:
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 energymomentum 4vector 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 : the cosine of the scattering angle , and an azimuthal angle (see Fig. 1). Then, the differential cross section is redundant and includes a deltafunction that ensures that the total energymomentum is conserved:
Here the latter condition is given by:
(7) 
where 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 KleinNishina
cross section does not depend on the azimuthal angle and reads (Heitler 1954; Rybicki & Lighman 1979):
where is the Thomson cross section, 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
angleaveraged cross section is now quite straightforward. However,
including the KleinNishina 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:
Integration is then performed over , and then . The latter integration is made with the variable . 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):
(10) 
with
where and are the Lorentz factors of the incoming and scattered particles respectively and is the momentum of the scattered particle. Here the following notations have been used:
(12)  
(13)  
(14)  
(15) 
The cross section is proportional to the difference of function F evaluated at the upper and lower boundaries of the last integration (b and a respectively). At these boundaries, the variables have the following expressions:
(16)  
(17) 
The critical photon frequencies for which the arguments of the and functions are equal define distinct domains of the scattered distribution. For the lower integration boundary a, both arguments are equal at a unique frequency: . For the higher boundary, both arguments are equal at two critical frequencies: and , where
(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):
where
When, , the photon can gain up to the entire electron energy. However, this absolute limit is not necessarily reached. When (as for Compton upscattering 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 roundoff 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 .
We have checked the accuracy of the scattered distribution for a large range of photon and lepton energies and identify the domain of the plane where Eq. (11) gives accurate results. Results are shown in Fig. 2. We find that Compton scattering of mild photons ( ) by midrelativistic particles () 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 upscattering of soft photons by high energy particles.
In the latter case, Jones (1968) derived asymptotic expansions in 1. However, a large number of orders must be kept when 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 Rewriting 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  F F, F
This happens typically when , as for example, at the distribution upper and lower photon energies and where the distribution vanishes, but where F^{a} and F^{b} can remain large.
Differences of hyperbolic functions between the upper and lower boundaries can be computed analytically by using the following trigonometric relations:
The difference of the two hyperbolic functions between the two integration boundaries a and b simply reads: if , or if , where
=  L_{+}(x_{+}^{b})L_{+}(x_{+}^{a})  (22)  
=  L_{}(x_{}^{b})L_{}(x_{}^{a})  (23)  
=  (24) 
Although the last term of Eq. (11) is not divergent it is numerically illbehaved when and it is best to use the following auxiliary function:
S(x) =  (25) 
which is evaluated as a series expansion for small argument:
Differences of all other terms in Eq. (11) between the two integration boundaries can also be computed analytically to factorize and get the following quite simple expression for the scattered distribution:
with
where and . The scattered distribution is proportional to . This factor holds as the difference between the two integration boundaries a and b and vanishes at the minimal and maximal photon energies and . For high particle energy, it is best computed as:
=  
(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 is not accurate enough and an alternative derivation must used:
(30) 
Equations (27)(28) are accurate over a larger domain than the original formula and can be used in many astrophysical situations.
Figure 2: Accuracy domains. Equation (11) is accurate for photons energy and electron momentum 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_{} z_{+}, z_{}
However, this new expression reveals other differences hidden in the original one and that result in cancelation issues, namely when . 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 .
By using the definition of z_{+} and z_{}, the differences in the first three terms of Eq. (28) can be computed as:
(31) 
(32)  
(33) 
respectively, where
=  (34)  
=  (35)  
(36) 
In the same manner, the last term can be computed as:
(37) 
where
=  (38)  
(39) 
The numerical cancelation of fails either when or when and these two cases must be dealt with separately. In the former case, first order Taylor developments of S around the mid value are sufficient whereas in the latter case (typically when ), series expansions of S in zero can be used (Eq. (26)). In both cases, the difference is proportional to which must be computed as:
(40) 
for low photon energy ( ). 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 upscattering 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:
(41)  
(42)  
(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.
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 p_{0}=10^{2},10^{2} and 10^{4}. 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 summedup in Fig. 2, the original cross section by Jones (1968) is inaccurate in typical astrophysical situations involving low energy photons upscattered 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 KleinNishina 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.
Figure 4: The scattered distribution for soft photons ( ) 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:
where is the energy of the scattered photon and is the differential cross section of the interaction . 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 FokkerPlanck 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 Xray binaries: the soft disc photons are upscattered 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 p_{0} of the incoming particle and the energy 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 energydependent. The accuracy of the numerical integration thus strongly depends on the region in the lepton momentumphoton energy space 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
is much larger than the local photon bin size
.
For linear grids, requiring that the scattered photon distribution spans at least n bins would simply read^{}:
.
For logarithmic grids, the same condition reads:
where
is the logarithmic increment of the photon grid. By denoting
the number of bins per decade energy, i.e. the resolution of the photon grid, and
,
the condition for the integral approach to be valid reads:
When the resolution is large enough: . However, as n is typically of several, 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 space. Using Eqs. (19) and (20) to compute the scattered distribution width, one can write an equation in 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:
(47) 
with
(48) 
=  
(49) 
(50) 
The inverse equation for this boundary would be more practical since it would be monovalued 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 and n=5. For comparison typical runs have a resolution . In the case = 13, only the scattering of low energy photons ( keV) off subrelativistic 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.
Figure 5: Photon equation: boundaries for the integral approach (solid lines) and for the FokkerPlanck approximation (dashed lines) for and respectively. The integral approach and the FP approximation are valid respectively right and left of the related boundaries. The dotted line shows (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:
where
are the largest and smallest momenta of the scattered particle distribution and
r_{p}=p^{i+1}/p^{i} is the logarithmic increment of the particle grid. By defining R_{p} and
as for the photons, it yields:
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 is computed easily since there is only one positive solution to the equation on and it happens to be always smaller than . A simple approximation for this boundary, accurate in all regimes to better than 30, is:
(52) 
When the resolution is good enough ( , ), condition 51 reduces to that used in Eq. (78) of Vurm & Poutanen (2009) for electron downscattering: for relativistic particles.
Figure 6 shows this boundary for various resolutions ( 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 scattering off high energy particles ( ) are accurately described by the integral approach only if R_{p} > 900, which is far beyond current desktop computer capabilities and the constraint is even more severe when mildly relativistic particles are involved.
Figure 6: Particle equation: boundaries for the integral approach (solid lines) and for the FokkerPlanck approximation (dashed lines) for and respectively. The integral approach and the FP approximation are valid respectively above and below the related boundaries. The dotted line shows . 

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 monoenergetic photons are upscattered 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 and comparing it to the energy gained by photons .
Figure 7: Numerical error on energy conservation, for a single time step, as a function of the incident photon energy , for particles of momentum p_{0}=10^{2},10^{3},10^{4}. In the case p_{0}=10^{3}, the boundaries between zones where the scattered distribution spans 2,3,4,5,6,7,8,9, and at least 10 energy bins are overplotted (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: . For this run the resolutions are: R_{p}=64/6 and . 

Open with DEXTER 
For low energy photons, the scattering takes place far below the boundary in the 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 FokkerPlanck approach
Alternatively, Compton scattering can be described by the FokkerPlanck 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 wellknown FokkerPlanck equations:
(53)  
(54) 
with
(55)  
(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
does not vary much with
in the width
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:
,
then the FP approach is found to give a correct description of the particle evolution if:
where is a small parameter that describes the accuracy of the FP approximation: the smaller the more valid the FP approach. The equation for the boundary of this region of the space reads:
(58) 
where
(59) 
In the limit , a simple approximation for all particle energy is:
(60) 
which reduces to in the ultrarelativistic regime (). Boundaries are shown in Fig. 6 for and 1. The FokkerPlanck approximation for particles is a good approximation when very low energy photons are upscattered by midrelativistic particles. It fails easily when the particle energy becomes large. For example, by setting , the scattering of photons off particles with energy p=10^{3} can only be described by the FP method when the photon energy is less than . However, when one does not need such a precision and sets , 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
does not vary much with
in the width
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 ,
the FP approach for photons is valid if:
When solved as a function of , 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:
(62) 
where
(63) 
=  
(64) 
(65) 
As for the integral approach, the inverse equation for this boundary is more practical to use and is best evaluated numerically.
Boundaries are plotted in Fig. 5 for , and 1. Contrary to the equation for particles, the photon evolution is only poorly described by the FP approximation and only low energy photons ( ) 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 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 , and , if:
for the equation on photons and particles respectively. Typically, this implies that and , which sets a constraint on the grid resolution:
(68) 
As good accuracy for each methods requires n>5 and , this puts strong constraints on the resolution. For example, by setting and n=5 this leads to R>40. For particle grids that typically span 5 orders of magnitude, it implies the need to use 200bin grids accessible to most desktop computers. However, for photon grids that typically span 15 orders of magnitude, it corresponds to 600bin 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:
(69) 
Finally, a combination of the two methods is achieved by solving the following equations:
=  
(70) 
=  (71)  
(72) 
where integrations in the integral part are performed only above the boundary energy and momentum and where the FP coefficients are computed as integrals of the cross section only below these boundaries:
(73)  
(74) 
Figure 8 shows the steady photon and particle distributions for a more realistic simulation. Particles are injected with a monoenergetic distribution at high energy ( ) at a constant rate (the compactness of which is set to , 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 and a flux of compactness (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 FokkerPlanck and the combined approaches. The distributions look very different.
Figure 8: Particle distribution and spectra in steady state when Compton scattering for particles is treated with the integral method for different resolutions: R_{p}=512/6, 256/6, 128/6, and 64/6 (solid lines), when it is treated with the combined approach with and n=10 (dashed line), and with the FokkerPlanck approximation (dotted line). For these runs, the power injected in particles and photons is: , and the source size is set to R=10^{7} cm. The temperature of the seed photons is 10^{4} and the energy of the injected particles is . For the combined and FP approaches, the particle resolution is: R_{p}=64/6. For all runs, the resolution of the photon grid is . 

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 upscatter 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 and through photons: . Energy is also lost both through particles and photons. The steady particle distribution has low density 10^{2} so that the energy lost through particle escapes is rather small ( ) 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: , and 11.4 from lower to higher grid resolution. The corresponding effective electron temperatures are , and (see Eq. (2.8), Coppi 1992).
The FokkerPlanck equation provides good energy conservation and the plasma luminosity balances the injected power: . However, energy conservation does not guarantee accurate results and significant deviation is observed in the low energy part of the particle distribution, where the FokkerPlanck 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: .
The combined approach produces the best results. With low grid resolution (R_{p}=64/6), it gives the correct photon spectrum and the correct particle distribution, it conserves energy to better than accuracy ( ) 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 . 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 problemdependent 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 FokkerPlanck 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). AcknowledgementsThe 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, AddisonWesley)
 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 (WileyVCH), 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 and Eq. (8.1.22) should read
 ... equations^{}
 A few misprints were corrected. In Eq. (3.3.5), should read , the last expression for Eq. (3.3.9) should be , and in Eq. (3.3.10), should read: .
 ... lose any^{}
 We note that the energy is not conserved only with nonlinear grids. For linear grids and identical resolution in the photon and electron energy grids ( and 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 , 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.
All Figures
Figure 1: Angles and notations. The incoming photon direction is defined by the incident angle and an azimuthal angle that orientates the global picture (not shown here). The scattered photon direction is defined by the scattering angle and the azimuthal angle . In coordinates , and , so that . 

Open with DEXTER  
In the text 
Figure 2: Accuracy domains. Equation (11) is accurate for photons energy and electron momentum 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 
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 p_{0}=10^{2},10^{2} and 10^{4}. For each panel, relative errors are also shown. 

Open with DEXTER  
In the text 
Figure 4: The scattered distribution for soft photons ( ) 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 
Figure 5: Photon equation: boundaries for the integral approach (solid lines) and for the FokkerPlanck approximation (dashed lines) for and respectively. The integral approach and the FP approximation are valid respectively right and left of the related boundaries. The dotted line shows (see Eq. (21)). 

Open with DEXTER  
In the text 
Figure 6: Particle equation: boundaries for the integral approach (solid lines) and for the FokkerPlanck approximation (dashed lines) for and respectively. The integral approach and the FP approximation are valid respectively above and below the related boundaries. The dotted line shows . 

Open with DEXTER  
In the text 
Figure 7: Numerical error on energy conservation, for a single time step, as a function of the incident photon energy , for particles of momentum p_{0}=10^{2},10^{3},10^{4}. In the case p_{0}=10^{3}, the boundaries between zones where the scattered distribution spans 2,3,4,5,6,7,8,9, and at least 10 energy bins are overplotted (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: . For this run the resolutions are: R_{p}=64/6 and . 

Open with DEXTER  
In the text 
Figure 8: Particle distribution and spectra in steady state when Compton scattering for particles is treated with the integral method for different resolutions: R_{p}=512/6, 256/6, 128/6, and 64/6 (solid lines), when it is treated with the combined approach with and n=10 (dashed line), and with the FokkerPlanck approximation (dotted line). For these runs, the power injected in particles and photons is: , and the source size is set to R=10^{7} cm. The temperature of the seed photons is 10^{4} and the energy of the injected particles is . For the combined and FP approaches, the particle resolution is: R_{p}=64/6. For all runs, the resolution of the photon grid is . 

Open with DEXTER  
In the text 
Copyright ESO 2009