Fluctuation-dissipation theorem and fundamental photon commutation relations in lossy nanostructures using quasinormal modes

We provide theory and formal insight on the Green function quantization method for absorptive and dispersive spatial-inhomogeneous media in the context of dielectric media. We show that a fundamental Green function identity, which appears, e.g., in the fundamental commutation relation of the electromagnetic fields, is also valid in the limit of non-absorbing media. We also demonstrate how the zero-point field fluctuations yields a non-vanishing surface term in configurations without absorption, when using a more formal procedure of the Green function quantization method. We then apply the presented method to a recently developed theory of photon quantization using quasinormal modes [Franke et al., Phys. Rev. Lett. 122, 213901 (2019)] for finite nanostructures embedded in a lossless background medium. We discuss the strict dielectric limit of the commutation relations of the quasinormal mode operators and present different methods to obtain them, connected to the radiative loss for non-absorptive but open resonators. We show exemplary calculations of a fully three-dimensional photonic crystal beam cavity, including the lossless limit, which supports a single quasinormal mode and discuss the limits of the commutation relation for vanishing damping (no material loss and no radiative loss).


I. INTRODUCTION
Cavity-QED phenomena in photonics nanostructures and nanolasers, such as metal nanoparticles 1-5 or semiconductor and dielectric microvavities [6][7][8][9][10] , have become an important and rising field in the research area of quantum optics and quantum plasmonics over the last few decades, since it provides a suitable platform to study, e.g., non-classical light effects 11,12 and quantum information processes 13,14 .A rigorous quantum optics theory for these systems is of great importance to describe the underlying mechanisms and applications of lightmatter interaction in these dissipative systems.Over two decades ago, a seminal phenomenological quantization approach for general absorbing and dispersive spatialinhomogeneous media 15 based on preliminary work of Huttner and Barnett 16,17 as well as Hopfield 18 was introduced.The theory has already been successfully applied to many technologically interesting quantum optical scenarios, e.g., input-output in multilayered absorbing structures 19 , active quantum emitters in the vicinity of a metal sphere 20,21 , the vacuum Casimir effect 22 , strong coupling effects in quantum plasmonics 23,24 , and non-Markovian dynamics in nonreciprocal environments 25 .
A few years after the introduction of the Green function quantization scheme, the method was further confirmed by approaches using a canonical quantization formalism in combination with a Fano-type diagonalization 26,27 .While it was confirmed that, for the case of absorptive bulk media, the theory is consistent with the fundamental axioms of quantum mechanics, e.g., the preservation of the fundamental commutation relations between the electromagnetic field operators 15,28 , it has been debated recently, whether the theory can be applied to non-absorbing media or finite-size absorptive media, e.g., dielectric nanostructures in a non-absorbing background medium 29,30 .One of the reasons for an apparent problem is that the electric field operator in the formulation in Refs.15 and 28 is proportional to the imaginary part of the dielectric permittivity which, at first sight, seems to vanish in the case of non-absorbing media, and would be inconsistent with the limiting case of quantization in free space.While it was shown explicitly for one-dimensional systems, that the Green function quantization is indeed consistent with the case of non-absorbing media 28 , there were only arguments for the general case, based on the fact that one has to include a small background absorption until the very end of the calculations 22,31 .
Recently, a fully three-dimensional quantization scheme for quasinormal modes (QNMs) was presented on the basis of the mentioned Green function quantization approach and successfully applied to metal resonators and metal-dielectric hybrid resonators coupled to a quantum emitter and embedded in a lossless background medium 32 .The QNMs are solutions to the Helmholtz equation with open boundary conditions and have complex eigenfrequencies, describing loss as an inherent mode property.The classical QNM approach has been used successfully to describe light-matter interaction on the semi-classical level [33][34][35][36][37][38] , and the introduction of a mode expansion in the Green function quantization approach has the benefit to describe general dissipative systems with few modes instead of the continua used in the seminal Green function quantization approach.At the fundamental level, the quantization of such modes is essential to study and understand the physics of few quanta light sources 39,40 and aspects of quantum fluctuations of quan-tum nonlinear optics.The implicit condition, that one leaves a small imaginary part of the permittivity in the equations until the end of the calculations, was also used in Ref. 32 and led to a contribution that describes the radiative loss of the QNMs and is finite in the case of pure dielectric structures.In this way, the final results can equally be applied to lossy or lossless resonators that can be described by few QNMs, or indeed a combination of both.
In this work, we shed fresh insights on a few of these apparent problems, and show explicitly that there are no problems with taking the lossless limit of the Green function quantization as long as the limits are performed carefully.Specifically we derive, using the limit of vanishing absorption, an alternative formalism of the Green function quantization approach, which helps to clarify the underlying physics and can be used to clearly see how to quantize open-cavity modes with and without material losses.The theory can thus be consistently used for a wide range of problems in quantum optics and quantum plasmonics, including finite nanostructures in a lossless background medium.To show the importance of our modified approach, we apply the procedure to a rigorous QNM quantization scheme for finite nanostructures.Using this quantized QNM approach, we will show, explicitly, why the commutation relations between different QNM operators are connected to two dissipation processes in general; nonradiative loss into the absorptive medium (material loss) and radiative loss into the far field (radiation loss), which renders the theory also suitable for dielectric lossy resonators.In the limit of no material loss, we obtain a physically meaningful result in terms of the normalized power flow from the QNM fields.We also recover well known results for a lossless infinite medium.
The rest of our paper is organized as follows: In Section II, we first introduce the "system" and construct a sequence of geometry configurations, that fully satisfies the outgoing boundary conditions of any open cavity system.With this model, we reformulate the quantization scheme from Refs. 15, 28, 41-43.In Section III, we then rederive a Green function identity, connected to the fundamental commutation relations of the electromagnetic fields, and explicitly show the appearance of a surface term in the zero-point vacuum fluctuations for the medium-assisted electromagnetic field operators in the case of the constructed permittivity sequences.In Section IV, we recapitulate the QNM approach and adopt the modified Green function quantization approach to the quantization scheme for QNMs from Ref. 32.In Section V, we derive the fundamental commutation relations of the QNM operators connected to the radiative dissipation.First, we apply the method from Section III directly as a special case to the QNM commutation relation.Second, we derive the far field limit of the commutation relation on the basis of two different methods (the Dyson approach and the field equivalence principle), to obtain the fields outside.Finally, we presents a practical example of a three-dimensional resonator, supporting a single QNM on a finite-loss and lossless photonic crystal beam.We discuss the numerically obtained QNM commutation relation in the limit of vanishing absorption, and show how loss impacts the various quantization factors.We then discuss the key equations in Section VI and summarize the main results of the work in Section VII.The main part of our paper is followed by five appendices, showing two more detailed derivations connected to a fundamental Green function identity, the derivation of the limit of the QNM commutation relation for vanishing dissipation, a discussion on the commutation relation of the electric and magnetic field operator using the QNM expansion, and details on the numerical calculations of the photonic crystal beam cavity.

