#  Low-symmetry polymorph of GaP upends bonding paradigms of metallic high-
pressure III-V compounds

Barbara Lavina<EMAIL_ADDRESS>University of Chicago, GeoSoilEnviro
Center for Advanced Radiation Sources, Chicago, IL, 60439, USA  X-ray Science
Division, Argonne National Laboratory, Lemont, IL, 60439, USA  Enrique Zanardi
Andrés Mujica<EMAIL_ADDRESS>Departamento de Física and Instituto
Universitario de Materiales y Nanotecnología, MALTA Consolider Team,
Universidad de La Laguna, San Cristóbal de La Laguna 38200, Tenerife, Spain
Hyunchae Cynn  High-Pressure Physics Group, Physics and Life Sciences,
Lawrence Livermore National Laboratory, Livermore, CA, 94550, USA  Yue Meng
Current contact<EMAIL_ADDRESS>X-ray Science Division, Argonne National
Laboratory, Lemont, IL, 60439, USA  Vitali Prakapenka  University of Chicago,
GeoSoilEnviro Center for Advanced Radiation Sources, Chicago, IL, 60439, USA
Jesse S. Smith  HPCAT, X-ray Science Division, Argonne National Laboratory,
Lemont, IL, 60439, USA

###### Abstract

The pressure-induced polymorphism of binary octect compounds has long been
considered a settled problem although the possible atomic disordering of some
phases remains a puzzling observation. Taking GaP as a case study, we
conclude, through x-ray microdiffraction and first-principles calculations,
that its high-pressure phase II (previously reported as being disordered)
adopts in fact an ordered base-centered monoclinic structure previously
unknown in this class of compounds. The formation of layered patterns with
variable degrees of interlayer dimerization, as observed in GaP, marks a
paradigm shift of our understanding of ordering in octect high-pressure phases
which calls for a more extensive re-examination. A rich polymorphism with fine
tuning of chemical and physical properties can be envisioned.

Under ambient conditions GaP crystallizes in the zincblende structure (zb)
common to several other III-V and II-VI binary compounds and is an indirect-
gap semiconductor. A pressure-induced structural transition around 22 GPa,
which is associated with a change towards a semimetallic state, was first
reported in the 1970s [1] but the structure of the high-pressure phase GaP-II
remained elusive for more than twenty years. Taking advantage of interim
experimental advances in high-pressure crystallography, Nelmes, McMahon, and
coworkers re-evaluated in the 1990s the structure of the high-pressure phases
of the III-V and II-VI families unveiling several underlying systematic
behaviors which have prevailed up to now [2]. For GaP, their angle-dispersive
powder diffraction experiments showed that the phase obtained upon ambient-
temperature compression could be described as a distorted NaCl arrangement
[3], named _Cmcm_ after its space group, ruling out the occurrence of the
previously proposed $\beta$-Sn-type structure [4]. Diffraction data are
consistent with intra-crystalline atomic disorder [3, 5] whereas x-ray
absorption spectroscopy data suggest that the phase is ordered in the short
range [5]. As per the established systematics, the _Cmcm_ phase appears
ubiquitously in the phase diagram of the III-V and II-VI families [2, 6].

Recently, a new polymorph of GaP was synthesized at $\sim$43 GPa and high
temperature [7]. It adopts the so-called super-_Cmcm_ structure (or _oS_ 24 in
Pearson notation) similar to that previously described in InSb [8], but in GaP
the phase forms single homonuclear P-P bonds in striking contrast with all
previous understanding of bonding in the high-pressure polymorphism of III-V
compounds. Finally, the sc16 modification of GaP was found to be stable in a
narrow range between 17-20 GPa [9]. Considering that _oS_ 24 was synthesized
at pressures significantly greater than the onset of metallization, and
exhibited instability upon decompression within its calculated stability field
[7], we suspected that yet another phase between sc16 and _oS_ 24 might occur.
If the dimerization observed in _oS_ 24 occurs only in a higher-pressure
modification of the first metallic phase, it would not be necessarily related
to the physical property change. Therefore we explored such pressure range to
establish the nature of the first metallic phase of GaP. In addition, we
performed extensive theoretical analysis to understand more deeply the
evolution of bonding across polymorphs.

High-pressure and temperature conditions were generated using a laser-heated
diamond anvil cell (DAC). Structures were probed _in situ_ with synchrotron
microdiffraction techniques. In order to sample a sufficiently large portion
of reciprocal space, we used a wide-opening DAC equipped with conical anvils,
providing 64∘ 4$\theta$ access. Between the anvils’ culets, 400 $\mu$m in
diameter, Re gaskets indented to a thickness of about 40 $\mu$m and with
laser-drilled holes of 230 $\mu$m diameter provided the sample chambers. A
chip or powdered sample (99.999% purity from Alfa Aesar) along with ruby and
platinum as pressure gauges [10, 11] were placed in the sample chamber before
filling it with pre-pressurized neon [12]. High temperatures were obtained
using online double-sided IR laser-heating systems. The temperature was
determined from the emitted radiation [13]. The structure was probed using a
microfocused hard x-ray beam at Sector 16 (16-ID-B), HPCAT, and at Sector 13
(13-ID-D), GSECARS, of the Advanced Photon Source, Argonne National
Laboratory. Diffracted beams were collected with either MAR165-CCD (16-ID-D)
or Pilatus 1M (13-ID-D) detector that were calibrated using the powder
diffraction patterns of the CeO2 and LaB6 standards [14]. SXD data, with a
resolution down to 0.6 Å, were collected with the rotation method and
processed with CrysAlis Pro [15] and Shelxl [16].

In order to assess the energetics and stability of the new polymorph relative
to the established phases in the phase diagram of GaP, we performed total
energy calculations for this material including the novel phase described in
this Letter and phases previously observed or considered for the III-V family
[6]. The results that we will show here were obtained within the _ab initio_
framework of the density functional theory (DFT), using a basis set of plane
waves, with the outermost $s$ and $p$ valence electrons of Ga and P, as well
as the Ga semicore $d$ electrons, explicitly dealt with in a projector
augmented wave scheme [17, *Kresse:1999aa] as implemented in the vasp code
[19, *Kresse:1996fk]; the PBE generalized gradient approximation (GGA) was
adopted for the exchange-correlation functional [21, *Perdew:1996_erratum],
and the cutoff in the kinetic energy of the plane waves basis was set at 370
eV [23, 24, *Perdew:2008aa_erratum, 26]. Dense k-point grids were used for the
Brillouin zone integrations, e.g. a 16$\times$16$\times$8 grid was used for
the novel _mS_ 16 phase, which is metallic with eight atoms within the
primitive unit cell (see later). All the structures were fully relaxed under
hydrostatic conditions through the calculation of the stress tensor and the
forces on the atoms, to residual stress anisotropies and forces less than 0.1
GPa and 0.005 eV/Å, respectively, which ensures a convergence in the relative
energies of the phases at the chosen level of approximation of the order of
the meV per formula unit (pfu). The phonons for the main phases under
discussion were studied using the small-displacements method, and the effect
of temperature was taken into account at the level of the quasi-harmonic
approximation [27, 28]. A fuller account will be given elsewhere [29].

