Betaboltz: a Monte-Carlo simulation tool for gas scattering processes

TL;DR

Betaboltz is an open-source C++17 simulation tool for modeling electron and ion transport in gas mixtures under electric and magnetic fields, validated against existing drift velocity data.

Contribution

It introduces a flexible, open-source Monte Carlo simulation framework for gas scattering processes with user-defined mixtures and external fields, integrated with LXCat cross-section data.

Findings

Validated against drift velocity tables with acceptable accuracy.

Supports user-defined gas mixtures and external fields.

Open-source and easily integrable as a shared library.

Abstract

We present an open-source code for the simulation of electron and ion transport for user-defined gas mixtures with static uniform electric and magnetic fields. The program provides microscopic interaction simulation and is interfaced with cross-section tables published by LXCat[1]. The framework was validated against drift velocity tables available in literature obtaining an acceptable match for atomic and non-polar molecular gases with spherical symmetry. The code is written in C++17 and is available as a shared library for easy integration into other simulation applications.

Tables

Table 1. Table 1 : This is the list of the classes extending the class BaseBulletLimiter . If you need a behavior not listed here, it is possible to create your own limiter extending the class BaseBulletLimiter and implementing the method isOver() .

Class Name	Description
TimeBulletLimiter	Limit the simulation to a certain duration.
DistanceBulletLimiter	Destroy any particle which goes too distant regards to a given point.
EnergyBulletLimiter	Destroy all the particles when the energy goes outside the given energy range.
ChildrenBulletLimiter	Limit the number of particles in the simulation
InteractionBulletLimiter	Limit the number of interaction a particle can have in a simulation
OutOfDetectorBulletLimiter	Destroy any particle which goes outside the detector. This is the default if no other limiter is specified.

Table 2. Table 2 : Table representing the base handlers methods which are called during a simulation.

Method	Description
onRunStart	This method is called when the simulation starts.
onRunEnd	This method is called when the simulation ends.
onEventStart	This method is called on the beginning of every event.
onEventEnd	This method is called at the end of every event.
onBulletCreate	This method is called when a new bullet is created (i.e. by an ionization).
onBulletStep	This method is called between two collisions. The ’from’ prefix specify the state just after the last collision and the ’to’ the state just before the current collision.
onBulletCollision	This method is called at every collision. The ’before’ prefix specify the bullet state just before the collision, the ’after’ just after.
onBulletDestroy	This method is called when a bullet is destroyed (i.e. by an attachment).

Table 3. Table 3 : List of commands used to generate the plots presented in this article. < F >,< Z >,< R1 >,< R2 > represents placeholders for the reduced field [ Td / divide townsend absent \mathrm{Td}\text{/} ], the overhead and two gas ratios, respectively. The flag -- disable - check - monoatomic was omitted for some entries to improve readability. The cross-section tables selected here were used for our simulations, but no serious analysis was performed to determine if this section produces the best fit against experimental data.

Figures	Command
6(a), 8(a), 7(a), 9(a)	drifter -g Ar \ -d ’ela:BSR/BSR-500|ine:Biagi8’ \ \ \ -e 100 -n 20 -i 250000 --reduced-field <F>
6(b)	drifter -g Kr \ -d ’ela:COP/Gr1|ine:Biagi8’ \ \ \ \ \ \ \ -e 100 -n 20 -i 250000 --reduced-field <F>
6(c), 9(b), 5	drifter -g CF4 -d ’Bordage’ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ -e 100 -n 20 -i 250000 --reduced-field <F>
6(d), 7(d), 10(a)	drifter -g CH4 -d ’Morgan’ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ -e 100 -n 20 -i 250000 --reduced-field <F>
6(e), 8(b), 7(c), 5	drifter -g CO2 -d ’ela:Phelps|ine:Biagi8’ \ \ \ \ \ \ \ \ -e 100 -n 20 -i 250000 --reduced-field <F>
6(f)	drifter -g NH3 -d ’Morgan’ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ -e 100 -n 20 -i 250000 --reduced-field <F>
7(b), 5	drifter -g Xe \ -d ’ela:cop/Gr1|ine:Biagi8’ \ \ \ \ \ \ \ -e 100 -n 20 -i 250000 --reduced-field <F>
4	drifter -g CO2 -d ’ela:Phelps|ine:Biagi8’ --reduced-field 15 --trial-algo "variable(overhead:<Z>)"
12	drifter -g Ar CO2 -x <R1> <R2> \ \ \ \ \ \ \ \ \ \ \ --reduced-field <F> -t 297.2
	\ \ \ \ \ \ \ \ -d ’ela.el:BSR/BSR-500|ela.mt:Biagi8|ine:Biagi8’ ’ela.el:Phelps|ela.mt:Biagi8|ine:Biagi8’

Equations

$$\ln \left(\frac{1}{1-R}\right)={\int }_0^{t_c}\nu \left(t\right)dt$$

