#DOI : https://doi.org/10.1143/JPSJ.12.570



import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import expm

# Define Pauli matrices
sigma_x = np.array([[0, 1], [1, 0]])
sigma_y = np.array([[0, -1j], [1j, 0]])
sigma_z = np.array([[1, 0], [0, -1]])

# System parameters
omega_0 = 1.0  # Transition frequency
Omega = 0.1    # Driving strength
omega = 1.0    # Driving frequency
hbar = 1.0

# Time array
t = np.linspace(0, 50, 10000)
dt = t[1] - t[0]

# Initial density matrix (ground state)
rho = np.array([[1, 0], [0, 0]], dtype=complex)

# Hamiltonian matrices
H0 = (hbar * omega_0 / 2) * sigma_z
H_prime = lambda t: hbar * Omega * np.cos(omega * t) * sigma_x

# Time evolution
rho_t = np.zeros((len(t), 2, 2), dtype=complex)
rho_t[0] = rho

for i in range(1, len(t)):
    H = H0 + H_prime(t[i-1])
    U = expm(-1j * H * dt / hbar)
    rho = U @ rho @ U.conj().T
    rho_t[i] = rho

# Calculate expectation values of Pauli z-matrix
expectation_z = [np.trace(rho @ sigma_z).real for rho in rho_t]

# Plot the results
plt.figure(figsize=(10, 6))
plt.plot(t, expectation_z, label=r'$\langle \sigma_z \rangle$')
plt.xlabel('Time $t$')
plt.ylabel(r'$\langle \sigma_z \rangle$')
plt.title('Response of a Two-Level System to a Sinusoidal Driving Field')
plt.legend()
plt.grid()
plt.show()


import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import expm

# Define Pauli matrices
sigma_x = np.array([[0, 1], [1, 0]])
sigma_y = np.array([[0, -1j], [1j, 0]])
sigma_z = np.array([[1, 0], [0, -1]])

# System parameters
omega_0 = 1.0  # Transition frequency
Omega = 0.1    # Driving strength
hbar = 1.0

# Time array
t = np.linspace(0, 50, 10000)
dt = t[1] - t[0]

# Initial density matrix (ground state)
rho = np.array([[1, 0], [0, 0]], dtype=complex)

# Hamiltonian matrices
H0 = (hbar * omega_0 / 2) * sigma_z
H_prime = lambda t: hbar * Omega * (t >= 0) * sigma_x  # Step function perturbation

# Time evolution
rho_t = np.zeros((len(t), 2, 2), dtype=complex)
rho_t[0] = rho

for i in range(1, len(t)):
    H = H0 + H_prime(t[i-1])
    U = expm(-1j * H * dt / hbar)
    rho = U @ rho @ U.conj().T
    rho_t[i] = rho

# Calculate expectation values of Pauli z-matrix
expectation_z = [np.trace(rho @ sigma_z).real for rho in rho_t]

# Plot the results
plt.figure(figsize=(10, 6))
plt.plot(t, expectation_z, label=r'$\langle \sigma_z \rangle$')
plt.xlabel('Time $t$')
plt.ylabel(r'$\langle \sigma_z \rangle$')
plt.title('Response of a Two-Level System to a Step Function Driving Field')
plt.legend()
plt.grid()
plt.show()
