Effective Hamiltonians for interacting superconducting qubits: local basis reduction and the Schrieffer–Wolff transformation

The usual approach to identifying the computational basis states for theory and simulations, which is extensively used in the literature (cf for example [13, 15, 17, 31]), is based on a series expansion of the circuit Hamiltonian around a fixed value of one of its bias parameters (voltage or magnetic flux bias). For clarity, let us consider the specific case of the rf-SQUID qubit, whose circuit is shown in figure 1. Following the method introduced in section 2, we can write its circuit Hamiltonian (up to an additive constant) as [14]:

$$\begin{equation}{\hat{H}}_{e.,m.}\left({f}_{z}\right)=\frac{{\hat{Q}}^{2}}{2{C}_{J}}+\frac{{\hat{{\Phi}}}^{2}}{2L}-{E}_{J}\mathrm{cos}\left[2\pi \left(\frac{\hat{{\Phi}}}{{{\Phi}}_{0}}+{f}_{z}\right)\right],\end{equation} \tag{ 19 }$$

where ${f}_{z}={{\Phi}}_{z}/{{\Phi}}_{0}{:=}{{\Phi}}_{01}^{\mathrm{e}\mathrm{x}\mathrm{t}}/{{\Phi}}_{0}$ is the magnetic flux applied externally to the rf-SQUID loop, in units of ${{\Phi}}_{0}$. When ${f}_{z}\simeq 0.5$, we can rewrite the previous equation as:

$$\begin{equation}{\hat{H}}_{e.,m.}\left({f}_{z}\right)\simeq {\hat{H}}_{0}+\delta \hat{H}{:=}={\hat{H}}_{e.,m.}\left(0.5\right)+\delta {f}_{z}\frac{\partial {\hat{H}}_{e.,m.}\left({f}_{z}\right)}{\partial {f}_{z}}{\left\vert \right.}_{{f}_{z}=0.5},\end{equation} \tag{ 20 }$$

where $\delta {f}_{z}={f}_{z}-0.5$ and

$$\begin{equation}\frac{\partial {\hat{H}}_{e.,m.}\left({f}_{z}\right)}{\partial {f}_{z}}={E}_{J}\cdot \mathrm{sin}\left[\frac{2\pi }{{{\Phi}}_{o}}\left(\hat{{\Phi}}-{{\Phi}}_{z}\right)\right]\equiv -{{\Phi}}_{0}\frac{\hat{{\Phi}}}{L}{:=}-{{\Phi}}_{0}\hat{I},\end{equation} \tag{ 21 }$$

with $\hat{I}$ the loop current operator, which will define our computational basis. Notice that we used Kirchhoff's current law to go from the second to the third term in the last chain of equations [1].

At this point, we can invoke stationary perturbation theory to write the $n$th eigenstate of ${\hat{H}}_{e.,m.}\left({f}_{z}\right)$, up to first order in $\delta {f}_{z}$ as [32]:

$$\begin{equation}\vert {E}_{n}\rangle =\vert {E}_{n}^{\left(0\right)}\rangle +\sum _{m\ne n}\frac{\langle {E}_{m}^{\left(0\right)}\vert \delta \hat{H}\vert {E}_{n}^{\left(0\right)}\rangle }{{E}_{n}^{\left(0\right)}-{E}_{m}^{\left(0\right)}}\vert {E}_{m}^{\left(0\right)}\rangle ,\end{equation} \tag{ 22 }$$

where $\vert {E}_{n}^{\left(0\right)}\rangle :{\hat{H}}_{0}\vert {E}_{n}^{\left(0\right)}\rangle ={E}_{n}^{\left(0\right)}\vert {E}_{n}^{\left(0\right)}\rangle$ is the $n$th eigenstate of ${\hat{H}}_{0}$. As we can see, the various terms of the first order correction $\vert {E}_{n}\rangle -\vert {E}_{n}^{\left(0\right)}\rangle$ scale with the inverse of the differences between the unperturbed energies. Then, since for the rf-SQUID the spectrum of ${\hat{H}}_{0}={\hat{H}}_{e.,m.}\left({f}_{z}=0.5\right)$ is largely anharmonic, i.e. ${E}_{2}^{\left(0\right)}\gg {E}_{1}^{\left(0\right)}$, the two lowest-energy perturbed eigenstates are approximately linear combinations of their unperturbed counterparts only.

Since we are only interested in the two lowest eigenstates of the system (the qubit subspace), we can now project equation (20) on $\vert {E}_{0}^{\left(0\right)}\rangle$ and $\vert {E}_{1}^{\left(0\right)}\rangle$ and use the fact that $\langle {E}_{0}^{\left(0\right)}\vert \hat{I}\vert {E}_{0}^{\left(0\right)}\rangle =\langle {E}_{1}^{\left(0\right)}\vert \hat{I}\vert {E}_{1}^{\left(0\right)}\rangle =0$ (due to the symmetry of the Hamiltonian under magnetic field inversion about the point $f={f}_{z}=0.5$) to get:

$$\begin{equation}{\mathbf{H}}_{q}\left({f}_{z}\right)\simeq \frac{{E}_{0}^{\left(0\right)}+{E}_{1}^{\left(0\right)}}{2}{\boldsymbol{\sigma }}_{I}+\delta {f}_{z}{{\Phi}}_{0}{I}_{p}{\boldsymbol{\sigma }}_{x}+-\frac{{E}_{0}^{\left(0\right)}-{E}_{1}^{\left(0\right)}}{2}{\boldsymbol{\sigma }}_{z},\end{equation} \tag{ 23 }$$

where ${\boldsymbol{\sigma }}_{x}$ and ${\boldsymbol{\sigma }}_{z}$ are two of the standard Pauli matrices, ${\boldsymbol{\sigma }}_{I}$ is the $2{\times}2$ identity matrix and ${I}_{p}{:=}\langle {E}_{0}^{\left(0\right)}\vert \hat{I}\vert {E}_{1}^{\left(0\right)}\rangle =\langle {E}_{1}^{\left(0\right)}\vert \hat{I}\vert {E}_{0}^{\left(0\right)}\rangle { >}0$ (notice that we can always ensure these two conditions by multiplying $\vert {E}_{0}^{\left(0\right)}\rangle$ and $\vert {E}_{1}^{\left(0\right)}\rangle$ by appropriate phase factors). At this point we can diagonalise the current operator part of the Hamiltonian simply by introducing the two following computational states:

definition Definition 3.2 (Computational basis states (perturbative)).

$$\begin{equation}\begin{aligned}\hfill \vert 0\rangle & =\frac{\vert {E}_{0}^{\left(0\right)}\rangle +\vert {E}_{1}^{\left(0\right)}\rangle }{\sqrt{2}},\hfill \\ \hfill \vert 1\rangle & =\frac{\vert {E}_{0}^{\left(0\right)}\rangle -\vert {E}_{1}^{\left(0\right)}\rangle }{\sqrt{2}}\hfill \end{aligned}\end{equation} \tag{ 24 }$$

Using the results above, it is trivial to show that these are actually eigenstates of $\hat{I}$ with opposite eigenvalues: $\langle 0\vert \hat{I}\vert 0\rangle =-\langle 1\vert \hat{I}\vert 1\rangle ={I}_{p}$ [31]. In this basis, the effective Hamiltonian reads

$$\begin{equation}{\mathbf{H}}_{q}\left({f}_{z}\right)\simeq \frac{{E}_{0}^{\left(0\right)}+{E}_{1}^{\left(0\right)}}{2}{\boldsymbol{\sigma }}_{I}-\frac{{E}_{0}^{\left(0\right)}-{E}_{1}^{\left(0\right)}}{2}{\boldsymbol{\sigma }}_{x}+\delta {f}_{z}{{\Phi}}_{0}{I}_{p}{\boldsymbol{\sigma }}_{z}.\end{equation} \tag{ 25 }$$

Notice that every single-qubit Hamiltonian can be written in the general form

$$\begin{equation}{\mathbf{H}}_{q}=\sum _{i=I,x,y,z}{h}_{i}{\boldsymbol{\sigma }}_{i},\end{equation} \tag{ 26 }$$

where ${\boldsymbol{\sigma }}_{I}\equiv {\mathbb{I}}_{2}$, ${\boldsymbol{\sigma }}_{i=x,y,z}$ are the three standard Pauli matrices and the ${h}_{i}$'s are real coefficients. In the following we will call these Pauli coefficients, with specific reference to their values in the computational basis.

Although equation (25) already contains the analytic expressions of the Pauli coefficients (which apply to the rf-SQUID qubit), it is useful to consider the following equivalent derivation, which has a straightforward extension to the interacting qubit case. Once we have found the computational states according to (24), we can use the homomorphism between ${\mathbb{C}}^{2}$ and qubit subspace (${\mathcal{H}}_{q}=\mathrm{S}\mathrm{p}\mathrm{a}\mathrm{n}\left\{\vert 0\rangle ,\vert 1\rangle \right\}$) to introduce the following four operators, which represent the action of the Pauli matrices on ${\mathcal{H}}_{q}$:

$$\begin{equation}\begin{gathered}{c}{\hat{\sigma }}_{I}=\vert 0\rangle \langle 0\vert +\vert 1\rangle \langle 1\vert ,\quad {\hat{\sigma }}_{x}=\vert 0\rangle \langle 1\vert +\vert 1\rangle \langle 0\vert ,\\ {\hat{\sigma }}_{y}=-i\vert 0\rangle \langle 1\vert +i\vert 1\rangle \langle 0\vert ,\quad {\hat{\sigma }}_{z}=\vert 0\rangle \langle 0\vert -\vert 1\rangle \langle 1\vert .\end{gathered}\end{equation} \tag{ 27 }$$

Then, using the following property of the Pauli matrices,

$$\begin{equation}\mathrm{Tr}\left({\boldsymbol{\sigma }}_{i}\cdot {\boldsymbol{\sigma }}_{j}\right)=2{\delta }_{ij},\end{equation} \tag{ 28 }$$

we find that:

$$\begin{equation}{h}_{i}=\frac{1}{2}\mathrm{Tr}\left({\mathbf{H}}_{q}\cdot {\boldsymbol{\sigma }}_{i}\right)=\frac{1}{2}\mathrm{Tr}\left({\hat{H}}_{q}\cdot {\hat{\sigma }}_{i}\right)\equiv \frac{1}{2}\mathrm{Tr}\left({\hat{H}}_{e.,m.}\cdot {\hat{\sigma }}_{i}\right).\end{equation} \tag{ 29 }$$

