Binary collisions

Section author: Pawel Ordyna

1 Introduction

The base PIC algorithm can’t resolve interactions between particles, such as coulomb collisions, on a scale smaller than one cell. Taking them into consideration requires additional simulation steps. In many PIC software solutions collisions are introduced as binary collisions. With that approach every macro particle collides only with one other macro particle so that the computation time scales linearly with the number of particles, and not quadratically.

We have used the binary collisions algorithm introduced in [1] and updated it according to the corrections suggested in [2]. In that way, we were recreating the implementation from smilei [3]. We have also tested our code against smilei by running the test cases from smilei’s documentation with PIConGPU. We describe the algorithm and the underlying equations in detail in section 3. Section 2 explains the usage of the collision feature. We show in section 4 the test simulations that we used to compare PIConGPU with smilei.

2 Usage

To enable collisions in your simulation you have to edit collision.param. There you can define the collision initialization pipeline. The collisions are set up with a Collider class. Each collider defines a list of species pairs that should collide with each other. A pair can consist of two different species, for collisions between two different particle groups, or two identical species, for collision within one group. The underlying collision algorithm is defined by a collision functor. At the moment, there are two collision functor classes available: RelativisticCollisionConstLog for running with a fixed coulomb logarithm, and RelativisticCollisionDynamicLog that calculates the coulomb log automatically. You can put in your initialization pipeline as many Collider objects as you need. Leaving the CollisionPipeline empty disables collisions.

2.1 Basics

Imagine you have two species in your simulation: Ions and Electrons. To enable collisions between the Ions and Electrons with a constant coulomb logarithm you would edit your param file in a following way:

namespace picongpu
{
    namespace particles
    {
        namespace collision
        {
            namespace precision
            {
                using float_COLL = float_64;
            } // namespace precision

            constexpr bool debugScreeningLength = false;
            using CollisionScreeningSpecies = MakeSeq_t<>;

            struct Params
            {
                static constexpr float_X coulombLog = 5.0_X;
            };
            using Pairs = MakeSeq_t<Pair<Electrons, Ions>>;

            /** CollisionPipeline defines in which order species interact with each other
             *
             * the functors are called in order (from first to last functor)
             */
            using CollisionPipeline = boost::mp11::
                mp_list<Collider<relativistic::RelativisticCollisionConstLog<Params>, Pairs> >;

            /** Chunk size used for cell list allocations.
             *
             * To reduce the fragmentation of the heap memory on accelerators the collision algorithm is allocating a
             * multiple of this value to store a cell list of particle IDs. The value must be non zero.
             */
            constexpr uint32_t cellListChunkSize = particles::TYPICAL_PARTICLES_PER_CELL;
        } // namespace collision
    } // namespace particles
} // namespace picongpu

Notice how the Coulomb logarithm is send to the functor class in a struct.

If you now would like to add internal collisions (electrons – electrons and ions – ions) you just need to extend the line 20 so that it looks like that:

using Pairs = MakeSeq_t<Pair<Electrons, Ions>, Pair<Electrons, Electrons>, Pair<Ions, Ions>>;

But what if you don’t want to have the same Coulomb logarithm for all collision types? For that you need more colliders in your pipeline. Here is an example with \(\Lambda = 5\) for electron-ion collisions and \(\Lambda=10\) for electron-electron and ion-ion collisions.

struct Params1
{
   static constexpr float_X coulombLog = 5.0_X;
};
struct Params2
{
   static constexpr float_X coulombLog = 10.0_X;
};
using Pairs1 = MakeSeq_t<Pair<Electrons, Ions>>;
using Pairs2 = MakeSeq_t<Pair<Electrons, Electrons>, Pair<Ions, Ions>>;
using CollisionPipeline =
    boost::mp11::mp_list<
        Collider<relativistic::RelativisticCollisionConstLog<Params1>, Pairs1>,
        Collider<relativistic::RelativisticCollisionConstLog<Params2>, Pairs2>
    >;

Automatic coulomb log calculation can be enabled for a collider by changing the collision functor. For example the previous setup with automatic calculation for the inter-species collisions would have the following CollisionPieline (and Params1 is not longer needed)

