Skip to content

Solution to PSet 3

provided and graded by Xuntao Wu

Download this notebook

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
# Import required libraries
import numpy as np
from math import ceil
import matplotlib.pyplot as plt
from IPython.display import display, Latex
from datetime import datetime
from pytz import timezone

# Install and import QuTiP if needed
try:
    from qutip import *
except ImportError:
    %pip install qutip
    from qutip import *

# Configure matplotlib for LaTeX rendering
plt.rcParams.update(
    {
        "text.usetex": True,
        "text.latex.preamble": r"\usepackage{amsmath} \usepackage{physics}",
        "font.family": "mathpazo",
        "figure.dpi": 300,
        "savefig.dpi": 300,
    }
)
%config InlineBackend.figure_format = 'svg'

# Print version and execution time info
print(f"QuTiP version: {qutip.__version__}")
print(
    "Time of execution: ",
    datetime.now(timezone("America/Chicago")).strftime("%Y-%m-%d %H:%M:%S"),
)

1
2
QuTiP version: 5.1.1
Time of execution:  2025-04-17 14:45:47

For this homework we define the fidelity between a pure state \(|\psi\rangle\) and a general density matrix \(\rho\) as

\[ F = \langle \psi | \rho | \psi \rangle. \]

This differs from the fidelity function in Qutip, but we use this for simplicity.

Problem 3-1 [30 points]

Let's do some error correction!

In this problem, we store one qubit of information, but the physical system suffers from errors. Suppose we model the dynamics of the system with the Lindblad master equation:

\[ \frac{d\rho}{dt} = -i\left(H\rho - \rho H\right) + \sum_j \kappa_j \left[ J_j \rho J_j - \frac{1}{2} \left( J_j^\dagger J_j \rho + \rho J_j^\dagger J_j \right) \right] \]

a. [5 points]

Suppose we have only one qubit initialized in \(\ket{0}\). Ideally we want to keep it unchanged, so \(\hat{H} = I\). However, the system may suffer from stochastic \(X\) errors, which can be modeled with a master equation by choosing \(J = \sigma_x\). We set the error rate \(\kappa = 1\). Use QuTiP to simulate the dynamics of the qubit and plot the fidelity

\[ F(t) = \bra{0} \rho(t) \ket{0} \]

as the qubit evolves.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
### 1.1) Single qubit, waiting

psi_1 = basis(2, 0)
psi_tar = basis(2, 0)

t_list = np.linspace(0, 3, 1000)

# Evolve following master equation
result_1 = mesolve(qeye(2), psi_1, t_list, c_ops=sigmax())

data_1 = [expect(ket2dm(psi_tar), result_1.states[m]) for m in range(len(t_list))]

plt.plot(t_list, data_1)
plt.xlabel("Time")
plt.ylabel("fidelity")

1
Text(0, 0.5, 'fidelity')

svg

b. [5 points]

Now suppose we have three qubits and we initialize them in \(\ket{000}\). Each of the qubits suffers from stochastic \(X\) errors independently, \(J_k = \sigma_{x,k}\), at the same error rate \(\kappa_k = 1\) for \(k = 1,2,3\). Use QuTiP to simulate the dynamics of the 3-qubit system and plot the fidelity

\[ F(t) = \bra{000} \rho(t) \ket{000} \]

as the qubit evolves.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
### 1.2) 3 qubits, waiting

psi_2 = tensor(basis(2, 0), basis(2, 0), basis(2, 0))
psi_tar_2 = psi_2

c_ops_3qb = [
    tensor(sigmax(), qeye(2), qeye(2)),
    tensor(qeye(2), sigmax(), qeye(2)),
    tensor(qeye(2), qeye(2), sigmax()),
]


# Evolve following master equation
result_2 = mesolve(tensor(qeye(2), qeye(2), qeye(2)), psi_2, t_list, c_ops=c_ops_3qb)

data_2 = [expect(ket2dm(psi_tar_2), result_2.states[m]) for m in range(len(t_list))]

plt.plot(t_list, data_1, label="1 qubit")
plt.plot(t_list, data_2, label="3 qubits")
plt.xlabel("Time")
plt.ylabel("fidelity")
plt.legend()
plt.show()

svg

c. [5 points]

After an evolution time \(t\), we would like to perform a "stabilizer measurement" by measuring \(Z_1Z_2\) and \(Z_2Z_3\). This allows us to decide in which of the four subspaces the state is:

  • \(\mathcal{V}_0\) spanned by \(\{\ket{000}, \ket{111}\}\)
  • \(\mathcal{V}_1\) spanned by \(\{\ket{100}, \ket{011}\}\)
  • \(\mathcal{V}_2\) spanned by \(\{\ket{010}, \ket{101}\}\)
  • \(\mathcal{V}_3\) spanned by \(\{\ket{001}, \ket{110}\}\)