Notice that here, both ${\hat{H}}_{e.,m.}$ and ${\hat{\sigma }}_{I,x,y,z}$ are conveniently expressed in whatever basis we initially choose for ${\hat{H}}_{e.,m.}$.

The perturbative reduction approach has a clear disadvantage: the effective Hamiltonian (25) reproduces the two lowest energy levels of the full circuit Hamiltonian only in the limit in which the first order perturbative expansion (22) holds. This entails two requirements. Firstly that the spectrum of the unperturbed Hamiltonian (in other words, the circuit Hamiltonian at the point of the expansion) is highly anharmonic, which is only true for some SC qubit designs and not for others (such as the capacitively-shunted flux qubit and the transmon) [12, 13]. Secondly, the perturbation to the bias parameter must be small, for instance $\vert \delta {f}_{z}\vert \ll 1$ for the rf-SQUID qubit [17].

A second common approach to the numerical calculation of the Pauli coefficients is the use of semi-classical theory. In this case the quantum state of the system is approximated by one that minimises its semi-classical potential, which is the part of the classical Hamiltonian depending on the coordinate variable (i.e. the flux in a flux qubit and the charge in a charge qubit). At the operational point the semi-classical potential of qubit circuits assumes a general double-well shape (or, more generally, that of a system of wells in more than one dimension), with two local minima very close in energy, such that quantum tunnelling can occur between them.

In this picture, the longitudinal Pauli coefficient ${h}_{z}$ is identified with the difference in energy between the two potential minima, whereas the effective transverse field ${h}_{x}$ corresponds to the tunnelling energy. This is calculated using the semi-classical instanton method (or equivalently the WKB approximation) [33]. These calculations are only accurate in the limit in which the tunnelling action across the potential barrier is very large, which implies that the tunnelling energy has to be exponentially small [34]. The instanton calculation of the transverse field for the rf-SQUID qubit is described in detail in appendix A.4.

In order to overcome the difficulties of the standard reduction approaches outlined above, we propose an alternative reduction method which relies on a local definition of the computational basis, i.e. one that explicitly depends on all of the circuit bias parameters. In other words, in this case the computational basis states are built as a linear combination of the two local circuit low-energy states:

$$\begin{equation}\begin{aligned}\hfill \vert 0\rangle ={u}_{00}\vert {E}_{0}\rangle +{u}_{01}\vert {E}_{1}\rangle ,\\ \hfill \vert 1\rangle ={u}_{10}\vert {E}_{0}\rangle +{u}_{11}\vert {E}_{1}\rangle ,\end{aligned}\end{equation} \tag{ 30 }$$

where ${\hat{H}}_{e.,m.}\vert {E}_{i}\rangle ={E}_{i}\vert {E}_{i}\rangle$ and ${\hat{H}}_{e.,m.}$ is the local circuit Hamiltonian. In order for these two states to be appropriately orthonormal, the ${u}_{ij}$'s have to be the elements of a unitary matrix,

$$\begin{equation}\mathbf{U}=\left(\begin{pmatrix}{cc}\hfill {u}_{00}\hfill & \hfill {u}_{10}\hfill \\ \hfill {u}_{01}\hfill & \hfill {u}_{11}\hfill \end{pmatrix}\right),\end{equation} \tag{ 31 }$$

which we will have to find. The unitarity condition ensures that when we transform from the energy eigenbasis $\left\{\vert {E}_{0}\rangle ,\vert {E}_{1}\rangle \right\}$ to the local computational basis $\left\{\vert 0\rangle ,\vert 1\rangle \right\}$ the spectrum of the effective qubit Hamiltonian (18) is unchanged and the two lowest-energy levels of the circuit Hamiltonian are preserved.

Owing to the orthonormality of $\mathbf{U}$ columns, we can always rewrite $\mathbf{U}$, up to an irrelevant global phase multiplication factor, as:

$$\begin{equation}\mathbf{U}=\left(\begin{pmatrix}{cc}{\mathrm{e}}^{i{\varphi }_{1}}\mathrm{cos}\theta \hfill & \hfill {\mathrm{e}}^{i{\varphi }_{2}}\mathrm{sin}\theta \hfill \\ \hfill -{\mathrm{e}}^{-i{\varphi }_{2}}\mathrm{sin}\theta \hfill & \hfill {\mathrm{e}}^{-i{\varphi }_{1}}\mathrm{cos}\theta \hfill \end{pmatrix}\right),\end{equation} \tag{ 32 }$$

where

$$\begin{align}\hfill \theta & =\mathrm{a}\mathrm{c}\mathrm{o}\mathrm{s}\vert {u}_{00}\vert =\mathrm{a}\mathrm{c}\mathrm{o}\mathrm{s}\vert {u}_{11}\vert ,\hfill \end{align} \tag{ 33 }$$

$$\begin{align}\hfill {\varphi }_{1}& =\frac{1}{2}\mathrm{a}\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{{u}_{00}{u}_{11}^{{\ast}}+{u}_{00}^{{\ast}}{u}_{11}}{2\;co{s}^{2}\theta }\right),\hfill \end{align} \tag{ 34 }$$

$$\begin{align}\hfill {\varphi }_{2}& =\frac{1}{2}\left[\pi -\mathrm{a}\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{{u}_{01}{u}_{10}^{{\ast}}+{u}_{01}^{{\ast}}{u}_{10}}{2\;si{n}^{2}\theta }\right)\right].\hfill \end{align} \tag{ 35 }$$

so that $\theta \in \left[0,\pi \right]$ and ${\varphi }_{1},{\varphi }_{2}\in \left[0,\pi /2\right]$.

Now we consider again the operational definition of the computational states. This specifies that these should be eigenstates of a certain observable $\hat{O}$. For a flux qubit $\hat{O}=\hat{I}$, the current operator associated with the qubit SC loop, whereas for a charge qubit $\hat{O}=\hat{Q}$ represents the charge on the qubit SC island. One can easily see that imposing this condition on the states (30) is equivalent to finding the two eigenstates of the operator

$$\begin{equation}{\hat{O}}_{p}={\hat{P}}_{0}\hat{O}{\hat{P}}_{0}=\langle {E}_{0}\vert \hat{O}\vert {E}_{0}\rangle \vert {E}_{0}\rangle \langle {E}_{0}\vert +\langle {E}_{0}\vert \hat{O}\vert {E}_{1}\rangle \vert {E}_{0}\rangle \langle {E}_{1}\vert +\langle {E}_{1}\vert \hat{O}\vert {E}_{0}\rangle \vert {E}_{1}\rangle \langle {E}_{0}\vert +\langle {E}_{1}\vert \hat{O}\vert {E}_{1}\rangle \vert {E}_{1}\rangle \langle {E}_{1}\vert ,\end{equation} \tag{ 36 }$$

associated with a non-zero eigenvalue¹, that is

definition Definition 3.3 (Computational basis states (local)). $\vert 0\rangle$ and $\vert 1\rangle$ such that

$$\begin{equation}\begin{gathered}{c}{\hat{O}}_{p}\vert 0\rangle ={u}_{0}\vert 0\rangle ,\\ {\hat{O}}_{p}\vert 1\rangle ={u}_{1}\vert 1\rangle ,\end{gathered}\end{equation} \tag{ 37 }$$

with ${u}_{0}\ne {u}_{1}$ and $\vert {u}_{0}\vert ,\vert {u}_{1}\vert { >}0$.

Notice that this definition coincides with the one used in the perturbative method at the specific bias point at which the Hamiltonian expansion is performed (for instance at ${f}_{z}=0.5$ in the the rf-SQUID qubit case).

By identifying $\vert {E}_{0}\rangle$ with the vector $\left(1,0\right)$ and $\vert {E}_{1}\rangle$ with $\left(0,1\right)$, we can rewrite ${\hat{O}}_{p}$ in the matrix form

$$\begin{equation}{\mathbf{O}}_{p}=\left(\begin{pmatrix}{cc}\hfill \langle {E}_{0}\vert \hat{O}\vert {E}_{0}\rangle \hfill & \hfill \langle {E}_{0}\vert \hat{O}\vert {E}_{1}\rangle \hfill \\ \hfill \langle {E}_{1}\vert \hat{O}\vert {E}_{0}\rangle \hfill & \hfill \langle {E}_{1}\vert \hat{O}\vert {E}_{1}\rangle \hfill \end{pmatrix}\right).\end{equation} \tag{ 38 }$$

Finding the eigenvalues and the eigenvectors of this $2{\times}2$ matrix is straightforward. In particular, for the eigenvalues, we have:

$$\begin{equation}{u}_{0,1}=\frac{t{\pm}\sqrt{{t}^{2}-4d}}{2},\end{equation} \tag{ 39 }$$

where $t=\mathrm{Tr}\left({\mathbf{O}}_{p}\right)$ and $d=\mathrm{det}\left({\mathbf{O}}_{p}\right)$. In accordance with the operational definitions given above, we need to enforce one condition on these eigenvalues. For a flux-qubit type circuit, we need to have ${u}_{1}{< }0{< }{u}_{0}$, which implies $\mathrm{det}\left({\mathbf{I}}_{p}\right){< }0$, or, more explicitly

$$\begin{equation}\langle {E}_{0}\vert \hat{I}\vert {E}_{0}\rangle \langle {E}_{1}\vert \hat{I}\vert {E}_{1}\rangle {< }\vert \langle {E}_{0}\vert \hat{I}\vert {E}_{1}\rangle {\vert }^{2}.\end{equation} \tag{ 40 }$$

For a qubit of the charge type, instead, we will require ${u}_{1}={u}_{0}{\pm}2e$ (up to some suitably small numerical error). If this condition is not satisfied, then the circuit cannot be operated as a qubit with the desired computational states and the reduction protocol fails. Note, however, that since we are not making use of a perturbative expansion or a semi-classical approximation here, the range of applicability of this local reduction method should be wider than that of the standard methods presented before.

If we now write the eigenvectors of ${\mathbf{O}}_{p}$ as ${\to {u}}_{0}=\left({u}_{00},{u}_{01}\right)$ and ${\to {u}}_{1}=\left({u}_{10},{u}_{11}\right)$, then equation (30) returns our desired computational basis states, which make ${\hat{O}}_{p}$ diagonal. Armed with ${\to {u}}_{0}$ and ${\to {u}}_{1}$, we can easily calculate the general expression of the effective qubit Hamiltonian in the computational basis. If we keep working with $2{\times}2$ matrices, the qubit Hamiltonian in the energy eigenbasis (18) takes the obvious diagonal form

