Journal of Nuclear and Particle Physics

p-ISSN: 2167-6895    e-ISSN: 2167-6909

2012;  2(3): 42-56


A Calculation Method of Nuclear Cross-Sections of Proton Beams by the Collective Model and the Extended Nuclear-Shell Theory with Applications to Radiotherapy and Technical Problems

W. Ulmer1, E. Matsinos2

1Dept. of Radiation physics, Klinikum Muenchen-Pasing and MPI of Physics, Goettingen, Germany

2Varian Medical Systems, Baden, Switzerland

Correspondence to: W. Ulmer, Dept. of Radiation physics, Klinikum Muenchen-Pasing and MPI of Physics, Goettingen, Germany.


Copyright © 2012 Scientific & Academic Publishing. All Rights Reserved.


An analysis of total nuclear cross-sections of various nuclei is presented, which yields detailed knowledge on the different physical processes such as potential/resonance scatter and nuclear reactions. The physical base for potential/resonance scatter and the threshold energy resulting from Coulomb repulsion of nuclei are collective/oscillator models. The part pertaining to the nuclear reactions can only be determined by the microscopic theory (Schrödinger equation and strong interactions). The physical impact is the fluence decrease of proton beams in different media, the stopping power of secondary particles, and a ‘translation’ of the results of the microscopic theory to the collective model.

Keywords: Nuclear Cross-sections of Protons, Collective Model of Nuclei, Extended Nuclear Shell Theory

Cite this paper: W. Ulmer, E. Matsinos, A Calculation Method of Nuclear Cross-Sections of Proton Beams by the Collective Model and the Extended Nuclear-Shell Theory with Applications to Radiotherapy and Technical Problems, Journal of Nuclear and Particle Physics, Vol. 2 No. 3, 2012, pp. 42-56. doi: 10.5923/j.jnpp.20120203.04.

1. Introduction

The accurate knowledge of the total nuclear cross-section Qtot resulting from proton – nuclei interactions is a decisive feature of proton-therapy planning, Monte-Carlo and technical applications. Nuclear cross-sections determine the fluence decrease of primary protons, different ranges and scatter angles of secondary particles (mainly protons and neutrons), passage of protons through materials, and the creation of heavy recoils, which usually undergo β+ decay and emission of γ quanta. Since water represents the usual reference medium for the measurement/calculation of Bragg curves, we consider, at first, the behaviour of Qtot for oxygen, which is the most important case in proton radiotherapy.
Figure 1 shows that, for protons, a threshold energy ETh = 7 MeV exists to surmount the Coulomb repulsion of oxygen. At E = 20.12 MeV, Qtot exhibits a resonance maximum and a Gaussian shape in the environment. For E > 50 MeV, Qtot decreases exponentially; at E ≈ 120 MeV, the asymptotic behaviour is reached. The fluence decrease of primary protons Φp[1] can be evaluated by the knowledge of Qtot(E) according to Figure 1:
Figure 1. Total nuclear cross-section Qtot of oxygen[1 - 5]
Formula (1) is based on a method given in[6], which has the following shape:
The fluence of secondary particles Φsp cannot be determined in an easy manner and will be one essential aspect of this communication. In equation (1) τ refers to different kinds of energy/range straggling expressed in terms of Gaussian convolutions for primary/secondary protons. Further parameters of equations (1 – 2) are: CSDA range of protons RCSDA, density ρ of the medium, charge number Z of a nucleus and AN (mass number), initial proton energy E0 at material entrance, and NAvogadro refers to the Avogadro number. Figure 2 clearly shows that the slope of the straight lines depends on the difference E0 – ETh, if E0≥ ETh. If E0 < ETh, the expression (E0 – ETh)f becomes complex and we have to impose the condition Φp = 1. There have been many attempts[7] to fit the fluence decrease of primary protons by an approximated form, which solely depends on RCSDA. Apart from the fact that these fits are unnecessarily complicated, they are rather inaccurate, since they do not involve the threshold energy ETh, which decides, whether a nuclear interaction is possible at all or not.
The goal of this communication is to determine some features of the Qtot in terms of Z, AN and the range of strong interactions RStrong (RStrong = 1.2·10-13·AN1/3 cm) by the toolkits of the collective model and extended nuclear-shell theory. In order to gain the necessary simplicity the latter method is presented in following section.
Figure 2. Decrease of proton fluence in water according to Figure 1

2. Theoretical calculation methods of nuclear interactions

2.1. Monte-Carlo Calculations of Energy Transfer, Fluctuations and Multiple Scatter

The Monte-Carlo code GEANT4 is described in detail in the reference manual [8]. The calculation of the stopping power is based on a numerical manipulation of the Bethe equation with the Bloch corrections (BBE, [9 – 11]). The cut-off amounts to 1 MeV. Since GEANT4 is an open system, all correction terms according to the BBE in accordance with [11] have been implemented. These corrections only imply the removal of singularities, if the total integration procedure is carried out analytically and not by numerical step-by-step calculations of ΔE/Δz; these calculations depend on the actual velocity v and have to be performed for each term of the BBE separately [1]. The scatter of protons is treated by the Molière multiple-scatter theory [12 – 15]. The code also contains a hadronic generator for the simulation of nuclear interaction processes. With respect to proton dose deposition, the basic theory is BBE and a numerical fit of the Vavilov distribution function, which takes account of the Landau tails. In the limit case of fluctuations of small transfer energies from protons to environmental electrons, a Vavilov distribution function assumes a Gaussian shape. In GEANT4, it is also possible to restrict these fluctuations to a Gaussian shape. This fact is of interest with regard to the role of the Landau tails in the initial plateau (entrance region) of Bragg curves. An analysis of Gaussian convolution and its generalizations in physics of protons is given in [5, 16].
The accurate knowledge of the total nuclear cross-section Qtot resulting from proton – nuclei interactions is a decisive feature of proton-therapy planning with advanced models and Monte-Carlo calculations. Nuclear cross-sections determine the fluence decrease of primary protons Φp, different ranges and scatter angles of secondary particles (mainly protons and neutrons), collimator scatter, and passage of protons through bones, implants, materials of technical interest, and the creation of heavy recoils, which usually undergo β+ decay and emission of γ quanta. Since water represents the usual reference medium for the measurement/calculation of Bragg curves, we consider, at first, the behaviour of Qtot for Oxygen.

2.2. Nuclear Interactions of Protons, Release of Secondary Protons and Heavy Recoils