using CollisionPipeline =
    boost::mp11::mp_list<
        Collider<relativistic::RelativisticCollisionDynamicLog<>, Pairs1>,
        Collider<relativistic::RelativisticCollisionConstLog<Params2>, Pairs2>
    >;

The dynamic logarithm implementation uses a Debye length that is pre-calculated once for all colliders on each time step. So, whenever there is at least one collider with the RelativisticCollisionDynamicLog present in the CollisionPipeline this precalculation needs to be enabled by adding Species to the CollisionScreeningSpecies sequence. To include all species just set using CollisionScreeningSpecies = VectorAllSpecies;. But, this can be an arbitrary list of (filtered) species, see the CollisionsBeamRelaxation test for reference.

Note

The Debye length calculation requires at least 2 FieldTmp slots enabled in memory.param when CollisionScreeningSpecies has only one element and at least 3 otherwise.

2.2 Particle filters

You can also use particle filters to further refine your setup. The Collider class can take one more, optional, template argument defining a pair of particle filters. Each filter is applied respectively to the first and the second species in a pair. You need to define your filters in particleFilters.param and than you can use them, for example, like that:

using Pairs1 = MakeSeq_t<Pair<Electrons, Ions>>;
using Pairs2 = MakeSeq_t<Pair<Electrons, Electrons>, Pair<Ions, Ions>>;
using CollisionPipeline =
    boost::mp11::mp_list<
        Collider<
            relativistic::RelativisticCollisionConstLog<Params1>,
            Pairs1,
            FilterPair<filter::FilterA, filter::FilterB>>,
        Collider<
            relativistic::RelativisticCollisionConstLog<Params2>,
            Pairs2,
            OneFilter<filter::FilterA>>
        >;

Here only the electrons passing the A-filter will collide with ions but only with the ions that pass the B-filter. If the filters are identical you can use OneFilter instead of FilterPair. For collisions within one species the filters in FilterPair have to be identical since there is only one particle group colliding.

A full functional example can be found in the CollisionsBeamRelaxation test, where particle filters are used to enable each of the three colliders only in a certain part of the simulation box.

2.3 Precision

Highly relativistic particles can cause numerical errors in the collision algorithm that result in NaN values. To avoid that, by default, all the kinematics of a single binary collision is calculated in the 64 bit precision, regardless of the chosen simulation precision. Until now, this has been enough to avoid NaNs but we are looking into better solutions to this problem. You can change this setting by editing the

using float_COLL = float_64;

line. You can set it to

using float_COLL = float_X;

to match the simulation precision or

using float_COLL = float_32;

for explicit single precision usage. If you use PIConGPU with the 32 bit precision, lowering the collision precision will speed up your simulation and is recommended for non–relativistic setups.

2.4 Debug output

It is possible to write the average coulomb logarithm and s parameter (see 3.4 Details on the s parameter) values (averaged over all collisions in a time-step) for each collider. This debug output can be enabled per collider by setting the optional template parameter of the collision functor to true:

using CollisionPipeline =
    boost::mp11::mp_list<
        Collider<relativistic::RelativisticCollisionDynamicLog<true>, Pairs1>,
        Collider<relativistic::RelativisticCollisionConstLog<Params2, true>, Pairs2>
    >;

The debug info is written to a text file debug_values_collider_<collider index in the pipeline>_species_pair_<pair index in the list of pairs used with the collider>.dat The output file has three columns: iteration, coulomb log, s param. It it also possible to write out the precalculated Debye length averaged over all simulation cells by setting constexpr bool debugScreeningLength = true; The output is written to a file with two columns: iteration, Debye length [m]. The file name is average_debye_length_for_collisions.dat.

3 Algorithm

3.1 Algorithm overview

A short summary of the important algorithm steps in the case of inter-species collisions. The case of intra-collisions is very similar. See figures Fig. 2, Fig. 4, Fig. 3, Fig. 5 for more details.

1. Sort particles from a super cell into particle lists, one list for each grid cell.

2. In each cell, shuffle the list with more particles.

