{'fraction_non_alphanumeric': 0.03542006156074597, 'fraction_numerical': 0.015322288611261996, 'mean_word_length': 4.264506136065769, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the scattering power of radiotherapy protons

Scattering power (T = d/dx of mean squared multiple Coulomb scattering (MCS) angle), as used in proton transport theory, is properly viewed as a differential description of the Gaussian approximation to MCS theories such as Moliere's. That is, we seek a function T which, when integrated over a finite slab, will recover the Moliere/Fano/Hanson angle for that slab. To be accurate, T must include a single scattering correction, which means mathematically it must be nonlocal, depending on how much MCS has taken place as well as the energy and scattering material at the POI. We review five formulas for T and introduce a sixth, testing each against the Moliere/Fano/Hanson prediction as well as experimental data. We discuss how sensitive some practical problems are to the choice of T. That choice is probably most important for general Monte Carlo codes, which are expected to address a wide variety of problems.


Introduction
Consider an initially monoenergetic, monodirectional proton beam slowing down in a finite slab of some homogeneous material. At any depth x we can write down the stopping power S ≡ −dE/dx, the rate of decrease of energy with depth. There exists a well-known theory [1] for computing S, which depends only on the mean proton energy at x and atomic properties of the stopping medium. Integrating S over x we obtain, accurately, the total energy loss in a finite slab.
At any x we can also calculate a scattering power T ≡ d < θ 2 > /dx (1) the rate of increase with x of the mean squared projected multiple Coulomb scattering (MCS) angle. 1 There are a number of formulas for T in the literature. According to some of these T , too, depends only on local variables at x (mean proton energy and atomic properties of the stopping medium). However if we integrate these over x we do not, in general, obtain the correct rms MCS angle. The analogy with S is flawed. This might seem irrelevant since there exist accurate and well tested theories, notably that of Molière [2,3,4], which directly give the total MCS angle, at arbitrary incident energy, for homogeneous slabs of any element, compound or mixture from very thin to near stopping thickness. Moreover, the derivation [2,3] of Molière theory does not rely on the concept 'stopping power'. It is never mentioned.
Seen in that light, the aim of this paper is to discuss differential approximations to Molière theory, functions T which, if integrated over x for any slab, will recover the Molière MCS angle more or less accurately.
Why do we want to do that? Because practical problems in proton transport go far beyond merely finding the MCS angle in a finite slab. We might for instance wish to compute transverse spreading of a pencil beam in a single slab, or in a stack of slabs, or the equivalent source point 2 in a stack, or the transverse penumbra of a proton therapy beam formed by such a stack. Such computations are done either by Monte Carlo methods or by deterministic (numerical) methods using the Gaussian approximation to proton transport, generalized Fermi-Eyges theory. 3 Either way, we need to know, at every depth, the rate of change of the mean squared MCS angle. That need is exemplified by the formulas for the Fermi-Eyges moments, but it is more fundamental than that. Any transport calculation ultimately needs a recipe for the local rate of change of whatever variables are involved. It should go without saying that, when we do integrals like (2) numerically, approximating them by sums as in the midpoint rule or Simpson's rule, the answer should not depend on the step size ∆x over some reasonable range. The same applies if ∆x happens to be the step size in a Monte Carlo calculation. We mention this only because, in the case of T , this issue has lead to some confusion. For instance, Li and Rogers [8] in an electron paper 4 discuss the 'slab thickness dependence' of T at length, while a proton paper by Russell et al. [9] uses, without much discussion, Molière theory with the step size (there conflated with slab thickness) fixed at an arbitrary value. 1 In the early literature T is the rate of increase of the mean squared space angle, which is greater by a factor 2. In transport calculations projected quantities are usually more convenient. 2 We use the term loosely here. Various 'equivalent sources' are defined in the literature cited below. 3 We assume some familiarity with generalized Fermi-Eyges theory. For reviews see, for instance, ICRU Report 35 [5] or more recent papers by Hollmark et al. [6] and Kanematsu [7]. 4 For historical reasons Fermi-Eyges theory and the concept of T were first developed (in the radiotherapy context) for electrons, but the conceptual issues are the same for protons. Indeed, protons are much better Fermi-Eyges particles. The small angle approximation comes naturally and the correlation of mean energy with depth holds to far greater depth because range straggling and the MCS detour factor are much smaller.
We shall see that there is, in fact, no problem provided we keep the concepts 'step size', 'slab thickness' and 'depth of the point of interest (POI)' carefully separated. All numerical results given later hold over a wide range of step size. We often use ∆x (cm) equivalent to 0.1 g/cm 2 but a much smaller or larger step gives essentially the same answer.
The present paper deals with T and its implications in mixed slab geometry, with the slabs assumed to be infinite transversely. In most practical problems, whether beam line design or dose distribution in the patient, transverse heterogeneity or beam limiting devices eventually come into play. We must go beyond Fermi-Eyges theory, which is then relegated to propagating the individual pencil beams which make up the final (non Gaussian) distributions. The present paper ignores all that. We treat only the essentially one dimensional mixed slab problem, in order to understand that fully before introducing complications. Nor do we consider non electromagnetic effects like the halo from non-elastic nuclear reaction secondaries [10,11]. This is strictly about the Gaussian approximation to MCS in mixed slab geometry.
This paper overlaps a recent one by Kanematsu [7], whose notation we use except as noted. That work discusses the entire transport problem in the framework of Fermi-Eyges theory, treats heavy ions as well as protons, introduces a new scattering power T dH and presents many useful analytic approximations. The present paper focuses entirely on scattering power, protons, and single or mixed slabs, largely glossing over computational issues. It is an in-depth look at a small part of the ground covered by Kanematsu. His work is somewhat biased towards dose reconstruction in the patient, ours towards beam line design, but the fundamental issues are the same and there is no significant disagreement.
For completeness, we repeat a number of formulas from the literature, some of which is quite old. There is, however, considerable new content: simplifying the formula we call T IC , deriving a new non-local formula T dM , and proposing, for a case of practical importance, an experimental test to discriminate between various formulas for T . In short, we combine original material with a critical review which, we hope, will be useful to the student and of some interest to the expert.

Mixed-Slab Notation
Take, for concreteness, an example from beam line design: a monoenergetic proton beam of known emittance 5 enters a Pb/Lexan/air stack ( Figure 1). That is precisely what happens at the upstream range modulator/first scatterer of a modern passive beam spreading system [12]. 6 (In the patient, the variation of material properties is usually less extreme, but still present.) We wish to transport the beam through the stack to an arbitrary measuring plane MP at a depth x (cm) measured from the entrance of the first slab. In particular let us assume we wish to find y rms of the Gaussian which approximately describes the transverse fluence distribution on an MP located, say, at the second scatterer, because the precise match of y rms to the second scatterer design [13] will determine the flatness of the dose distribution at the patient.
M i is the material of which the i th slab is composed and stands for a set of material properties e.g. density ρ i (g/cm 3 ), radiation length X 0i (cm) and (to be defined) scattering length X Si (cm). There are also atomic weights A, atomic numbers Z and fractions by weight w of the constituents of M i , all implicit stepwise constant functions of x. x i refers to the entrance of the i th slab and x 1 ≡ 0.
E i (MeV) is the proton kinetic energy entering the i th slab, p i v i (MeV) the corresponding product of momentum and speed, and R i (cm) the mean proton range corresponding to E i stopping in M i , that is to say the residual range in M i at the entrance to the i th slab. E, pv and R (the latter always computed for the current material) are corresponding quantities at the depth of interest x. Given a specified stack and any value of x within it, all these quantities can be computed using range-energy tables and kinematic relations. 7 This example is a rare practical case where Fermi-Eyges theory alone gives a useful answer, because there is in fact no transverse heterogeneity and the final fluence distribution is in fact very nearly Gaussian. In more complicated beam line design problems, with collimators as well as scatterers whose thickness varies with radius, Fermi-Eyges theory can only be a building block, as noted earlier.

Limits of Molière Theory
Radiotherapy protons have kinetic energies roughly in the range 3 ≤ E ≤ 300 MeV corresponding to 0.015 ≤ R 1 ≤ 51.5 cm in water. The interesting range of normalized depth in a single slab is roughly 0.001 ≤ x/R 1 ≤ 0.97. The lower limit might for instance apply to a vacuum window in a beam scanning system. The upper, near stopping depth, is where Molière theory fails [4] as straggling blurs the correlation between depth and energy. It is reached and exceeded when protons stop in the patient.
It is well known that Molière theory applies only to slabs of some minimum thickness. There must be enough atomic encounters to establish the regime of multiple (as distinct from single or plural) scattering. 8 Molière [2] gives B = 4.5 as the lowest allowed value of his slab thickness parameter. 9 Table 1 lists target parameters at which that limit is reached for some common materials for 158.6 MeV incident protons. It is clear that the lower Molière limit need not concern us in practical proton radiotherapy calculations.

Kinematics
The kinematic expression 1/pv appears in all multiple scattering formulas. 10 If we define the reduced kinetic energy of any particle as where mc 2 is the particle's rest energy, we find In the early literature x is frequently expressed in units of the radiation length X 0 (therefore dimensionless) and instead of T , the mass stopping power T /ρ (radian 2 /(g/cm 2 )) is used. This is unhelpful for mixed slabs. 8 One can speak of the (extremely small) average energy loss in an atomic monolayer, or of the (extremely small) single scattering probability, but it makes no sense at all to speak of multiple scattering since there is at most one collision. That is the fundamental difference between stopping and multiple scattering theories. 9 B ∝ logarithm of normalized slab thickness. The constant of proportionality varies with the material [4]. 10 Multiple scattering is derived from single scattering as outlined in Section 3.1. 1/v enters the derivation of the single scattering probability (14) because the impulse delivered in a single collision is proportional to the interaction time or inverse speed. 1/p comes from the fact that the angle of deflection equals ∆p/p. For 3 ≤ E ≤ 300 MeV protons the coefficient of E varies from 2.00 to 1.76. Other useful relations in the same vein are β 2 = τ + 2 (τ + 1) 2 τ (5) and (pc) 2 = (τ + 2) mc 2 E (6) These formulas avoid differences between large quantities that can arise in relativistic calculations, and their relativistic and non-relativistic limits are obvious by inspection.

The Øverås Approximation
Øverås [14] found a simple empirical relation between (pv) 2 and normalized residual range. For a single slab where k depends on the material. If the material is characterized by its radiation length, Schneider et al. [15] found the empirical expression k is fairly small compared to 1, and often the 'weak Øverås' approximation is adequate. We shall use it later as a guideline. Figure 2 tests (9) over a wide range of materials, normalized slab thicknesses and outgoing kinetic energies.

Molière/Fano/Hanson Procedure: θ Hanson , T Hanson
For a single slab, given the outgoing energy, material and thickness, we can find the rms angle in the Gaussian approximation by using the Molière/Fano/Hanson procedure [4] with We use this angle as our standard of comparison for all rms angles obtained by integrating scattering powers over single slabs. It is Kanematsu's θ MH [7]. We have given it the longer name to emphasize that it is not derived from a scattering power. Indeed, we now derive a scattering power from it by numerical differentiation 11 Though we call this T Hanson it is not a scattering power in the usual sense. Evaluating it requires a lengthy procedure rather than a simple formula. However it is the 'correct' T for single slabs and we will use it later to derive an improved scattering power. 12

The Single Scattering Correction
Rather than thinking of a proton entering a given slab at some energy, let us focus on the proton at a depth x or equivalently, a proton leaving a single slab of thickness x. Assume the material is Be and the energy at x is 20 MeV. A 'local' formula for T uses only those two facts. But there is another parameter: the amount of material overlying x. Our 20 Mev might for instance represent an incident 23.7 MeV proton passing through an overlying 1 mm of Be (x/R 1 = 0.264), or a 102 MeV proton passing through an overlying 5 cm (x/R 1 = 0.95). Obviously the mean squared MCS angle itself depends not only on E(x) and M (x) but on the quantity of 'MCS buildup' material x/R 1 traversed to get to x. Does T , the rate of change of mean squared angle with x, also depend on x/R 1 ? It does. Figure 3 shows mass scattering power T /ρ, computed three different ways for three materials at 20 MeV, as a function of normalized overlying material x/R 1 . 13 T Hanson (11) is the 'correct' scattering power. It is obviously nonlocal: it depends on x/R 1 . T FR and T IC (discussed later) are local: they do not depend on x/R 1 . Of the two, T IC has the more accurate material dependence suggesting that it, with some logarithmic function of x/R 1 for the single scattering correction, might lead to an accurate approximation to T Hanson .

Highland's Formula: θ Highland
Highland [16], in order to simplify MCS calculations for the high-energy physicist, parameterized Molière/Bethe/Hanson theory [4] and obtained the elegantly simple formula 14 in which the single scattering correction, in parentheses, is evident. He assumed a slab that is finite but sufficiently thin that pv does not decrease much (x ≪ R 1 ). Thus (12) is already an integral expression. It cannot be applied to a thicker slab by interpreting x as a step size ∆x, using (12) for each step, and adding in quadrature: the sum decreases indefinitely with the number of steps.
We circumvented that problem in [4] by arbitrarily removing the logarithmic term from the sum (or integral) and proposing a generalized Highland formula which reduces to (12) for thinner slabs and was found [4] to agree with measurement almost as well as Molière/Fano theory. It is Kanematsu's θ iH , which we have given a longer name to emphasize that it is not computed from a scattering power. Figure 4 compares θ Highland (13) with θ Hanson (10). The behavior for Be reflects the difference between the Molière/Bethe/Hanson and Molière/Fano/Hanson forms of the theory. Otherwise, θ Highland is better than the ±5 % advertised [16], and much easier to compute than the full theory.

Prior Formulas for T
At this writing, three local and two nonlocal formulas for T can be found in the literature. We review the first two at some length because Rossi's excellent book [17] is no longer easily available and because the second, T IC , is the basis for our improved T . 13 3.1 Fermi-Rossi: T FR , θ FR Following Rossi [17] the single scattering probability for a singly charged particle in small angle approximation is 15 r e is the classical electron radius, N is Avogadro's number, A, Z and ρ are the atomic weight, atomic number and density of the target material and m e c 2 is the electron rest energy.
In the derivation of (14), the target nucleus is modeled as an unscreened point charge. More realistically, the scattering law must depart from 1/χ 4 at small angles (distant collisions) of order due to screening of the nuclear charge by atomic electrons, and at large angles (close collisions) of order 16 due to the finite size of the nucleus. α is the fine structure constant. Rossi now assumes that the value of Θ 2 at x + dx equals its value at x plus the mean squared space angle of scattering in dx. This step is equivalent to assuming the process is exactly Gaussian and ignoring the single scattering correction. It leads to Rossi now defines Later, Brahme [18] used (T /ρ) and the term 'mass scattering power' for Θ 2 s . We remind the reader that this and other early uses of (T /ρ) refer to the rate of change of space angle unlike later treatments and the present paper.
To perform the integral in (18) analytically we must assume some simple behavior of Ξ(χ) below χ 1 and above χ 2 . Rossi does this two different ways. In the first, Rossi assumes Ξ is zero in both regions. After integrating, simplifying and introducing 17 he obtains X 0 (cm) is the radiation length of the material. For a sufficiently thin single slab pv does not change much, the integral We use Molière's notation χ for single scattering space angle and related parameters and, in order to follow Rossi more closely, briefly use Θ for the multiple scattering space angle to distinguish it from θ, the projected angle. 16 We have kept two changes from ICRU Report 35. The 1.13 comes from the Thomas-Fermi radius of the atom, ra = 0.885 a 0 Z −1/3 (a 0 = Bohr radius), where Rossi used 1 instead of 0.885. In Eq.(16) ICRU35 rounded Rossi's (1/0.49) to 2. 17 Here we switch back to projected angle θ.
is trivial and the well known Rossi formula. For thicker slabs, pv can be related to x ′ as long as we know the range-energy relation in the material, and the integral is performed numerically. 18 Figure 5 compares θ FR , so computed, to θ Hanson . θ FR is far too large for thin scatterers, by an amount which depends on material.

ICRU Report 35: T IC , θ IC
Rossi next cites a much better approximation, assuming that Ξ(χ) behaves as 1/(χ 2 + χ 2 1 ) 2 at small angles (leveling off at small χ rather than suddenly vanishing). Eq. (18) then gives This is in essence the scattering power given in ICRU Report 35 [5], except that the version given there takes advantage of a chance cancelation with m e c 2 and therefore, unlike (22), only applies to electrons. Nevertheless, we shall call it T IC to distinguish it from T FR . For protons, T IC can be simplified considerably. When χ 2 (16) comes out larger than 1 radian it should be truncated to 1 [17]. For radiotherapy protons that never happens. The worst case is 3 Mev protons in Be where χ 2 = 0.913 rad. We can therefore simplify χ 2 /χ 1 , canceling the p dependence. If we also ignore the rightmost term in χ 2 /χ 1 , which is always much less than 1, and introduce a scattering length X S defined by we find for radiotherapy protons, 3 to 300 MeV, identical in form to T FR (20) For compounds or mixtures, any scattering power obeys a Bragg rule. Atoms act independently, and the compound or mixture is equivalent to very thin sheets of each constituent in the correct proportion. 19 That picture leads to where w i is the fraction by weight of the i th constituent. X 0 obeys a similar formula. Table 2 compares X S with X 0 for a few materials. Figure 6 compares θ IC , obtained by integrating T IC , with θ Hanson . Material dependence is greatly improved over T FR , but the large error for thin scatterers remains.
[6] use T IC for protons and heavy ions. However, the formula they quote (their Eq. (24)) is only valid for electrons, and they introduce an effective charge factor Z P,eff without comment. 20 Their values of T /ρ for protons up to 200 MeV in water (their Table 2) are larger than ours by a factor 1.19, perhaps due to incorrect application of the Bragg rule. Our Figure 14 seems to confirm that discrepancy. A footnote to a more recent paper [11] by the same group corrects the headings of Table 2 but does not mention any numerical error. 18 See Appendix A for details. 19 On fundamental grounds one would expect the Bragg rule for scattering powers to be very much better than for stopping powers, where (through I) there is some sensitivity to molecular binding [1]. 20 If we use the Barkas [19,20] then Z P, eff = 1 for radiotherapy protons.  Table 2: Comparison of scattering length X S with radiation length X 0 for six materials.

Linear Displacement: T LD , θ LD
In a recent note [21] Kanematsu proposes a simple scattering power for protons and heavy ions in tissue-like matter. He uses water as a reference material. For protons, his formula reduces to In a single slab of tissue-like material M where W stands for water, S is stopping power (MeV/cm) and R W (cm) is the proton's residual range in water at depth x. In the interest of a more complete survey we will, in what follows, explore the behavior of T LD for non tissue-like materials, even though it may not work well for those. We call it the 'Linear Displacement' scattering power after Kanematsu's derivation. For a single slab it can be integrated analytically giving A plot of θ LD vs. x/R 1 (Figure 7) is reminiscent of θ FR except near end-of-range. 21 Indeed T LD is a variant of T FR with three changes: the weak Øverås approximation (9) is used for pv, E s is adjusted downward, and water is used as a reference material. To show this consider a single slab of water. Then . When 1.25 is reduced to 1.00 and (31) is generalized to other materials by introducing ratios to water, (26) follows. Kanematsu [21] freely admits that T LD is semi-empirical, but his derivation [7, 21] is complicated and it is somewhat unclear exactly where the downward adjustment of E s takes place. It is probably mostly present in his Eq. (6) [21] which flows from Highland's equation which in turn is an empirical fit to Molière/Bethe/Hanson theory [16].
There is precedent for adjusting E s downward to improve the performance of T FR over a limited range of normalized thickness. Soukop where M (x) is the current material. t is a discontinuous function of x. It is that normalized depth which would be obtained if the protons were degraded from E 1 to E(x) in the current material. They derive their function of t by fitting a large body of single slab experimental data for θ rms at 158.6 Mev [4] with a two parameter analytic function of t (their Eq. (10)). Differentiating, they find 22 with fitted constants k(X 0 ) (8) and  34) and (35) is the mass radiation length of the current material. Figure 8 shows considerable improvement over the local formulas for T (note change in scale). Oscillations and divergent behavior at the ends are characteristic of polynomial fits. Two general remarks. First, Schneider et al. fit experimental data rather than some form of multiple scattering theory. Since theory and measurement agree rather well [4] this should not have a large effect, but it biases their formula towards the data that happen to be available. Second, their major advance is the introduction of a nonlocal correction based on a generalized definition of normalized depth. The Øverås approximation, tightly woven into their formalism to obtain formulas in closed form, is somewhat of a distraction. Similar results could be obtained without it. 3.5 Differential Highland: T dH , θ dH Kanematsu [7] derives a nonlocal scattering power applicable to mixed slabs which by construction gives the same result as (12) for a single thin slab. In his case the nonlocality parameter is depth weighted by inverse radiation length: he generalizes x/X 0 to a dimensionless radiative path length and writes a new scattering power as T FR times a correction factor He now constructs f so that, for a single thin slab, the integral of T dH will equal θ Highland . That requires that the average of f equal the squared ratio of the Highland formula (12) to the Rossi formula (21) 23 1 22 We have corrected two typographic errors. 23 lg ≡ log 10 , ln ≡ log e .
Unlike T ØS the extension of T dH to mixed slabs is totally straightforward, requiring only the new path integral ℓ(x).
4 Improved Nonlocal Formula: T dM , θ dM Instead of T FR as a basis let us use T IC (24) which, for radiotherapy protons, is as simple and has better material dependence ( Figure 6). From Figure 3, the correction should be logarithmic in total material overlying the POI.
In the weak Øverås approximation (9) we found that 1 − (pv/p 1 v 1 ) 2 , which depends only on local energy E and incident energy E 1 , is a reasonably good proxy for normalized depth x/R 1 for all materials and energies of interest over three orders of magnitude ( Figure 2). Let us therefore compute and plot the ratio of T Hanson , the ideal scattering power, to T IC , our proposed basis, for 0.001 ≤ 1 − (pv/p 1 v 1 ) 2 ≤ 0.97 and E 1 ≤ 300 Mev, for several values of the energy E at the POI and several materials (our usual Be, Cu and Pb). For each point we compute the exact x/R 1 for that material, without relying on the Øverås approximation, then use (11). The result is shown in Figure 10.
An adequate fit to these data, also shown in Figure 10, is a linear polynomial in lg (1−(pv/p 1 v 1 ) 2 ) whose two coefficients are in turn linear in lg (pv). Defining a new 'differential Molière' scattering power and writing everything out we have Although the weak Øverås approximation suggested the form of (40) the final result is an independent fit. It does not depend on the accuracy of either (7) or (9), or on the scattering material. Figure 11 compares θ dM obtained by integrating T dM with θ Hanson (note the vertical scale). The reader will object that the coefficients of (40), here given to excessive precision, are arbitrary. We could have chosen different sets of materials or energies. That is perfectly true, but it is also true of T ØS and T dH . The former is a fit to a specific data set for specific materials at a specific energy. The latter stems from Highland's formula, also a fit albeit to theory. That said, we have tried other materials and energies without much change in the general appearance of Figure 11.
Since θ dM agrees well with θ Hanson it may be inferred [4] that it agrees well with experiment for many materials. We will show that directly only for one low-Z and one high-Z material, also using the opportunity for a head-to-head comparison of six scattering powers and the generalized Highland formula. Figure 12 shows the comparison for polystyrene and Figure 13 shows it for lead, both with data from [4]. All three nonlocal T 's are better than any local T but T dM is the best. It almost agrees with measurement within the experimental error. 24 To facilitate numerical checks we include short tables of θ XX (Table 4) and (T /ρ) XX (Table 5) for single slabs of various materials and normalized thicknesses. We use T 1 = 158.6 MeV to correspond to [4] particularly Table 1. θ Hanson is not given directly there but may be found using (10).
The generalization of T dM to mixed slabs is the easiest of all. It does not even require an additional path integral. The single scattering correction is a logarithmic function of the fractional decrease in pv, with no explicit material dependence, from the incident beam to the point of interest. It diverges, as any single scattering correction must, at pv/(p 1 v 1 ) = 1 (no overlying material).

