A new calculation formula of nuclear cross-section of therapeutic protons

We have previously developed for nuclear cross-sections of therapeutic protons a calculation model, which is founded on the collective model as well as a quantum mechanical many particle problem to derive the S matrix and transition probabilities. In this communication, we show that the resonances can be derived by shifted Gaussian functions, whereas the unspecific nuclear interaction compounds can be represented by an error function, which also provides the asymptotic behavior. The energy shifts can be interpreted in terms of necessary domains of energy to excite typical nuclear processes. Thus the necessary formulas referring to previous calculations of nuclear cross-sections will be represented. The mass number AN determines the strong interaction range, i.e. RStrong = 1.2×10-13·AN1/3cm. The threshold energy ETh of the energy barrier is determined by the condition Estrong = ECoulomb. A linear combination of Gaussians, which contain additional energy shifts, and an error function incorporate a possible representation of Fermi-Dirac statistics, which is applied here to nuclear excitations and reaction with release of secondary particles. The new calculation formula provides a better understanding of different types of resonances occurring in nuclear interactions with protons. The present study is mainly a continuation of published papers.


Introduction
The knowledge of the total nuclear cross-section Q tot of protons is an important impact with regard to sophisticated features of therapy planning, since Q tot provides essential information of the following aspects: Decrease of the fluence of primary protons Φpp and release of secondary particles and their transport (secondary protons, neutrons, clusters like H2 1 , H3 1 , He3 2 , He4 3 , heavy recoil nuclei, which usually undergo either a β + or βdecay with additional emission of a γ-quant). With regard to secondary protons we have to differ between two kinds, namely protons resulting from nuclear reactions with production of heavy recoils and those protons, which are, in reality, primary protons and have undergone elastic and inelastic scatter by strong interactions according to the Breit-Wigner formula. Elastic scatter by nuclear forces is only a deflection of the projectile protons with additional energy-momentum transfer to the whole target nucleus, whereas inelastic scatter is connected to excitations of nuclear vibrations, rotations and transitions to excited states without releasing other nuclear particles, i.e. some quantum number will be changed. The resonances due to Breit-Wigner formula represent the main part of the resonance domain Eres according to Figure 1. Nuclear reaction types cannot be regarded as simple resonances, they mainly occur for proton energies E > Eres.
In previous publications 1, 2, 3 , we have developed a calculation method based on a nonlinear and nonlocal Schrödinger equation with a Gaussian kernel and on an interacting many-body system containing strong interactions, spin-orbit coupling and Coulomb interaction with inclusion of various excited configurations. The results of these calculation methods, which provide excited nuclear states, virtual compounds and nuclear

Introduction
The knowledge of the total nuclear cross-section Q tot of protons is an important impact with regard to sophisticated features of therapy planning, since Q tot provides essential information of the following aspects: Decrease of the fluence of primary protons Φpp and release of secondary particles and their transport (secondary protons, neutrons, clusters like H2 1 , H3 1 , He3 2 , He4 3 , heavy recoil nuclei, which usually undergo either a β + or βdecay with additional emission of a γ-quant). With regard to secondary protons we have to differ between two kinds, namely protons resulting from nuclear reactions with production of heavy recoils and those protons, which are, in reality, primary protons and have undergone elastic and inelastic scatter by strong interactions according to the Breit-Wigner formula. Elastic scatter by nuclear forces is only a deflection of the projectile protons with additional energy-momentum transfer to the whole target nucleus, whereas inelastic scatter is connected to excitations of nuclear vibrations, rotations and transitions to excited states without releasing other nuclear particles, i.e. some quantum number will be changed. The resonances due to Breit-Wigner formula represent the main part of the resonance domain Eres according to Figure 1. Nuclear reaction types cannot be regarded as simple resonances, they mainly occur for proton energies E > Eres.
In previous publications 1, 2, 3 , we have developed a calculation method based on a nonlinear and nonlocal Schrödinger equation with a Gaussian kernel and on an interacting many-body system containing strong interactions, spin-orbit coupling and Coulomb interaction with inclusion of various excited configurations. The results of these calculation methods, which provide excited nuclear states, virtual compounds and nuclear

