pyEPR.calcs.back_box_numeric module#

Numerical diagonalization of quantum Hamiltonian and parameter extraction.

@author: Phil Reinhold, Zlatko Minev, Lysander Christakis

Original code on black_box_hamiltonian and make_dispersive functions by Phil Reinhold Revisions and updates by Zlatko Minev & Lysander Christakis

pyEPR.calcs.back_box_numeric.black_box_hamiltonian(fs, ljs, fzpfs, cos_trunc=5, fock_trunc=8, individual=False, non_linear_potential=None)[source]#

Build the EPR Hamiltonian matrix in Fock space.

Constructs

\[H = H_{\rm lin} + H_{\rm nl} = \sum_m f_m \hat n_m - \sum_j \frac{E_{J,j}}{h}\,{\rm nl}(\hat\varphi_j)\]

where \(\hat\varphi_j = \sum_m \phi_{\rm zpf,jm} (\hat a_m + \hat a_m^\dag)\) (dimensionless reduced phase), and the nonlinear correction nl defaults to cos_approx() (\(\cos\varphi - 1 + \varphi^2/2\), the Taylor series starting at \(\varphi^4\)).

Parameters:

fs (array-like) – Linearised normal-mode frequencies in Hz, length M.
ljs (array-like) – Junction linearised inductances in Henries, length J.
fzpfs (array-like) – Generalised (not reduced) zero-point flux fluctuations in Wb, shape M×J. Internally divided by fluxQ = Φ_0/(2π) to obtain reduced phases.
cos_trunc (int) – Cosine Taylor-series truncation order (default 5). Ignored when non_linear_potential is given.
fock_trunc (int) – Fock-space truncation per mode (default 8).
individual (bool) – If True return [H_lin, H_nl] separately instead of their sum.
non_linear_potential (callable or None) – Custom nonlinear potential: nl(φ_op: Qobj) -> Qobj. Must return the correction beyond the quadratic, i.e. \(V(\varphi) - V(0) - \tfrac12 V''(0)\varphi^2\) normalised by \(|V''(0)|\). Use make_nonlinear_potential() to construct this from any scalar Python function.

Returns:

Full Hamiltonian H/h in Hz, or [H_lin/h, H_nl/h] when individual=True.

Return type:

qutip.Qobj or [qutip.Qobj, qutip.Qobj]

pyEPR.calcs.back_box_numeric.black_box_hamiltonian_nq(freqs, zmat, ljs, cos_trunc=6, fock_trunc=8, show_fit=False)[source]#: N-Qubit version of bbq, based on the full Z-matrix Currently reproduces 1-qubit data, untested on n-qubit data Assume: Solve the model without loss in HFSS.

pyEPR.calcs.back_box_numeric.cos_full_correction(op_cos_arg: Qobj)#

Exact EPR nonlinear correction: cos(φ) - I + φ²/2 (no truncation).

This is the infinite-order equivalent of cos_approx(). In the EPR Hamiltonian the linear eigenfrequencies already account for the harmonic (φ²/2) Josephson energy, so the nonlinear potential that must be subtracted is the remainder:

\[H_{\mathrm{nl}} = -E_J \bigl[\cos(\varphi) - 1 + \tfrac{\varphi^2}{2}\bigr] = -E_J \sum_{n\ge 2} \frac{(-1)^n \varphi^{2n}}{(2n)!}\]

For weakly anharmonic circuits (transmon, φ_zpf ≲ 0.3) the truncated cos_approx() is equivalent and faster. For strongly anharmonic circuits (fluxonium, φ_zpf ≳ 1) use this function to avoid systematic errors from the truncated series.

Parameters:: op_cos_arg (qutip.Qobj) – Phase operator φ (Hermitian), in the full Fock-space tensor product.
Returns:: cos(φ) - I + φ²/2 as a matrix.
Return type:: qutip.Qobj

See also

cos_approx: Truncated Taylor-series version (faster for small φ_zpf).

References

arXiv:2411.15039 — EPR analysis for very anharmonic superconducting circuits.

pyEPR.calcs.back_box_numeric.epr_numerical_diagonalization(freqs, Ljs, φzpf, cos_trunc=8, fock_trunc=9, use_1st_order=False, return_H=False, non_linear_potential=None, use_full_cos=False)[source]#

Numerical diagonalization of the EPR Hamiltonian.

Parameters:

