# Circuit quantum electrodynamics (cQED) with modular quasi-lumped models

Zlatko K. Minev<EMAIL_ADDRESS>www.zlatko-minev.com Thomas G. McConkey
Maika Takita  Antonio D. Corcoles  Jay M. Gambetta  IBM Quantum, IBM T.J.
Watson Research Center, Yorktown Heights, US

(September 3, 2024)

###### Abstract

Extracting the Hamiltonian of interacting quantum-information processing
systems is a keystone problem in the realization of complex phenomena and
large-scale quantum computers. The remarkable growth of the field increasingly
requires precise, widely-applicable, and modular methods that can model the
quantum electrodynamics of the physical circuits, including their more-subtle
renormalization effects. Here, we present a computationally-efficient method
satisfying these criteria. The method partitions a quantum device into compact
lumped or quasi-distributed cells. Each is first simulated individually. The
composite system is then reduced and mapped to a set of simple subsystem
building blocks and their pairwise interactions. The method operates within
the quasi-lumped approximation and, with no further approximation,
systematically accounts for constraints, couplings, parameter
renormalizations, and non-perturbative loading effects. We experimentally
validate the method on large-scale, state-of-the-art superconducting quantum
processors. We find that the full method improves the experimental agreement
by a factor of two over taking standard coupling approximations when tested on
the most sensitive and dressed Hamiltonian parameters of the measured devices.

Quantum phenomena offer a distinct advantage for information processing,
assuming we can faithfully design and realize the physical systems underlying
them. A leading platform to accomplish this goal has emerged in the form of
superconducting quantum technology (Devoret and Schoelkopf, 2013; Gambetta
_et al._ , 2017; Krantz  _et al._ , 2019; Kjaergaard  _et al._ , 2020; Blais
_et al._ , 2020), which employs macroscopic, lithographically-defined, and
configurable devices. Their design versatility, however, comes with inherent
challenges—parameter variability, a complicated design space, and the
difficulty of engineering their non-linear interactions. These challenges pose
a concern of central importance for the growth of the field and have received
a strong and growing interest in the form of new design and quantization
methods (Nigg  _et al._ , 2012; Bourassa  _et al._ , 2012; Solgun  _et al._ ,
2014; Solgun and DiVincenzo, 2015; Smith  _et al._ , 2016; Malekakhlagh and
Türeci, 2016; Gely  _et al._ , 2017; Malekakhlagh  _et al._ , 2017; Pechal and
Safavi-Naeini, 2017; Minev  _et al._ , 2018; Parra-Rodriguez  _et al._ , 2019,
2018; Ansari, 2019; Krupko  _et al._ , 2018; Malekakhlagh  _et al._ , 2020a;
Solgun  _et al._ , 2019; Petrescu  _et al._ , 2019; You  _et al._ , 2019; Di
Paolo  _et al._ , 2019; Gely and Steele, 2020; Ding  _et al._ , 2020; Kerman,
2020; Minev  _et al._ , 2020; Kyaw  _et al._ , 2020; Malekakhlagh  _et al._ ,
2020b; Menke  _et al._ , 2021). The rapidly growing pace of diversity,
complexity, and scale of quantum hardware (Wallraff  _et al._ , 2004; Paik
_et al._ , 2011; Barends  _et al._ , 2013; Minev  _et al._ , 2013; Yan  _et
al._ , 2016; Brecht  _et al._ , 2016; Rosenberg  _et al._ , 2017; Versluis
_et al._ , 2017; Naik  _et al._ , 2017; Charles James Neill, 2017; Yan  _et
al._ , 2018; Minev  _et al._ , 2019; Antonio D.  _et al._ , 2021) urges for
methods that are increasingly modular, widely-applicable, yet ever-more
precise. As such, these methods must closely incorporate details of the
layout, materials, and electromagnetic environment of the physical device,
while keeping approximations to a minimum without sacrificing computational
efficiency.

Here, we develop such a precise, modular method that operates at the physical-
device level and that is suitable for a wide array of quantum devices. The
method builds on the quantization of lumped models (Yurke and Denker, 1984;
Devoret, 1995; Burkard  _et al._ , 2004), which is more computationally
efficient than full-wave methods (Nigg  _et al._ , 2012; Solgun and
DiVincenzo, 2015; Minev  _et al._ , 2020). The method also handles distributed
resonant structures, such as co-planar waveguide (CPW) resonators. Modularity
is achieved in two ways. First, the physical layout of the quantum device is
systematically partitioned into disjoint cells—physical blocks of the device.
Each cell can be independently simulated to extract its electromagnetic
parameters. Second, the effective circuit of the device is partitioned into
non-dynamical coupler elements and nodes and into subsystem building blocks. A
subsystem is selected based on two requirements: i) it should constitute a
well-understood, basic system in isolation and ii) it should have sufficient
parameter flexibility. The latter enables the faithful reduction and mapping
of the larger device into subsystems, without loss of information. The larger
device model is stitched together from the results of the cell simulations. It
is then reduced according to the subsystem partitions. In the process,
constraints and non-dynamical degrees of freedom associated with the coupling
structures are systematically eliminated; dressing of the systems and their
interactions because of this are accounted for. The reduction is precise in
that no approximations are made. All dressing and effective parameter
renormalizations of the subsystems are tracked, which are mediated by the
coupling structures.

We experimentally tested the method on two large-scale, superconducting quantum processors (Antonio D.  _et al._ , 2021), which employed transmon qubits (Koch  _et al._ , 2007) and CPW structures. Each qubit was connected to four or five neighboring structures, spanning a wide range of coupling strengths. We observed renormalizations on the order of 25% for subsystem Hamiltonian parameters due to coupling dressing and a series of smaller renormalizations due to unwanted, indirect couplings and higher-order effects. We provide a detailed budget to account for the extent to which each model parameter influenced the results. We compare the full method presented here to one that resorts to standard approximations, weak coupling and no dressing of spatial eigenfields. The full method presented here yields a factor of two improvement on the experimental agreement. The method is applicable to the analysis of a broad class of quantum processors. We have automated it in the open-source project Qiskit Metal | for quantum device design 111For an early version of the open-source code (Minev  _et al._ , 2021) developed by the authors of this manuscript, see qiskit.org/metal. The Qiskit Metal project builds on github.com/zlatko-minev/pyEPR..

##  I Partitioning a quantum device into interconnected component cells

Figure 1: Method overview. (a) Illustration of an example quantum processor
layout (not-to-scale; partial). Center: transmon qubit subsystem (purple),
described by a dressed Hamiltonian $\hat{H}_{0}$, connected to three co-planar
waveguide (CPW) subsystems—a readout (green, $\hat{H}_{1}$) and two bus
resonators (brown, $\hat{H}_{2}$, and blue, $\hat{H}_{3}$)—and a charge line
(bottom, gray). The device layout is partitioned into subsystems and cells.
(b) Example simulation model of $\mathrm{cell}_{0}$, overlaid with its
simulation mesh. The cell incorporates elements of multiple subsystems: qubit
pads $P_{0}$ and $P_{1}$, CPW coupler pads and segments $B_{1},B_{2},B_{3}$,
and a segment of a CPW charge line $C$. (c) Partial schematic of the
composite-system network depicting the nodes and elements of
$\mathrm{cell}_{0}$ and their connections to neighboring cells, such as those
of the CPWs. The nodes of the cell are capacitively coupled by a fully-
connected graph (thick line). A Josephson tunnel junction connects $P_{0}$ to
$P_{1}$. Coupler nodes $C$ and $P_{0}$ (not colored) are eliminated. (d)
Depiction of the reduced, dressed subsystems—acting as building blocks—and
their dressed interactions, described by the Hamiltonians
$\hat{H}_{01},\hat{H}_{02},\ldots$

The physical layout of a quantum information processor, such as the one
depicted in Fig. 1(a), can be thought of as the interconnected collection of
quasi-independent subsystems and inter-system couplers. Each subsystem is
identified with a domain of the physical layout and supports a set of
quantized, potentially-anharmonic modes. When considered in isolation, a
subsystem could be one that is well understood and whose parameters could be
readily obtainable, using analytical or numerical techniques. However, due to
its embedding in the larger device network and loading by system-system
coupling structures, its properties can be dressed and can significantly
depend on parts of the larger network. Generally, this dressing of the
Hamiltonian parameters is necessarily unavoidable and non-local. Even for
relatively moderate non-linear coupling strengths, the renormalizations of the
effective subsystem Hamiltonian and its spectrum can be significant; for an
experimental example of a large, 1 GHz mode dressing, see Sec. IV. Similarly,
the interaction between two subsystems is dressed in a potentially non-local
manner that can depend on more than the immediate coupling structure and can
be influenced by the physical layout more broadly. Faithfully and
systematically accounting for such effects in the increasingly complex
physical device layouts is important for improved experimental agreement, see
Fig. 4.

To capture these effects, we aim to systematically and modularly study the
physical device layout by partitioning it into disjoint cells, see Fig. 1(b).
Each cell is first analyzed quasi-analytically or using numerical quasi-static
electromagnetic methods to extract device parameters used to construct the
full system Hamiltonian. Importantly, cells do not correspond to subsystems,
see Sec. II. A semi-classical model of the connected system is constructed and
reduced to a simpler, dressed model of the device, see Fig. 1(b) and 1(c).

For a subsystem coupled to $K$ neighbors, we explicitly construct the
composite-system Hamiltonian $\hat{H}_{\mathrm{full}}$, within the quasi-
lumped approximation but with no further approximations, in the pair-wise
interaction form,

$\hat{H}_{\mathrm{full}}=\hat{H}_{0}+\sum_{n=1}^{K}\hat{H}_{n}+\sum_{n=0}^{K-1}\sum_{m=n+1}^{K}\hat{H}_{nm}\;,$
(1)

where $\hat{H}_{0}$ and $\hat{H}_{n}$ are the dressed Hamiltonians of the
subsystem and its $n$-th neighbor, respectively, and $\hat{H}_{nm}$ is the
dressed, bi-linear interaction Hamiltonian between the $n$-th and $m$-th
systems.

In the construction, the $n$-th and $m$-th system Hilbert spaces are disjoint;
hence, $\left[\hat{H}_{n},\hat{H}_{m}\right]=0$. Nonetheless, their dressed
Hamiltonians $\hat{H}_{n}$ and $\hat{H}_{m}$ are interdependent in terms of
their parameters and physical layouts—due to the coupling structures and
dressing by the larger network.

The dressed interaction Hamiltonian of the $n$-th and $m$-th systems is

$\hat{H}_{nm}=\hat{Q}_{n}\hat{Q}_{m}/\text{$C_{nm}^{\mathrm{eff}}$}+\hat{\Phi}_{n}\hat{\Phi}_{m}/\text{$L_{nm}^{\mathrm{eff}}$}\;,$
(2)

where the effective coupling capacitance and inductance are
$C_{nm}^{\mathrm{eff}}$ and $L_{nm}^{\mathrm{eff}}$, respectively, and the
generalized charge and magnetic flux operators of the $n$-th and $m$-th
systems involved in the coupling are $\hat{Q}_{n}$, $\hat{Q}_{m}$,
$\hat{\Phi}_{n}$, and $\hat{\Phi}_{m}$. For systems coupled by multiple
physical degrees of freedom, see the more general form of Eq. (2) obtained in
Eqs. (12) and (13).

In the case of purely capacitive inter-system coupling (i.e.,
$1/L_{nm}^{\mathrm{eff}}=0$), we conveniently reexpress the interaction in
terms of the dimensionless operators $\hat{A}_{n}\coloneqq\hat{Q}_{n}/A_{n}$
and $\hat{B}_{m}\coloneqq\hat{Q}_{m}/B_{m}$, scaled by the choice scaling
factors $A_{n}$ and $B_{m}$,

$\hat{H}_{nm}=\hat{Q}_{n}\hat{Q}_{m}/\text{$C_{nm}^{\mathrm{eff}}$}=\hbar
g_{nm}\hat{A}_{n}\hat{B}_{m}\;,$ (3)

where the linear coupling energy $\hbar g_{nm}\coloneqq
A_{n}B_{m}/\text{$C_{nm}^{\mathrm{eff}}$}$. When operating with continuous
quadrature variables, such as
$\hat{A}_{n}=i(\hat{a}_{n}^{\dagger}-\hat{a}_{n})$, and employing a harmonic
basis, one typically employs the zero-point quantum fluctuations of the
quadrature with respect to the linearized system; i.e.,
$\text{$A_{n}$}^{2}=\left(A_{n}^{\mathrm{ZPF}}\right){}^{2}\coloneqq\left\langle\hat{A}_{n}^{2}\right\rangle-\left\langle\hat{A}_{n}\right\rangle^{2}$.
Sometimes in the case of a transmon qubit (Koch  _et al._ , 2007), especially
when interested in charge effects, one operates with wrapped or discrete-
variable quadratures. In the case of the discrete-variable Cooper pair number
operator $\hat{A}_{n}=\hat{n}$, the scaling is $A_{n}=2e$, where $e$ is the
elementary electron charge.