$$\nu \left(\epsilon \right)=\sqrt{\frac{2\epsilon }{m}}\sum _{k=1}^N{n}_k\left({\sigma }_k\left(\epsilon \right)+\sum _{i=1}^{J_k}{\sigma }_{ki}\left(\epsilon \right)\right)$$

$$t_c=-\frac{\ln (1-R)}{{\nu }^{\prime }}$$

$$R>\nu (t)/{\nu }^{\prime }$$

$$\cos (\theta )=1-\frac{2R}{1+8\epsilon (1-R)}$$

$$\cos (\theta )=1-\frac{2R(1-\xi )}{1+\xi (1-2R)}$$

$$\frac{{\sigma }_{mt}(\epsilon )}{{\sigma }_{el}(\epsilon )}=\frac{1-\xi }{2{\xi }^2}\left((1+\xi )\ln \frac{1+\xi }{1-\xi }-2\xi \right)$$

$$\frac{{\epsilon }_f}{{\epsilon }_i}=1-\frac{2mM(1-\cos \theta )}{(m+M{)}^2}$$

$$\frac{{\epsilon }_f}{{\epsilon }_i}=1-\frac{M}{m+M}\frac{{\epsilon }_k}{{\epsilon }_i}+\frac{2mM}{(m+M{)}^2}\left(\sqrt{1-\frac{m+M}{M}\frac{{\epsilon }_k}{{\epsilon }_i}}\cos \theta -1\right)$$

$${\text{V}}_1^{\prime }$$

$${\text{V}}_2^{\prime }$$

$$K$$

$$A$$

$$\sigma (\epsilon )$$

$$W=\frac{1}{N}\sum _n^N\frac{z_n}{t_n}$$

$$⟨\alpha ⟩$$

$$⟨\eta ⟩$$

$$⟨\epsilon ⟩$$

$$⟨{\epsilon }^{\prime }⟩$$

$$W$$

$${\text{d}}_n$$

$$D_{}$$

$$D_L$$

$$D_T$$

$${\sigma }_{}$$

$${\sigma }_L$$

$${\sigma }_T$$

$$\eta =\frac{{\nu }_{real}}{{\nu }_{real}+{\nu }_{null}+{\nu }_{fail}}$$

$${\nu }^{\prime }=\left(1+\zeta \right)\cdot \nu (\epsilon )$$

Full text

††journal: Computer Physics Communications

PROGRAM SUMMARY

Program Title: Betaboltz

Licensing provisions: LGPL v3

Programming language: C++17

Nature of problem: Simulations of electron and ion transport in arbitrary gas mixture under static uniform electric and magnetic fields.

Solution method: Particle motion using classical and relativistic equation via interaction sampling using Monte-Carlo techniques.

Additional comments including Restrictions and Unusual features: At the time of writing only static uniform electromagnetic fields are supported. However, the implementation of arbitrary fields can be added given an analytical solution is available. A custom XML format for cross-section was developed, because full compatibility with LXCat [1] XML format was not possible. Cross-section databases in the new format are available in the download section of the LXCat site [2].

1 Introduction

A comprehensive understanding of the charged particle transport in low-temperature gas mixtures, has crucial importance for the design and simulation of gas-based radiation detectors. The modeling of these processes can be performed by numerically solving the Boltzmann equation [1, 2], or by sampling the motion of charged particles using Monte-Carlo techniques [3]. Both approaches require the knowledge of the possible interaction types for all the components of the gas mixture. This information can be obtained using theoretical quantum models or by measurements of the electron/ion swarm experiments.

When performing any calculation involving particle transport in low-temperature gases, we have to choose which cross-section databases to use in the calculations. While, in the past, these tables had to be searched within the extensive available literature, now it is possible to access publicly available database such as LXCat [4], collecting cross-section tables, drift velocities, and other swarm attributes, for an extensive set of gases.

Several Boltzmann and Monte-Carlo simulation frameworks exist, but none of them can perform full detector simulations or have features like multi-thread execution. A well-known and widespread solution is Magboltz [5], an electron drift solver which contains the cross-section tables calculated by Biagi et al. [3] and allows fast determination of the electron swarm properties for a reasonable set of gases. In the same class of solvers, we can mention BOLSIG+ [6], a well-known Boltzmann solver and METHES [7], which allows Monte-Carlo simulations for an arbitrary static electric field using cross-sections tables available from LXCat. In addition to the previous solutions, we should mention PyBoltz [8], a Cython porting of the original Magboltz code and pyMETHES [9], a Python port of the METHES package.

In this paper, we present a flexible and robust solution for the simulation of particle drift in a custom detector’s using Monte Carlo methods. We will present here a list of the main characteristics of the simulation tool that we have developed:

C++ shared library interface-able with existing codes.

• 2.

Multi-thread execution.

• 3.

Access to cross-section tables from LXCat.

• 4.

