Mean excitation energy determination for Monte Carlo simulations of boron carbide as degrader material for proton therapy.

Boron carbide is a material proposed as an alternative to graphite for use as an energy degrader in proton therapy facilities, and is favoured due to its mechanical robustness and promise to give lower lateral scattering for a given energy loss. However, the mean excitation energy of boron carbide has not yet been directly measured. Here we present a simple method to determine the mean excitation energy by comparison with the relative stopping power in a water phantom, and from a comparison between experimental data and simulations we derive a value for it of 83.1 ± 2.8 eV suitable for use in Monte-Carlo simulation. This is consistent with the existing ICRU estimate (84.7 eV with 10-15% uncertainty) that is based on indirect Bragg additivity calculation, but it has a substantially smaller uncertainty. The method described can be readily applied to predict the ionisation loss of other boron carbide materials in which the atomic constituent ratio may vary, and allows this material to be reliably used as an alternative to graphite, diamond or beryllium.


Introduction
Proton therapy is a form of radiation therapy whose benefit in particular is the fact that protons release most of their energy at a specific depth (the so-called Bragg peak) to ensure maximum conformity of the deposited dose distribution to the tumour volume. State-of-the-art facilities use the so-called pencil beam scanning (PBS) technique [1][2][3], in which a pristine pencil beam of a defined energy is directed at a specific location of the target volume using steering magnets. For cyclotron-based facilities (which correspond to more than two-thirds of the proton therapy facilities worldwide, as of February 2020), the beam exits the accelerator with a fixed high energy, typically between 230 MeV and 250 MeV. To adjust the proton range to that required for a given treatment layer depth, the proton beam passes through an energy modulation system, which consists of material placed on the beam path either immediately before the patient or upstream of that soon after the cyclotron extraction. For upstream energy variations, such a system is called a degrader. Proton energies down to 70 MeV are routinely used in clinical treatments; the required degrader thicknesses for such an energy reduction are several centimetres, and consequently induce significant lateral scattering and momentum spread, which is then reduced using sets of collimators and a so-called energy selection system (ESS) [4].
Additionally, inelastic scattering and nuclear interaction processes cause a further reduction of the flux of particles, for example converting some of the incident protons into neutrons. The whole process introduces substantial beam losses, and the transmitted beam to the patient at lower energies is often less than 1% of the extracted current at the cyclotron. The losses lead to significant requirements when shielding the cyclotron and ESS, may cause activation of material along those parts of the beamline and radiation damage to equipment, and pose strong limits in the possible current delivered to the patient; treatment times and consequently overall costs will be affected. As the interest in the field of radiation therapy is shifting towards using much higher dose rates for so-called FLASH irradiation [5,6], the ability to deliver larger beam currents to the patient may become a priority in the not too distant future. For this purpose, the choice of the degrader material is of the highest importance. Different materials result in different beam sizes and scattering angles, and especially materials with a low atomic number (Z) may offer greater transmission for a given exit energy from the degrader, because of the smaller lateral scattering for a given ionisation energy loss [7,8]. Graphite is the most common degrader material as it combines a reasonably-low atomic number Z = 6 with good mechanical robustness, economical cost and reasonable density. Beryllium (Z = 4) has been used but is unpopular due to the potential for toxicity when manufacturing and handling degraders. Boron carbide (B 4 C, Z = 5, 6) is mechanically similar enough to graphite but its lower average Z is advantageous at lowering lateral scattering. In addition, several innovative designs have been proposed with the goal of limiting the beam size growth due to the degrader [9,10].
Standard beam dynamics codes that simulate particle transport (such as TRANSPORT/TURTLE [11] or OPAL [12]) are not aimed at precise estimation of beam losses and optimisation of degrader materials considering particle-matter interactions [13]. Usually such codes lack precise scattering cross sections (particularly inelastic ones). Monte-Carlo particle interaction codessuch as those based on GEANT4 [14][15][16] can provide more precise predictions of degrader physics and other relevant phenomena when studying degraders [17,8]; however, they rely on the precise definition of the different characteristics of the materials under study, which are not always well-known and therefore can introduce significant uncertainties in any predictions [18]. Earlier studies of graphite and beryllium as degrader materials [17] have highlighted the importance of the mean excitation energy I Z in accurate predictions of the ionisation slowing, and matching of simulations with experimental data by tuning I Z is not always consistent with reference data such as those reported by the International Commission on Radiation Units and Measurements (ICRU) [19,20]. An incorrect use of I Z directly leads to incorrect definitions of other quantities such as the material thickness, that in turn influences the predicted lateral scattering; a precise determination of I Z matched to the Monte-Carlo model used is required to obtain reliable calculations of other simulated phenomena. On the other hand, a direct measurement of I Z is very challenging, as it requires precise knowledge of several correction factors, as recently highlighted in Seltzer et al. [20].
We are currently investigating the potential of B 4 C as a degrader material to be used in the Paul Scherrer Institute Centre for Proton Therapy (PSI) [21]. B 4 C is a material in which the B to C atomic ratio can vary significantly [22], and it is important to carefully quantify that ratio, its homogeneity and the mean density in any given sample used. In the context of Monte-Carlo degrader studies, we are proposing a method to determine the value of I Z from experimental measurements of a given degrader or other material sample. With our study, we aim at answering the following questions: • How can we determine the I Z parameter for a material, keeping the uncertainties under control? • Can such a method be reliably and easily applied in a clinical facility, to help different institutes setting up their own simulation environments?
In this paper, we demonstrate a simple method that for the first time derives an estimate of I Z for B 4 C from a comparison with experimental data. Such a value is suitable to use for Monte-Carlo predictions and design. The method described belowbased on a relative measurement versus a water phantommay be straightforwardly used for other material samples, and we have applied it here also to conventional graphite degrader materials.

