Qubit Protected by Two-Cooper-Pair Tunneling

In this notebook, we try to reproduce the result of the “Superconducting circuit protected by two-Cooper-pair tunneling” paper.

Introdcution

Smith2020 designed a qubit that is protected by two Cooper-pair tunneling. We reproduced the main results of the paper such as energy spectrum, wavefunctions, and matrix elements by use of SQcircuit. The diagram of the circuit is

The circuit consists of four equal inductors with \(E_L=2\text{GHz}\), two equal Josephson Junction with \(E_J=15\text{GHz}\) and \(E_{C_J}=2\text{GHz}\), and one shunt capacitor of \(E_{C}=0.04\text{GHz}\).

Circuit Description

Firstly, we import the SQcircuit and the relavant libraries

import SQcircuit as sq
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.colors import LinearSegmentedColormap

We define the single inductive loop of the circuit via Loop class

loop1 = sq.Loop()

The elements of the circuit can be defined via Capacitor, Inductor, and Junction classes in SQcircuit, and to define the circuit, we use the Circuit class.

# define the circuit elements
CJ = sq.Capacitor(2, 'GHz', Q=1e6)
Cshunt = sq.Capacitor(0.04, 'GHz')
L = sq.Inductor(2,'GHz', loops = [loop1])
JJ = sq.Junction(15,'GHz', loops = [loop1])

# define the circuit
elements = {
    (0, 1): [L],
    (0, 2): [L],
    (1, 3): [CJ, JJ],
    (2, 4): [CJ, JJ],
    (3, 5): [L],
    (4 ,5): [L],
    (0, 5): [Cshunt]
}

cr = sq.Circuit(elements)

By creating an object of Circuit class, SQcircuit systematically finds the correct set of transformations and basis to make the circuit ready to be diagonalized.

Before setting the truncation numbers for each mode and diagonalizing the Hamiltonian, we can gain more insight into our circuit by calling the description() method. This prints out the Hamiltonian and a listing of the modes, whether they are harmonic or charge modes, and the frequency for each harmonic in GHz (the default unit). Moreover, it shows the prefactors in the Josephson junction part of the Hamiltonian \(\tilde{\textbf{w}}_k\), which helps find the modes decoupled from the nonlinearity of the circuit.

cr.description()

$\hat{H} =~\omega_1\hat a^\dagger_1\hat a_1~+~\omega_2\hat a^\dagger_2\hat a_2~+~E_{C_{33}}(\hat{n}_3-n_{g_{3}})^2~~-~E_{J_{1}}\cos(\hat{\varphi}_1-\hat{\varphi}_2+\hat{\varphi}_3-0.5\varphi_{\text{ext}_{1}})~-~E_{J_{2}}\cos(\hat{\varphi}_1+\hat{\varphi}_2+\hat{\varphi}_3+0.5\varphi_{\text{ext}_{1}})$

$------------------------------------------------------------$

$\text{mode}~1:~~~~~~~~~~~\text{harmonic}~~~~~~~~~~~\hat{\varphi}_1~=~\varphi_{zp_{1}}(\hat a_1+\hat a^\dagger_1)~~~~~~~~~~~\omega_1/2\pi~=~4.07921~~~~~~~~~~~\varphi_{zp_{1}}~=~9.71e-01$

$\text{mode}~2:~~~~~~~~~~~\text{harmonic}~~~~~~~~~~~\hat{\varphi}_2~=~\varphi_{zp_{2}}(\hat a_2+\hat a^\dagger_2)~~~~~~~~~~~\omega_2/2\pi~=~4.0~~~~~~~~~~~\varphi_{zp_{2}}~=~1.00e+00$

$\text{mode}~3:~~~~~~~~~~~\text{charge}~~~~~~~~~~~~~~~~n_{g_{3}}~=~0$

$------------------------------------------------------------$

$\text{parameters}:~~~~~~~~~~~E_{C_{33}}~=~0.154~~~~~~~~~~~E_{J_{1}}~=~15.0~~~~~~~~~~~E_{J_{2}}~=~15.0~~~~~~~~~~~$

$\text{loops}:~~~~~~~~~~~~~~~~~~~~\varphi_{\text{ext}_{1}}/2\pi~=~0.0~~~~~~~~~~~$

The above output shows that cr circuit has two harmonic modes with one charge mode.

To determine the size of the Hilbert space, we specify the truncation number for each circuit mode via set_trunc_nums() method. Note that this is a necessary step before diagonalizing the circuit.

