Transmon Qubit Design with NetKet#

Required extras

This notebook needs the models and netket extras:

uv add "qpdk[models,netket]"
# or, from a checkout of this repository:
uv sync --extra models --extra netket
# or with pip:
pip install "qpdk[models,netket]"

See the extras reference for what each extra installs.

This notebook uses NetKet [CCH+19, VHS+22] to analyse the transmon qubit Hamiltonian numerically for superconducting qubit design, and translates the resulting quantum parameters into physical layout parameters using qpdk.

NetKet is a JAX-based library for many-body quantum physics that provides exact diagonalisation and variational solvers for custom Hamiltonians. For superconducting qubit design, it computes energy spectra, charge dispersion, and coupling parameters directly from the microscopic circuit Hamiltonian.

Each section takes its circuit parameters from a qpdk layout cell or an existing qpdk model, and hands NetKet a Hamiltonian built from them. Where qpdk has no model for a parameter yet, we say so and use a design choice or the model’s default:

  • a double-pad transmon, with \(E_\text{C}\) from the Palace capacitance lookup of its pads and \(E_\text{J}\) as a design choice, solved by exact diagonalisation, by a hand-written variational loop, and by NetKet’s own variational Monte Carlo (VMC) driver;

  • the transmon–resonator system, including its real-time dynamics with NetKet’s experimental time-dependent variational principle (TDVP);

  • a Josephson junction array, a many-body problem whose Hilbert space grows exponentially with its length, solved with Monte Carlo VMC and a custom Flax network, and fed back into the array’s SAX model;

  • the unimon and gatemon, whose existing qpdk Hamiltonians become NetKet operators with one call; the gatemon takes \(E_\text{C}\) from its cell, and the unimon uses its model’s default energies.

Background#

The transmon qubit [KYG+07] is an anharmonic oscillator formed by shunting a Josephson junction with a large capacitance.

THE TWO ENERGY SCALES OF A TRANSMON A Josephson junction shunted by capacitance forms a weakly anharmonic circuit. 𝐶 Σ junction 𝐶 Σ sets 𝐸 C = 𝑒 2 2 𝐶 Σ 𝐼 c sets 𝐸 J = Φ 0 𝐼 c 2 𝜋

The shunt capacitance sets \(E_\text{C}\), while the junction critical current sets \(E_\text{J}\). The transmon Hamiltonian in the charge basis reads:

(1)#\[\hat{H} = 4 E_\text{C} (\hat{n} - n_\text{g})^2 - E_\text{J} \cos\hat{\varphi},\]

where \(E_\text{C}\) is the charging energy, \(E_\text{J}\) the Josephson energy, \(\hat{n}\) the Cooper-pair number operator, \(n_\text{g}\) the offset charge, and \(\hat{\varphi}\) the superconducting phase.

In the truncated charge basis \(\{|n\rangle\}\) with \(n \in \{-n_\text{cut}, \ldots, +n_\text{cut}\}\):

  • Diagonal elements: \(\langle n | \hat{H} | n \rangle = 4 E_\text{C} (n - n_\text{g})^2\)

  • Off-diagonal elements: \(\langle n \pm 1 | \hat{H} | n \rangle = -E_\text{J} / 2\)

In the transmon regime (\(E_\text{J} / E_\text{C} \gg 1\)), the qubit frequency and anharmonicity are approximately [KYG+07]:

(2)#\[\omega_\text{q} \approx \sqrt{8 E_\text{J} E_\text{C}} - E_\text{C}, \qquad \alpha \approx -E_\text{C}.\]

Hide code cell source

 
import os
import sys

if "google.colab" in sys.modules:
    import subprocess

    print("Running in Google Colab. Installing QPDK...")
    subprocess.check_call([
        sys.executable,
        "-m",
        "pip",
        "install",
        "-q",
        "qpdk[models,netket] @ git+https://github.com/gdsfactory/quantum-rf-pdk.git",
    ])
    # On a Colab TPU runtime, installing qpdk replaces the TPU build of JAX
    # with the CPU one, so put the TPU runtime library back at the same JAX
    # version.
    if os.environ.get("TPU_ACCELERATOR_TYPE") or os.environ.get("COLAB_TPU_ADDR"):
        from importlib.metadata import version

        subprocess.check_call([
            sys.executable,
            "-m",
            "pip",
            "install",
            "-q",
            f"jax[tpu]=={version('jax')}",
        ])

Hide code cell source

 
import time
from dataclasses import dataclass
from functools import partial
from typing import Any

import flax.linen as nn
import jax
import jax.numpy as jnp
import netket as nk
import netket.experimental as nkx
import numpy as np
import optax
import polars as pl
import scipy
from IPython.display import Math, display
from matplotlib import pyplot as plt
from tqdm.auto import tqdm

from qpdk import PDK
from qpdk.helper import display_dataframe
from qpdk.models.constants import Φ_0, c_0, e, h
from qpdk.models.cpw import cpw_parameters

PDK.activate()

from qpdk.cells import gatemon, josephson_junction_array
from qpdk.cells.transmon import (
    double_pad_transmon,
    xmon_transmon,
)
from qpdk.models.gatemon import (
    gatemon_capacitances,
    gatemon_frequency_and_anharmonicity,
    gatemon_hamiltonian,
)
from qpdk.models.junction import (
    josephson_junction_array as josephson_junction_array_model,
)
from qpdk.models.perturbation import (
    dispersive_shift_to_coupling,
    ej_ec_to_frequency_and_anharmonicity,
    measurement_induced_dephasing,
    purcell_decay_rate,
    resonator_linewidth_from_q,
)
from qpdk.models.qubit import (
    coupling_strength_to_capacitance,
    double_pad_transmon_capacitance,
    ec_to_capacitance,
    ej_to_inductance,
)
from qpdk.models.resonator import resonator_frequency
from qpdk.models.unimon import unimon_frequency_and_anharmonicity, unimon_hamiltonian

Hide code cell output

 
∣NK⟩ Tip: With many Markov Chains (e.g GPUs), n_discard_per_chain>5 is often inefficient.

Building the Transmon Hamiltonian#

The charging energy comes from the layout. We read the pad size and gap of a double_pad_transmon() cell and look up the capacitance between its pads with double_pad_transmon_capacitance(), which interpolates a stored Palace electrostatic sweep. The Josephson energy is set by the junction, which we treat as a design choice here.

We then construct the transmon Hamiltonian as a matrix in the truncated charge basis and wrap it in a NetKet LocalOperator, ready for NetKet’s exact-diagonalisation routines.

transmon_cell = double_pad_transmon()
pad_width, pad_height = transmon_cell.settings["pad_size"]
C_pads = float(
    double_pad_transmon_capacitance(
        pad_width=pad_width,
        pad_height=pad_height,
        pad_gap=transmon_cell.settings["pad_gap"],
    )
)
if not np.isfinite(C_pads):
    msg = f"Pad size {transmon_cell.settings['pad_size']} is outside the capacitance lookup grid"
    raise ValueError(msg)
# The lookup excludes the junction; add its capacitance for your process here.
C_junction = 0.0

EJ = 20.0  # Josephson energy in GHz, set by the junction
EC = e**2 / (2 * (C_pads + C_junction)) / h / 1e9  # Charging energy in GHz
ng = 0.0  # Offset charge (gate charge)

display(
    Math(rf"""
C_\text{{pads}} = {C_pads * 1e15:.1f}\,\text{{fF}} \quad\Rightarrow\quad
E_\text{{C}} = e^2 / 2C = {EC:.4f}\,\text{{GHz}}
""")
)

ncut = 15  # Charge basis truncation: n ∈ {-ncut, ..., +ncut}
n_states = 2 * ncut + 1

# Build the Hamiltonian matrix in the charge basis
charges = jnp.arange(n_states) - ncut  # physical charge quantum numbers
H_mat = jnp.diag(4 * EC * (charges - ng) ** 2)
for n in range(n_states - 1):
    H_mat = H_mat.at[n, n + 1].set(-EJ / 2)  # Josephson tunnelling
    H_mat = H_mat.at[n + 1, n].set(-EJ / 2)

# Wrap as a NetKet operator on a Fock Hilbert space
hi = nk.hilbert.Fock(n_max=n_states - 1, N=1)
H_transmon = nk.operator.LocalOperator(hi, np.asarray(H_mat), acting_on=[0])
\[\displaystyle C_\text{pads} = 89.0\,\text{fF} \quad\Rightarrow\quad E_\text{C} = e^2 / 2C = 0.2177\,\text{GHz} \]

Transmon Spectrum via Exact Diagonalisation#

NetKet provides lanczos_ed for efficient sparse-matrix diagonalisation. We extract the lowest five eigenvalues and compute the qubit frequency \(\omega_\text{q} = E_1 - E_0\) and the anharmonicity \(\alpha = (E_2 - E_1) - (E_1 - E_0)\).

evals_exact = nk.exact.lanczos_ed(H_transmon, k=5, compute_eigenvectors=False)
e0_exact = evals_exact[0]
evals = evals_exact - e0_exact  # shift ground state to zero

f01 = evals[1] - evals[0]
f12 = evals[2] - evals[1]
alpha_nk = f12 - f01  # negative for transmon

# Analytical transmon approximation from qpdk
omega_t_approx, _ = ej_ec_to_frequency_and_anharmonicity(EJ, EC)

display(
    Math(rf"""
\textbf{{Transmon Spectrum (NetKet):}} \\
E_\text{{J}}/E_\text{{C}} = {EJ / EC:.0f} \\
f_{{01}} = {f01:.4f}\,\text{{GHz}} \quad
(\text{{approx.}}\;\sqrt{{8 E_\text{{J}} E_\text{{C}}}} - E_\text{{C}} = {omega_t_approx:.4f}\,\text{{GHz}}) \\
f_{{12}} = {f12:.4f}\,\text{{GHz}} \\
\alpha = f_{{12}} - f_{{01}} = {alpha_nk:.4f}\,\text{{GHz}} \quad
(\text{{approx.}}\;-E_\text{{C}} = -{EC:.4f}\,\text{{GHz}})
""")
)
\[\begin{split}\displaystyle \textbf{Transmon Spectrum (NetKet):} \\ E_\text{J}/E_\text{C} = 92 \\ f_{01} = 5.6751\,\text{GHz} \quad (\text{approx.}\;\sqrt{8 E_\text{J} E_\text{C}} - E_\text{C} = 5.6840\,\text{GHz}) \\ f_{12} = 5.4355\,\text{GHz} \\ \alpha = f_{12} - f_{01} = -0.2396\,\text{GHz} \quad (\text{approx.}\;-E_\text{C} = -0.2177\,\text{GHz}) \end{split}\]

