Tutorial 5: Generic Junction Potential and Fluxonium EPR#

This tutorial covers two advanced features of the pyEPR numerical diagonalization:

  1. Exact cosine (use_full_cos=True) — avoids truncation errors for strongly anharmonic circuits such as fluxonium (φ_zpf ≳ 1).

  2. Generic nonlinear potential (make_nonlinear_potential) — lets you supply any scalar potential V(φ) expanded around the operating bias point.

Physics background#

The EPR Hamiltonian is split as

\[H = H_{\rm lin} + H_{\rm nl},\]

where \(H_{\rm lin}\) is built from the HFSS linearised eigenfrequencies (which already encode the harmonic Josephson energy \(E_J\varphi^2/2\) via the linearised inductance \(L_j\)). The nonlinear correction must therefore subtract exactly the constant and quadratic parts of the junction potential:

\[H_{\rm nl} = -\sum_j \frac{E_{J,j}}{h}\, {\rm nl}(\hat\varphi_j), \qquad {\rm nl}(\delta\varphi) = \frac{V(\varphi_0+\delta\varphi)-V(\varphi_0)-\tfrac12 V''(\varphi_0)\,\delta\varphi^2}{|V''(\varphi_0)|}\]

For the standard Josephson junction \(V(\varphi) = \cos\varphi\):

\[V''(0) = -1 \implies {\rm nl}(\delta\varphi) = \cos(\delta\varphi) - 1 + \tfrac{\delta\varphi^2}{2} = -\frac{\delta\varphi^4}{4!} + \frac{\delta\varphi^6}{6!} - \cdots\]

The cos_approx function computes the truncated Taylor series of this correction. For weakly anharmonic circuits (transmon, \(\varphi_{\rm zpf} \lesssim 0.3\)) truncation at order \(\varphi^{16}\) is more than sufficient. For strongly anharmonic circuits (fluxonium, \(\varphi_{\rm zpf} \sim 1{-}3\)) the series may converge poorly or diverge — the exact matrix-exponential path is then required.

References

  • Z. K. Minev et al., arXiv:2010.00620 — Energy-participation quantisation of Josephson circuits (pyEPR original paper).

  • A. Petrescu, C. T. Hann, Z. K. Minev et al., arXiv:2411.15039 — EPR analysis for very anharmonic superconducting circuits (fluxonium full-cosine method).

  • G. Zhu, D. G. Ferguson, V. E. Manucharyan, and J. Koch, Phys. Rev. B 87, 024510 (2013) — Circuit QED with fluxonium qubits.

  • N. A. Masluk, I. M. Pop, A. Kamal, Z. K. Minev, M. H. Devoret, Phys. Rev. Lett. 109, 137002 (2012) — Microwave characterisation of Josephson junction arrays.

import numpy as np
import matplotlib.pyplot as plt

from pyEPR.calcs.back_box_numeric import (
    epr_numerical_diagonalization,
    make_nonlinear_potential,
    cos_full_correction,
)

Part 1 — Transmon: truncated vs. exact cosine#

For a transmon-like circuit the zero-point phase fluctuation is small (\(\varphi_{\rm zpf} \approx 0.15\)), so the truncated Taylor series converges quickly. We verify that cos_trunc=16 and use_full_cos=True agree to better than 0.1 %.

# Transmon-like parameters
freqs   = np.array([5.0])     # GHz
Ljs     = np.array([10e-9])   # H
phi_zpf = np.array([[0.15]])  # reduced ZPF

f_trunc, chi_trunc = epr_numerical_diagonalization(
    freqs, Ljs, phi_zpf, cos_trunc=16, fock_trunc=15
)
f_full, chi_full = epr_numerical_diagonalization(
    freqs, Ljs, phi_zpf, fock_trunc=15, use_full_cos=True
)

