Transfer line scattering model of therapeutic hadron beams and applications to nozzle and gantry optimization

The field of hadron therapy is growing rapidly with several facilities currently being planned, under construction or in commissioning worldwide. In the ‘‘active scanning’’ irradiation technique, the target is irradiated using a narrow pencil beam that is scanned transversally over the target while the penetration depth is altered with the beam energy. Together, the target dose can thereby be conformed in all three dimensions to the shape of the tumor. For applications where a sharp lateral beam penumbra is required in order to spare critical organs from unwanted dose, beam size blowup due to scattering in on-line beam diagnostic monitors, air gaps and passive elements like the ripple filter must be minimized. This paper presents a model for transverse scattering of therapeutic hadron beams along arbitrary multislab geometries. The conventional scattering formulation, which is only applicable to a drift space, is extended to not only take beam optics into account, but also non-Gaussian transverse beam profiles which are typically obtained from the slow resonant extraction from a synchrotron. This work has been carried out during the design phase of the beam delivery system for MedAustron, an Austrian hadron therapy facility with first patient treatment planned for the end of 2015. Irradiation will be performed using active scanning with proton and carbon ion beams. As a direct application of the scattering model, design choices for the MedAustron proton gantry and treatment nozzles are evaluated with respect to the transverse beam profile at the focal point; in air and at the Bragg peak.


I. INTRODUCTION
The goal of any radiotherapy, be it with hadron beams or photons, is to expose the tumor to a lethal dose, while minimizing the dose given to healthy tissue [1].Compared to photon beams, beams of ions are attractive due to the following two properties: (1) the finite, well-defined penetration depth, which implies that little or no dose is given beyond the tumor, and (2) the Bragg peak behavior of the depth-dose curve, where each particle deposits a high dose at the end of the range.Together, these advantages allow for higher dose in the target and/or lower dose at the entrance channel.
There are many cases where a narrow or sharp transverse beam profile is clinically motivated in order to spare critical organs or healthy tissue abutting the tumor from an unwanted dose: the sharper the beam profile, the better the dose can be conformed to the target.While there are intrinsic limitations to the minimal beam size which can be delivered by the accelerator, in the case of therapeutic low-energy protons, the main limitation is often scattering in elements upstream of the patient: vacuum windows, monitors, passive elements, and air gaps.In order to estimate the beam size at the focal point including these scattering effects, a semiempirical scattering model [2] has been extended to take both beam optics as well as non-Gaussian beam profiles into account.This model allows for rapid evaluation of different design options without the need for time-consuming tracking simulations and/or Monte Carlo simulations of the entire beam-line geometry.Since the model is not limited to thin scatterers, evaluation of the quality of the beam profile deep into the target is straightforward.
In this paper, the scattering model is described and then used to evaluate various different aspects of the nozzles and proton gantry of the Austrian hadron therapy facility MedAustron.This synchrotron based facility is currently under construction in Wiener Neustadt and will utilize beams of protons (60-250 MeV) and carbon ions (120-400 MeV=n) for tumor irradiation.Two fixed beam-line irradiation rooms (one with a horizontal beam line, one with a horizontal plus vertical beam line) and one proton gantry will be available for patient treatment [3] that is planned to commence in 2015.Irradiation will be performed using active scanning over a 20 Â 20 cm 2 area, where two scanning magnets scan the beam transversally over the target.The different depth layers are irradiated by varying the extraction energy from the synchrotron.Neglecting scattering, beam sizes between 4 and 10 mm FWHM at the focal point are available.
Other examples of recent work where such a model would be applicable are: (i) insertion of a thin scatterer in a rotating gantry to produce a symmetric phase space at the patient [4]; (ii) using a range shifter close to the patient to reduce upstream scattering and optimize the lateral penumbra of a scanned proton beam [5]; (iii) evaluating dose distribution inside the target after lateral beam profile modifications [6].

