# The superconducting circuit companion – an introduction with worked examples

S. E. Rasmussen<EMAIL_ADDRESS>Department of Physics and Astronomy, Aarhus
University, DK-8000 Aarhus C, Denmark  K. S. Christensen  Department of
Physics and Astronomy, Aarhus University, DK-8000 Aarhus C, Denmark  S. P.
Pedersen  Department of Physics and Astronomy, Aarhus University, DK-8000
Aarhus C, Denmark  L. B. Kristensen  Department of Physics and Astronomy,
Aarhus University, DK-8000 Aarhus C, Denmark  T. Bækkegaard  Department of
Physics and Astronomy, Aarhus University, DK-8000 Aarhus C, Denmark  N. J. S.
Loft  Department of Physics and Astronomy, Aarhus University, DK-8000 Aarhus
C, Denmark  N. T. Zinner<EMAIL_ADDRESS>Department of Physics and
Astronomy, Aarhus University, DK-8000 Aarhus C, Denmark  Aarhus Institute of
Advanced Studies, Aarhus University, DK-8000 Aarhus C, Denmark

###### Abstract

This tutorial aims at giving an introductory treatment of the circuit analysis
of superconducting qubits, i.e., two-level systems in superconducting
circuits. It also touches upon couplings between such qubits and how microwave
driving and these couplings can be used for single- and two-qubit gates, as
well as how to include noise when calculating the dynamics of the system. We
also discuss higher-dimensional superconducting qudits. The tutorial is
intended for new researchers with limited or no experience with the field but
should be accessible to anyone with a bachelor’s degree in physics. The
tutorial introduces the basic methods used in quantum circuit analysis,
starting from a circuit diagram and ending with a quantized Hamiltonian, that
may be truncated to the lowest levels. We provide examples of all the basic
techniques throughout the discussion, while in the last part of the tutorial
we discuss several of the most commonly used circuits for quantum information
applications. This includes both worked examples of single qubits and examples
of how to analyze the coupling methods that allow multiqubit operations. In
several detailed appendices, we provide the interested reader an introduction
to more advanced techniques for handling larger circuit designs.

###### Contents

  1. I Introduction
  2. II Lumped-element circuit diagrams
    1. II.1 Circuit variables
    2. II.2 Circuit components
      1. II.2.1 Capacitors
      2. II.2.2 Inductors
      3. II.2.3 Josephson junctions
      4. II.2.4 dc SQUID
      5. II.2.5 Voltage and current sources
  3. III Equations of motion
    1. III.1 Applying Kirchhoff’s laws directly
    2. III.2 Method of nodes
      1. III.2.1 Spanning tree
    3. III.3 Lagrangian approach
      1. III.3.1 Using matrices
    4. III.4 Hamiltonian approach
    5. III.5 Normal modes
    6. III.6 Change of basis
  4. IV Quantization and effective energies
    1. IV.1 Operators and commutators
    2. IV.2 Effective energies
  5. V Recasting to interacting harmonic oscillators
  6. VI Time-averaged dynamics
    1. VI.1 Interaction picture
    2. VI.2 Rotating-wave approximation
  7. VII Truncation
    1. VII.1 Two-level model (qubit)
    2. VII.2 Three-level model (qutrit)
  8. VIII Microwave driving
    1. VIII.1 Single-qubit gates
    2. VIII.2 Generalization to qudit driving
  9. IX Coupling of modes
    1. IX.1 Capacitive coupling
    2. IX.2 Two-qubit gates
    3. IX.3 Linear resonators: control and measurement
    4. IX.4 Inductive coupling
  10. X Noise and decoherence
    1. X.1 Bloch-Redfield model
      1. X.1.1 Longitudinal relaxation
      2. X.1.2 Transverse relaxation
      3. X.1.3 Pure dephasing
    2. X.2 Master equation
  11. XI Examples
    1. XI.1 Charge qubits
      1. XI.1.1 Single Cooper pair box
      2. XI.1.2 Transmon qubit
    2. XI.2 Flux qubits
      1. XI.2.1 $C$-shunted flux qubit
      2. XI.2.2 Fluxonium
      3. XI.2.3 0-$\pi$ qubit
    3. XI.3 Tunable couplers
      1. XI.3.1 Tunable capacitive coupler
      2. XI.3.2 Delft coupler
      3. XI.3.3 Gmon coupler
    4. XI.4 $XXZ$ coupling and qutrits
      1. XI.4.1 External driving of the modes
    5. XI.5 Multibody interactions
      1. XI.5.1 External coupling to eigenmodes
  12. XII Summary and outlook
  13. A Fundamental graph theory of electrical networks
    1. A.1 Circuit matrices
  14. B Method of electrical network graph theory
    1. B.1 Kirchhoff’s laws
      1. B.1.1 Kirchhoff’s current law
      2. B.1.2 Kirchhoff’s voltage law
      3. B.1.3 Reducing the number of coordinates
    2. B.2 Equations of motion
    3. B.3 Voltage and current sources
    4. B.4 Lagrangian and Hamiltonian
  15. C Exact truncation
    1. C.1 Exact four-level model

##  I Introduction

Since Richard Feynman first proposed using quantum simulators to simulate
physics [1, 2], an increasing amount of attention has been given to quantum
processors and quantum technology, something which is only expected to
increase further in the coming years [3, 4, 5]. This increase in attention has
led to swift progress within the field of quantum mechanics, taking it from
basic science research to engineering of multiqubit quantum systems capable of
performing actual calculations [6, 7, 8, 9, 10]. During this evolution, a new
discipline has emerged, coined quantum engineering, bridging the basic science
of quantum mechanics with areas traditionally considered engineering fields
[11]. It is expected that the advent of quantum engineering will lead to
computational speedups, making it possible to solve classically unsolvable
problems [12, 13, 14].

A particularly prominent platform for scalable quantum technology is
superconducting circuits used for implementing qubits or even higher-
dimensional qudits. Compared to other quantum technology schemes, such as
trapped ions [15, 16, 17, 18, 19, 20], ultracold atoms [21, 22, 23, 24, 25],
electron spins in silicon [26, 27, 28, 29, 30, 31] and quantum dots [32, 33,
34, 35, 36], nitrogen vacancies in diamonds [37, 38], or polarized photons
[39, 40, 41, 42], which all encode quantum information in microscopic systems,
such as ions, atoms, electrons, or photons, superconducting circuits are quite
different. They are macroscopic in size and printed lithographically on wafers
much similar to classical computer chips [43, 44, 45, 46, 47]. The fact that
these systems exhibit microscopic behavior, i.e., quantum-mechanical effects,
while being macroscopic in size has led to the notion of mesoscopic physics in
order to describe this intermediate scale [48, 49, 50]. A mesoscopic advantage
of superconducting circuits is the fact that microscopic features such as
energy spectra, coupling strengths, and coherence rates depend on macroscopic
circuit parameters. This means that one can design circuits such that the
properties of the resulting quantum-mechanical system, sometimes called an
artificial atom [51, 52, 53, 54, 55], can be more or less tailormade to
exhibit a particular behavior.

In this tutorial, we aim to give an introduction to circuit analysis of
superconducting qubits intended for new researchers in the field. With this,
we aim to give the tools needed for tailoring macroscopic circuits to a
desired qubit behavior. We refer to a (superconducting) qubit as the two
lowest energy levels of a superconducting circuit or subcircuit, denoted by
the Fock states $\ket{0}$ and $\ket{1}$. There are, however, several examples
of superconducting qubits which exploit higher-lying states for coupling [56]
or control [57, 58].

The field of superconducting circuits is rapidly evolving, and new theoretical
frameworks are emerging which make use of the larger Hilbert space of both
harmonic and anharmonic resonator modes, e.g., bosonic qubits [59, 60, 61, 62,
63, 64, 65, 66] or the Kerr-cat qubits [67, 68, 69], which employs the entire
circuit including drives to yield an effective potential where the two lowest
levels are coherent cat states. Such continuous variable [70, 71] qubits are
outside the scope of this tutorial, but an understanding of the fundamentals
presented in this tutorial can act as a stepping stone towards an increased
understanding of emerging superconducting circuit designs.

The present tutorial can be viewed as an introduction to more advanced reviews
of the field, such as Refs. [72, 73, 74, 51, 75, 76, 54, 77, 78, 79, 80, 50,
11], and is by no means a review of current state-of-the-art technology or
practices, but rather a detailed introduction to the theoretical methods
needed to analyze superconducting circuits in order to produce and manipulate
qubits. We do not discuss the actual experimental production of
superconducting circuits, but limit the tutorial to theoretical analysis of
such circuits. The tutorial assumes knowledge of undergraduate-level quantum
mechanics, electrodynamics, and analytical mechanics, meaning that the
tutorial should be accessible to readers with a bachelor’s degree in physics.

The tutorial is organized as follows: First, we present the basic circuit
variables and components used in the analysis in Section II. Then we present
the classical analysis used for finding the Hamiltonian of a given
superconducting circuit in Section III, where we use the method of nodes. In
Section IV we quantize the Hamiltonian and in Section V we recast the
Hamiltonian as interacting oscillators. In Section VI we discuss time-averaged
dynamics using the interaction picture. The truncation of anharmonic
oscillators is discussed in Section VII. The use of microwave driving for
control and single-qubit gates is presented in Section VIII, and the simple
coupling of modes is presented in Section IX, where two-qubit gates are
discussed as well. In Section X we introduce a method for treating noise in
open two-level quantum systems, and finally in Section XI we present a variety
of examples ranging from single qubit implementations to tunable couplers and
multibody interactions. In Section XII we present an overview of the methods
and give a perspective on where to go from here.

To students and researchers entirely new to the field of superconducting
qubits, who just want to start analyzing their first circuit, the amount of
information in this tutorial might seem extensive at first. To distill this
down to the essential information needed to get started we therefore recommend
reading Sections II, III, IV, V, and VII, which should be sufficient for
analyzing your first superconducting circuit.

##  II Lumped-element circuit diagrams

In this section we start by introducing the dynamical variables used when
analyzing superconducting circuits and then present the basic components of
the circuits.

Our analysis takes its starting point in the lumped-element model. This model
simplifies the description of a spatially distributed system (in our case a
superconducting electrical circuit) into a topology of discrete entities. We
assume that the attributes of the circuit (capacitance, inductance, and
resistance) are idealized into electrical components (capacitors, inductors,
and resistors) joined by a network of perfectly conducting wires. An example
of a lumped circuit can be seen in Fig. 1. We discuss the different components
in Section II.2.

We assume all the circuits discussed in this tutorial to be superconducting,
meaning that there is no electrical resistance in the circuit and all magnetic
fields are expelled from the wires (the Meissner effect). We therefore ignore
losses to the external environment in the following analysis. In other words,
we will consider closed quantum systems for most of this tutorial. However, a
realistic description of any quantum system should include some interactions
with the environment, as these can never be completely ignored in an
experiment. Notwithstanding, it is a good description to treat losses to the
external environment as a correction to the dynamics of the system, something
which we discuss in Section X.

Figure 1: Example of a lumped-element circuit consisting of a Josephson
junction and a capacitor in parallel connected by an inductor to another
Josephson junction and capacitor pair. Such a Josephson junction and capacitor
pair is considered a transmonlike qubit, see Section XI.1.2. An external flux
is threading the inductive loop of the circuit.

###  II.1 Circuit variables

Figure 2: Arbitrary two-terminal component on a branch, $b$, between two nodes
(dots). The voltage, $V_{b}(t)$, over the component is defined from the start
of the branch to the end of the branch. The current, $I_{b}(t)$, through the
branch is defined in the opposite direction.

Circuit analysis aims at finding the equations of motion of an electrical
circuit. Typically this means determining the current and voltage through all
components of the circuit. For simplicity, we consider only circuit networks
containing two-terminal components, i.e., components connected to two wires.
Each such component is said to lie on a _branch_ , $b$, and is characterized
by two variables at any given time $t$: The voltage, $V_{b}(t)$, across it and
the current, $I_{b}(t)$, through it. We define the orientation of the voltage
to be opposite to the direction of the current, see Fig. 2. Thus these two are
defined by the underlying electromagnetic field by

$\displaystyle V_{b}(t)$ $\displaystyle=\int_{\text{start of $b$}}^{\text{end
of $b$}}\bm{E}(t)\cdot\differential\bm{\ell},$ (1a) $\displaystyle I_{b}(t)$
$\displaystyle=\frac{1}{\mu_{0}}\oint_{b}\bm{B}(t)\cdot\differential\bm{\ell},$
(1b)

where $\mu_{0}$ is the vacuum permeability, and $\bm{E}$ and $\bm{B}$ are the
electric field inside the wire and the magnetic field outside the wire,
respectively. The closed loop in the second integral is done in vacuum
encircling the given element. As we describe the circuits in the lumped-
element model, the voltage and current are independent of the precise path the
fields are integrated along in the following sense. For the line integral of
the electric field in Eq. 1a we take the integration path to be well outside
the wire of the inductors, meaning that the magnetic field is zero along the
path. Similarly for the loop integral of the magnetic field in Eq. 1b, we take
the integration path to be well outside the dielectric of the capacitors,
meaning that the electric field is zero along the path. For more details on
the integration of electromagnetic fields see, e.g., Ref. [81].

We define the branch flux and branch charge variables as

$\displaystyle\Phi_{b}(t)$
$\displaystyle=\int_{-\infty}^{t}V_{b}(t^{\prime})\differential{t^{\prime}},$
(2a) $\displaystyle Q_{b}(t)$
$\displaystyle=\int_{-\infty}^{t}I_{b}(t^{\prime})\differential{t^{\prime}},$
(2b)

where it is assumed that the system is at rest at $t^{\prime}=-\infty$ with
zero voltages and currents. As there are less degrees of freedom in the
circuit than there are branches in the circuit, these are, just as the
currents and voltages, not completely independent but related through
Kirchhoff’s laws

$\displaystyle\sum\limits_{\begin{subarray}{c}\text{all $b$ arriving}\\\
\text{at $n$}\end{subarray}}Q_{b}=q_{n},$ (3a)
$\displaystyle\sum\limits_{\text{all $b$ around
$l$}}\Phi_{b}=\tilde{\Phi}_{l},$ (3b)

where $q_{n}$ is the charge accumulated at node $n$ and $\tilde{\Phi}_{l}$ is
the external magnetic flux through the loop $l$. A _node_ can be understood as
a point where components, or branches, converge, see Fig. 1, where we denote
nodes with a dot. We can define any circuit as a set of nodes and a set of
branches.

The notion of nodes and branches comes from graph theory, which is the natural
mathematical language for analyzing circuits. The interested reader can find
more details of fundamental graph theory and its application to electrical
circuits in Appendix A.

###  II.2 Circuit components

We consider primarily three different components of a superconducting circuit:
linear capacitors, linear inductors, and nonlinear Josephson junctions. The
two linear components should be well known to most readers, and we therefore
introduce them only briefly. The Josephson junction, on the other hand, is a
nonlinear component that is specific to superconducting circuits, and it is
the main component when working with superconducting qubits.

As we are considering superconducting circuits we do not consider resistors or
other losses. Such dissipative components are not easily included in the
Hamiltonian formalism presented in this tutorial due to their irreversible
nature. However, it can be done using, for instance, the Caldeira-Leggett
model [82, 74].

####  II.2.1 Capacitors

The first component we consider is the capacitor. For a general capacitor, the
charge on the capacitor is determined as a function of the voltage,
$q(t)=f[V(t)]$. In this tutorial, we consider only linear capacitors where the
voltage is proportional to the charge stored on the capacitor plates

$V(t)=\frac{q(t)}{C},$ (4)

where $C$ is the capacitance of the capacitor. This linear relationship is the
defining property of the linear capacitor. In reality, this is merely an
approximation, as there are small nonlinearities, which makes $C$ a function
of $q$ and $V$. These effects are usually small and therefore it is standard
to neglect them. Equation 4 can be rewritten to the flux-charge relation using
Eq. 2a as

$\dot{\Phi}(t)=V(t)=\frac{q(t)}{C},$ (5)

where the dot indicates differentiation with respect to $t$. The charge $q(t)$
is equal to the branch charge, and using Eq. 2b we find the branch current

$I(t)=C\ddot{\Phi}(t).$ (6)

The energy stored in the capacitor is found by integrating the power
$P=V(t)I(t)$ from $t=-\infty$ to $t$

$E=\frac{1}{2}C\dot{\Phi}^{2}(t).$ (7)

For superconducting circuits, typical values of the capacitances are of the
order $10\text{\,}\mathrm{fF}$. In lumped-circuit diagrams we denote the
capacitor as a pair of parallel lines, see Fig. 1.

####  II.2.2 Inductors

The time-dependent current flowing through a general inductor is a function of
the flux through it, $I(t)=f[\Phi(t)]$. For a linear inductor, the current is
proportional to the magnetic flux,

$I(t)=\dot{q}(t)=\frac{1}{L}\Phi(t),$ (8)

where $L$ is the inductance of the inductor. Integrating over the power as
before, the energy stored in the inductor is then

$E=\frac{1}{2L}\Phi^{2}(t).$ (9)

For superconducting qubits, typical values of linear inductances are of the
order $1\text{\,}\mathrm{nH}$. In lumped-circuit diagrams, we denote the
linear inductor as a coil, see Fig. 1.

