Initial recombination in the track of heavy charged particles: Numerical solution for air filled ionization chambers

Abstract

Introduction. Modern particle therapy facilities enable sub-millimeter precision in dose deposition. Here, also ionization chambers (ICs) are used, which requires knowledge of the recombination effects. Up to now, recombination is corrected using phenomenological approaches for practical reasons. In this study the effect of the underlying dose distribution on columnar recombination, a quantitative model for initial recombination, is investigated. Material and methods. Jaffé's theory, formulated in 1913 quantifies initial recombination by elemental processes, providing an analytical (closed) solution. Here, we investigate the effect of the underlying charged carrier distribution around a carbon ion track. The fundamental partial differential equation, formulated by Jaffé, is solved numerically taking into account more realistic charge carrier distributions by the use of a computer program (Gascoigne 3D). The investigated charge carrier distributions are based on track structure models, which follow a 1/r² behavior at larger radii and show a constant value at small radii. The results of the calculations are compared to the initial formulation and to data obtained in experiments using carbon ion beams. Results. The comparison between the experimental data and the calculations shows that the initial approach made by Jaffé is able to reproduce the effects of initial recombination. The amorphous track structure based charge carrier distribution does not reproduce the experimental data well. A small additional correction in the assessment of the saturation current or charge is suggested by the data. Conclusion. The established model of columnar recombination reproduces the experimental data well, whereas the extensions using track structure models do not show such an agreement. Additionally, the effect of initial recombination on the saturation curve (i.e. Jaffé plot) does not follow a linear behavior as suggested by current dosimetry protocols, therefore higher order corrections (such as the investigated ones) might be necessary.

For characterization of radiotherapy beams plane parallel ionization chambers (ICs) are widely used. IC dosimetry is accepted worldwide to be one of the most accurate dosimetry systems, which can be applied on a daily routine base. Correction factors have to be applied to correct for the effect of several influence quantities, such as the beam quality or perturbations introduced by the IC itself. A widely used protocol for this purpose is the TRS-398 [1] of the International Atomic Energy Agency, Vienna (IAEA). Besides that the electrometer reading (i.e. the charge collection) is decreased due to recombination losses within the IC. This recombination taking place in air filled ICs is currently corrected for by the use of the two voltage method or by determining recombination effects by plotting the reciprocal of the chamber reading against the reciprocal of the polarizing voltage (Jaffé-plot).

TRS-398 suggests using the two voltage method for correction of general recombination. For initial recombination in ion beams, in case of a negligible general recombination, TRS-398 recommends a Jaffé plot to determine the saturation charge or current, and hence the recombination correction factor (see page 129 in [1]). Initial recombination can, however, also be modeled by more fundamental theories, such as Onsager theory [2], the cluster theory [3], which are both valid for low linear energy transfer (LET) radiation, and the theory of columnar recombination elaborated by Jaffé [4] in 1913, which is the only theory of these describing the processes taking place in the track of a heavy charged particle (HCP, here: proton, helium, etc.) with a high LET. Jaffé's theory was previously used for estimation of recombination effects in air filled ICs exposed to heavy charged particles with high LET [5]. Since its original formulation the theory has not changed and no extensions were made.

The objective of this study is to evaluate possible extensions of Jaffé's theory. This is done by the use of different radial dose distributions (RDD) around ion tracks as initial value of the diffusion and recombination process. This process is modeled as partial differential equation (PDE). In the initial formulation of Jaffé's theory of columnar recombination only one of these RDDs was taken into account. In this paper the impact of the different RDDs on the result of the diffusion and recombination PDE is examined. Further, the results are compared to experimental data.

Material and methods

Jaffé's theory of columnar recombination

Columnar recombination occurs between charged particles created by one incident HCP. Here, the incident particle passes the detector and ionizes the sensitive material in its track. This cloud of positive and negative charges is assumed to be cylindrically symmetric around the center of the track. After ionization the produced charge carriers undergo diffusion and recombination. Therefore, the basic PDE is

where n is the charge carrier density, D the diffusion constant, α the recombination coefficient, and t the time.

Based on this equation Jaffé [4] calculates, for the case that the electric field is parallel to the trajectory of the incoming HCP, the collection efficiency with respect to initial recombination to be