Before proceeding to the general treatment and explicit construction of Eqs.
(1)–(5), for definitiveness of example, consider the device and associated
reduction steps depicted in Fig. 1. The dressed Hamiltonian of the $k=1$
system, a co-planar-waveguide (CPW) resonator, is
$\hat{H}_{1}=\sum_{m=1}^{\infty}\hbar\omega_{1m}\left(\text{$\hat{a}_{1m}^{\dagger}\hat{a}_{1m}$}+\frac{1}{2}\right)$,
see Sec. III. Due to the coupling and network dressings, the normal mode
frequencies $\omega_{1m}$, annihilation operators $\hat{a}_{1m}$, and
associated zero-point-quantum fluctuations of the physical degrees of freedom
are not dependent solely on the geometry and materials of the CPW structure.
They additionally include contributions arising from the coupler structure and
transmon qubit self-capacitance. Similarly, the $k=0$ qubit system frequency
$\omega_{q}$ of $\hat{H}_{0}$ is weakly interdependent with $\omega_{1m}$ due
to the semi-classical hybridization; i.e., when one tunes $\omega_{q}$, one
has to be careful to also retune $\omega_{1m}$ for maximum accuracy. The
coupling parameters $C_{nm}^{\mathrm{eff}}$ and $L_{nm}^{\mathrm{eff}}$ are
dressed not only by the $n$-th and $m$-th systems, but also by couplers
coupling the $n$-th system to its other neighbors; e.g., the qubit-bus
coupling can weakly influence the qubit-readout coupling.

In the following section, we capture these potentially non-perturbative
effects, see Sec. III. Such a detailed construction is required to treat
larger and higher-order coupling effects and to obtain improved experimental
agreement, see Sec. IV. Additionally, effects such as impurity scattering on
the distributed CPW structures are accounted for—this is the cavity-QED
equivalent of a gauge-dependent diamagnetic $A^{2}$ contribution (Malekakhlagh
and Türeci, 2016).

##  II Theory of the composite model

In this section, we present the theory of the general construction of Eqs. (1)
and (2) starting from the analysis of disjoint partitions of the physical
device layout. The model is constructed to allow the mapping of the larger
interconnected network to a set of small, independent, subsystem building
blocks—each potentially well understood in isolation. Each such cell partition
can be independently simulated, thus decoupling the complexity of the network
and providing modular simulations for improved computational efficiency. From
the results of these independent simulations, we construct the composite-
system Lagrangian and then the quantum Hamiltonian of the systems and
interactions. In the process, we eliminate holonomic constraints and non-
dynamical, coupling degrees of freedom, and account for their dressing of the
systems and couplings. The treatment in this section, for generality, is more
formal and abstract. However, in practice, significant simplifications occur
due to the typically diagonal, block, or sparse structure of the circuit
matrices.

We simultaneously employ two different types of partitions of the device. The
first partitions the physical layout of the device into $N_{\mathrm{cell}}+1$
strictly disjoint cell modules—each of which can be independently simulated or
analyzed to extract effective circuit parameters. The second partition
operates at the level of the effective quasi-lumped schematic of the device.
It partitions the device network nodes and elements into $K+2$ subsets of
each. One of the $K+2$ subsets is for all non-system couplers. The remaining
$K+1$ are for the $K+1$ subsystems. The cell and system partitions of the
composite circuit are distinct from each other. This is illustrated by the
qubit cell of Fig. 1(b). It incorporates the transmon qubit elements (metal
pads $P_{0}$ and $P_{1}$ in a cutout of the ground plane $G$) but also
elements from three neighboring non-qubit systems (structures
$B_{1},B_{2},B_{3,}$ and $C$).

Let $\mathcal{N}_{\mathrm{full}}$ and $\mathcal{B}_{\mathrm{full}}$ denote the
sets of all nodes and elements in the device network, respectively. To
partition into cells, the elements of $\mathcal{B}_{\mathrm{full}}$ are
disjointly distributed among $N_{\mathrm{cell}}+1$ sets. The nodes in
$\mathcal{N}_{\mathrm{full}}$ are however distributed among sets with
potential overlap. If we denote the set of all non-datum nodes assigned to the
$n$-th cell as $\mathcal{N}_{n}^{\mathrm{cell}}$, then
$\mathcal{N}_{\mathrm{full}}=\mathcal{N}_{g}\cup\bigcup_{n=0}^{N_{\mathrm{cell}}}\mathcal{N}_{n}^{\mathrm{cell}}$,
where $\mathcal{N}_{g}$ is the single-element set comprising just the circuit
datum (ground node). The intersection
$\mathcal{N}_{n}^{\mathrm{cell}}\cap\mathcal{N}_{m}^{\mathrm{cell}}$ can be
non-empty. Each cell can incorporate the datum. To partition into subsystems,
$\mathcal{B}_{\mathrm{full}}$ and $\mathcal{N}_{\mathrm{full}}$ are broken up
into strictly disjoint sub-sets. If $\mathcal{N}_{k}$ denotes the set of all
non-datum nodes assigned to the $k$-th system and
$\mathcal{N}_{\mathrm{couple}}$ denotes the set of all non-system, non-datum
nodes, then
$\mathcal{N}_{\mathrm{full}}=\mathcal{N}_{g}\cup\mathcal{N}_{\mathrm{couple}}\cup\bigcup_{k=0}^{K}\mathcal{N}_{k}$,
where
$\mathcal{N}_{k}\cap\mathcal{N}_{k^{\prime}}=\mathcal{N}_{k}\cap\mathcal{N}_{\mathrm{couple}}=\emptyset$
for $k\neq k^{\prime}$. In the case of a continuous subsystem,
$\mathcal{N}_{k}$ comprises a continuum of nodes.

The assignment of nodes and elements to a subsystem is determined by the
design intention to realize quasi-independent subsystems, incorporating
dynamical degrees of freedom that support quantized modes. A subsystem can
also incorporate non-dynamical nodes and zero-frequency modes; e.g., a
transmission line open at both ends. The system’s non-dynamical degrees will
be preserved in the following reduction. This feature provides helpful
flexibility that allows the mapping of the larger network to a set of known
building blocks.

A second, necessary condition for a subsystem is that it admits sufficient
parameter flexibility in its definition. It must be able to admit the
renormalization of the larger network. For example, the simple, ideal
transmission line is not a suitable subsystem candidate. It cannot admit in
its homogenous construction the inhomogeneity introduced by the impurity
scattering effect of a coupler structure. The line must be allowed to
incorporate the inhomogeneity in line parameters due to the coupler structure
in its definition; see Fig. 2 for an example.

Nodes and elements not part of any subsystem are assigned to the coupler sets.
The coupler nodes should all be non-dynamical—they do not support quantized
modes on their own. A sufficient condition to identify a non-dynamical node is
that it is touched by only inductive or capacitive elements, but not by
elements of both classes.

Inductive elements can be linear or non-linear. A non-linear inductive dipole
is a two-terminal purely-inductive sub-circuit, such as a Josephson tunnel
junction, flux-biased SQUID (Zimmerman and Silver, 1966; Clarke and Braginski,
2004), SNAIL (Frattini  _et al._ , 2018), or a more complicated composite sub-
circuit. We describe such dipoles using the formulation of Ref. Minev  _et
al._ , 2020. The $j$-th non-linear dipole is fully characterized by a known
energy function
$\mathcal{E}_{j}\text{$\left(\Phi_{j};\Phi_{j}^{\mathrm{ext}}\right)$}$ of the
generalized magnetic flux $\Phi_{j}$ across it and any external bias
parameters $\Phi_{j}^{\mathrm{ext}}$, such as the magnetic flux of a d.c.
voltage bias. The dipole can also have an intrinsic capacitance $C_{j}$ that
spans its two terminals. The energy function of the dipole intrinsic
capacitance is $\frac{1}{2}C_{j}\dot{\Phi}_{j}^{2}$.

With respect to the circuit operating point (Minev, 2019; Minev  _et al._ ,
2020), one can partition the inductive energy into strictly linear and non-
linear contributions

$\displaystyle\mathcal{E}_{j}^{\mathrm{lin}}\text{$\left(\Phi_{j};\Phi_{j}^{\mathrm{ext}}\right)$}$
$\displaystyle\coloneqq\frac{1}{2}L_{j}^{-1}\left(\Phi_{j}^{\mathrm{ext}}\right)\Phi_{j}^{2}\;,$
(4a)
$\displaystyle\mathcal{E}_{j}^{\mathrm{nl}}\text{$\left(\Phi_{j};\Phi_{j}^{\mathrm{ext}}\right)$}$
$\displaystyle\coloneqq\mathcal{E}_{j}\left(\Phi_{j};\Phi_{j}^{\mathrm{ext}}\right)-\frac{1}{2}L_{j}^{-1}\left(\Phi_{j}^{\mathrm{ext}}\right)\Phi_{j}^{2}\;,$
(4b)

respectively, where $L_{j}^{-1}\left(\Phi_{j}^{\mathrm{ext}}\right)$ is the
linear-response inductance of the dipole at the bias point. For example, for a
SQUID dipole,
$\mathcal{E}_{j}\left(\Phi_{j};\Phi_{j}^{\mathrm{ext}}\right)=-E_{j}\left(\Phi_{j}^{\mathrm{ext}}\right)\cos\left(\Phi_{j}/\phi_{0}\right)$,
where $E_{j}$ is the effective Josephson energy as a function of the bias and
$\phi_{0}\coloneqq\hbar/2e$. For a Josephson tunnel junction, $L_{j}$ is
simply the Josephson inductance, which can be computed from the Ambegaokar-
Baratoff expression and room-temperature resistance measurements of the
junction (Gloos  _et al._ , 2000). Henceforth, we make the bias argument
$\Phi_{j}^{\mathrm{ext}}$ implicit.

The energy functions of the $J$ non-linear dipoles and the structure and
topology of the $N_{\mathrm{cell}}+1$ cells are known. The circuit parameters
of the linear elements contained in $\mathcal{B}_{\mathrm{full}}$ can be
extracted using analytical results (Wolff, 2006) or numerical techniques, see
Sec. II.2. Information about the linear part of the circuit of the $n$-th cell
can be organized in the cell’s geometric capacitance and inverse inductance
matrices $\mathbf{C}_{n,\mathrm{cell}}^{\prime}$ and
${\mathbf{L}^{\prime}}_{n,\mathrm{cell}}^{-1}$, respectively. We define these
matrices with respect to the node-to-datum generalized magnetic fluxes of the
circuit (Devoret, 1995; Minev  _et al._ , 2020). In the case of a cell that
can be considered in the lumped regime,
$\mathbf{C}_{n,\mathrm{cell}}^{\prime}$ is the reduced Maxwell capacitance
matrix, see Eq. (15). In this case, $\mathbf{C}_{n,\mathrm{cell}}^{\prime}$
represents a fully connected graph, see Fig. 1(c).

We incorporate $C_{j}$ and $L_{j}$ of the non-linear dipoles into
$\mathbf{C}_{n,\mathrm{cell}}^{\prime}$ and
${\mathbf{L}^{\prime}}_{n,\mathrm{cell}}^{-1}$ to define the total cell linear
capacitance $\mathbf{C}_{n,\mathrm{cell}}$ and inverse inductance
$\mathbf{L}{}_{n,\mathrm{cell}}^{-1}$ matrices. In terms of these, the
capacitance and inverse inductance matrices of the composite-system are

$\mathbf{C_{n}}=\sum_{n=0}^{N_{\mathrm{cell}}}\mathbf{C}_{n,\mathrm{cell}}\quad\text{and}\quad\mathbf{L}_{\mathbf{n}}^{-1}=\sum_{n=0}^{N_{\mathrm{cell}}}\mathbf{L}_{n,\mathrm{cell}}^{-1}\;,$
(5)

respectively. They are real, symmetric, and nearly block diagonal. The matrix
$\mathbf{C_{n}}$ is guaranteed to be positive semi-definite, but due to the
inclusion of $L_{j}$, $\mathbf{L}_{\mathbf{n}}^{-1}$ does not have this
guarantee. Typically, $\mathbf{\mathbf{L}_{n}}^{-1}$ is sparse and rank
deficient. In terms of system, not cell, partitions, $\mathbf{C_{n}}$ and
$\mathbf{L}_{\mathbf{n}}^{-1}$ have a simple block structure. All subsystems
are described by blocks in the diagonals of $\mathbf{C_{n}}$ and
$\mathbf{L}_{\mathbf{n}}^{-1}$. The only off-diagonal blocks in the two
matrices are due to coupling elements spanning subsystems. Coupler nodes
$\mathcal{N}_{\mathrm{couple}}$ also introduce diagonal blocks in the two
matrices.

To describe the physical model, we introduce the column vector of all node-to-
datum fluxes $\mathbf{\Phi_{n}}$, each entry of which is associated with one
unique node of $\mathcal{N}_{\mathrm{full}}-\mathcal{N}_{g}$. The Lagrangian
of the composite system is

$\mathcal{L}_{\mathrm{n}}\left(\mathbf{\Phi_{n}},\mathbf{\dot{\Phi}_{n}}\right)=\frac{1}{2}\mathbf{\dot{\Phi}_{n}}^{\intercal}\mathbf{C_{n}}\mathbf{\dot{\Phi}_{n}}-\frac{1}{2}\mathbf{\Phi_{n}}^{\intercal}\mathbf{L_{n}}^{-1}\mathbf{\Phi_{n}}-\sum_{j=1}^{J}\mathcal{E}_{j}^{\mathrm{nl}}\left(\Phi_{j}\right)\;,$
(6)

where $J$ denotes the total number of non-linear dipoles in the network. The
Lagrangian is potentially singular due to the rank deficiency of
$\mathbf{L}_{\mathbf{n}}^{-1}$ and that of $\mathbf{C_{n}}$. The singularity
leads to dynamics on a constrained submanifold of phase space. In
quantization, it is this reduced, constrained phase space that provides the
required physical Poisson structure required for canonical quantization
(Dirac, 1982).

