UCRL-ID-144510
July 11, 2001
N Division
Physics and Advanced Technologies Directorate
Lawrence Livermore National Laboratory
Livermore, CA 94551
Table of Contents
2. The Maxwellian-averaged thermonuclear reaction rate _{}
3. Maxwellian-averaged mean energies and cosines for two-body final state reactions
5. Maxwellian-averaged two-body final state reaction particle spectra
6. Digression on two-body final state spectrum moments
7. Data, algorithmic bases and results of _{} calculations
8. Algorithmic bases and results for _{}, _{} and _{}
9. Algorithmic basis of the final state spectral probability distributions
10. The Thermonuclear Data File (TDF) generator
11. Structure and content of the TDF
12. The TDFLIB subroutine library
Appendix A. Example TDF format
Appendix B. Erratum in Williams' original paper
Appendix C. Evaluation of a d function integral over finite limits
The rate of thermonuclear reactions in hot plasmas as a function of local plasma temperature determines the way in which thermonuclear ignition and burning proceeds in the plasma. The conventional model approach to calculating these rates is to assume that the reacting nuclei in the plasma are in Maxwellian equilibrium at some well-defined plasma temperature, over which the statistical average of the reaction rate quantity sv is calculated, where s is the cross-section for the reaction to proceed at the relative velocity v between the reacting particles. This approach is well-understood and is the basis for much nuclear fusion and astrophysical nuclear reaction rate data. (See, e.g., Chiu 1968, Clayton 1968)
The Thermonuclear Data File (TDF) system developed at the Lawrence Livermore National Laboratory (Warshaw 1991), which is the topic of this report, contains data on the Maxwellian-averaged thermonuclear reaction rates for various light nuclear reactions and the correspondingly Maxwellian-averaged energy spectra of the particles in the final state of those reactions as well. This spectral information closely models the output particle and energy distributions in a burning plasma, and therefore leads to more accurate computational treatments of thermonuclear burn, output particle energy deposition and diagnostics, in various contexts. In this report we review and derive the theoretical basis for calculating Maxwellian-averaged thermonuclear reaction rates, mean particle energies, and output particle spectral energy distributions for these reactions in the TDF system. The treatment of the kinematics is non-relativistic. The current version of the TDF system provides exit particle energy spectrum distributions for two-body final state reactions only. In a future report we will discuss and describe how output particle energy spectra for three- and four-body final states can be developed for the TDF system.
We also include in this report a description of the algorithmic implementation of the TDF system, which is in two parts: the code which generates the data file itself from the nuclear reaction cross-sections, and the lookup utilities (bundled in a Fortran 77 subroutine library) which are used to retrieve the appropriate thermo-nuclear data from the file on demand. Unique approaches to generating, retrieving and looking up the data that were developed will be described in some detail. These developments resulted in an unusually compact data file which can be quickly generated, and from which an enormous amount of thermonuclear reaction rate, output particle energy and particle energy spectrum data can be rapidly retrieved. It is this compactness, speed and relative platform independence that distinguishes the LLNL TDF system from others.
2. The Maxwellian-averaged thermonuclear reaction rate _{}
Consider a plasma containing a mixture of two nuclear species which react to produce other particles and nuclei. We quantify the rate of thermonuclear reaction by first considering a "local" microscopic model in which one of the nuclear species is at rest, and the other is traveling toward the stationary nuclei in a beam moving at velocity v, as shown schematically in Figure 1.
Figure 1. Particle-beam particle-target geometry
With the geometric notions shown in the figure, the number of particles in the beam segment is
_{},
where A is the area of the beam cross-section (and therefore the area of the particle interaction volume). The total area presented by the target particles within the particle interaction volume to the particles in the beam is
_{}.
Thus the fraction of the particle beam that is intercepted by the target particles is
_{},
and the total number of beam particle - target particle interactions is therefore
_{}.
Since the interaction volume is ADx and the interaction time is Dt, we then have that the number of interactions per unit volume per unit time is
_{}.
This quantity is considered to be the instantaneous local reaction rate density for the binary plasma in question when it is understood that v is the scalar relative velocity
_{}
between the nuclei of the different species, the binary interaction cross-section s is strictly a function of v_{R}, and the particle densities n_{1} and n_{2} are understood to be locally constant and uniform.
In a plasma the velocities v_{1} and v_{2} of two interacting component species of particle are distributed in magnitude and direction, so that the number of particles per unit volume of species type 1 having a mean velocity v_{1} within the differential cell d^{3}v_{1} = dv_{1x}dv_{1y}dv_{1z} in velocity space is
_{},
and the number of particles per unit volume of species type 2, similarly, is
_{},
where f(v_{i}) is symbolic for the normalized distribution function for particle species i in the plasma. Then the total reaction rate density for the interacting species in the plasma is the integral of the local reaction rate density over all possible particle pair combinations and vector relative velocities, as follows:
_{}
Note that this is an integral over a joint distribution f(v_{1})f(v_{2}). Since the functions f(v_{i}) are assumed normalized, that is,
_{}
this quantity is also the statistical average of n_{1}n_{2}sv over the joint distribution, and may therefore be written as
_{}.
In the case where both reacting particles are identical (as in the d(d,n)^{3}He reaction, for example), the RRD needs to be calculated differently, because during the double integration over dn_{1}dn_{2} the same particle population is integrated over twice. This means that each reacting particle pair is counted twice, and the RRD in this case is double what it should be. The correct formulation of the RRD in the general case is therefore as follows:
_{}
where i and j denote the reacting particle species, and d_{ij} is the Kronecker delta, which is equal to 1 when i = j and 0 otherwise. (This multiplying factor issue is discussed more thoroughly by Clayton (1968)).
Because the particle densities n_{i} and n_{j} are particular to a given plasma configuration and so may themselves be distributed in time and space, we have formulated the TDF system to provide only the "bare" statistical average _{}. It therefore does not include the pair count multiplier n_{i}n_{j}/(1 + d_{ij}). In this report we treat only the statistical average _{} and ignore the multiplier, as well.
In our modeling we assume the reacting particle species each have an isotropic Maxwellian distribution at a common kinetic temperature q = kT. Then the particle distribution functions for each species also depend on particle mass and plasma temperature, and have the well-known form
_{}
where m_{i} is the particle mass, v_{i}^{2} is the square of the magnitude of the vector velocity v_{i}, and i = 1 or 2 identifies the particle species. These functions are normalized. Using these functions, the statistically averaged reaction rate density for the two species in thermal equilibrium at a plasma temperature q = kT is given by the following expression:
_{}.
This is the formal development of the Maxwellian-averaged thermonuclear reaction rate _{} in integral form. It is clear that _{} is a function of q. (Brysk (1973) worked out _{} for the case where the two reacting components each have a Maxwellian distribution but at different temperatures.)
This integral can be put in a simpler form by changing the velocity variables from v_{1} and v_{2} to the relative velocity
_{}
and the center-of-mass velocity
_{}.
The total kinetic energy then transforms as
_{}
where M = m_{1} + m_{2} is the total mass of the system andm = m_{1}m_{2}/M is its reduced mass. The term _{} is denoted by T_{cm} and is the kinetic energy of the center of mass of the particle system. The term _{} is denoted by T_{R} and is the kinetic energy of relative motion; it is equal to the sum of the kinetic energies of the two reacting particles in their center-of-mass (cm) frame. It is also proportional to the incident laboratory frame bombarding energy of one particle when the other (the target) is at rest, since the relative velocity v_{R} is the same in all frames of reference.
The Jacobian of this transformation is unity, so the differential elements transform directly as d^{3}v_{1}d^{3}v_{2} = d^{3}v_{R}d^{3}v_{cm}. On applying this transformation to the integral for _{}, we find that this integral separates into two parts, as follows:
_{}.
_{}
_{and}
_{}_{.}
_{On the assumption of particle direction isotropy within the plasma, the integrations over the polar angles can be carried out immediately, to yield the result that}
_{}.
This may be further simplified by making one more variable substitution in the form _{}, and carrying out the integration over dv_{cm}. Then we obtain the single integral for _{} as
_{}.
The integrand _{}is sometimes called the Gamow distribution function.
Note that the above reaction rate results hold for any initial pair of reacting nuclei whose total reaction cross-section is specified, and are independent of how many particles or additional channels result in the final state. Thus, if one of the final state particles breaks up into additional particles, or any exit nuclei are left in excited states, the total reaction rate is still the same. On the other hand, if the cross-section chosen is specific to a particular final state channel out of many channels for the same two initial reacting particles, then the reaction rate associated with this channel will necessarily be different.
We also note in passing that Slaughter (1983, 1985, 1986, 1989) describes Monte Carlo calculations of reaction rates and output particle spectra for other types of particle distributions in which the particle motions are not isotropic, and therefore non-Maxwellian. This occurs, for instance, when a directed beam of one species of particle with a narrow energy spread is incident on a quiescent plasma containing other species of particles.
3. Maxwellian-averaged mean energies and cosines for two-body final state reactions
The intrinsic properties of the thermonuclear reaction components during the plasma conditions are also of interest. These include the mean kinetic energies, directions and energy spectra of the incident and exit particles, which can be used to estimate, e.g., energy deposition or diagnostic particle outputs as thermonuclear burn proceeds (Slaughter 1989). In this section we develop the Maxwellian-averaged mean energies and direction cosines for the reacting and exit particles of two-body final state reactions. In the next section we develop model estimates for the mean energies of the exit particles of three-body final state reactions, specializing to the t(t,2n)a reaction. Corresponding exit particle energy spectra for two-body final state reactions are discussed in the section after that. Spectra for three-body final state reactions are more complex and will be treated in a separate report.
A two-body final state nuclear reaction taking place in a plasma is depicted schematically in Figure 2, where the two incident and two exit particles are labeled with the indices 1 through 4, and the corresponding masses, velocities and angles are appropriately labeled. We observe that the incident and exit vector pairs (v_{1},v_{2}) and (v_{3},v_{4}) each define planes containing the common vector v_{cm}, which is the velocity of the center of mass of the interacting and exit particle pairs. We assume in our development that the nuclear reaction is such that the orientations of these two planes are not correlated by the interaction, and are therefore arbitrary. The incident vectors are therefore not coplanar with the exit vectors, in general.
Figure 2. Two-body final state reaction kinematics
The kinetic energies of each particle are labeled in the plasma (i.e. "lab") frame as _{}, while in the "center of mass" (cm) frame they are written as _{}. The kinetic energies of the motion of the center of mass and the relative motion of the incident particles are related to these particle energies as follows:
_{}
and
_{},
where
_{}.
The kinetic energy associated with the relative motion of the exit particles is
_{}
where
_{}
and
_{}.
We assume that the quantities M (= m_{1} + m_{2} = m_{3} + m_{4}), T_{cm}and v_{cm} are the same before and after the interaction, within a very small error (of the order of Q/Mc^{2} where M is the total reaction particle mass), and that the relative velocities v_{R} and v_{R}^{new} are each independent of frame of reference.
We calculate Maxwellian-averaged mean properties of these reaction quantities by first expressing all energies in terms of T_{R} and T_{cm}, since, as we will show, this facilitates developing the separable expressions for their statistical averaging. This transformation, however, has a much deeper significance apart from this convenient consequence, because it effectively separates the plasma-related phenomenology and the nuclear interaction. This is because in any interaction in the plasma, the center of mass motion v_{cm} and the relative particle motion v_{R} are independent of one another, i.e. uncorrelated. Thus the center of mass kinetics are determined solely by the plasma properties, and none of the nuclear interaction properties depend on any of the center of mass quantities — including the direction of v_{cm}. And, s is a function of T_{R} and not of T_{cm}. This is why — as we will show later — the statistical average of T_{cm} depends on q, that for T_{R} depends on both q and s, and the angular correlations between relative and center of mass motion quantities average to zero.
We first express the initial particle velocities in terms of v_{R} and v_{cm}, as follows:
_{}
_{}.
From this we see that v_{1}¢ = m_{2}v_{R}/M and v_{2}¢ = - m_{1}v_{R}/M. It is then a fairly straightforward procedure to use the previous definitions of the various kinetic energies to obtain the initial state energies and relations in terms of T_{cm} and T_{R}:
_{}
_{}
and
_{},
where the q_{i}¢ are the angles that v_{1}¢ and v_{2}¢(and therefore v_{R}) make with v_{cm}, as shown in Figure 2. In other words,
_{}
where the caret ( ˆ ) denotes the unit vector.
The final state particle energies and relations are correspondingly rendered, when it is understood that the reaction Q-value (the energy released by the reaction process) is included only in the final particle energy and velocities, and leaves the center of mass quantities T_{cm}, v_{cm} and M unaffected with negligible error. This says that
_{}
is the only essential difference; thus the relations for the final state can be copied directly from those of the initial state by appropriately relabeling the particle indices. We then have that
_{}
_{}
and
_{}.
q_{3}¢ and q_{4}¢ are the angles that v_{3}¢ and v_{4}¢(and therefore _{}) make with v_{cm}, as shown in Figure 2, so that
_{}.
We can now address the statistical average of these reaction energies and cosines over the plasma distribution, and this requires some care. Basically, the only way to obtain the mean value of a reaction quantity in a plasma is to calculate its average over a reaction probability distribution, and the only distribution that qualifies is the reaction rate function that determines _{} – not just the Maxwellian distribution of the plasma. This function is the integrand used to calculate _{}, or
_{}.
which is a joint probability distribution in d^{3}v_{R}d^{3}v_{cm} space. It is important to note that this function identifies, or selects, in a statistical sense, the distribution of that part of the particle population of the plasma that participates in the reaction, and only that part. It is also important to notice that this distribution identification holds for any and all final state configurations that result from the reaction, which is to say that it is independent of whatever values and directions the final state particle velocities do take. As we will show shortly, this concept is the key element that helps resolve anisotropic differential cross-section issues (Perkins, 1990 and 1991) that arise in calculating the Maxwellian-averaged reaction quantities for the final state.
Because F is not normalized to 1 (its integral is equal to _{}), the statistical average of any quantity, say f, over this reaction probability distribution, which we denote by _{}, is calculated by the classic ratio of two integrals,
_{}
where the subscript "reac" on the angular brackets denotes the statistical average over the reaction distribution F, and the angular brackets without a subscript continue to denote the statistical average over the Maxwellian distribution. Thus the denominator in the above is, of course, _{}, while the numerator is (now obviously) the Maxwellian average _{}. We can write this more compactly as
_{}.
Notice that the exponents in the integrands are T_{R}/q and T_{cm}/q, which is the reason for rewriting all kinematic quantities in terms of T_{R} and T_{cm}. This means that if the quantity f is separable into functions of T_{R}, T_{cm} (or v_{R}, v_{cm}) or the angular directions, then individual integrations over these arguments can proceed separately. As it turns out this will indeed be so for our needs.
We developed T_{i} in terms of T_{R}, T_{cm} and the cosines given by and . It is a fairly easy matter to obtain the reaction rate averages of these latter quantities, since their integrals are now separable. The statistical reaction average for T_{R} reduces to a single integral much the same way as the statistical Maxwellian average for sv worked out, and the result is
_{}
Here, too, numerical integration is required, in an equally simple form as for_{}.
Zimmerman (2001) has pointed out that the slope of the _{} curve as a function of q is related to _{}. This is immediately evident by differentiation: one finds, by including the ∂ /∂q operation on the _{} integrand, that
_{}.
Thus _{} can be expressed in terms of _{}. The reaction statistical average for T_{cm} easily reduces to its Maxwellian average:
_{}.
Finally, we consider the reaction averages of the particle direction cosines, which contain some subtle issues. To calculate the reaction statistical average of the cosine term _{}, we select the direction of v_{cm} as the z-axis for the polar coordinate system for v_{R}. Then
_{}
and the angular integral part (in the v_{R} system) of the reaction average of this quantity vanishes:
_{}.
Thus
_{}.
This reiterates our earlier observation that v_{R} and v_{cm} are uncorrelated. The reaction statistical average of the cross-terms in the expressions for T_{1} and T_{2} therefore vanish.
Next, we calculate the reaction statistical average of the post-reaction cosine term _{}. However, the angular anisotropy in the emission of final state products needs to be addressed, because the angular distributions ds/dW of final state particles in nuclear reactions are neither uniform nor constant. As it turns out (Zimmerman 1990) this is of no consequence, because the angular distribution is a function of the angle that _{} makes with the initial relative velocity direction _{}, and the direction cosine of _{}, as we showed earlier, has a reaction statistical average of zero. We may put this in another, more graphical way. Suppose one of the final state reaction particles was emitted preferentially at some fixed angle a which _{} makes relative to the original interaction particle direction of approach, _{}. The statistical average would then be taken of the cosine of this direction relative to the direction of _{} (which is the polar axis for the integration). We are therefore in fact calculating the statistical reaction average of the quantity
_{},
and this, too, results in zero. Thus the reaction statistical average of the cross-terms in the expressions for T_{3} and T_{4} also vanish.
The reaction statistical average of the particle kinetic energies in the plasma are thus as follows:
_{}
_{}
_{}.
Energy conservation clearly still holds during reaction averaging:
_{}
We emphasize again that the statistical averages which are labeled by the subscript "reac" are over only those particles that actually participate in the reaction, and not the entire plasma population. The above conservation law is thus so restricted. Perkins and Howerton (1988) obtain the same results but with a priori initial assumptions of average energy conservation and isotropy of reaction angular distributions. Our treatment does not require these assumptions.
4. Exit particle energies in the t(t,2n)a reaction: a foray into three-body final state energy averaging
The Maxwellian-averaged mean energies of the three exit particles in the t(t,2n)a reaction are handled differently, because at any given initial particle total energy, the center-of-mass energies of the exit particles are distributed, rather than sharply clustered at a single value. This is because this reaction proceeds partly by instantaneous three-particle breakup, and partly by a sequence of two-body states. In the latter process ^{5}He is formed by the t(t,n)^{5}He reaction, and then breaks up into a neutron and an alpha while in flight.
In the current version of the TDF system, we estimate the mean energies of the exit particles in the center-of-mass frame for the t(t,2n)a reaction by calculating the expectations of their energies over the final state energy distributions that were experimentally measured at LLNL by Wong, Anderson and McClure (1965). They directed a 500 keV beam of tritons on a tritium target and measured the energy spectrum of the neutrons and alphas emitted perpendicular to the beam line. The resulting center-of-mass frame neutron and alpha energy distributions suggested that the t(t,2n)a reaction proceeded 70% by direct three-body breakup, and 30% by sequential two-body final states, first by the t(t,n)^{5}He reaction, then by breakup of ^{5}He in flight into an alpha and a neutron. The 30% sequential breakup distribution was further estimated to subdivide into two competing processes in which ^{5}He was either left in the ground state (20% of the time) or in a broad first excited state at 2 MeV (10% of the time). (This first excited state has subsequently been assigned a level energy of 4 MeV, see Ajzenberg-Selove (1988)). Each one of these three processes (channels) results in different mean energies for the two neutrons and alpha particle. The mean energies of these particles for all reaction processes are therefore their expectation energies, which are the sum of their mean energies in each process weighted by the relative frequency of each process.
If we make the simplifying assumption that these branching ratios hold for the t(t,2n)a reaction in a warm plasma, then the mean energies of the exit particles can be approximately modeled in a relatively straightforward manner. We first consider the reaction kinetics at zero plasma temperature, for then the center-of-mass frame is at rest and the emitted particle energies take on simple and well-defined distributions. In the simultaneous three-body breakup process, the exit particles each have a phase space spectral energy distribution of the form
_{}
where
_{},
M is the sum of the particle masses, and i = 3, 4 and 5. (This corresponds to the two neutrons and one alpha particle of the t(t,2n)a reaction). The average of the particle energies over this symmetric distribution is obviously
In the sequential two-body final state mode, the first final state is the "usual" two-body case, and the energies of the emitted particles are
_{}
for i = 3 and 4 (corresponding to the neutron and ^{5}He nucleus of the t(t,n)^{5}He reaction). Their energy distributions are very sharp, so _{} for i = 3, 4. The second final state occurs when the heavier particle (say, m_{4}) breaks up in flight into two particles, say m_{5} and m_{6} (corresponding to a neutron and an alpha particle), and releases breakup energy Q_{B} in the process. The exit particles in this stage have a uniform (boxcar) distribution of energies between two limits, which are given (for, e.g., particle 5) by the relations
The average particle energy over this distribution is evidently
_{}.
For the t(t,2n)a reaction, the Q value is 11.33 MeV, and there are two known competing sequential two-body final states corresponding to ^{5}He left in the ground state and in a broad first excited state at 4 ± 1 MeV (Ajzenberg-Selove, 1988). The breakup energy (Q_{B}) for ^{5}He in the ground state is .89 MeV, and for the first excited state it is approximately 4.89 MeV. The mean particle energies (i.e. the average over that particle energy distribution) for each channel of this reaction are calculated from the above expressions, and the results are presented in the following Table I.
channel and fraction (b_{i}) |
_{} |
_{} |
_{} |
E_{x}(^{5}He) |
three-body breakup (70%) |
4.718 |
4.718 |
1.906 |
- |
1^{st} sequential two-body (20%) |
8.694 |
1.066 |
1.589 |
0 |
2^{nd} sequential two-body (10%) |
5.363 |
4.124 |
1.854 |
4 ± 1 |
weighted mean, 3 channels |
5.577 |
3.928 |
1.838 |
- |
Table I. Mean center-of-mass frame energies (MeV) for
t(t,2n)a exit particle distributions in each channel
The last line in the table gives the weighted mean energies of each particle over all reaction channels. These are calculated from their expectation sum, which is
_{},
where the subscript i is the channel label, _{} is the particle energy for each channel averaged over the particle's distribution, and the branching ratio fractions b_{i} for each channel (also given in the table) sum to unity: b_{1} + b_{2 }+ b_{3} = 1.
We estimate the reaction-averaged energies of the three exit particles as a function of plasma temperature by the following model, which is used in the TDF system until a more comprehensive and correct treatment of three-body final state energetics in a warm plasma is developed. The three weighted mean center-of-mass energies in the table add up to the final state total kinetic energy, which is the Q value. This corresponds to a plasma at zero temperature. The partition fraction that each exit particle i has of this total energy is f_{i} = _{}/Q, so for the three particles we have f_{n1} = 0.492 ± .007, f_{n2} = 0.347 ±. 007 and f_{a} = 0.162 ± 001. The uncertainties are due to the ± 1 MeV uncertainty in the energy level of the first excited state of ^{5}He. We now assume in our approximation model that these same fractions f_{i} (of the reaction-averaged final state total energy) also determine the reaction-averaged energy of the exit particles for any plasma temperature (White, 1990). That is, for any given q we calculate _{} for each exit particle i, where
and where subscripts 3, 4 and 5 correspond to n_{1}, n_{2} and a. Then the reaction-averaged mean kinetic energies of each exit particle are assumed to have the values _{} = 0.492 _{}, _{} = 0.347 _{}, and _{} = 0.161 _{} at that plasma temperature.
The procedure is thus to first calculate _{} for a given q, and then use the above partition fractions to estimate the three final state particle _{} values. These fractions are "hard-wired" into the TDF system mean energy lookup routine for the t(t,2n)a reaction. It is expected that these averaged energy estimates — as derived from the partition fractions — will be revised (perhaps minimally) when a correct three-body final state reaction formalism is eventually developed for the TDF system.
The t(t,2n)a reaction partition fractions used in earlier versions of the TDF system differ slightly from the ones given above: they are f_{n1} = 0.473, f_{n2} = 0.364,and f_{a} = 0.163. The reason for this is that a final state consisting of three sequential two-body final state channels was used to calculate the partitioning, and the three-body breakup mode was not considered. The ^{5}He nucleus was assumed to have its first and second excited states at 4 and 8 MeV, so that three channels corresponded to the ground, first and second energy levels. These channels were each assumed to occur with an equal probability of 1/3. These older partition fractions appear close enough to the correct ones to have caused little error in their use.
5. Maxwellian-averaged two-body final state reaction particle spectra
In the center of mass system of two reacting nuclei the kinetic energies of each particle are given by T_{1}¢ = m_{2}T_{R}_{/}M and T_{2}¢ = m_{1}T_{R}_{/}M where T_{R} = T_{1} + T_{2} - T_{cm}. If the final state contains two particles, their center-of-mass energies are then, just as simply, T_{3}¢ = m_{4}(T_{R} +Q)_{/}M and T_{4}¢ = m_{3}(T_{R} +Q)_{/}M. In the limit of a very cold plasma (i.e., q @ 0) the kinetic energies of these exit nuclei would be T_{3}¢ = m_{4}Q_{/}M and T_{4}¢ = m_{3}Q_{/}M. Such energies are sharply defined, and diagnostic measurements would show them to have a very narrow energy spread. In progressively warmer plasmas, i.e., as q increases, the energy spread around these values and the centroid of the distribution also increases; each final state particle energy becomes distributed over some energy range, and thus develops into an energy spectrum that is characteristic of the reaction and the plasma.
In this section we outline a practical analytic integral approach to deterministically calculating this energy spectrum for a thermonuclear reaction in a Maxwellian plasma, which is used in the TDF system. This approach was reported by Williams, and much of the development which follows was taken from his paper (Williams 1971). Williams defines a "super-distribution" for two-body final state nuclear reactions in a plasma with a function f(v_{1},v_{2},v_{3},v_{4}), which measures the probability that the reaction has occurred in such a way that the initial two particles have vector velocities v_{1} and v_{2}, and the final two particles have vector velocities v_{3} and v_{4}. Then the probability that the reaction does occur at all — for any values of v_{3} and v_{4} — is necessarily given by the following integral result:
where s is the total reaction cross-section corresponding to the initial two particles having the velocities v_{1} and v_{2} (or relative velocity v_{R} = |v_{1} -v_{2}|). To ensure conservation of energy and momentum in the reaction process, Williams redefines f(v_{1},v_{2},v_{3},v_{4}) by implicitly defining a new function g(v_{1},v_{2}), as follows:
This form is substituted in the above integral and the integration carried out, making full use of the d function properties. This results in the following "normalization function":
_{}.
When the above expressions for s, f and g are combined, we have a redefined form for s to insert in the integral which formally defines _{}, which is
_{},
or, a little more explicitly,
_{}.
We now employ a useful short cut developed by Williams: we "unpeel" from this multiple integral the presumed rate of production of final state particle 3, by formally omitting the integration over d^{3}v_{3} to obtain the distribution
_{}.
This is formally equal to the number of particle 3's produced in the reaction with velocity v_{3} within cell d^{3}v_{3}in velocity space; the integral of n_{3}d^{3}v_{3} recovers the full statistically averaged reaction rate _{}. It is also clear that n_{3} contains the spectral information associated with the plasma and its distribution functions f(v_{i}). We bring this out more explicitly by noting that the directivity in the production of particle 3 in the plasma is isotropic; the basis for this assumption was presented in Section 3 and is discussed more fully at the end of this section. Then the angular part of the integral of can be carried out with the result that
_{}.
Because dT_{3} is the integration element, we can identify the last integrand as the spectral energy density of reaction product particle 3 in the plasma; it is the number of particle 3's produced with energy T_{3} in the interval dT_{3} per unit volume per unit time. Thus the particle 3 spectral rate density is
This is the formal expression of the energy spectrum of the emission of particle 3 in the isotropic plasma. The evaluation of the integral for n_{3} results in the explicit numerically integrable expression for the energy spectrum dN_{3}/dT_{3} of particle 3 as a function of its energy. The full formal expression of the integral for dN_{3}/dT_{3} is therefore as follows:
_{}
where we have replaced the general plasma particle distribution functions f(v_{i}) by the Maxwellian forms f(m_{i},v_{i},q), as given previously in Section 2.
We sketch here only briefly the procedures we used to reduce this formidable integral. We first express all velocities in terms of v_{3}, v_{cm}, and v_{R}, and then carry out the integration over d^{3}v_{4}, utilizing the momentum d function in the process. Two integrals result, with v_{cm} and v_{R} as the variables of integration. which are not totally separated because the energy d function argument contains both variables. However this argument does contain the cosine of the angle between v_{3} and v_{cm}, so the angular part of the v_{cm} integral is evaluated using v_{3} as the z-axis, while "using up" this d function. Then the integration is carried out over v_{cm}. The "final result", representing much algebra and evaluation of the energy d function over finite angle limits (see Appendix C), is a single integral over the variable T_{R}, as follows:
_{}
where the exponents are given by
_{}.
The functional dependence of dN_{3}/dT_{3} on T_{3} is now clearly evident. In fact, if we integrate dN_{3}/dT_{3} over dT_{3} we indeed recover _{}, as we should. (This integral was also presented by Talley and Hale (1988) but without elaboration.)
Because of particle symmetry in the equations, we obtain the spectral energy distribution for particle 4 from the above merely by switching labels 3 and 4 in all of the preceding derivation. dN_{4}/dT_{4} is obtained, as it were, "by inspection."
Note that the exponents in the above integral can be rewritten as
_{},
and since v_{cm} = v_{3} - v_{3}¢ we can have, alternatively,
_{}
corresponding to the cases where v_{3} and v_{3}¢ are parallel or antiparallel, respectively. Also, _{} and _{} as defined above are the limits for the integration over dT_{cm} , which in turn come from the angular integration of the energy d function over the direction cosine of v_{cm}, or cosq_{cm}. This angular integration is carried out by the steps shown in Appendix C.
We close this section with the resolution of a subtle and historically persistent issue. Our definition of the f and g functions given above was taken from Williams (1971) and contains the total cross-section s rather than the differential cross-section ds/dW. This implies a priori that the details of the reaction angular distribution do not affect particle spectra. Heinrichs et al (1989) included this angular distribution explicitly in their formulation, and obtained an integral for the particle spectra similar but not identical to Williams'. Their spectrum integral also contains s, and they assert that it depends in no way on the angular variation of ds/dW. Perkins (1990) addressed this angular anisotropy issue by carrying out Monte Carlo calculations of two-body final state energy distributions for reactions having markedly non-isotropic differential cross-sections ds/dW. He found empirically that the angular variation of ds/dWhad no effect on the final state particle spectra.
This anisotropy issue is actually moot. We observe that s rather than ds/dW occurs in all joint distribution functions treated here, because s is the probability for the thermonuclear reaction to occur at all. It therefore selects, as it were, only that subset of the particle population in the plasma that participates in the reaction, and it thereby defines the reaction distributions. This selection is independent of final state properties, that is, it holds for any final state configuration that can occur. At no point does the angular distribution of the final state particles become a factor. Both F in Section 3 and g(v_{1},v_{2}) in this section have this independence property. The angular properties of the final state in the Maxwellian plasma — which is inherently isotropic to begin with — are determined by the conservation laws, and these are put into effect by the delta functions in the integrand of the spectral distribution functions dN_{i}/dT_{i}. The statistical reaction average of the direction cosines as treated in Section 3 also vanishes for this same reason.
6. Digression on two-body final state spectrum moments
The integral defining the two-body final state spectrum distribution dN_{3}/dT_{3} can be used to work out the moments of this distribution directly without actually evaluating it, because the exit particle energy T_{3} appears only in the part of the integrand which is the difference of two exponentials. The general integral for the n^{th} moment of the spectral distribution is the spectral expectation average of _{}, or
_{}.
The actual integral over dT_{3} that is needed to work out these moments is
_{}
where a(±) is as given in the previous section. This integral I_{n} is used to evaluate both the numerator and denominator of the expression for _{} because, evidently,
_{},
and the case n = 0 applies to the denominator. Although the multipliers in front of the integrals can cancel out in the ratio, they are not discarded here.
There exist different algebraic approaches for evaluating these auxiliary moment integrals I_{n}. The one we use is to first make the substitution _{} and rewrite the integral in the form
_{}
where _{} and _{}. We reconsider this as the difference of two separate integrals I_{n} = I_{-} - I_{+} where the new ones are
_{}.
These integrals are simplified further by making the substitutions u = x + b in I_{+} and u = x - b in I_{-}. The results are
_{}
and
_{}.
When these two integrals are subtracted, the finite integrals combine to zero, leaving the result
_{}.
This reduces (for any given n) to well-known definite integrals when the difference term is expanded and simplified further.
We are interested particularly in I_{n} for the cases n = 0 and n = 1, which evaluate to
_{} and _{}.
Higher moment integrals may be evaluated similarly. On replacing a and b with their original meanings, carrying out algebraic reductions, and inserting the results appropriately in the moment integrals, we find – as we asserted earlier – that
_{}
which is the denominator of _{}. It is also the numerator for the zeroth moment _{} and the area under the spectrum distribution curve dN_{3}/dT_{3}. Straightforward algebra then produces the following result for the first moment of the spectrum distribution:
_{}
where _{} is as given in Section 3. In other words, the zeroth spectrum moment is obviously unity, and the first spectrum moment is the reaction statistical average of the final state particle kinetic energy T_{3}. This is to be expected: the first moment of a distribution is the mean. Thus we have evaluated the mean exit particle energy by explicit use of the exit particle spectrum distribution. These spectrum average results were also derived by Williams (1971).
7. Data, algorithmic bases and results of _{} calculations
The Thermonuclear Data File (TDF) system as currently implemented uses the cross-section data for five exothermic hydrogen and helium nuclear reactions, which are listed with their Q values as follows:
d(d,n)^{3}He Q = 3.26891 MeV
d(d,p)t 4.03271
t(t,2n)a 11.3321
t(d,n)a 17.5894
^{3}He(d,p)a 18.3532
In Figure 3 we show log-log plots of the total reaction cross-section s(T_{R}) for each of these reactions as a function of the relative kinetic energy T_{R} of the two initial particles in their center-of-mass system, from about 100 eV to 15 MeV. This laboratory frame cross-section data was provided by the LLNL Nuclear Data Group (White, 1990). These are the cross-section functions that are used in all integration procedures described in this report.
Because of the smoothness of these data in log-log space, all interpolations and lookups on these cross-sections are carried out by the cubic Lagrange technique, which provides extremely smooth and highly accurate results even though logarithmic and exponential operations are used. The cubic Lagrange interpolation procedure fits a cubic in Lagrange form through four points, and then the interpolation is carried out between the two inner points. The Lagrange interpolating cubic function has the following form:
f(x) = |
f_{1}(x - x_{2})(x - x_{3})(x - x_{4})/((x_{1} - x_{2})(x_{1}
- x_{3})(x_{1} - x_{4})) + f_{2}(x - x_{1})(x - x_{3})(x - x_{4})/((x_{2} - x_{1})(x_{2} - x_{3})(x_{2} - x_{4})) + f_{3}(x - x_{1})(x - x_{2})(x - x_{4})/((x_{3} - x_{1})(x_{3} - x_{2})(x_{3} - x_{4})) + f_{4}(x - x_{1})(x - x_{2})(x - x_{3})/((x_{4} - x_{1})(x_{4} - x_{2})(x_{4} - x_{3})) |
where f_{1}, f_{2}, f_{3} and f_{4} are the values of the interpolant data at the points x_{1}, x_{2}, x_{3}, and x_{4}. In the case of the cross-section data in the TDF system we use ln s(T_{R}) for the f_{i} and ln T_{R} for the x_{i}. The four adjacent ln T_{R} entries in this array used to construct the cubic are chosen so that the interpolating point lies between the middle two. The LLNL cross-section file energies are spaced so that linear interpolation of the cross-section as a function of energy is accurate to within 1 %. Using the cubic Lagrange procedure on this array in log-log space in the manner just described can provide an interpolation accuracy of the order of 10^{-}^{5}, i.e. .001 %.
Figure 3. Total cross-section for five TDF reactions
The Gamow function integrands for the Maxwellian-averaged reaction rates _{} are similarly smooth, and typical cases are shown in Figure 4 for — by way of example — the t(d,n)a reaction at three different Maxwellian plasma temperatures. The cross-section function for this reaction — which is part of the integrands — is also displayed for comparison.
These integrand functions are displayed again in Figure 5 in a linear-linear plot, but this time renormalized to the same peak value, in order to show off differences in form and width for different plasma temperatures. From Figures 4 and 5 it is seen that the maximum contribution to the _{} integral comes from progressively larger values of s and T_{R} as the plasma temperature qincreases.
Figure 4. Gamow integrand functions and cross-section data
Figure 5. Linear-linear plot of Gamow integrands.
The integration procedure used to calculate _{} is illustrated schematically in Figure 6. For a given temperature q, the energy T_{R} and value at which the Gamow integrand function peaks is determined by table lookup, and the two values of T_{R} at which the function is at 10^{-}^{4} of this peak are calculated by Newtonian root iteration. These two values (_{}and _{}_{)} determine three integration intervals. The integration over the middle (main) interval between _{} and _{} is carried out by one or two 16-point Gauss-Legendre quadratures, depending on the interval size. The integration from 0 to _{} is done by fitting a power law function of T_{R} to the integrand over this interval, and integrating that function. The integration from _{} to ∞ is done by fitting an exponential in T_{R} to the integrand over that range, and integrating that function. The integral for _{} is then the sum of these three integration results. The accuracy of this approach has been verified to six places by using it to evaluate integrals of t^{z}^{-}^{1}e^{-}^{t}dt and t^{-}^{n}e^{-}^{zt}dt, and comparing answers with published values of the Gamma function G(z) and exponential integrals E_{n}(z). (The integral of _{} was also accurately evaluated as a special case of E_{n}(x).)
Figure 6. Gamow function integration intervals
There are three other practical reasons for preferring this algorithmic approach to the integration. First, only the relevant range of the integration variable T_{R} is used, as opposed to doing the numerical integration over the entire data range; the integration is thus efficient and also minimizes computational noise. Second, this integration procedure is well-suited to the other Maxwellian-distribution type integrals required for the TDF system, and is also used for them. Third, this integration technique is extremely fast, and permits function lookups and iterative root findings by direct integration, rather than by interpolation over pre-calculated arrays. It thereby facilitates accurate numerical double integration required to generate the output particle energy spectrum quantities.
The results for the reaction rate integration for the five TDF system nuclear reactions are shown in Figure 7, in which _{} is plotted as a function of the plasma temperature from 1 keV to 1 MeV for each reaction. As with the cross-sections, the reaction rates also range over several orders of magnitude.
Figure 7. Reaction rate _{} for five TDF reactions as
a function of the plasma kinetic temperature.
Previous calculations of _{} have been carried out by many authors. The ones done at LLNL prior to the current TDF system have been reported by Greene (1967) and Howerton (1979). The generation and use of _{} in stellar nucleosynthetic and other astrophysical contexts is outlined, for example, by Caughlan and Fowler in their 1988 compilation and earlier papers.
8. Algorithmic bases and results for _{}, _{} and _{}
The Maxwell-averaging integrals used to calculate the mean particle energies and final state spectra are evaluated in the same manner as for the reaction rates; the integrands are morphologically identical. In Figure 8 we show — using the d(d,n)^{3}He reaction in a 100 keV plasma as an example — the functional forms of the integrands for the reaction rate _{}, the mean relative initial state kinetic energy _{} (which is used in obtaining all the _{}), and the Maxwellian-broadened energy spectrum for each final state particle. The integration procedure for the Gamow function described in Section 7 and illustrated in Figure 6 applies equally well to these other integrands.
Figure 8. Integrand function forms for the Maxwell-averaging procedures,
using the d(d,n)^{3}He reaction in a 100 keV plasma to create this example.
For the Maxwellian averaged mean particle energies only the averaged values of T_{R} and T_{cm} are required since the _{} are linear combinations of _{} and _{} = 3q/2, and _{} is the only average requiring an integration. Example results are shown in Figure 9 as plots of _{}, _{} and the _{} for i = 1, 2, 3, 4 versus plasma temperature q, for the d(d,p)t reaction, where q ranges from 0 to 1000 keV. The indicial correspondence here is that particles 1 and 2 are d's, particle 3 is a p and particle 4 is a t. For the most part these variations are fairly linear in q (and _{} is exactly linear). Indeed, one could approximate _{} ≈ 2q for this reaction. Note that the values of _{} and _{} at zero temperature correspond to m_{i}Q/(m_{3} + m_{4}) for i = 4 and 3 (where Q ≈ 4 MeV), as they should.
Figure 9. Maxwellian averaged mean particle and system energies
for the d(d,p)t reaction as a function of plasma kinetic temperature.
The integrations for the reaction rate spectral distributions dN_{i}/dT_{i} are likewise straightforward, and handled the same way as the other integrations just described. We show in Figures 10a and 10b the final state energy distributions calculated in the TDF system for the two exit particles in the d(d,p)t reaction. The log-log plots in 10a indicate the stability the TDF integration procedures are capable of, as they show how well the regularity and smoothness of the spectrum functions are maintained over several orders of magnitude.
Figure 10a. Final state reaction rate spectral density distribution
for the d(d,p)t reaction at a plasma temperature of 90 keV.
Figure 10b. Final state reaction rate spectral density distribution
for the d(d,p)t reaction at a plasma temperature of 90 keV.
We need to point out that the energy averaging and spectral calculation procedures just described apply only to those TDF system reactions that have two-body final states; the t(t,2n)a reaction, with three particles in the final state, is handled differently. For this reaction no spectral information is provided, as this requires development of three-body final state spectral distributions and their adaptation to thermonuclear plasmas. This will be the subject of a separate and forthcoming report. The reaction-averaged mean energies of the final state particles in the t(t,2n)a reaction, however, are calculated using the approximate multiple channel model developed earlier in Section 4.
9. Algorithmic basis of the final state spectral probability distributions
One of the main uses for final state spectral distribution information is to provide a probabilistic basis for final state configuration sampling in Monte Carlo simulations involving plasma reaction products. The integral over the final state spectral distributions, when treated as a function of the spectral energy, can be used for Monte Carlo sampling. Recall that dN_{i}/dT_{i} is the final state particle spectral distribution function for some particular output particle for a given reaction in a plasma at a given plasma temperature, viz. i = 3 or 4 (or even 5). Let us define a normalized spectral distribution by the function
Then the integral of this normalized distribution as a function of energy is given by
_{}.
Because the area under the spectral distribution dN_{i}/dT_{i} is indeed _{}, this monotone integral function is asymptotically bounded on [0,1], that is, I_{i}(0) = 0 and I_{i}(∞) = 1. It is therefore suitable for the random statistical sampling over [0,1].
It is very useful to define an auxiliary function u_{i}(T_{i}) by the relationship
_{}.
This is inverted to obtain the following defining functional form for u_{i}:
_{}.
We also find that an equally convenient relationship defines S_{i} in terms of u_{i} and its first derivative, as follows:
_{}.
We show in Figure 11 plots of these three functions for the third particle (the neutron) from the t(d,n)a reaction, at a plasma temperature of 90 keV. The display is in three parts, one for each function, with a common energy abscissa.
Figure 11. Spectral, probability and auxiliary functions for the
neutron output from the t(d,n)a reaction at 90 keV.
The reasons for using the (1 + tanh) function to reformulate I_{i}(T_{i}) are, first, that I_{i}(T_{i}) has roughly that form, and, second, the auxiliary variable u_{i} provides a very convenient basis from which to do extremely accurate interpolation over a small collection of points that define the values of S_{i} and I_{i} as a function of T_{i}. The way this works in the TDF system is as follows: when an input random number value for I_{i} is provided, the corresponding u_{i} is determined by calculating u_{i} as a function of I_{i}. Then we interpolate both du_{i}/dT_{i} and T_{i} at this u_{i} and obtain S_{i}. The convenience of this arrangement stems from the fact that T_{i} and u_{i} are very smooth functions of each other over the entire practical energy range, and so the other quantities can be very accurately interpolated — over many orders of magnitude in spectrum amplitude — with very few points. This smoothness appears to be present for all reactions and particle output types in the TDF system. To show the power of this procedure, we used only the 11 points indicated by small squares in each of the first two graphs of Figure 11 to interpolate from, and thus generate all three curves shown.
It is clear that generating the data for the cumulative probability distribution I_{i}(T_{i}) requires numerical double integration, and we detail how this is done in the TDF system. For each spectrum at a prearranged given plasma temperature, we chose 11 points on the spectrum curve — at the peak spectrum amplitude and at the spectrum amplitudes at 0.6, .1, .01, .001 and .0001 of the peak value — to be used for the interpolating point set. The particle energies for each of these lower spectrum amplitudes are determined by using Newton's iterative root finding method on each side of the spectrum peak. In using this method every lookup of a spectrum point was done by actually carrying out the integral for the spectrum value at each iteration (rather than by interpolation). We show in Figure 12 these 11 points for the t(d,n)a reaction at a plasma temperature of 90 keV, which are indicated in the figure by the small open squares lying on the spectrum function curve. (The 11 points shown by open squares in each of the three plots in Figure 11 also have the same significance.)
The values of the cumulative probability function I_{i}(T_{i}) at the interpolation energies T_{i} are obtained by numerically integrating the function S_{i}(T_{i}) between each pair of adjacent points, using Gauss-Legendre quadrature. The S_{i}(T_{i}) values needed for the quadrature are obtained by actually doing the integral for the spectrum value dN_{i}/dT_{i} at each quadrature point T_{i} (as opposed to doing an interpolation) and normalizing to _{}. The double integration is thus in effect a Gaussian quadrature on the spectrum distribution function with the quadrature function values determined by Gaussian quadrature on the spectrum integrand. In this way we obtain the 11 number triplets per spectrum, i.e., S, T and I, for interpolating the particle energy spectrum S and particle energy T from the cumulative probability I at a given plasma temperature. Because of the auxiliary function u, the results of the spectrum lookups are quite accurate over several orders of magnitude. This is suggested by the regularity in the semi-log plot of S_{3}(T_{3}) in Figure 12, which follows, over the range of particle energy and amplitude shown.
Figure 12. Calculated spectral distribution function for the
t(d,n)a reaction at a plasma temperature of 90 keV where
the interpolation array points are indicated by open squares.
10. The Thermonuclear Data File (TDF) generator
The Thermonuclear Data File (hereinafter called TDF) is the data base from which thermonuclear reaction rate and output particle spectral information are "looked up" with TDF-based subroutines, which an application code can use to read the file and generate the thermonuclear information it requires. These subroutines reside in an ASCII source code library file called TDFLIB, and are written in Fortran 77. The TDF is created by a Fortran 77 code called TDFMAKER, which carries out the mathematical and integration operations (as described in Sections 2, 3, 4 and 5) to generate data for each reaction and final state particle in the TDF system. This data is placed in the TDF in a particular standard format for lookup operations. In this section we outline briefly what TDFMAKER does, how it is used, and describe some of its key algorithms. The TDF and its lookup routines in TDFLIB will be described in the next two sections.
TDFMAKER requires as input a file identifying the nuclear reactions to be used, and a separate cross-section data file for each nuclear reaction. The content of an example input file is presented in Table II, which shows the identifying textual data for the five thermonuclear reactions d(d,n)^{3}He, d(d,p)t, t(t,2n)a, t(d,n)a and ^{3}He(d,p)a used in the TDF, respectively. Each line corresponds to a reaction. The first item on each line is the designated (unique but arbitrary) reaction identifier index. The second item is the count of the number of final state (exit) particles of the reaction. The next four (or five) items are the alphameric symbols for the incident particle, target particle and exit particles (in order of increasing A). The last item is the name of the file containing the lab bombarding energy and cross-section data for that reaction. TDFMAKER contains accurate nuclear masses of isotopes from A = 1 to A = 20, inclusive, and it can therefore generate the required kinematic data, including the Q-values.
Table II: Example reaction input for the TDFMAKER code.
After reading the input file with the reaction identifiers, TDFMAKER opens each cross-section file and reads in the lab frame beam energy (in Mev) and 4p cross-sections (in barns) for each reaction. The beam energy is converted to center-of-mass energy T_{R}. Both T_{R} and the cross-section s are stored internally as logarithms of their values so as to perform cubic Lagrange interpolation in log-log space. This interpolation technique has proven good to about 5 to 6 significant figures in this application.
TDFMAKER generates the data for the TDF one reaction at a time, and keeps the generated data for each reaction in separate blocks in the TDF. After the cross-section data is read in and processed, the value of T_{R} for the peak of the Gamow function (i.e., the _{} integrand) is calculated as a function of plasma temperature, since it is used to help determine the three integration intervals of the numerical integration procedure as depicted in Figure 6. In what follows we discuss how TDFMAKER is set up to calculate the integrand peak and the three integration intervals.
Gamow's model of Coulomb barrier penetration accounts very well for the non-resonant variation of thermonuclear reaction cross-sections at very low energies (Evans, 1955). It is used here to estimate the T_{R} at which the _{} integrand peaks, and to provide an algorithm for the initial estimate of the numerical integration intervals. A well-known form of these cross-sections which includes Gamow's model is written as
where S is a slowly varying function of T_{R} that is very nearly constant at keV energies,
_{}
is the Gamow constant, a = 2π e^{2}/hc = 1/137.03604 is the fine structure constant, c is the speed of light, h is Planck's constant, e is the electronic charge, and the Z's and m's are the charge numbers and masses of the reacting isotopes (Evans, 1955). The argument of the exponential is, of course, dimensionless. In the stellar physics field, S(T_{R}) is referred to as the astrophysical S-factor (Chiu 1968, Clayton 1968) and is generally tabulated in cross-section data for stellar nucleosynthesis. Here, for each reaction, we take S to be a constant which normalizes s_{G}(T_{R}) to the reaction s(T_{R}) at low energies. In Figure 13 we plot the cross-sections s(T_{R}) for each reaction, with the corresponding Gamow model cross-section s_{G}(T_{R}) plotted superposed.
Figure 13. Reaction cross-sections: actual data and
corresponding Gamow analytic model s_{G} for each reaction.
At low energies the fit of s_{G} to the various s's for each reaction is excellent, but deteriorates at the higher energies. To ensure accurate determination of the T_{R} of the peak over all data energies, we use a "variable Gamow constant" scheme in TDFMAKER in which g is calculated at each energy point T_{R}(i) of the cross-section array s(i) so as to make s_{G} fit the data in the neighborhood of that T_{R}. Thus we estimate the "local" value of g at the i^{th} point of the cross-section array as follows:
It is worth mentioning that for all five reactions currently used in the TDF system, g(i) as a function of T_{R} is approximately constant up to about T_{R} = 100 keV, after which it varies significantly with higher values of T_{R}.
When s_{G} is used in place of s, the _{} integrand takes the following form:
_{}.
To find the value of T_{R} at the peak of this function we set ∂G/∂T_{R} = 0 and regard S and g as constants. This results in the following power law relationship between T_{R} of the peak and the kinetic temperatureq:
This recipe (i.e., the power law formula using data-fitted "variable" g) is an excellent peak finder for the Gamow integrand, as shown in Figure 14 for the d(d,p)t reaction and in Figure 15 for the t(d,n)a reaction, where we present a log-log plot of q versus T_{R} at the peak for these reactions in three different ways. The curved line represents the power law with the "variable" recipe for g. The straight line represents the power law with g fixed at Gamow's constant. The symbol ¥ locates the value of T_{R} at which the integrand array calculated for the plasma temperature q actually peaks. The "variable g" recipe appears to provide the best peak predictions, while the "constant g" case does well only at lower values of T_{R} and q.
The data for the d(d,p)t reaction shown in Figure 14 follows the "constant g" power law relationship fairly well over most of its range, but the peak T_{R} values for the t(d,n)a reaction seem to approach an 80 keV asymptote, as shown in Figure 15. (This asymptotic behavior of the peak T_{R} is also evident in Figures 4 and 5.) The reason for this is that the cross-section for the t(d,n)a reaction has a pronounced resonance which carries over into the integrand behavior. Because the TDF system is intended for other reactions which may or may not have resonances, the variable g recipe as given above is included in the TDFMAKER algorithmic lexicon as a permanent feature. This feature and the subsequent integration limit determinations underlie most of the integrations that TDFMAKER carries out, which is the reason why a detailed description of it is included in this report.
Figure 14. See text for details.
Figure 15. See text for details.
Given, then, that _{} has been determined for a specified q, we use the Gamow function G(T_{R},q) of the _{} integrand to determine a first estimate of the two integration interval limits shown in Figure 6. Let e be the fraction of the integrand peak value at which these limits occur, and denote by T_{R}^{*} the values of T_{R} which we estimate by this model to occur at these limits. These are computed by determining the values of T_{R}^{*} for which the equality
holds. Taking logarithms of both sides gives the algebraic equation
_{}.
This is a cubic in the variable _{}. Rather than solve a cubic, we approximate the two solutions we need, because at low T_{R}^{*} the quantity g/_{}is dominant, while at high T_{R}^{*} the quantity T_{R}^{*}/q dominates. Let the right-hand side of this equation be symbolized by q. Then g/_{} = q defines _{} and T_{R}^{*}/q = q defines _{}, and we therefore have
_{}
and
_{}.
and _{} are as depicted in Figure 6, where we recall that the three integration intervals are 0 to _{}, _{} to _{}, and _{} to ∞. In TDFMAKER, e has the value 10^{-}^{4}, and the coding has checks to make sure that _{} and _{} indeed bracket _{}. These interval estimate algorithms are applied to all the Maxwell-averaging integrals. When this is done, the inner interval limits are further refined by iterative root finding to be at 10^{-}^{4} of the peak of the integrand.
Recall next that the integrals of the spectrum function dN_{i}/dT_{i} leading to the cumulative probability functions I_{i}(T_{i}) are carried out over adjacent particle energy intervals, as outlined in Section 9. Accurate first estimates of the energy at which the spectrum function peaks and the energies bounding each integration interval can be obtained by use of the exponential difference term in the dN_{i}/dT_{i} integrand, which is given by
_{},
where i is the index label (3 or 4) for the particle of interest, the exponents are
and j is the index label (4 or 3) for the other particle. At a given fixed plasma temperature the dependence of F_{i} on T_{i} is where the variation of the final state spectrum dN_{i}/dT_{i} with particle energy T_{i} originates, so F_{i} determines the shape of dN_{i}/dT_{i}. We also note that F_{i} is also the integrand for the auxiliary integral I_{n} discussed in Section 6 on spectrum moments.
One of the more remarkable insights that emerged during the writing of this report is the realization that F_{i} is a very close approximation to the shape of the spectrum, dN_{i}/dT_{i}, apart from normalization. This is because the variation of F_{i} with T_{R} is quite slow, and to a very good approximation F_{i} can be taken outside the integral for dN_{i}/dT_{i}. In Figures 16a and 16b (as in Figures 10a and 10b) we compare log-log and linear-linear plots of dN_{i}/dT_{i} and F_{i} as functions of T_{i} for the d(d,p)t reaction at q = 90 keV, for both exit particle types. For these F_{i} plots T_{R} is set equal to its mean value _{} at that temperature. The closeness of the shapes of dN_{i}/dT_{i} and F_{i} to each other, for each particle i, is excellent. We thus obtain very accurate first estimates of the integration intervals by determining those pairs of energies T_{i} that make F_{i} equal to the fractions .6, .1, .01, .001, 10^{-}^{4} and 10^{-}^{5}, since the peak of F_{i} is very nearly 1. However, because this approach postdates the TDF system, its implementation is planned for a future version of TDFMAKER.
Figure 16a. Spectral distribution dN_{i}/dT_{i} (curved lines) and function
F_{i} (dots) for the d(d,p)t reaction at a plasma temperature of 90 keV.
Figure 16b. Spectral distribution dN_{i}/dT_{i} (curved lines) and function
F_{i} (dots) for the d(d,p)t reaction at a plasma temperature of 90 keV.
The approach originally developed and currently used in TDFMAKER to determine the spectrum integration intervals is based on first inverting the preceding equation for the exponents a to obtain (for, say, i = 3)
T_{R} is set equal to its value T_{R}^{peak} at the peak of the integrand, because this is where the dN_{i}/dT_{i} integrand makes its maximal contribution to the integral. When a(-) is zero, a(+) is large, and F_{3} is almost 1. The resulting value of T_{3} estimates the energy of the spectrum peak, which we denote by _{}. Next, when a is a large positive number, F_{3} will be much less than 1, and two values of T_{3} (one on either side of the spectrum peak) result, corresponding to the plus/minus signs in the parentheses. We estimate the useful energy limits of dN_{3}/dT_{3} by setting a = ln 10^{4}, so that F_{3} = e^{-}^{a} = 10^{-}^{4}, and we denote the resulting two values of T_{3} by _{} and _{}.
In the current version of TDFMAKER, these three numbers — _{}, _{}, and _{}^{r} — initiate the Newtonian root-finding and peak finder procedures that determine the value of T_{3} at which the spectrum peaks, and the T_{3} pairs (one on each side of the peak) at which the spectrum is at .6, .1, 10^{-}^{2}, 10^{-}^{3}, 10^{-}^{4} and 10^{-}^{5} of the peak. These procedures yield the 11 point triples (S, T and I) in the spectrum data arrays for each exit particle and each plasma temperature that are used to carry out the interpolations done in connection with Figures 11 and 12 of Section 9. It is between these values of T_{3} that the spectrum functions S_{i}(T_{i}) are integrated to obtain the cumulative distribution function point values I_{i}(T_{i}).
The remainder of TDFMAKER consists of coding and algorithms to carry out the integrations, inverse lookups (root findings) and maximum findings described above and referred to in Sections 2 through 9. Two similar Fortran 77 routines in TDFMAKER were written to handle these integrations, one for the "Gamow peak integrand types" used to calculate _{}, _{} and the spectrum distribution dN_{i}/dT_{i}, and the other to help calculate the integrals of the dN_{i}/dT_{i} functions used for the spectrum probability distributions I_{i}(T_{i}). The first subroutine (intgrt) is a straight-forward coding of the integration scheme and limit determinations as described in the previous paragraphs. The second routine (sp3fn) is this subroutine intgrt recoded to be a function call and to carry out only the integration of the spectrum distribution dN_{i}/dT_{i}. This is so that dN_{i}/dT_{i} (in the guise of sp3fn) is used as a "callable" function within TDFMAKER, as needed. This "functionality" is what facilitates the high efficiency in carrying out the integrations over the spectrum distributions, as well as determining where the maximum of these distributions occur. Subroutine mxfnd locates the maximum of non-analytic functions by recursively fitting parabolas to point triplets and calculating the maximum of the parabola. This subroutine is used to more precisely locate the maximum of the Gamow-type integrands as well as that of the spectrum distributions.
Other Fortran 77 subroutines were developed to perform inverse lookup tasks, that is, input the value of a function and output the corresponding argument value by Newtonian root-finding techniques. In TDFMAKER the inverse lookup subroutines all have names beginning with inv. It is within these inverse routines that sp3fn is heavily used as the function value generator in determining energy arguments corresponding to final state particle spectrum distribution values.
11. Structure and content of the TDF
The information in the TDF is in alphameric format so that the TDF can be read as a plain text file. It is therefore platform-independent. The information for each reaction is arranged in separate reaction data blocks within the file, and each block has the same format, which is illustrated by the example of blocks of data from a working TDF that is shown in Appendix A. Every TDF (as of this writing) starts with a header block which is a copy of the input data for TDFMAKER (cf. Table II) that identifies the reactions and the text files containing the reaction cross-section data. Each subsequent reaction data block starts with a line labeling the version of TDFMAKER that was used to generate the block. Within each such block this line is followed by eight different types of sub-blocks of data, each starting with a header line having a unique numeric label in parentheses. The headings for each of these sub-blocks are listed in Table III.
Table III. Data block header lines in the TDF.
Heading line (1) lists the names of the files used by TDFMAKER; the last one is the name of the TDF itself. This line is followed by an alphameric identifier code which the TDFLIB subroutines check to insure that they are reading a compatible version of the TDF. Heading line (2) identifies the constants on the line which follow it, which are, respectively, the MeV energy equivalent of the atomic mass unit, the inverse of the fine structure constant, the speed of light, the Q of the reaction, and the total energies of the final state particles if those were left in excited states by the reaction. Heading line (3) identifies the nuclear masses for each particle in the reaction, in atomic mass units. Heading line (4) identifies two groups of integer data which follow it. The first five numbers are the isotope identifiers for each particle in the reaction; these identifiers have the integer value 1000Z + A where Z is the nuclear charge (number of protons) and A is the atomic number (total number of protons and neutrons). The last four numbers are respectively the number of temperature blocks for each final state particle spectrum (heading (8)), the number of lines of data in each _{} block (heading (5)), the number of final state particles (a 2 or a 3), and the numerical index/identifier used within the TDF for that reaction. Heading line (5) identifies the data block containing the values of q, _{} and _{} for this reaction. Currently (as shown in Appendix A) there are 29 lines of these data triplets, each for a different plasma temperature. Also, there are 13 temperature blocks for each final state particle spectrum, where each of the 13 blocks corresponds to a different plasma temperature. Thus if there are two particles in the final state for a reaction there will be 26 spectrum blocks in that reaction section of the TDF. The 13 plasma temperatures q_{i} are the same for all spectrum blocks, and are equal to .0001, .001, .01, .04, .09, .16, .25, .36, .49, .64, .81, 1.00, and 1.05 MeV. The first and last temperature entries are included to preserve accuracy of the cubic Lagrange interpolation within the temperature intervals .001 to .01 and .81 to 1.00.
The remaining data blocks after line (5) refer to the output particle spectra. Heading line (6) contains the number of fractions of the spectrum peak used to define the spectrum interpolation points, and these fractions are listed in the next line. Heading line (7) identifies (for the spectrum block which follows) the plasma temperature and the final state particle index to which the block applies (a 3 or a 4 or (later) a 5), the number of lines of spectrum data, the integrated area under the spectrum curve, and the corresponding value of _{}. (A comparison of these last two numbers gives computational diagnostic information.) Heading line (8) identifies the spectrum data block and lists, respectively, the line index j and the values of T_{j}, S_{j} and I_{j} for the given output particle at the temperature q_{i} specified. These are the actual interpolation data points used for the spectral information lookups. Currently there are 11 lines of these data quartets per block, corresponding to the peak value of the spectrum distribution and five pairs of fractional values (currently 0.6, 0.1, 10^{-}^{2}, 10^{-}^{3}, and 10^{-}^{4}) of this spectrum peak value. That is, S_{1} = S_{11} = 10^{-}^{4} S_{6}, S_{2} = S_{10} = 10^{-}^{3} S_{6}, S_{3} = S_{9} = 10^{-}^{2} S_{6}, S_{4} = S_{8} = 0.1 S_{6}, and S_{5} = S_{7} = 0.6 S_{6}, where S_{6} is the peak value.
As shown in Appendix A, the last five entries of I_{j} (i.e. index 7 through 11) in each spectrum data block are negative. This is because the numerical value of 1. was subtracted from those I_{j} in TDFMAKER before putting it out in the TDF, in order to retain the double precision numerical accuracy in the TDF. When these negative numbers are read in by the TDFLIB routines, the value of 1. is added to them to restore them to their original (and positive) double precision value. (TDFMAKER and TDFLIB currently declare double precision for use with 32-bit processors.)
12. The TDFLIB subroutine library
TDFLIB is the name of an ASCII Fortran 77 source program library file that contains 8 subroutines, one function subroutine, and specifications for 3 labeled common blocks. Four subroutines — tdfrdr, mxavlu, speclu, tdmass — and one function subroutine — svblu — are directly used by the applications code that needs the thermonuclear information. Subroutines tdclag, tdchlu, tdmchl and tdselu are specialized interpolation subroutines based on cubic Lagrange and Hermite methods that are intended to be used only by the other subroutines. The three labeled common blocks are tdfcon, tdfsvb and tdfspe; their memory allocation details are spelled out in comments in subroutine tdfrdr.
The main function of the TDFLIB subroutines is to perform monovariate and bivariate interpolations on the condensed data in the TDF, so as to retrieve Maxwellian-averaged and reaction-averaged thermonuclear energy and spectral information, as desired. In what follows, we present a brief description of these subroutines, their argument lists and their purpose, and then follow this by some spectrum data interpolation examples that demonstrate the calculationally robust nature of the TDF system over several orders of magnitude in the variables. Some of these routines were used to generate the plotted data shown in previous sections.
subroutine tdfrdr(iunit,ntr,ierr)
Reads data from the TDF into the three labeled common blocks in memory. It must be called no more than once by the application code before any other routine can be used. iunit is the logical unit number assigned to the TDF by the application code. ntr is the number of thermonuclear reactions read in. ierr is an error flag: zero means that the TDF was read and the memory allocated correctly. A non-zero return means the TDF was not read properly and cannot be used; this is a fatal error indicating that either the TDF was corrupted, or the version of TDFMAKER used to create the TDF, and the current TDFLIB, do not match.
In all the other subroutines, the quantity ierr is a returned error flag, in which zero means that the subroutine call was completed correctly. A non-zero error flag occurring for the main four subroutines and one function usually means that the input reaction or particle index, or one or more input numerical values, was improperly specified. In some cases, however, a more serious error could be indicated. The error flags are explained in the documentation for TDFLIB (White, Warshaw and Resler, 1990).
subroutine tdmass(ireac,bm1,bm2,bm3,bm4,bm5,q,ierr)
Returns purely nuclear masses (i.e. atomic mass minus electron masses and electron binding energies) and the Q-value q (in MeV) for a specified reaction ireac. The mass units are in amu where ^{12}C = 12 amu.
bm1 = mass of lighter reacting particle
bm2 = mass of heavier reacting particle
bm3 = mass of lightest product particle
bm4 = mass of next lightest product particle
bm5 = zero for two-body final state reactions, else is the mass of the remaining product particle in three-body final state reactions
subroutine mxavlu(ireac,tp,sv,e1,e2,e3,e4,e5,ierr)
Returns the Maxwellian-averaged reaction rate sv (cm^{3}/sec) and reaction-averaged kinetic energy (MeV) of the reacting and final particles for the specified reaction ireac and plasma temperature tp (keV).
e1 = kinetic energy of the first reacting particle
e2 = kinetic energy of the second reacting particle
e3 = kinetic energy of the first final particle
e4 = kinetic energy of the second final particle
e5 = zero for 2-body final state reaction or kinetic energy of third final particle
mxavlu identifies the data array index for the specified reaction, and determines the four sequential data temperatures which bracket the input plasma temperature tp. It calls tdclag to interpolate the reaction rate sv at this input temperature, and tdchlu to interpolate the mean relative kinetic energy er of the two reacting particles. From this and the plasma temperature it computes the reaction-averaged kinetic energies of all particles of the reaction, e1 through e5. mxavlu was used to generate the averaged kinetic energy data plotted in Figure 9, and was modified specifically to also output the average of T_{R} for the plot. (Note that er is the _{} of this report, and the ei arguments are _{}, correspondingly.)
function svblu(ireac,tp,ierr)
Returns the Maxwell-averaged reaction rate svblu (cm^{3}/sec) for the specified reaction ireac and plasma temperature tp (keV). svblu identifies the data array index for the reaction, and determines the four sequential data temperatures which bracket the input temperature tp. It calls tdclag to interpolate the reaction rate over this temperature set. svblu was used to generate the reaction rate plots in Figure 7.
subroutine speclu(ireac,tp,ipart,pin,eout,sout,ierr)
speclu is an inverse particle spectrum lookup routine. Given the input plasma temperature tp (keV) and cumulative probability of emission pin (a real number between 0. and 1.) of exit particle ipart (3 or 4), it returns that particle's energy eout (MeV) and normalized spectral amplitude (i.e., relative probability) sout (MeV^{-}^{1}). The actual spectral amplitude is obtained by multiplying sout by the reaction rate. It is currently implemented for two-body final state reactions only. It returns a non-zero error flag and zero arguments if a reaction with three or more final state particles is selected.
ipart = 3 for the first product particle
ipart = 4 for the second (heavier) product particle
ipart = 5 is intended for a third product particle (not currently implemented)
speclu first checks the range of the input value of pin and returns special values of eout and sout if pin is exactly 0 or 1. For input values of pin between 0 and 1, speclu determines the four sequential data temperatures (and corresponding spectral data blocks in memory) which bracket the user's requested temperature tp. At each of these four temperatures it calls tdselu to interpolate the energy (e) and spectrum amplitude (s) at the user's specified cumulative probability input (pin). The resulting four values each of e and s are then interpolated at the user's specified temperature tp, by calls to tdchlu or tdmchl, as appropriate, and output as eout and sout, respectively. speclu was used to generate the data plotted in figures 10a, 10b, 11, 12, 16a and 16b.
subroutine tdclag(x1,x2,x3,x4,f1,f2,f3,f4,x,fval,ierr)
Cubic Lagrange interpolation routine for mxavlu and svblu. Interpolant x lies between x2 and x3.
subroutine tdselu(nnode,pfnode,ennode,spnode,p,e,s)
Spectrum-energy lookup routine for use by speclu. Input is p, returns are e and s, assuming 0 < p < 1. The formulas developed in Section 9 are used to interpolate e and s over the 11-point data set. Fitted power law and exponential functions are used to extrapolate e and s below and above this data set, respectively.
subroutine tdchlu(x1,x2,x3,x4,f1,f2,f3,f4,x,fval,ierr)
4-point cubic Hermite-Lagrange slope-value interpolation routine for mxavlu and speclu. Interpolant x lies between x2 and x3.
subroutine tdmchl(x1,x2,x3,x4,f1,f2,f3,f4,x,fval,ierr)
Mixed 4-point cubic Hermite-Lagrange slope-value interpolation routine for speclu. Interpolation is in linear space on nodes 1, 2 and 3 and in logarithmic space on nodes 2, 3 and 4. Interpolant x lies between x2 and x3.
The interpolation subroutines in TDFLIB were carefully checked for consistency and robustness before release. One such test was to generate contour plots of the bivariate functions I_{i}(S_{i},q) and I_{i}(T_{i},q) for the two output particle spectra (corresponding to i = 3 and i = 4) of each two-body final state reaction, and inspect them for smoothness over the entire variable range in the TDF. The spectrum data (T_{i} or S_{i} as appropriate) were interpolated from the TDF for 15 probability values I_{i} and 50 logarithmically spaced values of the plasma temperature q between 1 keV and 1 MeV, thus generating 750 data quartets (q, T, S, I) for each reaction and each output particle. From these quartets for each reaction we chose to generate two types of plots consisting of 15 separate curves, one for each value of I, for each of the two output particles. They are T_{3} versus q^{1/2} and T_{4} versus q^{1/2} on semi-log scales, and S_{3} versus qand S_{4} versus qon log-log scales. These curves are interpreted as contour plots of I.
Four such plots are shown in Figures 17 through 20 for the d(d,p)t reaction, in each of which the 15 contour levels of the cumulative probability function I_{i} are shown for the two functions, for each output particle. Figures 17 and 18 indicate approximately how the spectrum peak width increases with plasma temperature, where the various curves of I_{i} define the points on each side of the curve of the peak that are being followed. Figures 19 and 20 are double-sheeted in I_{i} since at each q there are two values of I_{i} for every S_{i}.
The essential point of these plots is that none of the 750 interpolated quartets coincide with any of the 143 spectrum data quartets for each reaction output particle in the TDF. Therefore, all calculations of the I_{i}(S_{i},q) and I_{i}(T_{i},q) functions for the two output particle spectra are done by bivariate Lagrange and Hermite cubic interpolation over the more limited plasma temperature and spectrum level data points in the TDF. The results are visually excellent over many orders of magnitude. These types of plots were used to choose and "shake down" various bivariate interpolation schemes during the code development phase. The ragged and "noisy" contours that were produced in the early coding stages were the best diagnostic that the interpolation procedures were not performing properly.
Figure 17. Contours of I_{3}as a function of T_{3} and q for the d(d,p)t reaction.
Figure 18. Contours of I_{4}as a function of T_{4}and q for the d(d,p)t reaction.
Figure 19. Contours of I_{3}as a function of S_{3}and q for the d(d,p)t reaction.
Figure 20. Contours of I_{4}as a function of S_{4}and q for the d(d,p)t reaction.
A word on the nature of the bivariate interpolation used in subroutine speclu is necessary, because it is rather specialized. Given that a plasma temperature q and a cumulative probability value I is input to a spectrum data calculation for a given reaction and a given output particle, the two spectrum data blocks at the two plasma temperatures immediately above and immediately below the input value q is are used to generate, for that input I and those four plasma temperatures q_{j}, four values of the spectrum amplitude S_{j} and particle kinetic energy T_{j}. These values of S and T are obtained by cubic Lagrange interpolation over the auxiliary function u that was described earlier in Section 9. For definiteness, we denote these four plasma temperatures by q_{1}, q_{2}, q_{3}, and q_{4}, in monotone increasing value, such that the input value q lies between q_{2}, and q_{3}, and we index the corresponding values of S and T similarly. The required temperature-interpolated values of S and T are obtained by modified cubic Hermite interpolation over the indexed Sand T. The way this is done is to convert the four plasma temperatures to their square root value, q_{j}^{1/2} and then regard S_{j} and T_{j} as functions of q_{j}^{1/2}, where j = 1, 2, 3 and 4. Then the slopes of these functions are modeled at the two inner abscissae — q_{2}^{1/2} and q_{3}^{1/2} — by fitting parabolas to the three points at the abscissa triplets (q_{1}^{1/2}, q_{2}^{1/2}, q_{3}^{1/2}) and (q_{2}^{1/2}, q_{3}^{1/2}, q_{4}^{1/2}). This gives slopes and values for the functions S and T at the inner abscissa pairs q_{2}^{1/2} and q_{3}^{1/2}. The cubic Hermite interpolation procedure is then carried out over these slope-value pairs to obtain the values of S and T at the square root of the input q.
The author would like to thank certain people for their assistance and contributions to the development and implementation of the LLNL TDF system described in this report. First, the late Ted Perkins, for bringing the key paper written by Williams (1971) to the author's attention, and for providing his own annotated copy of that paper to the author, with much useful commentary scribbled in the margins. Roger White and David Resler of the Physics Department's former Nuclear Data Group, for providing evaluated nuclear cross-section files and nuclear data subroutines used by the TDFMAKER code to retrieve nuclear masses and reaction cross-sections. Roger White, who, as Nuclear Data Group leader, provided a hospitable environment for development of the TDF system, and Eugene Brooks, who, as the succeeding group leader for Nuclear Data and Codes, provided support for the writing of this report now ten years later. Eugene Canfield, for indicating that the idea of conjoining the nuclear cross-section with the Maxwellian plasma distribution as the required joint distribution for the energy and spectrum averaging may have originated with John von Neumann. And, finally, it is a pleasure to acknowledge the many helpful discussions with and comments made by George Zimmerman, who took a keen interest in the development of the TDF system and its exploratory extensions to include spectra for three-body final state reactions. This report is dedicated to the memory of Sterrett "Ted" Perkins.
References to thermonuclear reaction data and modeling listed here represent neither an extensive nor an exhaustive survey of the field. Rather, they are citations of external and internal reports that the author has found helpful and suggestive in the course of developing the modeling and calculational approaches described in this report. Not all references listed were specifically cited.
F. Ajzenberg-Selove, "Energy levels of light nuclei A = 5 - 10", Nuclear Physics A490, 1 - 225 (1988)
R. E. Brown and N. Jarmie, "Differential Cross-sections at Low Energies for ^{2}H(d,p)^{3}H and ^{2}H(d,n)^{3}He", Physical Review C41, 1391 - 1400 (1990)
R. E. Brown, N. Jarmie and G. M. Hale, "Fusion-Energy Reaction ^{3}H(d,a)n at Low Energies", Physical Review C35, 1999-2004 (1987)
H. Brysk, "Fusion Neutron Energies and Spectra", Plasma Physics 15, 611 - 617 (1973)
G. R. Caughlan and W. A. Fowler, "Thermonuclear Reaction Rates V", Atomic Data and Nuclear Data Tables 40, 283 - 334 (1988) and earlier references cited therein.
Hong-Yee Chiu, "Stellar Physics", Volume I, pp. 349 - 363, Blaisdell, Waltham MA, 1968
D. D. Clayton, "Principles of Stellar Evolution and Nucleosynthesis", Chapter 4, pp. 283 - 361, McGraw-Hill, New York, 1968
R. D. Evans, "The Atomic Nucleus", pp. 71 - 74 and pp. 874 - 878, McGraw-Hill, New York, 1955
S. L. Greene, Jr., "Maxwell Averaged Cross Sections for Some Thermonuclear Reactions on Light Isotopes", LLNL report UCRL-70522, May 1967. Preprint for paper in journal Nuclear Physics.
G. Heinrichs, et al. "Analytical Calculation of Source Energy Spectra from Maxwellian Fusion Plasmas of Arbitrary Temperature", Proceedings of the 12th International Conference on Plasma Physics and Controlled Nuclear Fusion Research, Nice, 12 - 19 October 1988, Volume 3, pp. 459 - 468, IAEA-CN-50/G-III-6, Vienna, 1989
R. J. Howerton, "Maxwell-Averaged Reaction Rates for Selected Reactions Between Ions with Atomic Mass £ 11", LLNL report UCRL-50400, Volume 21, Part A, February 1979
N. Jarmie, R. E. Brown and R. A. Hardekopf, "Fusion-Energy Reaction ^{2}H(t,a)n from E_{t} = 12.5 to 117 keV", Physical Review C29, 2031 - 2046 (1984)
S. T. Perkins, "TNSPEC, a code to calculate particle spectra and energy deposition terms from thermonuclear reactions", LLNL Nuclear Data Group internal memorandum PD-61, July 1981
S. T. Perkins, "Average energy of interacting particles for thermonuclear reactions", LLNL Nuclear Data Group internal memorandum PD-62, July 1981
S. T. Perkins, "TNSP4DP, a code to calculate particle spectra and energy deposition terms from continuum thermonuclear reactions", LLNL Nuclear Data Group internal memorandum PD-72, July 1982
S. T. Perkins, "Average interaction energy of a suprathermal particle in a Maxwellian plasma", LLNL Nuclear Data Group internal memorandum PD-158, April 1988
S. T. Perkins, "The influence of the center of mass angular distribution on two-body thermonuclear particle spectra", LLNL Monte Carlo Group internal memorandum MC-10, February 1990
S. T. Perkins, "Brief comments on the center of mass system collision cosine for two body reactions on moving targets", LLNL Monte Carlo Group internal memorandum MC-12, January 1991
S. T. Perkins and R. J. Howerton, "Average Energy of Interacting Particles for Thermonuclear Reactions", LLNL report UCRL-87543. Also Nuclear Science and Engineering 82, 388 - 392 (1982)
D. R. Slaughter, "Fusion Reactivities for Several Beam and Target Ion Distributions", LLNL report UCRL-87416, September 1982. Also Journal of Applied Physics 54, 1209 - 1217 (1983)
D. R. Slaughter, "Fusion-Product Spectral Line Width in D-D, D-T and ^{3}He-D Plasmas as an Ion Speed Distribution Diagnostic", LLNL report UCRL-97595, October 1987. Also Review of Scientific Instruments 60, 552 - 61 (1989)
D. R. Slaughter, "Effect of Ion Speed Distribution on the Spectral Shape of the 2.5 MeV Neutron Line Source Produced by DD Fusion", LLNL report UCRL-93895, February 1986. Also Review of Scientific Instruments 57, 1751 - 1753 (1986)
D. R. Slaughter, "Fusion Product Energy Spectra in Beam Heated D-D, D-T and D-^{3}He Plasmas", LLNL report UCRL-93946, March 1986
D. R. Slaughter, "LINE: a code which simulates spectral line shapes for fusion reaction products generated by various speed distributions", LLNL internal report UCID-20374, March 1985
T. L. Talley and G. M. Hale, "Fusion Product Spectra", pp. 299 - 302, proceedings of the 1988 MITO Conference on Nuclear Data for Science and Engineering, JAERI, Japan
S. I. Warshaw, "Burning Issues in Maxwell-Averaged Reaction Rates, Probability Distributions, Data Files and Lookup Codes", LLNL Nuclear Data Group internal memorandum PD-180, Revision I, 1991
R. M. White, LLNL Nuclear Data Group, N Division, private communication (1990)
R. M. White, D. A. Resler and S. I. Warshaw, "Evaluation of Charged Particle Reactions for Fusion Applications", LLNL report UCRL-JC-107158, 7 pp. (May, 1991); Also pp. 834 - 839, proceedings of the International Conference on Nuclear Data For Science and Technology, Julich, Germany, 13 - 17 May 1991.
R. M. White, S. I. Warshaw and D. A. Resler, "Documentation for the files TDF and TDFLIB", LLNL Nuclear Data Group internal memorandum PD-179, December 1990
M. M. R. Williams, "A Generalized Energy Exchange Kernel For Inelastic Neutron Scattering and Thermonuclear Reactions", Journal of Nuclear Energy 25, 489 - 501 (1971)
M. M. R. Williams, "The Particle Spectrum from a D-T Plasma", Nuclear Science and Engineering 73, 219 - 220 (1980)
C. Wong, J. D. Anderson and J. W. McClure, "Neutron Spectrum from the T + T Reaction", Nuclear Physics 71, 106 - 112 (1965)
G. B. Zimmerman, LLNL Defense and Nuclear Technology Directorate, A Division, private communications (1990, 2001).
Appendix A. Example TDF format
-------index/input-------
1 2 d d n he3 ddnevalcs
2 2 d d p t ddpevalcs
3 3 t t n n a tt2nevalcs
4 2 d t n a dtevalcs
5 2 d he3 p a dhe3evalcs
-------------------------
TDFMKRdp version 8a of 11/12/91
(1) ddnevalcs tdfinput tdfout
8A3E911111
(2) amumev fnstnv celer qvalue exval
9.31501E+02 1.37036E+02 2.99792E+10 3.26891E+00 0.00000E-01
(3) amu1 amu2 amu3 amu4 amu5
2.01355E+00 2.01355E+00 1.00866E+00 3.01493E+00 0.00000E-01
(4) iaz1 iaz2 iaz3 iaz4 iaz5 nspec nsigv npart idreac
1002 1002 1 2003 0 13 29 2 1
(5) temp sigvb
extr
1.00000E-04 1.59243E-07 1.43504E-03
1.25000E-04 2.51163E-06 1.67293E-03
1.80000E-04 1.46593E-04 2.15117E-03
2.50000E-04 3.74694E-03 2.70057E-03
3.60000E-04 9.01754E-02 3.48017E-03
5.00000E-04 1.13287E+00 4.37873E-03
7.20000E-04 1.35184E+01 5.65891E-03
1.00000E-03 9.68333E+01 7.14118E-03
1.40000E-03 5.77409E+02 9.07899E-03
2.00000E-03 3.05447E+03 1.17374E-02
2.80000E-03 1.21623E+04 1.49933E-02
4.00000E-03 4.40668E+04 1.94980E-02
5.60000E-03 1.28096E+05 2.50697E-02
8.00000E-03 3.46285E+05 3.28673E-02
1.10000E-02 7.58454E+05 4.20404E-02
1.60000E-02 1.71175E+06 5.64938E-02
2.20000E-02 3.15706E+06 7.30058E-02
3.20000E-02 5.98733E+06 9.93100E-02
4.50000E-02 1.00430E+07 1.32057E-01
6.40000E-02 1.61926E+07 1.77967E-01
9.00000E-02 2.44813E+07 2.38144E-01
1.25000E-01 3.49645E+07 3.15636E-01
1.80000E-01 4.96368E+07 4.31629E-01
2.50000E-01 6.54818E+07 5.72492E-01
3.60000E-01 8.55791E+07 7.84893E-01
4.90000E-01 1.04190E+08 1.02797E+00
6.40000E-01 1.21178E+08 1.30202E+00
8.10000E-01 1.36601E+08 1.60651E+00
1.00000E+00 1.50556E+08 1.94023E+00
(6) spectrum levels, n = 5
6.00000E-01 1.00000E-01 1.00000E-02
1.00000E-03 1.00000E-04
(7) temp particle npoints
area sigvb
1.00000E-04 3 11 1.59173E-07
1.59243E-07
(8) i energy spectrum
integral
1 2.4031589E+00 3.5982959E-03
8.8424660E-06
2 2.4094762E+00 3.5982801E-02
1.0000389E-04
3 2.4169803E+00 3.5982773E-01
1.1942821E-03
4 2.4267775E+00 3.5982847E+00
1.5847336E-02
5 2.4393182E+00 2.1589762E+01
1.5551904E-01
6 2.4505695E+00 3.5982939E+01
5.0117482E-01
7 2.4617310E+00 2.1589762E+01
-1.5660226E-01
8 2.4743616E+00 3.5982850E+00
-1.6027927E-02
9 2.4842743E+00 3.5982782E-01
-1.2123457E-03
10 2.4918941E+00 3.5982817E-02
-1.0181209E-04
11 2.4983269E+00 3.5982980E-03
-9.0249845E-06
(7) temp particle npoints
area sigvb
1.00000E-04 4 11 1.59173E-07
1.59243E-07
(8) i energy spectrum
integral
1 7.7295021E-01 3.5990631E-03
8.6594925E-06
2 7.7915141E-01 3.5990511E-02
9.8207865E-05
3 7.8653932E-01 3.5990483E-01
1.1763575E-03
4 7.9621976E-01 3.5990548E+00
1.5668177E-02
5 8.0866874E-01 2.1594376E+01
1.5444567E-01
6 8.1983442E-01 3.5990628E+01
4.9730239E-01
7 8.3107668E-01 2.1594376E+01
-1.5768205E-01
8 8.4379420E-01 3.5990547E+00
-1.6207574E-02
9 8.5381973E-01 3.5990480E-01
-1.2303020E-03
10 8.6155273E-01 3.5990500E-02
-1.0360734E-04
11 8.6809902E-01 3.5990632E-03
-9.2044831E-06
(7) temp particle npoints
area sigvb
1.00000E-03 3 11 9.67906E+01
9.68333E+01
(8) i energy spectrum
integral
1 2.3063630E+00 1.1358504E-03
8.6424999E-06
2 2.3259784E+00 1.1358504E-02
9.8041043E-05
3 2.3493540E+00 1.1358504E-01
1.1746919E-03
4 2.3799935E+00 1.1358477E+00
1.5651439E-02
5 2.4194130E+00 6.8151017E+00
1.5434509E-01
6 2.4547851E+00 1.1358504E+01
4.9713511E-01
7 2.4904148E+00 6.8151017E+00
-1.5778328E-01
8 2.5307383E+00 1.1358480E+00
-1.6224601E-02
9 2.5625397E+00 1.1358504E-01
-1.2320086E-03
10 2.5870770E+00 1.1358504E-02
-1.0377778E-04
11 2.6078542E+00 1.1358504E-03
-9.2215053E-06
(7) temp particle npoints
area sigvb
1.00000E-03 4 11 9.67908E+01
9.68333E+01
(8) i energy spectrum
integral
1 6.7757772E-01 1.1370148E-03
8.0650292E-06
2 6.9602711E-01 1.1370141E-02
9.2370902E-05
3 7.1823285E-01 1.1370141E-01
1.1180973E-03
4 7.4769523E-01 1.1370124E+00
1.5085648E-02
5 7.8618278E-01 6.8220841E+00
1.5095047E-01
6 8.2126448E-01 1.1370141E+01
4.9147704E-01
7 8.5711205E-01 6.8220841E+00
-1.6117732E-01
8 8.9828560E-01 1.1370123E+00
-1.6790136E-02
9 9.3119974E-01 1.1370141E-01
-1.2885624E-03
10 9.5685723E-01 1.1370141E-02
-1.0943324E-04
11 9.7875840E-01 1.1370151E-03
-9.7870434E-06
Appendix B. Erratum in Williams' original paper.
In the course of developing the equations for the two-body final state particle spectra (viz., dN_{i}/dT_{i}, i = 3, 4) that are presented in Section 5 of this report, we made use of Williams' 1971 paper, and uncovered a nontrivial error in it, which, fortunately, does not affect the final formula for the spectra. It does, however, cause the "general energy exchange kernel" (Williams' equations 11 and 33, which corresponds in part to our function g(v_{1},v_{2}) as given in Section 5) to have the wrong magnitude, to not be symmetric under exchange of the two final state particle indices 3 and 4, and to be dimensionally incorrect. The multiplicative factor in this kernel is M_{3}(M_{3}/M_{4})^{1/2}, which does not display the interchange symmetry. The correct factor is (M_{3}M_{4})^{3/2}. This is the same factor (m_{3}m_{4})^{3/2} in the numerator of our normalization function g(v_{1},v_{2}) as derived in Section 5.
The discrepancy occurs at Williams' improper evaluation of the integral of the momentum d function in his equation 4 to reach his equation 6. Williams treats this three-dimensional integral as if it were one-dimensional. The integral is, in Williams' notation,
_{}.
The correct integration over dV_{4} that uses up the momentum d function produces a multiplier of 1/_{}, not the 1/M_{4} that Williams has. This is because a d function with a vector argument and vector differential, e.g., d(MV)dV, is symbolic for the product of d functions each with a component argument and differential, viz. d(MV_{x})d(MV_{y})d(MV_{z})dV_{x}dV_{y}dV_{z}. Thus the momentum integral over dV_{4} is actually the product of three "component" integrals, each of which contributes a factor 1/M_{4} to the magnitude. The appropriate expressions in Williams' paper from equation 6 on need to be corrected for this multiplicative error.
Appendix C. Evaluation of a d function integral over finite limits
The derivation of the output particle energy spectrum integral dN_{3}/dT_{3} requires, at two places, the evaluation of an integral containing a d function over a finite interval. This results in an expression that contains the interval endpoints, and not necessarily just the zero of the d function argument which would be the case if the interval were of infinite extent. The author believes that this aspect of d function integration is not well known, and is providing a discussion of it here as a helpful reference.
A basic d function property is that if the argument of thed function in an integrand is zero in some integration interval, then the integral is equal to the integrand evaluated at that argument zero. Consider the following integral, where the argument of the d function is a linear function in the integration variable:
The d function argument is zero when ax = b, or x = b/a, and we expect that the integral has the value 1/a when b/a falls between x_{1} and x_{2}, and zero otherwise. To realize this analytically, we make use of the Heaviside step function H(x), which is equal to 1 for positive x and zero for negative x. It can be defined as the integral of the delta function:
Therefore it follows that
_{},
and so
_{}
The point of this discussion is that this difference of the step functions is a boxcar function of a and b which defines the interval over which the integral is not zero. If this integral is part of an integrand, and the integrand variable is contained in a and b, then the boxcar function limits defined by ax_{2} = b and ax_{1} = b will modify the integrand variable limits if their intervals overlap. We now consider two such cases in the development of the spectrum integral dN_{i}/dT_{i} in Section 5.
The first case occurs during the evaluation of g(v_{1},v_{2}) and the relevant integral is
which becomes, after integrating over dj_{3},
where m_{3} = cosq_{3}, a = p_{0}p_{3}/m_{4}, b = (p_{0}^{2} + p_{3}^{2})/(2m_{4}) + p_{3}^{2}/(2m_{3}) - (T_{1} + T_{2} + Q), p_{0} = p_{1} + p_{2}, and q_{3} is the angle between p_{0} and p_{3}. Carrying out the integration over dm_{3} while using p_{0} as the angular reference direction results in
What happens now is that the boxcar function in the brackets effectively changes the integration limits, because the endpoints of the boxcar interval redefine the limits of integration over dp_{3}. That is, the limits _{} and _{} are obtained by solving the conditions a ± b = 0 at the zeros of the Heaviside function arguments. Then the integral easily evaluates to
The second case occurs during the evaluation of dN_{3}/dT_{3}, and the relevant integral is
_{}
where the symbols have the same meaning as above, and the d function argument is the same, but the integration variable this time is v_{cm}. Note that _{}, and that p_{0} = Mv_{cm} by definition of v_{cm}, so the directions of p_{0} and v_{cm} are the same. This means that q_{3} in d^{3}p_{3} is the same as q_{cm} in d^{3}v_{cm}: both are the angle between p_{3} and p_{0} = Mv_{cm}. Thus the angular integration for I˝ is the same as for I´ and proceeds in the same way.
After integrating over dj_{cm}, the integral I˝ becomes
_{}
where a and b have the same meaning as before. Carrying out the integration over dm_{cm} with p_{3} as the angular reference direction gives the result
_{}
As in the previous case, the boxcar function in the brackets will further limit the range of integration to lie between the limits _{} and _{}. Those are obtained by solving the conditions _{} at the zeros of the Heaviside function arguments, but this time for v_{cm}, or, equivalently, p_{0}/M. When this is done we obtain
which is recognized as part of the integrand for dN_{3}/dT_{3}. The exponents in the form shown are discussed at the end of Section 5.