# Spectral Theory for Non-linear Superconducting Microwave Systems: Extracting
Relaxation Rates and Mode Hybridization

Dung N. Pham  Department of Electrical and Computer Engineering, Princeton
University, Princeton, NJ 08544, USA  Richard D. Li  Department of Physics,
Yale University, New Haven, CT 06511, USA  Hakan E. Türeci  Department of
Electrical and Computer Engineering, Princeton University, Princeton, NJ
08544, USA

(August 28, 2024)

###### Abstract

The accurate modeling of mode hybridization and calculation of radiative
relaxation rates have been crucial to the design and optimization of
superconducting quantum devices. In this work, we introduce a spectral theory
for the electrohydrodynamics of superconductors that enables the extraction of
the relaxation rates of excitations in a general three-dimensional
distribution of superconducting bodies. Our approach addresses the long-
standing problem of formulating a modal description of open systems that is
both efficient and allows for second quantization of the radiative hybridized
fields. This is achieved through the implementation of finite but transparent
boundaries through which radiation can propagate into and out of the
computational domain. The resulting spectral problem is defined within a
coarse-grained formulation of the electrohydrodynamical equations that is
suitable for the analysis of the non-equilibrium dynamics of multiscale
superconducting quantum systems.

## Introduction

Engineered superconducting systems employed in analog quantum simulation [1,
2], quantum sensing [3, 4] and quantum computing [5, 6] are described by the
electrohydrodynamic theory of superconducting systems (EHDS), first proposed
without justification by R.P. Feynman in a special lecture he delivered in the
1960s [7]. Few attempts have been made to either justify [8, 9] or solve [10]
this model’s equations of motion exactly, but simplified versions [11, 12, 13,
12] have been considered via various modifications to minimal coupling. This
low-energy theory can be credited as the underpinning of Josephson phenomena,
flux quantization, and, generally, the physics of vortices. The second-
quantized version of it constitutes the basis of the formulation known as
circuit quantum electrodynamics (cQED) [14], which has been the workhorse
behind the current understanding of the physics governing superconducting
quantum computers. As systems grow increasingly sophisticated and their
electromagnetic environments become more complex, efficient computational
strategies in both classical and quantum regimes are much needed to produce
accurate reduced models containing the degrees of freedom of interest. This
has given rise to an active research area at the intersection of computational
electromagnetism and quantum electrodynamics of superconducting devices [15,
16, 17, 18, 19, 20, 21], to analyze and model the physics of these circuits.

An indispensable component in the electrodynamic modeling of superconducting
circuits is the extraction of relaxation rates. Purcell modification of
radiative lifetimes of qubits is today an essential mechanism for the
protection of qubits and plays an important role in the electromagnetic design
of the entire processor [22, 23, 24, 25]. Additionally, in applications that
require strong RF or microwave excitations of individual oscillators, new
dissipative processes and inter-level transitions can be activated [26, 27,
28, 29, 30, 31, 32, 33, 34]. These processes, in turn, depend strongly on the
spectral characteristics of the electromagnetic environment in which the non-
linear elements are embedded [30]. Resource-intensive numerical simulations
are needed to extract such information for complex quantum-electrodynamic
systems. Such simulations are not only computationally demanding but also face
conceptual difficulties arising from the infinite degrees of freedom inherent
in quantum electrodynamics (QED) [35]. These conceptual issues [22, 36, 18,
37] can be effectively addressed by the introduction of a rigorous spectral
theory [38]. This approach allows the accurate quantification of limitations
brought about by noise [39, 40, 41], dissipation [42, 43] and undesired
interactions [44, 45] within the system and its electromagnetic environment.

In the present work, we introduce a spectral theory of EHDS from which the
relaxation rates of general three-dimensional superconducting circuits can be
obtained. This spectral theory is adapted to a coarse-grained description of
the EHDS equations, known as DEC-QED, previously derived and analyzed in Ref.
[10]. DEC-QED provides the following advantages: (1) through its structure-
preserving geometric discretization procedure, known as discrete exterior
calculus (DEC), this formulation enables stable long-time simulations. (2) By
virtue of the fundamental fields being hybridized gauge-invariant fields
rather than the standard electromagnetic potentials, different materials can
be handled in a uniform fashion. (3) At a superconducting-normal-
superconducting (SNS) junction, the hybridized field is identical to the
gauge-invariant Josephson phase across the junction. (4) For transmission
lines and lumped-element circuits and in the limit of a perfect superconductor
($\lambda_{L}$ = 0), DEC-QED equations reduce to the standard cQED equations.

The difficulty with developing a spectral theory stems from the desire to
solve the EHDS equations, a set of non-linear PDEs describing the evolution of
the order parameter of a charged fluid coupled to Maxwell’s equations, in a
finite spatial domain. Ideally, one would like to keep the computational
domain as small as possible, just as large as the superconducting system
itself. However, this is not possible, because we also need to keep the volume
as large as possible to allow for radiative relaxation of the excitations in
the superconductor. The optimal solution to avoid this trade-off is to use a
modal description for a finite but open domain subject to transparent boundary
conditions for electromagnetic (EM) radiation. This allows us to keep the
computational domain small enough to be computationally feasible, while still
allowing for the radiative relaxation of the excitations. While boundary
conditions, such as the perfectly matched layer (PML) [46, 47] or absorbing
boundary conditions [48, 49] are known and very useful in classical EM
problems, it is also desirable to implement a transparent boundary so that
solving the corresponding spectral problem would provide a modal decomposition
which can serve as the basis for second quantizing the EHDS equations. This,
however, has been a difficult problem that remains unaddressed since the full
formulation of quantum electrodynamics [50] because such boundary conditions
result in non-Hermitian modes [51, 52, 53]. Recently, it was shown that by
using the Heisenberg equations of motion for quantum field operators, this
issue can be circumvented in the context of one-dimensional transmission line
systems through the use of singular function expansion and a suitable spectral
problem [38]. Here we generalize the statement of this spectral problem to the
solution of the non-linear EHDS equations (the original spectral problem has
been defined [38] in the context of cQED which assumes fields do not penetrate
the superconductors), and to general three dimensional domains. We do not
address the issue of quantization but focus on the spectral problem that
arises through the linearization of the EHDS equations, the resulting non-
Hermitian modes, and the calculation of dissipation rates for all modes. The
dissipation rates for “qubit-like” modes are precisely their Purcell radiative
lifetimes [54], while for the more spatially extended modes, they represent
their losses due to their hybridization with radiative channels in the system
[55, 56, 57].

The main results presented in this paper are organized as follows: we first
discuss the electrohydrodynamic model used to describe the dynamics of the
superconducting condensate coupled to the EM environment and derive the
spectral problem to be solved. We then briefly discuss the formulation of DEC-
QED used in the rest of the article for coarse-grained calculations. Next, we
demonstrate the convergence in the calculations for the spectral profile of a
sample superconducting cavity. Both simplicial and cubical meshing strategies
are benchmarked. We then show how mode hybridization due to accidental
degeneracies can be encountered in the calculations of the modes for a
symmetric cavity and that our method can distinguish the degenerate modes
through the inclusion of a small perturbation to the cavity shape. Next, we
demonstrate the coarse-grained calculations of hybridized modes for multiscale
systems containing multiple components. Finally, we present two main
approaches for implementing open boundary conditions based on (1) Green’s
boundary integrals and (2) vector spherical harmonic expansions, along with
numerical examples, to demonstrate the accuracy and effectiveness of these
methods. We have also made the implementation of the formulation discussed
here available as part of an open-source DEC-QED repository [58].

## Results

### Electrohydrodynamic formulation of superconducting materials

Consider a multiply connected region composed of different dielectric and
superconducting materials. Well below the critical temperature, where most
superconducting devices operate, the superconducting material can be modeled
as a charged condensate [8, 9] trapped by a background of positive charge
($\rho_{\text{src}}$). The dynamics of the condensate are wholly determined
through minimal coupling to the dynamical electromagnetic potentials
($\bm{A}$, $V$) and the Coulomb attraction to the static positive background
($U$) defining the superconductor. We refer to the resulting equations as the
ElectroHydroDynamics of Superconductors (EHDS) [7]

$i\hbar\frac{\partial\Psi({\bf r},t)}{\partial
t}=\bigg{[}\frac{1}{2m}\Big{(}\\!\\!-i\hbar{\bf\nabla}-q{\bf
A}\Big{)}^{2}+qV({\bf r},t)+U({\bf r})\bigg{]}\Psi\left({\bf r},t\right)$ (1)

and Maxwell’s equations

$\displaystyle{\bf\nabla}\times{\bf\nabla}\times{\bf
A}+\mu_{0}\epsilon\ddot{{\bf A}}$ $\displaystyle=\mu_{0}({\bf J}_{s}+{\bf
J}_{\text{src}})-\mu_{0}\epsilon\nabla\dot{V},$ (2)
$\displaystyle\nabla^{2}V+\frac{\partial}{\partial
t}({\mathbf{\nabla}}\cdot{\mathbf{A}})$
$\displaystyle=-\frac{q}{\epsilon}(\rho+\rho_{\text{src}}).$ (3)

where $\rho$ and ${\bf J}_{s}$ are the condensate density and the
supercurrent, respectively, and $\rho_{\text{src}}$ and ${\bf J}_{\text{src}}$
are the external charge and current sources. Here, $q=2e$ and $m=2m_{e}$ are
the charge and mass of a cooper pair, respectively, which are twice the charge
$e$ and mass $m_{e}$ of an electron. Generally, we include the positive
background defining the superconducting regions in $\rho_{\text{src}}$.
Dielectric regions are defined by
$\epsilon(\bm{r})=\tilde{\epsilon}(\bm{r})\epsilon_{0}$, where
$\tilde{\epsilon}$ is the step-wise constant relative permittivity function
that is unity where dielectric is not present and is greater than one
otherwise. Using the Madelung representation for the condensate wavefunction,
$\Psi({\bf r},t)=\sqrt{\rho({\bf r},t)}e^{i\theta({\bf r},t)}$, and
introducing the gauge-invariant hybridized field $\bm{\mathcal{A}}={\bf
A}-\frac{\hbar}{q}\nabla\theta$, Eqs. (1)-(3) can be written in the form [10]

