Validation of nuclear reaction models for incident α-particles

2.1 Absorption cross sections

Common to all models based on level densities as those studied here, are reaction cross sections as input. They enter the calculations twofold: the exciton model as well as the compound nucleus evaporation model are based on the principle of detailed balance. So the absorption cross section of the α-particle for the absorption in the initial state is entering the calculations, and of the neutron or proton cross sections enter the inverse reactions. They can be derived in principle from optical model calculations. Nolte et al. [12] derived a generalized set of optical model parameters for α-particles. Limitations of this approach are first neglecting surface effects, thus results for small energies might be questionable. Second the use of the Born approximation is another reason for uncertainties. Thus the results will not account for data beyond a certain energy. Target nuclei in their ground state have a rather large transparency. However, the inverse cross section deals with a highly excited system and the transparency is therefore reduced [13]. We will ignore this effect in the calculations. A versatile method to have absorption cross section at hand is to fix the parameters of empirical expressions. Examples of such expressions are given in Refs [14], [15]. Both are incorporated in the code PHITS which we use here.

For the exciton model calculations of the proton energy spectra we used a different method. Bauhoff [16] had compiled all proton absorption cross sections up to 1 GeV. Carlson [17] fitted the cross sections for some selected energies by

where At is the target nucleus mass number, and r0 and R are energy dependent fit parameters. Machner and Razen [18] fitted smooth curves to the energy dependent quantities deriving r0(Ep) and R(Ep). In principle, one might use these results as input for the calculations. However, the Carlson fits are poor for light nuclei and for heavy nuclei at small proton energies. We, therefore, employed the tunneling approach of Wong [19] for energies in the vicinity of the Coulomb barrier. The so calculated absorption cross sections do not differ much from those calculated by the standard Coulomb correction 1−VC/Ep. Calculations using this approach are compared with experimentally derived cross sections in Figure 1 for some selected target nuclei spanning the Periodic Table of Elements.

Figure 1:

Absorption cross section for protons impinging on the indicated target nuclei as a function of the proton energy. The dots indicate experimental results (full dots) taken from the compilation by Carlson [17]. The solid curve is the prediction of Eq. (1) with energy dependencies given in Ref. [18] with additional barrier penetration. The dashed line in the case of the aluminum target is the prediction of Ref. [14].

The joining point of the two different models is somewhat arbitrary, resulting in a kink for the lighter nuclei. Unfortunately data in the vicinity of the Coulomb barrier are missing and it is therefore impossible to prove the quality of the model for small energies. In the case of aluminum another calculation is in the literature [14] which is also shown in the figure. It clearly extends the range further down in energy than the present calculation. There is one major difference between the two approaches: the present one uses a fixed Coulomb energy while the other one applies an energy dependent one. Since this model is used in the PHITS calculations we may have a closer view to it. The model prediction can be accounted for by the present model choosing a rather large Coulomb radius of 2.4 fm.

The derivation of the absorption cross sections is discussed in Ref. [18] considering geometry, asymmetry and transparency. The model calculation is compared with data in Figure 2

Figure 2:

Absorption cross section for α-particles impinging on the indicated target nuclei as function of the α energy. Lower panel: the sources of the data are: full dots [20], full squares [21], full triangle up [22], closed diamond [23], crossed triangle up [24], crossed triangle down [25], crossed square [26]. The solid line curve is a the model calculation [18], the dashed line curve (red) the optical model prediction [12] and the dotted curve the model taken from [14]. Upper panel: Same as lower panel.

Some of the data shown are measured directly, others were derived from elastic scattering via optical model analyzes.

The model agrees favorably with the data. Also shown are optical model predictions making use of a generalized set of model parameters [12]. This approach obviously fails for energies of a few 100 MeV. This is not surprising since at these energies the validity of Born approximation, on which the optical model is based, breaks down. At such energies the impulse approximation is a better choice. DeVries and Peng [24] performed a Glauber model like calculation for the case of α-particle absorption on lead. This calculation is also shown in Figure 3. It predicts much larger cross sections for high energies than the optical model calculation.

Figure 3:

Same as Figure 2, but for the two different target nuclei. In the lower panel additional data are added: from [12] shown as hour glass, [27] open circle, [28] open square, [29] open triangle up and [30] open star. A semi microscopic model calculation [24] is shown as long dashed curve.

TALYS makes use of the generalized optical model parameters for nucleons with energies below 200 MeV and normalization functions depending on the asymmetry (N−Z)/A [31]. It is worth mentioning that the inverse cross sections given by Dostrovsky et al. [32] are still in use by several evaporation codes following intra nuclear cascades, behave very similar to the Machner-Razen model up to 100 MeV. Higher energy particle emission in the evaporation phase is, however, rare and hence the precision of the model results in this energy range of no importance.

2.2 Exciton model

We will first discuss the exciton model (EM). As stated above states are classified by the number of excitons m and n, i. e. the number of particles above and the number of holes below the Fermi energy. For the sake of simplicity we restrict ourselves to an angle-independent formulation. P(n,t) is the occupation probability of such a state which is measured in terms of the reaction cross section σabs. Transitions between different states are denoted by λ(m→n) and the transition into the continuum by λn↑. The occupation probabilities follow the Pauli master equation

The first sum is a gain term and the second one a loss term. Assuming only two body collisions, the exciton number can change at most by ±2. The differential cross section for a particle of type x is given by

This equation contains three factors which we now discuss. The quantity f(n,x) is a correction factor which takes into account, that the model does not distinguish between proton gas and neutron gas. It is discussed in Ref. [2]. The emission probability WPE,x(n,ε,E) is calculated as in compound nuclear theory from detailed balance

The state densities ρ are those of particle-hole states, the so called Ericson level density,