A. System and introduction of permittivity sequences
To keep our model as general and realistic as possible, we investigate a spatial-inhomogeneous geometry described by the Kramers-Kronig permittivity (r, ω) = 1 + χ s (r, ω) (or complex dielectric constant), which can be separated into two main regions as depicted in Fig. 1.The volume V in is a spherical volume with radius R in , which shall include all scattering sources (e.g., metallic or dielectric nano resonators), described by the susceptibility χ s (r, ω), and is otherwise filled with vacuum, described by the spatial-homogeneous background permittivity B = 1.Furthermore, V in is bounded by an artificial fixed spherical surface S. The volume V out (λ) is a spherical shell with variable thickness λ filled with spatial-homogeneous background medium B = 1, surrounding V in and is bounded by S and S ∞ (λ).We also define V (λ) = V in + V out (λ) and take V → R 3 (all 3D space, and without any artificial surrounding surface in the far field) for the limit λ → ∞.
In particular, the sequence of spatial-homogeneous background permittivity functions reads

B. Green function quantization approach with permittivity sequences
In Refs.15, 28, and 42, the inclusion of an artificial noise term in the Maxwell equations was proposed to preserve the fundamental spatial commutation relation between the electromagnetic field operators in the case of absorbing background media.The associated inhomogeneous Helmholtz equation for the medium-assisted electric field operator with the modified permittivity model from above then reads where ĵN,α (r, ω) is a phenomenological introduced noise current density of the form ĵN,α (r, ω) = ω 0 π (α) and bα (r, ω), b † α (r, ω) are annihilation and creation operator of the medium-assisted electromagnetic field.Equa-tion (4) also has the formal solutions Êα (r, ω) = Êhom α (r, ω) where Êhom α (r, ω) is the solution to the homogeneous Helmholtz equation and the second term is the scattering solution.Here, G α (r, r , ω) is the relevant Green function, fulfilling As discussed in Ref. 42, in order to preserve the fundamental commutation relation between the electromagnetic field operators, the only allowed solution of the homogeneous problem is the trivial zero solution.This is true for the sequence of homogeneous solutions, i.e., Êhom α (r, ω) = 0. Now, in the limit of α → 0, the quantities associated to the original permittivity describing the physical problem of interest are preserved.
It is important to note that different orders of the two limits, namely (α → 0, λ → ∞), can yield a different results; also note the direct application of α → 0 (under the spatial integral) leads to a vanishing electric field operator in the limit of non-absorbing media in V (λ) (note the integral kernel is proportional to ( (α) I (r, ω)) 1/2 ).However, as we will show in the next section, the prior application of the limit of unbounded volume λ → ∞ leads to a physical contribution associated with the vacuum fluctuations on a finite boundary, that survives also in the limit of non-absorbing media.The first (direct) application of α → 0 corresponds to a point-wise convergence 45 , while the first application of λ → ∞ corresponds to a convergence in the sense of tempered distributions.Therefore, the exchange of the limits is non-trivial and we will later clarify how to correctly carry out these limits, and how to obtain physically meaningful results.
Note that with respect to the phenomenological quantization approach, other models have also adopted the approach of leaving a small imaginary part of the permittivity until the very end of the calculations 31,42 ; below we will present a more detailed mathematical treatment, and one which can easily be applied to quantized mode theories (using QNMs) in a rigorous and intuitive way.