A. Coordinate system
The forward direction of the beam coincides with the s axis, as shown in Fig. 1.The horizontal and vertical coordinates of a particle are denoted by ''x'' and ''y'', while the angular direction (derivative of x or y with respect to s) is denoted x 0 and y 0 .ðz; z 0 Þ are generic transverse coordinates that refer to either the horizontal or the vertical plane.

B. Transverse phase space representation
With help of the Twiss functions [7], , , , the particle density of the beam in the transverse phase space in one plane, z ðz; z 0 Þ, of a Gaussian beam can be described as where E z;1 is the 1emittance.When the beam is subject to scattering, the emittance increases and the Twiss functions, which are typically calculated assuming vacuum up to the focal point (FP), are modified.It is therefore more convenient to represent the beam by its angular variance hz 02 i, covariance hzz 0 i, and spatial variance hz 2 i. h i denotes an average over the entire phase space, i.e., hfðz; Carrying out the integrals for hz 02 i, hzz 0 i, and hz 2 i gives the following relations for the Twiss functions: For convenience, the three ensemble averages have been bundled into a single vector, z, which represents the transverse phase space of the beam in one plane.Note that as the beam is focused, defocused, and scattered, z will vary along s.
As ð1 þ 2 Þ=, the emittance of the beam can be calculated from The rms beam divergence and beam width in one plane are given by For completeness, we also specify the beam FWHM, W z , and 80%-20% lateral penumbra, z , of a Gaussian beam in relation to z : i.e. the lateral penumbra is roughly half the beam FWHM.

