Improved Recursive Computation of Clebsch-Gordan Coefficients

Fast, accurate, and stable computation of the Clebsch-Gordan (C-G) coefficients is always desirable, for example, in light scattering simulations, the translation of the multipole fields, quantum physics and chemistry. Current recursive methods for computing the C-G coefficients are often unstable for large quantum numbers due to numerical overflow or underflow. In this paper, we present an improved method, the so-called sign-exponent recurrence, for the recursive computation of C-G coefficients. The result shows that the proposed method can significantly improve the stability of the computation without losing its efficiency, producing accurate values for the C-G coefficients even with very large quantum numbers.


Introduction
The Clebsch-Gordan (C-G) coefficients arise whenever the coupling of two angular momenta is involved. The computation of C-G coefficients has a broad range of applications including particle light scattering simulations [1,2,3,4], fast multipole methods [5,6], spherical polar Fourier transform [7], quantum physics and chemistry [8]. In the context of light scattering by small particles, for instance, the C-G coefficients are needed for the realization of analytical random orientation average using T-matrix. Developed by M. Mishchenko [1], the analytical random orientation scheme is perhaps one of the greatest advantages of T-matrix method, which could save massive amount of computational time if the orientation averaged properties are needed. In addition, the C-G coefficient also arises when one need to compute the translation of multipole fields [4]. Consequently, the importance of obtaining fast, reliable and accurate computation of C-G coefficients should not be underestimated. Currently, the computation of C-G coefficients with large quantum numbers is often based on a modified recursion method, which was originally proposed by Schulten and Gordon [9,10] and later modified and implemented by M. Mishchenko [1] and Wielaard et al. [2] in T-matrix simulations. According to [11], the modified recursive method can compute the C-G coefficients with quantum number up to 150 in a stable and accurate manner. Although the size of T matrix up to this order is already quite large, it is always desirable to develop a highly stable and accurate method for calculating the C-G coefficients with quantum number as large as possible. In this paper, we present an improved method which could significantly extend this limitation without losing the accuracy and efficiency.
The paper is organised as follows. In section 2, we briefly introduce the widely applied recursion method for the computation of C-G coefficients and analyse its limitations. Section 3 presents our method. In section 4, we demonstrate the stability, accuracy, and efficiency of the proposed method by comparing the results with those from the previous methods. Section 5 concludes this study.