III. GREEN FUNCTION IDENTITY AND FLUCTATION-DISSIPATION RELATION FOR FINITE SCATTERING OBJECTS
Important physical quantities in the context of macroscopic QED are the so-called zero-point field fluctuations, described by 0| Ê(r , ω) Ê † (r, ω)|0 , where |0 is the vacuum state, implicitly defined within the Green function framework via bα (r, ω)|0 = 0.In fact, for r = r , it is well known that the zero-point field fluctuations are directly connected to the spontaneous emission rate of a quantum emitter at position r, which interacts with its electromagnetic medium 20,21 .In the following, we will first show the consistency of the modified approach with the formulas obtained in Refs.15, 20, and 21 and then calculate a form for finite scattering objects, which gives two contributions.The first contribution is a volume integral covering the scattering regions (which vanishes in the limit of no absorption), and the second contribution is a surface integral (which remains finite also in nonabsorptive cases).
Using the formal solution of the sequence of electric field operators from Eq. ( 6) and the the bosonic nature of bα (r, ω), b † α (r, ω), we obtain Equation ( 8) relates the fluctuation of the electric field (LHS) to the dissipation in terms of the absorption I .Note also that this expression is related to other important quantities in Green function quantization, e.g., the fundamental commutation relation [ Êα (r), Bα (r )] between the electric field and the magnetic field operator 41 .We will show that, in the limit α → 0, the Green function identity is obtained, where G(r, r , ω) solves the Helmholtz equation, Eq. ( 7), with permittivity (r, ω).
In a first step, we note that I = ( − * )/(2i) and use the Helmholtz equation Eq. ( 7) in combination with the dyadic-dyadic Green second identity (see App. A), to get with and we have dropped the explicit ω notation as an argument of the functions in the following to simplify the notation.We focus on the integral over the spherical surface S ∞ (λ), and change to spherical coordinates with respect to s.Using the Dyson equation in combination with analytical properties of the Green function in spatial-homogeneous media (cf.App.B), we arrive at lim λ→∞ S∞(λ) with k α,I = (k α − k * α )/(2i), and the radius R in of V in and the I g,α (r, r ) is a geometric λ-independent factor and implicitly defined via Eq.(B11), (B12), Eq. (B13) and Eq.(11).We further note, that k 2 α = ω 2 (α) B (ω)/c 2 .Performing the limit λ → ∞ on Eq. ( 12), we see that the integral on the LHS Eq. ( 12) vanishes for any α > 0, since k α,I > 0 for ω, α > 0 by construction of χ(ω) (cf.Eq. ( 2)).Also note that we have used here the positive root of the complex number The procedure and ordering of limits is connected to the weak convergence of tempered distributions 45 .To analyze this in more detail, in the case of finite scattering objects with lossless background medium, i.e., χ s (r, ω) = 0 only over a finite spatial domain, the exponential function in Eq. ( 12) can be defined as a sequence of functions of the form: where λ is the sequence index.Here, the limits cannot be exchanged, because the sequence of functions p λ (α) does not uniformly converge to lim λ→∞ p λ (α) ≡ p(α) = 0 as a function of α on any interval α ∈ [0, r] with r > 0, because where "sup" is the supremum maximum.Importantly, this is also true for the special case of vacuum, i.e., χ s (r, ω) = 0 for all r ∈ R 3 .
In the case of a lossy spatial-homogeneous medium, i.e., χ s (r, ω) = χ s (ω) = 0 for all r ∈ R 3 , the sequence reads where χ s (ω) is complex.Here, pλ (α) does uniformly converge to 0 for any α ∈ [0, r], since the supremum maximum is still located at α = 0, but The key difference to the former case is, that there is still absorption at α = 0 (k α=0,I = 0).The same holds also true in any spatial-inhomogeneous media with a lossy background permittivity B = 1.
In any of the discussed cases, the exponential function and RHS of Eq. ( 12) vanish.Since there is no λdependent term left, we can apply the limit α → 0 on the surviving term from Eq. ( 10) to obtain the relation which is a corrected version of the relation from Ref. 15.In fact, if we have started with Êα=0 (r, ω) instead of Êα (r, ω), i.e. setting α = 0 in Eq. ( 12), the integral over S ∞ does not vanish for vacuum or finite scattering objects with lossless background medium, although G vanishes.To achieve this, we look at the vacuum case, (r, ω) = 1, and thus G = G B with k α=0 = ω/c.The exponential function on the RHS of Eq. ( 12) would immediately turn to 1 and we are left with where ŝ is the radial basis vector in spherical coordinates.Looking at the special case r = r = r 0 , we can perform the angular integrals analytically: Putting this back into Eq.( 10), with α = 0 and G = G B , gives However, we also note that Im[G B (r 0 , r 0 )] = ω 3 /(6πc 3 ) and hence M α=0 (r 0 , r 0 , ω) = 0.This is not consistent with the fluctuation-dissipation theorem, since it implies 0| Ê † (r 0 , ω) Ê(r 0 , ω )|0 = 0 (cf.Eq. ( 8)), which contradicts with the well known quantum interaction of an atom with free-space environment.In contrast, using the permittivity sequences, we get in complete agreement with the fluctuation-dissipation theorem, signified by the relation in Eq. ( 8).Thus the introduction of α > 0 is essential here.
It follows from Eq. ( 17) that the zero point fluctuations at r and r is equal to the imaginary part of the propagator between points at r, r , which is a physical appealing result with respect to the fluctuation-dissipation theorem, as discussed in, e.g., Ref. 28.However, we already showed that, in the limit of non-absorbing media, the lhs of Eq. ( 18) is zero, and hence the zero-point fluctuations vanish, which means that the formal quantization approach 15,28 cannot describe the limiting case of vacuum fluctuations.This is not the case in Eq. ( 17) anymore, since the lhs must always be regarded as a limit, where λ → ∞ must be applied before α → 0.
However, for some application/variations of the Green function quantization approach (e.g., for the QNM quantization below), we cannot directly exploit the relation in Eq. ( 17) and need to calculate a similar form of the LHS of Eq. (17), which in its current form seems to be impractical, since the integral domain is clearly the whole space.Especially in the important case for applications of a finite nanostructure within a lossless spatial-homogeneous background, it would be desirable to calculate the zeropoint fluctuations with integrals over a finite region.To circumvent the integration over a infinite region, and to show that there is a contribution connected to the radiative damping in the LHS of Eq. ( 18), we also derive a different variant of Eq. ( 17) (and Eq. ( 18)).
We start again with the rhs of Eq. ( 8), but now split V (λ) into volume integrals over V in and V out (λ).Once more we use the Helmholtz equation (7) in combination with the dyadic-dyadic Green second identity, to arrive at (similar to the derivation in App.A): with the combined surface S (λ) = S + S ∞ (λ), and V in (and its radius R in ) is chosen, such that r, r ∈ V in .The integral over S ∞ vanishes as we have already shown earlier.Since there is no λ-dependent term left, we can apply the limit α → 0 on the surviving terms from Eq. ( 23) (volume integral over V in and surface integral over S), to obtain the relation which is one of the main results of this work.
It is important to note that we can now choose V in and its surrounding surface S arbitrarily, as long as these remain finite and contain all sources I (s), as well as r and r.For a practical evaluation involving finite absorptive nanostructures, the volume integral is always only performed over the absorptive regions.In the case of no absorption, only the surface integral remains.