Figure 3: Simple $LC$-oscillator circuit. A capacitor with capacitance $C$ is
connected in a closed circuit with an inductor of inductance $L$. The voltages
over the two components are $V_{C}$ and $V_{L}$, respectively, while the
currents are $I_{C}$ and $I_{L}$, respectively. The resulting equation of
motion is a harmonic oscillator.

As a short clarifying example we consider the classical $LC$ oscillator shown
in Fig. 3. From Kirchhoff’s current law in Eq. 3a we know that $I_{C}=I_{L}$,
where $I_{C}$ and $I_{L}$ are the currents through the capacitor and inductor,
respectively. Kirchhoff’s voltage law gives us $V_{C}=-V_{L}$, assuming no
fluctuating external flux. Using Eqs. 2a, 2b, 5, and 8 we can set up the
equations of motion for the system

$\ddot{\Phi}(t)=-\frac{1}{LC}\Phi(t),$ (10)

where we introduce $\Phi(t)=\Phi_{C}(t)=-\Phi_{L}(t)$ to get rid of the
subscripts. The system behaves as a simple harmonic oscillator in the flux.
This is analogous to a spring, where the flux is the position, and the mass
and spring constants are replaced by the capacitance and inverse inductance,
respectively.

####  II.2.3 Josephson junctions

So far we have considered only components with linear current-voltage
relations. For reasons that will become clear when we quantize the lumped
circuit, constructing a qubit from only linear components is by no means
straightforward. We therefore need nonlinear components which come in the form
of the Josephson junction. The Josephson junction plays a special role in
superconducting circuits, as it has no simple analog in a nonsuperconducting
circuit since it is related to charge quantization effects that occur in
superconductors. We start with a short introduction to superconductivity (see
Ref. [83] for more details).

When the temperature is decreased some materials undergo a phase transition
where the resistivity drops to zero. Together with the Meissner effect, i.e.,
that the material perfectly expels all magnetic fields, the perfect conduction
is the defining property of a superconductor.

The phase transition between the nonsuperconducting phase and the
superconducting phase of a material happens because the conduction electrons
condense into a so-called BCS ground state, which is characterized by an
amplitude and a phase. A priori it might seem impossible for electrons to
condense into a single quantum state since the Pauli exclusion principle
forbids this. However, as Cooper suggested, some attractive force between the
electrons leads to the formation of electron pairs [84], which have integer
spin and thus behave like bosons. This makes it possible for these so-called
Cooper pairs to condense into a single quantum ground state and in this state
the solid becomes superconducting.

Figure 4: Sketch of a Josephson junction. Two superconducting materials are
separated by a thin insulator, with a thickness of the order of
$30\text{\,}\mathrm{\SIUnitSymbolAngstrom}$. If a nonsuperconducting metal is
used as a separator, it can be several micrometers wide. Cooper pairs can
tunnel back and forth between the two superconducting materials.

A Josephson junction consists of two superconducting islands separated by a
thin insulator, a nonsuperconducting metal, or a narrow superconducting wire.
Cooper pairs can then tunnel through the barrier from one island to the other,
a phenomenon known as the Josephson effect [85, 86], see Fig. 4. The tunneling
rate (current) and the voltage between the two islands depends on the
superconducting phase difference, $\phi$, between the islands through [87]

$\displaystyle I(t)$ $\displaystyle=I_{c}\sin[\phi(t)],$ (11) $\displaystyle
V(t)$ $\displaystyle=\frac{\hbar}{2e}\dot{\phi},$ (12)

where $I_{c}$ is the critical current of the junction, which depends on the
junction geometry. Equation 12 allows us to relate the junction phase
difference to the generalized flux through $\Phi=\hbar\phi/2e$. The charge and
flux are thus related through

$\dot{q}(t)=I_{c}\sin\left(2\pi\dfrac{\Phi(t)}{\Phi_{0}}\right),$ (13)

where we define the magnetic flux quantum $\Phi_{0}=h/2e$. The Josephson
junction works as a flux-dependent inductor with inductance given by [73]

$L(\Phi)=\left(\partialderivative{I}{\Phi}\right)^{-1}=\dfrac{L_{J}}{\cos\left(2\pi\dfrac{\Phi}{\Phi_{0}}\right)},$
(14)

where we define the Josephson inductance $L_{J}=\Phi_{0}/2\pi I_{c}$. Since
the inductance is associated with the inertia of the Cooper pairs it is often
referred to as kinetic inductance. See Section XI.2.2 for details on the use
of large kinetic inductance. For superconducting qubits, typical values of
Josephson inductances are of the order $100\text{\,}\mathrm{nH}$. The energy
of a Josephson junction is also nonlinear. We have

$E=\frac{\Phi_{0}^{2}}{(2\pi)^{2}}\frac{1}{L_{J}}\left[1-\cos\left(2\pi\frac{\Phi}{\Phi_{0}}\right)\right],$
(15)

where we often neglect the constant term when dealing with the Lagrangian or
Hamiltonian, as it is irrelevant for the dynamics of the system. We define the
factor in front of the bracket to be the Josephson energy of the Josephson
junction, $E_{J}=\Phi_{0}^{2}/(2\pi)^{2}L_{J}=\Phi_{0}I_{c}/2\pi$. In this
tutorial we denote Josephson junctions as a boxed ”x” in lumped-circuit
diagrams, see Fig. 1. In the literature sometimes an ”x” without a box is
used.

It is conventional to simplify notation in a way such that charges and fluxes
become dimensionless. This is done by using units where

$\hbar=2e=1\qquad\text{and thus}\qquad\frac{\Phi_{0}}{2\pi}=1.$ (16)

This means that we get rid of the cumbersome factor of $2\pi/\Phi_{0}$ in the
sinusoidal Josephson junction terms. Note that in this convention the units of
capacitance and inductance become inverse energy. Moreover, with this choice
of units the junction phase differences are equal to the generalized flux
$\phi=\Phi$, and the energy of a Josephson junction becomes equal to the
critical current, $E_{J}=I_{c}$.

####  II.2.4 dc SQUID

Figure 5: A dc superconducting quantum interference device (dc SQUID). (a)
Implementation of a dc SQUID. (b) Corresponding circuit diagram.

It is often desirable to be able to tune the parameters of the circuit
externally. Therefore many circuits employ a direct current superconducting
quantum interference device, or dc SQUID, instead of a single Josephson
junction. A dc SQUID consists of two Josephson junctions on a ring, with an
external magnetic field, $\tilde{\Phi}$, through the ring [88], see Fig. 5(a).
While this does not change the form of the energy of the Josephson junction,
it has the advantage that it makes the front factor in Eq. 15 tunable. To see
this consider the circuit diagram in Fig. 5(b). The energy of this component
must be the sum of two Josephson junctions

$U=-E_{J}\cos\left(\Phi_{L}+\frac{\tilde{\Phi}}{2}\right)-E_{J}\cos\left(\Phi_{R}+\frac{\tilde{\Phi}}{2}\right),$
(17)

where $\Phi_{L/R}$ is the branch flux of the left and right branch,
respectively, and we divide the external flux equally between the two arms of
the dc SQUID following Kirchhoff’s voltage law in Eq. 3b. Note that here we
consider symmetrical junctions, but it is a neat exercise to extend it to
asymmetrical junctions.

Since we are considering the arms of a loop, we can write
$\Phi=\Phi_{L}=-\Phi_{R}$ in Eq. 17. Using the trigonometric identity
$2\cos\alpha\cos\beta=\cos(\alpha-\beta)+\cos(\alpha+\beta)$ with
$\alpha=\tilde{\Phi}/2$ and $\beta=\Phi$, we can rewrite Eq. 17 into the form

$U=-2E_{J}\cos\left(\frac{\tilde{\Phi}}{2}\right)\cos\Phi.$ (18)

The so-called fluxoid quantization condition states that the algebraic sum of
branch fluxes of all the inductive elements along the loop plus the externally
applied flux must equal an integer number of superconducting flux quanta [89,
90, 11], i.e.,

$\Phi+\tilde{\Phi}=2\pi k,$ (19)

where $k$ is an integer. Together with Kirchhoff’s voltage law in Eq. 3b this
means that we can remove a degree of freedom. This explains how one goes from
two branch fluxes, $\Phi_{L/R}$, to just one branch flux, $\Phi$, since the
branch fluxes are the system degrees of freedom. In other words we obtain, an
effective Josephson energy of
$E^{\prime}_{J}(\tilde{\Phi})=2E_{J}|\cos(\tilde{\Phi}/2)|$, where the
Josephson energy can be dynamically tuned through the external flux,
$\tilde{\Phi}$. This idea is often implemented in superconducting circuits
instead of a single Josephson junction so that the spacing of the energy
levels can be tuned dynamically by tuning $\tilde{\Phi}$. However, we usually
just place a single Josephson junction in a circuit diagram. Due to the
sensitivity of the dc SQUID it has many uses especially in clinical
applications such as magnetoencephalography [91, 92], magnetocardiography, and
magnetic resonance imaging (MRI), where they are used for detecting tiny
magnetic fields in living organisms [93, 94].

####  II.2.5 Voltage and current sources

We can treat constant voltage and current sources by representing them as
capacitors or inductors. Consider a constant voltage source $V$. This can be
represented by a very large but finite capacitor, in which an initially large
charge $Q$ is stored such that $V=Q/C$ in the limit where
$C\rightarrow\infty$. Similarly, a constant current source can be represented
by a very large but finite inductor, in which an initially large flux $\Phi$
is stored, such that $I=\Phi/L$ in the limit where $L\rightarrow\infty$.

##  III Equations of motion

In order to describe the dynamics of the lumped-circuit diagrams we presented
in the previous section, we now determine the equations of motion for the
systems. The equations of motion depend on the circuit components and can be
written in terms of the circuit variables using either the voltage and current
in Eq. 1 or equivalently using the flux and charge in Eq. 2. There are several
ways of finding the equations of motion, and we start from the simplest
approach; applying Kirchhoff’s laws directly to the circuit. From this
starting point, we then progress to the method of nodes and then to the
Lagrangian and Hamiltonian.

###  III.1 Applying Kirchhoff’s laws directly

The simplest way to find the equations of motion for a given circuit is to
apply Kirchhoff’s laws. We have already done this for the simple $LC$
oscillator example in Fig. 3, which yielded the harmonic oscillator equation
of motion in Eq. 10. To get a better feel for this procedure, let us consider
a few additional examples.

Figure 6: A Josephson junction and capacitor circuit. A capacitor with
capacitance $C$ is connected in a closed circuit with a Josephson junction,
$E_{J}$. The voltage over the two components is $V_{C}$ and $V_{J}$,
respectively, while the currents are $I_{C}$ and $I_{J}$, respectively. The
resulting equation of motion is a Duffing oscillator.

The next natural step is to exchange the linear inductor in Fig. 3 with a
nonlinear Josephson junction. This yields the circuit in Fig. 6. From
Kirchhoff’s current law in Eq. 3a we know that $I_{C}=I_{J}$, where $I_{C}$
and $I_{J}$ are the currents through the capacitor and Josephson junction,
respectively. Kirchhoff’s voltage law implies $V_{C}=-V_{J}$. Using Eqs. 2a,
2b, 5, and 8 we can set up the equations of motion for the system,

$\ddot{\Phi}(t)=-\frac{I_{c}}{C}\sin\Phi(t),$ (20)

where we introduce $\Phi(t)=\Phi_{C}(t)=-\Phi_{J}(t)$. Equation 20 is
identical to the equation of motion for a simple pendulum, with the critical
current, $I_{c}$, playing the role of the gravitational constant and the
capacitance, $C$, becoming the mass of the pendulum, similar to the case of
the $LC$ circuit, see Eq. 10, which is the lowest order approximation to Eq.
20. Contrary to Eq. 10 this is not linear in $\Phi$, which is an effect of the
introduction of the nonlinear Josephson junction.

Figure 7: Example circuit of Fig. 1 with explicit directions of the voltages
and currents shown. All loops are propagated counterclockwise.

We now continue to the more complicated example of Fig. 1. This time
Kirchhoff’s voltage law gives us three equations, one for each loop of the
circuit. We denote the left capacitor and Josephson junction $C_{1}$ and
$E_{J,1}$, respectively. Similarly, we have to the right $C_{2}$, $E_{J,2}$.
The connecting inductor is denoted by $L_{12}$. Defining the direction of the
current and voltages as in Fig. 7, we find the following equations from Eq. 3b

$\displaystyle-\Phi_{EJ1}-\Phi_{C1}$ $\displaystyle=0,$ (21a)
$\displaystyle\Phi_{EJ2}+\Phi_{C2}$ $\displaystyle=0,$ (21b)
$\displaystyle\Phi_{EJ1}-\Phi_{EJ2}+\Phi_{L}$ $\displaystyle=\tilde{\Phi},$
(21c)

where $\tilde{\Phi}$ is the external flux in the inductor loop. We propagate
all loops counterclockwise, which yields negative signs on the terms in Eq. 21
when the voltage of the given branch is in the opposite direction to the loop
direction. Note that we can also include external fluxes in the two capacitive
loops. However, as we will see in Section III.3, as long as we consider only
time-independent fluxes, the external fluxes will only be relevant in purely
inductive loops. From Eqs. 21a and 21a we define
$\Phi_{1}=\Phi_{EJ1}=-\Phi_{C1}$ and $\Phi_{2}=\Phi_{EJ2}=-\Phi_{C2}$. Using
this we can also express the flux through the inductor as
$\Phi_{L}=\Phi_{2}-\Phi_{1}+\tilde{\Phi}$, which significantly reduces the
number of variables.

From Kirchhoff’s current law we find the following equations

$-I_{Cn}+I_{EJn}=\mp I_{L},$ (22)

for $n=1,2$, where the minus is for $n=1$ and the plus is for $n=2$. Inserting
the current relations for the respective components, we find the following
equations of motion

$C_{n}\ddot{\Phi}_{n}=\mp\frac{1}{L_{12}}(\Phi_{2}-\Phi_{1}+\tilde{\Phi})-E_{J,n}\sin\Phi_{n},$
(23)

for $n=1,2$.

The end goal of our analysis is to quantize the circuit to treat it quantum
mechanically. When doing quantum mechanics we are usually interested in the
Hamiltonian of the system, as it is closely related to the energy spectrum and
time evolution of the system. It is possible to infer the system Hamiltonian
from the equations of motion. This is usually done by finding a Lagrangian
that yields the equation of motion using Lagrange’s equations [see Eq. 27] and
then performing a Legendre transformation.

While the approach of applying Kirchhoff’s law directly always yields the
correct equations of motion, it quickly becomes cumbersome as the circuits
increase in complexity. We, therefore, seek a method for determining the
Lagrangian directly. This can be achieved using the method of nodes.

###  III.2 Method of nodes

In this section, we present the method of nodes which solves most practical
problems involving Josephson junctions. The discussion follows the method
proposed by Devoret [95, 50].

Our main obstacle when determining the Lagrangian of a given circuit is to
remove superfluous degrees of freedom and determine how to include the
external fluxes. As we saw above, we can solve these problems by manipulating
Kirchhoff’s law. Here we present an alternative approach.

We have already defined a node as a point where one or more components
connect. We now further define a _ground node_ as a node connected to ground.
These nodes are inactive since the flux through them is zero and thus they do
not contribute to the dynamics of the system, and can thus be ignored. For the
remaining nodes, we distinguish between active and passive nodes. An _active
node_ is defined as a node where at least one capacitor and one inductor
(either linear and Josephson junction) meet. A _passive node_ is defined as a
node where only one type of component meet, either only capacitors or only
inductors. It turns out that passive nodes represent superfluous degrees of
freedom and therefore only yield constraints on the dynamics of the system.
This is similar to how one may determine an effective capacitance for a serial
or parallel collection of individual capacitances.

Considering the example circuit in Fig. 1, we can represent the circuit as a
set of branches, $\mathcal{B}$, and a set of nodes, $\mathcal{N}$. The set of
nodes consists of three nodes; two active nodes and a ground node. The set of
branches is equal to the set of components in the circuit, i.e., the example
circuit has five branches; two capacitor branches, two Josephson junction
branches, and a single linear inductor branch.

We call such a representation consisting of a set of nodes and a set of
branches a _network graph_ or simply a _graph_. With this notation, we can
divide the circuit into subgraphs. For a given circuit there are many possible
subgraphs, but we focus on the capacitive subgraph and the inductive subgraph.
The _capacitive subgraph_ contains only branches of capacitors and the nodes
connected to such branches. The _inductive subgraph_ contains only branches of
inductors and nodes connected to such branches. In the example circuit in Fig.
1 the two capacitor branches and all three nodes are in the capacitive
subgraph, while the inductor branch and the two Josephson junction branches
are in the inductive subgraph together with the three nodes. Notice how the
nodes can be in both subgraphs at the same time, see Fig. 8.

Figure 8: Highlight of (a) the capacitive subgraph and (b) the inductive
subgraph of the example circuit in Fig. 1.

The capacitive subgraph consists only of linear capacitors, and thus we can
express the energy of a capacitive branch in terms of the voltages, i.e., the
derivative of the flux, using Eq. 7. By doing this we have now broken the
symmetry between the charge and flux, and the flux can now be viewed as the
“position.” With this treatment, the capacitive energy becomes equivalent to
the kinetic energy, while the inductive energy becomes equivalent to the
potential energy.