3. Collide each particle from the first longer list with a particle from the shorter one (or equally long). When you run out of particles in the shorter list, start from the beginning of that list and collide some particles more than once.

   1. Determine how many times the second particle will be collided with some particle from the longer list (in the current simulation step).

   2. Read particle momenta.

   3. Change into the center of mass frame.

   4. Calculate the \(s\) parameter.

   5. Generate a random azimuthal collision angle \(\varphi \in (0, 2\pi]\).

   6. Get the cosine of the 2nd angle \(\theta\) from its probability distribution (depends on \(s\)).

   7. Use the angles to calculate the final momenta (in the COM frame).

   8. Get the new momenta into the lab frame.

   9. Apply the new momentum to the macro particle A (smaller weighting).

      Do the same for the macro particle B (bigger weighting) but with a probability equal to the weighting ratio of the particles A and B.

4. Free up the memory used for the particle lists.

Fig. 2 Flow chart showing the complete algorithm. For more detail on intra-collisions see fig. Fig. 3, for more details on inter-collisions see fig. Fig. 4. Numbers in brackets refer to equations other to sections.

Fig. 3 Flow chart showing the part of the collision algorithm that is unique for intra-collisions. For more details on collisions functor see fig. Fig. 5 . Numbers in brackets refer to equations other to sections.

Fig. 4 Flow chart showing the part of the collision algorithm that is unique for inter-collisions. Numbers in brackets refer to equations other to sections.

Fig. 5 Flow chart showing the RelativisticBinaryCollision collisions functor. Numbers in brackets refer to equations other to sections.

3.2 Details on macro particle duplication

First step that requires some more detailed explanation is the step 3.1 . In a situation where there are less macro particles, inside one cell, of one species than the other one not every macro particle has its collision partner. Similar problem emerges in a case of intra-collisions when the particle number is odd. We deal with that issue using an approach introduced in [2]. We collide, in such situation, some macro particles more than once. To account for that, we use corrected particle weights \(w_{0/1} =\frac{1}{\max\qty{d_0, d_1}}\), where \(d_{0/1}\) are the number of collisions for the colliding macro particles.

Let us consider the inter-collisions first. The i–th particle from the longer list is collided with the (\(i \mod m)\) –th particle in the shorter one (\(m\) is the length of the shorter list). All of the particles from the longer list will collide just once. So the correction for each binary collision is \(1/d\) of the particle from the shorter list. \(d\) is determined in the following way:

d = floor(n / m);
if (i % m ) < (n % m) d = d + 1;

\(i\) – particle index in the long list, \(n\) – long list length, \(m\) – short list length, \(d\) – times the particle from the shorter list is used in the current step.

In the intra-collisions, the i–th (\(i\) is odd) particle collides with the \(i+1\)–th one. When there is, in total, an odd number of particles to collide, the first particle on the list collides twice. At first it is collided with the second one and in the end with the last one. All other particles collide once. So \(d\) will be 2 for the first collision (1st with 2nd particle) and for the last one (n-th with 1st particle). For the other collisions it’s 1.

3.3 Details on the coordinate transform

A binary collision is calculated in this model in the center of mass frame. A star \(^*\) denotes a COMS variable.

We use the coordinate transform from [1]:

(1)\[\mathbf{p}^* = \mathbf{p}_{\text{lab}} + \qty( \frac{\gamma_C -1}{\left|\mathbf{v}_C\right|^2} \mathbf{v}_C \cdot\mathbf{v}_{\text{lab}} - \gamma_C) m\gamma \mathbf{v}_C \ ,\]

where \(\mathbf{v}_C\) is the velocity of the CMOS in the lab frame, \(\gamma\) is the [list::duplications] factor in the lab frame, \(m\) the particle mass and \(\gamma_C\) the gamma factor of the CMOS frame.

(2)\[\mathbf{v}_C = \frac{\mathbf{p}_{\text{lab},0} + \mathbf{p}_{\text{lab}, 1}}{m_0\gamma_0 + m_1\gamma_1}\]

The inverse transformation:

(3)\[\mathbf{p_{\text{lab}}} = \mathbf{p}^* + \qty( \frac{\gamma_C -1}{\left|\mathbf{v}_C\right|^2} \mathbf{v}_C \cdot\mathbf{p}^* + m\gamma^* \gamma_C) \mathbf{v}_C \ ,\]

where

(4)\[\gamma^* = \gamma_C \gamma( 1 - \frac{\mathbf{v}_C \cdot \mathbf{v}_{\text{lab}}}{c^2}) \ .\]

