5. Maxwell eigenmode extraction for an Xmon qubit: inductor ports and EPR analysis

5.1. Requirements

5.1.1. Software components

• QTCAD

• Gmsh

5.1.2. Geometry file

• qtcad/examples/tutorials/meshes/xmon.py

5.1.3. Python script

• qtcad/examples/tutorials/maxwell_xmon.py

5.1.4. References

• Maxwell eigenmode extraction of a coplanar waveguide resonator with adaptive meshing

• Capacitance matrix extraction for an Xmon qubit

• Inductor boundaries (API reference)

• Energy-participation-ratio method

5.2. Briefing

In this tutorial, we explain how to use inductor boundaries in the Maxwell eigenmode Solver. The device considered is the Xmon [BKM+13] from the tutorial Capacitance matrix extraction for an Xmon qubit. Here, an inductor boundary, described in Inductor boundaries (API reference), is associated to the surface representing the Josephson junction; that is, we approximate it as a linear inductor. In general, such an analysis is relevant because the eigenmode extraction involving linearized Josephson circuits is an important step in the energy-participation quantization method [MLM+21]. The last section of this tutorial shows how to post-process the resulting Maxwell eigenmode with the energy-participation-ratio (EPR) module.

5.3. Xmon qubit under consideration

The operation of an Xmon qubit [BKM+13] is based on that of a transmon qubit [KYG+07]. This type of qubit is formed by a capacitor and a Josephson junction in parallel, as illustrated on the left of Fig. 5.3.4.

Fig. 5.3.4 Transmon qubit and the linearized LC oscillator approximation.

The Josephson junction is governed by the following equations [BGGW21]:

(5.3.1)\[i_{\text{JJ}}(t) = I_c \sin\left(\dfrac{2 \pi \Phi(t)}{\Phi_0}\right)\]

and

(5.3.2)\[\dfrac{d\Phi}{dt} = v_{\text{JJ}}(t),\]

where \(i_{\text{JJ}}\) and \(v_{\text{JJ}}\) are the current and voltage of the junction respectively and \(I_c\) is the critical current. Also, \(\Phi_0=h/(2e)\) is the superconducting flux quantum and \(\Phi(t)\) is the so-called flux variable, analogous to the magnetic flux in an inductor.

The capacitor and the Josephson junction form a weakly anharmonic oscillator, where, in the classical description, the capacitor charge and the flux variable oscillate around zero. The anharmonicity of the oscillator comes from the nonlinear nature of the Josephson junction, namely the non-linear relation Eq. (5.3.1) between the current and the flux variable.

The transmon exhibits quantized energy levels, with the separation between the first two levels given by [BGGW21]

(5.3.3)\[E_{01} = \sqrt{8 E_C E_J} - E_C,\]

where \(E_C=e^2/(2C)\) is the charging energy [1] and \(E_J = \Phi_0 I_c/(2\pi)\) is the Josephson energy [2].

The anharmonicity of the oscillator manifests itself in the unequal separation between energy levels, and is quantified as \(E_{12}-E_{01} = -E_C\) [BGGW21, KYG+07]. Here, \(E_{ij} \equiv E_j - E_i\) is the energy separation between the energy levels \(i\) and \(j\). For a transmon, the anharmonicity is much smaller than the energy level separation. Nevertheless, it is essential in ensuring that the quantum oscillator stays within the first two energy levels, allowing it to operate as a qubit.

5.3.1. Modeling Josephson junction as an inductor

The Josephson junction can be thought of as a nonlinear inductor, with the inductance given by [BGGW21]

(5.3.4)\[L = \dfrac{\Phi_0^2}{(2 \pi)^2 E_J} \dfrac{1}{\cos(2 \pi \Phi/\Phi_0)}.\]

For small oscillations of \(\Phi\), the Josephson junction can be linearized by setting \(\Phi\) to zero in Eq. (5.3.4), yielding an inductance that is independent of the flux variable:

(5.3.5)\[L \approx \dfrac{\Phi_0^2}{(2 \pi)^2 E_J}.\]

The energy level separation of an LC harmonic quantum oscillator is given by [BGGW21]

(5.3.6)\[E_{01, LC} = \hbar \omega_{LC},\]

where