print(f"Dressed frequency  — truncated: {np.real(f_trunc[0]):.6f} GHz  |  exact: {np.real(f_full[0]):.6f} GHz")
print(f"Anharmonicity χ    — truncated: {np.real(chi_trunc[0,0]):.4f} MHz  |  exact: {np.real(chi_full[0,0]):.4f} MHz")
rel_diff = abs(np.real(f_full[0]) - np.real(f_trunc[0])) / abs(np.real(f_full[0]))
print(f"Relative difference in frequency: {rel_diff:.2e}")

Part 2 — Fluxonium: exact cosine is essential#

For a fluxonium-like circuit \(\varphi_{\rm zpf} \approx 2\). The truncated series at order \(\varphi^8\) gives qualitatively wrong results; the exact cosine (use_full_cos=True) is required. This is the main result of arXiv:2411.15039.

# Fluxonium-like parameters (large phi_zpf → strongly anharmonic)
freqs_fl   = np.array([0.5])     # GHz
Ljs_fl     = np.array([500e-9])  # H  (large superinductance)
phi_zpf_fl = np.array([[2.0]])   # phi_zpf >> 1

# Low truncation — INACCURATE
f_low, chi_low = epr_numerical_diagonalization(
    freqs_fl, Ljs_fl, phi_zpf_fl, cos_trunc=4, fock_trunc=25
)

# Exact cosine — CORRECT
f_exact, chi_exact = epr_numerical_diagonalization(
    freqs_fl, Ljs_fl, phi_zpf_fl, fock_trunc=25, use_full_cos=True
)

print("Fluxonium — low truncation vs. exact cosine:")
print(f"  f (trunc=4):  {np.real(f_low[0]):.4f} GHz")
print(f"  f (exact):    {np.real(f_exact[0]):.4f} GHz")
rel = abs(np.real(f_exact[0]) - np.real(f_low[0])) / abs(np.real(f_exact[0]))
print(f"  Relative error from low truncation: {rel:.1%}  ← significant!")
print()
print(f"  χ (trunc=4): {np.real(chi_low[0,0]):.2f} MHz")
print(f"  χ (exact):   {np.real(chi_exact[0,0]):.2f} MHz")

Part 3 — Generic potential: expanding around the bias minimum#

make_nonlinear_potential(V, phi_min) converts any scalar Python function \(V(\varphi)\) into the operator-valued EPR correction. The function:

  1. Evaluates \(V_0 = V(\varphi_0)\) and \(V'' = V''(\varphi_0)\) (numerically).

  2. Returns a callable nl(φ_op) that computes $\(\mathrm{nl}(\delta\varphi) = \frac{V(\varphi_0+\delta\varphi) - V_0 - \tfrac{V''}{2}\delta\varphi^2}{|V''|}\)$ as a qutip operator via eigendecomposition of the phase operator.

When to use this: when the junction potential is not a simple cosine centred at zero — e.g., a flux-biased junction, an asymmetric SQUID, or any exotic Josephson element.

3a. Standard junction — recovering the default#

Passing V = np.cos recovers cos_full_correction to numerical precision.

nl_from_cos = make_nonlinear_potential(np.cos)   # V(phi) = cos(phi), phi_min=0

f_gen, chi_gen = epr_numerical_diagonalization(
    freqs, Ljs, phi_zpf, fock_trunc=15,
    non_linear_potential=nl_from_cos,
)
print("Standard JJ via make_nonlinear_potential vs. use_full_cos=True:")
print(f"  f  — generic: {np.real(f_gen[0]):.6f} GHz  |  exact: {np.real(f_full[0]):.6f} GHz")
print(f"  χ  — generic: {np.real(chi_gen[0,0]):.4f} MHz  |  exact: {np.real(chi_full[0,0]):.4f} MHz")

3b. Flux-biased junction at \(\varphi_{\rm ext}\)#

A Josephson junction biased by external flux \(\varphi_{\rm ext}\) has potential

\[U(\varphi) = -E_J \cos(\varphi - \varphi_{\rm ext})\]