Before eliminating the coupler constraints on the phase space associated with
$\mathcal{L}_{\mathrm{n}}$, we first rotate the Lagrangian basis
$\mathbf{\Phi_{n}}$ to place all non-linear dipole fluxes explicitly in the
basis. The flux of the $j$-th dipole $\Phi_{j}=\Phi_{n_{2}}-\Phi_{n_{1}}$,
where $\Phi_{n_{2}}$ and $\Phi_{n_{1}}$ are the node-to-datum fluxes of its
two nodes. We construct the simple, linear transformation
$\mathbf{S_{n}}^{-1}$ with elements in $\left\\{-1,0,1\right\\}$ to rotate the
basis according to
$\mathbf{\Phi}=\mathbf{S_{n}}^{-1}\mathbf{\mathbf{\Phi_{n}}}$ such that
$\mathbf{\Phi}$ explicitly contains all junction fluxes $\Phi_{j}$. In the
basis $\mathbf{\Phi}$, the transformed capacitance and inverse inductance
matrices are
$\mathbf{C}=\mathbf{S_{n}}^{\intercal}\mathbf{C_{n}}\mathbf{S_{n}}$ and
$\mathbf{L}^{-1}=\mathbf{S_{n}}^{\intercal}\mathbf{L_{n}}^{-1}\mathbf{S_{n}}$,
respectively.

To eliminate constraints due to the singularity of $\mathbf{L}^{-1}$, we
select all $r$ coupler fluxes or linear combinations of fluxes in the kernel
space of $\mathbf{L}^{-1}$. For example, these include the fluxes associated
with all nodes in $\mathcal{N}_{\mathrm{couple}}$ that have only capacitive
elements touching them. We construct the $N\times r$ matrix formed by joining
the $r$ flux vectors $\mathbf{s}_{1},\mathbf{s}_{2},\ldots,\mathbf{s}_{r}$
defined in the $\mathbf{\Phi}$ basis as
$\mathbf{S}_{r}\coloneqq\left[\mathbf{s}_{1},\ldots,\mathbf{s}_{r}\right]$.
The flux-vector image space of $\mathbf{S}_{r}$ is a non-dynamical subspace of
the inductors; i.e., $\mathbf{L}^{-1}\mathbf{S}_{r}=\mathbf{0}_{Nr}$, where
$\mathbf{0}_{Nr}$ is the $N\times r$ matrix of all zeros. The image space of
$\mathbf{S}_{r}$ is not necessarily the full kernel space of
$\mathbf{L}^{-1}$; i.e.,
$\operatorname{span}\left(\mathbf{S}_{r}\right)\subseteq\ker\left(\mathbf{L}^{-1}\right)$.

We purposefully do not fully reduce the dynamics to only the image space of
$\mathbf{L}^{-1}$ at this stage. This allows the composite network model to be
reduced to known subsystems. These subsystems, treated as building blocks of
the device, can themselves incorporate a singular Lagrangian. A simple example
is that of the open-ended transmission line, which effectively has one less
inductor than capacitor and hence a kernel space of $\mathbf{L}^{-1}$ with
dimension one, which leads to one zero-frequency solution. To map to this
known problem, which itself handles its own singularity, we retain the
singularity in its definition by not including the associated subspace in the
span of $\mathbf{S}_{r}$.

We denote the complement of $\mathbf{S}_{r}$ as $\mathbf{S}_{k}$, an
$N\times\left(N-r\right)$ matrix whose image is the image space of
$\mathbf{L}^{-1}$ together with subspace of the kernel of $\mathbf{L}^{-1}$
that is not spanned by $\mathbf{S}_{r}$; i.e.,
$\mathbb{R}^{N}=\operatorname{span}\left(\mathbf{S}_{r}\right)+\operatorname{span}\left(\mathbf{S}_{k}\right)$.
Its columns can be selected mostly if not entirely from those of the identity
matrix.

We partition the degrees of freedom $\mathbf{\Phi}$ of the composite system
into the $r$ non-dynamical degrees of freedom
$\mathbf{\Phi_{r}}\coloneqq\mathbf{S}_{r}^{\intercal}\mathbf{\Phi}$ and the
remaining ones
$\mathbf{\Phi_{k}}\coloneqq\mathbf{S}_{k}^{\intercal}\mathbf{\Phi}$. For the
non-dynamical coordinates, we impose the standard zero initial condition
$\mathbf{\Phi_{r}}\left(t_{0}\right)=\left(\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{\Phi_{r}}\right)\left(t_{0}\right)=\mathbf{0}_{r}$,
where $t_{0}$ denotes the initial time of the circuit and $\mathbf{0}_{r}$ is
the length-$r$ column vector of all zeros. Using this, we solve for the $N$
degrees comprising $\mathbf{\Phi}$ in terms of the $N-r$ reduced dynamical
degrees $\mathbf{\Phi_{k}}$, and find
$\mathbf{\Phi}=\mathbf{S}_{k}\mathbf{\Phi_{k}}$. In the reduced basis, the
reduced capacitance and inverse inductance matrices are

$\displaystyle\mathbf{L_{k}}^{-1}$
$\displaystyle\coloneqq\mathbf{S}_{k}^{\intercal}\mathbf{L}^{-1}\mathbf{S}_{k}\;,$
(7a) $\displaystyle\mathbf{C_{k}}$
$\displaystyle\coloneqq\mathbf{S}_{k}^{\intercal}\left(\mathbf{C}-\mathbf{C}\mathbf{S}_{r}\left(\mathbf{S}_{r}^{\intercal}\mathbf{C}\mathbf{S}_{r}\right)^{-1}\mathbf{S}_{r}\mathbf{C}\right)\mathbf{S}_{k}\;,$
(7b)

respectively. The matrix $\mathbf{C_{k}}$ can be seen as a Schur complement.

In Eq. (7), we have eliminated the holonomic constraints presented by
inductive coupler nodes and contracted the problem from $N$ to $N-r$ degrees
of freedom $\mathbf{\Phi_{k}}$. The reduction of the coupler nodes
renormalizes system parameters subject to the transformation $\mathbf{S}_{k}$
and Eq. (7). In particular, the effective capacitance matrix can be strongly
dressed, due to the term containing
$\left(\mathbf{S}_{r}^{\intercal}\mathbf{C}\mathbf{S}_{r}\right)^{-1}$. This
is only a first dressing. The capacitance matrix can be dressed two more times
in the following.

####

The second renormalization step parallels the first, but this time with the
roles of the capacitance and inductance matrices flipped. Namely, if there are
coupler flux combinations in the kernel space
$\ker\left(\mathbf{C_{k}}\right)$, we repeat the above basis partition
protocol to eliminate these with respect to $\mathbf{C_{k}}$ (rather than
$\mathbf{L}^{-1}$, as above). An example of a coupler flux in
$\ker\left(\mathbf{C_{k}}\right)$ is the flux of a coupler node that is only
touched by inductors and no capacitors. Since the steps repeat, for brevity,
we do not explicitly detail the process again. Instead, we proceed by assuming
the step has been carried out, all coupler nodes in
$\ker\left(\mathbf{C_{k}}\right)$ and in $\ker\left(\mathbf{L}^{-1}\right)$
have been eliminated. To avoid new notation, we redefine $\mathbf{\Phi_{k}}$,
$\mathbf{L_{k}}^{-1}$ and $\mathbf{C_{k}}$ to denote in the following the
doubly reduced flux basis vector and circuit matrices and $N-r$ to denote the
length of $\mathbf{\Phi_{k}}$.

The generalized, reduced charge vector is
$\mathbf{Q_{k}}\coloneqq\mathbf{C_{k}}\mathbf{\Phi_{k}}$. We employ canonical
Dirac quantization (Dirac, 1982), explicated in App. C of Ref. Minev  _et al._
, 2020, with the canonical commutator
$\left[\mathbf{\hat{\Phi}_{k}},\mathbf{\hat{Q}_{k}}\right]=\mathbf{\hat{\Phi}_{k}}\mathbf{\hat{Q}_{k}}^{\intercal}-\mathbf{\hat{Q}_{k}}\mathbf{\hat{\Phi}_{k}}^{\intercal}=i\hbar\hat{I}\mathbf{I}$,
where $\mathbf{I}$ is the identity matrix of dimension $N-r$ and $\hat{I}$ is
the identity operator on the composite-system Hilbert space. The Hamiltonian
of the composite system is

$\hat{H}_{\mathrm{full}}=\frac{1}{2}\hat{\mathbf{Q}}_{\mathbf{k}}^{\intercal}\mathbf{C_{k}}^{-1}\hat{\mathbf{Q}}_{\mathbf{k}}+\frac{1}{2}\hat{\mathbf{\Phi}}_{\mathbf{k}}^{\intercal}\mathbf{L_{k}}^{\prime-1}\hat{\mathbf{\Phi}}_{\mathbf{k}}+\sum_{j=1}^{J}\mathcal{E}_{j}\left(\hat{\Phi}_{j}\right)\;,$
(8)

where $\mathbf{L_{k}}^{\prime-1}$ is the reduced effective capacitance matrix
$\mathbf{L_{k}}^{-1}$ with all non-linear dipole inductance $L_{j}$
contributions subtracted out. Since we included the non-linear dipole fluxes
$\Phi_{j}$ in the basis $\mathbf{\Phi}_{\mathbf{k}}$, this amounts to a simple
subtraction of the inductances from the diagonal of $\mathbf{L_{k}}^{-1}$. The
inversion of the composite-system capacitance matrix leads to the third
capacitive dressing of the subsystem parameters by the couplers and coupled
systems.

Equation (8) brings us to the desired composite-system Hamiltonian
$\hat{H}_{\mathrm{full}}$ defined by Eqs. (1) and (2). The $n$-th subsystem
Hamiltonian is a partition of the matrix equation given in Eq. (8),

$\hat{H}_{n}=\frac{1}{2}\hat{\mathbf{Q}}_{n}^{\intercal}\mathbf{C}_{n}^{-1}\hat{\mathbf{Q}}_{n}+\frac{1}{2}\hat{\mathbf{\Phi}}_{n}^{\intercal}\mathbf{L}_{n}^{\prime-1}\hat{\mathbf{\Phi}}_{n}+\sum_{j\in\mathcal{J}_{n}}\mathcal{E}_{j}\left(\hat{\Phi}_{j}\right)\;,$
(9)

where $\mathcal{J}_{n}$ is the set of non-linear dipoles belonging to the
$n$-th subsystem. The total number of junctions is
$J=\sum_{n=0}^{K}\left|\mathcal{J}_{n}\right|$. The vectors of subsystem
charge and flux operators are

$\hat{\mathbf{Q}}_{n}\coloneqq\operatorname{block}\left(\hat{\mathbf{Q}}_{\mathbf{k}},\mathcal{N}_{n}\right)\quad\text{and}\quad\hat{\mathbf{\Phi}}_{n}\coloneqq\operatorname{block}\left(\hat{\mathbf{\Phi}}_{\mathbf{k}},\mathcal{N}_{n}\right)\;,$
(10)

where the $\operatorname{block}$ function yields the matrix partition of its
first argument with respect to the set of nodes indicated in its second
argument. Similarly, the $n$-th subsystem matrices are the matrix partitions
with respect to the set of nodes $\mathcal{N}_{n}$ of the $n$-th system,

$\displaystyle\mathbf{C}_{n}^{-1}$
$\displaystyle\coloneqq\operatorname{block}\left(\mathbf{C_{k}}^{-1},\mathcal{N}_{n},\mathcal{N}_{n}\right)\;,$
(11a) $\displaystyle\mathbf{L}_{n}^{\prime-1}$
$\displaystyle\coloneqq\operatorname{block}\left(\mathbf{L}_{n}^{\prime-1},\mathcal{N}_{n},\mathcal{N}_{n}\right)\;,$
(11b)

where the second and third argument indicates the set of nodes associated with
the matrix rows and columns, respectively. The $K\left(K-1\right)/2$ pair-wise
interaction terms among the $K+1$ subsystems take the form

$\hat{H}_{n^{\prime}m^{\prime}}=\sum_{n\in\mathcal{N}_{n^{\prime}}}\sum_{m\in\mathcal{N}_{m^{\prime}}}\left(\hat{Q}_{n}\hat{Q}_{m}/\text{$C_{nm}^{\mathrm{eff}}$}+\hat{\Phi}_{n}\hat{\Phi}_{m}/\text{$L_{nm}^{\mathrm{eff}}$}\right)\;,$
(12)

where the indices $n^{\prime}$ and $m^{\prime}$ label the subsystems, for
$n^{\prime}\neq m^{\prime}$. The effective coupling capacitance
$C_{nm}^{\mathrm{eff}}$ and inductance $L_{nm}^{\mathrm{eff}}$ associated with
the $n$-th and $m$-th node of the $n^{\prime}$-th and $m^{\prime}$-th
subsystems, respectively, are obtained from the off-diagonal blocks of the
dressed system matrices,

$\displaystyle 1/\text{$C_{nm}^{\mathrm{eff}}$}$ $\displaystyle\coloneqq
2\operatorname{block}\left(\mathbf{C_{k}}^{-1},\mathcal{N}_{n^{\prime}}\left(n\right),\mathcal{N}_{m^{\prime}}\left(m\right)\right)\;,$
(13a) $\displaystyle 1/\text{\text{$L_{nm}^{\mathrm{eff}}$}}$
$\displaystyle\coloneqq
2\operatorname{block}\left(\mathbf{L_{k}}^{\prime-1},\mathcal{N}_{n^{\prime}}\left(n\right),\mathcal{N}_{m^{\prime}}\left(m\right)\right)\;,$
(13b)