Possibility to simulate detector composed by multiple chambers with different electromagnetic fields and gases.

• 5.

Supports arbitrary static uniform electromagnetic fields.

• 6.

Dimensional check at compile time.

Our program’s most important feature is that we do not perform any calculation of the swarm attribute: the software will calculate the motion of every electron/ion inside the gas and will call a set of handlers for every collision. The user will be able to create custom handlers and calculate the macroscopic quantities of interest (however, in section 7, we describe a command-line utility that can be used to calculate the most common swarm attributes). This characteristic allows our solution to be used in various applications, ranging from the simulation of avalanche amplifications in detectors to the determination of drift velocities in gases under strong electromagnetic fields.

An important aspect is that we do not provide any recommended cross-section data, but we leave the end-user responsible for choosing the proper set of tables. Selecting the right table set can be a difficult task in some cases, especially when there are several cross-section databases with non-negligible differences, as shown by [10]. In addition, no correction or parametric fitting is performed by our solution, providing a result based on clear and definite formulas with no arbitrary correction factors.

This paper begins with a brief description of the theoretical framework, in section 2. The library architecture is presented in sections 3 and 4, while section 6, will present and discuss a new algorithm that can be used to choose an efficient trial frequency. In section 7, we provide some brief examples of the usage of our framework. In section 8, we will test our solution against data available in literature, highlighting the performance and the limitations of our theoretical framework, which are further discussed in sections 9 and 10.

2 Theoretical framework

We will report here the theoretical framework used in our simulation tool. A charged particle, under the effect of an electric and magnetic field, will be subject to an acceleration. For non-relativist velocities, it is possible to calculate the evolution of the particle state using the classical equation of motions, as presented in our previous work [11]. For higher energies, relativistic formulas can be used, as shown by S. Chin [12], providing more accurate results with the cost of increased computation time.

The collision occurs after a random time $t_{c}$, depending on the particle speed and the gas mixture combined cross-section [13]:

[TABLE]

where $R\in\left[0,1\right)$ is a uniform random number and $\nu(t)$ is the interaction frequency of the particle in the gas. The value of the interaction frequency depends on the particle energy and the gas composition. For a gas mixture of $N$ components, where each component has an elastic interaction mode and $J$ inelastic ones, we have:

[TABLE]

where $\varepsilon$ and $m$ are, respectively, the particle energy and mass, $n_{k}$ is the molecule number density and $\sigma_{k}$ and $\sigma_{ki}$ the elastic and inelastic cross-sections.

Resolving eq. 1 for every interaction would require the numerical solution of the integral for each collision, thus requiring a substantial computing time. To avoid that, it was used the null collision technique presented by Skullerud [14]: this technique allows to bypass the numeric integration of $\nu(t)$ replacing it with a constant trial frequency $\nu^{\prime}\geqslant\nu(t)$. The eq. 1 now becomes:

[TABLE]

This substitution is possible only if we consider a fraction of the interaction as null collision, interactions that does not alter the direction and energy of the particle. To decide if an interaction should be considered null, we use a uniform random number $R\in[0,1)$ and mark as null all collisions satisfying the condition:

[TABLE]

where $\nu(t)$ is the real interaction frequency right before the next collision.

It is important to remark that determining a reasonable value of $\nu^{\prime}$ is quite important because it directly affects computing performance. A too high value will generate a high number of null collisions, decreasing the computing performance. A too low value will generate situations where $\nu(t)>\nu^{\prime}$, invalidating the result and requiring a re-computation of the step. Because several strategies are available [14, 15] to compute a proper value for the trial frequency, we decided to allow the user to choose the best strategy to determinate the $\nu^{\prime}$ value (details in section 6).

Further, we have to handle the collisions: for a correct calculation, we would need the differential cross-sections for each process, in the form of $\sigma(\varepsilon,\theta)$. Unfortunately, such cross-section tables are difficult to be obtained experimentally, and in the literature can be found such tables just for a limited set of gases. Instead, it is quite common to find integral cross-section for all gases used in gas detectors in the form of total elastic cross-section $\sigma_{el}(\varepsilon)$ or momentum transfer cross-sections $\sigma_{mt}(\varepsilon)$.

In our implementation we decided to use the approach presented by Okhrimovskyy et. al. [16] (however, we have to mention that there are other valid approaches, such as the one presented by Longo and Capitelli [17]). This approach allows using a pseudo-differential cross-section generated by the combination of both $\sigma_{el}$ and $\sigma_{mt}$. In this way, the scattering angle $\theta$ for atomic gases like $Ar$, $Ne$, $Xe$, etc. can be calculated using a formula based on the screened Coulomb potential:

[TABLE]

where $R\in[0,1)$ is a uniform random number and $\varepsilon$ is the dimension-less energy in atomic units. For non-polar molecular gases, like $CO_{2}$, $CH_{4}$ or $O_{2}$, a similar formula can be used [16]:

[TABLE]

where $\xi$ is a dimensionless value derived solving the equation:

[TABLE]

Finally, we calculate the energy exchange for a given interaction, knowing the deviation angle $\theta$. We use the model of Fraser and Mathieson [18], which gives an analytical expression for both elastic collisions:

[TABLE]

and inelastic collision:

[TABLE]

where $m$ is the incoming particle mass and $M$ is the target molecule mass, $\varepsilon_{i}$ and $\varepsilon_{f}$ is the energy before and after the collision and $\varepsilon_{k}$ is the threshold energy characteristic of the given inelastic process.

Equation 9 can be used to calculate the energy transfer for any inelastic process. However, there are two exceptions: for attachment processes, where the $\varepsilon_{k}$ can be zero, we can assume all the energy is transferred, considering them as plastic collisions. In the case of ionizations, the equations remain valid. However, another step is necessary: after the extraction of the electron, we have to calculate the energy transfer between the bullet and the extracted electron. We can model it as an elastic collision. Failing to do so, may cause the drift velocity to diverge for high fields, when the ionization processes overtake the elastic ones.

2.1 Effect of gas temperature

Equations 8 and 9 are very useful to calculate the energy transfer at each collision, both for elastic and inelastic collisions. However, analyzing these formulas, we can notice that, while we account for the mass of both bullet particle and target molecule, we do assume the target at rest. This is a very reasonable assumption in any electron/ion drift experiment: at normal conditions the mean energy of any gas molecule is in the scale of $\approx$25\text{,}\mathrm{meV}\text{/}$$. This is much lower of the mean energies of electrons/ions in a common drift experiment, which is in the order of tens electron-volts.

However, we want our solution to be suitable for a wide range of EM fields, with values ranging between $1\text{\times}{10}^{-3}\text{\,}\mathrm{Td}\text{/}$ and up to over $1\text{\times}{10}^{4}\text{\,}\mathrm{Td}\text{/}$. For values lower than $0.01\text{\,}\mathrm{Td}\text{/}$, we can not assume anymore the gas components at rest. Each molecule will have a speed distributed according to the Maxwell-Boltzmann distribution which will affect the charged particle drift in the gas.

We calculate, in [19], the effects of the temperature of the gas in the collision dynamics:

[TABLE]

where $\textbf{V}_{1}$, $\textbf{V}_{2}$ and $\textbf{V}_{cm}$ are, respectively the bullet particle, target and the center-of-momentum velocities in the laboratory frame, and $\hat{\textbf{u}}_{1}^{\prime}$ and $\hat{\textbf{u}}_{2}^{\prime}$ the new bullet and target velocities in the center-of-momentum frame.

2.2 Extension of the cross-section tables

The extended range of the EM fields that we decided to support, put stress on the cross-section tables we can use. Many tables have a limited energy range being limited, usually, to $\approx$1\text{,}\mathrm{keV}\text{/}$$. This could seem a reasonable value, but we have not to ignore the fact that, if we have a big number of bullets and/or interactions, a fraction of them will overpass this energy boundary. Fortunately, if we look at a log-log plot of the cross-section tables, we can realize that it is possible to easily extend these tables.

For elastic processes, we can extend to the left (lower energies), using a constant value based on the first values of the table. For inelastic tables, instead, we can truncate the value to zero for values below the threshold energy. The right extension, for higher energies, we found it is possible to fit a straight line, using two near points:

[TABLE]

where $K$ is a dimensionless number, $\sigma_{1}$ and $\sigma_{2}$ are the cross-sections in two points near the boundary, $\varepsilon_{1}$ and $\varepsilon_{2}$ are the respective energies and $\varepsilon$ is the requested energy.

2.3 Determination of swarm attributes

As stated in the introduction, our core library, Betaboltz, does not perform, by default, any computation of swarm attributes (drift velocity, ionization coefficients, diffusion, etc.), delegating to the final user the duty to create its own analysis using the provided handler mechanism. However, we can not ignore the importance of such attributes within the low temperature plasma community, so we provide a command-line utility, drifter, which allows their calculation without setting up a full analysis.

In this section, we will briefly present how such attributes are calculated. First, we use an infinite volume with a uniform static electric field. The axes are arranged so that the $z$-axis is aligned with the same direction of the electric field. A magnetic field may be specified at an arbitrary angle $\theta$ in the plane $xz$.

Then, $N$ primary electrons are placed at $(0,0,0)$ and they are left drifting for an $I$ number of real collisions. We then can repeat the simulation for $E$ events and, eventually, for $R$ runs (we adopt the division in events and runs as commonly used in the HEP community).

At the end of each event, we calculate the position of the individual particle along the z-axis. This information allows determining the drift velocity $W$, along the $z$-axis, using the following formula:

[TABLE]

where $z_{n}$ is the $z-$position of the $n$ particle when the particle is destroyed (because it reached the $I$ collisions count), and $t_{n}$ is the time when the destruction occurs.