If you find that the state ends up in \(\mathcal{V}_k\) instead of \(\mathcal{V}_0\), you can perform a \(\hat{\sigma}_{x,k}\) operation on the \(k\)-th qubit to return it to \(\mathcal{V}_0\). This correction operation can be represented as a quantum channel \(\mathcal{E}\) such that:

\[ \rho_R(t) = \mathcal{E}[\rho(t)] = \sum_{j=0}^{3} R_j \rho(t) R_j^\dagger \]

where \(R_0 = \ket{000}\bra{000} + \ket{111}\bra{111}\) and \(R_j = R_0 \sigma_{x,j}\) for \(j = 1,2,3\).

Plot the fidelity

\[ F(t) = \bra{000} \rho_R(t) \ket{000} \]

after recovery and compare to the result in (a) and (b). Is it better or worse for different \(t\)? Can you briefly explain the reason?

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
### 1.3) 3 qubits, error correction

z0 = basis(2, 0)
z1 = basis(2, 1)

# Define recovery operators
R0 = (
    tensor(z0, z0, z0) * tensor(z0, z0, z0).dag()
    + tensor(z1, z1, z1) * tensor(z1, z1, z1).dag()
)
R1 = (
    tensor(z0, z0, z0) * tensor(z1, z0, z0).dag()
    + tensor(z1, z1, z1) * tensor(z0, z1, z1).dag()
)
R2 = (
    tensor(z0, z0, z0) * tensor(z0, z1, z0).dag()
    + tensor(z1, z1, z1) * tensor(z1, z0, z1).dag()
)
R3 = (
    tensor(z0, z0, z0) * tensor(z0, z0, z1).dag()
    + tensor(z1, z1, z1) * tensor(z1, z1, z0).dag()
)

R_ops = [R0, R1, R2, R3]

# Recovery process
data_3 = []
for rho in result_2.states:
    rho_R = 0
    for k in range(len(R_ops)):
        rho_R += R_ops[k] * rho * R_ops[k].dag()
    data_3.append(expect(ket2dm(psi_tar_2), rho_R))

plt.plot(t_list, data_1, label="1 qubit")
plt.plot(t_list, data_2, label="3 qubits")
plt.plot(t_list, data_3, label="3 qubits recovery")
plt.xlabel("Time")
plt.ylabel("fidelity")
plt.legend()
plt.show()

svg

1
2
3
4
5
6
7
8
9
### In the short time regime, with log-log plot...
plt.loglog(t_list, 1 - np.array(data_1), label="1 qubit")
plt.loglog(t_list, 1 - np.array(data_2), label="3 qubits")
plt.loglog(t_list, 1 - np.array(data_3), label="3 qubits recovery")
plt.xlim(3e-3, 2e-1)
plt.xlabel("Time")
plt.ylabel("1-F")
plt.legend()
plt.show()

svg

We can see that with recovery the scaling (slope) is changed in the small error regime (\(\kappa t \ll 1\)).

d. [5 points]

Now we want to perform gate operations. Suppose first we only have one qubit initialized in \(\ket{+X} = (\ket{0} + \ket{1})/\sqrt{2}\). We want to do a \(Z\) gate with the Hamiltonian \(H = \sigma_z\). If there is no noise, after what time \(\Delta t\) will the state reach \(\ket{-X} = (\ket{0} - \ket{1})/\sqrt{2}\)?

If during the evolution, there is a stochastic \(X\) error with error rate \(\kappa\), how does the final fidelity

\[ F = \bra{-X} \rho_f \ket{-X} \]

change with \(\kappa\)?

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
### 1.4) 1 qubit, gate operation

psi_4 = (basis(2, 0) + basis(2, 1)).unit()
psi_tar_4 = (basis(2, 0) - basis(2, 1)).unit()

Omega = 1.0
H_4 = Omega * sigmaz()

T = np.pi / 2 / Omega
print("Gate time:", T)

kappa_list = np.linspace(0, 3, 50)
data_4 = []
for kappa in kappa_list:
    result_4 = mesolve(
        H_4, psi_4, np.linspace(0, T, 1000), c_ops=np.sqrt(kappa) * sigmax()
    )
    data_4.append(expect(ket2dm(psi_tar_4), result_4.states[-1]))

plt.plot(kappa_list, data_4)
plt.show()

1
Gate time: 1.5707963267948966

svg

e. [5 points]

We now encode one logical qubit in 3 physical qubits in the following way: \(\ket{0_L} = \ket{000}\) and \(\ket{1_L} = \ket{111}\). We initialize the logical qubit in

\[ \ket{(+X)_L} = \frac{1}{\sqrt{2}}(\ket{0_L} + \ket{1_L}) = \frac{1}{\sqrt{2}}(\ket{000} + \ket{111}), \]

and use the Hamiltonian