where $\mathcal{N}_{n^{\prime}}\left(n\right)$ and
$\mathcal{N}_{m^{\prime}}\left(m\right)$ denote the single-element subsets of
$\mathcal{N}_{n^{\prime}}$ and $\mathcal{N}_{m^{\prime}}$ that comprise just
the $n$-th and $m$-th node of the sets, respectively.

By construction, we have assumed that the subsystem Hamiltonians $\hat{H}_{n}$
are of familiar, known structure and can be diagonalized by known analytical
or numerical methods. For instance, each Hamiltonian can be expressed in a
harmonic-oscillator basis and fully diagonalized, as detailed in the
methodology of the energy-participation-ratio method (Minev  _et al._ , 2020).
A diagonalization in this second quantized basis can provide the $n$-th system
Hamiltonian in the diagonal form

$\hat{H}_{n}=\sum_{\mu=0}^{M_{n}}\sum_{\nu=0}^{\infty}E_{n\mu\nu}\left|n,\mu,\nu\vphantom{n,\mu,\nu}\right\rangle\left\langle\vphantom{n,\mu,\nu}n,\mu,\nu\right|\;,$
(14)

where $E_{n\mu\nu}$ is the $\nu$-th Hamiltonian energy level of the $\mu$-th
resonant eigenmode of the $n$-th subsystem, which supports a total of $M_{n}$
number of quantized resonant modes. Since we have partitioned the composite
into small, manageable subsystems, for an effectively lumped subsystem of
interest one can diagonalize $\hat{H}_{n}$ independently from the circuit
using analytical or numerical techniques. In the case of numerics, one can use
a large, truncated basis for the diagonalization, but subsequently can retain
only a handful of the low-energy levels. The flux $\hat{\Phi}_{n}$ and charge
$\hat{Q}_{n}$ operators can be constructed in this truncated eigenbasis; i.e.,
in terms of the $\left|n,\mu_{n},\nu_{n}\right\rangle$ and
$\left|m,\mu_{m},\nu_{m}\right\rangle$ eigenstates. Thus, the interaction
$\hat{H}_{nm}$ and composite Hamiltonian $\hat{H}_{\mathrm{full}}$ can be
expressed more efficiently 222This is the form employed in a recent package
entitled scqubits; see also Ref. Kerman, 2020..

In this section, we explicitly constructed the subsystem Hamiltonians and
their interactions and eliminated all coupler degrees of freedom and
associated constraints, and incorporated all circuit effects, such as the
$A^{2}$ diamagnetic contribution, in Eqs. (8)–(13). The parameters of each
subsystem Hamiltonian $\hat{H}_{n}$ were obtained by three sequential
dressings. Thus, the spectrum of $\hat{H}_{n}$ and $\hat{H}_{m}$ are not
circuit parameter independent in general. They can be dressed by their
couplers and by each other’s circuit parameters. The dressing can in principle
be rather non-local. This depends on the particulars of the circuits. In
typical situations, nodes in the network graph separated by a minimal path of
several edges can be coupled but potentially negligibly so. For example, in
the inversion step of $\mathbf{C_{k}}$, these distance direct couplings are
allowed but severely dampened by the weight of the inversion. However,
dressing of node parameters by proximal coupler edges can be large. For
example, for the direct capacitance between the pads of a transmon qubit can
be renormalized by a factor of two subject to a large pad coupler, such as
featured in our experiment, see Sec. IV. The coupling terms are also similarly
dressed by the system parameters. Proper elimination at the composite level,
rather than at the individual-system level, is essential for more general
circuits and accurate results, as observed from data for the theory versus
experiment comparison presented in Sec. IV.

###  II.1 Example of the experimentally measured devices

Let us briefly illustrate the method by returning to the example of Fig. 1.
The composite system is comprised of the qubit
($\mathcal{N}_{0}=\left\\{P_{1}\right\\}$), CPW readout resonator
($\mathcal{N}_{1}=\left\\{B_{1},\ldots\right\\}$), two CPW bus resonators
($\mathcal{N}_{2}=\left\\{B_{2},\ldots\right\\}$ and
$\mathcal{N}_{3}=\left\\{B_{3},\ldots\right\\}$), and coupler nodes
($\mathcal{N}_{\mathrm{coupler}}=\left\\{P_{0},C\right\\}$); the ellipses in
the sets indicate the continuum of nodes associated with the distributed
lines. Here, $K=3$; however, we only need one lumped simulation cell. This
cell contains the set of nodes
$\mathcal{N}_{0}^{\mathrm{cell}}=\left\\{P_{0},P_{1},B_{0},B_{1},B_{2},C\right\\}$;
the ground node is considered accessible by every cell. The remaining three
cells associated with the bodies of the CPW transmission lines need not be
simulated explicitly since they can be handled quasi-analytically, see Sec.
III. Following the procedure outlined in Sec. II.2, the cell capacitance
matrix $\mathbf{C}_{0,\mathrm{cell}}$ is extracted; its fully-connected graph
is depicted with green lines in Fig. 1(c). For our experiment, we employed the
_Ansys Q3D Extractor_ quasi-field solver. The creation of the qubit cell and
extraction of its Maxwell matrix was automated using Qiskit Metal (Minev  _et
al._ , 2021).

The inductance matrix of the qubit cell contains a single inductor due to the
one junction in the cell; i.e., $J=1$,
$\mathcal{J}_{0}=\left\\{\Phi_{J}\right\\}$ and $\mathcal{J}_{n}=\emptyset$
for $n>0$, where the junction flux is $\Phi_{J}=\Phi_{P_{1}}-\Phi_{P_{0}}$.
The CPW cells all have diagonal capacitance matrices
$\mathbf{C}_{n,\mathrm{cell}}$, where $n\in\left\\{1,2,3\right\\}$; their
inductance matrices $\mathbf{L}_{n,\mathrm{cell}}^{-1}$ each have a one-
dimensional kernel space. The qubit cell inductance matrix in isolation has
$6$ columns but has unity rank. Since the three bus nodes are touched by
inductors in other systems, these three nodes will not be eliminated; i.e.,
one should consider the kernel space of the composite inductance matrix
$\mathbf{L}_{\mathbf{n}}^{-1}$ and not just that of a given cell
$\mathbf{L}_{n,\mathrm{cell}}^{-1}$, see Eq. (5).

Since there are no coupling inductors, we can obtain the reduced Hamiltonian
matrices $\mathbf{C_{k}}^{-1}$ and $\mathbf{L_{k}}^{-1}$ directly from Eq.
(7), without having to perform a second elimination of inductive coupling
nodes. Since the transmission line systems have diagonal node-to-datum cell
matrices, it suffices to only invert and reduce $\mathbf{C}_{0,\mathrm{cell}}$
to obtain all needed information comprised in $\mathbf{C_{k}}^{-1}$. The
reduced effective circuit is depicted in Fig. 1(d); the third bus is omitted
from the illustration for the sake of visual simplicity. The dressed transmon
qubit Hamiltonian, associated circuit shaded in green in Fig. 1(d), is

$\hat{H}_{0}=\frac{1}{2}\left(\hat{Q}_{J}-Q_{\mathrm{ofs}}\right)^{2}/C_{J}^{\mathrm{eff}}-E_{J}\cos\left(\hat{\Phi}_{J}\right)\;,$

where we have included a potential charge offset $Q_{\mathrm{ofs}}$, which
could be due to a charge fluctuation on any of the nodes in
$\mathcal{N}_{0}^{\mathrm{cell}}$, where $\hat{Q}_{J}\coloneqq
C_{J}^{\mathrm{eff}}\frac{\mathrm{d}}{\mathrm{d}t}\hat{\Phi}_{J}$, and
$1/C_{J}^{\mathrm{eff}}\coloneqq\operatorname{block}\left(\mathbf{C_{k}}^{-1},\left\\{P_{1}\right\\},\left\\{P_{1}\right\\}\right)$.
Importantly, $C_{J}^{\mathrm{eff}}$ comprises the dressed and renormalization
effect of all capacitances in the qubit cell model; i.e., while we may refer
to it as the dressed qubit capacitance, its value depends on the coupling and
transmission-line end-loading capacitances. That is, the qubit and the
neighboring system Hamiltonians are dressed by each others’ elements and their
parameters are not strictly independent. This first dressing is purely
classical and due to linear effects. The dressed loading capacitance of the
$n$-th line, see Eq. (16), is

$1/C_{L,n}^{\mathrm{eff}}\coloneqq\operatorname{block}\left(\mathbf{C_{k}}^{-1},\left\\{B_{n}\right\\},\left\\{B_{n}\right\\}\right)\;.$

The interaction Hamiltonian between the $n$-th and $m$-th subsystems is

$\hat{H}_{nm}=\hat{Q}_{n}\hat{Q}_{m}/\text{$C_{nm}^{\mathrm{eff}}$}\;,$

where $\hat{Q}_{0}=\hat{Q}_{J}$ and $C_{nm}^{\mathrm{eff}}$ is the
appropriately scaled off-diagonal term of the reduced and inverted
$\mathbf{C}_{0,\mathrm{cell}}$, see Eq. (13). In Sec. IV, we compare the
parameters of interest found from this method to those experimentally measured
on a large-scale quantum processor (Antonio D.  _et al._ , 2021).

###  II.2 Lumped cell model

In the physical model of a lumped cell, see for example Fig. 1(b), each
galvanically-isolated conducting island is represented by a node in the
circuit. The effect of the geometry and all other physical aspects of the
cell, such as dielectric permittivities and material properties, on the
electrostatics of the cell are succinctly expressed in the Maxwell capacitance
matrix(Zangwill, 2012) $\mathbf{C_{M}}$, which is reviewed in App. A. The
$\left(i,j\right)$-th off-diagonal element of the Maxwell matrix is the mutual
capacitance between the $i$-th and $j$-th node in the cell times negative
unity. It is efficiently extracted from the cell model using a quasi-field
solver, such as _Ansys Q3D Extractor_ ; this process __ is automated __ in our
open-source project Qiskit Metal (Note, 1)_._ The cell capacitance matrix with
respect to the cell node-to-datum generalized magnetic fluxes, see Eqs. (5)
and (A1), is

$\mathbf{C}_{\mathrm{cell}}=\mathbf{S}_{\mathbf{M}}^{\intercal}\mathbf{C_{M}}\mathbf{S_{M}}\;,$
(15)

where
$\mathbf{S_{M}}\coloneqq\left[\begin{array}[]{cc}\mathbf{0}_{N}&\mathbf{I}_{N}\end{array}\right]^{\intercal}$,
$\mathbf{0}_{N}$ is the zero column vector of length $N$, and $\mathbf{I}_{N}$
is the identity matrix of dimension $N$, where $N=\dim\mathbf{C_{M}}-1$.

##  III Quantum physics of the loaded transmission-line cell

####

Figure 2: Transmission line capacitively loaded at one end. (a) Continuous
model of a line with total length $L$ and loading capacitance $C_{L}$. Length
along the line is $z$. For $b=0$ (resp., $b=1$), the right end of the line is
terminated in an open (resp., short). (b) Circuit schematic of an equivalent,
discrete-space model, with $N$ nodes. The spatial discretization length is
$\Delta z\coloneqq L/(N-1)$. The capacitance and inductance per unit length of
the line are $c$ and $l$, respectively. The operators of the node-to-datum
generalized magnetic fluxes of the circuit are $\hat{\Phi}_{n}$, for
$n=0,1,\ldots N-1$.

Because of its practical and central importance, we revisit the
diagonalization of an end-loaded transmission line, depicted in Fig. 2(a).
Quantizing its discrete model and taking the continuous limit, one finds its
quantum Hamiltonian

$\hat{H}_{\mathrm{LTL}}=\frac{1}{2}C_{L}^{-1}\left[\hat{Q}(0)\right]^{2}+\int_{0^{+}}^{L}\mathrm{d}z\,c^{-1}\left[\hat{q}\left(z\right)\right]^{2}\\\
+\int_{0}^{L}\mathrm{d}z\,l^{-1}\left[\frac{\partial\hat{\Phi}\left(z\right)}{\partial
z}\right]^{2}\;,$ (16)

where $C_{L}$ is the effective end-loading capacitance to ground at $z=0$, the
total charge operator at $z=0$ is $\hat{Q}(0)$, and
$L,c,l,\hat{q}\left(z\right),$ and $\hat{\Phi}\left(z\right)$ are the line
total length, capacitance and inductance per unit length, and the charge-
density and magnetic-flux-field operators at length $z$ along the line,
respectively. Due to the end-singularity in the capacitance per unit length of
the line, strictly speaking, there is no operator $\hat{q}\left(0\right)$;
hence, the integration from $0^{+}$ in the lower integral bound, which
excludes zero. The scalar, one-dimensional quantum fields obey
$\left[\hat{\Phi}(z),\hat{q}(z^{\prime})\right]=i\hbar\delta\left(z-z^{\prime}\right)\hat{I}$
and $\left[\frac{\partial\hat{\Phi}\left(z\right)}{\partial
z},\hat{q}(z^{\prime})\right]=i\hbar\frac{\partial}{\partial
z}\delta\left(z-z^{\prime}\right)\hat{I}$, where $\hat{I}$ is the identity
operator.

#### Boundary conditions.