Variational Energy via VMC (Full Summation)#

While exact diagonalisation works for this small Hilbert space, larger systems require variational methods. Here we demonstrate a simple Variational Monte Carlo (VMC) approach using JAX and Flax. Here the energy is summed over all basis states rather than sampled.

We define a feedforward network in Flax as the variational ansatz for the ground-state wavefunction. The network takes the charge configuration as input and outputs the logarithm of the wavefunction amplitude (log-psi). We then compute the energy expectation value by summing over all charge states, and use JAX’s automatic differentiation to compute gradients for optimisation.

The network is a compact multilayer perceptron [Ros58] preceded by a sinusoidal input encoding at geometrically spaced frequencies. This encoding lifts the one-dimensional Fock-state index into a higher-dimensional feature space of sines and cosines, enabling the subsequent \(\tanh\)-activated hidden layers to approximate the smooth, Gaussian-like ground-state wavefunction much more efficiently than piecewise-linear activations such as Leaky ReLU.

For details on VMC with NetKet, see the documentation: https://netket.readthedocs.io/en/latest/vmc-from-scratch/index.html

class VariationalTransmon(nn.Module):
    r"""A compact MLP with sinusoidal input encoding for the transmon wavefunction.

    Sinusoidal (Fourier) features at geometrically spaced frequencies lift
    the one-dimensional Fock-state index into a richer feature space, allowing
    the two :math:`\tanh`-activated hidden layers to approximate the smooth,
    Gaussian-like ground state much more efficiently than a plain MLP with
    piecewise-linear activations.
    """

    n_max: int = 30
    n_freqs: int = 4

    @nn.compact
    def __call__(self, x):
        """Evaluate the model on a set of input configurations.

        Returns:
            The model evaluation result as a squeezed array.
        """
        # x has shape (..., 1)
        x = x.astype(jnp.float32)
        # Sinusoidal feature encoding at geometrically spaced frequencies;
        # maps the scalar input to a 2*n_freqs-dimensional feature vector.
        freqs = 2.0 ** jnp.arange(self.n_freqs) * (jnp.pi / self.n_max)
        x = jnp.concatenate([jnp.sin(x * freqs), jnp.cos(x * freqs)], axis=-1)
        x = nn.Dense(features=32)(x)
        x = jnp.tanh(x)
        x = nn.Dense(features=32)(x)
        x = jnp.tanh(x)
        x = nn.Dense(features=1)(x)
        return jnp.squeeze(x, axis=-1)


model = VariationalTransmon(n_max=n_states - 1)
parameters = model.init(jax.random.key(0), jnp.ones((1, 1)))


def to_array(model, parameters):
    """Compute the normalized wavefunction as an array.

    Returns:
        The normalized wavefunction as a JAX array.
    """
    all_configs = hi.all_states()
    logpsi = model.apply(parameters, all_configs)
    psi = jnp.exp(logpsi)
    return psi / jnp.linalg.norm(psi)


def compute_energy(model, parameters, H_sparse):
    """Compute the variational energy of the state.

    Returns:
        The expected value of the energy (real part).
    """
    psi = to_array(model, parameters)
    psi = psi.astype(H_sparse.dtype)
    return jnp.real(jnp.vdot(psi, H_sparse @ psi))


@partial(jax.jit, static_argnames="model")
def train_step(model, parameters, opt_state, H_sparse):
    """Perform a single variational optimization step.

    Returns:
        A tuple containing (updated parameters, updated opt_state, energy).
    """
    energy, grad = jax.value_and_grad(compute_energy, argnums=1)(
        model, parameters, H_sparse
    )
    updates, opt_state = tx.update(grad, opt_state, parameters)
    parameters = optax.apply_updates(parameters, updates)
    return parameters, opt_state, energy


H_sparse = H_transmon.to_sparse(jax_=True)
tx = optax.adam(0.01)
opt_state = tx.init(parameters)

energies = []
for _ in tqdm(range(200), desc="Variational Optimization"):
    parameters, opt_state, e_val = train_step(model, parameters, opt_state, H_sparse)
    energies.append(e_val)

e_gs_variational = float(energies[-1])

display(
    Math(rf"""
\textbf{{Ground State Energy Comparison:}} \\
E_{{0, \text{{exact}}}} = {e0_exact:.6f}\,\text{{GHz}} \\
E_{{0, \text{{variational}}}} = {e_gs_variational:.6f}\,\text{{GHz}} \\
\text{{Error:}} \quad |E_\text{{exact}} - E_\text{{var}}| = {abs(e0_exact - e_gs_variational):.2e}\,\text{{GHz}}
""")
)
\[\begin{split}\displaystyle \textbf{Ground State Energy Comparison:} \\ E_{0, \text{exact}} = -17.104643\,\text{GHz} \\ E_{0, \text{variational}} = -17.103051\,\text{GHz} \\ \text{Error:} \quad |E_\text{exact} - E_\text{var}| = 1.59e-03\,\text{GHz} \end{split}\]

The same network with NetKet’s VMC driver#

The loop above is what NetKet automates. A variational state wraps any Flax module, so the same VariationalTransmon network plugs straight into nk.vqs.FullSumState, which still sums over every basis state. nk.driver.VMC then runs the optimisation and logs the energy. For larger systems the only change is to swap the full-summation state for a Monte Carlo MCState with a sampler, as in the junction-array section below.

vs_transmon = nk.vqs.FullSumState(hi, VariationalTransmon(n_max=n_states - 1), seed=0)
vmc_driver = nk.driver.VMC(H_transmon, optax.adam(0.01), variational_state=vs_transmon)
vmc_log = nk.logging.RuntimeLog()
vmc_driver.run(n_iter=200, out=vmc_log, show_progress=False)
e_gs_netket = float(vs_transmon.expect(H_transmon).mean.real)

fig, ax = plt.subplots(figsize=(8, 4))
ax.semilogy(
    np.abs(np.asarray(energies) - e0_exact), label="Hand-written loop (Flax + optax)"
)
ax.semilogy(
    np.abs(np.asarray(vmc_log["Energy"]["Mean"].real) - e0_exact),
    "--",
    label="NetKet VMC driver",
)
ax.set_xlabel("Iteration")
ax.set_ylabel(r"$|E - E_{0,\text{exact}}|$ (GHz)")
ax.set_title("Same Flax ansatz, two optimisation loops")
ax.legend()
ax.grid(True, which="both", alpha=0.3)
fig.tight_layout()
plt.show()

display(
    Math(rf"""
E_{{0, \text{{NetKet VMC}}}} = {e_gs_netket:.6f}\,\text{{GHz}} \quad
|E_\text{{exact}} - E_\text{{NetKet VMC}}| = {abs(e0_exact - e_gs_netket):.2e}\,\text{{GHz}}
""")
)
../_images/ed0d870f649ecce07bb85ece0ac94b507a22513e4c419a9cf9fff70f2b05f361.svg ../_images/ed0d870f649ecce07bb85ece0ac94b507a22513e4c419a9cf9fff70f2b05f361_dark.svg
\[\displaystyle E_{0, \text{NetKet VMC}} = -17.103062\,\text{GHz} \quad |E_\text{exact} - E_\text{NetKet VMC}| = 1.58e-03\,\text{GHz} \]

Charge Dispersion#

A key design criterion for the transmon is its exponential insensitivity to charge noise. As \(E_\text{J}/E_\text{C}\) increases, the energy levels flatten as a function of the offset charge \(n_\text{g}\), suppressing dephasing from charge fluctuations [KYG+07].

We sweep \(n_\text{g}\) across one period and compute the first three transition energies for several values of \(E_\text{J}/E_\text{C}\).

ng_sweep = jnp.linspace(-0.5, 0.5, 101)
ej_ec_ratios = [1, 5, 20, 50]

fig, axes = plt.subplots(1, len(ej_ec_ratios), figsize=(14, 3.5), sharey=True)

for ax, ratio in zip(axes, ej_ec_ratios):
    EJ_sweep = ratio * EC
    levels = jnp.zeros((len(ng_sweep), 4))

    for i, ng_val in enumerate(ng_sweep):
        H_sweep = jnp.diag(4 * EC * (charges - ng_val) ** 2)
        for n in range(n_states - 1):
            H_sweep = H_sweep.at[n, n + 1].set(-EJ_sweep / 2)
            H_sweep = H_sweep.at[n + 1, n].set(-EJ_sweep / 2)
        ev = jnp.linalg.eigvalsh(H_sweep)
        levels = levels.at[i].set(ev[:4] - ev[0])

    for j in range(1, 4):
        ax.plot(ng_sweep, levels[:, j], label=rf"$E_{j} - E_0$")
    ax.set_xlabel(r"$n_\text{g}$")
    ax.set_title(rf"$E_J/E_C = {ratio}$")
    ax.grid(True, alpha=0.3)

axes[0].set_ylabel("Energy (GHz)")
axes[-1].legend(fontsize=8)
fig.suptitle("Charge Dispersion: Transmon Regime", y=1.02)
fig.tight_layout()
plt.show()
../_images/9ef312f552013bf464f9de126c0c1fa875171e8e86405e88a7b17571743c82b3.svg ../_images/9ef312f552013bf464f9de126c0c1fa875171e8e86405e88a7b17571743c82b3_dark.svg

The energy levels flatten with respect to the offset charge \(n_\text{g}\) as the ratio \(E_\text{J}/E_\text{C}\) grows. At \(E_\text{J}/E_\text{C} \gtrsim 50\), well inside the transmon regime, the charge dispersion is exponentially suppressed [KYG+07].

