# Dispersive Shift of a Transmon–Resonator System with scQubits
This notebook demonstrates how to use
[scqubits](https://scqubits.readthedocs.io/en/latest/) [[GK21](#id27)]
to numerically compute the **dispersive shift** of a readout resonator coupled
to a transmon qubit, and how to translate the resulting Hamiltonian
parameters into physical layout parameters using **qpdk**.
This notebook covers the same design workflow as the companion notebook on
**pymablock-based dispersive shift calculation**, but uses full numerical
diagonalization instead of closed-form perturbation theory—making it easy to
compare the two approaches.
## Background
In circuit quantum electrodynamics (cQED), the state of a transmon qubit is
typically measured via the *dispersive readout* technique
[[BGGW21](#id12)]. The qubit is detuned from
a readout resonator, and the qubit-state-dependent frequency shift of the
resonator—the **dispersive shift** $\chi$—allows non-destructive measurement
of the qubit state.
The two response frequencies differ by $\chi$ under the convention used below;
the sketch shows the separation qualitatively, without assigning its sign.
The full transmon–resonator Hamiltonian (without the rotating-wave approximation)
reads:
$$
\mathcal{H} = \omega_\text{t}\, a_\text{t}^\dagger a_\text{t}
+ \frac{\alpha}{2}\, a_\text{t}^{\dagger 2} a_\text{t}^{2}
+ \omega_\text{r}\, a_\text{r}^\dagger a_\text{r}
+ g\,(a_\text{t} + a_\text{t}^\dagger)(a_\text{r} + a_\text{r}^\dagger),
$$
where $\omega_\text{t}$ is the transmon frequency, $\alpha$ its anharmonicity
($\alpha < 0$ for a transmon),
$\omega_\text{r}$ the resonator frequency, and $g$ the coupling strength.
scQubits constructs and diagonalizes this Hamiltonian numerically in the
transmon charge basis and Fock basis [[KYG+07](#id39)].
The dispersive shift is then extracted directly from the dressed eigenvalues:
$$
\chi = (E_{11} - E_{10}) - (E_{01} - E_{00}),
$$
where $E_{ij}$ is the dressed eigenenergy with $i$ transmon excitations and
$j$ resonator photons.
```ipython3
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,scqubits] @ git+https://github.com/gdsfactory/quantum-rf-pdk.git",
])
```
```ipython3
import numpy as np
# Monkeypatch for NumPy 2.0 compatibility with scqubits.
# TODO: Remove once the minimum required scqubits version officially supports
# NumPy 2.x without relying on np.complex_ / np.float_ aliases.
if not hasattr(np, "complex_"):
# pyrefly: ignore[missing-attribute] # alias removed in NumPy 2
np.complex_ = np.complex128 # ruff: ignore[numpy2-deprecation]
if not hasattr(np, "float_"):
# pyrefly: ignore[missing-attribute] # alias removed in NumPy 2
np.float_ = np.float64 # ruff: ignore[numpy2-deprecation]
import jax.numpy as jnp
import polars as pl
import scipy
import scqubits as scq
from IPython.display import Math, display
from matplotlib import pyplot as plt
from qpdk import PDK
from qpdk.helper import display_dataframe
from qpdk.models.constants import c_0
from qpdk.models.cpw import cpw_parameters
PDK.activate()
from qpdk.models.perturbation import (
dispersive_shift,
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,
ec_to_capacitance,
ej_to_inductance,
)
from qpdk.models.resonator import resonator_frequency
```
```myst-ansi
SyntaxWarning: invalid escape sequence '\p'
/home/runner/work/quantum-rf-pdk/quantum-rf-pdk/.venv/lib/python3.12/site-packages/scqubits/core/zeropi.py: 284
```
## Building the System
We use the same Hamiltonian parameters as the companion pymablock notebook to
allow direct numerical comparison. A transmon qubit in the transmon regime
($E_\text{J} \gg E_\text{C}$) satisfies:
$$
\omega_\text{t} \approx \sqrt{8 E_\text{J} E_\text{C}} - E_\text{C}, \qquad \alpha \approx -E_\text{C}.
$$
```ipython3
# Hamiltonian parameters (matching pymablock notebook for comparison)
EJ = 20.0 # Josephson energy in GHz
EC = 0.2 # Charging energy in GHz
omega_r_val = 7.0 # Resonator frequency in GHz
g_val = 0.1 # Coupling strength in GHz
```
```ipython3
transmon = scq.Transmon(
EJ=EJ,
EC=EC,
ng=0.0, # Offset charge
ncut=30, # Charge basis cutoff
truncated_dim=6, # Keep 6 transmon levels
)
resonator = scq.Oscillator(
E_osc=omega_r_val,
truncated_dim=6, # Keep 6 Fock levels
)
omega_t_val, alpha_val = ej_ec_to_frequency_and_anharmonicity(EJ, EC)
display(
Math(rf"""
\textbf{{Hamiltonian Parameters:}} \\
E_\text{{J}} = {EJ:.1f}\,\text{{GHz}}, \quad
E_\text{{C}} = {EC:.1f}\,\text{{GHz}} \\
\omega_\text{{t}} = \sqrt{{8 E_\text{{J}} E_\text{{C}}}} - E_\text{{C}} = {omega_t_val:.3f}\,\text{{GHz}} \\
\alpha \approx -E_\text{{C}} = {alpha_val:.1f}\,\text{{GHz}} \\
\omega_\text{{r}} = {omega_r_val:.1f}\,\text{{GHz}} \\
g = {g_val:.1f}\,\text{{GHz}} \\
\Delta = \omega_\text{{t}} - \omega_\text{{r}} = {omega_t_val - omega_r_val:.3f}\,\text{{GHz}}
""")
)
```
```myst-ansi
```
## Transmon Spectrum
scQubits numerically diagonalizes the transmon Hamiltonian in the charge
basis. The anharmonicity $\alpha$ is extracted from the first two transition
frequencies and should satisfy $\alpha \approx -E_\text{C}$ in the transmon regime.
```ipython3
eigenvals_t = transmon.eigenvals(evals_count=5)
f01 = eigenvals_t[1] - eigenvals_t[0]
f12 = eigenvals_t[2] - eigenvals_t[1]
anharmonicity = f12 - f01
display(
Math(rf"""
\textbf{{Transmon Spectrum (scQubits):}} \\
0\rightarrow 1\ \text{{frequency:}}\ {f01:.3f}\,\text{{GHz}} \\
1\rightarrow 2\ \text{{frequency:}}\ {f12:.3f}\,\text{{GHz}} \\
\text{{Anharmonicity,}}\ \alpha =\ {anharmonicity:.3f}\,\text{{GHz}}
""")
)
```
```myst-ansi
```
## Dispersive Shift via Numerical Diagonalization
We build the coupled Hilbert space and extract the dispersive shift from
the dressed eigenvalues using [(2)](#equation-eq-dispersive-shift-definition-scq).
In circuit QED the dominant coupling is between the transmon charge operator
$n_\text{t}$ and the resonator electric field $(a_\text{r} + a_\text{r}^\dagger)$
[[BGGW21](#id12)]:
$$
\mathcal{H}_\text{int} = g\, n_\text{t}\, (a_\text{r} + a_\text{r}^\dagger),
$$
where $n_\text{t} \propto (a_\text{t} - a_\text{t}^\dagger)/(2i)$ is the transmon charge operator.
This charge-coupling form is equivalent to [(1)](#equation-eq-transmon-resonator-hamiltonian-scq)
when working in the transmon energy eigenbasis.
The dressed states are identified by their maximum overlap with the
corresponding bare states $|i_\text{t}, j_\text{r}\rangle$.
```ipython3
interaction = scq.InteractionTerm(
g_strength=g_val,
operator_list=[
(0, transmon.n_operator), # (subsystem_index, operator)
(1, resonator.annihilation_operator() + resonator.creation_operator()),
],
)
hilbert_space = scq.HilbertSpace([transmon, resonator], interaction_list=[interaction])
# Diagonalize the full Hamiltonian
evals, evecs = hilbert_space.eigensys(evals_count=12)
# Stack into (evals_count, hilbert_dim) matrix where each row is an eigenvector
evecs = np.array([np.array(v.full()).flatten() for v in evecs])
dim_t = transmon.truncated_dim
dim_r = resonator.truncated_dim
def bare_state_vec(qt_idx: int, res_idx: int) -> np.ndarray:
"""Return the bare state |qt_idx, res_idx⟩ in the full Hilbert space."""
vec = np.zeros(dim_t * dim_r)
vec[qt_idx * dim_r + res_idx] = 1.0
return vec
def find_dressed_index(bare_vec: np.ndarray) -> int:
"""Return the dressed-state index with maximum overlap with bare_vec."""
return int(np.argmax(np.abs(evecs @ bare_vec) ** 2))
idx_00 = find_dressed_index(bare_state_vec(0, 0))
idx_10 = find_dressed_index(bare_state_vec(1, 0))
idx_01 = find_dressed_index(bare_state_vec(0, 1))
idx_11 = find_dressed_index(bare_state_vec(1, 1))
E_00, E_10, E_01, E_11 = evals[idx_00], evals[idx_10], evals[idx_01], evals[idx_11]
chi_scqubits = (E_11 - E_10) - (E_01 - E_00)
# Compare to the analytical formula from qpdk.models.perturbation
chi_analytical = dispersive_shift(omega_t_val, omega_r_val, alpha_val, g_val)
display(
Math(rf"""
\chi_{{\text{{scQubits}}}} = {chi_scqubits * 1e3:.3f}\,\text{{MHz}} \\
\chi_{{\text{{analytical}}}} = {chi_analytical * 1e3:.3f}\,\text{{MHz}} \\
\text{{Relative error: }}
{abs(chi_scqubits - chi_analytical) / abs(chi_analytical) * 100:.1f}\,\%
""")
)
```
```myst-ansi
```
The scQubits result agrees well with the analytical perturbation-theory
formula; any residual difference reflects higher-order corrections beyond
the leading $g^2/\Delta$ term.
## Dispersive Shift vs. Coupling Strength
We sweep the coupling strength $g$ to compare the scQubits numerical result
against the analytical formula. The two approaches agree in the dispersive
regime ($g \ll |\Delta|$) and diverge as the system approaches the
resonance condition.
From the quantum model, we can extract physical parameters relevant for PDK design.
The coupling strength $g$ (in the dispersive limit $g\ll \omega_\text{q}, \omega_\text{r}$) can be related to a coupling capacitance $C_\text{c}$ via [[Sav23](#id65)]:
$$
g &\approx \frac{1}{2} \frac{C_\text{c}}{\sqrt{C_{\Sigma} \left( C_\text{r} - \frac{C_\text{c}^2}{C_\text{q}} \right) }} \sqrt{\omega_\text{q}\omega_\text{r}} \\
&\approx \frac{1}{2} \frac{C_\text{c}}{\sqrt{C_{\Sigma} C_\text{r}}} \sqrt{\omega_\text{q}\omega_\text{r}}, \quad \text{for } C_\text{c} \ll C_\text{q}
$$
where $C_\Sigma$ is the total qubit capacitance, $C_{\text{q}}$ is the capacitance between the qubit pads, and $C_\text{r}$ is the total capacitance of the resonator.
```ipython3
_a_plus_adag = resonator.annihilation_operator() + resonator.creation_operator()
g_sweep = jnp.linspace(0.01, 0.3, 200)
chi_sweep_analytical = jnp.array([
dispersive_shift(omega_t_val, omega_r_val, alpha_val, gi) for gi in g_sweep
])
def chi_scqubits_at_g(g: float) -> float:
"""Numerically compute χ via scQubits at a given coupling g.
Returns:
The dispersive shift value.
"""
hs = scq.HilbertSpace(
[transmon, resonator],
interaction_list=[
scq.InteractionTerm(
g_strength=g,
operator_list=[(0, transmon.n_operator), (1, _a_plus_adag)],
)
],
)
ev, evc = hs.eigensys(evals_count=12)
# Stack into (evals_count, hilbert_dim) matrix where each row is an eigenvector
evc = np.array([np.array(v.full()).flatten() for v in evc])
def di(qi: int, ri: int) -> int:
v = np.zeros(dim_t * dim_r)
v[qi * dim_r + ri] = 1.0
return int(np.argmax(np.abs(evc @ v) ** 2))
i00, i10, i01, i11 = di(0, 0), di(1, 0), di(0, 1), di(1, 1)
return float((ev[i11] - ev[i10]) - (ev[i01] - ev[i00]))
g_sweep_coarse = np.linspace(0.01, 0.3, 20)
chi_sweep_scqubits = np.array([chi_scqubits_at_g(gi) for gi in g_sweep_coarse])
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(
g_sweep * 1e3,
chi_sweep_analytical * 1e3,
"-",
linewidth=2,
label="Analytical (qpdk)",
)
ax.plot(
g_sweep_coarse * 1e3,
chi_sweep_scqubits * 1e3,
"o",
markersize=6,
label="scQubits (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{{t}} = {omega_t_val:.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()
```
```myst-ansi
```
## Flux-Tunable Transmon
A key advantage of scQubits over the analytical pymablock approach is its
native support for **flux-tunable transmons** with a SQUID junction.
The effective Josephson energy varies with external flux $\Phi$:
$$
E_\text{J}(\Phi) = E_{\text{J,max}}
\left|\cos\!\left(\pi\frac{\Phi}{\Phi_0}\right)\right|
\sqrt{1 + d^2 \tan^2\!\left(\pi\frac{\Phi}{\Phi_0}\right)},
$$
where $d$ is the junction asymmetry parameter. The qubit frequency and
dispersive shift vary continuously with flux, which scQubits handles without
any change to the perturbation-theory derivation.
Key parameters of the tunable transmon:
- $E_{\text{J,max}}$: maximum Josephson energy at zero flux
- $E_\text{C}$: charging energy (sets the anharmonicity)
- $d$: SQUID junction asymmetry
See the [scqubits documentation](https://scqubits.readthedocs.io/en/latest/guide/qubits/tunable_transmon.html) for more details.
```ipython3
tunable_transmon = scq.TunableTransmon(
EJmax=40.0, # Maximum Josephson energy in GHz
EC=0.2, # Charging energy in GHz
d=0.1, # Junction asymmetry
flux=0.0, # External flux in units of Φ₀
ng=0.0,
ncut=30,
)
tunable_transmon.plot_evals_vs_paramvals(
"flux", np.linspace(-0.5, 0.5, 201), subtract_ground=True, evals_count=5
)
plt.show()
```
```myst-ansi
Spectral data: 0%| | 0/201 [00:00
```
## Design Workflow: From Hamiltonian to Layout Parameters
The purpose of computing $\chi$ is to **design** a qubit–resonator system
with a target dispersive shift. The workflow mirrors the pymablock notebook,
using qpdk analytical helpers for speed and scQubits for numerical verification:
1. **Choose target $\chi$** based on readout speed requirements
2. **Select $\omega_\text{t}$, $\omega_\text{r}$, $\alpha$** from qubit design goals
3. **Determine $g$** from target $\chi$
4. **Convert to circuit parameters** ($C_\Sigma$, $L_\text{J}$, $C_\text{c}$) using qpdk helpers
5. **Convert to layout parameters** (resonator length, capacitor geometry)
### Step 1–3: From target $\chi$ to coupling $g$
```ipython3
# 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
EJ_design = 20.0 # GHz
EC_design = 0.2 # GHz
omega_r_design = 7.0 # GHz
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
)
# Verify analytically and with scQubits
chi_verify_analytical = dispersive_shift(
omega_t_design, omega_r_design, alpha_design, float(g_design)
)
chi_verify_scq = chi_scqubits_at_g(float(g_design))
display(
Math(rf"""
\textbf{{Step 1–3: Hamiltonian Design}} \\
\text{{Target:}}\quad \chi = {chi_target_mhz:.1f}\,\text{{MHz}} \\
\text{{Required coupling:}}\quad g = {float(g_design) * 1e3:.1f}\,\text{{MHz}} \\
\text{{Verification (analytical):}}\quad \chi = {chi_verify_analytical * 1e3:.3f}\,\text{{MHz}} \\
\text{{Verification (scQubits):}}\quad \chi = {chi_verify_scq * 1e3:.3f}\,\text{{MHz}}
""")
)
```
```myst-ansi
```
### 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`
- Coupling capacitance $C_\text{c}$ from $g$ via `coupling_strength_to_capacitance`
```ipython3
# 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:
"""Minimize the squared frequency error.
Returns:
The squared error between the calculated and target 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 = 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}}
""")
)
```
```myst-ansi
```
### Step 5: Layout parameters
The circuit parameters map directly to layout dimensions:
```ipython3
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}}
""")
)
```
```myst-ansi
```
## Readout System Design Considerations
A complete readout design also requires considering the resonator
linewidth $\kappa$, the Purcell decay rate, and the measurement-induced
dephasing.
```ipython3
# 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_verify_analytical, kappa, n_bar)
# SNR estimate: the signal-to-noise ratio scales as χ/κ
chi_over_kappa = abs(chi_verify_analytical) / 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}})
""")
)
```
```myst-ansi
```
## 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.
```ipython3
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:.1f}", "GHz"),
("Hamiltonian", "$\\omega_\\text{t}$", f"{omega_t_design:.3f}", "GHz"),
("Hamiltonian", "$\\alpha$", f"{alpha_design:.1f}", "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)
```
```myst-ansi
._DualFormatTable at 0x7f00d0ad3260>
```
## References