Figure 1:  Example of diffraction data of coarse-grained _mS16_ -GaP.  Figure 2:  (Top) The _mS_ 16 structure viewed from the ${a_{0}}$ axis of the primitive unit cell used in the calculations. Ga1 and P1 sites (within the L1 layers, viewed here in profile) are respectively shown in dark blue and red, whereas a lighter shade is used for the Ga2 and P2 sites that make the L2 layers. (Middle) View along the $c$ axis showing two adjacent, dimerized L1 layers (left) and two adjacent, overlapped L2 layers (right). (Bottom) Same as (Top) but showing a 2$\times$2$\times$2 slab with a modicum of perspective so as to reveal the full connectivity among sites. In all these plots the profiles of the primitive-cell boundaries used in the calculations $({a_{0}},{b_{0}},{c_{0}})$ are shown in black whereas those of the conventional cell $({a},{b},{c})$ are shown in green. (See also further plots in the Supplemental Material.)  Table 1:  Experimental structural data for the _mS_ 16 phase, corresponding to $p$ = 30.7 GPa and $T$ = 300 K. The structural refinement was conducted against 180 measured structure factors, 103 of which unique. A total of 13 parameters were refined, including eight fractional coordinates, four isotropic displacement parameters ($U_{eq}$) and the scale factor. Structure factors statistics is satisfactory with $R_{int}$=0.044 and $R\sigma$=0.028. The single-crystal structural refinement converged with $R_{4\sigma}$=0.038, $R_{all}$=0.045, $wR_{2}$=0.11. For comparison, our calculated theoretical values of the structural parameters are given in square brackets.  _mS_ 16 |  | SG _C2/m_ , No. 12
---|---|---
|  $a$ (Å) |  $b$ (Å) |  $c$ (Å) |  |  $\beta$ (∘) |
| 4.931(3) | 4.7189(11) | 9.553(2) |  | 102.67(3) |
| [4.905] | [4.732] | [9.590] |  | [102.97] |
Sites |  | Wyckoff positions |  $U_{eq}$ (Å2)
Ga1 | 4 _i_ | 0.0443(4), 0, 0.36638(11) | 0.0119(9)
|  | [0.0408, 0, 0.3662] |
Ga2 | 4 _i_ | 0.2504(4), 0, 0.11899(10) | 0.0115(9)
|  | [0.2504, 0, 0.1190] |
P1 | 4 _i_ | 0.5549(8), 0, 0.3935(3) | 0.0142(8)
|  | [0.5504, 0, 0.3934] |
P2 | 4 _i_ | 0.7510(9), 0, 0.1347(3) | 0.0146(8)
|  | [0.7516, 0, 0.1342] |

We performed several synthesis experiments between 15 and 30 GPa. A new phase
was observed between about 20 and 30 GPa. Diffraction images showed that the
samples developed variable grain size after laser-heating (Fig. 1). Single-
crystal structural analysis was performed on the largest grains harvested from
diffraction maps of the sample chambers [30]. At about 31 GPa a monoclinic
lattice with $a$=4.931(3) Å, $b$=4.7189(11) Å, $c$=9.553(2) Å, and
$\beta$=102.67(2)∘ was readily identified. The lattice showed obvious
geometric relations with that of the recently characterized _oS_ 24 polymorph
[7], which measured $a$=4.621 Å, $b$=4.927 Å, $c$=13.335 Å at 43 GPa. Taking
into account the pressure difference, the lengths of the $a$ and $b$ axes of
the two lattices are close, and the projection of the $c$ parameter of the new
phase onto the perpendicular to its $ab$ plane is roughly 2/3 that of the _oS_
24 phase. We inferred that the two phases likely have a similar atomic
arrangement in the $ab$ plane, _viz._ a distorted NaCl-like pattern, and that
the monoclinic unit cell of the new phase contains four NaCl-like layers
(compared to six in _oS_ 24), resulting in a crystallographic unit-cell
content of 16 atoms. With such constraints, solving the crystal structure was
rather straightforward, nonetheless the correct space group, _C2/m_ , was
assigned after calculations. Single-crystal structural refinements were
performed against the structure factors of selected grains showing the best
statistical parameters for observed intensities, and a selected result is
shown in Table 1. The new polymorph will be hereafter referred to using the
Pearson notation as _mS_ 16\. To the best of our knowledge, this structural
type has not been previously identified among the high-pressure polymorphs of
semiconductors and is a new structure type. Plots of the _mS_ 16 structure are
shown in Figure 2.

Two distinct features appear prominently in both the _mS_ 16 and _oS_ 24
phases: (i) their underlying two-dimensional NaCl-like arrangements, and (ii)
the partial dimerization of the phosphorus atoms. The structures of both _mS_
16 and _oS_ 24 can be described as different stackings of layers corresponding
to distorted NaCl-(001) planes, of which two distinct types can be found in
both phases: L1 layers, hosting the dimerizing P1 phosphorus and the Ga1
atoms; and L2 layers, containing the non-dimerizing P2 phosphorus and the Ga2
atoms (see Fig. 2). The P1-P1 dimers lie across adjacent L1 layers in both
polymorphs. The stacking is different than in NaCl as like atoms nearly
overlap (anti-NaCl). The L1 layer buckling is the result of two nearest P1
atoms approaching to form the dimer, whereas the two nearby Ga1 atoms maintain
a greater distance. In _mS_ 16 there are also two adjacent L2 layers, and
their relative stacking, while distorted, is NaCl-like, with interlayer L2-L2
bonds directed along the $c$ direction, hence all closest neighbors are unlike
atoms both within and across the pair of L2 layers. The stacking of L1 and L2
adjacent layers is shifted with respect to the NaCl structure by roughly
1/4-unit cell along the $a$-axis for _mS_ 16 and along the $b$-axis for _oS_
24\. The resulting stacking sequence is L1-L1-L2-L2 in _mS_ 16 and
L1-L1-L2-L1-L1-L2 in _oS_ 24, so that half the P atoms are dimerized in _mS_
16 whereas 2/3 form monatomic bonds in _oS_ 24\. As an effect of the
interlayer dimerization, the L1 layers exhibit a slightly larger amount of
corrugation than the L2 layers, as can be seen in Fig. 2. In _mS_ 16 there is
no puckering of the atomic rows _within_ the NaCl-like layers (either L1 or
L2), only transverse (out of plane) buckling/corrugation. This is unlike _oS_
24 (and _Cmcm_ [2]) for which there is puckering within the L2 layers whereas,
conversely, these layers do not have any corrugation – they are perfectly
flat, see Supplemental Material. The intra-layer puckering of _oS_ 24 is
related to the reduced coordination of the P2 sites in this phase [31].

The likeness of _mS_ 16 and _oS_ 24 is further evidenced in their similar
densities (Fig. 3, top). The $b$-axis of _mS_ 16 and the $a$-axis of _oS_ 24
lie approximately on the same compression curve (Fig. 3, center), reflecting
the stacking correspondence in these directions between the two phases. The
$a$-axis of _mS_ 16 is systematically smaller than the $b$-axis of _oS_ 24\.
At any given pressure, _mS_ 16 shows a larger average thickness of the NaCl-
like layers, the greater the degree of dimerization the larger is the $ab$
plane and the shorter is the average distance between layers.