cr.set_trunc_nums([17,17,9])

Circuit Spectrum

To generate the spectrum of the circuit, firstly, we need to change and sweep the external flux of loop1 by the set_flux() method. Then, we need to find the eigenfrequencies of the circuit that correspond to that external flux via diag() method. The following lines of code find the spec a 2D NumPy array so that each column of it contains the eigenfrequencies with respect to its external flux.

# number of eigenvalues we aim for
n_eig=6

# array that contains the spectrum
phi = np.linspace(-0.1,0.6,50)

# array that contains the spectrum
spec = np.zeros((n_eig, len(phi)))

for i in range(len(phi)):
    # set the value of the flux external flux
    loop1.set_flux(phi[i])

    # diagonlize the circuit
    spec[:, i], _ = cr.diag(n_eig)

E_L = 2
E_C = 2
x=0.02
plt.figure()
for i in range(n_eig):
    plt.plot(phi, (spec[i, :] - spec[0, :])/np.sqrt(16*x*(E_L/2)*E_C))

plt.xlabel("$\Phi_{ext}/\Phi_0$", fontsize=13)
plt.ylabel(r"$f_i-f_0$[GHz]", fontsize=13)
plt.show()

The next cell shows the spectrum from the figure 3a of the paper, which is the same spectrum that SQcircuit calculated.

Eigenfunctions

We can get the phase space eigenfunction of a specific eigenvector of a circuit by using the eig_phase_coord() method. To calculate the eigenfunction at \(\Phi_{ext} = 0.5\Phi_0\) similar to paper, we set back the flux of our loop to \(0.5\) and diagonalize the cr again.

loop1.set_flux(0.5)
cr.set_trunc_nums([17,17,17])
_, _ = cr.diag(n_eig=5)

We generate the eigenfunction in the phase space by eig_phase_coord() method for the first four eigen state of the circuit similar to the paper.

# create a range for each mode
phi1 = 0
phi2 = np.linspace(-0.8,0.8,100)
phi3 = np.pi*np.linspace(-1.5,2.5,100)

# creat the grid list
grid = [phi1, phi2, phi3]

# the ground state
state0 = cr.eig_phase_coord(0, grid = grid)

# the first excited state
state1 = cr.eig_phase_coord(1, grid = grid)

# the second excited state
state2 = cr.eig_phase_coord(2, grid = grid)

# the third excited state
state3 = cr.eig_phase_coord(3, grid = grid)

fig, axs = plt.subplots(1, 4,figsize=(16,4), sharey='row')
axs[0].pcolor(phi3, phi2, np.abs(state0.T)**2,cmap="Purples",shading='auto',label='state0')
axs[1].pcolor(phi3, phi2, np.abs(state1.T)**2,cmap="Greens",shading='auto',label='state1')
axs[2].pcolor(phi3, phi2, np.abs(state2.T)**2,cmap="Oranges",shading='auto',label='state2')
axs[3].pcolor(phi3, phi2, np.abs(state3.T)**2,cmap="Reds",shading='auto',label='state3')
for i in range(4):
    axs[i].set_xlabel(r"$\varphi_2$",fontsize=13)
    axs[i].legend(handletextpad=-0.1, handlelength=0.0)
axs[0].set_ylabel(r"$\varphi_3$",fontsize=13)
plt.subplots_adjust(wspace=0)

The next cell shows the eigenfunctions from the figure 3c of the paper, which is the same as what SQcircuit calculated.

Matrix Elements

Smith2020 defined the normalized matrix element as a measure of degree of freedoms to which circuit can couple. For coupling operator \(\mathcal{O}\) the normalized matrix element is

where \(|\psi\rangle\) is the eigenstates of the circuit and \(|0o\rangle\) is the ground state. We calculate the normalized matrix element for inductive coupling to inductor between node 0 and node 2 using cr.coupling_op() method.

n_eig = 2
phi = np.linspace(0.46, 0.54, 1000)
# normalized matrix element
inductive_01 = np.zeros(len(phi))

for i in range(len(phi)):
#     print(i)
    loop1.set_flux(phi[i])
    evals, estates = cr.diag(n_eig)
    O = cr.coupling_op("inductive", nodes=(2,0))
    norm = (estates[0].dag()*(O.dag()*O)*estates[0])[0].real
    inductive_01[i] =  np.abs((estates[0].dag()*O*estates[1])[0])**2/norm

plt.semilogy(phi, inductive_01, 'green')

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

The next figure shows the same result from the figure 4a of the paper.