$$\begin{equation}{\mathbf{H}}_{q}^{\prime }=\left(\begin{pmatrix}{cc}\hfill {E}_{0}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {E}_{1}\hfill \end{pmatrix}\right).\end{equation} \tag{ 41 }$$

Going from this basis to the computational basis amounts to applying the unitary transformation $\mathbf{U}$ defined above; this gives the effective qubit Hamiltonian in the computational basis as:

$$\begin{align}\hfill {\mathbf{H}}_{q}={\mathbf{U}}^{{\dagger}}{\mathbf{H}}_{q}^{\prime }\mathbf{U}& =\frac{{E}_{0}+{E}_{1}}{2}{\mathbb{I}}_{2}+\frac{{E}_{0}-{E}_{1}}{2}\left(\begin{pmatrix}{cc}\hfill \vert {u}_{00}{\vert }^{2}-\vert {u}_{01}{\vert }^{2}\hfill & \hfill {u}_{00}^{{\ast}}{u}_{10}-{u}_{01}^{{\ast}}{u}_{11}\hfill \\ \hfill {u}_{10}^{{\ast}}{u}_{00}-{u}_{11}^{{\ast}}{u}_{01}\hfill & \hfill \vert {u}_{10}{\vert }^{2}-\vert {u}_{11}{\vert }^{2}\hfill \end{pmatrix}\right)\hfill \\ \hfill & =\frac{{E}_{0}+{E}_{1}}{2}{\mathbb{I}}_{2}+\frac{{E}_{0}-{E}_{1}}{2}\left(\begin{pmatrix}{cc}\hfill \mathrm{cos}2\theta \hfill & \hfill {\mathrm{e}}^{-i\varphi }\mathrm{sin}2\theta \hfill \\ \hfill {\mathrm{e}}^{i\varphi }\mathrm{sin}2\theta \hfill & \hfill -\mathrm{cos}2\theta \hfill \end{pmatrix}\right)\hfill \\ \hfill & =\frac{{E}_{0}+{E}_{1}}{2}{\mathbb{I}}_{2}+\frac{{E}_{0}-{E}_{1}}{2}\left[\mathrm{sin}2\theta \cdot \left(\mathrm{cos}\varphi {\boldsymbol{\sigma }}_{x}+\mathrm{sin}\varphi {\boldsymbol{\sigma }}_{y}\right)+\mathrm{cos}2\theta {\boldsymbol{\sigma }}_{z}\right],\hfill \end{align} \tag{ 42 }$$

where ${\mathbb{I}}_{2}$ is the $2{\times}2$ identity matrix and $\varphi ={\varphi }_{1}-{\varphi }_{2}\in \left[-\pi /2,\pi /2\right]$.

We observe that, by rescaling the computational states ${\to {u}}_{0}$ and ${\to {u}}_{1}$ by two phase factors, say ${\mathrm{e}}^{i{\phi }_{0}}$ and ${\mathrm{e}}^{i{\phi }_{1}}$, i.e. by applying some local gauge transformation in the qubit subspace, we can always remove the imaginary component ${h}_{y}{\boldsymbol{\sigma }}_{y}$ of ${\mathbf{H}}_{q}$. In fact such a gauge transformation $\mathbf{G}\left({\phi }_{0},{\phi }_{1}\right)$ corresponds to a spin rotation around the z axis, multiplied by a global phase:

$$\begin{equation}\mathbf{G}\left({\phi }_{0},{\phi }_{1}\right)=\left(\begin{pmatrix}{cc}\hfill {\mathrm{e}}^{i{\phi }_{0}}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\mathrm{e}}^{i{\phi }_{1}}\hfill \end{pmatrix}\right)={\mathrm{e}}^{i\frac{{\phi }_{0}+{\phi }_{1}}{2}}\left(\begin{pmatrix}{cc}\hfill {\mathrm{e}}^{i\frac{{\phi }_{0}-{\phi }_{1}}{2}}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\mathrm{e}}^{-i\frac{{\phi }_{0}-{\phi }_{1}}{2}}\hfill \end{pmatrix}\right)={\mathrm{e}}^{i\frac{{\phi }_{0}+{\phi }_{1}}{2}}\cdot {\mathrm{e}}^{i\frac{{\phi }_{0}-{\phi }_{1}}{2}{\boldsymbol{\sigma }}_{z}}.\end{equation} \tag{ 43 }$$

Hence $\mathbf{G}\left(-{\varphi }_{1},-{\varphi }_{2}\right)$, which represents a rotation around z by the angle ${\varphi }_{2}-{\varphi }_{1}=-\varphi$ (followed by a rescaling by ${\mathrm{e}}^{-i\left({\varphi }_{1}+{\varphi }_{2}\right)/2}$), transforms $\mathrm{cos}\varphi {\boldsymbol{\sigma }}_{x}+\mathrm{sin}\varphi {\boldsymbol{\sigma }}_{y}$ into ${\boldsymbol{\sigma }}_{x}$, and makes the effective qubit Hamiltonian real, that is:

$$\begin{equation}{\mathbf{H}}_{q}=\frac{{E}_{0}+{E}_{1}}{2}{\mathbb{I}}_{2}-\frac{{\Delta}}{2}{\boldsymbol{\sigma }}_{x}-\frac{\varepsilon }{2}{\boldsymbol{\sigma }}_{z},\end{equation} \tag{ 44 }$$

where ${\Delta}=\left({E}_{1}-{E}_{0}\right)\mathrm{sin}2\theta$ and $\varepsilon =\left({E}_{1}-{E}_{0}\right)\mathrm{cos}2\theta$. Notice that this gauge transformation can equivalently be written as:

$$\begin{equation}\begin{aligned}\hfill {\to {u}}_{0}& \to \frac{\vert {u}_{00}\vert }{{u}_{00}}\cdot {\to {u}}_{0},\hfill \\ \hfill {\to {u}}_{1}& \to \frac{\vert {u}_{10}\vert }{{u}_{10}}\cdot {\to {u}}_{1}.\hfill \end{aligned}\end{equation} \tag{ 45 }$$

Expression (44) for the effective qubit Hamiltonian, is the one adopted by most of the literature on SC qubits [17, 25, 35]. (Note that the coefficient ${\Delta}$ is usually further assumed to be positive, a condition which can also always be achieved with a $\pi$ rotation about z.)

An equivalent and more convenient way of calculating the four Pauli coefficients ${h}_{i}$, $i=I,x,y,z$ than using equations (33), (42) and (26) together is again to use the computational states to build the Pauli operators and then to apply equation (29).

In section 4 we will present numerical simulations which benchmark the performance of the local reduction method against the standard methods and demonstrate the increased accuracy of the former relative to the latter ones.

In this subsection we briefly review the reduction method outlined in a recent work by Ozfidan et al [18]. This method starts by writing the low-energy part of the total circuit Hamiltonian ${\hat{H}}_{e.,m.}$, i.e. the component associated with its lowest ${2}^{N}$ eigenvalues, in its diagonal form: ${\mathbf{H}}_{q}^{\prime }=\mathrm{diag}\left({E}_{0},\dots ,\;{E}_{{2}^{N}-1}\right)$. Then a sequence of two rotations, say ${\mathbf{R}}_{1},{\mathbf{R}}_{2}\in \mathrm{S}\mathrm{O}\left({2}^{N}\right)$, is applied to it, producing the effective qubit Hamiltonian ${\mathbf{H}}_{q}={\mathbf{R}}_{2}^{T}{\mathbf{R}}_{1}^{T}{\mathbf{H}}_{q}^{\prime }{\mathbf{R}}_{1}{\mathbf{R}}_{2}$. Since orthogonal operations do not change the spectrum of an operator, this protocol guarantees by construction that the spectrum of ${\mathbf{H}}_{q}$ matches the low-energy spectrum of the circuit Hamiltonian.

The first rotation applied in this protocol, ${\mathbf{R}}_{1}$, maps from the low-energy eigenbasis of the total Hamiltonian ${\hat{H}}_{e.,m.}$, $\left\{\vert {E}_{0}\rangle ,\dots ,\;\vert {E}_{{2}^{N}-1}\rangle \right\}$ to that of the unperturbed Hamiltonian ${\hat{H}}_{0}$, i.e. $\left\{\vert {E}_{0}^{\left(0\right)}\rangle ,\dots ,\;\vert {E}_{{2}^{N}-1}^{\left(0\right)}\rangle \right\}$, and is initially calculated as

$$\begin{equation}{\left({\mathbf{R}}_{1}\right)}_{ij}=\langle {E}_{i}\vert {E}_{j}^{\left(0\right)}\rangle .\end{equation} \tag{ 51 }$$

However, as we pointed out before, $\vert {E}_{i}\rangle$ also has components outside of the subspace $\mathrm{S}\mathrm{p}\mathrm{a}\mathrm{n}\left\{{E}_{0}^{\left(0\right)}\rangle ,\dots ,\;\vert {E}_{{2}^{N}-1}^{\left(0\right)}\rangle \right\}$, which implies that this matrix is not orthogonal. ${\mathbf{R}}_{1}$ must therefore be explicitly orthonormalised, for instance using the Gram–Schmidt procedure. This step is only justified if the columns of ${\mathbf{R}}_{1}$ are already approximately orthonormal [18]. Since in our case orthogonality follows from normalisation, it suffices to check that

$$\begin{equation}\sum _{i=0}^{{2}^{N}-1}{\left({\mathbf{R}}_{1}\right)}_{ij}^{2}\approx 1,\;\forall j=0,\dots ,\;{2}^{N}-1,\end{equation} \tag{ 52 }$$

before we apply the Gram–Schmidt procedure.

To obtain the qubit Hamiltonian we now need the second rotation ${\mathbf{R}}_{2}$ to map from the basis of the energy eigenstates $\vert {E}_{0}^{\left(0\right)}\rangle ,\dots \;$ to the computational basis. We then take

$$\begin{equation}{\left({\mathbf{R}}_{2}\right)}_{ij}=\langle i\vert {E}_{j}^{\left(0\right)}\rangle ,\end{equation} \tag{ 53 }$$

