In genral, Rabi frequency and detuning are not constants,
$$ \Omega = \Omega(t) \\ \Delta = \Delta(t) $$
Let’s make 2-pulse input
$$ \Omega_0 = (2\pi) 1~{\rm MHz} \\ \Delta_0 = (2\pi) 0.1~{\rm MHz} \\ \Omega(t) = \begin{cases}\Omega_0 & {\rm for}~0<t<T_1, & T_1+\tau<t<T_1+T_2+\tau \\ 0 & {\rm otherwise} \end{cases} \\ \Delta(t) = \begin{cases}\Delta_0 & {\rm for} ~T_1<t<T_1+\tau \\ 0 & {\rm otherwise} \end{cases} \\ $$
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
T_1 = 0.25e-6
T_2 = 0.25e-6
tau = 1e-6
# Pauli matrices
sigma_x = np.array([[0, 1], [1, 0]]) # X matrix
sigma_n = np.array([[0, 0], [0, 1]]) # Z matrix
# Parameters
Omega_0 = (2*np.pi)*1e+6 # 1 MHz
Delta_0 = (2*np.pi)*0.1e+6 # 0.1 MHz
Omega = np.zeros(t_size) # Rabi frequency
Omega[np.logical_and(t_array >= t_ini, t_array < t_ini + T_1)] = Omega_0
Omega[np.logical_and(t_array >= t_ini + T_1 + tau, t_array < t_ini + T_1 + T_2 + tau)] = Omega_0
Delta = np.zeros(t_size) # Detuning
Delta[np.logical_and(t_array >= t_ini + T_1, t_array < t_ini + T_1 + tau)] = Delta_0
# Initial state
c = np.zeros((2, t_size), dtype=np.complex128) # c-vector
c[0, 0] = 1 # Initial condition
# Time evolution using vectorized computation
for idx_time in range(t_size - 1):
# Hamiltonian (time-dependent)
H = (hbar/2) * (Omega[idx_time] * sigma_x + (-2 * Delta[idx_time]) * sigma_n)
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()