$\displaystyle{\bf\nabla}\\!\times\\!{\bf\nabla}\\!\times\\!\bm{\mathcal{A}}+\mu_{0}$
$\displaystyle\epsilon_{0}\frac{\partial^{2}\bm{\mathcal{A}}}{\partial
t^{2}}+\frac{\mu_{0}q^{2}}{m}\rho\bm{\mathcal{A}}-\frac{\mu_{0}\epsilon_{0}q}{2m}\frac{\partial}{\partial
t}\nabla\big{|}\bm{\mathcal{A}}\big{|}^{2}$
$\displaystyle+\frac{\mu_{0}\epsilon_{0}\hbar^{2}}{2mq}\frac{\partial}{\partial
t}\nabla\bigg{[}\frac{\nabla^{2}(\sqrt{\rho})}{\sqrt{\rho}}\bigg{]}=\mu_{0}{\bf
J}_{src},$ (4)

and

$\frac{\partial\rho}{\partial
t}={\bf\nabla}\cdot\Bigg{[}\frac{q}{m}\rho\bm{\mathcal{A}}-\frac{{\bf
J}_{src}}{q}\Bigg{]}-\frac{\partial\rho_{src}}{\partial t}.$ (5)

These equations were derived and the resulting real-time dynamics under
specific conditions were analyzed in Ref. [10]. Here, we are interested in the
spectral problem associated with these non-linear equations. Linearization can
be done by splitting the condensate density $\rho$ into the mean value
$\rho_{0}$ that exactly balances the positively charged ionic background and
the fluctuation $\delta\rho$ arising from interactions with the external EM
field. The linear sector of Eq. (Electrohydrodynamic formulation of
superconducting materials) that corresponds to the transverse excitations then
reproduces London theory [59] and is sourced by the fluctuation $\delta\rho$.
The resulting inhomogeneous source-field equations can be solved through the
spectral problem of a vector Helmholtz equation for the gauge-invariant
hybridized field $\bm{{\mathcal{A}}}$:

$\displaystyle{\mathbf{\nabla}}\times{\mathbf{\nabla}}\times\bm{{\mathcal{A}}}+\bigg{(}\frac{1}{\lambda_{L}^{2}({\mathbf{r}})}-n^{2}({\mathbf{r}})k^{2}\bigg{)}\bm{{\mathcal{A}}}=0,$
(6)

where the distribution of superconducting and dielectric materials is given by
the local London penetration depth function
$\lambda_{L}({\mathbf{r}})=\sqrt{\frac{m}{\mu_{0}q^{2}\rho_{0}({\mathbf{r}})}}$
and the refractive index function
$n({\mathbf{r}})=\sqrt{\tilde{\epsilon}({\mathbf{r}})}$, respectively, and $k$
is the wavenumber.

### DEC-QED formulation

DEC-QED provides a spatially coarse-grained description of a physical system
governed by nonlinear PDEs. The fundamental field variables are small
integrals of the original microscopic continuous fields over finite spatial
intervals. This results in a discretized model that is computationally
efficient and still accurate within the resolution of a given measurement
apparatus. In this section, we provide a minimal introduction to the geometric
constructions in DEC that are needed for computing electromagnetic modes in
systems that may contain an arbitrary distribution of superconducting and
dielectric materials, all within a finite computational domain. For a more
comprehensive discussion of the formulation, we refer to Ref. [10].

In DEC, the discretization of PDEs requires a dual-mesh construction, within
which the $d$-dimensional computational space is discretized by a primal mesh
$M$ that conforms to the boundaries of the enclosed physical domain and the
interfaces between materials. The fundamental building blocks of the primal
mesh can be simplices (i.e. triangles in 2D and tetrahedra in 3D) or cubical
elements. The vertices of the dual mesh $M^{\dagger}$ are then circumcenters
of the primal $d-$cells, and the edges of $M^{\dagger}$ are generated by
connecting the neighboring circumcenters. Throughout this paper, we will use
$\dagger$ to denote a dual quantity. For elements strictly inside the
computational space, the mappings from vertices ($v$), edges ($e$), faces
($f$), and cells ($c$) in $M$ to cells ($v^{\dagger}$), faces ($e^{\dagger}$),
edges ($f^{\dagger}$), vertices ($c^{\dagger}$) respectively in $M^{\dagger}$
are one-to-one. This bijective correspondence is not applicable at the
computational boundary, where the dual elements are truncated, which leads to
auxiliary dual nodes lying on the boundary of $M$ (See Figs. (1a) and (1b)).
During numerical calculations, this truncation of the dual mesh is taken care
of by the appropriate application of boundary conditions.

Figure 1: (a) An example primal mesh. The highlighted edges at the lower-left
corner form the boundary of an auxiliary dual volume. (b) Close-up view of a
node $v$ that lies on the boundary of the primal mesh. The neighboring dual
nodes are labeled in dark green, while the auxiliary boundary dual nodes are
labeled in bright green. The shaded volume is the truncated volume
$v^{\dagger}$ dual to the node $v$. (c) 2D example of a vertex that lies at
the interface among multiple material regions. The material properties
assigned to this node is the weighted average of the values in the surrounding
regions. (d) 3D example of an edge lying at a material interface composed of
multiple cells ${c_{i}}$. The darkened triangle is the intersection of the
dual face $e^{\dagger}$ with the cell $c_{i}$, and the shaded volume is the
portion of the support volume of $e$ that lies inside $c_{i}$.

In this DEC framework, scalar fields live on either primal vertices $v$ or
their dual volume $v^{\dagger}$, while vectorial quantities are projected onto
the discrete primal edges $e$ or assigned to the dual faces $e^{\dagger}$. For
example, given a vector field $\bm{{\mathcal{A}}}$ we construct the coarse-
grained edge field

$\Phi(e)=\int_{e}{\mathbf{d}\ell}\cdot\bm{{\mathcal{A}}},$ (7)

where the integral is done along the primal edge $e$. By solving the equations
governing these coarse-grained fields, we can probe the properties of the
system.

To properly account for the distribution of different materials within the
computational domain, material properties such as the London penetration depth
$\lambda_{L}$ or dielectric function $n$ can be assigned to objects living on
the dual mesh. For example, for every edge $e$ a value for $\lambda_{L}$ is
assigned to its dual $e^{\dagger}$. This procedure is particularly convenient
if $e$ lies at the interface between multiple materials. The dual
$e^{\dagger}$ in such cases would intersect with all the different material
domains that share this edge, and the effective penetration depth
$\overline{\lambda}_{L}(e^{\dagger})$ there will be the weighted average of
the values of $\lambda_{L}$ in the surrounding domains (See Fig. 1d). In a
three-dimensional setting, this is formally defined as followed: let
$\\{c_{i}\\}$ be the list of all the cells that share an edge $e$, then

$\overline{\lambda}_{L}(e^{\dagger})=\sum_{c_{i}}\frac{|e^{\dagger}\cap
c_{i}|}{\Delta A(e^{\dagger})}\lambda_{L}(c_{i}),$ (8)

where $|e^{\dagger}\cap c_{i}|$ in the area of the intersection between
$e^{\dagger}$ and $c_{i}$, and $\Delta A(e^{\dagger})$ is the area of
$e^{\dagger}$.

With this setup for DEC, in the following sections, we demonstrate with
specific examples how this coarse-grained formulation is utilized to compute
electromagnetic modes.

### Modes of closed superconducting systems

#### Performance benchmarks for simplicial- and cubical-DEC

Figure 2: (a) Convergence in the eigenvalues $E=k^{2}$ of the vector Helmholtz
equation applied to a perfect superconducting cavity as a function of the
number of edges $N_{e}$ used in the discretization. The solid lines correspond
to results obtained using simplicial meshing, while the dashed lines are
results obtained using cubical meshing. The relative error is calculated with
respect to the exact analytical solution. (b) The relative error in the edge
fields of the eigenmodes with respect to analytical solutions. The errors are
averaged over all edges and plotted as a function of $N_{e}$. The cavity
dimensions are $L_{x}=1\text{\,cm},L_{y}=1.5\text{\,cm}$, and
$L_{z}=2\text{\,cm}$.

We are interested in solving Eq. (6) for the electromagnetic modes of systems
composed of possibly spatially disjoint superconducting structures. To
numerically compute the modes of such systems, we derive the DEC equations
corresponding to Eq. (6)

$\displaystyle\sum_{e_{0}\in\partial(e^{\dagger})}$
$\displaystyle\sum_{e_{1}\in\partial(e_{0}^{\dagger})}\frac{\Delta\ell(e_{0})}{\Delta
A(e_{0}^{\dagger})}\Phi(e_{1})$ (9)
$\displaystyle+\bigg{(}\\!\frac{1}{\overline{\lambda}_{L}^{2}(e^{\dagger})}-\overline{n}^{2}(e^{\dagger})k^{2}\\!\bigg{)}\frac{\Delta
A(e^{\dagger})}{\Delta\ell(e)}\Phi(e)=0,$

where $\partial(e^{\dagger})$ is the boundary of the dual face $e^{\dagger}$,
and $\\{\overline{\lambda}_{L},\overline{n}\\}$ are defined as in Eq. (8).
This form is universal and is independent of the mesh elements used.

First, we investigate the numerical convergence of DEC equations using two
types of elements: simplicial and cubical. Simplicial elements are generally
the preferred choice for a meshing scheme that conforms to arbitrary
superconducting domains. To analyze the convergence of simplicial meshing, we
consider shapes of superconducting cavities for which analytical solutions are
available. A good choice is a rectangular cavity, which conforms well to
cubical elements. A comparison of the convergence of error of simplicial and
cubical meshing would illustrate the efficacy of simplicial meshing.

Consider a three-dimensional rectangular cavity with dimensions
$L_{x}=1\text{\,cm},L_{y}=1.5\text{\,cm}$, and $L_{z}=2\text{\,cm}$, where
$L_{x},L_{y},$ and $L_{z}$ are the sizes of the cavity along $x,y,$ and $z$
respectively. To allow for direct comparisons with analytical solutions, we
first assume the cavity walls are perfect superconductors
$(\lambda_{L}\rightarrow 0)$, so that the tangential component of
$\bm{{\mathcal{A}}}$ vanishes at the boundary. This is equivalent to the
boundary condition

$\Phi(e_{t})=0,$ (10)