For the calculation of the Townsend coefficients, $\alpha$ and $\eta$, we use the concept of dummy ionizations and attachments. We do calculate the energy exchange of such interactions, but we do not create/suppress any particle, limiting our-self to count the number of particles that would have been created or destroyed. This method is effective and allows to quickly calculate the coefficients in very high fields, where a full simulation would be hard to achieve::

[TABLE]

where $G_{n}$ is the number of new particles created by ionizations and $D_{n}$ the number of particles destroyed by attachments.

For the mean energy, we use a different approach. We perform a plain mean of the energies just before and after the collisions for all particles and events:

[TABLE]

where $n$ is the particle index, $e$ is the event index and $t_{i}$ is the time when the $i$ collision occurs.

For diffusion coefficients and magnitude, we calculate, at the event end, the mean velocity $W$ of the swarm along the $z$ axis:

[TABLE]

where $z_{n}$ is the position of the $n$ particle along the $z$ axis, and $t_{n}$ the time needed to move such distance. Then, for each particle, we calculate the displacement vector as:

[TABLE]

where $W$ is the mean velocity calculated above, $t_{n}$ the time when the particle is destroyed and $\hat{\textbf{z}}$ a dimensionless versor pointed along the $z$ direction. We can now calculate the transversal, $D_{T}$, and longitudinal, $D_{L}$, diffusion coefficients as:

[TABLE]

and their respective magnitude:

[TABLE]

3 Software architecture

In this section, we will present the architecture of our software. The project is divided into several modules:

Univec: a custom-created library which allow vector operations with compile time dimensional analysis [20].

• 2.

ZCross: a library that allows reading the cross-section tables in XML format [21, 22].

• 3.

ZCrossGPU: an optional library that allows performing the cross-section interpolations on GPU [23].

• 4.

Betaboltz: the main simulation library [24].

• 5.

Drifter: a command-line utility to calculate the swarm attributes of a gas mixture [25].

In fig. 1, we show the organization of the classes in the Betaboltz software library. Some classes are defined as abstract classes, and we encourage the end-user to implement them for their specific needs, although we provide a set of concrete implementations to cover the most common scenarios.

In this section, we will analyze each class, providing some details about their purpose and the implementations already existing by default. We want to remark that this is only a very brief introduction and more details and examples can be found in the software user guide [26].

It is worth mentioning that we heavily use the concept of static dimensional checking [27] aiming to reduce of the probability of software bugs and an improved type consistency.

3.1 Betaboltz Simple

This class is the starting point of the library and allows the simulation of an arbitrary number of electrons or ions for an indefinite amount of time. To be able to perform its job, we need to provide:

A BaseDetector, providing the geometry, fields, and gases of the detector.

• 2.

A list of BaseBulletLimiter, providing the conditions which will halt the simulation.

• 3.

A list of BaseHandler, providing the actions to perform when something notable happens in the simulation.

The only mandatory element is BaseDetector, which divides the space into volumes with specific gas conditions and EM fields. The other two elements are optional: if no BaseBulletLimiter is added, a default one will be added, which will dispose the particles once outside the detector boundaries. If no BaseHandler is provided, no handler will be executed during the simulation.

3.2 BaseBulletLimiter

This class will take care to remove from the simulation particles when a specific condition occurs. It can be used to limit the simulation up to a particular time or a specific energy range. In table 1, is presented a list of existing implementations of this class, covering the most common usage. If none of the existing limiters satisfies the user’s needs,it is possible to create a custom limiter implementing the method BaseBulletLimiter::isOver().

3.3 BaseHandler

This class contains the actions to perform when notable events take place during the simulation. In table 2, are listed the methods which are called during the simulation. Two implementations of this class already exist: PrintProgressHandler is used to print to the console the simulation progress while ExportCSVHandler is used to write into a CSV file the result of a simulation. A custom handler can be created overriding one or more of the methods listed in the table mentioned above.

3.4 BaseDetector

This class contains the detector geometry: the physical space is divided by volumes, identified by a non-negative integer. Each volume has its own GasMixture, BaseField and BaseTrialFrequency. Negative integers are reserved and used to indicate the out-of-detector state. Two implementations exist to perform a simulation without the need to define a detector: InfiniteDetector, a detector with a single infinite volume, and BoxDetector, a rectangular cuboid containing a single volume.

3.5 GasMixture

This class is used to define the gas mixtures. It is possible to define gases of any molecule, using common formula notations like Si(CH3)4: the software will compute the molecule mass and find a match in the cross-section database. Gas components can be specified by mass densities, molar densities or molecule number densities.

3.6 BaseField

This class specifies an electromagnetic field. For the moment, only uniform and static fields can be used. However, it can be easily expanded to non-uniform and non-static fields, given that an analytical solution can be implemented in the moveParticle() method. At this time, only two implementations are provided: UniformFieldClassic, which provides the classical motion equation for both E and B fields and the class UniformFieldRelativisticChin implementing the equation of motion for relativistic velocities as presented by S. A. Chin [12].

