I want to look at time evolution of the density matrices of some, very simple, spin systems, but I am having trouble with my approach.

I want to use a simple for-loop time-stepping method, not as ODE solver like ode45 as this is supposed to be a demonstration code.

I work in Matlab (numerical calculations), so I will include code snippets where relevant.

What I do

I construct the spin-only hamiltonian for my system in two parts. First part describes the Zeeman interactions of the spin system with an external magnetic field. $S_A$ and $S_B$ are two electrons. $S_C$ is a nucleus: $$ H_{magnetic} = \omega_0^e(S_{Az} + S_{Bz}) + \omega_{0}^nS_{Cz} $$ where $\omega_0^e = B_0*\gamma_e$ is the electron Larmor frequency, and $\omega_0^n$ the Nuclear Larmor frequency. (Arguably the nuclear Zeeman interaction could be ignored, but I've included it for the completeness' sake).

The second part of the Hamiltonian describes the hyperfine interaction of one of the electrons with the nucleus via the hyperfine interaction: $$ H_{hyperfine} = a\mathbf{S_A}\mathbf{S_C} = a*(S_{Ax}S_{Cx} + S_{Ay}S_{Cy} + S_{Az}A_{Cz}) $$ where $a$ is the hyperfine coupling constant. My spin-only hamiltonian is then: $$ H = H_{magnetic} + H_{hyperfine} $$

In my Matlab code this reads as:

% Pauli spin matrices
Ix = 0.5 * [0  1;  1   0];
Iy = 0.5 * [0 -1i; 1i  0];
Iz = 0.5 * [1  0;  0  -1];
Ie = eye(2);

% Successive kronecker product of three elements
tkron = @(x,y,z) kron(kron(x,y),z);

% Spin matrices
SAx = tkron(Ix,Ie,Ie); SAy = tkron(Iy,Ie,Ie); SAz = tkron(Iz,Ie,Ie);
SBx = tkron(Ie,Ix,Ie); SBy = tkron(Ie,Iy,Ie); SBz = tkron(Ie,Iz,Ie);
SCx = tkron(Ie,Ie,Ix); SCy = tkron(Ie,Ie,Iy); SCz = tkron(Ie,Ie,Iz);

% Input constants
gamma_n = 42.577e3; % Hz/mT
gamma_e = 1.76e8; % Hz/mT
hfc = 1.0; % Hyperfine coupling in mT
B0 = 0.01; % Magnetic field in mT

% Derived constants
a = gamma_e * hfc;
omega_0 = gamma_e*B0; % Electron Larmor frequency
omega0_n = -gamma_n*B0; % Nuclear Larmor frequency

% Hamiltonian
H_mag = omega_0.*(SAz + SBz) + omega0_n.*SCz;
H_hyperfine = a*(SAx*SCx + SAy*SCy + SAz*SCz);
H = H_mag + H_hyperfine;

Now, because there is no election-electron spin interaction I then define a similarity transform from xyz basis to electronic Singlet-Tiplet (ST) basis such that: $$ H_{ST} = M^{-1}HM $$ which in Matlab is:

s2 = 1/sqrt(2);
%       S  T0  T+  T-
Me = [ 0  0   1   0
        s2 s2  0   0
       -s2 s2  0   0
        0  0   0   1];
Mn = eye(2); % Nucleus stays in alpha/beta basis
M = kron(Me,Mn);

% Change to ST basis
Hst = M'*H*M; % M is unitary therefore use transpose

I then start my spin density in a pure singlet state and use the integrated form one the Liouville von Neumann equation $$ \hat\rho(t) = e^{-I\hat{H}t}\hat\rho(0)e^{I\hat{H}t} $$ to iterate through time points as follows

rh0 = zeros(8);
rh0(1,1) = 0.5;
rh0(2,2) = 0.5;
T = linspace(0,1e-6,1000); % 1 us of spin evolution
for i = 1:length(T);
    t = T(i);
    e = 1i*H*t;
    em = expm(-e);
    ep = expm(e)';
    rh_t = em*rh0*ep; % Propagate rh0 to time t
    Sd(i) = rh_t(1,1) + rh_t(2,2);


Well. It doesn't work. More specifically; I'd expect the singlet density Sd to oscillate between 1 and 0.25 due to S$\leftrightarrow$T interconversion. The results I get are indeed oscillatory, but between +1 and -1. Furthermore the trace of the density matrix is not always 1.


  • Can you address any of the problems mentioned above directly?
  • Can you reference me to a piece of literature (or, even better, code) which does this step-by-step propagation of the density matrix with the LvN equation?
  • If this is not a possible way of doing density matrix propagation, can you explain to me why? (Please remember that I realize this is not the optimal way of doing spin density simulations. I specifically want to do a simulation this very way)
  • 1
    $\begingroup$ Would Computational Science be a better home for this question? $\endgroup$
    – Qmechanic
    Oct 9, 2017 at 17:51
  • $\begingroup$ I didn't know about this stack exchange. I think it would be a better home indeed. Can I ask a moderator to move it please?? $\endgroup$ Oct 9, 2017 at 17:56

1 Answer 1


Ok, I figured it out. There were a couple of problems:

  • My similarity transform was bad (this was the major issue)
  • I should not have transposed the ep propagator
  • My nuclear gyromagnetic ratio was wrong (although that ended up not really mattering due to the nuclear Zeeman interaction being small in comparison to other parts of the Hamiltonian)

I have fixed the problem by doing the (correct) similarity transform first, and then working in ST basis throughout. What helped was to define the singlet and triplet projection operators (question 3d from this assignment) because even if I end up transforming to some cookie basis, as long as I'm consistent throughout (i.e. also transforming the projection operator), I end up with the correct result.

Anyway; HERE is my code, and here is why I wanted to understand that subject better. I hope this helps someone at some point.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.