The minimum is at \(\varphi_0 = \varphi_{\rm ext}\). After expansion around the minimum, the potential in the fluctuation variable \(\delta\varphi = \varphi - \varphi_0\) is simply \(-E_J^{\rm eff}\cos(\delta\varphi)\), where the effective Josephson energy \(E_J^{\rm eff} = E_J\cos(\varphi_{\rm ext})\) is encoded in the linearised inductance \(L_j\) from the HFSS EPR run.

This means the nonlinear correction is identical to the unbiased case — make_nonlinear_potential confirms this automatically.

phi_ext = np.pi / 6   # 30 degrees external flux

nl_biased = make_nonlinear_potential(
    lambda phi: np.cos(phi - phi_ext),
    phi_min=phi_ext,
)

f_bias, chi_bias = epr_numerical_diagonalization(
    freqs, Ljs, phi_zpf, fock_trunc=15,
    non_linear_potential=nl_biased,
)

print(f"Flux-biased junction at phi_ext = pi/6:")
print(f"  f  = {np.real(f_bias[0]):.6f} GHz")
print(f"  χ  = {np.real(chi_bias[0,0]):.4f} MHz")
print()
print("(Same as unbiased since Ej_eff is encoded in Lj passed to the EPR solver.")
print(" Only use a different Lj = (Φ₀/2π)² / Ej_eff for the biased case.)")

3c. Asymmetric SQUID#

For an asymmetric SQUID (two junctions with energies \(E_{J1}\), \(E_{J2}\)) at external flux \(\varphi_{\rm ext}\), the effective potential seen by the common mode phase is

\[V_{\rm SQUID}(\varphi) = \cos\varphi\cos\varphi_{\rm ext} + d\sin\varphi\sin\varphi_{\rm ext}\]

where \(d = (E_{J1}-E_{J2})/(E_{J1}+E_{J2})\) is the asymmetry parameter and the result is normalised by \(E_{J1}+E_{J2}\) (see arXiv:2010.00620, §IV). The potential minimum is at

\[\varphi_0 = \arctan(-d\tan\varphi_{\rm ext}).\]

make_nonlinear_potential handles the expansion around \(\varphi_0\) and normalisation automatically.

phi_ext = 0.4   # external flux (reduced units)
d       = 0.15  # asymmetry parameter

def V_squid(phi):
    """Asymmetric SQUID normalised potential."""
    return np.cos(phi) * np.cos(phi_ext) + d * np.sin(phi) * np.sin(phi_ext)

# Minimum of -Ej * V_squid (maximum of V_squid)
phi_min_squid = np.arctan(-d * np.tan(phi_ext)) % np.pi

print(f"SQUID bias phi_ext = {phi_ext:.2f} rad,  asymmetry d = {d}")
print(f"Potential minimum at phi_min = {phi_min_squid:.4f} rad")
print(f"V''(phi_min) = {-(np.cos(phi_min_squid)*np.cos(phi_ext) + d*np.sin(phi_min_squid)*np.sin(phi_ext)):.4f}")

nl_squid = make_nonlinear_potential(V_squid, phi_min=phi_min_squid)

f_sq, chi_sq = epr_numerical_diagonalization(
    freqs, Ljs, phi_zpf, fock_trunc=15,
    non_linear_potential=nl_squid,
)
print(f"\nAsymmetric SQUID:  f = {np.real(f_sq[0]):.4f} GHz,  χ = {np.real(chi_sq[0,0]):.4f} MHz")

3d. Visualising the potential and its nonlinear correction#

Understanding what make_nonlinear_potential subtracts is key to using it correctly.

phi = np.linspace(-np.pi, np.pi, 500)

# ── Standard JJ ─────────────────────────────────────────────────────────────
V = np.cos(phi)                        # full potential (normalised)
V0 = np.cos(0)                         # = 1
Vpp = -np.cos(0)                       # = -1  (V'' at the minimum)
V_harm = V0 + Vpp / 2 * phi**2         # quadratic approx (already in H_lin)
nl = (V - V0 - Vpp / 2 * phi**2) / abs(Vpp)  # nonlinear correction