This symmetry breaking also explains why passive nodes do not contribute to
the dynamics of the system. A passive node in between two inductors does not
have any kinetic energy and can therefore be considered stationary. On the
other hand, a node in between two capacitors does have kinetic energy, but no
potential energy, and can therefore be considered a free particle which does
not interact with the rest of the system.

However, any realistic inductor (both linear and nonlinear) will always
introduce some capacitance since a capacitance occurs whenever two conducting
materials are in close proximity to each other. Consider the Josephson
junction in Fig. 4, it quite closely resembles a linear plate capacitor, thus
it is expected that some parasitic capacitance will be present in parallel
with the inductor. Nonetheless, we can often make this parasitic capacitance
so small that it can be neglected in the lumped-element circuit. One should,
however, be aware of these capacitances when designing superconducting
circuits.

####  III.2.1 Spanning tree

We are now ready to consider the most important subgraph of the circuits: the
spanning tree. The _spanning tree_ is constructed by connecting every node in
the circuit to each of the other nodes by only one path. See Definition 3 in
Appendix A for a more mathematical definition using graph theory. Note that
there are often several choices for the spanning tree. This is not a problem
for the analysis and can be seen analogous to the choice of a particular gauge
in electromagnetic field theory or the choice of a coordinate system in
classical mechanics.

Choosing a spanning tree for a given circuit partitions the branches into two
sets: The set of branches on the spanning tree, $\mathcal{T}$, and its
complementary set, $\bar{\mathcal{T}}=\mathcal{B}\setminus\mathcal{T}$, i.e.,
the branches _not_ on the spanning tree. We call the latter set the set of
closure branches because its branches close the loop of the spanning tree.

We use the spanning tree to determine where to include the external fluxes of
the system. Following Kirchhoff’s laws, the flux $\phi_{n}$ of node $n$ can be
written as the sum of incoming and outgoing branch fluxes, with a suitable
sign depending on the direction of the flux. With this in mind, we can write
the branch fluxes in terms of the node fluxes

$\displaystyle\Phi_{b\in\mathcal{T}}$
$\displaystyle=\phi_{n}-\phi_{n^{\prime}},$ (24a)
$\displaystyle\Phi_{b\in\bar{\mathcal{T}}}$
$\displaystyle=\phi_{n}-\phi_{n^{\prime}}+\tilde{\Phi},$ (24b)

where $n$ and $n^{\prime}$ are the nodes at the start and end of the given
branch, respectively, and $\tilde{\Phi}$ is the external flux through the loop
closed by the branch. Note that the external flux occurs only if the branch is
a closure branch. The fact that external fluxes do not appear in every branch
is due to Kirchhoff’s law in Eq. 3b, which eliminates the external flux on
some of the branches. One can therefore choose onto which branches these
external fluxes should be included, as long as Eq. 3b is satisfied, which is
exactly the choice we make by choosing the spanning tree.

Substituting the node fluxes into the expressions for the energy of the
different components, i.e., into Eqs. 7, 9, and 15, we can express the
energies as a function of the node fluxes. The results can be seen in Table 1.

Table 1: Energies of different components on either the spanning tree or a closure branch of the circuit. The magnetic flux through the closure branch due to external fields is denoted $\tilde{\Phi}_{b}$. The time derivative of the magnetic flux is included for linear capacitors on closure branches for completeness. For the rest of this tutorial, we assume time-independent external fluxes, i.e., $\dot{\tilde{\Phi}}_{b}=0$. We refer to Refs. [96, 97] for a discussion of time-dependent fluxes. Element | Spanning tree | Closure branch
---|---|---
Linear capacitor | $\dfrac{C}{2}(\dot{\phi}_{n}-\dot{\phi}_{n^{\prime}})^{2}$ | $\dfrac{C}{2}(\dot{\phi}_{n}-\dot{\phi}_{n^{\prime}}+\dot{\tilde{\Phi}}_{b})^{2}$
Linear inductor | $\dfrac{1}{2L}(\phi_{n}-\phi_{n^{\prime}})^{2}$ | $\dfrac{1}{2L}(\phi_{n}-\phi_{n^{\prime}}+\tilde{\Phi}_{b})^{2}$
Josephson junction | $-E_{J}\cos\left(\phi_{n}-\phi_{n^{\prime}}\right)$ | $-E_{J}\cos(\phi_{n}-\phi_{n^{\prime}}+\tilde{\Phi}_{b})$

Note that if the circuit contains only time-independent external fluxes, it is
often an advantage to choose a spanning tree containing as few capacitors as
possible, such that the capacitors lie on the closure branches. The reason is
that a time-independent external flux disappears from capacitive terms since
$\dot{\tilde{\Phi}}=0$. When working with time-independent external fluxes
these are therefore only relevant in purely inductive loops. Time-dependent
external fluxes are beyond the scope of this tutorial, see Refs. [96, 97] for
a treatment of this case.

If we consider the example circuit in Fig. 1 we can choose the spanning tree
in many different ways. Since we consider only time-independent external
fluxes, a particularly nice choice of spanning would be over the two Josephson
junction (JJ) branches, which means that any external flux will appear only in
the linear inductor term. For this reason, we do not need to worry about any
external fluxes through the two capacitive loops.

###  III.3 Lagrangian approach

Having chosen a spanning tree for our circuit, we are now ready to determine
its Lagrangian. The Lagrangian is found by subtracting the potential
(inductive) energies from the kinetic (capacitive) energies

$\mathcal{L}=T-U=T_{\text{cap}}-U_{\text{ind}}-U_{\text{JJ}},$ (25)

where $T$ is the kinetic energy and $U$ is the potential energy. The
subscripts indicate the type of element each term refers to.

With the definition of the Lagrangian and the energies of Table 1, we can
write the Lagrangian for the example circuit in Fig. 1 as

$\displaystyle\mathcal{L}$
$\displaystyle=\frac{C_{1}}{2}\dot{\phi}_{1}^{2}+\frac{C_{2}}{2}\dot{\phi}_{2}^{2}-\frac{1}{2L_{12}}(\phi_{2}-\phi_{1}+\tilde{\Phi})^{2}$
(26) $\displaystyle\phantom{=}+E_{J,1}\cos\phi_{1}+E_{J,2}\cos\phi_{2},$

where $C_{n}$ and $E_{J,n}$ is the capacitance and Josephson energy of the
capacitor and Josephson junction, respectively. The index $n=1,2$ corresponds
to the left and right side, respectively. The inductance of the inductor is
denoted $L_{12}$. With the Lagrangian, one can obtain the equations of motion
from Lagrange’s equations

$\derivative{t}\partialderivative{\mathcal{L}}{\dot{\phi}_{n}}=\partialderivative{\mathcal{L}}{\phi_{n}}.$
(27)

Applying this to the example circuit, we find the equations of motion

$\ddot{\phi}_{n}=\mp\frac{1}{L_{12}C_{n}}(\phi_{2}-\phi_{1}+\tilde{\Phi})-\frac{E_{J,n}}{C_{n}}\sin\phi_{n},$
(28)

where the minus is for $n=1$ and the plus is for $n=2$. This is identical to
Eq. 23 written up with node fluxes instead of branch fluxes.

####  III.3.1 Using matrices

Writing the Lagrangian as in Eq. 26 can be rather tedious for larger circuits
since it includes a lot of sums. We, therefore, seek a more elegant way to
write the Lagrangian. This is achieved using matrix notation. First we list
all the nodes 1 to $N$ and define a flux column vector
$\bm{\phi}^{T}=(\phi_{1},\dots,\phi_{N})$, where $T$ indicates the transpose
of the vector. Note that for a grounded circuit we do not include the ground
node since its flux equals zero and it does not contribute to the true degrees
of freedom in any case. We can always choose a ground node in our circuits as
one mode will always decouple from the remaining modes for ungrounded
circuits, see Section III.6.

We are now ready to set up the capacitive matrix $\bm{C}$ of the system. The
nondiagonal matrix elements are minus the capacitance, $C_{jk}$, connecting
nodes $j$ and $k$. The diagonal elements consist of the sum of the nondiagonal
values in the corresponding row or column, multiplied by $-1$, i.e.,
$C_{jj}=\sum_{k\neq j}C_{jk}$. If a node is connected to ground via a
capacitor, this capacitance must also be added to the diagonal element. With
this $N\times N$ matrix we can write the kinetic energy term as

$T=\frac{1}{2}\dot{\bm{\phi}}^{T}\bm{C}\dot{\bm{\phi}}.$ (29)

In the case of the example circuit in Fig. 1, the flux column vector is
$\bm{\phi}^{T}=(\phi_{1},\phi_{2})$, and the capacitive matrix becomes

$\bm{C}=\begin{bmatrix}C_{1}&0\\\ 0&C_{2}\end{bmatrix}.$ (30)

We now consider the contribution from the linear inductors. We set up the
inductive matrix $\bm{L^{-1}}$ in the same way as the capacitive matrix. The
nondiagonal elements are $-1/L_{jk}$ if an inductance $L_{jk}$ connects nodes
$j$ and $k$, and zero otherwise, while the diagonal elements consist of the
sum of values in the corresponding row or column, multiplied by minus one,
$1/L_{jj}=\sum_{k\neq j}1/L_{jk}$. If a node is connected to the ground via an
inductor, this inductance must also be added to the diagonal element. Of
course, if no inductor is connecting two nodes, the element should be zero. We
must also include the external magnetic flux in this term. Thus the energy due
to linear inductors becomes

$U_{\text{ind}}=\frac{1}{2}\bm{\phi}^{T}\bm{L^{-1}}\bm{\phi}+\sum_{b\in\bar{\mathcal{T}}}\frac{1}{L_{b}}(\phi_{n}-\phi_{n^{\prime}})\tilde{\Phi}_{b},$
(31)

where we remove all irrelevant constant terms. The second term sums over all
the inductive closure branches of the circuit, where $n$ and $n^{\prime}$ are
the nodes connected by branch $b$.

If we consider the example circuit again, the inductive matrix is

$\bm{L^{-1}}=\begin{bmatrix}1/L_{12}&-1/L_{12}\\\
-1/L_{12}&1/L_{12}\end{bmatrix},$ (32)

where $L_{12}$ is the inductance of the linear inductor. With this the
inductive energy of the example circuit becomes

$U_{\text{ind}}=\frac{1}{2}\bm{\phi}^{T}\bm{L^{-1}}\bm{\phi}+\frac{1}{L_{12}}(\phi_{1}-\phi_{2})\tilde{\Phi},$
(33)

where $\tilde{\Phi}$ is the external flux through the inductive loop. When
there are only a few linear inductors, as in the example circuit, it might be
more straightforward to write the energy without the matrix notation. We do
not attempt to write the Josephson junction terms using matrix notation as
they are nonlinear functions of the node flux variables.

###  III.4 Hamiltonian approach

The Hamiltonian of the circuit can be found by a simple transformation of the
Lagrangian through what is commonly referred to as a Legendre transformation.
First, we define the conjugate momentum to each node flux by

$q_{n}=\partialderivative{\mathcal{L}}{\dot{\phi}_{n}},$ (34)

which in vector form becomes $\bm{q}=\bm{C}\dot{\bm{\phi}}$. If the
capacitance matrix is invertible we can express $\dot{\bm{\phi}}$ as a
function of $\bm{q}$. We denote the conjugate momenta as node charges since
they correspond to the algebraic sum of the charges on the capacitances
connected to node $n$.

The Hamiltonian can now be expressed in terms of the node charges, $q_{n}$,
for the kinetic energy and node fluxes, $\phi_{n}$, for the potential energy
through the Legendre transform

$\displaystyle\mathcal{H}$
$\displaystyle=\dot{\bm{\phi}}^{T}\bm{q}-\mathcal{L}$ (35)
$\displaystyle=\frac{1}{2}\bm{q}^{T}\bm{C}^{-1}\bm{q}+U(\bm{\phi}),$

where the potential energy is a nonlinear function of the node fluxes. Note
that the functional form of the Hamiltonian may differ depending on the choice
of spanning tree. This is because the choice of flux-node coordinates is not
unique, much like the electrodynamic potentials, which have a “gauge freedom”
in which certain functions can be added to the potentials without any change
to the physics, or more concretely; without changes to the electric and
magnetic fields [81]. Here a different choice of flux variables would
correspond to a change of gauge as well and a physical quantity like the total
energy should not change under such a transformation.

With the Hamiltonian, it is possible to find the equations of motion using
Hamilton’s equations

$\dot{\phi}_{n}=\partialderivative{\mathcal{H}}{q_{n}},\qquad\dot{q}_{n}=-\partialderivative{\mathcal{H}}{\phi_{n}},$
(36)

which yields results for the equations of motion that are equivalent to
Lagrange’s equations Eq. 27.

###  III.5 Normal modes

Lagrange’s equations tell us that for all passive nodes $\dot{q}_{n}=0$, since
for a passive note we have $\partialderivative{\mathcal{H}}{\phi_{n}}=0$. This
means that the circuit has at most the same number of true degrees of freedom
as the number of active nodes except the ground node. The number of true
degrees of freedom turn out to be identical to the number of normal modes of
the system. If all inductors can be approximated as linear inductors (and
external fluxes are ignored), the Lagrangian takes the form

$\mathcal{L}=\frac{1}{2}\bm{\dot{\phi}}^{T}\bm{C}\bm{\dot{\phi}}-\frac{1}{2}\bm{\phi}^{T}\bm{L^{-1}}\bm{\phi}.$
(37)

This simple form of the Lagrangian means that the equations of motion become

$\bm{C}\bm{\ddot{\phi}}=-\bm{L^{-1}}\bm{\phi},$ (38)

which is essentially Hooke’s law in matrix form where the capacitances play
the role of the masses and inductances play the role of the spring constants
[98]. The normal modes of the full systems can be found as the eigenvectors of
the matrix product $\bm{\Omega}^{2}=\bm{C}^{-1}\bm{L^{-1}}$ associated with
nonzero eigenvalues. These nonzero eigenvalues correspond to the squared
normal mode frequencies of the circuit. Note that $\bm{C}^{-1}$ and
$\bm{L^{-1}}$ can always be diagonalized simultaneously since they are both
positive definite matrices [98]. It can be advantageous to find these
eigenmodes and use them as a basis as it reduces the number of couplings
between modes.

###  III.6 Change of basis

Here we present a method for changing into the normal mode basis of a circuit.
Given a circuit with $N$ nodes and a matrix product
$\bm{\Omega}^{2}=\bm{C^{-1}}\bm{L^{-1}}$, let
$\bm{v}_{1},\bm{v}_{2},\dots,\bm{v}_{n}$ be the orthonormal eigenvectors of
$\bm{\Omega}^{2}$, with eigenvalues $\chi_{1},\chi_{2},\dots,\chi_{n}$. Let
$\bm{\phi}$ be the usual vector of the node fluxes of the circuit. We can then
introduce the normal modes $\bm{\psi}$ via

$\displaystyle\bm{\phi}=\bm{\mathcal{V}}\bm{\psi},$ (39)

where

$\displaystyle\bm{\mathcal{V}}=\begin{bmatrix}\mid&\mid&&\mid\\\
\bm{v}_{1}&\bm{v}_{2}&\cdots&\bm{v}_{N}\\\ \mid&\mid&&\mid\end{bmatrix},$ (40)

is a matrix whose columns are the eigenvectors of $\bm{\Omega}^{2}$. The
kinetic energy term in Eq. 29 can now be written

$\displaystyle T=\frac{1}{2}\dot{\bm{\psi}}^{T}\bm{K}\dot{\bm{\psi}},$ (41)

where we introduce the capacitance matrix in the transformed coordinates
$\bm{K}=\bm{\mathcal{V}}^{T}\bm{C}\bm{\mathcal{V}}$.

While we assume the columns of Eq. 40 to be the eigenvectors of
$\bm{\Omega}^{2}$, this is not a requirement, and one can rotate to any frame
using an orthonormal basis to construct $\bm{\mathcal{V}}$. However, only if
one uses the eigenvectors of $\bm{\Omega}^{2}$ will the transformed
capacitance matrix, $\bm{K}$, be diagonal with entries $\lambda_{i}$. In terms
of the canonical momenta $\bm{p}$ conjugate to $\bm{\psi}$, the kinetic energy
takes the usual form

$\displaystyle T=\frac{1}{2}\bm{p}^{T}\bm{K}^{-1}\bm{p},$ (42)

where the inverse of $\bm{K}$ is trivial to find if it is a diagonal matrix,
yielding the entries $1/\lambda_{i}$. In the above, we have assumed that
$\bm{C}$ is positive definite, which is usually the case. This means that
$\lambda_{i}\neq 0$. We comment on the case where $\bm{C}$ is not positive
definite below.

We must also consider how contributions from the higher-order terms of
inductors behave under this coordinate transformation. Even though we have
approximated all inductors as linear to find the normal modes, higher-order
terms from Josephson junctions still contribute as corrections, often leading
to couplings between the modes. Such terms transform the following way

$\displaystyle\phi_{k}-\phi_{l}\rightarrow\sum_{i}\left[(\bm{v}_{i})_{k}-(\bm{v}_{i})_{l}\right]\psi_{i},$
(43)

where $(\bm{v}_{i})_{k}$ is the $k$th entry of $\bm{v}_{i}$. Considering for
instance fourth-order terms in $\phi$, this can result in both two-body
interactions as well as interactions beyond two body. These multibody
interactions can complicate the equations of motion beyond what the change of
basis adds in terms of simplification. Coordinate transformations are
therefore often most useful in cases where the capacitors are symmetrically
distributed, which results in simple normal modes.