freqs (array-like) – Linearised normal-mode frequencies in GHz (not radians), length M.
Ljs (array-like) – Junction linearised inductances in Henries, length J.
ϕzpf (array-like) – Reduced zero-point flux fluctuations, shape M × J.
cos_trunc (int, optional) – Order of the Taylor-series truncation of cos(φ). Only used when use_full_cos=False (default 8). Ignored when use_full_cos=True.
fock_trunc (int, optional) – Fock-space truncation (number of levels per mode). Default 9.
use_1st_order (bool, optional) – Use first-order perturbation theory instead of full diagonalization. Default False.
return_H (bool, optional) – If True, also return the Hamiltonian Qobj. Default False.
non_linear_potential (callable, optional) –
Custom nonlinear potential function. Must accept a qutip Qobj φ (the reduced phase-fluctuation operator) and return a qutip Qobj equal to the correction beyond the quadratic of the junction potential:

\[\mathrm{nl}(\delta\varphi) = \frac{V(\varphi_0 + \delta\varphi) - V(\varphi_0) - \tfrac{1}{2}V''(\varphi_0)\,\delta\varphi^2} {|V''(\varphi_0)|}\]

Use make_nonlinear_potential() to build this callable from any scalar Python function. Overrides both cos_trunc and use_full_cos.
use_full_cos (bool, optional) – If True, use the exact matrix-exponential cosine cos(φ) = (e^{iφ} + e^{-iφ}) / 2 instead of the truncated Taylor series. Recommended for strongly anharmonic circuits such as fluxonium, where large zero-point phase fluctuations make the truncated expansion inaccurate. Default False (truncated series).

Returns:

f_ND (numpy.ndarray) – Dressed mode frequencies in GHz.
χ_ND (numpy.ndarray) – Cross-Kerr / anharmonicity matrix in MHz (positive = red shift).
Hs (qutip.Qobj) – Full Hamiltonian Qobj — only returned when return_H=True.

Notes

Choosing the right potential path:

Small φ_zpf (transmon, φ_zpf ≲ 0.3): default cos_trunc=8 is fast and accurate.
Large φ_zpf (fluxonium, φ_zpf ≳ 1): set use_full_cos=True to avoid truncation errors. See arXiv:2411.15039.
Non-cosine or flux-biased junction: use non_linear_potential=make_nonlinear_potential(V).

Examples

Fluxonium with exact cosine:

>>> f_ND, chi_ND = epr_numerical_diagonalization(
...     freqs, Ljs, phi_zpf, fock_trunc=25, use_full_cos=True
... )

Flux-biased junction at φ_ext = π/4:

>>> import numpy as np
>>> phi_ext = np.pi / 4
>>> nl = make_nonlinear_potential(lambda phi: np.cos(phi - phi_ext),
...                               phi_min=phi_ext)
>>> f_ND, chi_ND = epr_numerical_diagonalization(
...     freqs, Ljs, phi_zpf, fock_trunc=15, non_linear_potential=nl
... )

pyEPR.calcs.back_box_numeric.make_dispersive(H, fock_trunc, fzpfs=None, f0s=None, chi_prime=False, use_1st_order=False)[source]#

Input: Hamiltonian Matrix.: Optional: phi_zpfs and normal mode frequencies, f0s. use_1st_order : deprecated
Output:: Return dressed mode frequencies, chis, chi prime, phi_zpf flux (not reduced), and linear frequencies
Description:: Takes the Hamiltonian matrix H from bbq_hmt. It them finds the eigenvalues/eigenvectors and assigns quantum numbers to them — i.e., mode excitations, such as, for instance, for three mode, \(|0,0,0\rangle\) or \(|0,0,1\rangle\), which correspond to no excitations in any of the modes or one excitation in the 3rd mode, resp. The assignment is performed based on the maximum overlap between the eigenvectors of H_full and H_lin. Based on the assignment of the excitations, the function returns the dressed mode frequencies \(\omega_m^\prime\), and the cross-Kerr matrix (including anharmonicities) extracted from the numerical diagonalization, as well as from 1st order perturbation theory. Note, the diagonal of the CHI matrix is directly the anharmonicity term.

pyEPR.calcs.back_box_numeric.make_nonlinear_potential(V, phi_min=0.0, _eps=1e-05)[source]#

Build an EPR-compatible nonlinear potential from a user-supplied scalar function.