For an open (short) right-end termination, the right boundary condition is of
the standard Dirichlet (respectively, Neumann) type,
$\frac{\partial\hat{\Phi}\left(L\right)}{\partial z}=\hat{0}$ (resp.,
$\hat{\Phi}\left(L\right)=\hat{0}$). However, the left-end boundary condition
is non-standard. It is not of the type covered by standard Sturm–Liouville
theory, $\frac{\partial\hat{\Phi}\left(0\right)}{\partial
z}=lC_{L}\frac{\partial^{2}\hat{\Phi}\left(0\right)}{\partial t^{2}}$. This
condition makes the boundary condition eigenvalue-dependent leading to
transcendental eigensolutions. The equation of motion in the interior of the
line is the standard wave equation,
$v_{p}^{2}\frac{\partial^{2}\hat{\Phi}}{\partial
z^{2}}\left(z\right)=\frac{\partial^{2}\hat{\Phi}}{\partial
t^{2}}\left(z\right),$ valid on $0<z<L$, where the phase velocity is
$v_{p}=\sqrt{1/\left(lc\right)}$.

#### Characteristic eigenvalue equation.

The Hamiltonian is diagonalized using harmonic solutions and the superposition
principle for linear systems. It is easier to show by demoting the problem to
classical fields, then recovering the quantum field operators and their
quantum zero-point fluctuations using the energy-participation-ratio method
(Minev  _et al._ , 2020); this establishes the classical-quantum
correspondence of the problem and avoids technical detail regarding the proper
orthonormalization of the field eigensolution in the presence of the loading
irregularity. Using the harmonic ansatz
$\Phi_{m}\left(z,t\right)=u_{m}\left(z\right)\Phi_{m}\left(t\right)=A_{m}\sin\left(\omega_{m}t\right)\cos\left(k_{m}z+\phi_{m}\right)$
and eliminating $k_{m}$ and $\phi_{m}$, we find the characteristic eigenvalue
equation for $\omega_{m}$,

$\boxed{\omega_{m}\frac{L}{v_{p}}+\arctan\left(\frac{\omega_{m}}{\omega_{L}}\right)=m\pi+b\frac{\pi}{2}:m\in\mathbb{Z}_{\geq
0}\;,}$ (17)

where $b=0$ (resp., $b=1)$ for an open (resp., short) right-end boundary
condition and $\omega_{L}\coloneqq\frac{1}{lC_{L}v_{p}}=\frac{1}{C_{L}Z_{0}}$.
The non-negative integer $m$ labels the discrete, uncountable number of modes.
The associated spatial phase shift and wavenumber are $k_{m}=\omega_{m}/v_{p}$
and $\phi_{m}=\arctan\left(\omega_{m}/\omega_{L}\right)$, respectively.

The equation is transcendental. For a loading knee frequency
$\omega_{L}\gg\omega_{m}$ (resp., $\omega_{L}\ll\omega_{m}$) the loading
capacitor acts as an open (resp., short), while for $\omega_{m}$ of the same
order of magnitude as $\omega_{L}$, the loading exerts a significant
renormalization on the eigenfrequency and eigenfields away from those
obtainable from the unloaded transmission line solutions; see Fig. 3. In other
words, in this regime, the naive eigenfunctions of an unloaded transmission
line are not a good starting-point basis choice, since the phase shift
$\phi_{m}$ and renormalization of the mode frequencies can be large. It is
this latter regime that is experimentally relevant for the devices and
measurements presented in Sec. IV. For our devices, the renormalization of the
loading was on the order of 20–25%.

Figure 3: Dressing of the eigenfrequencies (a) and eigenfields (b) of a
transmission line resonator capacitively loaded at its left end. Line and
loading parameters correspond to those of measured devices. (a) Solutions
$\omega_{m}$ of the transcendental eigenvalue equation are found at the
intersection of the diagonal (red) line with the curved (resp., horizontal)
line for the loaded (resp., unloaded) line, see Eq. (17). The unloaded line
corresponds to a loading capacitance $C_{L}=0$. The point labeled $m=1$
denotes the fundamental loaded solution. (b) Corresponding eigenfield spatial
distributions $u_{m}\left(z\right)$, depicted for the first three modes,
indexed by $m$. Vertical dashed line is a guide to the eye, highlighting the
symmetry breaking due to the asymmetric pull on the fields toward the loading.

#### Quantum Hamiltonian.

Exploiting the eigensolutions $\omega_{m}$ of Eq. (17), in second
quantization, Eq. (16) becomes

$\hat{H}_{\mathrm{LTL}}=\sum_{m=m_{0}}^{\infty}\hbar\omega_{m}\left(\hat{a}_{m}^{\dagger}\hat{a}_{m}+\frac{1}{2}\right)\;,$

where $\hat{a}_{m}$ is the annihilation operator for the $m$-th mode. For the
case of an open termination $b=0$, a trivial zero-frequency (d.c.) mode
solution of Eq. (17) exists, and the lower bound of the sum is $m_{0}=1$;
otherwise, for $b=1$, $m_{0}=0$. By linearity, the fields are linear
superpositions of the modal operators,

$\displaystyle\hat{\Phi}\left(z\right)$
$\displaystyle=\sum_{m=m_{0}}^{\infty}\Phi_{m}^{\mathrm{ZPF}}\left(z\right)\left(\hat{a}_{m}^{\dagger}+\hat{a}_{m}\right)\;,$
$\displaystyle\hat{q}\left(z\right)$
$\displaystyle=\sum_{m=m_{0}}^{\infty}iq_{m}^{\mathrm{ZPF}}\left(z\right)\left(\hat{a}_{m}^{\dagger}-\hat{a}_{m}\right)\;,$
$\displaystyle\hat{Q}\left(0\right)$
$\displaystyle=\sum_{m=m_{0}}^{\infty}iQ_{m}^{\mathrm{ZPF}}\left(0\right)\left(\hat{a}_{m}^{\dagger}-\hat{a}_{m}\right)\;,$
(18a)

where we will use the EPR to find the values of the quantum zero-point
fluctuations $\Phi_{m}^{\mathrm{ZPF}}\left(z\right)$,
$q_{m}^{\mathrm{ZPF}}\left(z\right)$, and $Q_{m}^{\mathrm{ZPF}}\left(0\right)$
of the magnetic flux, charge density, and charge, respectively.

#### Quantizing the fields using the energy-participation ratio (EPR).

The EPR of the loading capacitor $p_{mL}$ in mode $m$ is the fraction of
energy stored in the loading capacitor $\mathcal{E}_{C_{L}}$ relative to the
total capacitive energy of the mode $\mathcal{E}_{\mathrm{cap}}$. As a
fraction, the participation is independent of the normalization of $u_{m}$; in
the classical setting,

$p_{mL}\coloneqq\frac{\mathcal{E}_{C_{L}}}{\mathcal{E}_{\mathrm{cap}}}=\frac{\frac{1}{2}C_{L}\left[u_{m}\left(0\right)\right]^{2}}{\frac{1}{2}C_{L}\left[u_{m}\left(0\right)\right]^{2}+\frac{1}{2}\int_{0}^{L}\mathrm{d}z\,c\left[u_{m}\left(z\right)\right]^{2}}\;,$
(19)

which is evaluated using the solutions of Eq. (17).

In the quantum setting,
$\left\langle\mathcal{E}_{C_{L}}\right\rangle=\left\langle\frac{1}{2C_{L}}\left[\hat{Q}\left(0\right)\right]^{2}\right\rangle$,
where we take the expectation value over a single photon state in mode $m$ and
disregarding zero-point energy contributions; substituting Eq. (18a),
$\left\langle\mathcal{E}_{C_{L}}\right\rangle=\frac{1}{C_{L}}Q_{\mathrm{ZPF}}^{2}$.
The total capacitive energy is
$\left\langle\hat{\mathcal{E}}_{\mathrm{cap}}\right\rangle=\frac{1}{2}\hbar\omega_{m}$;
hence, the quantum ZPF of the field charge operator at $z=0$ can be found in
terms of the EPR, which we can calculate from the classical solutions, using
Eq. (19),

$\left[Q_{m}^{\mathrm{ZPF}}\left(0\right)\right]^{2}=\text{$\frac{\hbar\omega_{m}}{2}C_{L}p_{mL}$}\;.$
(20)

Using the same line of EPR reasoning, one finds

$\displaystyle\left[q_{m}^{\mathrm{ZPF}}\left(z\right)\right]^{2}$
$\displaystyle=\frac{\hbar\omega_{m}}{2}cp_{mc}\left(z\right)\;,$ (21)
$\displaystyle\Phi_{m}^{\mathrm{ZPF}}\left(z\right)$
$\displaystyle=\frac{1}{c\omega_{m}}q_{m}^{\mathrm{ZPF}}\left(z\right)\;,$
(22)

where the EPR density
$p_{mc}\left(z\right)\coloneqq\frac{1}{2}c\left[u_{m}\left(z\right)\right]^{2}/\mathcal{E}_{\mathrm{cap}}$
is the density of the fraction of capacitive energy stored in the
infinitesimal capacitance at $z$. Equations (20) and (21) fully specify the
quantum fields, Eq. (18), and thus complete the solution of the loaded line.

##  IV Measurements and comparison between theory and experiment

Figure 4: Measured vs. predicted qubit-readout cross-Kerr coupling $\chi_{qr}$ for 10 qubit subsystems across two different 14-qubit processors. The parameter $\chi_{qr}$ is the most sensitive and dressed parameter of the composite-system Hamiltonian. The frequency of the readout CPW resonator is dressed down by 25% due to its embedding in the larger network. Circles: predictions from the full method of this paper. Triangles: predictions from a simpler model that performs an approximate capacitive reduction based on weak coupling and does not account for the dressing of the zero-point quantum fluctuations of the CPWs. The agreement between the experiment and the full (resp., naive) theory is -10.5% (resp., +19%). The simulation error bars are smaller than the size of the dots. Model feature / Parameter |  Magnitude of effect on experimental agreement
---|---
Including all coupling Hamiltonians $\hat{H}_{nm}$ |  $5\%$
Including qubit coupling to all bus resonators $\hat{H}_{qm}$ |  $1\%$
Transmission-line impedance $Z_{0}$ |  $1\%$ for $3\%$ variation on $Z_{0}=50\,\Omega$
Chip separation |  $1\%$ for $20\%$ variation on separation
Substrate permittivity $\epsilon_{r}$ |  $0.5\%$ for $2\%$ variation on $\epsilon_{r}$
Including readout first harmonic in analysis $\omega_{12}$ |  $0.3\%$
Cell partition bounding box length |  Less than $0.5\%$ and increasingly negligible over 100 $\mathrm{\mu m}$
Bus resonator frequency $\omega_{n1}$ |  Negligible for $\pm 5\%$ variation on $\omega_{n1}$
Cell bounding box padding distance |  Negligible over 100 $\mathrm{\mu m}$
Transmission-line phase velocity $v_{p}$ |  Negligible
Including charge line node in reduced model |  Negligible
Sample holder enclosure distance |  Negligible
Table 1: Influence and error budget for the experimental agreement due to
model features and parameters. Left column: model feature or parameter. Right
column: magnitude of the change in the average experimental agreement reported
in Fig. 4 for the dispersive coupling $\chi_{qr}$. The intrinsic capacitance
$C_{j}$ of a Josephson junction in the model is fixed, and its inductance
$L_{J}$ is determined based on the measured qubit frequency.

We applied the lumped-oscillator-model (LOM) method, presented in this work,
on two 14-qubit superconducting quantum processors, measured over multiple
cooldowns(Antonio D.  _et al._ , 2021). Ten of the most well measured
subsystems were chosen for analysis with the LOM method. The processor
architecture was based on floating transmon qubits (Koch  _et al._ , 2007)
interconnected by CPW bus resonators (Blais  _et al._ , 2004), controlled by
charge lines, and readout by CPW resonators. Qubit, bus, and resonator
frequencies were allocated in the 5.1–5.5 GHz, 6.5–6.7 GHz, and 7–7.1 GHz
bands, respectively. Qubit anharmonicities were measured in the 300–350 MHz
range. Readout coupling was strong so as to provide fast qubit readout
(Antonio D.  _et al._ , 2021); qubit-cavity dispersive shifts $\chi_{qr}$ were
in the 3–7 MHz range. Each qubit had either 3 or 4 couplers—one coupler for
the readout CPW resonator and the remainder for bus CPW resonators; see Fig.
1(b). The qubit connectivity implemented the topology of 3 lattice placates.
FPGAs controlled the experiment. Single-qubit gates were Gaussian-shaped
pulses of length 30 ns. The measurement readout time was 320 ns, performed in
reflection, through a Traveling Wave Parametric Amplifier (TWPA). The
experimental is further described in Antonio D.  _et al._ , 2021.

We measured the qubit and cavity system parameters—the qubit frequency
$\omega_{q}$ and anharmonicity $\alpha_{q}$ and the readout frequency
$\omega_{r}$—using standard spectroscopic and time-resolved protocols. Of
note, we used a more-sensitive Ramsey dressed-dephasing-based method (Gambetta
_et al._ , 2006, 2008; Antonio D.  _et al._ , 2021) to measure their
interaction $\chi_{qr}$. This much smaller parameter $\chi_{qr}$ was extracted
from dephasing experiment measurements simultaneously fitting the dephasing in
both the $\left\langle X\right\rangle$ and $\left\langle Y\right\rangle$
quadratures of the qubit Bloch vector as a function of the readout probe
frequency and the effective readout-strength photon number $\bar{n}$. This
collection of curves, which are very sensitive to the Hamiltonian parameters,
was simultaneously fit to extract the readout photon number $\bar{n}$ and
dispersive shift $\chi_{qr}$. We also accounted for the rise and fall times of
the readout pulse in the extraction by accounting for the time-domain
evolution of the reduced stochastic master equation of the qubit (Gambetta
_et al._ , 2008).

