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.

[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\).