Assessing Uncertainties of Theoretical Atomic Transition Probabilities with Monte Carlo Random Trials

This paper suggests a method of evaluation of uncertainties in calculated transition probabilities by randomly varying parameters of an atomic code and comparing the results. A control code has been written to randomly vary the input parameters with a normal statistical distribution around initial values with a certain standard deviation. For this particular implementation, Cowan’s suite of atomic codes (R.D. Cowan, The Theory of Atomic Structure and Spectra, Berkeley, CA: University of California Press, 1981) was used to calculate radiative rates of magnetic-dipole and electric-quadrupole transitions within the ground configuration of titanium-like iron, Fe V. The Slater parameters used in the calculations were adjusted to fit experimental energy levels with Cowan’s least-squares fitting program, RCE. The standard deviations of the fitted parameters were used as input of the control code providing the distribution widths of random trials for these parameters. Propagation of errors through the matrix diagonalization and summation of basis state expansions leads to significant variations in the resulting transition rates. These variations vastly differ in their magnitude for different transitions, depending on their sensitivity to errors in parameters. With this method, the rate uncertainty can be individually assessed for each calculated transition.


Introduction
With the rapid improvement in computer power and quality of atomic structure codes, calculations of atomic structure and transition properties now become widely used in largescale simulations of physical conditions in complex plasma environments, such as found in fusion devices. Such simulations have wide range of applications, from plasma diagnostics to prediction of technological characteristics of industrial devices. To assess the accuracy of these simulations, it is important to have well-defined uncertainties for all calculated atomic parameters. Thus, critical evaluation of atomic data, implying estimation of uncertainties, has recently become one of the top priorities in fusion research [ 1 ].
Currently existing methods of evaluation of uncertainties of calculated transition probabilities were summarized by Wiese [ 2 ] and Kramida [ 3 ]. These methods heavily rely on comparisons between different calculations and between calculations and experiments.
In this paper, I suggest a new method of evaluation of uncertainties. It also relies on comparisons; however, these comparisons do not use any external data, but only the data produced by the same computational procedure. The base for comparisons is built from data generated with varied parameters of the atomic code. Application of this method is illustrated for the case of magnetic-dipole (M1) and electric-quadrupole (E2) transitions within the ground configuration of titanium-like iron (Fe V). The calculations were made with the suite of atomic structure codes by Cowan [ 4 ]. Uncertainties evaluated with this method are compared with recent critical compilation [ 5 ], where radiative rates of these transitions were also calculated with the same Cowan codes, and their uncertainties were evaluated by comparison with calculations of Nahar et al. [ 6 ].