Figure 3:  Comparison of the pressure evolution of the unit-cell geometry of
_mS_ 16 and _oS_ 24: Volume per formula unit (top panel) and $a$, $b$, and $c$
lattice parameters (central panel). Blue symbols represent experimental _mS_
16 data, whereas red symbols show _oS_ 24 data, and the lines correspond to
calculated values. The $c$ parameters of _mS_ 16 and _oS_ 24 are shown
respectively divided by 2 and 3. Bottom panel: Calculated pressure evolution
of near-neighbour distances in _mS_ 16\. As a reference, the nearest-neighbor
Ga-P distance for the ideal NaCl phase is also provided (dashed black curve).
Figure 4:  Energy-volume relations (per formula unit) for the various
calculated structures. The energy is given with respect to the zero-pressure
state of the zb phase. In the insets, the enthalpies of the calculated phases
(given with respect to that of the NaCl phase) and the $p$-$T$ phase boundary
between _mS_ 16 and _oS_ 24 are provided.

Figure 4 shows the calculated energy-volume curves, $E(V)$, and enthalpy-
pressure curves, $H(p)$, for the different structural phases of GaP considered
in our _ab initio_ DFT calculations, including the novel _mS_ 16 phase as well
as the _Cmcm_ , $\beta$-Sn, _oS_ 24, and sc16 phases, among others [6]. These
results show _oS_ 24 and _mS_ 16 to be very close in energy/enthalpy (and
hence also density) in the whole range of compressions investigated, with _oS_
24 favored at the highest pressures and, viceversa, _mS_ 16 favored in the
lower pressure range, which is a qualitative trend consistent with the
experimental observations about the pressure conditions of the synthesis of
these two polymorphs and their intervals of existence. The calculated _mS_
16/_oS_ 24 coexistence pressure is around 32 GPa however this value should be
taken _cum grano salis_ given the closeness of their respective $E(V)$ curves
which results in a large uncertainty, with temperature favoring the highest-
pressure _oS_ 24 phase, see inset in Fig. 4 [32]. Notwithstanding this
unavoidable uncertainty related to the similarities of the two phases, all our
calculations consistently support the high-pressure stability of _mS_ 16
followed by that of _oS_ 24\. The calculated phonon frequencies of _mS_ 16
(see Fig. 5) show that it is dynamically stable at these high pressures.

Figure 5:  (Top) Calculated phonon bands along selected segments within the
Brillouin zone of the _mS_ 16 phase at $\sim$32 GPa. (Bottom) Idem for the
_oS_ 24 phase.

As recently reported [9], a narrow stability range for the sc16 polymorph is
found between zb and _mS_ 16\. On the other hand, around $\sim$30 GPa, the
enthalpies of the $\beta$-Sn and _Cmcm_ phases lie significantly higher
($\sim$70 meV pfu) than either _mS_ 16 or _oS_ 24, which rules $\beta$-Sn and
_Cmcm_ out as thermodynamically stable in the phase diagram of GaP, although
they have been previously considered in experiments [3, 5] and theoretical
studies [33, 34, 35]. The _Cmcm_ phase can be understood in terms of a simple
sliding pattern of alternate planes from the ideal NaCl structure [2], which
relates to an intrinsic TA($X$) zone-boundary mode softening of the NaCl-type
phases in III-Vs and the ensuing dynamical instability [36, 35]. Our
calculations are in agreement with this established view, however any possible
stability range for _Cmcm_ -GaP is completely obliterated when the dimerized
polymorphs are taken into consideration. As previously described, _mS_ 16 and
_oS_ 24 can also be viewed in terms of larger patterns of sliding NaCl planes
but while the structural flexibility of a four-plane (_mS_ 16) or six-plane
(_oS_ 24) stacking pattern allows the formation of interlayer dimers, in
various proportions, the simplest two-plane pattern of the _Cmcm_ structure
does not allow any such dimerization.

The _mS_ 16 phase is characterized by bonding anisotropy and distorted
coordination geometries (Fig. 2). The shortest bond of the structure is that
of the P-P dimer, measuring 2.217(6) Å at 30.7 GPa, a length typical of a
single covalent P-P bond [37]. P-Ga interatomic distances (Fig. 3, bottom) are
clustered between $\sim$2.35 and $\sim$2.5 Å, with a clear gap to longer
distances starting at about 2.7 Å, thus we consider bonding interactions all
the Ga-P distances up to about 2.5 Å. In addition to the dimer, P1 forms four
bonds with Ga1 atoms measuring between 2.36 and 2.48 Å, the coordination
geometry is strikingly similar to that of the corresponding atom in _oS_ 24
[7]. The non-dimerizing P2 atom is surrounded by six Ga atoms in a distorted
octahedral coordination with bonds lengths ranging between 2.36 and 2.50 Å.
The angles between P2 - Ga2 bonds within the L2 layer are close to 90∘,
whereas the angles involving off-plane bonds are as small as 67∘ and 72∘. In
_oS_ 24, the corresponding non bonding Ga atom lies at about 2.6 Å from P2,
and while all the angles between P2 and the four closest Ga atoms are close to
90∘, the bond angles formed with the apical Ga measure between 116∘ and 127∘,
defining a clear fivefold coordination. Therefore, unlike P1, the coordination
geometry of P2 shows some small but significant differences in the two
dimerized phases, resulting in a coordination number six for P2 in _mS_ 16 and
five for P2 in the higher-pressure phase _oS_ 24\. Hence the two layered GaP
phases show a rare exception to the general rule of increasing coordination
number with pressure in polymorphs.

The behavior under pressure of the P-P dimer that stems from the _ab initio_
calculations (see Fig.3, bottom) shows that it is remarkably rigid within the
range of observation of the _mS_ 16 phase, whereas its length _increases_
slightly at higher pressures, which is indicative of a slight weakening of the
homonuclear bond. On the other hand, the various intralayer Ga-P bonds follow,
unremarkably, approximately the same compressibility curve as that calculated
for the reference NaCl phase in the interval of interest.

Figure 6:  Calculated ELF isosurfaces for the valence $s$ and $p$ electrons of
the _mS_ 16 (top) and _oS_ 24 (bottom) phases, as seen along the $a$ axes of
their respective primitive cells (outlined in black). The value of the ELF was
set at 0.75.

Calculated ELF isosurfaces (valence only) for the two dimerized polymorphs are
shown in Fig. 6 [38]. Once again, these plots show that there are strong
similarities between the two phases and also subtle differences. The covalent
P-P bonding is clearly seen in both phases as well-defined large-ELF domains
between nearby P1 atoms. The ELF around P1 is nearly identical in the two
phases, the four Ga–P bonds produce a doughnut shaped ELF domain in the _ab_
plane. Fig. 6 highlights the subtle differences in ELF around the non-
dimerizing phosphorus in _mS_ 16 and _oS_ 24\. In the projection shown in Fig.
6 the ELF shows nearly perfect threefold symmetry around P2 in _oS_ 24, such
symmetry is broken in _mS_ 16, and the ELF domain at 0.75 protrudes in the
direction of the sixth Ga atom, supporting the occurrence of a weak bond.

In summary, two decades after the high-pressure polymorphism of the binary
semiconductors was considered a substantially settled matter, our findings
reveal the occurrence of a new structure type, adopted by the first high-
pressure metallic polymorph of gallium phosphide. We have unraveled a marked
complexity in the metallic phases, arising from strong bonding differentiation
and finely tuned pressure dependence of the degree of dimerization. It is
worth noting that the work presented here is not a mere redetermination of the
phases observed in previous studies and interpreted as _Cmcm_ [3, 5]. Without
heating, the complex ordering of _mS_ 16 and _oS_ 24 cannot be achieved (the
powder diffraction pattern would show more peaks than observed in the cited
studies). Small domains locally resembling _mS_ 16 or _oS_ 24 could however
form. Hence a distorted NaCl-like phase, locally ordered, might still be a
reasonable description for the metastable disordered material obtained upon
ambient-temperature compression of GaP.