Transmon–Resonator Coupling#

In circuit QED, the transmon is coupled to a readout resonator via the electric-field interaction [BGGW21]:

(3)#\[\hat{H}_\text{int} = g\, (\hat{a} + \hat{a}^\dagger)\, (\hat{b} + \hat{b}^\dagger),\]

where \(\hat{a}\) and \(\hat{b}\) are the transmon and resonator ladder operators respectively, and \(g\) is the coupling strength.

We build the full transmon–resonator Hamiltonian in the energy eigenbasis of the transmon (Duffing oscillator approximation) and diagonalise it using NetKet to extract the dispersive shift \(\chi\):

(4)#\[\chi = (E_{11} - E_{10}) - (E_{01} - E_{00}),\]

where \(E_{ij}\) is the dressed eigenenergy with \(i\) transmon excitations and \(j\) resonator photons.

omega_r_val = 7.0  # Resonator frequency in GHz
g_val = 0.1  # Coupling strength in GHz
n_levels = 8  # Truncation level for each mode

# Combined Hilbert space: site 0 = transmon, site 1 = resonator
hi_combined = nk.hilbert.Fock(n_max=n_levels - 1, N=2)

# Transmon as Duffing oscillator: f01 * n + (α/2) * n(n-1)
H_qubit = jnp.diag(
    jnp.array([f01 * n + (alpha_nk / 2) * n * (n - 1) for n in range(n_levels)])
)

# Resonator: ω_r * n
H_res = jnp.diag(jnp.array([omega_r_val * n for n in range(n_levels)]))

# (a + a†) ladder operator
a_plus_adag = jnp.zeros((n_levels, n_levels))
for n in range(n_levels - 1):
    a_plus_adag = a_plus_adag.at[n, n + 1].set(jnp.sqrt(n + 1))
    a_plus_adag = a_plus_adag.at[n + 1, n].set(jnp.sqrt(n + 1))

# Coupling operator: g * (a + a†) ⊗ (b + b†)
coupling = g_val * jnp.kron(a_plus_adag, a_plus_adag)

# Build full Hamiltonian as sum of NetKet operators
H_full = (
    nk.operator.LocalOperator(hi_combined, np.asarray(H_qubit), acting_on=[0])
    + nk.operator.LocalOperator(hi_combined, np.asarray(H_res), acting_on=[1])
    + nk.operator.LocalOperator(hi_combined, np.asarray(coupling), acting_on=[0, 1])
)

# Full exact diagonalisation
evals_full, evecs_full = nk.exact.full_ed(H_full, compute_eigenvectors=True)
sort_idx = jnp.argsort(evals_full)
evals_full = evals_full[sort_idx]
evecs_full = evecs_full[:, sort_idx]


# Identify dressed states by maximum overlap with bare states
def _find_dressed(qubit_idx: int, res_idx: int) -> int:
    """Return the index of the dressed state closest to |qubit_idx, res_idx⟩."""
    bare = jnp.zeros(n_levels**2)
    bare = bare.at[qubit_idx * n_levels + res_idx].set(1.0)
    return int(jnp.argmax(jnp.abs(evecs_full.T @ bare) ** 2))


idx_00 = _find_dressed(0, 0)
idx_10 = _find_dressed(1, 0)
idx_01 = _find_dressed(0, 1)
idx_11 = _find_dressed(1, 1)

chi_netket = (evals_full[idx_11] - evals_full[idx_10]) - (
    evals_full[idx_01] - evals_full[idx_00]
)

display(
    Math(rf"""
\textbf{{Dispersive Shift (NetKet):}} \\
g = {g_val * 1e3:.0f}\,\text{{MHz}}, \quad
\omega_\text{{r}} = {omega_r_val:.1f}\,\text{{GHz}} \\
\Delta = \omega_\text{{q}} - \omega_\text{{r}} = {f01 - omega_r_val:.3f}\,\text{{GHz}} \\
\chi = {chi_netket * 1e3:.3f}\,\text{{MHz}}
""")
)
\[\begin{split}\displaystyle \textbf{Dispersive Shift (NetKet):} \\ g = 100\,\text{MHz}, \quad \omega_\text{r} = 7.0\,\text{GHz} \\ \Delta = \omega_\text{q} - \omega_\text{r} = -1.325\,\text{GHz} \\ \chi = -2.298\,\text{MHz} \end{split}\]

Dispersive Shift vs. Coupling Strength#

We sweep the coupling strength \(g\) and compare the dispersive shift with its quadratic scaling, which holds only while \(g \ll |\Delta|\); compare the upper end of the sweep with the detuning \(\Delta\) printed above.

g_sweep = jnp.linspace(0.01, 0.25, 20)
chi_sweep = jnp.zeros_like(g_sweep)

for i, g_i in enumerate(g_sweep):
    H_i = (
        nk.operator.LocalOperator(hi_combined, np.asarray(H_qubit), acting_on=[0])
        + nk.operator.LocalOperator(hi_combined, np.asarray(H_res), acting_on=[1])
        + nk.operator.LocalOperator(
            hi_combined,
            np.asarray(g_i * jnp.kron(a_plus_adag, a_plus_adag)),
            acting_on=[0, 1],
        )
    )
    ev_i, evc_i = nk.exact.full_ed(H_i, compute_eigenvectors=True)
    si = jnp.argsort(ev_i)
    ev_i, evc_i = ev_i[si], evc_i[:, si]

    def _fd(qi: int, ri: int, _evc: np.ndarray = evc_i) -> int:
        bv = jnp.zeros(n_levels**2)
        bv = bv.at[qi * n_levels + ri].set(1.0)
        return int(jnp.argmax(jnp.abs(_evc.T @ bv) ** 2))

    j00, j10, j01, j11 = _fd(0, 0), _fd(1, 0), _fd(0, 1), _fd(1, 1)
    chi_sweep = chi_sweep.at[i].set((ev_i[j11] - ev_i[j10]) - (ev_i[j01] - ev_i[j00]))

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(g_sweep * 1e3, chi_sweep * 1e3, "o-", markersize=5, label="NetKet (numerical)")
ax.axhline(0, color="gray", linestyle="--", alpha=0.5)
ax.set_xlabel("Coupling strength $g$ (MHz)")
ax.set_ylabel(r"Dispersive shift $\chi$ (MHz)")
ax.set_title(
    rf"Dispersive shift vs. coupling ($\omega_\text{{q}} = {f01:.2f}\,\text{{GHz}}$, "
    rf"$\omega_\text{{r}} = {omega_r_val:.1f}\,\text{{GHz}}$)"
)
ax.legend()
ax.grid(True, alpha=0.3)
fig.tight_layout()
plt.show()
../_images/0aa106e4c0013c35da30162ebd7757cbd727c82f72d02ffeba366d9b607bbfdc.svg ../_images/0aa106e4c0013c35da30162ebd7757cbd727c82f72d02ffeba366d9b607bbfdc_dark.svg

Vacuum Rabi Oscillations with TDVP#

NetKet also evolves states in time. The experimental nkx.TDVP driver integrates the time-dependent variational principle [CT17]: at every step it solves \(S\,\dot{\theta} = -\mathrm{i}\,F\) for the parameter velocities, with \(S\) the quantum geometric tensor and \(F\) the energy gradient.

Here we tune the transmon from the spectrum above to within \(\delta = g\) of the readout resonator and start with one excitation in the qubit, \(|1, 0\rangle\). In the frame rotating at \(\omega_\text{r}\) and within the rotating-wave approximation, the excitation swaps partially into the resonator and back at the generalised vacuum Rabi frequency \(\sqrt{4g^2 + \delta^2}\) [BGGW21]. The ansatz is LogStateVector, one complex parameter \(\log\psi(\sigma)\) per basis state, so we can check TDVP against direct exponentiation. Two caveats come with the log parametrisation. First, no amplitude may pass through zero, which is why we detune: on resonance the qubit amplitude \(\cos(2\pi g t)\) hits zero and the evolution diverges. Second, every unoccupied state needs a small finite amplitude to start from. Swap in any other Flax network and the same driver evolves it.

n_levels_rabi = 4
hi_rabi = nk.hilbert.Fock(n_max=n_levels_rabi - 1, N=2)
a_op = np.diag(np.sqrt(np.arange(1, n_levels_rabi)), k=1)  # lowering operator
detuning_rabi = g_val  # qubit-resonator detuning, in GHz

H_rabi = (
    2
    * np.pi
    * (
        nk.operator.LocalOperator(
            hi_rabi,
            np.diag([
                detuning_rabi * n + float(alpha_nk) / 2 * n * (n - 1)
                for n in range(n_levels_rabi)
            ]),
            acting_on=[0],
        )
        + nk.operator.LocalOperator(
            hi_rabi,
            g_val * (np.kron(a_op, a_op.T) + np.kron(a_op.T, a_op)),
            acting_on=[0, 1],
        )
    )
)  # angular frequency in rad/ns, so time is in ns
n_qubit_op = nk.operator.LocalOperator(
    hi_rabi, np.diag(np.arange(n_levels_rabi, dtype=float)), acting_on=[0]
)

# Start in |1, 0>; a small floor keeps every log-amplitude finite.
psi_initial = np.full(hi_rabi.n_states, 1e-3, dtype=complex)
psi_initial[hi_rabi.states_to_numbers(np.array([[1, 0]]))[0]] = 1.0
psi_initial /= np.linalg.norm(psi_initial)

vs_rabi = nk.vqs.FullSumState(
    hi_rabi, nk.models.LogStateVector(hi_rabi, param_dtype=complex)
)
vs_rabi.parameters = {"logstate": jnp.log(jnp.asarray(psi_initial))}

t_final = 1.5 / g_val  # in ns
t_points = np.linspace(0.0, t_final, 121)
tdvp = nkx.TDVP(
    H_rabi,
    variational_state=vs_rabi,
    ode_solver=nkx.dynamics.RK4(dt=0.005),
    # log ψ is the parameter itself, so the ansatz is holomorphic
    # pyrefly: ignore[bad-argument-type]  # NetKet accepts a QGT constructor here
    qgt=partial(nk.optimizer.qgt.QGTJacobianDense, holomorphic=True),
)
tdvp_log = nk.logging.RuntimeLog()
tdvp.run(
    T=t_final,
    out=tdvp_log,
    obs={"n_qubit": n_qubit_op},
    tstops=t_points,
    show_progress=False,
)