where i_∞ is the saturation current or charge, k the mobility of the charge carriers, d is the length of the columns (or the electrode spacing), E the electric field, i the measured current or charge and N₀ the total number of charge carriers of one sign per unit length. The expression li(e^x) is the logarithmic integral of x. For Equation 2 it is assumed that the distribution of the charge carriers follows a Gaussian style charge carrier distribution (see Equation 3) with the width of b² = 2ς², centered at the trajectory of the incoming HCP. This parameter b can be taken as fitting parameter. More details can be found in Appendix A and in references [4] and [6]. For the energy interval of ions encountered in clinical applications the differences between the numerical calculation of the collection efficiencies based on Equation 1 and the closed solution (Equation 2) are insignificant [7]. Therefore, the closed solution may be used instead of a numerical solution of Equation 1, which is done in this paper.

Initial charge carrier distributions

Jaffé used for his theory a Gaussian type initial charge carrier distribution (CCD) which is a solution of Equation 7. Nevertheless, experiments by Varma et al. [8–10] employing gas filled detectors have shown that in the “penumbra” the decrease in dose is proportional to 1/r². These findings led to a variety of different models and theories of RDDs, which describe the dose distribution around the track (i.e. trajectory) of the HCP (see e.g. [11]). These RDDs can directly be translated into radial CCDs by multiplication with eρ/W, where the W is the mean energy to produce an ion pair (see ICRU report 31 [12]), ρ the density of the medium (in this case air), and e the elementary charge.

For our calculations we used two basic amorphous track structure models which describe the RDD around an ion track. The first model is based on the assumption of a constant core and a continuous transition between core and penumbra. This CCD is used for the prediction of the biological effects in radiotherapy using HCPs (see [13–17]). This RDD is piecewise defined by Equations 8 to 11 in Appendix B. In the following we refer to this CCD as “Scholz–Kraft (SK) CCD”. In Figure 1 a qualitative comparison between the SK CCD and the later presented CCDs can be seen.

Figure 1. Comparison between the different charge carrier distributions used. The Kiefer–Chatterjee (KC) and the Scholz–Kraft (SK) charge carrier distribution share the same r _min = 4 nm, the width of the Gaussian style charge carrier distribution is b = 20 μm.

The other RDD used is the Kiefer-Chatterjee parameterization (Equations 13 and 14 in Appendix B, see also [18]), which has been used for the modeling of the response of diamond detectors [19]. Here, the transition between core and penumbra is not continuous. Based on considerations from collision dynamics, the dose in the penumbra D_p is calculated and the dose in the core D_c is then normalized to the total energy loss. In the following we refer to this CCD as “Kiefer–Chatterjee (KC) CCD”.

The core radius r_min was taken in both cases as free parameter, as up to now, no experimental determination of the core radius r_min was found in the literature. Nevertheless, some models provide an equation or a value for the core radius r_min (see Equation 15 in Appendix B).

The last CCD used is a Gaussian-style initial charge carrier distribution.

Using this charge carrier distribution the collection efficiency with respect to initial recombination f_I can be calculated quickly for clinically relevant parameters using the relation presented in Equation 2 (see [7]). This CCD was already introduced by Jaffé in the original formulation of the theory. For the calculations the closed solution (Equation 2) was used, instead of the procedure presented in the next section, which was used for the KC and the SK CCDs.

Gascoigne 3D

Gascoigne 3D is flexible toolkit for numerical solving of PDEs [20]. It is based on a finite element analysis and is usually used for the modeling of combustion problems and stationary and dynamic viscous flows. It uses adaptive mesh refinement and multigrid solvers. The package is written in a modular structure in C++. Equation 6 was implemented in the program. A fixed grid of 440 nodes was chosen, placing the nodes equidistant on a log(r) scale. Besides that, the radius r and the charge carrier density n were regarded as log(r) and log(n) which required some changes of Equation 6 in its implementation. Additionally, the width of the time steps Δt can be arbitrarily chosen by the user. For these calculations a pseudo logarithmic (“stairway”) time step width was chosen.

Using this program the diffusion and recombination process perpendicular to the trajectory of one single HCP is modeled solving numerically Equation 6 for a given charge carrier density n based on the dose distributions presented in the previous section as initial condition. The total number of charge carriers is later calculated for every timestep. The effect of the electric field on the collection efficiency f_I with respect to initial recombination is then calculated by the formula

where Q₀ is the total number of charges initially produced by the incident particle and Q(t) or Q _n are the number of not recombined charges at time t or time-step n, which is calculated at every time-step according to Equation 1. The time to separate the columns of charge carriers is . Using this method only the recombination within one particle track is modeled, no inter track recombination or track overlap is taken into account.