B 4 C as a degrader material in the PSI proton therapy system
An optimal combination of low atomic number (to minimise multiple scattering angle) and sufficiently high density (to minimise the amount of material needed for the energy degradation) is the goal when choosing degrader material [7]. The present degrader material used at PSI is graphite, with density ρ = 1.870 ± 0.005 g cm − 3 , and is used to degrade the energy for treatments in the range 70-230 MeV. But due to the losses introduced in the degrader and ESS sections, in the PSI proton therapy beamline the overall transmission efficiency at the lowest beam energy of 70 MeVfrom the cyclotron extraction to the treatment nozzle -is only around 0.2%, comparable to other similar facilities. We are also investigating an alternative candidate materialboron carbide (B 4 C); the lower Z together with the higher density (ρ = 2.52 g cm − 3 ) should lead to lower emittance growth after the degrader in comparison to graphite, and thereby could enable treatments with proton energies below 70 MeV at an acceptable dose rate. B 4 C is the third-hardest material known (after diamond and cubic boron nitride), and as such is used in such applications as body and vehicle armour and to coat cutting tools. It can also absorb neutrons without forming long-lived radionuclides, and is widely used in nuclear reactor control rods, shielding materials and for neutron detection [23]; in particle physics it has also been investigated for use in muon experiments [24]. A first experiment with B 4 C at PSI showed promising results, indicating a transmission increase of more than 30% using B 4 C with respect to graphite, although with a simplified degrader geometry and not at the lowest 70 MeV energy [21]. A follow-up study is in progress that examines a wider range of energies and uses the PSI double-wedge degrader geometry.
The mean excitation energy of B 4 C is presently known only from an estimate by Seltzer and Berger [19], who use a 13% correction to the value predicted by the Bragg additivity rule; I Z was estimated in that work to be 84.7 eV. This estimate and the correction is however based upon experimental values from other compoundsand particularly compounds not even containing boronand has never been validated experimentally. Uncertainties on this value of I Z are conservatively estimated to be somewhere from 10% to 15%; this estimate is reported in the ICRU tables and in several simulation codes (including GEANT4). Variation of the B:C ratio from 3.5:1 to 10.5:1 gives a variation of the additivity-derived I Z value of 0.7 eV.
Several blocks of B 4 C were manufactured using hot powder pressing by an industrial supplier 1 ; the manufacturer estimate of B:C atomic ratio was 78.3%:21.7%, slightly less than the nominal 4:1 ratio. Porosity was estimated by the supplier to be below 0.5% by volume; we independently assessed the quality of the samples (manufacturing defects, uniformity and density variations) by carrying out computed tomography (CT) scans. The observed Hounsfield unit variation in the CT scans were at the level of the scanner noiseindicating good homogeneityand there were no observed voids or other significant defects. We used these samples in tests in the PSI Gantry 2 treatment room, to estimate the mean excitation energy of the material from range measurements. The experimental setup for the measurements is explained in Section 2.2, whilst the method used for the estimation of I Z is explained in Section 2.4.

