Numerical modeling of air-vented parallel plate ionization chambers for ultra-high dose rate applications

Purpose: Air-vented ionization chambers have been the secondary standard for radiation dosimetry since the origins of radiation metrology. However, the feasibility of their use in ultra-high dose rate pulsed beams has been a matter of discussion, as large losses are caused by ion recombinations and no suitable theoretical model is available


Introduction
FLASH radiotherapy has arisen as a new form of therapy based on the seminal work of researchers including Fauvaudon et al. [1] who have examined the radiobiological advantages of ultra-short radiation delivery times. Existing scientific evidence only confirms that research is needed to clarify the limits, benefits and potential risks of this novel form of therapy [2][3][4][5][6]. However, in the context of this new ultrahigh dose rate application [7], many established ionization chambers currently used as a secondary standard for absorbed dose to water exhibit severe saturation problems [8]. Thus, in order to define the parameters in which the FLASH effect is observed and to investigate its clinical feasibility, we require reliable, accurate and convenient dosimetry standards.
The European UHDpulse project has the fundamental objective of developing and improving dosimetry standards for FLASH radiotherapy, very high energy electron (VHEE) radiotherapy and laser-driven medical accelerators. Such a challenge cannot be addressed in a unidirectional way. In fact, there is a collection of possible paths to perform metrology in the ultra-high dose rate range. These include dosimetry like Fricke, alanine, calorimetry, ionization chambers and solid state detectors [9][10][11].
In this paper, we will focus on the use of ionization chambers in the ultra-high dose rate (UHDR) regime for pulsed beams. At low doses per pulse, the performance and, more specifically, the charge collection efficiency (CCE) of plane-parallel chambers is well known. The incomplete charge collection is mainly due to the so-called volume recombination, which occurs when two ions with opposite charges neutralize themselves during the drift of the charge. Other contributions to the incomplete charge collection such as initial recombination and back-diffusion recombination are independent of the dose per pulse. The ion recombination correction factor is usually determined using the two-voltage method based on Boag's formalism [12]. However, in the Intra-Operative Radiotherapy (IORT) and UHDR regimes, the corrections derived from the theories developed by Boag fail. When an ionization chamber is irradiated with a UHDR beam that has a large dose per pulse, the so-called space charge effect becomes more relevant [13]. Thus, a theory capable of reproducing the charge transport and collection in an ionization chamber must take into account this effects not considered in Boag's theory in order to be accurate.
Below, we present a computational model capable of handling these effects in the UHDR regime. By solving the transport equations, our aim is to obtain a more complete picture of the behavior of a parallel plate ionization chamber (PPIC). As we will describe in the following, accurate knowledge of transport parameters such as drift velocities, attachment to neutral molecules and volume recombination parameter is necessary for a precise description of the behavior of the ionization chambers. The results of the model presented are satisfying, bringing us closer to a deeper understanding of the processes taking place in PPICs at UHDR and to an accurate description of ion recombination effects in air ionization chambers under UHDR and conventional dose rates.