Applications
What kinds of computation are significantly affected by the choice of T ? If we are only interested in the Gaussian MCS angle itself in a single slab, Figures 4 through 9 and 11 through 13 already answer that question. However, in that case we would not use transport theory or T at all but simply evaluate either θ Hanson (10) or θ Highland (13). Let us examine some less trivial cases.

Pencil Beam Spreading in a Single Slab
The archetype for dose reconstruction in the patient is pencil beam spreading in a homogeneous water slab. An early paper by Preston and Koehler [23] derived a universal formula for beam spreading and compared theory with measurements. Hollmark et al.
[6] refer to additional measurements. Figure 5a of [4] already suggests that beam spreading in water is insensitive to T , comparing two very different models of MCS: Preston and Koehler's which is a local model similar to T IC , and a beam spreading model based on the nonlocal generalized Highland formula. Insensitivity of beam spreading in water to T implies that almost any T will work reasonably well for dose reconstruction, and is therefore worth a closer look. Figure 14 shows spreading of a 127 MeV pencil beam in water with y rms = (A 2 ) 1/2 (Eq. 2) computed according to all six formulas for T and according to Table 2 of Hollmark et al.
[6] (T HO ). Except T F R and T HO , all agree with experiment and are barely distinguishable from each other. Figure 15 is an expanded version where we plot the difference between each calculation and Preston and Koehler's formula (Appendix B). The spread, of order 0.1 mm, would be negligible in any dose reconstruction problem.
To gain some insight into this insensitivity to T , Figure 16 shows the integrand of A 2 when the POI is at 0.97 R 1 , near stopping depth. The near linearity of curves for the local formulas (T FR , T IC and T LD ) is a consequence of the Øverås approximation. The nonlocal formulas do in fact give a lower result because of the single scattering correction, important for small x, but not by much. Except for T FR , which is too high everywhere, the areas under the curves are nearly the same. That breaks down when we consider the integrand to x = 0.1 R 1 (Figure 17) but by then the absolute effect is so small as to be negligible.
It is instructive to look at the same problem in terms of the evolution of the beam phase space ellipse [5]. Figure 18 studies the same case as the preceding figures, dividing the water slab up to x = 0.97 R 1 into five sub-slabs. We also show bounding boxes for the final slab. 25 The spread of vertical bounds shows that the different T 's do give different answers for the final rms angle, but the effect on spatial spreading (horizontal bounds) is almost nil. The same study for Pb, Figure 19, shows that this is a fortuitous property of tissue-like materials, presumably due to a particular combination of drift and scatter. Beam spreading in near-stopping high-Z slabs is not entirely academic. In computing collimator scatter, we are basically asking how many protons leave the bore of a Cerrobend or other aperture that, except for MCS, would have stopped.