Ionization chambers

In our experiments we used a Bragg peak ionization chamber (BPIC) which is a plane parallel IC with a disk shaped sensitive volume with an electrode spacing of 2 mm. The electrode radius is 40.8 mm. The entrance window has a thickness of approximately 3.5 mm, whereas the exit window has an approximate thickness of 6 mm. The housing of this IC is made of PMMA, thus these thicknesses correspond to water equivalent path lengths of 4 mm and 7 mm, respectively. This ionization chamber is mainly used for dosimetry in proton and carbon ion beams. The recommended chamber voltage is 400 V which corresponds to an electric field of 2000 V/cm. Nevertheless, for this experiment the voltage was varied between 200 V and 800 V. More details on the Bragg Peak IC can be found in its manual [21].

Additionally a specially made advanced Roos Chamber (PTW, TM34073-1, 80-0001) was used. This chamber has an entrance window thickness of 1.1 mm and an exit window of approximately 1.5 mm. The electrodes are disk shaped with a spacing of approximately 2 mm and a radius of 19.8 mm. Here, the recommended voltage of 400 V was used. A more detailed description on the use of this chamber may be found in [22].

Both ICs were read out by remotely controlled Unidos electrometers (Unidos webline, PTW, Freiburg, Germany). The advanced Roos IC was connected directly to the 400V HV provided by the electrometer, whereas the BPIC was connected with an adapter to an external high voltage power supply SHQ224M (ISEG, Radeberg/Rossendorf, Germany).

Irradiation with carbon ion beams

The experiments were performed at the Heidelberg Ion Therapy Center (HIT) located in Heidelberg, Germany. There, carbon ions are accelerated by a synchrotron to energies between 88 MeV/u and 430 MeV/u. Six different spot size widths can be chosen by the user resulting in spot sizes from approximately 4 mm (FWHM) to approximately 20 mm. The intensity is also variable and 10 different intensity levels can be chosen by the user. For this experiment the smallest intensity was chosen to minimize any recombination effects caused by general recombination. For the same reason also a spot size of 10 mm was chosen (N.B. At the time when the experiments were performed only four different spot sizes were available. The used spot size was at that time the widest spot size available.). Here, we present the results of calculations for clinically used (nominal) energies of 88 to 430 MeV/u, which were degraded through the setup to 75 and 420 MeV/u.

Experiments

For these experiments the BPIC was placed at a distance of 1423 mm from the vacuum window in a horizontal beam. In upstream position, i.e. facing to the nozzle, the advanced Roos IC was placed next to the BPIC. This IC was used as beam monitor and to later normalize the electrometer readings of the BPIC. The advanced Roos IC was chosen as monitor to avoid too high changes in the mean stopping power (or unrestricted LET) caused by the thicker exit window of the BPIC.

Fitting

Saturation correction of experimental data. From the measured data the relative charge, i.e. ⁱField/ⁱMonitor was calculated dividing the electrometer reading of the field chamber (BPIC) by the reading of the monitor chamber (advanced Roos IC). The reciprocal of the relative saturation charge a was determined using Equation 5

where a and b are fit parameters and V the applied voltage. Later, the inverse collection efficiency 1/f_I was calculated, for each data point, dividing the reciprocal of the relative charge, i.e. by ⁱField/ⁱMonitor the value a. This was done as a could deviate from one. In the case of ideal conditions (e.g. monitor at the same position as field chamber with the same calibration factor and the same dimensions of the sensitive volume, etc) this step would not be necessary.

Fitting calculated data to experimental data. The numerically calculated data for the inverse collection efficiency with respect to initial recombination 1/f_I was fitted to the inverse collection efficiencies obtained in the experiment. Fit refers here to the best estimate for the KC and SK CCDs. This method was chosen due to the long computational times of several days per curve for the KC and SK curves. For Jaffé's model, i.e. Equation 2, the nls fitting algorithm of the language R [23] was used to determine the best fit.

Results

For several residual energies and electric fields the inverse collection efficiency with respect to initial recombination 1/f_I can be seen in Figure 2, and Figures 5 and 6 in the Appendix, for experimental and calculated data A summary of the values of the core radius r_min or width b can be found on table 1.

