5. Energy-participation-ratio method

This section introduces the energy-participation-ratio (EPR) method for quantizing superconducting circuits containing Josephson junctions [MLM+21]. The method provides a practical way to incorporate Josephson nonlinearities into the quantum Hamiltonian of a superconducting circuit, starting from linear electromagnetic simulations. These simulations yield the circuit’s linear normal modes, which are promoted to quantum harmonic modes described by bosonic creation and annihilation operators \(\hat a_m\) and \(\hat a_m^\dagger\). The nonlinear Josephson relation is expressed in terms of the reduced junction flux. In the EPR formalism, each reduced junction flux is written as a sum of contributions from the linear normal modes, with weights determined from classical electromagnetic data. These weights connect the simulated linear modes to the many-body quantum Hamiltonian of the circuit.

Linearized Josephson circuit

To obtain the linear modes used in the EPR formalism, each Josephson junction is first approximated by its linear inductive response.

For a Josephson junction \(j\), the potential energy may be written as [BGGW21]

(5.41)\[U_{J}^j(\varphi_j) = -E_{J}^j\cos\varphi_j, \qquad \varphi_j = \frac{\Phi_j}{\left[\Phi_0/(2\pi)\right]} = \frac{\Phi_j}{\left[\hbar/(2e)\right]},\]

where \(E_{J}^j\) is the Josephson energy, \(\varphi_j\) is the reduced junction flux, \(\Phi_j\) is the flux across the junction, and \(\Phi_0=h/(2e)\) is the superconducting flux quantum. Expanding the expression for \(U_{J}^j(\varphi_j)\) in Eq. (5.41) around \(\varphi_j=0\) gives

(5.42)\[-E_{J}^j\cos\varphi_j = -E_{J}^j + \frac{E_{J}^j}{2}\varphi_j^2 - \frac{E_{J}^j}{24}\varphi_j^4 + \cdots .\]

The quadratic term has the same form as the energy of a linear inductor. This allows each Josephson junction to be associated with an effective linear inductance \(L_{J}^j\) as follows:

(5.43)\[\frac{E_{J}^j}{2}\varphi_j^2 = \frac{\Phi_j^2}{2L_{J}^j}.\]

yielding

(5.44)\[L_{J}^j = \frac{\left[\Phi_0/(2\pi)\right]^2}{E_{J}^j} = \frac{\left[\hbar/(2e)\right]^2}{E_{J}^j}.\]

By replacing each Josephson junction \(j\) by a linear inductor \(L_J^j\), the device becomes a linear electromagnetic system. Quantizing this linear system yields the Hamiltonian

(5.45)\[\hat H_\mathrm{lin} = \sum_m \hbar\omega_m \left(\hat a_m^\dagger \hat a_m + \frac{1}{2}\right),\]

where \(\omega_m\) is the angular frequency of the \(m\)-th linearized normal mode, and \(\hat a_m^\dagger\) and \(\hat a_m\) are creation and annihilation operators for that mode. In QTCAD^®, the normal-mode angular frequencies \(\omega_m\) are computed with the Maxwell eigenmode Solver. In this linearized problem, the Josephson junctions are represented by inductor boundaries. For additional information, see the Maxwell Eigenmodes theory section.

Once the eigenmodes have been computed, the remaining task is to account for the higher-order terms in Eq. (5.42). This requires expressing each reduced junction flux operator \(\hat\varphi_j\) in terms of the creation and annihilation operators \(\hat a_m^\dagger\) and \(\hat a_m\) of the linearized normal modes, as described in the next section.

Hamiltonian expansion

The reduced flux operator of junction \(j\) can be expanded in the linearized-mode basis as

(5.46)\[\hat\varphi_j = \sum_m \varphi_{mj} \left(\hat a_m+\hat a_m^\dagger\right),\]

where each coefficient \(\varphi_{mj}\) is the zero-point fluctuation (ZPF) of the reduced junction flux associated with mode \(m\). These coefficients are the central local quantities needed for the nonlinear Hamiltonian [MLM+21, NPV+12, YD84].

