NonMaxwellian electron distributions in clusters of galaxies
J. S. Kaastra^{1,2}  A. M. Bykov^{3}  N. Werner^{4}
1  SRON Netherlands Institute for Space Research, Sorbonnelaan 2, 3584 CA Utrecht, The Netherlands
2  Sterrenkundig Instituut, Universiteit Utrecht, PO Box 80000, 3508 TA Utrecht, The Netherlands
3  A.F. Ioffe Institute of Physics and
Technology, St. Petersburg 194021, Russia
4  Kavli Institute for Particle Astrophysics and Cosmology,
Stanford University, 452 Lomita Mall/mc 4085, Stanford, CA 94305, USA
Received 14 May 2009 / Accepted 28 May 2009
Abstract
Context. Thermal Xray spectra of clusters of galaxies and other sources are commonly calculated assuming Maxwellian electron distributions. There are situations where this approximation is not valid, for instance near interfaces of hot and cold gas and near shocks.
Aims. The presence of nonthermal electrons affects the Xray spectrum. To study the role of these electrons in clusters and other environments, an efficient algorithm is needed to calculate the Xray spectra.
Methods. We approximate an arbitrary electron distribution by the sum of Maxwellian components. The decomposition uses either a genetic algorithm or an analytical approximation. The Xray spectrum is then evaluated with a linear combination of those Maxwellian components.
Results. Our method is fast and leads to an accurate evaluation of the spectrum. The Maxwellian components allow the use of the standard collisional rates that are available in plasma codes such as SPEX. We give an example of a spectrum for the suprathermal electron distribution behind a shock in a cluster of galaxies. The relative intensities of the satellite lines in such a spectrum are sensitive to suprathermal electrons. These lines can only be investigated with high spectral resolution. We show that the instruments on future missions like AstroH and IXO will be able to demonstrate the presence or absence of these suprathermal electrons.
Key words: acceleration of particles  radiation mechanisms: general  Xrays: galaxies: clusters  Xrays: general
1 Introduction
Most of the visible baryonic matter in clusters of galaxies is in the form of a hot, tenuous gas. Usually it is assumed that this gas is in collisional ionisation equilibrium. In the central parts, the density is high enough and the timescales long enough to fulfil the conditions for collisional ionisation equilibrium. In the outer parts this may not always be the case. Freshly accreted gas may still be ionising. For a review of equilibration processes in such tenuous plasmas, see Bykov et al. (2008). Nevertheless, usually a Maxwellian electron distribution is still assumed to be valid; however, there are situations where this is not an obvious constraint. For instance, in cooling cores of clusters, electrons from hot gas may penetrate the colder parts. Also, when shocks are present, deviations from a Maxwellian distribution may occur associated with the temperature gradients or caused by particles accelerated by the shock. Shocks can stem from merger activity, AGN jets, and accretion shocks may occur in cluster outskirts. As the Coulomb thermal relaxation times increases with electron energy E as E^{3/2}, suprathermal electrons are relatively longlived and may yield a pressure contribution that is potentially interesting for mass profiles. The presence of such suprathermal electrons can be revealed by excess emission in satellite lines.
In this paper we describe a way to model the emerging Xray emission spectra for the case of a nonMaxwellian electron distribution. Although our main focus here is on clusters of galaxies, the same procedure can be applied to different circumstances, for instance in stellar coronae, supernova remnants, and the Galactic ridge; for the last case see also Masai et al. (2002) and Tatischeff (2003).
2 Xray spectra for suprathermal electron distributions
To calculate an Xray spectrum, essentially two steps are needed: 1) determine the ionisation balance and; 2) calculate the corresponding spectrum. For more details see e.g. Kaastra et al. (2008).
We first consider the ionisation balance. To determine that, it is necessary to evaluate the collisional ionisation and recombination rates for all ions. These rates are derived by integrating the relevant cross sections over the electron distribution. For instance, the collisional ionisation rate
is given by
where n_{e} and n_{i} are the densities of electrons and the ion i considered, is the energydependent collisional ionisation crosssection, is the kinetic energy of an electron colliding with the ion and is the electron velocity distribution.
In the most simple approximation, the cross section for collisional ionisation can be written as (Lotz 1967):
where E is the kinetic energy of the free electron, n_{s} is the number of electrons in a given atomic shell and the normalisation m^{2} keV^{2}. The cross section (2) can be inserted into (1), and the integration can be done analytically for simple electron distribution functions like Maxwellians or power laws, resulting in an expression proportional to E_{1}(I/kT)/IT^{0.5}, where I is the ionisation potential of the relevant shell, and is the exponential integral. Equation (2) is obviously too simple in realistic cases. But more sophisticated approximations to the collisional ionisation cross section can be made, e.g. Younger (1981), that still allow analytical integration for Maxwellian or power law electron distributions:
with and A, B, C and D adjustable parameters. Analytical integration has the obvious advantage of enhanced computational speed, as during the evaluation of a spectrum many atomic shells of hundreds of ions need to be taken into account.
A similar treatment can be made for other relevant rates, like the collisional
excitation rates that are important for the line emission. However, some rates
cannot be written in a form that allows analytical integration over a Maxwellian electron distribution. A good example are the radiative recombination rates. We recall that the cross section for radiative recombination
can be expressed in terms of the photoionisation cross section through Milne's relation as
where is the photoionisation cross section at the photon energy . Although there exist analytical approximations to the photoionisation cross sections (for instance Verner & Yakovlev 1995), they are too complicated to allow for analytical integration of (4) over a Maxwellian electron distribution, in particular as they are more complicated than a simple power law near the ionisation edges (Cooper minima). Even for the hydrogenic case, where the cross section can be calculated analytically, its shape is too complex for analytical convolution with the electron distribution. Therefore, the integration must be done numerically. In practice, often the exact cross sections (including also resonances near the edge) are convolved numerically with a Maxwell distribution, and the resulting rates are then approximated by analytical functions.
It is clear that we need to follow an alternative approach for nonMaxwellian
electron distributions. There have been attempts to make a generalisation of the
Maxwell distribution. For instance, Porquet et al. (2001) combine a Maxwellian at
low energies with a powerlaw at high energies, but they only consider the
ionisation balance, not the resulting spectrum. Prokhorov et al. (2009) adopt a
similar approximation for the electron distribution, but only consider continuum
emission and the Fe XXV and Fe XXVI 1s2p blends.
Owocki & Scudder (1983) use the socalled kappadistribution (Olbert et al. 1967):
This distribution tends to a Maxwellian for and it has a power law tail proportional to for high energies. However, the above approach has also clear disadvantages. Not every nonthermal electron distribution can be cast into the shape of (5), and not all rates can be convolved analytically with (5), in particular again the recombination rates. Owocki & Scudder (1983) had to assume very simple analytical approximations to the cross sections, like (2), in order to be able to evaluate analytically the plasma rates, which are given in terms of hypergeometric functions. Therefore we propose a different approach.
3 The multiMaxwellian approach
Plasma codes generally may contain (ten)thousands of spectral lines, and at high spectral resolution, spectra contain many bins. Therefore, numerical integration for each rate over the electron distribution is a very time consuming task, and not of practical use, in particular when spectral fitting is performed and the model needs to be evaluated many times.
But all relevant rates in the case of the collisionally ionised plasmas that we
consider here, are proportional to the product
n_{e}n_{i}, as
for each of these rates the interaction (collision) of an electron with an ion
is the driving process. Therefore, all rates are linearly proportional to the
electron density
n_{e}. If we decompose the electron distribution into
a linear combination of elementary components, the total rate for any process is
simply the sum of the rates for these individual elementary components. Now, for
Maxwellian electron distributions all rates are well known and fast to evaluate
analytically or with analytical approximations. Thus, if we decompose an
arbitrary electron distribution f(E) into a linear combination of Maxwellians:
with A_{i} the normalisation constant for component i and g_{M}(E,T_{i}) the Maxwell distribution for temperature T_{i}, then the calculation of any rate (recombination, ionisation, excitation) is simple and straightforward. The problem is then reduced to obtaining the normalisations A_{i} and temperatures T_{i} for an arbitrary electron distribution f(E).
In principle, one might solve formally the integral equation
Using
(7) can be cast into the shape of a Laplace transform as
with
(10) 
The calculation of the inverse Laplace transform is not always trivial. For instance, when we consider a monoenergetic beam, f(E) is essentially a deltafunction, and the formal solution is a sinewave with infinite frequency, which is not very practical. Fortunately, in most practical cases f(E) is a smooth function of energy, spanning a broad range of energies.
There are various ways to solve (9). There is software available that can do the inverse Laplace transform numerically (for example d'Amore et al. 1999), but in practice this is difficult, because these algorithms only work if the lefthandside of (9) is an analytical function. However, in practice the electron distributions are given in tabular form, and the algorithms used to interpolate such tables do not produce a formal analytical function (the requirement that the function is infinitely differentiable is usually violated). There are two solutions to this problem that we elaborate in the next subsections: direct approximation of the electron distribution by an analytical function, or direct decomposition of the electron distribution as a sum of Maxwellians.
3.1 Approximation of the electron distribution by an analytical function
In practice, apart from the most pathological cases, electron distributions contain a core that is close to Maxwellian, plus a high energy tail. For an example see Sect. 4. After some experimentation, we found
that the electron distribution of that example can be approximated by the following series:
with x the dimensionless momentum defined by
(12) 
and where T_{0} is a characteristic temperature corresponding to the electron distribution, and the parameters c_{k} and a_{k} can be adjusted to give the best fit.
Equation (11) is an analytical function, and the inverse Laplace
transform in this case can be done analytically, yielding
where is the function and Dirac's deltafunction.
For our example, the normalisations c_{0}c_{4} are 0.755, 0.0609, 2.54, 13.3, and 17.58, respectively; the scale factors a_{0}a_{4} are 0.483, 152, 6843, 57.4 and 12.4. The fit is not perfect, but better than 0.5% for x<30, and better than 1% for x<700. The relative difference to the exact distribution are shown in Fig. 2.
3.2 Direct decomposition of the electron distributions as the sum of Maxwellians
We seek an approximation to f(E) that can be written as a linear combination of n Maxwellians, with n a given number. The solution can be represented as a set of pairs (T_{i},A_{i}) with T_{i} the temperatures and A_{i} the normalisations. We will choose for convenience a set of increasing temperatures ( T_{i+1}>T_{i}), and we will demand that (physically allowed components). Most realistic electron distributions can be represented as a Maxwellian with highenergy tails. Therefore, the first temperature T_{1} should be close to the temperature T_{0} of this dominant Maxwellian component.
This then leads to the following approach. We define a set of nlogarithmically spaced temperature intervals (T_{1,i},T_{2,i}) with T_{1,1} equal to T_{0} and for all i, we take T_{2,i}=sT_{1,i} and T_{1,i}=T_{2,i1}. For the parameter s>1 we adopt a value of 1.5. The total temperature range spanned by these n contiguous intervals is therefore T_{0}  s^{n} T_{0}. For each interval, we choose an arbitrary temperature T_{i} within this interval ( ) and an arbitrary normalisation . We scale the A_{i} in such a way that their sum corresponds to the total electron density. We then check how close the resulting electron distribution is to the electron distribution f(E) that we want to approximate. By varying all values T_{i} within their allowed ranges (the above mentioned intervals) and also the values for A_{i}, we try to get the best possible approximation.
We use the genetic algorithm PIKAIA developed by Charbonneau (1995) to find the best solution. This algorithm needs a formal function of the parameters (T_{i}, A_{i}) that needs to be maximised. For this function we choose the inverse of the sum of the squared relative deviations of the approximation to the true electron distribution f(E).
After some experimenting, we obtained satisfactory convergence with the following initial parameters: number of individuals 500, number of generations 500. All other parameters in PIKAIA were kept to their default values. The number of individuals represents here the number of sets (T_{i},A_{i}) from which the iteration starts; the number of generations is the number of iterations during which the individual sets evolve towards the optimum solution.
We have tested our algorithm on various different electron distributions, but give here only one illustrative example.
4 Example of decomposition of an electron distribution into Maxwellians
Figure 1: Electron distributions immediately after a shock with different Mach numbers, in the downstream region. The distributions are expressed here in terms of the dimensionless momentum (normalised to the thermal momentum in the far upstream region). The distributions are normalised to integral unity. Solid line: M = 2.2; dashed line: M=1.5; dashdotted line: M=1.2; dotted line: far upstream Maxwellian distribution. 

Open with DEXTER 
Figure 2: Relative differences of approximations to the electron distribution with M = 2.2 of Fig. 1. Upper panel: analytical approximation (11), see Sect. 3.1; lower panel: approximation using the sum of 32 Maxwellian components, see Sect. 3.2. 

Open with DEXTER 
Figure 3: Relative normalisations versus temperature for the 32 Maxwellian components used in the approximation of the electron distribution that is shown in Fig. 2. Histogram: the analytical approximation (13); Dots: the solution from the genetic algorithm. 

Open with DEXTER 
We have used electron distributions based on the models of Bykov & Uvarov (1999) for collisionless MHD shocks taking into account particle acceleration and a shock precursor. We consider here a case for a preshock temperature of T_{0}=10^{8} K, electron density 10^{3} m^{3}, magnetic field 10^{10} T, and shock size 100 kpc. In the simulation we assumed a pure Maxwellian distribution of electrons in the far upstream region of the flow and thus only the direct injection of electrons from the upstream thermal pool to Fermi type shock acceleration is responsible for the nonthermal electron population. The electron distribution in the immediate postshock region is shown in Fig. 1 for three different Mach numbers. The model distribution was calculated taking account of both Fermitype acceleration in a collisionless shock and Coulomb losses of the electrons in both the upstream and downstream regions. It is normalised to the momentum p_{0} corresponding to the typical energy kT_{0}. Note that the suprathermal electron distribution well away from the shock front differs from that shown in Fig. 1, because the efficient Coulomb losses washout the nonrelativistic suprathermal electrons.
There is a variant of the model where mildly relativistic electrons (with Lorentz factors 30) comprising a putative longlived cosmic ray electron population in the ICM are reaccelerated by the MHDshocks. The energy efficiency problem is alleviated in that model in comparison with the case of only direct particle injection from the thermal plasma. However, in this paper we concentrate mostly on the diagnostic of nonrelativistic electrons and thus we consider only a direct injection model.
For our example, we take the case of M=2.2. In our decomposition with the genetic algorithm we have taken n=32 Maxwellians. The relative differences of the approximation compared to the exact electron distribution is shown in Fig. 2. The temperatures and relative normalisations of the solution are shown in Fig. 3. We have also used the analytical approximation (13) and binned it into similar temperature bins as the solutions from the genetic algorithm.
5 Spectrum for the suprathermal electron distribution
We have adjusted all models in the spectral fitting package SPEX that involve emission or absorption from a hot plasma. All these models now include an option to account for suprathermal electrons. They have an additional parameter, which is the name of a file containing the temperatures and relative emission measures of the Maxwellian components. As SPEX does not use precalculated tables but calculates spectra on the fly, all relevant rates (ionisation, recombination, excitation) are simply calculated by adding the contributions from the Maxwellian components. Obviously, this process is done in two steps: first the composite multiMaxwellian electron distribution is applied to determine the ionisation balance, and with the resulting nonequilibrium ion concentrations, we calculate the Xray spectrum for the nonequilibrium electron distribution. It is tacitly assumed here that we consider a timeindependent, steady situation. For the example given in the previous section, only a few dozen Maxwellian components are needed. This allows fast and accurate evaluation of the spectrum, without the need to make simplifications to the atomic physics.
The most important effects of nonthermal electrons on the spectrum is then a shift of the ionisation balance towards higher ionisation, the production of a nonthermal Bremsstrahlung tail on the continuum spectrum, and enhanced satellite line emission. These satellite lines (for instance in the FeK band) can be detected easily with highresolution spectrometers that will fly on future missions such as AstroH and IXO. Their relevance as indicators for nonthermal electrons was already indicated by Gabriel & Phillips (1979).
Figure 4: Crosses: simulated 100 ks calorimeter spectrum for IXO ( top panel) and AstroH ( bottom panel) as described in the text, for the suprathermal electron distribution of Fig. 1. Solid line: bestfit model to a pure Maxwellian plasma, with temperature 16.99 keV. Note the excess emission of satellite lines in the data, in particular the Fe XXV jline. 

Open with DEXTER 
To illustrate the effects of such a suprathermal electron distribution on data, we simulated an XMMNewton EPIC/pn spectrum extracted from a circular region with a radius of 1 centred on the core of a bright cluster with a 0.310 keV luminosity of W within the extraction region, at an assumed redshift of z=0.055. In the simulation of the spectrum, we assumed the above mentioned postshock downstream electron distribution for the Mach number of M=2.2 and preshock temperature of kT=8.62 keV (10^{8} K). We assume a deep 100 ks observation. The resulting spectrum has very high statistics, which should in principle allow to detect any nonisothermality of the plasma. The best fit temperature of the simulated spectrum is keV ( K), consistent with the postshock temperature of the plasma given by the RankineHugoniot jump condition. An isothermal model fits the data extremely well (reduced ) and the nonthermal tail of the electron distribution cannot be detected in the spectrum.
We also simulated a spectrum with the same input parameters as observed during a deep 100 ks observation with the Xray microcalorimeter on the proposed International Xray Observatory (IXO). We fitted the simulated spectrum with an isothermal model and obtained a temperature of keV, about 1 keV lower than the expected postshock temperature. In Fig. 4 we show the 6.97.0 keV part of the spectrum (restframe energies) which shows the Fe XXVI Ly lines and the Fe XXV jsatellite line. Enhanced equivalent widths of satellite lines are good indicators of nonthermal electrons. The satellite line in the simulated spectrum in Fig. 4 is clearly stronger than that predicted by the thermal model with a Maxwellian electron distribution. This exercise illustrates, that in order to observationally reveal nonMaxwellian tails in the electron distributions, we will need highresolution spectra obtained by future satellites with a large effective area.
However, even before IXO, AstroH (expected launch 2014) will be able to detect suprathermal electrons. We simulated the same spectrum for the same extraction region as above for the main instruments of AstroH, again for 100 ks exposure time. The two hard Xray telescopes detect the source up to 75 keV, and the spectrum is well approximated ( for 223 degrees of freedom) by an isothermal model with measured temperature of keV. Thus, the presence of such a small amount of suprathermal electrons cannot be revealed as a hard tail, but it can be revealed in highresolution spectra. Figure 4 shows the simulated spectrum for the Soft Xray Spectrometer (SXS) of AstroH. The excess flux at the Fe XXV jsatellite has a significance. Obviously, longer exposure times will enhance the significance.
Interestingly, in all our simulations above, the bestfit iron abundance for the isothermal model is about 30% higher than the actual abundance that we have put into our model spectrum with suprathermal electrons. This holds also if we restrict our fit only to the Fe Lshell or Fe Kshell band. This is because the higherenergy electrons have less efficient line emissivity relative to the continuum, compared to lowerenergy electrons. Without knowing the amount of suprathermal electrons, which can be determined only from highresolution spectra, this abundance bias cannot be resolved and will result into incorrect interpretations.
6 Concluding remarks
In practice most effort goes into finding a good decomposition of an electron distribution into Maxwellians. We have indicated and illustrated in this paper two different methods: fits to simple analytical models that allow analytical inversion of the Laplace transform, and a direct decomposition with a genetic algorithm.
Finally we note that all relevant plasma rates that are used in the SPEX code are calculated with nonrelativistic approximations. For instance, ionisation cross sections are approximated by analytical functions like (3) that loose their validity for relativistic energies. Thus, for electron distributions containing a significant fraction of relativistic electrons, the results will be less accurate. Fortunately, in most situations this is not a problem. For instance, the electron distribution of Fig. 1 has a highenergy tail roughly proportional to p^{3}. Inserting this for instance into (1) and using (2) for the highenergy limit, shows that the integrand scales, apart from a logarithmic term, proportional to E^{3}, and therefore the highest energy electrons do not contribute much to the rates. Therefore, even though the approximations made to the decomposition of the electron distribution, in particular the genetic algorithm, are not always very accurate at high energies (see Fig. 2), this affects the final spectrum to a much lesser extent.
Only in astrophysical situations with a significant number of relativistic electrons our method will not apply. Bykov (2002) has given an example of this in the context of supernova remnants, considering only line fluorescence due to collisional ionisation. Good approximations to the relativistic collisional ionisation cross section of Kshell and Lshell electrons are available (see references in Bykov 2002), but for a full plasma model relativistic corrections to all rates would be needed, which is beyond the scope of the present paper.
Acknowledgements
The Netherlands Institute for Space Research is supported financially by NWO, the Netherlands Organisation for Scientific Research. A.M. Bykov was supported by RBRF grant 090212080. Support for this work was provided by the National Aeronautics and Space Administration through Chandra Postdoctoral Fellowship Award Number PF890056 issued by the Chandra Xray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics and Space Administration under contract NAS803060.
References
 Bykov, A. M. 2002, A&A, 390, 327 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Bykov, A. M., & Uvarov, Y. A. 1999, Sov. J. Exper. Theor. Phys., 88, 465 [NASA ADS] [CrossRef] (In the text)
 Bykov, A. M., Paerels, F. B. S., & Petrosian, V. 2008, Space Sci. Rev., 134, 141 [NASA ADS] [CrossRef] (In the text)
 Charbonneau, P. 1995, ApJS, 101, 309 [NASA ADS] [CrossRef] (In the text)
 d'Amore, L., Laccetti, G., & Murli, A. 1999, ACM Trans. Math. Softw., 25, 306 [CrossRef] (In the text)
 Gabriel, A. H., & Phillips, K. J. H. 1979, MNRAS, 189, 319 [NASA ADS] (In the text)
 Kaastra, J. S., Paerels, F. B. S., Durret, F., Schindler, S., & Richter, P. 2008, Space Sci. Rev., 134, 155 [NASA ADS] [CrossRef] (In the text)
 Lotz, W. 1967, Z. Phys., 206, 205 [NASA ADS] [CrossRef] (In the text)
 Masai, K., Dogiel, V. A., Inoue, H., Schönfelder, V., & Strong, A. W. 2002, ApJ, 581, 1071 [NASA ADS] [CrossRef] (In the text)
 Olbert, S., Egidi, A., Moreno, G., & Pai, L. G. 1967, Trans. AGU, 48, 177 (In the text)
 Owocki, S. P., & Scudder, J. D. 1983, ApJ, 270, 758 [NASA ADS] [CrossRef] (In the text)
 Porquet, D., Arnaud, M., & Decourchelle, A. 2001, A&A, 373, 1110 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Prokhorov, D. A., Durret, F., Dogiel, V., & Colafrancesco, S. 2009, A&A, 496, 25 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Tatischeff, V. 2003, in EAS Publications Series, ed. C. Motch & J.M. Hameury, EAS SP, 7, 79 (In the text)
 Verner, D. A., & Yakovlev, D. G. 1995, A&A, 109, 125 [NASA ADS] (In the text)
 Younger, S. M. 1981, J. Quant. Spectr. Rad. Trans., 26, 329 [NASA ADS] [CrossRef] (In the text)
All Figures
Figure 1: Electron distributions immediately after a shock with different Mach numbers, in the downstream region. The distributions are expressed here in terms of the dimensionless momentum (normalised to the thermal momentum in the far upstream region). The distributions are normalised to integral unity. Solid line: M = 2.2; dashed line: M=1.5; dashdotted line: M=1.2; dotted line: far upstream Maxwellian distribution. 

Open with DEXTER  
In the text 
Figure 2: Relative differences of approximations to the electron distribution with M = 2.2 of Fig. 1. Upper panel: analytical approximation (11), see Sect. 3.1; lower panel: approximation using the sum of 32 Maxwellian components, see Sect. 3.2. 

Open with DEXTER  
In the text 
Figure 3: Relative normalisations versus temperature for the 32 Maxwellian components used in the approximation of the electron distribution that is shown in Fig. 2. Histogram: the analytical approximation (13); Dots: the solution from the genetic algorithm. 

Open with DEXTER  
In the text 
Figure 4: Crosses: simulated 100 ks calorimeter spectrum for IXO ( top panel) and AstroH ( bottom panel) as described in the text, for the suprathermal electron distribution of Fig. 1. Solid line: bestfit model to a pure Maxwellian plasma, with temperature 16.99 keV. Note the excess emission of satellite lines in the data, in particular the Fe XXV jline. 

Open with DEXTER  
In the text 
Copyright ESO 2009