The most important aspect is the hadronic generator and the energy transport of secondary (and higher-order) particles (protons, neutrons, deuterons, etc.[17 – 21]). However, the default nuclear cross-section implemented in GEANT4 is very poor. Instead of using the default routine, which is based on data of [2], we have implemented the cross-section data of O16 of [4]. Furthermore, we have calculated this nuclear cross-section with the help of the extended nuclear shell theory containing, apart from the strong interaction spin-spin and spin-orbit couplings, the electrostatic interaction and the exchange interactions between the nucleons due to the Pauli principle. Thus, the wave-functions of ground and excited states can be calculated by a perturbed SU3, see [22 – 24]. The calculation of the cross-section due to energy transfer by an external proton is based on well-elaborated principles (determination of the transition probability and density of states). At first, let us return to Figures 1 and 2. The decrease of fluence of primary protons can also be calculated from the results of this figure. Protons with energy lower than the threshold energy ETh cannot surmount the potential wall of O16. Figure 1, which presents the total nuclear cross-section, shows that there is a threshold energy ETh = 7 MeV (more accurate: 6.997 MeV), which a proton should have to perform nuclear interactions with the O16. For proton energies lower than the resonance maximum at Eres = 20.12 MeV, the primary proton is preferably scattered by the nucleus (the secondary proton is now identical to the deflected primary one); the nucleus is excited, to undergo rotations/oscillations and emission of X-rays of very low energy (around 1 keV), leading to the release of Auger electrons. A complete classification of the total nuclear cross-section of the proton – nucleus (O) interaction for therapeutic protons is given by the following two types:
1. Potential scatter of protons by the strong interaction potential in the environment of the nucleus (note that Rstrong: ≈ 1.2∙AN1/3 ∙10-13 cm determines a typical distance where the strong-interaction and the Coulomb force ‘balance’ one another; AN = 16). For R > Rstrong, only the Coulomb part is present. Potential scatter accounts for most of the protons undergoing a nuclear interaction. Resonance scattering of the incident protons at the nucleus occurs by inducing transitions between different states of the nucleus (e.g., vibrations leading to intermediate deformations, rotation bands, excited states by changing the spin multiplicity).
2. Nuclear reactions, which produce heavy recoils (see the listing 3 below). These protons are sometimes referred to as reaction protons – or ‘secondary protons’[15 – 21]. For protons in the therapeutic energy domain, the amount of reaction protons is about 1 % - 4.5 %.
Within the total nuclear cross-section, the case 1 above plays the dominant role, if the residual proton energy E is lower than 150 MeV. However, for E > 150 MeV case 2 is dominant. In first order, case 1 is described by the Breit-Wigner formula; for a detailed representation of this formula, see[6]. The Breit-Wigner formula results from the exact scatter theory and the restriction of the general S-matrix to ‘S states’, i.e., to l = 0. This restriction is only valid for light nuclei with rotational symmetry. If the number of neutrons is significantly different from the atomic number Z, the related S-matrix has also to account for ‘higher states’, i.e., for l ≠ 0. The total nuclear cross-section obtained by the Breit-Wigner formula has an elastic and inelastic part. The restrictions of the Breit-Wigner formula are insufficient in our problem; an extension to more resonances has been given in [25]. By taking account for the above-mentioned collective vibrations (oscillations of the nucleus by deformations and distortions) and rotations, we have taken into account all these degrees of freedom. Therefore, the total nuclear cross-section shows, around its maximum at Eres, a Gaussian behaviour for E > ETh, before it exponentially decreases to reach the asymptotic behaviour.
It has to be added that, in all three cases, about 1 – 7 MeV of the proton energy (depending on the deflection angle) is transferred to whole nucleus to satisfy the energy and the momentum conservation in the center-of-mass system. This implies that for a neutron release, the proton energy has to be 21 – 27 MeV and not simply 20 MeV. Due to the potential barrier, the energy of the colliding proton has, at least, to be 30 MeV in order to release a secondary proton. These two processes have, however, an exception and may also occur at very low energies via exchange of mesons due to the Pauli principle, as pointed out in the discussion of listing (3).
The reaction protons (which always result from an inelastic process) are closely related to heavy recoil particles; the most probable heavy recoil elements resulting from the nuclear reactions of therapeutic protons are given by the types (1 – 5); types (6 – 7) result from types (1 – 2).
All types of β+-decay emit one γ quantum; its energy is around 0.6 MeV – 1 MeV. The β+-decay of F916 has a half-life of about 20 seconds, and γ quanta are produced by collisions of positrons with environmental electrons. Some of the remaining heavy recoil fragments have half-times up to ten minutes (N15). Since Figure 1 refers to the total nuclear cross-section in relation to the actual (residual) proton energy, we have to add some qualitative aspects to the five different types with regard to the required proton energy: if E < 50 MeV, the type (1) is the most probable case with rapid decreasing tendency between 50 MeV < E < 60 MeV, to vanish for E > 60 MeV. Type (2) also pushes out a neutron, but the incoming proton is not absorbed; the required energy amounts, at least, to 50 MeV. Type (3) is similar, but requires, at least, about 60 MeV, with probability increasing with the energy.
The release of α-particles, resulting from clusters in the nucleus, requires an energy E ≈100 MeV and the probability is increasing up to E ≈ 190 MeV; thereafter, it is decreasing rapidly, since higher-energy protons destroy these clusters by pushing out deuterons (type 5). Thus, case 5 is energetically possible for E > 60 MeV, but the significance is only increasing for E > 200 MeV. It has to be pointed out that, with regard to the nuclear reactions in listing (3), the exchange interactions between protons and neutrons, which result from the Pauli principle (applied to spin and iso-spin), play a dominant role. If the incident proton hits the nucleus, it is also included in these exchange interactions. According to the Pauli principle, a preferred interaction (resonance) of this proton with an environmental neutron occurred by virtue of the exchange of a π- meson (neutron n of the nucleus + p (incident) → p + n (outgoing)); this implies that a neutron of the nucleus converts to a proton and the incident proton leaves the nucleus as a neutron (reaction type (1) of the listing (3)). The inclusion of nuclear reaction processes to a generalized Breit-Wigner formula according to[25] requires results of the extended nuclear shell theory. It is obvious that secondary protons may again release tertiary protons by additional inelastic scatter or by the types (6 – 7) of the listing (3), where released neutrons are responsible for resonance effects, e.g., type (6) represents the reversal process of type (1). The incoming (secondary) neutron is converted to a proton via π+ exchange. The calculation procedure of the stopping power contribution of reaction protons Ssp,r is presented in the final section. A particular feature of the nucleus is the lack of a central force, unlike the case of atomic electrons. Each nuclear constituent (proton, neutron) is therefore moving in the field of all remaining constituents (self-consistent field). The average kinetic energy amounts to 24 MeV. This fact implies that, due to the spatial charge distribution of the nucleus, the potential barrier with the threshold energy ETh is not proportional to Z2; the exponent has to be somewhat lower, i.e., Z k (k < 2). An analysis of theoretical results and of the nuclear cross-section data of the Los Alamos library resulted in the following connection between ETh and Z (valid for those nuclei, for which the number of neutrons is approximately equal to that of protons):
The constant factor CF amounts to CF = 0.222265, and the form-factor function F(Z, AN) accounts for the total mass/charge distribution of a nucleus, which is defined by:
The parameters of the form-factor function are given in Table 1.
Table 1. Form-factor function F(Z, AN)
From relation (4 - 5), it follows that κ converges to 2, if Z → 1 (proton – proton repulsion); F(AN = 12.01, Z = 6) = 1.328917. With regard to ETh, Figure 3 shows an application of the formulas (7 - 9). Referring to point 1 (potential scatter) and point 2 (resonance scatter), we have evaluated the Breit-Wigner formula for S and P states, and only the nuclear reaction processes producing neutrons, secondary protons, α-particles, deuterons, etc., have partially accounted for by data of the Los Alamos library. An analysis of the presented total nuclear cross-sections Qtot (and of some further cases) available in the Los Alamos library suggests the following adaptation model (Em: the characteristic energy of the maximum value Qtotmax of Qtot, Qtotas: the asymptotic value of Qtot, approximately equal to the geometric cross-section; σres: the half-width of the resonance region; σas: a characteristic value used in the description before reaching the asymptotic behaviour).
The parameters Eres, Qtotmax, σres, Ic (Ic = Qtotc/Qtotmax), Qtotas, and σas are still not defined. A discussion on these parameters (as well as an overview of some theoretical aspects of nuclear physics, e.g., of the nuclear shell theory and extensions will be given in the later section). The parameters required in equations (7 – 9), referring to Figure 3, are given in Table 2. An inspection of Figure 3 and of equations (7 - 9) indicates that two Gaussian distributions are needed between ETh ≤ E ≤ Eres and Eres < E ≤ Ec, this being a result of the different interaction mechanisms of the proton with the nuclei. The parameter σas describes the asymptotic behaviour of Qtot for E> Ec and is determined by the condition that, at E = Ec, equations (7 - 9) have to be compatible (that is, the functions and their first derivatives must be equal). It is known from nuclear (and even particle) physics that inelastic cross-sections show an exponential decrease represented by the sum of some exponential functions before reaching the asymptotic region. We have verified that a single exponential function is not sufficient and that the hyperbolic-tangent (tanh) function provided more accurate results.
Figure 3. The total nuclear cross-section for C, O, Ca, Cu, and Zn (taken from [5])
Based on the results of this section we are now able to calculate the fluence of secondary protons, so far the origin refers to non-reaction protons.
Table 2. Numerical parameters for the evaluation of Figure 3 according to equations (7- 9)

