Precise measurement of the Bragg curve for 800 MeV / u 238 U ions stopping in polyethylene and its implications for calculation of heavy ion ranges

,


A
: Stopping power predictions in radiation transport codes are based on the Bethe-Bloch formula and different corrections. For very heavy ions at relativistic energies the available experimental data are scarce and therefore verification of stopping power predictions is only possible to a limited extent.
In this work, a full experimental Bragg curve for 800 MeV/u 238 U ions stopping in polyethylene is presented. The measurements were conducted at the experimental area Cave A at GSI Helmholtzzentrum für Schwerionenforschung in Darmstadt, Germany. The 800 MeV/u 238 U beam was provided by the SIS18 heavy ion synchrotron. The Bragg curve was measured with a setup consisting of a binary range shifter and two large area parallel plate ionization chambers. Complementary Monte Carlo simulations using the FLUKA code were performed and compared with the experimental Bragg curve. The mean ionization potential of polyethylene was fine-tuned to match the measured primary ion range with FLUKA simulations. The impact of the Bloch and Mott corrections to the stopping power calculation were studied by switching them off intentionally in separate simulations. A detailed description of the implementation of the stopping power formulae and the Mott correction in FLUKA is provided.
Supplementary data for this article is available online.

K
: Models and simulations; Simulation methods and programs; Detector modelling and simulations I (interaction of radiation with matter, interaction of photons with matter, interaction of hadrons with matter, etc); Gaseous detectors

Introduction
The accurate calculation of energy loss of heavy ions in media is important in many fields of application, such as heavy ion therapy [1], space radiation protection [2] and shielding design for high-energy accelerators [3]. Moreover, high energy uranium ions are a very useful tool to study single event effects in microelectronics as they provide a maximum linear energy transfer [4]. Modern radiation transport codes predict the energy loss and ranges of heavy ions with great precision. However, for very heavy ions like uranium at relativistic energies, where many different effects and corrections have to be taken into account for calculation of the stopping power [5], the available experimental data are scarce and therefore the verification of stopping power predictions is only possible to a limited extent.
The stopping power for heavy charged particles is commonly derived from the Bethe formula [6]. Various corrections to the Bethe formula have been proposed over the years to account for different effects: e.g. the effective projectile charge (re-capturing of electrons), shell corrections, the nuclear stopping power, the density effect [7], higher order Born approximation terms like 3 (Barkas-Andersen correction [8][9][10]) and 4 (Bloch correction [11]) as well as the Mott correction [12].
For high energy uranium beams, it has been pointed out that predicted stopping powers and ranges are particularly sensitive to the Bloch and Mott corrections [13][14][15]. The Bloch correction takes into account the saturation of energy transfer in collisions with small impact-parameter [14] and therefore decreases the predicted stopping power compared to the pure Bethe formula. The Mott correction describes the direct scattering of the projectile nuclei on electrons and therefore increases the predicted stopping power and also affects the production of high-energy delta electrons.
-1 -In this work, a full experimental Bragg curve for 800 MeV/u 238 U ions stopping in polyethylene is presented. The measurements were conducted at the experimental area Cave A at GSI Helmholtzzentrum für Schwerionenforschung in Darmstadt, Germany. The 800 MeV/u 238 U beam was provided by the SIS18 heavy ion synchrotron.
Complementary Monte Carlo simulations using the FLUKA code were performed and compared with the experimental Bragg curve. The mean ionization potential (I-value) of polyethylene used by FLUKA was fine-tuned to match the measured primary ion range. The impact of the Bloch and Mott corrections to the stopping power calculation were studied by switching them off intentionally in separate simulations. A detailed description of the implementation of the stopping power formulae and the Mott correction in FLUKA is provided.

Beam characteristics
The 238 U ions were accelerated to 800 MeV/u by the SIS18 synchrotron at GSI with the charge state 238 U 73+ . After passing through the vacuum exit window they are mostly fully stripped. The beam spot size was about 5 mm full width at half maximum. The absolute precision of the energy can be considered better than 0.04% for the SIS18 accelerator [14]. The intensity was about 5 × 10 6 particles per spill with about 5 s spill length. Two spills were requested per measurement point and for most thicknesses at least two independent measurements were taken to verify the robustness of the measurement.