(5.3.7)\[\omega_{LC} = \dfrac1{\sqrt {LC}}\]

is the resonant angular frequency of the classical LC harmonic oscillator.

5.4. Geometry of the problem

The geometry of the device is the same as in the tutorial Capacitance matrix extraction for an Xmon qubit. The layout and dimensions were kindly provided by Christopher Xu from Red Blue Quantum.

Fig. 5.4.9 Layout of the Xmon.

5.5. Setting up the device and finding the Maxwell eigenmodes

The procedure for finding the Maxwell eigenmodes is similar to that described in the tutorial Maxwell eigenmode extraction of a coplanar waveguide resonator with adaptive meshing, with the additional step of adding the inductor to the device.

5.5.1. Header, input parameters, and file paths

"""
Maxwell eigenmodes for an Xmon with the Josephson junction represented as a linear
inductor.

The layout and dimensions for this example were kindly provided by Christopher Xu from
Red Blue Quantum.

See the following article for more details on the operation of an Xmon qubit coupled to
a readout line, XY control line, and a quantum bus resonator.
    Barends, Rami, et al. "Coherent Josephson qubit suitable for scalable quantum
    integrated circuits." Phys. Rev. Lett., 111.8 (2013): 080502.
"""

from pathlib import Path
import os
from time import time
import numpy as np

# Import relevant modules of QTCAD.
from qtcad.device.maxwell_eigenmode import Solver
from qtcad.device.maxwell_eigenmode import SolverParams
from qtcad.device.device import Device
from qtcad.device.mesh3d import Mesh
from qtcad.device import materials as mt
from qtcad.device import constants as ct

The Solver class contains algorithms needed to find the Maxwell eigenmodes and the SolverParams class can be used to control parameters of the solver.

To proceed, let us define some general variables related to the output directory and files, as well as the mesh:

# Scale in the Gmsh files. That is, Gmsh file coordinates are in μm.
scale = 1e-6

# Directories and file paths.
script_dir = Path(__file__).parent.resolve()
# For mesh and raw geometry files.
input_dir = script_dir / "meshes"
# For results.
result_dir = script_dir / "output" / Path(__file__).stem
# Mesh file.
fpath_mesh = input_dir / "xmon.msh"
# Raw geometry file (needed for adaptive meshing).
fpath_xao = input_dir / "xmon.xao"

# Check if the mesh and raw geometry files exist.
if not os.path.isfile(fpath_mesh) or not os.path.isfile(fpath_xao):
    raise Exception(
        "Please run %s/xmon.py to generate the mesh and raw geometry files."
        % (input_dir)
    )

5.5.2. Loading the initial mesh and defining the device

Next, we parse the mesh file, initialize the device and assign the media to 3D physical groups and boundary conditions to 2D physical groups:

######################################################################################
# Setup the device.
######################################################################################
# Parse the mesh and initialize the device.
mesh = Mesh(scale, fpath_mesh)
dvc = Device(mesh)

material_sub = mt.Si
material_air = mt.vacuum
# Assign media to regions.
dvc.new_region("substrate", material_sub)
dvc.new_region("air", material_air)
# Assign perfect electric conductor boundary condition to all conductors.
dvc.new_pec_bnd("gnd")
dvc.new_pec_bnd("xmon_cross")
dvc.new_pec_bnd("xy_ctrl")
dvc.new_pec_bnd("readout")
dvc.new_pec_bnd("qbus")

Then, we add an inductor to represent the Josephson junction:

# Calculate inductance.
#
# For the expression of the inductance as a function of the flux, see
#   Alexandre Blais, Arne L. Grimsmo, S. M. Girvin, and Andreas Wallraff.
#   Circuit quantum electrodynamics. Rev. Mod. Phys., 93:025005, May 2021.
EJ = 23e9 * ct.h
Phi0 = ct.h / (2 * ct.e)
inductance = Phi0**2 / ((2 * np.pi) ** 2 * EJ)

# Add the inductor.
dvc.new_inductor("jj", inductance, dir="y", length=24e-6, width=24e-6)

