Self-consistent numerical evaluation of concrete shielding activation for proton therapy systems

Due to the advancement of proton therapy for cancer treatment, there has been a worldwide increase in the construction of treatment facilities. Therapy centres are often coupled with clinical, biological or material-science research programs. Research activities require proton beams at energies spanning an extensive range with higher beam currents and longer irradiation times than clinical conditions. Additionally, next-generation proton therapy systems are evolving towards more compact designs. In addition to the increased centres’ workloads, reducing the system in size produces a more significant number of secondary particles per unit volume and time. Therefore, the activation level of materials constituting those future proton therapy centres is expected to be higher, increasing the ambient dose and the amount of radioactive waste collected at the end of a centre’s lifetime. These operating conditions pose new challenges for the shielding design and the reduction of the concrete activation. To tackle them, we propose a novel approach to seamlessly simulate all the processes relevant for the evaluation of the concrete shielding activation using, as an illustration, the Ion Beam Applications Proteus®\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^\circledR $$\end{document} One system. A realistic model of the system is developed using Beam Delivery Simulation (BDSIM), a Geant4-based particle tracking code. It allows a single model to simulate primary and secondary particle tracking in the beamline, its surroundings, and all particle-matter interactions. The code system and library database FISPACT-II allows the computation of the shielding activation by solving the rate equations using ENDF-compliant group library data for nuclear reactions, particle-induced or spontaneous fission yields, and radioactive decay. As input, FISPACT-II is provided with the secondary particle fluences scored using the BDSIM Monte Carlo simulations. This approach is applied to the proton therapy research centre of Charleroi, Belgium. Results compare the evolution of the clearance level and the long-lived nuclide concentrations throughout the facility lifetime when using regular concrete or the newly developed Low Activation Concrete (LAC). A comparison with the initial shielding dimensioning has been performed for all the shielding walls to validate the methodology and highlight the clear benefits of integrating LAC inserts in the shielding design. The effectiveness of coupling BDSIM and FISPACT-II gives a glimpse of the possibility of a complete activation study following the actual workloads of the centre, allowing a better assessment of the shielding activation level at any time of the facility lifespan.


Introduction
The first human proton therapy treatment was performed at the Lawrence Berkeley Laboratory in 1954. Since then, the number of proton therapy facilities has experienced a sharp worldwide increase, with more than 90 new centres commissioned over the past two decades [1]. Indeed, cancer treatments by proton therapy have proven to be highly efficient for specific types of tumours such as paediatric cancers, skull-based, sino-nasal malignancies or brain tumours [2].
Recent designs of hadron therapy installations are evolving towards more compact systems for which cyclotron-based solutions are favoured. Unfortunately, cyclotron-based proton therapy systems have the drawback of producing a large number of secondary particles, mainly neutrons [3]. Those neutrons are generated during the acceleration, the energy degradation processes and, to a lesser extent, during the protons propagation along the beamline. Appropriate shielding is necessary to mitigate the propagation of those secondary particles in the treatment room or the vicinity of the machine to limit the extra dose received by the patients and the medical staff. The accelerator and the energy degradation system are the two most important radiation sources. Their isolation in a concrete vault separated from the patient room is essential. However, the interactions with the trace impurities present in the concrete produce short-and long-lived nuclides through nuclear reactions, mainly neutron capture and spallation. As proton therapy a e-mail: eliott.ramoisiaux@ulb.be (corresponding author)