where $\vert i\rangle =\vert {i}_{{2}^{N}-1}\rangle \otimes \cdots \otimes \vert {i}_{0}\rangle$ is an outer product of single qubit computational states, with ${i}_{{2}^{N}-1}{i}_{{2}^{N}-2}\cdots {i}_{1}{i}_{0}$ the $N$-digit binary representation of the integer $i\in \left\{0,\;1,\dots ,\;{2}^{N}-1\right\}$. These computational states are found from the reduction of the unperturbed single-qubit Hamiltonians. If the local reduction method is used for this, the rotation matrix ${\mathbf{R}}_{2}$ is guaranteed to be orthogonal.

Note that, although the effective Hamiltonian calculated with this method has the correct spectrum, the procedure is based on the approximate equality (52), which is not often satisfied, particularly in the case of relatively large interactions. (This can be seen by considering, once again, the perturbative expansion (22).) Additionally, the previous derivation implicitly assumes that the circuit Hamiltonian is real, so that all the eigenstates and computational states can be chosen to have only real components. This ensures that ${\mathbf{R}}_{1},{\mathbf{R}}_{2}\in \mathrm{S}\mathrm{O}\left({2}^{N}\right)$. Some circuits, however, may have an efficient matrix representation of the Hamiltonian which is complex. In this case the definition of the two rotations would lead to the presence of arbitrary complex phases in their elements, which would need to be somehow taken care of. (Notice that even in the real case the scalar products defining the elements of ${\mathbf{R}}_{1}$ and ${\mathbf{R}}_{2}$ are only defined up to an arbitrary sign.)

A second method of determining the effective Hamiltonian of a multi-qubit system is presented in a recent work by Melanson et al [37]. This method works under the more restrictive assumption that the effective Hamiltonian is diagonal in the computational basis. In this case the lowest ${2}^{N}$ eigenstates of the circuit are also eigenstates of the single-qubit operators ${\hat{O}}_{i}$ specifying the computational basis and the corresponding eigenvalues can be calculated numerically as the expectation values $\langle {E}_{n}\vert {\hat{O}}_{i}\vert {E}_{n}\rangle$. Additionally the ${2}^{N}$ non-zero Pauli coefficients of the system can be expressed as a linear combination of its low-energy eigenvalues [37]. For instance, in the two-qubit case one has:

$$\begin{equation}\left(\begin{pmatrix}{c}\hfill {E}_{00}\hfill \\ \hfill {E}_{01}\hfill \\ \hfill {E}_{10}\hfill \\ \hfill {E}_{11}\hfill \end{pmatrix}\right)=\left(\begin{pmatrix}{cccc}\hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill \\ \hfill 1\hfill & \hfill 1\hfill & \hfill -1\hfill & \hfill -1\hfill \\ \hfill 1\hfill & \hfill -1\hfill & \hfill 1\hfill & \hfill -1\hfill \\ \hfill 1\hfill & \hfill -1\hfill & \hfill -1\hfill & \hfill 1\hfill \end{pmatrix}\right)\cdot \left(\begin{pmatrix}{c}\hfill {h}_{II}\hfill \\ \hfill {h}_{zI}\hfill \\ \hfill {h}_{Iz}\hfill \\ \hfill {h}_{zz}\hfill \end{pmatrix}\right){:=}\mathbf{M}\cdot \left(\begin{pmatrix}{c}\hfill {h}_{II}\hfill \\ \hfill {h}_{zI}\hfill \\ \hfill {h}_{Iz}\hfill \\ \hfill {h}_{zz}\hfill \end{pmatrix}\right),\end{equation} \tag{ 54 }$$

where ${E}_{ij}$, $i,j\in \left\{0,1\right\}$ is the eigenvalue of the circuit Hamiltonian corresponding to the computational state $\vert i\rangle \vert j\rangle$. We can therefore determine the Pauli coefficients of the two-qubit system by finding the lowest four energy eigenvalues of its circuit Hamiltonian, calculating the expectation value of the operators ${\hat{O}}_{1,2}$ on the each eigenstate to identify its corresponding computational state and by inverting the previous equation to get

$$\begin{equation}\left(\begin{pmatrix}{c}\hfill {h}_{II}\hfill \\ \hfill {h}_{zI}\hfill \\ \hfill {h}_{Iz}\hfill \\ \hfill {h}_{zz}\hfill \end{pmatrix}\right)={\mathbf{M}}^{-1}\cdot \left(\begin{pmatrix}{c}\hfill {E}_{00}\hfill \\ \hfill {E}_{01}\hfill \\ \hfill {E}_{10}\hfill \\ \hfill {E}_{11}\hfill \end{pmatrix}\right).\end{equation} \tag{ 55 }$$

The same procedure can be applied to systems with three or more qubits (plus eventual additional couplers).

In practice, an effective Hamiltonian diagonal in the computational basis is verified when the qubit tunnelling barriers are high (negligible transverse field ${h}_{x}$) and the qubits are coupled only through their z degree of freedom (that is when the coupling is inductive between flux qubits or capacitive between charge qubits). A Hamiltonian of this form is however classical and cannot be sufficient for universal quantum computation [28]. This method can nevertheless still be useful when it is reasonable to assume that the different non-commuting terms of the qubit Hamiltonian can be turned on and off independently.

In this final subsection we introduce a new reduction protocol for multi-qubit systems which overcomes some of the limitations of the methods described above. In particular, this method does not require the mixing between the qubit subspace (equation (50)) and its complement, resulting from the interactions, to be negligible, which is a crucial assumption of the approximate rotation reduction. Secondly, unlike the approximate rotation reduction, it can be applied directly to circuit Hamiltonians with complex elements, since the arbitrary phase choices made when numerically evaluating the Hamiltonian eigenvectors cancel out in all the necessary expressions. Thirdly, the reduction method introduced here can be applied to find arbitrary non-diagonal effective Hamiltonians. This is all made possible by the Schrieffer–Wolff transformation (SWT), which by construction maps the total circuit Hamiltonian ${\hat{H}}_{e.,m.}$ to a new Hermitian operator acting on the qubit subspace ${\mathcal{H}}_{q}$ and whose spectrum matches the low-energy spectrum of ${\hat{H}}_{e.,m.}$, which is precisely what we expect from the effective qubit Hamiltonian [23].

The SWT relies on a single assumption regarding the form of the full system Hamiltonian, namely that the spectrum of the unperturbed part of the Hamiltonian (excluding the interactions) has a sufficiently large gap, as we will see below. For the purpose of this reduction method, we will replace this assumption with an equivalent pair of two distinct conditions. In order to state the first one, let us rewrite the unperturbed part of the N-qubit M-coupler system Hamiltonian (49) as

$$\begin{equation}{\hat{H}}_{0}={\hat{P}}_{0}{\hat{H}}_{0}{\hat{P}}_{0}+{\hat{Q}}_{0}{\hat{H}}_{0}{\hat{Q}}_{0},\end{equation} \tag{ 56 }$$

where

$$\begin{equation}{\hat{P}}_{0}=\sum _{i=0}^{{2}^{N}-1}\vert {E}_{i}^{\left(0\right)}\rangle \langle {E}_{i}^{\left(0\right)}\vert \end{equation} \tag{ 57 }$$

is the projector on the low-energy eigenspace ${\mathcal{H}}_{\mathrm{l}\mathrm{o}\mathrm{w}}^{\left(0\right)}$, spanned by the eigenstates corresponding to the lowest ${2}^{N}$ eigenvalues of ${\hat{H}}_{0}$, and ${\hat{Q}}_{0}=\hat{\mathbb{I}}-{\hat{P}}_{0}$ projects on the complementary subspace $\mathcal{H}{\backslash}{\mathcal{H}}_{\mathrm{l}\mathrm{o}\mathrm{w}}^{\left(0\right)}$. Notice that, given ${E}_{i}^{\left(0\right)}$, the ith eigenvalue of ${\hat{H}}_{0}$, the spectrum of ${\hat{P}}_{0}{\hat{H}}_{0}{\hat{P}}_{0}$ is by definition ${\mathcal{S}}_{\mathrm{l}\mathrm{o}\mathrm{w}}^{\left(0\right)}=\left\{{E}_{0}^{\left(0\right)},{E}_{1}^{\left(0\right)},\dots ,\;{E}_{{2}^{N}-1}^{\left(0\right)}\right\}$, whereas ${\hat{Q}}_{0}{\hat{H}}_{0}{\hat{Q}}_{0}$ has the set of eigenvalues ${\mathcal{S}}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}^{\left(0\right)}=\left\{{E}_{{2}^{N}}^{\left(0\right)},\dots \;\right\}$. The first assumption of our reduction is that ${\mathcal{H}}_{q}\equiv {\mathcal{H}}_{\mathrm{l}\mathrm{o}\mathrm{w}}^{\left(0\right)}$, i.e. that no additional excited state of the independent circuits is mixed in the low energy subspace of ${\hat{H}}_{0}$, and that the two sets ${\mathcal{S}}_{\mathrm{l}\mathrm{o}\mathrm{w}}^{\left(0\right)}$ and ${\mathcal{S}}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}^{\left(0\right)}$ are separated by at least ${\Delta}{ >}0$, that is $\vert {E}_{{2}^{N}}^{\left(0\right)}-{E}_{{2}^{N}-1}^{\left(0\right)}\vert {\geqslant}{\Delta}$. This composite condition can be written, more explicitly, in the following form:

$$\begin{equation}\begin{gathered}{c}\left\vert {\Delta},{E}_{i,2},-,\sum _{j=1}^{N},{\Delta},{E}_{j,1}\right\vert {\geqslant}{\Delta},\;\forall i=1,\dots ,\;N\\ \left\vert {\Delta},{E}_{ci,1},-,\sum _{j=1}^{N},{\Delta},{E}_{j,1}\right\vert {\geqslant}{\Delta},\;\forall i=1,\dots ,\;M,\end{gathered}\end{equation} \tag{ 58 }$$

where we have introduced the notation ${\Delta}{E}_{i,j}={E}_{i,j}-{E}_{i,0}$ and ${\Delta}{E}_{ci,j}={E}_{ci,j}-{E}_{ci,0}$ (with ${E}_{i,j}$ (${E}_{ci,j}$) again the jth eigenstate of the ith qubit (coupler) unperturbed Hamiltonian). Since the summations above grow linearly with the number of qubits in the system, this condition limits the size of the systems to which we can apply our reduction method. Intuitively this limit reflects the impossibility of finding any coherent description of the low-energy spectrum of a composite system in terms of interacting two-level subsystems, whenever the second excited state of one of these subsystems appears in the spectrum. Therefore if we are interested in characterising a very large circuit, we should first subdivide it into smaller connected subsystems for which the inequalities (58) hold.