Range scanner measurements at PSI Gantry 2
We used an in-house built water phantomthe so-called Range Scanner (RS) shown in Fig. 1 to measure the Bragg peak profile as a function of depth in water. The scanner box is made of Plexiglass and allows continuous movement of an ionisation chamber through the contained water over a range of 450 mm and with a resolution of ±0.025 mm, using a motor controlled by an in-house-programmed LabView interface. Several chamber holders may be mounted inside the scanner that allow for measurements using a variety of ionisation chambers; in our experiment we mounted as measurement chamber a large, plane-parallel ionisation chamber (80 mm sensitive area diameter), the same as used for the range commissioning and quality assurance (QA) tests for treatments at our centre.
The LabView interface allows the user to define the step resolution for the measurement, allowing the combination of regions of coarser resolution (for example, in the buildup region) with others of finer resolution (for example, around the Bragg peak); this gives flexibility and efficiency in taking the data. The normalization of the dose measured along the peak is done using a reference chamber. At PSI Gantry 2 we use the dose monitor in the gantry nozzle as a reference chamber; in the LabView interface it is possible to define a preset number of monitor units (MU) for the reference chamber, and the measurement system will accumulate the integrated charge (IC) collected in the measurement chamber while the reference chamber reaches that preset. The charge collected as a function of position of the measurement chamber in the RS is used to produce the Bragg peak profile. The fall-off of the curve is then fitted with a cubic spline interpolation and the range is extracted as the 80% point from the peak height (to make this better visible to the reader, we normalised all plots in the paper to the peak height). In our present measurement, we used the same settings for preset and step resolution as in our yearly QA tests, to ensure good quality of our measurements.
The measurements were performed in the Gantry 2 treatment room; therefore the beam characteristics are the same as described in Pedroni et al. [25], with typical spot sizes (σ) at the isocentre of around 2.3 mm at 230 MeV. The nominal beam energy is defined by range measurements performed with the RS and matched to the ICRU 49 tables [26]. The energy spread (around 0.7%) is estimated from the energy selection system acceptance and from matching to measured Bragg peak profiles.

