Ambient dose simulation of the ProtherWal proton therapy centre radioactive shielding decay using BDSIM and FISPACT-II

. Next-generation proton therapy centres couple treatment and research programs, leading to higher beam currents and longer irradiation times than in clinical conditions. Large ﬂuxes of energetic secondary particles are produced and long-and short-term radioactive nuclides are generated in the concrete shielding of the cyclotron vault. While the overall long-term activation of the centre is well known from the shielding design activation studies, the short-term activation peaks are still of importance when radiation protection studies are involved. The centre shielding design was validated using the BDSIM/FISPACT-II methodology combining particle tracking and Monte Carlo particle-matter interactions simulations using Beam Delivery Simulation (BDSIM) and the computation of the activation using FISPACT-II. We establish, as the next stage of our methodology, the simulation of the decay radiation of the activated concrete shielding and the accurate scoring of the related radiation protection quantities. A single BDSIM simulation per radioactive nuclide is performed based on the nuclide concentration obtained from the prior FISPACT-II activation computations at the start of a given cooling period. The evolution of the radiation protection quantities is obtained by scaling the results with the nuclides activity obtained at later times from fast FISPACT-II computations. We show the evolution of the ambient dose equivalent in the centre vault when considering regular concrete and Low Activation Concrete (LAC) as shielding material to demonstrate the eﬃciency of LAC mix in mitigating the shielding activation.