where $e_{t}$ is an edge lying tangentially on the cavity boundary. The
eigenmodes of the cavity are computed with two choices of meshing: tetrahedral
and cubical. The convergence with respect to analytical solutions is shown in
Fig. 2. As shown in Fig. 2a, the error in the eigenvalues $E=k^{2}$ of both
mesh choices converge at the same rate as the number of edges $N_{e}$ in the
discretized domain increases. The effectiveness of simplical-DEC is therefore
shown to be comparable to cubical-DEC, even when it is applied to a geometry
where cubical meshing holds an advantage due to its elements sharing the same
symmetry as the cavity. In Fig. 2b, we also study how simplicial-DEC produces
accurate solutions to the edge fields of the eigenmodes. The convergence to
analytical solutions of the edge fields is achieved.

These results help justify our shift to using simplicial meshing onwards, as
its flexibility allows us to apply DEC to complicated, realistic structures
where cubical symmetry is rarely present. Moreover, simplicial meshing allows
for the implementation of different spatial resolutions in different regions,
an important feature needed for the efficient application of DEC to systems
comprised of multiple spatial scales.

Figure 3: Comparison between the convergence of computed eigenvalues for a
cavity embedded in a superconducting shell and those for a perfect cavity with
the hard-wall boundary condition imposed. $N_{e}$ is the number of edges used
in discretizing the cavity. The dimensions of the vacuum regions within both
cavities are $L_{x}=1\text{\,cm},L_{y}=1.5\text{\,cm}$, and
$L_{z}=2\text{\,cm}$. A schematic of the embedded cavity is shown in the inset
at the center top, while the inset at the top-right corner shows schematically
how the effective penetration depth is assigned on the inner boundary of this
cavity via coarse-graining.

Next, we consider a rectangular cavity enclosed by a superconducting shell
that has a finite thickness. The penetration depth of this shell is set to be
short enough so that the field inside the cavity decays immediately at the
material interface, i.e. the inner walls of the cavity. Instead of directly
imposing the Dirichlet boundary condition at the cavity inner walls as in the
previous case, the values of the field are left floating there, and we only
impose a “hard-wall” boundary condition (Eq. (10)) at the outer boundary of
the shell, which is also the boundary of the computational domain. The purpose
of this numerical experiment is to investigate the validity of DEC when there
are sharp interfaces, such as the vacuum-superconductor boundary here, where
the field undergoes abrupt variations. The procedure for applying material
properties to the edges lying on the interface is given in Eq. (8). In a non-
uniform tetrahedral mesh, the partitioning of the dual faces $e^{\dagger}$
that lie on multi-material interfaces to the neighboring material domains is
also highly non-uniform. This leads to the effective penetration depth of each
edge on the same boundary being different from one another. One can imagine
the vacuum-material interface made of many small patches, each with a slightly
different material property, as shown in the top-right inset in Fig. 3. We
compute the eigenmodes of Eq. (9) for an embedded cavity that has the same
dimensional ratios as in the perfect cavity case, and the convergence is shown
in Fig. 3. The embedded cavity calculation achieves a similar order of
accuracy as the perfect cavity case, and the two converge at the same rate as
mesh density increases.

#### Detection and removal of hybridization between degenerate modes

Figure 4: Demonstration of the appearances and removals of hybridization
between degenerate modes in a system with hidden symmetries. (a) The
degenerate $3^{\text{rd}}$ and $4^{\text{th}}$ modes of a rectangular cavity
with perfect superconducting boundaries and dimension ratios of
$1\\!:\\!1.5\\!:\\!2$. (b) The $3^{\text{rd}}$ and $4^{\text{th}}$ modes of a
perturbed cavity, when the dimension ratios are now $1\\!:\\!1.5\\!:\\!2.01$.
The two modes are now well-characterized by their respective set of quantum
numbers, which are schematically shown below the field distributions.

In symmetric structures such as the rectangular cavity discussed earlier in
this article, there is a possibility of hybridization of degenerate modes in
the numerically obtained eigenspectrum. These degeneracies are sometimes
called “accidental” because they are not predicted by the symmetry group of
the Hamiltonian but originate from a hidden symmetry of the system [60]. In
such cases, numerical solvers tend to face difficulties in distinguishing
these degenerate modes. For a rectangular superconducting cavity, the
eigenvalues corresponding to each field component ${\mathcal{A}_{i}}$ of the
field $\bm{{\mathcal{A}}}$ reads

$E_{i}\sim\bigg{(}\frac{n_{x}^{2}}{L_{x}^{2}}+\frac{n_{y}^{2}}{L_{y}^{2}}+\frac{n_{z}^{2}}{L_{z}^{2}}\bigg{)},$
(11)

with the corresponding eigenfunction that is separable in cartesian
coordinates

${\mathcal{A}_{i}}=X_{i}(x)Y_{i}(y)Z_{i}(z).$ (12)

If the ratio between any two of the three lengths $L_{x},L_{y}$ or $L_{z}$ are
integers, then there can be degeneracies. To numerically lift these
degeneracies, the hidden symmetry needs to be eliminated. In this specific
case, we can do so by either extending the cavity boundary to have a finite
thickness, as was done in our earlier discussion on the embedded cavity, or by
introducing a small perturbation to the cavity shape. To demonstrate this, we
consider again the perfect cavity with dimensions
$L_{x}\\!\\!:\\!\\!L_{y}\\!\\!:\\!\\!L_{z}=1\\!:\\!1.5\\!:\\!2$. In Fig. 4a,
where the $3^{\text{rd}}$ and $4^{\text{th}}$ modes of the unperturbed cavity
are plotted, we can see that each field component of the two modes is not
characterized by any single set of quantum numbers $(n_{x},n_{y},n_{z})$, but
rather a linear combination of the two accidentally degenerate modes. In Fig.
4b, we plot the same modes, but the cavity is now slightly perturbed by an
amount $\Delta L_{z}=0.01$ along $z$, so that now $L^{\prime}_{z}=L_{z}+\Delta
L_{z}=2.01$. The degeneracy is then lifted, and each of the two modes is now
well-described by a distinct set of quantum numbers, which helps demonstrate
the sensitivity of our numerical scheme to small changes in the geometry of
the system simulated.

#### Calculations for closed multiscale systems

Figure 5: (a) Schematic of a system consisting of a three-dimensional
superconducting cavity containing two spatially separated chips. Each chip is
a dielectric substrate on top of which a Josephson qubit is mounted. The qubit
is made of two superconducting capacitor pads connected together by a bridge
containing a Josephson junction. A few exemplary modes of the system are
plotted, such as (b) the mode where only qubit 1 participates, (c) only qubit
2 participates, (d) both qubits are activated, and (e) a hybridization of the
two qubits and the cavity field. The units for all the axes in plots (b)-(e)
are in millimeters.

So far, in this article, we have validated the accuracy of DEC by applying the
method to simple systems with a large degree of symmetry and comparing the
results to analytical solutions. The strength of DEC, however, lies in its
ability to correctly capture the properties of systems containing multiple
spatial scales through coarse-graining. An example of a multiscale system is a
three-dimensional superconducting cavity containing one or a few
superconducting chips whose modes can be efficiently computed using DEC. The
dimensions of a 3D cavity are typically on the order of $\sim\\!1$ cm, while a
dielectric substrate holding the qubits is a few mm in size, and a qubit
itself can have its smallest components in the range from micrometers down to
tens of nanometers. To model such systems, one possible workaround to avoid
computational bottlenecks is dividing the problem into individual simulations
of separate parts. However, due to hybridization, the spectral characteristics
of the entire system composed of these devices being in the vicinity of each
other can be vastly different from that of the individual components. Hence,
it is imperative to be able to model the entire structure and directly extract
its modes.

A schematic of the system we shall consider is shown in Fig. 5a. It consists
of a 3D cavity containing two dielectric substrates. A transmon qubit is
mounted on each substrate. Each transmon qubit is composed of two
superconducting capacitor pads connected by a Josephson junction. It is a
well-known property of Josephson junctions that their dynamics are well
characterized by the coarse-grained phase across it [61]. In other words, to
an observer outside and away from the junction and who can only make
measurements with a limited precision, the detailed dynamics inside the
junction is unimportant to the dynamics that results from its interaction with
the surrounding electromagnetic environment. This property can be utilized in
coarse-grained calculations using DEC to reduce the complexity of the mesh
needed; this is done by modeling a junction by a single edge instead of finely
meshing the junction geometry. In the calculations of electromagnetic modes,
this is equivalent to the linearization of the junction across the edge. To
obtain modes whose wavelengths are orders of magnitude larger than the
longitudinal size of the junction, the fine details of the material
distribution within the junction are irrelevant and do not need to be resolved
during meshing. In Fig. 5b-5e, a few example modes of this multiscale system
are shown. In this calculation, we have set the dimensions of the cavity to be
$1160\\!\times\\!160\\!\times\\!550$ (mm), while the sizes of the two
identical substrates are both $25\\!\times\\!50\\!\times\\!3$ (mm) and are
separated by a distance of $600$ mm. The size of the capacitor pads of the
qubits is $6\\!\times\\!3\\!\times\\!0.1$ (mm), and the two pads belonging to
the same qubit are connected by a $0.2$ mm-long bridge. We show the results
for different types of modes that the system supports; in Figs. 5b and 5c,
single qubit modes are shown, where the hybridization with the other qubit and
the cavity field is suppressed. On the other hand, there can also be modes in
which both qubits participate, such as the one shown in Fig. 5d. Finally, in
Fig. 5e, we demonstrate a hybridized mode in which both the qubits and the
cavity participate.

### Calculations of open modes

So far in this article, we have considered only closed systems, i.e. systems
that only support modes that decay exponentially beyond their boundaries. We
would also like to extend our method to include the calculations of radiative
modes, which are crucial to the understanding of open systems. The correct
modeling of open systems, where fields can propagate into and out of a
confined domain, has been of great interest since it is directly related to
the quantification of qubit lifetimes and radiative losses. Beyond the realm
of superconducting microwave circuits, the problem has also been important to
the broader scope of electromagnetic devices. We are particularly interested
in the implementation of finite, open boundaries that are transparent so that
fields can propagate through without any reflection. Moreover, the formulation
needs to be compatible with the eventual second quantization of the
electromagnetic field everywhere within the system, including the boundary
itself. Here, we provide the derivations of two such implementations and
present the numerical demonstrations.

#### Green’s boundary integral formulation for scalar fields