The quadratic Josephson term has the same form as the other linear inductive terms and is included in the linearized normal modes. Under the weakly nonlinear assumption, the leading remaining contribution is therefore the fourth-order term in Eq. (5.42). For transmon-like circuits [KYG+07], higher-order terms are usually smaller, so the nonlinear Hamiltonian is approximated by

(5.47)\[\hat H_4 = -\sum_j \frac{E_{J}^j}{24} \left[ \sum_m \varphi_{mj} \left(\hat a_m+\hat a_m^\dagger\right) \right]^4 + \mathcal O(\hat\varphi_j^6).\]

Equation (5.47) is already a many-body Hamiltonian: expanding the fourth power produces products of operators from one mode, two modes, or more modes.

Kerr and cross-Kerr parameters

When the zero-point fluctuations are small and the modes are well separated in frequency, the quartic Hamiltonian may be reduced to the terms that conserve the occupation of each mode. This is the rotating-wave approximation (RWA). The names Kerr and cross-Kerr come from nonlinear optics, where the Kerr effect describes an intensity-dependent refractive index.

In a quantized circuit, the analogous effect is an occupation-dependent mode frequency: self-Kerr gives the frequency shift of a mode due to its own occupation, while cross-Kerr gives the frequency shift of one mode due to the occupation of another mode. The reduction can be carried out explicitly from Eq. (5.47). First write

(5.48)\[\hat X_m=\hat a_m+\hat a_m^\dagger, \qquad \hat n_m=\hat a_m^\dagger\hat a_m.\]

For one junction, the fourth power in Eq. (5.47) contains

(5.49)\[\left(\sum_m\varphi_{mj}\hat X_m\right)^4 = \sum_m \varphi_{mj}^4\hat X_m^4 + 6\sum_{m<n}\varphi_{mj}^2\varphi_{nj}^2 \hat X_m^2\hat X_n^2 +\cdots .\]

The omitted terms contain an odd power of at least one mode operator, or products of three or four different modes. Under the RWA, these terms oscillate rapidly and do not contribute to the leading Kerr Hamiltonian.

The number-conserving part of the expansion of (5.49) with a single mode index is

(5.50)\[\hat X_m^4 \rightarrow 6\hat n_m(\hat n_m-1)+12\hat n_m+3.\]

Similarly, for terms involving two different mode indices,

(5.51)\[\hat X_m^2\hat X_n^2 \rightarrow 4\hat n_m\hat n_n +2\hat n_m +2\hat n_n +1, \qquad m\ne n.\]

Here, the arrow means that the operator has been normal-ordered and only the number-conserving terms have been retained. Substituting Eqs. (5.50) and (5.51) into Eq. (5.47) gives

(5.52)\[\begin{split}\begin{aligned} \frac{\hat H_{4,\mathrm{RWA}}}{\hbar} &= -\sum_{j,m} \frac{E_{J}^j}{24\hbar} \varphi_{mj}^4 \left[ 6\hat n_m(\hat n_m-1)+12\hat n_m+3 \right] \\ &\quad - \sum_{j,m<n} \frac{E_{J}^j}{4\hbar} \varphi_{mj}^2\varphi_{nj}^2 \left[ 4\hat n_m\hat n_n+2\hat n_m+2\hat n_n+1 \right]. \end{aligned}\end{split}\]

Equation (5.52) contains three physically distinct operator structures. Terms proportional to \(\hat n_m(\hat n_m-1)\) give the self-Kerr contribution, from which the anharmonicity is derived. Terms proportional to \(\hat n_m\hat n_n\) give the cross-Kerr interaction. Terms proportional to a single \(\hat n_m\) shift the linear mode frequency. Constant terms only shift the zero of energy and are usually omitted.

It is convenient to collect the coefficients of these terms in the symmetric matrix

(5.53)\[\chi^{(0)}_{mn} = \sum_j \frac{E_{J}^j}{\hbar} \varphi_{mj}^2\varphi_{nj}^2,\]

which is the leading quartic Kerr matrix before separating diagonal and off-diagonal conventions. In terms of this matrix, Eq. (5.52) becomes