2.3. Secondary Protons Φsp

In this section, we separate the whole number of secondary protons by their origins, i.e., we differ between reaction protons Φsp,r and non-reaction protons Φsp,n. Due to the complexity the contributions of reaction protons will be determined in a later section.
It should be noted that σsp is somewhat different from τ of equation (1). The argument of the error function in equation (10) is slightly changed by the additional zshift, which results from an average energy loss of the secondary protons. The uncertainty intervals in the value v = 0.958 are + 0.40 % and – 0.42 %. Thus, v = 0.958 represents the value with the lowest mean standard deviation. The essential formulas for the determination of some parameters of equations (10) and (15) are given by the following formulas (11 – 14), which are also the result of calculations of the system ‘proton - nucleon’.

2.4. Recoil Protons/Neutrons Φrp

For the recoil factor η, we use η = 0.042. The remaining parameters are the same as in equations (10) and (15). The fluence decrease of primary protons, resulting from Figure 4, is analogous to that of water (Figure 3).
Figure 4. Decrease of the fluence of 100, 200, and 250 MeV primary protons in copper (solid lines) and in calcium (dashes). The RCSDA ranges are stated by perpendicular lines
A comparison of Figures 3 and 4 shows that the slope of the fluence decrease is more bending in the latter case, in particular for energies below 100 MeV. For copper, the threshold energy ETh now amounts to 8.29 MeV and for calcium 7.09 MeV. Thus, for proton energies below ETh, the fluence remains constant. Due to the rather significant fluctuations, we cannot verify the constancy of the fluence at the end of the track. Formula (15) can be applied to Figure 4. However, ETh = 8.29 MeV (Cu) and ETh = 7.09 MeV (Ca) has to be taken into account. The power f = 1.032, which is valid for water to compute the slope of the straight line in equations (10 – 15), has to be slightly modified (f = 0.755 for Cu and f = 0.86 for Ca), and the general formulas are given by:
After these modifications, formulas (16– 17) can be applied. Formula (17) is closely related to the nuclear collective model, where the parameters a, b, c and d are interpreted in terms of different contributions to the total cross section. These parameters have been determined by this model: a = -0.087660001, b = -6.379250217 , c = 5.401490050 and d = - 0.054279999. There are two applications, in which equations (16 - 17) are relevant: 1. Passage of protons through collimators. 2. Passage of protons through materials different from water. In case 2, only a small path length has to be corrected; however, Figure 4 shows that the fluence decrease has also to be corrected to fulfil continuity at the boundaries. With regard to Figure 4 and to the qualitative listing (3) valid for the proton - oxygen nucleus interaction, we should finally point out that by a suitable modification a similar listing for nuclear reactions will be obtained, if the oxygen nucleus is replaced by another one such as calcium or copper. All formulas necessary for the calculation of Qtot are applicable, as long as the rotational symmetry of the nuclei approximately holds.

2.5. Fitting Procedure for the Determination of ETh and Qtot

