Uncertainties of optical-model parameters for the study of the threshold anomaly

In the analysis of elastic-scattering experimental data, optical-model parameters (usually, depths of real and imaginary potentials) are fitted and conclusions are drawn analyzing their variations at bombardment energies close to the Coulomb barrier (threshold anomaly). The judgement about the shape of this variation (related to the physical processes producing this anomaly) depends on these fitted values but the robustness of the conclusions strongly depends on the uncertainties with which these parameters are derived. We will show that previous published studies have not used a common criterium for the evaluation of the parameter uncertainties. In this work, a study of these uncertainties is presented, using conventional statistic tools as well as bootstrapping techniques. As case studies, these procedures are applied to re-analyze detailed elastic-scattering data for the $^{12}$C + $^{208}$Pb and the $^6$Li + $^{80}$Se systems.


Introduction
In the study of elastic scattering of atomic nuclei at low energies there has been a long-lasting interest in the so-called threshold anomaly (TA [1]). For the case of tightly bound nuclei, this phenomenon, related to the closure of reaction channels, consists in a decrease of the depth of imaginary part of the optical potential at bombarding energies below the Coulomb barrier V B . Due to the causality-related dispersion relation (DR [1]) linking the real and imaginary parts of optical potentials, the depth of the real part also varies strongly, peaking around V B (see Refs. [2,3,1]). In the case of weakly bound projectiles, the coupling to nonelastic channels (e.g. breakup) generates a repulsive polarization potential [4] that can produce either the absence of TA (no TA [5,6,7,8]) or the so-called breakup threshold anomaly (BTA [7,9,8]), in which the imaginary potential increases while approaching V B and, conversely, there is a reduction on the real part of the potential [9,10].
For an unambiguous determination of the kind of anomaly, the variation of the opticalmodel parameter values as a function of the energy must be evaluated taking into account their uncertainties (hereafter called parameter uncertainties), which must be derived from the experimental angular distribution data and their own uncertainties (hereafter called experimental uncertainties). However, publications on the subject apply different, sometimes non-rigorous criteria for the estimation of the parameter uncertainties or do not consider uncertainties at all. Here, a study of the parameter uncertainties is presented, using conventional statistic tools as well as the bootstrap technique [11]. It is hoped that these methods will help to achieve a better determination of the parameters of interest in the characterization of the TA. In section 2 we discuss mathematical and physical issues that influence the determination of uncertainties when adjusting optical-model parameters using experimental angular distributions, in section 3 the 12 C + 208 Pb and the 6 Li + 80 Se systems will be reanalyzed and their uncertainties evaluated with different techniques. Finally, in section 4 we make our final remarks.

The determination of optical-model parameter uncertainties
In this section we will firstly consider the χ 2 test (2.1) since it is the most popular method to evaluate the optical-model parameter uncertainties, we will review the different χ 2 criteria found in the literature (2.2) and finally we will consider the two techniques that we will apply to the case studies in the next section, namely: covariance ellipses (2.3) and bootstrap method (2.4).

