Source term calculation and validation for ${}^{18}$F-production with a cyclotron for medical applications at HZDR

In this document we present the calculation and experimental validation of a source term for ${}^{18}$F-production with a cyclotron for medical applications operating at 18 MeV proton energy and 30 $\mu$A proton current. The Monte Carlo codes MCNP6 and FLUKA were used for the calculation of the source term. In addition, the radiation field around the $^{18}$O-enriched water target was simulated with the two codes. To validate the radiation field obtained in the simulation, an experimental program has been started using activation samples which are placed close to the water target during an ${}^{18}$F-production run of the cyclotron. After the irradiation, the samples are analyzed and the resulting activation is compared to Monte Carlo calculations of the expected sample activation. We find good agreement between simulations and experimental results, with most calculation to experiment (C/E) ratios well between 0.6 and 1.4.


Motivation
Positron emission tomography (PET) has developed into a standard tool for imaging methods in medicine. The required radionuclides are often produced with the aid of cyclotrons. Depending on the emitter to be produced, different nuclides are bombarded with protons that trigger nuclear reactions. In addition to the desired nuclide, neutron and gamma radiation is also produced during these nuclear reactions. These are the main source of the radioactive dose rate on the outside of the protection buildings and determine the shielding design. In addition, neutron radiation leads to the activation of the construction and building materials, which could be important for the decommissioning of the facilities. Therefore, the correct determination of neutron and gamma source terms is the imperative basic condition for a correct shielding calculation and thus for a sufficient protection of the employees.
Several approaches can be taken to obtain the source term needed for the shielding calculations. In one approach, the source term spectrum is determined using nuclear model programs such as ALICE-91 [1], and subsequently, the corresponding transport calculation (shielding calculation) is done using the obtained spectrum. The determination of the absolute number of emitted neutrons is then carried out on the basis of generated activities. For many reactions these are available in tabular form for different proton energies and a standard current [2]. This approach for the 18 F-production is applied in [3], assuming that the neutron source term originates exclusively from the desired reaction. In [4] the energy and angular distributions of the neutron source term were taken from the double differential data of the nearby reaction 14 N(p,n) 14 O, but also basing the absolute number of emitted neutrons on the production rate of 18 F. Often the source term can be traced back to a confidential information from the manufacturer of the cyclotron with little information on how it was obtained [5,6,7]. For the shielding calculations for the new cyclotron with proton beam energies of 24 MeV and 28 MeV at the HZDR [8], a different approach was used, namely the direct calculation of the full neutron and gamma source terms including all contributing reaction channels with the radiation transport codes themselves. Different nuclear models are integrated into the program MCNP6 [9] to calculate the generation of neutrons in the target. Likewise, the source term can also be determined with the help of corresponding cross section tables. Both possibilities were used. The source terms calculated in this approach show a large difference respect to the values that were obtained using the approach mentioned above on the basis of tabulated activities. A calculation with the FLUKA [10,11] program gave similar results like MCNP6. These results have been already published and discussed in [8].
To validate the results from independent radiation transport codes, in addition to source term calculations with MCNP6 and FLUKA for HZDR's 18 MeV cyclotron ‡, this work gives experimental results for neutron fluence measurements using activation sample monitors. For these measurements, existing experience from the field of reactor dosimetry was applied.