IV. FORMAL RESULTS OF A QUASINORMAL MODE QUANTIZATION SCHEME
A. Quasinormal mode approach Here, we recapitulate the definition and properties of QNMs, used for the quantization scheme (described in the next subsection).Similar to the construction of the electric field operator in Subsection II B, we also take into account the permittivity sequences (α) (r, ω).In this context, we define the QNM eigenfunctions fµ (r) (implicitly α-dependent) with mode index µ = {±1, ±2 . . .} as solutions to the Helmholtz equation: together with outgoing boundary conditions for the specific geometry in Fig. 1, i.e, the Silver-Müller radiation conditions 34 , which are asymptotic conditions for |r| → ∞.The radiation condition leads to complex QNM eigenvalues ωµ = ω µ − iγ µ , where ω µ and γ µ are the center frequency and width or the decay rate of the QNM resonance, respectively.We note, that although the QNM eigenfunctions and eigenvalues are implicitly α-dependent, αχ(ω µ ) is very small in the frequency interval of interest compared to (r, ωµ ), if ω 0 is chosen far away from the relevant QNM frequencies (cf.Eq. ( 2)).
The decay of the QNMs together with the Silver-Müller radiation condition leads to a divergent behaviour of the QNM eigenfunctions in far-field regions V out (λ).This makes the normalization of these modes a challenging task and it involves usually volume and surface integrals of the geometry of interest 33,35,46,47 .However, for many practical cases, e.g., for spherical structures with spatialhomogeneous background medium, it was shown 46 that, for positions inside the scattering geometry, the QNMs form a complete basis or (at least) can be used approximately to expand the full electric field.In particular, the full (transverse) Green function G α (r, r , ω) can be expanded in these basis functions for positions in the scattering geometry, using the representation for the Green dyad [33][34][35]46,48 , where fµ (r) are normalized QNMs and the coefficients A µ (ω) are given by We note that there are alternatives forms of A µ (ω) 46 , which can be converted into one another by using sum rules of the QNMs, but require additional terms in Eq. ( 27) (cf.Ref. 49).However, we emphasize that the QNM Green function from Eq. ( 27) with the above choice of A µ (ω) is consistent with the fundamental commutation relation in the quantization approach for lossy and absorptive media from Section II (cf.App.D).We further note, that while the QNM expansion approximates solely the transverse part of the total Green function, the longitudinal part can easily be obtained from the background Green function (cf.App.D), which however, is negligible in cavity-QED scenarios, where, e.g., a quantum emitter is placed near the scattering object.
To obtain the Green function also for positions outside the resonator geometry (i.e., in the background medium), one can exploit the Dyson equation.We do this, since the QNMs are not a good representation of the fields outside the scattering geometry, and they need a regularization to prevent spatial divergence 50,51 .Here, we again assume that G ff (r, r , ω) is an accurate approximation to the full Green function inside the resonator geometry region, to obtain the Green function for r, r outside the scattering geometry as 50 G α (r, r ) = G FF,α (r, r ) + G other,α (r, r ) with the QNM contribution and all other contributions associated to the background are summarized in G other,α (r, r ).The α-dependent fields Fµ (r, ω; α) are regularized QNM functions, (30)  which can be obtained from the QNMs inside the scattering geometry, or by using numerical calculations of the QNMs in real frequency space 51 .Here, B (ω) is the permittivity difference with respect to the background medium and we recall, that G B,α (r, r , ω) is the Green function of the spatialhomogeneous background medium, i.e., the solution to the Helmholtz equation Eq. ( 7) with (r, ω) = 1, together with suitable boundary conditions.We note, that for most practical cases, already a few QNMs are sufficient to accurately approximate the Green function on the range of frequencies of interest [51][52][53] .
In Ref. 54, it was shown how to construct an approximated regularized QNM via the field equivalence principle 55 , where the sources inside the scattering region are replaced by artificial sources on a virtual surface S around the sources (cf.Fig. 1).This yields the same field Fµ (r, ω; α) as Eq. ( 30) for positions outside the resonator, and within a frequency interval ∆ω, that covers the resonance of the QNM µ through: where are the sources on the boundary, and hµ (r ) = ∇ × fµ (r )/(iω µ µ 0 ) is the magnetic field of the associated QNM µ and n is the normal vector on the surface S .
Note that indeed within a small frequency interval around ω µ , where (r, ω) is approximately constant and Q µ > 10, the regularized expression from Eq. ( 30) and Eq. ( 31) are found to be practically the same in the far field (cf.Ref. 54).We highlight that the approximate expression from the near field to far field transformation can simplify the numerical calculations involving, e.g., far field propagation of QNMs, significantly compared to Eq. ( 30), since the internal volume integral is replaced by a surface integral (for numerical details, cf.Ref. 54).However, for single mode behaviour, one can also adopt a highly accurate regularization procedure 51 to obtain the field inside and outside the resonator, as we will discuss below for the numerical examples.