in a Fermi gas with equidistant spacing 1/g close to the Fermi surface [33]. Although σinv has to be taken on an excited nucleus it is usually assumed to take place on a nucleus in its ground state and σinv=σabs.

The integral in Eq. (3) is the product of life time and depletion factor [34] where usually teq≈∞ is assumed.

There are two limiting cases in the sum in Eq. (3): n0 and n‾. The latter indicates statistical equilibrium, for which creation and annihilation of a particle-hole pair have the same probability. Particle emission from such a state is usually compound nucleus emission. With n0 the initial state is denoted. It is a crucial model parameter since it defines the spectral shape at the high energy end. For a nucleon induced reaction the incident nucleon hits a nucleon from the Fermi sea. Then this state is simply n0=3=2p+1h. However, for complex projectile particles the situation is not that simple. Machner [35] had shown that for heavy ion induced reactions n0∝E with the energy E consisting of two parts: one depends on the binding energy and the other on the kinetic energy per projectile nucleons. Here we are interested in α-particle induced reactions, where the simple model assumptions made in Ref. [35] may not hold. We have therefore extracted n0 for α induced reactions from Griffin plots. The physics behind such a plot is the assumption that the initial state dominates the energy spectrum at forward angles. For U=E−ε−Bx and ε sufficiently larger than the Coulomb barrier, Eq. (3) gives for an Ericson level density

Taking the logarithm on both sides, one has

This relation will show up on a plot as a linear dependence with a slope s=n0−2. It is called a Griffin plot. An example is given in Figure 4. There is clearly an increase of the slope and hence the initial exciton number with beam energy. The deduced value of a slope 10 or initial exciton number 12, is in agreement with the assumption that every projectile nucleon hits a target nucleon and hence n0=8p0+4h0.

Figure 4:

Griffin plot for the indicated reactions. The data in the left panel are from Ref. [36], [37], those in the middle panel from Ref. [38] and those in the right panel from [39], [40].

Chevarier et al. [41] had measured proton spectra from α induced reactions on target nuclei spanning the range from ⁵¹V to ¹⁹⁷Au at a beam energy of 54.8 MeV. In Figure 5 we show the Griffin plots for two target nuclei not too different with respect to their masses. The slope parameters are smaller as one should expect from the small bombarding energy. The difference between the two slopes was explained by Chevarier et al. [41] by odd-even effects: odd- Z targets require n0 = 5 and even Z, n0 = 4.

Figure 5:

Griffin plot for the indicated reactions. The cross section data are from Ref. [41].

Although the Griffin plot points to n0 between four and five, a calculation with n0=5 reproduces nicely the data (see Figure 6). A second chance emission in the pre-equilibrium part is negligible. This is not true for the low energy evaporation part. The deviation between data and model calculation is due to multi particle emission from the equilibrated system, which is not considered in the present computer code. Unfortunately Chevarier et al. did not give angle-integrated cross sections for the ⁵⁹Fe(α,p) reaction, so no comparison can be made. The comparison of the data with TALYS and PHITS calculations will be in section 3.1.

Figure 6:

Proton spectrum from α-particles impinging on an even-even (Z=30,N=36) target nucleus: ⁶⁶Zn. The data (full circles) are from Chevarier et al. [41]. The exciton model calculation with the initial exciton number n0=4p+1h is shown as a solid line curve. The result for a TALYS calculation with n0=5p+1h is shown as broken line curve (see text). The histogram (red) is a Monte Carlo calculation performed with the code PHITS.

Back to the Griffin plot analyses. The finding of Chevarier et al. [41] contradicts the results from Machner et al. [36], who had measured charged particle spectra from α induced reactions at a beam energy of 100 MeV in s-d shell target nuclei. From the Griffin plot analysis they concluded an initial exciton number n0=5 (see Figure 7). A possible explanation might be that the Griffin plot is only one half of the truth. Let us insert the Ericson level density Eq. (5)ρ(p,h,E)=gnEn−1/p!h!(n−1)!, where n=p+h and g is the single particle state density into Eq. (4).

Figure 7:

Griffin plots for the indicated reactions on target nuclei from the s-d shell. The spectra were taken at 20⁰[36]. The slope parameter is s = 3.

In the Fermi gas model is g=A/c with c ≈ 13.3 MeV. The cross section from the initial state is then, when applying this level density

The cross section, therefore, does not depend only on the Griffin term but also on the inverse of the level density.

In Figure 8 we compare data from Machner et al. [36], [37] with different exciton model calculations. The Griffin plot indicates n0 = 5 and this is shown as solid line curve with n0=4p+1h. This curve overestimates the data. The agreement becomes better for n0=5p+1h (short dashed curve), but is in disagreement with the Griffin plot. The best agreement is achieved when n0=4p+1h but the level density is increased by a factor of f=3.2 (short dashed curve).

Figure 8:

Comparison of experimental proton cross sections from the bombardment of an aluminum target with 100 MeV α-particles [36], [37] with results of different exciton model calculations given as curves (see text).

In order to calculate the lifetime of the different exciton states, the internal transition rates λ+ and λ− have to be known. These transition rates can be derived by the golden rule by assuming the same matrix element M for all processes. By applying the Ericson state density with gp=gh=g closed formulas can be derived:

and

with parameters aj, derived from nucleon-nucleon scattering cross section, given in Ref. [34]. As usual, we applied a factor of four to correct the free nucleon -nucleon cross section for scattering in medium.

In TALYS a four component Fermi gas is used throughout. In the Ericson level density a correction to the energy is added to account for the finiteness of the potential well. Koning & Duijvestijn [7] derived a semi-empirical formula for the transition matrix element

The parameters αi were derived by fitting a high number of proton and neutron spectra from proton and neutron induced reactions in the range from 7 to 200 MeV.