Introduction
The shielding design studies of proton therapy systems are typically performed using general-purpose Monte Carlo codes such as MCNPX [1], Geant4 [2,3], FLUKA [4] or PHITS [5]. They use discrete neutron or proton sources in radiation transport simulations with limited modelling of the magnetic fields in the transport beamline. The source terms of these radiation transport simulations are obtained from primary beam tracking simulations using codes such as MAD-X [6] or Transport [7] from which beam loss patterns are extracted. The Monte Carlo radiation transport simulations track primary and secondary particle types and compute quantities required for shielding design: the secondary particle fluence, shielding activation, residual decay or ambient dose equivalent. However, recently designed proton therapy systems tend to be more compact, single-room systems combining clinical and research activities. This is the case of the * e-mail: eliott.ramoisiaux@ulb.be single-room Ion Beam Applications (IBA) Proteus One (P1) "ProtherWal" centre in Belgium, featuring a complete clinical, biological, material-science and beam delivery research program leading to higher currents, longer irradiations and lower degraded energy. The production of secondary particles is therefore expected to be much higher than in non-research-oriented facilities. The shielding walls activity of the ProtherWal centre is expected to increase by a factor of three. This imposes additional constraints on the predictive power of the radiation protection Monte Carlo simulation results. Our methodology using Beam Delivery SIMulation (BDSIM) [8] and FISPACT-II [9] fulfils these requirements and has been validated [10]. It challenges the paradigm of using discrete source terms and combines magnetic tracking of the primary proton beam and the particle-matter interactions within a single simulation.
BDSIM is a Geant4-based Monte Carlo tracking code that allows one model to simulate primary and secondary particle tracking in any beamline and surrounding geometries while also considering the particle-matter interactions. BDSIM Monte Carlo models with integrated beam tracking consider all losses and continuously simulate the secondary particle generation along the beamline, which inevitably leads to results with higher accuracy. Furthermore, the polyvalence of these models allows them to have predictive roles and drive an essential part of the research studies. Indeed, new beam optics, energy degradation techniques or materials that impact the losses and, subsequently, the secondary particle production of the system can be optimised to decrease the beamline and shielding activation and lower the residual decay distribution; in addition to allowing direct beamline optimisations.
FISPACT-II solves the rate equations using ENDFcompliant group library data for nuclear reactions, particleinduced or spontaneous fission yields, and radioactive decay [9]. The methodology was first applied to the concrete shielding activation of the to-be-built proton therapy facility in Belgium, which uses the newly developed Low Activation Concrete (LAC). Only the secondary neutrons were considered in the FISPACT-II computations, as the neutron-induced activation in the shielding walls dominates the contributions from other particle types. LAC was introduced to decrease the amount of radioactive waste produced during the centre lifetime. The activation results comparing the efficiency between the regular concrete and the LAC were validated against prior MCNPX shielding design studies [10]. Since then, the BDSIM/FISPACT-II methodology has been applied to improve knowledge of the centre activation [11,12] and to prepare the future experimental measurements campaign [13]. Those results confirm the efficiency of the BDSIM/FISPACT-II methodology to simulate the prompt irradiation and compute the induced activation.
To extend the capabilities of this workflow to the simulation of secondary fluences over arbitrary timescales and irradiation profiles including cooldown times, we implement a new feature in BDSIM to provide simulations of radioactive decays supporting all the capabilities of the Geant4 radioactive decay module [14]. This enables a BDSIM/FISPACT-II workflow where activation results are used as input for decay radiation studies. This allows to compute radiation protection quantities from arbitrary irradiation time profiles. This new feature is first applied to an example with simplified geometry and compared with FLUKA results.
The extended version of the BDSIM/FISPACT-II methodology is depicted in Figure 1. The method consists of several stages. The user generates the beamline and shielding models by combining the default accelerator component geometries provided by BDSIM and more complex geometries that can be imported from Geometry Description Markup Language (GDML) files. 3D meshes are positioned in the BDSIM model for the secondary neutrons differential fluence scoring. BDSIM simulations are run using a given input beam and an appropriate physics list for the particle-matter Monte Carlo simulations. Then, the differential fluences are extracted for each voxel and given to FISPACT-II, which computes the nuclide inventories and the related radiological outputs for given irradiation and cooling time sequences using a given cross-section database. Finally, the mesh nuclides total activity and its bins nuclides atomic concentration are used in BDSIM for the decay radiation simulations and the scoring of any desired radiation protection quantity.
This new workflow is then applied in full to the study of the ProtherWal shielding to compute the ambient dose resulting from the activated shielding. In particular, the efficiency of LAC to mitigate decay radiation fluence is evaluated in detail. We compare the evolution of ambient dose equivalent induced by the decay radiation of the activated shielding for two different configurations: entirely in regular concrete and with a fraction of LAC. Additionally, we examine the evolution of the shielding activation levels for each radioactive nuclide. The results show that LAC significantly decreases the ambient dose equivalent in the centre vault and is an asset for future compact proton therapy centre shielding design. Section 2 discusses the implementation of the radioactive decay feature in BDSIM and presents the new BDSIM/FISPACT-II workflow. The validation against FLUKA simulation results is presented in Section 3. The application to the ProtherWal model is introduced in Section 4, where all the simulation results are discussed in detail. We draw conclusions in Section 5.
2 Implementation of the Geant4 radioactive decay feature in BDSIM The radioactive decay feature of Geant4 [15], which relies on the Evaluated Nuclear Structure Datafile (ENSDF) data library [16] and accounts for the α, β − , β + , electron capture, isomeric transition, neutron emission, and proton emission, has been implemented and integrated within BDSIM as a new physics list labelled "radioactivation". This feature contains built-in variance reduction methods, which can be activated through biasing parameters. The Geant4 biasing parameters are NSplit to set the number of times a nucleus undergoes a decay, BRBias to give all branching ratios the same probability, DecayBias to define a decay time window and SourceTimeBias to define a source profile window. The option AnalogueMC Fig. 2. BDSIM/FISPACT-II radioactive decay workflow. For each event, a random number is generated to select a mesh based on the target nuclide activity and then a bin based on the mesh nuclide concentration distribution. In this example, the radioactive nuclide was chosen to be placed in Bin 1 of Mesh 2. must be activated for analogue Monte Carlo. The user also has the choice to simulate the whole decay chain by setting the option fullChain.
The BDSIM library was developed to propose a solution for simulating accelerator-based systems coupling particle tracking and particle-matter simulations. As BDSIM was developed with primary beam tracking in mind, new options had to be implemented to support the definition of radioactive nuclei at rest as particle sources. The option radioactiveDecay allows the primary beam to be defined as static particles with the particle kinetic energy E0 set to 0.
Punctual sources of radioactive nuclides can now be simulated inside a BDSIM simulation. However, more complex simulations, for example, decay radiation from activated materials, cannot be defined by hand. Such simulations require allowing prior activation study results to be used as sources for the decay simulations. The BDSIM/FISPACT-II methodology for activation study is thoroughly described in Ref. [10]. The methodology is divided into the differential secondary particle fluence scoring and the FISPACT-II activation computation. The information about the scoring meshes positioning in the BDSIM model (translation and rotation), the nuclide total activity in each mesh and the nuclide atomic concentration in each mesh and bins are extracted and written in a ROOT file. The new disrType beam distribution labelled "radioactivedecaysource" reads this file and follows the specific workflow described in Figure 2 for a target nuclide. A first random number is generated for each event to select a scoring mesh based on the target nuclide activity. Then, a second random number is generated to choose a bin in the selected mesh based on the mesh nuclide concentration distribution. Finally, the bin position and dimensions are extracted, and the radioactive nuclide is positioned randomly inside the bin volume to be decayed.
A BDSIM decay simulation, scoring any desired quantity Q, is performed for each target nuclide j to obtain the studied quantity distribution per target nuclide total activity Q Aj within the scored volume. While this part can be time-consuming, once the Q Aj distributions are obtained, fast FISPACT-II computations are used to provide snapshots of each nuclide total activity at any requested moment A j (t) and the evolution of Q is simply Fig. 3. BDSIM and FLUKA models used for comparison. The primary beam is represented in blue. The beryllium target is represented in green. The regular concrete shielding is represented in grey with the floor and the roof made transparent for visualisation purposes. The secondary particle production is not shown for clarity.
given by, The influence of time of the decay-emitted γ spectrum on the value Q is considered through the time-dependent linear combination of the Q Aj values.
This workflow simulates the decay radiation of activated material computed from a prior activation study. Radiation protection quantities related to the decay radiation can now be scored in BDSIM.