B. Quasinormal mode quantization scheme
Next, we generalize the QNM quantization procedure from Ref. 32 using permittivity sequences.Using the approximated form of the full Green function using QNMs, i.e., G ff,α (r, r , ω) and the regularization G FF,α (r, r , ω), we can rewrite the electric field operator from Eq. ( 6) in terms of QNMs.For positions r nearby the scattering geometry, we obtain the expression for the total electric field operator Êα (r) = ∞ 0 dω Êα (r, ω) + H.a. (from Eq. ( 6)), is a QNM operator, that depends implicitly on α.
The first term corresponds to the medium in the absorptive regions (V in ), while the second part corresponds to the contribution from the photons.
As shown in Ref. 32, applying a symmetrization transformation to the operators ãµ (ã † µ ), we obtain where a µ (a † µ ) are suitable annihilation (creation) operators acting on the symmetrized QNM Fock subspace and f s µ (r) = η ω η /ω µ S 1/2 ηµ fη (r) are the symmetrized QNM eigenfunctions.The radiative and material losses give rise to a loss induced symmetrization matrix, with matrix elements: At first glance, it seems that in the non-absorptive limit of I → 0, would yield S µη = 0 and hence all QNM operator related quantities vanish.However, as explained earlier, the exchange of the limits α → 0 and λ → ∞ is non-trivial.As we have shown in Section III, there is a general relation between fluctuation and dissipation including a boundary term, and with the help of the permittivity sequences, we will subsequently resolve this apparent issue in Section V.

V. COMMUTATION RELATION OF THE QNM OPERATORS IN THE DIELECTRIC LIMIT
Having discussed the general relation between fluctuation and dissipation, we will now apply the ideas from Section III to the QNM quantization scheme from Section IV.
Inspecting Eqs. ( 37)- (39) suggests that lim α→0 S µη,α seem to vanish in the limit of non-absorbing media, i.e., when I → 0. In the following, we show explicitly that S µη does not vanish in the limit without absorption, e.g., in the case of a dielectric cavity structure embedded in vacuum.In the non-absorptive limit, obviously Sin µη,α (ω) vanishes immediately, since it is independent of λ.This contribution reflects the absorptive part of the mode dissipation.Furthermore, we note that the expression for Sin µη,α (ω) for a single mode is similar to the classical absorptive part connected to the Poynting theorem 56 , and we will henceforth call it the non-radiative part Sin µη,α (ω) = Snrad µη,α (ω).We are left with the contribution connected to Sout µη (ω) = lim α→0 lim λ→∞ Sout µη,αλ (ω): Next we discuss and compare different approaches to obtain a non-vanishing contribution for Sout µη,αλ (ω) in the limits λ → ∞ and α → 0 applied in this order, which will yield the second term of the Poynting theorem, reflecting radiative loss.

A. Helmholtz equation and Green second identity
Here we apply directly a method similar to that described in Section III.We use the regularized modes Fµ (r, ω; α) from Eq. ( 30), and we note that the derivation is analogue to the case, where Fµ (r, ω; α) is obtained from the field equivalence principle (cf.Eq. ( 31)).We first use the definition of the regularized QNM fields Fµ (r, ω; α) from Eq. ( 30) and Eq. ( 40) to obtain with where ω is implicitely included (to simplify the notation).We remark that we can apply the limit α → 0 immediately to most parts on the rhs of Eq. ( 41), except for M B α,λ (r , r ), such that the volume integrals over V in reduces to the scattering region inside V in .Now using I = ( − * )/(2i), and noting that (α) (r) = (α) B for r ∈ V out (λ), we can exploit the Helmholtz equation of the background Green function from (Eq. ( 7) with (r, ω) = 1), together with the dyadic-dyadic Green second identity (Eq.(A8)), to arrive at with ) Equation ( 43) is a special case of Eq. ( 23), since r , r are always inside the scattering region, by construction of the regularized QNM fields Fµ (r, ω; α).Following a similar derivation as in Section III, we get in the limit λ → ∞ the remaining contribution Taking the limit λ → ∞ in Eq. ( 41), and inserting lim λ→∞ M B α,λ (r , r ) into the resulting equation yields Using the definition of C B,α (s, r , r ) from Eq. ( 44), and exploiting the definition of Fµ (r, ω; α), we obtain which does not vanish in the limit α → 0 and hence I → 0. Thus, after taking the limit α → 0, we arrive at Srad where Hµ (s) = ∇× Fµ (s)/(iωµ 0 ) is the regularized QNM magnetic field and n = −n points outwards of S. Note that for a single QNM the expression in Eq. ( 48) is indeed similar to the classical radiative output flow connected to the Poynting theorem 56 .We note that the surface S can be chosen arbitrarily as long as it is far away from the resonator region, such that G FF is an accurate approximation to the full Green function expansion outside the resonator.
In the limit of no radiative and no nonradiative loss, namely in the non-dissipative limit of γ µ = 0, in App.C we show that S µµ → 1.Thus the model fully recovers the well known result for normal mode quantization 57 .