Machner [34] had derived a formula for the matrix element, starting from nucleon-nucleon scattering in the range up to 1 GeV [34]. Both approaches are compared with each other in Figure 9 for an exciton number 3. Both approaches contain a dependence M2∝g−3 with g the single particle state density and hence a dependence M2∝A−3. The energy dependence is almost identical except for very small energies. However, in this region the pre-equilibrium fraction is very small. We can therefore conclude that the code of the author and TALYS will yield similar pre-equilibrium spectra. An example of a comparison with the two codes is already given in Figure 6. In both calculations standard input parameters were used. The TALYS spectrum taken from [9] applied n0=5p+1h.

Figure 9:

The energy dependence of the transition matrix element for n = 3. The semi-empirical approach of Koning & Duijvestijn [7] is shown as solid curve. The derivation from nucleon-nucleon scattering labelled Machner [34] is shown as dashed curve.

At high energies the residual system after particle emission is still sufficient highly excited to emit more particles before reaching statistical equilibrium. In this stage further pre-equilibrium particles may be emitted. This is treated in the same formalism as for the first particle, obeying energy, flux and exciton mumbler conservation [42].

At the end of a pre-equilibrium phase the system is in statistical equilibrium: a compound nucleus. Emission from such a system is called evaporation and will be discussed in section 2.4. The TALYS calculation contains always multi-particle emission in this phase which is not the case in the code of the author. The effect of this option on spectral height is discussed below.

Q-values for the reactions were taken from the mass tables [43].

2.3 Intranuclear cascade model

In the cascade codes, which have been constructed to date, it is assumed that the sequence of collisions can be described by classical trajectories, resulting from free nucleon-nucleon interactions. Two important underlying assumptions are: (1) the de Broglie wavelength of the incident particle is small compared to other relevant distances (e. g., the mean free path), so that a classical trajectory is appropriate, and (2) correlations among nucleons enter only through the Pauli principle, and average binding energies, permitting a free scattering description of the system. The direction of the trajectory can change only, if the particle hits another particle or scatters on the potential well, which defines the target nucleus boundary.

Let us first discuss the model assumptions of the Bertini cascade model [44]. The volume of the target nucleus is divided into concentric spheres. Each sphere has a different density and hence a different potential depth. The momentum distribution in each sphere is assumed as a Fermi distribution f(p) with zero temperature, i. e. step like functions. The proton density in each region is set equal to the average value of the charge distribution in that region. The neutron-to-proton ratio is the same in each region. It is equal to the ratio in the whole nucleus. Each region i thus has a Fermi momentum pF,i with

Here Ni is the number of protons or neutrons in the ith zone, ρ the density of protons or neutrons in this zone, respectively. In this model of the Fermi gas the Fermi momentum in the ith zone is related to the density by

The code as discussed so far represents a pure classical model. However, Bertini included one piece of quantum mechanics: a Pauli blocking factor. In a system of Fermions in its ground state all levels with energy below the Fermi energy are occupied, and that the nucleons after being scattered have energies above the Fermi energy. It is assumed that the interactions of projectile with target nucleons do not lead to a state which is seriously different from that.

This model was improved, among by other groups, by Cugnon and others from the University of Liege [46], [47]. Hence, this version is called the intranuclear cascade model Liege (INVL). Characteristics are that initial positions of target nucleons are taken at random in the spherical nuclear target volume with a sharp surface; the initial momenta are generated stochastically in a Fermi sphere, relativistic kinematics is used and inelastic collisions, pion production, and absorption are supposed to proceed via Δ excitation and subsequent decay. Isospin degrees of freedom are introduced for all types of particles, and isospin symmetry is respected; the Pauli principle is treated by statistical blocking factors. This is what in principle is called the Liège INCL model [47]. Different from the Bertini INC model is the density distribution of the target nucleus. A trapezoidal behavior of the nuclear density at the surface is assumed. It agrees much more with the Fermi distributions as obtained from electron scattering [45]. A new parametrization of the nucleon-nucleon scattering data was performed, based on newer experimental data.

By way of example experimental cross sections for elastic np scattering close to 650 MeV were compared with the older INC parametrization labeled Bertini and the new parametrization labelled Cugnon (see Figure 10). Obviously, the latter agrees much better with the experimental results than the former.

Figure 10:

Angular distribution of elastic np scattering. The symbols with asterisk are in the c.m. system. The experiments are from Evans et al. [48] and Terrien et al. [49]. The curves were taken from Refs. [44], [47]. The cross sections were transformed into the present coordinate system.

In order to extend the range of validity towards lower energies, improvements beyond the quasi-classical picture due to the quantum properties of the strongly interacting Fermion systems became necessary. Therefore, in the model calculations the nucleon-nucleon scattering cross sections were modified due to in-medium effects (see [47]). The Pauli blocking was improved compared to the Bertini INC model. Let fi denote the phase space occupation probability in the vicinity of the ith nucleon. It is evaluated by counting identical particles inside a reference volume in phase space around the point (ri,pi). The reference volume is the direct product of a sphere in r-space, of 2 fm radius, and of a sphere in p-space, with a radius of 200 MeV/e. The collision (or decay of a Δ) is blocked when

is larger than a random number chosen between 0 and 1.

The model became thoroughly tested and expanded by the Saclay-Liège collaboration [50], [51], [52] named INCL4.2 to INVL4.6. Here we will concentrate on improvements in the code which are of importance for the present study. First the density profile of the surface was assumed to be the one of the Fermi distribution as discussed above. Whereas INCL treats Pauli blocking within phase space with fixed radii in coordinate space as well as in momentum space, the radii were now calculated in a dynamic manner. If two nucleons i and j are going to suffer a collision at positions ri→ and rj→ leading to a final state with momenta pi→ and pj→, the phase-space occupation probabilities fi are calculated, by counting nearby nucleons in a small phase-space volume,

