Numerical calculation of the Clebsch-Gordan Coefficients

Examples – Quantum Mechanics

Last edited: January 2nd 2021


This notebook is actually not written by us from the NumFys-team, instead it is written by André Flakke, and after some guidance from us, we are happy to publish André's work! He will demonstrate how to calculate the Clebsch-Gordan coefficients, coefficients arising when dealing with addition of quantum mechanical angular momentum. If you have any questions to André, he can be reached at [email protected]


Motivation for Computing the Clebsch-Gordan Coefficients

The motivation for computing Clebsch-Gordan coefficients is mainly that they answer the question: If I pair two angular momentum states, how does the new angular momentum state relate to the original, individual states? The emergent two-particle state can be written as a linear combination (superposition) of states, and we are interested in the coefficients of this linear combination. These coefficients are known as the Clebsch-Gordan (CG) coefficients.

For concreteness, the term angular momentum states can be thought of as states describing particles with angular momentum. The angular momentum can be either orbital momentum, or the intrinsic orbital momentum (spin), carried by most elementary particles.

The CG coefficients are computable by hand, normally through the use of ladder operators, but the computation quickly becomes time consuming as the total angular momentum increases. Thus, it is more common to look them up in a table or calculate them numerically. In this notebook, we will demonstrate how to do the latter.

Online tables of CG coefficients are for example the Wikipedia.org web page [1] or Particle Data Group's [2]. We will use these tables to verify our results.

Addition of angular momenta and the CG Coefficients

We will not cover the details of addition of angular momentum in this notebook. We recommend [3] (English) or [4] (Norwegian) for a gentle introduction on the subject. However, we will cover the most basic properties needed to understand the CG coefficients. The total angular momentum vector of a particle, $\textbf{J}$, can be written as the vector sum $\textbf{J} = \textbf{L} + \textbf{S}$, where $\textbf{L}$ denotes the orbital angular momentum and $\textbf{S}$ denotes the intrinsic angular momentum. One can furthermore add the total angular momentum vectors of two particles $\textbf{J}_{1} + \textbf{J}_{2} = \textbf{J}$. This restricts the possible values for the total angular momentum for the new state. In particular, this enforces the quantum number $j$ (corresponding to $\textbf{J}$), to take the values $|j_{1} - j_{2}| \leq j \leq j_{1} + j_{2}$, with unit integer increments.

Before we define the CG coefficients, we need to make a short detour and introduce the total angular momentum projection quantum number $m$. It is instructive to first consider a single particle with orbital momentum $l_1$, spin $s_1$, and thus $j_1 \in \{|l_1 - s_1|, |l_1 - s_1| + 1, ..., l_1 + s_1\}$. The possible values for the magnetic quantum number $m_1$ is then given by $m_1 \in \{ -j_1, -j_1 + 1, ..., j_1 - 1, j_1 \}$. Physically, one measures $m_1$ by projecting the angular momentum on a certain axis, conventionally along the $z$-axis. These relations also hold for the two-particle example described above, such that for each possible $j$, we have a corresponding set of $m$-s given by $m\in \{ -j, -j + 1, ..., j - 1, j \}$. $m$ must also fulfill $m = m_{1} + m_{2}$. Note that quantum numbers with subscripts refer to a specific particle, whereas the quantum numbers without subscripts denote the new two-particle state.

We are now ready to write the expression for the CG-coefficients. In bra-ket notation the equation relating the two-particle state $|j m \rangle$ in terms of the initial angular momentum states, $| j_{1} m_{1} \rangle | j_{2} m_{2} \rangle$, is given by

\begin{equation} | j m \rangle = \sum_{m_{1},m_{2}}\langle j_{1} m_1 j_{2}m_{2} | j m \rangle | j_{1} m_{1} \rangle | j_{2} m_{2} \rangle, \label{CG_def} \end{equation}