# Reference: exact propagation of the same initial state
H_rabi_dense = H_rabi.to_dense()
n_qubit_dense = n_qubit_op.to_dense()
n_qubit_exact = [
    float(np.real(np.vdot(psi_t, n_qubit_dense @ psi_t)))
    for psi_t in (
        scipy.linalg.expm(-1j * H_rabi_dense * t) @ psi_initial for t in t_points
    )
]

n_qubit_tdvp = np.asarray(tdvp_log["n_qubit"]["Mean"].real)
print(
    f"Max |TDVP - exact| in qubit population: "
    f"{np.max(np.abs(n_qubit_tdvp - n_qubit_exact)):.1e}"
)

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(t_points, n_qubit_exact, label="Exact propagation")
ax.plot(
    np.asarray(tdvp_log["n_qubit"].iters),
    n_qubit_tdvp,
    "--",
    label="NetKet TDVP",
)
ax.set_xlabel("Time (ns)")
ax.set_ylabel(r"Qubit population $\langle \hat{n}_\text{q} \rangle$")
ax.set_title(
    rf"Vacuum Rabi oscillations, $g = \delta = {g_val * 1e3:.0f}\,\text{{MHz}}$"
)
ax.legend()
ax.grid(True, alpha=0.3)
fig.tight_layout()
plt.show()

Max |TDVP - exact| in qubit population: 6.8e-04
../_images/0859501c60e5b37cda647f33b0656e2e3d63bda105baef516682d62481a69346.svg ../_images/0859501c60e5b37cda647f33b0656e2e3d63bda105baef516682d62481a69346_dark.svg

Design Workflow: From Target \(\chi\) to Layout Parameters#

The purpose of computing \(\chi\) is to design a qubit–resonator system with a target dispersive shift. The workflow is:

  1. Choose target \(\chi\) based on readout speed requirements

  2. Select \(\omega_\text{q}\), \(\omega_\text{r}\), \(\alpha\) from qubit design goals

  3. Determine \(g\) from target \(\chi\) using qpdk helpers

  4. Convert to circuit parameters (\(C_\Sigma\), \(L_\text{J}\), \(C_\text{c}\))

  5. Convert to layout parameters (resonator length, capacitor geometry)

FROM A READOUT TARGET TO LAYOUT PARAMETERS Each stage fixes the quantities that the next stage needs. steps 1–2 TARGETS 𝜒 , 𝛼 with qubit and resonator frequencies step 3 COUPLING 𝑔 solved from the target 𝜒 step 4 CIRCUIT 𝐶 Σ , 𝐿 J , 𝐶 c from qpdk conversion helpers step 5 LAYOUT length of the readout resonator

Steps 1–3: From target \(\chi\) to coupling \(g\)#

# Design targets
chi_target_mhz = -1.0  # Target dispersive shift in MHz
chi_target = chi_target_mhz * 1e-3  # Convert to GHz

# Qubit and resonator parameters (from NetKet calculation above)
EJ_design = EJ
EC_design = EC
omega_r_design = omega_r_val

# Use the analytical transmon approximation for design
omega_t_design, alpha_design = ej_ec_to_frequency_and_anharmonicity(
    EJ_design, EC_design
)

# Find g using the qpdk analytical inversion
g_design = dispersive_shift_to_coupling(
    chi_target, omega_t_design, omega_r_design, alpha_design
)

display(
    Math(rf"""
\textbf{{Step 1–3: Hamiltonian Design}} \\
\text{{Target:}}\quad \chi = {chi_target_mhz:.1f}\,\text{{MHz}} \\
E_\text{{J}} = {EJ_design:.1f}\,\text{{GHz}}, \quad
E_\text{{C}} = {EC_design:.3f}\,\text{{GHz}} \\
\omega_\text{{q}} = {omega_t_design:.3f}\,\text{{GHz}}, \quad
\alpha \approx -E_\text{{C}} = -{EC_design:.3f}\,\text{{GHz}} \\
\omega_\text{{r}} = {omega_r_design:.1f}\,\text{{GHz}} \\
\text{{Required coupling:}}\quad g = {float(g_design) * 1e3:.1f}\,\text{{MHz}}
""")
)
\[\begin{split}\displaystyle \textbf{Step 1–3: Hamiltonian Design} \\ \text{Target:}\quad \chi = -1.0\,\text{MHz} \\ E_\text{J} = 20.0\,\text{GHz}, \quad E_\text{C} = 0.218\,\text{GHz} \\ \omega_\text{q} = 5.684\,\text{GHz}, \quad \alpha \approx -E_\text{C} = -0.218\,\text{GHz} \\ \omega_\text{r} = 7.0\,\text{GHz} \\ \text{Required coupling:}\quad g = 68.1\,\text{MHz} \end{split}\]

Step 4: Convert to circuit parameters#

Using the qpdk helper functions, we convert the Hamiltonian parameters to circuit parameters:

  • Total qubit capacitance \(C_\Sigma\) from \(E_\text{C}\) via ec_to_capacitance

  • Josephson inductance \(L_\text{J}\) from \(E_\text{J}\) via ej_to_inductance

  • Resonator length for \(\omega_\text{r}\) via resonator_frequency, and the resonator capacitance \(C_\text{r}\) from the CPW impedance and phase velocity

  • Coupling capacitance \(C_\text{c}\) from \(g\) via coupling_strength_to_capacitance

# Total qubit capacitance
C_sigma = float(ec_to_capacitance(EC_design))

# Josephson inductance
L_J = float(ej_to_inductance(EJ_design))

# Determine resonator length for the target frequency
ep_eff, z0 = cpw_parameters(width=10, gap=6)


def _resonator_objective(x: np.ndarray) -> float:
    """Minimise the squared frequency error.

    Returns:
        The squared error between the calculated and design frequency.
    """
    freq = resonator_frequency(
        length=float(x[0]),
        epsilon_eff=float(jnp.real(ep_eff)),
        is_quarter_wave=True,
    )
    return (freq - omega_r_design * 1e9) ** 2


result = scipy.optimize.minimize(_resonator_objective, 4000.0, bounds=[(1000, 20000)])
resonator_length = result.x[0]

# Total resonator capacitance from CPW impedance and phase velocity
# {cite:p}`gopplCoplanarWaveguideResonators2008a`:
# C_r = l / Re(Z_0 * v_p)
C_r = float(1 / jnp.real(z0 * (c_0 / jnp.sqrt(ep_eff))) * resonator_length * 1e-6)  # F

# Coupling capacitance
C_c = float(
    coupling_strength_to_capacitance(
        g_ghz=float(g_design),
        c_sigma=C_sigma,
        c_r=C_r,
        f_q_ghz=omega_t_design,
        f_r_ghz=omega_r_design,
    )
)

display(
    Math(rf"""
\textbf{{Step 4: Circuit Parameters}} \\
C_\Sigma = {C_sigma * 1e15:.1f}\,\text{{fF}} \\
L_\text{{J}} = {L_J * 1e9:.2f}\,\text{{nH}} \\
C_\text{{r}} = {C_r * 1e15:.1f}\,\text{{fF}} \\
C_\text{{c}} = {C_c * 1e15:.2f}\,\text{{fF}}
""")
)
\[\begin{split}\displaystyle \textbf{Step 4: Circuit Parameters} \\ C_\Sigma = 89.0\,\text{fF} \\ L_\text{J} = 8.17\,\text{nH} \\ C_\text{r} = 724.2\,\text{fF} \\ C_\text{c} = 5.48\,\text{fF} \end{split}\]

Step 5: Layout parameters#

The circuit parameters map directly to layout dimensions:

f_resonator_achieved = resonator_frequency(
    length=resonator_length,
    epsilon_eff=float(jnp.real(ep_eff)),
    is_quarter_wave=True,
)

display(
    Math(rf"""
\textbf{{Step 5: Layout Parameters}} \\
\text{{Resonator length:}}\quad l = {resonator_length:.0f}\,\text{{µm}} \\
\text{{Resonator CPW width:}}\quad w = 10\,\text{{µm}}, \quad
\text{{gap}} = 6\,\text{{µm}} \\
\text{{Resonator frequency achieved:}}\quad
f_\text{{r}} = {f_resonator_achieved / 1e9:.3f}\,\text{{GHz}}
""")
)
\[\begin{split}\displaystyle \textbf{Step 5: Layout Parameters} \\ \text{Resonator length:}\quad l = 4348\,\text{µm} \\ \text{Resonator CPW width:}\quad w = 10\,\text{µm}, \quad \text{gap} = 6\,\text{µm} \\ \text{Resonator frequency achieved:}\quad f_\text{r} = 7.000\,\text{GHz} \end{split}\]

Readout System Design Considerations#

A complete readout design must also account for the resonator linewidth \(\kappa\), the Purcell decay rate, and the measurement-induced dephasing.

# Resonator quality factor and linewidth
Q_ext = 10_000  # External quality factor
kappa = resonator_linewidth_from_q(omega_r_design, Q_ext)

# Purcell decay rate
gamma_purcell = purcell_decay_rate(
    float(g_design), omega_t_design, omega_r_design, kappa
)
T_purcell = 1 / (gamma_purcell * 1e9) if gamma_purcell > 0 else float("inf")

# Measurement-induced dephasing
n_bar = 5.0  # Mean photon number during readout
gamma_phi = measurement_induced_dephasing(chi_target, kappa, n_bar)

# SNR estimate: the signal-to-noise ratio scales as χ/κ
chi_over_kappa = abs(chi_target) / kappa