B. Far field expression of Srad µη (ω)
Equation ( 48) is already a significant result and confirms the derivations in Ref. 32, namely that there is indeed a non-vanishing contribution in the case of nonabsorptive cavities.However, it would be desirable to find an expression which only depends on integrals and quantities in the system region.In the following, we will use the fact that the approximation involving Fµ (s) as a function for positions outside the scattering region improves with increased distance from the resonator.We will then derive an exact limit of Srad µη (ω).First, we choose S = S far as a spherical surface in the far field region, such that n s = s/|s|.Next, we recapitulate the radiation condition of the background Green function: for |s| → ∞, which is an analogue of the Silver-Müller radiation condition from Eq. ( 26) for real frequencies ω.
Since the regularized QNM obtained from the Dyson approach as well as the field equivalence principle, involve only terms proportional to integrals over G B (s, r), we use this relation and apply it on the far field surface S far (Eq.( 48)), to obtain Since S far is located in the far field |s| |r| and only the far-field contributions of G B play a significant role, we can rewrite Srad µη (ω) by using the limiting case of G B for far field positions in Eq. ( B9) and (B10) to get (51)  with the two representations I µη (ω) = I vol µη (ω) = I sur µη (ω), yielding using the Dyson approach (Eq.( 30)), and using the field equivalence principle (Eq.( 31)).The dyadic forms K (i) (ŝ, r) are given implicitly in Eq. ( B9) and (B10).We see that Eq. ( 51) has no explicit appearance of the far field surface and only involves integrals over the system regions, which refines the QNM quantization; Srad/nrad µη (ω) and all related quantities in the QNM quantum model can be calculated from the QNMs within the system of interest and all integrals are performed over the scattering region or bounded surfaces.In addition, it was recently pointed out that relations such as in Eq. ( 48) are generally not positive definite, because, e.g., of optical backflow 58 .However, the (exact) far field limit explicitly shows, that S rad µη is a positive definite form, since it constitutes a scalar product between different entries µ and η.Indeed, one has to choose the surface S in the general expression from Eq. ( 48) sufficiently far away from the resonator, where Fµ is a suitable representation of the fields.In fact, as shown in App.E, one needs to choose the integration surface in the general expression, Eq. ( 48), in the intermediate to far-field region with respect to the resonator in order to achieve numerical convergence.

C.
Numerical example for three-dimensional photonic crystal beam cavities Next, we discuss a numerical example to study the QNM quantization as a function of material loss and also in the limit I → 0, using a photonic crystal (PC) beam cavity (see Fig. 2 (a)), whose real part of the dielectric constant is similar to silicon nitride 52 .The threedimensional beam supports a single QNM with index 'c' (cavity mode) in the spectral regime of interest.
The precise details of the structure and numerical implementation details are given in App.E.Here we just mention that the QNM function fc and complex frequency ωc = ω c − iγ c are calculated using the method from Ref. 59 in combination with the normalization method shown in Ref. 51.In addition, because of our chosen numerical implementation, these modes are regularized QNMs ( f r ), computed in real frequency space, so they have the correct behaviour in the far field, but we refer to them here as just the QNM for ease of notation.In the near field regime, then f r → f , while in the far field regime, f r → F. Thus, we construct the Green function through 50,60 with ] for all spatial locations and at frequencies close to the resonance frequency ω c 51 .The permittivity PC (ω) of the PC beam is described by a single Lorentz oscillator model, where ∞ = 1.0, s = 2.04 2 = 4.1616.A relatively large Lorentz resonance energy ω t = 12 eV is used to ensure that it is far away from the PC resonance.We also use γ p = γ p0 = 0.131 eV, and can easily cover low to high loss regions with the same model by changing the loss term; below we use γ p values of 0.5γ p0 , 0.3γ p0 , 0.2γ p0 , and 0.1γ p0 and γ p → 0 to analyze the limit behaviour of S rad cc ≡ S rad , S nrad cc ≡ S nrad and S cc = S, from finite to vanishing absorption.Outside the PC beam, we consider free space with permittivity B = 1.0.The z-polarized dipole is positioned at d = 5 nm above the PC cavity (Fig. 2 (a)).
The classical generalized Purcell factor for a emitter with dipole moment d (=dn d ), at location r 0 , can be obtained though Green function via 36,61 where n B (= √ B = 1) is the background refractive index, and the emitter is assumed to be outside the beam which accounts for the extra factor of 1 50 .
It is also useful to compare with full dipole solutions to Maxwell's equations (i.e., with no approximations) to  56)) from single QNM, quantum Purcell factor F quan P (Eq.( 64)) from quantized QNMs, and numerical Purcell factor F num P (Eq.( 57)) from full dipole method for γp = 0. check the accuracy of the QNM expansion and to confirm that we only need one QNM.The numerical Purcell factor is defined as follows where S dipole is a small spherical surface (with radius 1 nm) surrounding the dipole point and n is a unit vector normal to S dipole , pointing outward.The vector S(r, ω) is the Poynting vector at this small surface and the subscript terms 'total' and 'background' represent the case with and without resonator.
As shown in Fig. 2 (b), there is an excellent agreement between the classical Purcell factor F QNM P (Eq.( 56)) using single QNM theory and the full-dipole numerical Purcell factor F num P (Eq.( 57)); here we use γ p = 0.1γ p0 , but more comparisons are shown explicitly in App.E for different loss values.In addition, the numerical radiative and nonradiative beta factors are defined as where the surface S PML is the interface just before the PML, and the vector S PML,total (r PML , ω) is the Poynting vector at this interface.Note that the classical beta factors are frequency dependent in general, but (in terms of the mode of interest) are most important to define near ω c .In a single-QNM case, the form of the commutation relation from Eq. (37) (in the limit α → 0) together Eq. ( 38) and Eq. ( 48)/(50) simplify as with where we have chosen the surface representation I sur c (ω) from Eq. ( 53) with µ = η = c.Performing a pole approximation at ω = ω c , leads then to The calculation of S and its approximated forms involves a volume integral over the absorptive resonator region in S nrad and a surface integral in S rad , which are encoded in I spat in/out and I spat out .See App.E for further details.These quantum-derived S factors are unitless quantities, and for well isolated single QNMs in metal resonators, we have numerically found that 54,62 S = S rad + S nrad ≈1, also in this work.This is in line with the derivation of the non-dissipative limit of S in App.C, especially because we are investigating a case with high Q factor.
After obtaining the S parameters and assuming the bad cavity limit, one can obtain the quantum SE rate 62 : where ∆ ce = ω c − ω e is the frequency detuning between the emitter and the single QNM, and gc = ω c /(2 0 )d• fc (r 0 ) is the emitter-QNM coupling.Then the quantum Purcell factor is In the limit that S → 1, Eqs. ( 63)-( 64) recover the wellknown decay rate and Purcell factor from the dissipative Jaynes-Cummings model 63 , which is also in agreement with the limit of vanishing dissipation of the QNM quantum model, discussed in App. C. Furthermore, in the quantized QNM theory, the beta factors are defined from As a specific example example for γ p = 0.1γ p0 , Fig. 2  (b) shows the excellent agreement between the quantum Purcell factor F quan P (Eq.( 64)) and the numerical Purcell factor F num P (Eq.( 57)) from the full dipole method.Note that here we show the quantum Purcell factors using S rad obtained from Eq. (62b), and emphasize again, that the results using Eq.(62c) give the same value within the numerical precision.The corresponding QNM fields for γ p = 0.1γ p0 in the resonator region are shown in Fig. 2 (c) and (d).
Table I gives a complete summary of the various S factors and beta factors as a function of loss including the lossless limit.Within the numerical precision, we also find that the total S = 1 ± 0.03 for the PC beams (with and without material loss), although there is a very slight trend from S = 0.979 to S = 1.009 from the case γ p = γ p0 to γ p = 0. Since for all cases Q 1, it is clear from the last subsection that S remains close or equal to 1 for the inspected absorption regimes.Interestingly, as mentioned above, we also find the same trend for low loss metal resonators with a Drude model, where Q ≈ 10 − 20 54 .However, S rad and S nrad change drastically in the inspected absorption range towards the limit γ p → 0, as can also be seen in Table I.