Introduction
The knowledge of the total nuclear cross-section Q tot of protons is an important impact with regard to sophisticated features of therapy planning, since Q tot provides essential information of the following aspects: Decrease of the fluence of primary protons Φpp and release of secondary particles and their transport (secondary protons, neutrons, clusters like H2 1 , H3 1 , He3 2 , He4 3 , heavy recoil nuclei, which usually undergo either a β + or βdecay with additional emission of a γ-quant). With regard to secondary protons we have to differ between two kinds, namely protons resulting from nuclear reactions with production of heavy recoils and those protons, which are, in reality, primary protons and have undergone elastic and inelastic scatter by strong interactions according to the Breit-Wigner formula. Elastic scatter by nuclear forces is only a deflection of the projectile protons with additional energy-momentum transfer to the whole target nucleus, whereas inelastic scatter is connected to excitations of nuclear vibrations, rotations and transitions to excited states without releasing other nuclear particles, i.e. some quantum number will be changed. The resonances due to Breit-Wigner formula represent the main part of the resonance domain Eres according to Figure 1. Nuclear reaction types cannot be regarded as simple resonances, they mainly occur for proton energies E > Eres.
In previous publications 1, 2, 3 , we have developed a calculation method based on a nonlinear and nonlocal Schrödinger equation with a Gaussian kernel and on an interacting many-body system containing strong interactions, spin-orbit coupling and Coulomb interaction with inclusion of various excited configurations. The results of these calculation methods, which provide excited nuclear states, virtual compounds and nuclear reactions via S-matrix, transition probabilities and finally, total nuclear cross-sections Q tot (E), can be translated to the collective nuclear excitation model. This model only uses Z and AN as parameters and suitable analytic functions to describe all properties of Q tot (E). The complete contents of this figure have been previously discussed. For protons, a threshold energy ETh exists to surmount the potential barrier of the oxygen nucleus. Thus for proton energies less than ETh nuclear reactions cannot occur. At E = Eres = 20.12 MeV Q tot (E) shows a maximum value, and thereafter, it decreases exponentially to reach the asymptotic behavior at about E = 100 MeV. According to an integration procedure 2, 6 the analytic version of Figure 1 (equation (1)) provides the decrease of primary protons. The following integration procedure has to be carried out: By that, we obtain the following formula for the fluence decrease of protons:  reactions via S-matrix, transition probabilities and finally, total nuclear cross-sections Q tot (E), can be translated to the collective nuclear excitation model. This model only uses Z and AN as parameters and suitable analytic functions to describe all properties of Q tot (E). The complete contents of this figure have been previously discussed. For protons, a threshold energy ETh exists to surmount the potential barrier of the oxygen nucleus. Thus for proton energies less than ETh nuclear reactions cannot occur. At E = Eres = 20.12 MeV Q tot (E) shows a maximum value, and thereafter, it decreases exponentially to reach the asymptotic behavior at about E = 100 MeV. According to an integration procedure 2, 6 the analytic version of Figure 1 (equation (1)) provides the decrease of primary protons.
The following integration procedure has to be carried out: By that, we obtain the following formula for the fluence decrease of protons:  reactions via S-matrix, transition probabilities and finally, total nuclear cross-sections Q tot (E), can be translated to the collective nuclear excitation model. This model only uses Z and AN as parameters and suitable analytic functions to describe all properties of Q tot (E). The complete contents of this figure have been previously discussed. For protons, a threshold energy ETh exists to surmount the potential barrier of the oxygen nucleus. Thus for proton energies less than ETh nuclear reactions cannot occur. At E = Eres = 20.12 MeV Q tot (E) shows a maximum value, and thereafter, it decreases exponentially to reach the asymptotic behavior at about E = 100 MeV. According to an integration procedure 2, 6 the analytic version of Figure 1 (equation (1)) provides the decrease of primary protons.
The following integration procedure has to be carried out: By that, we obtain the following formula for the fluence decrease of protons:   A further result is that in the environment of a nucleus the effective potential can be calculated by a linear combination of two Gaussians. This is not true for longer ranges of Coulomb forces. However, this is not a simple task, since the shielding of the nuclear repulsion by the shell electrons has to be accounted for. The nuclear potential according to Figure 3 and the related parameters can be best calculated by equation (3), which assumes the shape: The units of this formula are stated in terms of units of Rstrong (range of strong interaction) according to Figure 3.