BDSIM and FISPACT-II integration for concrete shielding activation evaluation
Shielding activation studies need to consider the complex nuclear reactions experienced by the primary and secondary particles while keeping inventories of the short-and long-lived activated nuclides.
The primary beam interactions with the beamline elements generate secondary neutrons. The energy spectrum of the secondary neutrons is continuous from thermal energy up to the energy of the primary beam. Those neutrons propagate until they reach the concrete shielding, where they interact via nuclear reactions, mainly neutron capture and spallation, and produce radioactive elements of which only a few are long-lived and responsible for the long-term activation of the shielding.
Low energy and thermal neutrons interact predominantly by neutron captures and generate long-lived isotopes, of which the most important ones are produced by interaction with the trace amounts of europium, cobalt and caesium present in the concrete, as detailed in Table 1. Their concentration varies with the aggregates used during the concrete casting process but can be estimated from decommissioning measurements. As a baseline, we use the results of [9], obtained from measurements at a decommissioned cyclotron and Monte Carlo simulations using the TART98 code [10].
On the other hand, high-energy neutrons produce long-lived isotopes through spallation reactions. Table 2 lists the most important ones along with their respective half-lives. The spallation cross-sections are usually smaller than those of neutron capture and have a high characteristic energy threshold in the 10 MeV to 30 MeV range. However, the number of long-lived isotopes produced by spallation is non-negligible due to the significant production of high-energy spallation neutrons from proton losses in the beamline. Proton therapy shielding design studies are typically performed using general-purpose Monte Carlo codes such as MCNPX [12], GEANT4 [13,14], or FLUKA [15]. These studies use beam losses from tracking codes as discrete sources in passive radiation transport simulations (i.e. without magnetic fields). A first enhancement has been proposed by Tesse [16] with simulations directly coupling the magnetic tracking of primary and secondary particles in the beamline elements with the computation of the particlematter interactions. This enhancement significantly improves the accuracy of the secondary particles production, which is essential for activation studies of compact systems such as single-room proton therapy installations. This method is implemented using BDSIM, a program that combines standard high-energy physics codes such as GEANT4, ROOT [17], and CLHEP to create a computational model of a particle accelerator that uses accurate accelerator tracking routines with all the physics processes available in GEANT4.
BDSIM provides default geometries that allow the user to quickly model any particle accelerator-based systems. If a more complex geometry, either for the beamline or the surrounding elements, is required, Geometry Description Markup Language (GDML) files [18] of the geometries can be implemented in the BDSIM model of the proton therapy system. The protons lost during the beam propagation along the beamline interact with the elements of the studied system, thus directly playing the role of secondary particle sources. As a direct improvement to the model, the origin of the secondary particle fluence is therefore distributed in a continuous manner along the beamline.
In the work of Tesse [16], the long-term activation is computed with a BDSIM scorer that computes the reaction rate of any reaction whose cross-section has been provided. This method is highly time-consuming as for each reaction, a specific scorer must be defined. Furthermore, if a radionuclide is forgotten or added to the study, an entirely new simulation should be done again. The time evolution decay was implemented in an ad-hoc manner using cross-section data for the selected isotopes presented above.
Our methodology for proton therapy shielding design uses the enhancements provided by BDSIM while also improving the activation computation process by coupling with the FISPACT-II code for the activation computation [6]. FISPACT-II is a code system and nuclear database able to model activation-transmutation processes induced by nuclear reactions and decays and keep track of the nuclides in time-dependent inventories. FISPACT-II requires, as input, the incident particle energy-differential fluence. It uses the Livermore Solver for Ordinary Differential Equations (LSODE) to integrate the rate equation in time. In addition, a wide range of radiological output quantities such as activity, output heat, dose rate, clearance index or pathway analysis can also be computed. The shielding activation is characterised by the clearance index of each nuclide, defined as the ratio between the specific activity and the clearance level allowed in the Belgian legislation in Bq/g. The FISPACT-II output is processed using an in-house Python analysis package.
The new method consists of several stages and is described in Fig. 1. First, the user generates the beamline and the shielding models by combining the default accelerator component geometries provided by BDSIM and more complex geometries imported from GDML files. Second, 3D meshes are positioned in the BDSIM model for the secondary neutrons differential fluence scoring. Then, BDSIM simulations are run using a given input beam and an appropriate physics list for the particle-matter Monte Carlo simulations. Finally, the differential fluences are extracted for each scored 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-sections database. This method separates the Monte Carlo simulations and the activation computation, increasing flexibility and time efficiency. Fig. 1 Diagram of the BDSIM/FISPACT-II activation computation methodology. The user inputs are in red, the codes are in green and the outputs are in blue. BDSIM requires information about the input beam, the beamline and the shielding model with a chosen physics list. Then the secondary fluxes extracted from the BDSIM simulations are provided to FISPACT-II which, using a given set of cross section databases, provides nuclide inventories and related radiological outputs

BDSIM model of the IBA Proteus ® One and beam loss fluence simulation
The P1 system is a single-room proton therapy centre equipped with a 230 MeV SuperConducting Synchro-Cyclotron (S2C2), followed by an extraction beamline comprising an energy degrader mounted on a wheel and a rotating compact gantry (CGTR) with a 220 degrees rotation range. A complete and realistic model of the P1 has been implemented in BDSIM, as reported in detail in Ref. [19]. This model has been fully validated in terms of beam envelope in the transport beamline, dose deposition profile (Bragg peak) at isocenter and global transmission efficiency. The latter is crucial as it validates the beam loss pattern, a crucial input for the activation computation.
The specific shielding of the Charleroi centre and its different types of concrete has been implemented in BDSIM as described in Sect. 4. Figure 2 is a ParaView [20] representation of the P1 BDSIM model where the beamline is placed in its concrete shielding with the primary proton beam represented in blue. This visualisation is achieved using an in-house Python-based toolbox [21] that has been developed to provide conversion mechanisms, from the Geant4 geometry format and ROOT data to a 3D visualisation using the VTK framework [22]. It supports the complex geometries used by the BDSIM models, the transported beam tracks, and the scored quantities (see Figs. 16,21). BDSIM default geometries are used for the majority of the beamline elements, while for more complex geometries such as the accelerator, the energy degrader, and the last bending magnet B3G, detailed Computer-Aided Design (CAD) geometries have been implemented. The main improvement made on the previous BDSIM model from Tesse [16] was the implementation of the beamline components correct geometry, which improved considerably the simulation accuracy. To further improve the model, a detailed magnetic field map is used for B3G. A thorough study of the P1 system losses is the key to understanding the origin of the activation. Proteus One geometry as implemented in the BDSIM model. The GDML concrete shielding is represented around the P1 with regular concrete in light brown and the Low-Activation Concrete in light grey. The GDML SuperConducting Synchro-Cyclotron (S2C2), the extraction line, the energy degrader system and the rotating gantry are shown. The concrete shielding wall separating the cyclotron vault from the treatment area is visible with the cylindrical opening. A cubic water volume is placed at the isocenter for energy deposition scoring. The location of measurement devices -the beam profile monitors (P1E, P1-4G), the ionisation chambers (IC) and the isocenter-are also indicated (green)