The second requirement is simply that the strength of the interaction Hamiltonian should be small compared to the spectral gap of ${\hat{H}}_{0}$, ${\Delta}$. Namely:

$$\begin{equation}{{\Vert}{\hat{H}}_{\mathrm{int}}{\Vert}}_{op}{< }\frac{{\Delta}}{2},\end{equation} \tag{ 59 }$$

where ${\Vert}\cdot {{\Vert}}_{op}$ is the operator norm:

$$\begin{equation}{\Vert}\hat{O}{{\Vert}}_{op}=\mathrm{sup}\left\{{\Vert}\hat{O}\vert {\Psi}\rangle {\Vert}:{\Vert}\vert {\Psi}\rangle {\Vert}=1\right\},\end{equation} \tag{ 60 }$$

with ${\Vert}\cdot {\Vert}$ the two-norm $\sqrt{\langle \;\cdot \;\vert \;\cdot \;\rangle }$.

Since the addition of the interaction term ${\hat{H}}_{\mathrm{int}}$ can shift the eigenvalues of ${\hat{H}}_{0}$ by at most ${{\Vert}{\hat{H}}_{\mathrm{int}}{\Vert}}_{op}$, this second inequality implies that the spectrum of ${\hat{H}}_{e.,m.}$ remains gapped. This in turn allows us to rewrite the total Hamiltonian in the block-diagonal form ${\hat{H}}_{e.,m.}=\hat{P}{\hat{H}}_{e.,m.}\hat{P}+\hat{Q}{\hat{H}}_{e.,m.}\hat{Q}$, where $\hat{P}$ is the projector on the ${2}^{N}$-dimensional low-energy eigenspace of ${\hat{H}}_{e.,m.}$, ${\mathcal{H}}_{\mathrm{l}\mathrm{o}\mathrm{w}}$, and $\hat{Q}=\hat{\mathbb{I}}-\hat{P}$ [23].

Additionally, according to [23], since ${\mathcal{H}}_{\mathrm{l}\mathrm{o}\mathrm{w}}$ and ${\mathcal{H}}_{q}$ have the same dimension, they are connected by a direct rotation $\hat{U}$ such that

$$\begin{equation}\begin{gathered}{c}\hat{U}\hat{P}{\hat{U}}^{{\dagger}}={\hat{P}}_{0},\\ \hat{U}\hat{Q}{\hat{U}}^{{\dagger}}={\hat{Q}}_{0}.\end{gathered}\end{equation} \tag{ 61 }$$

$\hat{U}$ is called the Schrieffer–Wolff transformation and can be written, in terms of the projectors, as [23]:

$$\begin{equation}\hat{U}=\sqrt{\left(2{\hat{P}}_{0}-\hat{\mathbb{I}}\right)\left(2\hat{P}-\hat{\mathbb{I}}\right)}.\end{equation} \tag{ 62 }$$

The principal square root $\sqrt{\cdot }$ above is well-defined as long as

$$\begin{equation}{{\Vert}\hat{P}-{\hat{P}}_{0}{\Vert}}_{op}{< }1,\end{equation} \tag{ 63 }$$

which in our case can be shown to be equivalent to (59) [23].

Now the action of the SWT on ${\hat{H}}_{e.,m.}$ is given by

$$\begin{equation}\begin{aligned}\hfill \hat{U}{\hat{H}}_{e.,m.}{\hat{U}}^{{\dagger}}& =\hat{U}\hat{P}{\hat{H}}_{e.,m.}\hat{P}{\hat{U}}^{{\dagger}}+\hat{U}\hat{Q}{\hat{H}}_{e.,m.}\hat{Q}{\hat{U}}^{{\dagger}}\hfill \\ \hfill & ={\hat{P}}_{0}\hat{U}{\hat{H}}_{e.,m.}{\hat{U}}^{{\dagger}}{\hat{P}}_{0}+{\hat{Q}}_{0}\hat{U}{\hat{H}}_{e.,m.}{\hat{U}}^{{\dagger}}{\hat{Q}}_{0},\hfill \end{aligned}\end{equation} \tag{ 64 }$$

where we used the identities $\hat{U}\hat{P}={\hat{P}}_{0}\hat{U}$ and $\hat{U}\hat{Q}={\hat{Q}}_{0}\hat{U}$. According to equation (64), $\hat{U}{\hat{H}}_{e.,m.}{\hat{U}}^{{\dagger}}$ is block-diagonal with respect to ${\hat{P}}_{0}$ and ${\hat{Q}}_{0}$. This finally leads us to the conclusion that

$$\begin{equation}{\hat{H}}_{q}{:=}{\hat{P}}_{0}\hat{U}{\hat{H}}_{e.,m.}{\hat{U}}^{{\dagger}}{\hat{P}}_{0}\end{equation} \tag{ 65 }$$

is an Hermitian operator, acting on ${\mathcal{H}}_{q}$, whose ${2}^{N}$ non-zero eigenvalues are the same as the lowest eigenvalues of the original interacting Hamiltonian ${\hat{H}}_{e.,m.}$ (because the unitary $\hat{U}$ leaves the spectrum of $\hat{P}{\hat{H}}_{e.,m.}\hat{P}$ unchanged) [23]. ${\hat{H}}_{q}$ therefore represents our effective qubit Hamiltonian, from which we can directly extract the Pauli coefficients by rewriting equation (48) as

$$\begin{equation}{h}_{\to {\eta }}=\frac{1}{{2}^{N}}\mathrm{Tr}\left({\hat{H}}_{q}\cdot {\hat{\sigma }}_{\to {\eta }}\right).\end{equation} \tag{ 66 }$$

In this case, the Pauli operator ${\hat{\sigma }}_{\to {\eta }}={\hat{\sigma }}_{{\eta }_{1}}\otimes \cdots \otimes {\hat{\sigma }}_{{\eta }_{N}}\otimes {\hat{P}}_{c}$ is built from the single-qubit Pauli operators $\left\{{\hat{\sigma }}_{{\eta }_{i}}\right\}$, which, in turn, are obtained as in the single-qubit case, starting from the unperturbed Hamiltonian ${\hat{H}}_{i}$ of each qubit and the appropriate operator ${\hat{O}}_{p,i}$. The operator

$$\begin{equation}{\hat{P}}_{c}=\underset{i=1}{\overset{M}{\otimes }}\vert {E}_{ci,0}\rangle \langle {E}_{ci,0}\vert \end{equation} \tag{ 67 }$$

represents the required identities acting on each of the ground-state energy subspaces of the coupler circuits.

Finally note that since both the approximate rotation reduction and the SWT reduction method determine an effective qubit Hamiltonian with the correct spectrum, the two results must be equivalent up to a unitary transformation. However, as mentioned before, the SWT reduction extends the range of applicability of the method to Hamiltonians with complex elements and does not involve the restrictive assumption (52).

We start by considering the simplest example of a flux qubit, i.e. the rf-SQUID circuit. As shown in figure 1, this consists of a Josephson junction, with tunnelling energy ${E}_{J}$, shunted by a superconducting inductive loop with self-inductance L and in parallel with its intrinsic capacitance ${C}_{J}$ [4].

Figure 2 shows the lowest five energy eigenvalues of the circuit, calculated as a function of the dimensionless external magnetic flux ${f}_{z}={{\Phi}}_{z}/{{\Phi}}_{0}\equiv {{\Phi}}_{01}^{\mathrm{e}\mathrm{x}\mathrm{t}}/{{\Phi}}_{0}$. (Note that the constant offset ${E}_{0}\left({f}_{z}=0.49\right)$ has been subtracted from all the energies.) The parameters used for the simulations are ${E}_{J}=125$ GHz, ${C}_{J}=5$ fF and $L=2.5$ nH, which are typical for this type of device [14]. In this case, the Hamiltonian was represented in a basis of harmonic oscillator occupation number states, truncated at a maximum occupation number of 40, which ensured the convergence of the low energy spectrum (cf appendix A.3) [40].

Zoom In Zoom Out Reset image size

As we can see from the graph in figure 2, the lowest two energy levels of the system (i.e. the qubit states), vary approximately linearly with the flux ${f}_{z}$, except around the symmetry point ${f}_{z}=0.5$, where they show a characteristic avoided crossing. In fact, as we saw previously, for small values of $\vert \delta {f}_{z}\vert =\vert {f}_{z}-0.5\vert$ the rf-SQUID Hamiltonian is well approximated by its first order expansion in $\delta {f}_{z}$. This maps to an effective qubit Hamiltonian of the form (see equation (25))

$$\begin{equation}{\mathbf{H}}_{q}\left({f}_{z}\right)=\frac{{\Delta}}{2}{\boldsymbol{\sigma }}_{x}+\frac{\varepsilon \left({f}_{z}\right)}{2}{\boldsymbol{\sigma }}_{x},\end{equation} \tag{ 68 }$$

where we have neglected the term proportional to the identity, ${\Delta}={\left({E}_{1}\left({f}_{z}\right)-{E}_{0}\left({f}_{z}\right)\right)}_{{f}_{z}=0.5}$ and $\varepsilon \left({f}_{z}\right)=2{{\Phi}}_{0}{I}_{p}\delta {f}_{z}$. The lowest two energy levels of the circuit are therefore approximately ${E}_{0,1}=\text{const.}\mp \sqrt{{{\Delta}}^{2}/4\left({{\Delta}}^{2}/4\right)+{{\Phi}}_{0}^{2}{I}_{p}^{2}\vert \delta {f}_{z}{\vert }^{2}}$, which become linear in ${f}_{z}$ for larger values of $\vert \delta {f}_{z}\vert$.