Materials and Methods
In previous publications 2, 3 , we have performed the task to calculate all nuclear properties to determine finally with the help of the S-matrix and transition probability of all possible configuration states the total nuclear crosssection Q tot (E) in terms of Z and AN. The results are closely related to the Bethe-Weizsäcker formula for the nuclear binding energy. The old formula for the calculation of ETh and Q tot (E) according to a previous paper will be stated in section "Summary of previous investigations". In order to determine all these properties we have to know the threshold energy ETh, i.e. that energy necessary to surmount the nuclear potential mount according to Figure  3. There exists also the possibility to circumvent this threshold restriction by the quantum mechanical tunneling process with energy E < ETh, but this is only a rather small effect, which leads to a little roundness of Q tot (E) at E < ETh. This aspect will be considered in section "Calculation procedure of the energy levels of excited states and quantum mechanical tunneling effect for E < Eth".

Summary of previous investigations
According to the results presented in 2, 3 we need for the calculation of Q tot (E), at first, the threshold energy ETh as a function of Z and AN. This function can be obtained via the formulas (4 -6). The second step provides the determination of the total nuclear cross-section Q tot (E), which can be accounted for by preceding studies 2, 6 with the help of equations (7 -9). Thus the formulas (4 -6) are necessary to determine the threshold energy ETh by a balance equation of Coulomb repulsion and oscillator model for strong interactions.
In the domain of the resonance energy Eres we obtain:   In the domain of the transition to the asymptotic behavior of the total nuclear cross-section the following formula is applicable: In order to give a qualitative motivation of the new calculation method we consider again Figure 1. At first, we have to repeat that Q tot (E) is not only restricted to proper nuclear reactions of protons with release of secondary particles: Thus for proton energies E > ETh (the threshold energy ETh amounts to 7 MeV) we can verify a rapid increase of Q tot up to a resonance maximum Eres = 20.12 MeV. This behavior up to the environment of the resonance maximum can be described by a proper Gaussian distribution. What happens in this domain? Elastic scatter of proton at the nuclear potential is dominated by strong interaction and mediated mainly by mesons, if the quantum state of the nucleus is not changed. This implies that in order to satisfy energy-momentum relation only the whole nuclear adopts energy and momentum (kinetic energy), the impinging proton is slightly deflected. If the nucleus is also excited by vibrations, rotations or transitions to an excited configuration, then the whole process is inelastic and by emitting γ-quanta it is damped to finally return to the ground state. These effects are mainly described by the Breit-Wigner formula and its generalization. 6 -8 It has to be pointed out that the secondary proton under these conditions is still the primary proton, which is deflected by a slightly higher scatter angle compared to Molière multiple scatter theory. 9 Figure 1 presents the total nuclear cross-section of oxygen; we can verify that, besides the Gaussian distribution of the environment of the maximal value Emax, after a slower decrease of Q tot , the asymptotic behavior is reached (this is certainly valid for proton energies E < 300 MeV). The asymptotic branch can be represented by a suitable error function erf(E), which has to satisfy some further boundary conditions. This error function erf(E) results by an integration over Gaussian resonance distributions of the energy within finite boundaries. The transition from the domain Emax to the asymptotic domain can either be represented by a sum of exponential functions or by a Gaussian distribution with an additional energy shift. In every case, this consideration indicates that there exists an alternative representation of the total nuclear cross-section Q tot (E) besides the previously studied one according to equations (7 -9). Thus equation (10) provides the new formula for Q tot (E), which appears to lead to a better access to quantum mechanical resonance mechanisms expressed by harmonic oscillators. Equation (11) provides all terms necessary for the determination of boundary conditions of the whole problem, and equation (12) an alternative calculation procedure for ETh. Aboundary has the purpose to ensure that Q tot (E) assumes zero at E = ETh, since the Gaussians do not vanish at this position. A modification will be accounted for, when we shall include the tunneling effect.     Table 2 presents the corresponding parameter values to perform calculations of equation (10) with the help of equation (11). All necessary parameters of equations (10, 11) can be determined via Table 2 and the formula: The total nuclear cross-section requires also the threshold energy ETh, which can be calculated by the following formula (12), which is easier to handle than formulas (4 -6):

New calculation formulas for Q tot (E) and ETh
The parameters of formula (12)