Determination of the neutron source term
To calculate the neutron source term, simulation models were created with both MCNP6 and FLUKA which consist of a cylinder with radius 0.55 cm and a length of 4 cm, filled with water enriched with 97% 18 O. These dimensions correspond to typical target ‡ Cyclone 18/9 model by IBA bodies used at the IBA cyclotron. The protons' direction is along the cylinder axis, hitting the target on one of the circular base surfaces. The precise shape of the proton beam is not known. Therefore, two approaches were calculated. In the first case the proton beam was simulated as an infinitesimally small pointlike beam and in the second case a circular surface beam with a Gaussian distribution with a standard deviation of 0.125 cm cut off at the target radius was chosen. The two approaches gave identical results (see also [8]). In the following we will use the results obtained with a pointlike proton beam in the simulations. The emitted neutron spectrum is determined on the surface of a surrounding sphere with a radius of 10 m, large enough to minimize geometrical effects due to the target shape. The generation of neutrons in the target was carried out using nuclear physics models of reaction cross sections. In MCNP a cascade exciton model (CEM) [13] was used, while FLUKA uses a pre-equilibrium cascade model (PEANUT) [14] for the nuclear interactions. In addition, MCNP6 calculations were also carried out with evaluated nuclear data of the (p,n) reaction. Since 18 O-data are not included in the standard library of MCNP6, they were generated using the NJOY [15] program and imported into MCNP6. The required reaction crosssections were read from the nuclear data library TENDL, which is based on the nuclear core model code TALYS [16]. This possibility to use externally generated cross sections does not exist for FLUKA. Since FLUKA does not include neutron cross sections for 18 O, cross section data for 16 O was used instead for the interactions of neutrons. Further differences in the calculations consist of the used cross section data libraries. The data libraries ENDF/B-VI.8 [17] were used for the interactions of neutrons with energies below 20 MeV at FLUKA and ENDF/BVII.1 [18] for MCNP6. In fig. 1, both the proton and neutron fluences per primary proton are depicted as obtained with the FLUKA simulation code. The protons penetrate only about 0.5 cm into the water target before they are stopped. The water target absorbs almost all protons in the forward direction, leaving only the backscattered ones to the left of the target. Neutrons are produced along the trajectory of the proton beam in the water and are then emitted isotropically. Fig. 2 shows the differential neutron rate recorded across the surrounding sphere

Energy in MeV
FLUKA MCNP6 Figure 2: Differential neutron rate for an 18 MeV proton beam with 1µA hitting the water target. The spectra are available from [19].
for 1µA of proton beam current obtained with MCNP6 (version 6.1.1) and FLUKA (version 2011.2x). Integrating the spectrum over energy we find a total neutron yield of 3.21 × 10 10 n/s for 1 µA of proton current for the FLUKA calculation, and 2.99 × 10 10 n/s for 1 µA of proton current for MCNP6. The higher yield obtained with FLUKA respect to the MCNP6 calculations has already been observed for 24 and 28 MeV protons in [8], and is attributed to differences of the underlying nuclear physics models. The values are about a factor 3 higher than the value of 1.115×10 10 n/s for 1 µA of proton current obtained from [20] for the 18 O(p,n) 18 F channel. We attribute the difference to additional neutron-producing reaction channels opening at 18 MeV proton energy, as suggested in [21]. It should be noted however that measurements of the neutron yield rate reported in [22] and [23] give results which are close to the value of 1.115×10 10 n/s for 1 µA §.

Experimental validation of the radiation field around the target
To validate the calculation of the source terms in sec. 2, activation monitor foils where placed on top of the irradiation target during a routine run for 18 F production. After irradiation, the activation of the foils was measured and compared to predictions from the radiation transport and reaction codes MCNP6 (version 6.1.1) and FLUKA. For these activation studies, a special developer version of FLUKA was used which includes updated information on branching ratios to meta-stable states [24] that is not yet § It was confirmed by the authors of [23] that the values quoted need to be corrected by a factor of 10 and the resulting neutron production yield should therefore read (1.55×10 10 ±1.03×10 9 ) n/s for 1 µA of beam current. We thank M. Hagiwara for this information.
available in the official FLUKA version [25].

Experimental setup to measure the radiation field with sample activation
Figures 3a and 3b show the individual activation monitor samples as well as the sample packages and their position on the irradiation container at the cyclotron. The samples consist of different metal foils made of pure metals or alloys. Table 1 shows the monitor samples and the reactions under study with the generated nuclides, the reaction threshold and their half-life. The selected metals are standard monitor materials which are inserted for neutron flux and fluence measurements at fission reactors for power determinations as well as for the validation of the results in reactor dosimetry [12]. As can be seen from table 1, several of the materials have reactions starting at different threshold energies. This makes it possible to study different energy regions in the spectrum. The monitor packages were positioned directly on top of the irradiation target in order to achieve a high neutron flux and thus high reaction rates. The irradiation took place during a regular 18 F production run. The energy of the protons was 18 MeV with average beam current of 25 µA and the irradiation lasted for 50 minutes.  Two packages with the same stacks of activation monitors were irradiated simultaneously. This allowed independent measurements using two independent laboratories for the sample activation analysis. The activation measurements were carried out at the "Department of Environmental Monitoring" and at the "Laboratory for Environment and Radionuclide Analysis" of the "VKTA -Strahlenschutz, Analytik & Entsorgung Rossendorf e. V." . The activity was determined by gamma spectrometry using high-purity germanium detectors. In order to detect nuclides with a relatively short half-life, some of the activation monitors were already examined about one hour after the end of irradiation. For some nuclides with longer half-life, the measurements were repeated at longer cooling times and the activity value at the end of irradiation was extrapolated back using the known half-life values. However, it was found in previous studies that these extrapolations were not always reliable, especially when there is a delayed production of a nuclide from an excited state. Therefore both labs were asked to provide the measured values at the time of the measurement and also the values extrapolated back to the time at end of irradiation (EOI). The quoted uncertainties by both labs were in the range between 5% and 25%, depending on the channel.  Given the neutron flux rate, the activities A i (t meas ) for each produced nuclide at a time t meas after irradiation in an energy bin i can be determined using the relation