and similar for the particle j. The factor 1/2 is introduced, because spins are ignored. The sum in Eq. (2.3) is limited to particles with the same isospin component as particle i (or j). The collision between participant i and j is allowed or forbidden following the comparison of a random number with the product similar to Eq. (14)

The centroid of the phase space volume rPB and pPB are free model parameters. They should not be taken too small; otherwise fi(j) is going to be always vanishingly small, nor too large, otherwise the details of the phase-space occupation can be missed. There is no a priori criterion for the appropriate choice of these parameters.

In Ref. [50] extension to incident light clusters among other topics was included. These clusters were approximated in the initial phase by Gaussians in ordinate and momentum space. For ⁴He-particles, in which we are interested here, the radii are <r2>=1.63 fm and <p2>=1.53 MeV/c. The neglect of a nuclear mean field inside the cluster was corrected in the case of ⁴He by a decrease of the incident kinetic energy in order to have the correct total incident energy. In Ref. [51] the emission of light clusters was considered. Finally, in Ref. [52] a detailed procedure for the treatment of light-cluster-induced reactions was included taking care of the effects of binding energy of the nucleons inside the incident cluster. The deficiency of the cluster treatment in Ref. [50], which gives the right total energy but an incorrect too small momentum was corrected for.

It is this model which is incorporated in the computer code PHITS [6]. It is this code we apply to make model calculations.

2.4 Compound nucleus models

At the end of the fast first process, taken into consideration by either the exciton model or the intranuclear cascade model, a thermal equilibrated nucleus remains. Emission from this system, the compound nucleus, is treated by the statistical model. Such model formulations are contained in both approaches. In statistical equilibrium the decay from the compound nucleus is independent of its formation. Hence the cross section can be written as a product of both processes. In its simplest version, a formula similar to Eq. (2.2) is applied to calculate the emission rate per unit time for a particle of type x with channel energy ε

The exciton level density ρ, which depends on the number of particles and holes, is replaced by the compound nucleus level density ω. The cross section is then given by

This is the so called Weisskopf-Ewing formula [53].

Again the target nucleus and the compound nucleus are approximated by the Fermi gas model. In the case of an equidistant single Fermion level density, its level density is

The model, so far, can only produce isotropic emission. However, this is not in agreement for low energy reactions. So the main assumption of compound nuclear theory, that the composite nucleus lives rather long and has “forgotten” the way it was produced, breaks down. However, it seems to be still valid for each partial wave, and for each wave angular momentum and parity are conserved [55]. This implies that the exit partial wave is the same as the initial one, i. e. σinv=Σlσl.

It was already shown in 1937 by Bethe [56] that thermal energy and rotational energy are almost independent of each other. Their functional dependencies can therefore be factorized. We now follow the review by Ericson [33]. The level density for a Fermion system with excitation U and spin J, ignoring the magnetic quantum numen, is

It should be possible to associate a rotational energy and a moment of inertia Θ with the nuclear spin distribution. This comes via σ, the cut off parameter of the spin distribution, which is determined by

Here t is the thermodynamic temperature U=at2−t and Θ the moment of inertia.

There are several spin dependent level density formulae. Sarantites and Pate [57] derived from an approach by Lang [58] a formula with only a very small number of assumptions:

where we have inserted

with Θ the moment of inertia. Since J∝p, the rotational energy depends on E.

Another formula is for the density of states is

with U=E−εx−Bx. The factor 1/2 indicates the equality distribution of the parity. The factor

is the spin distribution. It is the level density of which the TALYS code makes use of [59]. The thermodynamic temperature can be associated with the nuclear temperature T

and U=aT2−3T. Replacing t by T we have the so called constant temperature model. This formula is used within the TALYS cade. For the spin dependent level density, we have then

Here one sees clearly that the portion of the energy Erot is not available for internal energy and thus not for particle evaporation.

The exciton model, as applied here to calculate particle spectra, depends on the Weisskopf-Ewing model Eq. (18). But TALYS applies the constant temperature model. For the PHITS calculations a code treating multiple evaporation is available, the GEM code by Furihata [61]. It makes also use of the Weisskopf-Ewing model [53] but it includes pairing effects in the exponential of the level density. This model, with a Bertini cascade [44] before evaporation, yielded nice results for proton induced reactions [61].

2.5 Moving source model

Starting from the Weisskopf-Ewing evaporation model Eqs. (17)–(19), a handsome equation for the evaporation cross section was derived by Blatt and Weisskopf [62]. Within this model the emission rate is

with T the temperature of the compound nucleus. Here spin and angular momentum dependencies are ignored. The cross section is then given by

In reactions induced by complex particles and not too low energies there are more sources of emission. One is emission from an excited projectile like particle, named fast, and one from an intermediate projectile plus target system. So one is dealing with typically three sources, each of the described by Eq. (30). These cross sections are given in the corresponding centimeter systems. Transformation into the laboratory system yields

In this equation one has three parameters for each source to be fitted to experimental data: an overall normalization constant N0, the source velocity vi and the source temperature Ti. Angle integration gives

3.1 Secondary particle spectra

We now compare experimental proton spectra from α induced reactions with those, calculated within the two non-equilibrium models discussed above. Contrary to the discussion in section 2.2 we do not vary model parameters in order to achieve best possible reproduction of data. However, we use standard values for model parameters to test the predictive power of the models., except the initial exciton number in the case of the exciton number. Here we restrict ourselves to only angle-integrated spectra. Unfortunately the data body for α-particle induced reactions is meager compared to proton induced reactions. We select data from ≈50 to 720 MeV. Not in all cases angle-integrated cross sections were published. In these cases the angle integration was performed by us. The applied methods are discussed in the text. Angle integrated cross sections were calculated within the intranuclear cascade plus evaporation model incorporated in the code PHITS [6].