Calculation procedure of the energy levels of excited states and quantum mechanical tunneling effect for E < ETh
The investigations referring to the tunneling effect are as old as the quantum mechanics itself, since this theory has been used to explain the α-decay of heavy nuclei by Gamow. In many textbooks of quantum mechanics the tunnel effect is treated in detail of a particle passing with energy E through a potential box V with E < V. In the meantime, some other important tunnel effects have been studied, e.g. Josephson junctions and, recently, the H bonds between the complementary DNA strands through a double minimum potential. 14 [Also, see reference list in Ulmer et al. 14 ]. After this digression, we intend to return to Figure 3, which refers to the nuclear potential of O8 16 , but all formulas developed here can be applied to other nuclei, which exhibit similar properties of their nuclear potential. In following, we need the deflection point ξ of this figure. If we reduce formula (3) to one Gaussian with V1 = 0, this point is easy to determine, since vanishing of the second derivation provides ξ 2 = σ0 2 /2. The deflection point of the nuclear potential according to formula (3) can only be obtained by an iteration procedure: In a first step, we restrict equation (13) to first-order terms; by that, we obtain a quadratic equation in terms of ξ0 2 : In the next step, we have to solve a quadratic equation, too, since we insert the solution with ξ0 into the Gaussians of equation (13) and determine ξ1 2 according this equation by putting ξ1 2 = x 2 , i.e. the following step is given by: In principle, this iterative procedure can be repeated, but stopping after step 2 is very sufficient. We now determine the energy levels of the potential type like that of formula (3). It offers to approximate the potential by a 3D harmonic oscillator. However, this is not sufficient, since the energy levels of the excited states are not equidistant. 2 Therefore we use a nonlinear and nonlocal Schrödinger equation, which incorporates nuclear interactions as a self-interacting field. With the help of equation 3 this generalized Schrödinger equations assumes the form: The coupling constant λ has to be chosen such that the dimension agrees on both sides, but it can be put λ = 1; the magnetic interaction, i.e. spin-orbit-coupling, can be added using principles previously given. 2,3 With the help of the generating functions of Hermite polynomials we are able to write equation (15) in the form: The equation above represents a highly inharmonic oscillator equation of a self-interacting field. Since the square of the wave-function is always positive definite, all terms with odd numbers of n1, n2, and n3 vanish due to the anti-symmetric properties of those Hermite polynomials. For rc ≤ ξ (domain with positive curvature), the whole equation is reduced to a harmonic oscillator with self-interaction; the higher-order terms are small perturbations: The solutions of this equation are those of a 3D harmonic oscillator; the classification of the states by SU3 and all previously developed statements with regard to the angular momentum are still valid. The only difference is that the energy levels are not equidistant; this property can easily be verified in one dimension. The usual ground state energy is ћω0/2. This energy level is lowered by the term ~Φ0,0,0, depending on the ground-state wavefunction. The energy difference between the ground and the first excited state amounts to ћω0; this is not true in the case above, since the energy levels of all excited states depend on the corresponding eigen-functions themselves (these are still the oscillator eigen-functions!). Next, we will include the terms of the next order, which are of the form ~ λ•( Φ 0,2,2, Φ2,2,0, Φ2,0,2): The additional term T represents tensor forces. The whole problem is still exact soluble. In further extensions of the nonlinear/nonlocal Schrödinger equation, we are able to account for spin, isotopic spin, and spin-orbit coupling.
The spin-orbit coupling, as an effect of an internal field with nonlocal self-interaction, is plausible, since the extended nucleonic particle has internal structure; consequently, we have to add Hso to the nonlinear term: With the help of the function set (21) all integrals can be evaluated analytically. By that, the task remains to determine the coefficients of the function set by an iterative procedure, since equation (23) yields a third order equation of all coefficients forming the wavefunction. However, this so-called self-consistent field method is very familiar in many-particle problems, and, therefore, a detailed description is superfluous. The performance of the above task provides, besides the ground state energy, excited states for E < ETh and virtually excited states, if E > ETh. However, we are interested in the role of excited states of a nucleus, when a proton can penetrate the potential wall for E < ETh and virtually excited states are only important with regard to nuclear 2 1 With respect to the Dyson expansion the terms according to equations (29, 29a) are very convenient, since iterations always yield terms of the same structure. Please note that the kinetic energy term of the proton assumes complex values, if the impinging energy is lower than the potential wall, i.e. the position probability within the positive, repulsive potential suffers damping. Some consequences will be given in the result section.

