Shell model calculation on 18O

This post is copy from the book Theory Of The Nuclear Shell Model by R. D. Lawson, chapter 1.2.1

---

The model space is only the 0d5/2 and 1s1/2, and the number of valence nucleon is 2. The angular coupling of the 2 neutrons in these 2 orbitals are

Note that for identical particle, the allowed J coupled in same orbital must be even due to anti-symmetry of Fermion system.

The spin 3, and 4 can only be formed by (0d5/2)(1s1/2) and (0d5/2)^2 respectively.

Since the Hamiltonian commute with total spin, i.e., the matrix is block diagonal in J that the cross J matrix element is zero,

or to say, there is no mixture between difference spin. The Hamiltonian in matrix form is like,

The metrix element of J=3 and J=4 is a 1 × 1 matrix or a scalar.

where is the single particle energy.

Suppose the residual interaction is an attractive delta interaction

---

Be fore we evaluate the general matrix element,

We have to for the wave function ,

where is the isospin, and the single particle wave function is

Since the residual interaction is a delta function, the integral is evaluated at , thus the radial function and spherical harmonic can be pulled out in the 2-particle wave function at is

Using the product of spherical Harmonic,

using the property of Clebsch-Gordon coefficient for spin half system

where

which is equal to

For

With some complicated calculation, the J-J coupling scheme go to L-S coupling scheme that

with

---

Return to the matrix element

Since the matrix element should not depends on , thus, we sum on M and divide by ,

with

( i give up, just copy the result ), for ,

---

The block matrix are

To solve the eigen systems, it is better to find the from experimental data. The single particle energy of the d and s-orbtial can be found from 17O, We set the reference energy to the binding energy of 16O,

the ground state of 18O is

Solving the , the eigen value are

Thus,

The solution for all status are

---

The 2nd 0+ state is missing in above calculation. This is due to core-excitation that 2 p-shell proton promotes to d-shell.

In the sd- shell, there are 2 protons and 2 neutrons coupled to the lowest state. which is the same s-d shell configuration as 20Ne. The energy is

In the p-shell, the configuration is same as 14C, the energy is

Thus, the energy for the 2-particle 2-hole of 18O is

,

where is the p-sd interaction, there are 4 particle in sd-shell and 2 hole in p-shell, thus, total of 8 particle-hole interaction.

The particle-hole can be estimate using 19F 1/2- state, This state is known to be a promotion of a p-shell proton into sd-shell.

In the sd-shell of 19F, the configuration is same as 20Ne. In the p-shell of 19F , the configuration is same as 15N, the energy is

Thus, the energy for the 1/2- state of 19F relative to 16O is

And this energy is also equal to

Thus,

Therefore, the 2nd 0+ energy of 18O is

Compare the experimental value of 3.63 MeV, this is a fair estimation.

---

It is interesting that, we did not really calculate the radial integral, and the angular part is calculated solely base on the algebra of J-coupling and the properties of delta interaction.

And since the single particle energies and residual interaction are extracted from experiment. Thus, we can think that the basis is the “realistic” orbital of d and s -shell.

The spectroscopic strength and the wave function of the 18O ground state is  . In here the basis wavefunctions are the “realistic” or “natural” basis.

Construction of many-body wavefunction

This post is an extraction from “Theory of the nuclear shell model” by R.D. Lawson, Chapter 1.1.2

---

In a simple construction, the single-particle wave function is

And the many body Slater determinants is

where is a set of  n states from the model space. Thus, for given , there are many possible choices of .

And the Hamiltonian is expanded into the space of .

I thought this is the end of the story, until I read the Chapter 1.1.2

---

To continuous further, we need 2 theorems. The 1st one is that.

If , then, and

This is kind of trivial, as the Ladder operator cannot promote the m-state any further, the m-state must be maximized and equal to J.

Before go to the 2nd theorem, let introduce the creation operator , and denote the Slater determinant as

The 2nd one is a bit complicated.