Bragg curve measurement
A schematic of the experimental setup is shown in figure 1.
The adjustable polyethylene thickness for the Bragg curve measurement was realized using a remote controllable binary range shifter plus additional polyethylene plates and foils that were placed by hand. The areal densities of the binary range shifter plates were 0.005923 g/cm 2 , 0.011855 g/cm 2 , 0.02228 g/cm 2 , 0.044897 g/cm 2 , 0.090054 g/cm 2 , 0.17999 g/cm 2 , 0.36 g/cm 2 , 0.7211 g/cm 2 , 1.4411 g/cm 2 and 2.8827 g/cm 2 . For thicknesses that were too thick or too thin to be realized using the binary range shifter alone, offset plates and foils with the following areal densities were placed by hand directly in front of IC2 in order to measure the delta electron build-up: 0.01 g/cm 2 , 0.04 g/cm 2 , 0.08 g/cm 2 , 0.19 g/cm 2 , 0.38 g/cm 2 , 0.56 g/cm 2 , 0.94 g/cm 2 , 1.39 g/cm 2 , 1.87 g/cm 2 , 4.68 g/cm 2 , 7.03 g/cm 2 and 9.36 g/cm 2 . The estimated uncertainty of the areal density is 2% for the foils and 0.3-0.5% for the thicker plates.
The relative ionization was measured by one large area (20 20 cm 2 ) parallel plate ionization chamber [16] before (IC1) and another one behind (IC2) the polyethylene target. The IC1 had an electrode gap of 1 cm and the IC2 a gap of 0.5 cm. Both were filled with a mixture of 80% Ar and 20% CO 2 . The charge released in both ionization chambers was read out by Keithley K6517A electrometers and the ratio of the charges measured with IC2 and IC1 is a measure of the relative ionization or laterally integrated dose, respectively. The closest polyethylene plate or foil was placed directly against the IC2 in order to avoid the loss of any large angle fragments or delta electrons. The ionization chambers were operated at 1800 V. As visible in figure 2, the signal of the measurement ionization chamber (IC2) was far in the saturation plateau at this voltage so that initial or volume recombination effects in the detector gas are certainly negligible. The ratio of IC2/IC1 is ∼ 0.5 because IC2 had about half the gap size of IC1 (see above).

FLUKA Monte Carlo simulations
Monte Carlo simulations were performed using the FLUKA Monte Carlo code (version 2021.2.5) [17][18][19]. For the simulations, the polyethylene plates were approximated as a single thick block target. Lateral loss of particles and other geometry effects caused by the gaps between the plates (see figure 1) can be considered negligible because of the large acceptance of the experimental setup combined with the high rigidity of the uranium beam. The offset materials were placed in front of the polyethylene block: the vacuum exit window (0.027 g/cm 2 aluminum), IC1 plus half IC2 (0.0035 g/cm 2 mylar plus 0.03015 g/cm 2 electrode material consisting of 36% nickel and 64% mylar) and the air gaps between the components (0.15 g/cm 2 ). These offset materials add up to ∼ 0.21 g/cm 2 (0.18 g/cm 2 polyethylene equivalent). The uncertainty of this offset material thickness is estimated to be in the order of 5%. The laterally integrated absorbed dose inside the polyethylene absorber was scored with the USRBIN estimator. The energy spread was tuned to 0.3% full width at half maximum to match the measured distal falloff behind the Bragg peak with the simulation. This value is in the order of what is expected from the SIS18 synchrotron. The mean ionization potential (I value) of polyethylene was adjusted to match the measured 80 range (thickness at which the dose decreased to 80% of the Bragg maximum). The delta electron production and electron transport thresholds were set to 100 keV. A dedicated sub-routine allowed us to switch off single corrections used in the stopping power calculation implemented in FLUKA. Special attention was paid to the Bloch and Mott corrections.