GATE/GEANT4 simulations
GEANT4 is a toolkit for Monte-Carlo simulation of particle interaction with matter. Its use case mainly originated in high-energy physics, but its models have been expanded to cover a broad range of particle types and energies, including applications in space radiation, cosmic ray modelling, nuclear physics, heavy-ion and radiation physics, and medical applications. For the simulations described in this paper we used the GATE application for GEANT4 [27,28]. Although initially conceived for imaging applications, GATE has been developed to allow radiotherapy applications, including proton and hadron therapy [29]. We use the term 'GEANT4' for brevity, meaning our simulations developed using GATE and GEANT4; we used GATE version 8.1.p01 and GEANT4 version 10.04 patch 2.
Results obtained using GEANT4 simulations are dependent on the physical models employed and by tuning parameters such as range cuts and energy limits. We followed the recommendations of the GATE collaborationsummarised in Grevillot et al. [29] and updated in Fuchs et al. [30] and Resch et al. [31] for scattering validationto fix the parameters and the physics processes used. GEANT4 provides different electromagnetic transport methods with varying precision; we used G4EmStandardPhysics_option4 (also called EMZ) as it provides the greatest precision [32] for the low-energy regime. As noted for example by Sabin et al. [33], a particular mean excitation energy I Z must be matched to the correct physics model in GEANT4 to give a reliable energy loss estimate. We used the settings of Grevillot et al. [29] for the Bethe-Bloch energy loss, which can be used down to 2 MeV; below that a parameterisation based on ICRU 49 is implemented. This list uses the highest precision low-energy models (G4KleinNishinaModel, G4Low-EPComptonModel, G4PenelopeGammaConversion and G4GoudsmitSaun-dersonModel). For proton multiple scattering simulations, it has been shown that the models available in the GEANT4 physics lists show different degrees of agreement with experimental data; Fuchs et al. [30] have shown that the WentzelVI model for multiple Coulomb scattering (MCS) results in a better agreement with data published in the literature; the WentzelVI model has been the default MCS model used in the EMZ option since version 10.02 of GEANT4, and is used also by us in the present work. For nuclear models, we also followed the recent recommendation of Resch et al. [31] and used the QGSP_BIC_HP physics list for our simulations; this list contains elastic scattering cross sections that have been shown to provide a better agreement with experimental data than other models (for example QBBC) for the energy range of interest of our simulations (between 60 MeV and 250 MeV). However, the difference between the two physics lists is mainly in the predicted transverse dose profiles, andsince we are integrating the dose over a large area transverse to the beam directionwe don't expect a large impact from the choice of the elastic scattering model. The additional option 'HP' used with this physics list ensures that high-precision neutron models and cross sections are used in the simulation.
GEANT4 simulations can also be affected by several other parameters; in particular, the production threshold for secondary particles (electrons/positrons and gammas) can be defined as a distance (or range cut-off) that is internally converted to an energy for individual materials. Production thresholds are defined for a given geometrical region. We fixed the production threshold at 1 mm for the region surrounding the samples and the water phantom, and to a lower value for the target materials and the water phantom where a better precision is needed. The impact of this lower production threshold was verified using simulations without degrader samples, by checking that the 80% proton range in water does not change as the threshold is lowered. A lower production threshold of course increases computing time, and therefore a higher production threshold is to be preferred for computational efficiency. We observed a 0.02 mm change in predicted proton range when varying the production threshold from 0.01 to 0.001 mm, for an 80 MeV simulated beam onto a water phantom; reducing the production threshold comes at a great additional expense in computing time (from 128 primaries/s for 0.01 mm production cut to 1.8 primaries/s for 0.001 mm production cut). We chose an intermediate 0.005 mm production cut for our simulations, to minimise the impact of the production cut on the resulting range but at the same time to keep the computing time reasonably lowaround 10 primaries/s simulated. We could not find a significant variation in the range from changing other cuts (e.g. tracking cuts and step limiters), and therefore for those cuts we used the standard GEANT4 parameters.

Determination of mean excitation energy
In our measurement we are attempting to estimate indirectly the I Z of a material from a comparison between simulations and experimental measurements. It is well known that the Bragg curve shape, and most of all the position of the peak, depends on the I Z of the material the particles are stopped in [34,18,35]. However, the variation with I Z is quite subtle, and can easily be masked by uncertainties both in the measurements as well as in the simulation parameters [18]. Therefore we decided to base our comparison of experimental data to simulation not on the single Bragg curves, but on the difference in range (ΔR) between the Bragg curve measured with a given sample and the one without the sample. In this way, we expect many sources of uncertainty to cancel out or be lowered, and in particular the incident proton energy does not need to be known precisely, in contrast to previously-described methods. We have described in 2.2 how we can extract the range from the RS measurements. For the simulations, we simulated the same setup (water phantom with or without sample at the entrance) and scored energy deposited in a fine grid of 0.01 mm along the water phantom depth, while integrating it on a large area to fully collect the energy deposited by large-angle scattered particles. We then interpolated the obtained curve on the falloff region and extracted the range as the 80% from the maximum.