VI. DISCUSSION OF OUR KEY RELATIONS AND FINDINGS
In this Section, we discuss key relations obtained from the Green function quantization in combination with the introduction of permittivity sequences.In the first step, it was shown, that the integral relation holds.This is a corrected version of the results (cf.Eq. ( 18)) from Refs. 15 and 28, that holds not only for absorptive bulk media, but also for, e.g., vacuum case or nanostructures in a lossless environment.In a second step, we have derived a different variant of the relation in Eq. ( 67) for finite nanostructures in a lossless background medium: This second relation permits to calculate the zero-point vacuum fluctuation and related quantities, such as the Purcell factor of an atom in dissipative environment, with a finite spatial domain.It also distinguishes the two fundamental dissipative processes more clearly: Absorptive loss through the imaginary part of the permittivity and radiation loss through the Poynting vector-like term on a surface surrounding the nanostructure.Interestingly, by substituting Eq. ( 67) into Eq.( 68), a third relation is obtained, where it is important to note, that r, r ∈ V in and S is strictly finite.Equation (69) can indeed be derived by just integrating the Helmholtz equation of the Green function over a finite volume V in with boundary S (without the necessity of introducing permittivity sequences), rendering the modified approach consistent with the limit of α → 0.
A technically interesting case, is the situation where a two-level quantum emitter at position r a interacts (weakly) with the surrounding photon field.It can be shown in the framework of the Green function quantization, that the spontaneous emission rate of that quantum emitter is proportional to lim α→0 lim λ→∞ M α,λ (r a , r a ).Ref. 30 pinpoints for Eq.(67), that it seems to not recover the limit of vanishing absorption.However, Eq. ( 68) does indeed recover this.To make this clear again, we set (r, ω) = 1 for all space points.It follows, that Using Eq. ( 69), this leads to which is indeed consistent with the relation in Eq. ( 68) and shows, that the limit is preserved.
Table I.Summary of the calculated classical and quantum QNM parameters for the PC beam cavities: QNM resonance energy ωc, quality factor Qc, quantum S parameters (including S nrad from Eq. (62a), S rad from Eq. ( 62b), (we also show S rad from Eq. (62c)), and S = S rad + S nrad ), quantum nonradiative beta factor β nrad quan from Eq. ( 66) and numerical nonradiative beta factor β nrad num (Eq.( 59)) from full dipole method for various material losses γp (as well as various Im[ PC(ωc)]).When calculating S rad , the surface we are choosing is a cuboid, which is 500 nm and 200 nm away from the PC beam surface along x, z direction and y direction (see App. E).In all cases, we find that S is close to 1, within numerical uncertainty.Furthermore, Eq. ( 68) and the other variants of this relation are the key to obtain the radiative part of the commutation relation of QNM operators: which is in line with Ref. 32.Since the surface S can be chosen arbitrarily (as long as it remains strictly finite and is far away from the resonator), we have found an exact limit of Srad µη (ω), which does not depend on the shape of S and yields a symmetric, positive and bilinear form (73) In Eq. (73), only integrals over the resonator region remain, which improves the modal quantization approach.

VII. CONCLUSIONS
In summary, we have reformulated a Green function quantization approach by introducing a sequence of permittivities including a fictitious (lossy) Lorentz oscillator with weight parameter α.In the limit α → 0, the sequence converges to the actual permittivity, allowing us to address explicitly the situation of a lossy or lossless resonator in a infinite lossless background medium.It was shown that the approach, in this form, rigorously recovers the connection of the fluctuation and dissipation, even in the case without absorption.Furthermore, we have explicitly derived a form of the zero-point field fluctuations that includes a volume integral over the absorptive region and a surface term, connected to the radiative dissipation into the lossless background medium, which does not depend on the imaginary part of the permittivity.We have then applied the method to a recent QNM quantization scheme and confirmed a contribution associated to the radiative loss.
In addition, we have studied a concrete numerical example of a lossy photonic-crystal cavity, supporting a single QNM.The model includes the dispersion and loss using a rigorous solution to the full 3D Maxwell equations, thus combining practical classical and quantum calculations on an equal footing.We then demonstrated how the radiative/non-radiative part of the QNM commutation relation is identical with the radiative/non-radiative β factor, obtained via the full Maxwell equation solution, within numerical precision.Calculations for various amounts of material loss and dispersion are presented (using a causual Lorentz model), including the practical example of a completely lossless material.These results shed light on recent controversies that have been raised in the literature about Green function quantization to finite size resonators, and on their own further demonstrate the power of using a QNM approach for rigorous quantization with open cavity modes.