Stopping power for "heavy" charged particles and Mott correction in FLUKA
The most important atomic processes undergone by charged particles when traversing media are related with Coulomb scattering with both atomic nuclei and electrons. Interactions with the latter give rise to the electronic stopping power which is the dominant process for particle energy losses down to the very low energies.
Energy losses of charged particles are commonly expressed as an average energy loss per unit path length, and by the associated fluctuations around the average value. The slowing down of protons and ions is governed by collisions with the atomic electrons and produces the characteristic shape of the depth dose profile for charged particles heavier than electrons, exhibiting the characteristic Bragg peak at the end.
A description of the implementation in FLUKA of the electronic stopping power and associated fluctuations for particles heavier than electrons can be found in [19]. Here a few reminders are given, together with a focus on the practical implementation of the Mott correction which was never described in detail before.
The formula for the average, unrestricted, energy loss of particles much heavier than electrons and with charge , can be expressed by:

JINST 17 P12019
for spin 0 particles and by: for spin 1/2 particles. FLUKA uses the spin 1/2 formulae for example for protons while for all ions heavier than 4 He the spin 0 formulae are used. is the projectile velocity relative to the speed of light, is the target material electron density ( = Av for an element), its mean ionization potential, is the projectile mass, , are the projectile kinetic energy and momentum, and max is the maximum energy transfer to a stationary electron, which is dictated by kinematics and given by ( √ is the center-of-mass energy): The terms , / , 1 , 2 and are all corrections to the Bethe-Bloch formalism. The "density correction" , and the shell correction , which are important at high and low energy, respectively, will not be discussed further. The reader can refer to the extensive literature about these corrections, and to ref. [19] for details about the FLUKA implementation. eff is the projectile "effective charge" which takes into account the partial neutralization of the projectile charge when its velocity is not much larger than those of the atomic electrons, the FLUKA implementation is briefly discussed in [19], as well as the implementation of the 3 (Barkas) [8], and 4 (Bloch) [11], corrections (indicated by 1 and 2 ).
The expression for the restricted energy loss of particles much heavier than electrons and charge , with energy transfers to atomic electrons restricted at is given by: for spin 0 particles and by a corresponding expression for spin 1/2 particles. The Mott correction, , is associated with the electron-ion Mott cross section [12], which is the Coulomb cross section for electron scattering on a charge computed beyond the first Born approximation. The expression for the stopping power of particles heavier than electrons is customarily derived making use of the distinction between "distant" and "close" collisions, where the names refer to the magnitude of the momentum transfer and hence to the collision impact parameter (even though in a quantum mechanical sense). In distant collisions the particle interacts with the atom as a whole, while for close collisions the interaction can be considered to be with free electrons, and atomic properties are not involved. The Bethe-Bloch equation, for the close collision part, is based on the electron-ion scattering cross sections computed in first Born approximation. Under this approximation the cross section for producing an electron of energy for an incident particle of velocity is given by: for spin 1/2 particles. is the classical electron radius. However whenever is no longer negligible ( is the fine structure constant), higher order corrections must be applied. The electron z-charge Mott cross section includes those corrections, however it is mathematically very complex, requiring the evaluation of several transcendental functions which are computationally cumbersome.
In FLUKA the Mott cross section parameterization proposed in [20] as further modified in [21] is used to compute the correction to the average stopping power, as well as the associated corrections to the secondary electron production cross section and to the energy loss fluctuations. Specifically, following refs. [20,21], the Mott correction to the Rutherford cross section is parameterized as a function of particle velocity and charge as: where is the electron scattering angle on the point charge . The ( , ) coefficients are in turn defined as: The Mott cross section is derived under the assumption of an infinitively heavy point-like target, hence one can safely confuse the scattering angle with the center-of-mass scattering angle * . For the derivation of the stopping power of heavy charged particles, the interest is in the reverse kinematics frame, where the electron is at rest and the ion of charge is the projectile. Going in reverse kinematics and using the relation among * and the target electron recoil energy ( √ being the center-of-mass energy): With a bit of mathematics, one can then derive the following expressions for the integrated cross section, and the first moment term which should be substituted into the stopping power formula, with -6 -1 = max for the unrestricted stopping power, while 1 = if secondary electrons are explicitly generated above : 14) The curly bracket terms are those coming from the adopted parameterizations for the Mott cross section and they should be compared with unity, the corresponding value for the Rutherford cross section. It should also be observed that the Mott cross sections reduces to the Rutherford one for → 0, and therefore the 0 ( ) coefficients must be construed in such a way to fulfill this condition. This constraint was not explicitly imposed in refs. [20,21], however their fit parameters are in substantial agreement with this condition. In FLUKA, the numerical parameters employed in ref. [21] have been adopted with a minimal re-normalization in order to assure 0 = 1 exactly. Under the usual assumption 0 , max , the Mott correction to the (un)restricted stopping power can be eventually expressed starting from eq. (2.14) as: It is important to note that the Mott correction, contrary to other ones, is dependent on the threshold for the explicit generation of secondary electrons. Obviously the correction for the unrestricted stopping power is obtained setting = max in eq. (2.15). The fluctuations associated with charged particle energy losses are also an important ingredient since they determine the shape and in part also the position of the Bragg peak. Indeed, its location does not correspond to the nominal range corresponding to the particle energy, but is situated slightly in front.
The approach used for FLUKA [22] makes use of very general statistical properties of the problem. This approach exploits the properties of the cumulants [23] of distributions, and in particular of the cumulants of the distribution of Poisson distributed variables. The approach can account for an arbitrary threshold for the explicit production of secondary electrons ("delta" rays), for arbitrary step-lengths, and for the contribution to the energy loss fluctuations of distant collisions, the latter with a formalism in part inspired by [24], while assuring the exact match of the average restricted stopping power. The effect of the Mott correction on energy loss fluctuations must also be accounted for both the spectra of the explicitly generated secondary electrons (see eq. (2.11), (2.12)), and the higher moments of the close collision energy losses. The adopted parameterization allows for an analytical calculations of all moments, < >, of energy losses due to close collisions. They are then used in FLUKA in order to reproduce the first 6 moments of the energy loss distribution. Figure 3 shows the full measured Bragg curve and a zoom into the Bragg peak region. A large dose build-up can be noticed in the first 500 mg/cm 2 . This dose build-up is caused by high-energy delta electrons and can also be observed to a lesser extent in Bragg curves of protons [25] and other ions [26]. After this initial build-up, the dose decreases due to nuclear fragmentation of the primary 238 U ions into lighter ions with smaller energy loss until the slowing down of the primary ions creates the Bragg peak shortly before 6 g/cm 2 . After the Bragg peak, a dose extension due to the lighter fragments with ranges higher than that of the primary ions can be observed. This fragment tail is also a well known characteristic of Bragg curves of carbon ions as used for radiotherapy [1] and is apparently very pronounced for 238 U ions.  The horizontal error bars represent the uncertainty of the absolute areal density including the uncertainty of the offset materials and the uncertainty of the individual plates and foils. The uncertainty in relative ionization is harder to estimate, but is believed to be small. At least two measurement points per thickness are shown, however, they can almost not be distinguished by eye. Due to the excellent signal to noise ratio of the ICs (which have originally been developed as beam monitors for ion beam therapy [16]) and the high accuracy electrometers, the fluctuations of repeated measurements were < 0.1%. Therefore the statistical uncertainties of the measurement can be considered negligible. However, possible systematic factors are energy dependencies of the stopping power ratio of polyethylene to ArCO 2 or the w-value of the detector gas (generated electron-ion pairs per dose). Those effects should affect the measurement not more than a few percent [27][28][29] and -8 -would be most relevant in the Bragg peak region and in the distal falloff where the primary particles have low energies. The vertical error bars represent these systematic uncertainties of the measured relative ionization. An uncertainty of ±2% was estimated for the entrance region and fragment tail and larger uncertainties of ±3% for the Bragg peak and ±5% for the distal falloff.