The discovery of _mS_ 16, along with the previously-described _oS_ 24
polymorph of GaP, challenges previous views concerning the unfavorable role
played by bonds between like atoms in the energetics of the high-pressure
phases of the III-V and II-VI binary compounds and emphasizes the role of
NaCl-derived layered structures _with_ partial interlayer dimerization as a
relevant type of ordering in these compounds. The role of complex ordering
derived from the NaCl basic arrangement has been minimized [39], but it is
clear that while the total-energy difference between such polymorphs is small,
they are distinct phases and their occurrence is hereby explained. The layer
stacking allows for the accommodation of variable degrees of phosphorus
dimerization. Distinct ordering sequences of such layers, and therefore new
structures, are likely to occur in other octet compounds as well, an aspect
relevant to material design, as phases and properties might be finely tunable
with pressure.

###### Acknowledgements.

This research was partially supported by COMPRES, the Consortium for Materials
Properties Research in Earth Sciences, under NSF Cooperative Agreement EAR
1606856. AM acknowledges financial support through Projects
MAT2016-75586-C4-3-P (MINECO, Spain), PID2019-106383GB-C43 and
PID2022-138076NB-C44 (MCIN/AEI/10.13039/501100011033). The MALTA Consolider
Team Research Network is supported by Project RED2022-134388-T (MCIN, Spain).
Portions of this work work were conducted at GeoSoilEnviroCARS (GSECARS),
Sector 13 and High Pressure Collaborative Access Team (HP-CAT), Sector 16,
Advanced Photon Source (APS), Argonne National Laboratory. GSECARS is
supported by the National Science Foundation - Earth Sciences (EAR-1634415)
and DOE-GeoSciences (DE-FG02-94ER14466). HPCAT operations are supported by
DOE-NNSA’s Office of Experimental Sciences. The Advanced Photon Source is a
U.S. Department of Energy (DOE) Office of Science User Facility operated for
the DOE Office of Science by Argonne National Laboratory under Contract No.
DE-AC02-06CH11357 Lawrence Livermore National Laboratory is operated by
Lawrence Livermore National Security, LLC, for the U.S. Department of Energy
(DOE), National Nuclear Security Administration under Contract DE-
AC52-07NA27344. Use of the COMPRES-GSECARS gas loading system was supported by
COMPRES under NSF Cooperative Agreement EAR-1606856 and by GSECARS through NSF
grant EAR-1634415 and DOE grant DE-FG02-94ER14466.

## References

  * Homan  _et al._ [1975] C. G. Homan, D. P. Kendall, T. E. Davidson, and J. Frankel, Solid State Commun. 17, 831 (1975).
  * Nelmes and McMahon [1998] R. J. Nelmes and M. I. McMahon, Semicond. Semimet. 54, 145 (1998).
  * Nelmes  _et al._ [1997] R. J. Nelmes, M. I. McMahon, and S. A. Belmonte, Phys. Rev. Lett. 79, 3668 (1997).
  * Baublitz and Ruoff [1982] M. Baublitz and A. L. Ruoff, J. Appl. Phys. 53, 6179 (1982).
  * Aquilanti  _et al._ [2007] G. Aquilanti, H. Libotte, W. Crichton, S. Pascarelli, A. Trapananti, and J.-P. Itié, Phys. Rev. B 76, 064103 (2007).
  * Mujica  _et al._ [2003] A. Mujica, A. Rubio, A. Muñoz, and R. J. Needs, Rev. Mod. Phys. 75, 863 (2003).
  * Lavina  _et al._ [2018] B. Lavina, E. Kim, H. Cynn, P. F. Weck, K. Seaborg, E. Siska, Y. Meng, and W. Evans, Inorg. Chem. 57, 2432 (2018).
  * Nelmes and McMahon [1995] R. J. Nelmes and M. I. McMahon, Phys. Rev. Lett. 74, 106 (1995).
  * Lavina  _et al._ [2022] B. Lavina, E. Zanardi, A. Mujica, H. Cynn, Y. Meng, J. Smith, M. Kong, and Y. Lee, High Pressure Research 42, 294 (2022), https://doi.org/10.1080/08957959.2022.2113874 .
  * Mao and Bell [1976] H. Mao and P. Bell, Science 191, 851 (1976).
  * Dewaele  _et al._ [2004] A. Dewaele, P. Loubeyre, and M. Mezouar, Phys. Rev. B 70, 094112 (2004).
  * Rivers  _et al._ [2008] M. Rivers, V. B. Prakapenka, A. Kubo, C. Pullins, C. M. Holl, and S. D. Jacobsen, High Press. Res. 28, 273 (2008).
  * Meng  _et al._ [2015] Y. Meng, R. Hrubiak, E. Rod, R. Boehler, and G. Shen, Rev. Sci. Instrum. 86, 072201 (2015).
  * Prescher and Prakapenka [2015] C. Prescher and V. B. Prakapenka, High Press. Res. 35, 223 (2015).
  * Diffraction [2018] R. O. Diffraction,  _CrysAlisPro Software (version 1.171.39-44), Rigaku Corporation, Oxford, UK_ (2018).
  * Sheldrick [2008] G. M. Sheldrick, Acta Crystallogr. A 64, 112 (2008).
  * Blöchl [1994] P. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Kresse and Hafner [1993] G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Perdew  _et al._ [1997] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 78, 1396 (1997).
  * Not [a] So as to check the effect of the XC functional, we did further calculations using the so-called PBEsol GGA functional [24, *Perdew:2008aa_erratum], and the effect of frozing the Ga semicore $d$ electrons was also checked. The vasp total-energy calculations for the PBE XC functional with explicit treatment of the Ga $d$ electrons were double-checked with Quantum espresso [26], using recommended potentials and cutoff. The intercomparison among the results of the different calculations will be given elsewhere [29] and only referred briefly here.
  * Perdew  _et al._ [2008] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Phys. Rev. Lett. 100, 136406 (2008).
  * Perdew  _et al._ [2009] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Phys. Rev. Lett. 102, 039902 (2009).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. Fabris, G. Fratesi, S. de Gironcoli, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, J. Phys.: Condens. Matter 21, 395502 (2009).
  * Togo and Tanaka [2015] A. Togo and I. Tanaka, Scripta Materialia 108, 1 (2015).
  * Alfè [2009] D. Alfè, Comput. Phys. Commun. 180, 2622 (2009).
  * [29] E. Zanardi  _et al._ , in preparation.
  * Lavina  _et al._ [2013] B. Lavina, P. Dera, and Y. Meng, JOVE , e50613 (2013).
  * Not [b] The intralayer L2 puckering is absent in the description of _oS_ 24-InSb by Nelmes and McMahon [2].
  * Not [c] Due to their closeness in energy, not only the temperature can have an effect on the _mS_ 16/_oS_ 24 stability boundary but also the level of approximation of the calculation, and in particular the choice of XC functional. For example, considering the PBEsol functional [24, *Perdew:2008aa_erratum] results in a calculated coexistence presssure of $\sim$26 GPa.
  * Mujica and Needs [1997] A. Mujica and R. J. Needs, Phys. Rev. B 55, 9659 (1997).
  * Mujica  _et al._ [1998] A. Mujica, A. Muñoz, and R. J. Needs, Phys. Rev. B 57, 1344 (1998).
  * Ozoliņš and Zunger [1999] V. Ozoliņš and A. Zunger, Phys. Rev. Lett. 82, 767 (1999).
  * Mujica and Needs [1996] A. Mujica and R. J. Needs, J. Phys.: Condens. Matter 8, L237 (1996).
  * Pöttgen  _et al._ [2011] R. Pöttgen, W. Hönle, and H. G. von Schnering, Phosphides: Solid-state chemistry, in  _Encyclopedia of Inorganic and Bioinorganic Chemistry_ (American Cancer Society, 2011).
  * Momma and Izumi [2011] K. Momma and F. Izumi, J. Appl. Crystallogr. 44, 1272 (2011).
  * Kelsey and Ackland [2000] A. A. Kelsey and G. J. Ackland, J. Phys.: Condens. Matter 12, 7161 (2000).


