1. Introduction
Computer simulations for kinetic plasma processes are of core interest for a variety of scenarios, ranging from astrophysics to laboratory experiments. The particle-in-cell (PIC) methodology (Evans & Harlow Reference Evans and Harlow1957; Dawson Reference Dawson1983; Hockney & Eastwood Reference Hockney and Eastwood1988; Bird Reference Bird1989; Birdsall & Langdon Reference Birdsall and Langdon1991) is one of the most popular and widely used techniques, which pioneered the study of collisionless plasmas. The standard PIC loop can be enriched with various quantum electrodynamics (QED) cross sections to investigate astrophysical environments and model laboratory experiments where quantum processes affect the plasma dynamics. These modules rely on Monte Carlo techniques by taking advantage of the inherent stochasticity of QED processes. The coupling of QED Monte Carlo modules to the PIC loop represents a unique numerical tool that allows such scenarios to be studied from first principles. For example, the inclusion of nonlinear Compton scattering (QED synchrotron) is essential to simulate the interaction of matter with ultra-intense electromagnetic fields (Nerush et al. Reference Nerush, Kostyukov, Fedotov, Narozhny, Elkina and Ruhl2011; Ridgers et al. Reference Ridgers, Brady, Duclous, Kirk, Bennett, Arber, Robinson and Bell2012; Blackburn et al. Reference Blackburn, Ridgers, Kirk and Bell2014; Vranic et al. Reference Vranic, Martins, Vieira, Fonseca and Silva2014; Gonoskov et al. Reference Gonoskov, Bastrakov, Efimenko, Ilderton, Marklund, Meyerov, Muraviev, Sergeev, Surmin and Wallin2015; Grismayer et al. Reference Grismayer, Vranic, Martins, Fonseca and Silva2016; Jirka et al. Reference Jirka, Klimo, Bulanov, Esirkepov, Gelfer, Bulanov, Weber and Korn2016; Lobet et al. Reference Lobet, d'Humiéres, Grech, Ruyer, Davoine and Gremillet2016; Vranic et al. Reference Vranic, Grismayer, Fonseca and Silva2016a,Reference Vranic, Grismayer, Fonseca and Silvab; Grismayer et al. Reference Grismayer, Vranic, Martins, Fonseca and Silva2017). Several other radiative energy loss channels can participate in the production of high-energy photons, such as curvature radiation, inverse Compton emission, and Bremsstrahlung. These photons are produced in astronomical sources such as active galactic nuclei, X-ray binaries, supernova remnants, pulsars and gamma-ray bursts can further interact with matter and, in particular, with the surrounding plasma. The wavelengths of high-energy photons are typically smaller than the average inter-particle distance of any tenuous plasma, implying only binary interaction between single photons and electrons. The leading photon–electron (positron) interaction mechanism is single Compton scattering (Compton Reference Compton1923).
The collision of high-energy photons with the plasma electrons is at the core of some fundamental scenarios: explains the saturation properties of cyclotron radiation masers (Dreicer Reference Dreicer1964), the relaxation to the thermal equilibrium of a photon–electron gas (Kompaneets Reference Kompaneets1957; Peyraud Reference Peyraud1968a,Reference Peyraudb,Reference Peyraudc), or the Comptonisation of the microwave background (Sunyaev & Zel'dovich Reference Sunyaev and Zel'dovich1980). These seminal studies approximate the plasma as a gas of free electrons and, thus, neglect its collective behaviour. Frederiksen, Haugbølle & Nordlund (Reference Frederiksen, Haugbølle and Nordlund2008) and, more recently, the present authors (Del Gaudio et al. Reference Del Gaudio, Fonseca, Silva and Grismayer2020) have shown that bursts of hard X-rays can couple to the collective plasma dynamics via incoherent Compton scattering events and drive plasma wakes. Such phenomena can be studied numerically by coupling a Monte Carlo Compton module to the PIC loop (Haugbølle Reference Haugbølle2005; Haugbølle, Frederiksen & Nordlund Reference Haugbølle, Frederiksen and Nordlund2013), in a binary collision module.
The implementation of binary collisions in PIC codes is discussed extensively in the literature, with the main focus on Coulomb collisions (Takizuka & Abe Reference Takizuka and Abe1977; Wilson, Horwitz & Lin Reference Wilson, Horwitz and Lin1992; Miller & Combi Reference Miller and Combi1994; Vahedi & Surendra Reference Vahedi and Surendra1995; Nanbu Reference Nanbu1997; Larson Reference Larson2003; Kawamura & Birdsall Reference Kawamura and Birdsall2005; Sentoku & Kemp Reference Sentoku and Kemp2008; Sherlock Reference Sherlock2008; Peano et al. Reference Peano, Marti, Silva and Coppa2009; Turrell, Sherlock & Rose Reference Turrell, Sherlock and Rose2015; Higginson Reference Higginson2017). The usual implementation relies on the approximation of small cumulative scattering angles (Takizuka & Abe Reference Takizuka and Abe1977; Miller & Combi Reference Miller and Combi1994; Nanbu Reference Nanbu1997), which allows the simulation time step that is not bound to resolve the collision frequency to be relaxed. Recently, Turrell et al. (Reference Turrell, Sherlock and Rose2015) and Higginson (Reference Higginson2017) included the effect of large angle deflections in Coulomb collision algorithms. The definition of a cut-off angle allows the occurrence of small-angle collisions or large-angle collisions to be identified, based on the impact parameter of the colliding particles. Coulomb collisions differ from Compton collisions in kinematics. In electron–ion collisions, the recoil on the massive ion can usually be neglected. Instead, in photon–electron collisions the electron recoil cannot be neglected unless in the Thomson limit. Notably, Goudsmit & Saunderson (Reference Goudsmit and Saunderson1940) have developed a electron multiple scattering theory in the case of elastic collisions, which can be applied to the Thomson regime. To the best of the authors’ knowledge, a multiple scattering theory for Compton scattering has not yet been developed. Thus, we employ a single scattering procedure where the collision frequency has to be properly resolved by the time discretisation. In general, the smallness of the cross section would result in single scattering events per simulation time step. For particularly high photon densities, the time resolution must be imposed by the collision routine rather than the usual Courant condition for the field solver.
In this article, we describe the implementation of a single Compton scattering collision module for PIC codes. It relies on first principles that the Klein–Nishina (Klein & Nishina Reference Klein and Nishina1923) cross section is employed with no approximations and allows a self-consistent treatment of the high-frequency radiation coupling with the plasma dynamics. In § 2, we review the basic theory for Compton scattering with particular attention paid to the Lorentz invariant quantities that the model must enforce for reproducing the correct scattering rates in the collision at relativistic energies (Peano et al. Reference Peano, Marti, Silva and Coppa2009). Section 3 is devoted to the implementation of our collision procedure. In § 4, we benchmark our code against problems for which exact analytical solution or formulation exist, namely the scattering photon spectrum of a relativistic charge (Blumenthal & Gould Reference Blumenthal and Gould1970), and the Kompaneets equation (Kompaneets Reference Kompaneets1957). Finally, in § 5, we comment on the computational cost that our module brings as compared with a standard PIC loop. Summary and conclusions are presented in § 6.
2. Compton scattering
Single Compton scattering is the inelastic collision between a photon and an electron (Compton Reference Compton1923). It is the generalisation of Thomson scattering (Thomson Reference Thomson1906), for any value of the incident photon of energy 
$\hbar \omega$ in the electron proper frame of reference. By applying energy and momentum conservation in the electron rest frame, later denoted reference frame (O),
\begin{gather} \hbar\omega + mc^2 = \hbar\omega' + \gamma'mc^2, \end{gather}
\begin{gather}\hbar\boldsymbol{k} = \hbar\boldsymbol{k}'+\boldsymbol{p}', \end{gather}
where 
$\gamma '=\sqrt {1+p'^2/m^2c^2}$ and 
$\omega =ck$, the photon frequency shift over one collision is
\begin{equation} \frac{\omega'}{\omega} = \frac{mc^2}{mc^2+\hbar\omega(1-\cos\theta)}, \end{equation}
where 
$\omega$ (
$\omega '$) is the absorbed (emitted) frequency, and 
$\theta$ is the scattering angle. For 
$\hbar \omega \ll mc^2$, the Thomson limit 
$\omega '\simeq \omega$ is recovered. However, when the incident photon energy approaches and exceeds the electron rest mass energy 
$\hbar \omega \gtrsim mc^2$, the energy transfer becomes relevant. For 
$\hbar \omega \gg mc^2$, at 
$\theta \simeq -{\rm \pi}$, the photon transfers up to half its energy 
$\hbar (\omega -\omega ')\simeq \hbar \omega /2$ to the electron. The classical theory of radiation explains Thomson scattering in terms of plane wave absorption and consequent dipole radiation from the oscillating charge (Landau & Lifshitz Reference Landau and Lifshitz1975; Jackson Reference Jackson1999), but does not predict Compton scattering, which is intrinsically a quantum process.
2.1. Klein–Nishina cross section
In the rest frame of an electron, the single Compton scattering probability is determined by the Klein–Nishina differential (in solid angle 
$\varOmega$) cross section (Klein & Nishina Reference Klein and Nishina1923), which, for unpolarised photons, reads
\begin{equation} \frac{\textrm{d}\sigma}{\textrm{d}\varOmega} = \frac{r_e^2}{2}\left(\frac{\omega'}{\omega}\right)^2\left(\frac{\omega'}{\omega}+\frac{\omega}{\omega'}-\sin^2\theta\right), \end{equation}
where 
$r_e=e^2/mc^2$ is the classical electron radius. By combining (2.3) and (2.4), and integrating over the solid angle 
$\textrm {d}\varOmega =\sin \theta \,\textrm {d}\theta \,\textrm {d}\phi$ (
$\phi$ is the symmetry angle around the direction of the incoming photon) the total cross section reads
\begin{equation} \sigma(\epsilon)= \frac{{\rm \pi} r_e^2}{\epsilon}\left[\left(1-\frac{2}{\epsilon}-\frac{2}{\epsilon^2}\right)\log(1+2\epsilon)+\frac{1}{2}+\frac{4}{\epsilon}-\frac{1}{2(1+2\epsilon)^2} \right], \end{equation}
where 
$\epsilon =\hbar \omega /mc^2$. In the limit for low photon energies
\begin{equation} \lim_{ \epsilon \to 0} \sigma(\epsilon)=\sigma_T \end{equation}
the Thomson cross section is recovered. For high photon energies 
$\epsilon \gg 1$, the cross section has the limiting expression
\begin{equation} \lim_{ \epsilon \gg 1} \sigma(\epsilon)=\frac{3}{8}\sigma_T\frac{\log(2\epsilon)}{\epsilon} \end{equation}and decreases with respect to the incident photon energy.
2.2. Relativistic kinematics and Lorentz invariants
We consider a relativistic electron which propagates along the 
$z$ coordinate at velocity 
$\beta c$ and scatters with a photon at an incident angle 
$\phi$ in the laboratory frame, see figure 1. In the electron proper frame of reference 
$_O$, the incident angle is modified by relativistic effects. In the frame 
$_O$ the incident photon is confined within a small cone (Blumenthal & Gould Reference Blumenthal and Gould1970)
\begin{equation} \tan\phi_O=\frac{\sin\phi}{\gamma(\cos\phi-\beta)} \end{equation}
of aperture 
$1/\gamma$. The photon energy in the frame 
$_O$ reads
\begin{equation} \epsilon_O=\gamma\epsilon(1-\beta\cos\phi). \end{equation}
It varies in the range 
$\epsilon _O\in [\epsilon /2\gamma ,2\gamma \epsilon ]$ according to the incident angle 
$\phi$. In the 
$_O$ frame, the photon energy after scattering obeys (2.3) and in the laboratory frame reads
\begin{equation} \epsilon'=\gamma\epsilon_O'[1+\beta\cos({\rm \pi}-\theta_O-\phi_O)]\simeq \gamma\epsilon_O'(1-\cos\theta_O ), \end{equation}
owing to the Lorentz transformation, where 
$\beta \simeq 1$ and 
$\phi _O\sim 1/\gamma$. In the Thomson regime, 
$\omega _O'\simeq \omega _O$ and the maximum energy achieved over one collision is 
$\epsilon '\simeq 4\gamma ^2\epsilon$, for 
$\phi \simeq {\rm \pi}$ and 
$\theta _O\simeq {\rm \pi}$. In the extreme Klein–Nishina limit, the maximum energy achieved over one collision can be obtained by combining (2.3), (2.9), and (2.10), and reads 
$\epsilon '\simeq \gamma$.
Schematic of the Compton scattering relativistic kinematics.
We now consider the scattering between photons, with distribution function 
$f_{\omega }$, and electrons, with distribution function 
$f_e$. Within a portion of space-time 
$\textrm {d}\boldsymbol {x}\,\textrm {d}t$, the number of collisions is a Lorentz invariant quantity (Groot, Leeuwen & van Weert Reference Groot, Leeuwen and van Weert1980) that is given by
\begin{equation} N = \sigma(\boldsymbol{p},\boldsymbol{k}) c f_{\omega}\,\textrm{d}\boldsymbol{k} f_{e}\,\textrm{d}\boldsymbol{p} \,\textrm{d}\boldsymbol{x}\textrm{d}t. \end{equation}
In general, the cross section 
$\sigma (\boldsymbol {p},\boldsymbol {k})$ depends on the electron momentum 
$\boldsymbol {p}$, and on the photon wavevector 
$\boldsymbol {k}$. As the space-time element 
$\textrm {d}\boldsymbol {x}\,\textrm {d}t$, the distribution functions 
$f_{\omega }$ and 
$f_{e}$, and the speed of light 
$c$ are Lorentz invariant, therefore 
$\sigma (\boldsymbol {p},\boldsymbol {k})\,\textrm {d}\boldsymbol {k}\,\textrm {d}\boldsymbol {p}$ is also Lorentz invariant (Landau & Lifshitz Reference Landau and Lifshitz1975). This invariance allows us to obtain the cross section in any inertial frame (
$\gamma =\sqrt {1+\boldsymbol {p}^2/m^2c^4}$, 
$\epsilon =\hbar |\boldsymbol {k}|/mc$). Knowing the cross section in the electron proper frame of reference (
$\gamma _O=1$, 
$\epsilon _O=\gamma \epsilon -\hbar \boldsymbol {p}\cdot \boldsymbol {k}/m^2c^2$)
\begin{equation} \sigma(\boldsymbol{p},\boldsymbol{k})\,\textrm{d}\boldsymbol{k}\,\textrm{d}\boldsymbol{p} = \sigma(\omega_O)\,\textrm{d}\boldsymbol{k_O}\,\textrm{d}\boldsymbol{p_O}, \end{equation}we finally obtain
\begin{equation} \sigma(\boldsymbol{p},\boldsymbol{k}) = \sigma(\epsilon_O)\frac{\epsilon_O}{\gamma\epsilon}, \end{equation}
because 
$\textrm {d}\boldsymbol {k}/\epsilon$ and 
$\textrm {d}\boldsymbol {p}/\gamma$ are Lorentz invariants (Landau & Lifshitz Reference Landau and Lifshitz1975).
3. Single Compton scattering algorithm
The implementation of single Compton scattering in a PIC code must not only recover the correct microphysics of the process (frequency shift, angle, momentum recoil) but must preserve the invariant number of collisions to obtain the correct scattering rates (Peano et al. Reference Peano, Marti, Silva and Coppa2009). The implementation follows naturally as each macro-particle represents an ensemble of real particles that are close neighbours in phase space. Each macro-particle has a weight 
$q$ that relates to the number of real particles it represents and, thus, samples a portion of the distribution function of real particles. Figure 2 outlines our implementation that follows three steps: (i) binning of the macro-particles into collision cells 
$\Delta \boldsymbol {x}$, a volume in configuration space, (ii) pairing of the colliding macro-particles according to their probability 
$P^{i,j}$ of interaction within 
$\Delta \boldsymbol {x}\Delta t$, (iii) update of the momenta of the scattering macro-particles.
Schematic of the Compton scattering algorithm. It follows three steps: (i) the macro-particles are binned into collision cells 
$\Delta \boldsymbol {x}$, (ii) the probability 
$P^{i,j}$ of interaction within 
$\Delta \boldsymbol {x}\Delta t$ is computed and scattering macro-particles are chosen using the no-time-counter method, (iii) the momenta of the scattering macro-particles are updated.
3.1. Macro-particles binning
The binning of macro-particles in collision cells naturally uses the single PIC cell as the smallest binning volume. The size of a PIC cell is also the smallest scale over which the self-consistent plasma collective fields are computed. For this reason, the collision cells are usually set equal to the PIC cells. Macro-photons and macro-electrons are binned in the collision cells and sorted such that we identify the indexes of macro-electrons and macro-photons within each collision cell.
3.2. Pairing
For each collision cell, we pair the scattering couples and add them to a scattering list using the no-time-counter (NTC) method (Bird Reference Bird1989; Abe Reference Abe1993). The NTC method is a popular Monte Carlo scheme for collision procedures involving single scattering events (not for cumulative scattering). The standard pairing routines for cumulative Coulomb collisions allow for a time step larger than the collision frequency, thus all macro-particles are involved in the scattering process each time step. Instead, the NTC method applies when the time step resolves the collision frequency such that the maximum possible number of macro-scatterings within a time step involves only a subset of all the macro-particles. Developed three decades ago (Bird Reference Bird1989), NTC provides a cost reduction for the sampling of a discrete probability distribution function. We detail now the NTC algorithm applied to single Compton scattering.
We consider a collision cell containing 
$N_{\omega }$ macro-photons and 
$N_e$ macro-electrons. A conservative upper-bound to the maximum probability of any macro-particle to collide within 
$\Delta t$ is
\begin{equation} P_{\textrm{max}} = 2 \sigma_T c \Delta t \max[q^i_e,q^j_{\omega}], \end{equation}
where 
$\max [q^i_e,q^j_{\omega }]$ is the largest weight with units of a density among all macro-particles in the collision cell (
$i\in [1,N_e]$ macro-electrons and 
$j\in [1,N_{\omega }]$ macro-photons). The factor 
$2$ appears conservatively as the upper bound in the relativistic transformation of the cross section 
$\sigma = \sigma _{T,O}\epsilon _O/\gamma \epsilon$, where 
$\max (\epsilon _O)=2\gamma \epsilon$. The maximum number of macro-particles that can scatter 
$N_{\textrm {max}}$ is given by the maximum probability 
$P_{\textrm {max}}$ times the number of all the possible unsorted pairing combinations 
$N_eN_{\omega }$ (potential scatterings) of the macro-photons with the macro-electrons. It reads
\begin{equation} N_{\textrm{max}}=P_{\textrm{max}} N_e N_{\omega}. \end{equation}
The number 
$N_{\textrm {max}}$ is usually not an integer and is rounded to the next or previous integer by a Monte Carlo sampling of the residue. This procedure preserves statistically the correct number of collisions within 
$\Delta \boldsymbol {x}\Delta t$. We randomly pair 
$N_{\textrm {max}}$ macro-photons and 
$N_{\textrm {max}}$ macro-electrons. This follows two steps: (i) the random sorting of the macro-photons and the macro-electrons, (ii) the selection of the first 
$N_{\textrm {max}}$ indexes. At this point, we have a shortlist of 
$N_{\textrm {max}}$ randomly paired macro-particles, which contains the maximum possible scatterings in the collision cell. For each pair in the short list, a random number 
$\textrm {rnd}\in [0,1]$ is rolled and compared with the joint probability
\begin{equation} P^{i,j}=\sigma(\boldsymbol{p}^i,\boldsymbol{k}^j) c \Delta t \max[q^i_e,q^j_{\omega}]/P_{\textrm{max}} \end{equation}
of scattering after having being selected within the 
$N_{\textrm {max}}$ pairs.
To compute the cross section 
$\sigma (\boldsymbol {p}^i,\boldsymbol {k}^j)$ we proceed as follows. The energy of the photon 
$\boldsymbol {k}^j$ is Lorentz boosted in the rest frame of the electron 
$\boldsymbol {p}^i$
\begin{equation} \epsilon_O^j=\gamma^i\epsilon^j-\hbar\boldsymbol{p}^i\cdot \boldsymbol{k}^j/m^2c^2. \end{equation}
Then, the cross section 
$\sigma _C(\epsilon _O^j)$ is computed in this frame and boosted back into the simulation frame using (2.13). A macro-electron/macro-photon pair from the shortlist is admitted to the scattering list based on a rejection method (accepted if 
$\textrm {rnd} < P^{i,j}$).
3.3. Momentum update
For each pair in the scattering list, the momenta are updated according to the Compton frequency shift and momentum recoil. The macro-photon four-wavevector 
$\mathcal {K}=(\epsilon , \hbar \boldsymbol {k}/mc)$ is Lorentz boosted in the rest frame of the electron, of momentum 
$\boldsymbol {p}$ in the simulation frame, as 
$\mathcal {K}_O=\boldsymbol{\mathsf{L}}(\boldsymbol {p})\mathcal {K}$, where the boost matrix is
\begin{equation} \boldsymbol{\mathsf{L}}(\boldsymbol{p}) = \left[ \begin{array}{cc} \gamma & -\boldsymbol{p}/mc \\ \displaystyle -\boldsymbol{p}^T/mc & \quad \boldsymbol{\mathsf{I}}+\boldsymbol{p}^T\boldsymbol{p}/m^2c^2(1+\gamma) \\ \end{array} \right] , \end{equation}
and 
$\boldsymbol{\mathsf{I}}$ the 
$3\times 3$ identity matrix. In the frame 
$O$, we identify the unit vector along the photon propagation direction 
$\hat {\boldsymbol {k}}_0$, which defines the symmetry axis for the scattering. We define an orthonormal unit vector base 
$\hat {\boldsymbol {e}}_1 = \hat {\boldsymbol {k}}_0, \hat {\boldsymbol {e}}_2\perp \hat {\boldsymbol {e}}_1, \hat {\boldsymbol {e}}_3=\hat {\boldsymbol {e}}_1\times \hat {\boldsymbol {e}}_2$. The two scattering angles 
$\theta$ and 
$\phi$ are then sampled, 
$\theta$ is the angle with respect to 
$\hat {\boldsymbol {e}}_1$ and 
$\phi$ is the angle on the plane 
$\hat {\boldsymbol {e}}_2,\hat {\boldsymbol {e}}_3$. This latter parameter, being the angle of rotational symmetry, is chosen randomly between 
$0$ and 
$2{\rm \pi}$. The angle 
$\theta$, or rather the parameter 
$\mu =\cos \theta$, is obtained by the inverse transform sampling method of the cumulative probability function given by the differential cross section of the process (see appendix A). We preferred this method rather than a rejection method, whose efficiency decreases for 
$\epsilon _O\gtrsim 1$ owing to the steepening of the probability density function close to 
$\mu \simeq -1$. The scattered photon energy is 
$\epsilon _O'$ given by (2.3)
\begin{equation} \epsilon_O' = \frac{\epsilon_O}{1+\epsilon_O(1-\mu)} \end{equation}and the scattered wavevector is
\begin{equation} \frac{\hbar\boldsymbol{k}_O'}{mc} = \epsilon_O'\left(\mu\hat{\boldsymbol{e}}_1+\sqrt{1-\mu^2}\cos\phi\hat{\boldsymbol{e}}_2 +\sqrt{1-\mu^2}\sin\phi\hat{\boldsymbol{e}}_3\right). \end{equation}
We transform back to the simulation frame 
$\mathcal {K}_O'=(\epsilon _O',\hbar \boldsymbol {k}_O'/mc)$ simply as 
$\mathcal {K}'=\boldsymbol{\mathsf{L}}(-\boldsymbol {p})\mathcal {K}_O'$, and by conservation of momentum the scattered electron has a new momentum 
$\boldsymbol {p}'=\boldsymbol {p}+\hbar (\boldsymbol {k}-\boldsymbol {k}')$.
3.4. Macro-particles with difference in weight
In PIC codes, it is unlikely that two scattering particles possess the same weight. Two main techniques to approach the problem have been discussed in Sentoku & Kemp (Reference Sentoku and Kemp2008). The first approach is based on a rejection method for which the scattering occurs with a probability based on the weights of the two-scattering macro-particles. This method does not reproduce the energy and momentum transfer of each collision but only on average, for a sufficiently high number of macro-particles in the collision cell. To preserve the energy and momentum transfer per collision, an alternative is first to spilt the scattering macro-particle of weight 
$q$ into a scattering fraction 
$q_s$ (
$\boldsymbol {p}_s$) and a non-scattering fraction
$q_{\textrm {ns}}$ (
$\boldsymbol {p}$). After the collision takes place, the two fractions are merged again into the macro-particle 
$q$, which has the average energy and momentum of 
$q_s$ and 
$q_{\textrm {ns}}$. This last method becomes inaccurate when 
$\boldsymbol {p}_s$ differs significantly from 
$\boldsymbol {p}$ such that the two fractions 
$q_s$ and 
$q_{\textrm {ns}}$ refer to two well-distinct portions of the phase space. This issue has already been addressed by Haugbølle (Reference Haugbølle2005) and relies on splitting and merging at two different steps. Here, we address this problem similarly but only the largest weight macro-particle is split before scattering. We briefly recall the main steps of the splitting procedure:
(i) identify the two scattering macro-particles of weight
$q_e^i$ and $q_{\omega }^j$, and select the largest weight between the two $\max [q_e^i,q_{\omega }^j]$;(ii) create a new particle of weight equal to
$\min [q_e^i,q_{\omega }^j]$;(iii) the two macro-particles of equal weight
$\min [q_e^i,q_{\omega }^j]$ are now paired and can be Compton scattered as described previously; and(iv) reassign to the split macro-particle the weight
$\max [q_e^i,q_{\omega }^j]-\min [q_e^i,q_{\omega }^j]$.
The splitting is performed within the scattering routine and can lead to a significant increase of macro-particles in the simulation. Merging algorithms (Vranic et al. Reference Vranic, Grismayer, Martins, Fonseca and Silva2015) can be used at a different step of the PIC loop to preserve the number of macro-particles in the simulation within a reasonable maximum, thus avoiding their exponential increase. The advantage of merging at a later stage is that only macro-particles close in phase space merge. Details can be found in Vranic et al. (Reference Vranic, Grismayer, Martins, Fonseca and Silva2015). Here we have not used the merging algorithm, but it can be straightforwardly combined with the algorithm described here.
4. Benchmarks
To benchmark our algorithm, we choose two problems that possess an exact analytical solution:
(i) the inverse Compton spectra produced by an electron scattering with an isotropic photon gas (Blumenthal & Gould Reference Blumenthal and Gould1970); and
(ii) the relaxation to the thermal equilibrium of a photon gas by Compton collisions with a thermal electron gas of fixed non-relativistic temperature described by the Kompaneets equation (Kompaneets Reference Kompaneets1957).
4.1. Inverse Compton spectra
Blumenthal & Gould (Reference Blumenthal and Gould1970) derived the inverse Compton spectra produced by the collision of a relativistic electron, 
$\gamma \gg 1$, with an isotropic gas of photons (see appendix B). The scattered photon distribution function reads
\begin{equation} f(\varGamma,\mathcal{E}') = 2q\log q +(1+2q)(1-q)+\frac{1}{2}\frac{\varGamma^2 q^2}{1+\varGamma q}(1-q), \end{equation}
where 
$q=\mathcal {E}'/[1+\varGamma (1-\mathcal {E}')]$, 
$ \, \mathcal {E}'=\epsilon '/\epsilon '_{\textrm {max}}$ is the scattered photon energy normalised to its maximum 
$\epsilon '_{\textrm {max}}=\gamma \varGamma /(1+\varGamma )$. The parameter 
$\varGamma =4\epsilon \gamma$ relates to the energy of the scattering photons in the electron rest frame and distinguishes two regimes: (i) Thomson limit 
$\varGamma \ll 1$ and (ii) extreme Klein–Nishina limit 
$\varGamma \gg 1$.
Figure 3 shows the excellent agreement between our simulations (dashed lines) and theory (4.1) (solid lines) for 
$\varGamma =0.1, 10, 100$. The scattered photon distribution function 
$f(\varGamma ,\mathcal {E}')$ is normalised 
$\int \,\textrm {d}\mathcal {E}' f(\varGamma ,\mathcal {E}') =1$. In our simulations, the photon gas is initialised with 
$1.5\times 10^7$ macro-photons, which mimic an emission line. All macro-photons have the same energy and are propagating in random directions, distributed uniformly on the surface of a sphere in momentum space. We considered the interaction at different photon energies 
$\epsilon =0.00025, 0.025, 0.25$. An equal number of macro-electrons is initialised at a Lorentz factor of 
$\gamma =100$, all collimated in one direction. To avoid that, each macro-photon scatters more than once the simulation runs for only a single time step where about 
$1\times 10^6$ macro-scatterings occur. The only constraint on 
$\Delta t$ is to be low enough such that 
$P_{\textrm {max}}<1$, to prevent multiple collisions of the photons to occur within a single time step.
Scattered photon distribution function 
$f(\varGamma ,\mathcal {E}')$ for 
$\varGamma =0.1, 10, 100$. The simulation results are shown by dashed lines and the theory (4.1) by solid lines (Blumenthal & Gould Reference Blumenthal and Gould1970).
4.2. Photon–electron gas equilibrium (Kompaneets equation)
Kompaneets addressed the thermodynamic equilibrium established between photons and free electrons if their interaction is only mediated by Compton scattering events (Kompaneets Reference Kompaneets1957). Kompaneets derived the partial differential equation that describes the temporal evolution of the photon occupation number 
$n$ resulting from the interaction with an electron gas of fixed non-relativistic temperature 
$k_BT\ll mc^2$, where 
$k_B$ is the Boltzmann constant. The full collision operator reads
\begin{equation} \frac{\partial n}{\partial t} = c\int\,\textrm{d}\boldsymbol{p} \frac{\textrm{d}\sigma}{\textrm{d}\varOmega}\left[f_e'n'(1+n)-f_en(1+n')\right], \end{equation}
where 
$f_e=f_e(\gamma )$ and 
$f_e'=f_e(\gamma ')$ refer to the electron energy distribution function evaluated at a Compton transition 
$\gamma mc^2 + \hbar \omega \rightleftharpoons \gamma 'mc^2+\hbar \omega '$. The evaluation of the photon occupation numbers 
$n=n(\omega )$ and 
$n'=n(\omega ')$ follow the same definition. The 
$n^2$ terms account for the photon Bose–Einstein statistics when phenomena such as stimulated scattering and superposition of states are considered. The full Boltzmann operator can be reduced to a Fokker–Planck form within the Thomson limit 
$\hbar \omega \ll mc^2$ (see appendix C). In regimes where the photon occupation number is small, 
$n\ll 1$, the photon electron gas interaction is mediated by single Compton scattering events and the linear Kompaneets equation in terms of the photon energy distribution function 
$f=\xi ^2 n$ reads
\begin{equation} \frac{\partial f}{\partial y} = \frac{\partial }{\partial \xi} \left[\xi^2 \frac{\partial f}{\partial \xi} +\left(\xi^2-2\xi\right) f \right], \end{equation}
where 
$y = t/t_C$, 
$t_C = mc/\sigma _Tn_ek_BT$ is the characteristic relaxation time and 
$\xi =\hbar \omega /k_BT$ is the photon energy normalised to the electron temperature.
Figure 4 shows the excellent agreement between our algorithm and the numerical solution of the linear Kompaneets equation (4.3) obtained with a finite-difference centred scheme. In our simulation, more than 
$10^5$ macro-photons are initialised to mimic an emission line at an average energy of 
$\bar {\xi }=\langle \xi \rangle =0.2$. The emission line has a small energy spread of 
$\sigma _{\xi }^2=\langle \xi ^2-\bar {\xi }^2\rangle =0.1$, and the initial distribution
\begin{equation} f(\xi, y=0) \propto \exp\left[-\frac{\left(\xi-\bar{\xi}\right)^2}{2\sigma_{\xi}^2}\right] \end{equation}
is Maxwellian. The same number of macro-electrons is sampled according to a Maxwellian distribution at a temperature of 
$k_BT=5\ \textrm {keV}$. To enforce a constant electron temperature during the simulation for a rigorous comparison with theory, we turn off the Lorentz force, which will arise from fluctuations in the electron density. We also omit the momentum, and energy ceded by the electron to the photons at each collision such that the electron population does not cool down. At 
$t\gtrsim 3\ t_C$, the photon energy distribution reaches equilibrium and converges towards the Wien's spectrum 
$f\propto \xi ^2\exp (-\xi )$, the correct equilibrium for the linear Kompaneets equation, as expected from the underlying hypothesis.
Time evolution of the photon distribution function 
$f(\xi )$ by the interaction with an electron gas of density 
$10^{18}\ \textrm {cm}^{-3}$ at 5 keV temperature for times 
$t=0, 0.5, 1, 3\ t_C$. After 
$t=3\ t_C$ the photon distribution function resembles the Wien spectrum and does not evolve significantly. Simulations in dashed lines and solution of the linear Kompaneets equation (4.3) in solid lines (Kompaneets Reference Kompaneets1957).
5. Considerations on the algorithm performance
In this section, we compare the computational cost of our Compton scattering algorithm with the standard PIC loop. The computational performance of our algorithm is usually dependent on the physical parameters of the particular simulated system. A thorough benchmark of its performance should then cover a variety of physical parameters of relevant case scenarios. In collisional plasmas, a typical benchmark of the performance of a collisional algorithm relies on the simulation of thermal plasma with and without collisions. Our choice for the comparison follows a similar criterion. We simulate in one dimension a thermal plasma in equilibrium with a photon gas, both at a temperature of 
$5\ \mathrm {keV}$. The plasma density is 
$n_p=10^{18}\ \mathrm {cm}^{-3}$ and the photon density is 
$n_{\omega }=3\times 10^{27}\ \mathrm {cm}^{-3}$, chosen such that the electron Compton collision frequency is a tenth of the plasma frequency 
$c \sigma _T n_{\omega }=\omega _p/10$. The electrons follow a Maxwell–Boltzmann distribution 
$f\propto \sqrt {W_k}\exp (-W_k/k_BT)$, and the photons follow a Wien distribution 
$f\propto W_k^2\exp (-W_k/k_BT)$. The computational domain is divided into 
$240$ cells. The time step is 
$\Delta t = 0.099\ \omega _p^{-1}$. Periodic boundary conditions are used. Figure 5 shows the time of a PIC loop per simulation particle with (
$\times$) and without (
$\circ$) Compton collisions as a function of the initial number of particles per cell (ppc). In this comparison macro-particle splitting occurs, but not particle merging. For a low number of ppc, the loop time is determined by the particle sorting routine (dependent also on the number of grid cells) used by the Compton module. For a high number of ppc, the scaling of our collision algorithm is proportional to the number of simulated particles. The computational cost of the sorting routine scales with both the number of simulated macro-particles and the number of cells in which they are sorted. Therefore, there is a trade-off, which for the set of parameters of these simulations occurs around 64 ppc and is highlighted by the dashed line in figure 5. Beyond the region delimited by the dashed line (which changed depending on the grid size and the number of ppc), the inclusion of the Compton algorithm does not affect significantly the standard PIC loop performance. This trade-off must be assessed for the different numerical parameters/configurations to determine the optimal performance conditions.
Time of a PIC loop per simulation particle with (
$\times$) and without (
$\circ$) Compton collisions as a function of the number of particles per cell (ppc). For a low number of ppc, the loop time is determined by the particle sorting routine (dependent also on the number of grid cells) used by the Compton module. For a sufficiently high number of ppc, the scaling of our collision algorithm is proportional to the amount of simulated particles.
6. Summary
We have presented a collision algorithm that incorporates the effect of single Compton scattering from high-frequency photons in PIC codes. This allows a self-consistent treatment of the high-frequency radiation coupling with the plasma dynamics from first principles. The algorithm shows excellent agreement with respect to the benchmarks: scattering photon spectrum from the collision with relativistic electrons (Blumenthal & Gould Reference Blumenthal and Gould1970) and the relaxation to thermal equilibrium of a photon population with an electron gas (Kompaneets Reference Kompaneets1957). This framework is at the forefront for the numerical modelling of photon–plasma interaction and opens new and exciting opportunities in the numerical investigation of plasma phenomena where a significant population of hard photons is present in the system.
Acknowledgements
This work was supported by the European Research Council (ERC-2015-AdG grant no. 695088), FCT (Portugal) grants SFRH/IF/01780/2013 and PD/BD/114323/2016 in the framework of the Advanced Program in Plasma Science and Engineering (APPLAuSE, FCT grant no. PD/00505/2012). Simulations were performed at IST cluster (Portugal) and at MareNostrum (Spain) under a PRACE award.
Editor Antoine C. Bret thanks the referees for their advice in evaluating this article.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Scattering angle by the Inverse Transform Sampling method
The probability distribution function (pdf) of the scattered macro-photon over the scattering angle 
$\mu =\cos \theta \in [1,-1]$ reads
\begin{equation} \mathrm{pdf}(\mu, \epsilon_O) = \frac{1}{\sigma(\epsilon_O)}\frac{\textrm{d}\sigma}{\textrm{d}\mu}, \quad \mathrm{and}\quad \int_1^{-1}\,\textrm{d}\mu \frac{\textrm{d}\sigma}{\textrm{d}\mu}=\sigma(\epsilon_O) \end{equation}with
\begin{equation} \frac{\textrm{d}\sigma}{\textrm{d}\mu}=-{\rm \pi} r_e^2\left(\frac{1}{1+\epsilon_O(1-\mu)}\right)^2\left(\frac{1}{1+\epsilon_O(1-\mu)}+\epsilon_O(1-\mu)+\mu^2\right). \end{equation}The cumulative distribution function (cdf) is
\begin{equation} \mathrm{cdf}(\mu, \epsilon_O) = \frac{1}{\sigma(\epsilon_O)}\int_1^{\mu}\,\textrm{d}\mu'\frac{\textrm{d}\sigma}{\textrm{d}\mu'} \end{equation}with
\begin{align} & \int_1^{\mu}\,\textrm{d}\mu'\frac{\textrm{d}\sigma}{\textrm{d}\mu'} = \frac{{\rm \pi} r_e^2}{\epsilon_O}\left\{\left( 1-\frac{2}{\epsilon_0}-\frac{2}{\epsilon_0^2}\right)\log\left[1+\epsilon_O(1-\mu)\right] \right.\nonumber\\ &\qquad \left. +\frac{1-\mu}{\epsilon_O}\left[1+\frac{1+2\epsilon_O}{1+\epsilon_O(1-\mu)}\right] +\frac{1}{2} -\frac{1}{[1+\epsilon_O(1-\mu)]^2} \right\} \end{align}
In the inverse transform sampling method a random number is generated in the range 
$\textrm {rnd}\in [0,1]$, then 
$\mu =\textrm {cdf}^{-1}(\textrm {rnd},\epsilon _O)$. Given the nonlinear dependence of the cdf on 
$\mu$, we use the bisection method to solve 
$\textrm {cdf}(\mu ,\epsilon _O)-\textrm {rnd}=0$.
Appendix B. Photon spectrum: single scattering with a relativistic electron
We briefly recall the main steps in the derivation of (4.1), see Blumenthal & Gould (Reference Blumenthal and Gould1970). If the photon gas is isotropic in the laboratory frame, it appears beamed at a small angle 
$\sim 1/\gamma$ in the proper frame 
$_O$ of reference of an incident relativistic electron 
$\gamma \gg 1$, as shown by (2.8). The Compton scattering differential rate in the laboratory frame reads
\begin{equation} \frac{\textrm{d}N_{\omega}}{\textrm{d}t\,\textrm{d}\epsilon'} = \int \textrm{d}\epsilon_O\ \int \textrm{d}\varOmega_O \frac{\textrm{d}N}{\textrm{d}t_O\,\textrm{d}\epsilon_O\,\textrm{d}\varOmega_O\,\textrm{d}\epsilon_O'} \frac{\textrm{d}t_O}{\textrm{d}t}\frac{\textrm{d}\epsilon_O'}{\textrm{d}\epsilon'}. \end{equation}
The time interval in the frame 
$_O$ is 
$\textrm {d}t_O= \textrm {d}t/\gamma$, and the energy transforms according to (2.10) as 
$\textrm {d}\epsilon ' \simeq \gamma (1-\cos \theta _O)\,\textrm {d}\epsilon _O'$. The Compton scattering differential rate in the frame 
$_O$ is
\begin{equation} \frac{\textrm{d}N}{\textrm{d}t_O\,\textrm{d}\epsilon_O\,\textrm{d}\varOmega_O\,\textrm{d}\epsilon_O'} =c \frac{\textrm{d}\sigma(\epsilon_O)}{\textrm{d}\varOmega_O}\delta(\epsilon_O'-\epsilon_O) \frac{\textrm{d}n_O}{\textrm{d}\epsilon_O}. \end{equation}
Here 
$\textrm {d}\sigma (\epsilon _O)/\textrm {d}\varOmega _O$ is the Klein–Nishina cross section. The photon density spectrum 
$\textrm {d}n_O/\textrm {d}\epsilon _O$ in the 
$_O$ frame can be related with the isotropic differential photon density in the laboratory frame by the Lorentz invariance of the ratio 
$\textrm {d}n/\epsilon$:
\begin{equation} \frac{1}{\epsilon_O}\frac{\textrm{d}n_O}{\textrm{d}\epsilon_O} = \frac{1}{\epsilon}\frac{\textrm{d}n}{\textrm{d}\epsilon}. \end{equation}
The isotropic differential photon density in the laboratory frame reads 
$\textrm {d}n=n(\epsilon ) \textrm {d}{\cos }\phi /2$, where 
$n(\epsilon )$ is the density of photons of a given energy 
$\epsilon$. According to (2.9), the incident angle in the laboratory frame results in a change in the photon energy in the 
$_O$ frame as 
$\vert \textrm {d}\epsilon _O/\textrm {d}{\cos }\phi \vert \simeq \gamma \epsilon$. One thus obtain from (B 3)
\begin{equation} \frac{\textrm{d}n_O}{\textrm{d}\epsilon_O} = \frac{\epsilon_O}{2\gamma \epsilon^2} n(\epsilon). \end{equation}By combining (B 4) and (B 2) with (B 1), the Compton scattering differential rate reads (Blumenthal & Gould Reference Blumenthal and Gould1970)
\begin{equation} \frac{\textrm{d}N}{\textrm{d}t\,\textrm{d}\mathcal{E}'} = \frac{3\sigma_T c}{4\gamma}\frac{n(\epsilon)}{\epsilon} f(\varGamma, \mathcal{E}'), \end{equation}
where 
$\mathcal {E}'=\epsilon '/\epsilon '_{\textrm {max}}$ is the scattered photon energy normalised to its maximum 
$\epsilon '_{\textrm {max}}=\gamma \varGamma /(1+\varGamma )$. The parameter 
$\varGamma =4\epsilon \gamma$ relates to the energy of the scattering photons in the electron rest frame and distinguishes two regimes: (i) Thomson limit 
$\varGamma \ll 1$ and (ii) extreme Klein–Nishina limit 
$\varGamma \gg 1$. The scattered photon distribution function reads
\begin{equation} f(\varGamma,\mathcal{E}') = 2q\log q +(1+2q)(1-q)+\frac{1}{2}\frac{\varGamma^2 q^2}{1+\varGamma q}(1-q), \end{equation}
where 
$q=\mathcal {E}'/[1+\varGamma (1-\mathcal {E}')]$.
Appendix C. Relaxation to thermal equilibrium of a photon gas: Kompaneets equation
We recall the main steps in the derivation of (4.3), see Kompaneets (Reference Kompaneets1957). In the Thomson limit 
$\hbar \omega \ll mc^2$, the energy exchange of one transition is small compared with the energy of the photon 
$\delta \omega =\vert \omega '-\omega \vert \ll \omega$. The energy exchange over one Compton event is
\begin{align} \hbar\delta\omega &= \hbar\omega\frac{c\boldsymbol{p}\cdot\left(\hat{\boldsymbol{k}}'-\hat{\boldsymbol{k}}\right)-\hbar\omega\left(1-\hat{\boldsymbol{k}}'\cdot\hat{\boldsymbol{k}}\right)}{\gamma mc^2+\hbar\omega\left(1-\hat{\boldsymbol{k}}'\cdot\hat{\boldsymbol{k}}\right)-c\boldsymbol{p}\cdot\hat{\boldsymbol{k}}} \end{align}
\begin{align} &\simeq \hbar\omega\left[\frac{\boldsymbol{p}}{mc}\cdot\left(\hat{\boldsymbol{k}}'-\hat{\boldsymbol{k}}\right)-\frac{\hbar\omega}{mc^2}\left(1-\hat{\boldsymbol{k}}'\cdot\hat{\boldsymbol{k}}\right)\right], \end{align}
where 
$\hat {\boldsymbol {k}}=\boldsymbol {k}/k$ and 
$\hat {\boldsymbol {k}}'=\boldsymbol {k}'/k'$ are the unit vectors that identify the photon propagation direction before and after scattering. In such regime, the functions 
$f_e'$ and 
$n'$ can be expanded to second order in the small parameter 
$\delta \omega$ allowing the reduction of the full collision operator (4.2) to a Fokker–Planck equation
\begin{align} \frac{\partial n}{c\partial t} &\simeq \left[\frac{\partial n}{\partial \xi} +n(1+n)\right] \int \textrm{d}\boldsymbol{p} \frac{\textrm{d}\sigma}{\textrm{d}\varOmega}f_e \frac{\hbar \delta\omega}{k_BT} \nonumber\\ &\quad +\left[\frac{\partial^2 n}{\partial \xi^2} +(1+n)\left(2\frac{\partial n}{\partial \xi} +n\right)\right] \int \textrm{d}\boldsymbol{p} \frac{\textrm{d}\sigma}{\textrm{d}\varOmega}f_e \left(\frac{\hbar \delta\omega}{k_BT}\right)^2, \end{align}
where the electron distribution function is assumed to be Maxwellian, and 
$\xi =\hbar \omega / k_BT$ is the energy of the photon normalised to the electron temperature. The expansion parameter 
$\delta \omega$ is small in the laboratory frame only if it is also small in the proper frame of each electron. This holds for non-relativistic electron temperatures 
$k_BT\ll mc^2$. The two integrals in 
$\delta \omega$ and in 
$\delta \omega ^2$ can be evaluated assuming the differential cross section in the Thomson limit
\begin{equation} \frac{\textrm{d}\sigma}{\textrm{d}\varOmega} = \frac{r_e^2}{2}\left( 1+ \cos^2\theta\right). \end{equation}
Then, the time evolution of the average occupation photon number 
$n$ reads
\begin{equation} \xi^2\frac{\partial n}{\partial y} = \frac{\partial}{\partial \xi} \left[\xi^4\left(\frac{\partial n}{\partial \xi} +n +n^2\right)\right], \end{equation}
where 
$y=t/t_C$ is the time normalised to 
$t_C=mc/\sigma _Tn_ek_BT$ and 
$n_e$ is the electron gas density. The time 
$t_C$ is the characteristic relaxation time of the process and the thermal equilibrium is reached when 
$y > 1$.
In regimes where the photon occupation number is small, 
$n\ll 1$, the photon electron gas interaction is mediated by single Compton scattering events and the equation reduces to its linear form
\begin{equation} \xi^2\frac{\partial n}{\partial y} = \frac{\partial}{\partial \xi} \left[\xi^4\left(\frac{\partial n}{\partial \xi} +n \right)\right]. \end{equation}
In terms of the photon energy distribution function 
$f=\xi ^2 n$, the linear Kompaneets equation reads
\begin{equation} \frac{\partial f}{\partial y} = \frac{\partial }{\partial \xi} \left[\xi^2 \frac{\partial f}{\partial \xi} +\left(\xi^2-2\xi\right) f \right]. \end{equation}
