Simulation of gain stability of THGEM gas-avalanche particle detectors

Charging-up processes affecting gain stability in Thick Gas Electron Multipliers (THGEM) were studied with a dedicated simulation toolkit. Integrated with Garfield++, it provides an effective platform for systematic phenomenological studies of charging-up processes in MPGD detectors. We describe the simulation tool and the fine-tuning of the step-size required for the algorithm convergence, in relation to physical parameters. Simulation results of gain stability over time in THGEM detectors are presented, exploring the role of electrode-thickness and applied voltage on its evolution. The results show that the total amount of irradiated charge through electrode's hole needed for reaching gain stabilization is in the range of tens to hundreds of pC, depending on the detector geometry and operational voltage. These results are in agreement with experimental observations presented previously.


Introduction
The time-evolution of the avalanche gain in Micro-Patterned Gaseous Detectors (MPGDs) has been a topic of major interest in the last years. Gain variations with time have been reported and are common in detectors containing dielectric materials in contact with the gas medium where multiplication occurs; the effects have been noticed particularly in detectors based on hole-type MPGDs with insulating surfaces in contact with active gases, e.g. Gas Electron Multipliers (GEMs) [1][2][3] and Thick-GEMs (THGEMs) [4][5][6][7][8][9][10].
Gain variation over time in GEM and THGEM-based detectors can arise from different processes, e.g. polarization of the insulator's volume under applied voltage, and charging-up of the insulator's surfaces exposed to the free drifting charges (electrons and ions) in the gas volume [4]. It has been also observed that detector's environmental conditions, such as gas impurities and moisture, can affect gain transients [9].
A simulation program has been developed for phenomenological studies of charging-up effects in GEM detectors [11,12]. The tool is based on the superposition principle according to which the field inside the holes is the sum of the field prior to any charge accumulating on the insulating surface and the field due to the charge accumulated on the insulators. This simulation tool predicted the time dependence of the gain in the GEM detectors with double-conical holes.
The superposition principle is also used to discuss experimental observation of gain variation over time [9,13]. In this work, the simulation tools were further improved and extended to THGEM-based detectors. We have focused on the development of a new toolkit that allows for calculations of charging-up effects of the detector's insulator surfaces in a more efficient way. The polarization of the insulator volume and its influence on gain evolution were not studied in this document. In the current simulations, the detector was assumed to be biased prior to irradiation.
The present toolkit allows for the simulation of charging-up effects, based on the superposition principle of electric fields. This toolkit, which can be included to Garfield++ [14], can be extended to virtually all MPGDs in a straightforward way. It has been applied in this work to the study of gain evolution in THGEM detectors, under typical experimental conditions, and can be downloaded from [15].

The simulation toolkit 2.1 Electric field calculations
To simulate charging-up effect in THGEMs, the method for the charging-up calculations described in [12] was chosen.
The total field E tot can be estimated using the superposition principle Here E uncharged is the electric field resulting from the voltage applied to the THGEM electrodes prior to any charging-up, E charges is the field due to the charges accumulating in the insulator surfaces (electrons and ions). This definition was also at the base of the analysis presented in [9]. While E uncharged remains constant, E charges varies due to the continuous accumulation of charges. Hence, the total field at equilibrium can be obtained in an iterative process, where in each iteration E charges is modified according to the cumulative surface charge. During the first iteration, that corresponds to the simulation of the multiplication process in a THGEM before irradiation, i.e., prior to charge accumulation on the insulator surfaces, a set of primary avalanches is simulated. A fraction of the free charges (drifting electrons and ions) end up accumulating on the insulator, as depicted in Fig. 1a. These charges will modify the electric field by some amount, resulting in a new field map, E tot , that can be obtained from the superposition of E uncharged and E charges , as depicted in Fig. 1b.