# Neutron Production in The Scattering S250 Mevion Proton-Therapy System

Nima Tatari
Washington University in St. Louis, Department of Physics
<EMAIL_ADDRESS>
##  1 Introduction

A study of the neutron production in a single-room Mevion S250 passive
scattering proton-therapy system is performed. The proton-therapy system has
been in operation at Washington University in St. Louis since 2013 and was
decommissioned in 2020. Geant4 Monte Carlo toolkit has been used to calculate
neutron fluence spectra from various nozzle components. Each fluence spectrum
is discussed based on the cross-sections of the materials contributing to
secondary neutron production.

##  2 Materials and methods

Several parameters in the simulation consist of the geometrical dimensions of
the room, synchrocyclotron, and nozzle components. The nozzle includes a first
scatterer, a range modulation wheel, the main range shifter, a second
scatterer, a beam collimation system, a brass aperture, and a range
compensator. A water phantom models the patient’s body. The geometry and
physics models in the simulation are discussed in the following section.

###  2.1 Geometry

Modeling a geometry in Geant4 is based on the dimensions, compositions, and
densities of the objects in the environment of radiation transport. The
treatment room, synchrocyclotron, beam-delivery system, and water phantom
define the geometry, each part contributing to the total neutron production
based on different physical processes. The neutrons from the room are created
by previously generated neutrons in the nozzle. The contribution of the
protons that are lost in the synchrocyclotron’s extraction stage to the
production of secondary neutrons is calculated. Figures 1 (a,b,c) show the
geometry of the room containing the proton-therapy system. The whole vault is
modeled with a 13.5$\times$15 m2 floor and with a 15 m height.

Figure 1: a) 3D view, b) top view, and c) side view of the treatment room
geometry design in the Geant4 model of the Mevion S250 scattering proton
system.

The components of the nozzle are shown in Figure 2. The first scatterer is
positioned after the beam extraction point from the synchrocyclotron.

Figure 2: Nozzle components in a Geant4 model of Mevion S250 scattering proton
system.

The scattering proton machine is modeled with a primary mono-energetic proton
beam of 250 MeV and is spread laterally to get an even dose throughout a
specific area and depth. Rotation of the RMW superposes the modulated proton
beams with different energies. The result is a spread-out Bragg peak (SOBP),
with a uniform proton dose curve, inside the water phantom.

Broadening of the proton beam in the model is achieved through the second
scatterer. The unwanted part of the beam is trimmed by the collimation system
which is modeled by highly attenuating materials such as stainless steel and
brass with large stopping power. The collimation system is modeled with a
stainless-steel cylinder with a 3 cm thickness, 9.75 cm inner radius, and 36
cm height. A similar structure with a thickness of 2 cm, 15 cm inner radius,
and 52.5 cm height was used for the collimator with the same material. In
addition to the outer cylinder, the collimator contains two inner rings at its
entrance and exit sides. Each ring has an inner radius of 12.5 cm with the
same material as the outer collimator. A brass aperture with a thickness of
6.84 cm is placed next to the collimator. The center of the brass aperture
contains a square opening.

During the actual treatment, the fine-tuning of the output beam energy is done
by a patient-specific lucite compensator that varies the depths of the beam at
different positions to shield the healthy tissues from the proton beam. The
compensator is modeled with a lucite cylinder of 5 cm height placed as the
last step of the nozzle.

###  2.2 Geant4 Monte Carlo simulation toolkit

Geant4 (GEometry ANd Tracking) has been developed by CERN as a general-purpose
Monte Carlo (MC) simulation software for particle transport through matter.
Two fundamental aspects of Geant4 are handling the geometrical shapes and
tracking various particles through objects made of different materials.
Detailed radiation transport calculations are performed by the toolkit for
that purpose. The building blocks of Geant4 simulations are G4Step objects
from which information can be stored. Particle transport involves a 3D
simulation of the interactions of particles irradiated through materials with
specific geometries. In this work, the primary particle source is protons.

##  3 Results and discussion

The proton PDDs and neutron spectra are benchmarked with measured data. The
neutron fluence spectra calculated by Geant4 form the basis of the presented
analysis. The properties of the spectra are described based on the neutron
cross-section data obtained from ENDF datasets [1].

Protons form a well-defined dose deposition peak as they travel through
matter, due to the stopping power. The simulation PDD is obtained by emitting
a primary proton beam with 250 MeV energy, 2 MeV energy bandwidth, and 1 mm
spot size.

Figure 3: Matching the Geant4-calculated pristine percentage depth dose to the
experimental data. Both PDDs correspond to the first step of the range
modulation wheel. The benchmarked parameters include the thickness of the
range shifter (affecting the range in water), the energy width of the primary
proton beam (affecting the range straggling), and the spot size of the initial
beam (affecting the entrance dose height).

Each range modulation step delivers a specific proton beam energy to the
isocenter with a corresponding pristine peak depth. Rotation of the range
modulation wheel leads to a simultaneous superposition of the proton beams
with various energies, and depending on the design of the angular extent of
each step that corresponds to their weights, leads to a spread-out Bragg peak
or SOBP.

The main characteristics of an SOBP are range and modulation. The range is
defined as the distance from the surface of the water phantom to the distal
90% proton dose. Modulation is defined as the distance from the proximal 95%
dose to the distal 90% [2]. Specifically, for the treatment setup of the whole
brain the range of 17 cm and modulation of 16 cm has been used in the work of
Howell et al. [2]. The same SOBP has been simulated in this work for
benchmarking the neutron calculations.

Simulation of an SOBP requires a linear combination of several pristine peaks
with different ranges. Each pristine peak is obtained by positioning a step of
the range modulation wheel in the primary proton beam path. The final SOBP has
a 17 cm range and 16 cm modulation. The first 15 pristine peaks delivered by
the 15 are used as the base functions in a singular value decomposition (SVD)
method to obtain the composition weights. Figures 4 and 5 show the normalized
pristine peaks and the SOBP with a range of 17 cm, respectively.

Figure 4: Normalized proton pristine peaks corresponding to 15 steps of RMW.
Figure 5: Constructing the SOBP from the 14 pristine Bragg peaks using the
singular value decomposition method.

To obtain the agreement between the transverse proton dose profiles, the
parameters that define the shape of the second scatterer are modified. The
second scatterer consists of bell-shaped, concave Lexan, and thin, convex lead
components. Contracting or widening these parts determines the uniformity of
the transverse proton dose profile.

Figure 6: Transverse matching of the simulated proton dose profile with the
measurement data (Zhao et al. 2016 [3]).