3.4 Details on the \(s\) parameter

(5)\[s = \frac{1}{2} N \left<\theta^{*2}\right>\]

\(N\) is the number of real collisions. It’s the number of small angle collisions of a test particle represented by one of the macro particles with all the potential collision partners in a cell (here real particles not macro particles) in the current time step assuming the relative velocity is the one of the two colliding macro particles. \(\left<\theta^{*2}\right>\) is the averaged squared scattering angle for a single collision (of real particles). According to [2] \(s\) is a normalized path length.

To calculate this parameter we use the relativistic formula from [1] and adjust it so it fits the new corrected algorithm from [2].

(6)\[\begin{split}\begin{split} s_{01} =& \frac{\Delta T \log \Lambda q_0^2 q_1^2}{4\pi \varepsilon_0^2 c^4 m_0 \gamma_0 m_1 \gamma_1} \\ & \times \frac{\gamma_C\left|\mathbf{p}_0^*\right|}{m_0\gamma_0 + m_1 \gamma_1} \qty( m_0 \gamma_0^* m_1 \gamma_1^* c^2 \left|\mathbf{p}_0^*\right|^{-2} +1 )^2 \\ & \times N_{\text{partners}} V_{\text{cell}}^{-1} \max\qty{\frac{w_0}{d}, \frac{w_1}{d}} \ . \end{split}\end{split}\]

Here: \(\Delta T\) – time step duration, \(\log \Lambda\) – Coulomb logarithm, \(q_0,q_1\) – particle charges, \(\gamma_0, \gamma_1\) particles gamma factors(lab frame), \(N_{\text{partners}}\) is the number of collision partners (macro particles), \(V_{\text{cell}}\) – cell volume, \(w_0, w_1\) particle weightings, \(d\) was defined in 3.2.

For inter-species collisions \(N_{\text{partners}}\) is equal to the size of the long particle list. For intra-species collisions \(N_{\text{partners}}\) = \(n - 1 + \qty(n \mod 2)\),where \(n\) is the number of macro particles to collide.

The fact that \(s_{01}\) depends only on the higher weighting is accounted for by the rejection method in the 3.9 step.

3.4.1 Low temperature limit

According to [1] equation (6) will provide non physical values for low temperatures. More specifically, it will result in \(s\) values corresponding to scattering lengths smaller than the average particle distance \(\qty(\frac{V}{n})^{\frac{1}{3}}\). [1] provides a maximal value for \(s_{01}\):

(7)\[\begin{split}\begin{split} s_{01}^{\max} =& \qty(\frac{4\pi}{3})^{1/3} \frac{\Delta T\qty(m_0 + m_1)}{\max \qty{m_0 n_0^{2/3}, m_1 n_1^{2/3}}} \mathbf{v}_{\text{rel}}^* \\ & \times N_{\text{partners}} V_{\text{cell}}^{-1} \max\qty{\frac{w_0}{d}, \frac{w_1}{d}} \ . \end{split}\end{split}\]

with

(8)\[\mathbf{v}^*_{\text{rel}} = \frac{\qty(m_1\gamma_1 + m_2\gamma_2)p_1^*}{m_1\gamma_1^*m_2\gamma_2^*\gamma_C} \ .\]

where the relativistic factor \(\qty(1 + v_1^*v_2^*/c^2)^{-1}\) has been left out.

For each binary collision both values are calculated and the smallest one is used later. The particle density is just the sum of all particle weightings from one grid cell divided by cell volume

(9)\[n = \frac{1}{V_{\text{cell}}} \sum_{i} w_i \ .\]

Note

It is not checked if the collision is really non-relativistic. If the low temp limit is smaller than \(s_{01}\) due to some other reason, e.g. an overflow in \(s_{01}\) calculation, the code will use this limit regardless of the particle being relativistic or not which could be physically incorrect.

3.5 Details on the scattering angle distribution

The distribution for the cumulative angle \(\chi\) as a function of \(s\) was introduced in [4]

(10)\[F(\chi) = \frac{A(s) \sin \chi}{2 \sinh A(s)} e^{A(s) \cos \chi} \ .\]

We obtain a random value for the cosine from \(F\) with