MCS Angle in a Double Slab (Range Modulator)
Unlike beam spreading in water, the choice of T is important in at least one practical problem in proton transport. The upstream range modulator in a passive beam spreading system is a sequence of high-Z/low-Z sandwiches designed to produce a Gaussian fluence distribution of constant width y rms either at the patient or, in a double scattering system, at the second scatterer [12,13]. Also, each sandwich is designed to pull back the pristine depth-dose by some fixed water equivalent amount so as to produce the desired spread-out Bragg peak.
The design procedure [24] is fairly complicated and the details need not concern us. We simply wish to define a sequence of Pb/Lexan/air stacks that might occur in a practical beam line, transport the beam through each such stack, and see how much difference T makes. We could close the loop with an experiment to see which T is best. Table 6 lists the parameters of a typical design. θ Hanson, 3 is the angle obtained by combining the Molière/Fano/Hanson angles for Pb and Lexan in quadrature. Scattering in air is ignored. More important, the single scattering correction for Lexan is wrong because its MCS angle is calculated for a beam of energy E 2 entering the Lexan de novo, ignoring the fact that some MCS has already taken place. For this and other reasons we do not expect the design and transport calculations to agree exactly for any choice of T . Figure 20 shows Fermi-Eyges transport results for each Pb/Lexan/air stack, with an ideal incident beam. 26 We used Kanematsu's finite increment form ([7] Eqs. [19][20][21], equivalent to integrating by the midpoint rule. The finite depth of the midpoint of the first step sidesteps the divergence of nonlocal T 's at x = 0. We discretize each slab separately with a single minimum step size parameter ρ∆x = 0.1 g/cm 2 for all materials, yielding typically tens of steps in Pb, a few hundred in Lexan, and one step for air. Increasing ρ∆x even by a factor 20 changes y rms less than a percent. We have not yet done an experiment to see what y rms these or similar setups actually produce, but Figure 20 is encouraging. The fact that our standard design procedure in fact produces a flat dose at isocenter (matches the second scatterer design) suggests that y rms = 3.5 cm is probably correct to a few percent. If, lacking a direct measurement of y rms , we assume that to be the case, we find that T dM indeed comes closest to the right answer whereas some of the other T 's are off by an amount much larger than the 1-2% a careful experiment could measure.
It is amusing that T LD also gives a very good result, even in a non tissue-like problem. That is consistent with Figure 7. θ rms for a moderately thick Pb slab is slightly low, whereas for a thick plastic slab it would be slightly high. It augurs well for 'corrected Rossi' T 's in general [22] even though the excellent result here may be somewhat accidental (or, depending on the outcome of a measurement, possibly wrong).
One could argue that, because we have an adequate 'slabwise-Hanson' design procedure, we again don't really need either beam transport or T . And indeed we don't, if we are merely designing beam spreading systems. However, the Pb/Lexan/air stack could be part of a larger Monte Carlo computation [25]. That brings us to the last topic.