The center-of-mass (CM) mode plays a special role in analytical mechanics, as
it often decouples from the dynamics of the system. The same is the case for
electrical circuits. The center-of-mass mode corresponds to
$\bm{v}_{\text{CM}}=(1,1,\dots,1)^{T}/\sqrt{N}$, which yields
$\psi_{\text{CM}}=(\phi_{1}+\phi_{2}+\cdots+\phi_{N})/\sqrt{N}$. This mode is
always present and it corresponds to charge flowing equally into every node of
the circuit from ground and oscillating back and forth between ground and the
nodes. Furthermore, since all its entries are identical it always disappears
in the linear combination of Eq. 43
[$(\bm{v}_{\text{CM}})_{k}-(\bm{v}_{\text{CM}})_{l}=0$]. Hence, this mode is
completely decoupled from the dynamics.

The decoupling of this mode is related to how we can arbitrarily choose a node
in our circuit as the ground node, whose node flux does not enter into our
equations, or rather is identically set to zero. For an ungrounded circuit
$\bm{C}$ is no longer positive definite and we have $\lambda_{\text{CM}}=0$,
making $\bm{K}$ singular. We therefore always assume the circuit is grounded
such that $\lambda_{\text{CM}}\neq 0$.

For an example of multibody interactions see Section XI.5. For other examples
of changes of basis see Sections XI.3.2 and XI.2.3.

##  IV Quantization and effective energies

###  IV.1 Operators and commutators

We now quantize the classical Hamiltonian to obtain a quantum-mechanical
description of the circuit. This is done through canonical quantization,
replacing all the variables and the Hamiltonian with operators

$\displaystyle\phi_{n}$ $\displaystyle\rightarrow\hat{\phi}_{n},$ (44)
$\displaystyle q_{n}$ $\displaystyle\rightarrow\hat{q}_{n},$
$\displaystyle\mathcal{H}$ $\displaystyle\rightarrow\hat{\mathcal{H}},$

where $\hat{\phi}_{n}$ is the node flux operator corresponding to position
coordinates, $\hat{q}_{n}$ is the conjugate momentum, and $\hat{\mathcal{H}}$
is the Hamiltonian operator. If the flux operator and the conjugate momentum
operator are not constants of motion they obey the canonical commutator
relation

$[\hat{\phi}_{n},\hat{q}_{m}]=\hat{\phi}_{n}\hat{q}_{m}-\hat{q}_{m}\hat{\phi}_{n}=i\hbar\delta_{nm},$
(45)

where $\delta_{nm}$ is the Kronecker delta. The commutator relation in Eq. 45
does not hold if a given node, $n$, is not a true degree of freedom. This can
happen in case the variable does not appear in $\hat{\mathcal{H}}$, and
therefore the commutator between the variable and the Hamiltonian will be
zero. This means that $\hat{\phi}_{n}$ or $\hat{q}_{n}$ will be constant of
motion according to Heisenberg’s equation of motion. This is of course also
true for the classical variables as seen in Hamilton’s equations in Eq. 36.

The commutator relation can be found using the value of the classical Poisson
bracket, which determines the value of the corresponding commutator up to a
factor of $i\hbar$, as Dirac argued [99]. Using this for the branch flux
operators and the charge operators, both defined in Eq. 2, we find that the
Poisson bracket is

$\\{\Phi_{b},Q_{b}\\}=\sum_{n}\left[\partialderivative{\Phi_{b}}{\phi_{n}}\partialderivative{Q_{b}}{q_{n}}-\partialderivative{Q_{b}}{\phi_{n}}\partialderivative{\Phi_{b}}{q_{n}}\right]=\pm
1,$ (46)

where the sign is plus for a capacitive branch and minus for an inductive
branch. Following Dirac’s approach, we arrive at the following commutator
relation

$[\hat{\Phi}_{b},\hat{Q}_{b}]=\pm i\hbar,$ (47)

which is equivalent to the commutator in Eq. 45. Note that in general, these
branch operators are not conjugate in the Hamiltonian. One must still find the
true degrees of freedom before quantization is applied.

###  IV.2 Effective energies

Consider the generalized momentum $\hat{\bm{q}}=\bm{C}\dot{\hat{\bm{\phi}}}$.
The time derivative of the generalized momentum is exactly the current through
the capacitors, $\hat{\bm{I}}=\bm{C}\ddot{\hat{\bm{\phi}}}$. Note that in the
limiting case of one node, this reduces to the current over a single parallel-
plate capacitor, as it should. For this reason, it makes sense to think of the
conjugate momentum as the sum of all charges on the capacitors attached to a
given node. We therefore define

$\hat{n}_{n}=-\frac{\hat{q}_{n}}{2e}$ (48)

as the net number of Cooper pairs stored on the $n$th node. If we consider the
kinetic energy of a circuit, we can write

$\hat{T}=\frac{1}{2}\hat{\bm{q}}^{T}\bm{C}^{-1}\hat{\bm{q}}=4\frac{e^{2}}{2}\hat{\bm{n}}^{T}\bm{C}^{-1}\hat{\bm{n}}.$
(49)

Now for each diagonal element, we have a contribution of
$4E_{C,n}\hat{n}_{n}^{2}$, where we define the effective capacitive energy of
the $n$th node as

$E_{C,n}=\frac{e^{2}}{2}(\bm{C}^{-1})_{(n,n)},$ (50)

which is equivalent to the energy required to store a single charge on the
capacitor. Note that in our dimensionless notation from Eq. 16 we have
$\hat{n}_{n}=-\hat{q}_{n}$, while the effective energy becomes
$E_{C,n}=(\bm{C}^{-1})_{(n,n)}/8$.

Similarly, we introduce the effective energies of the linear inductances and
Josephson junctions, $E_{L,n}$ and $E_{J,n}$, of each node. The effective
inductive energy is the diagonal elements of $\bm{L^{-1}}$, which is
equivalent to the sum of the inverse inductances of the inductors connected to
the given node. The effective Josephson energy is found as the sum of the
Josephson energies of the junctions connected to the given node.

Returning to our example circuit in Fig. 1, we can now write it using
operators and effective energies. It becomes

$\displaystyle\hat{\mathcal{H}}$
$\displaystyle=4\left(E_{C,1}\hat{n}_{1}^{2}+E_{C,2}\hat{n}_{2}^{2}\right)+E_{L,12}(\hat{\phi}_{1}-\hat{\phi}_{2}+\tilde{\Phi})^{2}$
(51)
$\displaystyle\phantom{=}-E_{J,1}\cos\hat{\phi}_{1}-E_{J,2}\cos\hat{\phi}_{2},$

where the coupling energy of the linear inductor is $E_{L,12}=1/2L_{12}$. The
effective energies of the Josephson junctions is the Josephson energies. Note
that since our example does not include any coupling capacitors, we do not
obtain any coupling term $\bm{(}\bm{C}^{-1})_{(1,2)}\hat{n}_{1}\hat{n}_{2}$
since $\bm{C}$ is diagonal. In reality, this is rarely the case.

##  V Recasting to interacting harmonic oscillators

We want to consider the low-energy limit of the superconducting circuit since
we want to create a qubit using the two lowest-lying states of the nonlinear
oscillator quantum system. This can be done by suppressing the kinetic energy
of the system, such that the ’position’ coordinate will be localized near the
minimum of the potential. We consider a single anharmonic oscillator (AHO) as
in Fig. 6 but with a possible linear inductor in parallel, which means that we
can omit subscripts in this section as there is only a single mode. The
Hamiltonian we thus consider is

$\hat{\mathcal{H}}_{\text{AHO}}=4E_{C}\hat{n}^{2}+E_{L}\hat{\phi}^{2}-E_{J}\cos\hat{\phi}.$
(52)

If the effective capacitive energy, $E_{C}$, of the mode is much smaller than
the effective Josephson energy, $E_{J}$, the flux will be well localized near
the bottom of the potential. This is equivalent to a heavy particle moving
near its equilibrium position. In this case, we can Taylor expand the
potential part of the Hamiltonian up to fourth order in $\phi$ such that the
Josephson-junction term takes the form

$E_{J}\cos\phi=E_{J}-\frac{1}{2}E_{J}\phi^{2}+\frac{1}{24}E_{J}\phi^{4}+O(\phi^{6}).$
(53)

Throwing away the irrelevant constant term, we are left with a Hamiltonian
consisting of second- and fourth-order terms. If we require the couplings
between different parts of the superconducting circuit to be small, we can
treat each mode individually as a harmonic oscillator perturbed by a quartic
anharmonicity and possibly some couplings to other modes of the system. For
each mode in our system, we have a simple harmonic oscillator (SHO) of the
form

$\hat{\mathcal{H}}_{\text{SHO}}=4E_{C}\hat{n}^{2}+\left(E_{L}+\frac{1}{2}E_{J}\right)\hat{\phi}^{2}.$
(54)

The simple harmonic oscillator is well understood quantum mechanically, and
using the algebraic approach [100] we define the annihilation and creation
operators

$\displaystyle\hat{b}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\frac{1}{\sqrt{\zeta}}\hat{\phi}-i\sqrt{\zeta}\hat{n}\right),$
(55a) $\displaystyle\hat{b}^{\dagger}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\frac{1}{\sqrt{\zeta}}\hat{\phi}+i\sqrt{\zeta}\hat{n}\right),$
(55b)

where we define the impedance

$\zeta=\sqrt{\frac{4E_{C}}{E_{L}+E_{J}/2}}.$ (56)

When restoring dimensions and going away from the dimensionless notation
defined in Eq. 16 the impedance in Eq. 56 must be multiplied with a factor of
$R_{Q}/2\pi$, where
$R_{Q}=h/(2e)^{2}\simeq$6.45\text{\,}\mathrm{k\SIUnitSymbolOhm}$$ is the
resistance quantum, which emerges in the quantum Hall effect. The annihilation
and creation operators fulfill the usual commutator relation

$[\hat{b},\hat{b}^{\dagger}]=1.$ (57)

Expressing the flux and conjugate momentum operators in terms of the
annihilation and creation operators,

$\displaystyle\hat{\phi}$
$\displaystyle=\sqrt{\frac{\zeta}{2}}(\hat{b}+\hat{b}^{\dagger}),$ (58a)
$\displaystyle\hat{n}$
$\displaystyle=\frac{i}{\sqrt{2\zeta}}(\hat{b}-\hat{b}^{\dagger}),$ (58b)

we can rewrite the oscillator part of the Hamiltonian as

$\hat{\mathcal{H}}_{\text{SHO}}=4\sqrt{E_{C}\left(E_{L}+\frac{1}{2}E_{J}\right)}\left(\hat{N}+\frac{1}{2}\right),$
(59)

where we introduce the usual number operator
$\hat{N}=\hat{b}^{\dagger}\hat{b}$.

Using the creation and annihilation operators, we can rewrite all quadratic
and quartic interaction terms. The results are in Table 2 for the most
commonly occurring terms.

Table 2: Overview of the different components and the operators they map to. Subscripts are included where appropriate and refer to different nodes. All constant terms are neglected. The impedance factor can be found in Eq. 56. Component | Hamiltonian term |  | Annihilation and
---
creation operators
All terms | $\hat{n}^{2}+\hat{\phi}^{2}$ | $4\sqrt{E_{C}\left(E_{L}+\frac{1}{2}E_{J}\right)}\hat{b}^{\dagger}\hat{b}$
Linear capacitors | $\hat{n}$ | $\frac{i}{\sqrt{2\zeta}}(\hat{b}-\hat{b}^{\dagger})$
| $\hat{n}^{2}$ | $-\frac{1}{2\zeta}(\hat{b}-\hat{b}^{\dagger})^{2}$
| $\hat{n}_{i}\hat{n}_{j}$ | $-\frac{1}{2\sqrt{\zeta_{i}\zeta_{j}}}(\hat{b}_{i}-\hat{b}_{i}^{\dagger})(\hat{b}_{j}-\hat{b}_{j}^{\dagger})$
Linear inductors | $\hat{\phi}$ | $\sqrt{\frac{\zeta}{2}}(\hat{b}^{\dagger}+\hat{b})$
| $\hat{\phi}^{2}$ | $\frac{\zeta}{2}(\hat{b}^{\dagger}+\hat{b})^{2}$
| $\hat{\phi}_{i}\hat{\phi}_{j}$ | $\frac{\sqrt{\zeta_{i}\zeta_{j}}}{2}(\hat{b}_{i}^{\dagger}+\hat{b}_{i})(\hat{b}_{j}^{\dagger}+\hat{b}_{j})$
Josephson junctions | $\hat{\phi}^{4}$ | $\frac{\zeta^{2}}{4}(\hat{b}^{\dagger}+\hat{b})^{4}$
| $\hat{\phi}^{3}$ | $\frac{\zeta^{3/2}}{2^{3/2}}(\hat{b}^{\dagger}+\hat{b})^{3}$
| $\hat{\phi}_{i}^{3}\hat{\phi}_{j}$ | $\frac{\zeta_{i}^{3/2}\zeta_{j}^{1/2}}{4}(\hat{b}_{i}^{\dagger}+\hat{b}_{i})^{3}(\hat{b}_{j}^{\dagger}+\hat{b}_{j})$
| $\hat{\phi}_{i}^{2}\hat{\phi}_{j}^{2}$ | $\frac{\zeta_{i}\zeta_{j}}{4}(\hat{b}_{i}^{\dagger}+\hat{b}_{i})^{2}(\hat{b}_{j}^{\dagger}+\hat{b}_{j})^{2}$

Returning to our example circuit in Fig. 1, we can write the Hamiltonian in
Eq. 51 using annihilation and creation operators as

$\displaystyle\hat{\mathcal{H}}$
$\displaystyle=\omega_{1}\hat{b}_{1}^{\dagger}\hat{b}_{1}+\omega_{2}\hat{b}_{2}^{\dagger}\hat{b}_{2}+\frac{\alpha_{1}}{12}(\hat{b}_{1}^{\dagger}+\hat{b}_{1})^{4}$
(60)
$\displaystyle\phantom{=}+\frac{\alpha_{2}}{12}(\hat{b}_{2}^{\dagger}+\hat{b}_{2})^{4}+g_{12}(\hat{b}_{1}^{\dagger}+\hat{b}_{1})(\hat{b}_{2}^{\dagger}+\hat{b}_{2})$
$\displaystyle\phantom{=}+\chi_{1}(\hat{b}_{1}+\hat{b}_{1}^{\dagger})-\chi_{2}(\hat{b}_{2}+\hat{b}_{2}^{\dagger}),$

where we omit all constant terms. We further define

$\displaystyle\omega_{n}$
$\displaystyle=4\sqrt{E_{C,n}\left(E_{L,12}+\frac{1}{2}E_{J,n}\right)},$ (61a)
$\displaystyle\alpha_{n}$ $\displaystyle=-\frac{\zeta_{n}^{2}}{8}E_{J,n},$
(61b) $\displaystyle g_{12}$
$\displaystyle=-\sqrt{\zeta_{1}\zeta_{2}}E_{L,12},$ (61c)
$\displaystyle\chi_{n}$ $\displaystyle=\sqrt{2\zeta_{n}}E_{L,12}\tilde{\Phi},$
(61d) $\displaystyle\zeta_{n}$
$\displaystyle=\sqrt{\frac{4E_{C,n}}{E_{L,12}+E_{J,n}/2}},$ (61e)

where we refer to $\omega_{n}$ as the frequency, $\alpha_{n}$ as the
anharmonicity, and $g_{12}$ the oscillator coupling strength. Note that if the
effective inductive energy is zero, $E_{L,n}=0$, then the anharmonicity in Eq.
61b becomes $\alpha_{n}=-E_{C,n}$, which is often the case.

Note that in the presence of an external flux, one should be careful in
identifying the minimum of the potential around which one can then perform the
expansion, as in Eq. 53.

##  VI Time-averaged dynamics

When analyzing the Hamiltonian of the circuit, it is often advantageous to
consider which terms dominate the time evolution and which terms only give
rise to minor corrections. The latter can often be neglected without changing
the overall behavior of the system. It can often be difficult to determine
which terms dominate, as different scales influence the dynamics of the
system. This stems from the fact that the frequencies, $\omega_{n}$, of the
oscillators are usually of the order $\mathrm{GHz}$ while the interactions
between the different oscillators are usually much smaller, on the order of
$\mathrm{MHz}$. We therefore employ separation of scales to remove the large
energy differences of the modes from the Hamiltonian. This makes it possible
to see the details of the interactions. In order to do this, we first
introduce the concept of the interaction picture, where the interacting part
of the Hamiltonian is in focus.

To summarize which terms we consider in the Hamiltonian, we divide the terms
into three categories.

  * •

Large trivial terms: Well understood energy difference terms, such as the
qubit frequencies, which we remove using separation of scales by transforming
into the interaction picture. Usually of the order $\mathrm{GHz}$.

  * •

Smaller but interesting terms: The dominant part of the interaction we are
interested in. Usually of the order $\mathrm{MHz}$

  * •

Small negligible terms: The suppressed part of the interaction, which does not
contribute significantly to the time evolution. These can be removed using the
rotating-wave approximation (RWA).