Exciton model calculations were performed with a code written by the author. In these calculations, transition rates from [34] have been used. Furthermore up to two particle emission during the pre-equilibrium phase has been considered. For some cases the different contributions to the cross section are shown separately. These calculations are indicated by the initial exciton number n0. They include always evaporation of one particle in the equilibrium state applying the compound nucleus model as described above. We ignore angular momentum effects such as angular momentum conservation as well as parity conservation. For the level density parameter the standard value a = A/8 MeV was used, δ=0 and κ=5/4 were assumed.

In the previous section we have stressed the point that multiple particle emission is important at low emission energies. This part is often suppressed in the data, either by a large Coulomb barrier, which exists in heavy target nuclei, or by deficiencies in the experiments, i. e. thick ΔE counters. In Figure 11 we show a case, where the low energy part is visible. The pre-equilibrium part nicely describes the spectrum for energies above 20 MeV. The cross section for the range below 20 MeV is accounted for by evaporation of one particle, secondary particle evaporation and even third particle evaporation.

Figure 11:

Angle integrated cross sections for the indicated reaction. The experimental data [63] are shown as full dots, the pre-equilibrium part as dashed curve, the first equilibrium proton as thin red dotted curve, the proton following the first neutron as blue dashed-dotted curve, the proton following the first proton as green long dashed-dotted curve. The sum of third chance proton emission is shown as thin solid curve (magenta) and the sum of all contributions as thick solid curve.

Chevarier et al. [41] studied proton emission from nuclear reactions at 54.8 MeV for target nuclei, spanning the range from vanadium to gold. They performed Griffin plot analyses, as discussed above, and derived slope parameters s, which showed dependencies of the odd-even structure of the residual nuclei. Since the exciton numbers have to be integer numbers, the initial exciton numbers derived from such analyses, have to be increased or decreased. They used a variation of the exciton model, i. e. the hybrid model [65], and compared the data with calculations with n0=4 and n0=5.

The spectrum for the ⁵⁶Fe(α, p) reaction is shown in Figure 12. The initial exciton number was chosen to be n0=4p+1h. The exciton model calculation agrees nicely with the data with respect to absolute height as well as to spectral form, except for energies below 10 MeV. It is this range which is dominated by multi-particle evaporation not included in the calculation. This finding is in stark contrast to the result of Chevarier et al. [41] who found best agreement with data for n0 = 4 while n0 = 5 underestimates the pre-equilibrium yield. This may be due to the application of a different model. The PHITS calculation overestimates the data by a factor of two in the energy range 10–20 MeV.

Figure 12:

Angle integrated spectrum of the indicated reaction. The data (full circles) are from Ref. [41], the PHITS calculation is shown as histogram (red), the exciton model calculation indicated by the initial particle and hole numbers as solid curve.

In the following figure more data from Chevarier et al. are compared with model calculations (see Figure 13). The exciton model calculation reproduces the shape and height of the experimental spectrum. Although in the previous cases the PHITS calculation tends to overestimate the experimental data, here it underestimate the data. The next heavier target studied is ⁶⁶Zn. Data and model calculations were already shown in Figure 6. The TALYS calculation employing n0=5p+1h underestimate the pre-equilibrium part of the spectrum similar as the PHITS calculation does. The exciton model with n0=4p+1h accounts well for this part.

Figure 13:

Comparison between data and calculations for the indicated target nucleus. The data are from Ref. [41]. the exciton model with the initial exciton number five is shown as solid curve, the PHITS calculation as histogram (red) (red).

We now turn to a lighter target nucleus and almost double bombarding energy. In Figure 14 we compare the data from Ref. [36] The contributions of the first and seconds chance pre-equilibrium emission are shown separately. The calculations agrees very well with the data, while the PHITS calculation vastly overestimates the experimental result.

Figure 14:

Comparison between experiment (data from [36], [37]) and calculations for ²⁵Mg. The contributions of the first (black dashed curve) and seconds chance (green and blue long dashed line curves) pre-equilibrium emission are separately shown. The sum of these is shown as thick solid curve. Also shown is the PHITS calculation (red histogram).

We then discuss the case of the lightest target nucleus in this study (see Figure 15). The data and the best exciton model case have already been shown in Figure 8.

Figure 15:

Same as Figure 8. Also shown is the PHITS calculation (histogram (red)).

Again PHITS overestimates the experimental data. The TALYS calculation underestimates the experimental yield especially in the high energy range. The calculation with the present exciton model code seems to reproduce best the data. However, in this case parameters were adjusted, as was discussed above.

We now proceed to higher beam energies. The next case are data from Wu et al. [64] at a beam energy of 140 MeV. Now the exciton model clearly requires n0 = 6 (5p+1h). The comparison is shown in Figure 16. While the exciton model accounts for the data the PHITS calculation is to large.

Figure 16:

Same as Figure 12, but wit data from Ref. [64].

Further data for nickel exist at an even higher energy. But while Wu et al. [64] report angle-integrated cross sections, this is not the case for the work at the higher energy of 172 MeV. The authors [38] report a quite complete set of spectra from forward to backward angles for a nickel target, but measurements at only three angles in the case of a gold target. In order to extract angle-integrated cross sections in case of the nickel target we have fitted three moving sources to the angle-dependent cross sections. These are Maxwell-Boltzmann distributions in moving frames.

