Study of prompt gamma and neutron emission for real-time range verification in proton and carbon-ion therapy.

INTRODUCTION

Hadrontherapy is an advanced radiotherapy technique that allows for the effective treatment of deep‐seated and radio‐resistant tumors, as well as tumors located near critical organs. These advantages are related to the use of charged particles, such as protons and heavier ions (i.e., carbon‐ions), to deliver the dose with higher spatial precision compared to conventional photon therapy.1,2This feature is due to the physical characteristics of charged particles, which deposit most of their energy at the end of their path, resulting in the so called Bragg Peak (BP). This confined distribution allows to give less dose to the surrounding healthy tissues. Additionally, carbon‐ions exhibit higher linear energy transfer (LET) and relative biological effectiveness (RBE), offering further benefits for radio‐resistant tumors.3

However, to fully exploit these advantages, it is essential to ensure accurate in vivo range verification, as even small uncertainties in beam range can result in underdosing the tumor or, more critically, overdosing adjacent healthy tissues. These range uncertainties arise from several sources, including limitations in CT imaging, inaccuracies in the calibration from Hounsfield units (HU) to stopping power ratios (SPR),4,5,6,7patient setup errors, anatomical changes during the treatment course, and organ motion.8Additional contributions may arise from dose calculation approximations and beam delivery imperfections.9To mitigate these uncertainties, clinical protocols typically incorporate both range‐dependent and absolute safety margins during treatment planning.10Although this strategy increases robustness, it may partially compromise the spatial conformity of ion therapy by expanding the high‐dose region into surrounding healthy tissues.11

To mitigate these uncertainties, several in vivo range verification techniques have been developed; these are based on detecting secondary radiation emitted during beam irradiation, which is correlated with the beam path and dose deposition. One of the most promising methods is positron emission tomography (PET) imaging, which relies on the detection of$\beta^{+}$emitters produced along the beam path. PET systems can be implemented in off‐line, in‐room, or in‐beam configurations. While in‐beam PET reduces the impact of biological washout, PET‐based techniques still suffer from a progressive reduction in tracer activity and from practical limitations related to background radiation and acquisition constraints, which restrict their real‐time applicability.12,13,14,15Alternative approaches include ionoacoustic and bremsstrahlung‐based methods, which determine the BP position from acoustic signals or secondary photons generated during irradiation.16,17Another promising technique is Interaction Vertex Imaging (IVI), which tracks secondary charged particles emerging from nuclear interactions to reconstruct their origin and estimate the beam range within the patient.4Finally, it has been proposed that a small percentage of helium ions could be added to a carbon‐ion treatment beam for online treatment monitoring.

Other approaches are based on prompt gamma (PG) techniques, which exploit the almost instantaneous emission of nuclear de‐excitation photons directly correlated with the beam interaction positions. prompt gamma Spectroscopy (PGS) focuses on measuring the energy spectra of these photons to estimate the beam range with high specificity although it requires fast and low‐noise detection systems.18,19,20,21prompt gamma Timing (PGT) derives range information from the time‐of‐flight distribution of prompt gammas.22,23Finally, prompt gamma Imaging (PGI) reconstructs the spatial origin of prompt gammas using collimated or Compton‐based detectors, providing a visual representation of the beam path and enabling real‐time range verification.24,25,26

PGI, in particular, has emerged as one of the most effective approaches for real‐time monitoring of beam range and dose deposition.11,18,19,20Since its initial proposal in 2003,27it has evolved significantly thanks to advances in detection technologies. Systems based on slit and knife‐edge collimator have achieved one and two dimensional imaging with clinical accuracy down to approximately 1–2 mm in proton therapy.24,28,29Compton cameras, on the other hand, offer three dimensional reconstruction of gamma emission points.25,30More recently, time‐of‐flight (TOF) based PGI systems employing Cherenkov detectors have demonstrated millimetric sensitivity in pulsed beam conditions.11,31

PGI systems have also been successfully translated into clinical environments, where prompt gamma imaging has been used for in‐vivo monitoring of the proton beam range during patient treatments, demonstrating the capability to detect range deviations and anatomical changes under realistic treatment conditions.32,33,34

Despite these advances, PGI continues to face significant challenges due to the intense background from neutrons and their induced photons, especially in the context of heavy ion beams, where nuclear fragmentation and the associated neutron field are significantly stronger and can dominate the measured detector signal. These limitations highlight the need for improvements in high‐resolution TOF detection,31,35spectral filtering,36and Monte Carlo‐informed detector design.28,37

Monte Carlo simulations are fundamental in addressing these challenges, as they enable accurate modeling of beam–matter interactions, secondary radiation generation, and detector response. Widely used codes such asfluka,geant4, andmcnpallow for the optimization of PGI system geometry, material selection, and data acquisition strategies.28,35,38,39For instance, FLUKA simulations have been used to evaluate knife‐edge collimator‐based cameras for carbon‐ion therapy, demonstrating sub‐centimetric range sensitivity even in the presence of complex radiation backgrounds.28Similarly, Wang et al.36used FLUKA to study a pixelated LYSO detector with shielding and energy selection, showing that spectral discrimination can effectively suppress background and improve range resolution. Finally, a recent benchmark study comparing FLUKA, GEANT4, and MCNP6 revealed discrepancies in prompt gamma yield predictions and provided guidance for cross‐code calibration.37,40It should be noted that, from this study FLUKA showed the more reliable results with respect to the other codes.

