%load_ext autoreload
%autoreload 2
%config IPCompleter.greedy = True
pyEPR Startup Example#
Single transmon qubit mode analysis#
Author: Zlatko Minev
Introduction: The Energy-Participation Ratio (EPR) Method#
The energy-participation ratio (EPR) method is a quantum circuit design framework that provides a rigorous and efficient bridge between a classical electromagnetic (EM) simulation and the quantum Hamiltonian of a superconducting circuit. Rather than deriving circuit parameters by hand from a schematic, the EPR method extracts them directly from the full-wave, distributed electromagnetic field solutions computed by a finite-element solver such as Ansys HFSS.
The method is described in detail in:
Z. K. Minev, Z. Leghtas, S. O. Mundhada, L. Christakis, I. M. Pop, and M. H. Devoret,
“Energy-participation quantization of Josephson circuits”,
npj Quantum Information 7, 131 (2021).
Preprint: arXiv:2010.00620
This tutorial walks through the complete analysis pipeline for a single transmon qubit coupled to a 3D microwave cavity. By the end you will have extracted the qubit frequency, anharmonicity, and dispersive coupling from an HFSS eigenmode simulation.
The three-stage EPR pipeline#
Stage 1 — Classical eigenmode simulation (HFSS)
The circuit is simulated in the linearized regime, meaning the Josephson junction (JJ) is replaced by its linear inductance \(L_J\). HFSS solves Maxwell’s equations to find the electromagnetic eigenmodes: resonant frequencies \(\omega_m / 2\pi\) and the associated spatial field distributions \(\mathbf{E}_m(\mathbf{r}),\, \mathbf{H}_m(\mathbf{r})\).
Stage 2 — EPR extraction (DistributedAnalysis)
For each mode \(m\) and each Josephson junction \(j\), pyEPR computes the energy-participation ratio: $\(p_{mj} = \frac{\omega_m}{E_J} \frac{\partial E_J}{\partial (1/L_j)} \cdot \frac{1}{2} = \frac{\text{inductive energy in junction } j \text{ for mode } m}{\text{total inductive energy of mode } m}\)$
In practice \(p_{mj}\) is computed from the zero-point voltage across the junction using a line integral of the electric field along the polyline that spans the junction rectangle: $\(p_{mj} = \frac{\omega_m L_j}{\hbar \omega_m / (2 e^2)} \cdot \frac{V_{zpf,j}^2}{2 \omega_m} \;\approx\; \frac{L_j \left|\int_{\ell_j} \mathbf{E}_m \cdot d\mathbf{\ell}\right|^2}{4 \omega_m^2 \mathcal{E}_m^{(E)}}\)\( where \)\mathcal{E}_m^{(E)}\( is the peak electric energy of mode \)m$. See Eq. (4) and Appendix B of arXiv:2010.00620 for the precise definition.
The zero-point flux fluctuation across junction \(j\) in mode \(m\) is: $\(\varphi_{zpf}^{(mj)} = \sqrt{\frac{p_{mj}\, \hbar \omega_m}{2 E_J}}\)$ This dimensionless quantity quantifies how strongly each mode couples to the nonlinearity of the junction.
Stage 3 — Quantum Hamiltonian analysis (QuantumAnalysis)
The EPR method maps the system to the quantum Hamiltonian (see Eq. (1) of the paper): $\(\hat{H} = \sum_m \hbar\omega_m \hat{a}_m^\dagger \hat{a}_m - \sum_j E_J \left[\cos\left(\sum_m \varphi_{zpf}^{(mj)}(\hat{a}_m + \hat{a}_m^\dagger)\right) - 1 + \frac{\varphi_j^2}{2}\right]\)$ Numerical diagonalization (or perturbation theory) of this Hamiltonian yields:
Dressed mode frequencies \(\tilde{\omega}_m\) (shifted from \(\omega_m\) by the Lamb shift due to junction nonlinearity)
Anharmonicities \(\alpha_m = \chi_{mm}\) (diagonal of the \(\chi\) matrix)
Dispersive couplings \(\chi_{mn}\) (off-diagonal; quantifies qubit-cavity interaction strength)
The three main pyEPR classes#
Class |
Role |
|---|---|
|
Configuration: HFSS project path, junction definitions, solver settings |
|
Runs field extraction in HFSS; computes EPRs (Stage 1 & 2) |
|
Loads EPR results; diagonalizes quantum Hamiltonian (Stage 3) |
Let us now work through each stage step by step.
Load pyEPR#
import pyEPR as epr
This loads the pyEPR package under the shorthand name epr.
We can access the main submodules as follows:
epr.calcs : calculation tools (EPR math, unit conversion, Josephson energy formulas)
epr.ansys : tools for controlling Ansys HFSS, Q3D, etc. via COM
epr.toolbox : useful toolbox for plotting, pandas helpers, and Python utilities
epr.core : the core analysis classes — ProjectInfo, DistributedAnalysis, and QuantumAnalysis
pyEPR also exposes several convenient utility functions at the top level. For example, parse_entry converts human-readable length strings to SI units, and the calcs.Convert submodule handles physical constant conversions such as Josephson energy from junction inductance:
print('Parsing unis: 1um =',
epr.parse_entry('1um', 'meters'), 'meters')
print(f"""For L_J = 11 nH, the Josephson junction energy is
E_J = {epr.calcs.Convert.Ej_from_Lj(11, 'nH', 'GHz'):.1f} GHz""")
The Josephson energy \(E_J\) is related to the junction inductance \(L_J\) through the Josephson relation:
$\(E_J = \frac{\Phi_0^2}{(2\pi)^2 L_J} = \frac{\hbar^2}{4 e^2 L_J}\)\(
where \)\Phi_0 = h / (2e)\( is the magnetic flux quantum. For a typical transmon with \)L_J \approx 11,\text{nH}\(, this gives \)E_J / h \approx 15,\text{GHz}\(. The transmon qubit frequency is approximately \)\omega_q \approx \sqrt{8 E_J E_C} - E_C\( (where \)E_C\( is the charging energy), so knowing \)E_J$ is essential for device design.
Load Ansys HFSS tutorial file#
Let us first determine where the example file is stored. The tutorial HFSS project ships with the pyEPR repository under _example_files/. We use pathlib.Path to construct the path relative to the installed package location.
# Load Path temporarily just to find where the tutorial folder is
# return path_to_project
from pathlib import Path
path_to_project = Path(epr.__file__).parent.parent / '_example_files'
print(f'We will the example project located in\n {path_to_project}')
Connect to Ansys HFSS#
Now, we will open Ansys Desktop and connect to a specific project.
ProjectInfo is the first of the three main pyEPR classes. It manages the connection to Ansys HFSS and stores all the configuration needed for the analysis — including junction definitions, dissipative element lists, and solver settings.
You can use the following creation keyword arguments to specify which design to look at:
project_path: Directory path to the HFSS project file. Should be the directory, not the file itself. None assumes the project is already open and gets the project based on project_name.
project_name: Name of the project within project_path. None will get the current active one.
design_name: Name of the design within the project. None will get the current active one.
setup_name: Name of the setup within the design. None will get the current active one.
Use ?epr.ProjectInfo to see all arguments.
pinfo = epr.ProjectInfo(project_path = path_to_project,
project_name = 'pyEPR_tutorial1',
design_name = '1. single_transmon')
Let’s see what the design type was using the design object. For the EPR analysis we always use an Eigenmode solution type, which finds the resonant frequencies and field distributions of the linearized circuit.
pinfo.design.solution_type
Let us see what all the objects in the design are. pinfo contains a convenience function to retrieve these. In the tutorial design you will find:
cavity_enclosure: the 3D metallic box forming the microwave cavity
substrate: the silicon or sapphire chip on which the qubit is fabricated
rect_jj1: the rectangle sheet that defines the Josephson junction lumped boundary
pad_lower, pad_upper: the capacitor pads of the transmon qubit
line_jj1: the polyline spanning the junction rectangle (defines current direction)
pinfo.get_all_object_names()
Let’s see what design variables are available. HFSS design variables parametrize the geometry and material properties of the model. The variable Lj_1 is of particular importance — it is the Josephson inductance that pyEPR will read from HFSS to compute the EPR.
pinfo.get_all_variables_names()
We can also get the names of the setups.
pinfo.design.get_setup_names()
We can also retrieve information such as mesh statistics and convergence (if the model has already been solved). The mesh quality determines the accuracy of the field integrals used in the EPR calculation. A well-converged eigenmode solution typically requires \(\Delta f / f < 0.5\%\) between successive adaptive passes.
pinfo.setup.get_mesh_stats()
We can also see how many modes we solved for.
pinfo.setup.n_modes
Josephson Tunnel Junctions (Non-Linear Elements) #
Physics of the Josephson junction in HFSS#
In an HFSS eigenmode simulation the Josephson junction is modeled as a lumped RLC boundary condition — a rectangular sheet with a specified sheet inductance \(L_J\). This is the correct linearized representation because, for small oscillations around the superconducting phase minimum, the JJ potential energy \(-E_J \cos(\varphi) \approx E_J \varphi^2 / 2 - E_J \varphi^4 / 24 + \ldots\) is dominated by the harmonic (quadratic) term, which corresponds exactly to an inductor with inductance:
$\(L_J = \frac{\Phi_0^2}{(2\pi)^2 E_J} = \frac{\hbar^2}{4 e^2 E_J}\)$
The anharmonicity — the deviation from a perfect harmonic oscillator that makes the qubit addressable — comes from the quartic and higher-order terms in the cosine expansion. The EPR method captures these nonlinear corrections in Stage 3 (quantum Hamiltonian analysis), after extracting the linear mode structure in Stages 1 and 2.
Each junction in the model is described by four pieces of information:
Lj_variable: the name of the HFSS design variable that specifies \(L_J\) on the lumped boundary condition. pyEPR reads this variable to compute \(E_J = \Phi_0^2 / (4\pi^2 L_J)\).
rect: the name of the HFSS rectangle on which the lumped boundary condition is applied. The area of this rectangle determines the current density normalization.
line: the name of the HFSS polyline that spans the length of the rectangle, parallel to the current flow direction. The zero-point voltage is computed by integrating the electric field \(\mathbf{E}\) along this line.
length: the physical length of the junction (used for normalization in the EPR line-voltage calculation).
An optional Cj_variable can be specified if a lumped capacitance is placed in parallel with the junction in the HFSS model.
To analyze the quantum part of the model, let us now register the junctions present in the design. We call the single junction j1:
pinfo.junctions['j1'] = {'Lj_variable' : 'Lj_1',
'rect' : 'rect_jj1',
'line' : 'line_jj1',
'length' : epr.parse_units('100um')}
# Check that valid names of variables and objects have been supplied.
# An error is raised with a message if something is wrong.
pinfo.validate_junction_info()
Microwave Analysis:
Run analysis on an eigenmode solution #
Assuming we have field solutions saved for the HFSS eigenmode design, we can now create the analysis object from the class DistributedAnalysis.
Run analysis (tutorial action): Run the analysis in the design model. This can be done manually in the HFSS GUI, or programmatically using pyEPR.
# Calling this function will run an analysis of the design in the background.
# The console will wait for HFSS to complete.
pinfo.setup.analyze()
Run optimetrics
If an optimetrics (parametric sweep) setup is defined in HFSS, the following call will execute it. Each parametric point generates a separate set of field solutions that pyEPR can subsequently analyze.
The following will get the frequencies and Qs from HFSS.
# Calling this function will run an analysis on all the optimetrics the design in the background.
# The console will wait for HFSS to complete.
# This requires optimetrics license by HFSS.
pinfo.design.optimetrics.solve_setup(pinfo.design.optimetrics.get_setup_names()[0])
EPR Analysis object
Let us introduce the core analysis object of pyEPR, which will be a class instance of epr.DistributedAnalysis. This is the classical-to-quantum bridge: it connects to HFSS, retrieves the electromagnetic field solutions, and computes the energy-participation ratios for each mode and junction.
Let us recover the eigenmode frequencies of the first variation we ran, which is indexed by '0'.
# This is the core object for interacting with HFSS
# and running analysis within HFSS.
eprh = epr.DistributedAnalysis(pinfo) # epr hfss analysis
eprh.get_ansys_variables()
eprh.get_ansys_frequencies_all()
eprh.hfss_report_full_convergence();
We can see that the convergence is reasonable for this tutorial example. In production designs, one typically targets \(\Delta f / f < 0.1\%\) between adaptive mesh passes and a monotonically decreasing convergence curve. Mesh refinement operations (e.g., seeding finer elements around the junction rectangle or the qubit capacitor pads) can improve convergence.
For this tutorial, let us press forward.
You can use these autogenerated reports to keep track of your simulations. This is still all at the classical level. Take a look at the other available methods:
hfss_report_full_convergence
hfss_report_f_convergence
get_convergence
get_mesh_statistics
get_freqs_bare_pd
…
For example, you can get tables of the frequencies and Q factors using:
Fs, Qs = eprh.get_freqs_bare_pd(variation='0')
import pandas as pd
for variation in eprh.variations[:2]: # just for the first 2
Fs, Qs = eprh.get_freqs_bare_pd(variation=variation, frame=False)
display(pd.DataFrame({'Freq. (GHz)':Fs, 'Quality Factor':Qs}))
The Q is infinite, since we have not included dissipation in this example — there are no lossy boundaries or bulk dielectrics defined in the HFSS model. The row index is the mode number:
Mode 0 (~5 GHz): the transmon qubit mode
Mode 1 (~9 GHz): the cavity mode
Dissipation can be modeled by specifying the corresponding HFSS objects for bulk dielectrics, surface losses, resistive surfaces, or seams when initializing ProjectInfo. For example:
pinfo = epr.ProjectInfo(project_path = path_to_project,
project_name = 'pyEPR_tutorial1',
design_name = '1. single_transmon',
dielectrics_bulk = ['si_substrate'],
dielectric_surfaces = ['interface'],
resistive_surfaces = None,
seams = None)
or equivalently, after the fact:
pinfo.dissipative['dielectrics_bulk'] = ['si_substrate']
pinfo.dissipative['dielectric_surfaces'] = ['interface']
See Tutorial 2 for a detailed treatment of dielectric loss calculation.
Full EPR Analysis#
What does do_EPR_analysis compute?#
This is the heart of the pyEPR workflow — the classical-to-quantum bridge. For each eigenmode and each parametric variation (e.g., each value of \(L_J\) in a sweep), do_EPR_analysis performs the following steps:
Field retrieval: Queries HFSS for the electric and magnetic field distributions \(\mathbf{E}_m(\mathbf{r})\), \(\mathbf{H}_m(\mathbf{r})\) of mode \(m\).
Energy balance check: Computes the total electric energy \(\mathcal{E}_m^{(E)}\) and magnetic energy \(\mathcal{E}_m^{(H)}\) by volume integration over the entire simulation domain:
$\(\mathcal{E}_m^{(E)} = \frac{1}{4}\mathrm{Re}\int_V dv\, \mathbf{E}_m^* \cdot \overleftrightarrow{\epsilon} \cdot \mathbf{E}_m\)\(
For a well-converged eigenmode, \)\mathcal{E}_m^{(E)} \approx \mathcal{E}_m^{(H)}$ (energy equipartition). A large imbalance signals poor mesh convergence.
Junction EPR calculation: For each junction \(j\), computes the participation ratio \(p_{mj}\) by integrating the electric field along the polyline \(\ell_j\):
$\(p_{mj} = \frac{\omega_m L_j \left| \int_{\ell_j} \mathbf{E}_m \cdot d\mathbf{\ell} \right|^2}{4 \omega_m^2 \mathcal{E}_m^{(E)}}\)$
See Eq. (4) and Appendix B of arXiv:2010.00620.
Sign determination: Computes the sign \(s_{mj} = \pm 1\) of the zero-point flux fluctuation, which determines whether modes are hybridized constructively or destructively at the junction.
Data storage: Saves all results (frequencies, participation matrices, sign matrices, and energy values) to an HDF5 file for subsequent quantum analysis.
The output printed during the run shows, for each mode, the energy balance \((\mathcal{E}_E - \mathcal{E}_H)/\mathcal{E}_E\) and the computed EPR \(p_{mj}\) for each junction.
# Do EPR analysis for all modes and all variations
# If you want to only analze a certain variation, you can change the
# default keyword argeumetns variations=None, modes=None
eprh.do_EPR_analysis();
Getting more info#
In general, to see more verbose output when running commands, set the logger level to DEBUG.
import logging
epr.logger.setLevel(logging.DEBUG)
Quantum Hamiltonian Analysis:
Load microwave analysis results and perform Hamiltonian analysis on the solutions #
The quantum Hamiltonian and its diagonalization#
Having computed the EPR participation matrix \(p_{mj}\) and zero-point flux fluctuations \(\varphi_{zpf}^{(mj)}\) for each mode–junction pair, we can now construct and diagonalize the full quantum Hamiltonian. This is done by the QuantumAnalysis class.
The system Hamiltonian (in the basis of the linearized HFSS eigenmodes) is (Eq. (1) of arXiv:2010.00620):
$\(\hat{H} = \sum_m \hbar\omega_m \hat{a}_m^\dagger \hat{a}_m - \sum_j E_J^{(j)} \left[\cos\!\left(\hat{\varphi}_j\right) - 1 + \frac{\hat{\varphi}_j^2}{2}\right]\)\(
where the phase operator for junction \)j\( is:
\)\(\hat{\varphi}_j = \sum_m s_{mj}\, \varphi_{zpf}^{(mj)} \left(\hat{a}_m + \hat{a}_m^\dagger\right)\)$
The cosine is expanded as a Taylor series (truncated at order cos_trunc) and the resulting polynomial in bosonic operators is diagonalized numerically in a truncated Fock space of dimension fock_trunc per mode. Two levels of analysis are provided:
First-order perturbation theory (PT): Analytic result valid when \(\varphi_{zpf} \ll 1\) (transmon regime). Fast, but less accurate for strongly anharmonic or hybridized modes.
Numerical diagonalization (ND): Full numerical diagonalization in the Fock space. More accurate, but computationally more expensive.
Key parameters:
cos_trunc: Number of terms in the Taylor expansion of \(\cos(\hat{\varphi})\) kept in the Hamiltonian. Typical value: 8. Higher values improve accuracy at the cost of computation time.
fock_trunc: Dimension of the truncated Fock space per mode. Typical value: 15 (i.e., states \(|0\rangle, |1\rangle, \ldots, |14\rangle\)). Must be large enough that truncation errors are negligible.
%matplotlib inline
# plot the figures in the notebook.
# This is the core object for working with post processing data from HFSS
# This iwll load the saved solutions from above
epra = epr.QuantumAnalysis(eprh.data_filename)
# Analyze
epra.analyze_all_variations(cos_trunc = 8, fock_trunc = 15)
epra.plot_hamiltonian_results();
Interpreting the Hamiltonian results#
The plot_hamiltonian_results() output and the printed tables give the following quantum quantities:
Chi matrix (\(\chi\)):
The \(\chi\) matrix encodes all second-order quantum corrections from the junction nonlinearity. In the standard dispersive approximation it takes the form \(\hat{H} \supset \sum_{m,n} \chi_{mn} \hat{a}_m^\dagger \hat{a}_m \hat{a}_n^\dagger \hat{a}_n\).
Diagonal elements \(\chi_{mm}\) (labeled anharmonicity \(\alpha_m\)): The self-Kerr coefficient, which sets how much the frequency of mode \(m\) shifts for each additional photon in that mode. For the transmon, \(\alpha_q \approx -E_C/h\) (where \(E_C = e^2 / (2C_\Sigma)\) is the charging energy). A transmon with \(\alpha \approx 200\,\text{MHz}\) can be selectively addressed with \(\sim 10\,\text{ns}\) pulses.
Off-diagonal elements \(\chi_{mn}\) (\(m \ne n\)): The cross-Kerr (dispersive coupling) between modes \(m\) and \(n\). In a qubit-cavity system, \(\chi_{qc}\) is the dispersive shift used for qubit readout: the cavity resonance shifts by \(\chi_{qc}\) depending on whether the qubit is in \(|0\rangle\) or \(|1\rangle\).
Dressed frequencies vs. bare frequencies:
The bare frequencies \(\omega_m\) are the HFSS eigenfrequencies of the linearized circuit. The dressed frequencies \(\tilde{\omega}_m\) include the Lamb shift — a downward frequency renormalization due to the quantum zero-point fluctuations coupling each mode to the junction nonlinearity. This shift is \(\tilde{\omega}_m = \omega_m + \chi_{mm}/2 + \ldots\) to leading order in perturbation theory.
The two columns labeled O1 PT and ND give the first-order perturbation theory and numerical diagonalization results respectively. For a typical transmon with small \(\varphi_{zpf}\), these agree to within a few percent.
See Section III and Table I of arXiv:2010.00620 for a systematic comparison of these methods against experiment.
Swept variable analysis driven through pyEPR:
Sweep variable in HFSS, analyze microwaves, then analyze quantum#
Let us now run a more complete example, where we sweep the Josephson inductance \(L_J\). This demonstrates the full power of the pyEPR workflow: by sweeping \(L_J\) (which is equivalent to varying the Josephson energy \(E_J = \Phi_0^2 / (4\pi^2 L_J)\)), we can trace out how the qubit frequency, anharmonicity, and dispersive coupling all depend on the junction parameter.
This is exactly how one would characterize the design space during the development of a new qubit device.
Classical microwave, distributed analysis
Sweep variable in HFSS, run the eigenmode solver for each value.
import numpy as np
### Configure setup
### Set the maximum number of passes on the analysis and a few other
pinfo.setup.passes = 15
pinfo.setup.min_freq = '3.2GHz'
pinfo.setup.modes = 2
pinfo.setup.delta_f = 0.02
pinfo.setup.basis_order = epr.ansys.BASIS_ORDER['Mixed Order']
### What to sweep
# Define sweep variables to go over
swp_var = 'Lj_1' # name
# how to convert sweep param to string variable in nanoHeries.
to_swp_val = lambda x: f'{x}nH'
### Run Sweep
for swp_param in np.linspace(6, 12, 6):
swp_val = to_swp_val(swp_param)
epr.logger.info(f'Setting sweep variable {swp_var}={swp_val}')
pinfo.design.set_variable(swp_var, swp_val)
pinfo.setup.analyze()
EPR microwave analysis
Classical-to-quantum bridge: compute participation ratios for all sweep points.
eprh = epr.DistributedAnalysis(pinfo) # epr hfss analysis
eprh.hfss_report_full_convergence()
# Do EPR analysis for all modes and all variations
eprh.do_EPR_analysis();
Quantum analysis
Perform the Hamiltonian diagonalization and plot the results as a function of the swept parameter \(L_J\).
# This is the core object for working with post processing data from HFSS
# This will load the saved solutions from above
epra = epr.QuantumAnalysis(eprh.data_filename)
# Analyze
epra.analyze_all_variations(cos_trunc = 8, fock_trunc = 15);
epra.plot_hamiltonian_results(swp_variable='Lj_1');
There is some jitter in the results since some of the runs were performed with only a few adaptive passes. In a production analysis one would ensure all variations have converged to the same \(\Delta f\) threshold before proceeding.
The sweep makes several physically important trends visible:
As \(L_J\) increases, \(E_J\) decreases, so the qubit frequency \(\omega_q \propto \sqrt{E_J}\) decreases.
The anharmonicity \(\alpha \approx -E_C\) is relatively insensitive to \(L_J\) in the transmon regime (since \(E_C\) is set by the capacitor geometry, not the junction).
The dispersive coupling \(\chi_{qc}\) decreases as \(L_J\) increases (equivalently, as \(E_J / E_C\) increases and the qubit moves deeper into the transmon regime).
See Figure 3 of arXiv:2010.00620 for an example of these trends in a real device.
epra.quick_plot_mode(0,0,1,numeric=True,swp_variable='Lj_1')
epra.report_results(swp_variable='Lj_1', numeric=True)
Addendum:
Snippets for ease #
An easy way to define many junctions is the following, using Python 3 f-strings for ease. This pattern is useful for multi-junction circuits such as flux qubits, SQUID arrays, or fluxonium qubits, where each junction must be registered individually but shares a common naming convention:
# Specify multiple junctions in HFSS model
n_junctions = 5
for i in range(1, 1+n_junctions):
pinfo.junctions[f'j{i}'] = {'Lj_variable' : f'Lj{i}',
'rect' : f'JJrect{i}',
'line' : f'JJline{i}',
'length' : epr.parse_units('100um')}
NEXT: Please see Tutorial 2 to continue. It covers the calculation of dielectric energy-participation ratios, which are needed to estimate qubit lifetimes \(T_1\) from dielectric loss tangents.
import pyEPR as epr
This loads the pyEPR package under the shorthand name epr.
We can access the main submodules as follows:
epr.calcs: calculation tools (EPR math, unit conversion, Josephson energy formulas)epr.ansys: tools for controlling Ansys HFSS, Q3D, etc. via COMepr.toolbox: useful toolbox for plotting, pandas helpers, and Python utilitiesepr.core: the core analysis classes —ProjectInfo,DistributedAnalysis, andQuantumAnalysis
pyEPR also exposes several convenient utility functions at the top level. For example, parse_entry converts human-readable length strings to SI units, and the calcs.Convert submodule handles physical constant conversions such as Josephson energy from junction inductance:
print('Parsing unis: 1um =',
epr.parse_entry('1um', 'meters'), 'meters')
print(f"""For L_J = 11 nH, the Josephson junction energy is
E_J = {epr.calcs.Convert.Ej_from_Lj(11, 'nH', 'GHz'):.1f} GHz""")
The Josephson energy \(E_J\) is related to the junction inductance \(L_J\) through the Josephson relation: $\(E_J = \frac{\Phi_0^2}{(2\pi)^2 L_J} = \frac{\hbar^2}{4 e^2 L_J}\)\( where \)\Phi_0 = h / (2e)\( is the magnetic flux quantum. For a typical transmon with \)L_J \approx 11,\text{nH}\(, this gives \)E_J / h \approx 15,\text{GHz}\(. The transmon qubit frequency is approximately \)\omega_q \approx \sqrt{8 E_J E_C} - E_C\( (where \)E_C\( is the charging energy), so knowing \)E_J$ is essential for device design.
Load Ansys HFSS tutorial file#
Let us first determine where the example file is stored. The tutorial HFSS project ships with the pyEPR repository under _example_files/. We use pathlib.Path to construct the path relative to the installed package location.
# Load Path temporarily just to find where the tutorial folder is
# return path_to_project
from pathlib import Path
path_to_project = Path(epr.__file__).parent.parent / '_example_files'
print(f'We will the example project located in\n {path_to_project}')
Connect to Ansys HFSS#
Now, we will open Ansys Desktop and connect to a specific project.
ProjectInfo is the first of the three main pyEPR classes. It manages the connection to Ansys HFSS and stores all the configuration needed for the analysis — including junction definitions, dissipative element lists, and solver settings.
You can use the following creation keyword arguments to specify which design to look at:
project_path: Directory path to the HFSS project file. Should be the directory, not the file itself. None assumes the project is already open and gets the project based on project_name.
project_name: Name of the project within project_path. None will get the current active one.
design_name: Name of the design within the project. None will get the current active one.
setup_name: Name of the setup within the design. None will get the current active one.
Use ?epr.ProjectInfo to see all arguments.
pinfo = epr.ProjectInfo(project_path = path_to_project,
project_name = 'pyEPR_tutorial1',
design_name = '1. single_transmon')
Let’s see what the design type was using the design object. For the EPR analysis we always use an Eigenmode solution type, which finds the resonant frequencies and field distributions of the linearized circuit.
pinfo.design.solution_type
Let us see what all the objects in the design are. pinfo contains a convenience function to retrieve these. In the tutorial design you will find:
cavity_enclosure: the 3D metallic box forming the microwave cavity
substrate: the silicon or sapphire chip on which the qubit is fabricated
rect_jj1: the rectangle sheet that defines the Josephson junction lumped boundary
pad_lower, pad_upper: the capacitor pads of the transmon qubit
line_jj1: the polyline spanning the junction rectangle (defines current direction)
pinfo.get_all_object_names()
Let’s see what design variables are available. HFSS design variables parametrize the geometry and material properties of the model. The variable Lj_1 is of particular importance — it is the Josephson inductance that pyEPR will read from HFSS to compute the EPR.
pinfo.get_all_variables_names()
We can also get the names of the setups.
pinfo.design.get_setup_names()
We can also retrieve information such as mesh statistics and convergence (if the model has already been solved). The mesh quality determines the accuracy of the field integrals used in the EPR calculation. A well-converged eigenmode solution typically requires \(\Delta f / f < 0.5\%\) between successive adaptive passes.
pinfo.setup.get_mesh_stats()
We can also see how many modes we solved for.
pinfo.setup.n_modes
Josephson Tunnel Junctions (Non-Linear Elements) #
Physics of the Josephson junction in HFSS#
In an HFSS eigenmode simulation the Josephson junction is modeled as a lumped RLC boundary condition — a rectangular sheet with a specified sheet inductance \(L_J\). This is the correct linearized representation because, for small oscillations around the superconducting phase minimum, the JJ potential energy \(-E_J \cos(\varphi) \approx E_J \varphi^2 / 2 - E_J \varphi^4 / 24 + \ldots\) is dominated by the harmonic (quadratic) term, which corresponds exactly to an inductor with inductance:
$\(L_J = \frac{\Phi_0^2}{(2\pi)^2 E_J} = \frac{\hbar^2}{4 e^2 E_J}\)$
The anharmonicity — the deviation from a perfect harmonic oscillator that makes the qubit addressable — comes from the quartic and higher-order terms in the cosine expansion. The EPR method captures these nonlinear corrections in Stage 3 (quantum Hamiltonian analysis), after extracting the linear mode structure in Stages 1 and 2.
Each junction in the model is described by four pieces of information:
Lj_variable: the name of the HFSS design variable that specifies \(L_J\) on the lumped boundary condition. pyEPR reads this variable to compute \(E_J = \Phi_0^2 / (4\pi^2 L_J)\).
rect: the name of the HFSS rectangle on which the lumped boundary condition is applied. The area of this rectangle determines the current density normalization.
line: the name of the HFSS polyline that spans the length of the rectangle, parallel to the current flow direction. The zero-point voltage is computed by integrating the electric field \(\mathbf{E}\) along this line.
length: the physical length of the junction (used for normalization in the EPR line-voltage calculation).
An optional Cj_variable can be specified if a lumped capacitance is placed in parallel with the junction in the HFSS model.
To analyze the quantum part of the model, let us now register the junctions present in the design. We call the single junction j1:
pinfo.junctions['j1'] = {'Lj_variable' : 'Lj_1',
'rect' : 'rect_jj1',
'line' : 'line_jj1',
'length' : epr.parse_units('100um')}
# Check that valid names of variables and objects have been supplied.
# An error is raised with a message if something is wrong.
pinfo.validate_junction_info()
Microwave Analysis:
Run analysis on an eigenmode solution #
Assuming we have field solutions saved for the HFSS eigenmode design, we can now create the analysis object from the class DistributedAnalysis.
Run analysis (tutorial action): Run the analysis in the design model. This can be done manually in the HFSS GUI, or programmatically using pyEPR.
# Calling this function will run an analysis of the design in the background.
# The console will wait for HFSS to complete.
pinfo.setup.analyze()
Run optimetrics
If an optimetrics (parametric sweep) setup is defined in HFSS, the following call will execute it. Each parametric point generates a separate set of field solutions that pyEPR can subsequently analyze.
The following will get the frequencies and Qs from HFSS.
# Calling this function will run an analysis on all the optimetrics the design in the background.
# The console will wait for HFSS to complete.
# This requires optimetrics license by HFSS.
pinfo.design.optimetrics.solve_setup(pinfo.design.optimetrics.get_setup_names()[0])
EPR Analysis object
Let us introduce the core analysis object of pyEPR, which will be a class instance of epr.DistributedAnalysis. This is the classical-to-quantum bridge: it connects to HFSS, retrieves the electromagnetic field solutions, and computes the energy-participation ratios for each mode and junction.
Let us recover the eigenmode frequencies of the first variation we ran, which is indexed by '0'.
# This is the core object for interacting with HFSS
# and running analysis within HFSS.
eprh = epr.DistributedAnalysis(pinfo) # epr hfss analysis
eprh.get_ansys_variables()
eprh.get_ansys_frequencies_all()
eprh.hfss_report_full_convergence();
We can see that the convergence is reasonable for this tutorial example. In production designs, one typically targets \(\Delta f / f < 0.1\%\) between adaptive mesh passes and a monotonically decreasing convergence curve. Mesh refinement operations (e.g., seeding finer elements around the junction rectangle or the qubit capacitor pads) can improve convergence.
For this tutorial, let us press forward.
You can use these autogenerated reports to keep track of your simulations. This is still all at the classical level. Take a look at the other available methods:
hfss_report_full_convergence
hfss_report_f_convergence
get_convergence
get_mesh_statistics
get_freqs_bare_pd
…
For example, you can get tables of the frequencies and Q factors using:
Fs, Qs = eprh.get_freqs_bare_pd(variation='0')
import pandas as pd
for variation in eprh.variations[:2]: # just for the first 2
Fs, Qs = eprh.get_freqs_bare_pd(variation=variation, frame=False)
display(pd.DataFrame({'Freq. (GHz)':Fs, 'Quality Factor':Qs}))
The Q is infinite, since we have not included dissipation in this example — there are no lossy boundaries or bulk dielectrics defined in the HFSS model. The row index is the mode number:
Mode 0 (~5 GHz): the transmon qubit mode
Mode 1 (~9 GHz): the cavity mode
Dissipation can be modeled by specifying the corresponding HFSS objects for bulk dielectrics, surface losses, resistive surfaces, or seams when initializing ProjectInfo. For example:
pinfo = epr.ProjectInfo(project_path = path_to_project,
project_name = 'pyEPR_tutorial1',
design_name = '1. single_transmon',
dielectrics_bulk = ['si_substrate'],
dielectric_surfaces = ['interface'],
resistive_surfaces = None,
seams = None)
or equivalently, after the fact:
pinfo.dissipative['dielectrics_bulk'] = ['si_substrate']
pinfo.dissipative['dielectric_surfaces'] = ['interface']
See Tutorial 2 for a detailed treatment of dielectric loss calculation.
Full EPR Analysis#
What does do_EPR_analysis compute?#
This is the heart of the pyEPR workflow — the classical-to-quantum bridge. For each eigenmode and each parametric variation (e.g., each value of \(L_J\) in a sweep), do_EPR_analysis performs the following steps:
Field retrieval: Queries HFSS for the electric and magnetic field distributions \(\mathbf{E}_m(\mathbf{r})\), \(\mathbf{H}_m(\mathbf{r})\) of mode \(m\).
Energy balance check: Computes the total electric energy \(\mathcal{E}_m^{(E)}\) and magnetic energy \(\mathcal{E}_m^{(H)}\) by volume integration over the entire simulation domain:
$\(\mathcal{E}_m^{(E)} = \frac{1}{4}\mathrm{Re}\int_V dv\, \mathbf{E}_m^* \cdot \overleftrightarrow{\epsilon} \cdot \mathbf{E}_m\)\(
For a well-converged eigenmode, \)\mathcal{E}_m^{(E)} \approx \mathcal{E}_m^{(H)}$ (energy equipartition). A large imbalance signals poor mesh convergence.
Junction EPR calculation: For each junction \(j\), computes the participation ratio \(p_{mj}\) by integrating the electric field along the polyline \(\ell_j\):
$\(p_{mj} = \frac{\omega_m L_j \left| \int_{\ell_j} \mathbf{E}_m \cdot d\mathbf{\ell} \right|^2}{4 \omega_m^2 \mathcal{E}_m^{(E)}}\)$
See Eq. (4) and Appendix B of arXiv:2010.00620.
Sign determination: Computes the sign \(s_{mj} = \pm 1\) of the zero-point flux fluctuation, which determines whether modes are hybridized constructively or destructively at the junction.
Data storage: Saves all results (frequencies, participation matrices, sign matrices, and energy values) to an HDF5 file for subsequent quantum analysis.
The output printed during the run shows, for each mode, the energy balance \((\mathcal{E}_E - \mathcal{E}_H)/\mathcal{E}_E\) and the computed EPR \(p_{mj}\) for each junction.
# Do EPR analysis for all modes and all variations
# If you want to only analze a certain variation, you can change the
# default keyword argeumetns variations=None, modes=None
eprh.do_EPR_analysis();
Getting more info#
In general, to see more verbose output when running commands, set the logger level to DEBUG.
import logging
epr.logger.setLevel(logging.DEBUG)
Quantum Hamiltonian Analysis:
Load microwave analysis results and perform Hamiltonian analysis on the solutions #
The quantum Hamiltonian and its diagonalization#
Having computed the EPR participation matrix \(p_{mj}\) and zero-point flux fluctuations \(\varphi_{zpf}^{(mj)}\) for each mode–junction pair, we can now construct and diagonalize the full quantum Hamiltonian. This is done by the QuantumAnalysis class.
The system Hamiltonian (in the basis of the linearized HFSS eigenmodes) is (Eq. (1) of arXiv:2010.00620):
$\(\hat{H} = \sum_m \hbar\omega_m \hat{a}_m^\dagger \hat{a}_m - \sum_j E_J^{(j)} \left[\cos\!\left(\hat{\varphi}_j\right) - 1 + \frac{\hat{\varphi}_j^2}{2}\right]\)\(
where the phase operator for junction \)j\( is:
\)\(\hat{\varphi}_j = \sum_m s_{mj}\, \varphi_{zpf}^{(mj)} \left(\hat{a}_m + \hat{a}_m^\dagger\right)\)$
The cosine is expanded as a Taylor series (truncated at order cos_trunc) and the resulting polynomial in bosonic operators is diagonalized numerically in a truncated Fock space of dimension fock_trunc per mode. Two levels of analysis are provided:
First-order perturbation theory (PT): Analytic result valid when \(\varphi_{zpf} \ll 1\) (transmon regime). Fast, but less accurate for strongly anharmonic or hybridized modes.
Numerical diagonalization (ND): Full numerical diagonalization in the Fock space. More accurate, but computationally more expensive.
Key parameters:
cos_trunc: Number of terms in the Taylor expansion of \(\cos(\hat{\varphi})\) kept in the Hamiltonian. Typical value: 8. Higher values improve accuracy at the cost of computation time.
fock_trunc: Dimension of the truncated Fock space per mode. Typical value: 15 (i.e., states \(|0\rangle, |1\rangle, \ldots, |14\rangle\)). Must be large enough that truncation errors are negligible.
%matplotlib inline
# plot the figures in the notebook.
# This is the core object for working with post processing data from HFSS
# This iwll load the saved solutions from above
epra = epr.QuantumAnalysis(eprh.data_filename)
# Analyze
epra.analyze_all_variations(cos_trunc = 8, fock_trunc = 15)
epra.plot_hamiltonian_results();
Interpreting the Hamiltonian results#
The plot_hamiltonian_results() output and the printed tables give the following quantum quantities:
Chi matrix (\(\chi\)):
The \(\chi\) matrix encodes all second-order quantum corrections from the junction nonlinearity. In the standard dispersive approximation it takes the form \(\hat{H} \supset \sum_{m,n} \chi_{mn} \hat{a}_m^\dagger \hat{a}_m \hat{a}_n^\dagger \hat{a}_n\).
Diagonal elements \(\chi_{mm}\) (labeled anharmonicity \(\alpha_m\)): The self-Kerr coefficient, which sets how much the frequency of mode \(m\) shifts for each additional photon in that mode. For the transmon, \(\alpha_q \approx -E_C/h\) (where \(E_C = e^2 / (2C_\Sigma)\) is the charging energy). A transmon with \(\alpha \approx 200\,\text{MHz}\) can be selectively addressed with \(\sim 10\,\text{ns}\) pulses.
Off-diagonal elements \(\chi_{mn}\) (\(m \ne n\)): The cross-Kerr (dispersive coupling) between modes \(m\) and \(n\). In a qubit-cavity system, \(\chi_{qc}\) is the dispersive shift used for qubit readout: the cavity resonance shifts by \(\chi_{qc}\) depending on whether the qubit is in \(|0\rangle\) or \(|1\rangle\).
Dressed frequencies vs. bare frequencies:
The bare frequencies \(\omega_m\) are the HFSS eigenfrequencies of the linearized circuit. The dressed frequencies \(\tilde{\omega}_m\) include the Lamb shift — a downward frequency renormalization due to the quantum zero-point fluctuations coupling each mode to the junction nonlinearity. This shift is \(\tilde{\omega}_m = \omega_m + \chi_{mm}/2 + \ldots\) to leading order in perturbation theory.
The two columns labeled O1 PT and ND give the first-order perturbation theory and numerical diagonalization results respectively. For a typical transmon with small \(\varphi_{zpf}\), these agree to within a few percent.
See Section III and Table I of arXiv:2010.00620 for a systematic comparison of these methods against experiment.
Swept variable analysis driven through pyEPR:
Sweep variable in HFSS, analyze microwaves, then analyze quantum#
Let us now run a more complete example, where we sweep the Josephson inductance \(L_J\). This demonstrates the full power of the pyEPR workflow: by sweeping \(L_J\) (which is equivalent to varying the Josephson energy \(E_J = \Phi_0^2 / (4\pi^2 L_J)\)), we can trace out how the qubit frequency, anharmonicity, and dispersive coupling all depend on the junction parameter.
This is exactly how one would characterize the design space during the development of a new qubit device.
Classical microwave, distributed analysis
Sweep variable in HFSS, run the eigenmode solver for each value.
import numpy as np
### Configure setup
### Set the maximum number of passes on the analysis and a few other
pinfo.setup.passes = 15
pinfo.setup.min_freq = '3.2GHz'
pinfo.setup.modes = 2
pinfo.setup.delta_f = 0.02
pinfo.setup.basis_order = epr.ansys.BASIS_ORDER['Mixed Order']
### What to sweep
# Define sweep variables to go over
swp_var = 'Lj_1' # name
# how to convert sweep param to string variable in nanoHeries.
to_swp_val = lambda x: f'{x}nH'
### Run Sweep
for swp_param in np.linspace(6, 12, 6):
swp_val = to_swp_val(swp_param)
epr.logger.info(f'Setting sweep variable {swp_var}={swp_val}')
pinfo.design.set_variable(swp_var, swp_val)
pinfo.setup.analyze()
EPR microwave analysis
Classical-to-quantum bridge: compute participation ratios for all sweep points.
eprh = epr.DistributedAnalysis(pinfo) # epr hfss analysis
eprh.hfss_report_full_convergence()
# Do EPR analysis for all modes and all variations
eprh.do_EPR_analysis();
Quantum analysis
Perform the Hamiltonian diagonalization and plot the results as a function of the swept parameter \(L_J\).
# This is the core object for working with post processing data from HFSS
# This will load the saved solutions from above
epra = epr.QuantumAnalysis(eprh.data_filename)
# Analyze
epra.analyze_all_variations(cos_trunc = 8, fock_trunc = 15);
epra.plot_hamiltonian_results(swp_variable='Lj_1');
There is some jitter in the results since some of the runs were performed with only a few adaptive passes. In a production analysis one would ensure all variations have converged to the same \(\Delta f\) threshold before proceeding.
The sweep makes several physically important trends visible:
As \(L_J\) increases, \(E_J\) decreases, so the qubit frequency \(\omega_q \propto \sqrt{E_J}\) decreases.
The anharmonicity \(\alpha \approx -E_C\) is relatively insensitive to \(L_J\) in the transmon regime (since \(E_C\) is set by the capacitor geometry, not the junction).
The dispersive coupling \(\chi_{qc}\) decreases as \(L_J\) increases (equivalently, as \(E_J / E_C\) increases and the qubit moves deeper into the transmon regime).
See Figure 3 of arXiv:2010.00620 for an example of these trends in a real device.
epra.quick_plot_mode(0,0,1,numeric=True,swp_variable='Lj_1')
epra.report_results(swp_variable='Lj_1', numeric=True)
Addendum:
Snippets for ease #
An easy way to define many junctions is the following, using Python 3 f-strings for ease. This pattern is useful for multi-junction circuits such as flux qubits, SQUID arrays, or fluxonium qubits, where each junction must be registered individually but shares a common naming convention:
# Specify multiple junctions in HFSS model
n_junctions = 5
for i in range(1, 1+n_junctions):
pinfo.junctions[f'j{i}'] = {'Lj_variable' : f'Lj{i}',
'rect' : f'JJrect{i}',
'line' : f'JJline{i}',
'length' : epr.parse_units('100um')}
NEXT: Please see Tutorial 2 to continue. It covers the calculation of dielectric energy-participation ratios, which are needed to estimate qubit lifetimes \(T_1\) from dielectric loss tangents.
Now, we will open Ansys Desktop and connect to a specific project.
ProjectInfo is the first of the three main pyEPR classes. It manages the connection to Ansys HFSS and stores all the configuration needed for the analysis — including junction definitions, dissipative element lists, and solver settings.
You can use the following creation keyword arguments to specify which design to look at:
project_path: Directory path to the HFSS project file. Should be the directory, not the file itself.Noneassumes the project is already open and gets the project based onproject_name.project_name: Name of the project withinproject_path.Nonewill get the current active one.design_name: Name of the design within the project.Nonewill get the current active one.setup_name: Name of the setup within the design.Nonewill get the current active one.
Use ?epr.ProjectInfo to see all arguments.
pinfo = epr.ProjectInfo(project_path = path_to_project,
project_name = 'pyEPR_tutorial1',
design_name = '1. single_transmon')
Let’s see what the design type was using the design object. For the EPR analysis we always use an Eigenmode solution type, which finds the resonant frequencies and field distributions of the linearized circuit.
pinfo.design.solution_type
Let us see what all the objects in the design are. pinfo contains a convenience function to retrieve these. In the tutorial design you will find:
cavity_enclosure: the 3D metallic box forming the microwave cavitysubstrate: the silicon or sapphire chip on which the qubit is fabricatedrect_jj1: the rectangle sheet that defines the Josephson junction lumped boundarypad_lower,pad_upper: the capacitor pads of the transmon qubitline_jj1: the polyline spanning the junction rectangle (defines current direction)
pinfo.get_all_object_names()
Let’s see what design variables are available. HFSS design variables parametrize the geometry and material properties of the model. The variable Lj_1 is of particular importance — it is the Josephson inductance that pyEPR will read from HFSS to compute the EPR.
pinfo.get_all_variables_names()
We can also get the names of the setups.
pinfo.design.get_setup_names()
We can also retrieve information such as mesh statistics and convergence (if the model has already been solved). The mesh quality determines the accuracy of the field integrals used in the EPR calculation. A well-converged eigenmode solution typically requires \(\Delta f / f < 0.5\%\) between successive adaptive passes.
pinfo.setup.get_mesh_stats()
We can also see how many modes we solved for.
pinfo.setup.n_modes
Josephson Tunnel Junctions (Non-Linear Elements) #
Physics of the Josephson junction in HFSS#
In an HFSS eigenmode simulation the Josephson junction is modeled as a lumped RLC boundary condition — a rectangular sheet with a specified sheet inductance \(L_J\). This is the correct linearized representation because, for small oscillations around the superconducting phase minimum, the JJ potential energy \(-E_J \cos(\varphi) \approx E_J \varphi^2 / 2 - E_J \varphi^4 / 24 + \ldots\) is dominated by the harmonic (quadratic) term, which corresponds exactly to an inductor with inductance: $\(L_J = \frac{\Phi_0^2}{(2\pi)^2 E_J} = \frac{\hbar^2}{4 e^2 E_J}\)$
The anharmonicity — the deviation from a perfect harmonic oscillator that makes the qubit addressable — comes from the quartic and higher-order terms in the cosine expansion. The EPR method captures these nonlinear corrections in Stage 3 (quantum Hamiltonian analysis), after extracting the linear mode structure in Stages 1 and 2.
Each junction in the model is described by four pieces of information:
Lj_variable: the name of the HFSS design variable that specifies \(L_J\) on the lumped boundary condition. pyEPR reads this variable to compute \(E_J = \Phi_0^2 / (4\pi^2 L_J)\).rect: the name of the HFSS rectangle on which the lumped boundary condition is applied. The area of this rectangle determines the current density normalization.line: the name of the HFSS polyline that spans the length of the rectangle, parallel to the current flow direction. The zero-point voltage is computed by integrating the electric field \(\mathbf{E}\) along this line.length: the physical length of the junction (used for normalization in the EPR line-voltage calculation).
An optional Cj_variable can be specified if a lumped capacitance is placed in parallel with the junction in the HFSS model.
To analyze the quantum part of the model, let us now register the junctions present in the design. We call the single junction j1:
pinfo.junctions['j1'] = {'Lj_variable' : 'Lj_1',
'rect' : 'rect_jj1',
'line' : 'line_jj1',
'length' : epr.parse_units('100um')}
# Check that valid names of variables and objects have been supplied.
# An error is raised with a message if something is wrong.
pinfo.validate_junction_info()
Microwave Analysis:
Run analysis on an eigenmode solution #
Assuming we have field solutions saved for the HFSS eigenmode design, we can now create the analysis object from the class DistributedAnalysis.
Run analysis (tutorial action): Run the analysis in the design model. This can be done manually in the HFSS GUI, or programmatically using pyEPR.
# Calling this function will run an analysis of the design in the background.
# The console will wait for HFSS to complete.
pinfo.setup.analyze()
Run optimetrics
If an optimetrics (parametric sweep) setup is defined in HFSS, the following call will execute it. Each parametric point generates a separate set of field solutions that pyEPR can subsequently analyze.
The following will get the frequencies and Qs from HFSS.
# Calling this function will run an analysis on all the optimetrics the design in the background.
# The console will wait for HFSS to complete.
# This requires optimetrics license by HFSS.
pinfo.design.optimetrics.solve_setup(pinfo.design.optimetrics.get_setup_names()[0])
EPR Analysis object
Let us introduce the core analysis object of pyEPR, which will be a class instance of epr.DistributedAnalysis. This is the classical-to-quantum bridge: it connects to HFSS, retrieves the electromagnetic field solutions, and computes the energy-participation ratios for each mode and junction.
Let us recover the eigenmode frequencies of the first variation we ran, which is indexed by '0'.
# This is the core object for interacting with HFSS
# and running analysis within HFSS.
eprh = epr.DistributedAnalysis(pinfo) # epr hfss analysis
eprh.get_ansys_variables()
eprh.get_ansys_frequencies_all()
eprh.hfss_report_full_convergence();
We can see that the convergence is reasonable for this tutorial example. In production designs, one typically targets \(\Delta f / f < 0.1\%\) between adaptive mesh passes and a monotonically decreasing convergence curve. Mesh refinement operations (e.g., seeding finer elements around the junction rectangle or the qubit capacitor pads) can improve convergence.
For this tutorial, let us press forward.
You can use these autogenerated reports to keep track of your simulations. This is still all at the classical level. Take a look at the other available methods:
hfss_report_full_convergence
hfss_report_f_convergence
get_convergence
get_mesh_statistics
get_freqs_bare_pd
…
For example, you can get tables of the frequencies and Q factors using:
Fs, Qs = eprh.get_freqs_bare_pd(variation='0')
import pandas as pd
for variation in eprh.variations[:2]: # just for the first 2
Fs, Qs = eprh.get_freqs_bare_pd(variation=variation, frame=False)
display(pd.DataFrame({'Freq. (GHz)':Fs, 'Quality Factor':Qs}))
The Q is infinite, since we have not included dissipation in this example — there are no lossy boundaries or bulk dielectrics defined in the HFSS model. The row index is the mode number:
Mode 0 (~5 GHz): the transmon qubit mode
Mode 1 (~9 GHz): the cavity mode
Dissipation can be modeled by specifying the corresponding HFSS objects for bulk dielectrics, surface losses, resistive surfaces, or seams when initializing ProjectInfo. For example:
pinfo = epr.ProjectInfo(project_path = path_to_project,
project_name = 'pyEPR_tutorial1',
design_name = '1. single_transmon',
dielectrics_bulk = ['si_substrate'],
dielectric_surfaces = ['interface'],
resistive_surfaces = None,
seams = None)
or equivalently, after the fact:
pinfo.dissipative['dielectrics_bulk'] = ['si_substrate']
pinfo.dissipative['dielectric_surfaces'] = ['interface']
See Tutorial 2 for a detailed treatment of dielectric loss calculation.
Full EPR Analysis#
What does do_EPR_analysis compute?#
This is the heart of the pyEPR workflow — the classical-to-quantum bridge. For each eigenmode and each parametric variation (e.g., each value of \(L_J\) in a sweep), do_EPR_analysis performs the following steps:
Field retrieval: Queries HFSS for the electric and magnetic field distributions \(\mathbf{E}_m(\mathbf{r})\), \(\mathbf{H}_m(\mathbf{r})\) of mode \(m\).
Energy balance check: Computes the total electric energy \(\mathcal{E}_m^{(E)}\) and magnetic energy \(\mathcal{E}_m^{(H)}\) by volume integration over the entire simulation domain:
$\(\mathcal{E}_m^{(E)} = \frac{1}{4}\mathrm{Re}\int_V dv\, \mathbf{E}_m^* \cdot \overleftrightarrow{\epsilon} \cdot \mathbf{E}_m\)\(
For a well-converged eigenmode, \)\mathcal{E}_m^{(E)} \approx \mathcal{E}_m^{(H)}$ (energy equipartition). A large imbalance signals poor mesh convergence.
Junction EPR calculation: For each junction \(j\), computes the participation ratio \(p_{mj}\) by integrating the electric field along the polyline \(\ell_j\):
$\(p_{mj} = \frac{\omega_m L_j \left| \int_{\ell_j} \mathbf{E}_m \cdot d\mathbf{\ell} \right|^2}{4 \omega_m^2 \mathcal{E}_m^{(E)}}\)$
See Eq. (4) and Appendix B of arXiv:2010.00620.
Sign determination: Computes the sign \(s_{mj} = \pm 1\) of the zero-point flux fluctuation, which determines whether modes are hybridized constructively or destructively at the junction.
Data storage: Saves all results (frequencies, participation matrices, sign matrices, and energy values) to an HDF5 file for subsequent quantum analysis.
The output printed during the run shows, for each mode, the energy balance \((\mathcal{E}_E - \mathcal{E}_H)/\mathcal{E}_E\) and the computed EPR \(p_{mj}\) for each junction.
# Do EPR analysis for all modes and all variations
# If you want to only analze a certain variation, you can change the
# default keyword argeumetns variations=None, modes=None
eprh.do_EPR_analysis();
Getting more info#
In general, to see more verbose output when running commands, set the logger level to DEBUG.
import logging
epr.logger.setLevel(logging.DEBUG)
Quantum Hamiltonian Analysis:
Load microwave analysis results and perform Hamiltonian analysis on the solutions #
The quantum Hamiltonian and its diagonalization#
Having computed the EPR participation matrix \(p_{mj}\) and zero-point flux fluctuations \(\varphi_{zpf}^{(mj)}\) for each mode–junction pair, we can now construct and diagonalize the full quantum Hamiltonian. This is done by the QuantumAnalysis class.
The system Hamiltonian (in the basis of the linearized HFSS eigenmodes) is (Eq. (1) of arXiv:2010.00620):
$\(\hat{H} = \sum_m \hbar\omega_m \hat{a}_m^\dagger \hat{a}_m - \sum_j E_J^{(j)} \left[\cos\!\left(\hat{\varphi}_j\right) - 1 + \frac{\hat{\varphi}_j^2}{2}\right]\)\(
where the phase operator for junction \)j\( is:
\)\(\hat{\varphi}_j = \sum_m s_{mj}\, \varphi_{zpf}^{(mj)} \left(\hat{a}_m + \hat{a}_m^\dagger\right)\)$
The cosine is expanded as a Taylor series (truncated at order cos_trunc) and the resulting polynomial in bosonic operators is diagonalized numerically in a truncated Fock space of dimension fock_trunc per mode. Two levels of analysis are provided:
First-order perturbation theory (PT): Analytic result valid when \(\varphi_{zpf} \ll 1\) (transmon regime). Fast, but less accurate for strongly anharmonic or hybridized modes.
Numerical diagonalization (ND): Full numerical diagonalization in the Fock space. More accurate, but computationally more expensive.
Key parameters:
cos_trunc: Number of terms in the Taylor expansion of \(\cos(\hat{\varphi})\) kept in the Hamiltonian. Typical value: 8. Higher values improve accuracy at the cost of computation time.
fock_trunc: Dimension of the truncated Fock space per mode. Typical value: 15 (i.e., states \(|0\rangle, |1\rangle, \ldots, |14\rangle\)). Must be large enough that truncation errors are negligible.
%matplotlib inline
# plot the figures in the notebook.
# This is the core object for working with post processing data from HFSS
# This iwll load the saved solutions from above
epra = epr.QuantumAnalysis(eprh.data_filename)
# Analyze
epra.analyze_all_variations(cos_trunc = 8, fock_trunc = 15)
epra.plot_hamiltonian_results();
Interpreting the Hamiltonian results#
The plot_hamiltonian_results() output and the printed tables give the following quantum quantities:
Chi matrix (\(\chi\)):
The \(\chi\) matrix encodes all second-order quantum corrections from the junction nonlinearity. In the standard dispersive approximation it takes the form \(\hat{H} \supset \sum_{m,n} \chi_{mn} \hat{a}_m^\dagger \hat{a}_m \hat{a}_n^\dagger \hat{a}_n\).
Diagonal elements \(\chi_{mm}\) (labeled anharmonicity \(\alpha_m\)): The self-Kerr coefficient, which sets how much the frequency of mode \(m\) shifts for each additional photon in that mode. For the transmon, \(\alpha_q \approx -E_C/h\) (where \(E_C = e^2 / (2C_\Sigma)\) is the charging energy). A transmon with \(\alpha \approx 200\,\text{MHz}\) can be selectively addressed with \(\sim 10\,\text{ns}\) pulses.
Off-diagonal elements \(\chi_{mn}\) (\(m \ne n\)): The cross-Kerr (dispersive coupling) between modes \(m\) and \(n\). In a qubit-cavity system, \(\chi_{qc}\) is the dispersive shift used for qubit readout: the cavity resonance shifts by \(\chi_{qc}\) depending on whether the qubit is in \(|0\rangle\) or \(|1\rangle\).
Dressed frequencies vs. bare frequencies:
The bare frequencies \(\omega_m\) are the HFSS eigenfrequencies of the linearized circuit. The dressed frequencies \(\tilde{\omega}_m\) include the Lamb shift — a downward frequency renormalization due to the quantum zero-point fluctuations coupling each mode to the junction nonlinearity. This shift is \(\tilde{\omega}_m = \omega_m + \chi_{mm}/2 + \ldots\) to leading order in perturbation theory.
The two columns labeled O1 PT and ND give the first-order perturbation theory and numerical diagonalization results respectively. For a typical transmon with small \(\varphi_{zpf}\), these agree to within a few percent.
See Section III and Table I of arXiv:2010.00620 for a systematic comparison of these methods against experiment.
Swept variable analysis driven through pyEPR:
Sweep variable in HFSS, analyze microwaves, then analyze quantum#
Let us now run a more complete example, where we sweep the Josephson inductance \(L_J\). This demonstrates the full power of the pyEPR workflow: by sweeping \(L_J\) (which is equivalent to varying the Josephson energy \(E_J = \Phi_0^2 / (4\pi^2 L_J)\)), we can trace out how the qubit frequency, anharmonicity, and dispersive coupling all depend on the junction parameter.
This is exactly how one would characterize the design space during the development of a new qubit device.
Classical microwave, distributed analysis
Sweep variable in HFSS, run the eigenmode solver for each value.
import numpy as np
### Configure setup
### Set the maximum number of passes on the analysis and a few other
pinfo.setup.passes = 15
pinfo.setup.min_freq = '3.2GHz'
pinfo.setup.modes = 2
pinfo.setup.delta_f = 0.02
pinfo.setup.basis_order = epr.ansys.BASIS_ORDER['Mixed Order']
### What to sweep
# Define sweep variables to go over
swp_var = 'Lj_1' # name
# how to convert sweep param to string variable in nanoHeries.
to_swp_val = lambda x: f'{x}nH'
### Run Sweep
for swp_param in np.linspace(6, 12, 6):
swp_val = to_swp_val(swp_param)
epr.logger.info(f'Setting sweep variable {swp_var}={swp_val}')
pinfo.design.set_variable(swp_var, swp_val)
pinfo.setup.analyze()
EPR microwave analysis
Classical-to-quantum bridge: compute participation ratios for all sweep points.
eprh = epr.DistributedAnalysis(pinfo) # epr hfss analysis
eprh.hfss_report_full_convergence()
# Do EPR analysis for all modes and all variations
eprh.do_EPR_analysis();
Quantum analysis
Perform the Hamiltonian diagonalization and plot the results as a function of the swept parameter \(L_J\).
# This is the core object for working with post processing data from HFSS
# This will load the saved solutions from above
epra = epr.QuantumAnalysis(eprh.data_filename)
# Analyze
epra.analyze_all_variations(cos_trunc = 8, fock_trunc = 15);
epra.plot_hamiltonian_results(swp_variable='Lj_1');
There is some jitter in the results since some of the runs were performed with only a few adaptive passes. In a production analysis one would ensure all variations have converged to the same \(\Delta f\) threshold before proceeding.
The sweep makes several physically important trends visible:
As \(L_J\) increases, \(E_J\) decreases, so the qubit frequency \(\omega_q \propto \sqrt{E_J}\) decreases.
The anharmonicity \(\alpha \approx -E_C\) is relatively insensitive to \(L_J\) in the transmon regime (since \(E_C\) is set by the capacitor geometry, not the junction).
The dispersive coupling \(\chi_{qc}\) decreases as \(L_J\) increases (equivalently, as \(E_J / E_C\) increases and the qubit moves deeper into the transmon regime).
See Figure 3 of arXiv:2010.00620 for an example of these trends in a real device.
epra.quick_plot_mode(0,0,1,numeric=True,swp_variable='Lj_1')
epra.report_results(swp_variable='Lj_1', numeric=True)
Addendum:
Snippets for ease #
An easy way to define many junctions is the following, using Python 3 f-strings for ease. This pattern is useful for multi-junction circuits such as flux qubits, SQUID arrays, or fluxonium qubits, where each junction must be registered individually but shares a common naming convention:
# Specify multiple junctions in HFSS model
n_junctions = 5
for i in range(1, 1+n_junctions):
pinfo.junctions[f'j{i}'] = {'Lj_variable' : f'Lj{i}',
'rect' : f'JJrect{i}',
'line' : f'JJline{i}',
'length' : epr.parse_units('100um')}
NEXT: Please see Tutorial 2 to continue. It covers the calculation of dielectric energy-participation ratios, which are needed to estimate qubit lifetimes \(T_1\) from dielectric loss tangents.
Run analysis on an eigenmode solution
Assuming we have field solutions saved for the HFSS eigenmode design, we can now create the analysis object from the class DistributedAnalysis.
Run analysis (tutorial action): Run the analysis in the design model. This can be done manually in the HFSS GUI, or programmatically using pyEPR.
# Calling this function will run an analysis of the design in the background.
# The console will wait for HFSS to complete.
pinfo.setup.analyze()
Run optimetrics
If an optimetrics (parametric sweep) setup is defined in HFSS, the following call will execute it. Each parametric point generates a separate set of field solutions that pyEPR can subsequently analyze.
The following will get the frequencies and Qs from HFSS.
# Calling this function will run an analysis on all the optimetrics the design in the background.
# The console will wait for HFSS to complete.
# This requires optimetrics license by HFSS.
pinfo.design.optimetrics.solve_setup(pinfo.design.optimetrics.get_setup_names()[0])
EPR Analysis object
Let us introduce the core analysis object of pyEPR, which will be a class instance of epr.DistributedAnalysis. This is the classical-to-quantum bridge: it connects to HFSS, retrieves the electromagnetic field solutions, and computes the energy-participation ratios for each mode and junction.
Let us recover the eigenmode frequencies of the first variation we ran, which is indexed by '0'.
# This is the core object for interacting with HFSS
# and running analysis within HFSS.
eprh = epr.DistributedAnalysis(pinfo) # epr hfss analysis
eprh.get_ansys_variables()
eprh.get_ansys_frequencies_all()
eprh.hfss_report_full_convergence();
We can see that the convergence is reasonable for this tutorial example. In production designs, one typically targets \(\Delta f / f < 0.1\%\) between adaptive mesh passes and a monotonically decreasing convergence curve. Mesh refinement operations (e.g., seeding finer elements around the junction rectangle or the qubit capacitor pads) can improve convergence.
For this tutorial, let us press forward.
You can use these autogenerated reports to keep track of your simulations. This is still all at the classical level. Take a look at the other available methods:
hfss_report_full_convergencehfss_report_f_convergenceget_convergenceget_mesh_statisticsget_freqs_bare_pd…
For example, you can get tables of the frequencies and Q factors using:
Fs, Qs = eprh.get_freqs_bare_pd(variation='0')
import pandas as pd
for variation in eprh.variations[:2]: # just for the first 2
Fs, Qs = eprh.get_freqs_bare_pd(variation=variation, frame=False)
display(pd.DataFrame({'Freq. (GHz)':Fs, 'Quality Factor':Qs}))
The Q is infinite, since we have not included dissipation in this example — there are no lossy boundaries or bulk dielectrics defined in the HFSS model. The row index is the mode number:
Mode 0 (~5 GHz): the transmon qubit mode
Mode 1 (~9 GHz): the cavity mode
Dissipation can be modeled by specifying the corresponding HFSS objects for bulk dielectrics, surface losses, resistive surfaces, or seams when initializing ProjectInfo. For example:
pinfo = epr.ProjectInfo(project_path = path_to_project,
project_name = 'pyEPR_tutorial1',
design_name = '1. single_transmon',
dielectrics_bulk = ['si_substrate'],
dielectric_surfaces = ['interface'],
resistive_surfaces = None,
seams = None)
or equivalently, after the fact:
pinfo.dissipative['dielectrics_bulk'] = ['si_substrate']
pinfo.dissipative['dielectric_surfaces'] = ['interface']
See Tutorial 2 for a detailed treatment of dielectric loss calculation.
Full EPR Analysis#
What does do_EPR_analysis compute?#
This is the heart of the pyEPR workflow — the classical-to-quantum bridge. For each eigenmode and each parametric variation (e.g., each value of \(L_J\) in a sweep), do_EPR_analysis performs the following steps:
Field retrieval: Queries HFSS for the electric and magnetic field distributions \(\mathbf{E}_m(\mathbf{r})\), \(\mathbf{H}_m(\mathbf{r})\) of mode \(m\).
Energy balance check: Computes the total electric energy \(\mathcal{E}_m^{(E)}\) and magnetic energy \(\mathcal{E}_m^{(H)}\) by volume integration over the entire simulation domain: $\(\mathcal{E}_m^{(E)} = \frac{1}{4}\mathrm{Re}\int_V dv\, \mathbf{E}_m^* \cdot \overleftrightarrow{\epsilon} \cdot \mathbf{E}_m\)\( For a well-converged eigenmode, \)\mathcal{E}_m^{(E)} \approx \mathcal{E}_m^{(H)}$ (energy equipartition). A large imbalance signals poor mesh convergence.
Junction EPR calculation: For each junction \(j\), computes the participation ratio \(p_{mj}\) by integrating the electric field along the polyline \(\ell_j\): $\(p_{mj} = \frac{\omega_m L_j \left| \int_{\ell_j} \mathbf{E}_m \cdot d\mathbf{\ell} \right|^2}{4 \omega_m^2 \mathcal{E}_m^{(E)}}\)$ See Eq. (4) and Appendix B of arXiv:2010.00620.
Sign determination: Computes the sign \(s_{mj} = \pm 1\) of the zero-point flux fluctuation, which determines whether modes are hybridized constructively or destructively at the junction.
Data storage: Saves all results (frequencies, participation matrices, sign matrices, and energy values) to an HDF5 file for subsequent quantum analysis.
The output printed during the run shows, for each mode, the energy balance \((\mathcal{E}_E - \mathcal{E}_H)/\mathcal{E}_E\) and the computed EPR \(p_{mj}\) for each junction.
# Do EPR analysis for all modes and all variations
# If you want to only analze a certain variation, you can change the
# default keyword argeumetns variations=None, modes=None
eprh.do_EPR_analysis();
Getting more info#
In general, to see more verbose output when running commands, set the logger level to DEBUG.
import logging
epr.logger.setLevel(logging.DEBUG)
Quantum Hamiltonian Analysis:
Load microwave analysis results and perform Hamiltonian analysis on the solutions #
The quantum Hamiltonian and its diagonalization#
Having computed the EPR participation matrix \(p_{mj}\) and zero-point flux fluctuations \(\varphi_{zpf}^{(mj)}\) for each mode–junction pair, we can now construct and diagonalize the full quantum Hamiltonian. This is done by the QuantumAnalysis class.
The system Hamiltonian (in the basis of the linearized HFSS eigenmodes) is (Eq. (1) of arXiv:2010.00620):
$\(\hat{H} = \sum_m \hbar\omega_m \hat{a}_m^\dagger \hat{a}_m - \sum_j E_J^{(j)} \left[\cos\!\left(\hat{\varphi}_j\right) - 1 + \frac{\hat{\varphi}_j^2}{2}\right]\)\(
where the phase operator for junction \)j\( is:
\)\(\hat{\varphi}_j = \sum_m s_{mj}\, \varphi_{zpf}^{(mj)} \left(\hat{a}_m + \hat{a}_m^\dagger\right)\)$
The cosine is expanded as a Taylor series (truncated at order cos_trunc) and the resulting polynomial in bosonic operators is diagonalized numerically in a truncated Fock space of dimension fock_trunc per mode. Two levels of analysis are provided:
First-order perturbation theory (PT): Analytic result valid when \(\varphi_{zpf} \ll 1\) (transmon regime). Fast, but less accurate for strongly anharmonic or hybridized modes.
Numerical diagonalization (ND): Full numerical diagonalization in the Fock space. More accurate, but computationally more expensive.
Key parameters:
cos_trunc: Number of terms in the Taylor expansion of \(\cos(\hat{\varphi})\) kept in the Hamiltonian. Typical value: 8. Higher values improve accuracy at the cost of computation time.
fock_trunc: Dimension of the truncated Fock space per mode. Typical value: 15 (i.e., states \(|0\rangle, |1\rangle, \ldots, |14\rangle\)). Must be large enough that truncation errors are negligible.
%matplotlib inline
# plot the figures in the notebook.
# This is the core object for working with post processing data from HFSS
# This iwll load the saved solutions from above
epra = epr.QuantumAnalysis(eprh.data_filename)
# Analyze
epra.analyze_all_variations(cos_trunc = 8, fock_trunc = 15)
epra.plot_hamiltonian_results();
Interpreting the Hamiltonian results#
The plot_hamiltonian_results() output and the printed tables give the following quantum quantities:
Chi matrix (\(\chi\)):
The \(\chi\) matrix encodes all second-order quantum corrections from the junction nonlinearity. In the standard dispersive approximation it takes the form \(\hat{H} \supset \sum_{m,n} \chi_{mn} \hat{a}_m^\dagger \hat{a}_m \hat{a}_n^\dagger \hat{a}_n\).
Diagonal elements \(\chi_{mm}\) (labeled anharmonicity \(\alpha_m\)): The self-Kerr coefficient, which sets how much the frequency of mode \(m\) shifts for each additional photon in that mode. For the transmon, \(\alpha_q \approx -E_C/h\) (where \(E_C = e^2 / (2C_\Sigma)\) is the charging energy). A transmon with \(\alpha \approx 200\,\text{MHz}\) can be selectively addressed with \(\sim 10\,\text{ns}\) pulses.
Off-diagonal elements \(\chi_{mn}\) (\(m \ne n\)): The cross-Kerr (dispersive coupling) between modes \(m\) and \(n\). In a qubit-cavity system, \(\chi_{qc}\) is the dispersive shift used for qubit readout: the cavity resonance shifts by \(\chi_{qc}\) depending on whether the qubit is in \(|0\rangle\) or \(|1\rangle\).
Dressed frequencies vs. bare frequencies:
The bare frequencies \(\omega_m\) are the HFSS eigenfrequencies of the linearized circuit. The dressed frequencies \(\tilde{\omega}_m\) include the Lamb shift — a downward frequency renormalization due to the quantum zero-point fluctuations coupling each mode to the junction nonlinearity. This shift is \(\tilde{\omega}_m = \omega_m + \chi_{mm}/2 + \ldots\) to leading order in perturbation theory.
The two columns labeled O1 PT and ND give the first-order perturbation theory and numerical diagonalization results respectively. For a typical transmon with small \(\varphi_{zpf}\), these agree to within a few percent.
See Section III and Table I of arXiv:2010.00620 for a systematic comparison of these methods against experiment.
Swept variable analysis driven through pyEPR:
Sweep variable in HFSS, analyze microwaves, then analyze quantum#
Let us now run a more complete example, where we sweep the Josephson inductance \(L_J\). This demonstrates the full power of the pyEPR workflow: by sweeping \(L_J\) (which is equivalent to varying the Josephson energy \(E_J = \Phi_0^2 / (4\pi^2 L_J)\)), we can trace out how the qubit frequency, anharmonicity, and dispersive coupling all depend on the junction parameter.
This is exactly how one would characterize the design space during the development of a new qubit device.
Classical microwave, distributed analysis
Sweep variable in HFSS, run the eigenmode solver for each value.
import numpy as np
### Configure setup
### Set the maximum number of passes on the analysis and a few other
pinfo.setup.passes = 15
pinfo.setup.min_freq = '3.2GHz'
pinfo.setup.modes = 2
pinfo.setup.delta_f = 0.02
pinfo.setup.basis_order = epr.ansys.BASIS_ORDER['Mixed Order']
### What to sweep
# Define sweep variables to go over
swp_var = 'Lj_1' # name
# how to convert sweep param to string variable in nanoHeries.
to_swp_val = lambda x: f'{x}nH'
### Run Sweep
for swp_param in np.linspace(6, 12, 6):
swp_val = to_swp_val(swp_param)
epr.logger.info(f'Setting sweep variable {swp_var}={swp_val}')
pinfo.design.set_variable(swp_var, swp_val)
pinfo.setup.analyze()
EPR microwave analysis
Classical-to-quantum bridge: compute participation ratios for all sweep points.
eprh = epr.DistributedAnalysis(pinfo) # epr hfss analysis
eprh.hfss_report_full_convergence()
# Do EPR analysis for all modes and all variations
eprh.do_EPR_analysis();
Quantum analysis
Perform the Hamiltonian diagonalization and plot the results as a function of the swept parameter \(L_J\).
# This is the core object for working with post processing data from HFSS
# This will load the saved solutions from above
epra = epr.QuantumAnalysis(eprh.data_filename)
# Analyze
epra.analyze_all_variations(cos_trunc = 8, fock_trunc = 15);
epra.plot_hamiltonian_results(swp_variable='Lj_1');
There is some jitter in the results since some of the runs were performed with only a few adaptive passes. In a production analysis one would ensure all variations have converged to the same \(\Delta f\) threshold before proceeding.
The sweep makes several physically important trends visible:
As \(L_J\) increases, \(E_J\) decreases, so the qubit frequency \(\omega_q \propto \sqrt{E_J}\) decreases.
The anharmonicity \(\alpha \approx -E_C\) is relatively insensitive to \(L_J\) in the transmon regime (since \(E_C\) is set by the capacitor geometry, not the junction).
The dispersive coupling \(\chi_{qc}\) decreases as \(L_J\) increases (equivalently, as \(E_J / E_C\) increases and the qubit moves deeper into the transmon regime).
See Figure 3 of arXiv:2010.00620 for an example of these trends in a real device.
epra.quick_plot_mode(0,0,1,numeric=True,swp_variable='Lj_1')
epra.report_results(swp_variable='Lj_1', numeric=True)
Addendum:
Snippets for ease #
An easy way to define many junctions is the following, using Python 3 f-strings for ease. This pattern is useful for multi-junction circuits such as flux qubits, SQUID arrays, or fluxonium qubits, where each junction must be registered individually but shares a common naming convention:
# Specify multiple junctions in HFSS model
n_junctions = 5
for i in range(1, 1+n_junctions):
pinfo.junctions[f'j{i}'] = {'Lj_variable' : f'Lj{i}',
'rect' : f'JJrect{i}',
'line' : f'JJline{i}',
'length' : epr.parse_units('100um')}
NEXT: Please see Tutorial 2 to continue. It covers the calculation of dielectric energy-participation ratios, which are needed to estimate qubit lifetimes \(T_1\) from dielectric loss tangents.
Load microwave analysis results and perform Hamiltonian analysis on the solutions
The quantum Hamiltonian and its diagonalization#
Having computed the EPR participation matrix \(p_{mj}\) and zero-point flux fluctuations \(\varphi_{zpf}^{(mj)}\) for each mode–junction pair, we can now construct and diagonalize the full quantum Hamiltonian. This is done by the QuantumAnalysis class.
The system Hamiltonian (in the basis of the linearized HFSS eigenmodes) is (Eq. (1) of arXiv:2010.00620): $\(\hat{H} = \sum_m \hbar\omega_m \hat{a}_m^\dagger \hat{a}_m - \sum_j E_J^{(j)} \left[\cos\!\left(\hat{\varphi}_j\right) - 1 + \frac{\hat{\varphi}_j^2}{2}\right]\)\( where the phase operator for junction \)j\( is: \)\(\hat{\varphi}_j = \sum_m s_{mj}\, \varphi_{zpf}^{(mj)} \left(\hat{a}_m + \hat{a}_m^\dagger\right)\)$
The cosine is expanded as a Taylor series (truncated at order cos_trunc) and the resulting polynomial in bosonic operators is diagonalized numerically in a truncated Fock space of dimension fock_trunc per mode. Two levels of analysis are provided:
First-order perturbation theory (PT): Analytic result valid when \(\varphi_{zpf} \ll 1\) (transmon regime). Fast, but less accurate for strongly anharmonic or hybridized modes.
Numerical diagonalization (ND): Full numerical diagonalization in the Fock space. More accurate, but computationally more expensive.
Key parameters:
cos_trunc: Number of terms in the Taylor expansion of \(\cos(\hat{\varphi})\) kept in the Hamiltonian. Typical value: 8. Higher values improve accuracy at the cost of computation time.fock_trunc: Dimension of the truncated Fock space per mode. Typical value: 15 (i.e., states \(|0\rangle, |1\rangle, \ldots, |14\rangle\)). Must be large enough that truncation errors are negligible.
%matplotlib inline
# plot the figures in the notebook.
# This is the core object for working with post processing data from HFSS
# This iwll load the saved solutions from above
epra = epr.QuantumAnalysis(eprh.data_filename)
# Analyze
epra.analyze_all_variations(cos_trunc = 8, fock_trunc = 15)
epra.plot_hamiltonian_results();
Interpreting the Hamiltonian results#
The plot_hamiltonian_results() output and the printed tables give the following quantum quantities:
Chi matrix (\(\chi\)): The \(\chi\) matrix encodes all second-order quantum corrections from the junction nonlinearity. In the standard dispersive approximation it takes the form \(\hat{H} \supset \sum_{m,n} \chi_{mn} \hat{a}_m^\dagger \hat{a}_m \hat{a}_n^\dagger \hat{a}_n\).
Diagonal elements \(\chi_{mm}\) (labeled
anharmonicity\(\alpha_m\)): The self-Kerr coefficient, which sets how much the frequency of mode \(m\) shifts for each additional photon in that mode. For the transmon, \(\alpha_q \approx -E_C/h\) (where \(E_C = e^2 / (2C_\Sigma)\) is the charging energy). A transmon with \(\alpha \approx 200\,\text{MHz}\) can be selectively addressed with \(\sim 10\,\text{ns}\) pulses.Off-diagonal elements \(\chi_{mn}\) (\(m \ne n\)): The cross-Kerr (dispersive coupling) between modes \(m\) and \(n\). In a qubit-cavity system, \(\chi_{qc}\) is the dispersive shift used for qubit readout: the cavity resonance shifts by \(\chi_{qc}\) depending on whether the qubit is in \(|0\rangle\) or \(|1\rangle\).
Dressed frequencies vs. bare frequencies: The bare frequencies \(\omega_m\) are the HFSS eigenfrequencies of the linearized circuit. The dressed frequencies \(\tilde{\omega}_m\) include the Lamb shift — a downward frequency renormalization due to the quantum zero-point fluctuations coupling each mode to the junction nonlinearity. This shift is \(\tilde{\omega}_m = \omega_m + \chi_{mm}/2 + \ldots\) to leading order in perturbation theory.
The two columns labeled O1 PT and ND give the first-order perturbation theory and numerical diagonalization results respectively. For a typical transmon with small \(\varphi_{zpf}\), these agree to within a few percent.
See Section III and Table I of arXiv:2010.00620 for a systematic comparison of these methods against experiment.
Swept variable analysis driven through pyEPR:
Sweep variable in HFSS, analyze microwaves, then analyze quantum#
Let us now run a more complete example, where we sweep the Josephson inductance \(L_J\). This demonstrates the full power of the pyEPR workflow: by sweeping \(L_J\) (which is equivalent to varying the Josephson energy \(E_J = \Phi_0^2 / (4\pi^2 L_J)\)), we can trace out how the qubit frequency, anharmonicity, and dispersive coupling all depend on the junction parameter.
This is exactly how one would characterize the design space during the development of a new qubit device.
Classical microwave, distributed analysis
Sweep variable in HFSS, run the eigenmode solver for each value.
import numpy as np
### Configure setup
### Set the maximum number of passes on the analysis and a few other
pinfo.setup.passes = 15
pinfo.setup.min_freq = '3.2GHz'
pinfo.setup.modes = 2
pinfo.setup.delta_f = 0.02
pinfo.setup.basis_order = epr.ansys.BASIS_ORDER['Mixed Order']
### What to sweep
# Define sweep variables to go over
swp_var = 'Lj_1' # name
# how to convert sweep param to string variable in nanoHeries.
to_swp_val = lambda x: f'{x}nH'
### Run Sweep
for swp_param in np.linspace(6, 12, 6):
swp_val = to_swp_val(swp_param)
epr.logger.info(f'Setting sweep variable {swp_var}={swp_val}')
pinfo.design.set_variable(swp_var, swp_val)
pinfo.setup.analyze()
EPR microwave analysis
Classical-to-quantum bridge: compute participation ratios for all sweep points.
eprh = epr.DistributedAnalysis(pinfo) # epr hfss analysis
eprh.hfss_report_full_convergence()
# Do EPR analysis for all modes and all variations
eprh.do_EPR_analysis();
Quantum analysis
Perform the Hamiltonian diagonalization and plot the results as a function of the swept parameter \(L_J\).
# This is the core object for working with post processing data from HFSS
# This will load the saved solutions from above
epra = epr.QuantumAnalysis(eprh.data_filename)
# Analyze
epra.analyze_all_variations(cos_trunc = 8, fock_trunc = 15);
epra.plot_hamiltonian_results(swp_variable='Lj_1');
There is some jitter in the results since some of the runs were performed with only a few adaptive passes. In a production analysis one would ensure all variations have converged to the same \(\Delta f\) threshold before proceeding.
The sweep makes several physically important trends visible:
As \(L_J\) increases, \(E_J\) decreases, so the qubit frequency \(\omega_q \propto \sqrt{E_J}\) decreases.
The anharmonicity \(\alpha \approx -E_C\) is relatively insensitive to \(L_J\) in the transmon regime (since \(E_C\) is set by the capacitor geometry, not the junction).
The dispersive coupling \(\chi_{qc}\) decreases as \(L_J\) increases (equivalently, as \(E_J / E_C\) increases and the qubit moves deeper into the transmon regime).
See Figure 3 of arXiv:2010.00620 for an example of these trends in a real device.
epra.quick_plot_mode(0,0,1,numeric=True,swp_variable='Lj_1')
epra.report_results(swp_variable='Lj_1', numeric=True)
Addendum:
Snippets for ease #
An easy way to define many junctions is the following, using Python 3 f-strings for ease. This pattern is useful for multi-junction circuits such as flux qubits, SQUID arrays, or fluxonium qubits, where each junction must be registered individually but shares a common naming convention:
# Specify multiple junctions in HFSS model
n_junctions = 5
for i in range(1, 1+n_junctions):
pinfo.junctions[f'j{i}'] = {'Lj_variable' : f'Lj{i}',
'rect' : f'JJrect{i}',
'line' : f'JJline{i}',
'length' : epr.parse_units('100um')}
NEXT: Please see Tutorial 2 to continue. It covers the calculation of dielectric energy-participation ratios, which are needed to estimate qubit lifetimes \(T_1\) from dielectric loss tangents.
Sweep variable in HFSS, analyze microwaves, then analyze quantum
Let us now run a more complete example, where we sweep the Josephson inductance \(L_J\). This demonstrates the full power of the pyEPR workflow: by sweeping \(L_J\) (which is equivalent to varying the Josephson energy \(E_J = \Phi_0^2 / (4\pi^2 L_J)\)), we can trace out how the qubit frequency, anharmonicity, and dispersive coupling all depend on the junction parameter.
This is exactly how one would characterize the design space during the development of a new qubit device.
Classical microwave, distributed analysis
Sweep variable in HFSS, run the eigenmode solver for each value.
import numpy as np
### Configure setup
### Set the maximum number of passes on the analysis and a few other
pinfo.setup.passes = 15
pinfo.setup.min_freq = '3.2GHz'
pinfo.setup.modes = 2
pinfo.setup.delta_f = 0.02
pinfo.setup.basis_order = epr.ansys.BASIS_ORDER['Mixed Order']
### What to sweep
# Define sweep variables to go over
swp_var = 'Lj_1' # name
# how to convert sweep param to string variable in nanoHeries.
to_swp_val = lambda x: f'{x}nH'
### Run Sweep
for swp_param in np.linspace(6, 12, 6):
swp_val = to_swp_val(swp_param)
epr.logger.info(f'Setting sweep variable {swp_var}={swp_val}')
pinfo.design.set_variable(swp_var, swp_val)
pinfo.setup.analyze()
EPR microwave analysis
Classical-to-quantum bridge: compute participation ratios for all sweep points.
eprh = epr.DistributedAnalysis(pinfo) # epr hfss analysis
eprh.hfss_report_full_convergence()
# Do EPR analysis for all modes and all variations
eprh.do_EPR_analysis();
Quantum analysis
Perform the Hamiltonian diagonalization and plot the results as a function of the swept parameter \(L_J\).
# This is the core object for working with post processing data from HFSS
# This will load the saved solutions from above
epra = epr.QuantumAnalysis(eprh.data_filename)
# Analyze
epra.analyze_all_variations(cos_trunc = 8, fock_trunc = 15);
epra.plot_hamiltonian_results(swp_variable='Lj_1');
There is some jitter in the results since some of the runs were performed with only a few adaptive passes. In a production analysis one would ensure all variations have converged to the same \(\Delta f\) threshold before proceeding.
The sweep makes several physically important trends visible:
As \(L_J\) increases, \(E_J\) decreases, so the qubit frequency \(\omega_q \propto \sqrt{E_J}\) decreases.
The anharmonicity \(\alpha \approx -E_C\) is relatively insensitive to \(L_J\) in the transmon regime (since \(E_C\) is set by the capacitor geometry, not the junction).
The dispersive coupling \(\chi_{qc}\) decreases as \(L_J\) increases (equivalently, as \(E_J / E_C\) increases and the qubit moves deeper into the transmon regime).
See Figure 3 of arXiv:2010.00620 for an example of these trends in a real device.
epra.quick_plot_mode(0,0,1,numeric=True,swp_variable='Lj_1')
epra.report_results(swp_variable='Lj_1', numeric=True)
Addendum:
Snippets for ease #
An easy way to define many junctions is the following, using Python 3 f-strings for ease. This pattern is useful for multi-junction circuits such as flux qubits, SQUID arrays, or fluxonium qubits, where each junction must be registered individually but shares a common naming convention:
# Specify multiple junctions in HFSS model
n_junctions = 5
for i in range(1, 1+n_junctions):
pinfo.junctions[f'j{i}'] = {'Lj_variable' : f'Lj{i}',
'rect' : f'JJrect{i}',
'line' : f'JJline{i}',
'length' : epr.parse_units('100um')}
NEXT: Please see Tutorial 2 to continue. It covers the calculation of dielectric energy-participation ratios, which are needed to estimate qubit lifetimes \(T_1\) from dielectric loss tangents.
Snippets for ease
An easy way to define many junctions is the following, using Python 3 f-strings for ease. This pattern is useful for multi-junction circuits such as flux qubits, SQUID arrays, or fluxonium qubits, where each junction must be registered individually but shares a common naming convention:
# Specify multiple junctions in HFSS model
n_junctions = 5
for i in range(1, 1+n_junctions):
pinfo.junctions[f'j{i}'] = {'Lj_variable' : f'Lj{i}',
'rect' : f'JJrect{i}',
'line' : f'JJline{i}',
'length' : epr.parse_units('100um')}
NEXT: Please see Tutorial 2 to continue. It covers the calculation of dielectric energy-participation ratios, which are needed to estimate qubit lifetimes \(T_1\) from dielectric loss tangents.