\[ \hat{H} = \frac{1}{3}(\hat{\sigma}_{z,1} + \hat{\sigma}_{z,2} + \hat{\sigma}_{z,3}) \]

that lasts for time \(\Delta t\) to drive the state toward

\[ \ket{(-X)_L} = \frac{1}{\sqrt{2}}(\ket{0_L} - \ket{1_L}). \]

However, as in (b), each physical qubit suffers from independent \(X\) errors at the same rate \(\kappa\). In this case:

  • How does the final fidelity

    \[ F = \bra{(-X)_L} \rho_f \ket{(-X)_L} \]

    change with \(\kappa\) (make a plot in QuTiP)? - How does it compare to the results in (d)? - How does this change if we perform the recovery channel described in (c) for the quantum state \(\rho_f\) we get in the end?

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
### 1.5) 3 qubits, gate operation

psi_5 = (
    tensor(basis(2, 0), basis(2, 0), basis(2, 0))
    + tensor(basis(2, 1), basis(2, 1), basis(2, 1))
).unit()
psi_tar_5 = (
    tensor(basis(2, 0), basis(2, 0), basis(2, 0))
    - tensor(basis(2, 1), basis(2, 1), basis(2, 1))
).unit()

Omega = 1.0
H_5 = (
    Omega
    * (
        tensor(sigmaz(), qeye(2), qeye(2))
        + tensor(qeye(2), sigmaz(), qeye(2))
        + tensor(qeye(2), qeye(2), sigmaz())
    )
    / 3
)
T = np.pi / 2 / Omega

data_51 = []
rho_f_5 = []  # store each final states achieved with corresponding kappa
for kappa in kappa_list:
    c_ops_3qb = [
        tensor(sigmax(), qeye(2), qeye(2)) * np.sqrt(kappa),
        tensor(qeye(2), sigmax(), qeye(2)) * np.sqrt(kappa),
        tensor(qeye(2), qeye(2), sigmax()) * np.sqrt(kappa),
    ]
    result_5 = mesolve(H_5, psi_5, np.linspace(0, T, 1000), c_ops=c_ops_3qb)
    rho_f_5.append(result_5.states[-1])
    data_51.append(expect(ket2dm(psi_tar_5), result_5.states[-1]))

# Recovery process
data_52 = []
for rho in rho_f_5:
    rho_R = 0
    for k in range(len(R_ops)):
        rho_R += R_ops[k] * rho * R_ops[k].dag()
    data_52.append(expect(ket2dm(psi_tar_5), rho_R))

plt.plot(kappa_list, data_4, label="1 qubit")
plt.plot(kappa_list, data_51, label="3 qubits")
plt.plot(kappa_list, data_52, label="3 qubits recovery")
plt.xlabel("kappa")
plt.ylabel("fidelity")
plt.legend()
plt.show()

svg

We can see that with bit-flip correction in the end cannot beat the 1-qubit case without protection. Of course one subtle thing is that we use the same gate time in both 1-qubit and 3-qubit cases. In fact, in this case the relative phase between \(|0\rangle\) and \(|1\rangle\) (or \(|0\rangle_L\) and \(|1\rangle_L\)) of the final state will depends on when the error happens. So, it is important for the next problem to know the time information of the errors.

Here is an easy way to understand that:

For \(H = (Z_1 + Z_2 + Z_3)/3\), if restricted in \(\{|000\rangle, |111\rangle\}\) subspace, the Hamiltonian will be diag(1,-1), while if a bit-flip \(X_1\) happens at some time, the computational basis will change into \(\{|100\rangle, |011\rangle\}\), and the Hamiltonian in that subspace will be diag(1/3,-1/3), which cannot give you the correct phase after evolution.

f. [Bonus! – 3 extra pts]

Suppose for the 3-qubit case, we split the gate time \(\Delta t\) into \(M\) intervals, each with equal length (say \(M = 100\)). At the end of each interval

\[ t_k = \frac{k \Delta t}{M}, \]

you quickly perform the stabilizer measurement to check whether an \(X\) error happened within \(t \in [t_{k-1}, t_k]\) and on which qubit the error occurred.

You may run into trouble when two or more qubits have \(X\) errors in this interval. For instance:

  • \(X_1 X_2\) will be misidentified as an \(X_3\) error
  • \(X_1 X_2 X_3\) will be misidentified as no error

Furthermore, you do not know the exact time an error occurred, so if you find an error in \([t_{k-1}, t_k]\), you may decide that it happened at

\[ t = \frac{t_{k-1} + t_k}{2}. \]

Now that you approximately know the time and place of each \(X_i\) error during the gate operation, can you design better recovery strategies on the final state in order to achieve higher final fidelity?

You may:

  • Use the idea from quantum trajectories
  • Use the mcsolve() function in QuTiP
  • For each trajectory, find out when (Result.col_times) and where (Result.col_which) an error happens
  • Assume in your recovery strategy that you do not have perfect knowledge in practice (e.g., place the error at midpoint of interval)