In Figure 17 we show data and fits for three selected angles: a forward angle, an intermediate angle and an almost backward angle. The fits are quite good, although not perfect. The choice of three sources is reasonable because of different reaction mechanisms. There might be break-up of the projectile particle, which will manifest itself in an enhancement of the cross section at forward angles at energies corresponding to the beam velocity. There is indeed such a component visible at 20⁰around 43 MeV which corresponds to the bombarding energy per nucleon. Then there should be an intermediate source from the pre-compound phase and finally evaporation from the compound nucleus. The assumed form of the fit function can be analytically integrated. The resulting parts to the angle-integrated cross sections are shown in Figure 18.

Figure 17:

Energy spectra for some selected angles. The data (full circles) are from [38], the fits with three moving sources are shown as solid line curves.

Figure 18:

The contributions of the three moving sources to the angle-integrated cross section are shown as curves as indicated in the figure (see text).

The three components show distinct spectral shapes: a bell shaped form for the fast component, intermediate a smooth decreasing shape with increasing proton energy and the compound part steeply decreasing with an exponential shape.

We now add the three components to yield the total angle-integrated cross section. This is shown in Figure 19. There the so derived data are compared with model calculations. Again an initial exciton number 5p+1h accounts very well for the data. It looks almost like a fit. The same is true for the PHITS calculation. Deviations are for both models only in the compound nucleus area. Also shown is a calculation with a different exciton model formulation, performed with the TALYS code [7], [8], [9]. While it nicely accounts for the evaporation dominated part of the spectrum, it completely underestimates the range from 20 to 80 MeV.

Figure 19:

Comparison between data (full dots) from the moving source fits performed here to spectra from [38] and model calculations. The exciton model (black dashed curve) employed n0=5p+1h, the PHITS calculation is shown as red histogram (red). In addition a TALYS calculation for E_α = 180 MeV is also shown (blue dotted curve).

We then proceed to the data from the gold target. There are data given for only three spectra: 20, 55 and 77⁰. In addition the lowest registered proton energy is 20 MeV, thus the compound part is missing. The data are, therefore, sensitive to only two sources. The fits are not so well as in the case of nickel. The so derived angle-integrated cross sections are shown in Figure 20. Also shown are the results from the theoretical calculations of the exciton model and the PHITS code. The former shows a quite good agreement with the data while the PHITS code overestimates data for energies below ≈100 MeV. The data cover the range from almost the Coulomb barrier on.

Figure 20:

Same as Figure 19, but a gold target. Angle integrated data are from moving source fits to spectra from Ref. [38].

We now proceed to higher beam energies. Unfortunately the data body above 172 MeV is meager. We have, to the best of our knowledge, only further data at 720 MeV α energy. In Ref. [39] energy spectra ranging from protons to ⁴He were reported for three angles: 60, 90 and 100⁰. In a second publication [40] the group added cross sections for 30 and 150⁰.

In order to extract differential energy spectra, angle-integration has to be performed. We used the same method as before: fitting three sources. The result of such a procedure in the case of the ²⁷Al is shown in Figure 21, together with the combined set of data. The fit parameters have rather large uncertainties. This might be an indication that the model function in not the best one to account for the underlying physics. The fits overestimate the data at the highest energy for all angles. We can therefore conclude that the angular integral will be too large. However, the discrepancy between data and fits are very large on a semi-logarithmic plot.

Figure 21:

Same as Figure 17, but data for the indicated reaction are from [39], [40].

Another method to obtain the differential cross sections from the double differential ones is to use angular distributions for constant proton energies. We then fitted simple formulae on a logarithmic scale to the distributions. From these fits we obtain differential cross sections by analytic integration of the fitted function. Results are also shown in Figure 22 for the ²⁷Al target, and in Figure 23 for the ¹⁸¹Ta target, as full dots. This method does not depend on a model for the underlying reaction mechanism. However, the results of both procedures are not too different. We then compare the differential energy spectra with model calculations.

Figure 22:

Differential cross sections for the indicated reaction, obtained by fitting angular distributions to data from [39], [40], are shown as full dots. The angle-integral from the moving source fit is shown as blue solid curve. The exciton model calculations with different initial exciton numbers are shown as dashed and dotted curves. The PHITS calculation is shown as red histogram (red).

Figure 23:

Angle integrated spectrum of the indicated reaction. The data marked with full circles are obtained from fitting angular distributions, those with full squares are from extrapolations of these fits. The proton being emitted as first particle in the exciton model is labelled as px, the proton following that proton is labelled as pp (green) and the proton following a first neutron as np (blue). Also shown is the sum of these contributions as solid curve. The PHITS calculation is shown as histogram (red).

The PHITS calculations agree to some extent with the data at small energies. At higher energies the calculations show a bell like structure which is not visible in the data. It can be that this structure is due to projectile break-up. This process may not be in the data since the most forward angle is 30⁰, while break-up yield is much more focussed to more forward angles. At the rather high beam energy the fragments may move into a cone with smaller opening angle.

For aluminum we present two exciton model calculations, one with n0=5p+1h, which was the best choice for most of the data at lower energy, and n0=8p+4h as is suggested by the Griffin plot. The latter is a better representation for the data, although not of the same quality as in the lower energy regime. The situation is almost the same in the case of the tantalum target.

In Figure 24, we show angular distributions for some selected proton energies. The slopes become steeper with increasing proton energies. The steep slope at 155 MeV can be understood as emission from a projectile like system. The beam velocity is 160 MeV/nucleon.

Figure 24:

Angular distribution of emitted protons at the indicated energies. The data are from Refs. [39], [40], the solid curves are fits. The dotted curves indicate the error band for 155⁰. For the other angles the uncertainty is in the range of the curve width.

The flat distribution at 12.5 MeV can be interpreted as emission from a compound nucleus. The exponential behavior for the energy in between would have its origin, if there is only one dominant emitting source. In order to examine this question further we study the energy range between the two extremes. This range is dictated by the availability of experimental data. We then fitted exponentials to the experimental angular distributions