Figure 2. Numerical solution of the diffusion and recombination partial differential equation (PDE) for 420 MeV/u carbon ions with a linear energy transfer of LET = 9.4 keV/μm using the Gaussian (JAFFE), the Kiefer–Chatterjee (KC), and the Scholz–Kraft (SK) charge carrier distribution for the best fitting b or r_min (indicated in the plot) as initial condition compared to measurements. The numerical calculations are displayed as lines whereas the experimental data is displayed as open circles.

Table I. Approximate fit parameters of the core radius r_min (KC and SK) or width b (Gaussian CCD, see Equation 3) and the correction in the relative saturation charge extracted from Figures 2, 5 and 6.

Energy [MeV/u]	420	240	75
LET [keV/μm]	9.4	13	29
KC fit to absolute values	6 μm	10 μm	25 μm
KC fit to curvature	0.1 μm	0.4 μm	20 μm
KC resulting correction	1.2%	1.4%	0.1%
SK fit to absolute values	30 nm	40 nm	700 nm
SK fit to curvature	0.1 nm	0.1 nm	40 nm
SK resulting correction	1%	1.6%	1.35%
Gaussian fit to absolute values	12 μm	15 μm	28 μm
Gaussian fit to curvature	6.8 μm	13 μm	27 μm
Gaussian resulting correction	0.43%	0.09%	0.02%

KC, Kiefer–Chatterjee; LET, linear energy transfer; SK, Scholz–Kraft.

The values for the best fitting core radius r_min or width b to the absolute values of the calculated 1/f_I curves are in the range from 6 to 28 μm for the KC and the Gaussian-style CCD, whereas the values for the SK CCD vary between 30 and 700 nm. For the fits to the curvature the values follow a similar pattern, but all the values are smaller: The values for the KC and the Gaussian-style CCD vary between 0.1 μm and 27 μm and between 0.1 nm and 40 nm for the SK CCD.

In Figures 3 and 4 an example for the evolution of the diffusion process can be seen. Here, the best fitting core radius or width was chosen. It can be observed that all CCDs change their shape with time. Besides that, the KC looks similar to a Gaussian CCD, whereas the SK CCD becomes more and more Gaussian with time. Additionally, the height is reduced; nevertheless the integral of the area under the curves stays equal.

Figure 3. Diffusion of a Scholz–Kraft (SK) and Gaussian style charge carrier distribution for a 240 MeV/u carbon ion. The values for the core radius or the width are taken from the fit to the absolute values.

Figure 4. Diffusion of a Kiefer–Chatterjee (KC) and Gaussian style charge carrier distribution for a 240 MeV/u carbon ion. The values for the core radius or the width are taken from the fit to the absolute values.

Discussion

Form Figures 2, 5, and 6 it can be observed that the Gaussian charge carrier distribution already introduced by Jaffé in 1913 reproduces the values measured in the experiments. Neither the KC nor the SK CCDs based calculations reproduce the experimental data. Looking at the core radius r_min of the KC CCD and the width b of the Gaussian CCD for the absolute values, it can be observed that these values are similar. Whereas, the values obtained using the SK CCD differ from the values obtained with the Gaussian CCD and the KC CCD.

Taking the slope or curvature as reference instead of the absolute values, changes the result especially for the SK and KC CCD. The best fit to the experimental data can be achieved by reducing the core radii r_min several order of magnitudes, whereas for the Gaussian CCD the width b changes at the maximum by a factor of 2. This is also reflected by the factor with which the experimental values have to be scaled/multiplied to fit the calculated curve best. In two of three cases this factor is higher for the SK CCD.

This multiplication factor could imply that either the true ⁱField/ⁱMonitor saturation charge or current, or the true at infinite voltage, is higher than those determined with the recombination correction proposed by TRS-398 [1] and other dosimetry protocols. Nevertheless, in the case of the Gaussian style CCD, these shifts are small compared to the accuracy required in clinical dosimetry and comparable to the uncertainty in the calculation of the saturation current or charge.

For all charge carrier distributions the radius increases with LET or in other words the radius increases with decreasing energy. This behavior is not in agreement with theories typically used (see e.g. [24]). The width of the core should decrease upon decreasing kinetic energy as it is, for example predicted by Equation 15. Nevertheless, the observed behavior of an increasing radius or width with decreasing energy was also observed by Kanai et al. in [5] applying Jaffé's theory.

The encountered radii r_min or widths b of the KC and Gaussian CCDs, using the fits to the absolute values and to the curvature, are in the region of approximately 1 μm. Density scaling these values by a factor 800 (ratio of water and air densities) to water a radius in water of approximately 1.3 nm is estimated.

