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

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.

I 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].

With reference to the prior work on modeling of superconducting devices, our approach is within the class of full-wave methods [58, 59]. With respect to existing full-wave methods for superconducting systems, DEC-QED has two important distinctions: 1. Calculation of Radiative Losses: The boundary conditions implemented in existing approaches to capture radiative losses often implement absorbing boundary conditions. Most frequently, radiative losses at a 3D cavity’s ports are modeled by adapting from the lumped-element models and terminating the ports with $50\,\Omega$ matched resistors [16] or by covering the holes on the cavity walls with resistive disks [20]. Although this approach is effective in obtaining the frequencies of the eigenemodes, it is known to fail to accurately capture the spatial structure of the open modes. More importantly, no known method exists to second quantize using such modes, and therefore such EM full wave simulations can not rigorously be used to synthesize a quantum mechanical Hamiltonian or Liouvillian. Often the quantum Hamiltonian is derived from simulations ignoring radiative losses, and the relaxation rates are introduced through a perturbative approach. In contrast, DEC-QED builds on prior work in classical electromagnetic systems for the implementation of correct open (transparent) boundary conditions [60, 61], which recently has been shown to serve as a mathematically consistent basis for second quantization of the modes of a transmission line cavity [38]. 2. Time-dependent Evolution: DEC-QED aims at numerically solving the set of non-linear PDEs describing the electrohydrodynamics of the superconducting condensate coupled to fields described by Maxwell equations. This allows DEC-QED to describe the evolution of gauge-invariant hybridized fields living throughout the free space and the superconductors, forgoing issues that arise in defining boundary conditions at superconductor/vacuum boundaries whose accuracy is not controlled. Compared to existing full-wave methods, this is a unique feature of our approach. The linearized spectral solutions of DEC-QED then serves as a basis for expanding this non-linear dynamics. In the current work, we do not discuss the non-linear aspects of evolution, which is discussed at length in Ref.  [10].

The main results presented in this paper are organized as follows: in Section II, we 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. Then in Section III we 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 in Section IV.1. In Section IV.2 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, in Section IV.3 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 in Sections V.1-V.2 and (2) vector spherical harmonic expansions in Section V.3, 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 [62].

II 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

$${\bf\nabla}\times{\bf\nabla}\times{\bf A}+\mu_{0}\epsilon\ddot{{\bf A}} =\mu_{0}({\bf J}_{s}+{\bf J}_{\text{src}})-\mu_{0}\epsilon\nabla\dot{V},$$ (2)

$$\nabla^{2}V+\frac{\partial}{\partial t}({\mathbf{\nabla}}\cdot{\mathbf{A}}) =-\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]

$${\bf\nabla}\!\times\!{\bf\nabla}\!\times\!\bm{\mathcal{A}}+\mu_{0} \epsilon\frac{\partial^{2}\bm{\mathcal{A}}}{\partial t^{2}}+\frac{\mu_{0}q^{2}}{m}\rho\bm{\mathcal{A}}-\frac{\mu_{0}\epsilon q}{2m}\frac{\partial}{\partial t}\nabla\big{|}\bm{\mathcal{A}}\big{|}^{2}$$

$$+\frac{\mu_{0}\epsilon\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. (II) that corresponds to the transverse excitations then reproduces London theory [63] and is sourced by the fluctuation in the source current. A derivation of the linearization is provided in Appendix A. 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}}}$:

$${\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.

III 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$, and given a scalar function $\rho$ we define the coarse-grained variable

$$Q(v^{\dagger})=\int_{v^{\dagger}}dV\rho,$$ (8)

where the integral is done over the dual volume $v^{\dagger}$. 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 dielectric function $n^{2}$ or $1/\lambda_{L}^{2}$ (which is proportional to the bulk condensate density) can be assigned to objects living on the dual mesh. For example, for every edge $e$ a value for $1/\lambda_{L}^{2}$ 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 value $\overline{1/\lambda_{L}^{2}}(e^{\dagger})$ there will be the weighted average of the values 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{\bigg{(}\frac{1}{\lambda_{L}^{2}}\bigg{)}}(e^{\dagger})=\sum_{c_{i}}\frac{|e^{\dagger}\cap c_{i}|}{\Delta A(e^{\dagger})}\frac{1}{\lambda_{L}^{2}(c_{i})},$$ (9)

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}$.

This DEC formulation of the electrohydrodynamics of coarse-grained quantities offers several benefits: (1) it provides a natural framework for describing gauge-invariant fields that are consistent with the standard flux-based description of Josephson dynamics. (2) Since the gauge-invariant fields $(\bm{{\mathcal{A}}},\rho)$ permeate all of space, the issues related to boundary conditions at material interfaces are simplified by introducing effective material properties. (3) Because the Josephson phase is, by definition, a gauge-invariant and coarse-grained generalized flux variable, our approach is a rigorous generalization of the coarse-grained formulation beyond the junction and apply to the entire discretized 3D(2D) computational space. It is therefore different from other fullwave approaches because here the fundamental variables of interest are the coarse-grained fields that hybridize light and matter. The quantization of this electromagnetic theory therefore follows a similar procedure as is typically done for 1D transmission line cQED, but now with fluxes living on the edges and charges living on the nodes of the fully discretized 3D(2D) system. This approach, in which all of the classical and quantum analysis is done at the EM field level, is distinct from circuit diagram-based methods such as blackbox quantization [15, 16] or energy-participation ratio [20]. We sacrifice the simplicity brought about by the reduction to the circuit picture to obtain accurate mode structure and the full transparency in the ensuing quantization.