Note, however, that the above categorization is only a guide, and one should
always consider each term in relation to the concrete system at hand.

###  VI.1 Interaction picture

Consider the state $\ket{\psi,t}_{S}$ at time $t$. This state satisfies the
Schrödinger equation,

$i\partialderivative{t}\ket{\psi,t}_{S}=\hat{\mathcal{H}}\ket{\psi,t}_{S},$
(62)

where $\hat{\mathcal{H}}$ is the Hamiltonian. The subscript $S$ refers to the
Schrödinger picture. In the Schrödinger picture operators are time independent
and states are time dependent. We wish to change into the interaction picture
by splitting the Hamiltonian in a way such that the dynamics are separated
from the noninteracting part,
$\hat{\mathcal{H}}=\hat{\mathcal{H}}_{0}+\hat{\mathcal{H}}_{I,S}$. There are
often several ways to make this splitting depending on what interaction we
want to highlight. This separation comes at the cost that both the operators
and states become time dependent. The advantage of using a specific splitting
of the full Hamiltonian is that we can highlight some desired physics while
ignoring other parts that are well understood. This is analogous to choosing a
reference frame rotating with the Earth when doing classical physics in a
reference frame fixed on the surface of the Earth.

States in the interaction picture are defined as

$|\psi,t\rangle_{I}=e^{i\hat{\mathcal{H}}_{0}t}|\psi,t\rangle_{S},$ (63)

where the subscript $I$ refers to the interaction picture. The operators in
the interaction picture are defined as

$\hat{\mathcal{O}}_{I}=e^{i\hat{\mathcal{H}}_{0}t}\hat{\mathcal{O}}_{S}e^{-i\hat{\mathcal{H}}_{0}t},$
(64)

where $\hat{\mathcal{O}}_{S}$ is an operator in the Schrödinger picture.

It is then possible to show that the state satisfies the following Schrödinger
equation

$i\partialderivative{t}|\psi,t\rangle_{I}=\hat{\mathcal{H}}_{I}|\psi,t\rangle_{I},$
(65)

where
$\hat{\mathcal{H}}_{I}=e^{i\hat{\mathcal{H}}_{0}t}\hat{\mathcal{H}}_{I,S}e^{-i\hat{\mathcal{H}}_{0}t}$
is the interaction part of the Hamiltonian in the interaction picture.

In general, one can transform a Hamiltonian to any so-called rotating frame
using the transformation rule

$\hat{\mathcal{H}}\rightarrow\hat{\mathcal{H}}_{R}=\hat{\mathcal{U}}(t)^{\dagger}\hat{\mathcal{H}}\,\hat{\mathcal{U}}(t)+i\derivative{\hat{\mathcal{U}}(t)^{\dagger}}{t}\hat{\mathcal{U}}(t),$
(66)

where $\hat{\mathcal{U}}(t)$ is a unitary transformation. This transformation
rule holds for any unitary transformation and is quite useful to keep in mind.
Note that Eq. 66 is equivalent to transforming the Hamiltonian into the
interaction picture $\hat{\mathcal{H}}\rightarrow\hat{\mathcal{H}}_{I}$ when
$\hat{\mathcal{U}}(t)=\exp(-i\hat{\mathcal{H}}_{0}t)$, and
$\hat{\mathcal{H}}_{0}$ is the noninteracting part of the Hamiltonian, as the
second term removes the noninteracting part of the Hamiltonian.

One can also show that the time evolution of the operators in the interaction
picture is governed by a Heisenberg equation of motion

$\derivative{t}\hat{\mathcal{O}}_{I}=i[\hat{\mathcal{H}}_{0},\hat{\mathcal{O}}_{I}],$
(67)

where we assume no explicit time dependence in $\hat{\mathcal{O}}_{S}$. Note
that this implies that the voltage of the $b$th branch can be calculated as
$\hat{V}_{b}=i[\hat{\mathcal{H}},\hat{\Phi}_{b}]$. For more information about
the interaction picture see e.g., Ref. [100].

###  VI.2 Rotating-wave approximation

Consider now the weakly anharmonic oscillator as seen in Fig. 9(c), which has
the quantized Hamiltonian

$\hat{\mathcal{H}}=\omega\hat{b}^{\dagger}\hat{b}+\frac{\alpha}{12}(\hat{b}^{\dagger}+\hat{b})^{4},$
(68)

where we remove all constant terms. The frequency is
$\omega=\sqrt{8E_{C}E_{J}}$ and the anharmonicity is $\alpha=-E_{C}$ where
$E_{C}$ and $E_{J}$ are the effective capacitive energy and Josephson energy,
respectively. Now we choose the first term as the noninteracting Hamiltonian,
$\hat{\mathcal{H}}_{0}=\omega\hat{b}^{\dagger}\hat{b}$. We want to figure out
how the annihilation and creation operators behave in the interaction picture,
i.e., we want to calculate Eq. 64 for the annihilation and creation operators.
First, we notice that
$\hat{\mathcal{H}}_{0}^{n}\hat{b}^{\dagger}=\hat{b}^{\dagger}(\hat{\mathcal{H}}_{0}+\omega)^{n}$.
Using this and expanding the exponential functions, we can prove that

$e^{i\hat{\mathcal{H}}_{0}t}\hat{b}^{\dagger}e^{-i\hat{\mathcal{H}}_{0}t}=\hat{b}^{\dagger}e^{i\omega
t}.$ (69)

By taking the complex conjugate, we find a similar expression for $\hat{b}$,
but with a minus in the exponential factor on the right-hand side.

We now wish to consider how different combinations of the annihilation and
creation operators transform in the interaction picture. Starting with the
number operator $\hat{N}=\hat{b}^{\dagger}\hat{b}$, we see the exponential
factor from $\hat{b}^{\dagger}$ cancels the exponential factor from $\hat{b}$,
meaning that the number operator is unaffected by the transformation. This is
not surprising as the noninteracting Hamiltonian is chosen exactly as the
number operator. However, if we consider terms like
$J\hat{b}^{\dagger}\hat{b}^{\dagger}$ we find that in the interaction picture
they take the form $J\hat{b}^{\dagger}\hat{b}^{\dagger}e^{2i\omega t}$. If
$\omega$ is sufficiently large compared to the factor, $J$, in front of the
term (which is often the case in superconducting circuit Hamiltonians, where
$\omega\sim$\mathrm{GHz}$$, while other terms are usually of the order
$J\sim$\mathrm{MHz}$$), these terms will oscillate very rapidly on the
timescale induced by $J$. The time average over such terms on a timescale of
$\tau\sim 1/J$ is zero, and we can therefore neglect them as they only give
rise to minor corrections. This is the rotating-wave approximation which is
widely used in atomic physics [101, 102]. The story is the same for
$\hat{b}\hat{b}$ terms. All terms that do not conserve the number of
excitations (or quanta) of the system, i.e., terms where the number of
annihilation operators is not equal to the number of creation operators, will
rotate rapidly and can therefore safely be neglected. Note that while these
individual terms are nonconserving, they always appear in conjugate pairs in
the Hamiltonian such that the full Hamiltonian conserves the excitations, as
it should.

It is important to point out that despite the naming, the ’conservation’ is
not related to a conservation law resulting from a symmetry, i.e., like in
Noether’s theorem. Rather, the statement here means that the excitation
conserving terms are much more important than the nonconserving terms as long
as the conditions for using the rotating-wave approximation are satisfied.

Now consider the anharmonicity term of Eq. 68. When only including excitation
conserving terms and removing irrelevant constants, the anharmonicity term
takes the form

$\displaystyle\frac{\alpha}{12}(\hat{b}^{\dagger}+\hat{b})^{4}$
$\displaystyle=\alpha\left(\frac{1}{2}\hat{b}^{\dagger}\hat{b}^{\dagger}\hat{b}\hat{b}+\hat{b}^{\dagger}\hat{b}\right)$
(70) $\displaystyle\phantom{=}+\text{Nonconserving terms}.$

The last term, $\hat{b}^{\dagger}\hat{b}$, is the number operator, and we can
therefore consider it a correction to the frequency, such that the dressed
frequency becomes $\tilde{\omega}=\omega+\alpha=\sqrt{8E_{C}E_{J}}-E_{C}$. The
remaining $(\hat{b}^{\dagger}\hat{b}^{\dagger}\hat{b}\hat{b})$ term makes the
oscillator anharmonic. For this reason, we call $\alpha$ the anharmonicity of
the anharmonic oscillator. If we remove terms that do not conserve the number
of excitation, the Hamiltonian takes the form (in the Schrödinger picture)

$\hat{\mathcal{H}}=\tilde{\omega}\hat{b}^{\dagger}\hat{b}+\frac{\alpha}{2}\hat{b}^{\dagger}\hat{b}^{\dagger}\hat{b}\hat{b}.$
(71)

Next, consider an interaction term like the one in Eq. 60

$(\hat{b}_{i}^{\dagger}+\hat{b}_{i}^{)}(\hat{b}_{j}^{\dagger}+\hat{b}_{j})=\hat{b}_{i}^{\dagger}\hat{b}_{j}+\hat{b}_{i}\hat{b}_{j}^{\dagger}+\hat{b}_{i}^{\dagger}\hat{b}_{j}^{\dagger}+\hat{b}_{i}\hat{b}_{j}.$
(72)

Changing into the interaction picture, we realize that the two last terms
obtain a phase of $\exp{[\pm i(\omega_{i}+\omega_{j})t]}$, which can be
considered a fast oscillating term if the frequencies $\omega_{i}+\omega_{j}$
are much larger than the interaction strength, which is usually the case. We
can therefore safely neglect these nonconserving terms. The two first terms on
the other hand obtain a phase of $\exp({\pm i\delta t})$, where
$\delta=\omega_{i}-\omega_{j}$ is called the detuning of the two oscillators.
It is therefore tempting to say that these terms only contribute if
$\omega_{i}\approx\omega_{j}$. This is, however, not the whole story. More
precisely, we find that

$\displaystyle\hat{b}^{\dagger}_{i}\hat{b}_{j}e^{i\delta
t}+\hat{b}_{i}\hat{b}_{j}^{\dagger}e^{-i\delta t}$
$\displaystyle=(\hat{b}^{\dagger}_{i}\hat{b}_{j}+\hat{b}_{i}\hat{b}_{j}^{\dagger})\cos\delta
t$ (73)
$\displaystyle+i(\hat{b}^{\dagger}_{i}\hat{b}_{j}-\hat{b}_{i}\hat{b}_{j}^{\dagger})\sin\delta
t,$

which can be useful in some situations, e.g., when driving qubits, see Section
VIII. However, as a general rule of thumb, one can neglect these terms unless
$\omega_{i}\simeq\omega_{j}$, i.e., $\delta=0$. For a more general discussion
on the validity of the time averaging dynamics see Ref. [103].

If we consider the example circuit in Fig. 1, under the assumption that
$|\alpha_{n}|\ll\omega_{n}$, we can time average its Hamiltonian in Eq. 60. We
choose the noninteracting Hamiltonian as
$\mathcal{H}_{0}=\omega_{1}\hat{b}_{1}^{\dagger}\hat{b}_{1}+\omega_{2}\hat{b}_{2}^{\dagger}\hat{b}_{2}$,
which means that the interacting part of the Hamiltonian in Eq. 60 becomes

$\displaystyle\hat{\mathcal{H}}_{I}$
$\displaystyle=\frac{\alpha_{1}}{2}(\hat{b}_{1}^{\dagger}\hat{b}_{1}^{\dagger}\hat{b}_{1}\hat{b}_{1}+2\hat{b}_{1}^{\dagger}\hat{b}_{1})+\frac{\alpha_{2}}{2}(\hat{b}_{2}^{\dagger}\hat{b}_{2}^{\dagger}\hat{b}_{2}\hat{b}_{2}+2\hat{b}_{2}^{\dagger}\hat{b}_{2})$
(74)
$\displaystyle\phantom{=}+g_{12}\left(\hat{b}_{1}^{\dagger}\hat{b}_{2}e^{i\delta
t}+\hat{b}_{1}\hat{b}_{2}^{\dagger}e^{-i\delta t}\right.$
$\displaystyle\phantom{=+g_{12}}+\left.\hat{b}_{1}\hat{b}_{2}e^{-i(\omega_{1}+\omega_{2})t}+\hat{b}_{1}^{\dagger}\hat{b}_{2}^{\dagger}e^{i(\omega_{1}+\omega_{2})t}\right)$
$\displaystyle\phantom{=}-\sum_{n=1}^{2}(-1)^{n}\chi_{n}(\hat{b}_{n}e^{-i\omega_{n}t}+\hat{b}_{n}^{\dagger}e^{i\omega_{n}t}),$

where we define the detuning $\delta=\omega_{1}-\omega_{2}$. Assuming
$\omega_{1}\simeq\omega_{2}$, i.e., $\delta\simeq 0$, and if we further assume
that $\omega_{n}\gg g_{12}$, then the coupling terms $\hat{b}_{1}\hat{b}_{2}$
and $\hat{b}_{1}^{\dagger}\hat{b}_{2}^{\dagger}$ are fast oscillating and can
thus be neglected.

We can also write the Hamiltonian in the Schrödinger picture, removing terms
that do not conserve excitations. This yields

$\displaystyle\hat{\mathcal{H}}$
$\displaystyle=\tilde{\omega}_{1}\hat{b}_{1}^{\dagger}\hat{b}_{1}+\tilde{\omega}_{2}\hat{b}_{2}^{\dagger}\hat{b}_{2}+\frac{\alpha_{1}}{2}\hat{b}_{1}^{\dagger}\hat{b}_{1}^{\dagger}\hat{b}_{1}\hat{b}_{1}$
(75)
$\displaystyle\phantom{=}+\frac{\alpha_{2}}{2}\hat{b}_{2}^{\dagger}\hat{b}_{2}^{\dagger}\hat{b}_{2}\hat{b}_{2}+g_{12}(\hat{b}_{1}^{\dagger}\hat{b}_{2}+\hat{b}_{1}\hat{b}_{2}^{\dagger}),$

where we introduce the revised frequencies
$\tilde{\omega}_{n}=\omega_{n}+\alpha_{n}$. Writing the Hamiltonian in this
frame without nonconserving terms reveals the effect of the anharmonicity. In
Eq. 75 we also assume that $\omega_{n}\gg\chi_{n}$, meaning that all terms
related to the external flux are neglected. However, since $\chi_{n}$ depends
on $\tilde{\Phi}$, which can be controlled externally, it is possible to tune
$\chi_{n}$ such that the terms involving $\chi_{n}$ are not suppressed. This
can be used to drive the modes, i.e., to add excitations to the two degrees of
freedom.

Figure 9: (a) Circuit of an $LC$ oscillator with inductance $L$ and
capacitance $C$. We denote the phase on the superconducting island $\phi$,
while the ground node has phase zero. (b) Energy potential of a quantum
harmonic oscillator, as can be obtained by an LC-circuit. Here the energy
levels are equidistantly spaced $\hbar\omega$ apart, where
$\omega=\sqrt{1/LC}$. (c) Josephson junction qubit circuit, where the linear
inductor is replaced by a nonlinear Josephson junction of energy $E_{J}$. (d)
The Josephson junction changes the harmonic potential (blue dashed) into a
sinusoidal potential (orange solid), yielding nonequidistant energy levels.

##  VII Truncation

A harmonic oscillator, as one gets from a regular LC-circuit, has a spectrum
consisting of an infinite number of equally spaced energy eigenstates [see
Fig. 9(b)]. This is not desirable, as we wish to consider only the lowest
states of the system in order to realize a qubit. However, when we introduce a
Josephson junction instead of a linear inductor, we introduce an
anharmonicity, compare Fig. 9(a) and (c). The anharmonicity stems from the
$(\hat{b}^{\dagger}\hat{b}^{\dagger}\hat{b}\hat{b})$ terms [see Eq. 71] and
can be viewed as perturbations to the harmonic oscillator Hamiltonian if
$|\alpha|\ll\omega$. This anharmonicity changes the spacing between the energy
levels of the harmonic oscillator, making it an anharmonic oscillator [see
Fig. 9(d)]. Formally, the anharmonicity is defined as the difference between
the first and second energy gap, while we define the relative anharmonicity as
the anharmonicity divided by the first energy gap

$\alpha=E_{12}-E_{01},\qquad\alpha_{r}=\frac{\alpha}{E_{01}}.$ (76)

Note that this anharmonicity is the same anharmonicity factor in front of the
$(\hat{b}^{\dagger}\hat{b}^{\dagger}\hat{b}\hat{b})$ terms mentioned in
previous sections.

To operate only on the two lowest levels of the oscillator, the anharmonicity
must be larger than the bandwidth of operations on the qubit. That is, if we
want to drive excitation between the two lowest levels of the anharmonic
oscillator, the anharmonicity must be larger than the amplitude of the driving
field (also known as the Rabi frequency, see Section VIII). If the
anharmonicity is smaller than the amplitude of the driving field, we cannot
distinguish between the energy gaps of the oscillator, and we end up driving
multiple transitions in the spectrum instead of just the lowest one.

