• April 1st, 2025

Computing angular momentum operators

  • Author: Quentin Stern
  • Editor: Joshua Straub

Hi reader! This post is the first of hopefully a long series, where we will show how to simulate magnetic resonance experiments, starting from the basics (eg: the liquid-state 1D NMR spectrum of a simple spin system), and building up to more complicated topics like multidimensional experiments or quantum gates.

The aim of these posts is not to obtain the most performant computation. If that’s what you’re searching for, you’d rather use simulation packages that have been well optimized by clever people. Our aim here is to learn how things work by coding them. For this reason, we don’t want to use a lot of external functions that would hide away some parts of the code. In all the posts that I will write, I will make use of a single predefined function, to compute angular momentum operators. Why? Because even with only two spins, you already have a set of 16 operators, and as you add more spins, the list gets very long. Yet, there’s not much to learn in this lengthy and repetitive process of constructing the operators. In fact, this approach is exactly what I use for my own research. I’m not making up something different for the purpose of these posts.

This first post is about how angular momentum operators are generated in MATLAB and Python. We will start by showing how single-spin operators are generated from scratch for arbitrary values of the the spin number I = 1/2, 1, 3/2, .... Then, we will use the Kronecker product to compute operators in systems containing multiple spins. Once we’ve learnt how to code the operators from scratch, we will see how to compute them using the function that I wrote. Finally, we will show some examples of Hamiltonians (eg: for the five 1H spins of ethanol in the liquid state).

Note that we will omit \hbar in the expression of the operators because all Hamiltonians will be expressed in terms of rad.s-1.

Writing from scratch

Single-spin systems

For a single spin 1/2, the three Cartesian operators expressed in terms of Pauli matrices in the Zeeman basis read

(1)   \begin{equation*}\begin{split} \hat{I}_x &= \frac{1}{2}\begin{pmatrix}0 & 1\\ 1 & 0 \end{pmatrix} \\ \hat{I}_y &= \frac{i}{2}\begin{pmatrix}0 & -1\\ +1 & 0 \end{pmatrix} \\ \hat{I}_z &= \frac{1}{2}\begin{pmatrix}+1 & 0\\ 0 & -1 \end{pmatrix},\end{split}\end{equation*}

with the identity

(2)   \begin{equation*} \hat{1} = \begin{pmatrix}1 & 0\\ 0 & 1 \end{pmatrix} .\end{equation*}

For a single spin 1, they read

(3)   \begin{equation*}\begin{split} \hat{I}_x &= \frac{1}{\sqrt{2}}\begin{pmatrix}0 & 1 & 0\\ 1 & 0 & 1 \\ 0 & 1 & 0\end{pmatrix} \\ \hat{I}_y &= \frac{i}{\sqrt{2}}\begin{pmatrix}0 & -1 & 0\\ +1 & 0 & -1 \\ 0 & +1 & 0\end{pmatrix} \\ \hat{I}_z &=\begin{pmatrix}+1 & 0 & 0\\ 0 & 0 & 0 \\ 0 & 0 & -1\end{pmatrix},\end{split}\end{equation*}

with the identity

(4)   \begin{equation*} \hat{1} = \begin{pmatrix}1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\end{pmatrix} .\end{equation*}

The dimension of the matrix representing the angular momentum operator for a spin I is (2I+1)\times(2I+1), yielding 2\times 2 and 3\times 3 for spins 1/2 and 1, respectively. The expressions for cases with I>1 are given below.

For I > 1

For any spin number I, there are 2I+1 possible projections of the angular momentum m, with m\in[-I,-I+1,..., +I-1, +I]. Writing the Zeeman states as \ket{m}, the matrix elements of the angular momentum operators read

    \begin{equation*}\begin{split} \braket{m'\vert\hat{I}_x\vert m}&= \left(\delta_{m'+1,m}+\delta_{m',m+1}\right)\frac{1}{2}\sqrt{I(I+1)-m'm} \\ \braket{m'\vert\hat{I}_y\vert m} &= \left(\delta_{m'+1,m}-\delta_{m',m+1}\right)\frac{i}{2}\sqrt{I(I+1)-m'm} \\ \braket{m'\vert\hat{I}_z\vert m} &= \delta_{m',m}I(I+1),\end{split}\end{equation*}