That being said, before the theory can be quantization-ready, a rigorous spectral theory for the coarse-grained quantities that can correctly account for radiative losses and mode hybridization, both in the eigenfrequencies and in the spatial profile of the modes (including and particularly on the open boundaries), is needed and is the focus of this paper. In the following sections, we demonstrate with specific examples how this is achieved in our approach.

IV Modes of closed superconducting systems

IV.1 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)

$$\sum_{e_{0}\in\partial(e^{\dagger})} \sum_{e_{1}\in\partial(e_{0}^{\dagger})}\frac{\Delta\ell(e_{0})}{\Delta A(e_{0}^{\dagger})}\Phi(e_{1})$$ (10)

$$+\bigg{(}\overline{\frac{1}{\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{1/\lambda^{2}_{L}},\overline{n^{2}}\}$ are defined as in Eq. (9). 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,$$ (11)

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. (11)) 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. (9). 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. (10) 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.

IV.2 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 [64]. 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{)},$$ (12)

with the corresponding eigenfunction that is separable in cartesian coordinates

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

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.

IV.3 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. The ability to simulate the full multiscale system is especially important for analyzing package modes and how they affect the on-chip operations [45] - a study that requires a computational mesh capable of encompassing multiple spatial scales to cover the volume of the entire package while simultaneously resolve the fine details on the quantum chip. Another example of a multiscale system is a three-dimensional superconducting cavity containing one or a few superconducting chips. 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 [65]. 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.

V 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 (e.g. a system surrounded by a superconducting enclosure much thicker than the relevant penetration depths). Here, we would like to extend DEC-QED to allow for a mathematically rigorous consideration of radiative losses. 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. A superconducting qubit placed in an open cavity acquires a relaxation rate through hybridization of the qubit with the electromagnetic modes of the cavity. Therefore its relaxation rate is given by the imaginary part of the complex-valued eigenfrequency of the qubit-like mode, while the difference between the real part in the frequency of a hybridized mode with the internal natural frequency of the qubit gives the Lamb shift [55]. In a similar vein, radiative loss of cavity modes is given by the imaginary parts of the cavity-like modes and is modified through hybridization with the qubit. 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 approaches and present numerical demonstrations.

V.1 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 Surface of Last Scattering (SoLS), the imaginary (transparent) boundary $\partial\mathcal{D}$ (the dashed contour) containing all material systems 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 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. The definition of radiative boundary conditions is analytically stated only at spatial infinity $(r\!\rightarrow\!\infty)$. However, our goal is to state a mathematically equivalent condition at a finite distance from the cavity to keep the computationally volume as small as possible. We refer to this imaginary boundary that includes all material systems within its volume as the “Surface of Last Scattering” (SoLS). We then employ the analytically known free-space frequency-domain Green’s function to ‘propagate’ back the boundary condition at infinity to the chosen SoLS.

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,$$ (14)

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

$$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}$$ (15)

where $H_{0}^{(1)}$ is the Hankel function of the first kind. At first glance, the Green’s functions in Eq. (15) and their derivatives diverging at ${\mathbf{r}}\!=\!{\mathbf{r}^{\prime}}$ might seem like a problem. This difficulty can be circumvented by casting the boundary condition in an integral form to regularize the Green’s function singularity [66, 67]. Upon applying Green’s second 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 [61]

$$\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},$$ (16)

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. (16) is given in Appendix B for completeness. In Eq. (16), the field at a point on the boundary is determined by the field and its gradient everywhere else on the same boundary. Eq. (16) 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. (16) 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. (16) 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. (16) the unit normal vector $\hat{n}$ on $\partial\mathcal{D}^{\prime}$ points into $\mathcal{D}$.

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

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

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 [68]. 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. (14) applied to a 2D 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. Note that since all modes are radiative, their eigenvalues are complex with the imaginary parts giving the corresponding relaxation rates. 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.

V.2 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

$$\bm{\mathcal{A}}({\mathbf{r}^{\prime}})=-2 \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{]}$$ (18)

$$-\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. (15), and the principal value integral is generalized to the vector case. A detailed derivation of Eq. (18) is given in Appendix B.

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. (18) in DEC by calculating the modes of Eq. (6) applied to the same 2D 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. (18), when applied to a 2D circular boundary, then becomes

$$\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{\}},$$ (19)

for the radial boundary edge, and

$$\Phi(e^{\theta}_{N_{R},n_{\theta}})=-\Phi(e^{\theta}_{N_{R},n_{\theta}-1})-\sum_{\theta} \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{)}$$ (20)

$$-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{]}$$

$$-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. (d) Spontaneous emission rate of a qubit-like mode as a function of the qubit frequency.

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. (10) 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. Finally, in Fig. 10d, we show how one can tune the internal resonance frequency of the qubit through a cavity resonance, exposing the Purcell enhancement of the qubit relaxation rate [54]. The relaxation rate due to radiative loss of a qubitlike mode, i.e. the spontaneous emission rate, is plotted as a function of the qubit frequency. We observe a finite peak as the qubit freqency is on-resonant with a cavitylike mode, signifiying an enhanced emission rate due to the hybridization of this qubit mode with the cavity field [18]. In this case, the Purcell enhancement factor with respect to its value away from resonance is not large, because we intentionally chose a very (radiatively) lossy cavity represented by placing two cylinders facing each other, a situation that corresponds to the regime of overlapping resonances (finesse $\ll 1$).

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.