Taking this into account we find that as a rule-of-thumb, the relative
anharmonicity should be at least a couple of percent for the system to make an
effective qubit. In actual numbers, this converts to an anharmonicity around
100-$300\text{\,}\mathrm{MHz}$ for a qubit frequency around
3-$6\text{\,}\mathrm{GHz}$ [104, 11]. It does not matter whether the
anharmonicity is positive or negative. For transmon-type qubits, it will be
negative, while it can be either positive or negative for flux-type qubits.
The relative anharmonicity is proportional to $\sqrt{E_{J}/E_{C}}$, which
means that this ratio must be of a certain size for the anharmonicity to have
an effect. This is in contrast to what was discussed at the beginning of
Section V, where we argued that we required this ratio to be as low as
possible to allow for the expansion of cosines. Thus we need to find a
suitable regime for the ratio, $E_{J}/E_{C}$. This regime is usually called
the transmon regime and is around 50-100.

In the following section we assume that we have a sufficiently large
anharmonicity to truncate the system into a two-level system. However, nothing
is stopping us from keeping more levels, as we do in Section VII.2.

As an alternative to the methods for truncation presented in this tutorial,
black-box quantization can be useful for determining the effective low-energy
spectrum of a weakly anharmonic Hamiltonian [105, 106, 107]. This approach is
especially useful when dealing with impedances in the circuit, but is beyond
the scope of this tutorial.

###  VII.1 Two-level model (qubit)

Table 3: Overview of the different combinations of the annihilation and creation operators and their truncation to two-dimensional Pauli operators. Subscripts are included for the interaction terms and refer to different nodes. All constant terms are ignored. | Annihilation and
---
creation operators
Pauli operators
$\hat{b}^{\dagger}-\hat{b}$ | $-i\sigma^{y}$
$\hat{b}^{\dagger}+\hat{b}$ | $\sigma^{x}$
$(\hat{b}^{\dagger}-\hat{b})^{2}$ | $-\sigma^{z}$
$(\hat{b}^{\dagger}+\hat{b})^{2}$ | $-\sigma^{z}$
$(\hat{b}^{\dagger}+\hat{b})^{3}$ | $3\sigma^{x}$
$(\hat{b}^{\dagger}+\hat{b})^{4}$ | $-6\sigma^{z}$
$(\hat{b}_{i}^{\dagger}-\hat{b}_{i})(\hat{b}_{j}^{\dagger}-\hat{b}_{j})$ | $-\sigma^{y}_{i}\sigma^{y}_{j}$
$(\hat{b}_{i}^{\dagger}+\hat{b}_{i})(\hat{b}_{j}^{\dagger}+\hat{b}_{j})$ | $\sigma^{x}_{i}\sigma^{x}_{j}$
$(\hat{b}_{i}^{\dagger}+\hat{b}_{i})^{3}(\hat{b}_{j}^{\dagger}+\hat{b}_{j})$ | $3\sigma^{x}_{i}\sigma^{x}_{j}$
$(\hat{b}_{i}^{\dagger}+\hat{b}_{i})^{2}(\hat{b}_{j}^{\dagger}+\hat{b}_{j})^{2}$ | $\sigma^{z}_{i}\sigma^{z}_{j}-2\sigma^{z}_{i}-2\sigma^{z}_{j}$

In a two-level system, which is equivalent to a qubit, we can represent the
state of the system with two-dimensional vectors

$|0\rangle\sim\begin{bmatrix}1\\\
0\end{bmatrix},\qquad|1\rangle\sim\begin{bmatrix}0\\\ 1\end{bmatrix}.$ (77)

In this reduced Hilbert space all operators can be expressed by the Pauli
matrices,

$\sigma^{x}=\begin{bmatrix}0&1\\\
1&0\end{bmatrix},\quad\sigma^{y}=\begin{bmatrix}0&-i\\\
i&0\end{bmatrix},\quad\sigma^{z}=\begin{bmatrix}1&0\\\ 0&-1\end{bmatrix},$
(78)

and the identity, since these four matrices span all $2\times 2$ Hermitian
matrices. If we view the unitary operations as rotations in the Hilbert space,
we can parameterize the superposition of the two states using a complex phase,
$\phi$, and a mixing angle, $\theta$

$|\psi\rangle=\alpha\ket{0}+\beta\ket{1}=\cos\left(\frac{\theta}{2}\right)|0\rangle+e^{i\phi}\sin\left(\frac{\theta}{2}\right)|1\rangle,$
(79)

where $0\leq\theta\leq\pi$ and $0\leq\phi<2\pi$ and
$|\alpha|^{2}+|\beta|^{2}=1$. With this, we can illustrate the qubit as a unit
vector on the Bloch sphere, see Fig. 10. It is conventional to let the north
pole represent the $\ket{0}$ state, while the south pole represents the
$\ket{1}$ state. These lie on the $z$ axis, which is called the longitudinal
axis as it represents the quantization axis for the states in the qubit. The
$x$ and $y$ axes are called the transverse axes.

Solving the Schrödinger equation in Eq. 62 for the state in Eq. 79 shows that
it precesses around the $z$ axis at the qubit frequency. However, changing
into a frame rotating with the frequency of the qubit, following the approach
in Section VI.1. makes the Bloch vector stationary.

Unitary operations can be seen as rotations on the Bloch sphere and the Pauli
matrices are thus the generators of rotations. Linear operators will then be
represented by $2\times 2$ matrices as

$\displaystyle
M_{2}[\hat{O}]=\begin{bmatrix}\matrixelement{0}{\hat{\mathcal{O}}}{0}&\matrixelement{0}{\hat{\mathcal{O}}}{1}\\\
\matrixelement{1}{\hat{\mathcal{O}}}{0}&\matrixelement{1}{\hat{\mathcal{O}}}{1}\end{bmatrix}.$
(80)

In general we denote the $n\times n$ matrix representation of an operator
$\hat{\mathcal{O}}$ with $M_{n}[\hat{\mathcal{O}}]$.

Figure 10: The Bloch sphere. Each point on the Bloch sphere corresponds to a
quantum state. Rotations around the sphere correspond to transformations of
the state.

In order to apply this mapping to the Hamiltonian, we must map each operator
in each term. As an example, we truncate the $(\hat{b}^{\dagger}+\hat{b})^{3}$
term from Table 3:

$\displaystyle(\hat{b}^{\dagger}+\hat{b})^{3}|0\rangle$
$\displaystyle=(\hat{b}^{\dagger}+\hat{b})^{2}|1\rangle$
$\displaystyle=(\hat{b}^{\dagger}+\hat{b})\left(\sqrt{2}|2\rangle+|0\rangle\right)$
$\displaystyle=\sqrt{6}|3\rangle+3|1\rangle,$
$\displaystyle(\hat{b}^{\dagger}+\hat{b})^{3}|1\rangle$
$\displaystyle=(\hat{b}^{\dagger}+\hat{b})^{2}\left(\sqrt{2}|2\rangle+|0\rangle\right)$
$\displaystyle=(\hat{b}^{\dagger}+\hat{b})\left(\sqrt{6}|3\rangle+3|1\rangle\right)$
$\displaystyle=\sqrt{24}|4\rangle+6\sqrt{2}|2\rangle+3|0\rangle.$

Using the orthonormality of the states we obtain the representation of the
operator

$\displaystyle M_{2}[(\hat{b}^{\dagger}+\hat{b})^{3}]$
$\displaystyle=\begin{bmatrix}\langle
0|(\hat{b}^{\dagger}+\hat{b})^{3}|0\rangle&\langle
0|(\hat{b}^{\dagger}+\hat{b})^{3}|1\rangle\\\ \langle
1|(\hat{b}^{\dagger}+\hat{b})^{3}|0\rangle&\langle
1|(\hat{b}^{\dagger}+\hat{b})^{3}|1\rangle\end{bmatrix}$
$\displaystyle=\begin{bmatrix}0&3\\\ 3&0\end{bmatrix}=3\sigma^{x}.$

Truncation of the remaining terms is presented in Table 3.

If we consider the example circuit in Fig. 1, after we remove nonconserving
terms as in Eq. 75 and assume an anharmonicity large enough for truncation to
a two-level system, we obtain the following Hamiltonian:

$\hat{\mathcal{H}}=-\frac{\tilde{\omega}_{1}}{2}\sigma_{1}^{z}-\frac{\tilde{\omega}_{2}}{2}\sigma_{2}^{z}+g_{12}(\sigma_{1}^{+}\sigma_{2}^{-}+\sigma_{1}^{-}\sigma_{2}^{+}),$
(81)

where we define $\sigma^{\pm}_{n}=(\sigma^{x}_{n}\mp i\sigma^{y}_{n})/2$. This
Hamiltonian represents two qubits that can interact by swapping excitation
between them, i.e., interacting via a swap coupling.

###  VII.2 Three-level model (qutrit)

It can be desirable to truncate to the three lowest levels of the anharmonic
oscillator, i.e., the three lowest states of Fig. 9(d). This can, e.g., be
useful if one wants to study qutrit systems [108, 109, 56], or the leakage
from the qubit states to higher states [110, 111]. In this case, the operators
will be represented as $3\times 3$ matrices. The matrix representation of the
annihilation and creation operators become

$M_{3}[\hat{b}^{\dagger}]=\begin{bmatrix}0&0&0\\\ 1&0&0\\\
0&\sqrt{2}&0\end{bmatrix},\quad M_{3}[\hat{b}]=\begin{bmatrix}0&1&0\\\
0&0&\sqrt{2}\\\ 0&0&0\end{bmatrix},$ (82a) while the number operator is
$M_{3}[\hat{b}^{\dagger}\hat{b}]=\begin{bmatrix}0&0&0\\\ 0&1&0\\\
0&0&2\end{bmatrix},$ (82b) and powers of $\hat{b}^{\dagger}+\hat{b}$ become
$\displaystyle M_{3}[(\hat{b}^{\dagger}+\hat{b})^{2}]$
$\displaystyle=\begin{bmatrix}1&0&\sqrt{2}\\\ 0&3&0\\\
\sqrt{2}&0&5\end{bmatrix},$ (82c) $\displaystyle
M_{3}[(\hat{b}^{\dagger}+\hat{b})^{3}]$ $\displaystyle=\begin{bmatrix}0&3&0\\\
3&0&6\sqrt{2}\\\ 0&6\sqrt{2}&0\end{bmatrix},$ (82d) $\displaystyle
M_{3}[(\hat{b}^{\dagger}+\hat{b})^{4}]=$
$\displaystyle\begin{bmatrix}3&0&6\sqrt{2}\\\ 0&15&0\\\
6\sqrt{2}&0&39\end{bmatrix}.$ (82e)

From Eq. 82e it is clear to see the varying size of the anharmonicity, as the
differences $15-3=12$ and $39-15=24$ between the levels changes. This pattern
continues for higher levels and means that we can distinguish between all the
levels in principle.

As we are dealing with $3\times 3$ matrices we can no longer use the Pauli
spin-1/2 matrices as a basis for the operators. In this case one can use the
Gell-Mann matrices as a basis. However, often it is more convenient to leave
the annihilation and creation operators as above. We are not limited by three
levels, and it is possible to truncate the system to an arbitrary number of
levels, thus creating a so-called qudit.

It is also possible to truncate the system before expanding the cosine
functions of the Josephson junctions. This approach is discussed in Appendix C
where we also truncate an anharmonic oscillator to the four lowest levels.

##  VIII Microwave driving

Figure 11: Circuit diagram of a single transmonlike superconducting qubit
capacitively coupled to a microwave drive line.

Single-qubit rotations in superconducting circuits can be achieved by
capacitive microwave driving. In this section, we go through the steps of
analyzing a microwave-controlled transmonlike qubit and then generalize to a
$d$-level qudit. To this end, we consider the superconducting qubit seen in
Fig. 11, which is capacitively coupled to a microwave source. Using the
approach presented in Section III.2 the Lagrangian of this circuit becomes

$\mathcal{L}=\frac{C}{2}\dot{\phi}^{2}+E_{J}\cos\phi+\frac{C_{\text{ext}}}{2}\left(V(t)-\dot{\phi}\right)^{2},$
(83)

where $\phi$ is the node flux. Expanding the last term, we obtain

$\mathcal{L}=\mathcal{L}_{0}+\frac{C_{\text{ext}}}{2}\left(V(t)^{2}+\dot{\phi}^{2}-2V(t)\dot{\phi}\right),$
(84)

where $\mathcal{L}_{0}$ is the static part of the Lagrangian, i.e., the two
first terms of Eq. 83. The first term in the parenthesis is an irrelevant
offset term, the second term is a change of the capacitance of the node, while
the last term is our driving term. We throw away the offset term and rewrite

$\mathcal{L}=\frac{C+C_{\text{ext}}}{2}\dot{\phi}^{2}+E_{J}\cos\phi-
C_{\text{ext}}V(t)\dot{\phi}.$ (85)

The conjugate momentum of the node flux, $\phi$, is then

$q=(C+C_{\text{ext}})\dot{\phi}-C_{\text{ext}}V(t).$ (86)

Doing the usual Legendre transformation, our Hamiltonian takes the form

$\mathcal{H}=\underbrace{\frac{1}{2}\frac{1}{C+C_{\text{ext}}}q^{2}-E_{J}\cos\phi}_{\mathcal{H}_{\text{AHO}}}+\underbrace{\frac{C_{\text{ext}}}{C+C_{\text{ext}}}V(t)q}_{\mathcal{H}_{\text{ext}}},$
(87)

where we denote the anharmonic oscillator part of the Hamiltonian
$\mathcal{H}_{\text{AHO}}$ and the external driving part
$\mathcal{H}_{\text{ext}}$. We are now ready to perform the quantization and
the driving part becomes

$\hat{\mathcal{H}}_{\text{ext}}=\frac{i}{\sqrt{2\zeta}}\frac{C_{\text{ext}}}{C+C_{\text{ext}}}V(t)(\hat{b}^{\dagger}-\hat{b}).$
(88)

Assuming a large enough anharmonicity, we can truncate the Hamiltonian into
the two lowest levels

$\hat{\mathcal{H}}=-\frac{1}{2}\omega\sigma^{z}+\Omega V(t)\sigma^{y},$ (89)

where $\omega$ is the qubit frequency and
$\Omega=C_{\text{ext}}/[\sqrt{2\zeta}(C+C_{\text{ext}})]$ is the Rabi
frequency of the transition between the ground state and the excited state.
Note that the size of the Rabi frequency is limited by the size of the
anharmonicity, as discussed in Section VII. The name Rabi frequency may cause
a bit of confusion at first as it is not the frequency of the driving
microwave but rather the amplitude. However, the Rabi frequency is named so
since it is equal to the frequency of oscillation between the two states in a
qubit when the driving frequency, $\omega_{\textup{ext}}$, is equal to the
qubit frequency, $\omega$, i.e. when we drive the qubit ’on resonance’ [102].

We now change into a frame rotating with the frequency of the qubit, also
known as the interaction frame as discussed in Section VI. In particular we
use $\mathcal{H}_{0}=-\omega\sigma_{z}/2$ for the transformation in Eq. 66. In
this frame the Hamiltonian becomes

$\hat{\mathcal{H}}_{R}=\Omega V(t)\left(\cos(\omega t)\sigma^{y}-\sin(\omega
t)\sigma^{x}\right),$ (90)

which is equivalent to the external driving part of the Hamiltonian in the
interaction picture, i.e., $\mathcal{H}_{R}=\mathcal{H}_{\text{ext}}^{I}$. We
assume that the driving voltage is sinusoidal

$\displaystyle V(t)$
$\displaystyle=V_{0}\eta(t)\sin(\omega_{\text{ext}}t+\varphi)$ (91)
$\displaystyle=V_{0}\eta(t)\left[\cos(\varphi)\sin(\omega_{\text{ext}}t)+\sin(\varphi)\cos(\omega_{\text{ext}}t)\right],$

where $V_{0}$ is the amplitude of the voltage, $\eta(t)$ is a dimensionless
envelope function, $\omega_{\text{ext}}$ is the external driving frequency,
and $\varphi$ is the phase of the driving. One usually defines the in-phase
component $I=\cos(\varphi)$ and the out-of-phase component $Q=\sin(\varphi)$
[11]. Inserting the voltage in Eq. 91 into the Hamiltonian in Eq. 90 and
rewriting we obtain

$\displaystyle\hat{\mathcal{H}}_{R}=\frac{1}{2}\Omega V_{0}\eta(t)\big{\\{}-$
$\displaystyle\left[Q\sin(\delta t)+I\cos(\delta t)\right]\sigma^{x}$ (92)
$\displaystyle+$ $\displaystyle\left[Q\cos(\delta t)-I\sin(\delta
t)\right]\sigma^{y}\big{\\}},$

where $\delta=\omega-\omega_{\text{ext}}$ is the difference between the qubit
frequency and the driving frequency and we neglect fast oscillating terms,
i.e., terms with $\omega+\omega_{\text{ext}}$, following the rotating-wave
approximation. This Hamiltonian can be written very simple in matrix form

$\hat{\mathcal{H}}_{R}=-\frac{1}{2}\Omega
V_{0}\eta(t)\begin{bmatrix}0&e^{-i(\delta t-\varphi)}\\\ e^{i(\delta
t-\varphi)}&0\end{bmatrix}.$ (93)

From this, we conclude that if we apply a pulse at the qubit frequency, i.e.,
$\omega_{\text{ext}}=\omega$, we can rotate the state of the qubit around the
Bloch sphere in Fig. 10. By setting $\varphi=0$, i.e., using only the $I$
component we rotate about the $x$ axis. By setting $\varphi=\pi/2$, i.e.,
using only the $Q$ component, we rotate about the $y$ axis.