Here, when adding the inductor with new_inductor, we specify dir="y" to indicate that the current of the inductor flows between the conductor of the Xmon superconducting island and the ground plane. Currently, we support the orthogonal directions "x", "y" and "z". Also, the dimensions length and width are the dimensions of the boundary of the inductor "jj". They must be passed and, in this example, can be found in the geometry file.

More information on inductor boundaries can be found in Inductor boundaries (API reference).

5.5.3. Creating and running the Maxwell eigenmode solver

######################################################################################
# Setup the solver.
######################################################################################
params = SolverParams()
# Number of modes to find.
params.num_modes = 1
# Adaptive meshing (relative) tolerance on the frequency.
params.tol_rel = 0.03
# Directory where output files will be stored
params.output_dir = result_dir

######################################################################################
# Run the solver (results are stored in the folder specified by params.output_dir).
######################################################################################
slv = Solver(dvc, params, geo_file=fpath_xao)
t0 = time()
slv.solve()
dt = time() - t0
print("Solution completed in %.2f s" % dt)

5.5.4. Extracting further solver results

Next, we can display the computed frequency, along with the analytical frequency of the the LC resonator formed by the Xmon capacitor and the linearized inductor. The latter is computed by dividing Eq. (5.3.7) by \(2 \pi\).

######################################################################################
# Extract additional information from the final results.
######################################################################################

# From the cap_xmon.py tutorial.
capacitance = 1.014e-13

freq_lc = 1 / np.sqrt(inductance * capacitance) / (2 * np.pi)
freq_lc_ghz = freq_lc / 1e9
print(f"Frequency of the LC resonator (analytical): {freq_lc_ghz:.3f} GHz")

freq_ghz = dvc.maxwell_freqs[0] / 1e9
print(f"Frequency of the fundamental mode (computed): {freq_ghz:.3f} GHz")

The resulting frequencies are

Frequency of the LC resonator (analytical): 5.929 GHz
Frequency of the fundamental mode (computed): 5.741 GHz

It should be noted that the LC resonator model is a good reference for comparison, as the size of the device is much smaller than the wavelength of electromagnetic waves at these resonant frequencies [3].

For completeness, we also show the frequency of the original qubit, computed from Eq. (5.3.3):

# Blais et al. (2021).
EC = ct.e**2 / (2 * capacitance)
freq_qbit = (np.sqrt(8 * EC * EJ) - EC) / ct.h
freq_qbit_ghz = freq_qbit / 1e9
print(f"Frequency of the qubit (analytical): {freq_qbit_ghz:.3f} GHz")

Frequency of the qubit (analytical): 5.738 GHz

The resulting fields from the xmon-fields.vtu file are shown in the figures below:

Fig. 5.5.4 Magnitude of the electric field viewed in ParaView. Clip filter at \(z=0\) is used to visualize the field at the level of the the ground plane. The values are displayed in logarithmic scale with the “Black, Blue and White” colour scheme.

Fig. 5.5.5 Magnitude of the magnetic flux density viewed in ParaView. Clip filter at \(z=0\) is used to visualize the field at the level of the the ground plane. The values are displayed in logarithmic scale with the “Black, Blue and White” colour scheme.

5.6. Energy-participation-ratio analysis

The above solution of the linearized-junction problem is a requirement for the use of the EPR method [MLM+21]. As explained in Energy-participation-ratio method, the Josephson potential can be expanded around zero reduced flux. The resulting quadratic term is associated with the linear inductor already used in the Maxwell eigenmode problem, while the leading nonlinear correction comes from the quartic term.

The EPR module uses the solved linear modes to compute the fraction of each mode’s inductive energy stored in each Josephson junction. For mode \(m\) and junction \(j\), this energy participation is

\[p_{mj} = \frac{ U_{Jj,m}^{\mathrm{lin}} }{ U_m^{\mathrm{ind}} },\]

where \(U_{Jj,m}^{\mathrm{lin}}\) is the linear inductive energy stored in the junction and \(U_m^{\mathrm{ind}}\) is the total inductive energy of the mode.

For a mode frequency \(f_m\), this energy participation determines the reduced-flux zero-point fluctuation through

\[\varphi_{mj} = s_{mj} \sqrt{ p_{mj}\frac{h f_m}{2E_{Jj}} },\]

where \(s_{mj}=\pm 1\) is set by the current orientation.

