Measurement of secondary neutron spectra induced by 480 MeV proton and 430 MeV/u 4He beams with a thick aluminum target

Knowledge of the characteristics of secondary neutrons produced by the interaction of Galactic Cosmic Radiation with spacecraft shielding materials is becoming increasingly important for predicting and mitigating biological risks of space explorers during deep-space travel. Hadron accelerators for medical applications are well suited to reproduce part of the conditions found in deep-space in terms of ion species and energies. The objectives of this work are to measure the secondary neutron spectra produced by proton and helium ion beams hitting an aluminum target with energies that correspond to the Galactic Cosmic Radiation peak during solar minimal activity and to validate and compare physical models of Monte Carlo simulations. Neutron spectra were measured with the extended-range Bonner sphere system NEMUS at two positions, 0° and 90° relative to the direction of the primary ion beam. The experimental setup consisted of 480 MeV proton and 430 MeV/u 4He beams colliding with a 30×30×63.5 cm3 aluminum target. The experimental neutron spectra were analyzed using the MAXED unfolding code and compared to several Monte Carlo simulation codes. The results show deviations in terms of the shape of the neutron energy distributions ranging between 1% and 14% and of the integral quantities of fluence and ambient dose equivalent ranging between 1% and 5.2%.


Introduction
Galactic Cosmic Radiation (GCR) poses a major health risk to spacecraft crews [1].It is composed of an isotropic low-intensity flux of baryons, from protons up to ions with high atomic number and energy (HZE), with energies from hundreds of MeV/u to TeV/u, and intensity modulated by the solar activity [2,3].Currently, the only technologically feasible strategy for protecting astronauts against GCR is passive shielding.It has been estimated that behind 20 g/cm 2 aluminum shielding, when the solar activity is minimal, protons, alpha particles and heavy ions (3 ≤  ≤ 28) contribute 58.1%, 12.2% and 29.7% of the relative effective dose, respectively [4].An important contribution of protons and alpha particles to the relative effective dose is caused by nuclear collisions of these primary ions through the interaction with the shield that produces a secondary radiation field composed mainly of high energy light charged fragments and neutrons.Many studies have shown that this secondary radiation field intensifies as the thickness of the shielding material increases [5,6].Neutron contributions are particularly important because of their high penetrative power and high Relative Biological Effectiveness (RBE).In fact, they may become the main cause of absorbed dose and equivalent environmental dose behind shielding [7].However, studies of the secondary radiation fields performed with stochastic and deterministic transport codes are mainly based on physical models and not on empirical data given by differential cross sections [8].In recent years, many studies of secondary neutrons produced by the interaction of high-energy ions with space-relevant materials are being conducted to fill this data gap [9][10][11].This scarcity of experimental data is mainly due to the low number of facilities capable of accelerating ions to energies relevant for the study of GCR effects.Medical hadron-therapy accelerators may be good candidates, although many of them can only accelerate protons to up to 250 MeV and/or 12 C ions to up to ∼400 MeV/u, these energies are present in the GCR spectrum [12].A second problem concerns the detectors (or tools) needed to characterize neutron spectra at energies above 20 MeV.Currently, the most used detectors are organic scintillators [13] which have good energy resolution but are limited in the measurable energy range and Bonner Sphere sets [14] which have low energy resolution but a large measurable energy range.
-1 -For this work, a Bonner sphere set was used to characterize neutron spectra at two measurement positions produced by two beams of 480 MeV protons and 430 MeV/u 4 He ions bombarding a thick aluminum target at the Heidelberg Ion Therapy (HIT) [15] experimental room.This configuration is of particular interest for two reasons, the first is that the contributions of protons and 4 He ions become relevant when thick shielding is used [8], and the second is that the chosen energies are on the peaks of the GCR spectra when solar activity is minimal [1].The secondary radiation field was simulated with four Monte Carlo codes, FLUKA, MCNP6, PHITS and Geant4, using the available physical models, and a quantitative analysis was done by comparing the simulated spectra with the experimental spectra obtained through an unfolding of the measured data.