###  VIII.1 Single-qubit gates

One of the objectives of using superconducting circuits is to be able to
perform high-quality gate operations on qubit degrees of freedom [12].
Microwave driving of the qubits can be used to perform single-qubit rotation
gates. To see how this works we consider the unitary time-evolution operator
of the driving Hamiltonian. At qubit frequency, i.e., $\delta=0$, it takes the
form

$\displaystyle\hat{\mathcal{U}}(t)$
$\displaystyle=\exp\left[-i\int_{0}^{t}\hat{\mathcal{H}}_{R}(t^{\prime})\differential{t^{\prime}}\right]$
(94)
$\displaystyle=\exp\left[\frac{i}{2}\Theta(t)(I\sigma^{x}-Q\sigma^{y})\right],$

where we take the Pauli operators outside the integral as there is no time
dependence other than on the envelope $\eta(t)$. Note that this holds only for
$\delta=0$, as here the Hamiltonian commutes with itself at different times.
For nonzero $\delta$, one needs to solve the full Dyson’s series in principle
[100]. Equation 94 is known as Rabi driving and can be used for engineering
efficient single-qubit gate operations. The angle of rotation is defined as

$\Theta(t)=\Omega V_{0}\int_{0}^{t}\eta(t^{\prime})\differential{t^{\prime}},$
(95)

which depends on the macroscopic design parameters of the circuit, via the
coupling $\Omega$, the envelope of the pulse, $\eta(t)$, and the amplitude of
the pulse, $V_{0}$. The latter two can be controlled using arbitrary wave
generators (AWGs). In case one wishes to implement a $\pi$ pulse one must
adjust these parameters such that $\Theta(\tau)=\pi$, where $\tau$ is the
length of the driving pulse.

Consider a $\pi$ pulse. For the in-phase case, i.e., $\varphi=0$, the time-
evolution operator takes the form

$\hat{\mathcal{U}}_{\textsc{x}}(\tau)=\exp\left[\frac{i}{2}\pi\sigma^{x}\right]=\begin{bmatrix}0&1\\\
1&0\end{bmatrix},$ (96)

which is a Pauli-x gate, also known as a not-gate, which maps $\ket{0}$ to
$\ket{1}$ and vice versa [112, 113, 114, 115]. This corresponds to a rotation
by $\pi$ radians around the $x$ axis of the Bloch sphere. By changing the
value of $\Theta(\tau)$ it is possible to change the angle of the rotation.
Had we instead considered the out-of-phase case, i.e., $\varphi=\pi/2$ then we
would have obtained a Pauli-y gate which maps $\ket{0}$ to $i\ket{1}$ and
$\ket{1}$ to $-i\ket{0}$, corresponding to a rotation around the $y$ axis of
the Bloch sphere.

A Pauli-z gate can be implemented in one of three ways:

  * •

By detuning the qubit frequency with respect to the driving field for some
finite amount of time. This introduces an amplified phase error, which can be
modeled as effective qubit rotations around the $z$ axis [116].

  * •

Driving with an off-resonance microwave pulse. This introduces a temporary
Stark shift, which causes a phase change, corresponding to a rotation around
the $z$ axis.

  * •

Virtual z gates where a composition of x and y gates rotates the qubit state
around the $x$ and $y$ axes, which is equivalent to a rotation around the $z$
axis [117]. This can be achieved very effectively simply by adjusting the
phases of subsequent microwave gates [118].

Finally, we note that the Hadamard gate can be performed as a combination of
two rotations: a $\pi$ rotation around the $z$ axis and a $\pi/2$ rotation
around the $y$ axis.

###  VIII.2 Generalization to qudit driving

Now let us generalize the discussion to a $d$-dimensional qudit. Quantizing
and truncating the anharmonic oscillator part of the Hamiltonian in Eq. 87 to
$d$ levels, the qudit Hamiltonian becomes

$\hat{\mathcal{H}}_{\text{osc}}=\sum_{n=0}^{d-1}\omega_{n}\outerproduct{n}{n},$
(97)

where $\omega_{n}$ is the energy of qudit state $\ket{n}$. This is a rewriting
of the $\omega\hat{b}^{\dagger}\hat{b}$ term and the anharmonicity term, where
the anharmonicity has been absorbed into the set of $\omega_{n}$. Starting
from Eq. 88 and for simplicity setting the phase in Eq. 91 to $\pi/2$ such
that $V(t)=V_{0}\cos(\omega_{\text{ext}}t)$, we can move to the rotating frame
as was also done above for the qubit using Eq. 66. We choose the frame
rotating with the external driving frequency

$\displaystyle\mathcal{H}_{0}=\sum_{n=0}^{d-1}n\,\omega_{\text{ext}}\outerproduct{n}{n},$
(98)

which is contrary to what we did for the qubit, where we rotated into a frame
equal to the qubit frequency. We see that for a qubit ($d=2$) we get
$\mathcal{H}_{0}=-\omega_{\text{ext}}\sigma_{z}/2$ up to a global constant,
which we could have also chosen to use above, instead of the qubit frame.

Applying Eq. 66 to the qudit Hamiltonian in Eq. 97, we get

$\displaystyle\hat{\mathcal{H}}_{\text{osc},R}=\sum_{n=0}^{d-1}\left(\omega_{n}-n\omega_{\text{ext}}\right)\outerproduct{n}{n}.$
(99)

The same transformation is performed on $\mathcal{H}_{\text{ext}}$ by using
the standard expansion of the bosonic operators. By expanding the cosine in
the voltage drive using Euler’s formula, the total Hamiltonian in the rotating
frame can be found. It becomes

$\displaystyle\hat{\mathcal{H}}_{R}=\sum_{n=0}^{d-1}\delta_{n}\outerproduct{n}{n}+i\Omega_{n}\left(\outerproduct{n+1}{n}-\outerproduct{n}{n+1}\right),$
(100)

where $\delta_{n}=\omega_{n}-n\omega_{\text{ext}}$ is the detuning of the
$n$th state relative to the ground state driven by the external field and

$\displaystyle\Omega_{n}=\sqrt{n+1}\Omega=\sqrt{n+1}\frac{C_{\text{ext}}}{C+C_{\text{ext}}}\frac{V_{\text{0}}}{\sqrt{2\zeta}}$
(101)

is the Rabi frequency of the $n$th transition. Thus, by using a single drive,
we achieve great control over this specific qudit transition. Transitions
between other neighboring qudit states can be performed simultaneously by
using a multimode driving field. Note that the $i$ in the second term of Eq.
100 comes from the choice of $\varphi=\pi/2$, which can of course be changed
if desired.

The external field enables transitions between two states in the qudit if the
effective detuning, $\Delta_{n,n+1}$, is small compared to the size of the
Rabi frequency, $\Omega_{n}$. The effective detuning between the $n$th and
$(n+1)$th states is given as the difference between the detuning of the two
states:

$\Delta_{n,n+1}=\delta_{n+1}-\delta_{n}=\omega_{n+1}-\omega_{n}-\omega_{\text{ext}},$
(102)

from which we see that the frequency of the external field,
$\omega_{\text{ext}}$, has to match the energy difference between the two
states, $\omega_{n+1}-\omega_{n}$, for the driving to be efficient.

As mentioned in Section VII, leakage to other states when driving between two
states depends on the size of the anharmonicity. This can be understood from
Eq. 102. For a small anharmonicity, $\Delta_{n,n+1}$ is approximately the same
for all $n$ since $\omega_{n}$ will be approximately the same for all $n$,
thus it becomes difficult to single out the desired transition we want to
drive since the driving frequency, $\omega_{\text{ext}}$, will overlap with
multiple transition frequencies. Luckily, tailored control pulse methods such
as derivative removal by adiabatic gate (DRAG) and its improvements [114, 119]
can reduce this leakage significantly, which allows for relative
anharmonicities of just a couple of percent. The topic of tailored control
pulses is beyond the scope of this discussion, and we refer to the cited
works.

##  IX Coupling of modes

In our central example of Fig. 1, we considered direct inductive coupling.
While this coupling is rather straightforward theoretically it is rather
difficult to implement experimentally. We therefore now consider simpler ways
to couple qubits. By coupling qubits we also open up the possibility of
implementing two-qubit gates. Examples of more sophisticated approaches to
coupling qubits are discussed in Section XI.3.

###  IX.1 Capacitive coupling

Figure 12: (a) Two transmonlike qubits coupled by a single capacitor with
capacitance $C_{g}$, resulting in a static coupling between the modes. (b) A
transmonlike qubit coupled to a linear resonator via a capacitor of
capacitance $C_{g}$.

The simplest form of coupling both experimentally and theoretically is
arguably capacitive coupling. Consider two transmonlike qubits coupled by a
single capacitor with capacitance $C_{g}$, as seen in Fig. 12(a). Note the
similarities between this coupling and the circuit in Fig. 1. As we see, the
resulting Hamiltonian of Fig. 12(a) is close to the Hamiltonian in Eq. 75.
However, capacitive coupling are much simpler to achieve experimentally.

The Hamiltonian is easily found following the approach in Section III.2

$\mathcal{H}=\frac{1}{2}\bm{q}^{T}\bm{C}^{-1}\bm{q}-E_{J,1}\cos\phi_{1}-E_{J,2}\cos\phi_{2},$
(103)

where $\bm{q}=(q_{1},q_{2})^{T}$ is the vector of conjugate momentum and the
capacitance matrix is

$\bm{C}=\begin{bmatrix}C_{1}+C_{g}&-C_{g}\\\ -C_{g}&C_{2}+C_{g}\end{bmatrix},$
(104)

which is invertible

$\bm{C}^{-1}=\frac{1}{C_{\Sigma}}\begin{bmatrix}C_{2}+C_{g}&C_{g}\\\
C_{g}&C_{1}+C_{g}\end{bmatrix}\simeq\begin{bmatrix}\frac{1}{C_{1}}&\frac{C_{g}}{C_{1}C_{2}}\\\
\frac{C_{g}}{C_{1}C_{2}}&\frac{1}{C_{2}}\end{bmatrix},$ (105)

where $C_{\Sigma}=\det(C)=C_{1}C_{2}+C_{1}C_{g}+C_{2}C_{g}$. In the
approximation of the second step above, we assume that the shunting
capacitances are larger than the coupling capacitance, $C_{n}\gg C_{g}$, as is
usually the case. After rewriting to interacting harmonic oscillators the
diagonal elements of $\bm{C}^{-1}$ contribute to the respective modes with the
frequencies

$\omega_{n}=\sqrt{E_{C,n}E_{J,n}}+\alpha_{n},$ (106)

where the effective capacitive energy is $E_{C,n}=(C_{n}+C_{g})/C_{\Sigma}$
and the anharmonicity is $\alpha_{n}=-E_{C,n}$. The off-diagonal elements on
the other hand contribute to the interaction. The interaction term of the
Hamiltonian is

$\mathcal{H}_{\text{int}}=\frac{C_{g}}{C_{\Sigma}}q_{1}q_{2}.$ (107)

Quantizing the Hamiltonian and changing into annihilation and creation
operators the interaction part takes the form

$\mathcal{H}_{\text{int}}=g_{12}\left(\hat{b}_{1}^{\dagger}\hat{b}_{2}+\hat{b}_{2}\hat{b}_{1}^{\dagger}\right),$
(108)

where we remove terms that do not conserve the total number of excitations by
using the RWA. The coupling strength is

$g_{12}=\frac{C_{g}}{2C_{\Sigma}\sqrt{\zeta_{1}\zeta_{2}}},$ (109)

where $\zeta_{n}$ is the impedance in Eq. 56. Note the similarity with Eq. 61c
if one defines $E_{C,12}=C_{g}/2C_{\Sigma}$. Such a coupling is called a
transverse coupling since the interaction Hamiltonian only has nonzero matrix
elements in off-diagonal entries. This is contrary to the longitudinal
coupling discussed in Section IX.4.

###  IX.2 Two-qubit gates

As with the single-qubit gates in Section VIII.1, we can calculate the time-
evolution operator, as in Eq. 94 of the interacting Hamiltonian in order to
determine the gate operation. However, contrary to microwave driving we cannot
turn the interaction on and off directly. Luckily there are several approaches
to this problem, the simplest being tuning the two qubits in and out of
resonance such that the interaction terms time average to zero due to the RWA
discussed in Section VI. Examples of more complex and tunable coupling schemes
are discussed Section XI.3.

Consider the interaction part of the Hamiltonian in Eq. 108, we calculate the
time-evolution operator of the two-level truncation of this

$\displaystyle\hat{\mathcal{U}}(t)$
$\displaystyle=\exp\left[i\int_{0}^{t}\eta(t^{\prime})\hat{\mathcal{H}}_{\text{int}}\differential
t^{\prime}\right]$ (110)
$\displaystyle=\exp\left[i\Theta(t)(\sigma_{1}^{+}\sigma_{2}^{-}+\sigma_{1}^{-}\sigma_{2}^{+})\right]$
$\displaystyle=\begin{bmatrix}1&0&0&0\\\ 0&\cos\Theta(t)&-i\sin\Theta(t)&0\\\
0&-i\sin\Theta(t)&\cos\Theta(t)&0\\\ 0&0&0&1\end{bmatrix},$

where $\eta(t)$ is the envelope constructed so that it correspond to tuning
the two qubits in and out of resonance, and we assume that this is the only
part of the integral with time dependence. We also assume that the Hamiltonian
commutes with itself at different times. The angle of the coupling is given as

$\Theta(t)=g\int_{0}^{t}\eta(t^{\prime})\differential t^{\prime},$ (111)

which depends on the coupling strength, $g$, and the envelope $\eta(t)$. By
setting $\Theta(\tau)=\pi/2$ we obtain the $i$swap gate from Eq. 110 and
taking $\Theta(\tau)=\pi/4$ we find the $\sqrt{i\textsc{swap}}$ gate.

Note that a similar procedure to the $i$swap gate can be used to create a cz
gate [120].

###  IX.3 Linear resonators: control and measurement

So far we have considered how to engineer anharmonic oscillators and truncate
them into qubits as well as how to drive the qubits. However, for a qubit to
be useful we must also be able to control it and perform measurements on it
[3]. These two things can be accomplished by coupling the qubit to a linear
resonator, which is a simple harmonic oscillator [121].

Consider therefore the circuit presented in Fig. 12(b) consisting of a
transmonlike qubit capacitively coupled to an $LC$ oscillator or linear
resonator. This circuit is similar to the example circuit presented in Fig.
12(a) and the analysis up until truncation is identical with $1\rightarrow q$,
$2\rightarrow r$, and only one anharmonicity meaning that we must change
$-E_{J,2}\cos\phi_{2}$ to $\phi_{r}^{2}/2L$ in Eq. 103. Thus, we can truncate
only the mode with the anharmonicity which results in the following
Hamiltonian

$\hat{\mathcal{H}}_{\text{JC}}=\omega_{r}\hat{b}^{\dagger}\hat{b}+\frac{1}{2}\omega_{q}\sigma^{z}+g\left(\sigma^{+}\hat{b}+\sigma^{-}\hat{b}^{\dagger}\right),$
(112)

where $\hat{b}^{\dagger}$ and $\hat{b}$ are the creation and annihilation
operators for the linear resonator, $\sigma^{z}$ is the $z$ Pauli operator of
the qubit, and $\sigma^{\pm}$ represents the process of exciting and de-
exciting the qubit. The qubit frequency is given as in Eq. 106, the resonator
frequency is given by

$\omega_{r}=\sqrt{E_{C,r}/L},$ (113)

and the coupling strength is given as in Eq. 109

The Hamiltonian in Eq. 112 is known as the Jaynes-Cummings (JC) Hamiltonian,
which was initially used in quantum optics to describe a two-level atom in a
cavity [122, 123, 124]. Since then, the model has found application in many
areas of physics, including superconducting electronic circuits, where a qubit
is typically coupled to a transmission line resonator [125, 126, 127, 128,
129, 130, 131, 132, 133]. Because the Jaynes-Cummings Hamiltonian comes from
quantum optics and cavity quantum electrodynamics (cavity QED), coupling
between superconducting circuits and linear resonators is often denoted
circuit QED.

Consider the limit where the qubit frequency is far detuned from the resonator
frequency compared to the coupling rate and resonator linewidth
$\kappa=\omega_{r}/Q$, where $Q$ is the quality factor of the resonator, i.e.,
$\Delta=|\omega_{r}-\omega_{q}|\gg g,\kappa$. This is known as the dispersive
limit since there is no direct exchange of energy between the two systems,
i.e., only dispersive interactions between the resonator and the qubit occur.
Using second-order perturbation theory we see that the qubit and the resonator
change each other’s frequencies [134, 135, 136].

In the dispersive regime the Jaynes-Cummings Hamiltonian can be approximately
diagonalized using the unitary transformation $e^{\hat{S}}$ where
$\hat{S}=\lambda(\sigma^{+}\hat{b}-\sigma^{-}\hat{b}^{\dagger})$ and
$\lambda=g/\Delta$ is a small parameter. This transformation is called a
Schrieffer-Wolff transformation [137]. Using the Baker-Campbell-Hausdorff
formula [100] to second order in $\lambda$ we find the Hamiltonian in the
dispersive regime becomes