Measurement uncertainties
Range measurements with the RS have a good resolution but are affected by several systematic uncertainties, particularly the estimation of the thickness of the water tank window, the thickness of the waterproof housing of the chamber, the position of the chamber reference point, and the calibration of the RS motor that gives the absolute position. Using ΔR measurements instead of absolute range measurements lowers the uncertainties since all thicknesses and calibration factors are common to both measurements and act on a measurement as a global shift of the range, therefore cancelling out when calculating ΔR. To limit other experimental variations arising from such things as the daily setup and positioning of the beam, we performed both range measurements (i. e. with and without the degrader sample) on the same day and with the same beam setup.
Density variations of the material samples will result in an uncertainty on the measured range which does not cancel when calculating ΔR. The density of our B 4 C was not well known from the data we obtained from the supplier, due to the specific method used to produce the samples. We determined the average density of our B 4 C samples directly by weighing each sample and measuring its dimensions with calipers; we measured an average density of 2.52 g cm − 3 from all samples, the precision of which measurement is limited by the caliper measurement of the sample thickness. We observed a maximum 1% mass variation across all samples, which could be explained either by density variation or by volume variation (due to manufacturing tolerances). CT scans of the degrader material confirmed the homogeneity of each sample but did not allow (from the Hounsfield units) a more accurate determination of the density than that obtained from direct size and mass measurements. Due to irregularities in the shapes (presence of grooves etc.), checking the consistency of the volume between all samples proved difficult. So, to estimate conservatively the impact of possible density variations on the measurements, we compared ranges measured with the beam impacting on a different position on a sample and by using two samples, and considered the resulting difference as a density uncertainty.
The Bragg curves and resultant proton range predicted from simulations depend on a variety of parameters [18], and uncertainties in variables or changes to the physical models used to determine the uncertainty of the range. A variation of the used value of the I Z for water, or the choice of physics list and production cuts, gives a systematic shift in the predicted range. However, when calculating ΔR these systematic shifts are expected to cancel out. The same is true when considering the uncertainty on the incident proton energy; we consider this briefly in Section 3.1, but again the effect of (small) variations of the incident energy upon ΔR is also small. In contrast, uncertainties such as the precision with which the thickness of the degrader sample is known, the density variation within that sample, and statistical uncertainties in the Monte-Carlo simulation itself do not cancel in ΔR and therefore contribute to the final uncertainty.
To validate our method and to investigate its limitations we applied it also to a graphite sample, here a spare degrader as used in the current PSI degrader; the results are reported in Section 3.4.

Benchmarking of GATE/GEANT4 simulations of the experimental setup
The experimental setup corresponds to the setup used during commissioning and yearly QA measurements of the clinical beam line of PSI Gantry 2; therefore, we benchmarked our simulations against the last yearly QA measurements, performed only three months before our experiment. Fig. 2 shows a comparison between four simulations of the PSI Gantry 2 beam in the water phantom and how they compare with the measured Bragg curves; we observe a good match in the peak width when using the known energy spread of 0.7%, whilst the range predicted by the simulation is slightly overestimated with respect to the data; the mismatch is corrected by using a mean beam energy in simulation that is 0.35 MeV lower than the nominal beam energy expected at the entrance to the experiment (note that proton energy loss in the intervening air is already taken account of). We repeated the Bragg curve measurements for a 230 MeV beam incident on the same day of the experimentwithout the degrader samplesand found a excellent match between those two curves; the measured difference in the fitted 80% range between data and simulations is 0.02 mm, similar to the RS resolution. This shows that variation in measured range is small. We investigated whether the range mismatch with the simulations arises from the scattering model used, and/or from the chosen value of the I Z for water. As shown in Fig. 3a, the scattering model does not impact the range measured in our experimental conditions. However, the choice of the mean excitation energy of water used in the simulation does. The nominal beam energy is defined by matching the measured ranges with the ICRU 49 report [26], which uses an I Z of 75 eV for water; GEANT4 uses 78 eV, so the discrepancy in range is due to this difference in I Z , as shown in Fig. 3b. However, since the method we use to determine the I Z of the B 4 C sample is based on the differential measurement ΔR, the additional uncertainty caused by the range mismatch will cancel out almost completely. For this reason, we decided to match the Bragg curves between measurements and simulations by adjusting the energy of the input beam in the GEANT4 simulation. Fig. 4 shows an example of the measured curves with the range scannerboth with and without the B 4 C degrader sampledemonstrating the shift in the Bragg peak which can be observed when the B 4 C sample is placed in the beam; measurements with either of our B 4 C samples gave almost the same range, and the difference was used in the density uncertainty estimate. We optimised the number of points collected along the Bragg curve for efficiency, using a coarser resolution in the plateau region and a finer resolution around the peak. Each point corresponds to a fixed number of monitor units (100, 000, about 9.2 × 10 8 protons), from which the dose is obtained using the dose monitor in the PSI Gantry 2 nozzle as a normalization signal. The measured curves are smooth, except for a few points where instabilities in the delivered current from the accelerator occurred (for example, in the plateau of the B 4 C data in Fig. 4). These instabilities however occurred only in a few points collected on the plateau regions of the measured curves; therefore they do not affect the analysis of the falloff of the measured curves. Fig. 5 shows how the variation of the I Z of B 4 C affects the position of the Bragg peak in the water phantom. We carried out simulations with the B 4 C sample for fourteen different values of I Z , using a finer resolution of 1 eV between 80 and 90 eV and a coarser resolution of 10 eV between 70-80 and 90-110 eV; a smaller range was used for the comparative test with graphite since we have a fairly precise reference value to compare to. We fitted the full distal edge of the curves with a cubic spline interpolation, and extracted the range as the 80% point from the fitted curve. The same fit was performed with the reference simulation that had no material sample. Fig. 6 shows the ΔR variations as function of I Z in the simulation, overlaid onto the measured value; the bands shown indicate the experimental and simulation uncertainties. As mentioned in Section 2.4, using ΔR removes most sources of systematic uncertainties. For the experimental data, the uncertainties left are the precision of the water scanner measurement, which is ±0.025 mm, and the density uncertainty, which is ±0.1 mm; when added in quadrature, the density  In simulations, we had to account for thickness and density variations of the sample, for the statistical uncertainty in the simulations, as well as for the uncertainty in the range extracted from the fit of the falloff. To estimate the uncertainties due to thickness and density variation, we repeated the simulations for a representative I Z value with a 1% variation in the density and a 0.002 cm variation of the thickness, and obtained the resulting change in range. As mentioned in Section 2.4, from mass measurements we observe a maximum 1% mass variation across all samples; since we could not estimate precisely volume variations for the different samples, we have directly translated this variation into a 1% density uncertainty. Combining the thickness and density uncertainties gives an overall uncertainty in the WER of ±0.1 mm, while the interpolation and the statistical uncertaintywhen obtaining the range from a fit to the Bragg peakgives an uncertainty in WER of ± 0.05 mm; since they are not correlated, we sum them in quadrature to give the overall uncertainty in WER.