From these quantities, QTCAD^® extracts dressed frequencies, Kerr coefficients and the anharmonicity of the qubit mode via EPRAnalysis, which must be used after the Maxwell eigenmode solver has been run. The usage of EPRAnalysis is demonstrated below.

5.6.1. Importing the EPR classes and identifying junctions

In addition to the modules imported earlier, let us import the relevant classes for EPR analysis: EPRAnalysis and JunctionSpec, as well as DielectricSpec and SurfaceDielectricSpec.

from qtcad.device.epr import EPRAnalysis
from qtcad.device.epr import JunctionSpec, DielectricSpec, SurfaceDielectricSpec

The JunctionSpec object identifies one Josephson junction boundary. Its boundary name must match the physical group’s name used to set up the inductor port in the eigenmode solver. In this example, we have a single junction whose boundary name is "jj". Let us then create a list with the junctions on which we would like to run the EPR analysis.

epr_junctions = [JunctionSpec("jj")]

Similarly, the DielectricSpec object represents a dielectric region on which we want to evaluate bulk losses. It requires the physical region name (for instance "substrate") and an optional loss tangent, which should be non-zero if one wants to study bulk dielectric losses (for more information, check Dielectric-loss-limited quality factors):

spec_dielectrics = [
    DielectricSpec("substrate", loss_tangent=5e-7),
]

To evaluate losses from thin oxide or other contaminant layers sandwiched between distinct components or regions [WAG+15, WBB+11], we can use the SurfaceDielectricSpec class. It allows the modelling of losses due to thin dielectric layers at different interfaces such as metal–air (MA), metal–substrate (MS) and substrate–air (SA). This class takes the following arguments:

• surface: The name (or list of names) of the boundary physical group(s).

• thickness: The thickness of the thin interface layer.

• relative_permittivity: The relative permittivity of the thin-film layer.

• loss_tangent: The loss tangent of the thin-film layer.

• interface_type: The physical model to utilize. The options are "ma", "ms" and "sa" for the MA, MS and SA interface models, respectively. For more information on the different models, please check Dielectric-loss-limited quality factors.

Here, we perform a simplified analysis and we define the different surface-loss specifications based on the appropriate physical groups in the mesh as follows:

spec_ms = SurfaceDielectricSpec(
    "xmon_cross",
    thickness=3e-9,
    relative_permittivity=10.0,
    loss_tangent=7e-4,
    interface_type="ms"
)
spec_ma = SurfaceDielectricSpec(
    "xmon_cross",
    thickness=3e-9,
    relative_permittivity=10.0,
    loss_tangent=4e-3,
    interface_type="ma"
)
spec_sa = SurfaceDielectricSpec(
    "substrate_top",
    thickness=3e-9,
    relative_permittivity=4,
    loss_tangent=6e-4,
    interface_type="sa"
)

5.6.2. Running the EPR analysis

We can now run the EPR analysis, passing the list of junctions and our bulk and surface loss specifications directly to the EPRAnalysis constructor.

Note

EPRAnalysis allows evaluating multiple surface-dieletric-loss models using the same physical group.

This is relevant for simulations where the superconducting sheet is approximated as a 2D surface and we perform analyses using the MA and MS models, where they may share the same physical boundary.

In the current example, this is the case for the physical group "xmon_cross", which is used in the set up of the MA and the MS models.

######################################################################################
# Solve EPR analysis.
######################################################################################
# The EPR analysis reads the solved Maxwell fields and the current through the
# inductor boundary named "jj", then converts the linear EPR data into
# first-order nonlinear quantities such as the dressed frequency, zero-point
# fluctuations and anharmonicity.
epr = EPRAnalysis(
    dvc,
    epr_junctions,
    dielectrics=spec_dielectrics,
    surface_dielectrics=[spec_ms, spec_ma, spec_sa],
).solve()

Here, the first positional argument is the solved device dvc. The second positional argument is the list of Josephson junctions to include in the EPR analysis.

Since we only solved for a single Maxwell eigenmode, the resulting arrays will have one mode row and one junction column. For example, epr.participation[0, 0] is the junction participation of the fundamental mode and epr.phi_zpf[0, 0] is the corresponding reduced-flux zero-point fluctuation.