The χ 2 test
Let us consider a series of experimental pairs x i , y i (x i ), where the y i values have an uncertainty ∆y i (statistical, background substraction, etc.) which is assumed to be normally distributed. To compare these experimental data with a theoretical model that adjust the data with a function y = f (x) taking into account the experimental uncertainties, we could use an adjusting program that minimizes the value of χ 2 defined as where N is the number of experimental points and R i are the so-called residuals. The adjusting program should minimize the effective distance between the experimental results y i and their theoretical estimation f (x i ). The simplest case to study is the repetition of a measurement of a physical quantity (say a life-time or the speed of light) by several laboratory groups. The theoretical model f (x i ) could be as simple as the construction of the best evaluated value of that physical quantity, e.g. the arithmetic or weighted average. In this simple case the expected χ 2 should be close to N − 1, where N is the number of measurements. If not, procedures such as the proposed by Birge [12,13,14] should be applied. In its original form, it applies for interlaboratory evaluations (external uncertainty) when χ 2 /(N − 1) > 1 and consists in multiplying the experimental (internal) uncertainties by χ 2 /(N − 1). With the increased uncertainties a new, normalized value χ 2′ is obtained, here using N instead of N − 1 for large number of points: We consider now cases in which the theoretical function f has k parameters α = (α 1 , ..., α k ), denoted from now on as f (x i |α), that should be adjusted to achieve the best fit of f to experimental data. It is a well-known rule [15,11] that the uncertainty with which a single, uncorrelated parameter α j is adjusted, can be obtained varying this parameter around its optimum value α j0 (while all others parameters are optimized) until the value of χ 2 increments in one unit, i.e.: Here χ 2 0 denotes the minimum value of χ 2 . For this rule to be valid the following conditions must apply: • The uncertainties of raw experimental data ∆y i are properly determined.
• There are no significant systematic errors. • The theoretical model is a true description of the data being studied.
As consequence of the previous conditions the residuals R i will follow a normal distribution with mean value R i = 0 and variance σ 2 R = 1 (we will denote this R i ∼ N (0, 1)) which implies χ 2 0 ∼ ν (although the reciprocal is not necessarily true). Here ν is the number of degrees of freedom ν = N − n f being n f the number of free parameters. In most of the cases n f = 2 and N ≥ 20, thus ν ≈ N and no significant difference arises from the use of N instead of ν.
As we will see in section 2.2 many data sets analyzed in the area of low-energy elastic scattering do not fulfil all the conditions previously mentioned, but the rule of eq. 3 is nevertheless applied. Hence the uncertainty of the obtained parameters is underestimated. For these cases, we will extend the procedure of Birge (eq. 2) for the experimental adjustment of one or several parameters α of a theoretical function f (x i |α) in the following way: For the optimum set of parameters α 0 we obtain Now, with the experimental uncertainties ∆y i scaled by the Birge factor χ 2 0 /ν the rule of eq. 3 can be applied: This does not imply that the experimental uncertainties ∆y i are effectively changed in fact the original ones should be reported. To avoid confusions we prefer to replace in eq. 6 with eq. 4 and 5 to obtain the equivalent rule In the case of adjusting several (even correlated) parameters the rule of eq. 1 should be extended as χ 2 = χ 2 0 + ∆χ 2 , where for two parameters and a 68% confidence level ∆χ 2 = 2.3. See Ref. [11] (p. 815) for a table of ∆χ 2 for different confidence levels and number of adjustable parameters.
In the more general case in which χ 2 0 /ν > 1 and several parameters are adjusted simultaneously eq. 7 is extended to In the study of scattering angular distributions, the experimental data (x i , y i ) consist in angles and differential cross-sections (θ i , dσ i /dΩ). The theoretical model f (x i |α) can be either the calculated dσ(θ i )/dΩ using the phenomenological Woods-Saxon optical model or the microscopic models involving folding of nuclear densities with appropriated nucleon-nucleon potentials. In this work the São Paulo global microscopic potential (SPP [16,17]) is used with α 1 = N R and α 2 = N I (depth of real and imaginary parts of the potential, respectively). Table 1 presents a non-exhaustive list of experiments studying the TA, either with weakly or tightly bound projectiles. It can be seen that most studies use the standard rule χ 2 ≤ χ 2 0 + 1 (eq. 3) to determine the uncertainties regardless of their χ 2 /ν value . We will show in the next section that due to parameter correlations, this rule gives too small uncertainties bars, even in the best case when χ 2 0 /ν ≈ 1, being even worse when χ 2 /ν ≫ 1. It is apparent that this problem has been recognized but not clearly acknowledged and, thus, many of the different criteria listed in Table 1 attempted to increase the χ 2 ≤ χ 2 0 + 1 limits. A good example of this can be seen in Ref. [18]. Since their angular distribution data for the 12 C + 208 Pb system have an average value χ 2 0 /ν = 6, the usual criterium χ 2 ≤ χ 2 0 + 1 would have clearly been an underestimation of the parameter uncertainties and an extreme criterium χ 2 ≤ 2χ 2 0 was used, which, as we will show in section 3.1, is an overestimation. In Fig. 1 we show the χ 2 vs. N R curve corresponding to our reanalysis using SPP of the 84.9 MeV angular distribution of Ref. [18]. The optimum fit is obtained applying a combination of grid, and steepest descent searching routines and gives: N R0 = 0.5738, N I0 = 0.5173, χ 2 0 = 198, χ 2 0 /ν = 9.9. The covariance matrix, including the variance of each of the two parameters N R and N I and their covariance, is calculated using second derivatives calculated numerically around the minimum [15]. The parabolic approximation of the curve describing χ 2 as a function of N R is shown in Fig. 1 with dashed lines. Its expression is given by

