Hamiltonian dynamic simulation¶
In this notebook, we will simulate the dynamics of the Hamiltonian system using the basic module, only numpy/scipy/matplotlib.
import numpy as np
import scipy.linalg
import matplotlib.pyplot as plt
time-inependent Hamiltonian¶
We choose single-qubit Hamiltonian as follow for simplicty
$$ H =\frac{1}{2} \omega \sigma_x + \delta\sigma_z$$
and the initial state is
$$ |\psi(0)\rangle = |0\rangle$$
The dynamics of the system is governed by the Schrodinger equation
$$ i\hbar \frac{d}{dt}|\psi(t)\rangle = H|\psi(t)\rangle$$
The solution of the Schrodinger equation is
$$ |\psi(t)\rangle = e^{-iHt/\hbar}|\psi(0)\rangle$$
# parameter
omega = 2*np.pi #radian/microsecond
delta = 0.8*np.pi #radian/microsecond
duration = 2 #microsecond
# time-independent hamiltonian
hamiltonian = np.array([[delta, omega/2], [omega/2, -delta]]) #0.5 omega X + delta Z
time_list = np.linspace(0, duration, 100)
unitary_list = scipy.linalg.expm(-1j*hamiltonian*time_list[:,None,None])
print(hamiltonian)
[[ 2.51327412 3.14159265] [ 3.14159265 -2.51327412]]
initial_state = np.array([1,0]) #|0>
state_list = unitary_list @ initial_state
probability_list = np.abs(state_list)**2
fig,ax = plt.subplots()
ax.plot(time_list, probability_list[:,0], label=r'$|0\rangle$')
ax.plot(time_list, probability_list[:,1], label=r'$|1\rangle$')
ax.set_xlabel('time (microsecond)')
ax.set_ylabel('probability')
ax.legend()
<matplotlib.legend.Legend at 0x7f6bd4262480>
time-dependent Hamiltonian¶
More general Hamiltonian is time-dependent, and to make general quantum computer, we do need time-dependent Hamiltonian. The basic block of a quantum circuit is quantum gate. In this part, we will make a $R_y(\frac{\pi}{2})$ quantum gate.
$$ R_y(\theta) = e^{-i\sigma_y\theta/2}=\begin{pmatrix} \cos\frac{\theta}{2} & -\sin\frac{\theta}{2} \\ \sin\frac{\theta}{2} & \cos\frac{\theta}{2} \end{pmatrix} $$
We choose the Hamiltonian as follow
$$ H=\frac{1}{2}\omega(t)\sigma_y $$
where the frequency $\omega$ is time-dependent. And you can prove that for that $\int_0^T \omega(t) dt=\frac{\pi}{2}$ for a time duration $T$ the evolution operator is $R_y(\frac{\pi}{2})$.
PauliY = np.array([[0,-1j],[1j,0]])
hf_ry = lambda x: scipy.linalg.expm(-1j*PauliY*x/2)
target_gate = hf_ry(np.pi/2)
print(target_gate)
[[ 0.70710678+0.j -0.70710678+0.j] [ 0.70710678+0.j 0.70710678+0.j]]
In numerical calculation, it's more common to slice the time interval into $N$ pieces with the time step is $\Delta t=T/N$ and assume that the Hamiltonian is constant in each time interval.
num_segment = 50
omega_list = 2*np.pi * np.exp(-np.linspace(-3, 3, num_segment)**2) #MHz
duration = 0.5*num_segment*np.pi/np.abs(omega_list).sum() #microsecond
fig,ax = plt.subplots()
tmp0 = np.linspace(0, duration, num_segment)
tmp1 = np.stack([tmp0[:-1],tmp0[1:]], axis=1).reshape(-1)
tmp2 = np.stack([omega_list[:-1],omega_list[:-1]], axis=1).reshape(-1)
ax.plot(tmp1, tmp2)
ax.set_xlabel('time (microsecond)')
ax.set_ylabel(r'frequency $\omega$ (MHz)')
ax.set_title('piecewise constant frequency')
Text(0.5, 1.0, 'piecewise constant frequency')
delta_t = duration/num_segment
unitary_list = [np.eye(2)]
for x in omega_list:
hamiltonian = x * PauliY/2
unitary_list.append(scipy.linalg.expm(-1j * hamiltonian * delta_t) @ unitary_list[-1])
unitary_list = np.stack(unitary_list)
print(unitary_list[-1])
[[ 0.70710678+0.j -0.70710678+0.j] [ 0.70710678+0.j 0.70710678+0.j]]
We can calculate the infedlity of the gate with the target gate
$$ 1 - \frac{\lvert \mathrm{Tr}[R_y(\pi/2)^\dagger U]\rvert^2}{d^2} $$
where $d=2$ is the dimension of the Hilbert space.
infedelity = 1 - abs(np.trace(target_gate.conj().T @ unitary_list[-1]))**2/4
print(infedelity)
4.440892098500626e-16
initial_state = np.array([1,0]) #|0>
state_list = unitary_list @ initial_state
probability_list = np.abs(state_list)**2
time_list = np.linspace(0, duration, num_segment+1)
fig,ax = plt.subplots()
ax.plot(time_list, probability_list[:,0], label=r'$|0\rangle$')
ax.plot(time_list, probability_list[:,1], label=r'$|1\rangle$')
ax.set_xlabel('time (microsecond)')
ax.set_ylabel('probability')
ax.legend()
<matplotlib.legend.Legend at 0x7f6bd42a1160>