Simulations steps
Previous simulation works related to charging-up effects in MPGDs relied on calculations of the electric field in the region of interest and simulation of the avalanche and detector's gain. The amount of charges produced on the avalanche that ended up trapped on the insulator surfaces was evaluated, a new electric-field map resulting from the addition of these charges was calculated, and the process was repeated, as shown in Fig. 2. Using this iterative method, the simulation results of the charging-up process were in agreement with experimental data obtained with GEM detectors [12].
The finite element method (FEM), as implemented in ANSYS® 1 , is used to calculate the electric potential in the gas volume, which then is used to infer the electric field at a given point of all the avalanche electrons' path [16]. Avalanche-size calculations and transport properties of (a) (b) Figure 1: (a) Schematic drawing of the charging-up process during an avalanche. (b) Charging-up superposition algorithm of the THGEM: the total electric field is a sum of E uncharged , calculated when voltages are applied on the electrodes and the detector was not exposed to ionizing radiation yet, superimposed with the electric field induced by accumulated charges on the insulator's surface of the detector E charges . the charges (electrons and ions) are calculated in Garfield++ [14] (which interfaces with Magboltz [17,18]), that simulates collisions of electrons in a given gas mixture and detector geometry. The previous implementation of the simulation toolkit was CPU and time consuming. It was based on an iterative transition between the FEM software and Garfield++ after each iteration, as depicted in Fig. 2. The update of the field-map, according to Garfield++ charge deposition results, required a new FEM calculation -to include the original electric field produced by the voltages applied in the detector as well as that arising from the accumulation of charges on the insulator surfaces. The simulations could take as long as several hundreds days of CPU time for a given gas mixture and amplification voltage. Even the use of parallelism tools (multiple processors, multithread, etc) can only reduce the computing time but not the algorithm efficiency, which sometimes even required manual intervention in the files between iterations.
Hence, in this work we have elaborated a faster solution based in a different algorithm. While in the previous method, the new filed-maps between each iterations were calculated in the FEM software, forcing a transition between the FEM software and Garfield++ during after each iteration, in the method now described, several field-maps, previously calculated in the FEM software, are imported to Garfield++, and further iterative calculations are performed without running new instances of the FEM software, speeding up the simulations. A schematic description is given in Fig. 3.
A new E charges is calculated at the end of each iteration and is used to updated E tot . The next iterations takes into consideration the updated E tot of the previous iterations. The different electric-field during iterations changes the avalanche process, which influences on the measured gain. Assuming a certain irradiation rate, a time-evolution of the gain can therefore be simulated. A flowchart of the superposition algorithm is depicted in Fig. 3. For a THGEM, the following procedure is applied for computing the field maps of interest: • The THGEM insulator surface is divided into 20 slices as shown in Fig. 4.
• The E uncharged is calculated from the voltages applied; • For each slice, a field map, considering only the electric field due to an unitary charge distributed uniformly over the exposed surface of the slice, is calculated -defined here as E slices .  Once these field maps are available, the Garfield++ simulation can start. To perform the simulation, Garfield++ imports the coordinates of the nodes and the associated electric potential as a list. In order to calculate E charges , an avalanche is calculated, and the position of all the electrons or ions deposited on an insulating surface is stored. The total number of accumulated charges N = N electrons − N ions (the number of accumulated electrons minus that of accumulated ions) is assigned to the corresponding slice. Since the expected variation of the field from N charges on a given slice for a single avalanche is very small, for fast algorithm convergence, the total number of accumulated charges attached to a single slice is usually multiplied by a constant value s (step-size). Therefore, in the case where N charges have stopped on the surface of a slice j (without other extra charges accumulated in the remaining surfaces), we need to calculate, for each node i, the electric potential V for that node which is given by Eq. 2.2: where V ( j, i) is the electric potential on node i due to the presence of a unitary positive charge in the surface of slice j, calculated from the E slices described previously. Depending on the accumulated charges, N can either be positive (more ions) or negative (more electrons).
Once E charges is calculated, the next iteration is launched, where the total electric field is given as a sum of E charges and E uncharged (Eq. 2.1). In the following iteration, a new E charges is calculated according to a new amount of accumulated charges, and this procedure is repeated until the total amount of accumulated charges N after a given number of iterations is equal to zero -E tot became constant with the increasing number of iterations.
3 Application of the method to THGEMs