+ ik
where R is the absolute value of the distance vector and R × 1 is a dyad, given as the cross product of R with every column of the unit matrix 1.In general the prefactor of G B,α (s, r) can be split into three parts scaling with R −3 , R −2 and R −1 , corresponding to near-field, intermediate field and far-field contributions, respectively.Since λ → ∞ is equivalent to |s| → ∞ in Eq. ( 10), we only take the far field contribution into account and expand where ŝ = s/|s| is a unit vector in the direction of s.In the denominator of g B α (R) (Eq.(B3)), only the first term in the Taylor expansion |s| plays a significant role.However, in the fast oscillating exponential function, also the second term −ŝ • r must be taken into account.The third term |r| 2 /(2|s|) can be neglected in both contributions, as long as holds for any finite r, which is the case here.This leads to e ikα|s| e −ikα(ŝ•r) .(B8) Furthermore, we note that R = s/R − r/R in the limit of |s| → ∞ yields the limit R → ŝ.The total far field contribution can therefore be written as and after the applying the curl as 2) α (ŝ, r), (B12) with L (i) α (ŝ, r) =K (i) α (ŝ, r) for i = 1, 2. Next, we insert Eq. (B11) and (B12) into Eq.(11).and the resulting relation then into Eq.(10).
Choosing spherical coordinates with respect to s finally leads to Eq. (12).
Appendix C: Non-dissipative limit of Sµη In the following, we look at the limit of S µη from Eq. (37) without damping in the system.We analyze the non-radiative and radiative contribution of S µη separately; the part describing ohmic loss (cf.Eq. ( 38 .

(C4)
We now inspect the case µ = η.Here the denominator i(ω µ − ω * η ) with ωµ = ω µ −iγ µ remains finite when taking the limit γ µ → 0, while the surface integral in the numerator vanishes, since the surface sources in the definition of Fµ (s, ω) from Eq. (31) vanish, and therefore S rad µη → 0 for µ = η.What remains to show is that S µµ → 1 for γ µ → 0. The diagonal elements read explicitly   56)) from single QNM, quantum Purcell factor F quan P (Eq.( 64)) from quantized QNMs, and numerical Purcell factor F num P (Eq.( 57)) from full dipole method for a z-polarized point dipole placed at 5 nm above the PC cavity with permitivity of (a) Re[ PC(1γp0)], (b) PC(0.2γp0),(c) PC(0.3γp0),(d) PC(0.5γp0) and (e) PC(1γp0). .Convergence behavior of S rad for the lossless PC beam case, as a function of the surface position.These surfaces are cuboids and we use hx, hy, and hz to represent the distance between these cuboid surfaces and PC beam surface along x, y and z direction.We define h ≡ hx ≡ hz.When h is not larger than 200 nm, hy = hx = hz.After that, we use fixed value hy = 200 nm.Within numerical precision, we found that the radiative S factor for the single QNM tends towards unity, as one might expect for a high Q resonator.
integration over both ϑ and ϕ.We found that S rad from Eq. (62c) converges when the angle step size is less than π/35, and the computed value is very close to those from Eq. (62b).This clearly indicates that the regularized field by this approach (normalization method in real frequency space 51 ) is indeed practically the same thing (as F), at least at the pole.
Finally, we subsequently obtain S = S rad + S nrad , and assuming the bad cavity limit, one can calculate the quantum Purcell factors from Eq. ( 64), which show excellent agreement with results from full dipole method (Figs.3(a)-(e)).Note that here the quantum Purcell factors are using S rad obtained from Eq. (62b), and we emphasize again, that the results using Eq.(62c) give the same value within the numerical precision.

Figure 1 .
Figure 1.Geometry of the complete system of interest.The scattering sources are contained in a spherical volume Vin with radius Rin and boundary S. Vin is surrounded by a spherical shell Vout(λ) of thickness λ, filled with a spatial-homogeneous medium and bounded by the surfaces S and S inf (λ).

Figure 2 .
Figure 2. (a) Schematic diagram of PC beam with permittivity PC in free space ( B = 1.0).The z-polarized dipole is put at d = 5 nm above the PC cavity.The width and the height of the PC beam is WPC = 376 nm and hPC = 200 nm.(b) Classical Purcell factor F QNM P (Eq.(56)) from single QNM, quantum Purcell factor F quan 1γp0.Regularized QNM field | f r z | (z-component) (c) at PC center surface and (d) at a surface 5 nm above PC surface for γp = 0.1γp0.The QNM calculations with different losses are not shown, as they look similar.

Figure 3 .
Figure 3. Classical Purcell factor F QNM P

Figure 4
Figure 4. Convergence behavior of S rad for the lossless PC beam case, as a function of the surface position.These surfaces are cuboids and we use hx, hy, and hz to represent the distance between these cuboid surfaces and PC beam surface along x, y and z direction.We define h ≡ hx ≡ hz.When h is not larger than 200 nm, hy = hx = hz.After that, we use fixed value hy = 200 nm.Within numerical precision, we found that the radiative S factor for the single QNM tends towards unity, as one might expect for a high Q resonator.
The quantum nonradiative beta factors β nrad quan are also very close to full-dipole classical results β nrad num .