display(
    Math(rf"""
\textbf{{Readout System Parameters:}} \\
Q_\text{{ext}} = {Q_ext:,} \\
\kappa = \omega_\text{{r}} / Q_\text{{ext}} = {kappa * 1e6:.1f}\,\text{{kHz}} \\
\gamma_\text{{Purcell}} = \kappa (g/\Delta)^2
  = {gamma_purcell * 1e6:.2f}\,\text{{kHz}}
  \quad \Rightarrow \quad T_\text{{Purcell}}
  = {T_purcell * 1e6:.0f}\,\text{{µs}} \\
\Gamma_\phi(\bar{{n}}={n_bar:.0f})
  = 8\chi^2 \bar{{n}} / \kappa
  = {gamma_phi * 1e6:.2f}\,\text{{kHz}} \\
|\chi| / \kappa = {chi_over_kappa:.1f}
  \quad (\text{{want}} \gtrsim 1 \text{{ for high-fidelity readout}})
""")
)
\[\begin{split}\displaystyle \textbf{Readout System Parameters:} \\ Q_\text{ext} = 10,000 \\ \kappa = \omega_\text{r} / Q_\text{ext} = 700.0\,\text{kHz} \\ \gamma_\text{Purcell} = \kappa (g/\Delta)^2 = 1.87\,\text{kHz} \quad \Rightarrow \quad T_\text{Purcell} = 534\,\text{µs} \\ \Gamma_\phi(\bar{n}=5) = 8\chi^2 \bar{n} / \kappa = 57142.86\,\text{kHz} \\ |\chi| / \kappa = 1.4 \quad (\text{want} \gtrsim 1 \text{ for high-fidelity readout}) \end{split}\]

qpdk Layout Components#

The design parameters computed above map directly to physical layout cells in qpdk. Below we show two transmon geometries that can be used for fabrication: the double-pad transmon and the Xmon transmon.

Each cell exposes design parameters (pad dimensions, gap widths, junction dimensions) that set the circuit’s capacitances and Josephson energy, and through them the Hamiltonian parameters \(E_\text{J}\), \(E_\text{C}\), and \(g\).

c_dpt = double_pad_transmon()
c_dpt.plot()
plt.title("Double-Pad Transmon")
plt.show()
../_images/23a53381c2bed4bc52ecc89fcffab3efb32f4ad5efb5014406e578eeb28371d6.svg ../_images/23a53381c2bed4bc52ecc89fcffab3efb32f4ad5efb5014406e578eeb28371d6_dark.svg 
c_xmon = xmon_transmon()
c_xmon.plot()
plt.title("Xmon Transmon")
plt.show()
../_images/a37364dfb6b8e5b191fcd01d450bdff0c98cd5cde078abecd216e801a0ebb890.svg ../_images/a37364dfb6b8e5b191fcd01d450bdff0c98cd5cde078abecd216e801a0ebb890_dark.svg

Design parameter mapping:

Layout parameter

Circuit parameter

Hamiltonian parameter

Pad dimensions, gap width

\(C_\Sigma\) (total capacitance)

\(E_\text{C} = e^2 / (2 C_\Sigma)\)

Junction area, critical current

\(L_\text{J}\) (Josephson inductance)

\(E_\text{J} = \Phi_0^2 / (4\pi^2 L_\text{J})\)

Coupling finger length / gap

\(C_\text{c}\) (coupling capacitance)

\(g \propto C_\text{c} / \sqrt{C_\Sigma C_\text{r}}\)

CPW length, width, gap

\(Z_0\), \(\varepsilon_\text{eff}\)

\(\omega_\text{r}\) (resonator frequency)

Summary: Complete Design Table#

The table below summarises the full design flow from Hamiltonian parameters to layout parameters, with the dispersive shift as the central design target.

data = [
    ("Target", "Target dispersive shift $\\chi$", f"{chi_target_mhz:.1f}", "MHz"),
    ("Hamiltonian", "$E_\\text{J}$", f"{EJ_design:.1f}", "GHz"),
    ("Hamiltonian", "$E_\\text{C}$", f"{EC_design:.3f}", "GHz"),
    ("Hamiltonian", "$\\omega_\\text{q}$ (NetKet)", f"{f01:.4f}", "GHz"),
    ("Hamiltonian", "$\\omega_\\text{q}$ (approx.)", f"{omega_t_design:.4f}", "GHz"),
    ("Hamiltonian", "$\\alpha$ (NetKet)", f"{alpha_nk:.4f}", "GHz"),
    ("Hamiltonian", "$\\omega_\\text{r}$", f"{omega_r_design:.1f}", "GHz"),
    ("Hamiltonian", "$g$", f"{float(g_design) * 1e3:.1f}", "MHz"),
    ("Circuit", "$C_\\Sigma$", f"{C_sigma * 1e15:.1f}", "fF"),
    ("Circuit", "$L_\\text{J}$", f"{L_J * 1e9:.2f}", "nH"),
    ("Circuit", "$C_\\text{r}$", f"{C_r * 1e15:.1f}", "fF"),
    ("Circuit", "$C_\\text{c}$", f"{C_c * 1e15:.2f}", "fF"),
    ("Layout", "Resonator length", f"{resonator_length:.0f}", "µm"),
    ("Layout", "Resonator CPW width", "10", "µm"),
    ("Layout", "Resonator CPW gap", "6", "µm"),
    ("Layout", "Resonator freq", f"{f_resonator_achieved / 1e9:.3f}", "GHz"),
    ("Readout", "$Q_{\\text{ext}}$", f"{Q_ext:,}", ""),
    ("Readout", "$\\kappa$", f"{kappa * 1e6:.1f}", "kHz"),
    ("Readout", "$T_{\\text{Purcell}}$", f"{T_purcell * 1e6:.0f}", "µs"),
    ("Readout", "$|\\chi|/\\kappa$", f"{chi_over_kappa:.1f}", ""),
]

df = pl.DataFrame(data, schema=["Category", "Parameter", "Value", "Unit"], orient="row")

display_dataframe(df)
Category Parameter Value Unit
Target Target dispersive shift χ -1.0 MHz
Hamiltonian EJ 20.0 GHz
Hamiltonian EC 0.218 GHz
Hamiltonian ωq (NetKet) 5.6751 GHz
Hamiltonian ωq (approx.) 5.6840 GHz
Hamiltonian α (NetKet) -0.2396 GHz
Hamiltonian ωr 7.0 GHz
Hamiltonian g 68.1 MHz
Circuit CΣ 89.0 fF
Circuit LJ 8.17 nH
Circuit Cr 724.2 fF
Circuit Cc 5.48 fF
Layout Resonator length 4348 µm
Layout Resonator CPW width 10 µm
Layout Resonator CPW gap 6 µm
Layout Resonator freq 7.000 GHz
Readout Qext 10,000
Readout κ 700.0 kHz
Readout TPurcell 534 µs
Readout |χ|/κ 1.4

Josephson Junction Arrays: a Many-Body Problem#

A transmon has one degree of freedom, so exact diagonalisation is all it needs. A chain of \(N\) junctions, such as the superinductor of a fluxonium [MKGD09, MPK+12], has one per island, and the Hilbert space grows as \((2n_\text{cut}+1)^{N-1}\). For the 43-junction josephson_junction_array() cell with \(n_\text{cut} = 8\) that is \(17^{42}\) states, far beyond diagonalisation. This is where NetKet’s Monte Carlo variational methods earn their place [CT17].

With both ends of the chain held at zero phase, the islands \(i = 1, \dots, N-1\) carry Cooper-pair numbers \(\hat{n}_i\) and phases \(\hat{\varphi}_i\). Quantising the circuit node by node [VD17] gives the standard Josephson-array Hamiltonian [FvandZant01]

(5)#\[\hat{H} = \hat{\boldsymbol{n}}^\mathsf{T} \mathbf{E}_\text{C}\, \hat{\boldsymbol{n}} - E_\text{J} \sum_{j=1}^{N} \cos\left(\hat{\varphi}_j - \hat{\varphi}_{j-1}\right), \qquad \mathbf{E}_\text{C} = 2 e^2 \mathbf{C}^{-1},\]

with \(\hat{\varphi}_0 = \hat{\varphi}_N = 0\) and the tridiagonal capacitance matrix \(C_{ii} = 2 C_\text{J} + C_\text{g}\), \(C_{i,i\pm1} = -C_\text{J}\) [MPK+12]. Because \(C_\text{g} \ll C_\text{J}\), \(\mathbf{C}^{-1}\) is dense: every island talks to every other one through the charging term. Each \(\cos(\hat{\varphi}_j - \hat{\varphi}_{j-1})\) moves one Cooper pair between neighbouring islands.

The helper below reads \(N\), the junction area and \(I_\text{c}\) straight from the cell’s info, so it works for any variant of the cell. The cell does not report capacitances, so we assume a specific junction capacitance of \(50\,\text{fF}/\mu\text{m}^2\), typical of Al/AlOx junctions, and take \(C_\text{g}\) from the default of the array’s SAX model josephson_junction_array(). Replace both with values for your own process.

JUNCTION_SPECIFIC_CAPACITANCE = 50e-15  # F/µm², assumed Al/AlOx value
ISLAND_GROUND_CAPACITANCE = 0.04e-15  # F, SAX model default


@dataclass(frozen=True)
class JunctionArrayProblem:
    """Circuit data and NetKet operators of a junction array, energies in GHz."""

    n_junctions: int
    ic: float  # critical current per junction, A
    c_j: float  # junction capacitance, F
    c_g: float  # island ground capacitance, F
    e_j: float
    charging: np.ndarray  # E_C matrix, H_C = nᵀ E_C n
    incidence: np.ndarray  # junction j measures φ_j - φ_{j-1}
    n_cut: int
    hilbert: nk.hilbert.Fock
    hamiltonian: nk.operator.LocalOperator
    cos_ops: list[nk.operator.LocalOperator]