fig, axes = plt.subplots(1, 2, figsize=(11, 4))

ax = axes[0]
ax.plot(phi, V,      label=r'$V(\varphi) = \cos\varphi$', lw=2)
ax.plot(phi, V_harm, label=r'Harmonic approx $V_0 + V\'\'/2\,\varphi^2$', lw=2, ls='--')
ax.set_xlabel(r'$\varphi$ (rad)')
ax.set_ylabel(r'$V / E_J$')
ax.set_title('Junction potential (standard JJ)')
ax.legend()
ax.set_xlim(-np.pi, np.pi)

ax = axes[1]
ax.plot(phi, nl, lw=2, color='C2',
        label=r'$\mathrm{nl}(\varphi) = \cos\varphi - 1 + \varphi^2/2$')
ax.axhline(0, color='k', lw=0.5)
ax.set_xlabel(r'$\varphi$ (rad)')
ax.set_ylabel(r'$\mathrm{nl}(\varphi)$')
ax.set_title('Nonlinear correction (what EPR adds)')
ax.legend()
ax.set_xlim(-np.pi, np.pi)

plt.tight_layout()
plt.show()
print("The harmonic approximation (dashed) is already captured by H_lin.")
print("The EPR nonlinear correction nl(φ) adds everything beyond the quadratic.")

Part 4 — Convergence comparison: transmon vs. fluxonium#

We sweep the Fock-space truncation and compare the dressed frequency as a function of truncation order for both circuit types.

fock_values = [8, 10, 12, 15, 20, 25]

f_transmon_full  = []
f_fluxonium_full = []

for ft in fock_values:
    f_t, _ = epr_numerical_diagonalization(freqs,    Ljs,    phi_zpf,    fock_trunc=ft, use_full_cos=True)
    f_f, _ = epr_numerical_diagonalization(freqs_fl, Ljs_fl, phi_zpf_fl, fock_trunc=ft, use_full_cos=True)
    f_transmon_full.append(np.real(f_t[0]))
    f_fluxonium_full.append(np.real(f_f[0]))

fig, axes = plt.subplots(1, 2, figsize=(11, 4))

for ax, label, values, ref in [
    (axes[0], 'Transmon (φ_zpf=0.15)', f_transmon_full,  f_transmon_full[-1]),
    (axes[1], 'Fluxonium (φ_zpf=2.0)', f_fluxonium_full, f_fluxonium_full[-1]),
]:
    rel_err = [abs(f - ref) / abs(ref) * 100 for f in values]
    ax.semilogy(fock_values, rel_err, 'o-', lw=2)
    ax.set_xlabel('Fock truncation')
    ax.set_ylabel('Relative error vs. fock_trunc=25 (%)')
    ax.set_title(label)
    ax.grid(True, which='both', alpha=0.3)

plt.tight_layout()
plt.show()

Summary and API cheat-sheet#

Scenario

Recommended call

Transmon (φ_zpf ≲ 0.3)

epr_numerical_diagonalization(..., cos_trunc=8) (default)

Fluxonium / strongly anharmonic (φ_zpf ≳ 1)

epr_numerical_diagonalization(..., use_full_cos=True)

Flux-biased JJ at φ_ext

make_nonlinear_potential(lambda φ: cos(φ-φ_ext), phi_min=φ_ext)

Asymmetric SQUID

make_nonlinear_potential(V_squid, phi_min=phi_min_squid)

Arbitrary junction potential V(φ)

make_nonlinear_potential(V, phi_min=phi_min_V)

The non_linear_potential argument of epr_numerical_diagonalization and black_box_hamiltonian accepts any callable nl(φ_op: Qobj) -> Qobj that returns the correction beyond the quadratic, normalised by \(|V''(\varphi_0)|\). Use make_nonlinear_potential to construct this automatically from a scalar Python function.

References#