where the CG-coefficients are defined as the real valued coefficient $\langle j_{1}m_{1}j_{2}m_{2} | j m \rangle$. From equation \eqref{CG_def}, the physical interpretation of the CG coefficients reemerges; They describe how much of each old state is contained in the new states. It is thus immediately clear that the square of the CG coefficients represent probability amplitudes.

Example: Addition of Angular Momentum of Two Electrons

Before we start with the numerical implementation, let's consider a simple example in order to familiarize us with the matrix notation that we will use later. Additionally, this will allow us to visualize how the bra-ket notation relates to the matrix notation.

Consider adding the angular momentum of two electrons with $l_{1} = l_{2} = 0$ and $s_{1} = s_{2} = 1/2$. Physically, these are electrons described by $s$-orbitals. We then have $j_{1} = j_{2} = 1/2$, and hence $m_{1},m_{2} \in \{-1/2, 1/2\}$. Thus, the total angular momentum for the new two-particle state, $j$, fulfills $0 \leq j \leq 1$. And since $j$ increases in integer steps, $j \in \{0,1\}$.

In the following, we again refer to [3] and [4] for a detailed derivation for the the two-particles states, here we will focus on the results. Using kets to denote quantum numbers, we have a singlet normalized product state with $j = 0$ by the combination

\begin{equation} | j=0, m=0 \rangle = {1 \over \sqrt{2}} \big ( | j_1=1/2, m_1=1/2 \rangle | j_2=1/2, m_2=-1/2 \rangle - | j_1=1/2, m_1=-1/2 \rangle | j_2=1/2, m_2=1/2 \rangle \big ), \end{equation}

or more compactly written as

\begin{equation} | 0, 0 \rangle = {1 \over \sqrt{2}} \big ( | 1/2, 1/2 \rangle | 1/2, -1/2 \rangle - | 1/2, -1/2 \rangle | 1/2, 1/2 \rangle \big ), \label{singlet} \end{equation}

and a triplet states with $j = 1$

\begin{equation} \begin{pmatrix} | 1, -1 \rangle \\ | 1, 0 \rangle \\ | 1, 1 \rangle \end{pmatrix} = \begin{pmatrix} | 1/2, -1/2 \rangle | 1/2, -1/2 \rangle \\ {1 \over \sqrt{2}} \big ( | 1/2, 1/2 \rangle | 1/2, -1/2 \rangle + | 1/2, -1/2 \rangle | 1/2, 1/2 \rangle \big ) \\ | 1/2, 1/2 \rangle | 1/2, 1/2 \rangle \end{pmatrix}, \label{triplet} \end{equation}

where we omitted to write all the labels for the quantum numbers in the interest of brevity.

By comparing equations \eqref{singlet} and \eqref{triplet} to the original expression for the CG coefficients in equation \eqref{CG_def} we see that the coefficients in front of each product of single particle states are indeed the CG coefficients. In matrix form, we can write this as

\begin{equation} \begin{pmatrix} | 0, 0 \rangle \\ | 1, -1 \rangle \\ | 1, 0 \rangle \\ | 1, 1 \rangle \end{pmatrix} = \begin{pmatrix} 0 & {1 \over \sqrt{2}} & -{1 \over \sqrt{2}} & 0 \\ 1 & 0 & 0 & 0 \\ 0 & {1 \over \sqrt{2}} & {1 \over \sqrt{2}} & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} | 1/2, -1/2 \rangle | 1/2, -1/2 \rangle \\ | 1/2, 1/2 \rangle | 1/2, -1/2 \rangle \\ | 1/2, -1/2 \rangle | 1/2, 1/2 \rangle \\ | 1/2, 1/2 \rangle | 1/2, 1/2 \rangle \end{pmatrix}, \label{CG_matrix} \end{equation}

where the CG-coefficients are contained in the matrix. In this notebook, we will numerically compute this matrix for general $j_1$ and $j_2$.