(11)\[\cos \chi = A^{-1} \ln(e^{-A} + 2U \sinh A) \ ,\]

where \(U\) is a random float between 0 and 1. The parameter \(A\) is obtained by solving

(12)\[\coth A - A^{-1} = e^{-s} \ \ .\]

Previously the algorithm was approximating \(A\) with a polynomial fit from [1]. Now the \(\cos \chi\) is obtained from a new fit that was introduced in smilei:

If \(\mathbf{ s < 4}\) then:

(13)\[\alpha = 0.37s - 0.005 s^2 - 0.0064 s^3 \ .\]

(14)\[\sin^2(x/2) = \frac{\alpha U}{\sqrt{(1-U) + \alpha^2 U}}\]

(15)\[\cos(x) = 1 - 2 \sin^2(x/2)\]

In the \(s\rightarrow \infty\) limit scattering becomes isotropic [4] so that we can take \(\cos \chi = 2U -1\) for \(s > 4\).

3.6 Details on final momentum calculation

The final particle momenta in the COMS frame are calculated with the following formula from [1]

(16)\[\begin{split}\mathbf{p}_{1f}^* = -\mathbf{p}_{2f}^* = \begin{pmatrix} \frac{p_{1x}^*p_{1z}^*}{p_{1\perp}^*} & \frac{p_{1y}^*p_{1}^*}{p_{1\perp}^*} & p_{1x}^* \\ \frac{p_{1y}^*p_{1z}^*}{p_{1\perp}^*} & \frac{p_{1x}^*p_{1}^*}{p_{1\perp}^*} & p_{1y}^* \\ -p_{1\perp}^* & 0 & p_{1z}^* \end{pmatrix} \cdot \begin{pmatrix} \sin\theta^*\cos\varphi^* \\ \sin\theta^*\sin\varphi^* \\ \cos\theta^* \end{pmatrix} \ .\end{split}\]

3.7 Dynamic Coulomb logarithm calculation

With the RelativisticCollisionDynamicLog functor the Coulomb logarithm is calculated individually for each collision following a formula from [1]:

\[\ln \Lambda = \max \qty[2, \frac{1}{2}\ln\qty(1 + \frac{\lambda_D^2}{b_\text{min}^2} ) ] \ ,\]

where \(b_\text{min}\) is a minimal impact parameter that depends on particle momenta, charges, and masses; and \(\lambda_D\) is the Debye length.

Please note, according to the smilei documentation, in the equation (22) in [1] for \(b_\text{min}\) the last factor should not be squared; we drop the square also in PIConGPU.

The Debye length is calculated once per time-step for each simulation cell using the formula:

\[\lambda_D^{-2} = \frac{1}{\epsilon_0} \sum_{\alpha} n_\alpha \left<q_\alpha\right>^2 / T_\alpha \ ,\]

where the sum goes over all charged particle species, \(n\) is the number density, \(\left<q\right>\) is the average charge, and \(T\) is the temperature. The temperature is assumed to be equal to \(\frac{2}{3}\left<E_\text{kin}\right>\).

In PIConGPU the contributions from each species are calculated as

\[\frac{2}{3 \epsilon_0} \rho^2 \varepsilon^{-1} \ ,\]

where \(\rho\) is the charge density and \(\varepsilon\) is the energy density. It can be shown that this is equal to \(\frac{1}{\epsilon_0} n <q>^2 / T\).

Additionally \(\lambda_D\) is cut-off at the mean interatomic distance of the species with the highest density:

\[\lambda_D \geq \qty(4\pi n_\text{max}/3)^{-1/3}\]

4 Tests