3.7 BaseScattering

This class is used to define the scattering algorithm, i.e. how the $\theta$ and $\phi$ angle should be computed for each collision. While in our solution only two implementations were provided, isotropic and Okhrimovskyy’s algorithm (described in section 2), we decide to use a modular architecture, which allows adding and testing new scattering theories via, implementing the BaseScattering::scatter() method.

3.8 BaseTrialFrequency

This class provides the strategy used to determine the trial frequency during the simulation. This parameter is quite important because it directly affects the simulation performances. Some concrete implementations are discussed in section 6.

4 Multi-threading implementation

As stated in section 3, we planned to have integrated support of multi-thread execution. Proper multi-thread implementation in C++ is not a trivial task: if not carefully implemented, multi-thread execution can lead to run-time data races and deadlocks. This kind of errors, can block the simulation or, even worse, lead to incorrect results.

We decided to base our implementation on a common and well-tested solution named OpenMP [28]. This helped us to focus only on the logical implementation of the parallel sections, letting the library manage the data synchronization between the different threads.

As shown by fig. 2, we decided to use a separate task for each charged particle. Every new particle, created by an ionization, has its separate thread and can be executed concurrently with the other tasks. The OpenMP framework will schedule each task’s execution and synchronize the shared data between the threads.

The drawback of this implementation is that only simulations containing multiple initial particles or presenting ionization processes can benefit from a higher number of execution threads.

Multi-threading support can be enabled or disabled via the ENABLE_OPENMP cmake variable (default value is ON). The number of maximum active threads can be set using the common OpenMP methods (omp_set_num_threads() function or the OMP_NUM_THREADS environment variable) or our BetaboltzBase::setNumThreads() function.

5 GPU acceleration

One of the most time-critical operations, during a simulation, is the solution of eq. 2. Solving this equation is computationally simple, requiring only basic arithmetic operations. However, if we analyze the real case usage of this equation, we realize it could become a critical bottle-neck, at least in some situations.

First, we have to sum all the components of the gas mixture (e.g. in common air, $\approx$10$$ components, to get a high precision simulation), when for each component, we have to sum all cross-section processes (between $10100$ tables) and then, for each table, we have to perform a linear interpolation. This operation must be repeated for each collision, including the null ones: we can realize how this operation can become very demanding when simulating complex gas mixtures.

However, due to the specific nature of this problem, we can take advantage of GPU computing to perform all the sum and interpolations using concurrent processing cores. Călin Banu implemented this feature as an optional external library [23], which can be compiled and linked to the main Betaboltz framework. It uses the CUDA library to offload the calculation of the real frequency on the GPU.

This implementation should be considered a proof of concept about the feasibility and the ability to provide results with the same accuracy of the non-GPU implementation. Additional work is required to improve the efficiency of the available solution and to make it competitive against the CPU only implementation.

6 Trial frequency strategies

In section 2, we highlighted the importance of choosing an optimized value for the trial frequency $\nu^{\prime}$ used by the null collision algorithm. In the original article by Skullerud [14], an algorithm is discussed on how a proper trial frequency can be determined. During our preliminary testing, we realized how vital a proper trial frequency is to get good performances. Several studies focus their attention on optimizing the execution of Monte-Carlo code using different algorithms built on lookup tables based on particle energy [15].

Let’s now introduce the concept of simulation efficiency. To do so, let’s briefly recall the null collision algorithm: during a simulation, after a time related to $1/\nu^{\prime}$ (see eq. 3) we will have a potential interaction, which can be real or null. Then, we use a uniform random number to decide is the collision is null or not, as shown by eq. 4. Also, we have to verify if, just before the collision, the criteria $\nu(\varepsilon)\leq\nu^{\prime}$ is respected: if not, we need to revert to the previous valid state and try again with an increased trial frequency. We want to remark that reverting to a previous state has a computational cost higher than a null collision, due to the need of restoring the previous state, discarding the null-collisions since the last real collision.

Observing a simulation, we can determinate three main figures related to the simulation frequency: the first one is the real frequency, the number of collisions per second altering the particle energy and direction, named $\nu_{real}$, the second one is the number of null collisions, $\nu_{null}$, and the last one, how many times the simulations had to be halted and restarted from a valid state, $\nu_{fail}$. We can define the efficiency of our simulation as:

[TABLE]

It is important to notice that the chosen trial frequency is not constant, but we can change it during the simulation to adapt to the current particle energy. However, changing it at each collision will invalidate the null-collision algorithm, so a trade-off must be accepted to ensure the correctness of the result. We implemented this trade-off in the form of a grace number of real collisions, expressed as a number of real collisions.

We query for the initial trial frequency, and we get a trial frequency and the number of grace collisions. 2. 2.