While most PGI studies have so far focused on proton beams,11the present work extends the analysis to carbon‐ion therapy, addressing the additional challenges associated with nuclear fragmentation, neutron‐induced background and simulated detector response. Although the feasibility of prompt gamma imaging for carbon‐ions has been experimentally demonstrated, only a limited number of studies provide a systematic comparison between proton and carbon‐ion beams across clinically relevant energies within the same Monte Carlo framework.11,35In particular, few investigations simultaneously combine a detailed characterization of prompt gamma production in the phantom, an assessment of the neutron field in the treatment room, and a component‐resolved analysis of the detector signal, distinguishing between primary prompt gammas, neutron‐induced secondary gammas, and direct neutron contributions. For these reasons, a prototype design of a PGI system is investigated, employing a LYSO scintillation detector coupled with a knife‐edge collimator. The choice of a LYSO scintillation detector is driven by its favorable physical and technological properties for prompt gamma detection in a clinical environment. LYSO offers high density and effective atomic number, ensuring good detection efficiency in the MeV energy range, together with fast scintillation decay time and high light yield, which are essential for real‐time applications.41The detector geometry investigated in this work is inspired by recently proposed LYSO‐based prompt gamma imaging prototypes employing collimated configurations,41but it is enlarged with respect to existing implementations in order to allow for a detailed investigation of the spatial profiles and relative contributions of primary gammas, secondary gammas, and neutrons, allowing for a clearer interpretation of the signal. This work aims at linking the treatment room neutron field to detector effects and quantifies how each radiation component influences fall‐off–based range observables, particularly in the more challenging carbon‐ion case. In this context, the study characterizes the energy spectra and spatial distributions of prompt gammas together with background contributions from secondary gammas and neutrons, with the aim of assessing the system performance for real‐time range verification in clinical settings.

MATERIALS AND METHODS

Detector and simulated geometry

The simulated detection system reported in Figure1is based on a prototype concept currently under development. The system consists of a gamma camera, composed of a knife‐edge collimator and a pixelated scintillator detector, both aligned with the BP. The BP position is used as the reference point for the entire geometry and is always placed at the origin of the reference system, specifically at$z = 0$along the beam direction. Consequently, both the collimator and the detector are centered with respect to the BP and aligned with the cylindrical phantom. The collimator is placed 25 cm from the center of the phantom, while the detector is positioned 25 cm from the collimator, resulting in a total distance of 50 cm between the detector and the phantom.28

Geometry of the simulated setup implemented in FLUKA.

The collimator consists of two tungsten alloy blocks, each 15 cm long and 4 cm thick, arranged to form a 63opening angle with a 6 mm gap between the knife edges. The simulated detector consist of a pixelated LYSO scintillator with an area of 20 × 10 cm2and a thickness of 3 cm, segmented into 1 × 1 cm2transverse pixels. For the purposes of the present analysis, the recorded signals were then grouped into 1 cm‐thick longitudinal slabs in order to extract the one‐dimensional profile along the beam axis. This choice allows for a more detailed study of the detector signal and of how the 1D Bragg peak profile, evaluated along the beam axis, is influenced by the different signal components (i.e., primary prompt gammas, neutrons, and secondary gammas).

The simulated target is a cylinder with a height of 30 cm and a radius of 15 cm, made of ICRP tissue‐equivalent material42already implemented in FLUKA. Both the target geometry and material were selected to create a rather general setup, aiming to approximate, as closely as possible, the response of real tissue.

The entire simulated setup is positioned at the center of a room at 1.5 m above the floor. The room is 7$\times$7$m^{2}$in area, 4 m high, and it is surrounded by 2 m thick walls made of Portland cement (see Figure2). The material composition is reported in Table1, along with the other simulated materials. This configuration was chosen to resemble an approximated hadrontherapy treatment room.43

Simulated treatment room geometry and positions of the neutron fluence scoring spheres.

Table: Elemental composition of the simulated materials.

Monte Carlo Physics

Monte Carlo simulations were performed using the FLUKA code44(version 4‐4.1) with the PRECISIO setting, enabling high‐accuracy transport and interaction modelling.

Hadronic interactions are treated using different models depending on the projectile type and energy. For nucleons and pions, the PEANUT (PreEquilibrium Approach to Nuclear Thermalization) model is used up to a few GeV, accounting for intranuclear cascades, pre‐equilibrium emission, and de‐excitation processes.45For ion‐induced reactions, three models are employed depending on energy: the Boltzmann Master Equation (BME) model below$\sim$150 MeV/u, a modified version of rQMD‐2.4 between$\sim$150 MeV/u and$\sim$5 GeV/u, and DPMJET‐III at higher energies. Electromagnetic dissociation of ions is also included via virtual photon emission coupled with the PEANUT framework.45Transport thresholds in the adoptedsetting are$10^{- 5}$eV for neutrons, 10 keV for photons, and 100 keV for electrons, positrons, muons, and charged hadrons across all regions.

Simulation workflow and scored quantities

This study simulated two types of primary beam particles, namely protons and carbon‐ions. For each particle type, three different energies were selected, resulting in a total of six simulation sets. In particular, mono‐energetic pencil beams with energies of 60, 150, and 210 MeV for protons, and 120, 280, and 400 MeV/u for carbon‐ions, were employed. These energies were chosen to cover the energy range typically used in hadrontherapy and to obtain comparable BP depths in soft tissue42for both particles: 31, 156 and 280 mm for protons, and 35, 153 and 273 mm for carbon‐ions. The Bragg Peak (BP) depth was determined from the simulated depth‐dose distributions and defined as the position of the maximum deposited dose along the beam axis. For the simulations,$10^{9}$primary protons and$10^{8}$primary carbon‐ions were simulated. Table2summarizes the beam parameters and associated Bragg peak depths used in the simulations.

Table: Beam energies and corresponding Bragg Peak (BP) depths in soft tissue for the simulated proton and carbon‐ion beams.42

Various FLUKA cards and user‐defined routines were exploited for scoring data throughout the simulations. The analysis starts from the study of particles generated by the beam interaction within the phantom. Then, the neutron field produced as a consequence of neutron interactions with the room's walls and objects is analyzed. Finally, the signal scored in the simulated detector is analyzed in all its components: primary gammas, secondary gammas, and neutrons.

For the analysis of the interactions of the beam with the phantom material, an MGDRAW routine was implemented. This routine allowed the identification of the particle from which each photon originated (e.g., primary beam particles and their fragments, or other particles such as neutrons). The following classification has been adopted:

The MGDRAW routine allowed to save the phase space (i.e. energies, positions and directions) of photons and neutrons generated by the nuclear interactions of the primary particle with the soft tissue.

These data enabled to obtain the energy spectra and spatial profiles of the secondary radiation generated inside the phantom. Specifically, from the phase space of the primary gammas, the energy of each photon was retrieved to plot the energy spectra. The resulting spectra were binned using 10 keV intervals and are presented separately for each particle type. In addition, the spatial distribution of primary gamma production was evaluated by correlating photon energies with their production along the path of the primary particle. The profiles were computed for each beam energy with a spatial binning of 1 mm. Two approaches were applied: the first consists in analyzing the profiles using different energy selections, while the second focuses on the contribution of specific energy peaks. To further investigate the differences between beam types, the production yields of both primary gammas and neutrons were calculated.