Figure 3(a) shows the values of the system Pauli coefficients as a function of ${f}_{z}$, calculated using equation (29). The solid lines correspond to values obtained by defining the Pauli operators according to the local reduction (LR) method introduced here (subsection 3.1.3). These are compared with the result of the perturbative (PR, empty circles) and instanton (crosses) methods. As we can see, the three reduction methods produce largely compatible results for this circuit. In particular, away from the symmetry point the LR method finds a 10% increase in the transverse field ${h}_{x}$ at the boundary of the flux interval considered, compared to its centre. The result of PR is instead independent of ${f}_{z}$, in agreement with equation (25). The values of ${h}_{z}$ and ${h}_{I}$ calculated with the LR and PR methods are compatible to 1% over the whole flux bias range. This implies that the definition of the computational basis in the LR method coincides, as it should, with that of the standard PR method in the limit in which the series expansion (20) and perturbation theory apply. As for the semi-classical calculations, these appear to over-estimate both the longitudinal field (by $\simeq 40$%) and the transverse field (by up to 6%), compared to the other two reduction methods.

Zoom In Zoom Out Reset image size

Since the semi-classical approximation applies in the limit where $\hslash$ is much smaller than the actions at play in the system, i.e. $S\gg \hslash$, and since the tunnelling energy ${h}_{x}$ decreases exponentially with the tunnelling action, ${h}_{x}\propto {\mathrm{e}}^{-S/\hslash }$ (see appendix A.4), we expect the result of the instanton calculations to be more accurate in the limit where ${h}_{x}$ is small [33, 34]. To verify this, we determined the qubit transverse field in the case of biasing at the symmetry point ${f}_{z}=0.5$ for increasing values of the loop inductance $L$. As we can see in figure 3(b), increasing $L$ causes the barrier between the two semi-classical potential wells (blue line) to rise, therefore suppressing the tunnelling ${h}_{x}$ (data in red). Since ${f}_{z}=0.5$, the perturbative and local reduction methods coincide, and they both determine the correct value of the tunnelling energy: ${h}_{x}=-{\Delta}E/2$, where ${\Delta}E$ is the energy separation between the ground and first excited state of the circuit (cf dots and solid line in figure 2). As expected, the instanton method result (crosses in figure 2) closely approaches that of the Hamiltonian reduction only as $L$ increases and $\vert {h}_{x}\vert$ becomes smaller.

At this point, as a consistency check, we can calculate the spectrum of the reduced qubit Hamiltonian simply as ${E}_{0,1}={h}_{I}\mp \sqrt{{h}_{x}^{2}+{h}_{y}^{2}+{h}_{z}^{2}}$ and compare it with that obtained from the full circuit model. PR and LR do a good job in reproducing the low-energy spectrum of the rf-SQUID qubit, as we can see from the plot in figure 3(c). This also shows the spectrum derived from the semi-classical model (crosses), which does not agree with the correct circuit spectrum as well.

Notice that LR is guaranteed to exactly reproduce the circuit levels as long as ${f}_{z}\simeq 0.5$. As mentioned in the previous section, the LR protocol only fails when, as $\vert {f}_{z}-0.5\vert$ increases, the two eigenvalues of ${\hat{I}}_{p}\left({f}_{z}\right)={\hat{P}}_{0}\left({f}_{z}\right)\hat{I}{\hat{P}}_{0}\left({f}_{z}\right)$ begin to have the same sign, meaning that no measurement distinguishing two qubit states with opposite persistent current is possible at the given bias. For the particular rf-SQUID circuit considered here, the local reduction method breaks down for $\vert {f}_{z}-0.5\vert \gtrsim 0.035$, as shown in figure 3(d) (region shaded in red). As we can see in this plot, as we approach this region the behaviour of the Pauli coefficients starts changing. In particular the transverse field increases considerably in magnitude, while the longitudinal field saturates. The green dotted lines in figure 3(d) show the circuit energy levels. We see that at the boundary of the unshaded region the second excited state starts mixing with the first, leading to an avoided crossing. This mixing means that, at this point, the two-level approximation does not hold any more, which leads to the failure of the LR.

Finally we might want to consider how well the reduced Hamiltonians are able to reproduce the correct expectation values of some circuit operator $\hat{O}$, i.e. whether the following relationship holds

$$\begin{equation}\langle {E}_{i}\vert \hat{O}\vert {E}_{j}\rangle =\langle {\mathcal{E}}_{i}\vert {\hat{O}}_{p}\vert {\mathcal{E}}_{j}\rangle ,\;\forall i,j\in \left\{0,1\right\},\end{equation} \tag{ 69 }$$

where ${\left\{\vert {E}_{i}\rangle \right\}}_{i=0,1}$ are energy eigenstates of ${\hat{H}}_{e.,m.}$ and ${\left\{\vert {\mathcal{E}}_{i}\rangle \right\}}_{i=0,1}$ are the eigenstates of the corresponding effective qubit Hamiltonian. ${\hat{O}}_{p}$ is defined locally as ${\hat{P}}_{0}\left({f}_{z}\right)\hat{O}{\hat{P}}_{0}\left({f}_{z}\right)$ (where ${\hat{P}}_{0}\left({f}_{z}\right)$ is the projector on the two-dimensional low-energy subspace of ${\hat{H}}_{e.,m.}\left({f}_{z}\right)$) in the LR method case, and is defined globally as ${\hat{P}}_{0}\left(0.5\right)\hat{O}{\hat{P}}_{0}\left(0.5\right)$ in the PR case. Figure 4 shows the matrix elements of the loop current operator $\hat{I}$ between qubit states, calculated with both the full and the reduced operators. We observe that LR is ensured to give the exact result, while PR produces a reasonable result.

Zoom In Zoom Out Reset image size

We have shown here that the approximations inherent in the perturbative reduction method are valid and sufficient for for determining the reduced Hamiltonian in the case of the simple rf-SQUID qubit of figure 1. We will see in the next subsection, however, that this is not true in the general case and that the local reduction method has a wider range of validity.

The accuracy of the perturbative reduction method deteriorates when we consider other flux qubit designs, particularly those with reduced anharmonicity like the capacitively-shunted flux qubit shown in figure 5. This consists of a superconducting loop interrupted by three Josephson junctions. The area of one junction is a factor $\alpha {< }1$ smaller than that of the other two and is shunted by a relatively large capacitor ${C}_{\text{sh}}\gg {C}_{JT}$. The capacitive shunt reduces the qubit sensitivity to charge noise, while improving the device reproducibility (by compensating for the fabrication variability of the junction size, which affects ${C}_{JT}$). At the same time, the effect of flux noise is mitigated by choosing small values of $\alpha$ (typically $0.125{< }\alpha {< }0.5$), which reduce the magnitude of the persistent current and therefore the magnetic dipole moment of the circuit [13]. The result is superconducting qubits with typical measured relaxation times ${T}_{1}$ in excess of 40 ms (three orders of magnitude longer than the standard rf-SQUID ${T}_{1}$) and decoherence times approaching the relaxation limit ${T}_{2}=2{T}_{1}$ [13].

Zoom In Zoom Out Reset image size

This substantial coherence enhancement comes at the cost of a decrease in the spectrum anharmonicity. We can see this by looking at figure 6(a), which shows the calculated low energy spectrum of a C-shunt qubit circuit as a function of ${f}_{z}={{\Phi}}_{z}/{{\Phi}}_{0}={{\Phi}}_{23}^{\mathrm{e}\mathrm{x}\mathrm{t}}/{{\Phi}}_{0}$, and comparing it with figure 2. The physical parameters used for the simulation are shown in table 1 (cf figure 5 for the meaning of the symbols). For the two lower junctions we used ${E}_{JL}={E}_{JR}={E}_{JT}/\alpha$ and ${C}_{JL}={C}_{JR}={C}_{JT}/\alpha$. These parameters are compatible with those reported in the experiments in reference [13].In this case, the Hamiltonian was represented numerically by projecting on a finite basis consisting of harmonic oscillator states for the mode associated with the circuit node 1 and charge number states for the modes associated with nodes 2 and 3 (cf appendix A.3) [35, 40, 41].

Zoom In Zoom Out Reset image size

Table 1. Physical parameters of the capacitively-shunted flux qubit used for the simulations in figure 6.

Parameter	Value
${E}_{JT}$	45 GHz
${C}_{JT}$	1.8 fF
$\alpha$	0.43
${C}_{\text{sh}}$	50 fF
$L$	100 pH

As we can see from figure 6(a), the two dispersion relations ${E}_{0,1}\left({f}_{z}\right)$ have first derivatives with the same sign everywhere. Since

$$\begin{equation}\langle {\hat{I}\rangle }_{0,1}{:=}-\langle \frac{\partial {\hat{H}}_{e.,m.}}{\partial {{\Phi}}_{z}}\rangle \simeq \frac{\partial {E}_{0,1}}{\partial {{\Phi}}_{z}},\end{equation} \tag{ 70 }$$

this means that the average persistent currents in the two energy eigenstates have equal sign (cf figure 6(d)). This is in contrast with the rf-SQUID flux qubit [13], but does not preclude the possibility to find two current eigenstates with opposite sign in the qubit subspace.

Figure 6(b) shows the Pauli coefficients obtained by the perturbative (circles) and local (lines) reduction methods. As anticipated, there is a clear discrepancy between the two results. In fact, owing to the much smaller anharmonicity of this circuit compared to the rf-SQUID, the two low-energy eigenstates of the circuit Hamiltonian at ${f}_{z}=0.5$ are not a good approximation for those away from ${f}_{z}=0.5$. This implies that projecting ${\hat{H}}_{e.,m.}\left({f}_{z}\right)$ on the states (24) does not preserve its low-energy spectrum and does not lead to the correct reduction. From the numerical results we see that the slope of ${h}_{z}\left({f}_{z}\right)$ in the local reduction case is smaller than in the perturbative reduction and further decreases away from ${f}_{z}=0.5$. Additionally, the transverse field ${h}_{x}\left({f}_{z}\right)$ shows a clear negative curvature in the LR results, whereas it is roughly constant in ${f}_{z}$ in the PR case (as in the rf-SQUID). The strong dependence of the transverse field on ${f}_{z}$ is a known distinguishing feature of the C-shunt flux qubit design when compared to more standard flux qubit circuits like the rf-SQUID [9, 13].

Calculating the spectra of the two reduced Hamiltonians leads to the result shown in figure 6(c). The local reduction result (filled dots) again reproduces the circuit ground and first excited states (lines) exactly, while the perturbative reduction fails to accurately predict the first excited state. Finally figure 6(d) shows the matrix elements of the current operator between the qubit energy eigenstates, calculated using the full circuit model (lines) and the two reduced two-level models (circles). The PR (empty circles) gives incorrect expectation values, which are opposite in sign for the two states.