Mathematical description of charge carrier transport
The charge carrier transport inside an air ionization chamber can be described via a system of coupled partial differential equations. In these equations, it is possible to allocate the different fundamental processes that take place in an ionization chamber during irradiation. The set of equations can be made more complex according to the level of detail that we want to assign to this description, since the actual number of species and processes is quite large [14]. Nevertheless, it can be assumed that a summarized approach of the physics could be sufficient for most of the dosimetry properties of interest. The following simplifications have been adopted for the building of the model: • The electric field is considered homogeneous in the ionization chamber and no edge effects are considered. • The ionization chamber is symmetric along the transverse dimensions and . This allows us to write the system of differential equations as a one-dimensional effective problem along the coordinate perpendicular to the electrode planes. Parallelplate chamber geometry is considered ideal (i.e. electrodes are considered parallel). • The release of the electron-ions pairs in the chamber occurs homogeneously in space. • The simulated signal produced in the ionization chamber is caused by the drift of the charges towards the electrodes (i.e., primary electrons contributing to the measured signal are not considered). • The only contribution to the recombination considered is the volume recombination. Initial recombination, back-diffusion recombination and electron-ion recombination are neglected. • The drift of ionic species is summarized in only one specie for positive ions and another for negative ions. Drift velocity is considered equal to the average velocity in a given electric field. • Ion and electron drift velocity are function of the local electric field.
where ( = +, −, ) are the positive ion, negative ion and electron densities in m −3 . We denote as the electric field across the chamber. In this work, we have not considered charge carrier multiplication in the chamber, taking the contribution from ( , ) term as negligible. The equations take all the carrier species as being charged with the elementary charge denoted as . The source term ( , ) represents the number of charge carriers of either sign escaping initial recombination produced in the sensitive medium per unit volume and time in m −3 s −1 . The parameter in m 3 s −1 accounts for the volume recombination effect between positive and negative ions in the same way that in m 3 s −1 accounts for the positive ion and electron recombination. For this work, the parameters and have been considered constant, although it could be argued that they may exhibit some dependence on the external electric field applied.
The electron attachment process is governed by in s −1 . The attachment coefficient (inverse of the attachment time ) has been taken as a function of the electric field. The coefficients ( = +, −, ) in m 2 s −1 model the diffusion of charge carriers inside the chamber. Finally, is the drift velocity of electrons in m s −1 and + and − in m 2 V −1 s −1 denote the mobility of positive and negative ions, respectively. In these equations, we have not included the de-attachment process and consider this effect to be negligible for the scope of this work.
The electric field dependence across the chamber can be computed by solving the following one-dimensional equation: where is the dielectric constant of the medium. The bias voltage provided to the chamber sets an additional boundary condition for the electric field. Assuming an ideal behavior of the high voltage supply and neglecting other additional electrical effects (i.e. cable, connector, electrometer and chamber impedance), the electric field between electrodes separated by a distance should fulfill the following boundary condition where is the high voltage applied to the ionization chamber. The instantaneous current density in A m −2 produced at any time may be calculated using the Shockley-Ramo [15] theorem in the following way: ] .
The numerical solution of this system of equations was realized via a Cython (version 0.29.28) code in which the equations were discretized in position and time. Similar approaches were proposed in order to reproduce computationally the behavior of the ionization chambers [16,17]. The discretization parameters were optimized to keep the error of the CCE below 0.1 % typically, with = and ≤ usually with = 3000 (see Fig. 1). This means that is dynamically set in each simulation iteration according to , which denotes the maximum speed of all carriers present in the area between the chamber electrodes at a given time.
For the purpose of analyzing the convergence of the numerical solution of the system of partial differential equations, a simulation were carried out by introducing the charge release in the medium instantaneously, turning off the electric field perturbation and considering only two carriers (positive ions and negative ions). In this situation, we can compare the charge collection efficiency via the numerical model against the evaluation using Boag's formula. Because this analysis does not involve the system as a whole, it allows us to obtain an absolute evaluation of the numerical error. Based on Fig. 1, 3000 spatial steps are sufficient to obtain a marginal numerical error compared to the uncertainties on the transport parameters.
The model can be used to estimate any characteristic of a PPIC such as the CCE or the free electron fraction for a pulsed or continuous beam.

Model parameters
One of the major problems we encounter when working with a computational model designed to obtain a satisfactory prediction of the experimental data is that the values of the physical variables involved must be known. Thus, in the model described above, we have diffusion, recombination and attachment coefficients together with mobility parameters for each charge carrier. The feasibility of the results of this model will be crucially dependent on the parametrization of these coefficients.
We used results from the work of Boissonnat et al. [18] to set the attachment parameter of the model. Average ion mobilities in air were taken from the work of Zhang et al. [19], in which full dependence on pressure, temperature and air humidity was considered. The Boissonnat mobilities were taken as an input parameter of the model to compare the simulation results against the experimental data. The electron drift velocity was simulated using the Magboltz code [20] (version 11.14).
Considering the conventional Nernst-Townsend formula, the diffusion coefficients can be written as = (273.2+ ) , where is the Boltzmann constant and is the mobility. Although the diffusion coefficients of ions are three orders of magnitude smaller than those of electrons, the standard deviation, due to the space diffusion after traveling a distance , would not be dependent on the mobility in a first order approximation and would amount to By contrast, the analogous numerical diffusion contribution [17] for a gap discretization in steps of the distance between electrodes has an approximate expression of Taking ℎ = , this leads to the relationship for the step distance where we have considered a reference numerical figure corresponding to a voltage bias of 300 V evaluated at 20 • C. Using a distance between electrodes of 1 mm, the number of steps is equal to As a result of the limit with equal numerical and physical diffusion, the number of steps needed to discretize a chamber does not depend on the distance between electrodes but only on its voltage bias . Normally, this discretization limit would yield long computing times for the simulation. For this reason, larger spatial step sizes with ≈ 3000 are applied, thus provoking numerical diffusion to dominate the overall diffusion contributions. In this scenario, we have neglected to first order the physical diffusion contribution in Eq. (1).
It should be noted that Nernst-Townsend formula is only valid when the charge carriers are in thermal equilibrium with the surrounding molecules. Considering the electric field regime in which established ionization chambers usually operate, this is true for positive and negative ions but not for electrons [21]. It is to be expected that the diffusion due to electrons is considerably greater than that predicted by the Nernst-Townsend relation and, therefore, that the number of steps necessary to equalize both physical and numerical contributions in the case of electrons is lower.
The volume recombination coefficient in air has values in the literature [17] from 0.9 × 10 −12 m 3 s −1 up to 1.3 × 10 −12 m 3 s −1 ; however, this coefficient was set to a value of 1.25 × 10 −12 m 3 s −1 for theoretical examples. The volume recombination coefficient may also exhibit some dependence on the electric field, temperature and humidity [22]; this dependence has not been included in this model. The recombination between ions and electrons is not considered here, as its contribution is estimated to be within the uncertainties of the transport parameters.
In the following simulations, unless not explicitly mentioned, standard temperature, pressure and humidity (20 • C, 101.325 kPa and 50 % relative humidity) were used.