NEMUS Bonner sphere set
The PTB uses for neutron field characterization a Bonner Sphere set called NEMUS [16] (NEutron MUltisphere Spectrometer), which consists of a set of 10 full polyethylene spheres with a radius between 3 and 12 inches and four "modified" spheres which contain a lead or copper shell contained between two polyethylene shells (designed specifically for high energy neutrons).When a neutron enters the sphere, it interacts through elastic scattering events with the hydrogen in the polyethylene, gradually reducing its energy.Neutrons that have not escaped or are captured can become thermal and reach the central thermal neutron sensor (CTNS) located in the center of the sphere, allowing neutron detection.The response or "reading"   of a Bonner sphere of diameter , can be calculated using a first-order Fredholm integral equation: where  is the energy spectrum of the neutrons incident on the sphere.  is the neutron response matrix of the sphere of diameter : it is the mathematical operator that contains all the responses of the sphere as a function of neutron energy.The spherical shape of the sphere guarantees the isotropy of the response.A key factor that characterizes sphere reading is the type of CTNS used.In this study, a spherical proportional counter (SP-9 type, produced by Centronic) filled with 3 He gas at 0.2 bar (∼20 kPa) was used.For more details on the data acquisition and CTNS spectrum processing system see reference [17].

HIT experimental setup
The results presented in this work were obtained at the Heidelberg Ion Therapy (HIT) Center [15].The facility has four ionic species, protons, 4 He, 12 C and 16 O, which can be accelerated between 48 MeV/u and 430 MeV/u (protons up to 480 MeV).The ion beam is delivered to a room used exclusively for physical measurements.The secondary radiation field measured in this study was obtained using two types of ion beam.The first consists of 480 MeV protons with an intensity of ∼ 4 • 10 8 ions/spill and a beam FWHM of 5 mm.The second is a 430 MeV/u 4 He ion beam with an intensity of ∼ 1•10 8 ions/spill and a beam FWHM of 10 mm.For both beams a spill duration of 10 s was chosen.A large area Ionization Chamber (IC), placed in front of the beam exit window, was used for primary particles monitoring [18].The two beams bombarded plates of the aluminum alloy AlMg3 (density = 2.67 g/cm 3 ) attached to -2 -each other in series for a total size of 30×30×63.5cm 3 .This target thickness was chosen to completely stop the primary ion beam within the target.Bonner spheres were placed at 0 • and 90 • to the beam direction and at 2.1 m to the beam exit window.Figure 1 shows a sketch of the experimental setup implemented in the HIT experimental room with labeled NEMUS measurement angles.In this work, the results obtained for the 90 • measurement position for both ion beams were analyzed, but they are not included in the current article for reasons of space.More emphasis is given instead to the results for the 0 • measurement position to focus on the results most relevant to the topic at hand.

Monte Carlo simulations
To get an estimate of the secondary neutron field in the HIT experimental room, four Monte Carlo codes, FLUKA, MCNP6, PHITS and Geant4, were selected.In all Monte Carlo simulations, a sufficiently realistic geometric model of the experimental room was used.The geometric model of the experimental room of HIT was built using the FLAIR [19] graphical user interface, figure 1 shows the geometric model built with FLAIR.The model is simplified but takes into account the actual size and material compositions (concrete walls, aluminum table, etc.) of the experimental room.Due to the scarcity and/or complete lack of neutron cross sections for energies above 20 MeV and ions, it is only possible to simulate particle transport and interactions through physical models.For each of the Monte Carlo codes selected for this study, the various available physical models were chosen to perform a quantitative analysis.
The first code is FLUKA [20,21] version 4.2.2.It uses the Dual Parton Model (PEANUT framework) for hadronic collisions (protons and neutrons) for energies above 20 MeV.For heavy ion interactions having energy above 0.1 GeV/u it uses the RQMD (Relativistic Quantum Molecular Dynamics) model, while below 0.1 GeV/u it uses the BME (Boltzmann Master Equation) model.For neutrons having energy below 20 MeV, FLUKA uses neutron cross sections from various Evaluated nuclear data libraries (see the manual for more information).The second code is MCNP6 [22] version 6.1.The CEM03.03 (Cascade-Exciton Model) was chosen for interactions at energies below 150 MeV, -3 -while for energies above this threshold the LAQGSM (Los Alamos Quark-Gluon String) model was activated.The ENDF/B-VII neutron cross-section library was used for neutrons having energies below 20 MeV.MCNP6 has a neutron induced cross-section library LA150N for some materials for energies up to 150 MeV.The geometric model used for MCNP6 was obtained through FLAIR which can export the geometric model implemented in FLUKA into one readable by MCNP6.The third code is PHITS [23,24] version 3.32.The physical models chosen for ion transport are INCL (Liege Intra-Nuclear Cascade) and JQMD (JAERI Quantum Molecular Dynamics), where for protons only the former was used while for 4 He particles both were used.PHITS uses the JENDL-4.0 library of neutron cross sections for energies below 20 MeV, while for energies below 200 MeV, it uses the JENDL-4.0/HElibrary.The geometric model used for MCNP6 can also be used in PHITS.The last code considered for this work is Geant4 [25] version 11.0.The physical models chosen for ion transport and collisions at high energies are BIC (Binary Intra-Nuclear Cascade), INCL and QMD (Quantum Molecular Dynamics).For protons only the first two models can be used, while for 4 He particles all three have been used.The physical list QGSP_BIC_HP was used for neutrons, which uses the JEFF-3.3neutron cross section library.The geometric model used in Geant4 was obtained by converting the geometric model used in FLUKA using the PyG4ometry [26] python library.
All neutron spectra simulated with the four codes for the measurement positions described in the previous section were extracted by an averaged cell fluence tally and they were normalized for simulated primary particle.

