Issue computing "inductively coupled Transmon" circuit

Question

Describe the Issue!

Dear all,

I have trouble using qutip to compute the eigenenergies of an "inductively coupled transmon" circuit:

First I tried to compute the eigenenergies of a Transmon using from the Koch et al. paper:

With:

• the charging energy
• the Josephson energy

From that I wrote the following function:

def spectrum(fluxes: np.ndarray,
             dim: int=11,
             c: float=10e-13,
             ic: float=10e-9,
             ng: float=0):

    n = qtp.charge(int(dim/2))
    phi = qtp.phase(dim)

    # Charging energy in GHz
    ec = (2*cst.e)**2/2/c/cst.h/1e9

    # Josephson energy in GHz
    ej = cst.phi_0/2/np.pi*ic/cst.h/1e9

    # Where we store the eigenenergies
    es = np.zeros((len(fluxes), dim-1))
    for i, flux in enumerate(fluxes):

        h = ec*(n + ng)**2 - ej*(phi.cosm() + (phi - 2*np.pi*flux).cosm())

        e = h.eigenenergies()
        e = (e - e[0])
        es[i] = np.abs(e[1:])

    return es

And got the correct Transmon transition frequency:

fluxes = np.linspace(-0.5, 1, 100)

fig, ax = plt.subplots(1, 1)
ys = spectrum(fluxes, ng=0., dim=21)
ax.plot(fluxes, ys[:,0])

ax.set_xlabel('fluxes')
ax.set_ylabel('transitions (GHz)')

We have the expected cos shape.

Then I went for the "inductively coupled transmon" circuit with the following hamiltonian

With:

• the charging energy
• the Josephson energy
• and
• and
• the inductance ratio

From that I wrote the following function:

def spectrumInductive(fluxes: np.ndarray,
             dim: int=11,
             c: float=10e-13,
             ic: float=10e-9,
             l: float=1000e-12,
             ng: float=0):

    n1 = qtp.tensor(qtp.charge(int(dim/2))+ng, qtp.qeye(dim))
    n2 = qtp.tensor(qtp.qeye(dim), qtp.charge(int(dim/2))+ng)

    phi1 = qtp.tensor(qtp.phase(dim), qtp.qeye(dim))
    phi2 = qtp.tensor(qtp.qeye(dim), qtp.phase(dim))

    nx = (n1+n2)/2
    ny = (n1-n2)/2

    phix = (phi1+phi2)/2
    phiy = (phi1-phi2)/2

    # Charging energy in GHz
    ec = (2*cst.e)**2/2/c/cst.h/1e9/2

    # Josephson energy in GHz
    ej = cst.phi_0/2/np.pi*ic/2/cst.h/1e9/2

    # The inductance ratio
    b = cst.phi_0/2/np.pi/ic/l

    h1 = ec*(nx**2 + ny**2)
    h2 = ej*(-phix.cosm()*phiy.cosm())

    es = np.zeros((len(fluxes), dim*dim-1))
    for i, flux in enumerate(fluxes):

        h = h1 + h2 + ej*b*(phiy - np.pi*flux)**2

        e = h.eigenenergies()
        e = (e - e[0])
        es[i] = np.abs(e[1:])

    return es

And got the incorrect "inductively coupled Transmon" transition frequency:

fluxes = np.linspace(-0.5, 0.5, 20)

fig, ax = plt.subplots(1, 1)
ys1 = spectrumInductive(fluxes, c=40e-15, ic=8e-9, l=7.5e-9, dim=21, ng=0.5)
ax.plot(fluxes, ys1[:,0])
ax.set_xlabel('fluxes')
ax.set_ylabel('transitions (GHz)')

We do not have the expected cos shape.

If anyone as any tips I would be very grateful!