Determination of the mean excitation energy of B 4 C
To extract I Z from a comparison between data and simulations, we then performed a polynomial fit of degree 3 through simulated ranges. The uncertainty in the estimated range is obtained as the 68% (1σ) confidence interval in the polynomial fit. By comparing the measurements to the simulations, we obtained a value for the I Z of our B 4 C samples to be 83.1 ± 2.8 eV.
To test the consistency of our result, we also estimated ΔR for an incident beam of 150 MeV and 83 eV mean excitation energy (it is worth reminding that ΔR is independent of the incident beam energy). The result was compatible with the one obtained in the 230 MeV simulation within uncertainties, further validating our simulation setup. Fig. 7 shows the application of the same method described above to a graphite sample, here a spare degrader for the PSI PROSCAN facility with a measured bulk density of 1.870 ± 0.005 g cm − 3 . The figure shows the ΔR changes as I Z is varied in simulation. We used again a 1% uncertainty in the density of the graphite as for B 4 C. We obtain an I Z for our graphite sample of 81 ± 5 eV. The relative error obtained is larger because of the larger uncertainties in the fit (lower number of points) as well as the lower density of graphite; this gives a smaller variation of ΔR with I Z with respect to B 4 C. The value obtained for I Z is in good agreement with the accepted reference value: the latest ICRU estimate gives 81.0 ± 1.8 eV [20]), thus supporting the validity of our method.