In a given model space ( j_1, j_2, …., j_k) with n particles.  For a p distinct Slater determinants with .  Suppose there are another p+q distinct Slater determinants , so that , thus, q states with angular momentum J= M-1 can be constructed$.

Lets use an example to illustrate, say, we have a model space of 0d5/2 with 2 particles, we can construct the state as,

And for , the only possible state is

Since the number of state for is 1, and so . The number of state for is also 1, thus . Therefore, according to the theorem 2, there is no state for . And we knew that, from a 2 particles in d5/2 shell can only form J = 0, 2, 4, due to the anti-symmetric of Fermions system.

we can continuous, , there are 2 possible states,

Thus, this time, (as M=3 has 1 state) and , therefore, there is 1 state of can be constructed.

For ,

Thus, the number of state is 0.

For ,

Therefore,  the number of state is 1. In summary,

The number of J = 4 state = 1.
The number of J = 3 state = 0.
The number of J = 2 state = 1.
The number of J = 1 state = 0.
The number of J = 0 state = 1.

Another way to do so is draw a table

m,

4

3

2

1

0

2

1

0

-1

0

-1

-2

-2

-3

-4

5/2

x

x

x

x

x

3/2

x

x

x

x

x

1/2

x

x

x

x

x

-1/2

x

x

x

x

x

-3/2

x

x

x

x

x

-5/2

x

x

x

x

x

Ignore the negative , and we can deduce the number of state for .

The prove of the 2nd theorem is skipped. And I like to add an auxiliary theorem that for n identical particle in the j-orbital, the maximum spin is

This theorem can be proved by using the 1st theorem. The reason that this theorem is included is that, the anti-symmetry property implicitly contains in this theorem, and I see the elegant and beauty in it. Also, this also an example that theorem is trivial but useful.

---

Back to the wave function construction, the wave function,

is only a component, or part of the total wave function of spin J. For example, the coupling of 0d5/2 and 1s1/2 states, we can construct

But hold-on, we missed the coupling coefficient, without the coupling coefficient, we actually did not know the total J!! The proper wave function should be

Up to now, we should separate two wave function, and . And more important, we drop the J in .

— single particle orbital
— m-state
— many-body wave function with spin J

when only 2 particles, the is the Clebsch-Gordon coefficient.

---

We now give an example on (0d5/2)^2 configuration. For J = 4, only 1 state is possible, this

Using the Ladder operator, we can find all .

For J = 2, there are 2 m-states,

In order to find the coefficient , we have the normalization condition , and we have to use the theorem 1, and apply the operator and require that

Thus,

solve

In fact,

In fact, for 2 particle wave function

---

The power of this construction method is the many-body wave function. we now demonstrate 4 identical particle in 0d5/2 orbital,

The maximum J is

.

the M-state are

For

For

For

For

For

or, in table

m,

4

3

2

2

1

0

1

0

-1

-2

0

-1

-2

-3

-4

5/2

x

x

x

x

x

x

x

x

x

x

3/2

x

x

x

x

x

x

x

x

x

x

1/2

x

x

x

x

x

x

x

x

x

x

-1/2

x

x

x

x

x

x

x

x

x

x

-3/2

x

x

x

x

x

x

x

x

x

x

-5/2

x

x

x

x

x

x

x

x

x

x

From the theorem 2, the number of m-state tell us the number of state for J are

To construct state,

, using the same method

---

I suspect the coefficient is deeply related to n-j symbol. And I am sure the Group theory can provide an elegant solution and display for that.

At last, in the construction of m-state, we restricted the order, therefore, state like

will be discarded.

But in general 2-body wave function,

there is no restriction.

Two Spin-1/2 system

I did it for very long time, unfortunately, cannot find the posts (here, here) show the detail, especially in Mathematica.

---

First, we have to distinguish a direct-sum and a direct-product. A direct sum is

and a direct-product is

where is KroneckerProduct in Mathematica, and is identity operator or matrix. And

Note that (the lower dot) and (the middle dot) are different. The $late \cdot $ is matrix dot-product.

A directly-sum acts on each spin separately, and a direct-product acts on each spin at the same time. The different can be seen using ladder operator.