Comparison and validation of the BDSIM implementation against FLUKA simulations
The new extension of the BDSIM/FISPACT-II methodology is compared to results obtained with the FLUKA Monte Carlo code [17,18] for validation. FLUKA has been extensively tested and validated for various radiation types and energies and is widely used in the radiation protection and dosimetry communities [19,20]. The model used for the comparison is of simple geometry: a beryllium target is placed inside a regular concrete (see composition in Tab. 3) shielding of 15 cm wall thickness as shown in Figure 3. A 230 MeV proton beam is sent through the target to be degraded to 100 MeV continuously for 20 years with a beam current of 5 nA.
The quantity chosen for the comparison is the ambient dose equivalent, H * (10), which is the operational quantity suggested for area monitoring of strongly penetrating radiation. H * (10) is defined as the dose produced in the International Commission on Radiation Units and Measurements (ICRU) [21] sphere at a depth of 10 mm on the radius opposing the direction of the irradiation field [22], where the ICRU sphere is defined as a 30 cm diameter tissue-equivalent sphere of unit density. BDSIM allows the computation of such operational quantity by making the product between the particle fluence φ (cm −2 ) and fluence-to-ambient dose equivalent conversion coefficients h 10 (Sv × cm 2 ). The resulting expression for the ambient dose equivalent is given by, The coefficients are taken from Ref. [22] and implemented in BDSIM as described in Ref. [23].
The H * (10) values are scored using a 5 cm×5 cm×5 cm scorer mesh for 29 of the most important γ-emitting radioactive nuclides in concrete, as detailed in Table 1. Those nuclides have been chosen based on the Γ power information provided in the FISPACT-II output files and defined as, with A i the activity of the nuclide i and E γ,i the γ decay energy for nuclide i.
It is important to note that the simulation of metastable states is possible and considered in this workflow. The physics lists used during the Monte Carlo simulations and the nuclear data used for the activation computation can impact the results [24]. The physics list and nuclear data configuration used in BDSIM are the same as the one used for the shielding design in Ref. [10]. The configuration used in the FLUKA simulations corresponds to the activation example on the FLUKA website [25]. This paper does not aim to fully investigate which physics list and nuclear data configuration offers similar results to FLUKA. This comparison aims to validate the new workflow that uses prior activation results with the BDSIM/FISPACT-II methodology as inputs for decay radiation simulations and related quantities scoring.
The H * (10) during the prompt irradiation obtained with BDSIM and FLUKA are respectively presented in  Figure 4a,b. Minor variations can be observed and are explained by the difference in physics list between the two Monte Carlo codes. Furthermore, the BDSIM simulation only scores the H * (10) in the air while FLUKA also computes it inside the whole volume leading to some discrepancies in the target volume. However, the correspondence between the results provides confidence in using the shielding activation for the decay-induced H * (10) comparison. The decay-related H * (10) values were computed in BDSIM following the workflow described in Section 2. The comparison is performed 12 hours after the end of the target irradiation. The results obtained with BDSIM and FLUKA are presented in Figure 5a,b. The differences between the methodologies are more pronounced than with prompt irradiation. BDSIM gives larger H * (10) values than FLUKA. However, these discrepancies can be explained by the accumulation of variations at each step of the code respective methods from the difference of physics list during the prompt and decay Monte Carlo simulations and the difference in nuclear data used for activation  In view of the results obtained, the extended BDSIM/FISPACT-II methodology for decay radiation and related quantities scoring is considered validated. Section 4 illustrates this new option with the realistic case of the to-be-built ProtherWal proton therapy centre in Belgium.