Figure 6: Schematic of an open system $D$ consisting of multiple material
regions with refractive indices $n_{i}$. The imaginary boundary
$\partial\mathcal{D}$ (the dashed contour) can have an arbitrary shape as long
as all regions of interest are enclosed within it. The secondary boundary
$\partial\mathcal{D}^{\prime}$ is created by contracting everywhere on the
original boundary $\partial\mathcal{D}$ by an amount $\Delta r$.

First, we discuss an implementation of open boundaries using Green’s boundary
integral formalism. Consider a domain $\mathcal{D}$ consisting of possibly
multiple disjoint regions all enclosed by an imaginary boundary surface.
Before considering the spectral problem associated with the vector Helmholtz
equation 6, we consider the warm-up problem of the Helmholtz equation for a
scalar field $\phi({\mathbf{r}})$:

$\nabla^{2}\phi({\mathbf{r}})+n_{i}^{2}({\mathbf{r}})k^{2}\phi({\mathbf{r}})=0,$
(13)

where $n({\mathbf{r}})$ is the dielectric function of the $i^{\text{th}}$
region $\Gamma_{i}$, and $k$ is again the wavenumber. The Green’s function of
the Helmholtz operator reads

$\displaystyle G({\bf
r,r^{\prime}},k)=\begin{cases}-\frac{i}{4}H_{0}^{(1)}(n_{i}k|{\mathbf{r}}-{\mathbf{r}^{\prime}}|)&\text{in
2D,}\\\
-\frac{e^{in_{i}k|{\mathbf{r}}-{\mathbf{r}^{\prime}}|}}{4\pi|{\mathbf{r}}-{\mathbf{r}^{\prime}}|}&\text{in
3D,}\end{cases}$ (14)

where $H_{0}^{(1)}$ is the Hankel function of the first kind. At first glance,
the Green’s functions in Eq. (14) and their derivatives diverging at
${\mathbf{r}}\\!=\\!{\mathbf{r}^{\prime}}$ might seem like a problem. However,
this difficulty can be circumvented by casting the boundary condition in an
integral form to regularize the Green’s function singularity. Upon applying
Green’s identity and taking the limit as
${\mathbf{r}^{\prime}}\rightarrow\partial\Gamma_{i}$, the field value at a
point ${\mathbf{r}^{\prime}}$ that lies on the boundary $\partial\Gamma_{i}$
of a domain is [62]

$\phi({\bf
r^{\prime}})=2\mathcal{P}^{d-1}\\!\\!\int_{\partial\Gamma_{i}}\\!\\!\big{[}\phi({\bf
r})\nabla G({\bf r,r^{\prime}},k)-G({\bf r,r^{\prime}},k)\nabla\phi({\bf
r})\big{]}\cdot{\mathbf{d}s},$ (15)

where $\mathcal{P}^{d-1}$ indicates the $(d\\!-\\!1)$-dimensional Cauchy
principal value integral, with $d$ being the dimension of the domain. A
detailed derivation of Eq. (15) is given in the Supplementary Material for
completeness. In Eq. (15), the field at a point on the boundary is determined
by the field and its gradient everywhere else on the same boundary. Eq. (15)
is applicable to smooth boundaries of any form and shape. They can also be
made up of different disjoint but closed segments.

To use Eq. (15) for applying boundary condition at the domain boundary
$\partial\mathcal{D}$, we need to numerically evaluate the normal gradient of
the boundary field in a way that is self-contained within the field living on
the boundary itself. To do so, consider a secondary boundary
$\partial\mathcal{D}^{\prime}$ that is formed by contracting
$\partial\mathcal{D}$ by an amount $\Delta r$ along the normal direction
everywhere on the surface, as shown schematically in Fig. 6. There is now a
thin layer bounded by $\partial\mathcal{D}\cup\partial\mathcal{D}^{\prime}$
over which Green’s identity can be performed to obtain Eq. (15) that
determines $\phi({\bf r^{\prime}})$ for every point
${\mathbf{r}^{\prime}}\in\partial\mathcal{D}$. The surface integral is now
done over $\partial\mathcal{D}\cup\partial\mathcal{D}^{\prime}$. Note that
during the evaluation of Eq. (15) the unit normal vector $\hat{n}$ on
$\partial\mathcal{D}^{\prime}$ points into $\mathcal{D}$.

The application of open boundaries using Eq. (15) results in a boundary
condition that depends parametrically on $k$. To solve for the eigenmodes of
Eq. (13), we rewrite it into the following matrix representation

$\big{[}\mathbb{H}_{s}(k)+k^{2}\big{]}\Phi=0,$ (16)

where $\mathbb{H}_{s}(k)$ is the scalar Helmholtz operator, and the vector
$\Phi$ contains the scalar field evaluated at all vertices in the
computational mesh. The task of finding the eigenvalues of $\mathbb{H}_{s}(k)$
can be cast into a singular value decomposition (SVD) problem [63]. The
eigenvalues $k_{\alpha}$ correspond to locations of local minima in the lowest
singular value of $\mathbb{M}_{s}(k)=\mathbb{H}_{s}(k)+k^{2}$.

Figure 7: (a) The lowest singular value of $\mathbb{M}_{s}(k)$ computed for a
sample range of $k$. The physical system is a dielectric disk with radius
$R_{d}=5$ mm and $n=1.5$. A schematic of the system is given in the top-left
corner. (b) The eigenvalues $k_{n}$ of $\mathbb{H}_{s}(k)$, are indicated by
the locations of the local minima of the lowest singular value of
$\mathbb{M}_{s}(k)$. The semi-analytical results are plotted in blue, while
the numerically computed values are plotted in orange. Figure 8: A selection
of open modes of the scalar Helmholtz equation on a dielectric disk whose
radius is $R_{d}=5$ mm and refractive index is $n=1.5$. In each panel: the
surface plot shows the distribution of $Re(\phi)$ in the entire domain inside
the boundary radius $R$. The blue curve shows the value of $Re(\phi)$ at a
fixed angle $\theta=4\pi/5$ as a function of the radial coordinate $r$. The
units for the horizontal axes are in mm. The orange vertical line indicates
the edge of the disk. The red curve shows the field at a fixed radius $r=4$ as
a function of the angular coordinate $\theta$.

We demonstrate the formulation presented here through the calculations of the
open modes of Eq. (13) applied to a dielectric disk placed in a vacuum. The
disk has radius $R_{d}\\!=\\!5$ mm and $n=1.5$, while the boundary is chosen
to be a concentric circle with radius $R\\!=\\!8$ mm. The results for the
lowest singular value of $\mathbb{M}_{s}(k)$ is shown in Fig. 7a for $k$
within the range $0\leq Re(kR_{d})\leq 4\pi$ and $-3.5\pi\leq Im(kR_{d})<0$.
Since the system supports incoming and outgoing modes equally, we only need to
focus on the outgoing modes, whose real part of $k$ is positive. The
distribution of the eigenvalues has a mirror symmetry across $Re(k)=0$, and
the incoming modes are only different from the outgoing ones by the sign of
$Re(k)$. The local minima in Fig. 7a, whose locations correspond to the
eigenvalues $k_{\alpha}$, are collected and plotted in Fig. 7b. We also
compare these numerically computed eigenvalues with the semi-analytical
solutions obtained by solving the transcendental equation arising from
matching the continuity condition of the field and its normal derivative at
the edge of the disk. As seen in Fig. 7b, the numerically obtained eigenvalues
(orange) agree with those obtained through the analytical equation (blue) over
a wide range of $k$.

The field distributions of a few eigenmodes are presented in Fig. 8, where the
real part of the fields, $Re(\phi)$, are plotted. The numerical labeling of
the modes are done first in ascending order of $Im(k_{\alpha}R_{d})$ within
the range $\\{-3.5\pi,0\\}$, then in ascending order of $Re(k_{\alpha}R_{d})$
in the range $\\{0,4\pi\\}$. The features exhibited by these distributions can
be understood through analytical considerations; in polar coordinates, one can
write the fundamental solution to the scalar Helmholtz equation as
$\phi(r,\theta)\\!=\\!R(r)\Theta(\theta)$. The angular dependence of the
fields is of the form $\Theta(\theta)\sim e^{\pm im\theta}$, which is
confirmed by the curves in red in Fig. 8, where the angular dependence of the
fields are plotted at a fixed radius. The radial component, on the other hand,
is given by $R(r)\sim Z_{m}(nkr)$, where $Z_{m}$ is either the Bessel function
of the first kind $J_{m}$, Bessel function of the second kind $Y_{m}$, or the
Hankel functions $H^{(1,2)}_{m}$. Inside the dielectric disk, the radial
component of the field is given by $J_{m}$, while the field outside the disk
is described by $H_{m}^{(1)}$ ($H_{m}^{(2)}\big{)}$ for an outgoing (incoming)
wave. This radial dependence is shown in the blue curves in Fig. 8.

#### Green’s boundary integral formulation for vector fields

A problem of greater interest is, however, the implementation of open
boundaries for vector fields. For this, we have extended Green’s boundary
integral method for scalar fields to address the vector Helmholtz problem as
well. Consider a divergence-free vector field
$\bm{\mathcal{A}}({\mathbf{r}^{\prime}})$ that satisfies Eq. (6) defined over
the domain $\mathcal{D}$. The dependence of the field value at a location
${\mathbf{r}^{\prime}}\in\partial D$ on the field everywhere else on the same
boundary is given by

$\displaystyle\bm{\mathcal{A}}({\mathbf{r}^{\prime}})=-2$
$\displaystyle\mathcal{P}^{d-1}\int_{\partial
D}\bigg{\\{}G({\mathbf{r}},{\mathbf{r}^{\prime}},k)\big{[}({\mathbf{\nabla}}\times\bm{\mathcal{A}})\times\hat{n}\big{]}$
(17) $\displaystyle-\nabla
G({\mathbf{r}},{\mathbf{r}^{\prime}},k)(\bm{\mathcal{A}}\cdot\hat{n})-\nabla
G\times({\bm{\mathcal{A}}\times\hat{n}})\bigg{\\}}ds,$

where $G({\mathbf{r}},{\mathbf{r}^{\prime}},k)$ is again the scalar Green’s
function as defined in Eq. (14), and the principal value integral is
generalized to the vector case. A detailed derivation of Eq. (17) is given in
the Supplementary Material.