For each null or real collision, we keep constant the trial frequency, and we decrease the grace collisions counter. 3. 3.

When the grace counter reaches zero, we calculate a new trial frequency (using the onNull or onReal method according to the last collision type). 4. 4.

When a collision fails, regardless of the grace counter value, we restore the state of the last real collision, and we query a new trial frequency using the onFail method, resetting the grace counter.

In the library we developed, we provide four concrete implementations of the BaseTrialFrequency strategy: two static strategies (the collision frequency does not change during the simulation) and two adaptive ones: we found a $5000$ collisions grace period to be a reasonable value for our adaptive trial strategy.

6.1 FixedTrialFrequency

This is the simplest one, the user can choose the trial frequency and it will remain fixed for the entire simulation. If the frequency is too low at any point of the simulation, the execution will halt and an exception will be thrown. This strategy is recommended when we know the nature of the problem well, and we know the maximum frequency we can get in the energy range of interest.

6.2 LazyTrialFrequency

This is the safest algorithm available to determine the trial frequency, but unfortunately, it has quite a poor performance when applied to real cases. Given a gas mixture, it scans the local maximum and minimum of the collision frequencies of the gas components. For each local peak, the algorithm calculates $\nu(\varepsilon)$ and at the end will keep the highest value from those calculated.

Using the value calculated in this way, we can assert that $\nu_{fail}=0$, so the simulation will never revert to a previous state. The drawback of this solution is that $\nu_{null}\gg\nu_{real}$, providing, in our testing, for several gases, has an efficiency low such as $\eta\approx 2.5\%$.

6.3 BinnedTrialFrequency

Conceptually similar to LazyTrialFrequency, it substantially improves its performance by dividing the energy range in bins. As we can see in fig. 3, the real interaction frequency can change abruptly within a small energy range. Using this class, the user may control the number of bins for each decade (we define $1\text{\,}\mathrm{eV}\text{/}10\text{\,}\mathrm{eV}\text{/}$ as a decade, $10\text{\,}\mathrm{eV}\text{/}100\text{\,}\mathrm{eV}\text{/}$ as another, and so on). For each bin, we calculate the maximum possible frequency achievable. When we have a fail, e.g. when we cross from a bin to another, we retry with the maximum frequency between the two bins. If it still fails, then we throw an exception.

6.4 VariableTrialFrequency

During the development of the previous strategies, we noticed that the real interaction frequency is strictly related to the particle energy and gas conditions. This algorithm will try to balance the $\nu_{null}$ and $\nu_{fail}$, to achieve a good efficiency $\eta$. If we have $\nu^{\prime}=\nu_{real}$, we will surely get $\nu_{null}=0$ but, from the other side, we can get a $\nu_{fail}\gg\nu_{real}$, with the related performance loss. In this algorithm, we set, at the end of each grace period, a new value for $\nu^{\prime}$ with a specified overhead $\zeta$ assigned by the user (default value $25\text{\,}\mathrm{\char 37\relax}\text{/}$):

[TABLE]

Using this approach we were able to obtain an efficiency $\eta$ up to $\approx$80\text{,}\mathrm{\char 37\relax}\text{/} in several simulations when using a $\zeta\approx$25\text{\,}\mathrm{\char 37\relax}\text{/}, as shown in figs. 4 and 5. It is important to remark that this is a state-less approach: we do not keep any state variable for the moving particles to reduce memory usage. State-full approaches could improve the efficiency to higher values, but it was not the main aim of our study.

6.5 Writing custom trial frequency strategies

It is possible to write a custom trial frequency strategy, to increase the efficiency in certain situations. A custom strategy can be written by extending the class BaseTrialFrequency, implementing these methods:

getInitialTrialFrequency

• 2.

getNextGrace

• 3.

getNextTrialFrequencyOnReal

• 4.

getNextTrialFrequencyOnNull

• 5.

getNextTrialFrequencyOnFail

It is important to notice that these methods are const, forcing to write a state-less algorithm. In the future, we will evaluate the opportunity to make these methods non-const, allowing a more efficient strategy at the cost of increased memory usage. During the simulation, it is possible to check the current efficiency using the method BetaboltzSimple::getStatsEfficency().

7 Installation and Usage

In this section, we will briefly present some helpful information about the installation and the usage of the framework we developed. First, we need to install ZCross. We can get it via its public repository:

⬇

$ git clone https://gitlab.com/micrenda/zcross.git

$ mkdir zcross/build && cd zcross/build

$ cmake -DCMAKE_INSTALL_PREFIX=/opt/zcross ..

$ make && sudo make install

Then we need to download the cross-section databases of interests. It is possible to download them from the LXCat download page [34]. We recommend saving into a personal directory and setting an environment variable pointing there:

⬇

$mkdir$\sim$/.zcross_data/

$**export** ZCROSS_DATA=$\sim$/.zcross_data/