Discussion
We have developed a simple method to determine the I Z value of a given material sample such that the I Z may be reliably used in Monte-Carlo simulations; this method overcomes the limitations of simulations that must today employ the Bragg additivity rule to estimate I Z . Using our method, we estimate the I Z of B 4 C to be 83.1 ± 2.8 eV; this uncertainty is larger than the uncertainty arising from variation B:C atomic ratio. Our method may readily be applied in any clinical proton treatment facility, since it utilises detectors that are already commonly used for beam commissioning and QA, and it can be modelled using a simple simulation setup.
Our measurement is the first experimental determination of the I Z of B 4 C which was previously based only on the calculations of Seltzer and Berger [19]. The previous B 4 C estimate of I Z = 84.7 eV was determined using the Bragg additivity rule and a 13% correction for compounds, and is the one commonly found in well-known databases such as ESTAR and the GEANT4 material database; the uncertainty on this value is poorly known, but estimated to be 10-15%. Our present result agrees with the reference value, but with a much smaller uncertainty of 3%. Our study indirectly confirms that the corrections proposed by Seltzer and Berger [19] for compounds indeed apply also to B 4 C although more precise experimental measurements with varying B:C atom ratios would be needed to investigate this point in detail.
It is important to improve in particular the determination of I Z for B 4 C since it is has been proposed by different groups [21,36] as a promising alternative material to graphite as a proton therapy degrader but whose slowing properties were hitherto not well known. We point out here that the atomic compositions of a given B 4 C sample may vary from those reported here such that the I Z of a given degrader sample may be different to the value we obtain; it is therefore important to test samples using the method we describe.
We also tested our method on graphite, a material whose characteristics are much better known than B 4 C. Albeit with a larger relative uncertainty, our method provides an estimate of I Z compatible with established reference values. The lower density of graphite means that ΔR has a smaller dependence upon I Z and therefore yields a poorer estimate of I Z for simulations. However, this will probably not be a very strong limitation when using simulations to design degraders, since both the small dependency of the range on I Z as well as the large degrader thickness needed (particularly for low energies) will probably together result in a small relative uncertainty on the thickness. For materials with a higher density and a larger dependency of the range on I Z the impact of these uncertainties on the predicted thickness is larger; therefore having a precise estimate is important, particularly when estimating quantities that strongly depend on a precise knowledge of the degrader thickness, such as emittance or transmission.
Previous studies of graphite degraders performed at PSI by van Goethem et al. [17] found that I Z = 95 eV; this resultobtained using a magnetic spectrometer to measure the degraded beam energydiffers significantly from the accepted ICRU value of 81 eV, and the difference was attributed by them to different models for the correction terms used in GEANT4 for Bethe-Bloch compared the ICRU calculations. It is common to 'tune' the I Z used in a particular simulation to obtain agreement with a measured proton range; however, such tuning can mask other discrepancies such as an incorrect density, inadequate measurement or simulation method, and can thereby lead to incorrect predictions of other phenomena such as scattering or beam current transmission. A limitation in Van Goethem's measurements may be the reliance on a knowledge of the incident proton energy, which can give systematic shifts in the simulated range that yield a different I Z . In the present work we obtain an I Z consistent with ICRU values [20] using a method that does not depend on the exact incident proton energy, and in which we have accounted carefully for the actual density, dimensions and composition of the particular degrader; this I Z value may be also used in simulationusing GEANT4 high precision electromagnetic models as well as the WentzelVI modeland therefore we conclude that this simulation method may be used reliably for predicting energy reductions from degraders.
We comment that the choice of MCS model does not significantly affect the predicted proton range; we compared in our simulations the Urban model used by Van Goethem with our preferred WentzelVI model, and we see little difference in the dose deposition. However, additional measurements of other phenomena such as scattering distributions would be needed to fully benchmark a model, and we have only reported a small part of that benchmarking (as discussed in Section 3.1) for the sake of brevity. Greater consideration of scattering models may be important if this method is applied to experiments using smaller ionisation chambers for the Bragg peak measurements, as highlighted in Grevillot et al. [37].
With Monte-Carlo simulation codes becoming more precise, faster and easier to set up, their use for proton therapy system design will only increase; this is already shown in the context of magnet design [38] and degrader design [8,17,21,13]. Our present work contributes to this by more carefully defining material properties for reliable simulation than has been previously described. With the present interest in high dose rates in the context of FLASH irradiation [5,39], optimisation of beamline design to achieve the best beam transmission will substantially increase in importance. The method reported here to obtain an improved determination in the I Z of B 4 C contributes to this, by allowing candidate degrader materials and configurations to be examined more reliably in proton therapy system design; this will be a first and essential step in designing future innovative degrader systems.