Results
The following Figures 4 -7 deal with the nuclei of carbon, oxygen, calcium and copper. They show that both calculation formulas provide equivalent results; we use the following abbreviations: FM (former method) and PM (present method). However, the chosen type of functions used in this communication offers new ways to analyze nuclear cross-sections due to the advantages of Gaussian functions and kernels.
The following Figures 8 -11 show with regard to the most important nuclei in radiotherapy with protons roundness in the environment of the threshold energy ETh and the previously assumed condition Q tot (E) for E < ETh does not hold due to the quantum mechanical tunneling effect. This roundness can be amplified by the energy spectrum of the impinging protons resulting from the beam-line and a jump of the fluence decrease of primary protons Φpp at E < ETh is prevented (see also Figure 2 with regard to the passage of protons through water). The passage of protons through other media, e.g. calcium, leads to a similar roundness. 3 By that, the quantum mechanical tunneling effect can be best studied by really mono-energetic protons with initial energy E < ETh.
Both Figures 8 and 9 also present a way to check ranges of protons in biologically significant media, since at the end of the particle track very specific nuclear interactions occur with the help of the quantum mechanical tunneling effect. The intermediary existence of the isotopes N7 13 and Sc21 41 provides a tool to study proton interactions for very low proton energies in that media. The quantum mechanical tunnel effect represents the necessary basis for these studies.

Discussion and Conclusion
Let us first consider the low energy domain, i.e. up to the resonance maximum Eres. Thus for E < ETh only due to the tunneling effect the projectile proton can enter the interior of the potential wall according to . This provides in the case of oxygen the isotope F9 17 as intermediary excited state (or N7 13 and Sc21 41 with respect to carbon or calcium). Since E < ETh the lifetime of the excited states depends on the transition probability to a lower energy state to finally produce a γ-quant and a neutron by an exchange interaction (Pauli principle) within the nuclei via a π + -meson to yield a neutron, which can easily escape the potential wall. The resulting processes are also described by equations (30 -32).

Discussion and Conclusion
Let us first consider the low energy domain, i.e. up to the resonance maximum Eres. Thus for E < ETh only due to the tunneling effect the projectile proton can enter the interior of the potential wall according to . This provides in the case of oxygen the isotope F9 17 as intermediary excited state (or N7 13 and Sc21 41 with respect to carbon or calcium). Since E < ETh the lifetime of the excited states depends on the transition probability to a lower energy state to finally produce a γ-quant and a neutron by an exchange interaction (Pauli principle) within the nuclei via a π + -meson to yield a neutron, which can easily escape the potential wall. The resulting processes are also described by equations (30 -32).