In the following, we compare the measured Hamiltonian parameters to those
obtained using the LOM method applied to the physical layout of the devices.
We focus on the most challenging parameter to obtain quantitative agreement
for, $\chi_{qr}$.

We model the qubit cell as depicted in Fig. 1(b) by including; a short segment
of charge control line, neighboring CPW structures, and all of the coupler
structures attached to the qubit in the cell model. We extract the qubit cell
Maxwell capacitance matrix $\mathbf{C_{M}}$ using _Ansys Q3D Extractor_ by
using the layout geometry and a nominal substrate relative permittivity
(Krupka  _et al._ , 2006) $\epsilon_{r}=11.45$. We estimate the nominal
Josephson junction intrinsic capacitance $C_{J}=2\,\mathrm{fF}$. The matrix is
treated and reduced as described in Sec. I, after which we find a notably
large capacitance to ground $C_{L}$, in excess of $320\,\mathrm{fF}$, loading
the readout coupling node $P_{0}$. The effective CPW lengths are extracted
from Eq. (17), where impedance $Z_{0}=53\,\Omega$ and phase velocity
$v_{p}=0.403c$, where $c$ is the speed of light, of the lines were found
analytically from the line geometry (Simons, 2001). Finally, the qubit
Josephson tunnel junction inductance $L_{j}$ is varied from its nominal
values, inferred from room-temperature resistance measurements (Gloos  _et
al._ , 2000), due to its unavoidable and inherent aging and cool-down
variability, to obtain a qubit frequency agreement at the 1% level.

While some parameters of the model are known with exceedingly high precision,
such as the geometry of the qubit device, we allow for small, reasonable
variation in several parameters that are less well known. The variation in
$\chi_{qr}$ caused by the roughly 5% uncertainty inherent in these parameters,
including $L_{J},C_{J},Z_{0}$ and $\epsilon_{r}$, is detailed in Table 1.
Furthermore, for achieving the highest level of agreement, we detail in Table
1 the relative importance of including finer-order effects in the Hamiltonian
model, such as the coupling of the qubit to the readout CPW first-harmonic
mode. While this particular contribution has only a single-percent level
effect on $\chi_{qr}$, including all direct CPW-CPW couplings originating from
the direct capacitive links in the qubit cell, described in the Hamiltonian as
terms of the form
$\beta_{nm}\hat{a}_{n}^{\dagger}\hat{a}_{n}\hat{a}_{m}^{\dagger}\hat{a}_{m}$,
significantly decreases the cross-Kerr $\chi_{qr}$ by approximately 5%. In
practice, since the qubit cell is only strongly coupled to a small handful of
modes, we numerically diagonalize each individual system Hamiltonian
$\hat{H}_{n}$ first, constructing the composite Hamiltonian
$\hat{H}_{\mathrm{full}}$, using Eqs. (1) and (2) and diagonalizing the full
system, from which we extract all final model parameters.

Accounting for these finer effects, we present the comparison between
experimentally measured and LOM-predicted values of $\chi_{qr}$ in Fig. 4.
Across the 10 devices, we find that the agreement between theory and
experiment is -10.5%. We compare this to a naive model, which disregards these
finer effects and, chiefly, does not account for the large renormalization of
the transmission line eigenmodes as detailed in Sec. III, illustrated in Fig.
3. Under these more conventional approximations, the average theory-experiment
agreement is +19%. A chief component in the difference is the large magnitude
of $C_{L}$, which dresses down the CPW readout mode frequency from
approximately 8.8 GHz to 7.0 GHz, a 25% renormalization effect.

In Table 1, we summarize the fluctuation in the agreement between simulated
and measured $\chi_{qr}$ from the modification of model features and
parameters . For example, the table reports the effect of varying the location
of the qubit cell partition and the length of neighboring CPWs included in the
qubit cell model. This is found to have a nearly negligible effect.
Determination of the consistently correct amount to include was not reached,
though for the transmon layouts this paper considered,
$\approx\leavevmode\nobreak\ 100\mu$m proved to be sound. The separation of
the device chips has a small effect when the gap is larger than the feature
dimensions (Minev  _et al._ , 2016), though becomes significant when said gap
is equivalent to or less than feature dimensions. Including the sample holder
packaging in the model has a negligible effect for the dimensions of our
device. It is not possible to achieve perfect experimental agreement by
varying the model parameters within experimental tolerances. In fact, the
model is fairly well constrained. We believe that better agreement requires a
more complete model description of distributed effects. These can be captured
by more computationally-expensive, but also more informationally-complete
treatments, such as impedance (Nigg  _et al._ , 2012) or energy-participation
ratio (EPR) quantization (Minev  _et al._ , 2020).

#### Conclusion.

We introduced the lumped-oscillator method (LOM) and studied its experimental
performance on state-of-the-art superconducting quantum processors (Antonio D.
_et al._ , 2021). We found the LOM method to be practical, due to its
systematic and modular analysis flow. It operates by partitioning the device
in a two-fold manner—subsystems and simulation cells. Cells are small domains
of the physical device layout, which can be independently analyzed using
analytical results or numerical simulations. Breaking up the device layout
into independent cells increases the computational modeling efficiency.
Subsystems serve as small, familiar building blocks. Within the lumped
approximation, the full Hamiltonian $\hat{H}_{\mathrm{full}}$ of the composite
system is extracted with no approximations on the strength or nature of the
inductive non-linear dipoles. Subsystem-subsystem couplings can be arbitrarily
large. This was experimentally important to account for the data. We observed
subsystem renormalizations on the order of 25% in the experiment due to
dressings induced by the coupler structures and the subsystem embedding in the
larger network. Overall, we found a two-times improvement in experimental
agreement using the LOM method over one that makes weak coupling
approximations and does not fully account for the non-perturbative dressing of
the distributed modes.

With the increased complexity and demands driven by the rapid development of quantum hardware, we believe fast, accurate, and systematic techniques such as the LOM method presented here will be an essential ingredient in the development of current and future quantum technology. For this reason, we contribute the automation and implementation of this method to the community, as part of our open-source project Qiskit Metal | for quantum device design (Note, 1).

#### Acknowledgements.

We thank M. Malekakhlagh, F. Solgun, D.C. McKay, D. Wang, and R. Gutiérrez-
Jáuregui for valuable discussions. We are grateful to all the early-access
participants of Qiskit Metal and the Metal team for discussions and stress
testing the open-source code. We thank G. Calusine and W. Oliver for providing
the traveling-wave parametric amplifier used in this work. We acknowledge
partial support for work on the simulations, and experimental bring-up and
characterization of the devices from the Intelligence Advanced Research
Projects Activity (IARPA) under Contract No. W911NF-16-1-0114.

#### Author contributions.

Z.K.M developed the theory and software, and wrote the manuscript. Z.K.M.
conceived of the project based on work by J.M.G., who oversaw the work. Z.K.M.
and T.G.M. simulated and analyzed the devices. M.T. and A.C. performed the
experimental bring-up and characterization of the devices and worked with
Z.K.M on the analysis of the experimental data. All authors discussed the
results and contributed to the manuscript.

#### Code availability.

The source code for Qiskit Metal | for quantum device design and for the LOM method implementation is open-sourced and can be found at github.com/qiskit/qiskit-metal.

##  Appendix A Maxwell capacitance matrix

The effective capacitances of a cell can be extracted from its physical
layout. A simulation incorporating the cell geometry, materials, and
electromagnetic boundary conditions can yield its Maxwell capacitance matrix
$\mathbf{C_{M}}$. For a cell with $N$ nodes (nets), $\mathbf{C_{M}}$ is an
$N+1$ square, full-rank matrix. Its off-diagonal entries are the negative of
the capacitance $C_{ij}$ between nodes $i$ and $j$, where $C_{ij}\geq 0$. Its
$i$-th diagonal is the self-capacitance of the $i$-th node to infinity
$C_{ii}$ plus the value of all coupled capacitances; i.e., the $i$-th row and
$j$-th column element of $\mathbf{C_{M}}$ is

$\left[\mathbf{C_{M}}\right]_{i,j}=\begin{cases}-C_{ij}\;,&i\neq j\;,\\\
\sum_{j^{\prime}=1}^{N}C_{ij^{\prime}}&i=j\;.\end{cases}$ (A1)

For example, the Maxwell matrix for the set of nodes
$\mathcal{N}=\left\\{n_{0},\ldots n_{5}\right\\}$ has the form

$\mathbf{C_{M}}={\NiceMatrixOptions{first-row,first-col,code-for-first-
col=\color[rgb]{.5,.5,.5},code-for-first-
row=\color[rgb]{.5,.5,.5}}\begin{pmatrix}&n_{0}&n_{1}&n_{2}&n_{3}&n_{4}&n_{5}\\\
n_{0}&C_{0\Sigma}&-C_{01}&-C_{02}&-C_{03}&-C_{04}&-C_{05}\\\
n_{1}&&C_{1\Sigma}&-C_{12}&-C_{13}&-C_{14}&-C_{15}\\\
n_{2}&&&C_{2\Sigma}&-C_{23}&-C_{24}&-C_{25}\\\
n_{3}&&&&C_{3\Sigma}&-C_{34}&-C_{35}\\\ n_{4}&&&&&C_{4\Sigma}&-C_{45}\\\
n_{5}&&&&&&C_{5\Sigma}\end{pmatrix}},$

where we introduce the shorthand $C_{i\Sigma}\coloneqq\sum_{j=1}^{N}C_{ij}$
and, for simplicity, we omit the lower-triangular block of the symmetric
matrix. The sum of the $n$-th row and identically the $n$-th column of
$\mathbf{C_{M}}$ is the self-capacitance (to infinity) $C_{nn}$ of the $n$-th
node. An abstract conductor at infinity serves as the effective datum of the
schematic corresponding to the Maxwell matrix. The matrix is expressed in the
basis of node fluxes $\mathbf{\Phi}_{\infty}$ referenced to that infinity
conductor, assumed at zero potential; note,
$\dim\mathbf{\Phi}_{\infty}=\left|\mathcal{N}\right|\times
1=\left(N+1\right)\times 1$.

By selecting $n_{0}$ as the physical ground node and datum of the circuit, we
set its respective flux $\Phi_{n_{0}}$ to zero. We need to thus eliminate
$\Phi_{n_{0}}$ as a degree of freedom and reduce the flux basis to an
$N$-length column vector of node-to-datum fluxes $\mathbf{\Phi_{n}}$; note,
$\dim\mathbf{\Phi_{n}}=N\times 1$. We perform the basis reduction
$\mathbf{\Phi}_{\infty}=\mathbf{S}_{N}\mathbf{\Phi_{n}}$ with the nearly-
trivial linear transformation
$\mathbf{S}_{N}=\left[\begin{array}[]{cc}\mathbf{0}_{N}&\mathbf{I}_{N}\end{array}\right]^{\intercal}$,
where $\mathbf{0}_{N}$ is the column vector of all zeros with length $N$ and
$\mathbf{I}_{N}$ is the identity matrix of dimension $N$. In the reduced basis
$\mathbf{\Phi_{n}}$, the reduced capacitance matrix of the cell is

$\mathbf{C}_{\mathrm{cell}}=\mathbf{S}_{N}^{\intercal}\mathbf{C_{M}}\mathbf{S}_{N}\;.$

Due to the simple form of $\mathbf{S}_{N}$, this is equivalent to simply
dropping the first row and column of $\mathbf{C_{M}}$; i.e.,

$\mathbf{C}_{\mathrm{cell}}=\operatorname{block}\left(\mathbf{C_{M}},\mathcal{N}_{\mathrm{cell}},\mathcal{N}_{\mathrm{cell}}\right)\;,$

where the set of non-datum nodes of the cell is
$\mathcal{N}_{\mathrm{cell}}=\mathcal{N}-\left\\{n_{0}\right\\}$.