Accelerator beam distribution and losses
The protons are accelerated to 230 MeV at a maximum current of 132 nA with a central field of 5.7 T, an RF frequency variation in the range of 90 to 60 MHz and a pulsed beam delivered at 1 kHz with a pulse duration of around 10 μs. The complex acceleration process has not been simulated in our model. To obtain the beam distribution at the exit of the S2C2, a separated three-stage simulation has been performed using an OPERA3D magnetic model of the accelerator and two IBA in-house simulation codes, Advanced Orbits Code (AOC) and "phase space motion" [23,24]. AOC integrates the equations of motion in a given magnetic field map with superimposed time variable or fixed electric field. In contrast, "phase space motion" tracks the orbit coordinates, the RF phase and the vertical motion. AOC is used to simulate in detail the injection process that occurs up to 5 MeV. Then, "phase space motion" performs the acceleration simulation from 5 MeV to 225 MeV. Finally, the extraction up to 230 MeV is simulated again by AOC based on the RF bucket and emittance evaluated at 225 MeV. The resulting phase space beam distribution is presented in Fig. 3. It can be seen that the horizontal distribution largely deviates from the usually assumed Gaussian distribution. This deviation originates from the specifics of the acceleration and extraction processes which strongly impact the horizontal plane. In the vertical plane, the distribution follows a Gaussian distribution with an emittance of 1.85 mm mrad. This distribution is used as input for our BDSIM model.
The S2C2 beam loss model used for our simulations follows the one used for the conservative shielding design [11] which considers that the worst-case scenario corresponds to an efficiency of 30% for the S2C2 extraction. The reamining of the beam is divided into three-loss sources. The first source corresponds to losses occurring close to the ion source. These losses are neglected in this paper as the protons involved are in an energy range close to the thermal energy (0.025 eV), which is too low to induce secondary neutrons production. Therefore, to follow a conservative approach, the remaining 70% of the beam will be distributed in the second and third loss sources. The second source corresponds to protons defocused in the vertical plane hitting the chamber. These losses occur when the proton beam is near its maximal energy, thus close to the extraction radius. They are assumed to be uniformly distributed along the S2C2 circumference and corresponds to 25% of the initial beam. The third and final source, corresponding to 45% of the initial beam, occurs when the accelerated beam interacts with the septum at the maximal energy of 230 MeV during the extraction phase. The different loss proportions and their respective spatial distribution in the stainless steel accelerator are considered in the BDSIM simulation of the beam losses, which will generate the secondary particles. The S2C2 loss distribution and the generated secondary neutrons are presented in Fig. 4.
The energy distribution of the secondary neutrons produced by the S2C2 losses, scored outside the accelerator, is presented in Fig. 5. The collisions generate high-energy neutrons (≥ 1 MeV) from spallation reactions. A spallation reaction is divided into two stages: intra-nuclear cascade and deexcitation. Different cascade and deexcitation models are implemented in Geant4 and can be chosen through the definition of the "physics lists" [25][26][27]. Ref. [28] recommends the use of the physics lists "g4QGSP_BIC_HP_EMZ" for global activation studies. The energy distribution of the neutrons shows that even if some high-energy neutrons still reach out of the S2C2, a majority of the emitted neutrons has thermalised in the accelerator structure itself and populates the low energy part of the spectrum. This distribution suggests that the S2C2 losses will primarily generate neutron capture-induced activation rather than spallation-induced (see Sect. 2).

Beamline losses
The beamline model is divided into two parts. First is the static extraction line, which consists in focusing the extracted beam on the degrader. The proton interactions with the degrader allow modulating the beam energy from the maximal energy of 230 MeV to any desired energy down to 70 MeV. The extraction line is entirely placed inside the vault as the energy degradation generates many secondary particles. Second, the Compact Gantry TRansport system (CGTR) transports the transmitted beam from the degrader to the isocenter in the treatment room. During its propagation, the beam inevitably interacts with the beamline elements. Those interactions are critical at the level of the extraction line as the S2C2 extracted beam has an inherent wide distribution along the horizontal plane (see Fig. 3). The extraction line elements are calibrated in such a way that the beam reaches a double waist in the degrader, which minimises the emittance increase [29]. Following the degrader, a circular tantalum collimator intercepts the deflected protons. The main CGTR elements which interact with the beam are the divergence (SL1G and SL2G) and momentum (SL3G) slits (presented in Fig. 2). The role of the divergence slits is to protect the following beamline elements from the diverging protons. The role of the momentum slits is to protect the machine and the target if the proton beam energy is too far from the target energy.

