A dynamic method for charging-up calculations: the case of GEM

The simulation of Micro Pattern Gaseous Detectors (MPGDs) signal response is an important and powerful tool for the design and optimization of such detectors. However, several attempts to simulate exactly the effective charge gain have not been completely successful. Namely, the gain stability over time has not been fully understood. Charging-up of the insulator surfaces have been pointed as one of the responsible for the difference between experimental and Monte Carlo results. This work describes two iterative methods to simulate the charging-up in one MPGD device, the Gas Electron Multiplier (GEM). The first method uses a constant step for avalanches time evolution, very detailed, but slower to compute. The second method uses a dynamic step that improves the computing time. Good agreement between both methods was reached. Despite of comparison with experimental results shows that charging-up plays an important role in detectors operation, should not be the only responsible for the difference between simulated and measured effective gain, but explains the time evolution in the effective gain.


Introduction
A considerable amount of work has been done over the last few years to improve the simulations of MPGDs. A full understanding of the micro-physics associated with MPGDs operation is vital for the improvement of such detectors, whose some behaviours are still not completely understood. [1][2][3].
Some works report that there is a transient period during which the effective gain changes, after voltages are applied and the detector irradiated [4,5]. The gain ends up stabilizing after minutes or even hours, depending on the MPGD and the rates of irradiation.
MPGDs were developed to detect radiation, and their main applications are for high energy physics, astrophysics, rare-event searches and medical imaging [6,7]. The Gas Electron Multiplier (GEM) [8] has been largely used in many of those applications.
The device consists of a thin polyimide (insulator) foil typically 50 µm thick. The foil is covered on both sides with 5 µm thick layers of a conductor and etched with an hexagonal pattern of holes. The device operates inside a gas medium and suitable electric potentials are applied between the upper and the lower electrodes of the structure. In this way, a very high electric field is created inside the holes. Electrons created in the drift region by interaction of external radiation, travel towards the micro-structure, being focused into the holes and accelerated. They acquire enough energy to ionize atoms/molecules of the gas, creating new ionizations. The secondary electrons undergo the same process while inside the hole and an avalanche ends up being produced, the Townsend avalanche.
The charges produced during multiplication have two possible destinations: they may be collected by conductor electrodes, both those of the GEM itself or of any other readout setup; and a fraction of them ends up accumulating in the insulator surfaces. The number of simulated avalanches is correlated with the number of primary electrons that undergo to the hole, since we assume that no charges will drift in the insulator.
Effective gain is defined as the number of secondary electrons, for each primary electron, that are collected in an electrode plane located below the GEM (figure 1b).
The electronic affinity of the polyimide usually used in these devices is high, 1.4 eV [9]. Once electrons are trapped, it is unlikely that they are able to leave the surface. The same process happens with ions, according to [10].
The effective gain of the detector strongly depends on the intensity of the electric field produced in the multiplication region. Charges accumulated in the insulator surfaces change locally the electric field, changing the amplification gain. This is known as the charging-up effect in the insulator.
Deposited charges can flow through the insulator surface and insulator bulk under the action of the electric field. Previous studies propose that the positive ions are not captured in insulators surfaces, instead they transfer their charge to intrinsic carriers of the insulator, and the conduction should be made by electrons and holes [10]. The time for charges evacuation is of the order of several hours to days, so we did not include this effect in our method, i.e. all deposited charges will remain in the same surface during the whole simulation time. This approach should remain valid if the charging-up process is much faster than the draining of charges.
Simulations of the charging-up influence in the GEM transparency were already reported [11]. In order to study the contribution of the charging-up for the effective gain variations, two methods to simulate the charge accumulation in the detector are presented. We also compare our results with available experimental data.

Calculations Geometry
An hexagonal hole pattern with a distance between holes of 140 µm was considered, the insulator thickness is 50 µm and the metal electrodes are 5 µm thick. Figure 1a shows the GEM geometry. The corresponding electric field configuration is depicted in figure 1b. The hole has a bi-conical shape, the outer diameter is 70 µm, while the narrower part has 50 µm of diameter.