Then:

  • Apply your recovery strategy (assuming operations are perfect) on each trajectory
  • Average all outcomes to get the final fidelity

Discuss the improvement by making use of your protocol.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
### 1.6) 3 qubits, with real-time syndrome monitoring
### Takes several minutes to run...

M = 100

kappa_list_6 = np.linspace(0, 3, 7)

psi_000 = tensor(z0, z0, z0)
psi_111 = tensor(z1, z1, z1)

X = [
    tensor(sigmax(), qeye(2), qeye(2)),
    tensor(qeye(2), sigmax(), qeye(2)),
    tensor(qeye(2), qeye(2), sigmax()),
]

data_6 = []

for kappa in kappa_list_6:
    c_ops_3qb = [
        tensor(sigmax(), qeye(2), qeye(2)) * np.sqrt(kappa),
        tensor(qeye(2), sigmax(), qeye(2)) * np.sqrt(kappa),
        tensor(qeye(2), qeye(2), sigmax()) * np.sqrt(kappa),
    ]

    # Use mcsolve to get the time information for errors in each trajectory
    result_6 = mcsolve(
        H_5,
        psi_5,
        np.linspace(0, T, 1000),
        c_ops=c_ops_3qb,
        ntraj=1000,
        options={"keep_runs_results": True, "progress_bar": False},
    )

    fidel = 0.0
    for k in range(result_6.num_trajectories):
        # Calculate the syndrome you will find in each mid-detection
        syndrome = np.zeros((M, 3))
        for count, time in enumerate(result_6.col_times[k]):
            piece_ind = int(time / (T / M))
            syndrome[piece_ind][result_6.col_which[k][count]] += 1
        syndrome = syndrome % 2

        # However, you may mis-identify that if 2 or 3 qubits suffer a quantum jump within same time bin
        for k2 in range(M):
            if sum(syndrome[k2]) > 1:
                syndrome[k2] = (syndrome[k2] + 1) % 2

        # Count the actual phase rotated due to the syndrome
        phase = 0.0
        total_sdrm = np.array([0, 0, 0])
        t_last_error = 0
        for k2 in range(M):
            if sum(syndrome[k2]) > 0:
                t_error = T / M * (k2 + 0.5)
                phase += (
                    2 * (1.0 - 2.0 / 3 * sum(total_sdrm)) * (t_error - t_last_error)
                )
                total_sdrm = (total_sdrm + syndrome[k2]) % 2
                t_last_error = t_error
        phase += 2 * (1.0 - 2.0 / 3 * sum(total_sdrm)) * (T - t_last_error)

        # Do the recovery
        #         psi_f = result_6.states[k][-1]
        #         psi_f = Qobj(psi_f)
        #         psi_f.dims = [[2, 2, 2], [1]]
        psi_f = result_6.states[k][-1]
        for ind in range(3):
            if total_sdrm[ind] == 1:
                psi_f = X[ind] * psi_f  # Do bit recovery

        rec_op = ket2dm(psi_000) + np.exp(1j * (np.pi - phase)) * ket2dm(psi_111)
        psi_f = rec_op * psi_f  # Do phase recovery

        fidel += expect(
            ket2dm(psi_tar_5), psi_f
        )  # calculate the fidelity after the recovery

    fidel /= result_6.num_trajectories  # Average over all trajectories
    data_6.append(fidel)

plt.plot(kappa_list, data_4, label="1 qubit")
plt.plot(kappa_list, data_51, label="3 qubits")
plt.plot(kappa_list, data_52, label="3 qubits final rec")
plt.plot(kappa_list_6, data_6, label="3 qubits mid-monitor")
plt.xlabel("kappa")
plt.ylabel("fidelity")
plt.legend()
plt.show()

svg

Problem 3-2 [25 points]

Long coherence quantum dots in silicon.

For this problem we will have a closer look at quantum dots realized in silicon, based on the UNSW Dzurak group's publication Veldhorst et al., "An addressable quantum dot qubit with fault-tolerant control-fidelity," Nature Nanotech. 9, 981–985 (2014). We will explore the device design and operation, and will look into gate fidelities and randomized benchmarking (Figure 4 in the paper) at a later stage in the class.

The device is based on a MOS (metal-oxide-semiconductor) silicon structure. Here, the electrons sit at the interface of silicon and silicon dioxide (which naturally forms on top of silicon). Unlike the 2DEGs that we discussed in class, without any gate voltage applied, there are no electrons at the interface. However, when a positive gate voltage is applied, electrons accumulate below the gate and form a 2DEG.

You can find more background on MOS silicon structures in: F. A. Zwanenburg et al., Rev. Mod. Phys. 85, 961–1019 (2013), especially on p. 980, and also in ref. 19 of Veldhorst et al. (it is not needed to solve the questions below).