Energy degrader
The energy degrader consists of a rotating wheel composed of blocks of variable thicknesses and materials (aluminium, graphite and beryllium). Aluminium is used to degrade energy between 230 MeV and 220 MeV, graphite between 220 MeV and 120 MeV and beryllium between 120 MeV and 70 MeV. These three materials have different physical properties and vary in transmission efficiency, emittance increase, and secondary particle yields. The accurate modelling of the degrader is of crucial importance for activation studies as it is the most interacting element of the P1 system [30]. The energy degrader geometry is directly imported in BDSIM using a GDML file.
The BDSIM model of the degrader is represented in Fig. 6 with the surrounding beamline elements. The proton beam (blue) exits the extraction line vacuum pipe from left to right and goes through an air gap. The beam interacts with the degrader and reaches its target energy with a particular transmission and emittance change. The diverging protons are then mainly stopped by the circular collimator (dark green). The vacuum in the beampipe is maintained using two 1 μm titanium foils.
Most of the activation is produced when the proton beam energy is degraded from 230 MeV towards 70 MeV. This critical beam condition has been simulated in BDSIM to provide an overview of the degrader effect on the beam in terms of shape, energy distribution and secondary particle production. The results upstream and downstream of the collimator for these three quantities are represented in Figs. 7 and 8.
At the degrader entry, the proton beam has not reached the double waist, and the beam distribution is still focusing along the X-axis while diverging along the Y-axis. The transmitted beam after the degrader has a Gaussian distribution whose width exceeds the collimator aperture. The collimator then catches the protons with high angular deviation, and only the core of the distribution continues through the beamline. This effect can be observed by comparing the beam size before and after the degrader.  6 Realistic model of the IBA energy degrader combining three degrader materials (aluminium (yellow), graphite (light green) and beryllium (light blue)) as well as the beam stop element (grey). The upstream quadrupole (red) and its drift, encapsulating a 1 μm titanium extraction foil keeping the void inside extraction line, are represented on the left. Following the degrader, another 1 μm titanium extraction foil protecting the CGTR void and the circular collimator (dark green) are placed before a drift, an H-kicker (dark blue) and a V-kicker (pink). The coils of the different magnets are also integrated in the geometry (orange). A sample of proton tracks is shown in blue The comparison of the energy proton beam distribution before and after the degrader is presented in Fig. 8a. The phenomenon called energy straggling [3] can be observed as, when the 230 MeV proton beam interacts with the degrader atoms, it experiences losses through a large number of individual interactions. A statistical error appears as each proton scatters differently in and out the degrader leading to an energy distribution centered around the targeted 70 MeV. The higher energies in the distribution are populated by diverging protons that are scattered out of the degrader before reaching its end. The impact of the collimator on the energy distribution is easily observed: the diverging protons of higher energy are stopped.
Spallation reactions occur during the collision of the proton beam with the degrader emitting high-energy neutrons. The resulting neutron energy spectrum is presented in Fig. 8b. In comparison with the neutron energy distribution of the S2C2 losses presented in Fig. 5, the amount of high-energy neutrons is greater as the produced neutrons could only interact with the degrader itself before reaching the air, which limits the neutron thermalisation. It is therefore expected that the beamline losses generate more spallation than neutron capture-induced activation (see Sect. 2).

Beamline
Proton beams of given energy and desired beam properties at the isocenter are obtained by calibrating the beamline elements for the target energy. These elements are adjusted to optimise the beamline efficiency, beam size at isocenter, and the beam energy spread in all energy conditions. The beam optics and loss pattern obtained along the beamline are thus expected to be energy-dependent. Figures 9 and 10 compare the beam envelope evolution along the beamline for a beam energy of, respectively, 70 MeV and 210 MeV. The loss patterns at these two energies are presented in Figs. 11 and 12.
The beam envelopes at 70 MeV and 210 MeV (Figs. 9, 10) show differences at characteristic points of the beamline. To compensate for the low transmission of the degrader at 70 MeV and keep a relative smooth transmission evolution with respect to the beam energy, the SL1E slit opening along the X-axis is wider at low energy than at high. The horizontal axis has been chosen to perform this transmission modulation as the beam distribution is wider in the horizontal plane (Fig. 3). The 70 MeV transmitted beam shows a greater divergence after the degrader, which leads to a larger envelope along the beamline and thus to higher losses: the loss patterns are then energy-dependent. The isocenter beam size also depends on the transmitted energy and is adjusted to be symmetrical in   Loss patterns along the beamline for a mean beam energy of 70 MeV on the left and for a target beam energy of 210 MeV on the right. Primary first hits represent the fraction of the primary beam that undergoes its first hit with a material. Primary loss represents the fraction of the primary beam that is lost along the beamline. The energy deposition represents the energy deposited by the primaries and the secondaries both planes. This symmetry is a design constraint imposed to facilitate Pencil Beam Scanning (PBS) beam models for Treatment Planning Systems (TPS) [31].
The losses and transmission evolution along the beamline at 70 MeV and 210 MeV (Fig. 11) shows expected behaviour related to the previous envelope analysis. The SL1E slit causes, as expected, more losses at higher energy due to its smaller opening, while the degrader and collimator generate more significant losses at low energy as the degradation and the beam divergence are more critical. The high divergence of the transmitted beam at 70 MeV leads to many protons colliding with the first elements of the CGTR before reaching the divergence slits. On the other hand, at 210 MeV, the divergent protons primarily interact with the SL1G and SL2G divergence slits.
The energy depositions along the beamline at 70 and 210 MeV of the primary and secondary particles are represented in Fig. 12. Additional information about the primary first hits, which represent the fraction of the primary beam that undergoes its first hit, and the primary loss, which describes the fraction of the primary beam that is lost, are also given. One can observe, as expected, that no primary first hit occurs after the degrader. Coherent results with the previous analyses are observed: more primary hits are found at the level of SL1E at 210 MeV than 70 MeV. Losses are distributed along the beamline with hotspots depending on the transmitted energy. The protons that do not continue until the isocenter deviate from the main trajectory and interact with the beamline elements or the external shielding, which generates important fluxes of secondary neutrons continuously along the beamline. It illustrates the model enhancement in the loss patterns and the secondary particle generation in comparison with the point source method previously used in shielding dimensioning studies [11].