In addition to the analysis of primary gammas, particular attention is given to the characterization of the neutron field, which represents one of the main challenges in a treatment room environment. Therefore, it is essential to characterize the neutron field generated within the treatment room during beam irradiation. To this purpose, several scoring volumes were placed at different positions in the room, and the neutron fluence (in neutrons per square centimeter) was evaluated using the FLUKAUSRTRACKcard.

Moreover, the signals scored in the LYSO were collected using theDETECTcard, which allows for the simulation of the behavior of a detector by scoring the energy depositions inside the sensitive volume. Thanks to the information provided by the MGDRAW routine on the parent particle responsible for each interaction, the detected signal could be subdivided into different components:

The scintillator detector was subdivided into slabs, and the signal was scored independently in each of them. By integrating the counts within each slab, a one‐dimensional profile of the detector signal as a function of position along the beam axis was reconstructed. This approach enables the evaluation of the spatial distribution of each radiation field component (primary gammas, secondary gammas, and neutrons) along the detector.

Fall‐off retrieval precision analysis

To evaluate the statistical precision in retrieving the distal fall‐off position, a procedure following the approach proposed byRoellinghoff et al.46was implemented. Starting from the simulated detector profiles expressed in counts per primary particle, a reference profile was obtained by scaling the signal to different numbers of primary particles$N_{prim}$. In this work, the analysis was performed for$N_{prim} = 10^{5} , 10^{6} , 10^{7} , 10^{8} , 10^{9}$and$10^{10}$.

For each value of$N_{prim}$, statistical fluctuations associated with counting statistics were simulated by generating 1000 noisy realizations of the signal through Poisson sampling of the expected counts in each bin. Each noisy profile was then compared with shifted versions of the reference profile. The reference curve was translated along the beam axis within a$\pm 2$cm interval using a shift step of 0.01 cm.

The similarity between the noisy profile and each shifted reference curve was quantified using the root mean square error (RMSE). The shift corresponding to the minimum RMSE was taken as the reconstructed fall‐off displacement. Repeating this procedure for many noisy realizations yields a distribution of reconstructed fall‐off shifts. The statistical uncertainty in the fall‐off determination is quantified from the width of this distribution, and the fall‐off retrieval precision (FRP) is defined as twice the standard deviation of the reconstructed shift distribution, corresponding to a$2 \sigma$confidence interval.

RESULTS

Radiation from Beam–Phantom interactions

At this stage, only the interactions between the beam particles and the phantom are presented; therefore, all the contributions coming from the room and the other objects are neglected. It should be noted that only the first photons generated by hadronic particle interactions were scored: photons resulting from the following interactions of primary photons or neutrons were not included in this initial analysis.

Energy spectra and depth profiles

The resulting spectra from the primary beam interaction are compared for each particle type in Figure3. The results are normalized to a primary proton or ion. The analysis was restricted to a 1–8 MeV window, since this interval is the most relevant for prompt gamma imaging. Above 8 MeV, the intensity of emitted gamma rays decreases significantly. Since the analysis is limited only to primary gamma, all the prompt gamma considered are originated from the de‐excitation of nuclei undergoing a nuclear inelastic collision. Therefore, the energy spectra exhibit peaks corresponding to the excited states of the nuclei present in the tissue‐equivalent material. Table3lists emission lines of the most abundant nuclei present in the phantom, providing an explanation for the main peaks observed in the energy spectra.

Primary gamma energy spectra for different proton and carbon‐ion energies.

Table: Main prompt gamma‐ray energies observed in the spectra of Figure3, associated with the de‐excitation of excited nuclei produced in beam–target interactions.20,21

Compared to protons, carbon‐ions exhibit a higher yield for gamma production near 2 MeV, in particular for the peak at 2.1 MeV. The most prominent peaks appear at 4.4 MeV, 5.2 MeV, and 6.1 MeV, corresponding to gamma lines derived from various nuclear reactions with12C,14N, and16O present in the phantom. In particular, the peaks correspond to the de‐excitation to the ground state of,,,,,, and.1

Figures4and5illustrate the depth profiles of primary gammas for different energy selections. The left panels show, for each simulation set, the primary gammas filtered into three intervals: no energy window, 1–8 MeV, and 3–7 MeV. The 3–7 MeV energy window was introduced to investigate a more realistic operating condition for a gamma detection system, characterized by reduced counting rates, while preserving the main prompt gamma features relevant for range verification. In addition, photons below approximately 3 MeV include a significant contribution from neutron‐induced processes and scattered radiation, which are weakly correlated with the Bragg peak. Restricting the analysis to the 3–7 MeV interval therefore improves the signal‐to‐background ratio by suppressing low‐energy neutron‐induced gamma contributions, as also discussed in literature byZarifi et al.35For clarity, the dose profile is scaled to match the maximum of the total profile. The fall‐off of the profiles of prompt gamma emission shows a clear correlation with the dose distribution and thus with the Bragg peak. The profiles for carbon‐ions appear narrower and with a sharper fall‐off than those for protons, which show a broader distribution. In the case of carbon‐ions, a tail beyond the Bragg peak is visible both in the dose distribution and the prompt gamma emission.

Gamma profiles and detailed distributions within the phantom for protons at various energies. The legend shown in (a) and (b) applies to all figures.

Gamma profiles and detailed distributions within the phantom for carbon‐ions at various energies. The legend shown in (a) and (b) applies to all figures.

On the right side of Figures4and5the 3–7 MeV interval is further subdivided into smaller intervals corresponding to the main peaks observed in the energy spectra. This enables a more detailed analysis of the spatial profiles associated with specific nuclear de‐excitation lines. When energy filtering is applied, the fall‐off distribution changes shape due to the different thresholds of the reactions involved in prompt gamma generation.

Table4provides the yields per primary particle of photons and neutrons. The values for prompt gamma are reported for three different energy windows: no selection, 1–8 MeV, and 3–7 MeV.

Table: Production yield (counts per primary particle) for different projectiles and energies.

Neutron field