(5.54)\[\frac{\hat H_{4,\mathrm{RWA}}}{\hbar} = - \sum_m \frac{\chi^{(0)}_{mm}}{4} \hat n_m(\hat n_m-1) - \sum_{m<n} \chi^{(0)}_{mn} \hat n_m\hat n_n - \frac{1}{2} \sum_{m,n} \chi^{(0)}_{mn}\hat n_m +\mathrm{const}.\]

Adding the linear Hamiltonian in Eq. (5.45) and dropping the constant energy offset gives the standard Kerr form

(5.55)\[\frac{\hat H_\mathrm{eff}}{\hbar} = \sum_m (\omega_m-\Delta_m)\hat n_m - \sum_m \frac{\alpha_m}{2}\hat n_m(\hat n_m-1) - \sum_{m<n} \chi_{mn}\hat n_m\hat n_n,\]

with

(5.56)\[\alpha_m = \frac{1}{2}\chi^{(0)}_{mm}, \qquad \chi_{mn} = \chi^{(0)}_{mn}\quad(m\ne n), \qquad \Delta_m = \frac{1}{2}\sum_n \chi^{(0)}_{mn}.\]

We denote:

• \(\alpha_m\) the self-Kerr coefficient, or anharmonicity magnitude, of mode \(m\),

• \(\chi_{mn}\) the cross-Kerr frequency shift between modes \(m\) and \(n\),

• \(\Delta_m\) the Lamb shift of mode \(m\),

• \(\omega_m-\Delta_m\) the dressed frequency of mode \(m\) when the other modes are in their ground states.

The remaining task is to determine the zero-point phase amplitudes \(\varphi_{mj}\). Directly extracting these ZPFs from an electromagnetic simulation is inconvenient, because the field amplitudes carry an arbitrary overall normalization. The EPR method removes this ambiguity by replacing the absolute field scale with a dimensionless energy-participation ratio, as described next.

Energy-participation ratio and zero-point fluctuations

The EPR method uses a classical energy ratio to compute the ZPFs \(\varphi_{mj}\) [MLM+21]. For mode \(m\) and junction \(j\), the Josephson energy-participation ratio is