We begin this subsection on coupled SC qubit systems by considering a simple two-qubit system, without any non-linear coupling element. As one such example we consider the system which Ozfidan et al characterised experimentally in [18]. This is composed of two compound-Josephson-junction rf-SQUID qubits (where the single Josephson junction is replaced by two junctions in parallel, forming a dc-SQUID) coupled both inductively and capacitively, as shown in figure 7. Assuming that the dc-SQUID loop is very small (such that its inductance is much smaller than both the main loop inductance and the Josephson inductance ${\left({{\Phi}}_{0}/2\pi \right)}^{2}/{E}_{J}$), we can effectively describe it as a single junction whose Josephson energy depends on the flux ${{\Phi}}_{x}$ threading the dc-SQUID [42]:

$$\begin{equation}{E}_{J}\left({{\Phi}}_{x}\right)={E}_{J0}\mathrm{cos}\left(\pi \frac{{{\Phi}}_{x}}{{{\Phi}}_{0}}\right).\end{equation} \tag{ 71 }$$

${E}_{J0}={E}_{J0,1}+{E}_{J0,2}$ here is the sum of the energies of the two junctions in parallel, which we are assuming to be equal.

Zoom In Zoom Out Reset image size

Within this approximation, the Hamiltonian describing our circuit is:

$$\begin{align}\hfill {\hat{H}}_{e.,m.}=\sum _{i=1}^{2}{\hat{H}}_{i}+{\hat{U}}_{C}+{\hat{U}}_{M}& =\sum _{i=1}^{2}\left[\frac{{\hat{Q}}_{i}^{2}}{2{\tilde {C}}_{i}}+\frac{{\hat{{\Phi}}}_{i}^{2}}{2{\tilde {L}}_{i}}\right.\left.-{E}_{J,i}\left({{\Phi}}_{x,i}\right)\mathrm{cos}\left(2\pi \frac{{\hat{{\Phi}}}_{i}-{{\Phi}}_{z,i}}{{{\Phi}}_{0}}\right)\right]\hfill \\ \hfill & +\frac{{C}_{12}{\hat{Q}}_{1}{\hat{Q}}_{2}}{{C}_{1}{C}_{2}+\left({C}_{1}+{C}_{2}\right){C}_{12}}+\frac{{M}_{12}{\hat{{\Phi}}}_{1}{\hat{{\Phi}}}_{2}}{{L}_{1}{L}_{2}-{M}_{12}^{2}},\hfill \end{align} \tag{ 72 }$$

where ${\tilde {C}}_{1\left(2\right)}={C}_{1\left(2\right)}+{C}_{12}{C}_{2\left(1\right)}/\left({C}_{2\left(1\right)}+{C}_{12}\right)$ and ${\tilde {L}}_{1\left(2\right)}={L}_{1\left(2\right)}-{M}_{12}^{2}/{L}_{2\left(1\right)}$ [18].

Using the physical parameters given in reference [18], i.e. ${C}_{12}=132$ fF and those in table 2, and calculating the lowest four eigenvalues of our Hamiltonian for different values of mutual inductance in the range $-2\mathrm{p}\mathrm{H}{< }{M}_{12}{< }2\mathrm{p}\mathrm{H}$, we obtained the graph shown in figure 8(a). This graph matches well with the corresponding one present in figure 3(c) of reference [18]. The avoided level-crossing at ${M}_{12}\simeq 0.7$ pH is proportional to the capacitive coupling ${C}_{12}$ and only occurs at finite longitudinal fields, i.e. ${{\Phi}}_{z,i}\ne 0$ [18]. (Notice that when $-1{< }{{\Phi}}_{x,i}/{{\Phi}}_{0}{< }0$, the effective Josephson energy ${E}_{J,i}\left({{\Phi}}_{x,i}\right)$ is negative and the symmetry point where ${h}_{z}=0$ is displaced from ${{\Phi}}_{z,i}={{\Phi}}_{0}/2$ to ${{\Phi}}_{z,i}=0$ [42]).

Table 2. Physical parameters of the two-coupled-qubit system used for the simulations in figures 8 and 9.

Qubit	${E}_{J0,i}$ (GHz)	${C}_{i}$ (fF)	${L}_{i}$ (pH)	${{\Phi}}_{x,i}/{{\Phi}}_{0}$	${{\Phi}}_{z,i}/{{\Phi}}_{0}$
Q1	$1.603\cdot 1{0}^{3}$	119.5	231.9	$-0.6538$	$1\cdot 1{0}^{-4}$
Q2	$1.568\cdot 1{0}^{3}$	116.4	239	$-0.6526$	$1\cdot 1{0}^{-4}$

Zoom In Zoom Out Reset image size

It is worth noting that, in order to efficiently represent a composite circuit Hamiltonian like (72) we cannot retain the representation of the circuit operators that we used for single circuits. In that case, the size of the total Hamiltonian matrix would equal the product of the sizes of all the individual circuit Hamiltonians, and would rapidly become unmanageable. Since we are, once again, only interested in the low-energy properties of the system, a good alternative basis choice is that of the outer products of some small number ${N}_{i}$ of low-energy eigenstates of each unperturbed (i.e. non-interacting) circuit Hamiltonian ${\hat{H}}_{i}/{\hat{H}}_{c,i}$. In this case, for example, we can write:

$$\begin{equation}{\hat{H}}_{e.,m.}\simeq \sum _{i,j=0}^{{N}_{1}-1}\sum _{k,l=0}^{{N}_{2}-1}\langle {E}_{1,i}{E}_{2,k}\vert {\hat{H}}_{e.,m.}\vert {E}_{1,j}{E}_{2,l}\rangle \vert {E}_{1,i}{E}_{2,k}\rangle \langle {E}_{1,j}{E}_{2,l}\vert ,\end{equation} \tag{ 73 }$$

with the meanings of the symbols introduced before. To ensure the convergence of our results, we first used 40 harmonic oscillator number states to represent the single qubit Hamiltonians and then projected onto their ${N}_{1}={N}_{2}=10$ lowest-energy eigenstates.

Now that we have determined the low energy spectrum of the system, we can apply some reduction method to calculate the effective qubit Hamiltonian. We begin with the Schrieffer–Wolff transformation method, introduced in section 3.2.3. After verifying that the hypotheses of its construction are satisfied, in particular observing that ${{\Vert}\hat{P}-{\hat{P}}_{0}{\Vert}}_{op}< sim 0.5$ in the whole range of ${M}_{12}$, we extracted the Pauli coefficients. These were calculated by defining the computational states and the Pauli operators locally for each qubit and then projecting the effective qubit Hamiltonian on them, as shown in section 3.2.3. The six one-local coefficients are shown in figure 8(b) by solid lines. (We do not consider the coefficient ${h}_{II}=\mathrm{Tr}\left({\hat{H}}_{q}\right)$ here since we are focussing on relative energies.) The dashed lines represent the same coefficients obtained by applying the SWT reduction to the non-interacting part of the circuit Hamiltonian, i.e. to the sum of the Hamiltonians of the isolated qubits (corrected for the static inductive and capacitive loading). Since for ${h}_{zI}$ and ${h}_{Iz}$ the solid and the dashed lines overlap, the values of the longitudinal fields of the coupled system are completely determined by the static loading of the unperturbed Hamiltonians. This effect appears approximately linear in ${M}_{12}$. The values of the transverse fields for the coupled system, instead, are $\sim 25\%$ lower in magnitude than those resulting from the loaded single-qubit Hamiltonians. The interaction with the other qubit, then, has an additional effect, which we call dynamic loading. The change in transverse field appears approximately quadratic in ${M}_{12}$ and is not centred around ${M}_{12}=0$ due to the presence of the capacitive coupling (as we verified by comparing against the case ${C}_{12}=0$). As usual, the components of the local field along the y direction have been removed by making the appropriate local gauge transformation. (Actually, the circuit Hamiltonian in this case is completely real, so that no imaginary terms can appear in the reduced Hamiltonian; the gauge transformation only ensures that the signs of different coefficients are consistent across the range of ${M}_{12}$.)

The empty circles in figure 8(b) are the one-local Pauli coefficients determined with the approximate rotation method, introduced in [18] and reviewed in section 3.2.1. Comparing with the previous results, we can see that we obtain qualitatively similar, but quantitatively different results. In particular the values for the transverse fields are close to those obtained with the SWT reduction, while the new longitudinal fields are everywhere smaller in magnitude, and, in this case, do not agree with their unperturbed values (dashed lines).

Figure 9(a) shows the coefficients of the nine effective qubit Hamiltonian two-local terms. According to the reduction based on the SWT (lines), the only non-negligible terms in the Hamiltonian are those proportional to ${\boldsymbol{\sigma }}_{z,1}{\boldsymbol{\sigma }}_{z,2}$, ${\boldsymbol{\sigma }}_{x,1}{\boldsymbol{\sigma }}_{x,2}$ and ${\boldsymbol{\sigma }}_{y,1}{\boldsymbol{\sigma }}_{y,2}$. The first term represents the inductive interaction, ${\hat{U}}_{M}\propto {M}_{12}{\hat{{\Phi}}}_{1}{\hat{{\Phi}}}_{2}$, the flux being our z degree of freedom, and it indeed scales linearly with ${M}_{12}$. Since we have chosen to identify a flux degree of freedom with the real operator ${\boldsymbol{\sigma }}_{z}$, the canonically conjugate charge operator must be complex (since $\left[\hat{{\Phi}},\hat{Q}\right]=i\hslash$), and therefore must be identified with ${\boldsymbol{\sigma }}_{y}$. The YY term, then, describes the capacitive interaction and, in fact, appears to be largely independent of ${M}_{12}$. Finally the $XX$ term is a result of the presence of the higher excited states of the system [18]. It is related to both the inductive and the capacitive Hamiltonian terms and appears to scale linearly with ${M}_{12}$.

Zoom In Zoom Out Reset image size

According to reference [43], a two-local two-qubit Hamiltonian of the form

$$\begin{equation}\mathbf{H}={h}_{xI}{\boldsymbol{\sigma }}_{x,1}+{h}_{Ix}{\boldsymbol{\sigma }}_{x,2}+{h}_{zI}{\boldsymbol{\sigma }}_{z,1}+{h}_{Iz}{\boldsymbol{\sigma }}_{z,2}+{h}_{xx}{\boldsymbol{\sigma }}_{x,1}{\boldsymbol{\sigma }}_{x,2}+{h}_{yy}{\boldsymbol{\sigma }}_{y,1}{\boldsymbol{\sigma }}_{y,2}+{h}_{zz}{\boldsymbol{\sigma }}_{z,1}{\boldsymbol{\sigma }}_{z,2}\end{equation} \tag{ 74 }$$