Simulated neutron field
The secondary radiation field inside the HIT experimental room induced by the 480 MeV proton beam and the 430 MeV/u 4 He particle beam bombarding the aluminum target was simulated with the Monte Carlo codes and the physical models mentioned in the previous paragraph.In addition to secondary neutrons, secondary protons are also of particular interest since they can interact with the metal shells of modified Bonner Spheres producing additional neutrons by nuclear spallation.The latter can thermalize and reach the CTNS creating parasitic counts [9,27] and consequently they must be subtracted to avoid overestimations.Simulations showed that secondary protons induced by the primary beam of protons and 4 He particles can be neglected since they contribute to the modified spheres responses at 0.24% and 1% for the 0 • measurement position, respectively, while for the 90 • position both contribute less than 0.01%.Figure 2 shows the comparison between simulated neutron spectra in lethargy representation at the 0 • measurement point obtained with 480 MeV proton beam (left) and 430 MeV/u 4 He beam (right).For the proton beam, in the energy range greater than 20 MeV, the direct peak has an average peak energy of (252 ± 128) MeV produced by nuclear spallation of the proton beam with the target, while for the 4 He beam, the direct peak has a mean energy of (350 ± 73) MeV produced by nuclear fragmentation.In both spectra, in the energy range between 100 keV and 20 MeV, called evaporation range, there are peaks of evaporation neutrons which are produced by the excited nuclei from the high-energy protons and neutrons of the materials in the experimental room.In the epithermal region, between 0.4 eV and 100 keV, neutrons interact through elastic scattering with hydrogen nuclei in the concrete walls until they thermalize, creating the thermal peak in the energy range below 0.4 eV.For this study, the component of interest is the direct peak, since the -4 -components constituting the neutron spectrum for energies lower than 20 MeV depend mostly on the room geometry implemented for the simulation, which is very simplified.
Table 1 lists neutron fluence values from Monte Carlo simulations using the proton and the 4 He beam for the 0 • measurement position; these values are integrated in the neutron energy ranges of interest and normalized for simulated primary particle.For the proton beam, the results showed that the direct peak contributes an average of 23% of the total neutron fluence, while the evaporation energy range contributes 41%.These results showed that FLUKA produces the highest fluence for the direct peak, while PHITS and Geant4, both using the INCL model produce the lowest fluences.For the 4 He beam, the results showed that the direct peak contributes an average of 55% of the total neutron fluence, while the evaporation energy range contributes 26%.MCNP-CEM, PHITS-ICNL and Geant4-QMD show the highest values of neutron fluence on the direct peak, while the lowest value is shown by Geant4-BIC.-5 -

