{'fraction_non_alphanumeric': 0.045283984308374554, 'fraction_numerical': 0.02012621524816647, 'mean_word_length': 4.2352678571428575, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonradiating and radiating modes excited by quantum emitters in open epsilon-near-zero cavities

Quantum emitters embedded in arbitrarily shaped epsilon-near-zero cavities can selectively excite both nonradiating and radiating modes.


Supplementary Note 1. Derivation of Eqs. 1 to 3.
In this section we derive equations (1)-(3) of the main text. To this end, we provide analytical expressions for the fields excited by a quantum emitter (QE) contained in a vacuum bubble immersed within an epsilon-near-zero (ENZ) background medium. First, we first introduce the general solution to Maxwell equations for the electromagnetic field excited by an arbitrary distribution of sources contained in a vacuum bubble immersed within an unbounded medium. Then, we observe the solution of a QE in ENZ media as a limiting case. We start by considering a distribution of sources ( ) immersed within a background medium characterized by relative permittivity , propagation constant and intrinsic medium impedance . In order to insulate the sources from the background medium, we assume that they are contained within a vacuum spherical bubble of radius 0 (see Fig. 1A). A − time-convention is assumed and omitted hereafter. Without loss of generality, the internal ( , ) and external fields ( , ) to the bubble can be written as a multipolar decomposition of Tesseral harmonics (56) where { } = { , , } is a multi-index defined so that the sum runs over all spherical multipoles and are the Stratton vector fields, which are defined as follows where the angular variation is defined by the Tesseral harmonics The functions ̂ ( ) are linear combinations of Schekunoff form of Bessel spherical functions, i.e., ̂( ) = √ /2 +1/2 ( ) with ( ) being any linear combination of the usual cylindrical Bessel functions of order (55,56). In our case, we select in order to describe the fields induced by the sources and cavity modes within the bubble, whereas outgoing waves outside the bubble. On the one hand, the coefficients { , }, corresponding to transversal magnetic (TM) and transverse electric (TE) modes, i.e., electric and magnetic multipoles, respectively, are defined by the source properties. On the other hand, the external and internal field coefficients are found by solving the boundary value problem imposed by the continuity of the tangential fields on the surface of the vacuum bubble. This exercise leads to the following solutions The above set of equations/coefficients represents the general solution to the problem. The fields excited by a quantum emitter immersed in an ENZ background medium can be found as a limiting case of such a solution. First, due to the symmetry of the problem the fields excited by a QE characterized by dipole moment =̂ correspond to those of the = 1, = 0, = multipole, with source coefficient 10 = 0 0 2 /(4 ) . Second, the impact of the ENZ background medium can be evaluated by taking the limits η → ∞ , k → 0 . In this manner, the external fields can be asymptotically written as follows =̂ 2cos +̂ sin 4 0 By comparing these fields with those excited by an electrostatic dipole it is clear that the spatial distribution of the fields excited by an insulated QE immersed in a ENZ medium corresponds to those of an electrostatic dipole (even though the dipole is dynamically oscillating with time with radian frequency ω), with effective dipole moment Note that according to (17) the effective electrostatic dipole moment is resonant at ̂1 ( 0 0 ) = 0.

Supplementary Note 2. Quasi-static solution to the problem.
Here we demonstrate that the 0 0 ≪ 1 limit of the general solution is consistent with the quasistatic solution to the problem. To this end, note that the quasi-static fields internal and external to the sphere containing the dipole can be written as Solving the boundary value problem at = 0 we determine the value of the and coefficients as a function of the source field . These are given by Therefore, it is clear than in the ENZ limit ( ℎ → 0) we find → 3 . Consequently, the external dipole is effectively three times larger than the internal dipole, in agreement with the full timeharmonic analysis derived in the previous section.

Supplementary Note 3. Derivation of Eq. 4.
In this note we derive equation (4) of the main text. Note that the expression for the effective electrostatic dipole moment (17) exhibits a resonance at ̂1 ( 0 0 ) = 0, where, in the absence of dissipation losses, eff becomes arbitrarily large. Here we introduce a correction expression for a finite amount of loss, characterized by the imaginary part of the relative permittivity ′′ of the background medium. To this end, we explicitly evaluate ̂1 ( 0 0 ) = 0 in the external field coefficient (11) and take the → 0 limit. In doing so, the external field coefficient can be written as follows Consequently, the external field at resonance is given by Thus, the effective electrostatic dipole at resonance can be approximated as follows

Supplementary Note 4. Magnetic dipole resonances.
Here we identify the radii of the vacuum spherical bubbles for which the magnetic dipole mode is at resonance. Inspecting the coefficients (13)- (14) it is clear that the magnetic dipole resonance of the general bubble-unbounded media system appears at the solutions of the following characteristic equation  Figure S3 depicts the l.h.s. of the characteristic equation as a function of 0 0 . It is apparent from the figure that the magnetic dipole resonance takes place approximately at 0 0 ≈ 3.14. For a vaccum bubble operating at = 10.32 this value corresponds to a radius of 0 = 5.17 , whereas for a silicon ( = 11.7) bubble it corresponds to a radius of 0 = 1.507 .

Supplementary Note 5. QEs shifted from the origin of the coordinates.
As the position of the QE is shifted from the origin of the coordinates, the fields excited by it can be written as a series of multipole sources centered at the origin of the coordinates. Specifically, the source coefficients for a QE with dipole moment placed at the position ′ are given by (61) Next, for a dipole positioned on the x-axis an oriented along ̂ the coefficients reduce to It is thus clear that the electric and magnetic dipole excitations oscillate following ̂′ 1 ( 0 ∆ )/ ( 0 ∆ ) and ̂1 ( 0 ∆ )/( 0 ∆ ), respectively, as the QE is shifted along the x-axis. Both functions are depicted in fig. S4. It is apparent from the figure that the optimal position for the excitation of the magnetic dipole response corresponds to the displacement 0 ∆ ≈ 2.08, or, equivalently, ∆ ≈ 2.08 0 /(2 √ ). Subsequently, at 0 = 10.32 the optimal displacements equal ∆ ≈ 3.42 in vacuum (i.e., = 1) and ∆ ≈ 1.00 in silicon (i.e., = 11.7).

Supplementary Note 6. Equivalent circuit model.
Here we introduce a simple equivalent circuit model to qualitatively explain the coupling dynamics for the configuration depicted in Fig. 4, i.e., for two QE embedded within resonant bubbles (for the coupled configuration ∆ = 1.00 immersed in an "arbitrarily shaped" open epsilon-near-zero (ENZ) body. The coupling parameters between two QEs can be determined by analyzing the fields excited by a classical dipole source in a system of two coupled resonators. This configuration can be modeled with the simple circuit model shown in fig. S6. In a recent article (67), we used a similar circuit model to predict the scattering from a dimer of high permittivity particles, which also corresponds to a system with two coupled resonators. In that case, the equivalent circuit model provided a quantitative agreement with the numerical simulations. In this case, a quantitative prediction is not possible, since closed-form expressions for the impedance terms of particles embedded in an arbitrarily shaped cavity are not available. However, the simple circuit model shown in the figure below is enough to capture the main physics behind the coupling of two identical resonators, and there is an excellent qualitative agreement with the numerical simulations.
In the circuit model, 0 represents the coupling of the dipole source to the resonator (i.e., the bubble containing the QE). is an impedance term that models the response of each individual resonator. Here, we adopt a time-convention to get the usual expressions for the impedance terms. Near the resonance frequency it can be simply modeled as a series RLC circuit On the other hand, is an impedance term that models the coupling between both resonators. For dipole-dipole interactions it can be assumed as a capacitive term Next, the currents excited in the equivalent circuits can be written as The coupling parameters studied in the main text are related to the currents excited in the circuits as follows: 11 ∝ Re[ 1 ], 21 ∝ Re[ 2 ], and ∆ 21 ∝ −Im[ 2 ], where the proportionality factor is determined by the mode overlap of the dipole source with the magnetic dipole mode supported by the bubbles. The minus sign in ∆ 21 arises from the different time convention − vs for time-harmonic fields in optics and electrical engineering communities.
In order to carry out a qualitative analysis we set, for the sake of simplicity, 0 = 1 , = 1 and 0 = 1 rad/s, and we choose = 125 and = 1/( 0 2 ) in order to fix the individual bubble resonance at frequency 0 , exhibiting a similar bandwidth than that observed in the numerical simulations. Moreover, we set = 0.1 and we vary the value of to model the changes in the strength of coupling that could occur, e.g., when changing the distance between the emitters or the amount of losses. Figures S7B, S7C and S7D, represent the currents excited in the equivalent circuits for three different coupling strengths characterized by = /0.001 (weak), = /0.01 (medium) and = /0.05 (strong), respectively. For weak coupling, = /0.001, we can approximate 1 ≈ 0 / and 2 ≈ 0 / 2 so that Re[ 1 ] peaks at the single resonator frequency 0 , while Re[ 2 ] changes its sign at 0 . On the other hand, for strong coupling, = /0.05, becomes comparable to . Consequently, the resonance of 1 splits into two peaks located at 2 = 0 2 ± 1/( ) corresponding to subrradiant and superradiant states according to the sign of Re[ 2 ]. In the intermediate case /0.01 we observe a mixture of both extreme cases. There is an excellent qualitative agreement with the numerical simulations in Fig. 4.   fig. S1. Sketch and dimensions of the system studied in Fig. 2. The cavity consists of an open epsilon-near-zero volume (shown as grey background) containing a few vacuum bubbles (shown in green) and a quantum emitter (shown as red arrow). The blue curves correspond to the second order polynomials ( ) = 0 + 1 + 2 2 that fit to the specified dimensions.  Fig. 4. The cavity consists of an open epsilon-near-zero volume (shown as grey background) containing two silicon (Si) spherical bubbles (shown in green) and a quantum emitter (shown as red arrow). The cavity is also pierced by a Si rod (shown in green). The blue curves correspond to the second order polynomials ( ) = 0 + 1 + 2 2 that fit to the specified dimensions.  S7. Coupling parameters in a cubic cavity. The simulation setup is identical to the one depicted in Fig. 4A, except here the "arbitrarily-shaped" cavity has been replaced by a cubic cavity of side = 20 µ . Couping parameters in cubic cavity. Individual decay rate, Г 11 , decay rate due to coupling, Г 21 , and cooperative Lamb shift, ∆ 21 , normalized with respect to the free-space decay rate Г 0 , as a function of frequency, in the coupled ( = 1 ) configuration. The ENZ host has been modeled with a dispersive Drude model ( ) = 1 − 2 / /( + ) with = 2 × 29.08 × 10 12 rad/s, and for three different amounts of loss: = 0.001 (first row), = 0.01 (second row) and = 0.1 (third row).