Application to the ProtherWal proton therapy centre
The extended BDSIM/FISPACT-II methodology is now applied to the future proton therapy centre of Charleroi, Belgium, to be built in the context of the ProtherWal project. References [23] and [10] describe the BDSIM modelling and experimental validation of the IBA Proteus One system and the ProtherWal complete shielding design using the BDSIM/FISPACT-II methodology. The critical result obtained during the shielding design was the efficiency of the new concrete mix, LAC, to decrease the amount and hazardousness of radioactive materials accumulated at the centre decommissioning. Indeed, LAC is produced with aggregates that allow reaching a concrete with lower concentrations in the trace elements Eu, Co and Cs at the origin of long-lived nuclides than in regular concrete, see Table 2. The LAC also features lower concentrations in Na, Al and Si responsible for the generation of radioactive nuclides through capture and spallation reactions such as 22 Na, 24 Na or 28 Al. Table 3 compares the atomic composition between regular concrete and the LAC. Figure 6 is a Paraview [27] representation of the cyclotron vault. The extraction of the P1 is composed of the 230 MeV SuperConducting Synchro-Cyclotron (S2C2), followed by four quadrupoles, two inside the S2C2 structure and two on the outside (Q1C and Q2C), the so-called "extraction slits" (SL1E), horizontal collimators protecting the downstream elements, and the rotating  Figure), the two quadrupoles Q1C and Q2C in red, the slit SL1E in green and the energy degrader are shown.
wheel energy degrader. The system losses during its operation mainly come from the S2C2 and the degrader. A complete system beam loss analysis can be found in Ref. [10].
The activation results obtained during the shielding design in Ref. [10] (considering a 20-year operation period for the centre) are used as input for the activated concrete decay radiation simulation following the workflow described in Section 2. The irradiation pattern used in the activation study followed the constraints of the IBA centre dimensioning study, which aimed to design a layout for an hourly dose rate of about 1 µSv/h in all public areas surrounding the vault and an annual dose below 1/3 mSv/year. Therefore, the irradiation conditions consist of a 150 nA S2C2 extracted current, an energy degradation to 70 MeV to maximise the secondary particle generation, and a maximum of 300 hours of operation per year distributed continuously throughout the year. This implies that the beam is still degraded to 70 MeV but the simulated S2C2 extracted current I S2C2 is given by, I S2C2 = 150 nA × 300 h/(365 * 24 h)) ≈ 5.137 nA. (4) The radiation protection quantity H * (10) is scored in the vault using a a 10 cm × 10 cm × 10 cm mesh for the same target nuclides as used in Section 3. The H * (10) evolution in the vault has been studied using the FISPACT-II activation results as scaling. The validity of this method depends on two conditions: either the assumption that the nuclide distribution inside the vault remains constant over time, or that any actual variations are below the targeted precision level. This hypothesis allows us to have only one time-consuming BDSIM simulation of the decay radiation per nuclide and shielding design for the entire study. Figures 7 and 8 show the H * (10) results integrated over 50 cm < y < 60 cm (with y being the vertical direction) for the concrete shielding made entirely in regular concrete and the optimised concrete shielding which uses a combination of regular concrete and LAC (labelled "realistic"). The results are presented four times after the system shutdown: 0 sec, 1 hour, 1 day and 7 days.
Several observations can be made on the results. First, the highest values are observed next to the degrader, which is expected behaviour as the degrader is the beamline element that generates most of the total secondary particles at the origin of the concrete activation. The shielding activation study showed that the walls next to the degrader are the most activated ones. The inner vault features large H * (10) values as a significant fraction of the S2C2 losses are uniformly distributed around the accelerator leading to the activation of all the vault inner walls [10,11]. The inner wall stops a fraction of the degrader and S2C2 losses leading to decreased H * (10) in the corridor to the treatment room. Comparing the shielding designs shows that the complete shielding in regular concrete delivers a higher ambient dose equivalent inside the vault. Figure 9 shows the evolution of the mean and maximum H * (10) values in the vault for both shielding configurations starting at the end of the irradiation period. The comparison between the two shielding configurations shows that LAC decreases by a factor of 3.5 the maximum H * (10) at the beginning of the cooling period. The evolution in time shows two bumps before reaching the long-term values after around 7 days of cooling. The first decrease is observed short times after the end of the irradiation in the range of minutes. A second decrease is observed several hours in the cooling period. Finally, the long-term values show a ratio of more than 12 between the maximum value obtained for the shielding in regular concrete and the one of the ProtherWal centre shielding. The long-term maximum H * (10) value obtained with the optimised shielding is comparable to the mean value obtained with the shielding in complete regular concrete. The results also show that using LAC in the shielding mix allows for a more significant gap between the maximum and the mean H * (10) values in the vault. These observations prove the beneficial impact of LAC in decreasing the hazardousness of the activated materials in the vault.
The BDSIM/FISPACT-II methodology allows access to information about the shielding activation that can explain how the ambient dose equivalent evolves in time for both shielding configurations. The most practical value, in this case, is the evolution of the activity in the vault shielding walls. The activity of a compound is determined by its clearance index, defined as the sum A i /CL i overall material radionuclides with A the specific activity and CL the clearance level allowed by the Belgian legislation. Figures 10 and 11 respectively show the evolution of the clearance index of the most radioactive part of the South left Wall (see Fig. 6) with time following the main radioactive nuclides in regular concrete and LAC. A first activity decrease is observed in both materials as extremely short-lived radioactive nuclides such as 28 Al, 38 K and 15 O decay. In regular concrete, 28 Al activity decreases until reaching a transient equilibrium with its parent nuclide 28 Mg. Hours after the end of the irradiation, short-lived nuclides as 28 Mg and 24 Na start to have their activity decline. After around 7 days, the activity of the materials is predominantly induced by long-lived nuclides as described in the activation studies results in Ref. [10] and the clearance index of both samples reaches the long-term values (delimited by the black dashed lines in the figures). The more significant decrease of the H * (10) values observed for the concrete mix using LAC between the beginning of the cooling period and 7 days later is explained by the greater activity of the extremely shortlived nuclides that decay in LAC compared to the longlived nuclides some minutes after the irradiation. This important aspect provides significant insights for shortterm radiation protection studies. The efficiency of LAC for mitigating the shielding activation and the related hazard effect is now fully explained.

