Secondary neutron spectrum from 250-MeV passively scattered proton therapy: Measurement with an extended-range Bonner sphere system

Purpose: Secondary neutrons are an unavoidable consequence of proton therapy. While the neutron dose is low compared to the primary proton dose, its presence and contribution to the patient dose is nonetheless important. The most detailed information on neutrons includes an evaluation of the neutron spectrum. However, the vast majority of the literature that has reported secondary neutron spectra in proton therapy is based on computational methods rather than measurements. This is largely due to the inherent limitations in the majority of neutron detectors, which are either not suitable for spectral measurements or have limited response at energies greater than 20 MeV. Therefore, the primary objective of the present study was to measure a secondary neutron spectrum from a proton therapy beam using a spectrometer that is sensitive to neutron energies over the entire neutron energy spectrum. Methods: The authors measured the secondary neutron spectrum from a 250-MeV passively scattered proton beam in air at a distance of 100 cm laterally from isocenter using an extended-range Bonner sphere (ERBS) measurement system. Ambient dose equivalent H*(10) was calculated using measured ﬂuence and ﬂuence-to-ambient dose equivalent conversion coefﬁcients. Results: The neutron ﬂuence spectrum had a high-energy direct neutron peak, an evaporation peak, a thermal peak, and an intermediate energy continuum between the thermal and evaporation peaks. The H*(10) was dominated by the neutrons in the evaporation peak because of both their high abundance and the large quality conversion coefﬁcients in that energy interval. The H*(10) 100 cm laterally from isocenter was 1.6 mSv per proton Gy (to isocenter). Approximately 35% of the dose equivalent was from neutrons with energies ≥ 20 MeV. Conclusions: The authors measured a neutron spectrum for external neutrons generated by a 250-MeV proton beam using an ERBS measurement system that was sensitive to neutrons over the entire energy range being measured, i.e., thermal to 250 MeV. The authors used the neutron ﬂuence spectrum to demonstrate experimentally the contribution of neutrons with different energies to the total dose equivalent and in particular the contribution of high-energy neutrons ( ≥ 20 MeV). These are valuable reference data that can be directly compared with Monte Carlo and experimental data in the literature. © 2014 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported License. [http://dx.doi.org/10.1118/1.4892929]


INTRODUCTION
Generation of secondary neutrons is an unavoidable consequence of proton therapy, as it is for high-energy photon therapy (>10 MV). For both photon and proton therapy, the secondary neutrons have a wide energy range, from thermal up to the maximum beam energy; for photon therapy, the maximum neutron is 18-20 MeV (most probable energy 1 MeV with an exceedingly small fraction >10 MeV), but for proton therapy the maximum neutron energy can be as high as the proton beam energy, which is as high as 250 MeV at some centers.
Secondary neutron spectra from proton therapy have three pronounced peaks. 1 The first peak is primarily comprised of direct neutrons and ranges in energy from approximately 10 MeV to the maximum proton energy, which may be as high as 250 MeV. The second is due to evaporation processes and ranges in energy from approximately 1 keV to approximately 15 MeV. The third is a thermal peak due to neutrons that have scattered around the room and returned to the measurement location, i.e., room return. Neutrons also exist between these peaks in the form of an intermediate energy continuum.
The secondary neutron spectra from high-energy photon therapy have been studied extensively using both Monte Carlo simulations and measurements. 2 In many of these studies, spectral measurements have been used to validate the computational models. [3][4][5][6][7][8] In recent years, as the number of proton centers in the United States and worldwide has grown, so too has the literature on secondary neutrons from proton therapy. However, in the vast majority of this literature, the spectra are based on computational methods rather than measurements. 1,[9][10][11][12][13][14] There are a few reports of neutron doseequivalent measurements made in proton therapy centers with either the wide-energy neutron detection instrument (WENDI) or smart WENDI (SWENDI), a type of Rem-meter that is sensitive over an energy range adequate for measuring neutrons from proton therapy. [14][15][16][17] Like standard Remmeters, however, these detectors provide only a single doseequivalent value and do not provide the neutron spectrum. There are also several reports of measurements with detectors such as track etch detectors 14,18 or thermoluminescent detectors 19 that have limited or no response to neutrons with energies greater than 10 MeV. There are a few reports that include neutron spectra measured with a Bonner sphere spectrometer (BSS), 11,20 which can be used to deduce neutron spectra. However, these studies used standard BSSs, which have a very low and nondifferentiating response to highenergy neutrons (>20 MeV). 21,22 While the BSS is ideal for measurements of secondary neutrons from high-energy photon therapy, which have energies from thermal to approximately 10 MeV, 23 it is much poorer for measurements of secondary neutrons from proton therapy because of the abundance of high-energy neutrons. Tissue-equivalent proportional counters (TEPC) have also been used to measure microdosimetric neutron spectra from proton therapy. [24][25][26] While TEPC are sensitive over the entire neutron energy range of interest in proton therapy, the published microdosimetric spectra reported frequency distributions as a function of lineal energy rather than neutron fluence as a function of energy. In short, measurements of secondary neutrons from proton therapy either have not included spectral data or spectral results did not provide neutron fluence as a function of energy or were measured using detectors with limited response to neutron energies greater than 20 MeV.
We developed an extended range Bonner sphere (ERBS) spectrometer with high-atomic-number (Z) shells placed within polyethylene spheres to increase high-energy response compared to a standard BSS. While various different ERBS systems [27][28][29][30][31] have been used to measure neutron spectra that extend up to very high energies, e.g., measurements of atmospheric radiation [32][33][34][35] and around nonmedical accelerators, 28,30,[36][37][38][39][40] this detection technique has not previously been used for measurement of neutron spectra associated with clinical proton therapy. The objectives of the present study were therefore to measure a secondary neutron spectrum using an ERBS, to report neutron fluence as a function of energy, and to calculate the ambient dose equivalent from proton therapy.

2.A. Description of Bonner sphere extension
Measurements were performed using the Bonner sphere extension (BSE) measurement system. The detector design is briefly described here; a more thorough description is provided by Burgett 28 and Howell et al. 30 The BSE is a modular multisphere system that expands upon the design of the exist-ing, commercially available BSSs by adding concentric shells of high-Z materials. The full system includes 6 standard Bonner spheres and 12 extended sphere combinations. The extended spheres comprise various combinations of two basic forms: the small assembly and the large assembly. Each uses a core Bonner sphere, the 7.6-cm [3-in.] sphere (small assembly) and the 12.7-cm [5-in.] sphere (large assembly), respectively. The core Bonner spheres are surrounded by a copper, tungsten, or lead shell. Both assemblies can be further encased in an outer polyethylene sphere (Fig. 1). The system can accommodate either the 6 LiI(Eu) scintillator or activation foils placed on a polyethylene holder (with dimensions similar to the Ludlum scintillator). In the present study we performed 18 separate consecutive measurements with all combinations of a 6 LiI(Eu) scintillator (Ludlum Measurements Inc., Sweetwater, TX) and each of the 6 standard Bonner spheres and 12 extended spheres.
This BSE system was previously tested in a high-energy neutron beam (thermal to 800 MeV) at Los Alamos Neutron Science Center. 28,30 Briefly, we measured a neutron spectrum with energies between thermal and 800 MeV with standard Bonner spheres and the BSE and compared those spectra with reference spectra, i.e., time of flight (TOF) and Monte Carlo (MC) simulated spectra. The agreement between both TOF and MC reference spectra was ∼6 times better for the spectrum measured with the BSE than for that measured with standard Bonner spheres. The BSE system was calibrated using monoenergetic neutron ISO 8529 reference beams at the Physikalisch-Technische Bundesanstalt (PTB), Braunschweig, Germany. 41

2.B. Facility-specific technical information
To enable a better comparison of our neutron data with that reported by other proton centers, we briefly review the technical features of our treatment system, i.e., The University of Texas MD Anderson Cancer Center Proton Therapy Center -Houston (PTCH). The proton treatment system (Probeat; Hitachi America Limited, Inc., Tarrytown, NJ) consists of a synchrotron (70-250 MeV) and a beam transport system that feeds into six treatment stations. The measurements described in this study were taken in a room with a rotating gantry with a passive-scattering nozzle; we have two such gantries at our institution. The nozzle includes a rotating range modulator wheel, a binary stack of range shifters, a contoured and compensated double scattering system, beam monitoring instrumentation, collimators, and a range compensator. Three distinct snouts are available: small, medium, and large, which provide different ranges of field sizes. This treatment equipment can produce beam energies (defined at the nozzle entrance) of 100, 120, 140, 160, 180, 200, 225, and 250 MeV and spread-out Bragg peaks ranging in size from 2 to 16 g/cm 2 . 42

2.C. Measurements
Measurements were performed using a 250-MeV passively scattered proton beam with a medium snout ( aperture opening) and closed brass aperture. All measurements were performed using the 6 standard Bonner spheres and the 12 extended sphere combinations in air with the 6 LiI(Eu) scintillator at a distance of 100 cm laterally (at 90 • ) from the isocenter [ Fig. 1(c)]. Measurements were repeated separately for each sphere. A total of 1 Gy was delivered for each measurement. The proton beam was calibrated by a board-certified medical physicist working at the PTCH facility, according to the method specified in International Atomic Energy Agency Technical Report 398, 43 to deliver 1 cGy per MU at isocenter under reference conditions. Sphere-specific calibration factors were applied to all measured data. 41

2.D. Spectrum unfolding, calculated neutron parameters, and uncertainty analysis
Data were unfolded using the MAXED MXD_FC33 unfolding algorithm part of the UMG software package, version 3.3. 44 The MAXED (MAXimum Entropy Deconvoloution) algorithm is fully described in the literature [44][45][46] and is briefly summarized here. MAXED is an iterative algorithm that utilizes the maximum entropy principle to solve the closed form solution to the Fredholm integral equation. The endpoint of the measurements with any Bonner sphere system is a count rate for each detector-sphere combination. These count rates are related to the neutron flux at the measurement location by the Fredholm integral equation where C i is the count rate for the ith detector-sphere combination (counts/s); R i (E) is the response function for the ith detector-sphere combination at energy E (counts cm 2 ) per unit flux, shown in Fig. 2 for the BSE; and φ(E) is the energydependent neutron flux (n/cm 2 s). The MAXED algorithm solves Eq. (1) to obtain φ(E) for a given set of measured C i data. Practically, this is accomplished by reducing Eq. (1) to a closed form matrix, Eq. (2): where φ g is the group flux between E g and E g+1 ; R ig is the multigroup form of the ith detector response; and G is the number of energy bins. This system of equations is underdetermined as there is a finite number of equations (one for each detector-sphere combination) and substantially more G group flux energy bins to determine. The group fluxes are determined by iterating on the predicted count rates (for each detector-sphere combination), which are obtained by folding the response functions with the group fluxes. Because the solution to Eq. (2) is underdetermined, the MAXED algorithm requires an a priori spectrum, e.g., an initial estimate of the spectrum for unfolding. This spectrum should be physically realistic and based on some knowledge of the spectrum being measured. The MAXED algorithm obtains the solution to the unfolding problem by maximization of the relative entropy subject to constraints imposed by the measurements. This approach permits the inclusion of prior information (an a priori spectrum) in a well-defined and mathematically consistent way, and it leads to a solution spectrum that is a non-negative function and can solve for a spectrum in closed matrix form. This last feature permits the use of standard methods for sensitivity analysis and propagation of uncertainties, which are implemented in the IQU_FC33 program (also part of the UMG software package). 46 The a priori spectrum and response functions used in this work are described in the following paragraphs. The a priori (i.e., starting) spectrum for neutron unfolding was determined using the Bayesian statistical estimator methodology. A statistical fit of the measured data was determined using a Bayesian estimation program written for the UMGjava package. 47 The realistic starting spectrum was estimated using physical features that are known to be present in a proton-induced spallation neutron field. The Bayesian estimator was constrained with a high-energy Gaussian shaped peak for the spallation neutrons, a Maxwellian evaporation peak, a thermal room return peak, and an intermediate energy continuum (with 1/υ dependence) between the thermal and evaporation peaks. The maximum energies of the high-energy Gaussian and evaporation peaks were set to 250 and 15 MeV, respectively. Each of the energy peaks was allowed freedom in amplitude, centroid energy location, full width at half maximum values, and skew. The intermediate energy region was allowed to be linearly continuous, varying its high-energy and low-energy end points in relation to the centroids of the thermal and evaporation peaks, its beginning and ending amplitude, and its slope. After reading in each of the sphere's independent responses and measured responses, it produced a "best fit" to the measured data using the physical parameters. The Bayesian estimation program then created a list of best fit parameters. These parameters were used to produce a single starting spectrum for use in the MAXED unfolding code.
We used response functions specifically calculated for the BSE. 28,30 Briefly, response functions were calculated using Monte Carlo N-Particle Code, eXtended (MCNPX) code, version 2.7a, 48 using 500 equally spaced energy bins in the loga-rithm of the neutron energy for each 6 LiI(Eu) detector-sphere combination. This high resolution response is needed to provide additional resolution in the 100 keV to 100 MeV energy range and integrates with other detector systems in use by the PIs which are used for in situ neutron spectral measurements which will be described in future work. Responses were modeled for the 6 standard Bonner spheres and the 12 extended spheres. Neutron transport was carried out using three sources of data: ENDF-B.VI.8 (Ref. 49) cross-section data were used for transport less than 20 MeV; LA150 cross-section data for transport between 20 and 150 MeV; and CEM03 physics models 50 for neutron transport greater than 150 MeV. The response function for each of the detector-sphere combinations was normalized to a unit neutron source, i.e., n,α reactions for the 6 LiI(Eu) scintillator per unit neutron fluence rate. The response functions for the extended spheres showed greater sensitivity to neutrons with energies greater than 10 MeV (Fig. 2) than the standard spheres.
The measured neutron spectra were used to calculate three parameters of interest: the average energy, total neutron fluence, and ambient dose equivalent [H* (10)]. The ambient dose equivalent was calculated for the solution spectrum using ICRP-74 conversion coefficients. 51,52 The uncertainty in the unfolding process was determined using the IQU_FC33 program. 44 This program uses the fluence output file to propagate the uncertainty in the measured data and the default spectrum. Reginatto et al. provided a full description of error propagation in the IQU_FC33 program. 46

RESULTS
The measured (unfolded) neutron fluence spectrum is shown in Fig. 3. The neutron fluence per proton Gy to isocenter per unit lethargy is shown on a linear scale versus the logarithm of neutron energy. This formatting of the data ensures that the area under the curve is conserved and the relative fluence and contribution from each energy bin is proportional to the area under the histogram. The total fluence for this measured spectrum at 100 cm laterally from the isocenter was 6.6 × 10 6 n/cm 2 per proton Gy (to isocenter); the mean energy was 12.3 MeV. As expected, the fluence spectrum (Fig. 3) has a high-energy direct neutron peak, an evaporation peak, a thermal peak (due to room return), and an intermediate energy continuum between the thermal and evaporation peaks.
The ambient dose-equivalent spectrum is shown in Fig. 4. The dose-equivalent data were plotted in a manner that was analogous to that used to plot fluence data (Fig. 3). It should be noted that there is no low-energy tail in the dose-equivalent spectrum. This is because the ambient dose-equivalent conversion factors for energies less than 1 keV are very low. The dose equivalent was dominated by the neutrons in the evaporation peak because of both the high fluence and the large fluence-to-ambient dose-equivalent conversion coefficients in that energy interval. Approximately 35% of the dose FIG. 4. Neutron ambient dose equivalent H*(10) spectrum for 250-MeV proton beam, medium snout, closed aperture, in air at a distance of 100 cm from isocenter. Dose-equivalent spectrum was calculated using measured spectrum (Fig. 3) and neutron fluence-to-ambient dose-equivalent conversion coefficients from ICRP Publication 74 (Refs. 51 and 52). equivalent was from neutrons with energies ≥ 20 MeV. The total dose equivalent 100 cm laterally from isocenter was 1.6 mSv per proton Gy (to isocenter).

DISCUSSION
We present here a secondary neutron spectrum (with fluence as a function of energy) from a 250-MeV passively scattered proton beam measured by using an extended-range Bonner sphere measurement system that is sensitive to neutron energies over the entire energy spectrum of interest, i.e., from thermal to the maximum energy of 250 MeV. We also report the neutron ambient dose equivalent and mean energy for this spectrum. The measurements were carried out using the highest energy available at our facility and for a closed aperture, thus representing the maximum external neutron production. The neutron spectrum reported here represents external neutrons created in the treatment head but does not include internal neutrons. It has been reported in the literature that, for passive scattering proton therapy, external neutrons are the dominant source of neutrons when treating small fields, as in prostate cancer. 53 In the spectrum reported here (Fig. 3), the fluence in the direct neutron peak (∼70 MeV) is lower than that in the evaporation peak (∼0.5 MeV). We attribute this to the measurement location being 100 cm from the isocenter. At isocenter, the direct neutron peak would be expected to be larger than the evaporation peak. As the distance from isocenter increases, the relative fluence under the direct neutron peak decreases and that under the evaporation peak increases. This general trend in spectra shape as a function of distance from isocenter was observed in Monte Carlo simulations of neutron spectra at increasing distances from the isocenter. 1,54 It is difficult to compare our data with that reported in the literature because of the inherent differences in beamlines at different proton facilities and in measurement conditions, distance from isocenter, aperture configuration, and similar parameters. Our dose equivalent result (1.6 mSv/Gy) is consistent with measured neutron dose equivalents reported in the literature, which vary from 0.25 to 10 mSv/Gy, 11, 15-17, 19, 20 the higher end of the scale occurring with high proton beam energies measured near the isocenter and more collimating material in the beamline.
Of the measured neutron dose equivalents reported in the literature, those reported by Yonai et al. had the experimental set-up most comparable to that used in our study (Table I). Yonai et al. reported an ambient dose equivalent of 1.15 mSv/Gy at 100 cm lateral to the isocenter for a 230 MeV proton beam, which is in good agreement (30% lower) with ours, 1.6 mSv/Gy. We attribute the difference in the reported dose equivalents to differences in the beamline, detector, and experimental set-up. For example, our final beam was completely closed, while Yonai et al. used a more moderate cone down (5 × 5 cm 2 ); we would expect the additional high-Z material (brass) in the primary proton field to increase the neutron production. Other differences in experimental set-up that would cause our neutron dose to be higher than that reported by Yonai et al. 17 are higher beam energy (250 versus 230 MeV) and shorter distance between the final brass aperture and the isocenter (20 versus 30 cm). Also, the dose equivalent reported by Yonai et al. 17 may be somewhat underestimated because the WENDI-II detector underresponds to neutrons with energies between 10 and 100 MeV. 55 One difference between the two experimental set-ups that would result in higher neutron production by the Yonai et al. set-up was that they performed their measurements with a water phantom at the beam isocenter, whereas we did not. They estimated that the contribution from internal neutrons generated in the phantom was relatively low (∼20%) because the incident field size was small (5 × 5 cm 2 ) and thus external neutrons were the dominant source of neutrons. 12,53 Given the various differences in the beamlines, detector types, and experimental set-ups between the two studies (most of which caused our dose to be higher), an agreement of 30% is very good.
We can also compare our data with those produced by Monte Carlo simulations. Data reported by Zhang et al. 1,12 are of particular relevance here, because these authors reported a neutron spectrum that was simulated for the MD Anderson Cancer Center PTCH under the same conditions with which we conducted our measurements: 250 MeV protons with medium snout and closed aperture in air at a distance of 100 cm laterally from the isocenter. Qualitatively, our measured spectrum is similar to the shape of their simulated spectrum, 1 with a few differences. Their evaporation peak was somewhat broader than ours and our low-energy tail is more  56 A strength of the present study is that we used a detector that is sensitive to the entire energy range of the neutron spectrum being measured. However, one limitation of this detector system is that measurements are time consuming because they must be repeated separately with each moderating sphere. Thus, we were only able to carry out measurements for a single experimental set-up, i.e., single nozzle configuration and measurement location. Another limitation of the present study is that we assessed only external neutrons and not internal ones. Nevertheless, this work provides detailed measurements that serve as a reference for the upper bound of externally generated neutrons at our facility and can also be directly compared with Monte Carlo data in the literature and with similar data from other institutions (when such data become available).
Another important strength of our work is that we reported a fluence spectrum compared to a single dose-equivalent value as would be recorded from a WENDI detector, which is also sensitive to neutrons over the entire energy spectrum of interest. The neutron fluence spectrum enabled us to calculate dose equivalent as a function of energy, thereby experimentally showing the contribution of neutrons of different energies to the total dose equivalent. In this study, we found that 35% of the dose equivalent resulted from neutrons of energies ≥20 MeV. This is the component of the dose equivalent that would not be measured or would be grossly underestimated if a neutron detector with limited response at energies ≥20 MeV was used for these measurements; we would expect the magnitude of this effect to be even greater closer to the isocenter where the direct neutron component is higher.

CONCLUSION
We measured a neutron spectrum for external neutrons generated by a 250-MeV proton beam using a BSE measurement system that was sensitive to neutrons over the entire energy range being measured, i.e., thermal to 250 MeV. We used the neutron fluence spectrum to experimentally demonstrate the contribution of neutrons with different energies to the total dose equivalent and in particular the contribution of high-energy neutrons (>20 MeV). These are valuable reference data that can be directly compared with Monte Carlo and experimental data in the literature.

ACKNOWLEDGMENTS
The authors acknowledge Ms. Kathryn Hale from The University of Texas MD Anderson Cancer Center Department of Scientific Publications for assistance in editing this paper. This research was supported in part by a subcontract of National Cancer Institute Award No. 1R01CA131463-01A1, Northern Illinois University (subcontract of Department of Defense Award No. W81XWH-08-1-0205), and by the Na-tional Institutes of Health through the MD Anderson Cancer Center Support Grant No. CA016672. The authors also thank Dr. Angelica Perez-Andujar and Dr. Stephen Kry for their feedback on various drafts of this paper.