The calculation of the total nuclear cross section requires some information acquired by fitting the Los-Alamos data, results of extended nuclear shell theory, and empirical rules. Let us at first consider equations (18 – 20) to compute ETh and k. We assume an iso-scalar nucleus, i.e., one in which the numbers of protons and neutrons are equal (AN = 2∙Z). This assumption holds for almost all light nuclei. For nuclei with spherical symmetry, the nuclear radius Rstrong is given by:
Thus for r > Rstrong, the contribution of strong interaction is negligible. The nuclear shell theory with oscillator potential gives for r = Rstrong:
The first term on the right-hand side of this equation represents the Coulomb repulsion of an incoming proton; the second one is the mutual Coulomb repulsion of Z protons in the nucleus. U0 is the depth of the potential and is put equal to AN∙EB∙Rstrong2; CF is a proportionality factor and EB the binding energy per nucleon. EB is equal to 8 MeV, if AN ≥ 12, smaller for AN < 12. We are able to rescale ω0 (if EB is constant) in such a way that U0 vanishes. Taking into consideration that AN = 2∙Z and multiplying both sides with (2∙Z)1/3, we finally can summarize the following results:
A least-squares fit of all available Los-Alamos data yielded k = 1.659, instead of k = 1 + 2/3. This might result from crude assumptions in the creation of our fitting model: we have assumed Mproton = Mneutron and, furthermore, neglected the spin-orbit coupling. The equation above may also provide the means for the calculation of the rescaled ω0. If the number of neutrons is slightly different from Z (i.e., AN = 2∙Z + εN and εN << Z), then a correction term is required, which is already taken into account in the form factor function F(Z, AN) according to equation (6); the rotational symmetry still has approximately to hold. With regard to the determination of the total nuclear cross section Qtot, we have borrowed some elements from the collective model of nuclear interactions. First of all, we have to know Qtotmax required for the calculation of other quantities. We start with the following ‘Ansatz’, and thereafter we give some explanations:
What is the physical interpretation? With the aid of formula (20), we obtain the following properties:
Term a: Connection of Qtotmax to the complete volume of the nucleus. It is important in the resonance domain; it includes resonance scatter via nuclear deformations (vibrations) of the whole nucleus, resonance excitations by changing the spin multiplicity (all effects are inelastic), and transformation of a nuclear neutron according the listing 3 (case 1, inelastic).
Term b: Proportional to the area of the geometric cross section. It contains potential scatter (major part, elastic), rotations induced by Coulomb repulsion/strong-interaction attraction (elastic and inelastic), and nuclear reactions by changing the iso-spin multiplicity (inelastic).
Term c: Proportional to the nuclear radius Rstrong. Excitations by spin-orbit coupling, when the whole nucleus changes its angular momentum, inelastic resonance effect, and elastic spin-spin scatter.
Term d: Proportional to Zk/Rstrong. Excitation of nuclear vibrations by Coulomb repulsion (resonance effect, inelastic) and elastic scatter.
It should be pointed out that formulas (19 – 22) stand in close relationship to the early mass formula of Bethe-Weizsäcker. This formula is thoroughly discussed in [6, 22 – 24].
The asymptotic behaviour Qtotas of Qtot is given by the relation:
This connection mainly contains the term b above. We have verified the validity of this property for nuclei up to Zn. The order of magnitude of the term b is a clear indication that elastic potential scatter of the nucleus via the strong interaction is the main contribution of the total nuclear cross section. However, equation (21) can also be used for the determination of Qtotc and Qtotas, and only the four coefficients are different. Therefore we write equation (21) in a modified form (the parameters are given in Table 4):
Table 4. Parameters a, b, c and d for some different types of Qtot
Eres is given by Eres = Em + ETh. Results obtained by using the extended nuclear theory and Los-Alamos data indicate the following connection:
The parameter σres, previously introduced, results from a fitting procedure of calculated and measured data. Ic is defined by the continuity condition for the Gaussian and the hyperbolic-tangent function. It is known that the nuclear cross-section can be described by a series of exponential functions, before the asymptotic behaviour is reached. We have verified that the hyperbolic-tangent function, which can be expanded in terms of exponential functions, provides optimal results, and it easily accommodates the continuity conditions at E = Ec and Qtotc = Qtotmax∙Ic. The Gaussian behaviour in the resonance domain is due to the numerous resonance excitations occurring at E ≈ Eres according the generalized Breit-Wigner formula[25].

3. Results of the Generalized Nuclear Shell Theory

3.1 Harmonic Oscillator Models