## References

  * Devoret and Schoelkopf (2013) M. H. Devoret and R. J. Schoelkopf, Science 339, 1169 (2013).
  * Gambetta  _et al._ (2017) J. M. Gambetta, J. M. Chow, and M. Steffen, npj Quantum Information 3, 2 (2017).
  * Krantz  _et al._ (2019) P. Krantz, M. Kjaergaard, F. Yan, T. P. Orlando, S. Gustavsson, and W. D. Oliver, Applied Physics Reviews 6, 021318 (2019).
  * Kjaergaard  _et al._ (2020) M. Kjaergaard, M. E. Schwartz, J. Braumüller, P. Krantz, J. I. Wang, S. Gustavsson, and W. D. Oliver, “Superconducting Qubits: Current State of Play,”  (2020), arXiv:1905.13641 .
  * Blais  _et al._ (2020) A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff,  (2020), arXiv:2005.12667 .
  * Nigg  _et al._ (2012) S. E. Nigg, H. Paik, B. Vlastakis, G. Kirchmair, S. Shankar, L. Frunzio, M. H. Devoret, R. J. Schoelkopf, and S. M. Girvin, Physical Review Letters 108, 240502 (2012).
  * Bourassa  _et al._ (2012) J. Bourassa, F. Beaudoin, J. M. Gambetta, and A. Blais, Physical Review A 86, 013814 (2012), arXiv:arXiv:1204.2237v2 .
  * Solgun  _et al._ (2014) F. Solgun, D. W. Abraham, and D. P. DiVincenzo, Physical Review B 90, 134504 (2014).
  * Solgun and DiVincenzo (2015) F. Solgun and D. P. DiVincenzo, Annals of Physics 361, 605 (2015), arXiv:1505.04116 .
  * Smith  _et al._ (2016) W. C. Smith, A. Kou, U. Vool, I. M. Pop, L. Frunzio, R. J. Schoelkopf, and M. H. Devoret, Physical Review B 94, 144507 (2016), arXiv:1602.01793 .
  * Malekakhlagh and Türeci (2016) M. Malekakhlagh and H. E. Türeci, Physical Review A 93, 012120 (2016), arXiv:1506.02773 .
  * Gely  _et al._ (2017) M. F. Gely, A. Parra-Rodriguez, D. Bothner, Y. M. Blanter, S. J. Bosman, E. Solano, and G. A. Steele, Physical Review B 95, 245115 (2017).
  * Malekakhlagh  _et al._ (2017) M. Malekakhlagh, A. Petrescu, and H. E. Türeci, Physical Review Letters 119, 073601 (2017).
  * Pechal and Safavi-Naeini (2017) M. Pechal and A. H. Safavi-Naeini, Physical Review A 96, 042305 (2017), arXiv:1706.05368 .
  * Minev  _et al._ (2018) Z. K. Minev, Z. Leghtas, S. O. Mudhada, P. Reinhold, A. Diringer, and M. H. Devoret, “pyEPR: The energy-participation-ratio (EPR) open-source framework for quantum device design,” (2018).
  * Parra-Rodriguez  _et al._ (2019) A. Parra-Rodriguez, I. L. Egusquiza, D. P. DiVincenzo, and E. Solano, Physical Review B 99, 014514 (2019), arXiv:1810.08471 .
  * Parra-Rodriguez  _et al._ (2018) A. Parra-Rodriguez, E. Rico, E. Solano, and I. L. Egusquiza, Quantum Science and Technology 3, 024012 (2018), arXiv:1711.08817 .
  * Ansari (2019) M. H. Ansari, Physical Review B 100, 024509 (2019), arXiv:1807.00792 .
  * Krupko  _et al._ (2018) Y. Krupko, V. D. Nguyen, T. Weißl, É. Dumur, J. Puertas, R. Dassonneville, C. Naud, F. W. J. Hekking, D. M. Basko, O. Buisson, N. Roch, and W. Hasch-Guichard, Physical Review B 98, 094516 (2018).
  * Malekakhlagh  _et al._ (2020a) M. Malekakhlagh, A. Petrescu, and H. E. Türeci, Physical Review B 101, 134509 (2020a), arXiv:1809.04667 .
  * Solgun  _et al._ (2019) F. Solgun, D. P. DiVincenzo, and J. M. Gambetta, IEEE Transactions on Microwave Theory and Techniques 67, 928 (2019), arXiv:1712.08154 .
  * Petrescu  _et al._ (2019) A. Petrescu, M. Malekakhlagh, and H. E. Türeci,  (2019), arXiv:1908.01240 .
  * You  _et al._ (2019) X. You, J. A. Sauls, and J. Koch, Physical Review B 99, 174512 (2019), arXiv:1902.04734 .
  * Di Paolo  _et al._ (2019) A. Di Paolo, T. E. Baker, A. Foley, D. Sénéchal, and A. Blais,  (2019), arXiv:1912.01018 .
  * Gely and Steele (2020) M. F. Gely and G. A. Steele, New Journal of Physics 22, 013025 (2020).
  * Ding  _et al._ (2020) D. Ding, H. S. Ku, Y. Shi, and H. H. Zhao, “Free mode removal and mode decoupling for simulating general superconducting quantum circuits,” (2020), arXiv:2011.10564 .
  * Kerman (2020) A. J. Kerman,  (2020), arXiv:2010.14929 .
  * Minev  _et al._ (2020) Z. K. Minev, Z. Leghtas, S. O. Mundhada, L. Christakis, I. M. Pop, and M. H. Devoret, arXiv (2020), arXiv:2010.00620 .
  * Kyaw  _et al._ (2020) T. H. Kyaw, T. Menke, S. Sim, N. P. D. Sawaya, W. D. Oliver, G. G. Guerreschi, and A. Aspuru-Guzik,  (2020), arXiv:2006.03070 .
  * Malekakhlagh  _et al._ (2020b) M. Malekakhlagh, E. Magesan, and D. C. McKay,  (2020b), arXiv:2005.00133 .
  * Menke  _et al._ (2021) T. Menke, F. Häse, S. Gustavsson, A. J. Kerman, W. D. Oliver, and A. Aspuru-Guzik, npj Quantum Information 7, 49 (2021), arXiv:1912.03322 .
  * Wallraff  _et al._ (2004) A. Wallraff, D. I. Schuster, A. Blais, L. Frunzio, R.-S. Huang, J. Majer, S. Kumar, S. M. Girvin, and R. J. Schoelkopf, Nature 431, 162 (2004), arXiv:0407325 [cond-mat] .
  * Paik  _et al._ (2011) H. Paik, D. I. Schuster, L. S. Bishop, G. Kirchmair, G. Catelani, A. P. Sears, B. R. Johnson, M. J. Reagor, L. Frunzio, L. I. Glazman, S. M. Girvin, M. H. Devoret, and R. J. Schoelkopf, Physical Review Letters 107, 240501 (2011), arXiv:1105.4652 .
  * Barends  _et al._ (2013) R. Barends, J. Kelly, A. Megrant, D. Sank, E. Jeffrey, Y. Chen, Y. Yin, B. Chiaro, J. Mutus, C. Neill, P. O’Malley, P. Roushan, J. Wenner, T. C. White, A. N. Cleland, and J. M. Martinis, Physical Review Letters 111, 080502 (2013), arXiv:1304.2322 .
  * Minev  _et al._ (2013) Z. Minev, I. M. Pop, and M. H. Devoret, Applied Physics Letters 103, 142604 (2013), arXiv:arXiv:1308.1743v1 .
  * Yan  _et al._ (2016) F. Yan, S. Gustavsson, A. Kamal, J. Birenbaum, A. P. Sears, D. Hover, T. J. Gudmundsen, D. Rosenberg, G. Samach, S. Weber, J. L. Yoder, T. P. Orlando, J. Clarke, A. J. Kerman, and W. D. Oliver, Nature Communications 7, 12964 (2016).
  * Brecht  _et al._ (2016) T. Brecht, W. Pfaff, C. Wang, Y. Chu, L. Frunzio, M. H. Devoret, and R. J. Schoelkopf, npj Quantum Information 2, 16002 (2016), arXiv:1509.01127 .
  * Rosenberg  _et al._ (2017) D. Rosenberg, D. Kim, R. Das, D. Yost, S. Gustavsson, D. Hover, P. Krantz, A. Melville, L. Racz, G. O. Samach, S. J. Weber, F. Yan, J. L. Yoder, A. J. Kerman, and W. D. Oliver, npj Quantum Information 3, 42 (2017), arXiv:1706.04116 .
  * Versluis  _et al._ (2017) R. Versluis, S. Poletto, N. Khammassi, B. Tarasinski, N. Haider, D. J. Michalak, A. Bruno, K. Bertels, and L. DiCarlo, Physical Review Applied 8, 034021 (2017).
  * Naik  _et al._ (2017) R. K. Naik, N. Leung, S. Chakram, P. Groszkowski, Y. Lu, N. Earnest, D. C. McKay, J. Koch, and D. I. Schuster, Nature Communications 8, 1904 (2017).
  * Charles James Neill (2017) Charles James Neill,  _Thesis_ , Tech. Rep. (2017).
  * Yan  _et al._ (2018) F. Yan, P. Krantz, Y. Sung, M. Kjaergaard, D. L. Campbell, T. P. Orlando, S. Gustavsson, and W. D. Oliver, Physical Review Applied 10, 054062 (2018), arXiv:1803.09813 .
  * Minev  _et al._ (2019) Z. K. Minev, S. O. Mundhada, S. Shankar, P. Reinhold, R. Gutiérrez-Jáuregui, R. J. Schoelkopf, M. Mirrahimi, H. J. Carmichael, and M. H. Devoret, Nature 570, 200 (2019), arXiv:1803.00545 .
  * Antonio D.  _et al._ (2021) Antonio D., M. Takita, K. Inoue, S. Lekuch, Z. K. Minev, J. M. Chow, and J. M. Gambetta,  (2021), arXiv:2102.01682 .
  * Yurke and Denker (1984) B. Yurke and J. S. Denker, Physical Review A 29, 1419 (1984).
  * Devoret (1995) M. H. Devoret, in  _Quantum Fluctuations, Les Houches, Sess. LXIII_ (1995).
  * Burkard  _et al._ (2004) G. Burkard, R. H. Koch, and D. P. DiVincenzo, Physical Review B 69, 064503 (2004).
  * Koch  _et al._ (2007) J. Koch, T. M. Yu, J. Gambetta, A. A. Houck, D. I. Schuster, J. Majer, A. Blais, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Physical Review A 76, 42319 (2007).
  * Note (1) For an early version of the open-source code (Minev  _et al._ , 2021) developed by the authors of this manuscript, see qiskit.org/metal. The Qiskit Metal project builds on github.com/zlatko-minev/pyEPR.
  * Zimmerman and Silver (1966) J. E. Zimmerman and A. H. Silver, Physical Review 141, 367 (1966).
  * Clarke and Braginski (2004) J. Clarke and A. I. Braginski, eds.,  _The SQUID Handbook_ (Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, FRG, 2004).
  * Frattini  _et al._ (2018) N. E. Frattini, V. V. Sivak, A. Lingenfelter, S. Shankar, and M. H. Devoret, Physical Review Applied (2018), 10.1103/PhysRevApplied.10.054020, arXiv:1806.06093 .
  * Minev (2019) Z. K. Minev, Ph.D. thesis, Yale University (2019), arXiv:1902.10355 .
  * Gloos  _et al._ (2000) K. Gloos, R. S. Poikolainen, and J. P. Pekola, Applied Physics Letters 77, 2915 (2000).
  * Wolff (2006) I. Wolff,  _Coplanar Microwave Integrated Circuits_ (John Wiley & Sons, Inc., Hoboken, NJ, USA, 2006).
  * Dirac (1982) P. A. M. Dirac,  _Principles of Quantum Mechanics_ (Oxford University Press, 1982).
  * Note (2) This is the form employed in a recent package entitled scqubits; see also Ref. Kerman2020.
  * Minev  _et al._ (2021) Z. K. Minev, T. G. McConkey, J. Drysdal, P. Shah, D. Wang, M. Facchini, G. Harper, J. Blair, H. Zhang, N. Lanzillo, S. Mukesh, W. Shanks, C. Warren, and J. M. Gambetta, “Qiskit Metal: An Open-Source Framework for Quantum Device Design & Analysis,” (2021).
  * Zangwill (2012) A. Zangwill,  _Modern Electrodynamics_ (Cambridge University Press, 2012).
  * Blais  _et al._ (2004) A. Blais, R.-S. Huang, A. Wallraff, S. M. Girvin, and R. J. Schoelkopf, Physical Review A 69, 062320 (2004), arXiv:0402216 [cond-mat] .
  * Gambetta  _et al._ (2006) J. Gambetta, A. Blais, D. I. Schuster, A. Wallraff, L. Frunzio, J. Majer, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Physical Review A - Atomic, Molecular, and Optical Physics (2006), 10.1103/PhysRevA.74.042318, arXiv:0602322 [cond-mat] .
  * Gambetta  _et al._ (2008) J. Gambetta, A. Blais, M. Boissonneault, A. A. Houck, D. I. Schuster, and S. M. Girvin, Physical Review A 77, 012112 (2008).
  * Krupka  _et al._ (2006) J. Krupka, J. Breeze, A. Centeno, N. Alford, T. Claussen, and L. Jensen, IEEE Transactions on Microwave Theory and Techniques 54, 3995 (2006).
  * Simons (2001) R. N. Simons,  _Coplanar Waveguide Circuits, Components, and Systems_ (2001).
  * Minev  _et al._ (2016) Z. Minev, K. Serniak, I. M. Pop, Z. Leghtas, K. Sliwa, M. Hatridge, L. Frunzio, R. J. Schoelkopf, and M. H. Devoret, Physical Review Applied 5, 044021 (2016), arXiv:1509.01619 .


# Hypersonic flow over a wedge in the detached shock range

H. G. Hornung
Graduate Aerospace Laboratories, California Institute of Technology
<EMAIL_ADDRESS>
## 1\. Introduction

One of the most important parameters in hypersonic flow is the inverse normal-
shock density ratio, which for a perfect gas is

$\varepsilon\,=\,{\rho_{\infty}\over{\rho_{s}}}\,=\,{\gamma-1+2/M_{\infty}^{2}\over{\gamma+1}},$
(1)

where $\rho$ is density, $\gamma$ is the ratio of specific heats, and $M$ is
the Mach number. The subscripts $\infty$ and $s$ refer to the free stream and
to the immediate post-normal-shock condition.

In the flow over a wedge with given free-stream conditions, the shock wave is
straight and attached to the wedge tip when the wedge half-angle $\theta$ is
sufficiently small. As $\theta$ is increased a point is reached at which the
flow downstream of the wedge is sonic, so that information about the length of
the wedge from tip to shoulder can be communicated to the tip. The shock
begins to curve, and, at a very slightly larger value of $\theta$, it detaches
from the tip. Results of Euler computations are shown for these three
conditions in Figure 1.