Free electron fraction
As shown above, the free electron fraction (FEF) parameter plays a major role in describing the behavior of PPICs. In addition, Boag derived a simple expression that allows the calculation of the FEF, , if the characteristic electron properties are known: where is the total free electron charge arriving to the positive electrode and 0 is the total charge of either sign escaping initial recombination produced by the ionizing radiation per pulse. However, due to the fact that, in the UHDR regime, the electric field disturbance is very important, the FEF is no longer constant. In fact, by varying the dose per pulse to which the PPIC is exposed, it is possible to observe significant changes in the FEF.
The FEF is defined with respect to the total charge released by ionization that escapes initial recombination. We can also define a fraction in relation to the charge actually collected, which we denote as ′ . The relationship between these two quantities is given by: where ′ denotes the ratio of the free electron charge that reaches the positive electrode to the total collected charge , with being the standard definition (see Eq. (9)) in Boag et al. formalism. The parameter corresponds to the inverse of the CCE. The behavior of these two fractions for different PPICs with different electrode distances (as PPC05: 0.6 mm, Advanced Markus: 1 mm, Roos: 2 mm) can be seen in Fig. 3A and B. Even for the smallest gap chamber of 0.6 mm, the simulated starts to decrease significantly when the dose per pulse exceeds 0.1 Gy. The Boag FEF monotonically reduces its value with the dose per pulse in contrast to the original Boag model, where is a constant (dashed lines).
If the free electron fraction ′ is considered, this fraction tends to increase at very high doses per pulse. This is due to the fact that, when large volume recombination is present, the ion collected charge decreases more dramatically than the corresponding free electron charge. For chambers with small gaps or large electron fractions, ′ reaches its minimum at a dose where ion recombination losses begin to be dominant (Fig. 3B).
The free electron fraction used in Boag's formalism does not correspond to the contribution to the total charge actually collected due to the presence of free electrons in the ionization chamber. The measured charge is the result of the electrically induced signal of the charge  [23,24] and Reid [25] and the parametrization of Laitano [26] are shown for comparison. carriers that can be calculated using the above-mentioned Shockley-Ramo theorem. In this case, the ratio of the signal induced by the electrons to the total charge released in the medium that escapes the initial recombination is given by: As in the previous case, it is convenient to define the alternative quantity ′ based on the total collected charge: The dependence of these two quantities, and ′ , on the dose per pulse can be seen in Fig. 3C and D. Solid lines represent the full simulation of the induced charge due to free electrons in the ionization chamber. While Boag's free electron fraction tends to unity (Fig. 3A) for a very small distance between electrodes, the signal induced due to electrons would contribute to half the chamber readout value for uniform irradiation (Fig. 3C). The quantity ′ represents the actual fractional contribution of the drifting free electrons to the electrometer readout.