This section is appropriate only for interested readers. At first, we consider the 3D harmonic oscillator, described by the Hamiltonian Hosc:
In this equation, iso-spin symmetry is assumed to hold, i.e., Mp = Mneutron = M. There are three ways to obtain the general solution of this equation, well-known from standard textbooks of quantum mechanics and nuclear physics; herein, we only present the results. 1. Use of creation and annihilation operators (algebraic method) based on the commutation relation:
2. Replacement of pk by -i∙ħ∙∂/∂qk in the Hamiltonian and solving the resulting Schrödinger equation by a Gaussian function multiplied with Hermite polynomials.
3. Solving the Schrödinger equation in terms of spherical harmonics and Laguerre polynomials.
First method
We rewrite the Hamiltonian of the 3D oscillator as:
These operators obey the commutation relations for bosons:
With the help of these operators, the Hamiltonian assumes the shape:
The angular-momentum operator commutes with the Hamiltonian and, therefore, it only connects degenerate states of the Hamiltonian H, by transforming a quantum state k to the state l and vice versa. The operator b+k (absorption operator) and bk (emission operator) modify (increase and decrease, respectively) the energy ħ∙ω0. There are nine independent types of bilinear products b+k∙bl (i.e., k = 1, …, 3 and l = 1, ..., 3), which implies that they can be the generators of SU3 in the configuration space. This means that there is a correspondence between SO3 (rotational symmetry in the configuration space) and SU3, in analogy to the one between the group SO2 (x/y – plane) and SU2 for the two-dimensional harmonic oscillator. In nuclear physics, the group SU2 is connected to the iso-spin, referring to both nucleons obeying anti-commutation rules. Although bilinear products of Fermion operators satisfy the above commutation rules, physical differences exist. Applied to a whole nucleus, the oscillator model is rather a collective description of physical properties as oscillations/vibrations via deformation or creation of rotational bands (quanta of the angular momentum of the whole nucleus) due to interactions with comparably low-energy protons (e.g., see resonance scatter of Qtot, in particular the first term of equation (32)). A further critical aspect is that the oscillator potential has a minimum for E = 0, not for E << 0; bound states exist for arbitrarily high energies. Nevertheless, with the help of some modifications this model will become suitable for practical problems.
Second method
This method is the configuration-space representation of the first method and implies the corresponding solution of the Schrödinger equation in Cartesian coordinates:
We make use of the substitutions
The complete solution is then given by:
N is a normalization factor; Hj, Hk, and Hl are Hermite polynomials as already previously defined[5, 26]. The advantage of this representation is the possible connection to other problems, involving Gaussian functions.
Third method
This method uses the framework of the H atom, i.e., the separation of the wave-function ψ in a radial function and spherical harmonics. The only difference is that the potential Z∙e02/r is replaced by the oscillator potential; this implies that the Laguerre polynomials have an argument different to that of the H atom:
The energy eigen-values are independent of the quantum number m (m = -l, -l+1, …, 0, l-1, l). The correspondence between Hermite polynomials and spherical harmonics may be found in[26].
The concept of spin can be introduced to the 3D harmonic oscillator by the commutation rule, see[27]:
In this case, σ represents the three Pauli spin matrices and the unit matrix. The result is the Pauli equation of a 3D harmonic oscillator. Similarly, we introduce the iso-spin τ by the substitution:
Figure 5. Total (effective nuclear potential plus Coulomb repulsion) for O
By that, we obtain the Pauli equation for a 3D oscillator for a proton and neutron. In nuclear many-particle theory, the Pauli principle is generalized, i.e., the total wave-function has to be anti-symmetric with regard to spin and iso-spin. The property gp/gn = -3/2 does not follow from the scope of the present theory and further principles have to be introduced[28], which we do not consider here. Even by extending equations (36) to a many-particle equation (including the spin-orbit coupling) and to a Slater determinant (Hartree-Fock: ground state), the problem of nuclear reactions due to the properties of the harmonic oscillator potential Vosc, cannot be solved. The problem of simple oscillator models can be verified in the following Figure 5, which shows the effective nuclear potential energy for O.
The abscissa is expressed in units of r = 1.2∙10-13∙AN1/3 cm (AN = 16). This type of Pauli equation reads:
We approximate the potential according to:
By that, the complete potential function can be expressed as a linear combination of two Gaussians:
A property of the Gaussian function is that its curvature changes sign at r = rc. For a single Gaussian, as the first one in equation (38), rc is given by:
Only for r ≤ rc, is a harmonic-oscillator approach useful, and the deviation to a Gaussian in this domain is small. This is, however, not true for r > rc. In the case of the linear combination of two Gaussians, rc is broader:
Figure 5 shows that, for r < 0.4982, strong interactions are dominant (all other interactions are negligible). If 0.4682 ≤ r < 1, strong interactions are still present with decreasing tendency, whereas Coulomb repulsion is increasing; finally, for r = 1, strong interactions are negligible. We will come back to these results in the next section.

3.2. Nonlinear/Nonlocal Schrödinger Equation, Inharmonic Oscillators with Self-interaction and Hartree-fock Method (Inclusive Configuration Interaction)