Figure 9: Open-mode calculations of the vector Helmholtz problem for a
dielectric disk placed in a vacuum. The disk has radius $R_{d}=5$ mm and
$n=1.5$. (a) A schematic of the polar grid used in the calculation. (b) The
surface plot of the lowest singular value of the operator $\mathbb{M}_{v}(k)$
as a function of $k$. The local minima indicate the eigenvalues of the vector
Helmholtz operator $\mathbb{H}_{v}(k)$. (c) The distributions of
$Re(\bm{\mathcal{A}})$ in a number of radiative modes (the units for all axes
are in mm).

We demonstrate the implementation of the boundary condition in Eq. (17) in DEC
by calculating the modes of Eq. (6) applied to the same dielectric disk
studied in the scalar case. Consider a discretization of the computational
domain into a polar grid that has $N_{R}$ and $N_{\theta}$ vertices along each
radial ray and circle, respectively. In the case of open BC, it is natural to
allow the mesh to also be “open” such that there are additional edges on the
physical boundary of the system that protrude outwards in the direction normal
to the surface, as shown in Fig. 9a. A primal edge $e$ in this setup is
identified by $n_{r}$ and $n_{\theta}$, its radial and angular indices,
accompanied by the superscript $\\{r,\theta\\}$ that indicates the direction
of the edge. In DEC, applying BC to a vector field translates to the
application of BC for the normal and tangential edge fields
$\Phi(e^{r}_{N_{R},n_{\theta}})$ and $\Phi(e^{\theta}_{N_{R},n_{\theta}})$
living on the boundary edges. The condition in Eq. (17) then becomes

$\displaystyle\Phi(e^{r}_{N_{R},n_{\theta}})=-\Phi(e^{r}_{N_{R}-1,n_{\theta}})+\sum_{\theta}2\bigg{\\{}\frac{\partial
G}{\partial
r}\Big{[}\Phi^{r}_{N_{R}+1}(\theta)+\Phi_{N_{R}}^{r}(\theta)\Big{]}R\Delta\theta-\frac{\partial
G}{\partial\theta}\Big{[}\Phi^{\theta}_{n_{\theta}-1}(R)+\Phi^{\theta}_{n_{\theta}}(R)\Big{]}\frac{\Delta
R}{R}\bigg{\\}},$ (18)

for the radial boundary edge, and

$\displaystyle\Phi(e^{\theta}_{N_{R},n_{\theta}})=-\Phi(e^{\theta}_{N_{R},n_{\theta}-1})-\sum_{\theta}$
$\displaystyle\bigg{\\{}2\Delta\theta\bigg{[}G({\mathbf{r}},{\mathbf{r}^{\prime}},k)\Big{(}1+\frac{R}{\Delta
R}\Big{)}-R\frac{\partial G}{\partial
r}\bigg{]}\Big{(}\Phi_{n_{\theta}}^{\theta}(R)+\Phi_{n_{\theta}-1}^{\theta}(R)\Big{)}$
(19) $\displaystyle-2G\frac{R\Delta\theta}{\Delta
R}\Big{(}\Phi_{n_{\theta}}^{\theta}(R-\Delta
R)+\Phi_{n_{\theta}-1}^{\theta}(R-\Delta
R)\Big{)}-2R\frac{\Delta\theta^{2}}{\Delta R}\frac{\partial
G}{\partial\theta}\Big{[}\Phi_{N_{R}+1}^{r}(\theta)+\Phi_{N_{R}}^{r}(\theta)\Big{]}$
$\displaystyle-G\frac{R\Delta\theta}{\Delta
R}\Big{[}\Phi_{N_{R}+1}^{r}(\theta+\Delta\theta)+\Phi_{N_{R}}^{r}(\theta+\Delta\theta)-\Phi_{N_{R}+1}^{r}(\theta-\Delta\theta)-\Phi_{N_{R}}^{r}(\theta-\Delta\theta)\Big{]}\bigg{\\}}$

for the angular boundary edges, where $\Delta R$ and $\Delta\theta$ are the
spacing between consecutive grid points on a radial ray and circle,
respectively. Similar to the scalar problem, to search for the eigenvalues, we
compute the lowest singular values of the operator
$\mathbb{M}_{v}(k)=\mathbb{H}_{v}(k)-k^{2}$, where $\mathbb{H}_{v}(k)$ is the
vector Helmholtz operator. The distribution of the local minima of
$\mathbb{M}_{v}(k)$ in $k$-space is shown in Fig. 9b. The field distributions
of a few randomly selected eigenmodes are presented in Fig. 9c, where the
plotted vector fields exhibit the correct behavior of how radiations permeate
and escape a dielectric disk.

Figure 10: Calculations of radiative modes of a transmon qubit sandwiched
between two cylindrical superconducting capacitors. (a) The schematic of the
system considered; the qubit is mounted on a circular dielectric disk that has
radius $R_{d}=4$ mm, refractive index $n=1.5$. The qubit is composed of two
superconducting islands that are $0.75$ mm away from each other and connected
by a bridge containing a junction in the middle. The width of both the bridge
and the junction is 50$\mu$m. (b) An example of a “dipole” mode; the qubit
acts like a dipole with the field lines starting from the surface of one
island and ending on the surface of another island. (c) An example of a mode
where the field lines flow around the individual capacitor islands with no
normal component at each island’s surface. The boundary of the qubit is
painted in yellow in the enlarged insets, and all the axes of the plots are in
units of mm.

To demonstrate the versatility of DEC in computing radiative modes of
arbitrarily shaped systems, we apply the formulation to the modeling of a
superconducting qubit chip placed in the region between two cylindrical
capacitors, as schematically shown in Fig. 10a. The qubit is made of two
superconducting islands that are separated by a distance of $0.75$ mm and
connected by a bridge that contains a Josephson junction. The width of both
the bridge and the junction is $50\mu$m, and the qubit is mounted on a
circular dielectric disk that has a radius $R_{d}=4$ mm with refractive index
$n=1.5$. We assume the two ideal superconducting cylindrical capacitors extend
indefinitely towards both ends, and the qubit is placed in the middle of the
slit separating them. The width $h$ of the vacuum spacing between the two
cylinders is taken to be very small compared to their radius $R_{c}$. This
allows us to effectively discretize the distance between the surfaces of the
two cylinders by two edges. Due to mirror symmetry at $z=0$, the amplitude of
the coarse-grained field on one $z$-edge is identical to that of its mirroring
edge. This allows us to decouple the in-plane field component from the out-of-
plane component, and we can focus on solving for the modes of Eq. (9) applied
to the 2D in-plane field on the $z=0$ slice where the substrate containing the
qubit is located. We choose the computational boundary to be a circle whose
radius $R=8$ mm is half that of the cylindrical capacitors $R_{c}$. There are
two notable types of modes, whose examples are shown in Figs. 10b-10c. Due to
the geometry of the qubit being two islands connected by a bridge, it can act
as a dipole. A mode that behaves this way is shown in Fig. 10b, where the
field lines exit in the normal direction from one island and enter the other
island. The field also decays in the bulks of the islands at the rate
determined by their penetration depths. In another scenario, due to their
separation, the two islands can behave as individual superconducting objects
around which the field lines flow with no normal component at the surfaces of
the islands. A demonstration of such a mode is shown in Fig. 10c.

As was mentioned earlier, the formulation for calculating radiative fields
using Green’s function method is flexible in terms of the shape and topology
of boundary surfaces it can be applied. The physical boundary of the system
can be made of multiple, possibly disjoint, but closed segments that bound a
non-simply connected structure. However, the vector field being calculated has
to be divergent-free as a prerequisite. Although this restriction is not a
concern in many useful cases, such as the calculation of the electric field in
a source-less region or of the magnetic vector potential in the Coulomb gauge,
it is sometimes helpful to have the freedom of not necessarily choosing a
divergence-less field. In the following section, we introduce an alternative
formulation that achieves this.

#### Vector spherical harmonics expansion

In this section, we discuss the implementation of open boundaries using vector
spherical harmonics (VSH) expansions. The VSH [64] is an extension for vector
fields of the perhaps more well-known scalar spherical harmonics. There are
three sub-classes of these vectors, defined in spherical coordinates
$(r,\theta,\varphi)$ as

$\displaystyle{\mathbf{Y}_{lm}}$ $\displaystyle=Y_{lm}(\theta,\varphi)\hat{r}$
(20) $\displaystyle{\mathbf{\Psi}_{lm}}$ $\displaystyle=r\nabla
Y_{lm}(\theta,\varphi)$ (21) $\displaystyle{\mathbf{\Phi}_{lm}}$
$\displaystyle=\mathbf{r}\times\nabla Y_{lm}(\theta,\varphi)$ (22)

that altogether form an orthonormal and complete basis. Any vector field can
then be expanded as follows

$\displaystyle\bm{\mathcal{A}}({\mathbf{r}})=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}\mathcal{A}^{r}_{lm}{\mathbf{Y}_{lm}}+\mathcal{A}_{lm}^{(1)}{\mathbf{\Psi}_{lm}}+\mathcal{A}_{lm}^{(2)}{\mathbf{\Phi}_{lm}},$
(23)

where coefficients $\mathcal{A}^{r}_{lm},\mathcal{A}_{lm}^{(1)}$ and
$\mathcal{A}_{lm}^{(2)}$ are functions of the radial coordinate $r$. In Eq.
(23) above, of the three mutually orthogonal terms, the first term on the
right-hand side is the radial component of the field, while the remaining two
terms are angular contributions. Note that neither of the two angular terms
aligns with azimuthal or polar directions but is a linear combination of both.
Using Eq. (23) as an ansatz for the field, Eq. (6) can then be split into
three decoupled ODEs, each for one of the three coefficients. Their solutions
are

$\displaystyle\mathcal{A}_{lm}^{r}=$
$\displaystyle\frac{Z_{l+1/2}(kr)}{r^{3/2}},$ (24)
$\displaystyle\mathcal{A}_{lm}^{(1)}=$
$\displaystyle\frac{1}{\l(l+1)}\bigg{\\{}\frac{Z_{l+1/2}(kr)}{r^{3/2}}+\frac{k}{\sqrt{r}}\Big{[}Z_{l-1/2}(kr)$
$\displaystyle\hskip 72.26999pt-Z_{l+3/2}(kr)\Big{]}\bigg{\\}},$ (25)
$\displaystyle\mathcal{A}_{lm}^{(2)}=$
$\displaystyle\frac{\pi}{2kr}Z_{l+1/2}(kr)$ (26)