Criteria used in previous works
where σ 2 1 is the variance of N R . This variance determines the value of |N R − N R0 | that produces a change in one unit in the χ 2 , as in eq. 10a. The solid lines in Fig. 1 represent the values corresponding to a grid on N R while N I is optimized. It is seen that the parabolic approximation is quite good even at several σ 1 away from the minimum. Uncertainty limits corresponding to different criteria found in the literature are: These equations are not meant to represent different confidence limits (CL) but different criteria.
In the particular example of Fig. 1, σ 3 represents the 68% CL for the adjustment of the single parameter N R , considered as uncorrelated with N I , while σ 4 represents the 68% CL for the simultaneous adjustment of the correlated N R and N I . In this case, since χ 2 0 /ν = 9.9, σ 1 and σ 2 do not represent 68% CL. In the parabolic approximation of eq. 9 it holds Perhaps, the great variety of criteria in Table 1 are due to the fact that in most cases the standard rule χ 2 ≤ χ 2 0 + 1 produces too small uncertainties. This turns to be evident when the optimum values of the parameters are in effect changed by their uncertainties: the corresponding change in the angular distributions is often indistinguishable from the optimum ones.

Error ellipses
One interesting example within low-energy nuclear physics of the use of error ellipses to estimate uncertainties for correlated parameters is presented in Ref. [39]. However, this has not been applied to take into account the correlation between the optical-model parameters. In fact, none of the works listed in Table 1 uses this technique. In our case the ellipses corresponding to the different limits previously mentioned are calculated as  where ∆χ 2′ is defined in eq.8, the covariance matrix C is calculated for each system and energy and α − α 0 represent deviations respect the optimum value of the parameters. Figure 2 shows three of such ellipses for the 12 C + 208 Pb angular distribution at 84.9 MeV corresponding to the σ 1 , σ 3 and σ 4 limits for the determination of the standard uncertainty in both parameters. The limit σ 2 should be used in cases where (χ 2 /ν) ≈ 1. In the cases were (χ 2 /ν) > 1, σ 4 should be used. We will show the error ellipses for all energies of both cases studied in Figs. 3 and 5.

Bootstrap
Bootstrap is one of the many resampling methods designed to go beyond regular statistic test. It creates a number N B of synthetic data sets, each of them consisting on N data points selected randomly (with reposition) from the original data set of N points (it is desirable to have N B ≫ N ). On the average a fraction 1/e (about 37%) of the N points will be repeated ones (thus having more weight in the fitting procedure) and, consequently, 37% of the elements will not be included. Each synthetic data set is used to fit the parameters of interest. The standard deviation of the N B values obtained for each parameter is a reliable estimation of its uncertainty. Hence, this procedure gauges the sensitivity of the fitted parameter to each individual data point by the simulation of N B experiments in which the experimenter could have chosen to repeat the measurement of some of the data points at the cost of missing some others. The bootstrap method has been one of the techniques applied to the evaluation of the half-life of 198 Au in Ref. [40]. To the best of our knowledge, bootstrap techniques have not been applied yet to the calculation of nuclear potential parameters and their uncertainties.