For example,

---

The total spin of two spin-1/2 system,

and are commute,

If using Mathematica, the calculation is straight forward, and the eigen-system will give out the spin and states in term of basics. Also, the Clebsch-Gordan coefficient. However, for human, we can rewrite using ladder operator.

The above expression is convenient in basis.

Here is a screenshot from Mathematica

Since the matrix of can be block diagonalized with with dimension 3 and 1, respectively. The matrix representation for arbitrary spin is done in this post.

 

The ladder operator for 2-spin system using basis, can be defined as,

It is because the ladder operator change the spin by flipping one at a time.

Nilsson Orbital using diagonalization method

Long time ago, I tried to tackle the Nilsson orbital by solving the Hamiltonian analytically. However, the Hamiltonian is without LS coupling. This times, I redo the calculation according to the reference B. E. Chi, Nuclear Phyiscs 83 (1966) 97-144.

---

The Hamiltonian is

using

The Hamiltonian becomes

Set , and

Set

and the perturbation is

with

For example, when A = 12, MeV, when A = 132, MeV

---

The wavefunction for the spherical harmonic is

---

The diagonal elements are

where

The off diagonal elements are

( I will evaluate this integral in future )

---

The rest is diagonalization the Hamiltonian

Here is the calculation for the 2nd harmonic for

The component of each orbital can be directly taken from the eigenvalue. Here is the [521]1/2 state.

Sum rule of Spectroscopic factors

For transfer reaction, the wave function of nucleus B = A + 1 can be written as

where the due to anti-symetrization between the single nucleon and the core nucleus A, the Clebsh-Gordon coefficient is come from the angular coupling, and is the square root of the spectroscopic factor for orbital between nuclei B and A.

We can see that, to find out the spectroscopic factor, we can integrate out the core nucleus A

In the last step, we write the Clebsh-Gordon coefficient in Winger 3j-symbol. Using the Wigner-Eckart theorem

Thus,

using the identity

we have

replacing the reduced matrix element in the Wigner-Eckart theorem with previous result,

multiple the adjoint

---

For A(d,p) B reaction, the nucleus A is in the ground state, we can sum all the and states, then sum all the , notice that

we have

recall

The last equality mean the number of nucleon in orbital in nucleus A. Thus, we obtain

---

For B(p,d)A reaction, the nucleus B is in the ground state and we can sum all ,

---

Thus, the adding A(d,p)B and removing A(p,d)C reaction from nucleus A, the sum of spectroscopic factors

 

 

3j, 6j, 9j-symbol for C++

Since I cannot find a odd c++ code ( < c++11 ) for the calculation of those j-symbol, I made one for myself!

The code just using the formula which you can find in wiki or here. So, the code is not optimized, but for simple calculation, I guess it does not matter. Also, I checked, in my pc, calculate 1000 9-j symbols only take 1 or 2 sec.

The program first build

1. factorial
2. Clebsch–Gordan coefficient
3. 3-j symbol
4. 6-j symbol
5. 9-j symbol

#include <stdlib.h>
#include <cmath>
#include <algorithm>

using namespace std;

double factorial(double n){
  if( n < 0 ) return -100.;
  return (n == 1. || n == 0.) ? 1. : factorial(n-1) * n ;
}