In this article, neutron dose-equivalent H* normalized by the treatment dose
at the isocenter D is analyzed. The dose D is defined as energy per unit mass,
and as a point quantity, it has a specific value at each geometric point. Dose
to the isocenter in a passive system per primary proton is in the range of
5.5$\times$10-13 Gy/pr to 6.5$\times$10-13 Gy/pr. This value remains roughly
stable regardless of the field size and delivered proton energy.

####  3.0.1 Benchmarking of the neutron fluence and neutron dose spectra

In a study by Howell et al [2] extended-range Bonner spheres (ERBS) have been
experimentally employed and the generated data unfolded to obtain the neutron
fluence spectra at 50 cm transverse to the isocenter in a scattering Mevion
proton-therapy system. The same geometry is set up with a beam angle in the
horizontal plane. Simulations corresponding to different energy outputs are
performed.

The neutron fluence spectra, corresponding to each range modulation step, are
weighted and summed according to the SOBP weights. Figure 7 indicates the
absolute values of the neutron fluence spectra per unit treatment (proton) Gy
to the isocenter per unit lethargy. Integration of the curves results in the
total fluence per proton Gy. The measured value by the Bonner sphere technique
has been reported 1.52$\times$107 n cm-2 Gy-1 [2]. The corresponding Geant4
total neutron fluence is 1.45$\times$107 n cm-2 Gy-1. The values agree within
4.6%.

Figure 7: Benchmarking of the neutron fluence per treatment proton Gy to the
isocenter per unit lethargy. The red curve shows the result of an ERBS
measurement by Howell et al. The blue curve shows the Geant4 simulation result
of the present work. Figure 8: Benchmarking of the neutron dose-equivalent per
treatment proton Gy to the isocenter per unit lethargy. The red curve shows
the result of an ERBS measurement by Howell et al. The blue curve shows the
Geant4 simulation result of the present work.

###  3.1 Analysis of the neutron sources using fluence ($\phi$(E)/D)

Characterizing the source of the neutrons incident onto the water phantom is
performed by decomposing the neutron fluence spectrum into several
contributions. The cross-section data are used to describe the behavior of the
spectra. The synchrocyclotron is made of borated polyethylene and stainless-
steel with radii of 74 to 80 cm and 80 to 90 cm respectively and 4 m height is
positioned horizontally behind the nozzle.

To account for the proton extraction step inside the synchrocyclotron, a cone
region with 90° opening is considered in its interior part. According to
Figure 9 (a), the neutrons produced in the cyclotron contain a thermal-to-
evaporation peak ratio close to 1. The borated-polyethylene layer has a large
thermalizing effect on the neutrons. The relatively large continuum of the
neutrons between the evaporation and thermal regions could be explained based
on the multiple objects that a neutron might interact with along its
trajectory.

(a) Neutrons from the primary protons lost inside the synchrocyclotron (a)
Neutrons from the primary protons sent through nozzle Figure 10: Neutron
fluence per unit treatment Gy incident onto the water phantom from neutrons
generated in the synchrocyclotron made of stainless-steel and borated-
polyethylene.

The first object that the primary protons hit is the first scatterer made of
lead. Neutrons are produced as a result of direct interactions of protons with
lead and move in a forward-peaked momentum distribution.

Figure 11: Neutron fluence incident onto the water phantom from neutrons
generated in the FS made of lead.

High-energy protons interact with the different steps of the range modulation
wheel containing lead. The calculated neutron spectrum at the water phantom
position is similar to the neutron spectrum generated in the first scatterer.

(a) RMW step 1 (a) RMW step 7 (a) RMW step 14 Figure 14: Neutron fluence
incident onto the water phantom from neutrons generated in the 1st, 7th, and
14th steps of the RMW

Using step 7 of the RMW increased the fast neutron peak by a factor of 5, as
can be seen in Figure 11 (b). The evaporation and thermal peaks remained in
the same order. The increase could be explained by the additional layer of
Lexan with 1.02 cm thickness and graphite with 2.22 cm thickness. Charge
exchange can happen frequently in light, symmetric nuclei due to direct
reactions with the protons. This type of reaction in carbon 12C(p,n)12N could
generate fast neutrons.

Figure 11 (c) shows the neutron fluence spectrum detected at the water phantom
generated by the 14th step of the RMW. This step contains 1.95 cm thick Lexan
and 6.9 cm thick graphite. The fast neutron peak increased by a factor of 1.5
compared to the results with the 7th step of the RMW.

####  3.1.1 Main range shifter (MRS)

MRS was modeled as a disk with 4 cm thickness made with carbon [4]. A trend
for the neutron fluence similar to that of the RMW steps 7 and 14 was observed
here, as shown in Figure 12. The protons passing through the RMW interacted
with carbon in the MRS and generated fast neutrons. The magnitude of the fast
neutron component was similar to that generated in the RMW step 14 with a 6.9
cm thickness of carbon.

Figure 15: Neutron fluence incident onto the water phantom from neutrons
generated in the MRS made of carbon.

####  3.1.2 Second scatterer (SS)

The SS is modeled by an outer brass ring and an inner part made of Lexan and
lead. The proton beam does not mainly interact with the brass ring of the SS.
Its Lexan components have a total mass of 0.52 kg. The lead has a mass of 0.06
kg. The large fast neutron peak in comparison to the evaporation peak could be
associated with the carbon elements in the Lexan. Similar to the RMW and MRS,
charge exchange (p,n) and knock-out (p, np) due to inelastic proton
interactions could be the main processes for fast neutron production. The
fluence of the water phantom is indicated in Figure 14.

Figure 16: Neutron fluence incident onto the water phantom from neutrons
generated in the second scatterer made of lead, Lexan, and brass.

####  3.1.3 Collimation system

The collimation system consists of a beam shielding cylinder and a collimator
tube both made of stainless steel, a mixture of 74% iron, 18% chromium, and 8%
nickel atoms as described in the NIST material database of Geant4 [5].

The main interactions in these components are initiated by the primary protons
that have been degraded down to the energy region of 160 MeV – 120 MeV
depending on the step used in the energy modulation wheel. Protons are
tangential to the inner edges of these tubes due to the transverse scattering
that is done by the first and the second scatterers. The fast neutrons are
generated within the collimation system, as indicated in Figure 14.

(a) Beam shield neutrons (a) Collimator neutrons Figure 18: Neutron fluence
incident onto the water phantom from neutrons generated in the collimation
process.

Fast neutrons mainly correspond to the charge exchange (p,n) and knock-out (p,
np) reactions. The average of the cross-section for (p, nx), x representing
any by-product nucleon or fragment aside from the secondary neutron, with a
stainless-steel target is 0.7 barn. For 12C target, the cross-section of
(p,nx) is of the order 0.3 barn. Given this value for 12C cross-section of (p,
nx), large values are observed for fast neutron peaks for the RMW and MRS that
contain carbon. The beam shield causes approximately 1/10 of fast neutrons
incident to the water phantom compared to the MRS, although the cross-section
of fast neutron production in the beam shield is more than twice larger than
that of the MRS. This could be explained by a reduced number of protons that
interact with the beam shield, in addition to the peripheral position of the
shield with respect to the proton beam, as compared to the central position of
the MRS in the beam path. The latter fact leads the fast neutrons from the MRS
to have a larger overlap with the surface of the water phantom.