Neutron spectra unfolding
To carry out a quantitative comparison between the simulated results obtained with the various physical models shown in the previous paragraph, it is necessary to compare them with the experimental results.
Figure 3 shows the comparison between experimental readings measured with the Bonner Sphere set for the measurement position at 0 • induced by the 480 MeV proton beam and those calculated with Monte Carlo codes by applying eq.(2.1) normalized per primary particle.The readings of the spheres numbered 22 to 24, indicate the results obtained with the modified spheres, see reference [16] for more information about the spheres.From this comparison it is difficult to assert which simulation best reproduces the neutron spectrum.It is necessary to reconstruct the measured neutron spectrum starting from the experimental data.This goal can be achieved through an unfolding procedure.It consists of the discretization and solution of the inverse problem of estimating the incident neutron spectrum ϕ that generated the sphere reading   .In this work, the code MAXED [28], based on the maximum entropy method, was chosen.The code to be initialized needs as input information: the response matrix of the NEMUS set, the experimental readings, their associated statistical uncertainties, and an initial default neutron spectrum.With the MAXED code, a set of solutions was obtained using all simulated spectra as input, since the solutions are similar but not coincident.This solution set represents the confidence interval in which the true experimental neutron spectrum most likely lies, taking into account among all unfolded spectra, bin by bin, the minimum and the maximum value.Figure 4 shows the comparison between the simulated neutron spectra (lines) and the MAXED solution set (purple area) obtained with the proton beam (left) and the 4 He beam (right) for the 0 • measurement position.Concerning the spectrum produced by the proton beam, the results showed that the simulated spectra that best reproduce the direct peak are Geant4 and PHITS with the INCL physical model, with a deviation of 1% and 1.5% respectively.Concerning the spectrum produced by the 4 He beam, the simulated spectra that best reproduce the experimental direct peak are MCNP-CEM, PHITS-INCL and Geant4-QMD with deviations of the order of 11%, 12% and 14%, respectively.
-6 -Analyzing the integral values of neutron fluence for the energy ranges described in section 3.1, the direct neutron peak of the proton beam-induced spectrum at the 0 • measurement position has an average fluence of 3.05•10 −6 cm −2 with a standard deviation of 5%.The average total neutron fluence amounts to 1.46•10 −5 cm −2 with a standard deviation of less than 1%.This means that regardless of the Monte Carlo code used and the physical model chosen, the integral value of the total neutron fluence is stable.Regarding the spectrum produced by the 4 He particle beam, the integral neutron fluence in the energy range of the direct peak is 1.09•10 −4 cm −2 with a standard deviation of 5.2%, and that of the total spectrum is 1.78•10 −4 cm −2 with a standard deviation of 2.37%.As with protons, for the 4 He beam, the total integral value calculated from the various unfolded neutron spectra is stable regardless of the Monte Carlo code and the corresponding physical model chosen.
The stability of the unfolded neutron fluence integrals is also reflected in the calculation of the neutron ambient equivalent dose  * (10), obtained by multiplying the unfolded neutron spectra with the fluence-to-dose conversion coefficients [29,30].The  * (10) value in the direct peak energy range is 3.06•10 −6 μSv and 4.14•10 −8 μSv both showing a standard deviation of 4.1% and contributing 30% and 72% of the total  * (10), for the proton and 4 He beam-induced spectra, respectively.

Conclusions
The NEMUS set was used to measure the neutron field produced by a 480 MeV proton and a 430 MeV/u 4 He beam bombarding a thick aluminum target, at two measurement positions (0 • and 90 • ) in the HIT facility experimental room.The measured results were unfolded with the MAXED code to calculate the experimental neutron spectra and compare them quantitively with those simulated with Monte Carlo codes.The results of the measurement position at 0 • showed that the simulations that best reproduce the direct neutron peak are PHITS-INCL and Geant4-INCL for the proton beam and MCNP-CEM, PHITS-INCL and Geant4-QMD for the 4 He beam.Finally, the total integral values of neutron fluence and  * (10) values in the unfolded spectra obtained using various Monte Carlo simulated spectra as initial estimates for the unfolding are similar to each other and exhibit low standard deviations.

Figure 1 .
Figure 1.Sketch of HIT experimental setup with the marked measurement positions.

Figure 2 .
Figure 2. Comparison between simulated neutron spectra in lethargy representation at the 0 • measurement point obtained with 480 MeV proton beam (left) and 430 MeV/u 4 He beam (right).

Figure 3 .
Figure 3.Comparison between measured and simulated NEMUS readings as a function of the BS number for the measurement position at 0 • induced by the 480 MeV proton beam, both data sets are normalized per primary particle.

Figure 4 .
Figure 4. Comparison between simulated neutron spectra (lines) and MAXED solution set (purple area) of the experimental neutron spectrum induced by protons (left) and 4 He beam (right).