Numerical scheme for treatment of Uehling-Uhlenbeck equation for two-particle interactions in relativistic plasma
We present a new efficient method to compute Uehling-Uhlenbeck collision integral for all two-particle interactions in relativistic plasma with drastic improvement in computation time with respect to existing methods. Plasma is assumed isotropic in momentum space. The set of reactions consists of: Moeller and Bhabha scattering, Compton scattering, two-photon pair annihilation, and two-photon pair production, which are described by QED matrix elements. In our method exact energy and particle number conservation laws are satisfied. Reaction rates are compared, where possible, with the corresponding analytical expressions and convergence of numerical rates is demonstrated.
keywords:Uehling-Uhlenbeck equations, collision integral, binary interactions, relativistic plasma.
Relativistic plasma, for which , where is Boltzmann constant, is speed of light, is electron mass, is temperature, is relevant in different branches of astrophysics. In the early universe ultrarelativistic electron-positron pairs contribute to the matter contents of the Universe weinberg2008cosmology (). X-ray and gamma-ray radiation from numerous astrophysical sources such as gamma-ray bursts 1999PhR…314..575P (); 2010PhR…487….1R (); 2015PhR…561….1K (), active galactic nuclei 2012AAT…27..557A (); blandford2013active (), and X-ray binaries 2006ARA&A..44..323F () points out to existence of relativistic electron-positron plasma in these objects. The upcoming high-energy laser facilities aiming at generation of femtosecond laser pulses with intensity more than aim generation of relativistic plasma by interacting laser pulses. At present relativistic electron-positron jets are generated by interaction of laser pulses with condensed matter 2015NatCo…6E6747S (); 0741-3335-53-1-015009 (); 1742-6596-454-1-012016 (); PhysRevLett.102.105001 ().
Solving the Boltzmann equations with collision integral containing a quantum cross-section represents the most general and complete method to describe a behavior of relativistic plasma vereshchagin2017relativistic (); cercignani2012relativistic (); groot1980relativistic (). The one-particle distribution function (DF) is defined on a seven dimensional space, three dimensions for the physical space and three dimensions for the momentum space, and one dimension for the time. Thus one has a multidimensional problem which is a real challenge from the computational point of view. Beside the dimensionality problem, there are other difficulties which are related to kinetic equations in general 1990JMP….31..245B (); 1991RvMaP…3..137B (). Our main goal in this paper is to tackle the challenge associated with the calculation of the collision integral, dealing with two key issues. First, the computational cost related to the evaluation of the collision operator involving multidimensional integrals which should be solved in each point of the coordinate space. Second, the presence of multiple scales requires the development of adapted numerical schemes capable of solving stiff dynamics. Different deterministic approaches are used to tackle collision integral from a numerical point of view: finite volume, semi-Lagrangian and spectral schemes 1991RvMaP…3..137B (); dimarco_pareschi_2014 (); 2006MaCom..75.1833M (); DIMARCO2017 (); WU201327 (). While the deterministic methods could normally reach high order of accuracy, the probabilistic ones, such as Monte-Carlo (MC) method, are often faster.
MC methods are traditionally used to model Coulomb interactions in non-relativistic plasmaSHERLOCK20082286 (); HUTHMACHER2016535 (); TURRELL2015144 (); BOBYLEV2013123 (). As a rule MC techniques are based on the random pairing of particles in close vicinity and the calculation of a scattering angle due to the interaction. Small-angle Coulomb collisions which allow small energy and momentum transfer are often described in diffusion approximation by the Fokker-Planck equation 1981phki.book…..L (). The principal feature of relativistic plasma is a presence of pair creation and pair annihilation processes, which are often included in MC based models 0741-3335-53-1-015009 (); 1742-6596-454-1-012016 (). However Fokker-Planck approximation is no longer valid in relativistic plasma 2009PhRvD..79d3008A ().
Classical Boltzmann equation does not take into account quantum statistics of particles. The generalization of classical Boltzmann equation including quantum corrections is Uehling-Uhlenbeck (U-U) equation, which contains additional Pauli blocking and Bose enhancement multipliers that give rise to equilibrium solution with Bose-Einstein and Fermi-Dirac distributions 1934PhRv…46..917U (); 1933PhRv…43..552U (). The main problem of the MC methods in application to U-U equations is that total reaction rate is unknown as distribution function is unknown too. Compensation methods include smoothing of the delta-function distribution of MC-particles over cells in the phase space, but it suffers from a large number of simulation particles and cells needed to reproduce Bose-Einstein steady state distribution. Spectral methods based on the Fourier transformation of the velocity distribution function require very dense computational grid to reach high accuracy 2017JCoPh.330.1010Y (); Hu2015 (); huying12 (); PhysRevE.68.056703 (); 2010arXiv1009.3352F (). Process-oriented approach to the U-U collision integral presented in this work allows one to get high accuracy results with low computational cost.
In this paper we further develop the method first used in the work 2004ApJ…609..363A (). This method was successfully applied to follow the thermalization of relativistic plasma 2007PhRvL..99l5003A (); 2009AIPC.1111..344A (); 2009PhRvD..79d3008A (); 2010AIPC.1205…11A () and to investigate thermalization timescales for an electron-positron plasma 2010PhRvE..81d6401A (). In section 2 we recall Boltzmann and UU equations and present usual scheme of their analytic treatment. Section 3 is devoted to the description of our numerical scheme, while section 4 shows comparison between our code results and known analytic formulae for non-degenerate case. Conclusion follows.
The Boltzmann equation governs an evolution of one-particle distribution function . We assume that plasma is homogeneous and isotropic in coordinate space and isotropic in momentum space, thus distribution function depends on absolute value of momentum (energy) and time. DF is normalized on particles concentration, so that .
Consider an interaction of two initial particles of type I and II which are in states 1 and 2, correspondingly, and creation of two final particles of type III and IV which are in states 3 and 4, correspondingly. Let us image the process by the following scheme:
The corresponding inverse process is:
If every particle has momentum , which lies in interval , then a number of interactions in unit time value and unit space volume is:
function is called a transition rate for a given reaction.
An effective cross-section of a given interaction is defined by the formula:
where is a relative velocity of initial particles.
In quantum field theory an expression for interaction cross-section is:
where are a matrix elements calculated with a methods of quantum field theory.
Comparing two last formulas one can derive the following expression for transition rate:
Let us write the Boltzmann equation for DF of particle I for a given process:
equations for particle DFs of remaining types can be derived by the corresponding replacement of indices.
Specifically, for a scattering with and the inverse process is the same as the direct one since pairs of indices and can be interchanged. The relation holds for all processes listed in Table1. The right hand side of the Boltzmann equation is a collision integral denoted as St. The first term in collision integral describes outgoing particles which leave a given element of the phase space (denoted ) and the second term describes incoming particles which arrive to a given element of the phase space (denoted ).
For the particle in the state U-U equation has the following form:
where is defined through
For numerical evaluation phase space is divided into zones, in calculations we approximate continuous DF by its average over each zone (see Eq. (22) below). For this purpose we add an integral over in U-U equation (8). The right-hand side of resulting equation will have the same form for each particle type differing only by sign and its integration limits:
where the upper sign corresponds to particles type I and II and the lower sign corresponds to particles type III and IV.
Particle in a given momentum state may participate in different reactions, which enter collision integrals multiple times. The essence of the reaction-oriented approach consists in evaluation of expression (10) to be used for calculation of collision integrals for all reactions simultaneously. In this way a given reaction with fixed particle state is considered only once, allowing considerable reduction in the computation time.
In this paper we deal with all two-particle QED processes in relativistic plasma, which are collected in Table 2. The exact QED matrix elements for given processes can be found in the standard textbooks, e.g. 2003spr..book…..G (); 1982els..book…..B ().
|Compton Scattering (CS)|
|Bhabha Scattering (BS)|
|Moeller Scattering (MS)|
|Pair Annihilation (PA)|
|Pair Creation (PC)|
Let us make a remark regarding conservation laws for interacting particles. Energy and momentum conservations read 2009PhRvD..79d3008A ()
In the integration over energy of particle it is necessary to take into account that is now a function of energy and angles of particles and , as well as angles of particle , so we have
where is the unit vector in the direction of particle momentum, is the absolute value of particle momentum, , and a dot denotes scalar product of 3-vectors.
We use spherical coordinates in momentum space: , , where is the particle energy, and and are polar and azimuthal angles, respectively. Then energy and angles of particle and energy of particle follow from energy and momentum conservations (11) and relativistic energy-momentum relation, namely
Then we introduce these relations into collision integral (10). We also use spherical symmetry in momentum space to fix angles of the particle : , and to perform the integration over azimuthal angle of particle : , setting in the remaining integrals. Then final expression for collision integral is
For numerical integration, however, another expression is proved useful
where the first term is expressed in the form ready for replacement by the sum over incoming particles and . In this term are given by relations (18) with indices exchange , , , .
This collision integral of any of two-particle processes is a four-dimensional integral in momentum space. In Sec. 3 we show how such integral is computed numerically on a finite grid.
Here we note that in the case of homogeneous and isotropic pair plasma one has to satisfy only two conservation laws, namely energy and particle conservation laws. Momentum conservation law should be added for anisotropic in momentum space DF, see e.g. BENEDETTI2013206 (). In our method electric charge is conserved due to conservation of particles because we use between cell interpolation for the same kind of particles.
3 Numerical Scheme
The phase space is divided in zones. The zone for particle specie corresponds to energy , cosine of polar angle and azimuthal angle , where indices run in the following ranges , , and . The zone boundaries are , , . The length of the -th energy zone is . On finite grid does not depend on and , and number density of particle in zone is
In this variables discretized U-U equation for particle and energy zone reads
where the sum is taken over all processes involving particle . Coefficients of particles income and outcome on the grid are obtained by integration of (21) for two-particle processes over the zone. The corresponding integrals are replaced by sums on the grid. For instance, coefficient of particles income in two-particle process (1) is
In integration of (21) over the zone one can integrate out the -function . However, when energies of incoming particles are fixed on the grid, the energies of outgoing particles are not on the grid. Hence an interpolation (26) is adopted, which enforces the exact number of particles and energy conservation in each two-particle process due to redistribution of outgoing particle with energy over two energy zones with .
The redistribution of final particles should also satisfy requirements of quantum statistics. Therefore if a process occurs, when final particle should be distributed over the quantum states which are fully occupied, such process should be forbidden. Thus we introduce the Bose enhancement/Pauli blocking coefficients in (24) and (25) as
The sum over angles can be found once and for all at the beginning of the calculations. We then store in the program for each set of the incoming and outgoing particles the corresponding terms and redistribution coefficients given by Eq. (26).
where constant coefficients are obtained from the summation over angles in the sums (24), (25). In the nondegenerate case of Boltzmann equation the indices in the first sum and in the second sum become dummy, equation (28) can be partially summed and takes the following form:
where . The last quantity is essentially reaction rate usually used for description of binary processes simply connected to the total cross section. The full U-U equation (23) contains similar sums for all processes from Table 2. Each individual term in these sums appears in the system of discretized equations four times in emission and absorption coefficients for each particle entering a given process. Then each term can be computed only once and added to all corresponding sums, that is the essence of our "reaction-oriented" approach ivanphd (); 2015AIPC.1693g0007S ().
We point out that unlike classical Boltzmann equation for binary interactions such as scattering, more general interactions are typically described by four collision integrals for each particle that appears both among incoming and outgoing particles.
4 Numerical results
The results of numerical calculations are presented below. As all known analytical expressions for reaction rates in relativistic plasma concern nondegenerate case, we present our results for nondegenerate plasma. We do not present numerical results for Coulomb scattering due to absence of the corresponding analytical results. Notice that for Coulomb scattering we have implemented a cutoff scheme based on minimal scattering angle 2009PhRvD..79d3008A (); 1988A&A…191..181H (). We consider mildly relativistic plasma with
where is particle kinetic energy divided by electron rest energy, this range contains both relativistic and non-relativistic domains. The upper limit is chosen to avoid thermal production of other particles such as neutrinos and muons, while the lower limit is required to have sufficient pair density.
We introduce logarithmic energy grid with nodes for all calculations and different homogeneous grids for angular variables, -grid is 2 time denser then -grid (typically -grid contains nodes). We use Coppi&Blandford 1990MNRAS.245..453C () analytical expressions (2.3),(3.2),(4.3) for , which corresponds to . Belmont 2009A&A…506..589B () formula (10), Svensson 1982ApJ…258..321S () formula (55), Peer&Waxman 2005ApJ…628..857P () formula (28) are used for quantity , which corresponds to .
To compare numerical results with analytical one, we introduce the following quantity for each given process
expressing relative deviations of numerical results from analytical one for all energy grid nodes, where is analytical expression for reaction rate. Table 2 presents values of for selected number of angular grid nodes. It is evident that the relative error decreases with increasing of number of angular grid nodes reaching about 1 % with 128 nodes. This demonstrates convergence of numerical results to the corresponding analytical ones.
Below we present several plots for the reaction rates of Compton scattering and creation/annihilation of pairs. Energy is measured in electron rest energy units. Presented results reproduce both nonrelativistic and relativistic energy cases.
Compton scattering presents well-known challenge for numerical treatment as all the analytical formulas for scattering rate behave badly numerically in different parameter areas, see e.g. 2005ApJ…628..857P (); 2009A&A…506..589B (). We easily bypass this difficulty as we numerically integrate well-behaved differential cross-section, as one can see for non-relativistic regime in Fig.1 and for relativistic regime in Fig.2. Figure 1 presents analytic photon spectrum for the reaction as solid line and our numerical results shown by dots. Overall there is good agreement between numerical and analytical results. Small deviations in high-energy of the spectrum arise from leakage of the particles to kinematically forbidden area at the boundary of energy zones.
Figure 2 shows the total reaction rate of the same process. Again there is good agreement between numerical and analytical results. Small discrepancy arises from truncation of reactions where final particles get out of the grid to higher or lower energies. As a result numerical reaction rates are systematically lower than analytic ones.
Annihilation photon spectrum for reaction is illustrated in Fig.3 and total reaction rate in this process in Fig.4. Figure 3 shows that the method is able to accurately reproduce the spectrum of annihilation photons in the range of more than two orders of magnitude. Reaction truncation errors, hardly seen at Fig.4, are much lower for annihilation as low-energy photons are rare in this process.
Balance between pair creation and annihilation represent an independent test for the numerical scheme, as it is not automatically satisfied due to different numerical treatment of incoming and outgoing particles in the reactions. Pair creation spectra for reaction are reproduced well, see Fig. 5, as well as total reaction rates, see Fig. 6. Numerical balance can be checked by the form of particle distributions in numerical equilibrium, that was verified to be within 5 % of corresponding Boltzmann distributions.
In this paper, we propose a new numerical method to accurately calculate Uehling–Uhlenbeck collision integral for two-particle interactions in relativistic plasma. After calculation of collision integral discretized Uehling–Uhlenbeck equations transforms into system of ODEs, which can be treated by various methods suitable to solve stiff ODEs. The method admits parallelization on GPU/CPU. Improvement in computation time with respect to previous work is achieved. Our reaction-oriented approach can be easily applied to any other types of particles and any other binary interactions, for instanse, weak interactions of neutrinos or electromagnetic ones of protons. Generalization of the proposed method for triple interactions is straightforward.
Our results show that reaction rates in relativistic plasma are well reproduced with moderate number of grid nodes in energy and angles (see Figures and Table 2) both for non-relativistic and relativistic particle energies. This allows development of an efficient method of solution for relativistic Uehling–Uhlenbeck equation.
- S. Weinberg. Cosmology. OUP Oxford, 2008.
- T. Piran. Gamma-ray bursts and the fireball model. Physics Reports, 314:575–667, June 1999.
- R. Ruffini, G. Vereshchagin, and S.-S. Xue. Electron-positron pairs in physics and astrophysics: From heavy nuclei to black holes. Physics Reports, 487:1–140, February 2010.
- P. Kumar and B. Zhang. The physics of gamma-ray bursts & relativistic jets. Physics Reports, 561:1–109, February 2015.
- R. Antonucci. A panchromatic review of thermal and nonthermal active galactic nuclei. Astronomical and Astrophysical Transactions, 27:557–602, 2012.
- P.R.D. Blandford, P.H. Netzer, P.L. Woltjer, T.J.L. Courvoisier, and P.M. Mayor. Active Galactic Nuclei. Saas-Fee Advanced Course. Springer Berlin Heidelberg, 2013.
- G. Fabbiano. Populations of X-Ray Sources in Galaxies. Annual Review of Astronomy and Astrophysics, 44:323–366, September 2006.
- G. Sarri, K. Poder, J. M. Cole, W. Schumaker, A. di Piazza, B. Reville, T. Dzelzainis, D. Doria, L. A. Gizzi, G. Grittani, S. Kar, C. H. Keitel, K. Krushelnick, S. Kuschel, S. P. D. Mangles, Z. Najmudin, N. Shukla, L. O. Silva, D. Symes, A. G. R. Thomas, M. Vargas, J. Vieira, and M. Zepf. Generation of neutral and high-density electron-positron pair plasmas in the laboratory. Nature Communications, 6:6747, April 2015.
- R Duclous, J G Kirk, and A R Bell. Monte carlo calculations of pair production in high-intensity laser–plasma interactions. Plasma Physics and Controlled Fusion, 53(1):015009, 2011.
- Toseo Moritaka, Luca Baiotti, An Lin, Li Weiwu, Youichi Sakawa, Yasuhiro Kuramitsu, Taichi Morita, and Hideaki Takabe. Plasma particle-in-cell simulations with qed reactions for pair production experiments using a high-z solid target. Journal of Physics: Conference Series, 454(1):012016, 2013.
- Hui Chen, Scott C. Wilks, James D. Bonlie, Edison P. Liang, Jason Myatt, Dwight F. Price, David D. Meyerhofer, and Peter Beiersdorfer. Relativistic positron creation using ultraintense short pulse lasers. Phys. Rev. Lett., 102:105001, Mar 2009.
- G.V. Vereshchagin and A.G. Aksenov. Relativistic Kinetic Theory: With Applications in Astrophysics and Cosmology. Cambridge University Press, 2017.
- C. Cercignani and G.M. Kremer. The Relativistic Boltzmann Equation: Theory and Applications. Progress in Mathematical Physics. Birkhäuser Basel, 2012.
- S.R. Groot, W.A. Leeuwen, and C.G. Weert. Relativistic kinetic theory: principles and applications. North-Holland Pub. Co., 1980.
- N. Bellomo and S. Kawashima. The discrete Boltzmann equation with multiple collisions: Global existence and stability for the initial value problem. Journal of Mathematical Physics, 31:245–253, January 1990.
- N. Bellomo and T. Gustafsson. The Discrete Boltzmann Equation:. a Review of the Mathematical Aspects of the Initial and Initial-Boundary Value Problems. Reviews in Mathematical Physics, 3:137–162, 1991.
- G. Dimarco and L. Pareschi. Numerical methods for kinetic equations. Acta Numerica, 23:369–520, 2014.
- C. Mouhot and L. Pareschi. Fast algorithms for computing the Boltzmann collision operator. Mathematics of Computation, 75:1833–1852, December 2006.
- Giacomo Dimarco, Raphaël Loubère, Jacek Narski, and Thomas Rey. An efficient numerical method for solving the boltzmann equation in multidimensions. Journal of Computational Physics, 2017.
- Lei Wu, Craig White, Thomas J. Scanlon, Jason M. Reese, and Yonghao Zhang. Deterministic numerical solutions of the boltzmann equation using the fast spectral method. Journal of Computational Physics, 250(Supplement C):27 – 52, 2013.
- M. Sherlock. A monte-carlo method for coulomb collisions in hybrid plasma models. Journal of Computational Physics, 227(4):2286 – 2292, 2008.
- Klaus Huthmacher, Andreas K. Molberg, Bärbel Rethfeld, and Jeremy R. Gulley. A split-step method to include electron–electron collisions via monte carlo in multiple rate equation simulations. Journal of Computational Physics, 322(Supplement C):535 – 546, 2016.
- A.E. Turrell, M. Sherlock, and S.J. Rose. Self-consistent inclusion of classical large-angle coulomb collisions in plasma monte carlo simulations. Journal of Computational Physics, 299(Supplement C):144 – 155, 2015.
- A.V. Bobylev and I.F. Potapenko. Monte carlo methods and their analysis for coulomb collisions in multicomponent plasmas. Journal of Computational Physics, 246(Supplement C):123 – 144, 2013.
- E. M. Lifshitz and L. P. Pitaevskii. Physical kinetics. 1981.
- A. G. Aksenov, R. Ruffini, and G. V. Vereshchagin. Thermalization of the mildly relativistic plasma. Phys. Rev. D, 79(4):043008, February 2009.
- E. A. Uehling. Transport Phenomena in Einstein-Bose and Fermi-Dirac Gases. II. Physical Review, 46:917–929, November 1934.
- E. A. Uehling and G. E. Uhlenbeck. Transport Phenomena in Einstein-Bose and Fermi-Dirac Gases. I. Physical Review, 43:552–561, April 1933.
- R. Yano. Fast and accurate calculation of dilute quantum gas using Uehling-Uhlenbeck model equation. Journal of Computational Physics, 330:1010–1021, February 2017.
- Jingwei Hu, Qin Li, and Lorenzo Pareschi. Asymptotic-preserving exponential methods for the quantum boltzmann equation with high-order accuracy. Journal of Scientific Computing, 62(2):555–574, Feb 2015.
- Jingwei Hu and Lexing Ying. A fast spectral algorithm for the quantum boltzmann collision operator. Commun. Math. Sci., 10(3):989–999, 2012.
- Alejandro L. Garcia and Wolfgang Wagner. Direct simulation monte carlo method for the uehling-uhlenbeck-boltzmann equation. Phys. Rev. E, 68:056703, Nov 2003.
- F. Filbet, J. Hu, and S. Jin. A Numerical Scheme for the Quantum Boltzmann Equation Efficient in the Fluid Regime. ArXiv e-prints, September 2010.
- A. G. Aksenov, M. Milgrom, and V. V. Usov. Structure of Pair Winds from Compact Objects with Application to Emission from Hot Bare Strange Stars. Astrophysical Journal, 609:363–377, July 2004.
- A. G. Aksenov, R. Ruffini, and G. V. Vereshchagin. Thermalization of Nonequilibrium Electron-Positron-Photon Plasmas. Physical Review Letters, 99(12):125003, September 2007.
- A. G. Aksenov, R. Ruffini, and G. V. Vereshchagin. Thermalization of pair plasma with proton loading. In G. Giobbi, A. Tornambe, G. Raimondo, M. Limongi, L. A. Antonelli, N. Menci, and E. Brocato, editors, American Institute of Physics Conference Series, volume 1111 of American Institute of Physics Conference Series, pages 344–350, May 2009.
- A. G. Aksenov, R. Ruffini, and G. V. Vereshchagin. Kinetics of the mildly relativistic plasma and GRBs. In R. Ruffini and G. Vereshchagin, editors, American Institute of Physics Conference Series, volume 1205 of American Institute of Physics Conference Series, pages 11–16, March 2010.
- A. G. Aksenov, R. Ruffini, and G. V. Vereshchagin. Pair plasma relaxation time scales. Phys. Rev. E, 81(4):046401, April 2010.
- J. Ehlers. Survey of general relativity theory. In Relativity, Astrophysics and Cosmology, pages 1–125, 1973.
- W. Greiner and J. Reinhardt. Quantum Electrodynamics. Berlin, Springer, 2003.
- V. B. Berestetskii, E. M. Lifshitz, and V. B. Pitaevskii. Quantum Electrodynamics. Elsevier, edition, 1982.
- A. Benedetti, R. Ruffini, and G.V. Vereshchagin. Phase space evolution of pairs created in strong electric fields. Physics Letters A, 377(3):206 – 215, 2013.
- Siutsou I. PhD thesis, University of Rome, Sapienza, 2013.
- I. A. Siutsou, A. G. Aksenov, and G. V. Vereshchagin. On thermalization of electron-positron-photon plasma. In American Institute of Physics Conference Series, volume 1693 of American Institute of Physics Conference Series, page 070007, December 2015.
- E. Haug. Energy loss and mean free path of electrons in a hot thermal plasma. Astronomy and Astrophysics, 191:181–185, February 1988.
- P. S. Coppi and R. D. Blandford. Reaction rates and energy distributions for elementary processes in relativistic pair plasmas. MNRAS, 245:453–507, August 1990.
- R. Belmont. Numerical computation of isotropic Compton scattering. Astronomy and Astrophysics, 506:589–599, November 2009.
- R. Svensson. The pair annihilation process in relativistic plasmas. Astrophysical Journal, 258:321–334, July 1982.
- A. Pe’er and E. Waxman. Time-dependent Numerical Model for the Emission of Radiation from Relativistic Plasma. Astrophysical Journal, 628:857–866, August 2005.