Pulse time structure
In many cases, the ion recombination correction factor is considered as an overall function of the dose per pulse delivered by the accelerator. One of the underlying questions in evaluating recombination effects is the extent to which the recombination is affected by the pulse-time structure. The dose rate variation during the pulse delivery may have a significant effect on the actual collected charge. Due to the possibility large variations of the dose rate within the pulse duration, we have chosen a simple test template. We considered the fact that the pulse of duration is divided into two subintervals with different dose rate as (see right part of Fig. 4): The first subinterval delivers a dose ratė1 in a time interval 1 and the seconḋ2 in 2 . This allows us to quantify the impact of the pulse structure in a simple way and show its impact in the final collection efficiency. Thus, the non-dimensional delta parameter controls the shape of the pulse. When = 0.5 the dose rate is constant over the J. Paz-Martín et al.  entire pulse duration while having a short and ultra high dose rate when → 0 and → 1 at the beginning or end of the pulse, respectively. The results were also analyzed in terms of different pulse durations .
The results for a 2 mm electrode distance operated at 300 V with a dose per pulse of 0.5 Gy can be seen in Fig. 4. For the pulse duration of conventional medical accelerators, which is typically less than 5 μs, the variation of the charge collection efficiency is below 3 % when the pulse structure is changed from an instantaneous dose delivery to a constant dose rate along the pulse. However, when using pulses with longer durations, relative variation of the ion recombination correction factor of the chamber can be dramatic (i.e., around 20 % for 100 μs).
It is important to realize that Boag's models do not take into account either pulse duration or pulse structure. Thus, it is important to verify the variation in charge collection efficiency of a chamber when the pulse duration/structure is altered in a specific irradiation facility. Due to the relatively low ion mobilities, variations within the pulse duration over periods of time shorter than the ion drift time do not significantly affect the chamber performance. However, when the duration of the pulse approaches the ion drift time , the recombination effect can exhibit a large variation.