Method Description
Calculation of transition probabilities involves several stages. These stages differ in different methods. In the non-relativistic or quasi-relativistic Hartree-Fock method, at first, the radial parts of wavefunctions are computed for each configuration in the single-configuration approximation. Then the Slater parameters and configuration-interaction (CI) parameters, as well as multipole transition integrals, are calculated from these radial functions. The parameters mentioned above are called hereafter the input parameters. The Hamiltonian matrix is built from these input parameters and diagonalized. The resulting eigenvalues and eigenvectors represent the initial calculated energy levels. Then these parameters are varied in the least-squares fitting (LSF) to find the set of parameters that best fits the experimental energy levels. Then these LSF parameters are used as input for the matrix diagonalization procedure, producing an improved set of eigenvalues and eigenvectors. These are used to compute the multipole transition probabilities. In most cases, because of the limited accuracy of the approximate atomic model, the LSF procedure does not exactly reproduce the experimental energy levels. Since the derivatives of the eigenvalues over the Slater and CI parameters can readily be calculated, the uncertainties (standard deviations, hereafter denoted as σ) of the LSF parameters can be computed from the differences between the fitted and experimental energies.
The idea of the method is to vary the LSF parameters in a random fashion (using a normal statistical distribution centered at the LSF values with a width equal to the σ of the LSF) and to evaluate the standard (root-mean-square) deviation of the line strengths calculated with these varied parameters.
distributions are uniform, meaning that all values of parameters within certain limits are equally probable. However, a normal statistical distribution seems to better reflect the common observation that the results of the LSF provide the best (most probable) values for the parameters insofar as the eigenvalues obtained with the fitted parameters exhibit symmetrical distributions around experimental values with the smallest standard deviations. The standard distribution has well-known statistical properties such as its symmetry around the central value (mean), its width at half-maximum (equal to 1.35 times the standard deviation around the mean) and the probabilities of occurrence of values deviating by more than one, two, or three standard deviations (32%, 5%, and 0.3%, respectively). Although other forms of distributions are possible, the normal distribution seems to be a good test case.
Cowan's codes [ 4 ] already have all necessary routines such as the LSF. Thus, this suite of codes was chosen as a test platform for the suggested method.
A few preliminary remarks should be made. The first one concerns the statistical distributions of the calculated quantities. In several recent papers, including [ 3 ], I made a statement that the logarithm of the calculated line strength, log(S), has much better statistical properties than the line strength itself, in the sense that the statistical distribution of log(S/ S*), where S* is a true value, is much closer to a normal distribution than S/S*. It turns out that this statement is generally incorrect. As discussed in Section 3, for the considered M1 and E2 transitions of Fe V, statistical distributions of both log(S/S*) and S/S* are asymmetrical. (In the context of this article, the "true value" S* is the one obtained with the LSF parameter values and ab initio values of the E2 transition moments). The distribution of log(S/S*) tends to be skewed to the negative side, while S/S* has a positive skew, meaning that there are more highly deviating values with S > S* than with S < S*. A somewhat better symmetry is observed for (S/S*) 1/3 . However, I have found that no single function of S can produce a normal statistical distribution of results for every transition. For any single function f(S), a few percent of all calculated transitions have extremely high volatility in the sense that small random variations of the input parameters with normal distributions around the initial values lead to a large number of strongly deviating results. For some transitions, there are too many trials in which the calculated quantity deviates from the initial value f(S*) by more than 3σ, as compared to the normal distribution. For such highly volatile transitions, which always exist if the same function f(S*) is used for all transitions, it is impossible to provide a definitive value of uncertainty in its normal sense, because there is a significant probability that the true value differs from the calculated one by ≥5σ. It is usually assumed that such highly volatile transitions are those that are affected by cancellations. This also turns out wrong. Some of the highly volatile transitions indeed have strong cancellation effects, but a significant fraction of them have weak cancellation. Thus, it is necessary to investigate the shape of the distribution function for each calculated transition and choose a function f(S) that has a statistical distribution closest to normal.
The second note concerns the identification of transitions. Generally, when the Slater parameters are changed, diagonalization of the Hamiltonian results in changes in both eigenvalues and eigenvectors, as well as in the predicted transition wavelengths. Identification of transitions produced by such different calculations poses a serious technical 6. The control code continues reading the part of the outg11 file from the main directory containing the transition data (initial and final energy levels, wavelengths, A-values, and cancellation factors).

7.
In each trial subdirectory, the control code randomly varies the Slater and CI parameters using the same procedure as for the E2 matrix elements, except that the widths of normal statistical distributions are set to be equal to the corresponding σ of the LSF. The parameters that were linked together at a fixed ratio in the LSF are varied in the same linked manner. Namely, their scaling factor is varied, but the ratios within each linked group remain fixed. For parameters that were fixed (not varied) in the LSF, the width of the normal distribution of the varied values was arbitrarily set to 2% of the parameter value. The varied parameters are substituted into the ing11 file (input file for RCG) prepared earlier in each subdirectory in step 3.

8.
In each trial subdirectory, the control code runs RCG, reads the resulting outg11 file containing new eigenvectors and sets of transition data, identifies the new eigenvectors with the old ones, kept in memory from step 5, and, for each transition, rescales the new A-values to the experimental (Ritz) wavelengths, and appends the statistics data.

9.
The accumulated statistics data are processed and results are printed to an output file.
The control code was written in Perl programming language. It has about 1,000 lines of code and uses several external Perl utility codes developed earlier and included with the Cowan code package [ 4 ]. These utilities include, for example, the eigenvector recognition routine and have about 3,000 lines of code. They were optimized for speed for the present purpose. With 10,000 random trials, the code execution takes about a few hours on a moderately powerful personal computer with 16 Gb of memory. Most of this time is taken by execution of RCG on each trial.

Results and Discussion
The 3d 4 configuration of Fe V consists of 34 energy levels spanning the range from zero to 94,000 cm −1 , all of which are experimentally known with uncertainties ranging from 0.3 cm −1 to 1.5 cm −1 [ 5 ]. My calculations with Cowan's codes yield 232 M1 transitions and 358 E2 transitions between these levels. For each of these transitions, the Monte Carlo method described in Section 2 produced a set of up to 10,000 A-values, providing a basis for statistical analysis. Results of this analysis are described below.