where $Z$ in general can be the Bessel function of the first kind $J$, the
Bessel function of the second kind $Y$, or the Hankel functions $H^{(1,2)}$.
For the outgoing(incoming) waves, the Hankel functions of the first(second)
kind provide the appropriate description. After the VSH expansion in Eq. (23),
the implementation of the open boundary condition on the field
$\bm{\mathcal{A}}$ now reduces to the BC on each of the three field components
$\bm{\mathcal{A}^{r}},\bm{\mathcal{A}^{(1)}}$ and $\bm{\mathcal{A}^{(2)}}$.
Here we discuss the BC for $\mathcal{A}^{r}$ as an example, while the details
about BCs for $\mathcal{A}^{(1)}$ and $\mathcal{A}^{(2)}$ are given in the
Supplementary Material. Consider a spherical boundary of radius $R$. The
series expansion for the radial component of the radiating field on the sphere
is given by

$\displaystyle\bm{\mathcal{A}}^{r}(R,\theta,\varphi)$
$\displaystyle=\sum_{l,m}\tilde{a}_{lm}^{r}A_{lm}^{r}(R){\mathbf{Y}_{lm}}(R,\theta,\varphi)$
$\displaystyle=\sum_{l,m}\tilde{a}_{lm}^{r}\frac{H_{l+1/2}(kR)}{R^{-3/2}}Y_{lm}\hat{r},$
(27)

where $\tilde{a}_{lm}^{r}$ are the expansion coefficients. Performing an
integral over the surface of the sphere and utilizing the orthonormality
condition of VSH, we arrive at an expression for the expansion coefficients as
follows

$\displaystyle a^{r}_{lm}$
$\displaystyle\equiv\tilde{a}_{lm}^{r}H_{l+1/2}(kR)$ (28)
$\displaystyle=R^{3/2}\int
d\Omega\bm{\mathcal{A}}^{r}\cdot{\mathbf{Y}^{*}_{lm}}$
$\displaystyle=R^{-1/2}\sum_{i,j}\mathcal{A}^{r}(R,\theta_{i},\varphi_{j})Y_{lm}^{*}(\theta_{i},\varphi_{j})\Delta
A(f),$

where $\int d\Omega$ is the solid angle integral that covers the entire
sphere, and the last line of Eq. (28) above is the discretized version of the
integral written as a sum over all the triangular patches $f$ whose centers
are located at $(\theta_{i},\varphi_{j})$ on the spherical boundary. To apply
the BCs, consider the sphere beneath the boundary that has radius
$R\\!-\\!\Delta r$ (i.e. the second-to-last layer). The radial derivative of
$\mathcal{A}^{r}$ is given by

$\displaystyle\frac{\mathcal{A}^{r}(R,\theta,\varphi)-\mathcal{A}^{r}(R-\Delta
R,\theta,\varphi)}{\Delta R}=$ (29)
$\displaystyle\sum_{l,m}\tilde{a}^{r}_{lm}\bigg{[}\frac{-3}{2R^{5/2}}H_{l+1/2}(kR)+\frac{k}{R^{3/2}}H^{\prime}_{l+1/2}(kR)\bigg{]}Y_{l,m},$

which allows us to write the field at the boundary
$\mathcal{A}^{r}(R,\theta,\varphi)$ in terms of the expansion coefficients
$\tilde{a}^{r}_{lm}$ and the field at the layer with radius $R\\!-\\!\Delta
R$. Using Eq. (29) and the expression for $\tilde{a}^{r}_{lm}$ derived in Eq.
(28), the boundary condition for $\mathcal{A}^{r}$ is then

$\displaystyle\mathcal{A}^{r}$
$\displaystyle(R,\theta,\varphi)\\!=\\!\sum_{l.m}\\!\sum_{i,j}\\!\mathcal{A}^{r}(R\\!-\\!\Delta
R,\theta_{i},\varphi_{j})Y^{*}_{lm}(\theta_{i},\varphi_{j})\frac{\Delta
A(f)}{R^{2}}$ $\displaystyle\times\bigg{\\{}1+\Delta
R\Big{[}-\frac{3}{2R}+k\frac{H_{l-1/2}(kR)-H_{l+3/2}(kR)}{2H_{l+1/2}(kR)}\Big{]}\bigg{\\}}$
$\displaystyle\hskip 86.72377pt\times Y_{lm}(\theta,\varphi).$ (30)

The expression in Eq. (Vector spherical harmonics expansion) is for the radial
component of the field, which in DEC is a scalar living on vertices of the
primal mesh. To arrive at the appropriate BC for the edge field $\Phi(e)$
governed by the discrete vector Helmholtz equation given in Eq. (9), we need
to rewrite $\mathcal{A}^{r}(R,\theta,\varphi)$ and
$\mathcal{A}^{r}(R\\!-\\!\Delta R,\theta_{i},\varphi_{j})$ in terms of the
edges belonging to the vertex located at $(R,\theta,\varphi)$. This means two
consecutive spherical layers beneath the boundary are needed to compute
$\mathcal{A}^{r}(R\\!-\\!\Delta R,\theta_{i},\varphi_{j})$. The procedure as a
whole can be substantially simplified by designing the mesh so that the
boundary surface and the two layers beneath it have identical triangulation
patterns. This ensures that for any vertex $v_{b}$ located at
$(R,\theta,\varphi)$ on the boundary surface, there is an edge $e^{r}$
oriented radially that connects $v_{b}$ with the vertex $v_{b}^{\prime}$ at
$(R\\!-\\!\Delta R,\theta,\varphi)$ that has the same angular coordinates as
$v_{b}$. Similarly, there is always an edge connecting $v_{b}^{\prime}$ with a
vertex $v_{b}^{\prime\prime}$ at $(R\\!-\\!2\Delta R,\theta,\varphi)$. Note
that we only need the three outmost layers of the mesh to be spherical and
have identical triangulation patterns, while the rest of the internal domain
can be freely and randomly discretized with tetrahedra. Similar to the
discussion on Green’s function approach in the previous section, here, for
open boundaries, it is convenient to use an open mesh such that there are
radial edges that protrude out of the boundary surface. The BC applied to
these boundary radial edges $e^{r}_{b}$ is then

$\displaystyle\Phi(e^{r}_{b})=-\Phi(e^{r^{\prime}}_{b})+\sum_{e^{r}}\sum_{f\supset
v|v\subset e^{r}}\sum_{l,m}\frac{\Delta
A(f)}{3R^{2}}\Phi(e^{r})\Bigg{\\{}1+\Delta
R\bigg{[}-\frac{3}{2R}+k\frac{H_{l-1/2}(kR)-H_{l+3/2}(kR)}{2H_{l+1/2}(kR)}\bigg{]}Y^{*}_{lm}(f)Y_{lm}(\theta,\varphi)\Bigg{\\}},$
(31)

where $e^{r^{\prime}}_{b}$ is the other radial edge that shares the boundary
vertex with $e_{b}^{r}$. In Eq. (31), the first sum is done over all radial
edges $e^{r}$ living in the two layers beneath the boundary surface, while the
second sum is done over all faces $f$ that share a vertex with $e^{r}$. This
concludes our derivation for the open BC applied to a radial edge in the DEC
framework. The boundary conditions for edges lying tangential to the surface
are discussed in the Supplementary Material.

Figure 11: Results for the open-boundary modes of a dielectric microsphere.
(a) Schematic of the system studied: a sphere of radius $R_{d}=12\mu$m. (b)
The lowest singular value of $\mathbb{M}_{v}(k)$ as a function of $k$. (c)
Error of the lowest eigenvalue as a function of the largest ordering of Bessel
terms. (d) Error of the lowest eigenvalue as a function of the boundary radius
$R$.

We apply the formulation developed here to calculate the modes for a three-
dimensional dielectric microsphere surrounded by a vacuum. The radius of this
dielectric sphere is kept fixed at $R_{d}=12\mu$m throughout all calculations
discussed in this section. The computational boundary is an “imaginary”
spherical surface at a finite radius $R>R_{d}$, as schematically shown in Fig.
11a. The results of the eigenvalue search found through singular value
calculations are shown in Fig. 11b, where we have chosen $R=24\mu$m. With a
primal mesh of $\approx 2900$ vertices and $\approx 18000$ edges, the local
minima of the smallest singular value are of order $10^{-7}$, proving them to
be reliable indicators of the eigenvalues of the system. The dielectric
microsphere is chosen to demonstrate our formulation because there exist semi-
analytical solutions that can be obtained by solving the transcendental
equation arising from enforcing the continuity of the field and its normal
derivative at the dielectric-vacuum interface. We utilize these solutions to
investigate the efficacy of our method in a number of ways. First, the
numerical implementation of the series expansion in Eq. (Vector spherical
harmonics expansion) requires a truncation in the number of terms considered.
We, therefore, investigate the dependence of convergence rate on the number of
Bessel functions included while keeping the mesh size fixed. The results can
be seen in Fig. 11c, where we compute the error of the numerically obtained
lowest eigenvalue with the semi-analytical solution. With a manageable number
of $L_{max}=10$, where $L_{max}$ is the largest ordering of Bessel terms in
the expansion, the error reduces to below $4\%$. Another important factor is
how large the computational boundary needs to be in order to achieve precise
solutions since accuracy is expected to improve as the boundary is moved
further away from where the devices are concentrated, but doing so also
increases mesh complexity, and hence the computational burden is exacerbated.
In Fig. 11d, the convergence as a function of boundary radius $R$ of the
lowest eigenvalue is plotted. We see that even with a tight boundary when
$R=16\mu$m (while the dielectric sphere is kept at $12\mu$m) the error is
relatively low at $15\%$, and for $R=2R_{d}$ it is reduced to $<4\%$. These
results verify that our coarse-grained implementation of open boundaries using
VSH is able to produce accurate solutions using a boundary of reasonable size.

