Evolution Under a Static Hamiltonian¶

In [1]:

from tins import spinsys
import numpy as np
from scipy.linalg import expm
np.set_printoptions(precision=2, suppress=True)

from tins import spinsys import numpy as np from scipy.linalg import expm np.set_printoptions(precision=2, suppress=True)

Start and Detect Operators¶

In [2]:

I = spinsys(0.5)

start = I.z
detect = I.x + 1j * I.y

I = spinsys(0.5) start = I.z detect = I.x + 1j * I.y

Hamiltonian¶

In [3]:

ωrf = 2 * np.pi * 10e3
Ham = ωrf * I.x
time = 25e-6

ωrf = 2 * np.pi * 10e3 Ham = ωrf * I.x time = 25e-6

Propagator¶

In [4]:

propagator = expm(-1j * Ham * time)

propagator = expm(-1j * Ham * time)

Evolution (Single Step)¶

In [5]:

end = propagator @ start @ propagator.conj().transpose()
print(end)

end = propagator @ start @ propagator.conj().transpose() print(end)

[[0.+0.j  0.+0.5j]
 [0.-0.5j 0.+0.j ]]

Detect¶

In [6]:

print(np.trace(end @ detect))

print(np.trace(end @ detect))

-0.49999999999999994j

In [7]:

print(np.trace(end @ detect))

print(np.trace(end @ detect))

-0.49999999999999994j

In [8]:

np.sum(end * detect)

np.sum(end * detect)

Out[8]:

np.complex128(0.49999999999999994j)

In [9]:

detect @ end

detect @ end

Out[9]:

array([[0.-0.5j, 0.+0.j ],
       [0.+0.j , 0.+0.j ]])

Nutation¶

In [10]:

start = I.x
detect = I.p
ωrf = 2 * np.pi * 10e3
Ham = ωrf * I.y

total_time = 5e-3
npoints = 1000
dwell_time = total_time / (npoints - 1)

prop = expm(-1j * Ham * dwell_time)
prop_dag = prop.conjugate().transpose()

start = I.x detect = I.p ωrf = 2 * np.pi * 10e3 Ham = ωrf * I.y total_time = 5e-3 npoints = 1000 dwell_time = total_time / (npoints - 1) prop = expm(-1j * Ham * dwell_time) prop_dag = prop.conjugate().transpose()

In [11]:

out = []
for t in range(npoints):
    out.append(start)
    end = prop @ start @ prop_dag
    start = end

out = [] for t in range(npoints): out.append(start) end = prop @ start @ prop_dag start = end

In [12]:

fid = np.array([np.trace(state @ detect) for state in out])

fid = np.array([np.trace(state @ detect) for state in out])

In [13]:

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(5, 2))
ax.plot(fid.imag)
ax.plot(fid.real)

import matplotlib.pyplot as plt fig, ax = plt.subplots(figsize=(5, 2)) ax.plot(fid.imag) ax.plot(fid.real)

Out[13]:

[<matplotlib.lines.Line2D at 0x7d8e5bf362a0>]

Fourier Transform¶

In [14]:

ft = np.fft.fftshift(np.fft.fft(fid))


ft_broadened = fid * np.exp(-np.linspace(0, total_time, npoints)/1e-3)
ft_broadened = np.fft.fftshift(np.fft.fft(ft_broadened))

ft = np.fft.fftshift(np.fft.fft(fid)) ft_broadened = fid * np.exp(-np.linspace(0, total_time, npoints)/1e-3) ft_broadened = np.fft.fftshift(np.fft.fft(ft_broadened))

In [15]:

fig, ax = plt.subplots(figsize=(5, 3), constrained_layout=True)
ax.plot(ft_broadened.real)

plt.show()

fig, ax = plt.subplots(figsize=(5, 3), constrained_layout=True) ax.plot(ft_broadened.real) plt.show()

Q1¶

• Remake this plot with the correct x-axis (in Hz)
• Do this calculation for a different rf amplitude.
• Modify the rf hamiltonian so that it can take an arbitrary phase.
• How do you modify the code if all you need is the final state?

In [16]:

# your answers here

# your answers here

For Later Use¶

In [17]:

import numpy as np
from scipy.linalg import expm


def mesolve_single(start, detect, ham, time):
    prop = expm(-1j * ham * time)
    prop_dag = prop.conjugate().transpose()

    end = prop @ start @ prop_dag

    if detect is not None:
        return np.trace(end @ detect)
    else:
        return end



def mesolve(start, detect, ham, time):
    dt = time[1] - time[0]
    prop = expm(-1j * ham * dt)
    prop_dag = prop.conjugate().transpose()

    out = []
    for t in time:
        out.append(start)
        end = prop @ start @ prop_dag
        start = end

    if detect is not None:
        return np.array([np.trace(state @ detect) for state in out])
    else:
        return np.array(out)

import numpy as np from scipy.linalg import expm def mesolve_single(start, detect, ham, time): prop = expm(-1j * ham * time) prop_dag = prop.conjugate().transpose() end = prop @ start @ prop_dag if detect is not None: return np.trace(end @ detect) else: return end def mesolve(start, detect, ham, time): dt = time[1] - time[0] prop = expm(-1j * ham * dt) prop_dag = prop.conjugate().transpose() out = [] for t in time: out.append(start) end = prop @ start @ prop_dag start = end if detect is not None: return np.array([np.trace(state @ detect) for state in out]) else: return np.array(out)