a. [5 points]

Let's learn about the device (Fig 1a in Veldhorst et al.):

(i) Specify the signs of the voltages applied to the gates R, G1, G2, G3, G4 and C. Why do these define a dot?

(ii) Which gate sets the chemical potential of the dot?

(iii) Sketch a cross-section through the device. Indicate the gates, magnetic fields, electrostatic potentials, and the location of the electrons.

(iv) What sets the tunneling rate from the dot to the reservoir?

(v) Estimate the value of \(B_1\) at the dot based on the Rabi frequency that the authors measure in Fig. 2b. What current \(I_{\mathrm{ESR}}\) does this correspond to? Does this match the reported \(P_{\mathrm{ESR}}\) in the paper?

(vi) Assume that the dot has a vertical confinement of 20 nm. Estimate how many nuclear spins the electron spin of the dot interacts with. How does this compare to dots defined in GaAs?

(i) \(R\) is grounded. \(G_1, G_2>0\), because in Supp.1 they mention "In the single quantum dot mode, G2 and G1 are set to a high potential, such that a continuous 2-DEG is formed underneath these gates". \(G_3<0\), as the quantum dot tunnel-coupled to \(R\) via \(G_3\), meaning \(G_3\) forms a electron barrier. \(G_4>0\) and \(C<0\) to form and confine a dot.

(ii) \(G_4\) sets the chemical potential of the dot.

(iii)PS3_21c_RyanWhite.png

(iv) \(G_3\) sets the tunneling rate.

(v) There are \(2.5\) periods during \(\tau_p\simeq8.56~\mu\mathrm{s}\), so the Rabi frequency \(\Omega\) satisfies \(\Omega\cdot\frac{8.56\mu\mathrm{s}}{2.5}=2\pi\), which means \(\Omega=2\pi\times292~\mathrm{kHz}\). Assuming the magnetic field takes the form of \(B_1\cos{\omega t}\), we have \(\frac{1}{2}\frac{\mu_B g B_1}{\hbar}=\Omega\), therefore \(B_1\simeq21.1~\mu\mathrm{T}\). Then, we have \(\frac{\mu_0 I}{2\pi r}=B_1\), from which we obtain \(I=31.6~\mu\mathrm{A}\) assuming \(r\sim300~\mathrm{nm}\). From the fitting in Fig. 2c, \(f_{\mathrm{Rabi}}=292~\mathrm{kHz}\) corresponds to \(P_{\mathrm{ESR}}\simeq4~\mathrm{mW}\), which is similar to \(P_{\mathrm{ESR}}=5~\mathrm{dBm}=3.2~\mathrm{mW}\) reported in the paper.

(vi) Let's say the planar dimension of the dot is \(100~\mathrm{nm}\times100~\mathrm{nm}\), it will roughly contain

\[ \frac{100~\mathrm{nm}\times100~\mathrm{nm}\times20~\mathrm{nm}\times\rho_{Si}}{28~\mathrm{g}/N_A} \simeq 1.0\times10^7 \]

\({}^{28}\mathrm{Si}\) atoms. However, only \(800\) ppm of them are \({}^{29}\mathrm{Si}\), which have nuclear spin \(1/2\). So the number of \({}^{29}\mathrm{Si}\) is roughly

\[ 1.0\times10^{7}\times800/10^{6}=8000 \]

b. [5 points]

Initialization:

(i) Given the energy scales set by the temperature and by the magnetic field, what is the probability to find a single, completely isolated electron spin in the up state?

(ii) How does this compare to the initialization fidelity of the dot's spin in the paper? What could be the reasons for the difference (if there is a difference)?

(i) \(p_{\uparrow}=e^{-\Delta E/k_B T}=4.75\times10^{-17}\)

(ii) The paper reports initialization fidelity of only \(95\%\) and claimed to be limited by thermal broadening of the electron reservoir. There are chances that \(E_e>\mu\), therefore leads to random tunneling events between the reservoir and the dot.

c. [5 points]

Readout:

(i) Figure 1d and Supplementary Figure 1 characterize the single-shot readout of the qubit. Qualitatively sketch the pulses that you would apply to G4 for the "load", "readout", and "empty" phases.

(ii) Based on Supplementary Figure 1, what would be a good "read level" for the readout?

PS3_23.png

d. [5 points]

Rabi oscillations:

(i) Take parameters from the paper to simulate Rabi oscillations as a function of detuning using QuTiP (Figure 2d). Can you reproduce the behavior shown in Figure 2d?

(ii) How would that behavior change if you assume a 10× faster dephasing rate than the authors measure? (Please illustrate by simulation in QuTiP.)

(i) The Hamiltonian is the same as in HW 2. (Rabi driven qubit with detuning)

1
2
3
4
5
6
7
8
T = 8.56e-6 / 2.5  # From the caption of figure 2.b