Conclusion
An extension of the BDSIM/FISPACT-II methodology has been developed by implementing a Geant4 feature that allows the simulation of decay radiation in BDSIM. The BDSIM/FISPACT-II methodology computes the activation of any accelerator-based system by obtaining the secondary particles differential fluence at any point of interest using BDSIM Monte Carlo simulations and solving the rate equations with FISPACT-II for any irradiation and cooling periods patterns. The extension of the methodology allows using the results of these activation studies to simulate the decay radiation of the system and its surroundings at any time in the system lifetime.
The extended methodology is validated against the Monte Carlo code FLUKA by comparing the prompt and decay ambient dose equivalent, H * (10), values. The results show good agreement between the codes, with minor discrepancies explained by the differences in the implementation of the nuclear processes and the nuclear data.
The complete methodology is applied to the evolution of the ambient dose equivalent in the cyclotron vault of the to-be-built ProtherWal centre in Belgium. The results of the shielding activation obtained with the BDSIM/FISPACT-II methodology [10] during the design phase are used as input for the simulations discussed in this work. The optimised concrete shielding of the ProtherWal centre, composed of a mix of regular concrete and LAC, is compared in terms of its efficiency to mitigate the ambient dose equivalent in the cyclotron vault  to shielding that would be entirely made of regular concrete. The results show a net decrease of the ambient dose equivalent in the vault when using the optimised shielding, illustrating the efficiency of LAC in mitigating the hazardousness of activated concrete. Comparing the long and short-term activation between regular concrete and LAC gives a complete explanation of the origin of the ambient dose equivalent evolution. The BDSIM/FISPACT-II methodology now allows a single Monte Carlo model to be used throughout all design phases, from beamline optimisation and shielding design to activation studies of both the shielding and Fig. 11. Evolution of the clearance index of the most radioactive part of the South left Wall (see Fig. 6) with time following the main radioactive nuclides when the shielding is made of the optimised mixed between regular concrete and LAC. The clearance index obtained when studying the long-term activation of the spot is represented in dashed black line. The days are represented by vertical dashed gray lines.
beamline and ultimately to estimate relevant quantities for radiation protection studies at any point in time during the lifetime of the system from first irradiation to the decommissioning period.