Simulation setup
With the method described in the previous section, we are able to simulate the charging-up effect using only Garfield++, even though more than one field map has to be introduced as input at the beginning of the simulation.
The THGEM-detector geometry simulated in this work has an electrode with an hexagonal pattern of hole-pitch a = 1 mm, hole diameter d = 0.5 mm, with no rim, and FR4 substrate having a thickness t = 0.4 mm, 0.8 mm and 1.2 mm. Electric fields of 0.2 and 0.5 kV cm −1 were set across the drift and induction gaps respectively.
Voltages ∆V THGEM were applied across the THGEM electrode. Several gas mixtures were used in the simulations. In most of the cases Ne and Ne/CH 4 (5%), however to obtain similar electron multiplication at higher voltages, the simulation work was extended to Ar-based gas mixtures: Ar, Ar/CH 4 (5%) and Ar/CO 2 (5,7%). Since some of these are believed to be Penning-mixtures [19], Penning factors of 0.40 for Ne/CH 4 (5%), 0.18 for Ar/CH 4 (5%), 0.35 for Ar/CO 2 (5%) and 0.4 for Ar/CO 2 (7%) [20][21][22] were used. Standard room-temperature conditions were considered (T=293 K and P=1 bar). A constant irradiation rate of 10 Hz mm −2 , 8 keV x-rays, was assumed. In our simulations, a random location for primary charges originated by the incoming x-rays is assumed and the charges are drifted towards the THGEM holes from different positions.
At each iteration, the gain for a single-electron avalanche multiplication was calculated. Electron multiplication is a stochastic process, and the gain is determined as the average number of electrons generated during the electron avalanche process. At each iteration a given number of avalanches between 100 and 1000 was simulated (denoted in the text as n AV ). In all simulations, the statistical error on the mean gain of the detector was taken to be In each avalanche, the position of all the electrons or ions deposited on an insulating surface is stored, and assigned to the corresponding slice. Electrons that do not end up on the insulator are either collected on the THGEM bottom electrode, on the induction plane, while non-trapped ions end up on the top THGEM electrode or on the drift plane. Some of the electrons are attached to the gas molecules (CO 2 or CH 4 ) and are not taken into account in the gain calculation.

Physical interpretation
The number of simulated iterations (n iter ) can be translated to a physical irradiation time (t) as follows: given the number of primary charges induced in the gas by the ionizing event 2 (n p ), the irradiation rate (R) and the step-size (s) can be related to a physical irradiation time: An example of a simulation result of the gain evolution in a THGEM detector, in Ne/CH 4 (5%) gas mixture with the parameters: t=0.8 mm, a=1 mm, d=0.5 mm, and no-rim, is shown in Fig. 5. Gain evolution curves over time were fitted using: where G(t) is the gain at a given instant, G 0 is the initial gain, τ is the relaxation time constant (or characteristic time), and ∆G is the gain variation. The final gain, G F , can be calculated in terms of G 0 and ∆G: The quality of the fit of Eq. 3.3 indicates upon a good agreement (χ 2 /n.d.f. = 0.71) to the simulation (n.d.f., number of degrees of freedom, being the number of iterations). The exponential model is driven by the assumption that the amount of the accumulated charges is proportionally decreasing with the amount of charges attached to the insulating surfaces. In Fig. 5, after about 100 min, the gain reaches a plateau, indicating that the ongoing charge accumulation is no longer capable of modifying the electric field in a way that could affect the gain because the electric fields due to new incoming electrons and ions cancel each other.
Using the formalism above, one can define a characteristic charge, Q tot -the total number of charges, which pass through the THGEM holes during the relaxation period (τ): Replacing the characteristic time by characteristic number of iteration (τ iter ), using Eq. 3.3, one can rewrite Eq. 3.5 in terms of the step-size and τ iter : The choice of the step-size s is an important parameter for simulation convergence. An example of simulation results of THGEM detector, in Ne/CH 4 (5%) gas mixture with the parameters: t=0.4 mm, a=1 mm, d=0.5 mm, and no-rim, with different step-size values (from 0.5 × 10 6 to 40 × 10 6 ) is shown in Fig. 6. The behaviour of the gain G(t) should remain unchanged, regardless the step-size. In Fig. 6a and 6b the fitted τ is similar, while in Fig. 6b and 6c clearly show that higher step − size lead to different trends.