rabi = 2 * np.pi / T  # Hz


# Same Hamiltonian as HW2
def H_rot(delta):
    return -delta / 2 * sigmaz() - rabi / 2 * sigmay()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
t_list = np.linspace(0, 30e-6, 1000)  # seconds

# definition of basis states
psi_up = basis(2, 0)
psi_down = basis(2, 1)
rho_up = psi_up * psi_up.dag()
rho_down = psi_down * psi_down.dag()

# Initial state, with fidelity of 95%
rho_0 = ((np.sqrt(0.05) * psi_up + np.sqrt(0.95) * psi_down).unit()) * (
    (np.sqrt(0.05) * psi_up + np.sqrt(0.95) * psi_down).unit()
).dag()  # Initial state with 95% fidelity

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
# Detuning list, from -1 to 1 MHz
detuning_list = 2 * np.pi * np.linspace(-1000000, 1000000, 250)

# Dephasing time as introduced in figure 3a
T2star = 120 * 1e-6  # seconds
for delta in detuning_list:
    result_list = mesolve(
        -delta / 2 * sigmaz() - rabi / 2 * sigmay(),
        rho_0,
        t_list,
        [np.sqrt(1 / (2 * T2star)) * sigmaz()],
    ).states  # result list, for a given delta
    result_delta = np.array(
        []
    )  # This is a list that will contain the fidelities of the state with rho_0 for a given delta
    for i in range(len(t_list)):
        result_delta = np.append(
            result_delta, fidelity(result_list[i], psi_up) ** 2
        )  # spin up fraction is given by the fidelity between state and psi_up

    # then we add those values to a final result array
    if delta == detuning_list[0]:
        result_array = result_delta
    else:
        result_array = np.vstack((result_array, result_delta))

# So we end up with a 2d array. x axis for time dependance, y axis for detuning dependance

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
# Plot of the result

A = 0.7  # Visibility from the paper

plt.title(
    r"Spin up fraction during a rabi pulse as a function of the pulse duration $\tau_p$ and the detuning $\nu_{ESR}-\nu_0$"
)
plt.ylabel(r"$\nu_{ESR} - \nu_0$ (MHz)")
plt.xlabel(r"$\tau_p$ ($\mu$s) ")
plt.contourf(
    t_list * 1e6, detuning_list * 1e-6 / (2 * np.pi), A * result_array, 20, cmap="hot"
)
cbar = plt.colorbar(label=r"$f_{\uparrow}$", location="right", shrink=0.5)
plt.show()

svg

(ii) 10 x faster dephasing. We can see in the plot that the effect is to drag state towards the middle of the Bloch sphere.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
# exact same thing but with dephasing time 10x smaller.

result_array_dephasing = np.array([])
detuning_list_dephasing = 2 * np.pi * np.linspace(-1000000, 1000000, 500)
T2star_smaller = 0.1 * 120 * 1e-6  # seconds
c_ops = [np.sqrt(1 / (2 * T2star_smaller)) * sigmaz()]
for delta in detuning_list_dephasing:
    result_list_dephasing = mesolve(
        -delta / 2 * sigmaz() - rabi / 2 * sigmay(),
        rho_0,
        t_list,
        c_ops,
        options={"max_step": 1e-6, "nsteps": 1000},
    ).states
    result_delta_dephasing = np.array([])
    for i in range(len(t_list)):
        result_delta_dephasing = np.append(
            result_delta_dephasing, fidelity(result_list_dephasing[i], psi_up) ** 2
        )
    if delta == detuning_list_dephasing[0]:
        result_array_dephasing = result_delta_dephasing
    else:
        result_array_dephasing = np.vstack(
            (result_array_dephasing, result_delta_dephasing)
        )

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
A = 0.7  # Visibility from the paper

plt.title(
    r"Spin up fraction during a rabi pulse as a function of the pulse duration $\tau_p$ and the detuning $\nu_{ESR}-\nu_0$"
)
plt.ylabel(r"$\nu_{ESR} - \nu_0$ (MHz)")
plt.xlabel(r"$\tau_p$ ($\mu$s) ")
plt.contourf(
    t_list * 1e6,
    detuning_list_dephasing * 1e-6 / (2 * np.pi),
    A * result_array_dephasing,
    20,
    cmap="hot",
)
cbar = plt.colorbar(label=r"$f_{\uparrow}$", location="right", shrink=0.5)
plt.show()

svg

e. [5 points]

Coherence times and extending coherence times (Fig. 3):

(i) In Figure 3a the authors perform a Ramsey experiment (similar to Pset 2, Problem 2) to determine the dephasing time \(T_2^*\). A slight detuning between the microwave source and the spin causes the qubit to pick up a phase in the rotating frame of the source. Verify this behavior by performing a Ramsey experiment in QuTiP with the microwave on-resonance and slightly detuned from resonance. What sets the oscillation frequency?

