我想看看一些非常简单的自旋系统的密度矩阵的时间演化，但我的方法遇到了麻烦。

我想使用一个简单的for-loop 时间步进方法，而不是像ode45应该是演示代码那样的 ODE 求解器。

我在 Matlab 中工作（数值计算），所以我会在相关的地方包含代码片段。

我做的事

我为我的系统构建了仅自旋的哈密顿量，分为两部分。第一部分描述了自旋系统与外部磁场的塞曼相互作用。$S_A$和$S_B$是两个电子。$S_C$是一个核：

$$H_{magnetic} = \omega_0^e(S_{Az} + S_{Bz}) + \omega_{0}^nS_{Cz}$$ 在哪里$\omega_0^e = B_0*\gamma_e$是电子拉莫尔频率，并且$\omega_0^n$核拉莫尔频率。（可以说，核塞曼相互作用可以忽略不计，但为了完整起见，我将其包括在内）。

哈密​​顿量的第二部分描述了一个电子与原子核通过超精细相互作用的超精细相互作用：

$$H_{hyperfine} = a\mathbf{S_A}\mathbf{S_C} = a*(S_{Ax}S_{Cx} + S_{Ay}S_{Cy} + S_{Az}A_{Cz})$$ 在哪里$a$是超精细耦合常数。我的纯自旋哈密顿是： $$H = H_{magnetic} + H_{hyperfine}$$

在我的 Matlab 代码中，它的内容如下：

% 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;

---

现在，因为没有电子-电子自旋相互作用，所以我定义了从 xyz 基到电子 Singlet-Tiplet (ST) 基的相似变换，使得：

$$H_{ST} = M^{-1}HM$$ 在 Matlab 中是：

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

---

然后我以纯单重态开始我的自旋密度，并使用 Liouville von Neumann 方程的积分形式

$$\hat\rho(t) = e^{-I\hat{H}t}\hat\rho(0)e^{I\hat{H}t}$$ 迭代时间点如下

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);

问题

好吧。它不起作用。进一步来说; Sd由于 S，我预计单线态密度会在 1 到 0.25 之间波动$\leftrightarrow$T 相互转换。我得到的结果确实是振荡的，但在 +1 和 -1 之间。此外，密度矩阵的迹并不总是 1。

问题

• 你能直接解决上面提到的任何问题吗？
• 您能否参考我的一篇文献（或者，甚至更好的代码），它使用 LvN 方程逐步传播密度矩阵？
• 如果这不是进行密度矩阵传播的可能方式，您能向我解释为什么吗？（请记住，我意识到这不是进行自旋密度模拟的最佳方式。我特别想以这种方式进行模拟）