Let us first consider the usual Schrödinger equation for a bound system:
During the past decades, this type has been encountered in many fields of physics, such as superconductivity, nuclear and plasma physics, e.g. [29 – 31] and various other references. A nonlinear Schrödinger equation is obtained by introducing the potential φ, proportional to the density of solutions:
The coupling constant λ is negative (in which case, the solutions are bound states with E < 0); equation (42) can be interpreted as a contact interaction. It is known from many-particle problems (e.g., quantum electrodynamics, Hartree/Hartree-Fock method, etc.), that the mutual interactions between the particles in configuration space lead to nonlinear equations in quantum mechanics. However, in these cases, there are not at all contact interactions; the nonlinear Schrödinger equation above is an idealistic case. By taking ε → 0, the Gaussian kernel is transformed into a δ kernel:
The nonlinear/nonlocal Schrödinger equation can be interpreted as a self-interaction of a many-particle system with internal structure, and it is possible to generalize this type by incorporation of additional internal symmetries (e.g., the introduction of the spin to obtain spin-orbit coupling, SU2, SU3, and also discrete-point groups).
According to the principles developed in [5], we are able to write equation (43) in the form of an operator equation (the Gaussian kernel is Green’s function):
Expanding this operator in terms of a Lie series and keeping only the terms up to Δ, equation (44) becomes a stationary Klein-Gordon equation, which describes the interaction between the particles obeying the Ψ-field:
By rescaling the Klein-Gordon equation, we obtain the more familiar form: 1 + 0.25 ε2 Δ → k2 + Δ; Green’s function is of the form:
By setting k → 0, the Poisson equation of electrostatics is obtained, if is interpreted as a charge density. The Gaussian kernel K also represents the exchange of virtual particles between the nucleons. In view of this fact, we point out that we have incorporated a many-particle system from the beginning. Which information now does this nonlinear/nonlocal Schrödinger equation provide? In order to obtain a connection of the combined equations (44 - 47) with the oscillator model of nuclear shell theory, we analyse the kernel K in detail. In the Feynman-propagator method [32 – 33] the expansion of K in terms of generating functions is an important tool:
Inserting this expression into the nonlinear/nonlocal Schrödinger equation provides:
The equation above represents a highly anharmonic 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 (domain with positive curvature), the whole equation is reduced to a harmonic oscillator with self-interaction; the higher-order terms are small perturbations. We summarize the results and refer to previous publications[28 – 31]):
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 amounts to: . This energy level is lowered by the term ~λ ∙Φ0,0,0, depending on the ground-state wave-function. The energy difference between the ground and the first excited state amounts to ; this is not true in the case above, since the energy level of the excited states depends on the corresponding eigen-functions (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, iso-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:
Ψ is now (at least) a Pauli spinor (i.e., a two-component wave-function), and together with Hso the SU3 symmetry is broken. We should like to point out that the operation acts on the Gaussian kernel K:
The expression in the bracket of the previous equation represents a vector, and p (p → -i) acts on the wave-function. Since the neutron is not a charged particle, the spin-orbit coupling of a neutron can only involve the angular momentum of a proton. In nuclear physics, these nonlinear fields are adequate for the analysis of clusters (deuteron, He, etc.), and has been extended the theory to describe nuclei with odd spin [31]. The complete wave-function Ψc is now given by the product of a function in configuration space Ψ multiplied with the total spin and iso-spin functions.
We should like to add that an extended harmonic oscillator model with tensor forces has been regarded in [22]. The application of oscillator models in nuclear physics goes back to [34]. In [32 – 33] it has been verified that the use of Gaussians in the description of meson fields provides many advantages compared to the Yukawa potential (Green’s function according to equation (43)).
In a final step, we consider the generalized Hartree-Fock method to solve the many-particle problem. In order to derive all required formulas, it is convenient to use second quantization. The method of second quantization is only suitable to derive the calculation procedure: extension of the Pauli principle to iso-spin besides spin, inclusion of spin-orbit coupling, and exchange interactions. This is the consequence of dealing with identical particles, in which case every state can only occupy one quantum number. In order to get numerical results (i.e., the minimum of the total energy of an ensemble of nucleons, the extraction of the excited states, the scatter amplitudes, etc.), we have to use representations of the wave-function by at least one determinant in the configuration space. In the ‘language’ of second quantization of fermions, we would have to regard expressions like:
The operators of the form ak+ and ak (k being a set of quantum numbers) are creation and destruction operators in the state space. The nonlinear/nonlocal Schrödinger equation with Gaussian kernel for the description of the strong interaction, including the spin-orbit coupling, can be written by these operators, leading from an extended particle with internal structure to a many-particle theory. Before we start to explain the calculations by including one or more configurations, we recall that, according to Figure 5, we have an increasing contribution of the Coulomb repulsion for r > rc,, though in the domain r < rc, the contributions of the Coulomb interactions are negligible. Since all basis elements of the calculation procedures, i.e., the calculation of eigen-functions in the configuration space, two-point kernels of strong interactions between nucleons, and the spin-orbit coupling can be expressed in terms of Gaussians and Hermite polynomials, we want to proceed in the same fashion with regard to the Coulomb part. According to results of elementary-particle models[28], the charge of the proton is located in an extremely small sphere with radius rp = 10-14 cm, not at one ‘point’. Therefore, we write the decrease of the proton Coulomb potential by 1/(r+rp); for r = 0, we then obtain 1014 cm-1, not infinite. In a sufficiently small distance of r = 2.4∙10-13 cm, we can approximate the Coulomb potential with high precision by:
The mean standard deviation amounts to 10-5, if the parameters of formula (55) are chosen as:
If necessary, it is possible to rescale r0, r1, and r2 by dividing by (AN)1/3. The contribution with c2 incorporates a long-range correction.
In the absence of an external electromagnetic field, the Hamiltonian reads as:
Note that it is possible to distinguish between the proton and the neutron masses by indexing M; the ε, previously used in equation (57), has been replaced by σs. The coupling constant of gs is 1, if the Coulomb interaction is scaled to:
Thus, in theoretical units with e0 = c = h/2π = 1, the coupling constant gs assumes 137. This relation can be best seen in the Dirac equation containing a Coulomb repulsion potential ~ e02 and a strong interaction term ~ -gs. The aforementioned relation is obtained by dividing the kinetic-energy operator c∙α∙p → - c∙α∙ħ∙v and β∙mc2 by (c∙ħ). In the calculations for deuteron, He3, and He, we have assumed the range length σs:
This assumption turned out to be not sufficient; a replacement of σs was justified to distinguish between the range length σsp (π-mesons) and σsk (K-mesons):
The range length σsk is proportional to 1/mk (mk: mass of the K-meson). The Hartree-Fock method provides the best one-particle approximation of the closed-shell case.
The one-particle functions φk1(1), …, φkN(N) contain all variables (configuration space of position coordinates, spin, and iso-spin). By using a complete system of trial functions, e.g., a Gaussian multiplied with Hermite polynomials, the Hartree-Fock limit is obtained. In view of our question to calculate the S-matrix and the cross-section of the proton-nucleon interactions (elastic, inelastic, resonance scatter, and nuclear reactions), this restriction is insufficient. In particular, we have to add excited configurations and virtually-excited configurations. The role of excited states is clear. As an example, we regard the O nucleus, where the total spin is 0. If a proton or neutron of the highest-occupied shell is excited, then the spin may change, and both, highest-occupied and lowest-unoccupied shell, are occupied by one nucleon. The emitted nucleon may be regarded as a ‘hole’. This procedure can be repeated to higher-unoccupied states and to linear combinations of configurations with different nucleon numbers. A virtually-excited state is produced, if the configuration of the excited state only formally exists for the calculation procedure, but cannot be reached physically. An example of this case is already the deuteron with iso-spin 0 and spin 1. An excited state with spin 1 or 0, where proton and neutron occupy different energy levels (shells), does not exist. In spite of this situation, the Hartree-Fock method does not provide the correct ground state, and linear combinations of determinants with different spin states (S = 1, -1, 0) and ‘holes, i.e. virtual states’ have to be included. These virtual states also enter the calculation of the S-matrix and of the cross-section.
We have previously performed[5] Hartree-Fock - configuration-interaction calculations (HF - CI) for the nuclei: deuteron, He3, He, Be, C, Si, O, Al, Cu, and Zn. The set of basic functions comprises 2∙(AN + 13) functions with the following properties:
Both α1-functions and α2-functions are chosen such that the number of functions is AN + 13. The different range parameters α1 and α2 are useful, since different ranges can be accounted for. If α >> β, the related wave-functions decrease much more rapidly (central part of the nucleus), whereas the β-contributions preferably describe the behaviour in the domain r ≥ rc. With the help of this set of trial functions (Ritz’s variation principle), we obtain the best approximation of the total energy E by Eapp and the nuclear shell energies (for occupied and unoccupied shells). For bound states, Eapp > E is always fulfilled. It should be noted that for computational reasons it is useful to replace the set of functions (63) by the non-orthogonal set:
By forming arbitrary linear combinations depending on α1 and α2 we obtain the same results as by the expansion (63). The exploding coefficients of the Hermite polynomials are
an obstacle in numerical calculations and can be avoided by the expansion (63). The minimal basis set for the calculation of deuteron would be one single trial function, i.e. a Gaussian without further polynomials. This is, however, a crude approximation and already far from the HF limit. Using this simple approximation, we obtain the result that the ground state Eg depends solely on α1. The best approximation exceeds the HF limit by about 15 %. Various tasks, such as resonance scatter, nuclear reactions, and spin-orbit coupling cannot be described; the cross section of the pure potential scatter is also 12 % too low.
Using 13 + AN α1-dependent and 13 + AN α2-dependent functions, we have obtained the HF limit and virtually-excited states (a bound excited state does not exist). The HF wave-function had to be subjected to virtually-excited configurations, i.e., all possible singlet and triplet states. This calculation had to be completed by introducing a further proton (interaction proton) and including all virtual configurations (besides a configuration with three independent nucleons, a configuration of a virtual He3 state). Thus, for low proton energies (slightly above ETh), the He3 formation is possible. The exceeding energy can be transferred to the total system and/or to rotations/vibrations of He3. In the same fashion, we have to proceed to the calculations for other nuclei: the configurations of all possible fragments have also to be taken into account. (The cases, corresponding to the O nucleus, are given in listing 3). In order to keep these considerations short, we now only give a skeleton of the calculation procedures, which are necessary to evaluate the cross sections. When – besides the ground state – all excited states (including virtually-excited states and configurations of fragments) are determined (wave-functions and related energy levels), then Green’s function is readily determined by taking the sum over all states. This function contains all coordinates in the configuration space (including the spin), quantum numbers of oscillations, and rotational bands:
The S-matrix is given by:
The transition matrix Tkl is defined by all transitions with k ≠ l:
In order to determine the differential cross section, we need the transition probability. For this purpose, we assume that, before the interaction of the proton with the nucleus, this nucleus is in the ground state. Thus, it might be possible that a proton produces excited states of the nucleus by resonance scatter (inelastic), and a second proton hits the excited nucleus before the transition to the ground state (by emission of a γ-quant) has occurred. The second proton would require a lower energy to release either a nucleon or to induce a much higher excited state of the nucleus. However, due to the nuclear cross-section, the probability for an inelastic nuclear reaction is very small and would require a very high proton density to yield a noteworthy effect. Therefore, we have calculated the transition probability using the assumption that the occupation probability of the ground state P0 is 1, i.e., P0 = 1 and Pk = 0 (k > 1). (This is very special case of the Pauli master equation). The differential cross section is obtained by the transition probability divided by the incoming proton flux:
At lower energies, this flux could be calculated by the current given by the Schrödinger equation. To be consistent, we have always used the Dirac equation, since proton energies E > 200 MeV show a significant relativistic effect. With regard to the incoming proton current, we have to point out an important feature:
The Breit-Wigner formula only considers S states and the incoming current is along the z direction. The generalization of this formula[25] includes P states, but the incoming beam is also restricted to the z direction.
Since for our purpose it is necessary to take account for the x/y/z direction by kx, ky, kz in the Dirac equation, we have not yet succeeded in obtaining a compact and simple analytical form.

4. Applications to Reaction Protons of the Inelastic Cross-section Ssp,r

We have already pointed out that the main purpose for calculations with the extended nuclear shell theory incorporate nuclear reaction contributions of protons, neutrons and further small nuclei to the total nuclear cross sections of nuclei discussed in this presentation. We should also mention that the default calculation procedure of nuclear reactions in GEANT4 is an evaporation/cascade model, which has been developed on the basis of statistical thermodynamics.
Figures 2 - 4 do not yet provide final information about the contributions Ssp,n and Ssp,r. The first case of non-reaction protons has already been treated. According to Figure 6 the contribution of reaction protons is particular important for E > 150 MeV with increasing energy. We now present the calculation formulas for this case. Thus Ssp,r is proportional to Φ0∙2∙υ∙Cheavy and a function Fr, depending on some further parameters. It should be mentioned that the parameters of equations (21 - 23) require modifications to be valid for various types of nuclei. However, from Figure 5 the corresponding parameters of some further nuclei can be verified, e.g., ETh, Eres, and some necessary information on the total cross-section. We use the following definitions and abbreviations:
Formulas (69 - 70) can only partially derived, and the adaptation to computed data with the help of the extended nuclear-shell theory is also needed.
The transport of secondary reaction protons resulting from the spectral distribution of these protons has to be taken into account; and the spectral distributions rather obey a Landau than a Gaussian distribution.
The result is the following connection:
The tails at z ≥ RCSDA result from tertiary protons induced by neutrons and the resonance interaction via meson exchange as pointed out in a previous section.
Some consequences of these contributions with regard to build-up have been thoroughly discussed in[16], since the role skew symmetric energy transfer (Landau distributions) and energy transfer from reaction protons along the proton track represents a principal question in understanding the physical foundation of Bragg curves. In Figure 6 we present the stopping power of the ‘reaction particles’ created by Oxygen in water as a function of the initial energy E0. The corresponding results with regard to some further nuclei of interest can be verified in Figures 7 – 10.
Figure 6. Stopping power of reaction protons (medium: water) in dependence of the initial energy E0
With regard to other media equations (69 - 70) are still applicable, if the parameters in equation (70) have been determined. This can be verified successfully by Figures 7 – 9.
Figure 7. Nuclear reaction part of the cross-sections of the nuclei C, O, Ca, Cu and Zn
Figure 8. Calculated cross-sections for Cu and Zn and smooth curves
Figure 9. Stopping power (proton - copper and comparison with water) of secondary/tertiary protons for the initial proton energies 100 MeV, 160 MeV, 200 MeV and 250 MeV
A final example with some technical importance is the nuclear cross-section Qtot of Cs137 according to Figure 10, which corresponds, in principle, to the previously shown cases. ETh amounts to 8.42 MeV (Cs236: 8.45 MeV).
All qualitative explanations of listing (3) can be transferred, inclusive the resonances obtained by nuclear interactions without creating secondary particles (e.g. Breit-Wigner formula). The principle difference to the previous cases like Copper or Oxygen results from surmount of neutrons, since in both cases of Figure10 Z = 55 is valid. By that, the excitation energy of a neutron of the highest occupied shell to the continuum is reduced to 9.2 MeV. Therefore the nuclear reaction ‘proton + Cs137 → Ba137 + neutron’ is able to start at a proton energy of about 10 MeV (please note that the proton energy is partially transferred the deformations of the nucleus with emission of γ-rays). This type of nuclear reaction corresponds to the case 1 of the listing 3, i.e. the incoming proton is converted to an outgoing neutron due to the exchange of a π- meson with the nucleus resulting from the exchange interaction due to the Pauli principle. Because of the general interest of this nuclear reaction we state some analogous cases of the listing 3 (nuclear reactions, which require a higher amount of energy such as the release of clusters, are not presented here):
Figure 10. Nuclear cross-section of the both nuclei Cs137 (dashes) and Cs136 (dots)
Since the lifetime of Cs137 amounts to ca. 30 years, and the corresponding value of some isotopes with less neutrons is comparably small (about 20 days), we are interested in the nuclear reactions ‘proton + Cs137 → Cs136 + proton + neutron’ or ‘proton + Cs137 → Cs135 + proton + two neutrons’. These reactions are very probable at proton energies between 25 MeV and 50 MeV. Since Cs135 exhibits an extremely long lifetime, only by further collisions with protons can lead to a removal of this isotope. Thus we can conclude that the reaction types (71a) or (71b) are very desirable with regard to an artificial change of the nuclide Cs137, which appears in many technical problems.
Thus it appears that the presented methods, which represent a synthesis between extended microscopic calculations and collective model, can provide a significant tool in various domains.

5. Conclusions

The presented results show that a suitable combination of the collective model with extended nuclear shell theory can be adequate to solve problems, which are rather outstanding in many practical problems. Besides the radiotherapy with protons it should be mentioned that the cross-sections of those nuclei/isotopes are important to reduce the half-times of the corresponding isotopes significantly. The storage of long-existing isotopes should be avoided. However, a discussion of technical details of these processes goes beyond the scope of this study.


[1]  Ulmer, W. ‘’Theoretical aspects of energy range relations, stopping power and energy straggling of protons’’ Rad. Phys. & Chem. 76, 1089 – 1107, 2007.
[2]  Berger, M.J., Coursey, J.S, Zucker, M.A. ESTAR, PSTAR and ASTAR: Computer Programs for calculating Stopping-Power and Range Tables for Electrons, Protons and α-particles (version 1.2.2) National Institute of Standards and Technology, Gaithersburg, MD 2000.
[3]  Berger, M.J. A Proton Monte Carlo transport program PTRAN National Institute of Standards and Technology, Report NISTIR-513, 1993.
[4]  Chadwick, M.B., Young, P.G. Los Alamos National Laboratory Report LAUR-96-1649, 1996.
[5]  Ulmer W. and Matsinos E. ‘’Theoretical methods for the calculation of Bragg curves and 3D distribution of proton beams’’ Eur. Phys. J. Special Topics 190, 1 – 81, 2010.
[6]  Segrè, E. Nuclei and Particles (Benjamin, New York) 1964.
[7]  Janni, J.F. Atomic Data & Nuclear Data Tables 27, 147 - 200, 1982.
[8]  GEANT4 Documents 2005; access available from the internet address: /Overview/ html/.
[9]  Bethe, H.A. ‚‘‘Zur Stossbremsung geladener Teilchen beim Durchgang durch Materie‘‘ Ann. Physik 5, 325 – 348, 1930.
[10]  Bloch, F. ‚‘Bremsung rasch bewegten Teilchen beim Durchgang durch Materie‘‘ Ann. Physik 16, 285 - 292, 1933.
[11]  ICRU Stopping powers and ranges for protons and α-particles ICRU Report 49 (Bethesda, MD), 1993.
[12]  Bethe, H.A. ‘’Molière’s theory of multiple scattering’’ Phys. Rev. 89, 1256 - 1262, 1953.
[13]  Molière, G. Z. ‚‘Theorie der Streuung schneller geladener Teilchen III‘‘ Naturforschung 10a, 177 – 200, 1955.
[14]  Ciangaru, G., Polf, J., Bues, M., Smith, A. ‘’Benchmarking analytical calculations of proton doses in heterogeneous matter’’ Med. Phys. 32, 3511 – 3519, 2005.
[15]  Boon, S.N. ‘’Dosimetry and quality control of scanning proton beams’’ PhD Thesis, Rijks University Groningen 1998.
[16]  Ulmer W. and Schaffner B. ‘’Foundation of an analytical proton beamlet model for inclusion in a general proton dose calculation system’’Radiation Physics and Chemistry 80, 378 – 402, 2011.
[17]  Rahu, M.R. Heavy Particle Radiotherapy (Academic Press, New York), 1980.
[18]  Paganetti, H. ‘’Nuclear interactions in proton therapy: dose and relative biological effect distributions originating from primary and secondary particles’’ Phys. Med. Biol. 47, 747 - 762, 2002.
[19]  Jiang, H., Paganetti, H. ‘’Adaptation of GEANT4 to Monte Carlo dose calculations based on CT data’’ Med. Phys. 31, 2811 – 2823, 2004.
[20]  Yan, X., Titt, U., Koehler, A.M., Newhauser, W.D. ‘’Measurement of neutron dose equivalent to proton therapy patients outside of the proton radiation field’’ Nucl. Instr. Meth. A 476, 429 - 434, 2002.
[21]  Polf, J.C., Newhauser, W.D. ‘’Calculations of neutron dose equivalent exposures from range-modulated proton therapy beams’’ Phy,s. Med. Biol. 50, 3859 - 3873, 2005 8/2005. e-Pub 8/2005.
[22]  Elliott, J.P The Nuclear Shell Theory and its relation with other models. Wien, International Atomic Energy Agency, 1963.
[23]  Evans, R.D. Atomic Nucleus (Robert E. Krieger, Malabar, FL, 1982), 1962.
[24]  Meyer-Goeppert, M., Jensen J. Nuclear Shell Theory (Wiley, New York), 1970.
[25]  Flügge, S. ‘’Zur Erweiterung der Breit-Wigner Formel auf höhere Quantenzustände‘‘ Z. Naturforschung 3a, 97 – 110, 1948.
[26]  Abramowitz, M., Stegun, I. Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables Natural Bureau of Standards, 1970.
[27]  Ulmer, W., Hartmann, H.’’ On a geometric interpretation of Heisenberg’s uncertainty relation and the algebraic structure of the Pauli equation’’ Nuovo Cimento 47A, 59 – 73. 1978.
[28]  Feynman, R.P. Photon-Hadron-Interaction (Benjamin, New York), 1972.
[29]  Ulmer, W. ‘’On the application of a Gauss transformation in nonlinear quantum mechanics II: Problems with spin’’ Nuovo Cimento 51A, 309 - 323, 1979.
[30]  Ulmer, W. ‘’On the representation of atoms and molecules as a self-interacting field with internal structure’’ Theoretica Chimica Acta 55, 179 - 201, 1980.
[31]  Milner, R.G., ‘’New physics with polarized HeIII targets’’ Nucl. Phys. A588, 599C – 615C, 1990.
[32]  Feynman, R.P., Hibbs, A.R. Quantum mechanics and path integrals (Mac Graw Hill, New York). 1965
[33]  Feynman, R.P. Quantum electrodynamics - Lecture Notes and Reprints (Benjamin, New York. 1962
[34]  Heisenberg, W. ‘’Zum Aufbau der Atomkerne und der starken Wechselwirkung‘‘ Z. Physik 96, 473 - 488, 1935.