5.6.3. Comparing EPR and analytical quantities

Let us then print a compact comparison table of some relevant quantities of interest, including the relaxation time \(T_1\) obtained from the dielectric-loss-limited quality factors and the dressed transmon frequency. Note that the EPR-specific rows are obtained from the EPRResult returned by solve.

######################################################################################
# EPR results and comparisons with simple analytical estimates.
######################################################################################
# Calculate relaxation time T1 from the total quality factor.
freq_dressed = epr.dressed_frequencies[0]
q_total = epr.mode_quality_factors[0]
t1 = q_total / (2 * np.pi * freq_dressed)

rows = [
    ("Frequency of the LC resonator (analytical)", f"{freq_lc / 1e9:.3f} GHz"),
    ("Frequency of the qubit (analytical)", f"{freq_qbit / 1e9:.3f} GHz"),
    ("Bare Maxwell frequency (computed)", f"{epr.bare_frequencies[0] / 1e9:.3f} GHz"),
    ("Dressed qubit frequency (EPR)", f"{freq_dressed / 1e9:.3f} GHz"),
    ("Junction participation", f"{epr.participation[0, 0]:.6f}"),
    ("Reduced flux ZPF", f"{epr.phi_zpf[0, 0]:.6f}"),
    ("Anharmonicity", f"{epr.anharmonicity[0] / 1e6:.3f} MHz"),
    ("Bulk substrate quality factor", f"{epr.dielectric_quality_factors[0, 0]:.3e}"),
    ("Surface MS quality factor", f"{epr.surface_dielectric_quality_factors[0, 0]:.3e}"),
    ("Surface MA quality factor", f"{epr.surface_dielectric_quality_factors[0, 1]:.3e}"),
    ("Surface SA quality factor", f"{epr.surface_dielectric_quality_factors[0, 2]:.3e}"),
    ("Total quality factor", f"{q_total:.3e}"),
    ("Relaxation time T1", f"{t1 * 1e6:.3f} μs"),
]

# Display data as a table.
header = f"{'Quantity':<45} | {'Value':>15}"
print("\n" + header)
print("-" * len(header))
for label, value in rows:
    print(f"{label:<45} | {value:>15}")

We obtain:

Quantity                                      |           Value
---------------------------------------------------------------
Frequency of the LC resonator (analytical)    |       5.929 GHz
Frequency of the qubit (analytical)           |       5.738 GHz
Bare Maxwell frequency (computed)             |       5.720 GHz
Dressed qubit frequency (EPR)                 |       5.547 GHz
Junction participation                        |        0.985526
Reduced flux ZPF                              |        0.350073
Anharmonicity                                 |    -172.717 MHz
Bulk substrate quality factor                 |       2.172e+06
Surface MS quality factor                     |       1.091e+07
Surface MA quality factor                     |       2.620e+08
Surface SA quality factor                     |       2.005e+07
Total quality factor                          |       1.651e+06
Relaxation time T1                            |       47.369 μs

Here,

• The bare frequency is the linearized Maxwell eigenfrequency before the quartic Josephson correction.

• The dressed frequency is the first-order EPR frequency after subtracting the Lamb shift associated with the junction nonlinearity.

• The participation is the normalized fraction of the mode’s inductive energy stored in the Josephson junction.

• The reduced-flux zero-point fluctuation (ZPF) is the coefficient relating the junction reduced-flux operator to the linear mode operators

• The anharmonicity is the resulting leading-order shift of the qubit transition spacing.

• The bulk and surface quality factors represent the individual loss contributions from the respective region (Si substrate) and thin layers of contaminants.

• The total quality factor \(Q_{\text{total}}\), epr.mode_quality_factors[0], represents the combined dielectric-loss-limited quality factor, satisfying:

  \[\frac{1}{Q_{\text{total}}} = \sum_{d} \frac{1}{Q_{\text{dielectric}, d}} + \sum_{s} \frac{1}{Q_{\text{surface}, s}}.\]

• The relaxation time \(T_1\) is the qubit’s relaxation time limit due to dielectric loss, computed using

  \[T_1 = \frac{Q_{\text{total}}}{2 \pi \times (\text{dressed qubit frequency})}.\]