double CGcoeff(double J, double m, double J1, double m1, double J2, double m2){
  // (J1,m1) + (J2, m2) = (J, m)

  if( m != m1 + m2 ) return 0;

  double Jmin = abs(J1 - J2);
  double Jmax = J1+J2;

  if( J < Jmin || Jmax < J ) return 0;

  double s0 = (2*J+1.) * factorial(J+J1-J2) * factorial(J-J1+J2) * factorial(J1+J2-J) / factorial(J+J1+J2 + 1.);
  s0 = sqrt(s0);

  double s = factorial(J +m ) * factorial(J -m);
  double s1 = factorial(J1+m1) * factorial(J1-m1);
  double s2 = factorial(J2+m2) * factorial(J2-m2);
  s = sqrt(s * s1 * s2);

  //printf(" s0, s = %f , %f \n", s0, s);

  int kMax = min( min( J1+J2-J, J1 - m1), J2 + m2);

  double CG = 0.;
  for( int k = 0; k <= kMax; k++){
    double k1 = factorial(J1+J2-J-k);
    double k2 = factorial(J1-m1-k);
    double k3 = factorial(J2+m2-k);
    double k4 = factorial(J - J2 + m1 +k);
    double k5 = factorial(J - J1 - m2 +k);
    double temp = pow(-1, k) / (factorial(k) * k1 * k2 * k3 * k4 * k5);
    if( k1 == -100. || k2 == -100. || k3 == -100. || k4 == -100. || k5 == -100. ) temp = 0.;

    //printf(" %d | %f \n", k, temp);
    CG += temp;
  }

  return s0*s*CG;

}

double ThreeJSymbol(double J1, double m1, double J2, double m2, double J3, double m3){

  // ( J1 J2 J3 ) = (-1)^(J1-J2 - m3)/ sqrt(2*J3+1) * CGcoeff(J3, -m3, J1, m1, J2, m2) 
  // ( m1 m2 m3 )

  return pow(-1, J1 - J2 - m3)/sqrt(2*J3+1) * CGcoeff(J3, -m3, J1, m1, J2, m2);

}

double SixJSymbol(double J1, double J2, double J3, double J4, double J5, double J6){

  // coupling of j1, j2, j3 to J-J1
  // J1 = j1
  // J2 = j2
  // J3 = j12 = j1 + j2
  // J4 = j3
  // J5 = J = j1 + j2 + j3
  // J6 = j23 = j2 + j3

  //check triangle condition
  if( J3 < abs(J1 - J2 ) || J1 + J2 < J3 ) return 0; 
  if( J6 < abs(J2 - J4 ) || J2 + J4 < J6 ) return 0;
  if( J5 < abs(J1 - J6 ) || J1 + J6 < J5 ) return 0;
  if( J5 < abs(J3 - J4 ) || J3 + J4 < J5 ) return 0;

  double sixJ = 0;

  for( float m1 = -J1; m1 <= J1 ; m1 = m1 + 1){
    for( float m2 = -J2; m2 <= J2 ; m2 = m2 + 1){
      for( float m3 = -J3; m3 <= J3 ; m3 = m3 + 1){
        for( float m4 = -J4; m4 <= J4 ; m4 = m4 + 1){
          for( float m5 = -J5; m5 <= J5 ; m5 = m5 + 1){
            for( float m6 = -J6; m6 <= J6 ; m6 = m6 + 1){

              double f = (J1 - m1) + (J2 - m2) + (J3 - m3) + (J4 - m4) + (J5 - m5) + (J6 - m6);

              double a1 = ThreeJSymbol( J1, -m1, J2, -m2, J3, -m3); // J3 = j12 
              double a2 = ThreeJSymbol( J1, m1, J5, -m5, J6, m6); // J5 = j1 + (J6 = j23)
              double a3 = ThreeJSymbol( J4, m4, J2, m2, J6, -m6); // J6 = j23
              double a4 = ThreeJSymbol( J4, -m4, J5, m5, J3, m3); // J5 = j3 + j12

              double a = a1 * a2 * a3 * a4;
              //if( a != 0 ) printf("%4.1f %4.1f %4.1f %4.1f %4.1f %4.1f | %f \n", m1, m2, m3, m4, m5, m6, a);

              sixJ += pow(-1, f) * a1 * a2 * a3 * a4;

            }
          }
        }
      }
    }
  }

  return sixJ;
}