Case studies
In the first step to study the 12 C + 208 Pb and the 6 Li + 80 Se systems, the optimum N R , N I values are calculated as well as the corresponding residuals, covariance matrices and error ellipses. Since a condition of the analysis is the normality of the residuals, the Kolmogorov-Smirnov test has been applied. This statistical test evaluates the maximum difference between the empirical and reference (theoretical) cumulative distribution functions and yields the acceptance or rejection of a null hypothesis at a given confidence level [41]. In our case, we compare to two different null hypothesis: i) the residuals follows a normal distribution with mean value R = 0 and variance σ 2 R = 1, i.e. H 01 : R i ∼ N (0, 1) and ii) the residuals follows a normal distribution but in this case with variance σ 2 3.1. The 12 C + 208 Pb system As expected, the residuals of the angular distributions of each individual energy, as well as the residuals corresponding to all data taken together, follow N (0, χ 2 0 /ν) but do not follow N (0, 1). Fig. 3 shows the error ellipses corresponding to the limits σ 1 , σ 3 and σ 4 of eq. 10 (σ 2 not shown for clarity) and points (N R , N I ) representing the results from fitting bootstrapped data sets as explained in section 2.4. For all energies, N B ≥ 300. The numbers n 1 , n 2 , n 3 and n 4 indicate the number of these bootstrapped points inside each ellipse (as it has been shown in Fig. 2,  expressed as a percentage of N B ). Except for the lowest and highest energies more than 50% of the bootstrapped points are enclosed inside the σ 4 ellipse. In general they extend beyond the last ellipse but keeping the same correlation. If cases c) and h) of this figure a group of points diverge  from the main trend. The effect was traced to the influence of one single experimental point in both cases, been the divergent group part of a population in which that point is omitted. In Fig. 4, the parameters and their uncertainties given by the bootstrap method is compared with the ones given by Ref. [18]. It is seen that even though our assigned uncertainties are smaller the character of a regular TA is maintained.

6
Li + 80 Se system Since in this experiment χ 2 0 /ν ≈ 1 both null hypothesis almost coincide. The Kolmogorov-Smironov test indicates that residuals corresponding to all data taken together are normally distributed, and so are individual angular distributions for each energy, except for the data corresponding to 23 and 14.5 MeV. This reveals that for these two energies there is a subtle difference between experimental data and the theoretical model. Even though this difference is statistically significant, it would have remained unnoticed without this analysis. Normalization factors of around 1.6% applied to these angular distributions are enough to make them pass the test, but this may require further analysis. Figure 5 shows the error ellipses and bootstrapped points as in Fig. 2. Here, more than 43% of the bootstraped points are enclosed inside the σ 4 ellipse. The behavior of the parameters and their uncertainties given by the bootstrap limits is compared with the ones given by in Ref. [24] in Fig. 4.  Figure 4. Comparison between the dispersion relations reported in the original works and those obtained in the present analysis with the bootstrap technique a) for the 12 C + 208 Pb system and b) for the 6 Li + 80 Se system. The bootstrap results are slightly displaced in energy for clarity.  Figure 5. Distributions of N R and N I parameters that best adjust synthetic data sets generated by the bootstrap technique for the 6 Li + 80 Se system. The analysis based on the ellipses is described in section 3.3. The uncertainty bars indicate the medium values of the parameter distribution and their standard deviations.

Conclusions
We have shown that most studies use the standard rule from curve fitting: χ 2 ≤ χ 2 0 + 1 (eq. 3), to determine the uncertainty in fitting parameters (usually N R and N I ) from experimental data, regardless of their experimental value of χ 2 0 /ν. This is not satisfactory even from the point of view of conventional statistic: For the cases where χ 2 0 /ν > 1 that prescription is incorrect and should be replaced by χ 2 ≤ χ 2 0 + 2.3χ 2 0 /ν of eq. 10d. Parameters uncertainties which has been already derived using the χ 2 ≤ χ 2 0 + 1 recipe , can be easily scaled up (under the parabolic approximation of eq. 11) by a factor of 2.3χ 2 0 /ν. The bootstrap resampling method applied in this work produced a more realistic distribution of the parameters as shown in Figs. 2 b) and c). The resulting uncertainties are similar or somewhat larger than the ones obtained with the σ 4 ellipse (eq. 10d).
To judge the kind of TA it is very helpful to look at the series of 2-dim N R vs. N I plots with their respective uncertainty ellipses. In this way, the correlation between the parameters can be considered. In Fig. 4 our results, plotted in the conventional way, show a more reliable estimation of the parameter uncertainties than the quoted in the original works (smaller for 12 C + 208 Pb, similar or slightly larger for 6 Li + 80 Se). In these cases, however, the conclusions about the kind of TA was not modified. Further work is in progress to extend the present analysis to other systems, using bootstrap as well as other resampling techniques.