For each simulation set, the neutron spectra collected in various locations of the room are plotted in Figure6. For visualization purposes, the spectrum corresponding to position E is scaled differently, with its corresponding scale displayed on the right side of the plot. This is because this specific scoring volume is positioned directly in front of the phantom, facing the incoming beam. Due to the strong forward‐peaked distribution of fast neutrons, the fluence rate in this location is significantly higher, making it difficult to visualize the spectra in other regions of the room on the same scale.

Neutron spectra scored in various locations around the room (see Figure2). Due to the significantly higher fluence measured at location E, this spectrum is shown with a separate scale reported on the right.

The neutron spectra reveal the presence of thermal ($<$0.5 eV), epithermal (0.5 eV–0.1 MeV), and fast ($>$0.1 MeV) components. The spectra extend to higher energies with increasing beam energy, and the maximum neutron energies are higher for carbon‐ions compared to protons. While the neutron component below 10–20 MeV remains relatively constant throughout the room, the fluence of neutrons at higher energies is strongly dependent on the scoring position. Above this energy threshold, the spectra recorded at positions B, C, and E exhibit additional peaks correlated to the intranuclear cascade, with fluence increasing as the position approaches the beam axis in the distal region.

Detector signals

The detector signals for each simulation, subdivided into the radiation field components, are presented in Figure7. The plots are shown within the 3–7 MeV energy window to better visualize the data and focus on the most relevant energy region. The signal component due to neutrons increases with beam energy and becomes dominant in the carbon‐ion cases. Additional peaks due only to neutron interactions with the detector are visible at 5.8 MeV, 6.3 MeV, 6.5 MeV and 6.6 MeV.

Energy spectra detected in the crystal region for different proton and carbon‐ion energies.The legend shown in (a) applies to all figures.

1‐D profiles

Figure8presents the resulting 1‐D profile obtained by including only prompt gammas in the 3–7 MeV energy window. The fall‐off related to the Bragg peak is more clearly visible in proton profiles than in carbon‐ion profiles, especially at lower energies. For carbon‐ions, the total profile exhibits a less pronounced fall‐off and a higher contribution from secondary components. The neutron and secondary gamma components form a relatively flat background. While secondary gammas are almost uniformly distributed, the neutron component shows an increase in the proximal region in respect to the beam.

1‐D profiles of the detected signal in the crystal region for different proton and carbon‐ion energies (3–7 MeV energy window). The legend shown in (a) applies to all figures.

The profiles vary across different orders of magnitude, particularly due to the different yields of the radiation field components in the case of carbon‐ion beams. Therefore, Figure9presents the six profiles of primary gamma radiation (Figure9a) and the six profiles of the total detected signal (Figure9b) scaled for better comparison. As the absolute quantities have already been shown in the previous figures, the profiles in Figure9are displayed after subtracting their minimum value, aligning all baselines to zero and thereby facilitating shape comparison. These plots are useful to highlight the slope induced by the BP, which would otherwise be less evident because of the large differences in component yields.

Primary gamma and total profiles detected in the crystal region.

To quantify the statistical precision in retrieving the distal fall‐off position, the fall‐off retrieval precision (FRP) was evaluated as a function of the number of primary particles for the different beam energies. The results are shown in Figure10for proton and carbon‐ion beams.

Fall‐off retrieval precision (FRP) as a function of the number of primary particles at different energies.

In all cases, the FRP decreases with increasing number of primary particles, reflecting the improvement in statistical precision with increasing counting statistics. For proton beams, the FRP shown in Figure10ais systematically lower for primary gammas than for the total detected profile. In addition, the FRP tends to increase with beam energy, indicating a reduced precision in determining the fall‐off position at higher proton energies.

A similar trend is observed for carbon‐ion beams, as shown in Figure10b, where the FRP decreases with increasing number of primaries and generally increases with beam energy. However, for the highest energy configuration (400 MeV/u), the FRP obtained from the total signal is smaller than that obtained from primary gammas. The origin of this behavior is discussed in the following section.

DISCUSSION

The results of the Monte Carlo simulations provide different insights into the behavior of prompt and secondary radiation generated during hadrontherapy treatments, as well as their effects on a system based on a LYSO crystal.

Energy spectra and spatial profiles

The energy spectra of the primary gammas (Figure3) exhibit distinct peaks corresponding to the de‐excitation process of the nuclei present in the tissue‐equivalent material. The most relevant lines for prompt gamma imaging are the 4.4 MeV line from12C, the 5.2 MeV line from15O, and the 6.1 MeV line from15O and16O.

Comparing the two particle types, carbon‐ions produce a higher prompt gamma yield than protons. This is primarily because carbon‐ions are heavier and more highly charged, increasing their probability to induce nuclear inelastic reactions along their path. Additionally, carbon‐ion spectra show extra peaks around 2 MeV due to prompt gamma emissions due to nuclear fragments. Indeed, the emission lines at approximately 2.0, 2.1, and 2.3 MeV correspond to11C,11B, and15O, respectively. In particular, a large increase is observed for the 2.1 MeV line associated with11B compared to the proton spectra.

After the energy spectra were analyzed, the spatial distribution of prompt gammas was investigated. Profiles presented in Figures4and5highlight the close correlation between prompt gamma emission and dose deposition for protons and carbon‐ions, respectively. As expected, the prompt gamma profile follows the shape of the dose distribution along the beam path and presents a fall‐off corresponding to the BP, confirming its suitability for range verification. While the overall shape of the profiles is comparable, some differences between protons and carbon‐ions can be noticed. In particular, the prompt gamma profile is broader for protons than for carbon‐ions, due to their lower mass and lower LET. Furthermore, while the proton profile ends after the fall‐off, the carbon‐ions profile exhibits a tail beyond the BP due to prompt gamma induced by nuclear fragments.

The filtered profiles in Figures4,5(left) are essential to select the prompt gamma energy range with the highest yield, while optimizing the signal‐to‐background ratio. Here, energy filtering refers to the application of energy thresholds or energy windows at the data acquisition or analysis level. The results show that applying energy filtering does not compromise the spatial correlation with the BP: for all three energy windows, indeed, the fall‐off remains unchanged within the spatial and statistical resolution of the analysis. Notably, for the 400 MeV/u beam, energy filtering improves the profile by reducing the contribution of low‐energy photons in the initial region of the phantom. The main drawback is the reduction in absolute yield, which, however, can be acceptable in order to improve the acquisition system.