double NineJSymbol( double J1, double J2, double J3, double J4, double J5, double J6, double J7, double J8, double J9){

  double gMin = min( min (min( abs(J1 - J2 ), abs(J4 - J5)) , abs( J4 - J6 )) , abs(J7 - J8));
  double gMax = max( max (max( J1+J2, J4+J5), J3+J6), J7+J8);

  //printf(" gMin, gMax = %f %f \n", gMin, gMax);

  double nineJ = 0;
  for( float g = gMin; g <= gMax ; g = g + 1){
    double f = pow(-1, 2*g) * (2*g+1);
    double s1 = SixJSymbol(J1, J4, J7, J8, J9, g);
    if( s1 == 0 ) continue;
    double s2 = SixJSymbol(J2, J5, J8, J4, g, J6);
    if( s2 == 0 ) continue;
    double s3 = SixJSymbol(J3, J6, J9, g, J1, J2);
    if( s3 == 0 ) continue;
    nineJ += f* s1*s2*s3;
  }

  return nineJ;
}

Sum rules of Clebsch-Gordon Coefficient

Since the CG coefficient is already normalized.Thus

Since the number of is , as . Thus,

At last, the number of dimension of the coupled space or (tensor product space) is equation to , i.e.

Thus,

Winger 6-j and 9-j symbol

The meaning of 3-j symbol is same as Clebsch-Gordan coefficient. So, we skip in here.

I am not going to construct the 6-j symbol from 3-j symbol. In here, I just state the meaning and usage in Mathematica.

---

The 6-j symbol is the coupling between 3 angular momenta, .

There are 2 ways to couple these 3 angular momenta. First,

the other way is

The 6-j symbol is

We can see that there are 4 vector-sum must satisfy.

If we draw a line to connect these 4 vector-sum, we have:

In Mathematica, there is a build in function

---

The 9-j symbol is the coupling between 4 angular momenta, .

The 9-j symbol can be used in coupling 2 nucleons, .

The 9-j symbol is

We can see, each row and column must satisfy the vector-sum.

Unfortunately, there is no build in function in Mathematica. The formula for 9-j symbol is

Where sum all possible value, which can be calculate using the 6 couplings inside the 3 6-j symbols.To check your result, the coupling between and to from a state is

 

 

 

Hartree method for Helium ground state

After long preparation, I am ready to do this problem.

The two electron in the helium ground state occupy same spacial orbital but difference spin. Thus, the total wavefunction is

Since the Coulomb potential is spin-independent, the Hartree-Fock method reduce to Hartree method. The Hartree operator is

where the single-particle Hamiltonian and mutual interaction are

In the last step, we use atomic unit, such that . And the energy is in unit of Hartree, .

---

We are going to use Hydrogen-like orbital as a basis set.

I like the left the , because in the integration , the can be cancelled. Also, the is a compact index of the orbital.

Using basis set expansion, we need to calculate the matrix elements of

---

Now, we will concentrate on evaluate the mutual interaction integral.

Using the well-known expansion,

The angular integral

where the integral .

From this post, the triplet integral of spherical harmonic is easy to compute.

The Clebsch-Gordon coefficient imposed a restriction on .

The radial part,

The algebraic calculation of the integral is complicated, but after the restriction of from the Clebsch-Gordon coefficient, only a few terms need to be calculated.

---

The general consideration is done. now, we use the first 2 even states as a basis set.

These are both s-state orbital. Thus, the Clebsch-Gordon coefficient

The radial sum only has 1 term. And the mutual interaction becomes

The angular part

Thus, the mutual interaction energy is

The radial part

We can easy to see that . Thus, if we flatten the matrix of matrixes, it is Hermitian, or symmetric.

We also notice that is a matrix of matrixes. The inner matrix runs for indexes, while the bigger matrix runs for index. This will simplify our minds.

---

Now, we can start doing the Hartree method.

The general solution of the wave function is

The Hartree matrix is

In this simple 2-states example, index runs from 1 to 2. And remember the is a matrix of matrixes, using the identify the sum becomes:

,

where the represent the bigger matrix in .

The first trial wave function is the Hydrogen-like orbital,

Solve for eigen system, we have the energy after 1st trial,

After 13th trial,

Thus, the mixing of the 2s state is only 3.7%.

Since the eigen energy contains the 1-body energy and 2-body energy. So, the total energy for 2 electrons is

In which ,

So the energies for

From He to He++.  
From He+ to He++, .
From He to He+, is