Shielding design and neutron fluence at the Charleroi proton therapy centre
The future proton therapy facility of Charleroi will have a workload quite different from the conventional Proteus One centres around the globe. Indeed, the clinical activities will be initially limited, and research activities will require intense proton beams for long and continuous periods (3-4h/day) of irradiating materials and biological tissue. A shielding design compliant with the legal radiation protection limits in Belgium has been developed by IBA using MCNPX 2.7.0 [11]. The shielding study was divided into two cases corresponding to the most critical irradiation conditions of the S2C2 vault and the treatment room. The vault design was realised using the simulation condition of the maximal energy degradation from 230 MeV to 70 MeV for which the highest secondary neutron production occurs at the degrader level.
On the other hand, the shielding of the treatment room was designed with a maximal energy beam of 230 MeV reaching a water phantom at the isocenter. A safety margin has been taken at the level of the S2C2 as the extracted beam current was considered to be 150 nA compared to the maximal extracted current of 132 nA. The IBA dimensioning studies aimed to design a layout for an hourly dose rate of about 1 µSv/h in all public areas surrounding the vault. The Belgian annual dose limit corresponding to 1/3 mSv/year limits the system to run for about 300 hours a year.
The concrete activation was studied in these two most critical conditions for a period of 20 years. This assessment has an essential role in the shielding design as it directly impacts the centre decommissioning costs. Indeed, concrete volumes activated with longlived nuclides whose specific activity (A) exceeds the clearance level (CL) allowed by the Belgian legislation are considered low-level nuclear waste. The clearance index, defined as the sum A i /C L i over all the material radionuclides, characterises the activation level and must stay below 1 to be considered as normal waste. The major isotopes produced in concrete are listed in Table 3 with their corresponding clearance level.
IBA performed a first activation study for the future centre concrete shielding with MCNPX and determined the amount of nuclear waste that would be generated inside the facility after 20 years of operation with the expected workload of Charleroi. The physics list and the nuclear library data used for the simulations were, respectively, "MNCPX_BERT" and "ENDFB.VII". The results showed an increase by a factor of 3 of the clearance level values obtained in the vault walls compared to the usual clinical centres. This activity increase corresponds to an additional 10 cm thickness of activated concrete in the walls. With the idea to let the activated concrete "cool down" and thus decrease the amount of activated concrete, a 20-year cooling period was studied. Because the vault clearance level decreased slowly and globally remained above a value of one, a cooling down period could not be the solution needed to cope with the significant amount of nuclear waste. The solution proposed was to use the newly developed Low-Activation Concrete (LAC) in the S2C2 vault [11]. The LAC was developed based on neutron activation analysis of different aggregates, sand and types of cement using the BR1 nuclear reactor from SCK-CEN in Belgium. In that study, the concentrations in Eu, Co and Cs were measured two months after the irradiations using a high-purity germanium spectrometer. The activity generated inside the LAC is lower compared to regular concrete as it exhibits low concentration in the trace elements at the origin of long-lived nuclides. The reduction factor for the concentration of the trace elements Eu, Co and Cs are, respectively, 30, 100 and 30 compared to regular concrete (see Table 4).
The LAC presents also lower concentration in Na, Al and Si elements responsible for the generation of radioactive nuclides through spallation reaction. Table 5 presents the comparison of the atomic composition between the regular concrete and the LAC.
An improved concrete shielding design was therefore obtained by replacing a certain layer of regular concrete (2.37 g/cm 3 ) with LAC (2.18 g/cm 3 ) in the vault walls. An activation study was performed by IBA using the new concrete shielding. As expected, the clearance level of almost of the walls stayed well below the transition value of one except for the floor and the walls next to the degrader. These areas require a cooling period of more than 10 years to be considered as regular waste which is too long for the small volumes concerned and will therefore be considered as potential nuclear waste at the dismantling phase. Even though a complete activation-free shielding could not be designed, the implementation of LAC avoids the production of a large amount of nuclear waste, and the price of the LAC (∼ 240 e/m 3 with respect to regular concrete at ∼ 100 e/m 3 ) is an investment much lower than the possible decommissioning costs. The Eu, Cs and Co concentrations used by IBA in their study were slightly different from the values presented in Table 5. These discrepancies are minor compared to the concentration variation between regular concrete and LAC and our calculations have showns that they do not impact the conclusions of our study.
The new shielding design was modelled in BDSIM using pyg4ometry, a Python library that enables users to rapidly create GDML-based geometry [32,33]. The regular and Low-Activation concrete have been implemented following the atomic composition presented in Table 5. The BDSIM model of the Charleroi shielding is presented in Fig. 13.
The "North Wall" and the "South left Wall" concrete shielding walls of the vault, presented in Fig. 14, will be thoroughly studied. These walls are close to the two main secondary particle sources, the S2C2 and the degrader. Studying the activation for these two walls allows assessing separately the activation coming from either the acceleration or the energy degradation processes.
The inventory code FISPACT-II requires the secondary neutrons differential fluence following some predefined energy group structure (ex: CCFE-709, UKAEA-1102, ...) as input for its inventory calculation. The so-called "4D scoring" option has been  Fig. 13 Realistic model of the Charleroi Concrete Shielding geometry as implemented in BDSIM. Regular concrete is shown in light brown while the Low-Activation concrete is shown in light grey. Some walls have been made transparent to view inside the shielding implemented in BDSIM. It defines a 3D spatial scoring mesh and allows to score a differential quantity, in this case, the energy differential fluence. In addition, the binning for the differential quantity can be chosen to follow a linear or logarithmic scale or an explicitly defined list of groups, as required for the FISPACT-II neutron group structure [21]. This new feature makes use of the generic boost::histogram library [34] and allows an event-by-event serialisation and storage in the ROOT data format. The nuclear reaction cross-sections handled by FISPACT-II are defined under different possible binning methods, including the CCFE-709 energy group structure, which is chosen in our simulations for the secondary neutrons differential fluence scoring. A 20-years irradiation period of the Charleroi proton therapy shielding is chosen to apply the proposed activation computation methodology. The worst irradiation case scenario for the vault shielding is considered [11]. These irradiation conditions have been implemented in BDSIM to simulate the corresponding secondary neutron generation. The physics list used in the BDSIM simulation is the "g4QGSP_BIC_HP_EMZ" as recommended in Ref. [28].The binary-cascade intranuclear model ("BIC") has been chosen over the Bertini model ("BERT") that was used in the MCNPX simulations. A comparative study showed that the BIC model generated more high-energy neutrons (≥ 1 MeV), inducing greater spallation reactions and thus improving the activation computation, especially for short-term studies. The differential fluence has been scored in all the walls of the vault using importance sampling to reach acceptable statistics for all depths (see "Appendix"). Figure 15a, b display the scoring results at the entry and the exit of the North Wall and South left Wall scorer meshes designated by red rectangles in Fig. 14. The neutron differential fluence generated from either the S2C2 or the beamline losses is presented. As expected by their respective position in the vault, the neutron fluence originating from the S2C2 losses is more important for the North Wall, while the beamline losses fluence is predominant for the South left Wall. The high-energy neutrons generated in the degrader (see Fig. 8b) are observed in the South left Wall beamline losses differential fluence in contrast to the North Wall. The high-energy neutrons disappearance is explained by the position of the North Wall far away from the degrader and separated by the beamline elements and the S2C2. Those neutrons undergo interactions while travelling toward the North Wall and thus populate the lower energies of the spectrum. This difference in the high-energy neutron spectrum explains why spallation reactions are less expected in the North Wall than in the South left Wall, leading to a clear difference in the monitored activation in each wall.

