Coverage for qpdk / models / perturbation.py: 86%
36 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-14 10:27 +0000
1r"""Perturbation theory helpers for superconducting quantum circuits.
3This module provides helper functions for computing effective Hamiltonians
4and energy corrections using quasi-degenerate perturbation theory via
5`pymablock <https://pymablock.readthedocs.io/>`_
6:cite:`arayaDayPymablockAlgorithmPackage2025`.
8The primary application is computing the **dispersive shift** :math:`\chi`
9of a resonator coupled to a transmon qubit, which is the key quantity for
10dispersive readout of superconducting qubits
11:cite:`blaisCircuitQuantumElectrodynamics2021,kochChargeinsensitiveQubitDesign2007a`.
13The transmon-resonator Hamiltonian (without the rotating wave approximation)
14is:
16.. math::
18 \mathcal{H} = -\omega_t\, a_t^\dagger a_t
19 + \frac{\alpha}{2}\, a_t^{\dagger 2} a_t^2
20 + \omega_r\, a_r^\dagger a_r
21 - g\,(a_t^\dagger - a_t)(a_r^\dagger - a_r)
23where :math:`\omega_t` is the transmon frequency, :math:`\omega_r` the resonator
24frequency, :math:`\alpha` the anharmonicity, and :math:`g` the coupling strength.
26References:
27 - :cite:`arayaDayPymablockAlgorithmPackage2025`
28 - :cite:`blaisCircuitQuantumElectrodynamics2021`
29 - :cite:`kochChargeinsensitiveQubitDesign2007a`
30"""
32from __future__ import annotations
34from functools import partial
36import jax
37import jax.numpy as jnp
38import sympy
39from jax.typing import ArrayLike
40from sympy.physics.quantum import Dagger
41from sympy.physics.quantum.boson import BosonOp
44def transmon_resonator_hamiltonian() -> tuple[
45 sympy.Expr, sympy.Expr, tuple[sympy.Symbol, ...]
46]:
47 r"""Build the symbolic transmon-resonator Hamiltonian.
49 Constructs the unperturbed (:math:`H_0`) and perturbation (:math:`H_p`)
50 parts of the transmon-resonator Hamiltonian using bosonic operators.
52 The Hamiltonian is split as :math:`H = H_0 + H_p` where:
54 .. math::
56 \begin{aligned}
57 H_0 &= -\omega_t\, a_t^\dagger a_t
58 + \frac{\alpha}{2}\, a_t^{\dagger 2} a_t^2
59 + \omega_r\, a_r^\dagger a_r \\
60 H_p &= -g\,(a_t^\dagger - a_t)(a_r^\dagger - a_r)
61 \end{aligned}
63 Returns:
64 A tuple ``(H_0, H_p, symbols)`` where ``symbols`` is
65 ``(omega_t, omega_r, alpha, g)``.
67 Example:
68 >>> H_0, H_p, (omega_t, omega_r, alpha, g) = transmon_resonator_hamiltonian()
69 """
70 omega_t, omega_r, alpha, g = sympy.symbols(
71 r"\omega_{t} \omega_{r} \alpha g", real=True, positive=True
72 )
74 a_t, a_r = BosonOp("a_t"), BosonOp("a_r")
76 H_0 = (
77 -omega_t * Dagger(a_t) * a_t
78 + omega_r * Dagger(a_r) * a_r
79 + alpha * Dagger(a_t) ** 2 * a_t**2 / 2
80 )
82 H_p = -g * (Dagger(a_t) - a_t) * (Dagger(a_r) - a_r)
84 return H_0, H_p, (omega_t, omega_r, alpha, g)
87@partial(jax.jit, inline=True)
88def dispersive_shift(
89 ω_t_ghz: float | ArrayLike,
90 ω_r_ghz: float | ArrayLike,
91 α_ghz: float | ArrayLike,
92 g_ghz: float | ArrayLike,
93) -> float | jax.Array:
94 r"""Compute the dispersive shift numerically.
96 Evaluates the second-order dispersive shift for a transmon coupled
97 to a resonator. Uses the analytical formula derived from perturbation
98 theory (without the rotating wave approximation)
99 :cite:`kochChargeinsensitiveQubitDesign2007a,blaisCircuitQuantumElectrodynamics2021`:
101 .. math::
103 \chi = \frac{2g^2}{\Delta - \alpha}
104 - \frac{2g^2}{\Delta}
105 - \frac{2g^2}{\omega_t + \omega_r + \alpha}
106 + \frac{2g^2}{\omega_t + \omega_r}
108 where :math:`\Delta = \omega_t - \omega_r`. The first two terms give
109 the rotating-wave-approximation (RWA) contribution
111 .. math::
113 \chi_\text{RWA} = \frac{2 \alpha g^2}{\Delta(\Delta - \alpha)}
115 and the last two are corrections from the counter-rotating terms.
117 All parameters are in GHz, and the returned value is also in GHz.
119 Args:
120 ω_t_ghz: Transmon frequency in GHz.
121 ω_r_ghz: Resonator frequency in GHz.
122 α_ghz: Transmon anharmonicity in GHz (positive value).
123 g_ghz: Coupling strength in GHz.
125 Returns:
126 Dispersive shift :math:`\chi` in GHz.
128 Example:
129 >>> χ = dispersive_shift(5.0, 7.0, 0.2, 0.1)
130 >>> print(f"χ = {χ * 1e3:.2f} MHz")
131 """
132 Δ = ω_t_ghz - ω_r_ghz
134 # Full expression including counter-rotating terms
135 return (
136 2 * g_ghz**2 / (Δ - α_ghz)
137 - 2 * g_ghz**2 / Δ
138 - 2 * g_ghz**2 / (ω_t_ghz + ω_r_ghz + α_ghz)
139 + 2 * g_ghz**2 / (ω_t_ghz + ω_r_ghz)
140 )
143@partial(jax.jit, inline=True)
144def dispersive_shift_to_coupling(
145 χ_ghz: float | ArrayLike,
146 ω_t_ghz: float | ArrayLike,
147 ω_r_ghz: float | ArrayLike,
148 α_ghz: float | ArrayLike,
149) -> float | jax.Array:
150 r"""Compute the coupling strength from a target dispersive shift.
152 Inverts the dispersive shift relation to find the coupling strength
153 :math:`g` required to achieve a desired :math:`\chi`. Uses only the
154 dominant rotating-wave term
155 :cite:`kochChargeinsensitiveQubitDesign2007a`:
157 .. math::
159 g \approx \sqrt{\frac{-\chi\,\Delta\,(\Delta - \alpha)}{2\alpha}}
161 where :math:`\Delta = \omega_t - \omega_r`.
163 Note:
164 The expression under the square root may be negative when the
165 sign of the target :math:`\chi` is inconsistent with the detuning
166 and anharmonicity (e.g., positive :math:`\chi` with
167 :math:`\Delta < 0`). In that case the absolute value is taken
168 so that the returned coupling strength is always real and
169 non-negative, but the caller should verify self-consistency
170 via :func:`dispersive_shift`.
172 Args:
173 χ_ghz: Target dispersive shift in GHz (typically negative).
174 ω_t_ghz: Transmon frequency in GHz.
175 ω_r_ghz: Resonator frequency in GHz.
176 α_ghz: Transmon anharmonicity in GHz (positive value;
177 the physical anharmonicity of a transmon is negative, but
178 following the Hamiltonian convention used throughout this
179 module, :math:`\alpha` is taken as positive).
181 Returns:
182 Coupling strength :math:`g` in GHz.
184 Example:
185 >>> g = dispersive_shift_to_coupling(-0.001, 5.0, 7.0, 0.2)
186 >>> print(f"g = {g * 1e3:.1f} MHz")
187 """
188 Δ = ω_t_ghz - ω_r_ghz
189 g_squared = -χ_ghz * Δ * (Δ - α_ghz) / (2 * α_ghz)
190 return jnp.sqrt(jnp.abs(g_squared))
193@partial(jax.jit, inline=True)
194def ej_ec_to_frequency_and_anharmonicity(
195 ej_ghz: float | ArrayLike, ec_ghz: float | ArrayLike
196) -> tuple[float | jax.Array, float | jax.Array]:
197 r"""Convert :math:`E_J` and :math:`E_C` to qubit frequency and anharmonicity.
199 Uses the standard transmon approximations
200 :cite:`kochChargeinsensitiveQubitDesign2007a`:
202 .. math::
204 \begin{aligned}
205 \omega_q &\approx \sqrt{8 E_J E_C} - E_C \\
206 \alpha &\approx E_C
207 \end{aligned}
209 Note:
210 The physical anharmonicity of a transmon is *negative*
211 (:math:`\alpha = -E_C`), but the Hamiltonian convention used
212 in this module and in pymablock takes :math:`\alpha` as positive.
214 Args:
215 ej_ghz: Josephson energy in GHz.
216 ec_ghz: Charging energy in GHz.
218 Returns:
219 Tuple of ``(ω_q_ghz, α_ghz)``.
221 Example:
222 >>> ω_q, α = ej_ec_to_frequency_and_anharmonicity(20.0, 0.2)
223 >>> print(f"ω_q = {ω_q:.2f} GHz, α = {α:.1f} GHz")
224 """
225 return jnp.sqrt(8 * ej_ghz * ec_ghz) - ec_ghz, ec_ghz
228@partial(jax.jit, inline=True)
229def purcell_decay_rate(
230 g_ghz: float | ArrayLike,
231 ω_t_ghz: float | ArrayLike,
232 ω_r_ghz: float | ArrayLike,
233 κ_ghz: float | ArrayLike,
234) -> float | jax.Array:
235 r"""Estimate the Purcell decay rate of a transmon through a resonator.
237 The Purcell effect limits qubit lifetime when coupled to a lossy
238 resonator. In the dispersive regime
239 :cite:`houckControllingSpontaneousEmission2008,blaisCircuitQuantumElectrodynamics2021`:
241 .. math::
243 \gamma_\text{Purcell} = \kappa \left(\frac{g}{\Delta}\right)^2
245 where :math:`\kappa` is the resonator decay rate and
246 :math:`\Delta = \omega_t - \omega_r`.
248 Args:
249 g_ghz: Coupling strength in GHz.
250 ω_t_ghz: Transmon frequency in GHz.
251 ω_r_ghz: Resonator frequency in GHz.
252 κ_ghz: Resonator linewidth (decay rate) in GHz.
254 Returns:
255 Purcell decay rate in GHz.
257 Example:
258 >>> γ = purcell_decay_rate(0.1, 5.0, 7.0, 0.001)
259 >>> T_purcell = 1 / (γ * 1e9)
260 >>> print(f"T_Purcell = {T_purcell * 1e6:.0f} µs")
261 """
262 Δ = ω_t_ghz - ω_r_ghz
263 return κ_ghz * (g_ghz / Δ) ** 2
266@partial(jax.jit, inline=True)
267def resonator_linewidth_from_q(
268 ω_r_ghz: float | ArrayLike, q_ext: float | ArrayLike
269) -> float | jax.Array:
270 r"""Compute resonator linewidth from external quality factor.
272 Converts external quality factor to linewidth (half-linewidth at half-maximum)
273 :cite:`gopplCoplanarWaveguideResonators2008a,m.pozarMicrowaveEngineering2012`:
275 .. math::
277 \kappa = \frac{\omega_r}{Q_\text{ext}}
279 Args:
280 ω_r_ghz: Resonator frequency in GHz.
281 q_ext: External quality factor.
283 Returns:
284 Resonator linewidth :math:`\kappa` in GHz.
286 Example:
287 >>> κ = resonator_linewidth_from_q(7.0, 10_000)
288 >>> print(f"κ = {κ * 1e6:.1f} kHz")
289 """
290 return ω_r_ghz / q_ext
293@partial(jax.jit, inline=True)
294def measurement_induced_dephasing(
295 χ_ghz: float | ArrayLike, κ_ghz: float | ArrayLike, n_bar: float | ArrayLike
296) -> float | jax.Array:
297 r"""Estimate measurement-induced dephasing rate.
299 During dispersive readout, photons in the resonator cause
300 additional dephasing of the qubit
301 :cite:`gambettaQubitphotonInteractionsCavity2006,blaisCircuitQuantumElectrodynamics2021`:
303 .. math::
305 \Gamma_\phi = \frac{8 \chi^2 \bar{n}}{\kappa}
307 where :math:`\bar{n}` is the mean photon number in the resonator.
309 Args:
310 χ_ghz: Dispersive shift in GHz.
311 κ_ghz: Resonator linewidth in GHz.
312 n_bar: Mean photon number in the resonator during measurement.
314 Returns:
315 Measurement-induced dephasing rate in GHz.
317 Example:
318 >>> Γ_φ = measurement_induced_dephasing(-0.001, 0.001, 5.0)
319 >>> print(f"Γ_φ = {Γ_φ * 1e6:.1f} kHz")
320 """
321 return 8 * χ_ghz**2 * n_bar / κ_ghz