The formulation based on VSH expansions introduced here allows for effective
calculations of open modes while also removing the requirement of divergence-
less fields. However, unlike Green’s boundary integral formulation, which is
adaptable to arbitrarily-shaped boundaries, this approach needs to be
implemented on a spherical surface. This makes the formulation most ideal when
applied to three-dimensional structures that have comparable lengths along
three directions, such as 3D cavities. In some other specialized cases, this
requirement may limit us from drawing the tightest boundary possible, which in
turn might cause the computational domain to be larger than it needs to be. An
example could be a coplanar waveguide whose length is much larger than its
planar width, which in turn is much larger than the thickness along the cross-
section of the waveguide. In such a case, replacing a spherical boundary with
an ellipsoidal one is preferable since an ellipsoid has three tuning
parameters that allow squeezing and stretching in three directions.
Fortunately, the series expansion-based approach discussed here is extendable
to ellipsoidal coordinates. This relies on the fact that the Helmholtz
equation is separable in ellipsoidal coordinates as well, which allows for the
expansion of functions, now in terms of ellipsoidal waves. Although the
mathematically intensive discussion on how to generalize our method to
ellipsoidal coordinates is beyond the scope of this article, in the
Supplementary Material, we provide a preamble by discussing a brief proof of
the separability of the Helmholtz equation.

## Discussion

In this article, we have introduced a spectral theory for the modeling of mode
hybridization and relaxation rates in general systems composed of
superconducting and dielectric materials. In doing so, we have demonstrated
the following important points:

(1) We show that the coarse-grained formulation of electrodynamics using
discrete exterior calculus is adaptable to both simplicial and complex meshing
strategies with equal precision. Although based on the calculations of coarse-
grained quantities, the numerical scheme is highly sensitive to small changes
in the geometry and topology of the system, as was demonstrated in the
categorization of degenerate modes arising from hidden symmetries when a small
perturbation is introduced.

(2) The method, which maps the problem of modeling vector fields to the
calculation of averaged projections living on discrete edges, is particularly
effective when applied to systems composed of objects of multiple spatial
scales, such as those used in superconducting electronic circuits. Since the
method keeps track of the average variables, it allows for more efficient
meshing strategies that take into account the sizes of the Josephson junctions
as well as the finite resolutions of the measurement apparatus. This is an
important advantage that this formulation offers because, here, the mesh does
not need to be finer than the size of JJs or the measurement apparatus.

(3) We introduce two implementations of open boundaries for the vector wave
equation that are applicable to a wide range of electromagnetic systems. By
drawing “imaginary,” transparent boundaries that are reasonably sized, we are
able to faithfully produce the open modes that agree well with analytical
solutions. We also demonstrate the flexibility of the formulation by applying
it to the calculations of radiative modes of a superconducting qubit.

The spectral theory presented here is suitable for the second quantization of
the electromagnetic field and the subsequent studies of open quantum systems
with arbitrary complexity. The ability to accurately and efficiently extract
relaxation rates, as well as the quantification of mode hybridization
demonstrated here, is useful to the modeling and optimization of
superconducting circuits.

## Methods

As mentioned in the Results section, the spectral problem of both closed and
open systems are computed through the coarse-graining of Eq. (6).
Operationally, this is obtained by evaluating the integral
$\int_{e}\bm{\mathcal{A}}\\!\cdot\\!{\mathbf{d}\ell}$ along an edge $e$ in the
discretized grid of the domain, which produces Eq. (9).

For systems enclosed by a cavity with boundaries made of either perfect
($\lambda_{L}=0$) or non-perfect superconductors ($\lambda_{L}>0$), the
application of the Dirichlet boundary condition on the field results in the
eigenvalue problem of type $(\mathcal{H}-k^{2})\Phi=0$ (where $\mathcal{H}$ is
independent of $k$) that can be solved using standard eigensolvers.

For the open-boundary problem, both implementations introduced here result in
an equation of the form $\big{[}\mathbb{H}(k)-k^{2}\big{]}\Phi=0$ where
$\mathbb{H}(k)$ depends non-trivially on $k$. To solve for eigenvalues
$k_{\alpha}$, we cast the equation into a singular value problem, in which
$k_{\alpha}$ ideally corresponds to when the lowest singular value of
$\mathbb{M}(k)=\mathbb{H}(k)-k^{2}$ is zero. In reality, due to numerical
artifacts, the lowest singular value of $\mathbb{M}$ never actually goes to
zero. We, therefore, look for local minima during the sweep over the complex
plane of $k$. As shown in the Result section, these local minima are sharp and
can be identified with high precision.

The source code for the calculations in this article is available in the open-
source DEC-QED repository [58].

## Acknowledgements

We gratefully acknowledge discussions with Nicholas Bronn, Thomas G. McConkey,
Anil N. Hirani, Thomas Maldonado, Zoe Zager, and Haley M. Cole. We are
grateful to Benjamin Lienhard and Wentao Fan for giving the paper their
critical read and for the insightful comments. We acknowledge support from the
US Department of Energy, Office of Basic Energy Sciences, Division of
Materials Sciences and Engineering, under Award No. DESC0016011. The
simulations presented in this article were performed on computational
resources managed and supported by Princeton Research Computing, a consortium
of groups including the Princeton Institute for Computational Science and
Engineering (PICSciE) and the Office of Information Technology’s High-
Performance Computing Center and Visualization Laboratory at Princeton
University.