$\hat{\mathcal{H}}_{\text{disp}}=e^{\hat{S}}\hat{\mathcal{H}}_{\text{JC}}e^{-\hat{S}}=(\omega_{r}+\chi\sigma^{z})\hat{b}^{\dagger}\hat{b}+\frac{1}{2}\tilde{\omega}_{q}\sigma^{z},$
(114)

where we define $\chi=g^{2}/\Delta$ as the qubit dependent frequency shift or
dispersive shift. The qubit frequency is Lamb shifted to
$\tilde{\omega}_{q}=\omega_{q}+\chi$, induced by the vacuum fluctuations in
the resonator. Note that Eq. 114 is derived for a two-level atom/qubit. Taking
the second excited state into account modifies the expression for the shift
into

$\chi=-\frac{g_{01}^{2}}{\Delta}\left(\frac{1}{1+\Delta/\alpha}\right),$ (115)

where $g_{01}$ is the coupling rate between the 0 and 1 state of the qubit and
$\alpha$ is the anharmonicity of the qubit, Eq. 76.

One can interpret the dispersive qubit-resonator interaction in two ways.
Either as a shift of the qubit frequency by a quantity proportional to the
photon population of the resonator
$2\chi\langle\hat{b}^{\dagger}\hat{b}\rangle$ or as a qubit-dependent pull of
the resonator frequency, $\omega_{r}\rightarrow\omega_{r}\pm\chi$.

In the first interpretation the bare qubit frequency is modified by a Lamb
shift and by an additional amount proportional to the number of photons
populating the resonator. This is known as the ac Stark shift. It has the
consequence that fluctuations in the photon number of the resonator induce
small shifts in the qubit frequency, which brings it slightly out of its
rotating frame and cause dephasing [138, 139, 140, 141, 142, 143, 144]. In an
experiment, it is therefore important to reduce photon-number fluctuations of
the resonator, e.g., by keeping the process properly thermalized.

In the second interpretation, the resonator frequency is dependent on the
state of the qubit. This means that it is possible to make a quantum
nondemolition (QND) measurement of the qubit by shinning microwaves into the
resonator at a frequency close to $\omega_{r}$ and then measuring the
transmitted signal using standard homodyne techniques [145]. However, the
approximation in Eq. 114 is only valid in the small-photon limit, i.e., when
when the resonator photon number, $\hat{N}=\hat{b}^{\dagger}\hat{b}$ is less
than the critical photon number $N_{c}=\Delta^{2}/4g^{2}$. This sets an upper
limit to the power of the resonator as a probe while maintaining the
conditions for a QND measurement. However, this is not the whole story; Ref.
[132] has shown that level crossings with other states of the qubit-resonator
system induce state transitions which can be explained by the Jaynes-Cummings
model. This is beyond the scope of our discussion, and we refer to the cited
work for more information.

In the other limit, when the detuning between the qubit and the resonator
frequency is small compared to the coupling rate, i.e., $\Delta\ll g$ we
obtain a hybridization of the energy levels of the two systems. This opens up
for a Rabi mode splitting, where each transition between the qubit and the
resonator splits into two states with distance $\sqrt{N}g/\pi$ where $N$
denotes the resonator mode, i.e., the photon number. Thus excitation is
coherently swapped between the two systems. While this cannot be used to
perform measurements on the qubit it can be used to mediate couplings between
two qubits by coupling another qubit to the resonator [146, 147]. We do not
dive deeper into the details of measurements and couplings to linear
resonators. For an experimental-minded review see, e.g., Ref. [11].

###  IX.4 Inductive coupling

Figure 13: Mutual inductive coupling between two modes.

So far we have considered only direct inductance as a way of coupling two
qubits. In this section we consider the mutual inductance of two modes as a
means of coupling the modes. Consider therefore the two circuits in Fig. 13,
consisting of a Josephson junction, a capacitor, and a linear inductor. Such
circuits are known as rf SQUIDs [148]. Each of the circuits have the following
Hamiltonian

$\hat{\mathcal{H}}_{j}=4E_{C,j}\hat{n}_{j}^{2}+\frac{1}{2L_{j}}\hat{\phi}_{j}^{2}-E_{J,j}\cos(\hat{\phi}_{j}+\tilde{\Phi}_{j}).$
(116)

If two such circuits are brought into proximity of each other they will share
a mutual inductance, yielding an interaction Hamiltonian

$\displaystyle\hat{\mathcal{H}}_{\text{int}}$
$\displaystyle=M_{12}\hat{I}_{1}\hat{I}_{2}$ (117)
$\displaystyle=M_{12}I_{c1}\sin(\hat{\phi}_{1}+\tilde{\Phi}_{j})I_{2c}\sin(\hat{\phi}_{2}+\tilde{\Phi}_{j}),$

where $\hat{I}_{j}$ is the current operator of the Josephson junction, see Eq.
11. The mutual inductance $M_{12}$ between the two circuits depends on the
relative geometrical placement of the circuits. This can be increased, e.g.,
by overlapping the circuits [149] or by letting them share the same wire or
Josephson-junction inductor [150, 151, 152, 153].

Consider now the case of no external flux, i.e., $\tilde{\Phi}=0$. If we
expand the potential to fourth order the interaction Hamiltonian takes the
form

$\hat{\mathcal{H}}=M_{12}\left[\hat{\phi}_{1}\hat{\phi}_{2}-\frac{1}{36}(\hat{\phi}_{1}\hat{\phi}_{2}^{3}+\hat{\phi}_{1}^{3}\hat{\phi}_{2})\right].$
(118)

Truncating into a two-level model using Tables 2 and 3 we find that the
coupling becomes transverse

$\hat{\mathcal{H}}_{\text{int}}=g_{x}\sigma_{1}^{x}\sigma_{2}^{x},$ (119)

where the coupling constant is

$g_{x}=\frac{1}{2}M_{12}\sqrt{\zeta_{1}\zeta_{2}}\left[1-\frac{1}{24}(\zeta_{1}+\zeta_{2})\right],$
(120)

with impedances given by Eq. 56.

Consider now an external flux of $\tilde{\Phi}=\pi/2$. In this case, the sine
terms obtain a phase effectively changing the terms into cosines. Expanding
these to second-order yields

$\hat{\mathcal{H}}_{\text{int}}=\frac{M_{12}}{4}\left[-(\hat{\phi}_{1}^{2}+\hat{\phi}_{2}^{2})+\frac{1}{144}(\hat{\phi}_{1}^{4}+\hat{\phi}_{2}^{4})+\hat{\phi}_{1}^{2}\hat{\phi}_{2}^{2}\right],$
(121)

where we recognize the two first terms as corrections to the qubit
frequencies, the next two terms as corrections to the anharmonicities, and
finally the last term is the interaction term. Considering only the last term
and truncating into a two-level model we find (see Tables 2 and 3)

$\hat{\mathcal{H}}^{\prime}_{\text{int}}=g_{z}(\sigma_{1}^{z}\sigma_{2}^{z}-2\sigma_{1}^{z}-2\sigma_{2}^{z}).$
(122)

The first term is a longitudinal coupling between the two qubits with coupling
constant

$g_{z}=\frac{1}{16}M_{12}\zeta_{1}\zeta_{2}.$ (123)

It is called longitudinal because all off-diagonal matrix elements are zero,
contrary to transverse coupling. Longitudinal coupling can be used to create
entanglement without exchanging energy between the modes by enabling a so-
called phase gate [146, 154, 11]. From the last two terms in Eq. 122, we see
that we obtain further corrections to the qubit frequencies.

##  X Noise and decoherence

So far we have considered only closed quantum systems, i.e., systems without
interaction with the environment. This is usually a good approximation as we
are dealing with cryogenic and thus isolated superconducting circuits.
However, even in the best experimental setups, random and uncontrollable
processes in the environment surrounding the system do occur. These are
sources of noise and lead to decoherence of the quantum system. It is,
therefore, necessary to develop a formalism to treat this theoretically as
well. We assume that the Hamiltonian of the system and the environment is
separable and has the form

$\hat{\mathcal{H}}=\hat{\mathcal{H}}_{\text{sys}}+\hat{\mathcal{H}}_{\text{env}}+\nu\hat{S}\cdot\hat{\lambda},$
(124)

where $\hat{\mathcal{H}}_{\text{sys}}$ is the Hamiltonian of the system,
$\hat{\mathcal{H}}_{\text{env}}$ is the Hamiltonian of the environment, and
the interaction strength between the system and the environment is given by
$\nu$, while $\hat{S}$ is an operator within the system Hamiltonian
$\hat{\mathcal{H}}_{\text{sys}}$ and $\hat{\lambda}$ represents the noisy
environment which produces fluctuations $\delta\lambda$.

The treatment of open quantum systems is a whole subject on its own and a
complete treatment is beyond the scope of this tutorial. We therefore present
only a method for modeling noise in qubit systems. For a more extensive
treatment of open quantum systems see e.g. Ref. [155], and for an introduction
on how to treat noise in an experiment see, e.g., Ref. [11].

###  X.1 Bloch-Redfield model

Consider an arbitrary state on the Bloch sphere as in Eq. 79. The density
matrix for such a pure state is [12]

$\rho=\outerproduct{\psi}{\psi}=\frac{1}{2}(I+\bm{a}\cdot\bm{\sigma})=\begin{bmatrix}|\alpha|^{2}&\alpha\beta^{*}\\\
\alpha^{*}\beta&|\beta|^{2}\end{bmatrix},$ (125)

where $I$ is the identity matrix, $\bm{a}$ is the Bloch vector, and
$\bm{\sigma}=(\sigma^{x},\sigma^{y},\sigma^{z})$ is the vector of Pauli spin
matrices. If $\rho$ represents a pure state, $\psi$, then $\tr\rho^{2}=1$ and
the Block vector becomes a unit vector,
$\bm{a}=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$, where $\theta$
and $\phi$ are the angles of the Bloch vector. If, on the other hand,
$|\bm{a}|<1$ the density matrix $\rho$ represents a mixed state with
$\tr\rho^{2}<1$. In this case the Bloch vector terminates at points inside the
unit sphere.

In the Bloch-Redfield formulation of two-level systems, sources of noise are
weakly coupled to the system with short correlation times compared to the
system dynamics [156, 157, 158, 155]. The noise in this formulation is
determined by two rates: The longitudinal relaxation rate and the transverse
relaxation rate.

Figure 14: Bloch-sphere representation of noise. (a) Longitudinal relaxation
is the result of energy exchange between the qubit and the environment.
Transverse noise couples to the qubit and drives a rotation (transition)
around an axis in the $x$-$y$ plane. Longitudinal relaxation is driven both by
emission of energy to the environment, $\Gamma_{-}$ and absorption of energy
from the environment, $\Gamma_{+}$. For a typical superconducting qubit, the
temperature is much lower than the frequency of the qubit,
$k_{B}T\ll\hbar\omega$, which suppresses the absorption rate, such that
$\Gamma_{1}\simeq\Gamma_{-}$. (b) Pure dephasing is the result of longitudinal
noise that drives a rotation around the $z$ axis. Due to stochastic frequency
fluctuations, a Bloch vector will diffuse both clockwise and counterclockwise
around the $z$ axis parallel to the equator.

####  X.1.1 Longitudinal relaxation

The longitudinal relaxation rate, $\Gamma_{1}=1/T_{1}$, describes
depolarization along the qubit quantization axis, often referred to as “energy
decay” or “energy relaxation”, which is why it is often referred to as the
relaxation time. Longitudinal relaxation is caused by transverse noise, via
the $x$ or $y$ axis on the Bloch sphere, see Fig. 14(a). Depolarization of the
superconducting circuit occurs due to exchange of energy with the environment,
leading both to excitation and relaxation of the qubits, meaning that one can
write

$\Gamma_{1}=\Gamma_{+}+\Gamma_{-}.$ (126)

Due to Boltzmann statistics and the fact that superconducting qubits are
operated at low temperatures ($T\lesssim$20\text{\,}\mathrm{mK}$$) and with a
qubit frequency in the $\mathrm{GHz}$ regime, the qubits generally lose energy
to the environment, meaning that the excitation rate $\Gamma_{+}$ is
suppressed exponentially as

$\frac{\Gamma_{+}}{\Gamma_{-}}=e^{-\beta\omega},$ (127)

where $\beta=1/k_{B}T$ is the inverse of the Boltzmann constant multiplied
with the temperature. Note, however, that empirically we often see stray
population of the excited state much higher than we would expect from this
theory. From Maxwell-Boltzmann statistics, we would expect a thermal
population of the excited state of $P_{\ket{1}}\sim 10^{-5}$%, but the
measured excited-state population is often orders of magnitudes higher at
around 1% [159].

The longitudinal relaxation rate can be determined using Fermi’s golden rule

$\Gamma_{1}=\frac{1}{\hbar^{2}}|\matrixelement{0}{\hat{S}}{1}|^{2}S_{\lambda}(\omega_{q}),$
(128)

where $\hat{S}$ is the transverse coupling of the qubit to the environment,
i.e., a coupling of the type $\sigma^{x}$ or $\sigma^{y}$. The qubit frequency
is denoted $\omega_{q}$. The noise power spectral density

$S_{\lambda}(\omega)=\int_{\infty}^{-\infty}\differential
t\langle\hat{\lambda}(t)\hat{\lambda}(0)\rangle e^{-i\omega t},$ (129)

characterizes the frequency distribution of the noise power for a stationary
noise process $\hat{\lambda}$. Note that the Wiener-Khintchine theorem states
that $S_{\lambda}(\omega)$ is the Fourier transform of the autocorrelation
function $c_{\lambda}(t)=\langle\hat{\lambda}(t)\hat{\lambda}(0)\rangle$ of
the noise source $\lambda$ [160, 161].

The longitudinal relaxation rate can be measured by preparing the qubit in
state $\ket{1}$ (e.g. using a $\pi$ pulse as in Section VIII) and then making
multiple measurements of the qubit excited-state population at a set of
subsequent times [11].

####  X.1.2 Transverse relaxation

The transverse relaxation time

$\Gamma_{2}=\frac{1}{T_{2}}=\frac{\Gamma_{1}}{2}+\Gamma_{\phi}$ (130)

describes the loss of coherence of a superposition. As seen in Eq. 130 it is
caused both by transverse noise, which leads to energy (longitudinal)
relaxation of the excited-state component of the superposition state, and by
longitudinal noise, which cause fluctuations of the qubit frequency and leads
to pure dephasing, see Section X.1.3 below. Note that the sum in Eq. 130 is
only valid for weak noise that also is only correlated at short times [162].

We introduce transverse relaxation as it is a measurable quantity, contrary to
pure dephasing which can only be inferred using Eq. 130. Transverse relaxation
can be measured using Ramsey interferometry [163, 11]. In Ramsey
interferometry, a $\pi/2$ pulse rotates the Bloch vector from $\ket{0}$ to the
equator of the Bloch sphere. If we know the qubit frequency perfectly, it
should remain stationary at the equator, and if we apply another $\pi/2$ pulse
at some time later, we should measure $\ket{1}$. However, if our knowledge of
the qubit frequency and our assumed frame does not match the qubit’s actual
rotation frame, then the state will not remain stationary at the equator of
the Bloch sphere after the first $\pi/2$ pulse is applied. Instead, it will
precess around the equator at a frequency equal to the difference between the
assumed frame and the actual qubit frequency. This means that if we perform
two $\pi/2$ pulses with variable delay in between, we should observe
oscillations in the measured state of the qubit. In reality, one often chooses
a frame that is intentionally detuned from the qubit frequency so that these
oscillations are observed even for the perfectly calibrated qubit. This means
that an error in qubit frequency will result in a difference from the expected
oscillation frequency [164]. For simple Markovian noise, these oscillations
are exponentially damped with characteristic time $T_{2}$ [11].

####  X.1.3 Pure dephasing

The pure dephasing rate $\Gamma_{\phi}$ describes depolarization in the
$x$-$y$ plane of the Bloch sphere. It is referred to as “pure dephasing,” to
distinguish it from other phase-breaking processes such as energy excitation
or relaxation. Pure dephasing is caused by longitudinal noise that couples to
the qubit via the $z$ axis. This longitudinal noise causes the qubit
frequency, $\omega$, to fluctuate such that it is no longer equal to the
interaction frame frequency, causing the Bloch vector to precess forward or
backward in the interacting frame as seen in Fig. 14(b).

To lowest order, the pure dephasing rate is orthogonal to the difference
between the two diagonal matrix elements [155, 74]

$\Gamma_{\phi}=\frac{1}{\hbar^{2}}\left(|\matrixelement{0}{\hat{S}}{0}-\matrixelement{1}{\hat{S}}{1}|^{2}\right)S_{\lambda}(0),$
(131)

where $\hat{S}$ is the longitudinal coupling of the qubit to the environment,
i.e., a coupling of the type $\sigma^{z}$. This means that pure dephasing
disappears if
$\matrixelement{0}{\hat{\mathcal{H}}_{\text{env}}}{0}=\matrixelement{1}{\hat{\mathcal{H}}_{\text{env}}}{1}$.
For superconducting circuits, this can often be realized by tuning the system
to the so-called “sweet spot” using external flux biasing. This means that the
transverse relaxation becomes approximately half the longitudinal relaxation
rate as in Eq. 130. Thus decreasing the longitudinal relaxation rate becomes
the main focus when developing qubits with the goal of increasing the
lifetime. Note, however, in reality, pure dephasing will never disappear