The experimental 1 electron ionization energy for Helium atom is

The difference with experimental value is 2.175 eV. The following plot shows the Coulomb potential, the screening due to the existence of the other electron, the resultant mean field, the energy, and

---

Usually, the Hartree method will under estimate the energy, because it neglected the correlation, for example, pairing and spin dependence. In our calculation, the energy is under estimated.

From the , we can see, the mutual interaction between 1s and 2s state is attractive. While the interaction between 1s-1s and 2s-2s states are repulsive. The repulsive can be easily understood. But I am not sure how to explain the attractive between 1s-2s state.

Since the mass correction and the fine structure correction is in order of , so the missing 0.2 eV should be due to something else, for example, the incomplete basis set.

If the basis set only contain the 1s orbit, the mutual interaction is 1.25 Hartree = 34.014 eV. Thus, the mixing reduce the interaction by 5.07 eV, just for 3.7% mixing

I included the 3s state,

The mutual energy is further reduced to 1.05415 Hartree = 28.6848 eV. The . If 4s orbital included, the . We can expect, if more orbital in included, the will approach to .

A little summery of what I am doing

My current position is developing a computational program that can measure the system of polarized target automatically and repeatedly. The program needs to connect with the microwave generator, voltage supply, power meter, and oscilloscope. That is purely technical and nothing special.

Later, after the program started to collecting data, another program is needed to analysis the data. Although there are many program that can do analysis, but those program are not so easy to use, in both control and display. So, I make an analysis program. The program is simply read the 2-D time-domain signal, and then determine some parameters of a specific function that fit the signal. However, the signal could be very noisy, so FFTW and wavelet analysis were implemented. That’s why the wavelet analysis appeared in this blog.

After that, the sample of polarized target is made from various pentacene derivatives, that no known energy level exist. One way to find out the energy levels of the singlet excited states is to measure the absorption spectrum. However, the energy of the triplet state is difficult to measure. And the energy level of the triplet state is critical for Dynamic Nuclear Polarization to be happened. Thus, one way to find out is doing computation chemistry.

My last chemistry class was like 15 years ago, when I was junior high school. But the basic of computation chemistry is solving the Schrodinger equation. That is what theoretical nuclear physicists do! The variation method, the Hartree-fock method, I heard it and somehow know it, but never do it with my own hand and computer. That is why I revisit the Hartree-fock method, and found out my previous understanding is so naive.

In the course of studying Hartree-Fock method, and one problem is evaluating the overlap integral

For a 3-D system, the integral involves product of multiple spherical harmonics. That is really troublesome. Therefore, I move to study the spherical harmonics, and the related rotational invariant, Wigner D-matrix, Clebcsh-Gordon series and Fourier series.  The spherical harmonics arises from solving the Laplace equation in spherical coordinate. A general theory of the solution of Lapalace equation involves Legendre polynomial, which is a special case for Hypergeometric function.  And a very interesting connection is that the elliptical function of the 1st and 2nd kind are also two special cases of hypergeometric function, that, the solution of Laplace equation for elliptical boundary condition is elliptical functions. That connects spherical harmonics and elliptical function! Wow!

I am now a bit off-track, that I am very interesting on the function of all functions. We know that there are many elementary functions, such as sin, cos, and Log. And even more kind of special functions, such as

• Hermite — solving 1-D harmonic oscillator
• Laguerre — the radial function of hydrogen
• Legendre — the solution of the “” of  Laplace equation
• Gamma — a continuation of factorial
• Bessel — solution of 3-D infinite square well
• Elliptic — magnetic field of a solenoid
• Dirac delta
• Gaussian or Error function

For discrete argument

• Clebsch-Gordon
• Factoral
• Binomial

As far as I know, the Hypergemetric function is like “mother of functions”, although not all special functions can be expressed as it. I am driven by curiosity, so, I am not sure where I will go. For instance, the transformation of Hypergeometric function is very interesting.

So, for now, as an ending, I found one article is very interesting. The History and Future of Special Functions, by Stephen Wolfram.