To quantify the stability of the fall‐off position, the depth at which the prompt gamma profile decreases to 50% of its maximum value ($z_{50}$) was extracted from the profiles obtained for the unfiltered spectrum and for the different energy windows. The differences between the corresponding$z_{50}$values were found to be below 0.4 mm for all beam configurations. This value is smaller than the adopted spatial binning (1 mm) and therefore within the resolution of the present analysis. Consequently, the position of the fall‐off can be considered unchanged when applying the different energy selections.

The production yields reported in Table4quantitatively confirm the qualitative trends observed in the spectral and depth profile analyses. For completeness, neutron yields are also reported. For proton beams, the neutron yield per primary particle is of the same order of magnitude as the prompt gamma yield. In contrast, for carbon‐ions the neutron production becomes dominant. In particular, for the 400 MeV/u beam, the neutron yield per primary particle is approximately seven times higher than the corresponding prompt gamma yield. This marked increase in the neutron‐to‐prompt gamma ratio highlights the progressively more unfavorable signal‐to‐background conditions in carbon‐ion therapy and underlines the importance of properly accounting for neutron‐induced contributions in range monitoring systems.

Possible mitigation strategies include the use of neutron absorbers, such as highly borated materials, to reduce the thermal and epithermal neutron field, which is expected to contribute significantly to the detected background signal through neutron capture processes in the detector materials. Preliminary observations suggest that thermal neutrons play a relevant role and this aspect will be investigated more systematically in future studies. In addition, discrimination techniques based on pulse‐shape analysis (PSD) may allow the separation of neutron‐induced signals from prompt gamma interactions when using adequate scintillation crystals.

In addition to the considerations discussed above, the choice of the 3–7 MeV energy window is further supported by the presence of a high neutron field in the treatment room. Indeed, this window excludes lower‐energy gammas that are more likely to originate from scattering and secondary processes occurring within the room walls and surrounding materials. Specifically, this energy window avoids the 2.2 MeV gamma line produced by neutron capture on hydrogen, which contributes significantly to the counts but does not provide any spatial correlation with the BP.

A quantitative evaluation confirms this behavior. In the 0–3 MeV interval, secondary gammas account for approximately 14% –28% of the total detected signal depending on the beam configuration. This fraction decreases to about 4–8% in the 1–3 MeV region and further drops to approximately 1–3% in the 3–7 MeV window. Therefore, selecting the 3–7 MeV interval significantly suppresses the contribution of secondary gammas while preserving the prompt gamma signal correlated with the Bragg peak.

Neutron field

As previously mentioned, the analysis of the neutron field is crucial when considering the implementation of a detection system inside a treatment room. The energy spectra of the neutron field, shown in Figure6, show the presence of thermal, epithermal, and fast neutron components. Thermal and epithermal neutrons are typically present in all neutron spectra, as they are produced through moderation processes in the surrounding materials, especially the treatment room walls. In contrast, fast neutrons, particularly those with energies above 10–20 MeV, are characteristic of interactions induced by high‐energy particle beams. These neutrons are produced from intranuclear cascade reactions, which are triggered when the projectile interacts inelastically with the nuclei in the target material. The spectrum of high‐energy neutrons is directly related to the primary beam, particularly to its energy per nucleon. As an example, for the most energetic beam simulated, the 400 MeV/u carbon‐ion beam, the neutron spectrum extends up to nearly 1 GeV.

As expected from the yields listed in Table4, the neutron fluence in the energy spectra increases with the energy of the beam and is significantly higher for carbon‐ions than for protons. In fact, for the same penetration depth, carbon‐ions produce at least one order of magnitude more neutrons than protons. This difference is related to the higher nuclear cross‐section and fragmentation probability.

However, while the thermal and epithermal neutron components remain relatively uniform across the room, as shown in Figure6, the fluence of fast neutrons varies considerably depending on the scoring position. Neutrons with energies above 10–20 MeV are strongly forward‐peaked along the beam direction, and this effect is clearly visible at positions B, C, and E, which show progressively increasing fluence. As shown in Figure2, positions A and D are behind the fast neutron emission direction and thus record limited high‐energy neutron fluence.

This analysis underlines the importance of characterizing the various components of the radiation field. Understanding the energy composition of the neutron field is crucial, particularly when evaluating the signal‐to‐background ratio of detectors. Knowledge of the neutron energy spectrum allows for the design of optimized shielding to reduce the background and improve prompt gamma detection. Although the very fast component of the neutron field remains challenging to suppress due to its high energy and penetration power, significant improvements could be achieved by specifically addressing the thermal and epithermal components, which still represent a high fraction of the total field.

Detector signal

Figure7provides the analysis of the signal recorded in the scintillator detector, separated into its main components, offering a valuable insight into the influence of the radiation field in the treatment room. The neutron contribution increases with beam energy, as expected.

Consequently, secondary gammas, mainly produced by neutron interactions, increase along with them, although not at the same rate. In fact, neutron‐induced gammas are often emitted at lower energies, most notably around 2.2 MeV, outside the 3–7 MeV window considered. Additionally, secondary gammas that reach the detector after scattering within the room typically undergo further energy degradation. Moreover, the secondary gammas detected within this energy window exhibit the same spectral peaks as the primary prompt gammas discussed in Section3.1, since high‐energy neutrons can trigger the intranuclear cascade as the primary particles.

On the other hand, neutrons mainly produce a continuous signal in the detector due to various scattering processes. However, additional peaks exclusively due to neutrons are visible above 5 MeV (i.e., 5.8, 6.3, 6.5, 6.6 MeV). These peaks result from gamma emissions following neutron capture on lutetium, one of the components of the detector material. Lutetium consists primarily of two isotopes:175Lu, with a natural abundance of 97.4% and a neutron capture cross section of 23 barns, and176Lu, with 2.6% abundance and a significantly higher cross section of 2090 barns.

Even if the energy deposited comes from gamma rays emitted during the capture process, the signal is attributed to neutrons in the signal plots, as the reaction takes place inside the detector itself. As previously noted, the signal associated only to neutrons in Figure7can therefore be largely attributed to the thermal component, as suggested by these distinct peaks, which are solely associated with neutron capture events inside the detector.

The realistic energy resolution of the scintillator (i.e., 10%) was not applied in order to preserve the visibility of narrow spectral structures (e.g., gamma lines from neutron capture in lutetium) that would otherwise merge with broader prompt gamma features.