V.3 Vector spherical harmonics expansion

In this section, we discuss the implementation of open boundaries using vector spherical harmonics (VSH) expansions. The VSH [69] 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

$${\mathbf{Y}_{lm}} =Y_{lm}(\theta,\varphi)\hat{r}$$ (21)

$${\mathbf{\Psi}_{lm}} =r\nabla Y_{lm}(\theta,\varphi)$$ (22)

$${\mathbf{\Phi}_{lm}} =\mathbf{r}\times\nabla Y_{lm}(\theta,\varphi)$$ (23)

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

$$\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}},$$ (24)

where coefficients $\mathcal{A}^{r}_{lm},\mathcal{A}_{lm}^{(1)}$ and $\mathcal{A}_{lm}^{(2)}$ are functions of the radial coordinate $r$. In Eq. (24) 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. (24) 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

$$\mathcal{A}_{lm}^{r}= \frac{Z_{l+1/2}(kr)}{r^{3/2}},$$ (25)

$$\mathcal{A}_{lm}^{(1)}= \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)$$

$$\hskip 72.26999pt-Z_{l+3/2}(kr)\Big{]}\bigg{\}},$$ (26)

$$\mathcal{A}_{lm}^{(2)}= \frac{\pi}{2kr}Z_{l+1/2}(kr)$$ (27)

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. (24), 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 Appendix C. Consider a spherical boundary of radius $R$. The series expansion for the radial component of the radiating field on the sphere is given by

$$\bm{\mathcal{A}}^{r}(R,\theta,\varphi) =\sum_{l,m}\tilde{a}_{lm}^{r}A_{lm}^{r}(R){\mathbf{Y}_{lm}}(R,\theta,\varphi)$$

$$=\sum_{l,m}\tilde{a}_{lm}^{r}\frac{H_{l+1/2}(kR)}{R^{-3/2}}Y_{lm}\hat{r},$$ (28)

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

$$a^{r}_{lm} \equiv\tilde{a}_{lm}^{r}H_{l+1/2}(kR)$$ (29)

$$=R^{3/2}\int d\Omega\bm{\mathcal{A}}^{r}\cdot{\mathbf{Y}^{*}_{lm}}$$

$$=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. (29) 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

$$\frac{\mathcal{A}^{r}(R,\theta,\varphi)-\mathcal{A}^{r}(R-\Delta R,\theta,\varphi)}{\Delta R}=$$ (30)

$$\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. (30) and the expression for $\tilde{a}^{r}_{lm}$ derived in Eq. (29), the boundary condition for $\mathcal{A}^{r}$ is then

$$\mathcal{A}^{r} (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}}$$

$$\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{\}}$$

$$\hskip 86.72377pt\times Y_{lm}(\theta,\varphi).$$ (31)

The expression in Eq. (V.3) 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. (10), 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

$$\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{\}},$$ (32)

where $e^{r^{\prime}}_{b}$ is the other radial edge that shares the boundary vertex with $e_{b}^{r}$. In Eq. (32), 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 Appendix C.

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. (V.3) 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 Appendix D, we provide a preamble by discussing a brief proof of the separability of the Helmholtz equation.

VI Conclusion

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 can still display the high sensitivity to small changes in the geometry and topology of the system that is warranted for instance 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 resolution imposed by a given 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. Future work will focus on adaptive meshers that adapt to the finite spatio-temporal resolution provided by the measurement apparatus interrogating a particular dynamics, allowing for the efficient and accurate simulation of only what is measurable.