where \delta is the Kronecker delta function

    \begin{equation*} \delta_{nm} = \biggl\{ \begin{matrix} 1 & \textrm{if} & n=m \\ 0 & \textrm{if} & n\ne m .\end{matrix}\end{equation*}

The elements of the corresponding identity matrix are given by

    \begin{equation*} \braket{m'\vert\hat{1}\vert m}=\delta_{m',m} .\end{equation*}

It is useful to define shift operators as

(5)   \begin{equation*} \hat{I}_\pm = \hat{I}_x\pm i\hat{I}_y.\end{equation*}

In the case of I=1/2, the shift operators read

(6)   \begin{equation*}\begin{split} \hat{I}_+ &= \begin{pmatrix}0 & 1\\ 0 & 0 \end{pmatrix} \\ \hat{I}_- &= \begin{pmatrix}0 & 0\\ 1 & 0 \end{pmatrix}.\end{split}\end{equation*}

The code snippet below shows how some of these operators are generated programatically.

  • MATLAB
  • Python
% Definition of the angular momentum operators for I = 1/2
Ix2x2   = [0  +1;+1  0]/2;
Iy2x2   = [0  -1;+1  0]*1i/2;
Iz2x2   = [+1  0; 0 -1]/2;
E2x2    = eye(2);
Ip2x2   = Ix2x2 + 1i*Iy2x2;
Im2x2   = Ix2x2 - 1i*Iy2x2;

% Definition of the angular momentum operators for I = 1
Ix3x3   = [0  +1 0;+1  0 +1; 0  +1  0]/sqrt(2);
Iy3x3   = [0  -1 0;-1  0 +1; 0  +1  0]*1i/sqrt(2);
Iz3x3   = [+1  0 0; 0  0  0; 0  0  -1];
E3x3    = eye(3);
Ip3x3   = Ix3x3 + 1i*Iy3x3;
Im3x3   = Ix3x3 - 1i*Iy3x3;
# Required package for numerical calculations
import numpy as np

# Definition of the angular momentum operators for I = 1/2
Ix2x2 = np.array([[0,  +1/2],[ +1/2, 0]])
Iy2x2 = np.array([[0, -1j/2],[+1j/2, 0]])
Iz2x2 = np.array([[+1/2,  0],[0,  -1/2]])
E2x2  = np.eye(2)
Ip2x2 = Ix2x2 + 1j*Iy2x2
Im2x2 = Ix2x2 - 1j*Iy2x2

# Definition of the angular momentum operators for I = 1
Ix3x3 = np.array([[ 0, +1, 0],[+1, 0, +1],[0, +1,  0]])/np.sqrt(2) 
Iy3x3 = np.array([[ 0, -1, 0],[-1, 0, +1],[0, +1,  0]])*1j/np.sqrt(2)
Iz3x3 = np.array([[+1,  0, 0],[ 0, 0,  0],[0,  0, -1]])
E3x3  = np.eye(3) 
Ip3x3 = Ix3x3 + 1j*Iy3x3
Im3x3 = Ix3x3 - 1j*Iy3x3

Multiple-spin systems

To represent multiple spins in the same system, we need to take Kronecker products of the individual spin operators. If the system contains N spins with spin numbers I_k, the operator \hat{I}_{k\mu} of spin k with \mu = x,y,z,+,- represented in the space of the N spins is given by

(7)   \begin{equation*}\begin{split} \hat{I}_{k\mu} &= \otimes_{j=1}^N \hat{u}_\mu^{(j)}\\&=\hat{1}^{(1)}\otimes\hat{1}^{(2)}\otimes ... \otimes \hat{I}_\mu^{(k)} \otimes... \otimes\hat{1}^{(N)},\end{split}\end{equation*}

where the operators \hat{u}_\mu^{(j)} are defined in the subspace of the j spin (i.e., one of the operators that we have defined in the previous section) as

(8)   \begin{equation*}\hat{u}_\mu^{(j)}=\biggl\{ \begin{matrix} \hat{1}^{(j)} & \textrm{if} & j \ne k \\ \hat{I}_\mu^{(j)} & \textrm{if} & j = k, \end{matrix}\end{equation*}

and have dimensions of (2I_j+1)\times(2I_j+1). What Eq. 7 tells us is that to obtain an angular momentum operator of spin k in the larger space of the N spins, we need to take the Kronecker product over all spins, using the identity operator for each spin in its own subspace, except for k, where we use the operator of interest in its own subspace (eg, \hat{I}_x^{(k)} or \hat{I}_+^{(k)}).

Kronecker product

The Kronecker product is an operation that takes two matrices with dimensions n \times m and n' \times m' and returns a matrix with dimension nn' \times mm'. For example, for two matrices of dimension 2\times 2, the Kronecker product is

    \begin{equation*}\begin{split}&\begin{pmatrix}a & b\\ c & d\end{pmatrix}\otimes\begin{pmatrix}\alpha & \beta\\ \gamma & \delta\end{pmatrix}\\=&\begin{pmatrix}a\begin{pmatrix}\alpha & \beta\\ \gamma & \delta\end{pmatrix} & b\begin{pmatrix}\alpha & \beta\\ \gamma & \delta\end{pmatrix}\\ c\begin{pmatrix}\alpha & \beta\\ \gamma & \delta\end{pmatrix} & d\begin{pmatrix}\alpha & \beta\\ \gamma & \delta\end{pmatrix}\end{pmatrix}\\=&\begin{pmatrix}a\alpha & a\beta & b\alpha & b\beta\\ a\gamma & a\delta & b\gamma & b\delta \\ c\alpha & c\beta & d\alpha & d\beta\\ c\gamma & c\delta & d\gamma & d\delta\end{pmatrix}\end{split}\end{equation*}

For example, if we want to compute the z angular momentum for a pair of spins 1/2, we have

(9)   \begin{equation*}\begin{split}\hat{I}_{1z} &= \hat{I}_z^{(1)} \otimes \hat{1}^{(2)}\\&= \left( \frac{1}{2}\begin{pmatrix}+1 & 0\\ 0 & -1 \end{pmatrix} \right)\otimes \begin{pmatrix}1 & 0\\ 0 & 1 \end{pmatrix} \\&= \frac{1}{2}\begin{pmatrix}+1 & 0 & 0 & 0\\ 0 & +1 & 0 & 0 \\ 0 & 0 & -1 & 0\\ 0 & 0 & 0 & -1\end{pmatrix},\end{split}\end{equation*}

for the first spin, and

(10)   \begin{equation*}\begin{split}\hat{I}_{2z} &= \hat{1}^{(1)} \otimes \hat{I}_z^{(2)}\\&= \begin{pmatrix}1 & 0\\ 0 & 1 \end{pmatrix} \otimes\left( \frac{1}{2}\begin{pmatrix}+1 & 0\\ 0 & -1 \end{pmatrix} \right)\\&= \frac{1}{2}\begin{pmatrix}+1 & 0 & 0 & 0\\ 0 & -1 & 0 & 0 \\ 0 & 0 & +1 & 0\\ 0 & 0 & 0 & -1\end{pmatrix},\end{split}\end{equation*}

for the second. We only show this simple case because, as you have probably noticed, the dimensions rapidly increase when we add more spins!

It is sometimes useful to define operators that correspond to the sum of the operators of all spins for a certain \mu

(11)   \begin{equation*}\hat{I}_\mu = \sum_{j=1}^N \hat{I}_{j\mu}.\end{equation*}

For example, the sum of the z angular momentum operators in a system consisting of a pair of spins 1/2 is

(12)   \begin{equation*}\begin{split}\hat{I}_z &= \hat{I}_{1z}+\hat{I}_{2z}\\&= \begin{pmatrix}+1 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0\\ 0 & 0 & 0 & -1\end{pmatrix}.\end{split}\end{equation*}

The code snippet below shows how the operators for a system of two spins 1/2 and one spin 1 are constructed, using the single-spin operators defined in the previous code snippet.

  • MATLAB
  • Python
% Definition of the angular momentum operators 
% of individual spins in a system of 2 spins 1/2 
% and 1 spins 1 
I1x     = kron(kron(Ix2x2,E2x2),E3x3);
I1y     = kron(kron(Iy2x2,E2x2),E3x3);
I1z     = kron(kron(Iz2x2,E2x2),E3x3);

I2x     = kron(kron(E2x2,Ix2x2),E3x3);
I2y     = kron(kron(E2x2,Iy2x2),E3x3);
I2z     = kron(kron(E2x2,Iz2x2),E3x3);

I3x     = kron(kron(E2x2,E2x2),Ix3x3);
I3y     = kron(kron(E2x2,E2x2),Iy3x3);
I3z     = kron(kron(E2x2,E2x2),Iz3x3);

% Sum of the angular mommentum operators
Ix      = I1x + I2x + I3x;
Iy      = I1y + I2y + I3y;
Iz      = I1z + I2z + I3z;
# Definition of the angular momentum operators 
# of individual spins in a system of 2 spins 1/2
# and 1 spin 1
I1x = np.kron(np.kron(Ix2x2, E2x2),E2x2)
I1y = np.kron(np.kron(Iy2x2, E2x2),E2x2)
I1z = np.kron(np.kron(Iz2x2, E2x2),E2x2)

I2x = np.kron(np.kron(E2x2, Ix2x2),E2x2)
I2y = np.kron(np.kron(E2x2, Iy2x2),E2x2)
I2z = np.kron(np.kron(E2x2, Iz2x2),E2x2)

I3x = np.kron(np.kron(E2x2, E2x2),Ix2x2)
I3y = np.kron(np.kron(E2x2, E2x2),Iy2x2)
I3z = np.kron(np.kron(E2x2, E2x2),Iz2x2)

# Sum of the angular mommentum operators
Ix  = I1x + I2x + I3x
Iy  = I1y + I2y + I3y
Iz  = I1z + I2z + I3z

Using a compact function

Single-spin operators and their sum

Now we move on to the function I that I promised in the introduction of this post. This function generates angular momentum operators for any spin system using a minimum of two arguments. The first argument is always a list containing the spin number of each spin in the system. If only a second argument is given, then it must be a string containing one of the \mu=x,y,z,+,-, and the output is the sum of all operators corresponding to \mu, as defined in Eq. 11. If three inputs are given, the second is the index or the list of indices corresponding to the spins for which the operator is to be computed. The third argument is a string containing one of the \mu=x,y,z,+,-. If a list is given as second input, then the output is the angular momentum operator for \mu summed over all spins in the list. The function can be used with more than three inputs to generate product operators. This is treated in the next section.

The code snippet below shows examples.

  • MATLAB
  • Python
% Definition of the spin system
SS  = [1/2 1/2 1];

% Definition of the angular momentum operators 
% of individual spins in a system of 3 spins 1/2
I1x = I(SS,1,'x');
I1y = I(SS,1,'y');
I1z = I(SS,1,'z');
I1p = I(SS,1,'+');
I1m = I(SS,1,'-');

I2x = I(SS,2,'x');
I2y = I(SS,2,'y');
I2z = I(SS,2,'z');
I2p = I(SS,2,'+');
I2m = I(SS,2,'-');

I3x = I(SS,3,'x');
I3y = I(SS,3,'y');
I3z = I(SS,3,'z');
I3p = I(SS,3,'+');
I3m = I(SS,3,'-');
# Definition of the spin system
SS = [1/2, 1/2, 1]

# Definition of the angular momentum operators 
# of individual spins in a system of 3 spins 1/2
I1x = I(SS, 0, 'x')
I1y = I(SS, 0, 'y')
I1z = I(SS, 0, 'z')  
I1p = I(SS, 0, '+')
I1m = I(SS, 0, '-')  

I2x = I(SS, 1, 'x')
I2y = I(SS, 1, 'y')
I2z = I(SS, 1, 'z')  
I2p = I(SS, 1, '+')   
I2m = I(SS, 1, '-')   

I3x = I(SS, 2, 'x')
I3y = I(SS, 2, 'y')
I3z = I(SS, 2, 'z')  
I3p = I(SS, 2, '+') 
I3m = I(SS, 2, '-')

The sum of the angular momentum operators can be computed in three different ways. First, using the function with only two arguments; second, using three arguments where the second argument is the list of all spin indices; last, writing the sum explicitly. The code snippet below shows these three cases for a system of three spins.

  • MATLAB
  • Python
Ix  = I(SS, 'x');
Iy  = I(SS, 'y');
Iz  = I(SS, 'z');
Ip  = I(SS, '+');
Im  = I(SS, '-');

Ix  = I(SS, 1:3, 'x');
Iy  = I(SS, 1:3, 'y');
Iz  = I(SS, 1:3, 'z');
Ip  = I(SS, 1:3, '+');
Im  = I(SS, 1:3, '-');

Ix  = I(SS,1,'x') + I(SS,2,'x') + I(SS,3,'x');
Iy  = I(SS,1,'y') + I(SS,2,'y') + I(SS,3,'y');
Iz  = I(SS,1,'z') + I(SS,2,'z') + I(SS,3,'z');
Ip  = I(SS,1,'+') + I(SS,2,'+') + I(SS,3,'+');
Im  = I(SS,1,'-') + I(SS,2,'-') + I(SS,3,'-');
Ix  = I(SS, 'x')
Iy  = I(SS, 'y')
Iz  = I(SS, 'z')
Ip  = I(SS, '+')
Im  = I(SS, '-')

Ix  = I(SS, [0, 1, 2], 'x')
Iy  = I(SS, [0, 1, 2], 'y')
Iz  = I(SS, [0, 1, 2], 'z')
Ip  = I(SS, [0, 1, 2], '+')
Im  = I(SS, [0, 1, 2], '-')

Ix  = I(SS,0,'x') + I(SS,1,'x') + I(SS,2,'x')
Iy  = I(SS,0,'y') + I(SS,1,'y') + I(SS,2,'y')
Iz  = I(SS,0,'z') + I(SS,1,'z') + I(SS,2,'z')
Ip  = I(SS,0,'+') + I(SS,1,'+') + I(SS,2,'+')
Im  = I(SS,0,'-') + I(SS,1,'-') + I(SS,2,'-')

Mulitple-spin operators

Another thing that the function I can do for you is to compute mutliple-spin operators, also called product operators. These operators are nothing more than the product of two or more of the operators defined above. However, in the interest of making the scripts as compact and readable as possible, I found it useful to make the function directly output products of operators. In this case, more than three arguments are used. The first argument is as always the list of spin numbers. The next arguments come in pairs. For each pair, the first and second of each pair are the index of the spin and the value of \mu defined as a string, respectively. The resulting operator is the product of the operators corresponding to each pair, Below is a list of possible uses of this option for a 2-spin and a 3-spin operators in the same spin system as above. In both cases, two ways to compute the same operators are shown.

  • MATLAB
  • Python
% Example of 2-spin operator
I1x2z   = I(SS,1,'x')*I(SS,2,'z');
I1x2z   = I(SS,1,'x',2,'z');

% Example of 3-spin operator
I1z2z3z = I(SS,1,'z')*I(SS,2,'z')*I(SS,3,'z');
I1z2z3z = I(SS,1,'z',2,'z',3,'z');
# Example of 2-spin operator
I1x2z   = np.matmul(I(SS,0,'x'),I(SS,1,'z'))
I1x2z   = I(SS,0,'x',1,'z')

# Example of 3-spin operator
I1z2z3z = np.matmul(np.matmul(I(SS,0,'z'),I(SS,1,'z')),I(SS,2,'z'))
I1z2z3z = I(SS,0,'z',1,'z',2,'z')

The last option of the function is to compute dot products. The dot product between two spins k and j can be written in several forms

(13)   \begin{equation*}\begin{split}\boldsymbol{\hat{I}}_k\cdot \boldsymbol{\hat{I}}_j & = \begin{pmatrix}\hat{I}_{kx} & \hat{I}_{ky} & \hat{I}_{kz} \end{pmatrix}\begin{pmatrix}\hat{I}_{jx} \\ \hat{I}_{jy} \\ \hat{I}_{jz} \end{pmatrix} \\&=\hat{I}_{kx}\hat{I}_{jx}+\hat{I}_{ky}\hat{I}_{jy}+\hat{I}_{kz}\hat{I}_{jz}. \end{split}\end{equation*}

Two of these forms can be computed using the I function, the most compact being the one corresponding to the left hand side of Eq. 13. In this case, the second, third, and fourth arguments are the index of the first spin, the string ‘dot’, and the index of the second spin, respectively. The function also supports lists of indices in the second and fourth arguments, which can be useful when dealing with equivalent spins, as in the example of the next section.

The code snippet below shows how to compute the scalar product using different approaches, the last one being the most compact. It is applied to two cases: the first is the dot product of spin 1 and spin 2; the second is the dot product of spins 1 and 2 with spin 3, in other words, the sum of the dot products of spin 1 and spin 2 and of spin 2 and spin 3.

  • MATLAB
  • Python
% Scalar product operator ot spin 1 with spin 2
I1dot2  = I(SS,1,'x')*I(SS,2,'x') + I(SS,1,'y')*I(SS,2,'y') + I(SS,1,'z')*I(SS,2,'z');
I1dot2  = I(SS,1,'dot',2);

% Scalar product operator ot spin 1 and 2 with spin 3 
I12dot3 = I(SS,1,'x')*I(SS,3,'x') + I(SS,1,'y')*I(SS,3,'y') + I(SS,1,'z')*I(SS,3,'z')+ ...
          I(SS,2,'x')*I(SS,3,'x') + I(SS,2,'y')*I(SS,3,'y') + I(SS,2,'z')*I(SS,3,'z');
I12dot3 = I(SS,1,'dot',3) + I(SS,2,'dot',3);
I12dot3 = I(SS,[1 2],'dot',3);
# Scalar product operator of spin 1 with spin 2 
I1dot2  = np.matmul(I(SS,0,'x'),I(SS,1,'x')) + np.matmul(I(SS,0,'y'),I(SS,1,'y')) \
          + np.matmul(I(SS,0,'z'),I(SS,1,'z'))
I1dot2  = I(SS,0,'x',1,'x') + I(SS,0,'y',1,'y') + I(SS,0,'z',1,'z')
I1dot2  = I(SS,0,'dot',1)

# Scalar product operator ot spin 1 and 2 with spin 3 
I12dot3  = np.matmul(I(SS,0,'x'),I(SS,2,'x')) + np.matmul(I(SS,0,'y'),I(SS,2,'y')) \
          + np.matmul(I(SS,0,'z'),I(SS,2,'z')) \
          + np.matmul(I(SS,1,'x'),I(SS,2,'x')) + np.matmul(I(SS,1,'y'),I(SS,2,'y')) \
          + np.matmul(I(SS,1,'z'),I(SS,2,'z'))
I12dot3 = I(SS,0,'x',2,'x') + I(SS,0,'y',2,'y') + I(SS,0,'z',2,'z') \
          + I(SS,1,'x',2,'x') + I(SS,1,'y',2,'y') + I(SS,1,'z',2,'z')
I12dot3 = I(SS,0,'dot',2) + I(SS,1,'dot',2)
I12dot3 = I(SS,[0,1],'dot',2)

Hamiltonian examples

The function I comes in handy when dealing with large numbers of spins and their couplings. Take for example the five 1H spins in the ethyl chain of the molecule of ethanol in the liquid state. The two methylene 1H spins (labelled 1 and 2) are J coupled with the three methyl 1H spins (labelled 3, 4, and 5). The J-Hamiltonian therefore reads

(14)   \begin{equation*}\begin{split}\hat{H}_J =& 2\pi J\left(\left( \boldsymbol{\hat{I}}_1+\boldsymbol{\hat{I}}_2 \right) \cdot \left( \boldsymbol{\hat{I}}_3+\boldsymbol{\hat{I}}_4+\boldsymbol{\hat{I}}_5 \right) \right)\\=&2\pi J(\boldsymbol{\hat{I}}_1\cdot\boldsymbol{\hat{I}}_3+\boldsymbol{\hat{I}}_1\cdot\boldsymbol{\hat{I}}_4+\boldsymbol{\hat{I}}_1\cdot\boldsymbol{\hat{I}}_5\\&+\boldsymbol{\hat{I}}_2\cdot\boldsymbol{\hat{I}}_3+\boldsymbol{\hat{I}}_2\cdot\boldsymbol{\hat{I}}_4+\boldsymbol{\hat{I}}_2\cdot\boldsymbol{\hat{I}}_5),\end{split}\end{equation*}

where the dot products can be expanded in Cartesian coordinates, using Eq. 13. This Hamiltonian would be quite long to write if we were to compute each operator manually and assemble them to make all the dot products. Using the I function, it takes the form shown in the snippet below.

  • MATLAB
  • Python
% Spin sytem - 5 spins 1/2
SS  = [1/2 1/2 1/2 1/2 1/2];

% J coupling, Hz
J   = 10;

% J-coupling Hamitlonian, rad.s-1
HJ  = 2*pi*I(SS,[1 2],'dot',[3 4 5]);
# Spin sytem - 5 spins 1/2
SS  = [1/2, 1/2, 1/2, 1/2, 1/2]

# J coupling, Hz
J   = 10

# J-coupling Hamitlonian, rad.s-1
HJ  = 2*np.pi*J*I(SS,[0, 1],'dot',[2, 3, 4])

Another example is the dipolar coupling Hamiltonian, which can be written in different forms

(15)   \begin{equation*}\begin{split}\hat{H}_\textrm{D} &= D_0\left(3\hat{I}_{1z}\hat{I}_{2z}- \boldsymbol{\hat{I}}_1\cdot\boldsymbol{\hat{I}}_1 \right)\\&= D_0\left(2\hat{I}_{1z}\hat{I}_{2z}- \frac{1}{2}(\hat{I}_{1+}\hat{I}_{2-}+\hat{I}_{1-}\hat{I}_{2+}) \right)\\&= D_0\left(2\hat{I}_{1z}\hat{I}_{2z}- (\hat{I}_{1x}\hat{I}_{2x}+\hat{I}_{1y}\hat{I}_{2y}) \right).\end{split}\end{equation*}

The code snippet below shows all three forms.

  • MATLAB
  • Python
% Spin sytem - 2 spins 1/2
SS  = [1/2, 1/2];

% Dipolar coupling, rad.s-1
D0  = 2*pi*10e3;

% Dipolar Hamitlonian, rad.s-1
HD  = D0*(3*I(SS,1,'z',2,'z') - I(SS,1,'dot',2));
HD  = D0*(2*I(SS,1,'z',2,'z') - 1/2*(I(SS,1,'+',2,'-') + I(SS,1,'-',2,'+')));
HD  = D0*(2*I(SS,1,'z',2,'z') - (I(SS,1,'x',2,'x') + I(SS,1,'y',2,'y')));
# Homonuclear dipolar coupling Hamiltonian
# Spin sytem - 2 spins 1/2
SS  = [1/2, 1/2]

# Dipolar coupling, rad.s-1
D0  = 2*np.pi*10e3

# Dipolar Hamitlonian, rad.s-1
HD  = D0*(3*I(SS,0,'z',1,'z') - I(SS,0,'dot',1))
HD  = D0*(2*I(SS,0,'z',1,'z') - 1/2*(I(SS,0,'+',1,'-') + I(SS,0,'-',1,'+')))
HD  = D0*(2*I(SS,0,'z',1,'z') - (I(SS,0,'x',1,'x') + I(SS,0,'y',1,'y')))

Conclusion

This post was perhaps not the most exciting. Very technical and not a lot of “concrete” results. But it gives a good basis to move on to spin physics simulations. The next post will be on simulating the NMR spectrum of simple spin systems. Stay tuned!

Leave a Reply
    Artist Credit:
    Yu SciVis & Art LLC (Dr. Chung-Jui Yu)
    Contact:
    raymondyu93@gmail.com
    Website designed and developed by:
    NetzOptimize Inc.
    © COPYRIGHT . QUANTUM-RESONANCE.ORG

    Quantum Insights