Self-consistent shielding activation computation
The concrete shielding activation is characterised by the clearance level computed using FISPACT-II based on the secondary particle differential fluences extracted from the BDSIM simulation of the vault walls following the irradiation conditions presented in Sect. 4. Internally, FISPACT-II uses the differential fluence to compute "collapsed" cross-sections that will allow the control of the irradiation rate through the total fluence value. The "collapsed" cross-sectionX is computed based on the cross-section X i and the projectile flux φ i in the energy group i ∈ [1, N ] followinḡ This methodology is efficient as it allows flexibility in the activation study: once the secondary particle spectrum is computed with an acceptable statistic, any irradiation configuration can be studied. The FISPACT-II simulation was run using the differential fluence extracted from the BDSIM simulation as input for a 20-year period of constant irradiation. The nuclear data library "TENDL17" was chosen for this work over "ENDFB.VII", which was used in the MCNPX simulations. The ENDFB.VII library, being a transport purpose library, does not have all the reaction channels defined. The isomer production is neglected, leading to an over-evaluation of some nuclides clearance index (ex: Eu 152 ). Finally, the data for most nuclides are defined as up to 20 MeV, not covering the incident neutron spectrum. A complete activation study was realised for the vault, and the results obtained for the South left Wall and the North Wall are presented in detail in the following two sections. Then, a summary of the results obtained for the totality of the walls is presented. The BDSIM/FISPACT-II workflow results are compared to those previously obtained by IBA for the centre dimensioning. A discrepancy between the two workflows is expected due to the difference in the Monte Carlo library used (Geant4 vs MCNPX), in the physics lists chosen (see Sect. 4), in the nuclear library chosen (TENDL17 vs ENDFB.VII), in the LAC composition and in the P1 system losses model accuracy. Considering that, a comparison of the main results is still possible to validate the new methodology workflow.