The evaporation peak is due to neutrons that are emitted from excited nuclei.
Fast neutron production processes, including (p,n) reactions, could form
compound nuclei. A compound nucleus is “an unstable nucleus formed by the
coalescence of an atomic nucleus with a captured particle” [6]. Such excited
systems will undergo subsequent neutron evaporation, with an exception of the
clustered systems such as 12C or 16O. Clustered systems usually break up into
smaller clusters such as the break-up of a compound 12C into three alphas.

In addition to proton reactions, compound systems could be formed by neutron
nonelastic reactions that can lead to further neutron evaporation. Neutrons
are thermalized due to elastic reactions with light targets. Hydrogen provides
a high thermalizing ratio per interaction. Thermal neutrons are captured by
nuclei. Neutron capture excites the nuclei and causes gamma radiation. Due to
radiative capture, neutron fields are mixed with gamma. According to the
measurements performed with a WENDI-2 chamber, the absorbed dose equivalent
from gamma is less than 3% of that from neutrons.

####  3.1.4 Brass aperture (BA)

The neutron spectrum incident on the water phantom (WP) from those neutrons
that have been initially generated in the BA shows a low thermal-to-
evaporation peak ratio. The integral of the neutron fluence incident onto the
WP from the BA is significantly larger than the other components. This is due
to the small distance between the BA and the WP. The forward-peaked fast
neutrons have a large probability to hit the WP surface because of the short
distance.

Also, the low thermal-to-evaporation peak ratio in the spectrum shown in
Figure 15 could be explained by the small amount of material between the BA to
WP which reduces the probability for the evaporation neutrons to undergo
elastic reactions with light nuclei and thermalize.

Brass is made of 67.4% copper atoms, 32.1% zinc atoms, and less than 0.5%
mixture of iron, lead, and tin atoms. The primary proton incident on the brass
aperture is mostly in this energy range. These protons have been degraded in
energy by the first scatterer, the range modulation wheel, and the second
scatterer.

Fast neutron production rate in brass per high-energy protons in the range of
120 to 160 MeV impinging on the brass aperture is expected to be in the same
order of magnitude as that in 12C which was discussed in the subsection on the
RMW. The elements in the brass aperture and the graphite in the RMW have a
$\sigma$ of about 0.4 and 0.7 barns for nonelastic (p, nx) reactions,
respectively.

Figure 19: Neutron fluence incident onto the water phantom from neutrons
generated in the brass aperture.

The nonelastic (n,nx) cross-section values for various neutron energies in the
MeV range, could be used to explain the evaporation of neutrons produced in
the brass aperture. Values of the cross-sections for (n,nx) reactions in brass
between incident neutron energies of 1 MeV to 10 MeV are comparable in
magnitude to the cross-section of (p,n) reactions in 12C with $\sigma$=0.7
barn. It was observed that (p,n) reaction in 12C leads to a large fast neutron
fluence incident on the water phantom. With a similar analogy, we might
conclude a large probability for the evaporation process based on the $\sigma$
values for the (n, nx) process in brass. That could explain the large
evaporation peak observed in Figure 15.

####  3.1.5 Ceiling and floor

The ceiling and floor were simulated by 15 cm and 20 cm thick concrete
respectively, with a 150 cm vertical distance from the isocenter each.

About half of the atoms in concrete are oxygen. The most probable reactions
that could happen on the floor would be from incident neutrons with various
energies, lower than the maximum primary proton energy. Neutrons that reach
the floor have interacted several times elastically or inelastically and have
attenuated in energy. Based on this fact, the absence of the fast neutron peak
in figures 16 (a) and (b) could be explained.

Neutrons with energies below 1 MeV have a large cross-section of interaction
with hydrogen nuclei. Hydrogen forms 30% of the total nuclei in concrete. The
attenuation ratio of neutron energy by elastic interactions with hydrogen per
collision is 0.5, on average. The thermal neutron peaks, as shown in Figure
16, were close to the evaporation peak in magnitude. Also, a relatively large
continuum of neutrons between thermal and evaporation energies was produced by
the simulations.

(a) Ceiling neutrons (a) Floor neutrons Figure 21: Neutron fluence incident
onto the water phantom from neutrons generated in the treatment room.

####  3.1.6 Internal neutrons

The calculated fluence for the internal neutrons in the water phantom has been
done by recording the energy of each neutron once produced inside the water
phantom. Similar to the fluence diagrams from other components, the resulting
internal neutron histogram has been divided by the average cross-sectional
area of the water phantom.

According to Figure 17, no neutrons were generated in the thermal region. This
is due to the fact that the neutrons are recorded at the initial step of their
products before being thermalized. Also, the magnitude of the fast and
evaporation neutrons are close, with the fast neutron peak being slightly
larger.

Figure 22: Internal neutron fluence from neutrons generated in the water
phantom.

####  3.1.7 Comparison of neutron fluence and neutron fluence per unit mass of
the nozzle items

Figure 18 shows the neutron fluence spectra per unit treatment Gy to the
isocenter per unit lethargy for all the nozzle components. According to Figure
18, the largest contribution to the external fluence received at the water
phantom comes from the brass aperture. Neutron fluence from other items in the
nozzle is about an order of magnitude less than the brass aperture’s
contribution. The abbreviations in the key in Figures 18 and 2.19 read as: FS=
first scatterer, RMW= range modulation wheel, MRS= main range shifter, SS=
second scatterer, SHLD= beam shield, CLMT= beam collimator, BRSAP= brass
aperture, CMPNS= Lucite compensator.

Figure 23: Neutron fluence spectra incident onto the water phantom from
neutrons generated in the nozzle’s items. The contribution from the brass
aperture dominates.

Figure 19 shows the neutron fluence spectra per unit treatment Gy to the
isocenter per unit lethargy for all the nozzle components per their unit mass.
According to Figure 19, the largest contribution to the external fluence
received at the water phantom per unit mass of the items comes from the first
scatterer and the first step of the range modulation wheel. This fact is due
to the lead component in these items that contributes largely to fast neutron
production. Neutron fluence per unit mass from other items in the nozzle is
about an order of magnitude less than the brass aperture contribution.

Figure 24: Neutron fluence spectra incident onto the water phantom from
neutrons generated in the nozzle’s items per unit mass of the corresponding
item. The first scatterer and the range modulation wheel dominate.

###  3.2 Neutron profiles inside the treatment room with and without the water
phantom

For this section, the profiles of the neutron dose-equivalent per treatment Gy
have been plotted and the general trends of the neutron dose behaviors at
various positions inside the treatment room have been obtained. The
simulations for calculating the profiles in Figures 22 and 2.23 were done with
a water phantom with dimensions 35$\times$35$\times$35 cm3 placed at the
isocenter. The first RMW step was used and the brass aperture opening was set
to 20$\times$20 cm2.

####  3.2.1 Neutron dose profile perpendicular to the beam direction

Figure 20 (a) plots the profiles for the absorbed neutron dose-equivalent per
treatment Gy to the isocenter for different heights measured from the
isocenter. The profiles show the behavior of H*/D as a function of distance
along the x-axis, which is the axis parallel to the synchrocyclotron cylinder.
As can be seen from Figure 20 (a), an approximately symmetric behavior was
observed for H*/D. As seen for the heights 80 cm and 120 cm, some deviations
from the symmetric behavior could be due to the asymmetry inside the room. The
entrance of the maze into the treatment room introduces a source of asymmetry
that might affect the neutron dose-equivalent profiles. The maximum value
observed for H*/D was 9.2 mSv/Gy, which was due to the neutrons counted at the
isocenter. A comparison between the two cases where a water phantom was placed
at the isocenter and where it was not there has been provided in Figure 20
(b).