Gases
In order to compare simulation results with our measurements, we simulated a gas mixture of Ar 70% / CO 2 30%. This is a penning-mixture, due to the presence of the quencher molecule CO 2 , that opens new de-excitation ways for the previously excited noble gas molecules. If this excess of energy is above the ionization threshold of the quencher molecule, an ionization may occur, whit a probability called Penning probability. We used a Penning probability of 0.7 in this simulations based on previous calculations for Ar 70% / CO 2 30% mixture [2,[12][13][14].
The drift and induction fields (electric fields applied up and above the GEM, respectively) were 0.2 and 0.3 kVcm −1 . Taking into account that the computational time strongly depends on the gain, because a higher number of electrons need to be tracked within the avalanche, a potential of 400 V between electrodes was used for the firsts simulations tests, corresponding to gains of ∼10 2 . All simulations were performed considering a temperature of 293 K and a pressure of 760 Torr.

Software platforms
The Monte Carlo calculations involved three programs. Due to the complex shape of the GEM structure, an analytic solution for the electric field in the interest region is not possible to obtain. To overcome this problem, the electric field is computed with Finit Element Methods (FEM) software, that is used to calculate the electric potential in discrete nodes of a mesh, using boundary conditions. We used ANSYS 1 to produce potential maps, to which we call generally as field maps, selecting the curved tetrahedral elements as our mesh elements, because they easily fit in sharp curved surfaces present in GEMs.
To simulate the drift and transport properties of electrons and ions in the MPGD gas medium, we used Garfield++ [15]. As input, this software requires the electric field configuration in the MPGD, the gas mixture, temperature, pressure and initial conditions of the primary charges (position, velocity and energy).
In what regards the electric field configuration, we used Garfield++ to read the potential maps calculated with ANSYS and to calculate the electric field in any point of the space by interpolation between nodes.
A microscopic approach is used to simulate the drift of the charges. This uses Monte-Carlo methods to calculate the probability to occur each type of collision during the drift (elastic, excitation or ionization). The cross sections associated with each collision type are obtained from Magboltz [16,17].
Primary electrons starts with assigned r start = (x s , y s , z s ), velocity v start = (v x,s , v y,s , v z,s ) and kinetic energy E start , drifting through the gas and producing secondary charges as it passes the multiplication region. The final position of each secondary charge, r end and the effective gain are the observables of interested that are recorded for further analysis of the charging-up effect.

Initial attempts
To start our simulations, we randomly distributed 10 4 primary electrons in the surface of a plane parallel to the GEM, located 100 µm above the GEM, indicated as the start plane in figure 1a.
In order to determine the number of collected and deposited electrons and ions, the final position of each electron and ion from avalanches are analysed: • Electrons are collected if the final z is -100 µm below the GEM (the collection plane represented in figure 1b. • Ions are collected if the final z coordinate is in the top (positive) electrode of the GEM.
• Electrons and ions are deposited in the insulator surface if after the drift the z coordinate is between -25 µm to +25 µm.
GEMs previous to charging-up, i.e. without deposited charges, are called uncharged GEMs, being the charged GEMs those who have already deposited charges due to charging-up. The deposition distributions of charges (electrons and ions separately) in the insulator, for the case previous to charging-up (figure 2a) shows that the charges are not deposited uniformly on the hole surface. In addition, the number of deposited electrons is higher than the number of deposited ions. The reason is related with the mass of each particle. Ions are heavier than electrons, and tend to follow very well the field lines, in the direction of the electrodes. Electrons has a much more chaotic movement, due to the lower mass, having more probability of ending in the insulator surfaces. This originates variations on the local electric field. After some avalanches, the distribution of new electrons and ions that reaches the insulator tend to compensate each other, due to Coulomb attraction between previous and future deposited charges (figure 2b). The local variation in the electric field will therefore vanish and a stable configuration will be achieved.
In order to simulate the effective gain variation as avalanches happen, we needed to iteratively include this charge deposition in the potential maps computed with ANSYS . The software does not provide the option to put single charges in their exact deposition position in the insulator surface. In addiction, this scenario would lead to discontinuities and numerical issues. Instead, we created small slice surfaces in the insulator foil and add the correspondent density charge to each surface. Due to the shape of the deposition, and to computational limitations of field maps files for very small finite elements, we used 24 different slices in the insulator, achieving in this way a good balance between the detail of the calculations and the needed computing power. The slices are not regularly distributed, as shown in figure 3, trying to match the z profile of the charge deposition

Constant Step Method
The flow-chart of the first iterative algorithm used to simulate charging-up iterations is in figure 4.  At the first iteration, we compute the electric field map assuming no charges in the insulator surface. Then, we import that field map into Garfield++, simulate 10 4 primary avalanches and determine the density charge deposited in each insulator slice surface. A new field map is created, with the contribution of previously deposited charges. The density charge in each slice is calculated taking into account the contribution of both the ions and the electrons ending up in the insulator surface. A new set of 10 4 primary avalanches is simulated and the process is repeated iteratively.
It was found that statistical fluctuations in the calculated gain depends on the number of simulated avalanches per step, but the number of deposited charges per avalanche seems to be less sensitive to fluctuations. A small step-size of 10 4 primary electrons was chosen in order to obtain good detail in the time evolution of charging-up. However, this small step implies hundreds of iterations until stabilization, which lead to a very heavy computation.
Since the number of deposited charges is the responsible for the local variation in the electric field, we use that observable as our control function for the iterative simulation, i.e. we stop our iterations when its value stabilizes over iterations (corresponding also to a gain stabilization).

Dynamic Step Method
In order to accelerate the simulation process, we developed an extended method that uses a dynamic step-size in each iteration. This step-size is smaller when the number of deposited charges per avalanche changes quickly, and is larger when this quantity is more constant, i.e. the deposition stabilizes.
To constrain the size of the step, we defined that the maximum total charge (sum of signs of ions and electrons) that can be added to the new field map should not be larger than 2×10 4 q e , (where q e is the elemental charge 1.6×10 −19 C). This way, the maximum allowed number of avalanches per step is equal to 2×10 4 q e G tot , where G tot is the absolute gain in each iteration. The output of this calculation give us an upper limit for the step-size, considering a maximum charge that can be added to new potential maps in each iteration. Our attempts show that this upper limit value an acceptable value in terms of the convergence and speed of the method, but other limits can be defined. The dynamic method is briefly described in the flow-chart in figure 5.
The method starts with an uncharged ANSYS field map of the GEM. In each iteration we simulate 10 3 primary avalanches, which is a good compromise between statistical fluctuations and computational time.
The number of deposited charges per avalanche, in each slice of the insulator surface, is multiplied by the variable step. For the firsts iterations, steps between 0.5 × 10 3 and 10 3 primary avalanches were used.
After the first 5 iterations (sufficient number of iterations that allows a reasonable fit), we fit the number of deposited charges per avalanche to a first order polynomial, and calculate for a given step, what the value of that function should be for the new iteration.
We then simulate iteration number 6, and compare with the predicted value from the fit: • If the difference between simulated and fitted value is larger than the maximum defined step, discard the iteration, the step is reduced to its half, and repeat the iteration.
• If the difference between simulated and fitted value is smaller than the maximum defined step, the iteration is saved, we increase the step to the double. A new iteration is calculated and new fit considers only the last 5 valid iterations

Comparison between methods
The sum of all electric charges (the integral of the deposition histograms in figures 2a and 2b) deposited in the insulator surface, per primary avalanche, is shown in figure 6a for both methods, as function of the charge produced by each avalanche, per hole (is simply the number of primary simulated electrons in each hole multiplied by the total gain).  The agreement between both methods is clear. However, the dynamic-step method saves computational resources, using about one tenth of iterations. Figure 6b represents the total gain evolution for the two methods. We observe an increase in effective gain, followed by a stabilization plateau, reached in both methods. Due to previous results, from now on we will only consider the dynamic-step method for calculations.

Charging-up effect in the GEM transmission
Primary electrons produced by incident radiation and drifting towards the GEM holes can be collected in the top electrodes, ending up not producing avalanches. The ratio between the number of primary electrons that enter the holes, producing avalanches and the total number of primary electrons simulated is defined as the electron transmission, shown in figure 7a, for several voltages applied to the GEM electrodes. The contribution of the charging-up effect in the electron transmission is more important when low voltages (<400 V) are used and negligible when higher electrical potentials are used.

Effective gain with and without charging-up
The dependence of effective gain on the voltage applied between electrodes in the GEM detector, is shown in figure 7b. The gain, after charging-up stabilization, is 10-15% higher than the situation without charging-up.

Electric field intensity variation
A 2D representation of the electric field in the GEM is shown in figure 8. Each plot is obtained by calculation of the intensity of the electric field along a plane corresponding to a vertical cross section of the GEM hole, at four different stages of the charging-up process.
We can observe that the biggest change in the electric field occurs near the electrodes. While the intensity of the electric field near the top (negative polarized) electrode decreases, it increases near the center of the hole and the bottom (positive polarized) electrode. The development of an avalanche inside the hole follows a nearly exponential model. The bigger fraction of secondary electrons is produced at the exit of the hole, in the last stages of the avalanches. There, the electric field is higher due to the charging-up effect, and thus, the effective gain increases as a result of this process.

Comparison with experimental results
Experimental measurements were performed at CERN. Physical parameters of the GEM, and the gas mixture used for measurements correspond to the simulation settings. X-ray photons were used as ionizing radiation. K α and K β photons corresponding to energies 8.0 keV and 8.9 keV respectively were emitted by the X-ray tube which employed a copper target.
Collimators to control the photon flux were used in order to regulate the rate of charging up. The gain was measured over the time, for a constant irradiation flux. The GEM structure was housed inside an air-tight chamber, shown in figure 9, that had a constant circulation of mixture, with a gas flow rate of 6 l.h −1 . Chamber pressure was maintained at 760 Torr. Figure 10 shows the schematics setup for gain calibration (left) and gain measurement (right). In order to compare with simulations, the experimental results were normalized from the time scale to charges per hole. This is done multiplying the number of primary electrons produced per second by the absolute gain, divided by the number of irradiated holes of the GEM.
For a detailed description of the measurement procedure refer to [18]. A comparison between Monte Carlo calculations and measurements is shown in figure 11a. The total gain increases as the GEM starts to be irradiated, in both cases, achieving then a plateau. A normalization to the plateau of both data was done to directly compare the time evolution of Monte Carlo simulations and experimental values, as shown in figure 11b. From this observation, Monte Carlo simulations reproduce the time evolution of the gain.
On the other hand, the value of the total gain still do not match. This can be related with the mobility of the charges in the insulator surfaces and bulk, numerical issues associated with the finite elements method, computed with ANSYS , impurities in gas and imperfections in GEMs dimensions introduced during the production.

Conclusion and future work
In this work we have presented two iterative methods for the simulation of the insulator surface charging-up in GEMs, allowing a better understanding of their response. Both methods agree between each other. However, the dynamic-step method saves computational resources. The Monte Carlo functional time behaviour of the gain as the GEM is irradiated reproduces that observed experimentally. However, the absolute scale still not agree.
Primary electrons transmission should be affected by the charging-up at lower voltages between electrodes of the GEM, but for higher voltages, used in regular applications, it does not play an important role.
Future work will include the application of the presented charging-up simulation methods to other MPGDs (e.g. THGEM) and the study of new geometries and detectors that could take advantages of this effect or minimize it.
The simulation of the mobility of deposited charges in the insulator surfaces could contribute to obtain more precise values. Refine the calculations of the electric field calculations (with ANSYS or other method) can also be important to get agreement between absolute simulated and measured   Figure 11. a) Absolute gain comparison between measuremments (red) and simulated (green) results. Same situation as figure 6b but with V GEM =380 V. b) Same plot as figure 11a but with gain normalized. Experimental data taken by Mythra Varun Nemallapudi at RD51 facilities, CERN.