To study the spatial distribution of the detected signal, the profiles detected are presented in Figure8. This comparison highlights that the slope related to the BP is more noticeable in the profiles associated to proton beams than to carbon‐ions, and it is more evident for lower energies. This is consistent with all the observations made in the previous sections on the relation and distribution of the various components.

For carbon‐ions, the higher initial dose and extended tail due to fragmentation reduce the contrast between the BP and the non‐relevant regions of the profile, making the fall‐off less pronounced with respect to protons. Moreover, the tail resulting from nuclear fragmentation causes a forward shift of the 1D profile, as shown in Figure9a. Furthermore, the neutron and secondary gamma components contribute a relatively uniform background that further flattens the signal. This behavior is consistent with the ratio between neutron and primary gamma yields reported in Table4. The higher neutron production leads to a flatter profile, making the fall‐off slope expected at the BP less evident. A small difference can be noted between the neutron and secondary gamma contributions to the profile: while secondary gammas tend to form a nearly uniform background, the neutron component shows a slight increase in the entrance region of the phantom, likely due to the fact that neutrons are primarily produced before the BP, where the beam energy is higher. These differences in the shape of the detected profiles also affect the statistical precision with which the fall‐off position can be determined. To quantify this effect, the fall‐off retrieval precision (FRP) was evaluated as a function of the number of primary particles.

The comparison between primary gammas and the total detected profiles highlights the impact of secondary radiation on the precision of the fall‐off determination. For proton beams, as shown in Figure10a, the FRP obtained using only primary gammas is systematically lower than that obtained from the total detected signal. This behavior is expected, since the total signal includes contributions from neutrons and secondary gamma radiation, which introduce additional fluctuations that are not directly correlated with the Bragg peak position.

The FRP also shows a dependence on beam energy. In particular, the precision tends to degrade with increasing proton energy. This trend can be related to the broader spatial distribution of prompt gamma emission at higher energies, which results in a less steep distal fall‐off of the detected profiles, as previously observed in Figure8. As the fall‐off becomes less pronounced, the determination of its position becomes more sensitive to statistical fluctuations in the detected signal.For carbon‐ion beams, a similar trend is observed (Figure10b), where the FRP decreases with increasing number of primaries and generally increases with beam energy.

It should be noted that the FRP values obtained for proton and carbon‐ion beams are of the same order of magnitude when evaluated for a fixed number of primary particles. However, this comparison does not correspond to equivalent clinical irradiation conditions. Due to the higher relative biological effectiveness of carbon ions and to the dependence of energy loss on$1 / Z^{2}$, resulting in a higher stopping power and a higher energy deposition, a smaller number of primary particles is required to deliver the same dose to the patient. Consequently, the number of detected prompt gammas per treatment spot is expected to be lower in carbon‐ion therapy. In typical clinical conditions, the number of protons per spot can reach approximately$10^{8}$–$10^{9}$, whereas the number of carbon ions is usually of the order of$10^{6}$–$10^{7}$. As a consequence, when comparing the FRP at clinically relevant intensities, proton beams achieve sub‐millimetric precision (approximately 0.2–0.7 mm), while carbon‐ion beams exhibit a reduced precision, typically in the millimetric range (approximately 2–6 mm), due to the lower number of detected prompt gammas per spot. Therefore, although the FRP obtained in the simulations is comparable for the same number of primaries, the achievable statistical precision in realistic treatment conditions may be more challenging for carbon‐ion beams due to the reduced number of incident particles.

A comparison between the lowest‐energy configurations further highlights the role of secondary radiation. In the case of the 60 MeV proton beam, the FRP obtained from primary and total signals is almost identical, indicating a negligible contribution from secondary components. In contrast, for the 120 MeV/u carbon‐ion beam the difference between the primary and total FRP is significantly larger, reflecting the stronger influence of neutron‐induced and secondary gamma radiation on the detected profiles. However, for the highest energy configuration (400 MeV/u), the FRP obtained from the total signal is slightly smaller than that obtained using only primary gammas. This effect is attributed to the geometrical limitation of the phantom used in the simulations. In this configuration, the Bragg peak is located close to the distal edge of the phantom, so that particles and nuclear fragments traveling beyond this region are not fully contained within the phantom volume. As a consequence, a fraction of forward‐going radiation, including fragments, neutrons, and the corresponding secondary gammas, is lost. This effect artificially sharpens the distal fall‐off of the total detected profile, producing an apparently improved FRP that does not correspond to a real improvement in the range sensitivity. It should be noted that this effect would be present also in a real patient depending on the irradiation geometry.

It should be underlined that this study adopts a simplified geometry to isolate the physical drivers of prompt‐gamma and neutron contributions. In fact, the phantom is homogeneous soft tissue and does not represent patient heterogeneities. Moreover, the irradiation conditions employ mono‐energetic pencil beams and do not model full treatment plans. Finally, an intrinsic scintillator energy resolution was not applied to preserve narrow spectral structures for analyzing its main components; in clinical detectors, these features may be broadened and partially merge with prompt‐gamma lines, although the underlying neutron‐capture mechanism remains relevant for background formation.

Nevertheless, the component detector analysis provides practical guidance for any future PGI system optimization. In particular, especially for carbon‐ion beams, restricting the analysis to the 3–7 MeV window improves signal‐to‐background conditions by suppressing neutron‐induced contributions while preserving the Bragg‐peak‐correlated prompt‐gamma component. However, the strong neutron field and neutron‐capture processes in the scintillation crystal introduce additional background and flatten the range‐sensitive fall‐off, emphasizing the need for moderated‐neutron mitigation and Monte Carlo informed interpretation strategies in carbon‐ion PGI. It should be noted that, all the analysis performed in this work can be extended to any scintillation crystal.

CONCLUSION

This work presented a Monte Carlo investigation of prompt gamma emission and neutron background in proton and carbon‐ion therapy, performed using the FLUKA code and a simulated prompt gamma imaging system based on a knife‐edge collimator coupled to a LYSO scintillator detector. The simulations aimed at performing a systematic comparison between proton and carbon‐ion beams at clinically relevant energies. The results confirm the spatial correlation between prompt gamma emission and the Bragg peak for both particle types, supporting PGI for range verification. However, carbon‐ion beams produce a stronger neutron background, which increases with energy and reduces the contrast of the distal fall‐off. The analysis of the simulated detector supports the choice of the 3–7 MeV energy window, where secondary gamma contributions are strongly suppressed while preserving the Bragg‐peak‐correlated prompt gamma signal, thus improving the signal‐to‐background conditions.