C. Scattering power
The immediate effect of scattering is an increase in the beam divergence (''thin lens'') which manifests itself only further downstream in the form of a broader beam.The scattering power SðsÞ is defined as the rate of increase of the angular variance hz 02 i with longitudinal coordinate s: In this paper the scattering power formulation of Kanematsu [2] is used.While there are different formulations of the scattering power (see e.g.[8]), most have in common that they are simplifications of the more accurate scattering theory of Molie `re [9,10], or modifications of the semiempirical Highland formula [11] in order to account for the energy loss in thick scatterers.This paper does not attempt to evaluate the accuracy of any given scattering power formulation: the presented extension from scattering power to include beam optics and non-Gaussian beam profiles could be made with any scattering power formulation.
Kanematsu defines the scattering power as where E s is a constant (15 MeV), p is the beam momentum (MeV=c), c the beam velocity, Q the beam charge number, X 0 the radiation length of the material (g=cm 2 ), and the material density (g=cm 3 ).fðlÞ is a correction function of the parameter lðsÞ-the integral path length in radiation lengths: Note that with the integral definition of lðsÞ, SðsÞ is nonlocal, i.e., it depends not only on the beam energy and material at s, but also on the amount of material the beam has already traversed.Other formulations contain instead the initial beam energy as a parameter, which is also nonlocal.
The scattering angle, s , is the 1divergence that a point beam (zero emittance) will have after passing through some multislab geometry of total length L. s is given by the square root of the integrated scattering power over the entire geometry: Kanematsu's correction function is chosen such that the scattering angle over a homogeneous material coincides with the Highland formula.

D. Energy loss
For thick scatterers, the energy loss of the beam must be taken into account in the scattering integral due to the 1=ðpcÞ 2 term.For the simulations presented in this paper, this has been done using energy-range tables from SRIM [12,13] which gives the beam range R in a wide range of materials and compounds as a function of initial kinetic beam energy T 0 , i.e.RðT 0 Þ, as well as the inverse-kinetic energy vs beam range.
Figure 2 shows the typical energy-loss curve for a hadron beam with initial energy T 0 and range R 0 in a homogeneous material, with a rapid energy loss toward the end of the beam range.At a depth d, the residual beam range is R 0 À d, and the energy at that depth is given by the inverse R À1 ðR 0 À dÞ.In a multislab geometry, the beam energy in the downstream slabs is given by first calculating the energy at the upstream interface.
In this paper, spline interpolation of the SRIM energyrange tables has been used to determine the residual beam energy as a function of s.

E. Transfer line scattering model
In order to formulate a scattering model that is applicable to a transfer line, we start with the three conventional drift space scattering integrals [14] which gives hz 02 i, hzz 0 i, and hz 2 i for an initial beam z a that is scattered along a drift space, as indicated in Fig. 3: More compactly, in matrix notation, where the drift space transport matrix, M Drift s!b is defined as FIG. 2. Typical energy-loss curve in a homogeneous material.At a depth d, the residual beam range is R 0 À d and the beam energy is given by T d ¼ R À1 ðR 0 À dÞ.
FIG. 3. The initial beam phase space z a will evolve to z b after scattering through N slabs of thickness, density, and radiation length d i , i , and X 0;i .
The interpretation of the scattering integrals is straightforward: the initial phase space z a is simply transported from a to b by the transport matrix, while the integral transports small divergence increments SðsÞds from s into phase space increments at b.At b, the contributions from initial phase space and scattering are added.The scattering part is completely independent of the initial phase space.The coefficients of the transport matrix are easily resolved by considering the motion of a single particle in a vacuum drift space: From the single particle motion matrix, the evolution of z from a to b in vacuum can be calculated from which gives exactly the coefficients of M Drift a!b .To generalize to a transfer line with uncoupled motion, we modify the elements of the transport matrix correspondingly.The Twiss functions define the single particle motion from s ¼ s 1 to s ¼ s 2 as (subscript 1 and 2 denotes value at s 1 and s 2 ): where and Á is the phase advance from s 1 to s 2 : With these values, the general transport matrix M s 1 !s 2 becomes and the general scattering equation, For beam profile considerations, the beam size is naturally the most interesting quantity.Writing it out explicitly gives where (33) 0 is the 1width of the unscattered beam, projected onto the z axis, at s ¼ b, while s is denoted the scattering term.

F. Non-Gaussian beam profile
If the initial beam distribution is Gaussian, the width of the resulting beam further downstream is given by a quadratic addition of the unscattered beam width and the scattering term [Eq.( 31)].The quadratic addition is equivalent to a convolution between two Gaussian profiles of 1widths 0 and s .In the case of the MedAustron synchrotron, with a slow, third-order resonance extraction in the horizontal plane, the horizontal beam profile is more trapezoidal, with sharp edges, rather than Gaussian, as shown in Fig. 4. In this case, the quadratic addition of two rms widths should be replaced by an actual convolution between the unscattered beam profile density, 0 ðzÞ, at the point of interest (e.g. the focal point) and a Gaussian function with rms width s , G s ðzÞ: The lateral penumbra and FWHM of the scattered beam profile is given through inspection of ðzÞ.

III. BEAM-LINE DESIGN OPTIONS A. Scanning nozzles
The MedAustron nozzles will contain monitors for online verification of the position and profile of the delivered beam and therefore the dose delivered to the patient.It also contains passive elements, such as ridge filters for longitudinal widening of the sharp Bragg peak (not considered in this paper) and a range shifter (RS).The range shifter is a rectangular slab inserted into the beam path in order to reduce the energy of the beam below the minimal extraction energy of the synchrotron, thereby shifting the Bragg peak closer to the skin.Typically, the range shifter is made of a low-Z material, like polymethyl methacrylate (PMMA), which features a favorable (i.e.low) ratio between transverse scattering and energy loss [15].Even so, the range shifter will significantly scatter the beam and it is therefore important to minimize the drift space from range shifter to target in order to maintain an acceptable beam size and penumbra.
Three of the considered nozzle versions are shown schematically in Fig. 5.The nozzle is movable with a minimum distance between the focal point and the last nozzle element of 10 cm (reducing risk of collision for patient safety), and a maximum distance of 65 cm, which is considered sufficient space for the patient and imaging equipment.
Nozzle 1 includes monitors which are movable in the longitudinal direction to reduce the drift space between monitors and focal point.
Nozzle 2 is identical to Nozzle 1, but with a helium bellow inserted between the vacuum window and the monitors, thus reducing the amount of air in the beam path when the nozzle is extended.
Nozzle 3 has fixed monitors, but a range shifter movable in the beam direction at the end of the nozzle (the range shifter is remotely removable when not used).
Note that the design of the MedAustron nozzles is currently ongoing, and that the examples shown in Fig. 5 are merely used to demonstrate useful applications of the implemented scattering model.
A 100 m exit window of PMMA prevents the patient from reaching into the nozzle.Additionally, it protects the interior equipment of the nozzle from dirt, dust, or liquid (water from, e.g., a leaking water phantom could hit the proton gantry nozzle when it is positioned directly under the patient table).
All relevant properties of the materials considered are summarized in Tables I and II.

B. Proton gantry
The proton gantry is a 15.58 m long (beam path) rotating construction that allows for irradiation of the patient from any direction (see Fig. 6).To separate the rotating vacuum pipe of the proton gantry from the fixed, upstream, vacuum pipe, double vacuum windows (2 Â 50 m kapton) were originally foreseen at the coupling point (CP), as this is the simplest solution.Note that as the cross section of the nozzle vacuum pipe is larger than the vacuum pipe at the gantry CP, the vacuum windows at the CP can be made thinner.Another possibility would be to use a windowless rotating joint to eliminate any beam growth due to scattering at the CP.The 90 bending dipole at the end of the gantry contains no dedicated vacuum chamber.Instead, the magnet yoke serves as vacuum chamber in order to increase the aperture, thereby maximizing the scanning area at the patient.In the event of a vacuum leakage, however, reparation of the 70 ton dipole could be time consuming.By filling the entire dipole with helium (which has a lower scattering power than air) under atmospheric pressure, it would be robust against minor leakages.The windows containing the helium are assumed to be kapton coated on both sides with copper and aluminum (see Table I), as in the proton therapy facility M. D. Anderson [17].
The beam growth due to scattering in gantry windows and helium in the dipole will be evaluated in Sec.IV, using the foreseen gantry optics shown in Fig. 7 and the four different slab geometries along the gantry shown in Fig. 6: 1. the only scattering element taken into account is the vacuum window at the nozzle; 2. double vacuum windows at the CP, plus the nozzle vacuum window; 3. no CP vacuum window, but helium in the dipole (since the helium extends all the way to the nozzle, the last vacuum window is moved in front of the dipole); 4. CP vacuum windows as well as helium in the dipole.

IV. RESULTS AND DISCUSSIONS A. Proton gantry
To evaluate the different gantry options, we consider only the scattering term, s (multiplied by 2.35 for conversion to FWHM) along the gantry.Application of Eq. ( 30) to the four different slab geometries gives the transverse FWHM of an incoming point beam if z ¼ 0. As a ''worst case,'' we consider a low-energy proton beam (60 MeV).
Important observations from Fig. 8 are: 1.With only the nozzle vacuum window as scattering element (geometry 1), the beam FWHM would be at about 5 mm in x and y at the focal point for an incoming point beam.2. Adding vacuum windows at the coupling point (geometry 2) blows up the horizontal FWHM to 13 mm, while the vertical beam size is unaffected.This can be understood from the different optics in the horizontal and vertical plane.As there is a phase advance of 2 rad from coupling point to focal point in the vertical plane, scattering at the coupling point will have a negligible effect on the focal point beam size [see Eq. ( 33)]. 3. Removing the CP vacuum window but filling the dipole with helium (geometry 3) will increase the  horizontal FWHM to about 15 mm and the vertical to 18 mm.4. With both CP vacuum windows and helium in the dipole, horizontal and vertical beam widths will be close to 20 mm.To summarize, CP vacuum windows alone would blow up the horizontal beam size, while a dipole filled with helium would significantly increase the beam size in both planes.In all cases except geometry 1, the scattering term is completely dominant compared to the pristine beam profile (4 mm FHWM).Thus, in order to produce small and symmetric beam sizes in both planes, the vacuum should be kept all the way to the nozzle.For evaluation of the nozzle design, geometry 1, with a single vacuum window at the nozzle, is therefore used.

B. Nozzle
Having fixed that the vacuum is kept all the way up to the nozzle, only scattering in the nozzle needs to be taken into account.Using the implemented scattering model, different nozzle design options can quickly be evaluated and compared: a few important results are presented below.

Beam growth through nozzle
Figure 9 shows the beam FWHM and 1divergence along Nozzle 1 and Nozzle 2 (Fig. 5) for a 60 MeV point beam.In Nozzle 1, the divergence growth in the air gap between monitors and vacuum window is around 2 mrad, while, as expected, the divergence growth in the helium chamber of Nozzle 2 is negligible.This results in a 2.5 mm smaller beam at the focal point.

Penumbra at focal point, in air
The lateral penumbra in air at the focal point is of interest as this defines the quality of the beam as it enters the patient.Shown in Fig. 10  FIG. 8. Horizontal (upper half) and vertical (lower half) FWHM along gantry for geometries 1-4 (Fig. 6).The incoming beam is a 60 MeV proton point beam.
total beam path in air is constant: at the lowest extraction energy (shortest beam range), the lateral penumbra is reduced by about 1.5 mm by moving the monitors from 65 to 10 cm away from the focal point.With a helium chamber, the gain is twice as large at low energies (about 3 mm difference), and penumbrae of less than 4 mm could be achieved with the lowest proton energy.
The figure also demonstrates that scattering is mainly a low-energy issue.Already at about 8 cm residual beam range, the ''worst case'' penumbra (Nozzle 1, 45 cm air gap) is smaller than the ''best case'' penumbra (Nozzle 2, 10 cm air gap) at lowest extraction energy.

Penumbra at focal point, at Bragg peak
Perhaps of more clinical interest than the focal point penumbra in air is the penumbra at the Bragg peak, which is always located in the target.To evaluate this, a water slab has been added just upstream of the focal point, as shown schematically in Fig. 11.The amount of water upstream of the focal point is adjusted to the energy of the incoming beam such that the beam stops exactly at the focal point.The (horizontal) penumbra at the Bragg peak has then been evaluated over the available energy range for Nozzle 1 and Nozzle 2. Results are shown for protons and carbon ions in Fig. 12, as a function of Bragg peak depth in water.
At low energies, the Bragg peak penumbra is similar to the penumbra in air.As the extraction energy increases, scattering through the nozzle and air gap is reduced, and the Bragg peak penumbra decreases.However, at higher energies, scattering in the target is dominant and, from 5-8 cm depth (depending on nozzle and air gap), the beam profile widens inevitably again.Any efforts to optimize the nozzle geometry would thus only have an impact on superficial tumors: the penumbra in deep-seated targets is the largest and independent of nozzle geometry.
For comparison, Fig. 12 also shows carbon ions penumbra, which are much less affected by scattering.

Range shifter
The penetration depth of the lowest energy proton or carbon ion beam in water is approximately 3 or 3.6 cm, respectively.In order to position the Bragg peak even closer to the surface, a range shifter is required.For a given penetration depth in the patient, there is an optimum range shifter thickness that minimizes the lateral penumbra of the beam at the patient.The optimum thickness depends on how close to the patient the range shifter can be positioned, as illustrated in Fig. 13.With a thin range shifter, a low extraction energy is required in order to position the Bragg peak correctly.A thicker range shifter requires a higher extraction energy.This reduces the scattering power through the nozzle [Eq.(10)] and the lateral penumbra of the beam at the range shifter will be smaller.On the other hand, the divergence of the beam leaving the thicker range shifter will be higher, since it has passed through more matter compared to the case with a thinner range shifter.The beam will therefore grow faster and a thick range shifter is therefore beneficial only if the air gap can be made small.
Figure 14 shows a contour plot of the beam FWHM at the focal point (in air) as a function of range shifter thickness and residual beam range for four different air gaps (0, 5, 10, and 20 cm), using Nozzle 3.With zero air gap, a thicker range shifter is beneficial: for a given Bragg peak R res , the FWHM at the focal point decreases with increasing range shifter thickness.However, at larger air gaps (20 cm), the trend is the opposite: the lower scattering in the nozzle is not enough to compensate for the high beam divergence after the range shifter.At about 10 cm air gap, the exact thickness of the range shifter would not have a great importance.
Ideally, one would like the thickness of the range shifter to be optimal for the achievable air gap in any given situation, according to Fig. 14: thin range shifter for large air gaps; thick range shifter for smaller air gaps.For practical reasons, however, the nozzles will only be equipped with a single range shifter of fixed thickness.The cross section of the range shifter is slightly larger than the scanning region (about 22 Â 22 cm 2 ).With these dimensions, it is not uncommon that the patient's shoulder will limit how close to the patient one can position an external range shifter during irradiation of head and neck tumors.This constraint gives an estimated air gap in the range of 10-20 cm (with large individual variations), which discourages to use a ''thick'' range shifter.
To conclude, with Nozzle 3, and foreseen average air gaps around 10-20 cm, a PMMA range shifter with a water equivalent thickness of about 3.6 cm is preferable.A thicker range shifter would only be beneficial-in terms of beam FWHM and lateral penumbra-if the air gap can be reduced.

V. DISCUSSION OF MODEL
Several useful applications of the transfer line scattering model have been demonstrated.The model should however not be stretched beyond its applicability.Effects not taken into account are here addressed.
First, the single-particle motion matrix (which defines the transport matrix M a!b ) is calculated via the Twiss functions, which are given by the field strength in the magnetic elements and the beam energy.If a thick scatterer is inserted along the transfer line, the energy loss will be significant and the magnetic rigidity of the beam will no longer match the magnetic field of the downstream magnets.This affects the single-particle motion and distorts M a!b .By matching the strength of the downstream magnets to the magnetic rigidity of the beam, the transport matrix can be recovered and the Twiss functions maintained.
Second, range straggling in scattering elements will increase the energy spread of the beam.In a dispersion region (such as the 90 gantry dipole), this would have an effect on the resulting beam profile at the focal point, coupling energy to position.However, in the example with the helium filled gantry dipole considered in this paper, the scattering effect alone is enough to dismiss filling it with helium.
Other effects of beam line material on the patient dose should also be considered.Apart from multiple Coulomb scattering increasing the penumbra, an unwanted dose is also caused by neutron production (protons and carbon ions) and fragmentation (carbon ions) in the beam-line material.With active scanning techniques, most secondary neutrons are produced in the patient tissue [15] and therefore unavoidable.The integral neutron dose to tissue surrounding the tumor is in any case low compared to the dose sparing that can be achieved with actively scanned proton beams compared to photons [18].Nuclear fragmentation reduces the number of primary ions and adds a dose beyond the Bragg peak.

VI. CONCLUSIONS
A general model for scattering of therapeutic hadron beams with non-Gaussian transverse beam profiles along a transfer line has been presented.While it has been implemented using one particular formulation of the scattering power, any scattering power formulation could be used.
The scattering model has been used to evaluate various different design options of the MedAustron proton gantry and nozzles with respect to the beam profile at the patient.The model has proven to be an important tool to quickly identify elements that are critical to the beam quality at the patient and how much a modification of, e.g., the nozzle design will actually improve the beam quality at the patient.FIG.14. FWHM (vertical) at focal point, in air, of a proton beam for four different air gaps: 0, 5, 10, and 20 cm, Nozzle 3. The horizontal axis is the thickness of the PMMA range shifter, the vertical axis the residual beam range in water as the beam passes through the focal point.

FIG. 6 .
FIG.6.Gantry and four different slab geometries along gantry which were considered.
FIG. 9. FWHM and divergence of a 60 MeV proton point beam impending on Nozzle 1 (top) and Nozzle 2 (bottom), 10 cm air gap.The focal point is indicated with a dashed vertical line.The background colors indicate material.

FIG. 11 .
FIG. 11.Insertion of water slab at the focal point.The thickness of the water column is matched to the energy of the beam such that it always stops at the focal point.
FIG.13.Thin range shifter: increased scattering and beam growth through nozzle.Thick range shifter: less scattering in nozzle, but larger beam divergence after range shifter.

TABLE II .
Nominal density, radiation length, and their ratio for materials used in this paper.