(ii) Now add dephasing to your Ramsey experiment. Can you reproduce Figure 3a?

(iii) The coherence time of a qubit is not necessarily limited by \(T_2^*\) but can often be extended by an echo sequence. Describe briefly and only qualitatively what the difference between Figure 3b and 3c is and why the coherence time in 3c is longer.

(i) We can see that if the detuning is 0, then the Ramsey sequence doesn't do anything to the state. For detuning, we can see that it oscillates like expected, furthermore, the more we increase the detuning, the greater the frequency gets. Also from my two points, this relation seems to be proportional.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
T_half_pi = 0.25 * 2 * np.pi / rabi


# Envelope function as in HW2
def ramsey_envelope(t, args):
    if (t > T_half_pi and t < args["tau"] + T_half_pi) or t > 2 * T_half_pi + args[
        "tau"
    ]:
        envelope = 0
    else:
        envelope = 1
    return envelope


# List of tau values
tau_list = np.linspace(0, 120e-6, 250)

# For these, I manually put a *0 in the c_ops to remove dephasing
# On resonance
delta = 0
H = [-delta / 2 * sigmaz(), [-rabi / 2 * sigmay(), ramsey_envelope]]
for i in range(len(tau_list)):
    t_list_ramsey = np.linspace(0, 2 * T_half_pi + tau_list[i], len(tau_list))
    result = mesolve(
        H,
        rho_0,
        t_list_ramsey,
        [np.sqrt(1 / (2 * T2star)) * 0 * sigmaz()],
        args={"tau": tau_list[i]},
        options={"max_step": 1e-7, "nsteps": 10000},
    ).states
    if i == 0:
        result_ramsey_list_on_res = fidelity(result[-1], psi_up) ** 2
    else:
        result_ramsey_list_on_res = np.append(
            result_ramsey_list_on_res, fidelity(result[-1], psi_up) ** 2
        )

# a little bit of resonance: 314 kHZ
delta = 2 * np.pi * 50000
H = [-delta / 2 * sigmaz(), [-rabi / 2 * sigmay(), ramsey_envelope]]
for i in range(len(tau_list)):
    t_list_ramsey = np.linspace(0, 2 * T_half_pi + tau_list[i], len(tau_list))
    result = mesolve(
        H,
        rho_0,
        t_list_ramsey,
        [np.sqrt(1 / (2 * T2star)) * 0 * sigmaz()],
        args={"tau": tau_list[i]},
        options={"max_step": 1e-7, "nsteps": 10000},
    ).states
    if i == 0:
        result_ramsey_list_detuned1 = fidelity(result[-1], psi_up) ** 2
    else:
        result_ramsey_list_detuned1 = np.append(
            result_ramsey_list_detuned1, fidelity(result[-1], psi_up) ** 2
        )

# a little bit of resonance: 628 kHZ
delta = 2 * np.pi * 100000
H = [-delta / 2 * sigmaz(), [-rabi / 2 * sigmay(), ramsey_envelope]]
for i in range(len(tau_list)):
    t_list_ramsey = np.linspace(0, 2 * T_half_pi + tau_list[i], len(tau_list))
    result = mesolve(
        H,
        rho_0,
        t_list_ramsey,
        [np.sqrt(1 / (2 * T2star)) * 0 * sigmaz()],
        args={"tau": tau_list[i]},
        options={"max_step": 1e-7, "nsteps": 10000},
    ).states
    if i == 0:
        result_ramsey_list_detuned2 = fidelity(result[-1], psi_up) ** 2
    else:
        result_ramsey_list_detuned2 = np.append(
            result_ramsey_list_detuned2, fidelity(result[-1], psi_up) ** 2
        )

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# plot of the whole thing
plt.title(
    r"Probability of being in spin up state at the end of a"
    "\n"
    r"Ramsay sequence $\frac{\pi}{2} - \frac{\tau}{2} - \frac{\pi}{2}$, as a function of $\tau$"
)
plt.xlabel(r"$\tau$ (ms)")
plt.ylabel(r"$f_{\uparrow}$")
plt.plot(tau_list * 1000, result_ramsey_list_on_res, color="blue", label="On resonance")
plt.plot(
    tau_list * 1000,
    result_ramsey_list_detuned1,
    color="lightgreen",
    label=r"$\nu - \nu_0$ = 314 kHz",
)
plt.plot(
    tau_list * 1000,
    result_ramsey_list_detuned2,
    color="hotpink",
    label=r"$\nu - \nu_0$ = 628 kHz",
)
plt.legend()
plt.show()

svg

(ii) Now add dephasing

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
# same thing, but by adding manually the T2star = 120microseconds

