Topic: How to fit astrophysical reaction rates
How to prepare fits of astrophysical reaction rates:
This is a brief outline of how to proceed, including some general warnings of what can go wrong. Some of the remarks were written with fitting rates to the 7-parameter format used in REACLIB in mind. However, they should also apply to fitting any other format.
Part A: General remarks
1) For details on the definition of rates and the 7-parameter REACLIB format, see ADNDT 75, 1 (2000).
2) It is essential to fit rates in the direction of POSITIVE Q-value unless absolutely no information is available for that reaction. Otherwise, any fit errors will be enhanced when computing the rate in the other direction by detailed balance. The fit parameters for the reverse direction can easily be computed from the relations given in Eq. (17) of ADNDT 75, 1 (2000). The fit parameters for both forward and backward reaction have to be included in REACLIB. It is also essential that both sets are derived from the same fit of ONE direction and application of Eq. (17) for the other direction! Separate fits of forward and backward reactions can NEVER be used, they are not consistent and may lead to numerically unstable network equations! This applies for all reaction networks, not just those using REACLIB.
3) Usually, fits are only valid in some temperature region. In REACLIB, this is T9=0.01-10. Nevertheless, it is extremely important that the fits behave well also outside this range. They must not diverge and must not show unphysical oscillation patterns (as sometimes obtained with polynomial fits). Specifically, charged particle rates should go to Zero for T9->0. Neutron capture reactions should become constant or go to zero for T9->0, depending on whether they are s-wave dominated or not.
Part B: Creating the required data set
1) Make sure the rate values (either from experiment or calculations) cover at least the temperature range of validity (see above) for the fit. Obviously, cross section data have to be converted to rates first. For the conversion and to check whether the covered temperature range is sufficient, my program exp2rate.f can be used (see also "How to compute astrophysical reaction rates from experimental data" in this forum). For interpolation (to obtain a more narrow energy/temperature grid for improved fitting) use, e.g. my small program interpolate.f (see also "How to obtain values at arbitrary energies/temperatures" in this forum).
2) If the obtained temperature range is not sufficient, the data have to be supplemented to cover the missing part(s). This has to be performed with extreme caution. Knowledge about reaction mechanisms, additional information, and a lot of experience are usually necessary in order to make an educated guess about the behavior of the cross section and rate. Renormalization of theoretical predictions (both for the Hauser-Feshbach and direct mechanism), inclusion of resonance information, simple extrapolations: Any combination of these may be necessary. Further details are beyond the scope of this description.
3) If distinct resonances are contributing, it may be helpful to split the data sets into the different contributions or to prepare analytic functions for the resonance contributions (e.g., as provided in Caughlan & Fowler, ADNDT 40 (1988) 283). However, limit the number of separate entries to as few as possible (see also C.8 below).
Part C: Fitting the data set
1) In the following, "data" refers to all sets of reaction rates prepared as in Part B, regardless of whether they stem from actual experimental data or theoretical calculations.
2) Usually, satisfactory fits of the astrophysical rates can only be achieved with fit routines that allow for different weights of the data points. When writing your own fitting program, allow for automatic setting of the weights according to some or all of the following rules. If the program cannot compute the weights automatically, iterations with modified weights have to be run manually.
3) The fit criterion should always include relative deviations of the rates as these can vary across many orders of magnitudes, i.e. (abs(fit-rate)/fit)**2.
4) A first run can use equal weights at all temperatures (w=1.0) or consider the experimental errors (as, e.g., included when having used the code exp2rate.f).
5) After the run, check the quality of the fit by inspecting the average deviation (or other fit criterion) AND by testing the behavior of the fit outside the fitted temperature range.
6) The truly complicated game starts when this first attempt does not yield good results. Then it becomes necessary to play around with the weights. Two cases should definitely get lower weighting:
a) Very slow rates, i.e. rates below 1.e-18 or 1.e-20 or so. Even at high density, those rates would be too slow to change anything within the time scale of astrophysical burning (or the lifetime of the Universe). Therefore it is not necessary to reach any accuracy in such a case and deviations even of several orders of magnitude are acceptable! (Unless they make the rate too fast.) The reason why such rates should not be discarded completely is because they serve as constraints to suppress divergences in the fit. It is possible to allow weights that are exponentially suppressed when the rates become smaller, e.g. w=(rate/1.e-18)**0.03. For the same reason, whenever comparing fits to data, it is of no importance when there are large deviations for such slow rates!!
b) Rates at high temperatures, i.e. rates at about T9=5 and above. At high temperatures, all rates attain equilibrium (with the exception of neutrino rates and decays). The rates then cancel out from the equilibrium abundance distribution. The abundances are just determined from spins and Q-values. Nevertheless, the rates have to be included in the network to compute the proper abundances from the relation of forward and backward reactions and to follow the freeze-out (or heat-up) properly. Due to the latter, the fit accuracy should be higher than for the slow rates, e.g. w=0.1.
7) At the same time (with every new fit and change of weights) it has to be assured that the fit still behaves properly outside the fitted range and does not diverge. Divergences and oscillations again may be treated by modified weights. Another trick is to supply additional data points at very low and very high temperature and assign low weights to them. The values of these additional data points are obviously only crude extrapolations or educated guesses but because of their low weights, they should not change the fit in the valid region too much. They merely serve as constraints to prevent divergence.
a) Neutron capture: Additional values at low and high energy can be obtained by assuming s-waves (constant reaction rates) or p-waves (rate proportional to T).
b) Charged particles in the entrance or exit channel: For positive Q-values, the behavior will be dominated by the entrance channel because of the higher energy of the ejectile in the exit channel. Then the temperature dependence of the rate can be estimated from the penetration through the Coulomb barrier. However, when the projectile is a neutron, the Coulomb barrier has to be accounted for in the exit channel at the proper energy (neutron energy + Q-value). Barrier penetrations can either be estimated from a simple WKB approximation or by employing a formula similar to what was used in the fits by Woosley et al. in ADNDT 18 (1976) 305 and 22 (1978) 321.
8) It may prove difficult to fit single resonance contributions. In that case, the fit may be split up in two or more functions for the different contributions (see B.3) which are fitted separately. In REACLIB, additional entries (with the same 7-parameter format) appearing for the same reaction are simply added up. Therefore, make sure that the total rate can be constructed from the sum of the different contributions. Usually, Hauser-Feshbach rates (obtained from calculations in the compound nucleus reaction model) are easy to fit and do not require such splitting. Generally, the number of REACLIB entries should be as few as possible because it is computationally expensive to evaluate a large number of entries. For large reaction networks it is preferrable to have only one or (at most) two entries, even if this leads to a lower fit accuracy. The true art of fitting is to get the best compromise between accuracy and the least number of separate entries for a rate! The only exception is the fitting of nonresonant and resonant contributions to a rate; these should always be fitted separately because they are treated differently when applying screening effects. Ideally, even in this case one limits the fits to one entry each for the resonant and the nonresonant part.
9) After modifying weights and/or supplying additional data points, the fit has to be redone (go to C.5 above). The above procedure can be automatized to a certain extent. However, manual interference may be necessary for some cases.
D. Partition functions
1) Finally, partition functions (see Eq. (12) in ADNDT 75 (2000) 1) have to be supplied for relating forward and backward reactions. They are not included in the fits and have to be supplied separately. For REACLIB, they are included in the file WINVN. It includes partition functions up to T9=10. For higher temperatures, you may check out Ap. J. Suppl. 147 (2003) 403.
Good luck! 