Finally, it is possible to test the ZCross installation:

⬇

$ /opt/zcross/bin/zcross list

It is now possible to install Betaboltz:

⬇

$ git clone https://gitlab.com/micrenda/betaboltz.git

$ mkdir betaboltz/build && cd betaboltz/build

$ git submodule init

$ git submodule update

$ cmake .. \

-DCMAKE_INSTALL_PREFIX=/opt/betaboltz \

-DZCross_DIR=/opt/zcross/share/zcross/cmake

$ make && sudo make install

and the Drifter package:

⬇

$ git clone https://gitlab.com/micrenda/drifter.git

$ mkdir drifter/build && cd drifter/build

$ git submodule init

$ git submodule update

$ cmake .. \

-DCMAKE_INSTALL_PREFIX=/opt/drifter \

-DBetaboltz_DIR=/opt/betaboltz/share/betaboltz/cmake

$ make && sudo make install

The last step is to make the executable accessible from command-line, adding the new paths to the environment variables:

⬇

$**export** PATH=$PATH:/opt/zcross/bin:/opt/drifter/bin

It is now possible to perform the calculation of the drift velocity using the drifter command-line utility:

⬇

$ drifter -g Ar CO2 -d ’ela:BSR/BSR-500|ine:Biagi8’ ’Phelps’ -x 93 7 --field 0.6 --disable-check-monoatomic

FIELD SYMBOL VALUE

run_id 0

events 1

particles 25

interactions 200000

gas pressure P 101.325 kPa

gas temperature T 293.15 K

gas density N 2.50348e+19 /cm3

electric field E 0.6 kV/cm

electric field (red) E/N 2.39667 Td

drift velocity W 4.56483 cm/$\mu$s

                               $\pm$   0.0464798 cm/$\mu$s

mobility $\mu$ 7608.05 cm2/V$\cdot$s

                               $\pm$     77.4664 cm2/V$\cdot$s

mobility (red) $\mu$0 7088.99 cm2/V$\cdot$s

                               $\pm$     72.1813 cm2/V$\cdot$s

mobility (norm) $\mu$N 1.90466e+23 /cm$\cdot$V$\cdot$s

                               $\pm$ 1.93935e+21 /cm$\cdot$V$\cdot$s

energy $\varepsilon$ 0.451974 eV

                               $\pm$ 4.94475e-05 eV

townsend alpha coeff $\alpha$ 0 /cm

                               $\pm$           0 /cm

townsend alpha coeff (red) $\alpha$/N 0 cm2

                               $\pm$           0 cm2

townsend eta coeff $\eta$ 0 /cm

                               $\pm$           0 /cm

townsend eta coeff (red) $\eta$/N 0 cm2

                               $\pm$           0 cm2

diffusion coeff D 12042 cm2/s

                               $\pm$     10482.5 cm2/s

diffusion coeff (norm) ND 3.01468e+23 /cm$\cdot$s

                               $\pm$ 2.62428e+23 /cm$\cdot$s

diffusion coeff (ratio) D/$\mu$ 1.58279 V

                               $\pm$     1.37792 V

diffusion coeff long DL 854.425 cm2/s

                               $\pm$     1070.07 cm2/s

diffusion coeff long (norm) NDL 2.13903e+22 /cm$\cdot$s

                               $\pm$  2.6789e+22 /cm$\cdot$s

diffusion coeff long (ratio) DL/$\mu$ 0.112305 V

                               $\pm$    0.140655 V

diffusion coeff tran DT 2583.28 cm2/s

                               $\pm$     2589.45 cm2/s

diffusion coeff tran (norm) NDT 6.46718e+22 /cm$\cdot$s

                               $\pm$ 6.48263e+22 /cm$\cdot$s

diffusion coeff tran (ratio) DT/$\mu$ 0.339546 V

                               $\pm$    0.340374 V

diffusion magnitude $\sigma$ 0.179933 cm

                               $\pm$   0.0762089 cm

diffusion magnitude long $\sigma$L 0.0609852 cm

                               $\pm$   0.0410937 cm

diffusion magnitude tran $\sigma$T 0.162378 cm

                               $\pm$   0.0800588 cm

real collisions 5e+06

null collisions 1.34285e+07

fail collisions 11958

efficiency 27.1142 %

elapsed time 264 s

We want to highlight the syntax we specify the cross-section tables to use for this simulation. In the previous example, we can see we use the Phelps database [32] as the $CO_{2}$ cross-section data. However, for the $Ar$ molecule, we specify both the BSR [33] and the Biagi [30] databases, separated by the vertical line. In this case, we combine the non momentum-transfer elastic cross-section tables from BSR database with the complete inelastic tables from the Biagi database, which provides exclusively momentum-transfer tables.

Finally, in section 7, we can find the list of commands used to generate the plots in this article, showing how it is possible to perform several complex analyses using basic shell scripting techniques.