(5.57)\[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 junction \(j\) when mode \(m\) is excited, and \(U_m^{\mathrm{ind}}\) is the total inductive energy of that same mode. The total inductive energy includes the magnetic-field energy and the linear inductive energy stored in all lumped inductive elements.

For example, with a peak-current convention [1], we have that

(5.58)\[U_{Jj,m}^{\mathrm{lin}} = \frac{1}{4}L_{Jj}I_{mj}^2,\]

where \(I_{mj}\) is the junction peak current output by QTCAD^® for mode \(m\).

The relation between EPR and the ZPFs follows from applying the same energy ratio to the quantized oscillator. Suppose only mode \(m\) is considered in Eq. (5.46), so that \(\hat\varphi_j=\varphi_{mj}\hat X_m\). The linear inductive energy stored in junction \(j\) is then

(5.59)\[\left\langle \hat U_{Jj,m}^{\mathrm{lin}}\right\rangle = \frac{E_{J}^j}{2} \varphi_{mj}^2 \left\langle \hat X_m^2 \right\rangle.\]

For a Fock state \(|n_m\rangle\), \(\left\langle \hat X_m^2 \right\rangle=2n_m+1\). The total inductive energy of a harmonic oscillator is one half of its total energy,

(5.60)\[\left\langle \hat U_m^\mathrm{ind}\right\rangle = \frac{1}{2}\hbar\omega_m \left(n_m+\frac{1}{2}\right) = \frac{\hbar\omega_m}{4}(2n_m+1).\]

Therefore,

(5.61)\[p_{mj} = \frac{ \left\langle \hat U_{Jj,m}^{\mathrm{lin}}\right\rangle }{ \left\langle \hat U_m^\mathrm{ind}\right\rangle } = \frac{2E_{J}^j\varphi_{mj}^2}{\hbar\omega_m}.\]

Solving Eq. (5.61) for \(\varphi_{mj}\) gives the EPR-to-ZPF relation

(5.62)\[\varphi_{mj} = s_{mj} \sqrt{ p_{mj}\frac{\hbar\omega_m}{2E_{J}^j} },\]

where \(s_{mj}=\pm 1\) is a sign set by the relative orientation of the mode current through the junction. The sign is not relevant for ordinary quartic Kerr magnitudes, which depend on even powers of \(\varphi_{mj}\), but it is part of the full flux expansion in Eq. (5.46).

From EPR to nonlinear parameters

EPR becomes useful because Eq. (5.62) expresses each zero-point fluctuation in terms of three linear quantities: the mode frequency, the Josephson energy, and the classical participation ratio.

Substituting Eq. (5.62) into Eq. (5.53) gives

(5.63)\[\chi^{(0)}_{mn} = \sum_j \frac{\hbar\omega_m\omega_n}{4E_{J}^j} p_{mj}p_{nj}.\]

This expression shows the central physical intuition of EPR. A mode has a large self-Kerr when it stores a large fraction of its inductive energy in nonlinear Josephson elements. Two modes have a large cross-Kerr interaction only when they both participate in the same Josephson junctions:

(5.64)\[\chi_{mn} \propto \sum_j p_{mj}p_{nj}.\]

Thus, EPR translates the spatial distribution of classical inductive energy into the nonlinear couplings of the quantum Hamiltonian.

Dielectric-loss-limited quality factors

In addition to accounting for nonlinearities arising from the quantum nature of superconducting circuits containing Josephson junctions via the EPR method, we can use the electric-field distribution obtained from electromagnetic simulations to estimate the contributions of different energy-loss mechanisms and compute dielectric-loss-limited quality factors.

Typically, dielectric losses are parameterized by the loss tangent \(\tan(\delta_i)\) of a given material \(i\), where \(\delta_i\) is the associated loss angle. Assuming independent dissipation channels, the total quality factor \(Q_m\) of mode \(m\) is determined by the sum of the contributions from all lossy regions:

(5.65)\[\frac{1}{Q_m} = \sum_i \frac{1}{Q_{m,i}} = \sum_i p_{m,i} \tan(\delta_i),\]

where \(p_{m,i}\) is the dielectric participation ratio of region \(i\) for mode \(m\). It is defined as the ratio of the time-averaged electric energy stored in region \(i\) to the total time-averaged electric energy of the mode:

(5.66)\[p_{m,i} = \frac{U_{m,i}^{\text{E}}}{U_{m,\text{total}}^{\text{E}}}.\]

Bulk dielectric losses

For a bulk dielectric volume \(V_i\), a physical group in the mesh, with relative permittivity \(\varepsilon_i\), the stored energy is computed as:

(5.67)\[U_{m,i}^{\text{E}} = \frac{1}{4} \int_{V_i} \varepsilon_0 \varepsilon_i \mathbf{E}_m^{*} \cdot \mathbf{E}_m \mathrm{d}V,\]

where \(\mathbf{E}_m\) is the peak electric field of mode \(m\).

To account for bulk dielectric losses in QTCAD^®, users should instantiate DielectricSpec for the different dielectric regions and use these objects when setting up EPRAnalysis.

Surface dielectric losses

For highly coherent transmon qubits, bulk dielectric losses are often negligible compared to losses originating from thin amorphous defect (oxide) layers at the interfaces [WAG+15, WBB+11] between the different domains (regions) of a superconducting-circuit chip.

Because these layers are usually a few nanometres thick, much thinner than the typical scale 3D finite-element meshes used for EM simulations, their participation can be calculated using the thin-film perturbation approximation described by Wenner et al. [WBB+11].

Here, \(U_{m,i}^{\text{E}}\) becomes \(U_{m,i}^{\text{surf}}\), the energy stored in a thin surface layer of thickness \(t_i\) and relative permittivity \(\varepsilon_i\) resting on a surface \(S_i\). It is given by

(5.68)\[U_{m,i}^{\text{surf}} = \frac{1}{4} t_j \int_{S_i} \varepsilon_0 \varepsilon_i \mathbf{E}_{m,i}^{*} \cdot \mathbf{E}_{m,i} \mathrm{d}S,\]

where \(\mathbf{E}_{m,i}\) is the electric field inside the thin film. In fact, this energy term can be further specialized by imposing the relevant boundary conditions at the different types of interfaces [WBB+11].

Extracting this internal field from the macroscopic simulation depends heavily on the type of physical interface being modelled and QTCAD^® supports three interface models to effectively model these contaminant layers:

• Metal-air (MA) model:

  At a boundary between a perfect electric conductor (PEC) and air, the tangential electric field is zero. Let us consider a sufficiently thin metal-oxide (contaminant) layer resting on a surface \(S_{\text{ma}}\) and with relative permittivity \(\varepsilon_{\text{ma}}\). The field inside it can be considered normal, with a negligible tangential component, and is determined by the continuity of the normal displacement field at the metal-oxide and air interface. This leads to \(\mathbf{E}_{\text{ma}} \approx \mathbf{E}_{\text{air},\perp}/\varepsilon_{\text{ma}}\). Hence, for a contaminant layer of thickness \(t_{\text{ma}}\), the stored energy becomes

  (5.69)\[ U_{\text{ma}} = \frac{1}{4} \varepsilon_0 t_{\text{ma}} \int_{S_{\text{ma}}} \frac{1}{\varepsilon_{\text{ma}}} |\mathbf{E}_{\text{air},\perp}|^2 \mathrm{d}S.\]

• Metal-substrate (MS) model:

  In the MS model, the tangential electric field is also zero at the boundary of the contaminant (oxide) layer resting on a surface \(S_{\text{ms}}\) and with relative permittivity \(\varepsilon_{\text{ms}}\), between a PEC and a substrate with relative permittivity \(\varepsilon_\text{substrate}\). The continuity of the normal displacement field leads to the following relation between the normal electric field inside the oxide and the substrate electric field: \(\mathbf{E}_{\text{ms}} \approx \varepsilon_\text{substrate} \mathbf{E}_{\text{substrate},\perp}/\varepsilon_{\text{ms}}\). Thus, for a contaminant layer of thickness \(t_{\text{ms}}\), the stored energy becomes

  (5.70)\[ U_{\text{ms}} = \frac{1}{4} \varepsilon_0 t_{\text{ms}} \int_{S_{\text{ms}}} \frac{\varepsilon_\text{substrate}^2}{\varepsilon_{\text{ms}}} |\mathbf{E}_{\text{substrate},\perp}|^2 \mathrm{d}S.\]

• Substrate-air (SA) model:

  At the interface between a dielectric substrate and air formed by a contaminant thin layer resting on a surface \(S_{\text{sa}}\) and with relative permittivity \(\varepsilon_{\text{sa}}\), both normal and tangential electric fields are non-zero. Here, the tangential field is continuous across the boundary, with \(\mathbf{E}_{\text{sa},\parallel} = \mathbf{E}_{\text{air},\parallel} = \mathbf{E}_{\text{substrate},\parallel}\). The normal field is discontinuous such that \(\mathbf{E}_{\text{sa}} \approx \mathbf{E}_{\text{air},\perp}/\varepsilon_{\text{sa}}\). We thus evaluate the field on the air side and, therefore, for a contaminant layer of thickness \(t_{\text{sa}}\), the stored energy becomes

  (5.71)\[ U_{\text{sa}} = \frac{1}{4} \varepsilon_0 t_{\text{sa}} \int_{S_{\text{sa}}} \left( \varepsilon_{\text{sa}} |\mathbf{E}_{\text{air},\parallel}|^2 + \frac{1}{\varepsilon_{\text{sa}}} |\mathbf{E}_{\text{air},\perp}|^2 \right) \mathrm{d}S.\]

To account for surface dielectric losses in QTCAD^®, users should instantiate SurfaceDielectricSpec for the different dielectric surfaces and use these objects when setting up EPRAnalysis.

[1]

The factor \(1/4\), rather than the \(1/2\) used in the EPR paper, comes from using averaged energy together with a peak sinusoidal current. For \(I(t)=I_\mathrm{peak}\cos(\omega t)\), the averaged inductor energy is \(\langle L I^2(t)/2\rangle_t = L I_\mathrm{peak}^2/4\). With an RMS current, the same averaged energy is \(L I_\mathrm{rms}^2/2\).