The fitted parameters α and −β are shown in Figure 25. Both parameters show the same dependence: a steep decrease in the compound area and a flat distribution in the beam energy per nucleon region. In between there is no distinct source visible but a smooth interpolation between the two extremes. The systems studied are compiled in Table 1.

Figure 25:

The fitted parameters to the angular distributions Eq. (33) as function of the proton energy (see text).

3.2 Excitation functions

The cross sections for a given reaction increases from threshold due to a larger phase space with increasing beam energy. When the next channel opens, it starts to decrease. This leads to an almost bell shape at small energies. This part is dominated by compound nucleus emission. A long tail at higher energies is dominated by non-equilibrium processes, which will be studied here. The sum of cross sections summed over all open channels is the total reaction cross section.

In the following we compare model calculations with data. There is a wealth of data at rather small energies, which were measured at low energy accelerators for application purposes. The bulk of data is below 25 MeV beam energy, and therefore dominated by evaporation. We are interested in the model’s ability to predict the non-equilibrium cross section. In order to confront the model results with experiments over a broad energy range, we have selected four target nuclei spanning the periodic table: ²⁷Al, ¹⁶⁵Ho, ¹⁸¹Ta, and ¹⁹⁷Au. In order to avoid complications due to Coulomb barrier effects or to complex particle condensation, we restrict the study mainly to neutron emission.

3.2.1 Reactions on ²⁷Al

The reactions on ²⁷Al are shown in Figure 26 and Figure 27.

Figure 26:

Cross sections for ²⁷Al(α, 1n) and ²⁷Al(α, 2n) reactions as indicated in the figures. Experiments are indicated : from [66] (full dots), [67] (red squares), [68] (blue crosses), [69] (magenta triangle up), [70] (green triangle down). TALYS calculations are shown as solid curves, PHITS calculations as broken curves.

Figure 27:

Cross sections for ²⁷Al(α, 2p) and ²⁷Al(α, 3p) reactions as indicated in the figures. Experiments are are from [71] (full dots) and [72] (open red dots). TALYS calculations are shown as solid curves, PHITS calculations as broken curves.

The upper frame of Figure 26 compares data and calculations for the ²⁷Al(α,1n) reaction. For the 1n reaction the data range up to 30 MeV. The two calculations are for Eα>10 MeV slightly below the experimental results. In the steep rising part with Eα<10 MeV there is a remarkable agreement between data and calculations. The lower frame compares data and calculations for the ²⁷Al(α,2n) reaction. We found only three data points in the literature for this case. It seems therefore useless to discuss agreement or disagreement between data and theories. More data exist for proton emission in ²⁷Al(α,2p) reaction (upper frame) and ²⁷Al(α,3p) reaction (lower frame in Figure 27). The PHITS calculation reproduces nicely the excitation function for the ²⁷Al(α,2p) reaction, while the TALYS calculation vastly underestimates the experimental results. This is also true for the ²⁷Al(α,3p) reaction, where the PHITS calculation is favorably, although the shape is different than the one of the data.

3.2.2 Reactions on ¹⁶⁵Ho

The next isotope studied is ¹⁶⁵Ho. Here the data body is much more rich than in the case of aluminum. The data for the ¹⁶⁵Ho(α,xn)169−xTm for x=1−4 are shown in Figure 28.

Figure 28:

Same as Figure 26, but for ¹⁶⁵Ho(α,xn) with x=1−4. Exp. data are from [73] (full dots), [74] (full red squares), [75] (blue triangle u), [76] (magenta closed diamond), [77] (green triangle down), [78] (brown star), [79] (purple hourglass), [80] (royal crossed square) and [81] (aqua blue open stars).

The data reach up to E_α ≈ 70 MeV for x=1 and up to 110 MeV for the larger x-values. The data were reproduced by the calculations with respect to shape and absolute height. The PHITS calculation does slightly better than the TALYS calculation. Unfortunately we could not find data with x≥5, and hence no further comparison could be made.

3.2.3 Reactions on ¹⁸¹Ta

The next heavy target nucleus for which a large data body exists, spanning a large energy range is ¹⁸¹Ta. The excitation functions for neutron emission are shown in Figures 29–32 for x=1−8.

Figure 29:

Same as Figure 28, but for ¹⁸¹Ta(α,xn)(185−x)Re with x=1−2. The experimental cross sections are from [82] (full dots). [83] (red full squares), [84] (blue full triangle up), [86] (magenta full diamond), [87] (green full triangle down), [85] (green hourglass), [88], (brown star), [90] (purple crossed square),[89] (cyan crossed triangle up), [91], (black crossed dot).

Figure 30:

Same as Figure 29, but for x=3−48.

Figure 31:

Same as Figure 29, but for x=5−6. More data are from [92] (magenta crossed diamond).

Figure 32:

Same as Figure 29, but for x=7−8.

Both calculations reproduce the data for x odd and x=1−3. For the even cases the reproductions are poor with TALYS, worse than PHITS. This is surprising since TALYS applies a pairing energy, both in the pre-equilibrium as in the equilibrium phase as correction to the level density. The INC code PHITS has the pairing energy only in the evaporation part. For x=6 both models predict cross sections at energies, which are below the lowest energy in the experiments. Since both model calculations behave almost identical in this range, the question arises: do the models have a deficiency or the data? From the similarity of the model calculations one might tend to attribute this discrepancy to the data. However, in the calculations by Hermes et al. [91] such an effect was not observed. Ernst et al. [93] claimed that such an effect will occur, when a shell is crossed. However, the present models make use of experimental masses [43], [94], which of course contain shell effects and in addition, there is no shell close by. Hermes et al. have used different separation energies and found the best agreement between data and calculations for the energies taken from the tables by Garvey et al. [95].