Convergence rules
Since the step-size should not affect the gain over time, a plot of the gain stabilization over time, simulated for different step-size values from 0.1 × 10 6 to 0.5 × 10 6 , for the same detector geometry, is shown in the Fig. 7a, only for those step-size values that do not vary the gain trend. According to the results obtained in Fig. 6a, for a step size of 0.5 × 10 6 the gain stabilizes after about 3-4 iterations, therefore the step-size was reduces by a factor of 5 to inspect the gain stabilization over 50 simulated iterations. From these curves, the average fitted τ is (13 ± 2) min. The plot of 1/τ iter as a function of s × (G 0 − ∆G e ), depicted in the Fig. 7b, allows to calculate the value of Q tot = 3.57 × 10 8 charges (≈ 57 pC).

Voltage effect on gain stabilization
Simulations were performed at different THGEM operation voltages, in the range ∆V THGEM = 500 -1100 V, aiming at investigating their influence on gain evolution. To compare the gain evolution for a broader range of ∆V THGEM values, the simulation was performed for different gas mixtures (with the highest operation voltages in Ar mixtures).
In Fig. 8, the calculated Q tot value (obtained from Eq. 3.6) increases with voltage, indicating that the higher the voltage, the higher the irradiation rate (or time) is needed to reach gain stabilization.

Electrode-thickness effect on the gain stabilization
Gain stabilization was simulated for THGEM electrodes of different thickness, in Ne/CH 4 (5%) and Ar/CO 2 (5%) (Ar-based mixtures have about two-fold higher operational voltages, and lower transverse electron diffusion for electric fields below 10 kV cm −1 [17]). The Q tot value calculated from the charging-up of THGEM hole walls and the fraction of the avalanche charges attached to the THGEM walls for an uncharged electrode are shown in Table 1. An increase of the Q tot is observed with the decrease of the electrode thickness; it is found to be more pronounced than the dependence on the THGEM voltages or electron transport parameters, e.g. for thicker structures, and higher applied voltages, a lower total-charge value would be needed for modifying the gain. For thicker electrodes, the fraction of accumulated charges on the THGEM holes is larger, due to a broader transverse expansion of the avalanche.

Conclusions
The present study introduced a new toolkit for charging-up calculations in MPGDs that can be included in a future version of Garfield++; it is expected to ease simulations of phenomena occurring in gas-avalanche detectors. A case-study is shown, where the new toolkit was used to model the gain stabilization due to charging-up effects in THGEM detectors; a quantitative study was conducted to characterize various aspects of these effects. The gain stabilization typically occurs in a few minutes to a few tens of minutes, depending on several parameters of the THGEM-electrodes. Higher voltages applied across the THGEM-electrode require a larger amount of charge to modify the electric field within the detector holes, resulting in longer gain-stabilization characteristic times. On the other hand, in thicker electrodes, the avalanche is broader due to the expansion of the avalanche-electrons within the holes, and the charging-up of the electrode's insulating surface occurs faster. Most of the simulation studies were carried at relatively low operation voltages (low gains). Differently to GEM detectors with single-conical or double-conical geometry, where the gain stabilization exhibits an increase over time, in THGEM electrodes with cylindrical holes, the effect is a decrease of the detector gain due to a different charge distribution along the hole's surface. A typical characteristic charging-up time in THGEM electrodes is found to be equivalent to a total irradiated charge of hundreds of pC similarly to the values reported for the GEM electrodes with double-conical geometry. The simulation results of the initial decrease of the gain are found to be in a good agreement with previously reported experimental data [5,6]. Further experimental studies are in course to validate our charging-up model.