Discussion and Conclusion
Let us first consider the low energy domain, i.e. up to the resonance maximum Eres. Thus for E < ETh only due to the tunneling effect the projectile proton can enter the interior of the potential wall according to . This provides in the case of oxygen the isotope F9 17 as intermediary excited state (or N7 13 and Sc21 41 with respect to carbon or calcium). Since E < ETh the lifetime of the excited states depends on the transition probability to a lower energy state to finally produce a γ-quant and a neutron by an exchange interaction (Pauli principle) within the nuclei via a π + -meson to yield a neutron, which can easily escape the potential wall. The resulting processes are also described by equations (30 -32). These equations are also valid for comparably low proton energies, but with E > ETh. However, we are now located in the energy domain E > ETh up to the resonance domain Eres. This domain, which is referred to as Breit-Wigner resonances, is characterized by various excitations, e.g. rotations of the whole nuclei, vibrations via proper deformations of the nuclei and excitations to excited configurations. All these processes are damped by emission of photons with proper energies. It should also be pointed out that these resonances belong to the category of inelastic scatter, but the resulting secondary proton is identical with the primary proton. Only the lateral scatter and energy losses are different to Molière scatter processes. If these scattering processes occur without nuclear excitations, rotations and vibrations, then they are referred to as elastic nuclear scatter of the protons, which are also contained in the Breit-Wigner resonance formula and in its generalization 7 , given by Flügge. (The original Breit-Wigner formula 6 is restricted to scatter processes induced by so-called 's-states' of the Smatrix; the inclusion of 'p-states' and states of higher order have later been accounted for 7, 8 ). The very often used way of notation 'inelastic nuclear scatter of protons' is not quite correct, since elastic scatter processes of protons at nuclei even occur in the asymptotic domain of Q tot (E). A very important nuclear process occurring in the same energy domain must be mentioned again, namely the exchange interaction of the proton with the mesons of the nucleus, i.e. we consider the reaction: The recoil nucleus F9 16 undergoes ß + -decay of electron capture (EC) to become finally again O8 16 . With regard to nuclear reactions of therapeutic protons the cross-section of oxygen is certainly most important. However, the interaction with calcium (bone) and proteins (carbon) seems also to be worthy of interest. Below we present lists of the most important reactions with carbon and calcium. For brevity, we do not state the decay reactions of the isotopes (detailed listings can be found in web), but for most of them electron capture (EC) and ß + -decay with emission of γ-quanta are preferred reaction channels. For these reasons we have state the symbols X in equations (31 -31g) and equations (32 -32h).
As already verified 2-3 in the case of oxygen, the nuclear reactions presented in the listings 4.1.1 and 4.1.2 depend on the available or residual proton energy. For this purpose we consider the domain E > Eres up to the beginning of the asymptotic behavior, i.e. the proton energy amounts to ca. 100 MeV. We refer to this energy branch of Q tot (E) as Etemperate. (In the case of copper, the asymptotic behavior is reached at a somewhat higher energy (ca. 150 MeV), but this element is mainly of importance with regard to the determination of beam-line properties). Thus in the energy domain Etemperate the nuclear reactions according to relations (31 -31c) and (32 -32d) preferably occur, whereas for E > Etemperate or E >> Etemperate the release of clusters according to relations (31d -31g) and (32e -32h) are more probable. A typical case of a cluster release is the α-particle as a secondary particle. According to 3 this threshold energy amounts to 100 MeV.
In the case of carbon it is increased to be 101 MeV, for calcium it amounts to 98 MeV and for copper to 97 MeV. At about E = 190 MeV the probability of this release of αparticles as outcome of nuclear reactions vanishes for all cases considered here. The main difference, however, exists in the yield of this secondary particle. We normalize total yield of α-particles to '1' for the reaction of proton with oxygen. Then in the case of carbon we obtain a total yield of 0.97. For calcium it amounts to 1.68 and for copper to 1.93. Therefore it is generally correct that mainly the yield of clusters is increasing, if the nuclear charge Z and mass number AN correspondingly increase. A previous check of GEANT4 10 in the papers 1-3, 13, 15 revealed that this Monte-Carlo code provides default reaction channels, which underrate cluster formations. However, since this Monte-Carlo code is an open system, it can be improved by the user. A further lack of this Monte-Carlo code is that the release of neutrons is not sufficiently accounted for, as we could verify at a glance at the papers. 11, 12 According to the present results the release of low energy neutrons should be much more in focus in dose calculations due to the rather high RBE.
The calculation method for the total nuclear cross-section Q tot (E) provides a significant advantage compared to the previously published method 2,3 : Each Gaussian distribution and the error-function distribution with a specific energy shift can be associated to probability distributions for the occurrence of some nuclear reaction types. Thus the energy shifts now refer to threshold energies with regard to the corresponding maxima and the half-breadths represent a measure for the yield of some reaction types. Furthermore, it is possible to analyze measured curves of Q tot (E) by suitable deconvolution procedures as worked out in. 16 These convolution methods are developed with regard to linear combinations of Gaussian kernels and shifts. Since there exists a connection between the number of kernels and the underlying statistics, which is strictly a non-relativistic Boltzmann distribution in the case of one Gaussian kernel and a linear combination of different Gaussian kernels with further shifts in the case of the Fermi-Dirac statistics, we have also obtained a way to interpret nuclear reactions by the operator formulation of Fermi-Dirac statistics. 2, 16 Such statistical models can help to calculate nuclear interaction processes in a simpler way compared to manybody-problems of relativistic quantum mechanics. The evaluation of measurement data in clinical practice of proton therapy, e.g. the contributions of secondary particles, requires relatively simple but reliable methods.