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);
end
plot(T,Sd);
Problem
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.
Questions
- 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)