Computing the CG Coefficients

As mentioned in the first section, computing the CG-coefficients numerically might seem strange at first, but as the complexity of the calculations grow, a numerical routine will be a useful aid up to the representation error ([5], 15.1) of the programming language.

There is a wide variety of ways to compute and represent the CG coefficients in exact arithmetic (i.e without rounding error), see e.g. [6]. We will use the expression \begin{equation} \langle j_{1}m_{1}j_{2}m_{2}|jm \rangle = { \delta_{m,m_{1} + m_{2}} \over \Delta(j_{1},j_{2},j) } \sqrt{ (2j + 1) } \sqrt{ { (j_{1} + m_{1})!(j_{1} - m_{1})!(j + m)!(j - m)! \over (j_{2} + m_{2})!(j_{2} - m_{2})! } } \times \sum_{z = 0}^{\min(j - m, j_{1} - m_{1})}\mathcal{X}(z), \label{CG_num} \end{equation}

where

$$ \mathcal{X}(z) = { (-1)^{j_{1} - m_{1} - z} \over z!} { (j_{1} + j_{2} - m - z)!(j_{2} + j - m_{1} - z)! \over (j_{1} - m_{1} - z)!(j - m - z)!(j_{1} + j_{2} + j + 1 - z)! }, $$

and $$\Delta(j_{1},j_{2},j) = \sqrt{ { (j_{1} + j_{2} - j)!(j_{1} - j_{2} + j)!(-j_{1} + j_{2} + j)! \over (j_{1} + j_{2} + j + 1)! } }.$$

For convenience, in our numerical routines we will label the prefactor in equation \eqref{CG_num} as

$$c := { 1 \over \Delta(j_{1},j_{2},j) } \sqrt{ (2j + 1) } \sqrt{ { (j_{1} + m_{1})!(j_{1} - m_{1})!(j + m)!(j - m)! \over (j_{2} + m_{2})!(j_{2} - m_{2})! } }.$$

Numerical implementation

The Python code contains four functions. Three of them are auxiliary functions which are called from the main function. Before we start going through the numerical routine, it is useful to remind ourselves that the angular momenta are either integer or an integer divided by two (half-integer). The routine is also only defined for half-integer or integer. There is however, nothing wrong with having a combination of half-integer and integer quantum numbers as input to the program code.

The first code block imports libraries and mathematical functions:

In [1]:
import numpy as np
import pandas
from math import factorial as fac
from math import sqrt 

The use of pandas is solely cosmetic, and no familiarity with the library is required to understand the code. Now we can look at the first auxiliary function.

In [2]:
def get_levels(j_i): 
    
    m_i = np.arange(-j_i,j_i + 1,1)  #quantum numers ms between -j_i and j_i in unit integer steps 
    n_i = len(m_i)
    
    J_i = np.zeros((n_i,2))         #allocate a zero-filled array
    
    for i in range(0,n_i):         #fill the array in a loop
        J_i[i][0] = j_i
        J_i[i][1] = m_i[i]
    
    return J_i,n_i;

This function computes and returns the $n_i \times 2$ table $\textbf{J}_i$, containing the different $m_i$-values one can have for the input $j_i$. The index $i$ referes to one of the inputs $j_1$ or $j_2$. $n_i$ is given by $n_i := (2j_i + 1)$. In $\textbf{J}_i$, the first column is simply $j_i$, and the second column of the array contains integer values $-m_i \leq j_i \leq m_i$ in unit integer increments (i.e. 1).

Let's illustrate with the example $j = 3/2$

In [3]:
print(get_levels(3/2))
(array([[ 1.5, -1.5],
       [ 1.5, -0.5],
       [ 1.5,  0.5],
       [ 1.5,  1.5]]), 4)

In the second auxiliary function, we use pandas to increase the readability of the output