Background — EPR Hamiltonian split

In the EPR framework the Hamiltonian is decomposed as

\[H = H_{\rm lin} - \sum_j \frac{E_{J,j}}{h}\, \underbrace{\Bigl[V(\varphi_{j,0} + \delta\varphi_j) - V(\varphi_{j,0}) - \tfrac12 V''(\varphi_{j,0})\,\delta\varphi_j^2 \Bigr] / |V''(\varphi_{j,0})|}_{\displaystyle \mathrm{nl}(\delta\varphi_j)}\]

\(H_{\rm lin}\) is built from the HFSS eigenmode frequencies, which already encode the harmonic (quadratic) Josephson energy via the linearised inductance \(L_j\). The nonlinear correction nl(δφ) must therefore subtract exactly the constant and quadratic parts of V that are already captured.

The function returned by make_nonlinear_potential() does this automatically: it evaluates V as an operator (via eigendecomposition of the phase operator), subtracts the constant \(V(\varphi_0)\) and the quadratic term \(\tfrac12 V''(\varphi_0)\,\delta\varphi^2\), and normalises by \(|V''(\varphi_0)|\) for consistency with the effective Josephson energy \(E_{J,\rm eff} = (\Phi_0/2\pi)^2 / L_j\).

Convention: V is the normalised potential appearing in \(U_J = -E_J V(\varphi)\). For a standard Josephson junction \(V(\varphi) = \cos\varphi\), giving \(V''(0) = -1\) and recovering the default cos_full_correction().

Parameters:

V (callable) –
Scalar potential function V(phi: float) -> float.

phi is the reduced phase in the absolute frame (not a fluctuation). Common choices:
- np.cos — standard Josephson junction (reproduces cos_full_correction())
- lambda phi: np.cos(phi - phi_ext) — flux-biased junction at phi_ext
- lambda phi: d_J * np.cos(phi - phi_ext) — asymmetric SQUID (scaled)
phi_min (float, optional) – Phase at the potential minimum (equilibrium / bias point) in the same frame as V. The EPR phase-fluctuation operator δφ is centred here. Default 0 (unbiased junction at its minimum).
_eps (float, optional) – Finite-difference step for estimating \(V''(\varphi_0)\). Default 1e-5.

Returns:

nl(phi_op: qutip.Qobj) -> qutip.Qobj, suitable as the non_linear_potential argument of epr_numerical_diagonalization() or black_box_hamiltonian().

Return type:

callable

Raises:

ValueError – If \(|V''(\varphi_0)| < 10^{-12}\), indicating phi_min is not a proper quadratic extremum of V.

Examples

Standard junction — equivalent to the built-in default:

>>> import numpy as np
>>> nl = make_nonlinear_potential(np.cos)   # V(phi) = cos(phi), phi_min=0
>>> # nl(phi_op) ≡ cos_full_correction(phi_op)

Flux-biased junction at half a flux quantum (φ_ext = π/2):

>>> phi_ext = np.pi / 2
>>> nl = make_nonlinear_potential(lambda phi: np.cos(phi - phi_ext),
...                               phi_min=phi_ext)
>>> f_ND, chi_ND = epr_numerical_diagonalization(
...     freqs, Ljs, phi_zpf, fock_trunc=20,
...     non_linear_potential=nl,
... )

Asymmetric SQUID tuned to flux bias phi_ext, asymmetry d = (Ej1-Ej2)/(Ej1+Ej2):

>>> import numpy as np
>>> phi_ext = 0.4  # external flux in reduced units
>>> d = 0.1       # junction asymmetry
>>> # Effective potential (see arXiv:2010.00620, §IV)
>>> def V_squid(phi):
...     return np.cos(phi) * np.cos(phi_ext) + d * np.sin(phi) * np.sin(phi_ext)
>>> phi_min = np.arctan(-d * np.tan(phi_ext))   # minimum of -Ej*V_squid
>>> nl = make_nonlinear_potential(V_squid, phi_min=phi_min)

Notes

Uses eigendecomposition of φ_op to evaluate V as a matrix function — exact for any analytic V but slightly slower than the Taylor-series cosine default for small φ_zpf. For production runs with standard junctions, prefer cos_full_correction() (use_full_cos=True) or the default cos_trunc path; use make_nonlinear_potential() when the junction potential is not a simple cosine centred at zero.