Analytical formula for numerical evaluations of the Wigner rotation matrices at high spins

0.4 MB·English

by  Naoki Tajima

Analytical formula for numerical evaluations of the Wigner rotation matrices at high
spins
∗
Naoki Tajima
Department of Applied Physics, University of Fukui, Bunkyo 3-9-1, Fukui, 910-8507, Japan
(Dated: January 27, 2015)
The Wigner d function, which is the essential part of an irreducible representation of SU(2) and
SO(3) parameterized with Euler angles, has been know to suffer from a serious numerical errors at
high spins, if it is calculated by means of the Wigner formula as a polynomial of cos and sin of
half of the second Euler angle. This paper shows a way to avoid this problem by expressing the d
functions as theFourier series of the half angle. A precise numerical table of the coefficients of the
series is provided as SupplementalMaterial.
5
1
0 I. INTRODUCTION where nmin and nmax are zero or positive integers,
2
n The matrix elements of the rotation operator between nmin =max(0,k m), (4)
−
a angular-momentum eigenstates are called the Wigner D n =min(j m,j+k), (5)
max
J function[1]. Whentherotationisspecifiedwiththethree −

6 Euler angles (φ, θ, ψ), and the eigenstates are labeled and
2 with the magnitude (j = 0,1,1,3, ) and the z com-
h] pmoennetnutm(mvecotrork,=the−Dj,−fujn+cti1o2,n··c·a,2nj)·b·oe·fdtehceomanpgousleadr imntoo- Wnjmk(θ)=wnjmk · cos2θ 2j+k−m−2n −sin2θ m−k+2n(6,)
t three factors, (cid:0) (cid:1) (cid:0) (cid:1)
-
cl Dj (φ,θ,ψ)= jme−iφˆjze−iθˆjye−iψˆjz jk (j+m)!(j m)!(j+k)!(j k)!
u mk h | | i wjmk = − − . (7)
n =e−i(mφ+kψ)djmk(θ), (1) n (jp−m−n)!(j+k−n)!(n+m−k)!n!
[
where However, this formula suffers from a serious loss of pre-

1
v djmk(θ)=hjm|e−iθˆjy|jki (2) cbiosriohnooadthofigθh=sp0in,πs.(i.e., for largej) exceptin the neigh-
7
4 is the nontrivial part which needs some method to eval- Forexample,assumingthatj isapositiveinteger,θ =
3 uate. It is called the Wigner (small) d function. In the π, and m=k =0, one obtains
6 standardphaseconventionforangular-momentumeigen- 2
0 states, the matrix elements of ˆjy are purely imaginary 2
π 1 j!
1. and thus djmk(θ) takes on real numbers. Wnj00 2 = 2j (j n)!n! , (8)
0 One may be anxious about the fact that Bohr and (cid:16) (cid:17) (cid:20) − (cid:21)
5 Mottelson[2] define the rotation matrix Dj as the
1 complex conjugate of the right-hand side (mr.kh.s.) of which, if j is even, becomes maximum at n= j,
2
v: Eq.(1). However,theirdefinitionfordj isidenticalwith
mk
i Wigner’s and the definition of the d function is unique. π 1 j! 2
X Wj00 = 2j. (9)
The Wigner D function is used in various fields of j/2 2 2j [(j/2)!]2 ≈ πj
r r
a physics. In some applications, those for large values of j (cid:16) (cid:17)
(say, > 50) are necessary. An example in nuclear struc-
(TheStirling’sformulaisusedinthelastapproximation.)
ture physics is the projection from spatially deformed
The absolute value of the d function is not greater than
solutions of modern realistic mean-field models to eigen-
one because it is a matrix element of a unitary operator
states of large angular momentum. For toy models, too,
between normalized states. Wigner’s formula expresses
one occasionally needs d functions for very large j to
the d function as a result of cancellation among terms
confirm the validity of one’s picture under extreme con-
of possibly huge size, Wjmk 2j. For j 54, the pre-
ditions. (e.g., a two-rotormodel of Refs. [3] and [4]). n ∼ ∼
cision of double-precision floating-point numbers (53-bit
The explicit form of the d function is given by the
mantissa) is lost completely, and even quadruple preci-
Wigner formula[1], which can be written as,
sion float numbers (113-bit mantissa) is lost completely
nmax for j 114.
djmk(θ)= (−1)nWnjmk(θ), (3) A f∼ew remedies have been proposed [5, 6] but the re-
n=Xnmin sultsarenotcompletelyfree ofthe precisionloss. Inthis
paper,Iinvestigatethedetailsofthis lossofsignificance,
and then present a perfect remedy to avoid such numer-
∗ [email protected] ical difficulty.
2
II. THE FOURIER-SERIES EXPRESSION OF d half-integervalues),one canexpressthe coefficients tjmk
ν
FUNCTIONS by an integral
1 4π
Onecanseeeasilythattermslikecosλ θ sinµ θ appear- tjmk = dj (θ)f(νθ)dθ, (15)
2 2 ν 2π(1+δ ) mk
inginther.h.s.ofEq.(6),whereλandµarezeroorpos- ν0 Z0
itive integers such that λ+µ 2j, can be expressed as
linear combinations of terms si≤nκθ (when µ is odd) and where δν0 =1 for ν =0 and δν0 =0 for ν = 21,1,32,···.
2 BysubstitutingthedfunctioninEq.(15)withEqs.(3-
cosκ2θ (when µ is even) with integers κ in 0≤κ≤λ+µ 7), I have derived a more useful expression for tjνmk con-
and κ λ+µ (mod 2), by means of repeated applica-
taining only four elementary operations of arithmetic,
≡
tionsoftheelementarytrigonometricidentitiescalledthe
product-to-sumidentitiesorthe prosthaphaeresisformu- 2( 1)m−k nmax [ν−12p]
las. tjmk = − ( 1)nwjmk ( 1)r 2ν
Becausethepowerofsinθ2 isµ=m−k+2ninEq.(6), ν 1+δν0 n=Xnmin − n Xr=0 − (cid:18)2r+p(cid:19)
onecanseethatµ m k (mod 2)andthatthe dfunc- 1
tionisaneven(odd≡)fun−ctionifm kiseven(odd),which × 2πI2(j+ν−n−r)−m+k−p,2(n+r)+m−k+p (16)
−
can be expanded only with cos (sin) function. This may
where n and n are those already defined by
also be deduced from Eq. (A3), one of the properties of min max
Eqs. (4) and (5), the square brackets are the floor func-
the d function which I have enumerated in the appendix
tion, i.e., [l+x]=l for integer l and real x in [0,1),
A.
From these considerations, one can conclude that the p m k (mod2), (17)
≡| − |
Fourier expansion of the d function has the following
form: i.e.,
p=0 for k =m,m 2,m 4, , (18)
dj (θ)= tjmkf(νθ), (10) ± ± ···
mk ν p=1 for k =m 1,m 3, , (19)
Xν ± ± ···
and
where
2π
cos even Iλµ = cosλxsinµxdx (20)
f = for m k, (11)
sin odd − Z0
(cid:26) (cid:27) (cid:26) (cid:27)
with zero or positive integers for λ and µ. If both λ and
and the summation runs over µ are even,
ν =ν ,ν +1, ,j, (12) 2π(λ 1)!!(µ 1)!!
min min ··· Iλµ = − − , (21)
(λ+µ)!!
with the values of ν given in Table I.
min
while I =0 otherwise.
λµ
Unlike the r.h.s. of Eq. (3), where the terms can have
TABLE I. The valuesof νmin appearing in Eq. (10) hugesizesandthusthenumericalerrorisaseriousprob-
lem, the r.h.s. of Eq. (10) is a summation of terms of
Even m−k Odd m−k order one or less and hence the problem is expected to
disappear. This can be seen by calculating the integrals
Even 2j 0 1
of the squares of the both sides of Eq. (10),
Odd 2j 1 1
2 2
4π 4π
dj (θ)2dθ = tjmktjmk f(νθ)f(µθ)dθ
mk ν µ
For example, for j = 7 and m = k = 1, the Wigner Z0 Xν Xµ Z0
formula (3) gives an exp2ression, − 2 = 4π tjmk 2. (22)
2 δ ν
ν0
d27 (θ)=sin7 θ 12cos2 θ sin5 θ +18cos4 θ sin3 θ Xν − (cid:0) (cid:1)
12,−12 2 − 2 2 2 2 Because the absolute values of d functions are 1, the
−4cos6 2θ sinθ2, (13) l2eft-δhan.dsideis≤4πand,consequently,itholds≤|tjνmk|≤
ν0
−
which can be rewritten in the form (10) as A further study from the numerical point of view has
indicatedthatthemaximum(amongallthepossiblecom-
7 35 sin7θ 5 sin5θ +15 sin3θ 9 sinθ binations of m, k, and ν) value of tjmk is 1 for j 1,
d2 (θ)= 2 − 2 2 − 2. | ν | ≤
12,−12 − 64 decreases as j increases from an integer to the next half
(14) integer, and does not change as j increases from a half
Byutilizing theorthogonalityofcosνθ andsinνθ over integertothenextinteger. Fortheinterval50 j 100,
≤ ≤
0 θ 4π (consideringthatν cantakebothintegerand the maximum value for integer j behaves as 1.13/√j.
≤ ≤ ≈
3
III. COMPUTATION OF THE NUMERICAL tj,−k,−m =tjmk. (26)
ν ν
VALUES OF THE COEFFICIENTS
Hence, the numerical data for tjmk in the Supplemental
ν
Material are given only for m 0 and k m.
Unfortunately, Eq. (16) also suffers froma seriousloss There are 1(j + 3≥)4 coeffici≤en|ts |for j =
of significant digits in ordinary floating-point numerical ≈ 2 max 2
0,1,1, ,j . Thesizeofthememorytostorethemas
calculations. Indeed, for even integer j, the term having 2 ··· max
64-bit floating-point numbers amounts to 27 (404) MiB
the maximum magnitude among those in the r.h.s. of
forj =50(100). Mydataaregivenastextfilesforthe
Eq. (16) occurs at m = k = 0, n = r = 1j, ν = j, with max
2 sake of compatibility, whose sizes are not very different
the maximum value of
from the above memory sizes after they are compressed
1 22(j+1) (with the software GZIP).
πwjj/020 2jj I2j,2j ≈ (πj)2 , (23) Themagnitudeofthe coefficienttjνmk becomessmaller
(cid:18) (cid:19) for larger values of j. However, even at j = 100, 80%
which is roughly 2j times as large as the value given by (90%) of the coefficients are larger than 10−5 (10−10).
Eq. (9). Some of the coefficients vanish according to rules un-
To avoid this problem, I first evaluated the r.h.s. of mentionedsofar. Forexample,tjmk =0ifjisaninteger,
ν
Eq. (16) rigorously as rational numbers or square root m = 0 and/or k = 0, and j ν 1 (mod 2). One can
− ≡
of rationalnumbers by means of a formula-manipulation prove this vanishment by using Eq. (A8). I do not use
software MAXIMA. Numerical values to be used in pro- these additional rules but simply give zero values in the
grams coded in Fortran, C, etc., can be calculated from data files.
such rigorous numbers to the full precision of the 64-bit The evaluations of cos νθ and sin νθ should be calcu-
floating-point number. lated by means of a recursion relation,
However,thecomputationtimeturnedouttobeexces-
sivelylongforlargevaluesofj. Tospeedupthecomputa- cos(ν+1)θ cos θ, sin θ cos νθ
= − ,
tion,IhavechangedthemethodtoevaluateEq.(16)not sin(ν+1)θ ! sin θ, cos θ ! sin νθ !
rigorously but in terms of high-precision floating-point (27)
numbers (of MAXIMA). This does notseem to be a ma- which is nothing but the trigonometric (angle) addition
jor drawback because numerical values are sufficient for theorem. For integer j, the initial values are cos 0 = 1
most of practical purposes. and sin 0 = 0. For half integer j, one has to calculate,
I change the precision of floating-point numbers de- first,the initialvaluescos θ andsin θ and,second,cos θ
2 2
pending on j in such a way that the number of digits and sin θ using identities cos θ = cos2 θ sin2 θ and
equals the common logarithm of the value of Eq. (23) 2 − 2
divided by 10−18. It increases with j, reaching 74 digits sin θ =2sin 2θ cos 2θ.
forj =100. Precise64-bitfloating-pointnumberscanbe
obtainedsimplybytruncatingthe high-precisionresults. straightforward, max
Empirically, the time necessary to compute all the co- 12 recursion, max
straightforward, rms
efficientstjmk foreachjincreasesasj4. Thenewmethod ) recursion, rms
ν 14
takes 43 h for j = 100 with a personal computer with a −
0
CPU Intel core-i7 3960X running at 3.3GHz, using one (18
physical core. r (cid:12)(cid:12)
o
I use the obtained 64-bit floating-point number coeffi- r
r
e
cients to evaluate the Fourier-series formula for d func- 4(cid:12)(cid:12)
tions. I provide data files of the numerical values of the
coefficients, together with a sample FORTRAN90 pro-
gram to read the data and calculate the values of the d 0
0 20 40 60 80 100
function, as Supplemental Material [7].
ν
Foreachvalueofj, there are(2j+1)2 possiblecombi-
nations of the values of m and k (because j m j,
− ≤ ≤ FIG. 1. Maximum and root-mean-square errors in the val-
j k j). Only about a quarter of them are in-
− ≤ ≤ ues of cosνθ and sinνθ, numerically calculated with 64-bit
dependent, however, because of the properties of the d
floating-point numbers, over 0 ≤ θ ≤ π. The values of θ are
function expressed by Eqs. (A4), (A6), and (A7) in the sampledwithauniformspacingof10−5 degree. Theabscissa
appendix A. Therefore, I consider only such combina- isν =0,1,1,···,100. RecursionmeanstheusageofEq.(27)
2
tions as m 0 and k m in the following analysis of togetherwiththeinitialvaluesdescribedinthefollowingsen-
≥ ≤ | |
the numerical precision. tences.
Inotherwords,thecoefficientstjmk havethefollowing
ν
symmetry, Straightforward evaluation of cos νθ and sin νθ (i.e.,
passing the value of νθ to the internal functions cos and
tjkm =( 1)m−ktjmk, (24)
ν − ν sin)requiresaboutj callstothefunctions andiscompu-
tj,−m,−k =( 1)m−ktjmk, (25) tationally very inefficient. Moreover,as shown in Fig. 1,
ν − ν
4
suchstraightforwardevaluationcausesslightlylargernu- log error ofd40 (60◦)
10 mk
mericalerrorsprobablyduetothelossofsignificantdigits 40 -6
35 (cid:12) (cid:12)
inreducingthevalueofνθ toanintervalsuchas[ π,π], (cid:12) (cid:12) -8
especiallywhenthevalueof νθ islarge. Thereas−on4w4hy 2350 -10
the recursion formula (27) d|oes|not suffer from large er- m 20 -12
15 -14
rors even after hundred steps may be attributed to the
10 -16
fact that the magnitudes of cos and sin functions are al- 5 (a) Wignerformula -18
ways not greater than one. 0 -20
40
-6
35
-8
30
IV. PRECISION OF THE NUMERICAL 25 -10
VALUES OF THE d FUNCTION m 20 -12
15 -14
10 -16
5 (b) Fourierseries -18
log error ofd40 (30◦) 0 -20
10 mk
40 -40 -30 -20 -10 0 10 20 30 40
-6
35 (cid:12) (cid:12)
(cid:12) (cid:12) -8 k
30
-10
25
m 20 -12 FIG. 3. The same as in Fig. 2 but for θ=60◦.
15 -14
10 -16
5 (a) Wignerformula -18 log error ofd40 (90◦)
0 -20 40 10 mk
-6
3450 -6 3305 (cid:12)(cid:12) (cid:12)(cid:12) -8
30 -8 25 -10
25 -10 m 20 -12
m 20 -12 15 -14
15 -14 10 -16
10 -16 5 (a) Wignerformula -18
5 (b) Fourierseries -18 0 -20
0 -20 40
-6
-40 -30 -20 -10 0 10 20 30 40 35
-8
k 30 -10
25
m 20 -12
FIG. 2. The common logarithm of the error in the 64-bit 15 -14
floating-point numerical value of the d function dj (θ) as a 10 -16
functionofmandk forj =40andθ=30◦. Thevmalkueofthe 5 (b) Fourierseries -18
0 -20
d function is calculated by means of the Wigner formula in
-40 -30 -20 -10 0 10 20 30 40
part(a) andtheFourierseries expression inpart(b). Onlya
quarter of the possible combinations of m and k are plotted k
because the error is symmetric in the lines m = ±k. Errors
smaller than 10−20 are painted with the same color as that FIG. 4. The same as in Fig. 2 but for θ=90◦.
for the error of 10−20.
In this section, I compare the errors of the values of dently recommend the Fourier series expression over the
dj (θ) calculatedwith 64-bitfloating-pointnumbers ac- Winger formula already at j 40.
mk ∼
cording to the Wigner formula (3) andthe Fourierseries For a special purpose, however, the Wigner formula
expression (10). The errors have been calculated as the still has an advantage. For regions of the arguments θ
◦ ◦ ∼
differencesfromtheexactvaluescalculatedbyapplyinga 0 (180 ) and m + k 0 (2j), the Wigner formula
formula-manipulation software MAXIMA to the Wigner has smaller erro|rs than| 1∼0−15, i.e., than the level of the
formula. almostconstanterrorofthe Fourierseriesexpression. In
Figures 2, 3, and 4 show the errors for θ = 30◦, 60◦, suchregionsofthearguments,themagnitudeofdj (θ)is
mk
and90◦,respectively. Theerrorsareexpressedasafunc- very small and a small number of terms dominate in the
tion of (m,k) while j is fixed at 40. I have found that summation of the Wigner formula, while many terms of
similar plots for θ > 90◦ look almost indistinguishable order 1 cancel among themselves to give the small value
from those for 180◦ θ except that the sign of k is re- in the Fourier series expression. Therefore, if one needs
−
versed, as could be foreseen from Eq. (A9). such high precision at those regions, it would be a good
One can see that the Wigner formula results in very idea to develop a program which switches between the
large errors, up to 10−5 for θ 90◦ and m k 0, two formulas depending on the values of θ, j, m, and, k.
while the Fourierse∼riesexpression∼givesprecisel∼y15∼dig- In Fig. 5, I compare the maximum errors of dj (θ)
mk
its irrespective of the values of θ, j, m, and k. I confi- as functions of j. The maximum is taken over the val-
5
15 Theirnumericalvalues,whicharepreciseasfaras64-bit
Wigner formula
θ) 10 Fourier series floating-pointnumberscanexpress,areprovidedas elec-
( tric files in Supplemental Material. Sample programs in
k
m
jd 5 FORTRAN90 to use the data files are also provided.
n
i
r 0(cid:12)(cid:12)
o
r
er -5 Appendix A: Properties of the d function utilized in
(cid:12)(cid:12)
10 this paper
og-10
l
x
ma-15 I enumerate the symmetries of the d function to be
referred to in this paper.
-20
0 20 40 60 80 100 First, the unitarity of rotations (i.e, the Hermite con-
j jugate operator is the inverse operator) means
dj (θ)=dj ( θ). (A1)
FIG.5. Thecommon logarithm ofthemaximumerrorinthe mk km −
numerical value of the d function dj (θ) versus j calculated Second,thecompositionoftworotationscanberewritten
mk
with the Wigner formula (blue solid line) and the Fourier as multiplication of matrices to represent them,
series expression (red dashed line). The maximum is taken
overthe valuesof m,k, and θ for each j. j
dj (θ +θ )= dj (θ )dj (θ ). (A2)
mk 1 2 mν 1 νk 2
ν=−j
X
◦ ◦ ◦
ues of θ from 0 to 180 with an increment of 5 and all
Ineedthreemorerelations,whoseeasiestderivationmay
possible combinations of m and k. For both formulas,
be to use the Wigner formula (3) as in Ref. [1],
the maximum error increasesas exponential functions of
j. The Wigner formula increases its error, however, far dj (θ)=( 1)m−kdj ( θ), (A3)
fasterthantheFourierseriesexpression. Fortheinterval mk − mk −
20 j 100, the error of the Wigner formula can be djmk(θ)=dj−k,−m(θ), (A4)
≤ ≤
aopfpthroexFimouartieedrbseyrileosg1e0x|perrersosri|o≈n c0a.2n9b4je−ap1p7r.o2xwimhailteedthbayt djmk(π)=(−1)j+mδm,−k. (A5)
log error 0.006j 14.8. In other words, by increas- From Eqs. (A1) and (A3), one can prove
10| |≈ −
ing j by one,the errorofthe Wignerformulais doubled,
while that of the Fourier series expressionincreases only dj (θ)=( 1)m−kdj (θ), (A6)
mk − km
by 1.4%.
from Eqs. (A4) and (A6),
V. SUMMARY dj (θ)=( 1)m−kdj (θ), (A7)
mk − −m,−k
and from Eqs. (A2), (A5), and (A6),
I have shown that the Wigner formula for the d func-
tionresults inintolerablelargenumericalerrorsfor large
dj (π+θ)=( 1)j+mdj (θ), (A8)
valuesoftheangularmomentumquantumnumberj. On mk − −m,k
the other hand, the Fourier series expression for the d
function is shown to be free of such errors, providing
precision of ∼ 10−14 even at j = 100. An analytic ex- djmk(π−θ)=(−1)j+mdjm,−k(θ). (A9)
pressionfor the coefficients of the Fourier series is given.
[1] M.Rose,ElementaryTheoryofAngularMomentum(John J. Chem. Phys.111, 8825 (1999).
Wiley and Sons, New York,1957). [6] H. Dachsel, J. Chem. Phys.124, 144115 (2006).
[2] A.BohrandB.Mottelson, Nuclear Structure vol. I (Ben- [7] See SupplementalMaterial at http://link.aps.org/
jamin, New York,1969). supplemental/10.1103/PhysRevC.91.014320 for the nu-
[3] N.Tajima,Journal of Physics: Conference Series 445, 012014 (20m13e)r.ical values of the coefficients of the Fourier series of
[4] N. Tajima and T. Otsuka, thed functions and sample programs to use thedata.
Phys. Rev.C 84, 064316 (2011).
[5] C. Choi, J. Ivanic, M. Gordon, and K. Ruedenberg,