In [4]:
def make_product_state_labels(J_1,J_2,n_1,n_2):
        
    r_1 = np.kron(np.ones((n_2,1)),J_1)  #kronecker tensor product
    r_2 = np.kron(J_2,np.ones((n_1,1)))  #kronecker tensor product
    
    J_12 = np.concatenate((r_1,r_2),axis = 1) #stacks column vectors r_1, r_2 along the column dimension 
    
    state_list = ['| j_1,','m_1 >','| j_2,','m_2 >']    #labels for the individual quantum numberinputs 
    J_12 = pandas.DataFrame(J_12,None, state_list)    #make a table for readability
    
    return J_12;    

As the name suggests, this function makes labels for the product state described by the two initial sets of quantum numbers. Hence it is a matrix containing labels for the possible products of $ |j_{1},m_{1} \rangle |j_{2},m_{2} \rangle $ for all possible values of $m_1$ and $m_2$. $j_{1}, j_{2}$ are fixed input parameters to the function.

The inputs $n_{1}, n_{2}$ to this function are used to construct the $(n_{1}n_{2}) \times 2$ matrices $\textbf{r}_{1} := \textbf{1}_{2} \otimes \textbf{J}_{1}$ and $\textbf{r}_{2} := \textbf{J}_{2} \otimes \textbf{1}_{1}$ in terms of the Kronecker tensor product $\otimes$ (which you can read up on here [7]).The vectors $\textbf{1}_{1}$ and $\textbf{1}_{2}$ consists of ones, with length $n_{1}$ and $n_{2}$, respectively. This ensures that the columns of the matrix $\textbf{J}_{12} = (\textbf{r}_{1}, \textbf{r}_{2})$ have equal length. The matrix $\textbf{J}_{12}$ will have column-dimension equal to $4$, where the two first columns label the states $ |j_{1},m_{1} \rangle $ and the two last columns label the states $|j_{2},m_{2} \rangle$. The matrix is turned into a table (or technically a dataframe) before it is returned from the function.

Hence the output from this function is only used to label states and hence keep track of Clebsch-Gordan coefficients that are not zero. Let's see how it works with the two-electron example $j_{1} = j_{2} = 1/2$

In [5]:
J_1,n_1 = get_levels(1/2)
J_2,n_2 = get_levels(1/2)
print(make_product_state_labels(J_1,J_2,n_1,n_2))
   | j_1,  m_1 >  | j_2,  m_2 >
0     0.5   -0.5     0.5   -0.5
1     0.5    0.5     0.5   -0.5
2     0.5   -0.5     0.5    0.5
3     0.5    0.5     0.5    0.5

Thus, we reproduce the rightmost vector in equation \eqref{CG_matrix}, as desired.

The next auxiliary function we will look at is the one that computes a particular CG-coefficient using equation \eqref{CG_num}.

In [6]:
def cg_coef(j_1,j_2,j,m_1,m_2):
    # This function computes a particular matrix entry, 
    # corresponding to a single Clebsch-Gordan coefficient
    
    m = m_1 + m_2 # delta function constraint 
    
    Delta = sqrt( fac(j_1 + j_2 - j)*fac(j_1 - j_2 + j)*fac(-j_1 + j_2 + j)/fac(j_1 + j_2 + j + 1)  )
    
    #c becomes a numerical prefactor in the formula  
    c = sqrt( fac(j_1 + m_1)*fac(j_1 - m_1)*fac(j + m)*fac(j - m) / (fac(j_2 + m_2)
                                                                     *fac(j_2 - m_2)) ) * sqrt(2*j + 1) / Delta
    
    def chi(z):
        numerator = fac(j_1 + j_2 - m - z)*fac(j_2 + j - m_1 - z)
        denominator = (fac(j_1 - m_1 - z)*fac(j - m - z)*fac(j_1 + j_2 + j + 1 - z))
        return numerator/denominator
    
    temp = 0
    zmax = min(j - m,j_1 - m_1)
    
    #sum over values consistent with 'zmax'
    
    for z in range(0,int(zmax + 1)):
        temp += (-1)**(j_1 - m_1 + z)*chi(z)/fac(z)
    
    return c*temp