The above EPR analysis provides an efficient way of bridging classical electromagnetic simulations and circuit QED, where we have quantum non-linear inductive components (the Josephson junctions). From the results of linear Maxwell-eigenmode simulations, EPR offers a scalable approach to quantum circuit quantization.

5.7. Full code

__copyright__ = "Copyright 2022-2026, Nanoacademic Technologies Inc."

"""
Maxwell eigenmodes for an Xmon with the Josephson junction represented as a linear
inductor.
The layout and dimensions for this example were kindly provided by Christopher Xu from
Red Blue Quantum.
See the following article for more details on the operation of an Xmon qubit coupled to
a readout line, XY control line, and a quantum bus resonator.
    Barends, Rami, et al. "Coherent Josephson qubit suitable for scalable quantum
    integrated circuits." Phys. Rev. Lett., 111.8 (2013): 080502.
"""
from pathlib import Path
import os
from time import time
import numpy as np

# Import relevant modules of QTCAD.
from qtcad.device.maxwell_eigenmode import Solver
from qtcad.device.maxwell_eigenmode import SolverParams
from qtcad.device.device import Device
from qtcad.device.mesh3d import Mesh
from qtcad.device import materials as mt
from qtcad.device import constants as ct

# Scale in the Gmsh files. That is, Gmsh file coordinates are in μm.
scale = 1e-6

# Directories and file paths.
script_dir = Path(__file__).parent.resolve()
# For mesh and raw geometry files.
input_dir = script_dir / "meshes"
# For results.
result_dir = script_dir / "output" / Path(__file__).stem
# Mesh file.
fpath_mesh = input_dir / "xmon.msh"
# Raw geometry file (needed for adaptive meshing).
fpath_xao = input_dir / "xmon.xao"

# Check if the mesh and raw geometry files exist.
if not os.path.isfile(fpath_mesh) or not os.path.isfile(fpath_xao):
    raise Exception(
        "Please run %s/xmon.py to generate the mesh and raw geometry files."
        % (input_dir)
    )

######################################################################################
# Setup the device.
######################################################################################
# Parse the mesh and initialize the device.
mesh = Mesh(scale, fpath_mesh)
dvc = Device(mesh)

material_sub = mt.Si
material_air = mt.vacuum
# Assign media to regions.
dvc.new_region("substrate", material_sub)
dvc.new_region("air", material_air)
# Assign perfect electric conductor boundary condition to all conductors.
dvc.new_pec_bnd("gnd")
dvc.new_pec_bnd("xmon_cross")
dvc.new_pec_bnd("xy_ctrl")
dvc.new_pec_bnd("readout")
dvc.new_pec_bnd("qbus")


# Calculate inductance.
#
# For the expression of the inductance as a function of the flux, see
#   Alexandre Blais, Arne L. Grimsmo, S. M. Girvin, and Andreas Wallraff.
#   Circuit quantum electrodynamics. Rev. Mod. Phys., 93:025005, May 2021.
EJ = 23e9 * ct.h
Phi0 = ct.h / (2 * ct.e)
inductance = Phi0**2 / ((2 * np.pi) ** 2 * EJ)

# Add the inductor.
dvc.new_inductor("jj", inductance, dir="y", length=24e-6, width=24e-6)

######################################################################################
# Setup the solver.
######################################################################################
params = SolverParams()
# Number of modes to find.
params.num_modes = 1
# Adaptive meshing (relative) tolerance on the frequency.
params.tol_rel = 0.03
# Directory where output files will be stored
params.output_dir = result_dir

######################################################################################
# Run the solver (results are stored in the folder specified by params.output_dir).
######################################################################################
slv = Solver(dvc, params, geo_file=fpath_xao)
t0 = time()
slv.solve()
dt = time() - t0
print("Solution completed in %.2f s" % dt)


######################################################################################
# Extract additional information from the final results.
######################################################################################

# From the cap_xmon.py tutorial.
capacitance = 1.014e-13

freq_lc = 1 / np.sqrt(inductance * capacitance) / (2 * np.pi)
freq_lc_ghz = freq_lc / 1e9
print(f"Frequency of the LC resonator (analytical): {freq_lc_ghz:.3f} GHz")