Figure 1: Flow over a wedge at $M_{\infty}=10$ and $\gamma=1.3$ with
$\theta=48$, 49 and 50∘, showing (l. to r.) entirely supersonic flow with a
straight attached shock, partly subsonic flow with a curved, attached shock
and flow with a detached shock. The white line is the sonic line.

An example of the power of the parameter $\varepsilon$ is the approximation
given by [3] for the values of $\theta$ and the shock angle $\beta$ at
detachment.

$\beta_{d}\,=\,{\rm
arctan}\sqrt{1\over{\varepsilon}},\qquad\theta_{d}\,=\,2\beta_{d}-{\pi\over{2}}.$
(2)

Here the subscript $d$ refers to the detachment condition. The exact values of
these detachment angles may be determined from

$\beta_{de}\,=\,{\rm
arcsin}\sqrt{{{(\gamma+1)M_{\infty}^{2}/4-1+\sqrt{\gamma+1}\sqrt{(\gamma+1)M_{\infty}^{4}/16+1+(\gamma-1)M_{\infty}^{2}/2}}\over{\gamma
M_{\infty}^{2}}}}$ (3) $\theta_{de}\,=\,{\rm
arctan}\left({(M_{\infty}^{2}\sin^{2}\beta_{de}-1)/\tan\beta_{de}\over{1+[(\gamma+1)/2-\sin^{2}\beta_{de}]M_{\infty}^{2}}}\right)$
(4)

Figure 2 shows the quality of the approximation of [3] by plotting exact
values for Mach numbers between 4 and 10 and $\gamma$ between 1.05 and 1.4
together with the approximation vs. $\varepsilon$. As $\varepsilon$ increases
the approximate value of $\theta_{d}$ falls above the exact curves, especially
at the lower Mach numbers.

Figure 2: Exact detachment angles for $4<M_{\infty}<10$ and $1.05<\gamma<1.4$,
plotted with black dashed lines, compared with the approximation of [3]
(equations 3 and 4) in red.

The range of $\theta$ for which the shock is detached is of particular
interest here. As in [5], we introduce the variable

$\eta\,=\,{\theta-\theta_{d}\over{\pi/2-\theta_{d}}},$ (5)

such that $\eta=0$ at detachment and $\eta=1$ at $\theta=\pi/2$. Since
analytical formulas for the exact detachment angles exists in the case of
wedge flow, we use the variable

$\eta_{e}\,=\,{\theta-\theta_{de}\over{\pi/2-\theta_{de}}}$ (6)

in place of $\eta$.

For the shock stand-off distance, $\Delta$, which is of special interest, we
again make the hypothesis that it follows the functional form

${\Delta\over{H}}\,=\,g(\varepsilon)f(\eta_{e})$ (7)

as in [5]. Here $H$ is the height of the wedge measured from the symmetry
plane to the shoulder. In order to test the hypothesis, we make a large number
of computations covering the parameter space $(M_{\infty},\gamma,\theta)$. If
the hypothesis is true in the case of flow over a wedge as it was for cone
flow, the results can be used to determine the functional forms of $g$ and
$f$.

Figure 3: Pseudo-schlieren images of flow over a wedge at $M_{\infty}=5$. The
grey-shading in these images is proportional to a monotonic function of the
magnitude of the density gradient. The white line is the sonic line. In rows
from left to right, $\theta=60$, 70, 80 and 90 deg. In columns from top to
bottom, $\gamma=1.05$, 1.1, 1.2 and 1.3. Similar sets of computations were
made for $M_{\infty}=4$, 7 and 10.

## 2\. Discussion of results of computations

### 2.1. Shock stand-off distance

The parameter space was explored by computing the flow over wedges using the
Euler equations. Details about the computational technique are given in the
Appendix. An example of the results is presented in Figure 3 for the case of
$M_{\infty}=5$. Similar results were also obtained for $M_{\infty}=4$, 7 and
10. Plotting the dimensionless shock stand-off distance against $\varepsilon$
in the case $\theta=90^{\circ}$, (square slab), i. e., for $\eta_{e}=1$, where
$f(\eta_{e})=f(1)$ is a constant, provides a partial test of the hypothesis
for the function $g(\varepsilon)$.

Figure 4: LEFT: Dimensionless shock stand-off distance for flow over a 90∘
wedge with $M_{\infty}=4$, 5, 7 and 10 and $\gamma=1.05$, 1.1, 1.2 and 1.3. In
the case of $M_{\infty}=10$, two cases of $\gamma=1.01$ and 1.02, and with
$M_{\infty}=7$, one case with $\gamma=1.02$ were added. RIGHT: Corresponding
plot for flow over a circular cylinder of radius $R$. In this case the results
are fitted well by $\Delta/R=2.14\,\varepsilon(1+\varepsilon/2)$. In [4] the
linear form $\Delta/R=2.32\,\varepsilon$ was found, which is only very
slightly different.

Figure 5: Four plots of $f(\eta_{e})$ vs. $\eta_{e}$. Top, left to right:
$M_{\infty}=10$ and 7. Bottom, left to right $M_{\infty}=5$ and 4. In all four
plots the curve is the same and given by equation 9.

This has been done in Figure 4. A fit to the points in Figure 4 yields the
interesting result that all points fall on a unique curve given by

$g(\varepsilon)\,=\,\sqrt{\varepsilon}\left(1+{3\over{2}}\varepsilon\right),$
(8)

thus confirming the first part of the hypothesis, i. e., that a unique
function $g(\varepsilon)$ exists. Again, as in the case of flow over cones,
the leading term is proportional to $\sqrt{\varepsilon}$ and the fact that all
the results fall on the same curve confirms the first part of the hypothesis.
It is interesting that in the corresponding function of $\varepsilon$ for the
90∘ cone, the factor 3/2 that appears in equation 8 is 1/2.

To test the second part, four plots of
$f(\eta_{e})=\Delta/[H\,g(\varepsilon)]$ vs. $\eta_{e}$ are shown in Figure 5
for the four Mach numbers and the four $\gamma$’s. All the computational
results fall on the same curve given by the unique function

$f(\eta_{e})\,=\,2.2\,\eta_{e}-0.3\,\eta_{e}^{2},$ (9)

thus confirming the second part of the hypothesis. Note however, that the fit
is not so good at $M_{\infty}=4$ as for the others. Had we used $\eta$ instead
of $\eta_{e}$, this deterioration of the fit at the lower Mach numbers would
have been slightly larger.

### 2.2. Drag coefficient

In [5] it was shown that for flow over cones, the drag coefficient could also
be expressed in the form of equation 5. Write the drag coefficient for a wedge
as

$C_{D}\,=\,{2D\over{\gamma p_{\infty}M_{\infty}^{2}HL}},$ (10)

where the drag force

$D\,=\,2L\int_{0}^{H}(p-p_{\infty})dy,$

$L$ is the transverse length of the wedge, and $y$ is the distance measured
from the symmetry plane of the wedge. Then, if the form of equation 5 holds
for $C_{D}$, expect that

$C_{D}\,=\,g_{1}(\varepsilon)f_{1}(\eta_{e}).$ (11)

The same set of computational results can now be used to check whether this is
correct. Again we use the case $\eta_{e}=1$ where $f_{1}(1)$ is a constant to
check if $g_{1}(\varepsilon)$ is unique. To this end, Figure 6 shows a plot of
$C_{D}$ vs. $\varepsilon$ for the 90∘ wedge and for the circular cylinder. In
both cases all the results collapse onto a single line. for the 90∘ wedge,

$C_{d}\,=\,g_{1}(\varepsilon)\,=\,2-1.4\,\varepsilon,$ (12)

and for the circular cylinder

$C_{D}\,=\,1.3-5(\varepsilon-0.085)^{2}.$ (13)

Figure 6: LEFT: Drag coefficient for flow over a 90∘ wedge with
$M_{\infty}=4$, 5, 7 and 10 and $\gamma=1.05$, 1.1, 1.2 and 1.3. For
$M_{\infty}=10$, two cases of $\gamma=1.01$ and 1.02, and with $M_{\infty}=7$,
one case with $\gamma=1.02$ were added. RIGHT: Drag coefficient for flow over
a circular cylinder with $M_{\infty}=4$, 5, 7 and 10 and $\gamma=1.05$, 1.1,
1.2 and 1.3.

Figure 7: Four plots of $f_{1}(\eta_{e})$ vs. $\eta_{e}$. Top, left to right:
$M_{\infty}=10$ and 7. Bottom, left to right $M_{\infty}=5$ and 4. In all four
plots the line is the same and given by equation 14.

In order to determine the form of the function $f_{1}(\eta_{e})$ Figure 7
shows four plots of $\Delta/[Hg_{1}(\varepsilon)]$ vs.$\eta_{e}$. While the
results agree approximately with

$0.85+0.15\,\eta_{e},$ (14)

they scatter fairly broadly around the line, so that the validity of the
functional form is not as convincing as in the case of the shock stand-off
distance.

## 3\. High-enthalpy effects

Although the results were obtained from perfect-gas computations, they may be
expected to apply also to flows at high enthalpy, where vibrational excitation
and dissociation may occur and non-equilibrium effects become important. As
has been shown by [11] and by [10], forming the density ratio with the average
density along the stagnation streamline instead of with the post-normal-shock
density causes the results to carry over to the high-enthalpy regime.

## 4\. Conclusions

It was shown that the reduction from three independent parameters to two, in
the parameter space defining the dimensionless shock stand-off distance and
drag coefficient, that was previously found for axisymmetric hypersonic flow
over cones, applies also to flow over wedges in the detached-shock range.
Useful analytical forms were found for these relations by performing a large
number of Euler computations and plotting the results in the appropriate form.
The results are also compared with new results for flow over a circular
cylinder.

Acknowledgement
This work was funded by AFOSR FA9550-19-1-0219, PI: J. M. Austin, Contract
Officer: I. A. Leyva.

Appendix
The software system Amrita, constructed by James Quirk, see [8], was used. A
detailed description of the features and phenomena encountered with some of
the algorithms used for Riemann solvers, including the one used here, has also
been given by [7]. An example of a test of the software against experiment may
be found in [9].

Amrita is a system that automates and packages computational tasks in such a
way that the packages can be combined (dynamically linked) according to
instructions written in a high-level scripting language. The present
application uses features of Amrita that include the automatic construction of
the Euler solver, documentation of the code, adaptive mesh refinement
according to simply chosen criteria, and scripting-language-driven
computation, archiving and post-processing of the results. The automation of
the assembly and sequencing of the tasks makes for dramatically reduced
possibility of hidden errors. It also makes computational investigations
transparent and testable by others. The ability to change one package at a
time, without changing the rest of the scheme, facilitates detection of
sources of error. In most of the work, the Euler solver generated was an
operator-split scheme with HLLE flux (after Harten et al.[2] and Einfeldt [1])
and kappa-MUSCL reconstruction. In some cases with $\gamma$ close to 1, the
carbuncle problem arose, and the more robust equilibrium flux method of [6]
was used. The ($x,y$) plane was discretized by a Cartesian grid of
300$\times$300 coarse–grid cells that are adaptively refined by a factor of 3
to make an effective grid of 900$\times$900 cells. The criterion for
adaptation was a chosen threshold of the magnitude of the density gradient.
Solid boundaries are represented by a level set defined as the smallest
distance of a field point from the solid boundary. The grey-shading of the
visualizations is a monotonic function of the magnitude of the density
gradient.

## References

  * [1] B. Einfeldt.  On Godunov-type methods for gas dynamics.  SIAM J. Numer. Anal., 25:294–328, 1988.
  * [2] A. Harten, P. D. Lax, and B. van Leer.  On upstream differencing and Godunov-type schemes for hyperbolic conservation laws.  SIAM Rev., 25:35–61, 1983.
  * [3] W. D. Hayes and R. F. Probstein.  Hypersonic flow theory.  Academic Press, 1959.
  * [4] H. G. Hornung.  Non-equilibrium flow of nitrogen over spheres and circular cylinders.  J. Fluid Mech., 53:149–176, 1972.
  * [5] H. G. Hornung, J. Martinez Schramm, and K. Hannemann.  Hypersonic flow over sphericcally blunted cone capsules for atmospheric entry. Part 1, the sharp cone and the sphere.  J. Fluid Mech., 871:1097–1116, 2019.
  * [6] D. I. Pullin.  Direct simulation methods for compressible inviscid ideal-gas flows.  J. Comput. Phys, 34:231–240, 1980.
  * [7] J. J. Quirk.  A contribution to the great Riemann solver debate.  Int. J. for Num. Methods in Fluids, 18:555–574, 1994.
  * [8] J. J. Quirk.  Amrita — a computational facility (for CFD modelling).  In VKI CFD Lecture Series, volume 29. von Karman Institute, 1998\.
  * [9] J. J. Quirk and S. Karni.  On the dynamics of a shock bubble interaction.  J. Fluid Mech., 318:129–163, 1996.
  * [10] V. P. Stulov.  Similarity law for supersonic flow past blunt bodies.  Izv. AN SSSR, Mechanika Zhidkosti i Gaza, 4:142–146, 1969.
  * [11] C.-Y. Wen and H. G. Hornung.  Non-equilibrium dissociating flow over spheres.  J. Fluid Mech., 299:389–405, 1995.