The recursive computation of C-G coefficients
The C-G coefficients are related to the Wigner's 3j symbol in accordance with where the j 1,2,3 and m 1,2,3 are the so-called principle and magnetic angularmomentum quantum numbers respectively. The C-G coefficient will be zero unless the following conditions are satisfied simultaneously, To be consistent with Schulten and Gordon [9] on notation, we focus on the discussion of Wigner's 3j symbol, while the C-G coefficients can be easily derived from Eq. 1. For small quantum numbers, the 3j symbols can be computed conveniently using Racah's formula [12]. For large quantum numbers, however, the computation using Racah's formula becomes very expansive and suffers numerical instabilities due to the high order factorials. It is believed that the 3j symbols with large quantum numbers can be evaluated most efficiently with the recurrence relations. In light scattering simulations, one needs to compute the 3j symbols with a range of principle quantum number j or magnetic quantum number m. Because the symmetry properties of the 3j symbols (see Appendix D of [11]), the recurrence relations for one of the three quantum numbers shall be enough to compute the recursion of all others. Without loss of generality, let us focus on the computation of j 3 and m 2 . The recurrence relations for the two quantum numbers are [9]: where For recursive computation using Eq.(4), j 3 lies in the following range For recursive computation using Eq.(5), m 2 lies in the following range − min(j 2 , j 3 + m 1 ) ≤ m 2 ≤ min(j 2 , j 3 − m 1 ).
In the original method proposed by Schulten and Gordon [9], the recursions can start from arbitrary real number, e.g., unity, and go both forward and backward from the minimum and maximum quantum number respectively. The method therefore requires the computation of a scaling factor such that the forward and backward recursions give the same number at the intermediate quantum number. The coefficients can then be determined by applying the unitary properties: In addition to iteration direction [9], there are two more sources that could cause inaccuracy or numerical instability via recursion. Firstly, because the recursion starts with arbitrary number, the computed values need to be scaled twice, first by the scaling factor and then by the normalization factor. The errors of the factors, possibly arising from the inaccuracy of particular coefficients or their ratios, could propagate to the whole group of the computed values. To illustrate this, we compare the exact values with those computed by scaling and normalization as displayed in Fig. 1. It can be seen that most of the computed values can be shifted by certain magnitude due to the multiplication of the factors. The second source of errors in the existing methods is the lack of a mechanism to avoid numerical outflows or underflows. Note that in the case of high order quantum numbers, the magnitude of quantities (or their ratios) involved in the computations could be extremely small or large.
In [13], Luscombe and Luban propose to iterate the ratio of two successive 3j symbols to avoid numerical overflows. Nevertheless, their method is not without drawbacks. Firstly, it has to perform normailzation, which could cause unnecessary shifting of values. Secondly, in their method, the values of C-G coefficients are to be obtained by multiplication of many ratios, and this could still induce numerical overflow/underflow, even though the iteration of the ratios may have no such risks. Thirdly, the necessity of identifying classical and nonclassical regions complicates the algorithm.
One way to remove the necessity of scaling and normalization is to start with an exact value of C-G coefficient at the minimum or maximum quantum number. As described in Appendix A, there are four different cases to be considered for j 3 = j 3 min , while there is only one case for j 3 = j 3max . The usage of exact starting values exclude the necessity of scaling and normalization, which improves the computational stability and accuracy. However, it does not exclude the possibility of numerical underflow/overflow, because for large quantum numbers, the factorial computations for the starting values will likely exceed the precision of the arithmetic.

The sign-exponent recurrence method
If the quantum numbers satisfy the conditions of Eq. 2 and Eq. 3, the 3j symbols are generally non-zero. However, it is well-known that some of the coefficients can be "accidentally" zero even if the conditions are fulfilled [14]. Such zeros are called non-trivial zeros. But the non-trivial zeros are quite rarely encountered, for the moment, let us assume that all the coefficients involved are non-zero. In this case, one can write arbitrary 3j symbols as where k(j 3 ) and h(m 2 ) are real functions of j 3 and m 3 respectively, and The sign function is defined as By introducing Eq. 14 into Eq. 4 , we obtain From Eq.18, one can obtain the recurrence relations for both s(j 3 ) and k(j 3 ), i.e., Similarly, recurrence relation for m 2 can be obtained by introducing Eq. 15 into Eq. 5, i.e., At the starting minimum quantum number, we have Therefore the iteration becomes It can be seen that we have separated the original recurrence relation into sign-recursion and exponent-recursion. Note that the above relations are valid for the forward direction which increases the quantum number. Similarly, it would be straightforward to derive the backward recurrence relation that reduces the quantum number. The basic idea behind such recurrence relations is that we could focus on computing the sign and exponent of the coefficient, rather than the coefficient itself. In principle, as long as the sign and exponent are computed accurately, the coefficient can always be calculated accurately via Eq. 14 or Eq. 15. Because the 3j symbols can vary many orders of magnitude, the sign-exponent recurrence can significantly reduce the risk of numerical underflows/overflows. To apply the derived relations, we just need to compute the starting exponent and sign of the coefficient at the minimum or maximum quantum numbers. According to our numerical tests, the method would not induce numerical underflow/overflow even with quantum number larger than 10 million. For now, we shall consider the problem of encountering non-trivial zeros. Once zero-value is encountered, the sign value will become 0 and the exponent will become −∞, making the computation of ∆(m 2 ) or ∆(j 3 ) meaningless. Therefore, the above recurrence relations must be avoided and the original three-term linear recurrence relations should be applied. The condition for using the sign-exponent recursions shall be that both s(j 3 − 1) and s(j 3 ) are non-zero (same for m 2 case), otherwise the three-term linear relations should be invoked. For large quantum numbers, the non-trivial zeros are rarely encountered. For most of the cases, only the sign-exponent iteration is invoked in the computation.

