- Email: [email protected]
Refinement of thermonuclear reaction rates
Fusion Engineering and Design 141 (2019) 51–58
Contents lists available at ScienceDirect
Fusion Engineering and Design journal homepage: www.elsevier.com/locate/fusengdes
Refinement of thermonuclear reaction rates A.A. Belov a b c
a,b,⁎
c
, N.N. Kalitkin , I.A. Kozlitin
T
c
Faculty of Physics, Lomonosov Moscow State University, Leninskie Gory 1, bld. 2, Moscow 119991, Russian Federation Peoples’ Friendship University of Russia (RUDN University), Miklukho-Maklaya St., 6, Moscow 117198, Russian Federation Keldysh Institute of Applied Mathematics of Russian Academy of Sciences, Miusskaya Sq., 4, 125047 Moscow, Russian Federation
A R T I C LE I N FO
A B S T R A C T
PACS: 00.00 20.00 42.10
We propose a new mathematical technique for processing of experimental data measured with large errors. The method is applied to all available experiments on 4 major thermonuclear reactions taken into account in simulations of deuterium and deuterium–tritium fusion targets. We present new approximations for the dependence of the cross sections on energy and for the dependence of the reaction rates on temperature. Along with these approximations, we propose a procedure allowing to estimate their confidence belts. Such estimations were not known before. New approximations provide error ∼0.3% for the cross sections and ∼4% for the reaction rates. The present data are up to ∼5 times more accurate than reported in literature.
Keywords: Thermonuclear reactions Cross sections reaction rates experiments processing
1. Introduction 1.1. Problem Since the mid-1960s, a great deal of interest is paid to the controlled fusion problem [1], [2]. In the projects under construction, deuterium and deuterium-tritium targets are mostly considered. The main points of interest are the ignition conditions [3–8], instabilities occurrence and ways of their suppression (see, e.g., [9–12] and references therein). An important part of these simulations concerns implied reaction kinetics models. The energetic outcome is determined by the following major thermonuclear reactions:
D + D → p + T,
(1)
3
D + D → n + He,
(2)
D + T → n + 4 He,
(3)
3
4
D + He → p + He.
(4)
However, the fusion behavior is also influenced by a number of other inelastic processes: ionization, recombination, charge transfer, excitation, deactivation, etc. Known experimental data on the mentioned processes are collected in large compendia (e.g. [13,14]). Starting from 1950s, experiments have been carried out in order to determine the dependencies of the thermonuclear reaction cross
⁎
sections σ(E) on the energy E in the center-of-mass system. The dependencies of the reaction rates K(T) on temperature T are commonly obtained by integration of the cross sections σ(E) with the Maxwellian distribution. The latter is quite reasonable to be implied because in fusion targets, fuel density increases up to 100 g/cm3. Under such conditions, local thermodynamical equilibrium is rapidly established. In later years, dozens of experimental studies have been performed. In the open-access segment, database [13] contains ∼2000 measurements for the reactions (1)–(4). These experiments are extremely difficult and sometimes possess large errors. For example, Fig. 1 shows all available experimental data on the S-factor of the reaction (1). We would like to remind that S-factor equals to the cross section divided by the Gamow multiplier (i.e., quasi-classical penetrability of the Coulomb barrier [15]). One can see that the results reported in different works diverge significantly, especially at low energies (up to 6 times). However, the required level of accuracy is about several percent. Commonly, the measurements which seem to be unreliable are simply discarded. However, there is considerable ambiguity in such discarding because the totality of all measurements looks like a spread “cloud” without any obvious kernel and it is simply unclear which points should be excluded from consideration. This ambiguity may sufficiently affect the result. Since the reaction rates represent sharp change as temperature varies (i.e., the problem is “stiff” in mathematical sense), slight alteration of approximation parameters in K(T) leads to considerable changes in the shape of K(T). Thereby, the inaccuracies of K(T) determination may spoil the result of overall simulation.
Corresponding author at: Faculty of Physics, Lomonosov Moscow State University, Leninskie Gory 1, bld. 2, Moscow 119991, Russian Federation E-mail addresses: [email protected] (A.A. Belov), [email protected] (N.N. Kalitkin), [email protected] (I.A. Kozlitin).
https://doi.org/10.1016/j.fusengdes.2019.02.082 Received 15 April 2018; Received in revised form 2 February 2019; Accepted 20 February 2019 0920-3796/ © 2019 Elsevier B.V. All rights reserved.
Fusion Engineering and Design 141 (2019) 51–58
A.A. Belov, et al.
validity for cross sections with broad maximum is questionable. Also, (7) ignores the Gamow law at E → 0. In order to provide reasonable approximation both at low energies and near the maximum, Davidenko proposed [24] to take the product of the formulae (5) and (7). Sometimes, the S-factor is approximated by polynomial dependencies, e.g., linear [25] or quadratic [26]. For the reactions (1) and (2), the S-factors have rather simple form, and polynomial expressions permit to obtain a good accuracy within the given interval. However, extrapolation of such formulae beyond the initial boundaries is unreliable. Thereby, their applicability range is limited by the one available in the experiments. Consequently, it is necessary to revise the formulae each time new data appear. Also, to calculate the reaction rates, the polynomials are “stitched” to some exponential or similar law which is selected for each type of process individually [14]. It is rather inconvenient. Duane represented the S˜ -function from (6) via a sum of the BreitWigner formula and a constant term [27]
S˜ (E ) =
Fig. 1. S-factor for the reaction D + D → p + T; points are experimental data [13]; lines stand for different approximations, numbers near the curves indicate the corresponding references; error bars are revised statistical weights.
A + E {A2 + E [A3 + E (A 4 + EA5 )]} S˜ (E ) = 1 . 1 + E {B2 + E [B3 + E (B4 + E )]}
1.2.1. Cross section Let us point out some known formulae for cross section dependencies. Let E be the energy in the center-of-mass system. At low energies, the cross sections are described by the Gamow formula [16,17]
A, B = const.
(5)
However, the domain of applicability for (5) is limited by E ⩽ 30–50 keV. Mott proposed generalization of (5) implying another expression for the Coulomb barrier penetrability [18]
σ (E ) =
S˜ (E ) , E [exp(B / E ) − 1]
S˜ ≈ A = const.
(6)
If the cross section exhibits a maximum, the Breit-Wigner formula [19] is often implied [20–23]
σ (E ) =
Γ12 , (E − E0)2 + Γ22
E0, Γ1,2 = const.
(9)
This formula was later used in Los Alamos [30] for processing of data on the reaction (3). Bosch and Hale implied [28] the Padé approximation (9) for the Sfactor instead of the S˜ -function from (6). For reactions (1) and (2), the denominator was equal 1; i.e., the S-factor was in fact approximated with a polynomial of the fourth order. For resonant reactions (3) and (4), the denominator coefficients Bi ≠ 0; however, the approximations did not describe high energy ranges E ⩾ 600 keV and E ⩾ 1 MeV, respectively. For these energies, special formulae were constructed which were smoothly conjuncted with the low-and-moderate energy approximations. Similar approach was implied by Nocente et al. [33] for the reaction (4). In his work, he used ratio of the fifth order polynomials. This allowed to improve the accuracy of the approximation (namely, in the high energy range). Kozlov proposed formulae containing up to 6 adjusting parameters [34]. The formulae were constructed individually for each reaction and their form was not the same. For the sake of laconism, we do not present them in the text. In recent years, a number of theoretical formulae were proposed (see, e.g., [35–38] and references therein). These approximations were constructed from the quasi-classical penetrability and attempted to account for the nuclear forces. However, accurate expressions for nuclear potential have not been proposed yet. This deteriorates accuracy of such models. For the reaction (1), several approximations proposed earlier are depicted in Figure 1. They agree reasonably at moderate energies, however, when E > 300–500 keV, they diverge from each other and from the measurements rather strongly. Obviously, extrapolation of these approximations beyond the initial energy segment leads to enormous error. The Bosch-Hale formula coincides with the experimental data within graphical accuracy.
1.2. Known approximations
A B ⎫ exp ⎧− , ⎨ E E⎬ ⎩ ⎭
(8)
Note that the exponent B was also treated as a fitting parameter along with A, E0, Γ1, Γ2. Since B did not coincide with its theoretical value, this formula did not permit extrapolation to low energies. This range is crucial for the controlled fusion problem. In [28], it was noted that extrapolation lower even E = 20 keV was almost impossible because it contradicted known experimental measurements. The nest step was performed by Peres who proposed the Padé approximation for the S˜ -function [29]
Considering the cross section formulae, one should pay special attention to the low energy range. The experimental data start from E ⩾ 1.5 keV. The temperatures achieved in laser fusion experiments are up to T ∼ 10–100 eV. Thereby, the really significant energy range is E ⩽ 200–300 eV. This means far extrapolation of the cross section dependencies. To construct approximation which allows such extrapolation is a non-trivial mathematical problem. In the present work, we have performed a refined processing of the cross sections and the reaction rates for the reactions (1)–(4). It provides not only the approximations but also their confidence belt estimations. The obtained accuracy for S(E) is about 0.3% in enormous energy range from 10 eV to ∼10 MeV and, for K(T), it is better than 4% also in large temperature range from 10 eV to 2 MeV. Such accuracy is considerably higher and the energy and temperature ranges are sufficiently wider than presented earlier in literature. However, we are not sure if one can restore the nuclear forces with confidence basing on such approximations. As it is well-known, this problem is mathematically ill-posed and satisfactory solution can be obtained only with additional theoretical suppositions. The latter exceeds the limits of our competence.
σ (E ) =
Γ12 + A. (E − E0)2 − Γ22
1.2.2. Reaction rates Hively [39,40] and Brunelli [41] derived the reaction rate expressions using the Duane formulae (8) for the cross sections. However,
(7)
However, this formula was derived for narrow resonances and its 52
Fusion Engineering and Design 141 (2019) 51–58
A.A. Belov, et al.
Here, δj are errors assigned to experimental arguments xj. Note that to obtain an, one has to solve a system of linear algebraic equations. The above described method behaves well if the values δj are small enough. However, if they are large then the approximation curve “attempts” to reproduce the saw-like behaviour of experimental points and high-frequency oscillations appear on it. In order to fix this issue, we have developed regularization of the double period method which damps high-frequency oscillations [45,46]. The idea is to add a stabilizer
since the latter suffer from poor low energy extrapolation, one cannot expect these reaction rates to be highly accurate at low temperatures. Peres proposed a sophisticated formula for the reaction rate [29]
K (T ) = C1 θ ξ (mc 2T 3)−1 e−3ξ , θ = T ⎡1 − ⎣
ξ = B2/3 (4θ)−1/3, −1
T (C2 + T (C4 + TC6)) ⎤ 1 + T (C3 + T (C5 + TC7)) ⎦
.
(10)
Bosch and Hale used the same functional form for all reactions (1)–(4) and revised the coefficients [28]. This allowed to improve the accuracy considerably. Kozlov performed numerical integration of his formulae for cross sections and approximated the results with some analytical dependencies. They were also constructed individually for each reaction. To conclude the literature survey, we would like to mention that the Bosch–Hale formulae for the S-factors and reactivities are widely used by European and American experimenters and theoreticians (see, e.g., [31,32]). The Kozlov formulae for these values are most utilized by Russian investigators (see, e.g., [2,42] and references therein). Therefore, we compare our results with the mentioned works. This comparison is presented below.
α
2.1. Mathematical technique It is convenient to process the experiments in such variables that provide the most simple form of the curve. These are lgS versus lgE. Sfactor is determined by nuclear forces only. Consequently, S ≈ const for E → 0 [43]. This allows to extrapolate S(E) to low energies. In the experiments, the targets are “cold”, i.e., they contain atoms but not bare nuclei. Strictly speaking, the atomic electrons slightly screen the nucleus. However, the energy of bonded electrons in D and T equals 13.6 eV and, in He, it is 54.5 eV. These values are sufficiently smaller than the energies of the incident projectiles E ⩾ 2 keV. Thereby, the screening can be neglected. As it was mentioned above, the S-factor approximations must permit accurate extrapolation beyond the lower boundary of the energy segment. If experimental data are highly accurate then good extrapolation can be obtained using the double period method [44]. It implies a specific decomposition of non-periodic function into Fourier series in which regular harmonics are added with several (∼3) odd harmonics with twice less frequency. The notation of the method is as follows. For convenience, let us linearly transform the argument E ∈ [Emin, Emax] into x ∈ [−π/2, π/2]
E − Emin , Emax − Emin
−π ≤ x ≤ π .
(16)
to the sum (15). Here, α > 0 is the regularization parameter and ωn stand for the harmonic frequencies. Commonly, terms with ωn0 and ωn2 are implied in the stabilizer. However, these terms provide weaker damping of high-frequency harmonics and also may distort the approximation curve and its first derivative because they introduce penalties for large values of u and u′. Stabilizer with ωn4 does not affect the values of u and u′ and thereby, provides sufficiently better results. The most difficult point is how to select the regularization parameter and the number of the fundamental period harmonics N. Increasing of α provides stronger damping of non-physical oscillations of the approximation curve. However, simultaneously, the residual R increases. Increasing of N decreases the residual but, on the other hand, the round-off errors in the solution of the linear system in an increase rapidly. Thereby, this is an optimization problem with contradictory criteria, and one has to seek a compromise. The procedure of obtaining the compromise solution of the problem for the reaction (1) is presented in Fig. 2. On lgN − lgα plane, we depict isolines of number m of decimal digits lost due to round-off errors in the solution of the linear system in an and isolines of the residual lgR. The result is acceptable if lgR < 0 and m < 810 for computations with 16 decimal digits (i.e., with 64-bit software). From Fig. 2, one can see that for the reaction (1), good results are provided by α = 0.3, N = 50. The reactions (2)–(4) are processed in the similar manner.
2. Cross section approximations
x = −π + 2π
π /2
∫−π /2 [u″ (x )]2 dx ∼ α ∑ an2 ωn4
2.2. Results The computed approximation for reaction D + D → p + T is shown in Fig. 1. One can see that it covers the entire experimental range 2 keV ⩽ E ⩽ 14 MeV and reasonably extrapolates to lower energies E < 2 keV as a constant. For reactions (2)–(4), the results are similar. They are presented in Fig. 3.
(11)
The ordinary Fourier harmonics are called the fundamental period
φ0 = 1, φ1 = sin x , φ2 = cos x , …, φ2N − 1 = sin 2Nx, φ2N = cos 2Nx.
(12)
They are supplied with M odd Fourier harmonics of the doubled period
φ2N + 1 = sin x /2, φ2N + 2 = cos x /2, φ2N + 3 = sin 3x /2, …
(13)
Careful analysis shows that M = 3 is sufficient, so, further, we shall use this value. The approximation takes the form 2N + M
u (x ) ≈
∑
an φn (x ).
(14)
n=0
Fourier coefficients an are obtained by the least squares method
R=
1 J−N
J
Fig. 2. Selection of N and α for the reaction D + D → p + T; R is the relative residual (15), m stands for the number of decimal digits lost because of roundoff error.
2
u (x j ) − Sj ⎞ ⎟ → min. δj ⎝ ⎠
∑ ⎜⎛ j=1
(15) 53
Fusion Engineering and Design 141 (2019) 51–58
A.A. Belov, et al.
2.2.1. Confidence interval The inaccuracy of each experimental point includes the random error of the concrete measurement and the systematic error of the selected laboratory. Let us consider the totality of measurements of all laboratories together, then the systematic error of each of them can be treated as a random value because the number of laboratories is rather large. Therefore, for the entire set of experimental data, all errors can be considered random. For each reaction, a sufficiently large number of points is available (300–700). Consequently, the obtained approximation is a reasonable averaging of errors and possesses sufficient accuracy. Thereby, deviation of the point from the approximation curve can be treated as the total error of this point. According to the law of large numbers, the distribution of this random value is close to the Gaussian one, whatever are the distributions for each single point and each single factor affecting its error. To construct the error estimate, let us calculate all deviations u (xj) − Sj and take them as amplitudes of the corresponding random values. Then, perform their independent statistical sampling and calculate new approximation u1(x) from this instant realization. Next, manifoldly repeat this procedure (∼30 times) and obtain a set of curves {uk(x)}. They are randomly distributed with respect to the first curve. At each x, we can calculate local dispersion and standard deviation of this set of curves. Obtained standards d, % for all reactions are depicted on Fig. 4. The confidence belt is smaller in the middle of the energy interval and slightly increases near its boundaries but even there equals ∼0.3%. Such an excellent accuracy could be achieved because we included a large number of points. According to the laws of statistics, averaging over J points decreases the error ∼ J times.
Fig. 3. S-factor curves from Table 1. Table 1 S-factors lgS, keV mbn. lgE, keV
DD → pT
DD → n3He
DT → n4 He
D3He → p4 He
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.1 3.2 3.3 3.4
4.774 4.773 4.771 4.770 4.769 4.768 4.768 4.767 4.768 4.770 4.773 4.776 4.781 4.789 4.799 4.813 4.830 4.853 4.880 4.913 4.952 4.998 5.048 5.102 5.157 5.212 5.267 5.321 5.373 5.426 5.481 5.540 5.600
4.773 4.773 4.773 4.772 4.769 4.767 4.768 4.772 4.776 4.781 4.786 4.794 4.803 4.812 4.824 4.841 4.862 4.887 4.919 4.958 5.003 5.053 5.106 5.160 5.215 5.272 5.330 5.389 5.448 5.504 5.557 5.607 5.658
7.086 7.086 7.086 7.086 7.090 7.098 7.108 7.120 7.137 7.161 7.196 7.241 7.297 7.360 7.414 7.428 7.368 7.227 7.029 6.807 6.590 6.395 6.226 6.082 5.954 5.844 5.754 5.684 5.635 5.605 5.596 5.605 5.627
6.852 6.852 6.852 6.852 6.852 6.852 6.851 6.848 6.844 6.840 6.836 6.835 6.836 6.842 6.859 6.892 6.940 7.003 7.076 7.153 7.216 7.231 7.161 7.004 6.803 6.610 6.435 6.283 6.156 6.052 5.970 5.908 5.863
2.2.2. Extrapolation to E → 0 According to the Gamow law, the S-factor should be extrapolated as a constant. The curves in Fig. 3 possess almost horizontal slope at low energies, thereby, their constant extrapolation is admissible. We emphasize, that the accuracy of the entire extrapolation to E < Emin equals the accuracy at Emin, i.e. ∼0.3%. 2.3. Comparison with known approximations Fig. 5 shows the discrepancy (in percent) between the present data for the S-factor and the Kozlov formulae [34]. One can see that the latter are applicable for E ⩽ 200 keV. Within this range, the discrepancy for the reaction (1) is 4–13%. For reaction (2), it is up to 18%. For reaction (3), the difference is merely ∼3–5%; while for (4), it is ∼14%.
The number of the Fourier coefficients required for the approximation turned out to be large (up to 150–200 for each reaction). It is unreasonable to list them in the text. Therefore, for the sake of convenience, we present the S-factors in Table 1 with step Δ lgE = 0.1. Note that, in the table, logarithmic scale is implied. It is especially handy for data at high energies. Earlier, thorough mathematical techniques for accuracy estimates of the obtained curves were not reported. In the present paper, we propose a method of constructing such estimates. Fig. 4. Standard deviations of the S-factor curves from Fig. 1. 54
Fusion Engineering and Design 141 (2019) 51–58
A.A. Belov, et al.
and Bosch-Hale formulae in the entire energy range 10 eV ⩽ E ⩽ 8 MeV. 2.4. Experimental errors revision An important part of work concerned selection of the statistical weights in the approximation. They should be equal to experimental errors. Quantitative error values given in the original contributions are quite moderate. However, the discrepancy between data reported by different authors is rather large and often exceeds their accuracy estimates. Consequently, it is risky to take them straightforwardly. Thereby, the error estimates should be revised. We have analysed the experimental setups and have revised the statistical weights. The errors are assigned not only to laboratories in general, but also to individual points. The error values for lgS varied from 0.004 to 0.78 (i.e. from 1% to 600%). The primary factors taken into account are the following. 2.4.1. Random errors If data of a particular experimental work follow some smooth line, the role of random errors is not large. But in a number of works, one can observe saw-like arrangement of points which indicates considerable random errors. Mostly, such behaviour takes place at low energies. The largest jitter of the curve occurs for heavy metal targets. Obviously, the random error is approximately equal to the amplitude of this “saw”.
Fig. 5. Divergence of the Kozlov formulae [34] for S(E) from Table 1.
At high energies E ⩾ 200 keV, the discrepancy rapidly increases up to 100%, i.e., the Kozlov formulae loose applicability. In Fig. 6, one can see similar comparison of the present data with the Bosch-Hale approximations [28] within the validity boundaries of the latter. For the reaction (1), the Bosch–Hale formula are valid in a very broad range 0.5 keV ⩽ E ⩽ 5 MeV. The difference between this formula and our data is ∼4–6%. For the reaction (2), the validity range of the Bosch–Hale formulae is also very broad 0.5 keV ⩽E ⩽ 4.9 MeV and the discrepancy between them and the present data is ∼3–6%. For the reaction (3), the Bosch–Hale formulae are valid for 0.5 ⩽ E ⩽ 550 keV and the divergence is ∼4–7%. Finally, for the reaction (4), the Bosch–Hale formulae are applicable in the range 0.5 ⩽ E ⩽ 900 keV. The difference is up to 20% at low energies E ⩽ 20 keV; for moderate and high energies 20 ⩽ E ⩽ 900 keV, it is ∼10%. The discrepancies between the present data and the Bosch–Hale formulae agree reasonably with the accuracy estimates of the R-matrix calculations outlined in [28]. The accuracy estimations of our formulae are ∼0.3%. This value is one standard deviation, it corresponds to the confidence probability 0.68. If one takes two standard deviations (i.e., the confidence probability 0.95) the error estimations should be multiplied by the factor 2. These estimations are obtained via mathematically rigorous procedure. Thereby, the present data sufficiently improve accuracy of the Kozlov
2.4.2. Systematic errors They are caused by deceleration of the ion beam in the targets. The experimenters introduce corrections to these effects; however, calculation of these corrections is complicated and is not reliable quite often. Therefore, the systematic error is, in fact, a subjective accuracy estimation of the implied corrections. We have performed such estimates via analysis of the target type and the processes in the target. If the target is filled with deuterium gas or thin deuterium ice (up to several nm) then the results are the most accurate. If the deuterium ice is thick, one should introduce corrections on absorption in the target ice so accuracy decreases. If a gas target is covered with windows made of thin foils, the deceleration in them sufficiently increases the correction and the inaccuracy of the result. On the contrary, these corrections are small for thin films made of CD2 (deuterated polyethylene). The largest error is given by heavy metals saturated by deuterium (Zn, Ti, Ta, W) because there is strong deceleration of deuterium nuclei on the highly energetic inner electrons of the metals. We do not list the assigned error values due to large amount of the material. For the reaction (1), the revised statistical weights for several points are shown in Fig. 1. The largest error corresponds to deuterium saturated Pb target. On the contrary, cryogenically pumped gas and CD4 targets provide rather small errors for both high and low energies. Naturally, our estimates are quite subjective and require verification. 2.4.3. Verification of the errors There is a simple way to perform a plausible verification of the assigned errors. Let us compare the deviation of the experimental points from the obtained regression with the statistical weights. The data from world-recognized research centers (Los Alamos, Oak Ridge, Sarov, Aberdeen, etc.) sustain such verification. The approximation curve passes through the error bars. To obtain quantitative estimation for the whole data array, let us consider the relative residual (15). If the weights δj agree with the deviation of the corresponding points from the approximation curve then the relative residual is ∼1. If the selected weights are too small (i.e., the assigned errors are too large), this value is small compared with 1. In contrast, if selected weights are too large (i.e., the assigned errors are too small), R is sufficiently greater than 1. In our case, for reaction (1), R is between 1 and ∼0.7. So, the assigned weights are fairly close to the real experimental inaccuracies. For the reactions (2)–(4), the results are similar.
Fig. 6. Divergence of the Bosch–Hale formulae [28] for S(E) from Table 1. 55
Fusion Engineering and Design 141 (2019) 51–58
A.A. Belov, et al.
Table 3 Approximation coefficients of lgK(T) with Fourier series (18).
ξ0 ξ1 η1 ξ2 η2 ξ3 η3 ξ4 η4 η5 ξ6 η6
DD → pT
DD → n3He
DT → n4 He
D3He → p4 He
−41.898 30.017 24.313 −9.150 −16.995 2.140 7.572 −0.316 −2.540 0.580 0.012 −0.068
−42.145 30.484 24.459 −9.451 −17.170 2.265 7.744 −0.346 −2.655 0.630 0.015 −0.079
−42.559 33.312 24.394 −9.876 −16.909 2.501 7.406 −0.414 −2.148 0.411 0 0
−60.504 53.039 43.244 −16.233 −31.497 3.980 15.707 −0.749 −6.055 1.734 0.004 −0.319
Fig. 7. Reaction rate curves (18). Note that rates of the reactions D + D → p + T and D + D → n + 3He almost coincide.
3. Reaction rates 3.1. Calculations The reaction rates are obtained by averaging of σv over the Maxwellian distribution
K (T ) =
∫0
∞
σ (E ) fMaxw (v )v dv, E =
Mv 2 . 2
(17)
We performed the integration numerically. The obtained curves are presented in Fig. 7; note that the ones for the reactions (1) and (2) almost coincide. In Table 2, we present the data on lgK(T) on a rare mesh with step Δ lgT = 0.5. These results are valid in knowingly sufficient range 10 keV ⩽ T ⩽ 2 MeV. For convenience, the tables of K(T) were approximated by the double period method with the following formula
lg K (T ) = t=
6 ∑k = 0
π (lg T 5.30
Fig. 8. Standard deviations of the reaction rate curves from Fig. 7.
statistics to be reasonable, such variation is performed 30 times for each reaction. Next, at each T, we calculate dispersion and standard deviation of the obtained K(T) set (i.e. the confidence belt). The values of the confidence belts (in percent) are presented in Fig. 8. One can see that at low temperatures, the error is maximal. For the reactions (1)–(4), it equals from 1% to 4%. As temperature increases T > 1 keV, all errors quickly decrease, and at T > 3 keV, they do not exceed 1%. Such accuracy is quite sufficient for gas-dynamic codes. We would like to emphasize that such estimations for the reaction rates and the cross sections were not known before. The achieved accuracy is sufficient for fusion target simulations. Note that these values are also one standard deviation (confidence 0.68).
(ξk cos kt + ηk sin kt),
− 3.65), 1.0 ⩽ lg T ⩽ 6.3.
(18)
The approximation (18) is given in physical units and temperature is in keV. For all 4 reactions, η0 = 0, ξ5 = 0; the rest coefficients are listed in Table 3. 3.2. Error estimates The confidence belts for K(T) are calculated in the same manner as mentioned above. For each random variation of the “experimental” points, S(E) is approximated and K(T) is calculated via it. For the
3.3. Comparison with previously reported formulae Table 2 Reaction rates lgK, cm3 s−1. lgT, keV
DD → pT
DD → n He
DT → n He
D He → p He
−2.0 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0
−50.402 −38.650 −30.750 −25.476 −21.988 −19.708 −18.237 −17.283 −16.635 −16.188 −15.888
−50.402 −38.651 −30.750 −25.475 −21.985 −19.695 −18.213 −17.246 −16.582 −16.124 −15.826
−50.492 −37.984 −29.567 −23.936 −20.177 −17.643 −15.947 −15.162 −15.082 −15.318 −15.589
−74.292 −54.241 −40.687 −31.559 −25.454 −21.396 −18.644 −16.736 −15.762 −15.552 −15.662
3
4
3
Fig. 9 shows divergence of the Kozlov formulae for reactivities from our data in percent for all reactions. At low temperatures T < 100 keV, the inaccuracy of the Kozlov formulae reaches ∼10–20% (whereas the error of the proposed approximations are ∼4%, i.e. ∼5 times less). In the vicinity of T = 100 keV, the difference decreases significantly; however, at T > 100 keV, the Kozlov formulae rapidly loose applicability. Fig. 10 presents similar comparison of the present data with the Bosch–Hale formulae. The applicability boundaries of these formulae are T = 0.2–100 keV for the reactions (1)–(3) and T = 0.5–190 keV for the reaction (4) [28]. For the reactions (1) and (3), the discrepancy between the present data and the Bosch–Hale formula is ∼3%, which approximately equals the accuracy of our data. For the reaction (2), the
4
56
Fusion Engineering and Design 141 (2019) 51–58
A.A. Belov, et al.
and Hale [28] is statistically insignificant from the point of view of the work [28] within the limits of their accuracy. From point of view of our approximations (providing higher accuracy) the discrepancy is statistically significant. To conclude, the present approximations improve the previously proposed formulae and provide accuracy and applicability range which exceed the requirements of contemporary applied calculations. 4. Conclusion In the present paper, we propose a new mathematical method for processing of experimental curves measured with large errors. For the first time, a procedure of the confidence interval estimation is proposed. These methods are implied for processing of the whole accessible totality of experimental data on 4 major thermonuclear reactions responsible for the energetic outcome in fusion reactors. The cross sections of these reactions are calculated with accuracy exceeding 0.3% within energy range 10 eV < E < 10 MeV. The reaction rates are computed for 10 eV < T < 2000 keV with accuracy 1–4%. The obtained results sufficiently improve the ones reported earlier. Currently, in the Keldysh Institute of Applied Mathematics of Russian Academy of Sciences, the group including the authors of the present paper constructs the TEFIS database [47] on thermophysical properties of matter under extreme conditions. We intend to supply it with detailed tables of the cross sections and the reaction rates for the considered reactions. The demo version of the tables can be accessed at [48]. Also, we have developed a program package in the Matlab/Octave language which realizes the regularized double period method. The package is open-access [49].
Fig. 9. Divergence of the Kozlov formulae [34] for K(T) from the proposed approximations (18).
Conflict of interests None. Acknowledgements This work was supported by Russian Science Foundation grant, Project No. 16-11-10001. We would like to express our gratitude to anonymous reviewer for valuable critical comments which promoted improvement of the paper. Fig. 10. Divergence of the Bosch-Hale formulae [28] for K(T) from the proposed approximations (18).
References
discrepancy is ∼7% for T ⩽ 100 keV, i.e., 2 times bigger than the accuracy of our data. Finally, for the reaction (4), the divergence if up to ∼1for T ⩽ 190 keV. These values exceed the error estimations outlined in [28]. The explanation is as follows. The cross section is obtained by multiplying the S-factor found with some error by the Gamow factor which is known exactly. To evaluate the reactivity, one averages σv over the Maxwellian distribution which is also exact. Multiplication and averaging can easily be performed via high-precision mathematical procedures. Therefore, the accuracy of the reactivity cannot be better than the error of the S-factor in the energy range from ∼E0/(2–3) to ∼(2–3) E0 where E0 is the maximum of the Maxwellian distribution. Consider, for instance, the reaction (4). In the energy range E ⩽ 600 keV, the authors of [28] state that the cross section accuracy equals ∼10% (since this range includes the resonant maximum; in its vicinity the uncertainty is maximal). Consequently, in the temperature range T < 200 keV, they can guarantee the accuracy of the K(T) expressions also equal ∼10% (instead of ∼2% as stated in [28]). The ∼2% estimate pointed out in [28] can be the auxiliary error caused by the mathematical procedure performed in this work. For the rest reactions, the arguments are similar. Thereby, the deviation of the present formulae from those of Bosch
[1] V.B. Rozanov, Phys.-Usp 47 (2004) 359. [2] V.B. Rozanov, Quantum Electron 27 (1997) 1063. [3] F. Mitarai, et al., Fusion Eng. Des. 136 A (2018) 82, https://doi.org/10.1016/j. fusengdes.2017.12.037. [4] S. Atzeni, et al., Nucl. Fusion 54 (2014) 054008. [5] H. Takabe, et al., Fusion Eng. Des. 44 (1999) 105. [6] A.I. Mahdy, et al., Fusion Eng. Des. 44 (1999) 255. [7] A.J. Kemp, et al., Nucl. Fusion 41 (2001) 235. [8] M.M. Basko, Nucl. Fusion 30 (1990) 2443. [9] P. Maget, et al., Nucl. Fusion 56 (2016) 086004. [10] V. Toigo, et al., Fusion Eng. Des. 86 (2011) 565. [11] S. Gordeev, et al., Fusion Eng. Des. 83 (2008) 1524. [12] T.C. Hender, et al., Nucl. Fusion 47 (2007) 128. [13] NEA Data Bank – Nuclear Data Services. http://www.oecd-nea.org/janisweb/ search/exfor. [14] EIRENE Atomic & Molecular Database. http://www.eirene.de/html/a-m_data.html. [15] A.A. Gamow, Phys.-Usp 10 (1930) 531. [16] W.R. Arnold, et al., Phys. Rev. 93 (1954) 483. [17] W. Wenzel, W.A. Whaling, Phys. Rev. 88 (1952) 1149. [18] N.F. Mott, H.S.W. Massey, The Theory of Atomic Collisions, (1949) Oxford. [19] G. Breit, E. Wigner, Phys. Rev. 49 (1936) 519. [20] J.P. Conner, et al., Phys. Rev. 88 (1952) 468. [21] T.W. Bonner, et al., Phys. Rev. 88 (1952) 473. [22] N. Jarmie, et al., Phys. Rev. C 29 (1984) 2031. [23] R.E. Brown, et al., Phys. Rev. C 35 (1987) 1999. [24] V.A. Davidenko, et al., Soviet J. Atomic Energy 2 (1957) 474. [25] R.E. Brown, N. Jarmie, Phys. Rev. C 41 (1990) 1391. [26] A. Krauss, et al., Nucl. Phys. A 465 (1987) 150. [27] B.H. Duane, Annual Report on CTR Technology 1972 (Wolkenhauer W C, ed.), Rep.
57
Fusion Engineering and Design 141 (2019) 51–58
A.A. Belov, et al.
[41] B. Brunelli, Nuovo Cim 55 (1980) 264. [42] G.V. Dolgoleva, E.A. Zabrodina, KIAM Preprint (2014) 2014-68. [43] L.D. Landau, E.M. Lifshitz, Quantum Mechanics (Non-relativistic Theory), 3rd edition, Butterworth-Heinemann, Oxford, 1981. [44] N.N. Kalitkin, L.V. Kuzmina, Dokl. Math. 62 (2000) 306. [45] A.A. Belov, N.N. Kalitkin, Dokl. Math. 94 (2016) 602. [46] A.A. Belov, N.N. Kalitkin, Comp. Math. Math. Phys. 57 (2017) 1741. [47] N.N. Kalitkin, et al., TEFIS Database, Keldysh IAM RAS, Moscow, http://tefis.ru, http://tefis.keldysh.ru. [48] A.A. Belov, N.N. Kalitkin, RaCheTheR Database on the Rates Of Chemical And Thermonuclear Reactions, Keldysh IAM RAS, Lomonosov MSU, Moscow, https:// github.com/ABelov91/RaCheTheR. [49] A.A. Belov, RDPM Package for the Regularized Double Period Method, Lomonosov MSU, Moscow, https://github.com/ABelov91/RMDP.
BNWL-1685, Battelle Pacific Northwest Laboratory, Richland, WA, 1972. [28] H.S. Bosch, G.M. Hale, Nucl. Fusion 32 (1992) 611. [29] A.J. Peres, J. Nucl. Mater. 50 (1979) 5569. [30] G. Sadler, P. Van Belle, Tech. Rep. JET-IR(87)08, Jet Joint Undertaking, Abingdon, Oxfordshire, 1987. [31] A.R. Christopherson, et al., Phys. Plasmas 25 (2018) 012703. [32] W. Biel, et al., Fusion Eng. Des. 123 (2017) 206. [33] M. Nocente, et al., Nucl. Fusion 50 (2010) 055001. [34] B.N. Kozlov, Soviet J. Atomic Energy 12 (1962) 247. [35] T. Koohrokhi, et al., J. Fusion Energy 35 (2016) 816. [36] V.I. Kukulin, et al., J. Phys. G: Nucl. Phys. 10 (1984) L213. [37] V.T. Voronchev, et al., Mem. Fac. Eng. Kyushu Univ. 50 (1990) 517. [38] V.T. Voronchev, Y. Nakao, J. Phys. G: Nucl. Phys. 29 (2003) 431. [39] L.M. Hively, Nucl. Fusion 17 (1977) 873. [40] L.M. Hively, Nucl. Technol. Fusion 3 (1983) 199.
58
Contents lists available at ScienceDirect
Fusion Engineering and Design journal homepage: www.elsevier.com/locate/fusengdes
Refinement of thermonuclear reaction rates A.A. Belov a b c
a,b,⁎
c
, N.N. Kalitkin , I.A. Kozlitin
T
c
Faculty of Physics, Lomonosov Moscow State University, Leninskie Gory 1, bld. 2, Moscow 119991, Russian Federation Peoples’ Friendship University of Russia (RUDN University), Miklukho-Maklaya St., 6, Moscow 117198, Russian Federation Keldysh Institute of Applied Mathematics of Russian Academy of Sciences, Miusskaya Sq., 4, 125047 Moscow, Russian Federation
A R T I C LE I N FO
A B S T R A C T
PACS: 00.00 20.00 42.10
We propose a new mathematical technique for processing of experimental data measured with large errors. The method is applied to all available experiments on 4 major thermonuclear reactions taken into account in simulations of deuterium and deuterium–tritium fusion targets. We present new approximations for the dependence of the cross sections on energy and for the dependence of the reaction rates on temperature. Along with these approximations, we propose a procedure allowing to estimate their confidence belts. Such estimations were not known before. New approximations provide error ∼0.3% for the cross sections and ∼4% for the reaction rates. The present data are up to ∼5 times more accurate than reported in literature.
Keywords: Thermonuclear reactions Cross sections reaction rates experiments processing
1. Introduction 1.1. Problem Since the mid-1960s, a great deal of interest is paid to the controlled fusion problem [1], [2]. In the projects under construction, deuterium and deuterium-tritium targets are mostly considered. The main points of interest are the ignition conditions [3–8], instabilities occurrence and ways of their suppression (see, e.g., [9–12] and references therein). An important part of these simulations concerns implied reaction kinetics models. The energetic outcome is determined by the following major thermonuclear reactions:
D + D → p + T,
(1)
3
D + D → n + He,
(2)
D + T → n + 4 He,
(3)
3
4
D + He → p + He.
(4)
However, the fusion behavior is also influenced by a number of other inelastic processes: ionization, recombination, charge transfer, excitation, deactivation, etc. Known experimental data on the mentioned processes are collected in large compendia (e.g. [13,14]). Starting from 1950s, experiments have been carried out in order to determine the dependencies of the thermonuclear reaction cross
⁎
sections σ(E) on the energy E in the center-of-mass system. The dependencies of the reaction rates K(T) on temperature T are commonly obtained by integration of the cross sections σ(E) with the Maxwellian distribution. The latter is quite reasonable to be implied because in fusion targets, fuel density increases up to 100 g/cm3. Under such conditions, local thermodynamical equilibrium is rapidly established. In later years, dozens of experimental studies have been performed. In the open-access segment, database [13] contains ∼2000 measurements for the reactions (1)–(4). These experiments are extremely difficult and sometimes possess large errors. For example, Fig. 1 shows all available experimental data on the S-factor of the reaction (1). We would like to remind that S-factor equals to the cross section divided by the Gamow multiplier (i.e., quasi-classical penetrability of the Coulomb barrier [15]). One can see that the results reported in different works diverge significantly, especially at low energies (up to 6 times). However, the required level of accuracy is about several percent. Commonly, the measurements which seem to be unreliable are simply discarded. However, there is considerable ambiguity in such discarding because the totality of all measurements looks like a spread “cloud” without any obvious kernel and it is simply unclear which points should be excluded from consideration. This ambiguity may sufficiently affect the result. Since the reaction rates represent sharp change as temperature varies (i.e., the problem is “stiff” in mathematical sense), slight alteration of approximation parameters in K(T) leads to considerable changes in the shape of K(T). Thereby, the inaccuracies of K(T) determination may spoil the result of overall simulation.
Corresponding author at: Faculty of Physics, Lomonosov Moscow State University, Leninskie Gory 1, bld. 2, Moscow 119991, Russian Federation E-mail addresses: [email protected] (A.A. Belov), [email protected] (N.N. Kalitkin), [email protected] (I.A. Kozlitin).
https://doi.org/10.1016/j.fusengdes.2019.02.082 Received 15 April 2018; Received in revised form 2 February 2019; Accepted 20 February 2019 0920-3796/ © 2019 Elsevier B.V. All rights reserved.
Fusion Engineering and Design 141 (2019) 51–58
A.A. Belov, et al.
validity for cross sections with broad maximum is questionable. Also, (7) ignores the Gamow law at E → 0. In order to provide reasonable approximation both at low energies and near the maximum, Davidenko proposed [24] to take the product of the formulae (5) and (7). Sometimes, the S-factor is approximated by polynomial dependencies, e.g., linear [25] or quadratic [26]. For the reactions (1) and (2), the S-factors have rather simple form, and polynomial expressions permit to obtain a good accuracy within the given interval. However, extrapolation of such formulae beyond the initial boundaries is unreliable. Thereby, their applicability range is limited by the one available in the experiments. Consequently, it is necessary to revise the formulae each time new data appear. Also, to calculate the reaction rates, the polynomials are “stitched” to some exponential or similar law which is selected for each type of process individually [14]. It is rather inconvenient. Duane represented the S˜ -function from (6) via a sum of the BreitWigner formula and a constant term [27]
S˜ (E ) =
Fig. 1. S-factor for the reaction D + D → p + T; points are experimental data [13]; lines stand for different approximations, numbers near the curves indicate the corresponding references; error bars are revised statistical weights.
A + E {A2 + E [A3 + E (A 4 + EA5 )]} S˜ (E ) = 1 . 1 + E {B2 + E [B3 + E (B4 + E )]}
1.2.1. Cross section Let us point out some known formulae for cross section dependencies. Let E be the energy in the center-of-mass system. At low energies, the cross sections are described by the Gamow formula [16,17]
A, B = const.
(5)
However, the domain of applicability for (5) is limited by E ⩽ 30–50 keV. Mott proposed generalization of (5) implying another expression for the Coulomb barrier penetrability [18]
σ (E ) =
S˜ (E ) , E [exp(B / E ) − 1]
S˜ ≈ A = const.
(6)
If the cross section exhibits a maximum, the Breit-Wigner formula [19] is often implied [20–23]
σ (E ) =
Γ12 , (E − E0)2 + Γ22
E0, Γ1,2 = const.
(9)
This formula was later used in Los Alamos [30] for processing of data on the reaction (3). Bosch and Hale implied [28] the Padé approximation (9) for the Sfactor instead of the S˜ -function from (6). For reactions (1) and (2), the denominator was equal 1; i.e., the S-factor was in fact approximated with a polynomial of the fourth order. For resonant reactions (3) and (4), the denominator coefficients Bi ≠ 0; however, the approximations did not describe high energy ranges E ⩾ 600 keV and E ⩾ 1 MeV, respectively. For these energies, special formulae were constructed which were smoothly conjuncted with the low-and-moderate energy approximations. Similar approach was implied by Nocente et al. [33] for the reaction (4). In his work, he used ratio of the fifth order polynomials. This allowed to improve the accuracy of the approximation (namely, in the high energy range). Kozlov proposed formulae containing up to 6 adjusting parameters [34]. The formulae were constructed individually for each reaction and their form was not the same. For the sake of laconism, we do not present them in the text. In recent years, a number of theoretical formulae were proposed (see, e.g., [35–38] and references therein). These approximations were constructed from the quasi-classical penetrability and attempted to account for the nuclear forces. However, accurate expressions for nuclear potential have not been proposed yet. This deteriorates accuracy of such models. For the reaction (1), several approximations proposed earlier are depicted in Figure 1. They agree reasonably at moderate energies, however, when E > 300–500 keV, they diverge from each other and from the measurements rather strongly. Obviously, extrapolation of these approximations beyond the initial energy segment leads to enormous error. The Bosch-Hale formula coincides with the experimental data within graphical accuracy.
1.2. Known approximations
A B ⎫ exp ⎧− , ⎨ E E⎬ ⎩ ⎭
(8)
Note that the exponent B was also treated as a fitting parameter along with A, E0, Γ1, Γ2. Since B did not coincide with its theoretical value, this formula did not permit extrapolation to low energies. This range is crucial for the controlled fusion problem. In [28], it was noted that extrapolation lower even E = 20 keV was almost impossible because it contradicted known experimental measurements. The nest step was performed by Peres who proposed the Padé approximation for the S˜ -function [29]
Considering the cross section formulae, one should pay special attention to the low energy range. The experimental data start from E ⩾ 1.5 keV. The temperatures achieved in laser fusion experiments are up to T ∼ 10–100 eV. Thereby, the really significant energy range is E ⩽ 200–300 eV. This means far extrapolation of the cross section dependencies. To construct approximation which allows such extrapolation is a non-trivial mathematical problem. In the present work, we have performed a refined processing of the cross sections and the reaction rates for the reactions (1)–(4). It provides not only the approximations but also their confidence belt estimations. The obtained accuracy for S(E) is about 0.3% in enormous energy range from 10 eV to ∼10 MeV and, for K(T), it is better than 4% also in large temperature range from 10 eV to 2 MeV. Such accuracy is considerably higher and the energy and temperature ranges are sufficiently wider than presented earlier in literature. However, we are not sure if one can restore the nuclear forces with confidence basing on such approximations. As it is well-known, this problem is mathematically ill-posed and satisfactory solution can be obtained only with additional theoretical suppositions. The latter exceeds the limits of our competence.
σ (E ) =
Γ12 + A. (E − E0)2 − Γ22
1.2.2. Reaction rates Hively [39,40] and Brunelli [41] derived the reaction rate expressions using the Duane formulae (8) for the cross sections. However,
(7)
However, this formula was derived for narrow resonances and its 52
Fusion Engineering and Design 141 (2019) 51–58
A.A. Belov, et al.
Here, δj are errors assigned to experimental arguments xj. Note that to obtain an, one has to solve a system of linear algebraic equations. The above described method behaves well if the values δj are small enough. However, if they are large then the approximation curve “attempts” to reproduce the saw-like behaviour of experimental points and high-frequency oscillations appear on it. In order to fix this issue, we have developed regularization of the double period method which damps high-frequency oscillations [45,46]. The idea is to add a stabilizer
since the latter suffer from poor low energy extrapolation, one cannot expect these reaction rates to be highly accurate at low temperatures. Peres proposed a sophisticated formula for the reaction rate [29]
K (T ) = C1 θ ξ (mc 2T 3)−1 e−3ξ , θ = T ⎡1 − ⎣
ξ = B2/3 (4θ)−1/3, −1
T (C2 + T (C4 + TC6)) ⎤ 1 + T (C3 + T (C5 + TC7)) ⎦
.
(10)
Bosch and Hale used the same functional form for all reactions (1)–(4) and revised the coefficients [28]. This allowed to improve the accuracy considerably. Kozlov performed numerical integration of his formulae for cross sections and approximated the results with some analytical dependencies. They were also constructed individually for each reaction. To conclude the literature survey, we would like to mention that the Bosch–Hale formulae for the S-factors and reactivities are widely used by European and American experimenters and theoreticians (see, e.g., [31,32]). The Kozlov formulae for these values are most utilized by Russian investigators (see, e.g., [2,42] and references therein). Therefore, we compare our results with the mentioned works. This comparison is presented below.
α
2.1. Mathematical technique It is convenient to process the experiments in such variables that provide the most simple form of the curve. These are lgS versus lgE. Sfactor is determined by nuclear forces only. Consequently, S ≈ const for E → 0 [43]. This allows to extrapolate S(E) to low energies. In the experiments, the targets are “cold”, i.e., they contain atoms but not bare nuclei. Strictly speaking, the atomic electrons slightly screen the nucleus. However, the energy of bonded electrons in D and T equals 13.6 eV and, in He, it is 54.5 eV. These values are sufficiently smaller than the energies of the incident projectiles E ⩾ 2 keV. Thereby, the screening can be neglected. As it was mentioned above, the S-factor approximations must permit accurate extrapolation beyond the lower boundary of the energy segment. If experimental data are highly accurate then good extrapolation can be obtained using the double period method [44]. It implies a specific decomposition of non-periodic function into Fourier series in which regular harmonics are added with several (∼3) odd harmonics with twice less frequency. The notation of the method is as follows. For convenience, let us linearly transform the argument E ∈ [Emin, Emax] into x ∈ [−π/2, π/2]
E − Emin , Emax − Emin
−π ≤ x ≤ π .
(16)
to the sum (15). Here, α > 0 is the regularization parameter and ωn stand for the harmonic frequencies. Commonly, terms with ωn0 and ωn2 are implied in the stabilizer. However, these terms provide weaker damping of high-frequency harmonics and also may distort the approximation curve and its first derivative because they introduce penalties for large values of u and u′. Stabilizer with ωn4 does not affect the values of u and u′ and thereby, provides sufficiently better results. The most difficult point is how to select the regularization parameter and the number of the fundamental period harmonics N. Increasing of α provides stronger damping of non-physical oscillations of the approximation curve. However, simultaneously, the residual R increases. Increasing of N decreases the residual but, on the other hand, the round-off errors in the solution of the linear system in an increase rapidly. Thereby, this is an optimization problem with contradictory criteria, and one has to seek a compromise. The procedure of obtaining the compromise solution of the problem for the reaction (1) is presented in Fig. 2. On lgN − lgα plane, we depict isolines of number m of decimal digits lost due to round-off errors in the solution of the linear system in an and isolines of the residual lgR. The result is acceptable if lgR < 0 and m < 810 for computations with 16 decimal digits (i.e., with 64-bit software). From Fig. 2, one can see that for the reaction (1), good results are provided by α = 0.3, N = 50. The reactions (2)–(4) are processed in the similar manner.
2. Cross section approximations
x = −π + 2π
π /2
∫−π /2 [u″ (x )]2 dx ∼ α ∑ an2 ωn4
2.2. Results The computed approximation for reaction D + D → p + T is shown in Fig. 1. One can see that it covers the entire experimental range 2 keV ⩽ E ⩽ 14 MeV and reasonably extrapolates to lower energies E < 2 keV as a constant. For reactions (2)–(4), the results are similar. They are presented in Fig. 3.
(11)
The ordinary Fourier harmonics are called the fundamental period
φ0 = 1, φ1 = sin x , φ2 = cos x , …, φ2N − 1 = sin 2Nx, φ2N = cos 2Nx.
(12)
They are supplied with M odd Fourier harmonics of the doubled period
φ2N + 1 = sin x /2, φ2N + 2 = cos x /2, φ2N + 3 = sin 3x /2, …
(13)
Careful analysis shows that M = 3 is sufficient, so, further, we shall use this value. The approximation takes the form 2N + M
u (x ) ≈
∑
an φn (x ).
(14)
n=0
Fourier coefficients an are obtained by the least squares method
R=
1 J−N
J
Fig. 2. Selection of N and α for the reaction D + D → p + T; R is the relative residual (15), m stands for the number of decimal digits lost because of roundoff error.
2
u (x j ) − Sj ⎞ ⎟ → min. δj ⎝ ⎠
∑ ⎜⎛ j=1
(15) 53
Fusion Engineering and Design 141 (2019) 51–58
A.A. Belov, et al.
2.2.1. Confidence interval The inaccuracy of each experimental point includes the random error of the concrete measurement and the systematic error of the selected laboratory. Let us consider the totality of measurements of all laboratories together, then the systematic error of each of them can be treated as a random value because the number of laboratories is rather large. Therefore, for the entire set of experimental data, all errors can be considered random. For each reaction, a sufficiently large number of points is available (300–700). Consequently, the obtained approximation is a reasonable averaging of errors and possesses sufficient accuracy. Thereby, deviation of the point from the approximation curve can be treated as the total error of this point. According to the law of large numbers, the distribution of this random value is close to the Gaussian one, whatever are the distributions for each single point and each single factor affecting its error. To construct the error estimate, let us calculate all deviations u (xj) − Sj and take them as amplitudes of the corresponding random values. Then, perform their independent statistical sampling and calculate new approximation u1(x) from this instant realization. Next, manifoldly repeat this procedure (∼30 times) and obtain a set of curves {uk(x)}. They are randomly distributed with respect to the first curve. At each x, we can calculate local dispersion and standard deviation of this set of curves. Obtained standards d, % for all reactions are depicted on Fig. 4. The confidence belt is smaller in the middle of the energy interval and slightly increases near its boundaries but even there equals ∼0.3%. Such an excellent accuracy could be achieved because we included a large number of points. According to the laws of statistics, averaging over J points decreases the error ∼ J times.
Fig. 3. S-factor curves from Table 1. Table 1 S-factors lgS, keV mbn. lgE, keV
DD → pT
DD → n3He
DT → n4 He
D3He → p4 He
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.1 3.2 3.3 3.4
4.774 4.773 4.771 4.770 4.769 4.768 4.768 4.767 4.768 4.770 4.773 4.776 4.781 4.789 4.799 4.813 4.830 4.853 4.880 4.913 4.952 4.998 5.048 5.102 5.157 5.212 5.267 5.321 5.373 5.426 5.481 5.540 5.600
4.773 4.773 4.773 4.772 4.769 4.767 4.768 4.772 4.776 4.781 4.786 4.794 4.803 4.812 4.824 4.841 4.862 4.887 4.919 4.958 5.003 5.053 5.106 5.160 5.215 5.272 5.330 5.389 5.448 5.504 5.557 5.607 5.658
7.086 7.086 7.086 7.086 7.090 7.098 7.108 7.120 7.137 7.161 7.196 7.241 7.297 7.360 7.414 7.428 7.368 7.227 7.029 6.807 6.590 6.395 6.226 6.082 5.954 5.844 5.754 5.684 5.635 5.605 5.596 5.605 5.627
6.852 6.852 6.852 6.852 6.852 6.852 6.851 6.848 6.844 6.840 6.836 6.835 6.836 6.842 6.859 6.892 6.940 7.003 7.076 7.153 7.216 7.231 7.161 7.004 6.803 6.610 6.435 6.283 6.156 6.052 5.970 5.908 5.863
2.2.2. Extrapolation to E → 0 According to the Gamow law, the S-factor should be extrapolated as a constant. The curves in Fig. 3 possess almost horizontal slope at low energies, thereby, their constant extrapolation is admissible. We emphasize, that the accuracy of the entire extrapolation to E < Emin equals the accuracy at Emin, i.e. ∼0.3%. 2.3. Comparison with known approximations Fig. 5 shows the discrepancy (in percent) between the present data for the S-factor and the Kozlov formulae [34]. One can see that the latter are applicable for E ⩽ 200 keV. Within this range, the discrepancy for the reaction (1) is 4–13%. For reaction (2), it is up to 18%. For reaction (3), the difference is merely ∼3–5%; while for (4), it is ∼14%.
The number of the Fourier coefficients required for the approximation turned out to be large (up to 150–200 for each reaction). It is unreasonable to list them in the text. Therefore, for the sake of convenience, we present the S-factors in Table 1 with step Δ lgE = 0.1. Note that, in the table, logarithmic scale is implied. It is especially handy for data at high energies. Earlier, thorough mathematical techniques for accuracy estimates of the obtained curves were not reported. In the present paper, we propose a method of constructing such estimates. Fig. 4. Standard deviations of the S-factor curves from Fig. 1. 54
Fusion Engineering and Design 141 (2019) 51–58
A.A. Belov, et al.
and Bosch-Hale formulae in the entire energy range 10 eV ⩽ E ⩽ 8 MeV. 2.4. Experimental errors revision An important part of work concerned selection of the statistical weights in the approximation. They should be equal to experimental errors. Quantitative error values given in the original contributions are quite moderate. However, the discrepancy between data reported by different authors is rather large and often exceeds their accuracy estimates. Consequently, it is risky to take them straightforwardly. Thereby, the error estimates should be revised. We have analysed the experimental setups and have revised the statistical weights. The errors are assigned not only to laboratories in general, but also to individual points. The error values for lgS varied from 0.004 to 0.78 (i.e. from 1% to 600%). The primary factors taken into account are the following. 2.4.1. Random errors If data of a particular experimental work follow some smooth line, the role of random errors is not large. But in a number of works, one can observe saw-like arrangement of points which indicates considerable random errors. Mostly, such behaviour takes place at low energies. The largest jitter of the curve occurs for heavy metal targets. Obviously, the random error is approximately equal to the amplitude of this “saw”.
Fig. 5. Divergence of the Kozlov formulae [34] for S(E) from Table 1.
At high energies E ⩾ 200 keV, the discrepancy rapidly increases up to 100%, i.e., the Kozlov formulae loose applicability. In Fig. 6, one can see similar comparison of the present data with the Bosch-Hale approximations [28] within the validity boundaries of the latter. For the reaction (1), the Bosch–Hale formula are valid in a very broad range 0.5 keV ⩽ E ⩽ 5 MeV. The difference between this formula and our data is ∼4–6%. For the reaction (2), the validity range of the Bosch–Hale formulae is also very broad 0.5 keV ⩽E ⩽ 4.9 MeV and the discrepancy between them and the present data is ∼3–6%. For the reaction (3), the Bosch–Hale formulae are valid for 0.5 ⩽ E ⩽ 550 keV and the divergence is ∼4–7%. Finally, for the reaction (4), the Bosch–Hale formulae are applicable in the range 0.5 ⩽ E ⩽ 900 keV. The difference is up to 20% at low energies E ⩽ 20 keV; for moderate and high energies 20 ⩽ E ⩽ 900 keV, it is ∼10%. The discrepancies between the present data and the Bosch–Hale formulae agree reasonably with the accuracy estimates of the R-matrix calculations outlined in [28]. The accuracy estimations of our formulae are ∼0.3%. This value is one standard deviation, it corresponds to the confidence probability 0.68. If one takes two standard deviations (i.e., the confidence probability 0.95) the error estimations should be multiplied by the factor 2. These estimations are obtained via mathematically rigorous procedure. Thereby, the present data sufficiently improve accuracy of the Kozlov
2.4.2. Systematic errors They are caused by deceleration of the ion beam in the targets. The experimenters introduce corrections to these effects; however, calculation of these corrections is complicated and is not reliable quite often. Therefore, the systematic error is, in fact, a subjective accuracy estimation of the implied corrections. We have performed such estimates via analysis of the target type and the processes in the target. If the target is filled with deuterium gas or thin deuterium ice (up to several nm) then the results are the most accurate. If the deuterium ice is thick, one should introduce corrections on absorption in the target ice so accuracy decreases. If a gas target is covered with windows made of thin foils, the deceleration in them sufficiently increases the correction and the inaccuracy of the result. On the contrary, these corrections are small for thin films made of CD2 (deuterated polyethylene). The largest error is given by heavy metals saturated by deuterium (Zn, Ti, Ta, W) because there is strong deceleration of deuterium nuclei on the highly energetic inner electrons of the metals. We do not list the assigned error values due to large amount of the material. For the reaction (1), the revised statistical weights for several points are shown in Fig. 1. The largest error corresponds to deuterium saturated Pb target. On the contrary, cryogenically pumped gas and CD4 targets provide rather small errors for both high and low energies. Naturally, our estimates are quite subjective and require verification. 2.4.3. Verification of the errors There is a simple way to perform a plausible verification of the assigned errors. Let us compare the deviation of the experimental points from the obtained regression with the statistical weights. The data from world-recognized research centers (Los Alamos, Oak Ridge, Sarov, Aberdeen, etc.) sustain such verification. The approximation curve passes through the error bars. To obtain quantitative estimation for the whole data array, let us consider the relative residual (15). If the weights δj agree with the deviation of the corresponding points from the approximation curve then the relative residual is ∼1. If the selected weights are too small (i.e., the assigned errors are too large), this value is small compared with 1. In contrast, if selected weights are too large (i.e., the assigned errors are too small), R is sufficiently greater than 1. In our case, for reaction (1), R is between 1 and ∼0.7. So, the assigned weights are fairly close to the real experimental inaccuracies. For the reactions (2)–(4), the results are similar.
Fig. 6. Divergence of the Bosch–Hale formulae [28] for S(E) from Table 1. 55
Fusion Engineering and Design 141 (2019) 51–58
A.A. Belov, et al.
Table 3 Approximation coefficients of lgK(T) with Fourier series (18).
ξ0 ξ1 η1 ξ2 η2 ξ3 η3 ξ4 η4 η5 ξ6 η6
DD → pT
DD → n3He
DT → n4 He
D3He → p4 He
−41.898 30.017 24.313 −9.150 −16.995 2.140 7.572 −0.316 −2.540 0.580 0.012 −0.068
−42.145 30.484 24.459 −9.451 −17.170 2.265 7.744 −0.346 −2.655 0.630 0.015 −0.079
−42.559 33.312 24.394 −9.876 −16.909 2.501 7.406 −0.414 −2.148 0.411 0 0
−60.504 53.039 43.244 −16.233 −31.497 3.980 15.707 −0.749 −6.055 1.734 0.004 −0.319
Fig. 7. Reaction rate curves (18). Note that rates of the reactions D + D → p + T and D + D → n + 3He almost coincide.
3. Reaction rates 3.1. Calculations The reaction rates are obtained by averaging of σv over the Maxwellian distribution
K (T ) =
∫0
∞
σ (E ) fMaxw (v )v dv, E =
Mv 2 . 2
(17)
We performed the integration numerically. The obtained curves are presented in Fig. 7; note that the ones for the reactions (1) and (2) almost coincide. In Table 2, we present the data on lgK(T) on a rare mesh with step Δ lgT = 0.5. These results are valid in knowingly sufficient range 10 keV ⩽ T ⩽ 2 MeV. For convenience, the tables of K(T) were approximated by the double period method with the following formula
lg K (T ) = t=
6 ∑k = 0
π (lg T 5.30
Fig. 8. Standard deviations of the reaction rate curves from Fig. 7.
statistics to be reasonable, such variation is performed 30 times for each reaction. Next, at each T, we calculate dispersion and standard deviation of the obtained K(T) set (i.e. the confidence belt). The values of the confidence belts (in percent) are presented in Fig. 8. One can see that at low temperatures, the error is maximal. For the reactions (1)–(4), it equals from 1% to 4%. As temperature increases T > 1 keV, all errors quickly decrease, and at T > 3 keV, they do not exceed 1%. Such accuracy is quite sufficient for gas-dynamic codes. We would like to emphasize that such estimations for the reaction rates and the cross sections were not known before. The achieved accuracy is sufficient for fusion target simulations. Note that these values are also one standard deviation (confidence 0.68).
(ξk cos kt + ηk sin kt),
− 3.65), 1.0 ⩽ lg T ⩽ 6.3.
(18)
The approximation (18) is given in physical units and temperature is in keV. For all 4 reactions, η0 = 0, ξ5 = 0; the rest coefficients are listed in Table 3. 3.2. Error estimates The confidence belts for K(T) are calculated in the same manner as mentioned above. For each random variation of the “experimental” points, S(E) is approximated and K(T) is calculated via it. For the
3.3. Comparison with previously reported formulae Table 2 Reaction rates lgK, cm3 s−1. lgT, keV
DD → pT
DD → n He
DT → n He
D He → p He
−2.0 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0
−50.402 −38.650 −30.750 −25.476 −21.988 −19.708 −18.237 −17.283 −16.635 −16.188 −15.888
−50.402 −38.651 −30.750 −25.475 −21.985 −19.695 −18.213 −17.246 −16.582 −16.124 −15.826
−50.492 −37.984 −29.567 −23.936 −20.177 −17.643 −15.947 −15.162 −15.082 −15.318 −15.589
−74.292 −54.241 −40.687 −31.559 −25.454 −21.396 −18.644 −16.736 −15.762 −15.552 −15.662
3
4
3
Fig. 9 shows divergence of the Kozlov formulae for reactivities from our data in percent for all reactions. At low temperatures T < 100 keV, the inaccuracy of the Kozlov formulae reaches ∼10–20% (whereas the error of the proposed approximations are ∼4%, i.e. ∼5 times less). In the vicinity of T = 100 keV, the difference decreases significantly; however, at T > 100 keV, the Kozlov formulae rapidly loose applicability. Fig. 10 presents similar comparison of the present data with the Bosch–Hale formulae. The applicability boundaries of these formulae are T = 0.2–100 keV for the reactions (1)–(3) and T = 0.5–190 keV for the reaction (4) [28]. For the reactions (1) and (3), the discrepancy between the present data and the Bosch–Hale formula is ∼3%, which approximately equals the accuracy of our data. For the reaction (2), the
4
56
Fusion Engineering and Design 141 (2019) 51–58
A.A. Belov, et al.
and Hale [28] is statistically insignificant from the point of view of the work [28] within the limits of their accuracy. From point of view of our approximations (providing higher accuracy) the discrepancy is statistically significant. To conclude, the present approximations improve the previously proposed formulae and provide accuracy and applicability range which exceed the requirements of contemporary applied calculations. 4. Conclusion In the present paper, we propose a new mathematical method for processing of experimental curves measured with large errors. For the first time, a procedure of the confidence interval estimation is proposed. These methods are implied for processing of the whole accessible totality of experimental data on 4 major thermonuclear reactions responsible for the energetic outcome in fusion reactors. The cross sections of these reactions are calculated with accuracy exceeding 0.3% within energy range 10 eV < E < 10 MeV. The reaction rates are computed for 10 eV < T < 2000 keV with accuracy 1–4%. The obtained results sufficiently improve the ones reported earlier. Currently, in the Keldysh Institute of Applied Mathematics of Russian Academy of Sciences, the group including the authors of the present paper constructs the TEFIS database [47] on thermophysical properties of matter under extreme conditions. We intend to supply it with detailed tables of the cross sections and the reaction rates for the considered reactions. The demo version of the tables can be accessed at [48]. Also, we have developed a program package in the Matlab/Octave language which realizes the regularized double period method. The package is open-access [49].
Fig. 9. Divergence of the Kozlov formulae [34] for K(T) from the proposed approximations (18).
Conflict of interests None. Acknowledgements This work was supported by Russian Science Foundation grant, Project No. 16-11-10001. We would like to express our gratitude to anonymous reviewer for valuable critical comments which promoted improvement of the paper. Fig. 10. Divergence of the Bosch-Hale formulae [28] for K(T) from the proposed approximations (18).
References
discrepancy is ∼7% for T ⩽ 100 keV, i.e., 2 times bigger than the accuracy of our data. Finally, for the reaction (4), the divergence if up to ∼1for T ⩽ 190 keV. These values exceed the error estimations outlined in [28]. The explanation is as follows. The cross section is obtained by multiplying the S-factor found with some error by the Gamow factor which is known exactly. To evaluate the reactivity, one averages σv over the Maxwellian distribution which is also exact. Multiplication and averaging can easily be performed via high-precision mathematical procedures. Therefore, the accuracy of the reactivity cannot be better than the error of the S-factor in the energy range from ∼E0/(2–3) to ∼(2–3) E0 where E0 is the maximum of the Maxwellian distribution. Consider, for instance, the reaction (4). In the energy range E ⩽ 600 keV, the authors of [28] state that the cross section accuracy equals ∼10% (since this range includes the resonant maximum; in its vicinity the uncertainty is maximal). Consequently, in the temperature range T < 200 keV, they can guarantee the accuracy of the K(T) expressions also equal ∼10% (instead of ∼2% as stated in [28]). The ∼2% estimate pointed out in [28] can be the auxiliary error caused by the mathematical procedure performed in this work. For the rest reactions, the arguments are similar. Thereby, the deviation of the present formulae from those of Bosch
[1] V.B. Rozanov, Phys.-Usp 47 (2004) 359. [2] V.B. Rozanov, Quantum Electron 27 (1997) 1063. [3] F. Mitarai, et al., Fusion Eng. Des. 136 A (2018) 82, https://doi.org/10.1016/j. fusengdes.2017.12.037. [4] S. Atzeni, et al., Nucl. Fusion 54 (2014) 054008. [5] H. Takabe, et al., Fusion Eng. Des. 44 (1999) 105. [6] A.I. Mahdy, et al., Fusion Eng. Des. 44 (1999) 255. [7] A.J. Kemp, et al., Nucl. Fusion 41 (2001) 235. [8] M.M. Basko, Nucl. Fusion 30 (1990) 2443. [9] P. Maget, et al., Nucl. Fusion 56 (2016) 086004. [10] V. Toigo, et al., Fusion Eng. Des. 86 (2011) 565. [11] S. Gordeev, et al., Fusion Eng. Des. 83 (2008) 1524. [12] T.C. Hender, et al., Nucl. Fusion 47 (2007) 128. [13] NEA Data Bank – Nuclear Data Services. http://www.oecd-nea.org/janisweb/ search/exfor. [14] EIRENE Atomic & Molecular Database. http://www.eirene.de/html/a-m_data.html. [15] A.A. Gamow, Phys.-Usp 10 (1930) 531. [16] W.R. Arnold, et al., Phys. Rev. 93 (1954) 483. [17] W. Wenzel, W.A. Whaling, Phys. Rev. 88 (1952) 1149. [18] N.F. Mott, H.S.W. Massey, The Theory of Atomic Collisions, (1949) Oxford. [19] G. Breit, E. Wigner, Phys. Rev. 49 (1936) 519. [20] J.P. Conner, et al., Phys. Rev. 88 (1952) 468. [21] T.W. Bonner, et al., Phys. Rev. 88 (1952) 473. [22] N. Jarmie, et al., Phys. Rev. C 29 (1984) 2031. [23] R.E. Brown, et al., Phys. Rev. C 35 (1987) 1999. [24] V.A. Davidenko, et al., Soviet J. Atomic Energy 2 (1957) 474. [25] R.E. Brown, N. Jarmie, Phys. Rev. C 41 (1990) 1391. [26] A. Krauss, et al., Nucl. Phys. A 465 (1987) 150. [27] B.H. Duane, Annual Report on CTR Technology 1972 (Wolkenhauer W C, ed.), Rep.
57
Fusion Engineering and Design 141 (2019) 51–58
A.A. Belov, et al.
[41] B. Brunelli, Nuovo Cim 55 (1980) 264. [42] G.V. Dolgoleva, E.A. Zabrodina, KIAM Preprint (2014) 2014-68. [43] L.D. Landau, E.M. Lifshitz, Quantum Mechanics (Non-relativistic Theory), 3rd edition, Butterworth-Heinemann, Oxford, 1981. [44] N.N. Kalitkin, L.V. Kuzmina, Dokl. Math. 62 (2000) 306. [45] A.A. Belov, N.N. Kalitkin, Dokl. Math. 94 (2016) 602. [46] A.A. Belov, N.N. Kalitkin, Comp. Math. Math. Phys. 57 (2017) 1741. [47] N.N. Kalitkin, et al., TEFIS Database, Keldysh IAM RAS, Moscow, http://tefis.ru, http://tefis.keldysh.ru. [48] A.A. Belov, N.N. Kalitkin, RaCheTheR Database on the Rates Of Chemical And Thermonuclear Reactions, Keldysh IAM RAS, Lomonosov MSU, Moscow, https:// github.com/ABelov91/RaCheTheR. [49] A.A. Belov, RDPM Package for the Regularized Double Period Method, Lomonosov MSU, Moscow, https://github.com/ABelov91/RMDP.
BNWL-1685, Battelle Pacific Northwest Laboratory, Richland, WA, 1972. [28] H.S. Bosch, G.M. Hale, Nucl. Fusion 32 (1992) 611. [29] A.J. Peres, J. Nucl. Mater. 50 (1979) 5569. [30] G. Sadler, P. Van Belle, Tech. Rep. JET-IR(87)08, Jet Joint Undertaking, Abingdon, Oxfordshire, 1987. [31] A.R. Christopherson, et al., Phys. Plasmas 25 (2018) 012703. [32] W. Biel, et al., Fusion Eng. Des. 123 (2017) 206. [33] M. Nocente, et al., Nucl. Fusion 50 (2010) 055001. [34] B.N. Kozlov, Soviet J. Atomic Energy 12 (1962) 247. [35] T. Koohrokhi, et al., J. Fusion Energy 35 (2016) 816. [36] V.I. Kukulin, et al., J. Phys. G: Nucl. Phys. 10 (1984) L213. [37] V.T. Voronchev, et al., Mem. Fac. Eng. Kyushu Univ. 50 (1990) 517. [38] V.T. Voronchev, Y. Nakao, J. Phys. G: Nucl. Phys. 29 (2003) 431. [39] L.M. Hively, Nucl. Fusion 17 (1977) 873. [40] L.M. Hively, Nucl. Technol. Fusion 3 (1983) 199.
58