def junction_array_hamiltonian(
    component,
    *,
    specific_capacitance: float = JUNCTION_SPECIFIC_CAPACITANCE,
    ground_capacitance: float = ISLAND_GROUND_CAPACITANCE,
    n_cut: int = 8,
) -> JunctionArrayProblem:
    """Build the charge-basis Hamiltonian of a ``josephson_junction_array`` cell.

    Energies are in GHz.  Each island is a ``Fock`` site whose occupation
    ``k`` encodes the Cooper-pair number ``n = k - n_cut``.  The ground-state
    charge spread per island is about 1.4 Cooper pairs for the default cell,
    so ``n_cut = 8`` keeps the truncation error well below the VMC error bars.

    Returns:
        The array parameters, NetKet Hilbert space, Hamiltonian and per-junction
        ``cos δφ`` operators.
    """
    n_junctions = int(component.info["n_junctions"])
    ic = component.info["critical_current_nA"] * 1e-9
    c_j = specific_capacitance * component.info["junction_area"]
    n_islands = n_junctions - 1

    e_j = ic * Φ_0 / (2 * np.pi) / h / 1e9
    capacitance = (2 * c_j + ground_capacitance) * np.eye(n_islands) - c_j * (
        np.eye(n_islands, k=1) + np.eye(n_islands, k=-1)
    )
    charging = 2 * e**2 * np.linalg.inv(capacitance) / h / 1e9

    hilbert = nk.hilbert.Fock(n_max=2 * n_cut, N=n_islands)
    n_local = np.diag(np.arange(-n_cut, n_cut + 1, dtype=float))
    raise_local = np.eye(2 * n_cut + 1, k=-1)  # |n+1><n|, i.e. exp(i φ)
    cos_ground = (raise_local + raise_local.T) / 2
    cos_pair = (
        np.kron(raise_local, raise_local.T) + np.kron(raise_local.T, raise_local)
    ) / 2

    hamiltonian = nk.operator.LocalOperator(hilbert, dtype=float)
    for i in range(n_islands):
        hamiltonian += nk.operator.LocalOperator(
            hilbert, charging[i, i] * n_local @ n_local, [i]
        )
        for j in range(i + 1, n_islands):
            hamiltonian += nk.operator.LocalOperator(
                hilbert, 2 * charging[i, j] * np.kron(n_local, n_local), [i, j]
            )
    # cos(φ_j - φ_{j-1}) for every junction j, the two end ones to ground
    cos_ops = [nk.operator.LocalOperator(hilbert, cos_ground, [0])]
    cos_ops += [
        nk.operator.LocalOperator(hilbert, cos_pair, [i, i + 1])
        for i in range(n_islands - 1)
    ]
    cos_ops += [nk.operator.LocalOperator(hilbert, cos_ground, [n_islands - 1])]
    for cos_op in cos_ops:
        hamiltonian -= e_j * cos_op

    incidence = np.eye(n_junctions, n_islands) - np.eye(n_junctions, n_islands, k=-1)
    return JunctionArrayProblem(
        n_junctions=n_junctions,
        ic=ic,
        c_j=c_j,
        c_g=ground_capacitance,
        e_j=e_j,
        charging=charging,
        incidence=incidence,
        n_cut=n_cut,
        hilbert=hilbert,
        hamiltonian=hamiltonian,
        cos_ops=cos_ops,
    )


def linear_modes(problem: JunctionArrayProblem, cos_mean=1.0):
    """Normal-mode frequencies and Gaussian ground state of the linearised array.

    Linearising ``-E_J cos δφ`` around each junction's ``<cos δφ>`` gives
    ``H = nᵀ E_C n + ½ φᵀ B φ`` with ``B = Dᵀ diag(E_J <cos>) D``.

    Returns:
        Mode frequencies in GHz and the charge-space precision matrix ``K`` of
        the ground state ``ψ(n) ∝ exp(-½ nᵀ K n)``.
    """
    stiffness = (
        problem.incidence.T
        @ np.diag(problem.e_j * np.broadcast_to(cos_mean, problem.n_junctions))
        @ problem.incidence
    )
    charging_sqrt = scipy.linalg.sqrtm(2 * problem.charging).real
    charging_sqrt_inv = np.linalg.inv(charging_sqrt)
    frequencies = np.sqrt(np.linalg.eigvalsh(charging_sqrt @ stiffness @ charging_sqrt))
    # Phase-space precision G solves G (2 E_C) G = B; the charge-space one is G⁻¹
    phase_precision = (
        charging_sqrt_inv
        @ scipy.linalg.sqrtm(charging_sqrt @ stiffness @ charging_sqrt).real
        @ charging_sqrt_inv
    )
    return frequencies, np.linalg.inv(phase_precision)

Checking the Hamiltonian on a short array#

A four-junction variant of the cell has \(17^3 = 4913\) states, which NetKet diagonalises exactly. Its gap is the lowest normal mode of the chain. Linearising every junction as an inductor \(L_\text{J} = \Phi_0 / 2\pi I_\text{c}\), as the SAX model does, overestimates that gap. The self-consistent harmonic correction replaces \(E_\text{J}\) by \(E_\text{J} \langle \cos\delta\hat{\varphi} \rangle\), the junction’s stiffness averaged over its zero-point phase fluctuations, and brings the linear model to within about 15 MHz of the full Hamiltonian.

short_array = junction_array_hamiltonian(josephson_junction_array(n_junctions=4))
ed_values, ed_vectors = nk.exact.lanczos_ed(
    short_array.hamiltonian, k=2, compute_eigenvectors=True
)
ed_ground = ed_vectors[:, 0]
cos_exact = np.array([
    ed_ground @ (cos_op.to_sparse() @ ed_ground) for cos_op in short_array.cos_ops
])
f_linear, _ = linear_modes(short_array)
f_renormalised, _ = linear_modes(short_array, cos_exact)

display(
    Math(rf"""
\textbf{{Four-junction array}} \quad (E_\text{{J}} = {short_array.e_j:.1f}\,\text{{GHz}},
\ C_\text{{J}} = {short_array.c_j * 1e15:.0f}\,\text{{fF}}) \\
\text{{Exact gap (NetKet ED):}} \quad {ed_values[1] - ed_values[0]:.3f}\,\text{{GHz}} \\
\text{{Lowest linear mode, }} E_\text{{J}}: \quad {f_linear[0]:.3f}\,\text{{GHz}} \\
\text{{Lowest linear mode, }} E_\text{{J}} \langle\cos\delta\varphi\rangle: \quad
  {f_renormalised[0]:.3f}\,\text{{GHz}}
  \quad (\langle\cos\delta\varphi\rangle = {cos_exact.mean():.3f})
""")
)
\[\begin{split}\displaystyle \textbf{Four-junction array} \quad (E_\text{J} = 22.4\,\text{GHz}, \ C_\text{J} = 25\,\text{fF}) \\ \text{Exact gap (NetKet ED):} \quad 11.146\,\text{GHz} \\ \text{Lowest linear mode, } E_\text{J}: \quad 11.754\,\text{GHz} \\ \text{Lowest linear mode, } E_\text{J} \langle\cos\delta\varphi\rangle: \quad 11.159\,\text{GHz} \quad (\langle\cos\delta\varphi\rangle = 0.901) \end{split}\]

Variational Monte Carlo with a custom Flax network#

For the full cell we sample instead of summing. The ansatz is a Flax module with two parts:

  • a Gaussian \(-\tfrac{1}{2}\boldsymbol{n}^\mathsf{T} \mathbf{K} \boldsymbol{n}\) whose precision matrix \(\mathbf{K}\) starts from the value linear_modes returns, the exact ground state of the linearised circuit, so the optimisation begins from the physics the SAX model already contains;

  • a one-hidden-layer network that learns the anharmonic correction, starting from zero.

A cold start from random parameters puts amplitude on states with large net island charge, whose energies run to thousands of GHz, and the optimisation diverges.

MCState draws configurations with a Metropolis sampler [MRR+53]. We give it a custom transition rule, ChargeHopRule, that proposes the moves the junctions themselves make: one Cooper pair onto or off an island, or across to the next one. VMC_SR preconditions the gradient with stochastic reconfiguration [Sor98], and for more parameters than samples it switches to the neural-tangent-kernel (NTK) form automatically [VHS+22]. We first check it against the exact energy of the four-junction array.

class ChargeGaussianNet(nn.Module):
    """log ψ(n) = -½ nᵀ K n + MLP(n), with K initialised to ``precision``."""

    n_cut: int
    precision: tuple[tuple[float, ...], ...]
    hidden: int = 16
    param_dtype: Any = float

    @nn.compact
    def __call__(self, x):
        """Return ``log ψ`` for a batch of island occupations ``x``."""
        n = x - self.n_cut
        k_matrix = self.param(
            "K", lambda _: jnp.asarray(self.precision, dtype=self.param_dtype)
        )
        k_matrix = (k_matrix + k_matrix.T) / 2
        gaussian = -0.5 * jnp.einsum("...i,ij,...j->...", n, k_matrix, n)
        y = nn.Dense(
            self.hidden,
            kernel_init=nn.initializers.normal(1e-2),
            param_dtype=self.param_dtype,
        )(n)
        y = nn.Dense(
            1, kernel_init=nn.initializers.zeros, param_dtype=self.param_dtype
        )(jnp.tanh(y))
        return gaussian + y[..., 0]


class ChargeHopRule(nk.sampler.rules.MetropolisRule):
    """Propose one Cooper pair added to or removed from an island, or moved to its neighbour.

    These are the moves the junctions make.  NetKet's default ``LocalRule``
    redraws a site uniformly from all ``2 n_cut + 1`` charges, which this
    strongly correlated chain almost always rejects.
    """

    def transition(self, sampler, machine, params, sampler_state, key, σ):
        """Return one proposed move per chain; out-of-range moves stay put.

        The kind of move is drawn first and then a site uniformly among those
        it can act on, so every move and its reverse are equally likely and
        the Metropolis acceptance needs no proposal correction.
        """
        key_site, key_step, key_kind = jax.random.split(key, 3)
        n_chains, n_sites = σ.shape
        rows = jnp.arange(n_chains)
        hop = jax.random.bernoulli(key_kind, 0.5, (n_chains,)) & (n_sites > 1)
        # N islands take on or lose a pair, N - 1 bonds pass one along
        n_targets = jnp.where(hop, max(n_sites - 1, 1), n_sites)
        sites = jax.random.randint(key_site, (n_chains,), 0, n_targets)
        steps = jax.random.choice(key_step, jnp.array([-1, 1]), (n_chains,)).astype(
            σ.dtype
        )
        σ_new = σ.at[rows, sites].add(steps)
        σ_new = σ_new.at[rows, jnp.minimum(sites + 1, n_sites - 1)].add(-steps * hop)
        # Proposals that leave the truncated charge range stay put
        n_max = sampler.hilbert.shape[0] - 1
        in_range = jnp.all((σ_new >= 0) & (σ_new <= n_max), axis=1, keepdims=True)
        return jnp.where(in_range, σ_new, σ), None