Experimental validation
Some of the existing experimental data on the recombination effect concern only ion recombination correction factors or CCE [8,27]. A detailed numerical model is also able to provide a dynamic solution for the electrically induced current in air ionization chambers. Our study of the chamber electrical current description is much more demanding than a study of the overall balance of charge. The instantaneous current (i.e., the time resolved current during and shortly after the irradiation pulse), depends not only on the density of electrons, negative and positive ions but also on their mobility and on the accurate description of the electric field during the carrier drift. For this reason, the experimental measurement of the instantaneous current constitutes a benchmark for numerical models.
For our experiments, we used one unit of a Roos chamber (Serial No. 2981) irradiated using PTB's ultra-high pulse dose rate reference electron beam [28,29] with a nominal electron energy of 20 MeV, doses per pulse ranging from 0.2 Gy to 5.8 Gy and a fixed pulse duration of 2.5 μs. The chamber was operated at either a negative and a positive polarity from 100 V up to 450 V. The time-resolved beam current during each electron beam pulse was measured by means of an integrating current transformer (ICT from Bergoz instrumentation) connected to a transient recorder (Spectrum M3i.4142). The instantaneous current from the Roos chamber was recorded via a time-resolved measurement of the voltage drop across a 50 Ω resistor connected directly at the chamber output in combination with a preamplifier (FEMTO HVA-200M-40-F) connected to the second input channel of the same transient recorder as used for beam current measurement. The dose rate (̇( )) at the position of the ionization chamber was obtained from the instantaneous beam current ( ( )). This was achieved via a previous cross calibration of the charge per beam pulse (∫ 0 ( ) ) against the dose per pulse measured by means of alanine pellets at the position of the ionization chamber [30]: where denotes the dose per pulse to water at the reference depth of the ionization chamber and denotes the pulse duration. This dose rate profile was introduced in the numerical simulation in order to predict the chamber current via the Shockley-Ramo theorem. The nominal collecting electrode area of the chamber, whose dimension is 1.91 × 10 −4 m 2 , was used to obtain the absolute value of the instantaneous current in order to compare it to the measured timeresolved current. For the conversion of the dose rate into a liberated charge inside the medium per unit of time, the absorbed dose-to-water calibration coefficient , was measured using PTB's calibrated Co-60 source [31] and the corresponding beam quality correction factors were calculated using a Monte Carlo model of the beam used [28]. Fig. 5 shows examples of the experimental data. Graph A shows an example of the beam current measured by means of the ICT and the simultaneously measured instantaneous current, which is produced in the ionization chamber. Graph B shows the instantaneous currents at different doses per pulse. Graph C shows the effect of the operating voltage on the ion collection time.
The instantaneous current can be divided into two distinct parts. The fast current is almost able to follow the intra pulse beam current variations, aside from the fact that the fast current always higher than the beam current at the beginning of the pulse. The existence of the fast current part is due to the high mobility of the free electrons in the medium. Subsequently, after the end of the beam pulse, a small and long-lasting tail is observed. This is due to the slow ions, whose mobility is several orders of magnitude lower than the mobility of the where denotes the induced current due to electrons, denotes the current induced by ions (positive and negative) and the total induced current. The EFEF is the integral of the instantaneous current up to the end of the 2.5 μs beam pulse over the integral of the total induced current. The same procedure was followed for the simulated current curves. The relative contribution of the instantaneous positive and negative ion current to is always lower than 10 %. For the simulation using the Zhang et al. ion mobilities, a volume recombination parameter of 1.15 × 10 −12 m 3 s −1 was chosen as it showed good agreement between the simulated and experimental data. Fig. 6 compares the experimental instantaneous current and the simulated current. The contribution of each charge carrier to the signal is also plotted. On the right, the integral of the simulated and experimental signals is shown. The dashed area shows the part of the integral used for the calculation of the EFEF.
The comparison between the experimental and simulated EFEFs can be seen in Fig. 7. The definition of the effective free electron fraction represents a figure of merit for the electron current during the pulse delivery. Experimental results indicate that the free electron fraction is independent of the bias voltage for the dose per pulse over 2 Gy. Fig. 7 includes experimental results that can be compared to Fig. 3D for the 2 mm electrode distance that corresponds to a Roos chamber. The results for in Fig. 7 differ from ′ shown in Fig. 3D by less than 10 % due to the integration of the ion-induced current during the radiation pulse. These results show an unequivocal variation of the EFEF with the dose per pulse. This is also supported by a numerical simulation that shows qualitatively good agreement with the experimental data. The differences may be attributed to the uncertainties of the transport parameters used.
The ion collection time was defined here as the time required in order to collect 99.95 % of the total charge (dashed line in Fig. 5C). The results of the simulation using the Boissonnat and the Zhang et al. ion mobility measurements are shown in Fig. 8 compared to the experimental values.
It can be observed that the collection times obtained from the simulations differ by up to 20 % from those measured experimentally, suggesting that possibly a value between both used parametrizations may be more adequate than the effective positive ion mobility (i.e., the ions with lower mobility). A particularly remarkable finding is that the ion collection times show no significant variation with the dose per pulse (each point in Fig. 8 displays the average and uncertainty bar corresponding to the dose per pulse from 0.2 Gy up to 5.8 Gy), although they are affected differently by the dose-dependent electric field perturbation inside the ionization chamber and entangled with the volume recombination parameter.
Many of the numerical models considered previously have only considered the recombination effect [13,16,17]. In the previous paragraphs, we have seen that the correct prediction of the instantaneous current is a much more demanding task and requires the proper calculation of the charge carrier densities, the electric field perturbation and drift speed. Alternatively, the results comparing the modeled and  experimental CCEs can be seen as a function of the dose per pulse in Fig. 9 for several positive operational voltages. The CCE was calculated using the following expression: where is the beam quality correction factor obtained from [32], 60 Co , is the absorbed dose to water calibration coefficient in 60 Co and is the field size correction factor for the non uniformity of the beam profile. The agreement shown in Fig. 9 is, in general, better . This suggest that the ion recombination correction may be calculated for a Roos chamber using the numerical approach described in this work.