Calculation of sample activation
In eq. 1, is the density of nuclei in the sample (in nuclei/(barn·cm)), V is the sample volume in cm 3 , σ i is the corresponding reaction cross section in barn for energy bin i,Φ i is the corresponding neutron flux rate obtained from the simulation in neutrons/cm 2 /s, t irr is the irradiation time in seconds and λ is the decay constant of the reaction product (in 1/s). The total activity is then the sum of the A i over all energy bins i. Given   an irradiation time profile, FLUKA conveniently gives the resulting nuclide activities in Bq/cm 3 for selected geometry regions at desired times directly in a tabular output. The required cross section data is hard-coded into the FLUKA program and cannot be changed by the user. For MCNP6, eq. 1 needs to be applied externally to the simulated neutron flux rates. In this case, the required cross section data had to be generated with the NJOY program. This procedure has the advantage that the neutron flux can be folded over cross sections from different nuclear data libraries. This allows to estimate systematic uncertainties coming from differences between the available cross section data sets. If more than one reaction channel contributed to a measured final state isotope, the resulting activities were added to obtain the final result.

Discussion of the results
In tables 2 to 5 the results for the measured and simulated activities for the different monitors at the corresponding time of measurement are presented. Measurements obtained by the "Department of Environmental Monitoring" are reported as "Analysis A" and the ones by the "Laboratory for Environment and Radionuclide Analysis" are reported as "Analysis B". We have only kept results for reactions for which both laboratories reported a significant value and for which the statistical uncertainty of the simulations was 15% or better. This e.g. excludes the reaction 55 Mn(n,g) 56 Mn in table 1, for which only Analysis A gave a measured result. In total 11 measurements for different nuclides remain. Numbers in parenthesis correspond to uncertainties on the last digits. The MCNP6 results in the tables use cross section data from the JEFF3.1A libraries for the Indium, Zinc and Tin monitors, while for the multi-component monitor ENDF/B-VII.1 were used (except for the 58 Ni(n,p) 58 Co reaction, for which also JEFF3.1A libraries were used). Table 2 gives the results for the multi-component monitor. The two analyses A and B were done with a time difference of 24 hours in the two different labs. Except for the 58 Ni(n,np) 57 Co reaction, the FLUKA code reproduces the experimental results reasonably well, with calculation-over-experiment (C/E) values between 0.77 and 1.20.   . An interesting fact is that the activity due to 58 Co is larger for analysis B. This is due to the fact that there is a delayed production of 58 Co due to the presence of the meta-stable state 58m Co which decays with a half-life of 9 hours to 58 Co. This is taken care of in the simulations, but not in the extrapolation of the measurements back to the end of irradiation. For this reason it was decided to compare the simulations with the experimental results at the time of measurement. For the results of the indium monitor in table 3, we find that both codes predict a factor 3 to 4 less activity for the production of 114m In, but are remarkably close to the measurements for the channel 115 In(n,n') 115m In.
The results for the production of 117m Sn with the tin monitor are given in table 4. The results with the FLUKA code give C/E ratios of 0.80 and 0.68, with the MCNP6 results being 17% lower in both cases. Table 5 gives the results for the zinc monitor. For the channel 64 Zn(n,p) 64 Cu, MCNP6 gives results which have a C/E ratio of 0.74 for analysis A and 0.97 for analysis B, with FLUKA results being consistently 7% lower. The activity due to the reaction 64 Zn(n,g) 65 Zn is simulated by the FLUKA code with C/E ratios of 0.71 for analysis A and 0.80 for analysis B, while for the channel 68 Zn(n,g) 69m Zn FLUKA predicts values significantly lower than the experimental results, namely a C/E ratio of 0.55 for analysis A and 0.38 for analysis B. For both reactions, MCNP6 predictions are about 25% lower than the FLUKA predictions.
In summary, for most of the reactions, the FLUKA code gives results with C/E ratios between 0.68 and 1.20, with the MCNP6 calculations giving in general results which are 10 to 40% lower (with the exception of the 64 Zn(n,p) 64 Cu and the 58 Ni(n,p) 58 Co reaction, for which MCNP6 results are 7 to 8% higher than the FLUKA results).
Both codes consistently give lower results for the production of 114m In and the channel 68 Zn(n,g) 69m Zn, with C/E ratios between 0.34 and 0.55 for FLUKA results and 0.24 and 0.42 for MCNP6 results. An especially large deviation between simulation and measured values is found for the reaction 58 Ni(n,np) 57 Co channel, with C/E ratios 0.16 and 0.18 for the FLUKA results and 0.12 and 0.14 for results obtained with MCNP6. Again, for these three channels, the MCNP6 results are lower than the FLUKA results by 20 to 30%.

Uncertainties in the calculations
To estimate the uncertainties on the calculations, we need to address the different terms in eq. (1). The primary sources of uncertainties in the simulations are the neutron flux rate at the sample position and the reaction cross sections. The neutron flux rate depends on the calculation of the source term (and therefore the underlying model for proton-induced neutron production in the water target), the proton beam current and to some extend the modeling of the target geometry. As mentioned in section 3.2, the geometry of the system has been implemented with great care. The two geometric models differ only in minor details (see fig. 4), and contain both identical material compositions and densities. We therefore consider systematic effects from the geometry negligible when comparing the two simulations, and also they are thought to have a minor effect when comparing simulation results to the measurement, even given the fact that the influence of possible backscattering of neutrons from surrounding walls was not considered. The neutron flux rate scales linearly with the proton beam current, and a deviation of the current from the nominal value of 25 µA will reflect on the neutron flux rate. However, the beam parameters can be determined very well and their uncertainties are very small. The dominant effect comes from the uncertainty of the source termas can be seen from fig. 2 this can reach up to 50% below energies of 1 MeV, and is certainly the reason why the MCNP6 results for most channels are significantly lower than the ones from FLUKA. A hint at the size of the uncertainty (which is energy dependent) is given by the C M /C F -values in tables 2 to 5. As can be seen, it can reach up to 40% and more in some cases.
Fundamental for the determination of the reaction rates are the data libraries used in the evaluations. Different evaluated nuclear data libraries are available for the majority of nuclides at e.g. NEA [26]. For some nuclides only data calculated from theoretical models is available. Using MCNP6's capability to include different cross section data sets, we have studied the effects of different cross sections on the simulation calculations. The cross section data was generated with the use of the NJOY program. Among the libraries tested were ENDF/B-VII, JEFF3.1A, JENDL33, EAF2010 and ROSFOND2010. If branching fractions for isomeric states were given, these were considered. While for some of the considered reactions (like 115 In(n,n') 115m In, 117m Sn production and 197 Au(n,g) 198 Au), the different libraries gave consistent results, in some cases discrepancies on the order of 30 to 40% could be observed ( 114 In production and the 68 Zn(n,g) 69m Zn reaction). As an example, in fig. 7 the evaluated cross section for the 68 Zn(n,g) 69 Zn/ 69m Zn reaction is shown obtained from two different data libraries (JEFF3.1A and the ENDF/B-VII). It can be seen that while for thermal energies, the cross sections are in good agreement, at high energies and around the first large resonance, the ENDF/B-VII data gives higher values, resulting in an activation result which is about 30% higher than the one calculated with the JEFF3.1A data. For the remaining reactions, the differences in the activation results due to the use of different cross section data sets is on the order of 10%.
Another point is the influence of the energy group structure. While the FLUKA code uses a fixed 260 energy groups structure, one can have a larger number of groups for the cross section data sets generated with NJOPY for MCNP6. For the reaction 58 Ni(n,p) 58 Co activities were calculated based on ENDF/B-VII.1 data determined with a 260 and 640 energy group structure. One finds an effect of 10% towards  a better agreement with the experiments using the high-resolution structure. This shows a general problem in the calculation of group-wise cross-sections. These may be underestimated in some cases, especially for threshold reactions. Generally, the flux decreases in the high energy region but cross-sections increase very strongly. Therefore it is very important to resolve the upper energy range well.
The volumes of the monitors in the two sample stacks are not fully identical, the simulation results were calculated using the average volume of two monitors of the same material. While this effect is at maximum 1% for the Indium, Zinc and the multicomponent monitor, it reaches about 5.5% for the Tin monitors (and therefore the simulation results on the 117m Sn production). Of course, this is negligible compared to the potentially large uncertainties on the neutron flux rate and the cross section spectra.
After the irradiation, we were notified by the operator of the cyclotron that the irradiation had to be stopped for about 10 minutes due to a vacuum problem. Once the problem was fixed, irradiation resumed to complete the 50 minutes of irradiation time. Due to the fact that the lifetimes of the produced isotopes in the reactions in table 1 are quite long, we do not expect a large effect due to this interruption. An additional simulation with the FLUKA package using an irradiation time profile of 30 minutes of beam, followed by 10 minutes of no beam and finally additional 20 minutes of beam gave indeed no significant differences within the statistical uncertainties respect to the calculation with a full uninterrupted beam for 50 minutes.
Finally, the values for material densities and half-life times λ were taken from the literature or are included in the simulation codes and the corresponding uncertainties are considered negligible for the present study.