Once this function has been called, the return value will be the particular instance of $\langle j_{1}m_{1}j_{2}m_{2}|jm \rangle$.

Now it is time to look at the main function which computes a matrix of CG-coefficients, i.e. the matrix in the matrix equation from the theory. We denote this matrix as $\textbf{C}$. The input parameters to this function are the quantum numbers $j_{1}, j_{2}$. For two electrons, these quantum numbers would be $j_{1} = 1/2, j_{2} = 1/2$.

In [7]:
def clebsch_gordan(j_1,j_2):
    # This is the main function computing the matrix of Clebsch-Gordan coefficients
    
    jmin = abs(j_1 - j_2)    # min. imum allowed total angular momentum
    jmax = j_1 + j_2         # maximum allowed total angular momentum
    
    increment = np.arange(jmin, jmax + 1,1)       #increment in integer steps
    
    # loop over allowed quantum numbers
    for k in range(0,len(increment)):
        
        if(k == 0):
            J,n = get_levels(increment[k])
            N_list = [n]
        else:
            jtemp,n = get_levels(increment[k])
            N_list.append(n)
            J = np.concatenate((J,jtemp))
        
    J_1,n_1 = get_levels(j_1)        
    J_2,n_2 = get_levels(j_2)

    N = sum(N_list)
    C = np.zeros((N,n_1*n_2))  # initialize matrix to hold the CG-coefficients to zero
    
    # loop over all possible quantum numbers 
    for m in range(0,N):
        for m_1 in range(0,n_1):
            for m_2 in range(0,n_2):
                if(J_1[m_1][1] + J_2[m_2][1] == J[m][1]):
                    C[m][m_1 + n_1*m_2] = ( 
                    cg_coef(J_1[m_1][0],J_2[m_2][0],J[m][0],J_1[m_1][1],J_2[m_2][1]))
    
    # the matrix C now contains the CG-coefficients, also the ones that are zero 
    
    final_states = ['| j,',' m >'] #labels the final states inside the main function
    J = pandas.DataFrame(J,None, final_states)
    
    J_12 = make_product_state_labels(J_1,J_2,n_1,n_2) 
    return J, J_12, C;

Note that the function returns all three quantities involved in the matrix equation from the theory. Firstly the vector on the left hand side is denoted $\textbf{J}$, and contains the labels for the states $| j m_{i} \rangle$ such that each row of $\textbf{J}$ is a vector with entries $[j, m_{i}]$. Additionally the right hand side quantities in the matrix equation, namely the matrix of Clebsch-Gordan coefficients and the original basis vector, are returned.

It should be mentioned at this point that only the returned matrix $\textbf{C}$ contains numerical values. The other return values are only used as labels to find a particular CG-coefficient.

Let's compare the numerical results for the two-electron case with the triplet and singlet found previously.

In [8]:
np.set_printoptions(precision=3) #truncate the display of decimals to three digits

J,J_12,C = clebsch_gordan(1/2,1/2) # two-electron model
print(J, "\n") # J is a vector of possible angular momentum states resulting from combining J_{1},J_{2} 
print(C, "\n") # C is the matrix of Clebsch-Gordan coefficients 
print(J_12, "\n") # J_12 is the product states that combine with coefficients taken from matrix C 
   | j,   m >
0   0.0  -0.0
1   1.0  -1.0
2   1.0   0.0
3   1.0   1.0 

[[ 0.     0.707 -0.707  0.   ]
 [ 1.     0.     0.     0.   ]
 [ 0.     0.707  0.707  0.   ]
 [ 0.     0.     0.     1.   ]] 

   | j_1,  m_1 >  | j_2,  m_2 >