Input Parameters
As described in [ 5 ], the calculation of the even parity complex of Fe V included eight configurations, 3d 4 , 3d 3 (4s + 5s + 4d + 5d), and 3d 2 (4s 2 + 4s4d + 4d 2 ). Thus, there were 38 non-zero E2 matrix elements for transitions between these configurations, 86 Slater parameters (average energy E av , spin-orbit parameters ς 3d and ς 4d , direct and exchange Coulomb interaction parameters F 2,4 (nd,n′d) and G 0,2,4 (nl,n′l′), respectively, and effective parameters α 3d , β 3d , and T 3d ; those are coefficients of Casimir operators representing manybody effects in shells with equivalent electrons, see Cowan's book [ 4 ], section 16-7), and 61 CI parameters. In the LSF, many of these parameters were grouped with fixed ratios within each group. In particular, all CI parameters were linked in one group. In total, there were 15 such groups, which included 121 parameters, nine parameters were allowed to vary independently, and seven parameters were fixed in the LSF (note that these were allowed to vary with a relative σ of 2% in the Monte Carlo trials uncertainties in the E2 transition integrals, which depend on the accuracy of the radial wavefunctions. In the present investigation, these uncertainties were arbitrarily assumed to be equal to 1%, just so that transitions strongly sensitive to these uncertainties would be detected. The statistical distribution of the random trials was preliminary tested with a million trials and was found to be indeed very close to normal with the given variance (described above in Section 2).
The LSF parameters used as initial input for Monte Carlo trials are listed in Table 1. They were obtained in the work reported in [ 5 ]. The acronym 'HFR' appearing in the last two columns of Table 1 means 'Hartree-Fock-Relativistic' and is commonly used to denote the approximation used in ab initio calculations with Cowan's codes, the Hartree-Fock method with relativistic corrections. It should be noted that the standard deviation of the eigenvalues produced by the LSF from experimental energies is 117 cm −1 for all 220 experimentally known even levels and 41 cm −1 for the 34 levels of 3d 4 . Thus, it is a fairly good fitting with well-defined parameters. As can be seen in Table 1, the fractional uncertainties of P LSF are in the range between a small fraction of a percent for E av and about 10% for α 3d and β 3d , except for β 3d (3d 5 5s), which had a larger uncertainty.

Statistical Distributions of A-Values Obtained with 1000 Random Trials
Let us start with something familiar to atomic physicists, at least to some extent. Namely, let us consider how the standard deviations σ of A values behave depending on the line strength S. This dependence is presented, separately for M1 and E2 transitions, in the upper half of Figure 1. At first glance, the qualitative behavior of these plots confirms the general trend observed by many researchers and traditionally used to estimate the uncertainties of A values in different ranges of line strength (see, for example, [ 3 , 5 ]). Namely, for strong transitions the A values are well defined, while for weaker transitions uncertainties grow rapidly with decreasing S. However, these plots differ from traditional comparison plots [ 3 , 5 ] in that each data point represents not a single comparison between two calculations, but a standard deviation over 1,000 different calculations.
For strong M1 transitions with S > 1, the calculations are especially stable, giving σ over 1,000 trials well below 1%. I was tempted to use the word "accurate" in this statement.
However, one should remember that this discussion is about statistical distributions of results of one particular computational model. The accuracy of this model is necessarily limited by the approximations used. Thus, the uncertainties of its results must be greater than the standard deviations discussed here. This point will be further discussed in the following subsections.
If we look closer at the layout of the data points in the upper part of Figure 1, we can see that, in the same region of line strengths, the σ of resulting A values strongly differ from each other. For example, for M1 transitions with S = (0.001-0.01), most of the data points are clustered around σ ≈ 10%. However, there are many data points having σ greater or smaller than that by up to an order of magnitude. In the standard method of evaluation of uncertainties by comparing a few different calculations with each other, the apparent outliers with too large deviations would be assigned uncertainties much greater than the average for a given line strength. However, because of low statistics in such comparisons, there is a big chance of occasional coincidences, so some of such highly volatile transitions would not be detected. On the other hand, for some transitions that happen to be much more stable than the others with a similar S, uncertainties would be greatly overestimated. However, for M1 transitions it is not so. There are some transitions with a rather small degree of cancellation D c having δA/A greater or close to 100%. At the same time, many transitions with truly strong cancellation (D c > 0.5) have unexpectedly small deviations δA/A. Thus, even this statistically augmented approach does not allow one to definitively differentiate between "good" and "bad" results. Now let us take a look at statistical distributions of results for individual transitions, something that was never considered before in the literature. These distributions are shown in Figure 2 for three typical transitions. In this figure, the quantity on the vertical axis is the relative deviation δf/f (unlike Figure  One can see from Figure 2 that statistical distributions of results obtained for all presented transitions are asymmetrical and have significant numbers of far outliers, the number and location of which depends on the function used. (Hereafter, "outlier" means a data point deviating from the mean by more than one σ, and "far outlier" means a datum deviating from the mean by more than 3σ). For each transition, some functions appear to give "better" statistical distributions than others, both in terms of symmetry and the number of outliers. Thus, for transition (b), the function f = A 1/3 appears to give a fairly symmetrical distribution with a relatively small number of high deviations, while for transition (c), the best of the three functions is f = ln(A). However, for transition (a), although the "straight A" distribution appears to have the lowest number and size of high deviations, none of the three functions have a symmetrical distribution.
To investigate how the different functions f(A) fair in regards to the closeness of their statistical distributions to normal, I counted the numbers of trials producing relative deviations δf/f > nσ, where n is an integer number from 0 to +5 and similar numbers for negative δf/f less than −nσ. The quantity defining the fractional numbers of these counts is denoted N(n) below. These counts were co-added for all 590 investigated transitions and compared to the counts predicted with the cumulative normal distribution function, N norm (n). For reference purposes, the values of N norm (n) are given in Table 2.
One can see from Table 2 that, with 1,000 trials made for each transition, with the normal distribution, the probability of obtaining a result deviating from the mean by >5σ is negligible. However, if counts for all transitions are co-added, this gives a total number of trials of 590,000, and the probability of having at least one outlier deviating by more than ±5σ should be about 0.3 for the normal statistical distribution.
The ratios of N(n)/N norm (n) obtained with different functions f(A) are depicted in Figure 3. However, for all these functions the number of far outliers (those with |n| > 4) is much greater than one would expect for a normal distribution.
I attempted to divide all considered transitions in three groups based on the "best" function f(A) out of the three ones considered above. The "straight A" function turned out to be the best (of the three) for about 14% of all transitions (10% of M1 transitions and 17% of E2 transitions). The function f = A 1/3 is the best one for 47% of all transitions (51% of M1 and 45% of E2), and f = ln(A) is the best for 39% of all transitions (40% of M1 and 38% of E2). Thus, neither of these functions provides the best statistics for a significant majority of transitions.
Furthermore, even when the best of the three functions is used for each transition, the number of far outliers is still much greater than for the normal distribution. No apparent connection was found between the choice of the best function and physical characteristics of a transition, such as the type (M1 or E2), line strength, wavelength, energy levels involved, and the degree of cancellation.
Thus, the logarithmic function, which I earlier assumed (in [ 3 , 5 ] and several other papers) to provide statistics close to normal, in fact does so only for less than half of all transitions (at least it is so for the M1 and E2 transitions considered here).
In retrospective, the failure to find one or a few simple functions of A having good statistics for all transitions should not be surprising. The A values are numerically calculated in a procedure involving a matrix diagonalization, so it is obvious that they are complicated functions of the many input parameters. There is no reason to believe that there exists a single simple function of A providing a normal statistical distribution for all transitions.
However, it turned out possible to find the best function of a general type described by the Box-Cox transformation function [ 7 ]: (2) Although the Box-Cox transformation is defined piecewise, it is a smooth function of p, i.e., it does not have any singularity at p = 0. In applied statistics, there is a well-known technique for finding the optimal values of the parameter p of the Box-Cox transformation providing statistics closest to normal. This technique utilizes the so-called normal probability plots [ 8 ]. In such plots, the trial data are plotted against a theoretical normal distribution in such a way that the points should form an approximate straight line. Departures from this straight line indicate departures from normality. Examples of such normality plots are given in Figure 4 for the 5 D 1 -1 S1 0 M1 transition at 826.53 Å. Panel (a) was drawn for the "straight A" transformation function (p = 1 means no transformation), while panel (b) presents the transformed trial data with the optimal value of p = 0.133. It can be seen from panel (a) that the raw A data significantly deviate from the normal statistical distribution, while the near linearity of the plot in panel (b) shows that the data transformed with p = 0.133 are nearly normally distributed. The quantity on the vertical axis of Figure 4 is the normalized response value (f − f*)/σ, where σ is the standard deviation of f and the starred symbols mean the reference values obtained in the original LSF fitting. The trial data are ordered in the order of increasing response value. For each data point i (1 ≤ i ≤ n, where n is the number of trials, 1,000 in this case), the corresponding value of the normal order statistic median (given on the horizontal axis) is calculated with the function N i = G(U i ), where G is the inverse of the cumulative normal distribution function (probability that its argument is less than or equal to some value), and U i are the uniform order statistic medians defined as follows: (3) The optimal value of p is found by maximizing the correlation coefficient of the normal probability plot. This is illustrated in Figure 5a, where the correlation coefficient C(p) is plotted against p (this plot uses the trial data for the same transition as in Figure 4). The shape of the function C(p) turned out to look similar for all studied transitions, but the position of its maximum varies in a wide range of p.
An alternate, somewhat simpler method of finding the optimal Box-Cox transformation is to find the value of p at which the response data have zero skewness, which is defined as follows: (4) For a normal distribution, the skewness is zero. However, in general the zero value of skewness does not guarantee that the distribution is close to normal; it only ensures that it is symmetric. Another parameter characterizing the shape of a distribution function is kurtosis.
For this quantity, I use the definition that is usually referred to as "excess kurtosis": (5) For a normal distribution, thus defined kurtosis is zero. Positive values indicate a "peaked" distribution, while negative values indicate a "flat" distribution.
For transitions considered in this work, it turned out that, for the optimal values of p that minimize skewness, the corresponding values of kurtosis were in the range −0.7 to +3.3 with an average of −0.04 and a standard deviation of 0.23. For all considered transitions, thus found optimal values of p are close to those found by maximizing the correlation coefficient of the normal probability plot. Both methods of optimizing the parameter p are illustrated in Figure 5.
As noted above, the two methods of optimizing the Box-Cox transformation work almost equally well for the considered transitions. However, to avoid ambiguity, hereafter, when mentioning the optimal parameter p, I will always mean the results obtained with the first method, i.e., maximizing the correlation coefficient of the normal probability plot. The overall statistic obtained with parameters p optimized individually for each transition is shown in panel (d) of Figure 3. Considering the leftmost and rightmost data points in this plot, one should keep in mind that the trial counts for these points are very small. In particular, the n = 5 point represents the only case (out of 590,000) in which δf/f exceeded 5.
Thus, the apparent deviations from unity for |n| ≥ 4 seem to be explained by statistical noise. However, further analysis with larger statistics indicates that this is not true (see Section 3.3).
The optimal values of p found for the considered transitions of Fe V span a wide range from about −8 to +34 (a few transitions had extremely large values p > 100). The histogram showing the distribution of optimal values of p (excluding the abnormal ones with extremely large p) is shown in Figure 6. In this figure, the counts plotted on the vertical axis correspond to the number of transitions having the value of p between the one shown below each column and the one below an adjacent column. For example, for positive p, the value shown for p = 1 is the number of transitions with 1 ≤ p < 2; the value for p = 2 is the number of transitions with 2 ≤ p < 3, etc. Similarly, for negative p, the value shown for p = −1 is the number of transitions with −2 < p ≤ 1; the value for p = −2 is the number of transitions with −3 < p ≤ 2, etc. The point labeled p = 1/128 is an exception, because the corresponding interval encompasses the p = 0 value.
The distribution of optimal powers p has two distinct peaks, a greater one at p ≈ 0.3 and a smaller one near p ≈ −0.7. However, there is no peak near p = 0, indicating that there are very few transitions for which the logarithmic transformation produces statistics close to normal.

Results with Larger Statistics
After making several runs with 1,000 trials each, I noticed that for many transitions the optimal p values differ strongly from one run to another. This indicated that 1,000 trials give insufficient statistics for accurate determination of p. This is due to the fact that the optimal transformation strongly depends on the far outliers (i.e., those values of A that strongly deviate from the initial ones). As discussed above, 1,000 trials result in a relatively small number of results deviating from the initial values by more than 4σ. Therefore, I made an extended run with 10,000 trials. Combining the results of this run with those from four 1,000-trials runs allowed me to derive more accurate values of p and estimate their uncertainties as standard deviations of the mean over several runs. The distribution of these final values of p among all transitions is shown in Figure 7. One can see that it differs significantly from Figure 6 obtained from one 1,000-trials run. Now the histogram has only one major peak near p = 0.5.
Larger statistics revealed that for a few transitions the optimal Box-Cox transformation is far from normal. Examples of normal probability plots for some of such "abnormal" transitions are given in Figure 8. All such abnormal transitions are extremely weak. A few transitions required very large values of p for an optimal transformation. Thus, five transitions have p = (50-75), seven have p = (100-600), and one has p = 2,621. The latter is the 3 P2 1 -3 P2 2 M1 transition at λ = 66,872 Å with A = 4.541 × 10 −2 s −1 .
When optimal Box-Cox transformations are applied to all transitions, the overall comparison of statistical distributions with normal (similar to the plot given in Figure 3d) looks as shown in Figure 9.
One can see from Figure 9 that the Box-Cox transformation is unable to rectify statistical distributions for all transitions. In the particular case of Fe V transitions considered here, all far outliers accounted for in Figure 9 in the statistics with n = ±5 (more than 5σ) correspond to extremely weak transitions.

Uncertainties of A-Values
Strictly speaking, the standard deviations of optimally transformed A values obtained with the described method do not yield dependable estimates of uncertainties. Rather than that, they represent lower limits of those uncertainties. Additional contributions to the uncertainties come from the approximations made in the atomic model used. When calculations are made with one theoretical atomic code, they should include a study of convergence effects in regards to systematical improvement of the atomic model, such as a growing number of included layers of orbitals with increasing quantum numbers. Furthermore, a realistic estimate of uncertainties should involve both internal estimations based on the sensitivity of each transition to the variation of input atomic parameters and from comparisons with results obtained with other theoretical models.
None of these uncertainty contributions are considered in the present work. Nevertheless, the obtained optimal transformation parameters p and relative standard deviations σ f of transformed A values f(A, p) may be useful for further investigations. Therefore, they are given in Table 3 for each transition included previously in [ 5 ]. Most of these transitions have a branching fraction greater than 0.001. Transitions having significant contributions of both M1 and E2 type (greater than 0.01) are considered to be of a mixed type M1 + E2. For such mixed transitions, the shapes of the statistical distributions were investigated for the sum A M1 + A E2 , and the optimal values of p were found for these sums.
A rough estimate of relative σ for A values (σ A , standard deviation of δA/A) is given in Table 3 along with the standard deviations of transformed A-values, σ f . For the vast majority of transitions, the values of σ A and σ f are very close to each other. For only 9 transitions out of total 229 listed in Table 3, σ A differs from σ f by more than 10%. Most of these transitions are very weak and have branching fractions smaller than 4%. The largest differences between σ A and σ f (by a factor of more than 3) occur for two transitions. One of them is the 5 D 2 -3 P2 2 M1+E2 transition at 3,837.58 Å having a tiny branching fraction of 0.005% and an optimal p value of 0.096. The large difference of σ A from σ f for this transition is explained by the peculiarity of its statistical distribution. The normal probability plot for this transition looks very similar to the one shown in Figure 8a for another transition with a different optimal value of p = 2.23. Due to smallness of the branching fraction, peculiarity of this transition hardly matters in practice. However, another peculiar transition having σ A greater than σ f by a factor of 5 deserves a closer look. It is the 3 P2 1 -3 P2 2 M1 transition at 66,872 Å. Its peculiarity is caused by the extremely large value of p = 2,621, which implies that the statistical distribution of trial data for this transition is extremely asymmetrical. Ninety five percent of trial A-values for this transition are distributed symmetrically around the nominal value of 4.5406 × 10 −2 s −1 with a relative standard deviation of 0.013% (equal to σ f ). However, the remaining 5% of trials produce much larger deviations, up to 3%, all on the negative side (i.e., these trial A-values are smaller than the nominal value from the LSF).
Since both σ A and σ f are very small for this transition (0.06% and 0.013%, respectively), this statistical peculiarity may not be of practical importance. However, one cannot rule out the presence of such peculiarities in other atomic problems, where they may have significant consequences.
The closeness of σ A and σ f , if it is confirmed for other types of transitions and for other atomic systems, may lead to a conclusion that investigation of shapes of statistical distributions is not really necessary. However, one should keep in mind that the standard deviations alone do not fully describe statistical properties of the calculated quantities. If derived from a sufficiently large statistical base, they do allow one to say that the "true" value of the calculated quantity is within the limits of plus or minus one standard deviation from the mean value with a probability of 68%. However, if one needs to know the limits of possible deviations more strictly, say, with a probability of 95%, the values of σ are of no use without additional knowledge of the shape of the distribution. Sufficient data about these shapes is provided in Table 3 by the values of p describing the optimal Box-Cox transformation. For the particular atomic problem considered here, most of them are clustered around 0.3 (see Figure 7). This means that statistical distributions of most considered transitions are skewed to the positive side, i.e., there are greater number of trials with A > A* than with A < A*. Values of p > 1 result in negative skew, meaning that it is more probable to obtain A < A*. This information cannot be obtained from standard deviations alone. For the particular 3 P2 1 -3 P2 2 M1 transition at 66,872 Å discussed above, knowledge of σ A = 0.06% allows one to say that the calculated A-value is (4.5406 ± 0.0027) × 10 −2 s −1 at the 68% confidence level. However, if one needs error limits at the 95% confidence level, the result can be obtained only if p and σ f are known.
As noted above, the standard deviations for straight A values are not of much use for statistical considerations. Nevertheless, they can be compared with similar estimates obtained earlier in [ 5 ] from comparisons with the calculations of Nahar et al. [ 6 ]. The latter authors used the Breit-Pauli R-matrix method to compute the E1, M1, and E2 transition probabilities of Fe V. The quality of their calculations was largely determined by the wavefunctions of the Fe VI ionic core calculated with the Breit-Pauli version of the SUPERSTRUCTURE code [ 9 ]. The accuracy of such calculations is usually similar to that of ab initio calculations with Cowan's codes. The SUPERSTRUCTURE code, as well as Cowan's codes, uses a non-relativistic approximation with relativistic corrections introduced as perturbations. The principal difference is that, to compute one-electron orbitals, SUPERSTRUCTURE uses a scaled statistical model of the transitions the uncertainties given in [ 5 ] were underestimated, and the corresponding category of accuracy should be degraded by one category (e.g., "C+" should be changed to "C", and "B" to "C+"). These few transitions are marked with an asterisk in the σ A column.
For the rest of the transitions, the uncertainties estimated in [ 5 ] are greater (in most cases much greater) than those found here by random trials. This confirms the general validity of the traditional method of critical assessment of uncertainties, but indicates that a small statistical basis of comparisons sometimes leads to underestimated uncertainties.
As noted in Section 3.3, for a few transitions the optimal Box-Cox transformation produces statistical distributions that are rather far from normal, meaning that there are too many trials in which transformed A-values deviate from the initial (LSF) ones by more than 5σ. Such abnormal transitions are marked in Table 3 with an asterisk in the last column.

Required Number of Comparisons
The number of comparisons required for a reliable estimation of uncertainty of calculated Avalues can be estimated by comparing the values of σ A obtained with a progressively increasing number of trials. Such comparison shows that, with just ten random trials, the resulting σ A differs from the "true" value (obtained with 10,000 trials) by more than 20% for

Further Considerations
Statistical distributions of A values considered here were obtained using several arbitrary assumptions described in the Introduction. First of all, the variations of the input atomic parameters such as Slater and CI parameters were assumed to have normal statistical distributions. There is no physical or mathematical reason behind this assumption. A further investigation is needed in order to find actual shapes of these distributions and their influence on the final A values. The same is true for the E2 transition integrals. In addition, the widths of statistical distributions of the latter were arbitrarily assumed to be equal to a fixed value of 1%. Estimation of the actual variances of the E2 transition matrix elements requires a separate study. Furthermore, the various atomic parameters are strongly correlated. These correlations were accounted in the present work only partially, via linking some of the parameters in groups having fixed parameter ratios. The actual correlations are much more complex and should be investigated.
The method suggested here was relatively easy to implement with Cowan's suite of codes, which include the LSF procedure with well-defined uncertainties of the fitted parameters.
Other codes also have adjustable parameters. For example, SUPERSTRUCTURE [ 9 ] uses a set of adjustable scaling parameters λ of the Thomas-Fermi-Dirac potential. Some theorists adjust these parameters to produce energy levels most close to observed ones. This procedure is, in principle, similar to Cowan's LSF. To make it a true LSF, one needs a matrix of derivatives [∂λ/∂E] similar to [∂P/∂E] in Equation (1). Then the uncertainties of adjusted λ could be estimated with a similar formula. I am not aware if such procedure is implemented in SUPERSTRUCTURE, and I never saw uncertainties of λ reported in publications. However, even if the adjustment of λ is manual, one can obtain some idea of their uncertainties by trial and error, and the resulting uncertainties in A-values can be estimated by a method similar to the one described here.
Adjustable parameters are also used by Hibbert in the CIV3 code [ 10 ]. Namely, in the process that he calls "fine tuning", he adjusts the diagonal Hamiltonian matrix elements to achieve the best agreement of calculated energies with experiment [ 11 ]. This adjustment is not much different from adjusting E av in Cowan's LSF. Thus, it must be possible to investigate the uncertainty of this fine tuning and its effect on A-values.
Other sophisticated ab initio atomic codes implementing multi-configuration Hartree-Fock and Dirac-Fock methods [ 12 , 13 ] also compute the Slater and CI parameters. However, in the current implementation of these codes these parameters cannot be adjusted, which is quite unfortunate. To implement the described method with these codes, they must be modified. Estimation of uncertainties of ab initio calculations without the use of experimental data would require different methods. However, these methods must necessarily involve an analysis of error propagation through matrix equations.
In principle, a clever statistician should be able to predict the expected transformations of distribution functions of calculated A-values on the basis of known matrix-diagonalization procedures and formulas for summation of various contributions to the line strength used in atomic codes, avoiding the computationally expensive empirical method used here. So far, I could not find any methods for that. In-depth investigation of statistical properties of atomic calculations requires much more than basic training in applied statistics. The present work shows only a first glimpse of these properties and gives more questions than answers. Some aspects of the problems touched upon appear to be on the cutting edge of statistical theory.
Solving these problems calls for development of a relatively new field of atomic physics, statistical atomic physics, the main subjects of which should be investigation of statistical

Conclusions
The present work suggests a method for evaluation of sensitivity of calculated transition probabilities (A values) to small variations of input atomic parameters by Monte Carlo random trials. This method provides an insight into shapes of statistical distributions of the A values. These distributions were found to be far from normal for most M1 and E2 transitions within the ground configuration of Fe V, which were used as a test problem. However, for almost all transitions, it turned out possible to find an optimal Box-Cox transformation that results in near normal statistical distribution. If the range of possible random variations of the input parameters is known, e.g., from a least-squares fitting of Slater parameters, this method leads to statistically sound estimates of internal uncertainties of the particular computational atomic model used.
In order to determine the total uncertainties of the calculation, these internal uncertainty estimates should be combined with studies of effects of approximations made in the atomic model, for example, by comparing A values obtained with increasing number of included basis states, estimating convergence trends, and extrapolating these trends to an infinite basis size. Additionally, the calculated A values should be compared with other independent ones based on different atomic models of comparable quality (for example, non-relativistic calculations with perturbative account for relativistic effects can be compared with fully relativistic calculations). Neither of these additional uncertainty estimates was made in the present work, so the numerical results given here represent the lower bounds of uncertainties. Even so, they helped detect and correct uncertainties underestimated earlier by other methods.
One conclusion following from the present results is that uncertainties should be estimated not for the straight A values, but for the transformed ones that have a statistical distribution close to normal. This requires that, for each transition, the optimized transformation parameter p should be given along with the A value and an estimate of σ for the transformed A value.
When results of two different theoretical models are compared to each other, one can expect that the shapes of statistical distributions of the A value for the same transition should be similar for different calculations, at least if cancellation is weak. This statement is not substantiated by any factual data, but is based on a general consideration of the method of computation of line strengths in all currently available atomic codes. Expanding the basis set should only add some minor terms in the summation series of these calculations. Thus, a significant change in the distribution shape should be expected only for transitions with strong cancellation.
There is an important implication for Monte Carlo simulations of plasma kinetics, namely, the distributions of the input atomic parameters (such as the A values and collision rates) should not be assumed normal. Instead, they must be skewed according to a Box-Cox transformation optimized individually for each transition.         Similar to Figure 3d, with larger statistics (10,000 trials) and with optimal Box-Cox transformations applied to all transitions.    Fractional number of outliers deviating from the mean by more than n standard deviations for the normal statistical distribution, N norm (n).         (7) 0.018 a The energy levels and Ritz wavelengths calculated from them are taken from [ 5 ]. Wavelengths between 2,000 Å and 20,000 Å are in standard air; shorter and longer wavelengths are in vacuum; Atoms. Author manuscript; available in PMC 2016 June 01. e Transition type is specified as mixed M1 + E2 for transitions having the fraction of the minor contribution to the total A value greater than 1%; f Fraction of E2 transition in the total A value was calculated for each trial. The given value is the result of the initial LSF calculation. The quantity in parentheses is the standard deviation of the E2 fraction over 10,000 trial calculations (in the units of the last digit of the value); g The parameter p of the optimal Box-Cox transformation is determined as a weighted mean over five runs, four with 1,000 trials each and one with 10,000 trials. The quantity in parentheses is the weighted standard deviation of the mean over five runs; i Relative standard deviation of transformed A values using the optimal Box-Cox transformation with the given parameter p, over 10,000 trial calculations (percent). The starred values denote transitions for which the optimal Box-Cox transformation yields statistical distributions that are far from normal.
Atoms. Author manuscript; available in PMC 2016 June 01.