freq_ghz = dvc.maxwell_freqs[0] / 1e9
print(f"Frequency of the fundamental mode (computed): {freq_ghz:.3f} GHz")

# Blais et al. (2021).
EC = ct.e**2 / (2 * capacitance)
freq_qbit = (np.sqrt(8 * EC * EJ) - EC) / ct.h
freq_qbit_ghz = freq_qbit / 1e9
print(f"Frequency of the qubit (analytical): {freq_qbit_ghz:.3f} GHz")

from qtcad.device.epr import EPRAnalysis
from qtcad.device.epr import JunctionSpec, DielectricSpec, SurfaceDielectricSpec

epr_junctions = [JunctionSpec("jj")]

spec_dielectrics = [
    DielectricSpec("substrate", loss_tangent=5e-7),
]

spec_ms = SurfaceDielectricSpec(
    "xmon_cross",
    thickness=3e-9,
    relative_permittivity=10.0,
    loss_tangent=7e-4,
    interface_type="ms",
)
spec_ma = SurfaceDielectricSpec(
    "xmon_cross",
    thickness=3e-9,
    relative_permittivity=10.0,
    loss_tangent=4e-3,
    interface_type="ma",
)
spec_sa = SurfaceDielectricSpec(
    "substrate_top",
    thickness=3e-9,
    relative_permittivity=4,
    loss_tangent=6e-4,
    interface_type="sa",
)

######################################################################################
# Solve EPR analysis.
######################################################################################
# The EPR analysis reads the solved Maxwell fields and the current through the
# inductor boundary named "jj", then converts the linear EPR data into
# first-order nonlinear quantities such as the dressed frequency, zero-point
# fluctuations and anharmonicity.
epr = EPRAnalysis(
    dvc,
    epr_junctions,
    dielectrics=spec_dielectrics,
    surface_dielectrics=[spec_ms, spec_ma, spec_sa],
).solve()

######################################################################################
# EPR results and comparisons with simple analytical estimates.
######################################################################################
# Calculate relaxation time T1 from the total quality factor.
freq_dressed = epr.dressed_frequencies[0]
q_total = epr.mode_quality_factors[0]
t1 = q_total / (2 * np.pi * freq_dressed)

rows = [
    ("Frequency of the LC resonator (analytical)", f"{freq_lc / 1e9:.3f} GHz"),
    ("Frequency of the qubit (analytical)", f"{freq_qbit / 1e9:.3f} GHz"),
    ("Bare Maxwell frequency (computed)", f"{epr.bare_frequencies[0] / 1e9:.3f} GHz"),
    ("Dressed qubit frequency (EPR)", f"{freq_dressed / 1e9:.3f} GHz"),
    ("Junction participation", f"{epr.participation[0, 0]:.6f}"),
    ("Reduced flux ZPF", f"{epr.phi_zpf[0, 0]:.6f}"),
    ("Anharmonicity", f"{epr.anharmonicity[0] / 1e6:.3f} MHz"),
    ("Bulk substrate quality factor", f"{epr.dielectric_quality_factors[0, 0]:.3e}"),
    (
        "Surface MS quality factor",
        f"{epr.surface_dielectric_quality_factors[0, 0]:.3e}",
    ),
    (
        "Surface MA quality factor",
        f"{epr.surface_dielectric_quality_factors[0, 1]:.3e}",
    ),
    (
        "Surface SA quality factor",
        f"{epr.surface_dielectric_quality_factors[0, 2]:.3e}",
    ),
    ("Total quality factor", f"{q_total:.3e}"),
    ("Relaxation time T1", f"{t1 * 1e6:.3f} μs"),
]

# Display data as a table.
header = f"{'Quantity':<45} | {'Value':>15}"
print("\n" + header)
print("-" * len(header))
for label, value in rows:
    print(f"{label:<45} | {value:>15}")

[1]

In order to account for the capacitance of the Josephson junction, one would add it to \(C\) [BGGW21].

[2]

It should be noted that in [BKM+13], there are two parallel Josephson junctions forming a superconducting quantum interference device (SQUID), which permits the tunability of \(E_J\). For simplicity, we only consider the case of a single Josephson junction.

[3]

At 5.929 GHz, the wavelength is 5.1 cm in vacuum and 1.5 cm in silicon.