Discussion
The results obtained via the numerical model showed good agreement with the experimental data from a Roos ionization chamber in UHDR; the numerical solution of the transport equations is able to qualitatively reproduce the effects observed on the instantaneous current as seen in Fig. 6. We have defined an alternative effective free electron fraction based on the contribution of free electron induction to the total charge collected in the ionization chamber. This quantity differs from the original Boag definition, which amounts to the fraction of free electrons produced in the gas reaching the positive electrode. To evaluate the effective free electron fraction, we integrated When evaluating the ion collection times, we were able to benchmark the average ion mobilities used; there were discrepancies between the experimental and calculated ion times of around 20 %. In general, the average ion mobilities measured by Zhang et al. provided slightly better agreement with the experimental determined ion collection times and with the instantaneous current of the ionization chamber than the Boissonnat parametrization as shown in Fig. 8. Of particular interest is the fact that the ion collection times follow a linear relationship with the inverse of the operational voltage, even when the electric field is strongly perturbed. For large doses per pulse greater than 2 Gy, ion collection times tend to be independent of the dose per pulse. In general, the discrepancies observed between the simulation and the experimental data can be attributed to the uncertainty of the charge transport parameters used in the model such as the average ion mobilities in air, the uncertainties of the nominal value of the area and distance between electrodes and the simplifications assumed in the model.
The polarity effect has a well know contribution to the charge readout of parallel plate ionization chambers [33] either as a Compton current or voltage dependent contribution. Unfortunately the version of the numerical model used in this work does not accommodate this effect and the data was compared with the simulation separately for each polarity without further correction. Small polarity differences can be seen in the comparison with experimental data both for the EFEF (Fig. 7) and CCE (Fig. 10). In general terms, the average experimental relative discrepancy between chamber readout at opposite polarity for similar dose per pulse is around 1.7 % for our data set. Considering the uncertainties of the experimental data it is not possible to establish a significant voltage dependence of such polarity effect. Nevertheless, due to the large dose per pulse delivered, we cannot discard that, depending on the high voltage supply and electronic readout chain, some perturbation of the drift electrode, collecting electrode and guard ring voltages may appear.
The evaluation of the contribution of ion and free electrons to the ionization chamber readout via the numerical model has been performed using the Shockley-Ramo theorem. To this end, we have defined the contribution of free electrons to the chamber signal with respect to the saturation charge and also ′ with respect to the collected charge. Predictions for these quantities are shown in this work as a function of the dose per pulse and the distance between electrodes for a PPIC. Our simulation could potentially be used to study the effect of the pulse time structure on the chamber response. To this end, we have taken a simple approach dividing a given dose per pulse into in two time subintervals with different instantaneous dose rates. Whenever the delivery time is much smaller than the ion collection time (i.e., below 3 μs for a Roos chamber), the pulse time structure effect in the collected charge is negligible. Alternatively, if the duration of the radiation pulse is comparable with the charge carrier drift time the effect of the pulse time structure can be significant.

Conclusions
In this paper, the design and results of a computational model for evaluating the performance of a plane-parallel plate ionization chamber under pulsed beams have been presented. We have used an electron UHDR beam for the experimental validation of this numerical model. Instead of only testing the charge collection efficiency as a summarized way to benchmark the model [17], the electrical instantaneous current due to the charge carrier drift in a Roos chamber was also measured up to full the ion collection time. The instantaneous dose rate in the chamber, which was needed for the numerical simulation, was taken from the non-destructive time-resolved beam current measurement along with the alanine calibration of the dose to water. This methodology allowed us to obtain the expected time-dependent instantaneous current from numerical simulations and to compare it with the experimental data. This method proved to be a more demanding test for this type of models being sensitive not only to the produced charge carrier densities and recombination coefficient, but also to the drift velocities of ions and electrons in the local electric field inside the chamber. The total ion collection time has small variation with the total dose per pulse (from 0.2 Gy up to 5.8 Gy) and exhibits a linear dependence with the inverse of the chamber bias voltage in spite of the large perturbation of the electric field inside the chamber. Moreover, this numerical approach can be used to study the effect of the time structure of the radiation delivery on the chamber readout.
We have also reviewed the Boag definition of the free electron fraction which does not correspond to the fraction of the total collected charge due to the free electron drift inside the chamber that we have denominated . To verify experimentally this contribution we have measured the integrated charge during a pulse of 2.5 μs, dominated by the free electron component, produced at the PTB electron accelerator facility MELAF and compared it to the total charge collected in a Roos chamber. This ratio, that we denominate the free electron fraction, depends on the dose per pulse and the agreement between the numerical model and measured data is in general better than 5 %.
This work demonstrates the great advantage of numerical tools for updating the current understanding of ionization chambers and describing their behavior in exhaustive detail, even in extreme regimes such as those of UHDRs. Specifically, this work supports the possibility of calculating the ion recombination correction factors in regimes where Boag's theories are not valid using the numerical approximation described. Unfortunately, the poor knowledge of the transport coefficient currently limits the accuracy of these models. Thus, further investigation of such quantities is required in order to produce a high precision tool for the medical physics community.