is non-stoquastic, and remains such after arbitrary local rotations, as long as ${h}_{xI},{h}_{Ix},{h}_{zI},{h}_{Iz}\ne 0$ and $\vert {h}_{yy}\vert { >}\vert {h}_{xx}\vert ,\vert {h}_{zz}\vert$. The region where this condition is satisfied is highlighted in green in figure 9(a). Non-stoquastic two-local catalyst Hamiltonians are know to provide an exponential speed-up to the convergence of quantum adiabatic optimisation, at least with specific problem classes, including the ferromagnetic p-spin model [44]. For this reason, they might be key to establish a quantum advantage over classical optimisation routines such as quantum Monte Carlo [18, 45].

Again, our implementation of the approximate rotation reduction produces qualitatively similar results to the SWT reduction for the two-local Pauli coefficients (see hollow circles in figure 9(a)), except for ${h}_{xz}\simeq {h}_{zx}$ (purple circles), which are now of the same order of magnitude as the other coefficients. As we mentioned in section 3.2.3, the approximate rotation and the SWT reduction methods actually find equivalent effective qubit Hamiltonians, modulo a unitary. This was in fact verified by showing that both sets of coefficients lead to qubit Hamiltonians with the same spectrum.

Notice that reference [18] actually reports the two ${h}_{xz}\simeq {h}_{zx}$ coefficients to be negligible, which we ascribe to the fact that the authors used a different form for the circuit Hamiltonian, and potentially a different definition of the computational basis, and hence of ${\mathbf{R}}_{2}$, as defined in section 3.2.1 [18]. (In our case the computational basis was defined locally as shown in section 3.1.3.) In fact, any mixed two-local term, involving different Pauli operators acting on the two qubits, can be eliminated from a two-qubit Hamiltonian by performing a local change of basis [43]. Applying this transformation produces a new set of coefficients which are within 5% of those found by the SWT reduction method (see figure 9(b)). In this case, then, the unitary mapping between the two is a local transformation.

We conclude the subsection on this two-qubit system by briefly considering how, in analogy to what we had in the single-qubit case, the reduced Hamiltonian not only contains information about the system low-energy spectrum, but also about state probabilities (as well as operator matrix elements). For instance, when we set ${M}_{12}=2$ pH, the SWT reduction produces the following effective qubit Hamiltonian:

$$\begin{align}\hfill {\mathbf{H}}_{q}=& -0.125{\boldsymbol{\sigma }}_{z,1}-0.121{\boldsymbol{\sigma }}_{z,2}-0.516{\boldsymbol{\sigma }}_{x,1}-0.509{\boldsymbol{\sigma }}_{x,2}\hfill \\ \hfill & -0.459{\boldsymbol{\sigma }}_{x,1}{\boldsymbol{\sigma }}_{x,2}+0.500{\boldsymbol{\sigma }}_{y,1}{\boldsymbol{\sigma }}_{y,2}+1.079{\boldsymbol{\sigma }}_{z,1}{\boldsymbol{\sigma }}_{z,2}.\hfill \end{align} \tag{ 75 }$$

One can easily find that the first excited state of this Hamiltonian is $\vert {\mathcal{E}}_{1}\rangle =-0.0046\vert 00\rangle +0.7041\vert 01\rangle -0.7101\vert 10\rangle -0.0012\vert 11\rangle$, i.e. an entangled state where the two qubits are in opposite computational states with probability approximately one (i.e. $p\left({q}_{1}=0\vert {q}_{2}=1\right)=\cdots \simeq 1$). This should translate to the fact that, at the circuit level, there is a high probability of measuring currents of opposite sign on the two qubits, when the system is in its first excited state. In other words, if the persistent current of one of the qubits is measured to be positive, the other qubit is projected on its negative persistent current state, and vice versa. We can verify that this is actually the case by using the projectors on the positive and negative subspaces of the qubit current operators and calculating their expectation value on the first excited state $\vert {E}_{1}\rangle$ of the circuit Hamiltonian. Table 3 gives the probabilities of measuring the different computational states on $\vert {\mathcal{E}}_{1}\rangle$ and different current sign combinations on $\vert {E}_{1}\rangle$. The two results are in good agreement.

Table 3. Probabilities of measuring different computational states on the two-coupled-qubit system determined from the effective qubit model and the full circuit model.

	Probabilities
Model	${p}_{00}/{p}_{++}$	${p}_{01}/{p}_{+-}$	${p}_{10}/{p}_{-+}$	${p}_{11}/{p}_{--}$
Qubit	$2\cdot 1{0}^{-5}$	0.50	0.50	$1.3\cdot 1{0}^{-6}$
Circuit	0.06	0.44	0.45	0.05

As the final example we consider a proposed circuit implementing a three-local $ZZZ$ interaction between three flux qubits, presented in [37]. The circuit diagram is shown in figure 10(a) and consists of the three flux qubits (in this case rf-SQUID qubits) and two compound-Josephson-junction rf-SQUID couplers. The main loops of the two couplers, one of which contains a twist, mediate a magnetic interaction between the superconducting loops of qubits ${q}_{1}$ and ${q}_{2}$ (see figure 10(a)). If the flux applied to the coupler main loop, ${{\Phi}}_{z,ci}$ is kept constant, the flux applied to its dc-SQUID loop, ${{\Phi}}_{x,ci}$, controls the effective mutual inductance between the qubits and therefore the magnitude and sign of the effective $ZZ$ interaction [42]. By magnetically coupling the current loop of qubit ${q}_{3}$ to the coupler dc-SQUID loop, one can control the two local interaction between ${q}_{1}$ and ${q}_{2}$ with the current state of ${q}_{3}$, therefore obtaining a three-local ${h}_{zzz}{\hat{\sigma }}_{z1}{\hat{\sigma }}_{z2}{\hat{\sigma }}_{z3}$ interaction [37].

Zoom In Zoom Out Reset image size

The solid lines in figure 10(b) show the effective Hamiltonian coefficients for the system consisting of the three flux qubits and the single coupler ${c}_{1}$, extracted using the SWT reduction method. The main loop of the coupler and those of the three qubits are all biased at ${{\Phi}}_{z,c1}={{\Phi}}_{z,i}={{\Phi}}_{0}/2$, such that the qubit longitudinal fields are all zero. The transverse fields are also zero with the physical parameters considered (which are given below). As expected, we find a three-local interaction term $\propto {h}_{zzz}$, in addition to a residual two-local interaction between qubits ${q}_{1}$ and ${q}_{2}$, $\propto {h}_{zzI}$ and a large longitudinal field ${h}_{IIz}$ on qubit ${q}_{3}$.

The parameters used in the simulations are as follows: all qubits ($i=1,2,3$) have ${E}_{J,i}=99.3$ GHz, ${L}_{q,i}=4.5$ nH and a large shunting capacitance ${C}_{\text{sh},i}=45$ fF; the two coupler junction Josephson energies are ${E}_{J1,c1}={E}_{J2,c1}=233.4$ GHz, the coupler main loop inductance is ${L}_{z,c1}=550$ pH, while the small loop has an inductance of ${L}_{x,c1}=170$ pH and is shunted by a capacitance ${C}_{\text{sh},c1}=10$ fF; all mutual inductances are 50 pH. As in the previous simulations, the rf-SQUID qubit Hamiltonians have been expressed in a basis of 40 occupation number states. The three degrees of freedom of the coupler are expressed using 20 occupation number states for the small plasma frequency mode and 7 for the higher plasma frequency modes. The total Hamiltonian is projected on the lowest 8 unperturbed eigenstates of each qubit and on the lowest 5 unperturbed coupler eigenstates. Since the effective Hamiltonian here is diagonal in the computational basis, its coefficients can also be calculated with the method used in [37] and reviewed in section 3.2.2. The result of this reduction is represented by the filled dots in figure 10b and matches very well with the result of the SWT reduction.

Introducing a twist in the coupler, for instance changing the mutual inductance between the coupler and qubit ${q}_{2}$ from 50pH to $-50$ pH (as in coupler ${c}_{2}$), and changing the sign of the coupler x-bias, reverses the sign not only of the two-local coefficient ${h}_{zzI}$, but also of ${h}_{IIz}$. The three-local interaction coefficient, however, remains of the same sign. Therefore attaching both couplers ${c}_{1}$ and ${c}_{2}$ to the qubits leaves us with a purely three-local Hamiltonian. The numerical simulation of the full system agrees with this picture. The Pauli coefficients extracted, as a function of ${f}_{x,c1}=-{f}_{x,c2}$, are shown in figure 10(c), with the solid lines and the dots being the result of the SWT and the diagonal Hamiltonian reduction method, respectively. Coupler ${c}_{2}$ shares the same physical parameters as ${c}_{1}$ and is also biased at ${f}_{z,c2}=0.5$. Its lowest 5 unperturbed eigenstates are kept for representing the full system Hamiltonian. As we can see, the size of the three-local $ZZZ$ interaction can be changed from zero to as much as 700 MHz in the range of fluxes considered. Its sign can also be changed to negative by biasing at ${f}_{x,c1}^{\prime }=-{f}_{x,c2}^{\prime }=2-{f}_{x,c1}$ [37].

Finally we can check that the reduced Hamiltonian has the correct spectrum. This is shown in figure 10(d), where the filled dots represent the effective qubit Hamiltonian transition energies and the solid lines those of the circuit Hamiltonian. The levels are grouped in two manifolds each of four degenerate levels, separated by an energy of $2\vert {h}_{zzz}\vert$. In the ground state manifold the expectation value of the product of the qubit currents, $\langle {\hat{I}}_{1}{\hat{I}}_{2}{\hat{I}}_{3}\rangle$, and therefore $\langle {\hat{\sigma }}_{z1}{\hat{\sigma }}_{z2}{\hat{\sigma }}_{z3}\rangle$ in the reduced model, is negative, while it is positive in the excited manifold states. At energies above 8 GHz we see the additional states of the system, specifically the first excited states of the couplers. As we can see, the interaction does not close the spectral gap of the Hamiltonian, which allows us to use the Schrieffer–Wolff transformation reduction method.