Besides that, the core radius r_min using the fits to the SK CCD data is in the order of 100 nm for the fit to the absolute values, density scaling such a value a core radius in water of 0.13 nm is estimated.

These core radii or widths obtained using the Gaussian and the KC CCD, as well as the SK CCD data, taking here only the fits to the absolute values, are in agreement with the data found in the literature. In for example [14–17] core radii varying from 0.3 to 40 nm are used for the calculation of biological effects. Taking the values of the SK CCD for the fit to the curvature, values in the order of 1 nm can be observed, density scaling these values would result in a core radius of 1.3 pm, which is less than the radius of an atom. This would imply a high charge carrier density of positive and negative charge carriers at the same region, with a radius smaller than an atom, which is not realistic.

The similarity of the values obtained using the KC CCD and a Gaussian charge carrier can be explained by the similarity between these two charge carrier distributions. The step or hat of the KC charge carrier distribution could be approximated by a Gaussian distribution as visible in Figure 4, thus the diffusion behavior is similar to those of a Gaussian charge carrier distribution neglecting the 1/r² decrease at wider radii (penumbra). In contrast to that, the shape of the SK CCD is not as similar to the Gaussian CCD, therefore a different diffusion behavior could be expected.

For the formulation of Jaffé's theory and the here presented extension the following assumptions were made, which could explain some of the observed effects:

i) Unique charge carrier mobility k for positive and negative charge carriers. This assumption was already made in the initial formulation of Jaffé's theory, however different values for positive and negative charge carriers (in air) are not known to the authors. The value of the charge carrier mobility used was taken from [5].

ii) Unique diffusion coefficient D for positive and negative charge carriers. This assumption was also made in the initial formulation. In the case of different diffusion coefficients D different values for the different ions and their charge (e.g. N⁺ , N²⁺ , O⁺ , O²⁺ , etc.) have to be taken into account. The value of the diffusion coefficient used was also taken from [5].

iii) The dose distributions in air are the dose distributions in water scaled by the ratio of the densities. The models of RDDs in water are based on measurements made in air of the penumbra scaled to water (see [8–10]). Models of dose distributions in water were chosen as input parameters as here the models are more established (see [11] or [14]). Nevertheless, measurements were only made in the penumbra region and the width and the form of the core region are only approximations to the real shape of the core [14].

iv) Charge carrier distributions are dose distributions multiplied with eρ/W. This assumption can be justified in the context of Amorphous track structure models [11], nevertheless it introduces an uncertainty.

At the same time the initial formulation of Jaffé's theory of columnar recombination, the best fitting theory in this investigation, is based on a Gaussian charge carrier distribution which has not been observed in this form in experiments. Experimental data found for example in [8–10] shows a 1/r² decrease in the penumbra, which lead to the development of more realistic radial dose distributions, such as the KC and SK radial dose distributions also investigated in this paper.

In conclusion, this study shows that the extension of Jaffé's theory by numerical solution of the underlying PDE, using radial dose distributions available from the amorphous track structure formalism, is not able to describe the recombination taking place in the track of a charged carbon ion with sufficient accuracy. Compared to this, the initial formulation of the theory of columnar recombination is able to describe the experimental data. For a final statement, the open questions presented in the previous paragraphs have to be investigated further and the simplifications made in the formalism have to be removed. Having such a description the response of air filled ICs may be predicted with higher accuracy if based on a theoretical basis and not only on phenomenological approaches as it is done today.

Acknowledgments

The authors want to thank Benjamin Ackermann and Stephan Brons for their help at the HIT facility. Acknowledgment is given to Thomas Richter for adjusting the code of Gascoigne to our needs, enabling us to perform the presented calculations. We also want to thank Steffen Greilich for his valuable remarks on this manuscript, helping us to improve this manuscript. Besides that, we also want to thank the reviewers and editors for their comments and suggestions to improve and complete this manuscript. This work is supported by the ”Kompetenzverbund Strahlenforschung” of the German Federal Ministry of Education and Research and CIRRO - The Lundbeck Foundation Center for Interventional Research in Radiation Oncology and The Danish Council for Strategic Research. Niels Bassler is supported by the Danish Cancer Society.

Declaration of interest: The authors report no conflicts of interest. The authors alone are responsible for the content and writing of the paper.