(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 calculation of radiative modes of a superconducting qubit.

The spectral theory presented here is suitable for the second quantization of the electromagnetic field and the time-dependent numerical simulation of the non-linear dynamics of open superconducting systems with arbitrary complexity. We expect the ability to accurately and efficiently extract relaxation rates, as well as the quantification of mode hybridization demonstrated here to be useful in the modeling and optimization of superconducting circuits.

VII 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.

Appendix A Linearization of the EHDS equations

In this Section, we provide a derivation for the linearization of the EHDS equations which leads to the spectral problem analyzed in this article. We start with the gauge-invariant form of the EHDS equations (Eq. (II) and (5) in the main text) and imposing that charge conservation is individually respected by both the charged condensate $\rho$ and the source $\rho_{src}$, meaning that the last two terms in Eq. (5) cancel each other. Consider a time-harmonic source current ${\mathbf{J}}_{src}={\mathbf{J}_{0}}e^{i\omega t}$ that, to first order in frequency, results in fluctuations in the gauge-invariant field $\bm{\mathcal{A}}=\bm{\mathcal{A}_{0}}e^{i\omega t}$. According to Eq. (5), this in turns leads to a fluctuating condensate

$$\rho=\rho_{0}+\delta\rho_{0}e^{i\omega t}$$ (33)

where $\delta\rho_{0}\ll\rho_{0}$. While the quadratic term $\partial_{t}\nabla|\bm{\mathcal{A}}|^{2}$ does not have a first-order contribution and can be neglected, the last term on the lhs of Eq. (II), which is related to the rate of change of the quantum pressure, needs careful consideration

$$\frac{\mu_{0}\epsilon_{0}\hbar^{2}}{2mq}\frac{\partial}{\partial t}\nabla\bigg{[}\frac{\nabla^{2}(\sqrt{\rho})}{\sqrt{\rho}}\bigg{]} =\frac{\mu_{0}\epsilon_{0}\hbar^{2}}{4mq}\nabla\left[\left(\frac{1}{\rho}\nabla^{2}+\frac{|\nabla\rho|^{2}}{\rho^{3}}-\frac{\nabla^{2}\rho}{\rho^{2}}-\frac{\nabla\rho\cdot\nabla}{\rho}\right)\frac{\partial\rho}{\partial t}\right]$$

$$=\frac{\mu_{0}\epsilon_{0}\hbar^{2}}{4m^{2}}\nabla\nabla^{2}\nabla\!\cdot\!\bm{\mathcal{A}},$$ (34)

where in the last line of Eq. (A) we have used the ansatz (33) for $\rho$ with $\rho_{0}$ being a step-wise constant function, used Eq. (5) to replace the time derivative $\partial_{t}\rho$, and only kept the leading-order terms. For transverse excitations, the operation on $\bm{\mathcal{A}}$ in Eq. (A) vanishes, resulting in the inhomogeneous source-field equation

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

where $\lambda_{L}({\mathbf{r}})=\sqrt{\frac{m}{\mu_{0}q^{2}\rho_{0}({\mathbf{r}})}}$ is the London penetration depth, $n({\mathbf{r}})=\sqrt{\tilde{\epsilon}({\mathbf{r}})}$, and $k=\mu_{0}\epsilon_{0}\omega^{2}$.

Appendix B Derivation of Green’s boundary integrals

In this Section, we show the derivations for the integral forms of the boundary fields that obey the scalar and vector Helmholtz equations as were shown in Eqs. (16) and (18), respectively.

Figure 12: Sketches of how deformations are done to evaluate diverging integrals containing the Green’s function. (a) A deformation on the boundary to include the boundary point. (b) The insertion of a hollow region internal to $\Gamma_{i}$ surrounding the point of interest. (c) Deformation of the boundary to exclude the boundary point of interest.

Consider a homogeneous domain $\Gamma_{i}$ bounded by the boundary $\partial\Gamma_{i}$, the Green’s function is the solution of the impulse response of scalar Helmholtz operator and is defined by

$$\big{[}\nabla^{2}+n_{i}^{2}k^{2}\big{]}G({\mathbf{r}},{\mathbf{r}^{\prime}},k)=\delta({\mathbf{r}}-{\mathbf{r}^{\prime}}),$$ (36)

where $\delta({\mathbf{r}}\!-\!{\mathbf{r}^{\prime}})$ is a $d-$dimensional Dirac $\delta-$function. Multiplying Eq. (13) by $G({\mathbf{r}},{\mathbf{r}^{\prime}},k)$ and Eq. (36) by $\phi({\mathbf{r}})$, then applying Green’s theorem gives

$$\phi({\bf r^{\prime}})=\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{ds},$$ (37)

for ${\mathbf{r}^{\prime}}$ strictly inside the domain boundary $\partial\Gamma_{i}$. As we take the limit ${\mathbf{r}^{\prime}}\rightarrow\partial\Gamma_{i}$, both $G({\bf r,r^{\prime}},k)$ and its normal gradient diverges. However, we can show for the 3D case that Eq. (37) is still integrable. A proof for 2D can be found in Ref. 66. We first consider an infinitesimal half-spherical deformation $C$ of the boundary at ${\mathbf{r}^{\prime}}$ such that the surface now avoids and encloses the singular point inside the domain $\Gamma_{i}$ as seen Fig. (12a). Let the radius of the deformation be $\epsilon$, then at the limit $\epsilon\rightarrow 0$ the integral of second term in Eq. (37) over $C$ vanishes

$$\lim_{\epsilon\rightarrow 0} \int_{C}\!\!G({\bf r,r^{\prime}},k)\nabla\phi({\bf r})\cdot\mathbf{ds}$$

$$=\lim_{\epsilon\rightarrow 0}-\frac{\partial\phi}{\partial n}\bigg{\rvert}_{\mathbf{r}^{\prime}}\frac{1}{4\pi}\bigg{(}\frac{1}{\epsilon}+in_{i}k+...\bigg{)}2\pi\epsilon^{2}=0,$$ (38)

while the integral of the first term over $C$ gives

$$\lim_{\epsilon\rightarrow 0}\int_{C}\!\!\phi\frac{\partial G}{\partial n}ds =\lim_{\epsilon\rightarrow 0}-\phi({\mathbf{r}^{\prime}})\frac{e^{in_{i}k\epsilon}}{4\pi}\bigg{(}\frac{in_{i}k}{\epsilon}-\frac{1}{\epsilon^{2}}\bigg{)}2\pi\epsilon^{2}$$

$$=\frac{1}{2}\phi({\mathbf{r}^{\prime}}),$$ (39)

where in the maneuvers above we have used the series expansion of the Green’s function, which reads

$$G({\bf r,r^{\prime}},k)=-\frac{e^{in_{i}k|{\mathbf{r}}-{\mathbf{r}^{\prime}}|}}{4\pi|{\mathbf{r}}-{\mathbf{r}^{\prime}}|}=-\bigg{(}\frac{1}{|{\mathbf{r}}-{\mathbf{r}^{\prime}}|}+in_{i}k+...\bigg{)}.$$ (40)

From the results in (B) and (B) it follows that the integral in Eq. (37) when ${\mathbf{r}^{\prime}}\in\partial\Gamma_{i}$ evaluated over the entire boundary surface is indeed the form shown in Eq. (16).

For the vector case, we start by introducing the vector analog of Green’s identity

$$\int_{\Gamma_{i}}\big{(}{\mathbf{u}}\cdot{\mathbf{\nabla}}\!\times\!{\mathbf{\nabla}}\!\times\!{\mathbf{v}}-{\mathbf{v}}\cdot{\mathbf{\nabla}}\!\times\!{\mathbf{\nabla}}\!\times\!{\mathbf{u}}\big{)}dV=\int_{\partial\Gamma_{i}}\big{(}{\mathbf{v}}\!\times\!\nabla\!\times\!{\mathbf{u}}-{\mathbf{u}}\!\times\!\nabla\!\times\!{\mathbf{v}}\big{)}\cdot\mathbf{ds},$$ (41)

where $\mathbf{u}$ and $\mathbf{v}$ are vector fields. We now choose ${\mathbf{v}}=\bm{\mathcal{A}}(\mathbf{r})$, our field of interest, and ${\mathbf{u}}={\mathbf{G}}\equiv G(\mathbf{r},\mathbf{r^{\prime}},k)\hat{\nu}$, where $\hat{\nu}$ is a unit vector that has a randomly selected, but constant, orientation. For a point ${\mathbf{r}^{\prime}}$ outside of $\Gamma_{i}$, Eq. (41) reduces to

$$\int_{\Gamma_{i}}\!(\hat{\nu}\!\cdot\!\nabla G)(\nabla\!\cdot\!\bm{\mathcal{A}})dV =\int_{\partial\Gamma_{i}}\!\!\Big{\{}(\hat{\nu}\!\cdot\!\nabla G)(\bm{\mathcal{A}}\!\cdot\!\hat{n})+\bm{\mathcal{A}}\!\times\!(\nabla G\!\times\!\hat{\nu})-\mathbf{G}\!\times\!{\mathbf{\nabla}}\!\times\!\bm{\mathcal{A}}\Big{\}}\cdot\mathbf{ds}$$ (42)

$$=\int_{\partial\Gamma_{i}}\!\!\Big{\{}(\hat{\nu}\!\cdot\!\nabla G)(\bm{\mathcal{A}}\!\cdot\!\hat{n})+(\hat{\nu}\!\cdot\!\bm{\mathcal{A}})(\nabla G\!\cdot\!\hat{n})-(\hat{\nu}\!\cdot\!\hat{n})(\nabla G\!\cdot\!\bm{\mathcal{A}})-\mathbf{G}\!\cdot\!\big{[}({\mathbf{\nabla}}\!\times\!\bm{\mathcal{A}})\!\times\!\hat{n}\big{]}\Big{\}}ds,$$

where $\hat{n}$ is the unit normal vector on $\partial\Gamma_{i}$. Since $\hat{\nu}$ is a common factor in all the terms, it can be dropped and we arrive at

$$\int_{\Gamma_{i}}\!\nabla G(\nabla\!\cdot\!\bm{\mathcal{A}})dV=\int_{\partial\Gamma_{i}}\!\!\Big{\{}\nabla G(\bm{\mathcal{A}}\cdot\hat{n})+\nabla G\!\times\!(\bm{\mathcal{A}}\!\times\!\hat{n})-G\big{[}({\mathbf{\nabla}}\!\times\!\bm{\mathcal{A}})\!\times\!\hat{n}\big{]}\Big{\}}ds.$$ (43)

We are interested in, however, the evaluation of the field inside $\Gamma_{i}$. To proceed, consider when ${\mathbf{r}^{\prime}}$ is inside the material region enclosed by $\partial\Gamma_{i}$, but there exists an infinitesimally small spherical hole $S$ of radius $\epsilon$ centered around ${\mathbf{r}^{\prime}}$. The boundary of $\Gamma_{i}$ now consists of the original outer surface $\partial\Gamma_{i}$ and the boundary of the sphere, as seen in Fig. (12b). At the limit $\epsilon\rightarrow 0$, each term in Eq. (43) evaluated on the sphere gives

$$\lim_{\epsilon\rightarrow 0}\int_{S}G\big{[}({\mathbf{\nabla}}\!\times\!\bm{\mathcal{A}})\!\times\!\hat{n}\big{]}ds =\lim_{\epsilon\rightarrow 0}-\frac{e^{in_{j}k\epsilon}}{4\pi\epsilon}\big{[}(\nabla\!\times\!\bm{\mathcal{A}})\!\times\!\hat{n}\big{]}\bigg{\rvert}_{\mathbf{r}^{\prime}}\!4\pi\epsilon^{2}$$

$$=0,$$ (44)

also

$$\lim_{\epsilon\rightarrow 0}\int_{S}\!\!\Big{\{}\nabla G(\bm{\mathcal{A}}\cdot\hat{n})+\nabla G\!\times\!(\bm{\mathcal{A}}\!\times\!\hat{n})\Big{\}}ds =\lim_{\epsilon\rightarrow 0}\int_{S}\Big{(}\frac{1}{\epsilon}-in_{i}k\Big{)}\frac{e^{in_{i}k\epsilon}}{4\pi\epsilon}\bm{\mathcal{A}}ds$$

$$=\bm{\mathcal{A}}{(\mathbf{r}^{\prime}}).$$ (45)

Therefore, the integral form of the field at ${\mathbf{r}^{\prime}}$ inside $\Gamma_{i}$ is given by

$$\bm{\mathcal{A}}({\mathbf{r}^{\prime}})=\int_{\partial\Gamma_{i}}\!\!\bigg{\{}G\big{[}({\mathbf{\nabla}}\times\bm{\mathcal{A}})\times\hat{n}\big{]}-\nabla G(\bm{\mathcal{A}}\cdot\hat{n})-\nabla G\times({\mathbf{A}\times\hat{n}})\bigg{\}}+\int_{\Gamma_{i}}\!\nabla G({\mathbf{\nabla}}\cdot\bm{\mathcal{A}})dV,$$ (46)

where the volume term on the right hand side is an improper integral. However, this term vanishes if the field is divergence-free, in which case the field at anywhere strictly inside $\Gamma_{i}$ is determined by the field on the boundary $\partial\Gamma_{i}$.

Finally, to find the integral form of the field on the boundary, at ${\mathbf{r}^{\prime}}\in\partial\Gamma_{i}$, we consider a half-spherical deformation of $\partial\Gamma_{i}$ such that the boundary avoids ${\mathbf{r}^{\prime}}$ and leaves it outside of the enclosed domain $\Gamma_{i}$, as shown in Fig. (12c). Note that this deformation is in the opposite direction as in the scalar case discussed earlier, where ${\mathbf{r}^{\prime}}$ was included into $\Gamma_{i}$. By following a similar limiting procedure as was laid out above for ${\mathbf{r}^{\prime}}$ strictly inside $\Gamma_{i}$ but now applied to a half-sphere, we arrive at the self-consistent integral form for the boundary field shown in Eq. (18).

Appendix C Coarse-grained open boundary conditions using VSH

In this Section, we derive the boundary conditions for the angular components of the field $\bm{\mathcal{A}}({\mathbf{r}})$ satisfying the vector Helmholtz equation. We also discuss the DEC implementation of such conditions on the tangential boundary edges.

The series expansion of the $\bm{\mathcal{A}}^{(1)}$ component is

$$\bm{\mathcal{A}}^{(1)}(r,\theta,\varphi) =\sum_{l,m}\tilde{a}_{lm}^{1}\mathcal{A}_{lm}^{(1)}{\mathbf{\Psi}}_{lm}(r,\theta,\varphi)$$

$$=\tilde{a}_{lm}\frac{1}{2l(l+1)}\bigg{[}\frac{H_{l+1/2}(kr)}{r^{3/2}}+\frac{k}{\sqrt{r}}\big{(}H_{l-1/2}(kr)-H_{l+3/2}(kr)\big{)}\bigg{]}{\mathbf{\Psi}}_{lm}(r,\theta,\varphi).$$

The coefficients $a^{1}_{lm}(R)$ in the expansion of the field on the boundary surface are then given by

$$a^{1}_{lm}(R) =\tilde{a}_{lm}^{1}\mathcal{A}_{lm}^{(1)}(R)\frac{1}{2l(l+1)}\bigg{[}\frac{H_{l+1/2}(kr)}{r^{3/2}}+\frac{k}{\sqrt{r}}\big{(}H_{l-1/2}(kr)-H_{l+3/2}(kr)\big{)}\bigg{]}$$

$$=\frac{1}{l(l+1)}\int d\Omega\bm{\mathcal{A}}^{(1)}\cdot{\mathbf{\Psi}}^{*}_{lm}$$

$$=\frac{1}{l(l+1)}\sum_{f}\frac{\Delta A(f)}{R^{2}}\bm{\mathcal{A}}_{t}(f)\cdot{\mathbf{\Psi}}_{lm}^{*}(f),$$ (47)

where $\bm{\mathcal{A}}_{t}$ is the component of $\bm{\mathcal{A}}({\mathbf{r}})$ tangential to the boundary surface, and the sum is perform over all the triangular faces $f$ on the boundary. In the last line of Eq. (C) we have utilized the orthogonality between $\bm{\mathcal{A}}^{(1)}$ and ${\mathbf{\Psi}}_{lm}^{*}$ to switch to using $\bm{\mathcal{A}}_{t}$ instead of $\bm{\mathcal{A}}^{(1)}$ in the computation of $a^{1}_{lm}(R)$, since the former is more readily accessible during numerical implementation. The field $\bm{\mathcal{A}}^{(1)}$ can be decomposed into azimuthal and polar contributions, for each of which we derive the boundary condition. The radial derivative of the polar component of $\bm{\mathcal{A}}^{(1)}$ at the boundary is given by

$$\frac{\mathcal{A}_{\theta}^{(1)}(R,\theta,\varphi)-\mathcal{A}_{\theta}^{(1)}(R\!-\!\Delta R,\theta,\varphi)}{\Delta R}=\sum_{l,m}\tilde{a}_{lm}^{(1)}\frac{1}{2l(l+1)}\bigg{[}\frac{H_{l+1/2}(kr)}{r^{1/2}}+k\sqrt{r}\big{(}H_{l-1/2}(kr)-H_{l+3/2}(kr)\big{)}\bigg{]}^{\prime}\bigg{\rvert}_{R}(\nabla Y_{lm})_{\theta}$$ (48)

from which the expression for $\bm{\mathcal{A}}_{\theta}^{(1)}$ on the boundary surface can be extracted

$$\mathcal{A}_{\theta}^{(1)}(R,\theta,\varphi)=\sum_{l,m}a^{(1)}_{lm}(R)\left(R+\Delta R\left\{1+R\frac{\Big{[}\frac{H_{l+1/2}(kr)}{r^{3/2}}+\frac{k}{\sqrt{r}}\big{(}H_{l-1/2}(kr)-H_{l+3/2}(kr)\big{)}\Big{]}^{\prime}\Big{\rvert}_{R}}{\Big{[}\frac{H_{l+1/2}(kR)}{R^{3/2}}+\frac{k}{\sqrt{R}}\big{(}H_{l-1/2}(kR)-H_{l+3/2}(kR)\big{)}\Big{]}}\right\}\right)(\nabla Y_{lm})_{\theta},$$ (49)

where $(\nabla Y_{lm})_{\theta}$ is the polar component of $\nabla Y_{lm}$, and $a^{1}_{lm}(R)$ given as in Eq. (C). The expression for the azimuthal component $\mathcal{A}_{\varphi}^{(1)}(R,\theta,\varphi)$ is exactly similar, with a replacement of $\theta\rightarrow\varphi$ in the subscripts in Eq. (49). Similar to the steps for $\mathcal{A}^{(1)}$, the derivation of BC for $\mathcal{A}^{(2)}$ starts with the expansion

$$\bm{\mathcal{A}}^{(2)}=\sum_{l.m}\tilde{a}^{(2)}_{lm}\frac{\pi}{2kR}H_{l+1/2}(kR){\mathbf{\Phi}}_{l.m}(r,\theta,\varphi)$$ (50)

from which the coefficients $a^{(2)}_{lm}$ can be computed

$$a_{lm}^{(2)}(R) =\tilde{a}^{(2)}_{lm}\frac{\pi}{2kR}H_{l+1/2}(kR)$$

$$=\frac{1}{l(l+1)}\int d\Omega\bm{\mathcal{A}}^{(2)}\cdot{\mathbf{\Phi}}_{lm}^{*}$$

$$=\frac{1}{l(l+1)}\sum_{f}\frac{\Delta A(f)}{R^{2}}\bm{\mathcal{A}}_{t}(f)\cdot{\mathbf{\Phi}}_{lm}^{*}(f).$$ (51)

We also write the boundary conditions for $\mathcal{A}_{\theta}^{(2)}$ and $\mathcal{A}_{\varphi}^{(2)}$ separately using their discrete form of radial derivative at the boundary. This leads to

$$\mathcal{A}_{\theta}^{(2)}(R,\theta,\varphi)=\sum_{l,m}a_{lm}^{(2)}(R)\left\{R+\frac{\Delta R}{2}\left[1+kR\frac{H_{l-1/2}(kR)-H_{l+3/2}(kR)}{H_{l+1/2}(kR)}\right]\right\}(\hat{r}\times\nabla Y_{lm})_{\theta},$$ (52)

and by replacing $\theta$ with $\varphi$ in the subscripts, we obtain the expression for $\mathcal{A}_{\varphi}^{(2)}(R,\theta,\varphi)$.

The boundary conditions given in Eqs. (49) and (52) are for vector individual components, which are scalars living on vertices of the primal computational mesh. We would like to translate these BCs to fields living the edges that lie on the boundary surface. First of all, the practical computation of the expansion coefficients in Eqs. (C) and (C) requires knowledge about the tangential field $\bm{\mathcal{A}}_{t}$ defined on the centers of the triangular faces $f$ that discretize the boundary surface. In DEC, however, vectors are replaced by their projections onto the discrete edges, which are quantities that we have access to. Therefore, at each triangular face, given the edge fields that are on the three sides of the triangle, we need to determine the value of the vector $\bm{\mathcal{A}}_{t}$ at the circumcenter $f^{\dagger}$. This mapping from primal edge field to dual vector field is done through a sharp (#) operator, whose formal definition is given in Ref. 70. Now given that the coefficient expansions are obtained, for every vertex on the boundary, the polar(azimuthal) component of $\mathcal{A}_{t}$ is the sum of projections from $\mathcal{A}^{(1)}$ and $\mathcal{A}^{(2)}$

$$\mathcal{A}_{\theta(\varphi)}(R,\theta,\varphi)=\mathcal{A}_{\theta(\varphi)}^{(1)}(R,\theta,\varphi)+\mathcal{A}_{\theta(\varphi)}^{(2)}(R,\theta,\varphi).$$ (53)

For an edge $e[v_{1},v_{2}]$ lying on the boundary, where $v_{1}$ and $v_{2}$ are the starting and ending vertices of the edge, the edge field is given by

$$\Phi(e[v_{1},v_{2}])=\frac{\big{[}\bm{\mathcal{A}}_{t}(v_{1})+\bm{\mathcal{A}}_{t}(v_{2})\big{]}}{2}\cdot({\mathbf{v}}_{2}-{\mathbf{v}}_{1}),$$ (54)

where ${\mathbf{v}}_{1}$ and ${\mathbf{v}}_{2}$ are the locations of the vertices $v_{1}$ and $v_{2}$, respectively, with the components of $\bm{\mathcal{A}}_{t}$ on the boundary given in Eq. (53). Eq. (54) is the boundary condition to be imposed on the edges lying tangentially to the boundary surface.

Appendix D Separability of Helmholtz equation in ellipsoidal coordinates

The ellipsoidal coordinate system is based on the equations

$$\frac{x^{2}}{\xi^{2}_{i}-a^{2}}+\frac{y^{2}}{\xi^{2}-b^{2}}+\frac{z^{2}}{\xi^{2}-c^{2}}=1,$$ (55)

where $a\geq b\geq c$, with $i=1,2,3$ such that

$$\xi_{1}>a>\xi_{2}>b>\xi_{3}>c.$$ (56)

Eqs. (55) represent three families of confocal quadric surfaces sharing the same foci. In the discussion here, to simplify the algebraic manipulations we consider the case where $c=0$. The relationship between the ellipsoidal coordinates $(\xi_{1},\xi_{2},\xi_{3})$ and the Cartesian coordinates are given by

$$x =\sqrt{\frac{(\xi_{1}^{2}-a^{2})(\xi_{2}^{2}-a^{2})(\xi_{3}^{2}-a^{2})}{a^{2}(a^{2}-b^{2})}},$$

$$y =\sqrt{\frac{(\xi_{1}^{2}-b^{2})(\xi_{2}^{2}-b^{2})(\xi_{3}^{2}-b^{2})}{b^{2}(b^{2}-a^{2})}},$$ (57)

$$z =\frac{\xi_{1}\xi_{2}\xi_{3}}{ab},$$

with the scale factors being

$$h_{1} =\sqrt{\frac{(\xi_{1}^{2}-\xi_{2}^{2})(\xi_{1}^{2}-\xi_{3}^{2})}{(\xi_{1}^{2}-a^{2})(\xi_{1}^{2}-b^{2})}}$$

$$h_{2} =\sqrt{\frac{(\xi_{2}^{2}-\xi_{1}^{2})(\xi_{2}^{2}-\xi_{3}^{2})}{(\xi_{2}^{2}-a^{2})(\xi_{2}^{2}-b^{2})}}.$$ (58)

It is known that the Helmholtz equation (Eq. (13) in the main text) is separable in eleven three-dimensional coordinate systems, with the ellipsoidal coordinates being the most general of them and the remaining ten - including the spherical coordinates used in the main text - are derived from it through limiting processes [71]. The separability in ellipsoidal coordinates can be shown by using the ansatz

$$\phi(\xi_{1},\xi_{2},\xi_{3})=\psi_{1}(\xi_{1})\psi_{2}(\xi_{2})\psi_{3}(\xi_{3})$$ (59)

for the field in Eq. (13), resulting in the separated ODEs of the form

$$4\sqrt{f(\kappa_{i})}\frac{d}{d\kappa_{i}}\!\left(\!\sqrt{f(\kappa_{i})}\frac{d\psi_{i}}{d\kappa_{i}}+(\lambda_{1}+\lambda_{2}\kappa_{i}+k^{2}\kappa_{i}^{2})\!\right)\psi_{i}=0,$$ (60)

where $f(\kappa_{i})=\sqrt{(\kappa_{i}-a^{2})(\kappa_{i}^{2}-b^{2})\kappa_{i}}$, with $\kappa_{i}=\xi^{2}$. Upon performing a change of variables

$$\kappa_{i}=b^{2}sn^{2}\left(\alpha,\frac{b^{2}}{a^{2}}\right),$$ (61)

where $sn(u,m)$ is a Jacobi elliptic function, we can write

$$f(\kappa_{i})=b^{4}a^{2}sn^{2}\left(\alpha,\frac{b^{2}}{a^{2}}\right)dn^{2}\left(\alpha,\frac{b^{2}}{a^{2}}\right)cn^{2}\left(\alpha,\frac{b^{2}}{a^{2}}\right),$$ (62)

with $dn(u,m)$ and $cn(u,m)$ also Jacobi elliptic functions. Eq. (60) can then be written in terms of the new variable as

$$\frac{d^{2}\psi_{i}}{d\alpha^{2}}+\left[\frac{\lambda_{1}}{a^{2}}+\frac{\lambda_{2}}{a^{2}}sn^{2}\left(\alpha,\frac{b^{2}}{a^{2}}\right)+\frac{k^{2}}{a^{2}}sn^{4}\left(\alpha,\frac{b^{2}}{a^{2}}\right)\right]\psi_{i}=0.$$ (63)

Eq. (63) is the ellipsoidal wave equation whose solutions are known as ellipsoidal wave functions. With the Helmholtz equation being separable in ellipsoidal coordinates with known solutions, the generalization to this coordinate system of our spherical harmonics method for open systems is therefore straightforward.