tau_list = np.linspace(0, 120e-6, 250)
delta = 0
H = [-delta / 2 * sigmaz(), [-rabi / 2 * sigmay(), ramsey_envelope]]
for i in range(len(tau_list)):
    t_list_ramsey = np.linspace(0, 2 * T_half_pi + tau_list[i], len(tau_list))
    result = mesolve(
        H,
        rho_0,
        t_list_ramsey,
        [np.sqrt(1 / (2 * T2star)) * sigmaz()],
        args={"tau": tau_list[i]},
        options={"max_step": 1e-7, "nsteps": 10000},
    ).states
    if i == 0:
        result_ramsey_list_on_res = fidelity(result[-1], psi_up) ** 2
    else:
        result_ramsey_list_on_res = np.append(
            result_ramsey_list_on_res, fidelity(result[-1], psi_up) ** 2
        )


delta = 2 * np.pi * 50000
H = [-delta / 2 * sigmaz(), [-rabi / 2 * sigmay(), ramsey_envelope]]
for i in range(len(tau_list)):
    t_list_ramsey = np.linspace(0, 2 * T_half_pi + tau_list[i], len(tau_list))
    result = mesolve(
        H,
        rho_0,
        t_list_ramsey,
        [np.sqrt(1 / (2 * T2star)) * sigmaz()],
        args={"tau": tau_list[i]},
        options={"max_step": 1e-7, "nsteps": 10000},
    ).states
    if i == 0:
        result_ramsey_list_detuned1 = fidelity(result[-1], psi_up) ** 2
    else:
        result_ramsey_list_detuned1 = np.append(
            result_ramsey_list_detuned1, fidelity(result[-1], psi_up) ** 2
        )

delta = 2 * np.pi * 100000
H = [-delta / 2 * sigmaz(), [-rabi / 2 * sigmay(), ramsey_envelope]]
for i in range(len(tau_list)):
    t_list_ramsey = np.linspace(0, 2 * T_half_pi + tau_list[i], len(tau_list))
    result = mesolve(
        H,
        rho_0,
        t_list_ramsey,
        [np.sqrt(1 / (2 * T2star)) * sigmaz()],
        args={"tau": tau_list[i]},
        options={"max_step": 1e-7, "nsteps": 10000},
    ).states
    if i == 0:
        result_ramsey_list_detuned2 = fidelity(result[-1], psi_up) ** 2
    else:
        result_ramsey_list_detuned2 = np.append(
            result_ramsey_list_detuned2, fidelity(result[-1], psi_up) ** 2
        )

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
plt.title(
    r"Probability of being in spin up state at the end of a"
    "\n"
    r"Ramsay sequence $\frac{\pi}{2} - \frac{\tau}{2} - \frac{\pi}{2}$, as a function of $\tau$ with added dephasing of $T_2^* = 0.12$ms"
)
plt.xlabel(r"$\tau$ (ms)")
plt.ylabel(r"$f_{\uparrow}$")
plt.plot(tau_list * 1000, result_ramsey_list_on_res, color="blue", label="On resonance")
plt.plot(
    tau_list * 1000,
    result_ramsey_list_detuned1,
    color="lightgreen",
    label=r"$\nu - \nu_0$ = 314 kHz",
)
plt.plot(
    tau_list * 1000,
    result_ramsey_list_detuned2,
    color="hotpink",
    label=r"$\nu - \nu_0$ = 628 kHz",
)
plt.legend()
plt.show()

svg

Trying to get the same figure as in the paper, figure 3a

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
# same thing, but by adding manually the T2star = 120microseconds

tau_list = np.linspace(0, 120e-6, 250)
delta = 2 * np.pi * 45000
H = [-delta / 2 * sigmaz(), [-rabi / 2 * sigmay(), ramsey_envelope]]
for i in range(len(tau_list)):
    t_list_ramsey = np.linspace(0, 2 * T_half_pi + tau_list[i], len(tau_list))
    result = mesolve(
        H,
        rho_0,
        t_list_ramsey,
        [np.sqrt(1 / (2 * T2star)) * sigmaz()],
        args={"tau": tau_list[i]},
        options={"max_step": 1e-7, "nsteps": 10000},
    ).states
    if i == 0:
        result_ramsey_list = fidelity(result[-1], psi_up) ** 2
    else:
        result_ramsey_list = np.append(
            result_ramsey_list, fidelity(result[-1], psi_up) ** 2
        )

plt.title(
    r"Reproduction of figure 3.a in the paper. With a manually added $T_2^* = 0.12ms$ and a detuning on $45kHz$"
)
plt.xlabel(r"$\tau$ (ms)")
plt.ylabel(r"$f_{\uparrow}$")
plt.plot(tau_list * 1000, result_ramsey_list, color="blue")

# plt.legend()
plt.show()

svg

(iii) The difference is how many \(\pi\) pulses are applied between the two \(\pi/2\) pulses. The reason why figure c yields longer coherence is that when more \(\pi\) pulses are applied, it samples only the high frequency component from the noise environment, which typically has smaller amplitudes.