Activation of the South Left Wall
Activation of the South left Wall has been simulated, and the clearance index in the wall has been computed and visualised with ParaView in Fig. 16. The clearance index of the different cells in the wall has been highlighted in blue when below the limit value of 1 and in red for parts of the wall that have a higher clearance index, therefore reaching the status of radioactive waste. The P1 BDSIM model, its shielding and its primary proton beam are also represented. The clearance index values along the thickness of the most radioactive part of the South left Wall are presented in Fig. 17 and compared to the IBA design results realised with MCNPX. As expected from the localisation of the wall, the beamline losses, which are mostly generated at the level of the degrader, govern the wall activation. The results show that the South left Wall will be completely activated at the end of the 20 years of the centre operation. This observation agrees with the results obtained by IBA Fig. 16 ParaView visualisation of the P1 BDSIM model with the clearance index highlighted for the South left Wall. The colour bar has been chosen to have the blue colour associated to clearance index values below 1 while the red colour has been assigned with clearance index values over 1 Fig. 17 Evolution of the clearance index along the thickness of the most radioactive part of the South left Wall while using regular concrete. The clearance index related to the S2C2 losses and beamline losses are, respectively, represented in green and blue while the total is represented in red. The clearance index limit characterising radioactive waste is represented in black (a) (b) Fig. 18 Evolution of the clearance index of the most radioactive part of the South left Wall with time following the main radioactive nuclides (left) and the activation reactions (right), using regular concrete during their dimensioning. However, a non-negligible difference in the clearance level amplitude can be observed between the two methodologies. The evolution with time of the clearance index of the most radioactive part of the South left Wall has been studied for a 20 year-cooling period filtered either by the main radioactive nuclides (Fig. 18a) or by the activation reactions (Fig. 18b). The analysis of the results shows that the clearance level does not reach a value below 1, implying that a 20-year period would not be enough for the wall to be considered normal waste. The clearance level evolution by nuclides clearly shows that the overall remaining activity is due to the presence of Eu 152 generated by the neutron capture of low-energy neutrons. The Eu 152 half-life above 13 years explains why the cooling down method is not suitable. The same conclusion was obtained by IBA during their shielding design and started the development of the LAC to prevent as much as possible the activation of the wall from the start. At shorter times after the start of the cooling period, the main radioactive nuclide is Na 22 , which comes from the spallation reactions in the concrete induced by the high-energy neutrons generated from the degrader (see Fig. 8b). Na 22 has a shorter half-life (≈ 2.6 years) and becomes the second most radioactive nuclide. Na 22 should have an important role in short-term activation studies of the vault.
The same study is performed by replacing regular concrete with LAC. The clearance level evolution along the thickness of the Wall using LAC is presented in Fig. 19. The clearance index values reached are much lower than those obtained when using regular concrete. Clearance level values below 1 are observed after around 1.1 m depth in the wall without any cooling time needed. This observation supports the fact that using LAC is a suitable solution to reduce considerably the amount of radioactive waste generated. The MNCPX results start at similar values, though they quickly decrease to an activated depth of around 0.8 m. This difference is directly linked to the choice of the nuclear library for the activation computation as comparative studies showed that using the ENDFB.VII library gives similar activated depth results.
A cooling period on the South left Wall when using LAC can also be studied. The results are presented in Fig. 20a, b and show that a cooling period of approximately 10 years would be required to reach a normal level of activity in the entire wall. It clearly shows that the use of LAC significantly decreases the hazards of the waste amassed as a shorter time would be required to recover a normal activity by comparison when regular concrete is used as a shielding material. As observed when regular concrete was used, the Na 22 nuclide is the most radioactive nuclide in the wall but quickly decreases, while in the long term, it is still the Eu 152 that gives the total activity.

Activation of the North Wall
The same comparative activation study considering the use of regular concrete and LAC as shielding material has been performed on the North Wall. Figure 21 presents a ParaView visualisation of the North Wall activation results when regular concrete is used in the walls after 20 years of operation. Figure 22 shows the evolution of the clearance index along with the thickness of the most radioactive part of the North Wall. The S2C2 losses govern the wall activation due to its proximity. After 20 years of operation and without any cooling time, the wall is activated until a depth of 0.5 m. The MCNPX results have a higher starting value but reach the same activated depths.
The North Wall activation is also studied during a cooling period of 20 years. The results are shown in Fig. 23a, b, respectively, representing the evolution with time of the clearance index of the most radioactive part of the North Wall following either the main radioactive nuclides or the activation reactions at the origin of the radioactive nuclides. As observed for the South left Wall, a 20 year-cooling period is not enough for the clearance index to reach a value below 1. The most critical radioactive nuclides in the wall are the Eu 152 and the Co 60 with respective half-lives of about 13 and 5 years. They are both induced by neutron captures of low-energy neutrons in the concrete. Spallation-induced activation is low in this wall as the main source of secondary particles reaching the wall is the S2C2 characterised by secondary neutron spectrum principally populated by low-energy neutrons as shown in Fig. 5.  Evolution of the clearance index along the thickness of the most radioactive part of the North Wall while using LAC. The clearance index related to the S2C2 losses and beamline losses are, respectively, represented in green and blue while the total is represented in red. Only the clearance index related to the total losses obatined with MCNPX was available for comparison. The clearance index limit characterising radioactive waste is represented in black Figure 24 shows the evolution of the clearance index along with the thickness of the most radioactive part of the North Wall when using LAC. The clearance index of the wall rapidly drops under 1, revealing that the use of LAC for the North Wall allows to eliminate completely the activation problem without any cooling period.