def run_array_vmc(
    problem: JunctionArrayProblem,
    *,
    n_iter: int = 100,
    n_samples: int = 1024,
    n_chains: int = 16,
    hidden: int = 16,
    learning_rate: float = 0.01,
    seed: int = 0,
    param_dtype: Any = None,
):
    """Optimise ``ChargeGaussianNet`` with VMC_SR and return the state and log.

    ``param_dtype`` defaults to float32 on a TPU, which computes it natively,
    and float64 elsewhere.  It has to be set here: importing qpdk's SAX models
    turns on JAX's float64 mode for the whole session.

    Returns:
        The optimised variational state and the run log.
    """
    if param_dtype is None:
        param_dtype = jnp.float32 if jax.default_backend() == "tpu" else jnp.float64
    _, precision = linear_modes(problem)
    model = ChargeGaussianNet(
        n_cut=problem.n_cut,
        precision=tuple(map(tuple, precision.tolist())),
        hidden=hidden,
        param_dtype=param_dtype,
    )
    sampler = nk.sampler.MetropolisSampler(
        problem.hilbert,
        ChargeHopRule(),
        n_chains=n_chains,
        sweep_size=4 * problem.hilbert.size,
    )
    vstate = nk.vqs.MCState(sampler, model, n_samples=n_samples, seed=seed)
    driver = nk.driver.VMC_SR(
        problem.hamiltonian.to_jax_operator(),
        optax.sgd(learning_rate),
        variational_state=vstate,
        diag_shift=1e-2,
    )
    log = nk.logging.RuntimeLog()
    driver.run(n_iter, out=log, show_progress=False)
    return vstate, log


vs_short, _ = run_array_vmc(short_array)
vs_short.n_samples = 8192
e_short_vmc = vs_short.expect(short_array.hamiltonian.to_jax_operator())

display(
    Math(rf"""
\textbf{{Four-junction array, ground-state energy}} \\
\text{{NetKet ED:}} \quad {ed_values[0]:.3f}\,\text{{GHz}} \\
\text{{VMC\_SR + Flax:}} \quad {e_short_vmc.mean.real:.3f}
  \pm {e_short_vmc.error_of_mean:.3f}\,\text{{GHz}}
""")
)

Automatic SR implementation choice:  QGT
\[\begin{split}\displaystyle \textbf{Four-junction array, ground-state energy} \\ \text{NetKet ED:} \quad -72.205\,\text{GHz} \\ \text{VMC\_SR + Flax:} \quad -72.206 \pm 0.005\,\text{GHz} \end{split}\]

Now the default 43-junction cell, with 42 coupled islands.

array_cell = josephson_junction_array()
array_problem = junction_array_hamiltonian(array_cell)
vs_array, array_log = run_array_vmc(array_problem)

vs_array.n_samples = 8192
e_array_vmc = vs_array.expect(array_problem.hamiltonian.to_jax_operator())
cos_array = np.array([
    vs_array.expect(cos_op.to_jax_operator()).mean.real
    for cos_op in array_problem.cos_ops
])
f_array_linear, _ = linear_modes(array_problem)
# The network's output layer starts at zero, so the untrained state is the
# harmonic Gaussian; evaluate it on the full cosine Hamiltonian as a baseline.
vs_gaussian, _ = run_array_vmc(array_problem, n_iter=0)
vs_gaussian.n_samples = 8192
e_array_gaussian = vs_gaussian.expect(array_problem.hamiltonian.to_jax_operator())

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
ax1.plot(np.asarray(array_log["Energy"]["Mean"].real), lw=1)
ax1.axhline(
    e_array_gaussian.mean.real, color="k", ls="--", label="Harmonic Gaussian start"
)
ax1.set_xlabel("VMC_SR iteration")
ax1.set_ylabel("Energy (GHz)")
ax1.set_title(
    f"{array_problem.n_junctions}-junction array, {array_problem.n_junctions - 1} islands"
)
ax1.legend()
ax1.grid(True, alpha=0.3)

ax2.plot(np.arange(1, array_problem.n_junctions + 1), cos_array, "o-", ms=3)
ax2.set_xlabel("Junction index $j$")
ax2.set_ylabel(r"$\langle \cos\delta\hat{\varphi}_j \rangle$")
ax2.set_title("Phase fluctuations along the chain")
ax2.grid(True, alpha=0.3)
fig.tight_layout()
plt.show()

display(
    Math(rf"""
E_0^\text{{VMC}} = {e_array_vmc.mean.real:.2f} \pm {e_array_vmc.error_of_mean:.2f}\,\text{{GHz}},
\qquad E_\text{{Gaussian}} = {e_array_gaussian.mean.real:.2f}
  \pm {e_array_gaussian.error_of_mean:.2f}\,\text{{GHz}},
\qquad \overline{{\langle\cos\delta\varphi\rangle}} = {cos_array.mean():.3f}
""")
)

Automatic SR implementation choice:  NTK

Automatic SR implementation choice:  NTK
../_images/526aa97127bd98a238f4333aa1dd7fa82d20c32875003931e53c27d30a994b45.svg ../_images/526aa97127bd98a238f4333aa1dd7fa82d20c32875003931e53c27d30a994b45_dark.svg
\[\displaystyle E_0^\text{VMC} = -723.08 \pm 0.04\,\text{GHz}, \qquad E_\text{Gaussian} = -722.75 \pm 0.05\,\text{GHz}, \qquad \overline{\langle\cos\delta\varphi\rangle} = 0.872 \]

VMC is variational, so the true ground state lies at or below its energy. The optimisation starts from the harmonic ground state, and most of the work is already done there: the Gaussian alone captures the bulk of the energy, and the network and the tuned \(K\) lower it by a further fraction of a GHz. Without exact diagonalisation at this size, that drop and the four-junction comparison above are the evidence the ansatz is adequate. The expectation value \(\langle\cos\delta\hat{\varphi}\rangle\) is nearly uniform along the chain, so we feed its average into the SAX model below.

Feeding the result back into the SAX model#

The array’s SAX model josephson_junction_array() treats each junction as a linear inductor set by ic. Passing \(I_\text{c} \langle\cos\delta\hat{\varphi}\rangle\) instead gives the self-consistent harmonic model of the same layout: the junctions are softer, the chain’s inductance rises and its self-resonances move down.

array_freqs = np.linspace(9e9, 12.5e9, 2001)
s21_bare = josephson_junction_array_model(
    f=array_freqs,
    n_junctions=array_problem.n_junctions,
    ic=array_problem.ic,
    capacitance=array_problem.c_j,
    ground_capacitance=array_problem.c_g,
)["o2", "o1"]
s21_dressed = josephson_junction_array_model(
    f=array_freqs,
    n_junctions=array_problem.n_junctions,
    ic=array_problem.ic * cos_array.mean(),
    capacitance=array_problem.c_j,
    ground_capacitance=array_problem.c_g,
)["o2", "o1"]
f_array_dressed, _ = linear_modes(array_problem, cos_array.mean())

fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(
    array_freqs / 1e9,
    20 * np.log10(np.abs(s21_bare)),
    label=r"$I_\text{c}$ (SAX default)",
)
ax.plot(
    array_freqs / 1e9,
    20 * np.log10(np.abs(s21_dressed)),
    label=r"$I_\text{c}\,\langle\cos\delta\varphi\rangle$ (from VMC)",
)
for f_mode in f_array_dressed[:3]:
    ax.axvline(f_mode, color="C1", ls=":", lw=1)
ax.set_xlabel("Frequency (GHz)")
ax.set_ylabel(r"$|S_{21}|$ (dB)")
ax.set_ylim(-60, 5)  # the lossless chain has transmission zeros at -∞ dB
ax.set_title(f"{array_problem.n_junctions}-junction array between 50 Ω ports")
ax.legend()
ax.grid(True, alpha=0.3)
fig.tight_layout()
plt.show()

display(
    Math(rf"""
\text{{Lowest chain mode:}} \quad {f_array_linear[0]:.2f}\,\text{{GHz}}
\ \to\ {f_array_dressed[0]:.2f}\,\text{{GHz}}, \qquad
L_\text{{array}} = {array_cell.info["array_inductance_nH"]:.0f}\,\text{{nH}}
\ \to\ {array_cell.info["array_inductance_nH"] / cos_array.mean():.0f}\,\text{{nH}}
""")
)
../_images/c64f281c58e34ab288b9b80e0740334b674c24f21253b261015de12321a97803.svg ../_images/c64f281c58e34ab288b9b80e0740334b674c24f21253b261015de12321a97803_dark.svg
\[\displaystyle \text{Lowest chain mode:} \quad 10.32\,\text{GHz} \ \to\ 9.64\,\text{GHz}, \qquad L_\text{array} = 314\,\text{nH} \ \to\ 360\,\text{nH} \]

The dotted lines are the first three renormalised normal modes of the chain with grounded ends, from linear_modes. Between 50 Ω ports the chain is loaded, so the \(|S_{21}|\) features sit near, not exactly on, these lines.

Scaling up on a TPU#

Running this section on a TPU in Google Colab

The variational calculation above is plain JAX, so NetKet runs it on whatever accelerator JAX finds, without code changes. To try it on a TPU:

  1. Open the notebook in Colab with the button at the top of the page.

  2. Choose Runtime → Change runtime type, select a TPU hardware accelerator (v2-8 or v5e-1 on the free tier) and save.

  3. Run all cells. The install cell detects the TPU and reinstalls the TPU build of JAX that matches qpdk’s pinned version.

TPUs are available on Colab’s free tier, but not guaranteed: access depends on demand and on your recent usage, and sessions are time-limited. If no TPU is available, a GPU runtime runs the same full calculation; on the CPU runtime the cell below runs a short version of it.