Results and Discussion
In this section, we shall discuss the stability, accuracy, and efficiency of the proposed method by comparing with the widely applied three-term linear recurrence with exact starting values. For large quantum numbers, in general, all recursion-based methods are much faster than those using the direct definition or formula. The efficiency of the proposed method is almost the same as the previous recursion methods. To quantify the accuracy of the computation, we define the following error term, which is consistent with [4].
Please note that all the results are obtained by applying the double-precision arithmetic. We compute the starting values for m 2 recurrence by firstly using the recurrence relation for j 3 , therefore the comparison on m 2 recurrence is preferred. The following figure displays a comparison between the three-term linear recurrence and sign-exponent recurrence method. It can be seen that the results of the two are indistinguishable from the figure. In fact, for the three-term linear recursion, R = 1.0334 × 10 −12 , while for the sign-exponent recursion, R = 1.0214 × 10 −12 . The values of R suggest that the sign-exponent recurrence method has at least the same performance as the three-term linear recurrence on numerical accuracy. The biggest advantage of the proposed method perhaps is that it can avoid numerical underflow when computing the starting values and this is crucial for dealing with large quantum numbers because the starting values could be extremely small. Figure 3 demonstrates the comparison for the case of very large quantum numbers. The three-term linear recurrence suffers numerical underflow and becomes zero throughout the iteration. On the contrary, the sign-exponent recurrence method excludes the risk of obtaining zero starting values. Consequently, for the three-term linear recurrence, R = 1, while for the sign-exponent recursion, R = 5.9769 × 10 −10 , which still maintains high accuracy. Based on our extensive tests, the proposed method is generally much more stable than the original linear recursion method, while having the same level of numerical accuracy and efficiency. To further demonstrate the accuracy and stability of our method, in Appendix B, we provide a link to our Matlab code and more test examples in comparisons with the most accurate package Python SymPy.

Conclusion
In this paper, an improved recursion method for computing the C-G coefficients is introduced. Specifically, the method separates the recursion process into sign-recursion and exponent-recursion, while the C-G value itself is not involved in the recursion except that a non-trivial zero occurs. The C-G values can be obtained after the computation of their signs and exponents.
By using the sign-exponent recursion, the method removes the risk of generating numerical overflows or underflows. The results presented in this paper, together with our extensive tests, have shown that the sign-exponent recurrence method is in general more stable than the original three-term linear recurrence method, while having the same level of accuracy and efficiency.

Appendix A. The computation of starting values
The recursive computation without using the scaling and normalization relies on the exact computation of the starting values. For backward recursion, the starting value is unique, i.e., For the forward iteration, we have four possibilities, depending on the values of j 3 min , i.e., In the computation of the above values, we can compute the logarithm of the factorials to avoid overflow. To be more specific, one can write arbitrary starting value as Obviously, the starting values for s and k functions are k(j start ) = 1 2 [ln(a 1 !) + ln(a 2 !) + ln(a 3 !) + ln(a 4 !) To compute the logarithm of the factorials, we may use the following formula to avoid numerical overflow,

Appendix B. Comparison with Python library SymPy
To demonstrate the accuracy of the proposed method, we compile some of the results from Python library SymPy and compare with the proposed method in this study. The code for our method can be accessed from [15]. The code for SymPy can be obtained from [16]. The SymPy is based on symbolic manipulation, which could be considered as the most accurate method. The computational time by the recursion methods ranges between 0.01 to 0.03 seconds for most of the cases with quantum number smaller than 1000, while the SymPy package could take up to roughly a few seconds. It can be seen that the two recursion methods are both very accurate comparing to SymPy and their levels of accuracy are pretty much the same. However, as we increase the quantum numbers, the three-term linear recursion becomes unstable and generates zeros due to numerical underflows. The Python SymPy simply produces errors at very large quantum numbers, while the sign-exponent method is still stable and produces reasonable results.