Euler method is the simplest way to solve differential equation numerically.
$$ y(t_{n+1}) = y(t_n) + \left.\frac{dy}{dt}\right|{t=t_n} (t{n+1}-t_{n}) $$
By adapting Euler method to Schrodinger equation,
$$ i\hbar \dot{\vec{c}} = \hat{H}\vec{c} \\
\vec{c}(t_{n+1}) = \vec{c}(t_{n}) + \left. \frac{\hat{H}}{i\hbar} \right|{t=t_n} (t{n+1}-t_{n}) $$
As a simplest example, let’s solve the dynamics of the two-level system.
The RWA Hamiltonian is
$$ \hat{H}_{\rm RWA} = \frac{\hbar}{2}\begin{bmatrix} 0 && \Omega \\ \Omega^* && -2\Delta\end{bmatrix} $$
and initial condition is
$$ \vec{c}(t=0) = \begin{bmatrix} c_0(t=0) \\ c_1(t=0)\end{bmatrix} = \begin{bmatrix} 1 \\ 0\end{bmatrix} $$
For simplicity, let
$$ \Omega = (2\pi) 1~{\rm MHZ} \\ \Delta = (2\pi) 0~{\rm MHZ} $$
import numpy as np
import matplotlib.pyplot as plt
# hbar
#hbar = 1.05e-34
hbar = 1/(2*np.pi) # natural unit: h=1
# time
t_ini = 0 # initial time
t_fin = 5e-6 # final time
t_size = 50000 # number of samples
t_array = np.linspace(t_ini, t_fin, t_size) # time array
dt = t_array[1] - t_array[0] # infinitesimal time
# Pauli matrices
sigma_x = np.array([[0, 1], [1, 0]]) # X matrix
sigma_n = np.array([[0, 0], [0, 1]]) # Z matrix
# Parameters
Omega = (2*np.pi)*1e+6 # 1 MHz
Delta = (2*np.pi)*0e+6 # 0 MHz
# Initial state
c = np.zeros((2, t_size), dtype=np.complex128) # c-vector
c[0, 0] = 1 # Initial condition
# Hamiltonian (time-independent)
H = (hbar/2) * (Omega * sigma_x + (-2 * Delta) * sigma_n)
# Time evolution using vectorized computation
for idx_time in range(t_size - 1):
c[:, idx_time + 1] = c[:, idx_time] + (H / (1j * hbar)) * dt @ c[:, idx_time]
# Populations
pop_0 = np.abs(c[0, :])**2 # Population of state |0>
pop_1 = np.abs(c[1, :])**2 # Population of state |1>
# Plotting results
plt.plot(t_array/1e-6, pop_0, label='|c_0|^2')
plt.plot(t_array/1e-6, pop_1, label='|c_1|^2')
plt.ylim([0, 1])
plt.xlim([t_ini/1e-6, t_fin/1e-6])
plt.xlabel('Time (us)')
plt.ylabel('Population')
plt.legend()
plt.show()
This is Rabi oscillation.
Let
$$ \Omega = (2\pi) 1~{\rm MHZ} \\ \Delta = (2\pi) 0.5~{\rm MHZ} $$
By just changing the variable ‘Delta’,