TPUs compute natively in single precision, so on a TPU run_array_vmc switches the network parameters to float32. The rest of the notebook keeps float64, which the resonator optimisation and the exact diagonalisations need. qpdk’s continuous integration only has CPUs, so it does not test the TPU path.

On a TPU the same driver can afford far more samples and a wider network, which reduces the Monte Carlo error and the variational bias in the 43-junction calculation.

backend = jax.default_backend()
on_accelerator = backend in {"tpu", "gpu"}
print(f"JAX backend: {backend}, devices: {jax.devices()}")

tpu_settings = (
    {"n_iter": 300, "n_samples": 16384, "n_chains": 1024, "hidden": 64}
    if on_accelerator
    else {"n_iter": 20, "n_samples": 1024, "n_chains": 16, "hidden": 16}
)
start = time.perf_counter()
vs_large, large_log = run_array_vmc(array_problem, **tpu_settings)
elapsed = time.perf_counter() - start

vs_large.n_samples = 65536 if on_accelerator else 8192
e_large = vs_large.expect(array_problem.hamiltonian.to_jax_operator())
print(
    f"{tpu_settings['n_iter']} VMC_SR iterations with {tpu_settings['n_samples']} samples "
    f"on {backend}: {elapsed:.1f} s ({elapsed / tpu_settings['n_iter'] * 1e3:.0f} ms/iteration)"
)
print(f"E0 = {e_large.mean.real:.2f} ± {e_large.error_of_mean:.2f} GHz")
if not on_accelerator:
    print("No accelerator found: ran a short CPU version. See the admonition above.")

JAX backend: cpu, devices: [CpuDevice(id=0)]
Automatic SR implementation choice:  NTK

20 VMC_SR iterations with 1024 samples on cpu: 6.8 s (340 ms/iteration)

E0 = -723.09 ± 0.03 GHz
No accelerator found: ran a short CPU version. See the admonition above.

Unimon and Gatemon Hamiltonians as NetKet Operators#

qpdk’s qubit models already build their Hamiltonian matrices: unimon_hamiltonian() in the harmonic-oscillator basis and gatemon_hamiltonian() in the charge basis. Any such matrix becomes a NetKet operator on a single Fock site, so it can be diagonalised, coupled to a resonator mode as in the transmon–resonator section, or evolved in time with TDVP.

def matrix_to_netket(matrix) -> nk.operator.LocalOperator:
    """Wrap a qpdk Hamiltonian matrix as a NetKet operator on one Fock site."""
    matrix = np.asarray(matrix)
    hilbert_1 = nk.hilbert.Fock(n_max=matrix.shape[0] - 1, N=1)
    return nk.operator.LocalOperator(hilbert_1, matrix, acting_on=[0])


def netket_frequency_and_anharmonicity(operator) -> tuple[float, float]:
    """f01 and α from NetKet's exact diagonalisation."""
    levels = nk.exact.lanczos_ed(operator, k=3)
    return float(levels[1] - levels[0]), float(levels[2] - 2 * levels[1] + levels[0])

For the gatemon, \(C_\Sigma\) comes from the Palace capacitance lookup gatemon_capacitances() evaluated at the default gatemon() cell’s dimensions, and the junction is the model’s default carbon-nanotube weak link. qpdk has no model yet that maps the unimon layout to \(E_\text{C}\) and \(E_\text{L}\), so it uses the model’s default energies.

gatemon_cell = gatemon()
gatemon_settings = gatemon_cell.settings
C_sigma_gatemon, _, _ = gatemon_capacitances(
    island_size=gatemon_settings["island_size"][0],
    pocket_gap=gatemon_settings["pocket_gap"],
    pocket_gap_top=gatemon_settings["pocket_gap_top"],
    gate_tip_gap=gatemon_settings["gate_tip_gap"],
)
if not np.isfinite(C_sigma_gatemon):
    msg = "Gatemon dimensions are outside the capacitance lookup grid"
    raise ValueError(msg)
EC_gatemon = float(e**2 / (2 * C_sigma_gatemon) / h / 1e9)

qubit_checks = {
    "Gatemon": (
        matrix_to_netket(gatemon_hamiltonian(ec_ghz=EC_gatemon)),
        gatemon_frequency_and_anharmonicity(ec_ghz=EC_gatemon),
    ),
    "Unimon": (
        matrix_to_netket(unimon_hamiltonian()),
        unimon_frequency_and_anharmonicity(),
    ),
}
rows = []
for name, (operator, (f01_model, alpha_model)) in qubit_checks.items():
    f01_netket, alpha_netket = netket_frequency_and_anharmonicity(operator)
    rows.append((
        name,
        f"{float(f01_model):.4f}",
        f"{f01_netket:.4f}",
        f"{float(alpha_model) * 1e3:.1f}",
        f"{alpha_netket * 1e3:.1f}",
    ))

display(
    Math(
        rf"\text{{Gatemon: }} C_\Sigma = {float(C_sigma_gatemon) * 1e15:.1f}\,\text{{fF}}"
        rf" \Rightarrow E_\text{{C}} = {EC_gatemon * 1e3:.0f}\,\text{{MHz}}"
    )
)
display_dataframe(
    pl.DataFrame(
        rows,
        schema=[
            "Qubit",
            "f01 model (GHz)",
            "f01 NetKet (GHz)",
            "α model (MHz)",
            "α NetKet (MHz)",
        ],
        orient="row",
    )
)
\[\displaystyle \text{Gatemon: } C_\Sigma = 51.9\,\text{fF} \Rightarrow E_\text{C} = 373\,\text{MHz}\]
Qubit f01 model (GHz) f01 NetKet (GHz) α model (MHz) α NetKet (MHz)
Gatemon 5.2959 5.2959 -195.5 -195.5
Unimon 1.4108 1.4108 3439.5 3439.5

Reusing NetKet in Other qpdk Flows#

The pattern above is the same for every component:

  1. Layout cell → circuit parameters. Read them from the cell’s settings and info, or from a qpdk model of the cell: capacitance lookups such as double_pad_transmon_capacitance() and gatemon_capacitances(), or the critical current the junction cells report.

  2. Circuit parameters → NetKet operator. Either an existing qpdk Hamiltonian matrix through matrix_to_netket, or a many-body LocalOperator such as junction_array_hamiltonian.

  3. Solve. nk.exact.lanczos_ed for small spaces, FullSumState with VMC or TDVP when the whole space fits in memory, MCState with VMC_SR and a Flax network when it does not.

  4. Back into the SAX models. Renormalised parameters, such as \(I_\text{c}\langle\cos\delta\hat{\varphi}\rangle\) here, plug into the existing SAX models and from there into any qpdk circuit simulation.

References#

[BGGW21] (1,2)

Alexandre Blais, Arne L. Grimsmo, S. M. Girvin, and Andreas Wallraff. Circuit quantum electrodynamics. Reviews of Modern Physics, 93(2):025005, May 2021. doi:10.1103/RevModPhys.93.025005.

[CCH+19]

Giuseppe Carleo, Kenny Choo, Damian Hofmann, James ET Smith, Tom Westerhout, Fabien Alet, Emily J Davis, Stavros Efthymiou, Ivan Glasser, Sheng-Hsuan Lin, Marta Mauri, Guglielmo Mazzola, Christian B Pereira, and Filippo Vicentini. Netket: a machine learning toolkit for many-body quantum systems. SoftwareX, 10:100311, 2019. doi:10.1016/j.softx.2019.100311.

[CT17] (1,2)

Giuseppe Carleo and Matthias Troyer. Solving the quantum many-body problem with artificial neural networks. Science, 355(6325):602–606, February 2017. doi:10.1126/science.aag2302.

[FvandZant01]

Rosario Fazio and Herre van der Zant. Quantum phase transitions and vortex dynamics in superconducting networks. Physics Reports, 355(4):235–334, December 2001. doi:10.1016/S0370-1573(01)00022-9.

[KYG+07] (1,2,3,4)

Jens Koch, Terri M. Yu, Jay Gambetta, A. A. Houck, D. I. Schuster, J. Majer, Alexandre Blais, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf. Charge-insensitive qubit design derived from the Cooper pair box. Physical Review A, 76(4):042319, October 2007. doi:10.1103/PhysRevA.76.042319.

[MKGD09]

Vladimir E. Manucharyan, Jens Koch, Leonid I. Glazman, and Michel H. Devoret. Fluxonium: Single Cooper-Pair Circuit Free of Charge Offsets. Science, 326(5949):113–116, October 2009. doi:10.1126/science.1175552.

[MPK+12] (1,2)

Nicholas A. Masluk, Ioan M. Pop, Archana Kamal, Zlatko K. Minev, and Michel H. Devoret. Microwave Characterization of Josephson Junction Arrays: Implementing a Low Loss Superinductance. Physical Review Letters, 109(13):137002, September 2012. doi:10.1103/PhysRevLett.109.137002.

[MRR+53]

Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller, and Edward Teller. Equation of state calculations by fast computing machines. The Journal of Chemical Physics, 21(6):1087–1092, June 1953. doi:10.1063/1.1699114.

[Ros58]

F. Rosenblatt. The perceptron: A probabilistic model for information storage and organization in the brain. Psychological Review, 65(6):386–408, 1958. doi:10.1037/h0042519.

[Sor98]

Sandro Sorella. Green function monte carlo with stochastic reconfiguration. Physical Review Letters, 80(20):4558–4561, May 1998. doi:10.1103/PhysRevLett.80.4558.

[VHS+22] (1,2)

Filippo Vicentini, Damian Hofmann, Attila Szabó, Dian Wu, Christopher Roth, Clemens Giuliani, Gabriel Pescia, Jannes Nys, Vladimir Vargas-Calderón, Nikita Astrakhantsev, and Giuseppe Carleo. Netket 3: machine learning toolbox for many-body quantum systems. SciPost Phys. Codebases, pages 7, 2022. doi:10.21468/SciPostPhysCodeb.7.

[VD17]

Uri Vool and Michel Devoret. Introduction to quantum electromagnetic circuits. International Journal of Circuit Theory and Applications, 45(7):897–934, July 2017. doi:10.1002/cta.2359.