Bragg curve measurement
The measured 80 range of 800 MeV/u 238 U ions in polyethylene obtained from the measurement is 5.9525 ± 0.018 g/cm 2 without the offset materials (0.18 ± 0.009 g/cm 2 , see above) and 6.1325 ± 0.027 g/cm 2 including the offset materials.
Golubev et al. [30] reported a comparable data set: Bragg curves for high-energy 238 U ions stopping in copper and steel targets, also measured at GSI using a calorimeter setup, however with much lower resolution than the Bragg curve in polyethylene presented in this work.

Adjustment of mean ionization potential
The only free parameter for the stopping power prediction is normally the mean ionization potential (I-value). In the ICRU report 37 [31], an I-value of 57.4 eV was suggested for polyethylene and is also the standard value used by FLUKA. However, this value has no major significance because it is simply the average of two available measurements at the time of the report (52.5 eV [32] and 62.2 eV [33]). Therefore, in this work, the I-value was adjusted manually to match the measured range. Figure 4 shows the relative range difference between simulation and experiment (( sim 80 −   The simulated and measured range fit for an I-value of 66.5 eV, which is the value used for the following simulations. The dotted lines represent the uncertainty of the determined range. From this plot, the uncertainty of the adjusted I-value is estimated to be approximately ±2.5 eV under the -9 -assumption that all other stopping power corrections for 238 U ions in FLUKA are correct. However, another important influencing factor of calculated heavy ion ranges beside the I value is the evolution of the effective projectile charge during passage through the absorber and therefore lighter ions or protons would be better suited for a precise measurement of the I-value for polyethylene.

Impact of Bloch and Mott corrections
In figure 5 the experimental Bragg curve is compared to FLUKA simulations using the I-value of 66.5 eV, on the left panel the full Bragg curve and on the right panel a zoom into the Bragg peak region. The solid line shows the result of the standard FLUKA version including all stopping power corrections. If the Bloch correction is switched off (dotted line), the stopping power prediction is too high and consequently the calculated range is 4.2% too short. If the Mott correction is switched off (dashed line) the impact is even higher: the stopping power prediction is far too low and results in an over-estimation of the range by 13.3%. Another interesting detail about the Bragg curve calculated without Mott correction is that the delta electron build-up in the first 500 mg/cm 2 is largely under-estimated. Therefore, it can be stated that Mott scattering is the mechanism that is mainly responsible for the production of the high-energy delta electrons that create the dose build-up in the 238 U Bragg curve presented in this work. The shape of the calculated and measured Bragg curves matches reasonably well. However, differences can be observed in the initial attenuation and the Bragg peak height. The FLUKA simulations under-estimate the initial attenuation and over-estimate the Bragg peak dose. The fragment tail is slightly under-estimated by FLUKA but its shape is very well reproduced. Those discrepancies can possibly be explained by a too low total reaction cross section for 238 U on polyethylene (CH 2 ) in the FLUKA physics model because it would predict too many primary 238 U -10 -ions reaching the Bragg peak depth. Similar effects were recently observed also for 4 He reaction cross sections and Bragg curves [34][35][36]. The FLUKA models predict a total reaction cross section of ∼ 3400 mb for 800 MeV/u 238 U on 12 C targets. It is possible that an increase of this cross section in the FLUKA models would improve the agreement of the simulated curve with the measured Bragg peak shape. The cross section for electromagnetic dissociation of the 238 U projectiles predicted by FLUKA, on the other hand, is only 60 mb and therefore is expected to have only minor impact on the Bragg peak shape. Based on the current data, an optimization of the nuclear models is difficult because experimental cross sections for uranium projectiles are scarce [37].

Conclusion
The measurement of a full Bragg curve for 800 MeV/u 238 U ions stopping in polyethylene is presented in this work. The measured Bragg curve was compared with FLUKA Monte Carlo simulations. The FLUKA code was able to reproduce the measured range after fine-tuning of the I-value. The impact of the Bloch and Mott corrections in the stopping power model implemented in FLUKA was studied by switching them off intentionally and the relevance of these two corrections for correct prediction of the ranges of highly charged heavy ions like 238 U was demonstrated. The shape of the Bragg curve was reproduced reasonably well by the FLUKA simulations, with slight deviations in the Bragg peak height. Further refinements in the total reaction cross section model for 238 U projectiles on polyethylene targets might resolve those discrepancies. The presented data set can be used for validation of any heavy ion transport code for the case of very heavy ions at high energies.
Supplementary data. The data points of the measured Bragg curve are provided as supplementary data attached to this paper.