Complete activation study of the vault walls
The activation study of the South left Wall, and the North Wall of the vault have been presented thoroughly. It showed that the use of LAC could either completely prevent walls from becoming waste or decrease the activation of the remaining activated volumes and thus considerably limit their hazards. The same activation study has been performed on all the walls inside the vault in order to quantify precisely the beneficial impact of LAC as a shielding material. The activated volumes amassed after 20 years of operation in each wall when using regular concrete or LAC is compared in Fig. 25. The efficiency of LAC is observed: the activated volumes drop from about 109.5 m 3 with the regular concrete to about 5.4 m 3 with the LAC. This result proves the benefit of using LAC for the shielding of the future Charleroi proton therapy centre.
On the other hand, the entire concrete shielding of the vault cannot be entirely built using LAC as it is 4 times more expensive than regular concrete. The solution is to design a specific shielding combining regular concrete and LAC based on insights extracted from the vault walls activation studies. The required information for such design is the maximally activated thickness reached if only regular concrete were used as a shielding material. Figure 26 presents the maximally activated thickness for the vault walls when only regular concrete or LAC is used as a shielding material. It directly characterises the exact thickness of LAC needed for each wall to minimise waste production and the shielding cost.
In Fig. 27, we compare the proportion of LAC in each wall obtained from the IBA shielding design and the results obtained with the new BDSIM/FISPACT-II method. The two workflows give identical results except for the Roof and the East Wall. These minor discrepancies would correspond to an increase of 1.5 m 3 in the activated volumes at the centre decommissioning. This increase is negligible compared to the amount of radioactive waste avoided through the use of LAC. Therefore, the new methodology is

Conclusion
The results and their comparison with prior MNCPX studies show that a self-consistent numerical evaluation of the concrete shielding activation can be obtained by coupling the multi-purpose simulation code BDSIM and the inventory code FISPACT-II. BDSIM allows particle tracking inside all the components of a hadron therapy system while also simulating the particle-matter interactions using the Geant4 library, and provides a significantly more accurate tracking of primary and secondary particles. FISPACT-II provides, in addition to nuclides inventory calculations, a large set of derived radiological quantities such as the dose rate, ingestion doses, inhalation doses and clearance index.
The dimensioning of the concrete shielding of the future proton therapy of Charleroi has been presented in detail from the activation in concrete to the specific development of Low Activation Concrete. The activation study of the future proton therapy centre of Charleroi has been realised using a complete BDSIM model of an IBA Proteus ® One system and its specific shielding. The simulations of the proton beam losses during the S2C2 acceleration process and during the beam transportation to the isocenter have been used to obtain reliable secondary particles generation patterns. The general beamline losses have been studied with particular attention to the degrader, as it generates the most significant number of secondary particles.
The differential fluence of the secondary neutrons has been extracted for all the vault walls by using scoring meshes in both position and energy. The spectra corresponding to the South left Wall and the North Wall, walls positioned close to the two main secondary particle sources, the S2C2 and the degrader, have been thoroughly analysed. Their long-term activation has been studied with and without the Low-Activation Concrete. The results have been validated against prior MCNPX simulations in terms of clearance level, quantity of concrete above a clearance level of 1, and comparison between the behavior of LAC and regular concrete. The evolution of the clearance index of the main radioactive nuclides was also analysed through a 20-year cooling period to characterise the evolution of the remaining waste activity level.
It is shown that the LAC considerably limits radioactive waste generation at the end of life. It reduces the cooling period required for the activated materials to reach acceptable activity levels. Therefore, the implementation of LAC in the shielding contributes to limiting the environmental impact and the long-term cost.
The proposed methodology will also allow precise tracking of the shielding activation over the lifespan of the centre using the actual irradiation history. Only the FISPACT-II simulation is affected as it uses the beam load as input and an integration between our tools and control system logging data is foreseen.
The results, in terms of simulated neutron fluences and activation levels in the vault are meant to be validated with measurements. A detailed experimental program featuring removable regular and LAC concrete cores to be installed at specific locations in the vault is part of the research program to be carried out at the Charleroi centre.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix: Importance sampling
Monte-Carlo simulations are usually combined with variance reduction techniques to improve the results statistics while maintaining adequate simulation time [35]. It is especially required when complex and large geometries are involved in the studied applications. While a large set of variance reduction techniques exists, the most popular is the Geometric Splitting/Russian roulette. For this technique, the region of interest is divided into cells of different importance values, and each particle starts with a weight of one. The Geometric Splitting is applied to a particle of weight w when it crosses the boundary from a cell of lower importance I l to higher importance I h . It consists of splitting the particle into n I h /I l identical particles of weight w/n (if n is not an integer, the number of split particles is chosen probabilistically to maintain the importance ratio statistically). On the other hand, if a particle goes the other way, the Russian roulette is applied, and the particle has a probability of I l /I h to continue with a modified weight of w × I l /I h . This technique increases the number of particles simulated in the region of interest for the same simulation time. Figure 28 illustrates the Geometric Splitting/Russian roulette variance reduction technique.
The variance reduction technique Geometry Splitting/Russian roulette was used in the BDSIM simulations to improve the secondary neutrons differential fluence computation accuracy. The scored volume in each wall of the vault was divided into eight cells of increasing importance values following the 2 k formula with k the cell index from one to eight. The importance value outside the cells is one. Figure 29 presents the importance cells in the BDSIM model.