Monte Carlo Calculations
Monte Carlo calculations are the gold standard in radiotherapy. A condensed history Monte Carlo should embody a differential model of multiple scattering analogous to scattering power, even if it is not called that. In each (finite but small) step of known material, one needs to compute the increase in the width parameter of some distribution (Gaussian, Molière or other) from which a random deflection angle is then drawn. Just as in deterministic calculations, one would like the final result to be independent of step size over some reasonable range. Table 3 shows for a simple case that results from Geant 4, whose MCS model is based on a variant of Highland's formula [26], depend somewhat on step size. Monte Carlos based on a 'corrected Rossi' formula [22], similar to T LD , should not have this problem, but may not always give the right answer (Figure 7) CERN. 158.6 MeV protons scatter in 20.196 g/cm 2 Pb (θ rms,expt = 108 ± 1 mrad, θ Hanson = 111.1 mrad). STEPMAX governs the step size and θ rms is obtained by fitting the projected angular distribution with a Gaussian.
lack the single scattering correction. T dM , easily generalized to other charged particles, might be a good compromise between accuracy and step size independence.

Summary
Unlike stopping theory, accurate theories of multiple Coulomb scattering such as Molière's do not flow from a differential form. If that is needed, for deterministic or Monte Carlo transport calculations, it must be devised retroactively as an approximation to the more exact theory. To be accurate for thin scatterers it must include a single scattering correction. That implies nonlocality: in addition to energy and properties of the material at the POI, it must depend in some way on how far the multiple scattering process has advanced. That does not conflict with being numerically integrable in the usual sense that an approximating sum approaches a limit as the step size decreases. Nonlocality may be characterized in different ways. Schneider et al. use a generalized definition of the normalized depth of the POI, Kanematsu uses a radiative pathlength integral up to the POI, and we use the diminution of pv from its incident value to the POI. All three can be tested against Molière theory for uniform slabs, and all three generalize to mixed slabs, as they must to be useful.
We have reviewed three local and two nonlocal formulas for T , comparing θ rms from each one graphically with the 'correct' MCS angle θ Hanson . One of them, T IC , can be simplified for protons by introducing a new material parameter, the scattering length X S , which is similar in form to radiation length X 0 . T IC provides a basis for a new nonlocal scattering power T dM (Eqs. 23,25,39 and 40) which, for single slabs with 0.001 ≤ x/R 1 ≤ 0.97, reproduces θ Hanson as well as the 'correct' scattering power T Hanson to a few percent.
In practical problems, the choice of T is frequently not critical. We have shown, for instance, that all T 's described here except T FR work well for beam spreading in water and, presumably, water-like materials. In particular T LD , typical of 'corrected Rossi' formulas, works well. Beam spreading in Pb and, presumably, other high-Z materials, is more sensitive to T .
Turning to a mixed slab problem, we have shown that Pb/Lexan/air combinations typical of an upstream range modulator yield very different answers for different T 's, so that a direct experimental test should be easy. Again in this case, T LD (even though used outside its supposed range of validity) yields almost the same answer as T dM , probably because of a lucky cancelation.
In the end, perhaps the strongest case for an accurate scattering power such as T dM can be made for general purpose Monte Carlo codes which are supposed to do everything well, rather than beam design or patient dose computations where special purpose workarounds can be devised.