Conclusion and Outlook
Inspired by the calculations of a shielding assessment for a new cyclotron bunker, investigations on the neutron source term for the 18 F production at a IBA Cyclone 18/9 cyclotron were carried out using the Monte Carlo transport and reaction codes MCNP6 and FLUKA. It was found that below 1 MeV, the MCNP6 code gave a differential neutron rate which is smaller than the one by FLUKA by up to 50%. The total neutron production yield for both codes was about 3 times larger than the value obtained from [20] for the exclusive 18 F production channel. To validate the results of the Monte Carlo codes, a more realistic model of the target geometry for the Cyclone 18/9 cyclotron was created with the two Monte Carlo codes which was used to calculate the activation of small monitor sample foils made of different metals and alloys during a typical run of 18 F production. These results were then compared to the actual activation of the sample foils after a 18 F run which was obtained using gamma spectroscopy with HPGe detectors at two independent labs. In total, 11 reactions were investigated, with C/E ratios between 0.6 and 1.4 for most cases. As a general trend, results calculated using the MCNP6 codes were 20 to 40% lower than the ones obtained with FLUKA. This may be (partially) explained by the fact that the source term obtained with MCNP6 is lower than the one from FLUKA below 1 MeV. For three reactions, the Monte Carlo simulations were consistently giving much lower results than the measured data (C/E values as low as 0.12 were observed). This was the case for the 58 Ni(n,np) 57 Co reaction, the production of 114m In and the reaction 68 Zn(n,g) 69m Zn. For these three reactions, the uncertainties discussed can not accommodate the discrepancies, and it is most likely that the underlying cross section data for these reactions is responsible for the results, and eventually this document may help to improve the cross section data base in the Monte Carlo programs in the future.
Despite the uncertainties of the measurements, the obtained results show that the calculation of the neutron source terms with the help of the methods and models which are implemented in radiation transport and reaction codes like MCNP6 and FLUKA should work better for a proton beam of 18 MeV than a calculation based solely on the 18 F yield. This is consistent with observations in [21], which reports significantly higher neutron yields for proton energies above 12 MeV for evaluations using a full ALICE-91 calculation respect to evaluations with tabulated yield values for the 18 O(p,n) 18 F reaction only. However, our observations seem to contradict the experimental results in [22,23] which find reasonable agreement with the yield of 1.115×10 10 n/s for 1 µA of proton current at 18 MeV energy obtained from [20] for the 18 O(p,n) 18 F channel.
In order to consolidate our results, further experiments will be planned at the new cyclotron of the HZDR. In these experiments, both the target material and the proton energy will to be varied. The aim is to provide validated absolute neutron fluence spectra for shielding calculations at medical cyclotrons.