To compare the masses from these inputs we have deduced separation energies from the mass tables of Audi and Wapstra. This comparison is shown in Figure 33. There is no serious deviation between the two mass tables.

Figure 33:

Separation energies for neutrons in the reaction chain Ta(α, xn)^185−xRe from Audi and Wapstra [43], [94] and Garvey et al. [95].

Another possibility is a fault in the energy determination of the experiments. The impinging α-particle energy is degraded within a foil stack. The energies are derived by making use of the stopping power tables by Williamson et al. [96]. We have compared the result from this table with those from Ziegler et al. [97] and Hubert et al. [98]. No deviations for the ranges of α-particles in tantalum could be found.

The next possibility for the differences between experiment and theoretical calculations may in the underlying models. Since both models behave similarly, there may be one common reason for this behavior.

In order to get a quantitative view we have fitted a Gaussian plus a constant to the theoretical excitation functions obtained with the PHITS code and the TALYS code and to the experimentally gained cross sections.

The fitted values of the maxima are shown in Figure 34. Obviously there seems to be a proportionality between the maxima and hence a shift in the energy scales represented by them. We have fitted a straight line to these relations. It is also shown in the figure together with the standard deviation. The result is

(34)Eexp=α∗EPHITS.

Figure 34:

The maxima in the experimental excitation functions [91] as a function of the maxima from PHITS calculations for the (α,xn) reactions on tantalum (full dots). A linear fit is shown as solid line together with the error band.

The slope parameter α is given in Table 2. The evaporation part in the PHITS code is the GEM code from Furihata [61]. Although it allows for a lot of open channels, it is a pure Weisskopf-Ewing approach, neglecting angular momentum effects. We can ,therefore, expect that the calculation gives cross section for small energies because Erot is not subtracted (see Eq. (27)).

Table 2:

Deduced rotational energies Qexp=Erotexp/EPHITS, Qr=Erotr/EPHITS and momenta of inertia.

Target	Α	Qexp	Qr	Θexp/Θr
181Ta	1.200 ± 0.011	0.200 ± 0.011	0.0552	0.27 ± 0.02

From this we find for the rotational energy.

(35)Erotexp=Eexp−EPHITS

(36)=(α−1)EPHITS.

The so derived rotational energy is also given in the table.

We may compare this value with the one of a rigid rotator. This is given by Eq. (22) and

(37)Θr=25mAR2

with m the nucleon mass. For the sake of simplicity we have assumed J to depend on the maximal impact parameter. The so derived rotational energy for a rigid rotator is also given in Table 2. We now get the ratio of the momenta of inertia from

(38)ΘexpΘr=ErotrErotexp.

These ratios are also given in Table 2. The ratio is always smaller than one. This is in agreement with findings by Bohr and Mottelson [100].

For the TALYS calculations we find almost the same slope parameter. However, TALYS deals with the constant temperature model. The energy shift between experiment and TALYS calculation remains therefore an unsolved problem. One anonymous referee conformed the present finding. However, he found out that this shift is not there for calculations with TALYS prior to version 1.6. It can therefore be speculated that from version 1.6 onwards the constant temperature model is not more correctly installed or addressed in the code.

In order to study the effect of rotational energy on the excitation functions, we compare the so modified calculations for the two models with experimental data. This is done for the tantalum target in Figure 35. Since all data at higher energies are from the work of Hermes et al. [91], we restrict the comparison to only these data. The procedure does not only bring the maxima of the theoretical excitation functions in agreement with those of the experiment, but improves also the agreement for the tails. The PHITS calculations predict smaller tails than the TALYS calculation. On the other hand TALYS gives too large cross sections in the maxima for the larger x values. In general, the shift vastly improved the quality of agreement.

Figure 35:

Upper frame: Excitation functions for neutron removal from Ta with x = even. The experimental data (curves with points) are from [91], the solid curves are the TALYS calculations and the broken line curves the PHITS calculations with shifted energy scale. Lower frame: Same as upper frame, bur for x = odd.

3.2.4 The Reactions ¹⁹⁷Au(α, xn)^(201-x)Tl

We then proceed to the next heavier target nucleus: ¹⁹⁷Au. The data cover the range up to 160 MeV with a wealth of data in the range where compound emission dominates.

Data comparison is made in Figures 36–38 for cross sections for reactions ¹⁹⁷Au(α, xn)^(201−x)Tl and x=1−8. The data agree with each other, except for one set for the ¹⁹⁷Au(α,1n) reaction. It is this reaction where both model calculations agree nicely with the experiments up to E_α = 70 MeV. For larger values of x the calculations underestimate the experimental yield with increasing deviation with increasing x. We perform the same procedure as in the case of tantalum. A linear fit yields

(39)Eexp=(1.08±0.01)*EPHITS.

Figure 36:

Similar as Fig. Same as Figure 29, but for ¹⁹⁷Au(α, xn)^(201−x)Tl and x=1−4. The experimental data are from [101] (open triangle up), [102] (red full diamond), [103] (blue full triangle down), [85] (magenta star), [104] (brown cross), [105] (purple crossed dot), [106] (green crossed triangle up), [107] (dark blue full triangle up), [108] (brown crossed triangle down), [109] (cyan crossed square), [110] (black hourglass), [111] (red open diamond), [112] (η dark blue) [113] (magenta open diamond), [114] (brown full square),[115] (green x), [117] (purple full triangle down), [118] (black full dot).

Figure 37:

Same as Figure 36, but for x=5−6.

Figure 38:

Same as Figure 36, but for x=7−8.

This value is smaller than in the case of tantalum. It yields a moment of inertia Θexp/Θr≈0.7, a value, almost a standard in the literature [60].