## References

  * Houck  _et al._ [2012] A. A. Houck, H. E. Türeci, and J. Koch, On-chip quantum simulation with superconducting circuits, Nature Physics 8, 292 (2012).
  * Carusotto  _et al._ [2020] I. Carusotto, A. A. Houck, A. J. Kollár, P. Roushan, D. I. Schuster, and J. Simon, Photonic materials in circuit quantum electrodynamics, Nature Physics 16, 268 (2020).
  * Hatridge  _et al._ [2011] M. Hatridge, R. Vijay, D. Slichter, J. Clarke, and I. Siddiqi, Dispersive magnetometry with a quantum limited squid parametric amplifier, Physical Review B 83, 134501 (2011).
  * Macklin  _et al._ [2015] C. Macklin, K. O’brien, D. Hover, M. Schwartz, V. Bolkhovsky, X. Zhang, W. Oliver, and I. Siddiqi, A near–quantum-limited josephson traveling-wave parametric amplifier, Science 350, 307 (2015).
  * Blais  _et al._ [2004] A. Blais, R.-S. Huang, A. Wallraff, S. M. Girvin, and R. J. Schoelkopf, Cavity quantum electrodynamics for superconducting electrical circuits: An architecture for quantum computation, Physical Review A 69, 062320 (2004).
  * Devoret and Schoelkopf [2013] M. H. Devoret and R. J. Schoelkopf, Superconducting circuits for quantum information: an outlook, Science 339, 1169 (2013).
  * Feynman  _et al._ [1965] R. P. Feynman, R. B. Leighton, and M. Sands, The Schrödinger Equation in a Classical Context: A Seminar on Superconductivity, in  _The Feynman Lectures on Physics_ , Vol. III (Addison–Wesley, 1965) Chap. 21.
  * Ao  _et al._ [1995] P. Ao, D. J. Thouless, and X.-M. Zhu, Nonlinear Schrödinger equation for superconductors, Modern Physics Letters B 9, 755 (1995).
  * Aitchison  _et al._ [1995] I. J. Aitchison, P. Ao, D. J. Thouless, and X.-M. Zhu, Effective Lagrangians for BCS superconductors at T=0, Physical Review B 51, 6531 (1995).
  * Pham  _et al._ [2023] D. N. Pham, W. Fan, M. G. Scheer, and H. E. Tureci, Flux-based three-dimensional electrodynamic modeling approach to superconducting circuits and materials, Phys. Rev. A 107, 053704 (2023).
  * Greiter  _et al._ [1989] M. Greiter, F. Wilczek, and E. Witten, Hydrodynamic relations in superconductivity, Modern Physics Letters B 3, 903 (1989).
  * Fisher  _et al._ [1990] M. P. A. Fisher, G. Grinstein, and S. M. Girvin, Presence of quantum diffusion in two dimensions: Universal resistance at the superconductor-insulator transition, Phys. Rev. Lett. 64, 587 (1990).
  * Salasnich [2009] L. Salasnich, Hydrodynamics of Bose and Fermi superfluids at zero temperature: the superfluid nonlinear Schrödinger equation, Laser physics 19, 642 (2009).
  * Blais  _et al._ [2021] A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff, Circuit quantum electrodynamics, Rev. Mod. Phys. 93, 025005 (2021).
  * 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, Black-box superconducting circuit quantization, Phys. Revs. Letts. 108, 240502 (2012).
  * Solgun  _et al._ [2014] F. Solgun, D. W. Abraham, and D. P. DiVincenzo, Blackbox quantization of superconducting circuits using exact impedance synthesis, Physical Review B 90, 134504 (2014).
  * Smith  _et al._ [2016] W. C. Smith, A. Kou, U. Vool, I. M. Pop, L. Frunzio, R. J. Schoelkopf, and M. H. Devoret, Quantization of inductively shunted superconducting circuits, Phys. Rev. B 94, 144507 (2016).
  * Malekakhlagh  _et al._ [2017] M. Malekakhlagh, A. Petrescu, and H. E. Türeci, Cutoff-free circuit quantum electrodynamics, Phys. Rev. Lett. 119, 073601 (2017).
  * Parra-Rodriguez  _et al._ [2019] A. Parra-Rodriguez, I. Egusquiza, D. DiVincenzo, and E. Solano, Canonical circuit quantization with linear nonreciprocal devices, Physical Review B 99, 014514 (2019).
  * Minev  _et al._ [2021] Z. K. Minev, Z. Leghtas, S. O. Mundhada, L. Christakis, I. M. Pop, and M. H. Devoret, Energy-participation quantization of josephson circuits, npj Quantum Information 7, 131 (2021).
  * [21] The time-coarse graining toolbox, https://github.com/leonbello/QuantumGraining.jl.
  * Houck  _et al._ [2008] A. Houck, J. Schreier, B. Johnson, J. Chow, J. Koch, J. Gambetta, D. Schuster, L. Frunzio, M. Devoret, S. Girvin, _et al._ , Controlling the spontaneous emission of a superconducting transmon qubit, Physical review letters 101, 080502 (2008).
  * Reed  _et al._ [2010] M. D. Reed, B. R. Johnson, A. A. Houck, L. DiCarlo, J. M. Chow, D. I. Schuster, L. Frunzio, and R. J. Schoelkopf, Fast reset and suppressing spontaneous emission of a superconducting qubit, Applied Physics Letters 96 (2010).
  * Jeffrey  _et al._ [2014] E. Jeffrey, D. Sank, J. Y. Mutus, T. C. White, J. Kelly, R. Barends, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. Megrant, P. J. J. O’Malley, C. Neill, P. Roushan, A. Vainsencher, J. Wenner, A. N. Cleland, and J. M. Martinis, Fast accurate state measurement with superconducting qubits, Phys. Rev. Lett. 112, 190504 (2014).
  * Bronn  _et al._ [2015] N. T. Bronn, Y. Liu, J. B. Hertzberg, A. D. Córcoles, A. A. Houck, J. M. Gambetta, and J. M. Chow, Broadband filters for abatement of spontaneous emission in circuit quantum electrodynamics, Applied Physics Letters 107 (2015).
  * Boissonneault  _et al._ [2009] M. Boissonneault, J. M. Gambetta, and A. Blais, Dispersive regime of circuit qed: Photon-dependent qubit dephasing and relaxation rates, Phys. Rev. A 79, 013819 (2009).
  * Slichter  _et al._ [2012] D. H. Slichter, R. Vijay, S. J. Weber, S. Boutin, M. Boissonneault, J. M. Gambetta, A. Blais, and I. Siddiqi, Measurement-induced qubit state mixing in circuit qed from up-converted dephasing noise, Phys. Rev. Lett. 109, 153601 (2012).
  * Mundhada  _et al._ [2016] S. Mundhada, S. Shankar, A. Narla, E. Zalys-Geller, S. Girvin, and M. Devoret, Dependence of transmon qubit relaxation rate on readout drive power, APS March Meeting V48 (2016).
  * Minev  _et al._ [2019] Z. Minev, S. Mundhada, S. Shankar, P. Reinhold, R. Gutiérrez-Jáuregui, R. Schoelkopf, M. Mirrahimi, H. Carmichael, and M. Devoret, To catch and reverse a quantum jump mid-flight, Nature 570, 10.1038/s41586-019-1287-z (2019).
  * Petrescu  _et al._ [2020] A. Petrescu, M. Malekakhlagh, and H. E. Türeci, Lifetime renormalization of driven weakly anharmonic superconducting qubits. ii. the readout problem, Phys. Rev. B 101, 134510 (2020).
  * Hanai  _et al._ [2021] R. Hanai, A. McDonald, and A. Clerk, Intrinsic mechanisms for drive-dependent purcell decay in superconducting quantum circuits, Phys. Rev. Res. 3, 043228 (2021).
  * Gusenkova  _et al._ [2021] D. Gusenkova, M. Spiecker, R. Gebauer, M. Willsch, D. Willsch, F. Valenti, N. Karcher, L. Grünhaupt, I. Takmakov, P. Winkel, D. Rieger, A. V. Ustinov, N. Roch, W. Wernsdorfer, K. Michielsen, O. Sander, and I. M. Pop, Quantum nondemolition dispersive readout of a superconducting artificial atom using large photon numbers, Phys. Rev. Appl. 15, 064030 (2021).
  * Shillito  _et al._ [2022] R. Shillito, A. Petrescu, J. Cohen, J. Beall, M. Hauru, M. Ganahl, A. G. Lewis, G. Vidal, and A. Blais, Dynamics of transmon ionization, Phys. Rev. Appl. 18, 034031 (2022).
  * Cohen  _et al._ [2023] J. Cohen, A. Petrescu, R. Shillito, and A. Blais, Reminiscence of classical chaos in driven transmons, PRX Quantum 4, 020312 (2023).
  * Caldeira and Leggett [1983] A. Caldeira and A. Leggett, Quantum tunnelling in a dissipative system, Annals of Physics 149, 374 (1983).
  * Filipp  _et al._ [2011] S. Filipp, M. Göppl, J. M. Fink, M. Baur, R. Bianchetti, L. Steffen, and A. Wallraff, Multimode mediated qubit-qubit coupling and dark-state symmetries in circuit quantum electrodynamics, Phys. Rev. A 83, 063827 (2011).
  * Gely  _et al._ [2017] M. F. Gely, A. Parra-Rodriguez, D. Bothner, Y. M. Blanter, S. J. Bosman, E. Solano, and G. A. Steele, Convergence of the multimode quantum rabi model of circuit quantum electrodynamics, Phys. Rev. B 95, 245115 (2017).
  * Sinha  _et al._ [2022] K. Sinha, S. A. Khan, E. Cüce, and H. E. Türeci, Radiative properties of an artificial atom coupled to a josephson-junction array, Phys. Rev. A 106, 033714 (2022).
  * Yoshihara  _et al._ [2014] F. Yoshihara, Y. Nakamura, F. Yan, S. Gustavsson, J. Bylander, W. D. Oliver, and J.-S. Tsai, Flux qubit noise spectroscopy using rabi oscillations under strong driving conditions, Physical Review B 89, 020503 (2014).
  * Kumar  _et al._ [2016] P. Kumar, S. Sendelbach, M. Beck, J. Freeland, Z. Wang, H. Wang, C. Y. Clare, R. Wu, D. Pappas, and R. McDermott, Origin and reduction of 1/f magnetic flux noise in superconducting devices, Physical Review Applied 6, 041001 (2016).
  * Wilen  _et al._ [2021] C. D. Wilen, S. Abdullah, N. Kurinsky, C. Stanford, L. Cardani, G. d’Imperio, C. Tomei, L. Faoro, L. Ioffe, C. Liu, _et al._ , Correlated charge noise and relaxation errors in superconducting qubits, Nature 594, 369 (2021).
  * Pop  _et al._ [2014] I. M. Pop, K. Geerlings, G. Catelani, R. J. Schoelkopf, L. I. Glazman, and M. H. Devoret, Coherent suppression of electromagnetic dissipation due to superconducting quasiparticles, Nature 508, 369 (2014).
  * Frattini  _et al._ [2018] N. Frattini, V. Sivak, A. Lingenfelter, S. Shankar, and M. Devoret, Optimizing the nonlinearity and dissipation of a snail parametric amplifier for dynamic range, Physical Review Applied 10, 054020 (2018).
  * Wenner  _et al._ [2011] J. Wenner, M. Neeley, R. C. Bialczak, M. Lenander, E. Lucero, A. D. O’Connell, D. Sank, H. Wang, M. Weides, A. N. Cleland, and J. M. Martinis, Wirebond crosstalk and cavity modes in large chip mounts for superconducting qubits, Supercond. Sci. Technol. 24, 065001 (2011).
  * Huang  _et al._ [2021] S. Huang, B. Lienhard, G. Calusine, A. Vepsäläinen, J. Braumüller, D. K. Kim, A. J. Melville, B. M. Niedzielski, J. L. Yoder, B. Kannan, _et al._ , Microwave package design for superconducting quantum processors, PRX Quantum 2, 020306 (2021).
  * Berenger [1994] J.-P. Berenger, A perfectly matched layer for the absorption of electromagnetic waves, Journal of computational physics 114, 185 (1994).
  * Berenger [1996] J.-P. Berenger, Three-dimensional perfectly matched layer for the absorption of electromagnetic waves, Journal of computational physics 127, 363 (1996).
  * Clayton and Engquist [1977] R. Clayton and B. Engquist, Absorbing boundary conditions for acoustic and elastic wave equations, Bulletin of the seismological society of America 67, 1529 (1977).
  * Higdon [1987] R. L. Higdon, Numerical absorbing boundary conditions for the wave equation, Mathematics of Computation 49, 65 (1987).
  * Gell-Mann and E. [1954] M. Gell-Mann and L. F. E., Quantum electrodynamics at small distances, Phys. Rev. 95, 1300 (1954).
  * Qi  _et al._ [2018] B. Qi, L. Zhang, and L. Ge, Defect states emerging from a non-hermitian flatband of photonic zero modes, Physical review letters 120, 093901 (2018).
  * El-Ganainy  _et al._ [2018] R. El-Ganainy, K. G. Makris, M. Khajavikhan, Z. H. Musslimani, S. Rotter, and D. N. Christodoulides, Non-hermitian physics and pt symmetry, Nature Physics 14, 11 (2018).
  * Gigli  _et al._ [2020] C. Gigli, T. Wu, G. Marino, A. Borne, G. Leo, and P. Lalanne, Quasinormal-mode non-hermitian modeling and design in nonlinear nano-optics, ACS photonics 7, 1197 (2020).
  * Purcell  _et al._ [1946] E. M. Purcell, H. C. Torrey, and R. V. Pound, Resonance absorption by nuclear magnetic moments in a solid, Physical review 69, 37 (1946).
  * Malekakhlagh  _et al._ [2016] M. Malekakhlagh, A. Petrescu, and H. E. Türeci, Non-markovian dynamics of a superconducting qubit in an open multimode resonator, Physical Review A 94, 063848 (2016).
  * Bosman  _et al._ [2017] S. J. Bosman, M. F. Gely, V. Singh, A. Bruno, D. Bothner, and G. A. Steele, Multi-mode ultra-strong coupling in circuit quantum electrodynamics, npj Quantum Information 3, 46 (2017).
  * Puertas Martínez  _et al._ [2019] J. Puertas Martínez, S. Léger, N. Gheeraert, R. Dassonneville, L. Planat, F. Foroughi, Y. Krupko, O. Buisson, C. Naud, W. Hasch-Guichard, _et al._ , A tunable josephson platform to explore many-body quantum optics in circuit-qed, npj Quantum Information 5, 19 (2019).
  * [58] DEC-QED toolbox, https://github.com/dnpham23/DEC-QED.
  * London and London [1935] F. London and H. London, The electromagnetic equations of the supraconductor, Proceedings of the Royal Society of London. Series A - Mathematical and Physical Sciences 149, 71 (1935).
  * Shaw [1974] G. B. Shaw, Degeneracy in the particle-in-a-box problem, Journal of Physics A: Mathematical, Nuclear and General 7, 1537 (1974).
  * Josephson [1965] B. D. Josephson, Supercurrents through barriers, Advances in Physics 14, 419 (1965).
  * Wiersig [2002] J. Wiersig, Boundary element method for resonances in dielectric microcavities, Journal of Optics A: Pure and Applied Optics 5, 53 (2002).
  * Bäcker [2003] A. Bäcker, Numerical aspects of eigenvalue and eigenfunction computations for chaotic quantum systems, in  _The Mathematical Aspects of Quantum Maps_ , edited by M. D. Esposti and S. Graffi (Springer Berlin Heidelberg, Berlin, Heidelberg, 2003) pp. 91–144.
  * Barrera  _et al._ [1985] R. G. Barrera, G. Estevez, and J. Giraldo, Vector spherical harmonics and their application to magnetostatics, European Journal of Physics 6, 287 (1985).