0     0.5   -0.5     0.5   -0.5
1     0.5    0.5     0.5   -0.5
2     0.5   -0.5     0.5    0.5
3     0.5    0.5     0.5    0.5 

The result is identical to equation \eqref{CG_matrix}, indicating that the implementation is correct. In order to further confirm the validity of the implementation, we can try different values for $j_1$ and $j_2$.

For example using $j_{1} = 3/2, j_{2} = 1/2$ yields the output

In [9]:
J,J_12,C = clebsch_gordan(3/2,1/2)

print(J, "\n") 
print(C, "\n") 
print(J_12, "\n") 
   | j,   m >
0   1.0  -1.0
1   1.0   0.0
2   1.0   1.0
3   2.0  -2.0
4   2.0  -1.0
5   2.0   0.0
6   2.0   1.0
7   2.0   2.0 

[[ 0.     0.5    0.     0.    -0.866  0.     0.     0.   ]
 [ 0.     0.     0.707  0.     0.    -0.707  0.     0.   ]
 [ 0.     0.     0.     0.866  0.     0.    -0.5    0.   ]
 [ 1.     0.     0.     0.     0.     0.     0.     0.   ]
 [ 0.     0.866  0.     0.     0.5    0.     0.     0.   ]
 [ 0.     0.     0.707  0.     0.     0.707  0.     0.   ]
 [ 0.     0.     0.     0.5    0.     0.     0.866  0.   ]
 [ 0.     0.     0.     0.     0.     0.     0.     1.   ]] 

   | j_1,  m_1 >  | j_2,  m_2 >
0     1.5   -1.5     0.5   -0.5
1     1.5   -0.5     0.5   -0.5
2     1.5    0.5     0.5   -0.5
3     1.5    1.5     0.5   -0.5
4     1.5   -1.5     0.5    0.5
5     1.5   -0.5     0.5    0.5
6     1.5    0.5     0.5    0.5
7     1.5    1.5     0.5    0.5 

By cross-checking with [1], we confirm that the function calculates the correct CG coefficients. We invite the reader to try other examples.

Conclusion

We have made a program which allows us to calculate the CG coefficients for any pair of $j_1$ and $j_2$. With the functionality presented here, a generalization to systems of more than two particles is straightforward (see that you understand why!). For larger quantum numbers, an automated routine will save us a great deal of time when calculating the CG coefficients.

It is also worth noting the physical significance of these states. For fermions, we know that the total wave function for a system of fermions should be anti-symmetric with respect to interchanging coordinates (both spin and spatial) of two identical particles. Thus, if the spin part of the wave function is described by the singlet in equation \eqref{singlet}, it is clear that an interchange in the spin coordinates would pick up a minus sign, and hence the spatial part of wave function must be symmetric with respect to interchange. The reverse is true for triplets, requiring triplets to have an anti-symmetric spatial wave function. This kind of symmetry often has profound consequences in physics, with examples from e.g. particle physics or condensed matter physics. In other words, what we have studied in this notebook has important, real-world implications!

References

[1]: https://en.wikipedia.org/wiki/Table_of_Clebsch%E2%80%93Gordan_coefficients (Wikipedia; Online: Accessed 26/11/2020)
[2]: https://pdg.lbl.gov/index.html (Particle Data Group) [3]: Bransden, B. H., & Joachain, C. J. (2000). Quantum mechanics (2nd ed., pp. XIV, 803). Pearson/Prentice Hall.
[4]: Hemmer, P. (2005). Kvantemekanikk (5. utg. ed.). Trondheim: Tapir akademisk forl.
[5]: https://docs.python.org/3/tutorial/floatingpoint.html (Online: Accessed 18/11/2020)
[6]: https://functions.wolfram.com/HypergeometricFunctions/ClebschGordan/06/01/
[7]: https://en.wikipedia.org/wiki/Kronecker_product (Wikipedia; Online: Accessed 19/11/2020)