C++ Clebsch-Gordon Coefficients

Clebsch-Gordon Coefficients.

About

This program computes Clebsch-Gordon Coefficients in a self-contained, minimal C++ program.

Wikipedia (currently, as of 201606) describes Clebsch-Gordon Coefficients as:

... [appearing] as the expansion coefficients of total angular momentum eigenstates in an uncoupled tensor product basis. In more mathematical terms, [they] are used in representation theory, particularly of compact Lie groups, to perform the explicit direct sum decomposition of the tensor product of two irreducible representations (i.e., a reducible representation) into irreducible representations, in cases where the numbers and types of irreducible components are already known abstractly.

Results are output in exact form (square roots of rational numbers) and in either an on-screen (human parseable, computer non-parseable) format, or in LaTeX code (computer parseable, human non-parseable.)

Rationale

As mentioned in the README, "This program was written (from scratch, no cheating, no parents helping me glue glitter or anything)." I created this program because I was bored with a Quantum Mechanics assignment in 2011 for Phys 500 at the University of British Columbia. Quick backstory: I had already taken several courses on Quantum Mechanics and Field Theory (during my BSc and first year of MSc) and both the logical successor courses to Phys 500 at UBC (Phys 526 - Quantum Electrodynamics and Phys 508 - Quantum Field Theory). I was initially able to skip the Phys 500 prerequisite when I started my MSc (in theory), but I ended up having to take it anyways after later switching to Medical Physics. Quantum Mechanics is my favourite subject, so I spent much of the course trying to challenge myself. Ok -- back to the assignment.
The Phys 500 assignment asked for some specific Clebsch-Gordon coefficients to be laboriously computed. As a general theme in my Quantum Mechanics and Field Theory courses, instead of just doing the assignments I liked to write software to provide more general results, and then use them to solve the specifically assigned problems. At this point, I was toying with the idea of writing a general Computer Algebra System (CAS), and so decided it would be neat to write a very basic CAS to generate Clebsch-Gordon Coefficients.

This program is not very efficient, and can be quite slow for large coefficients. There are faster ways to compute Clebsch-Gordon coefficients, but I've never needed any that this program could not compute easily. Furthermore, Clebsch-Gordon Coefficients are often not needed explicitly, because one can work out the decomposition indirectly. It is nice, however, to spend some time in pure mathematics every once in a while.

The code should be understandable. There are inline comments (examples below) showing what is happening. I have no immediate plans to continue development, and am not currently using it for research.

Download

The source is available here, and is released under GPLv3, and FDLv3 licenses. Please send questions or comments to . Or, even better, send a pull request on GitHub ☺.

Examples

Here are a listing of some coefficients for up to spin-4 particles:


Here are some of the code comments:

    ...

    //Step 1) Determine the maximum m1 value possible.
    frac m1_max = j1.abs();

    //Step 2) Determine the maximal m2 value with m1 at maximum.
    frac m2_max = (j - m1_max);

    //Step 3) Determine the minimum m1 value possible.
    frac m1_min = -j1.abs();

    //Step 4) Determine the minimal m2 value with m1 at minimum.
    frac m2_min = (-j - m1_min);

    //Step 5) Set the value of the (m1_max,m2_max) CGC coefficient
    //  to be 1. No shenanigans here with m (m = m1 + m2.)
    std::vector<frac> arb;
    arb.push_back( frac(1,1) );
    TheCoeffs[CGCcoeff(j1, j2, m1_max, m2_max, j, m1_max+m2_max)] = arb;

    //Step 6) Starting at (m1_max,m2_max-1), use the J- relation to
    // determine the coefficient (m1_max,m2_max-1) in terms of 
    // (m1_max,m2_max) and a known (or zero) term to the right.
    // If one cannot be determined, indicate so.
    for(frac m2 = (m2_max-1); m2 >= -j2; --m2){
        std::vector<frac> PreFactorsA;
        //   J- equation: What are A,B,C factors?
        //
        //    m2|   C 
        //      |    o (m1, m2+1)  x            x
        //      |    |
        //      |    |
        //      |    |
        //      |    |
        //      |    |
        //      |    |A             B
        //      |    o-------------o            x
        //      | (m1,m2)       (m1+1,m2)
        //      |    
        //   ---|----------------------------------
        //      |                          m1
        //      |
        frac m = m1_max + m2 + 1; //Comes from the use of J- on the state.
        frac A = (j  + m )*(j  - m  + 1);
        frac B(0,1); // = (j1 - m1_max)*(j1 + m1_max + 1); //Always zero here.
        frac C = (j2 - m2)*(j2 + m2 + 1);
        if(A == 0){
            //If A is ever zero, then we have an issue. Report it and then sepuku!
            std::cerr << "Found a troublesome zero in step 6. Tell the human to fix me." << std::endl;
            exit(0);
        }
        //Get the prefactors from the point above (the point denoted C.)
        std::vector<frac> PreFactorsC = TheCoeffs[ CGCcoeff(j1,j2,m1_max,m2+1,j,m) ];
        PreFactorsA.push_back( (C/A)*PreFactorsC[0] );  //Only one number to worry about here.
        if(PreFactorsA.size() != 0)  TheCoeffs[ CGCcoeff(j1,j2,m1_max,m2,j,m-1) ] = PreFactorsA;
    }

    //Step 7) Use the J+ relation to determine (m1_max-1,m2_max) in
    // terms of (m1_max,m2_max) and (m1_max,m2_max-1).
    //Step 8) Perform the previous step but with m1_max decremented. 
    // Continue this until m1_max-i == m1_min.
    for(frac m1 = (m1_max);      m1 > m1_min; --m1)
    for(frac m2 = (m2_max); (m2).abs() <= j2; --m2)
    if( (-j <= (m1-1+m2).abs()) && ( (m1-1+m2).abs() <= j ) ){
        std::vector<frac> PreFactorsB;
        //   J+ equation: What are A,B,C factors?
        //
        //    m2|   B              A
        //      |    o-------------o (m1,m2)    x
        //      |  (m1-1,m2)       |
        //      |                  |
        //      |                  |
        //      |                  |
        //      |                  |
        //      |                C |
        //      |    x             o            x
        //      |              (m1,m2-1)
        //      |    
        //   ---|----------------------------------
        //      |                          m1
        //      |
        frac m = m1 + m2 - 1; //Comes from the use of J+ on the state.
        frac A = (j - m)*(j + m + 1);
        frac B = (j1 + m1)*(j1 - m1 + 1);
        frac C = (j2 + m2)*(j2 - m2 + 1);
        if(B == 0){
            //If B is ever zero, then we have an issue. Report it and then sepuku!
            std::cerr << "Found a troublesome zero in step 7. Tell the human to fix me." << std::endl;
            exit(0);
        }
        //Get the prefactors from the local point (point A) and the one below (point C)
        // so we can determine the prefactor for the point to the left (point B.)
        std::vector<frac> PreFactorsA = TheCoeffs[ CGCcoeff(j1,j2,m1,m2,j,m+1) ];
        std::vector<frac> PreFactorsC;
        if( (m2-1).abs() <= j2 ) PreFactorsC = TheCoeffs[ CGCcoeff(j1,j2,m1,m2-1,j,m) ];

        for(size_t i=0; i<PreFactorsA.size(); ++i)  PreFactorsB.push_back( (A/B)*PreFactorsA[i] );
        for(size_t i=0; i<PreFactorsC.size(); ++i)  PreFactorsB.push_back( -(C/B)*PreFactorsC[i] );
        if(PreFactorsB.size() != 0)  TheCoeffs[ CGCcoeff(j1,j2,m1-1,m2,j,m) ] = PreFactorsB;
    }
    ...