The statistical precision in retrieving the distal fall‐off position was evaluated through FRP as a function of the number of primary particles. The comparison between primary and total detected profiles shows that secondary radiation, mainly neutrons and secondary gammas from, for example, the concrete walls, degrades the precision of the fall‐off determination by introducing fluctuations not correlated with the Bragg peak. Overall, these results provide a characterization of the radiation components affecting prompt gamma detection in hadrontherapy and highlight the main physical mechanisms that influence the accuracy of range verification. These findings provide guidance for the design and optimization of prompt gamma detection systems for real‐time range monitoring, especially in the more challenging radiation environment associated with carbon‐ion therapy.

CONFLICT OF INTEREST STATEMENT

The authors declare no conflicts of interest.

ACKNOWLEDGMENTS

This work was supported by the REALPATH project (Italian MUR), grant n. FISA‐2022‐00779.

Open access publishing facilitated by Politecnico di Milano, as part of the Wiley ‐ CRUI‐CARE agreement.

References

1. Durante M, Paganetti H. Nuclear physics in particle therapy: a review. Rep Prog Phys. 2016;79(9):096702. doi.org/10.1088/0034-4885/79/9/096702
2. Schardt D, Elsässer T, Schulz‐Ertner D. Heavy‐ion tumor therapy: Physical and radiobiological benefits. Rev Mod Phys. 2010;82:383–425.
3. Tommasino F, Scifoni E, Durante M. New ions for therapy. Int J Part Ther. 2015;2(3):428–438. doi.org/10.14338/IJPT-15-00027.1
4. Henriquet P, Testa E, Chevallier M, et al. Interaction vertex imaging (IVI) for carbon ion therapy monitoring: a feasibility study. Phys Med Biol. 2012;57(14):4655. doi.org/10.1088/0031-9155/57/14/4655
5. Lomax AJ. Intensity modulated proton therapy and its sensitivity to treatment uncertainties 1: the potential effects of calculational uncertainties. Phys Med Biol. 2008;53(4):1027. doi.org/10.1088/0031-9155/53/4/014
6. Paganetti H, Beltran C, Both S, et al. Roadmap: proton therapy physics and biology. Phys Med Biol. 2021;66(5):05RM01. doi.org/10.1088/1361-6560/abcd16
7. Yang M, Zhu XR, Park PC, et al. Comprehensive analysis of proton range uncertainties related to patient stopping‐power‐ratio estimation using the stoichiometric calibration. Phys Med Biol. 2012;57(13):4095. doi.org/10.1088/0031-9155/57/13/4095
8. Knopf AC, Lomax A. In vivo proton range verification: a review. Phys Med Biol. 2013;58(15):R131. doi.org/10.1088/0031-9155/58/15/R131
9. Paganetti H. Range uncertainties in proton therapy and the role of Monte Carlo simulations. Phys Med Biol. 2012;57(11):R99. doi.org/10.1088/0031-9155/57/11/R99
10. Unkelbach J, Chan TCY, Bortfeld T. Accounting for range uncertainties in the optimization of intensity modulated proton therapy. Phys Med Biol. 2007;52(10):2755. doi.org/10.1088/0031-9155/52/10/009
11. Pinto M. Prompt‐gamma imaging in particle therapy. Eur Phys J Plus. 2024;139(10):884.
12. Parodi K, Yamaya T, Moskal P. Experience and new prospects of PET imaging for ion beam therapy monitoring. Z Med Phys. 2023;33(1):22–34. doi.org/10.1016/j.zemedi.2022.11.001
13. Oelfke U, Lam GKY, Atkins MS. Proton dose monitoring with PET: quantitative studies in Lucite. Phys Med Biol. 1996;41(1):177. doi.org/10.1088/0031-9155/41/1/013
14. Shakirin G, Braess H, Fiedler F, et al. Implementation and workflow for PET monitoring of therapeutic ion irradiation: a comparison of in‐beam, in‐room, and off‐line techniques. Phys Med Biol. 2011;56(5):1281. doi.org/10.1088/0031-9155/56/5/004
15. Ferrero V, Fiorina E, Morrocchi M, et al. Online proton therapy monitoring: clinical test of a Silicon‐photodetector‐based in‐beam PET. Sci Rep. 2018;8(1):4100. doi.org/10.1038/s41598-018-22325-6
16. Assmann W, Kellnberger S, Reinhardt S, et al. Ionoacoustic characterization of the proton Bragg peak with submillimeter accuracy. J Med Phys. 2015;42(2):567‐574. doi.org/10.1118/1.4905047
17. Yamaguchi M, Nagao Y, Ando K, et al. Secondary‐electron‐bremsstrahlung imaging for proton therapy. Nucl Instrum Methods Phys Res Sect A Accel Spectrom Detect Assoc Equip. 2016;833:199‐207.
18. Testa E, Bajard M, Chevallier M, et al. Monitoring the Bragg peak location of 73MeV/u carbon ions by means of prompt gamma ray measurements. Appl Phys Lett. 2008;93(9):093506.
19. Min CH, Kim CH, Youn MY, Kim JW. Prompt gamma measurements for locating the dose falloff region in the proton therapy. Appl Phys Lett. 2006;89(18):183517.
20. Dal Bello R, Magalhaes Martins P, Brons S, et al. Prompt gamma spectroscopy for absolute range verification of 12C ions at synchrotron‐based facilities. Phys Med Biol. 2020;65(9):095010. doi.org/10.1088/1361-6560/ab7973
21. Verburg JM, Seco J. Proton range verification through prompt gamma‐ray spectroscopy. Phys Med Biol. 2014;59(23):7089. doi.org/10.1088/0031-9155/59/23/7089
22. Golnik C, Hueso‐González F, Müller A, et al. Range assessment in particle therapy based on prompt gamma‐ray timing measurements. Phys Med Biol. 2014;59(18):5399. doi.org/10.1088/0031-9155/59/18/5399
23. Werner T, Berthold J, Hueso‐González F, et al. Processing of prompt gamma‐ray timing data for proton range measurements at a clinical beam delivery. Phys Med Biol. 2019;64(10):105023. doi.org/10.1088/1361-6560/ab176d
24. Smeets J, Roellinghoff F, Prieels D, et al. Prompt gamma imaging with a slit camera for real‐time range control in proton therapy. Phys Med Biol. 2012;57(11):3371. doi.org/10.1088/0031-9155/57/11/3371
25. Peterson SW, Robertson D, Polf J. Optimizing a three‐stage Compton camera for measuring prompt gamma rays emitted during proton radiotherapy. Phys Med Biol. 2010;55(22):6841. doi.org/10.1088/0031-9155/55/22/015
26. Krimmer J, Dauvergne D, Letang J, Testa E. Prompt‐gamma monitoring in hadrontherapy: a review. Nucl Instrum Methods Phys Res Sect A Accel Spectrom Detect Assoc Equip. 2018;878:58–73.
27. Stichelbaut F, Jongen Y. Verification of the proton beam position in the patient by the detection of prompt gamma‐rays emission. In: Proceedings of 39th meeting of the particle therapy co‐operative group Vol 16; 2003.
28. Missaglia A, Bourkadi‐Idrissi A, Casamichiela F, et al. Prompt‐gamma fall‐off estimation with C‐ion irradiation at clinical energies, using a knife‐edge slit camera: a Monte Carlo study. Physica Med. 2023;107:102554. doi.org/10.1016/j.ejmp.2023.102554
29. Malekzadeh E, Rajabi H, Tajik‐Mansoury MA, Sabouri P, Fiorina E, Kalantari F. Design and performance evaluation of a slit‐slat camera for 2D prompt gamma imaging in proton therapy monitoring: a Monte Carlo simulation study. Med Phys. 2023;50(6):3701–3718. doi.org/10.1002/mp.16259
30. Richard MH, Chevallier M, Dauvergne D, et al. Design guidelines for a double scattering compton camera for prompt‐γimaging during ion beam therapy: a Monte Carlo simulation study. IEEE Trans Nucl Sci. 2011;58(1):87–94.
31. Jacquet M, Ansari S, Gallin‐Martel ML, et al. A high sensitivity Cherenkov detector for prompt gamma timing and time imaging. Sci Rep. 2023;13(1):3609. doi.org/10.1038/s41598-023-30712-x
32. Richter C, Pausch G, Barczyk S, et al. First clinical application of a prompt gamma based in vivo proton range verification system. Radiother Oncol. 2016;118(2):232–237. doi.org/10.1016/j.radonc.2016.01.004
33. Berthold J, Khamfongkhruea C, Petzoldt J, et al. First‐in‐human validation of CT‐based proton range prediction using prompt gamma imaging in prostate cancer treatments. Int J Radiat Oncol Biol Phys. 2021;111(4):1033–1043. doi.org/10.1016/j.ijrobp.2021.06.036
34. Berthold J, Pietsch J, Piplack N, et al. Detectability of anatomical changes with prompt‐gamma imaging: first systematic evaluation of clinical application during prostate‐cancer proton therapy. Int J Radiat Oncol Biol Phys. 2023;117(3):718–729. doi.org/10.1016/j.ijrobp.2023.05.002
35. Zarifi M, Guatelli S, Qi Y, Bolst D, Prokopovich D, Rosenfeld A. Characterization of prompt gamma ray emission for in vivo range verification in particle therapy: a simulation study. Physica Med. 2019;62:20–32. doi.org/10.1016/j.ejmp.2019.04.023
36. Wang JL. Pixelated prompt gamma imaging detector for online measurement of proton beam: Monte Carlo feasibility study by FLUKA. Radiat Detect Technol Methods. 2018;2(1):4.
37. Freitas H, Nobakht E, Grüner F, Seco J. A comparative analysis of GEANT4, MCNP6 and FLUKA on proton‐induced gamma‐ray simulation. Med Phys. 2025;52(6):4862–4870. doi.org/10.1002/mp.17754
38. El Kanawati W, Létang JM, Dauvergne D, et al. Monte Carlo simulation of prompt gamma‐ray emission in proton therapy using a specific track length estimator. Phys Med Biol. 2015;60(20):8067. doi.org/10.1088/0031-9155/60/20/8067
39. Skjerdal K, Kögler T, Lionheart WR, Ytre‐Hauge KS, Meric I. Range verification in proton therapy with computational prompt gamma‐ray spectroscopy: a Monte Carlo study. Radiat Phys Chem. 2026;238:113096.
40. Wrońska A, Kasper J, Ahmed AA, et al. Prompt‐gamma emission in GEANT4 revisited and confronted with experiment. Physica Med. 2021;88:250–261. doi.org/10.1016/j.ejmp.2021.07.018
41. Idrissi AB, Borghi G, Caracciolo A, et al. First experimental verification of prompt gamma imaging with carbon ion irradiation. Sci Rep. 2024;14(1):25750. doi.org/10.1038/s41598-024-72870-6
42. ICRP . Conversion coefficients for radiological protection quantities for external radiation exposures . Technical report, ICRP Publication 116, Ann. ICRP 2010;40(2‐5). doi.org/10.1016/j.icrp.2011.10.001
43. Clasie BM, Letourneau D, Schwarz M, Seuntjens J, Maughan RL. Proton therapy equipment installation, upgrades, and building design. Pract Radiat Oncol. 2024;14(3):e249‐e254. doi.org/10.1016/j.prro.2023.09.011
44. FLUKA collaboration . FLUKA 4‐4.4 Manual . CERN; 2024.https://fluka.cern
45. Battistoni G, Muraro S, Sala PR, et al. The FLUKA code: description and benchmarking. In: Hadronic Shower Simulition Workshop . Vol 896. American Institute of Physics Conference Proceedings; 2007:31‐49.
46. Roellinghoff F, Benilov A, Dauvergne D, et al. Real‐time proton beam range monitoring by means of prompt‐gamma detection with a collimated camera. Phys Med Biol. 2014;59:1327–38. doi.org/10.1088/0031-9155/59/5/1327

---

Republished from the open web under CC-BY. Authors: Bellotti ES, Mazzucconi D, Bortot D, Pola A, Fiorini C, Urban K, Agosteo S. Read the original.