Acknowledgments
We are indebted to Drs. Kanematsu, Schneider and Hollmark for correspondence regarding their work and to Dr. L. Urban of CERN for generating the data for Table 3. We particularly thank the Harvard Physics Department and the Laboratory for Particle Physics and Cosmology (LPPC) for ongoing support.

A Computational Details
Examples, figures and tables were computed with Fortran programs which may be downloaded free: see \BGware.zip at http://physics.harvard.edu/˜gottschalk. Source code is in \BGware\source. All calculations are single precision.
Formulas for variants of Molière theory are given in [4], implemented by module THETA0.FOR. The only significant change is that we now use cubic spline interpolation of ln R(ln T ) 27 , rather than a polynomial fit, to interpolate range-energy tables (RANGE.FOR). 28 Subroutine TOUT computes the energy out of a stack of slabs, as well as the outgoing rms projected MCS angle using what we have called the slabwise-Hanson or slabwise-Highland procedure, not Fermi-Eyges transport.
Figures and tables specifically for this paper were computed with various branches of SPWR.FOR. ProjScatPower(x) computes the projected scattering power at x in a stack according to the formula selected by software switch spMode. WhatsHere(x...) returns material properties and such quantities as Kanematsu's ℓ and Schneider's t at x. The body of each table set in L A T E X was produced as a text file to avoid errors of transcription. Data matrices required for graphics were also imported from text files.
We computed isolated examples with our 'proton desk calculator' LOOKUP, a WinXP executable distributed with BGware. LOOKUP is a convenient driver for some of the subroutines mentioned, offering a choice of 'tasks'. A useful one in the present context is STACK, which computes various quantities as a proton beam proceeds down a stack of slabs.
Results depend somewhat on the choice of range-energy tables. Common tables differ by 1-2% for a given material, presumably because of different choices of I, the mean excitation energy. We used ICRU Report 49 [1] throughout (our table MIXED.RET in \BGware\data) despite some experimental evidence [27] that Janni's 1982 tables [28] are better for water. Reference [4] used Janni's 1966 tables [29].
Finally, some comments on integration. We need to evaluate A 0 (Eq. 2) at many values of x/R 1 for the numerous graphs of θ XX compared with θ Hanson . By far the most efficient way is to divide the slab by equal ratios so that the contributions of each step to the sum are nearly equal instead of very different. 29 If, in addition, Simpson's rule is used rather than the midpoint rule, extremely fast and accurate integration is achieved. The math (subroutine SimpRat in module SPWRSUBS.FOR) is slightly confusing. Suppose we wish to divide a slab of thickness x into n steps (∆x) 1 . . . (∆x) n such that the ratio of successive steps is a constant r. Let us assume some provisional value for n. The formula for the sum of a geometric progression 1 + r + r 2 + · · · + r n−1 = (r n − 1)/(r − 1) leads directly to formulas for the last and first steps In our code T is kinetic energy, E is total energy, R is in g/cm 2 and the longitudinal coordinate is usually z. 28 In general, we interpolate published tables for range-energy relations but compute MCS quantities directly from one of the applicable theories. We use few analytical approximations, but many lookup tables. 29 There is no advantage in dividing the higher Fermi-Eyges integrals A 1 and A 2 by ratios. and we simply divide by r on each iteration. That works for any integer n ≥ 1 and any r > 1. It is physically reasonable to let because the need to subdivide at all (variation of pv) occurs only when the residual range (R 1 − x) is significantly smaller than the range R 1 . By the same argument it is foolish to pick n at the outset, since that leads to unnecessary subdivision of thin slabs. Therefore we input some desired maximum value of r and let For R 1 − x = 0.03 R 1 and r max = 1.6 we find n = 8, the most steps needed in practice.
Unlike the midpoint rule, Simpson's rule uses T (0) where any nonlocal formula for T diverges because the single scattering correction diverges. That is easily fixed by substituting a very small positive value, a given fraction of the first step size, for 0. Final results are quite insensitive to that fraction, which we have adjusted to suppress discontinuities, as n takes on successive values, in the graphs of θ XX vs. x/R 1 .

B Preston and Koehler's Formula
Preston and Koehler [23] show that, to a good approximation, the rms radius σ = √ 2 y rms at any normalized depth t in any single slab is related to its maximum value σ 0 by where t ≡ x/R 1 is the normalized depth. σ 0 in water is and β is given by (5) with E corresponding to the kinetic energy at a depth R 1 /2. Their derivation of (41) involves the Øverås approximation which they apparently discovered independently. They use a quantity similar to T IC based on the multiple scattering formalism of Bethe and Ashkin [30]. Kanematsu [21] gives a more modern derivation of (41) based on T LD .

List of Tables
A test using Geant4 v9.1, courtesy L. Urban, CERN. 158.6 MeV protons scatter in 20.196