(a) At various heights, with water phantom. (a) With and without water
phantom. Figure 26: Neutron dose-equivalent profiles per treatment proton dose
at the isocenter in mSv/Gy unit. The values have been plotted as a function of
transverse distance from the isocenter across the x-axis for various heights.
A water phantom with dimensions 35$\times$35$\times$35 cm3 has been placed at
the isocenter for the results shown in figure (a). The profiles in figure (b)
correspond to the H*/D values across the x-axis at height y=0 from the
isocenter with a water phantom (blue) and without a water phantom (red). The
first RMW step was used and the brass aperture opening was set to 20$\times$20
cm2.

####  3.2.2 Neutron dose profiles along the beam direction

Figure 21 plots the profiles for the absorbed neutron dose-equivalent per
treatment Gy to the isocenter along the z-axis at a height 50 cm below the
isocenter and parallel to the proton beam. The blue curve corresponds to a
simulation with a water phantom placed at the isocenter and the red curve
corresponds to a simulation without a water phantom. As can be seen from
Figure 21, the general trend of H*/D was a decrease in the neutron dose with
increasing distance from the nozzle. An exception to this trend was found at
the isocenter, where the water phantom was positioned. The internal neutrons
generated inside the water phantom shifted the profiles upward at the region
around z=0. In Figure 21, the point z=0 indicates the isocenter, and the
primary protons are generated at the point z=2 m toward the negative
z-direction. The calculated value of H*/D at the isocenter was 4 mSv/Gy when
the water phantom front side was positioned and 3.5 mSv/Gy when it was
removed. The largest difference between the two cases was observed at z=0.5 m,
close to the backside of the water phantom. At that position, 2.5 times larger
H*/D was scored for the configuration with the water phantom. At z=2 m, the
extraction point where the primaries are generated, the value of H*/D for both
cases was about 19 mSv/Gy. Finally, closer to the front wall at z=-6.5 m, a
slight increase in the H*/D was observed for the case without the water
phantom. That could be associated with the evaporation neutrons that back-
scatter from the front wall due to the high energy protons that reach the wall
because of the absence of blocking by the water phantom.

Figure 27: Neutron dose-equivalent profile per treatment proton dose at the
isocenter in mSv/Gy unit. The values have been plotted as a function of
distance along the beam direction on an axis 50 cm below the isocenter.
Isocenter is at z=0, positive z is toward the nozzle, and the proton beam is
directed toward the negative z-direction. A water phantom with dimensions
35$\times$35$\times$35 cm3 has been placed at the isocenter. The first RMW
step was used and the brass aperture opening was set to 20$\times$20 cm2.

##  4 Conclusions

A study of neutron production in the scattering Mevion proton-therapy system
is performed. The influence of several parameters on neutron production in the
scattering Mevion machine was considered and analyzed. The Geant4 code was
calibrated using previously reported data. The results for the total fluence
and total neutron dose-equivalent matched the extended-range Bonner sphere
(ERBS) measurements within 4.6% and 0.25%, respectively. In addition, the
neutron dose-equivalent calculations at 5 transverse points with comparable
beam and nozzle parameters to a previous TOPAS study were fit with a
difference of an order of 10%.

The simulations with Geant4, however, did not agree with the WENDI-2
measurements in the same study that provided the TOPAS results. Moreover,
other reported values of H*/D exist in the literature that my data has not
agreed with. These differences could be associated with the various types of
machines, the different proton source average energies, energy spreads, and
spot sizes, uncertainties regarding the proton extraction ratios in the
accelerators, the geometrical uncertainties in the items’ dimensions and
variations in the material densities, and to the uncertainties in the cross-
section libraries.

The analysis of neutron production in the scattering machine consists of
showing several diagrams of the neutron fluence spectra per treatment Gy to
the isocenter received from various items in the nozzle and those in the room
onto the water phantom, which models a human body. Data on the percentages of
thermal, intermediate, evaporation, and fast neutron production were provided
for each contributing item to the neutron dose. The neutron fluence spectra
per unit mass of each item are shown. Correspondingly, the largest neutron
production per mass of the objects is from the lead in the first scatterer and
in the RMW.

The pie charts of the total H*/D from different nozzle components and room
items are obtained and plotted. An analysis is done by changing the RMW step
at a fixed field size and changing the field size at a fixed RMW step. The
corresponding effects on the values of H*/D are observed and studied. It is
concluded from the pie charts that about 60% of the neutron H*/D is caused due
to the internal neutrons when the 1st RMW step is used. This ratio decreased
to 50% for 7th RMW step. Using the 14th step of the RMW further decreases this
ratio to about 25%.

A decrease of H*/D from the neutrons generated in the brass aperture is
observed as a result of rotating the RMW from the first step upward. For the
1st RMW step the contribution of the brass aperture is 2.7 mSv/Gy, for the 7th
RMW step it is 1.4 mSv/Gy, and for the 14th RMW step it is 0.28 mSv/Gy. In
this analysis, the RMW output beam energy is lowered from 227 MeV to 165 MeV,
while keeping the brass aperture opening constant with a 10 cm side length. As
a result of this shift in the RMW output beam energy, the absorbed neutron
H*/D to the water phantom decreases by about a factor of 10.

The profiles of the total, i.e. integrated, H*/D values as functions of
position inside the treatment room is obtained. A general decreasing trend of
H*/D is observed from the proton extraction point toward the isocenter and
beyond. The water phantom is found to cause an increase in the H*/D profile
around the isocenter. This fact is discussed to be due to the internal
neutrons that move outward from the water phantom. In conclusion, a set of the
expected values is presented for the absorbed neutron dose-equivalent per unit
treatment Gy to the isocenter for the scattering Mevion proton machine at
various positions, and as a function of various parameters. The H*/D values
could be used to obtain a better understanding of the difference between the
neutron production in the scattering Mevion S250 system and the scanning
Mevion S250i system that was studied in Chapter 3.

## References

  * [1] “ENDF Neutron Cross-Sections”  URL: https://www-nds.iaea.org/exfor/endf.htm
  * [2] Rebecca M. Howell et al.  “Measured Neutron Spectra and Dose Equivalents from a Mevion Single-Room, Passively Scattered Proton System Used for Craniospinal Irradiation”  In _Int. J. Radiat. Oncol. Biol. Phys._ 95, (2016), pp. pp 249–257  DOI: 10.1016/j.ijrobp.2015.12.356
  * [3] Tianyu Zhao et al.  “Commissioning and Initial Experience with The First Clinical Gantry-Mounted Proton Therapy System”  In _J. Appl. Clin. Med. Phys._ 17, (2016), pp. pp 24–40  DOI: 10.1120/jacmp.v17i2.5868
  * [4] Milad Baradaran-Ghahfarokhi et al.  “A Monte Carlo-based Analytic Model of Neutron Dose Equivalent for a Mevion Gantry-Mounted Passively Scattered Proton System for Craniospinal Irradiation”  In _Med. Phys._ 47, (2020), pp. pp 4509–4521  DOI: 10.1002/mp.14299
  * [5] “NIST Materials”  URL: https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/Appendix/materialNames.html
  * [6] “Merriam Webster”  URL: https://www.merriam-webster.com/dictionary/compound%20nucleus