For testing we plan to reproduce all the test cases from smilei’s documentation( https://smileipic.github.io/Smilei/collisions.html). For now we have done the thermalization and the beam relaxation tests. The simulations that we used are available under share/picongpu/tests.

4.1 Thermalization

In this example there are two particle populations — electrons and ions. They are thermally initialized with different temperatures and their temperatures get closer to each other with time. The usual PIC steps are disabled (there is no field solver and no pusher). The thermalization happens solely due to the binary collisions. We enable inter-collisions for ions and electrons as well as collisions between the two species. Simulation parameters are listed in table Table 2. The species temperatures are calculated in post processing from an openPMD output, a python script used for the analysis can be found in the tests directory. We extended the tests by running the simulation for longer and also by including runs using the automatic coulomb log calculation instead of a constant value.

The results from the original setup with constant logarithm are shown in fig. Fig. 6, Fig. 7, Fig. 8 for three different macro particle weight ratios. The figures Fig. 9, Fig. 10, Fig. 11 show the results from a set-up with dynamic coulomb logarithm calculation for the electron-ion collisions, and fig. Fig. 12, Fig. 13, and Fig. 14 come from a set-up with automatic logarithm for all three collision pairs. The theoretical curves are obtained from the same formula that was used by smilei’s developers and originates from the NRL plasma formulary [5].

\[\frac{\mathrm{d}T_\alpha}{\mathrm{d}t} = \nu_\epsilon(T_\beta -T_\alpha)\]

\[\nu_\epsilon = \frac{2}{3}\sqrt\frac{2}{\pi} \frac{e^4\,Z^{\star 2} \sqrt{m_em_i}\,n_i\,\ln\Lambda }{ 4 \pi\varepsilon_0^2 \,\qty(m_eT_e+m_iT_i)^{3/2} }\]

Since the collisions in different cells are independent of each other, one can treat each cells as an individual randomized run. The simulation values are obtained by averaging over the individual simulation cells. The upper right panel shows the average values together with \(\pm 2\) standard deviation of the distribution, while the left panel is showing the same values with \(\pm 2\) standard deviation of the mean. The bottom panel is just a zoom in on the upper left panel. The inputs for the smilei runs can be found in smilei_thermalization_inputs.zip.

Additionally, figures Fig. 6, Fig. 7, Fig. 8 show the average values of the coulomb logarithm, the s paramter, and the debye length. The coulomb logarithm theoretical value (blue curve) is calculated with the following formula from [5]:

\[\lambda_{ei} = 24 - \ln(n_e^{1/2}T_e^{-1})\]

Fig. 6 Electron (blue) and ion (red) temperature over time in the thermalization test with both \(\ln\Lambda\) constant. The electron to ion weight ratio in the simulation is 1:5. Black lines are the the theoretical curves.

Fig. 7 Electron (blue) and ion (red) temperature over time in the thermalization test with both \(\ln\Lambda\) constant. The electron to ion weight ratio in the simulation is 1:1. Black lines are the the theoretical curves.

Fig. 8 Electron (blue) and ion (red) temperature over time in the thermalization test with both \(\ln\Lambda\) constant. The electron to ion weight ratio in the simulation is 5:1. Black lines are the the theoretical curves.

Fig. 9 Electron (blue) and ion (red) temperature over time in the thermalization test with a dynamic \(\ln\Lambda\) for the electron-ion collisions The electron to ion weight ratio in the simulation is 1:5. Black lines are the the theoretical curves.

Fig. 10 Electron (blue) and ion (red) temperature over time in the thermalization test with a dynamic \(\ln\Lambda\) for the electron-ion collisions. The electron to ion weight ratio in the simulation is 1:1. Black lines are the the theoretical curves.

Fig. 11 Electron (blue) and ion (red) temperature over time in the thermalization test with a dynamic \(\ln\Lambda\) for the electron-ion collisions. The electron to ion weight ratio in the simulation is 5:1. Black lines are the the theoretical curves.

Fig. 12 Electron (blue) and ion (red) temperature over time in the thermalization test with a dynamic \(\ln\Lambda\) for all collisions. The electron to ion weight ratio in the simulation is 1:5. Black lines are the the theoretical curves.

Fig. 13 Electron (blue) and ion (red) temperature over time in the thermalization test with a dynamic \(\ln\Lambda\) for all collisions. The electron to ion weight ratio in the simulation is 1:1. Black lines are the the theoretical curves.

Fig. 14 Electron (blue) and ion (red) temperature over time in the thermalization test with a dynamic \(\ln\Lambda\) for all collisions. The electron to ion weight ratio in the simulation is 5:1. Black lines are the the theoretical curves.

Fig. 15 Collision debug values from the thermalization test with both \(\ln\Lambda\) constant. The electron to ion weight ratio in the simulation is 1:5.

Fig. 16 Collision debug values from the thermalization test with both \(\ln\Lambda\) constant. The electron to ion weight ratio in the simulation is 1:1.

Fig. 17 Collision debug values from the thermalization test with both \(\ln\Lambda\) constant. The electron to ion weight ratio in the simulation is 5:1.

Fig. 18 Collision debug values from the thermalization test with a dynamic \(\ln\Lambda\) for the electron-ion collisions The electron to ion weight ratio in the simulation is 1:5.

Fig. 19 Collision debug values from the thermalization test with a dynamic \(\ln\Lambda\) for the electron-ion collisions. The electron to ion weight ratio in the simulation is 1:1.

Fig. 20 Collision debug values from the thermalization test with a dynamic \(\ln\Lambda\) for the electron-ion collisions. The electron to ion weight ratio in the simulation is 5:1.

Fig. 21 Collision debug values from the thermalization test with a dynamic \(\ln\Lambda\) for all collisions. The electron to ion weight ratio in the simulation is 1:5.

Fig. 22 Collision debug values from the thermalization test with a dynamic \(\ln\Lambda\) for all collisions. The electron to ion weight ratio in the simulation is 1:1.

Fig. 23 Collision debug values from the thermalization test with a dynamic \(\ln\Lambda\) for all collisions. The electron to ion weight ratio in the simulation is 5:1.

Table 2 Simulation parameters in the thermalization test

parameter or setting	value
time step duration	2/3 fs
time steps in the simulation	100
density profile	homogeneous
density	1.1 × 1028 m-3
cell side length	\(\frac{1}{3}c \cdot 10^{-13} \approx 10 \mu m\)
ion mass	\(10 \ m_e\)
ion charge	+1
initial ion temperature	1.8 × 10−4 \(m_e c^2\)
initial electron temperature	2.0 × 10−4 \(m_e c^2\)
Coulomb logarithm (inter–collisions)	5 (if not dynamic)
Coulomb logarithm (intra–collisions)	1000 (if not dynamic)
geometry	2D
grid	12x12
super cell size	4x4
macro particles per cell (ions) setups 1, 2, 3	5000, 1000, 5000
macro pearticles per cell (electrons) setups 1, 2, 3	5000, 5000, 1000

4.2 Beam relaxation

A population of electrons with a very small temperature and a drift velocity (the beam) is colliding with ions. Due to the collisions the velocity distribution of electrons is changing and the drift momentum is transferred into the electron transversal momentum and partially into ion momenta. In this test only the inter-collisions (between ions and electrons) are enabled.

There are three slightly different setups with varying electron drift velocity, ion charge and time step duration. Additionally each setup performs the collisions with three different electron to ion weight ratios: 1:1, 5:1, 1:5. This is achieved by dividing the simulation box into three parts and enabling collisions only for one ratio in each part. All important simulation parameters can be found in tables Table 3 and Table 4. This test was also extended with runs utilizing the automatic Coulomb logarithm calculation.

The following figures show the electron and ion drift velocities \(\left<v_x\right>\), electron transversal velocity \(\sqrt{\left< v_\perp^2\right>}\), as well as the ion drift velocity, developing over time. All figures include a comparison values from smilei simulations, the smilei inputs can be found in smilei_beam_relaxation_inputs.zip.

Fig. 24 Electron drift velocity \(\left<v_x\right>\), electron transversal velocity \(\sqrt{\left< v_\perp^2\right>}\), and ion drift velocities from the beam equilibration example setup 1 (with a constant Coulomb logarithm).

Fig. 25 Electron drift velocity \(\left<v_x\right>\), electron transversal velocity \(\sqrt{\left< v_\perp^2\right>}\), and ion drift velocities from the beam equilibration example setup 2 (with a constant Coulomb logarithm).

Fig. 26 Electron drift velocity \(\left<v_x\right>\), electron transversal velocity \(\sqrt{\left< v_\perp^2\right>}\), and ion drift velocities from the beam equilibration example setup 3 (with a constant Coulomb logarithm).

Fig. 27 Electron drift velocity \(\left<v_x\right>\), electron transversal velocity \(\sqrt{\left< v_\perp^2\right>}\), and ion drift velocities from the beam equilibration example setup 1 (with dynamic Coulomb logarithm).

Fig. 28 Electron drift velocity \(\left<v_x\right>\), electron transversal velocity \(\sqrt{\left< v_\perp^2\right>}\), and ion drift velocities from the beam equilibration example setup 2 (with dynamic Coulomb logarithm).

Fig. 29 Electron drift velocity \(\left<v_x\right>\), electron transversal velocity \(\sqrt{\left< v_\perp^2\right>}\), and ion drift velocities from the beam equilibration example setup 3 (with dynamic Coulomb logarithm).

The debug average quantities for \(\ln\Lambda, s, \lambda_D\) are shown in figures … .

Fig. 30 Average Coulomb logarithm \(\ln\Lambda\), \(s\) (proportional to collision frequency), Debye length \(\lambda_D\) from the beam equilibration example setup 1 (with a constant Coulomb logarithm).

Fig. 31 Average Coulomb logarithm \(\ln\Lambda\), \(s\) (proportional to collision frequency), Debye length \(\lambda_D\) from the beam equilibration example setup 2 (with a constant Coulomb logarithm).

Fig. 32 Average Coulomb logarithm \(\ln\Lambda\), \(s\) (proportional to collision frequency), Debye length \(\lambda_D\) from the beam equilibration example setup 3 (with a constant Coulomb logarithm).

Fig. 33 Average Coulomb logarithm \(\ln\Lambda\), \(s\) (proportional to collision frequency), Debye length \(\lambda_D\) from the beam equilibration example setup 1 (with dynamic Coulomb logarithm).

Fig. 34 Average Coulomb logarithm \(\ln\Lambda\), \(s\) (proportional to collision frequency), Debye length \(\lambda_D\) from the beam equilibration example setup 2 (with dynamic Coulomb logarithm).

Fig. 35 Average Coulomb logarithm \(\ln\Lambda\), \(s\) (proportional to collision frequency), Debye length \(\lambda_D\) from the beam equilibration example setup 3 (with dynamic Coulomb logarithm).

Table 3 Collisions in the 3 parts of the simulation box in the beam relaxation example

parameter	upper part	middle part	lower part
macro particles per cell (ions)	1000	1000	100
macro particles per cell (electrons)	1000	100	1000

Table 4 Simulation parameters in beam the relaxation test

parameter or setting	value
	setup 1	setup 2	setup 3
time step duration	\(\frac{2}{3}\) fs	\(\frac{0.01}{3}\) fs	\(\frac{0.002}{3}\) fs
time steps in the simulation	200
density profile	homogeneous
density electrons	1.1 × 1028 m-3
density ions	1.1 × 1028 m-3	1.1 × 1028 m-3	3.7 × 1027 m-3
cell side length	\(\frac{1}{15}c \cdot 10^{-13} \approx 2\mu m\)
ion mass	\(10 \ m_e\)
ion charge	+1	+1	+3
initial electron drift	\(0.05c\)	\(0.01c\)	\(0.01c\)
initial ion temperature	0.00002 \(m_e c^2\)
initial electron temperature	0.0000002 \(m_e c^2\)
Coulomb logarithm	5
geometry	2D
grid	12x12
super cell size	4x4

References

[1]F. Pérez, L. Gremillet, A. Decoster, M. Drouin, and E. Lefebvre, Improved modeling of relativistic collisions and collisional ionization in particle-in-cell codes, Physics of Plasmas 19, 083104 (2012).

[2]D. P. Higginson, I. Holod, and A. Link, A corrected method for Coulomb scattering in arbitrarily weighted particle-in-cell plasma simulations, Journal of Computational Physics 413, 109450 (2020).

[3]J. Derouillat, A. Beck, F. Pérez, T. Vinci, M. Chiaramello, A. Grassi, M. Flé, G. Bouchard, I. Plotnikov, N. Aunai, J. Dargent, C. Riconda, and M. Grech, SMILEI: A collaborative, open-source, multi-purpose particle-in-cell code for plasma simulation, Computer Physics Communications 222, 351 (2018).

[4]K. Nanbu, Theory of cumulative small-angle collisions in plasmas, Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 55, 4642 (1997).

[5]A. S. Richardson, NRL Plasma Formulary, (2019).