Modal approximation for plasmonic resonators in the time domain: the scalar case
MModal approximation for plasmonic resonators in the timedomain: the scalar case
Lorenzo Baldassari ∗ Pierre Millien † Alice L. Vanel ∗ , ‡ February 11, 2021
Abstract
We study the electromagnetic field scattered by a metallic nanoparticle with dispersive materialparameters in a resonant regime. We consider the particle placed in a homogeneous medium ina low-frequency regime. We define modes for the non-Hermitian problem as perturbations ofelectrostatic modes, and obtain a modal approximation of the scattered field in the frequencydomain. The poles of the expansion correspond to the eigenvalues of a singular boundary integraloperator and are shown to lie in a bounded region near the origin of the lower-half complex plane.Finally, we show that this modal representation gives a very good approximation of the field inthe time domain. We present numerical simulations in two dimensions to corroborate our results.
Mathematics Subject Classification (MSC2000).
Keywords. plasmonic resonance, time-domain modal expansion, subwavelength resonators, quasi-normal modes
When describing the interaction of light with a resonating particle, summing the natural resonantmodes of the system is an intuitive and attractive approach. The modes are easily computed as theyare eigenmode solutions to a source-free problem. They are intrinsic quantities of the system and giveinsights to understand the underlying physics. Once they are calculated, the response of the systemto any given excitation can be computed at a low computational cost. A bounded, lossless systemis Hermitian and admits a basis of orthonormal eigenmodes associated to real eigenvalues. But for asystem that exhibits loss (by absorption or radiation), the classical spectral theorem cannot be usedto diagonalise the non-Hermitian operator and the eigenvalues become complex [25, 41, 30].Several authors have obtained modal expansions for non-Hermitian systems [13, 19, 22, 26, 31, 34,35, 39, 45]. Their use in nanophotonics is quite recent and is studied by many research groups inthe physics community (see the review paper [25] and references therein). Nevertheless, a number oftheoretical and numerical issues arise [15]. Modes of non-Hermitian systems are not orthogonal, usingclassical inner products. In order to satisfy the outgoing boundary conditions, these generalised modeshave complex frequencies with negative imaginary parts and, if they decay exponentially in time as t → ∞ , they grow far away from the resonating systems. This is known in the literature as Lamb’sexponential catastrophe [36]. Recently, frameworks for the computation and normalisation of thesegeneralised modes have been established in different settings [18, 23, 38, 39, 43, 37]. ∗ Department of Mathematics, ETH Zürich, Rämistrasse 101, CH-8092 Zürich, Switzerland. † Institut Langevin, 1 Rue Jussieu, 75005 Paris, France. ‡ Corresponding author: [email protected] a r X i v : . [ m a t h - ph ] F e b .2 Scope of the paper In this paper we consider the scattering of a scalar wave by an obstacle with dispersive parameters(described by a
Drude-Lorentz model ). This is a good model for the scattering of light by a dispersiveobstacle in the transverse magnetic polarisation (see [29, remark 2.1]). We work in a low-frequencyregime corresponding to relevant physical applications, such as the scattering of light in the visi-ble/infrared domain by a metallic nanoparticle whose characteristic size is a few tens of nanometers.The goal of this paper is to obtain an approximation of the low-frequency part of the scattered fieldby a dispersive obstacle in the time domain as a finite sum of modes oscillating at complex frequencies.The tools used are singular boundary integral equations and elementary functional analysis. In thispaper we do not deal with the high frequency part of the field that is usually studied with micro-localanalysis tools.
It has been shown in [1, 5, 7] that using boundary integral representation and layer potential analysis,one can define the resonant frequencies as solutions of a non-linear eigenvalue problem on the boundaryof the particle. In a low-frequency regime, i.e. at frequencies corresponding to wavelengths that areorders of magnitude larger than the particle’s size, asymptotic analysis techniques, as in [5], yield ahierarchy of boundary integral equations. The asymptotic small parameter is δωc − , where δ is the sizeof the particle, ω the frequency and c the velocity. At leading order the well-known Neumann-Poincaréoperator appears [44]. Using the Plemelj symmetrisation principle and the spectral theory of compactself-adjoint operators, the latter can be diagonalised in the appropriate functional spaces [21, 32],which allows the scattered field to be decomposed in a basis of orthogonal modes in the static case[9]. The properties of the eigenvalues of the Neumann-Poincaré operator have been extensively studiedin the literature, see the review paper [12] and references therein. For a smooth enough boundary,say C ,α for some α > , the operator is compact and its eigenvalues are real numbers convergingto zero. The eigenvalues of the Neumann-Poincaré operator in the two- and three-dimensional casesare intrinsically different. In two dimensions, the spectrum is symmetric with respect to the origin(except for the eigenvalue 1/2), so there are as many positive eigenvalues as negative. The decay rateof the eigenvalues depends strongly on the regularity of the boundary. For an analytic boundary, theeigenvalues have an exponential decay rate [10]. In three dimensions, very few surfaces are known tohave negative eigenvalues [20]. For a strictly convex C ∞ domain, there are infinitely many positiveeigenvalues and a finite number of negatives ones [11]. The eigenvalues rate of decay is much slowerthan in two dimensions: λ j = O ( j − / ) as j → ∞ [28] and zero is not in the essential spectrum [11]. We begin by describing the problem geometry and we formulate the governing equations in section 2.We introduce the layer potential and boundary integral formulation and recall the modal decomposi-tion of the static ( ω = 0 ) solution. In section 3, we prove that in three dimensions, for a strictly convexparticle, the modal expansion can be truncated due to the super-polynomial decay of the expansion’scoefficients. With a perturbation argument, we deduce from the static ( ω = 0 ) result a modal approxi-mation in the dynamic case (for a small non-zero frequency). The perturbation analysis yields size andfrequency dependent dynamic complex resonant frequencies. We show that all the resonant frequen-cies have a negative imaginary part and lie in a bounded region near the origin. Finally, in section 5,using only elementary complex analysis techniques, we give an approximation for the low-frequencypart of the scattered field in the time domain as a finite sum of modes oscillating at complex resonantfrequencies. We also show with a simple causality argument that the exponential catastrophe is notproblematic in practice. In section 6 we implement this expansion in the two-dimensional setting andillustrate the validity of our approach with numerical simulations.2 Problem geometry and formulation
We are interested in the scattering problem of an incident wave illuminating a plasmonic nanoparticlein R d , d = 2 , . The homogeneous medium is characterised by electric permittivity ε m and magneticpermeability µ m . Let D be a smooth bounded domain in R d , of class C ∞ , characterised by electricpermittivity ε c . We assume the particle to be non-magnetic, i.e., µ c = µ m . Let D = z + δB where B isthe reference domain and contains the origin, and D is located at z ∈ R d and has a characteristic size δ (cid:28) . We define the wavenumbers k c = ω √ ε c µ c and k m = ω √ ε m µ m . Let ε = ε c χ ( D ) + ε m χ ( R d \ ¯ D ) ,where χ denotes the characteristic function. We denote by c the speed of light in vacuum c = 1 / √ ε µ and by c the speed of light in the medium c = 1 / √ ε m µ m .Hereafter we use the Drude model [33] to express the electric permittivity of the particle: ε c ( ω ) = ε (cid:32) − ω p ω + iω T − (cid:33) , (1)where the positive constants ω p and T − are the plasma frequency and the collision frequency ordamping factor, respectively. Condition 1.
In two dimensions, we assume the domain D to be an algebraic domain of class Q , i.e.a quadrature domain. An algebraic domain is a domain enclosed by a real algebraic curve, namely thezero level set of a bivariate polynomial. A quadrature domain is the conformal image of the unit discby a rational function. Remark 2.1.
Algebraic domains are dense among all planar domains, so every smooth curve can bedescribed as a sequence of algebraic curves [6].
Condition 2.
In three dimensions, we assume the domain D to be strictly convex: for any two pointsin D , the line segment joining them is contained in D \ ∂D . Throughout the rest of the paper, D is assumed to satisfy conditions 1 or 2. Given an incident wave u in solution to the Helmholtz equation, the scattering problem in the frequencydomain can be modelled by ∇ · ε ( x ) ∇ u ( x ) + ω µ m u ( x ) = 0 , x ∈ R d , (2)subject to the Sommerfeld radiation condition (cid:12)(cid:12)(cid:12)(cid:12) ∂ ( u − u in ) ∂ | x | − ik m ( u − u in ) (cid:12)(cid:12)(cid:12)(cid:12) = O (cid:16) | x | − ( d +1) / (cid:17) , as | x | → ∞ , uniformly in x/ | x | , for (cid:60) k m > . The transmission conditions are given by u ( x ) | + = u ( x ) | − , x ∈ ∂D, ε m ∂u ( x ) ∂ν (cid:12)(cid:12)(cid:12)(cid:12) + = 1 ε c ∂u ( x ) ∂ν (cid:12)(cid:12)(cid:12)(cid:12) − , x ∈ ∂D. Here, ∂ · /∂ν denotes the normal derivative on ∂D , and the + and − subscripts indicate the limitsfrom outside and inside D , respectively. 3 efinition 2.1. We denote the contrast λ by λ ( ω ) = ε m + ε c ε m − ε c ) . Definition 2.2 (Resonant frequency, mode) . We say ω is a resonant frequency if there is a non-trivialsolution to equation (2) with u in = 0 . We call the solution a mode. A subwavelength resonance occurswhen a resonant frequency ω satisfies ωδc − < . Let H / ( ∂D ) be the usual Sobolev space and let H − / ( ∂D ) be its dual space with respect to theduality pairing (cid:104)· , ·(cid:105) , − . The field u can be represented using the single layer potentials S k c D and S k m D ,introduced in definition A.2, as follows: u ( x ) = (cid:40) S k c D [Φ]( x ) , x ∈ D,u in ( x ) + S k m D [Ψ]( x ) , x ∈ R d \ D, (3)where the pair (Φ , Ψ) ∈ H − ( ∂D ) × H − ( ∂D ) is the unique solution to S k m D [Ψ]( x ) − S k c D [Φ]( x ) = F , x ∈ ∂D, ε m (cid:18) I + K k m , ∗ D (cid:19) [Ψ]( x ) + 1 ε c (cid:18) I − K k c , ∗ D (cid:19) [Φ]( x ) = F , x ∈ ∂D, (4)and F = − u in ( x ) , F = − ε m ∂u in ( x ) ∂ν , x ∈ ∂D, where K k m , ∗ D is the Neumann-Poincaré operator defined in definition A.2. The trace relations for thesingle layer potential are given in lemma A.2. The goal of this section is to establish an equivalent formulation for (4) in the form A ωδc [Ψ] = F (proposition 2.1), in order to write an asymptotic expansion of the operator A ωδc (lemma 2.1) anda spectral decomposition for the limiting operator A (proposition 2.2). The scaling is new in thiscontext, but the asymptotic expansion and the spectral decomposition were first obtained in [5]. Werecall them here for the sake of completeness. The proofs are quite lengthy and technical, so they areincluded in the appendix.Recall that z is the centre of the resonator and δ its radius. We introduce the scaling x = z + δX .For each function Ξ defined on ∂D , we define a corresponding function on ∂B by (cid:101) Ξ( X ) := Ξ( z + δX ) , X ∈ ∂B . The scaling properties of the integral operators are given in appendix B. The solution (cid:101) u becomes (cid:101) u ( X ) = (cid:40) δ S k c δB [ (cid:101) Φ]( X ) , X ∈ B,u in ( z + δX ) + δ S k m δB [ (cid:101) Ψ]( X ) , X ∈ R d \ B, (5)where the single-layer potential S kδB and Neumann-Poincaré operator K kδ, ∗ B are defined by the funda-mental solution Γ kδ . The density pair ( (cid:101) Φ , (cid:101) Ψ) ∈ H − ( ∂B ) × H − ( ∂B ) is the unique solution to S k m δB [ (cid:101) Ψ]( X ) − S k c δB [ (cid:101) Φ]( X ) = 1 δ (cid:101) F , X ∈ ∂B, ε m (cid:18) I + K k m δ, ∗ B (cid:19) [ (cid:101) Ψ]( X ) + 1 ε c (cid:18) I − K k c δ, ∗ B (cid:19) [ (cid:101) Φ]( X ) = (cid:101) F , X ∈ ∂B, (cid:101) F = − u in ( z + δX ) , (cid:101) F = − δε m ∂u in ( z + δX ) ∂ν X , X ∈ ∂B. Since S k c δB : H − / ( ∂B ) → H / ( ∂B ) is invertible for k c δ small enough (see lemmas A.3 and A.6), thefollowing proposition holds. Proposition 2.1.
For d = 2 , , the following equation holds for (cid:101) Ψ : A ωδ/cB [ (cid:101) Ψ] = (cid:101)
F , (6) where A ωδ/cB = 1 ε m (cid:18) I + K k m δ, ∗ B (cid:19) + 1 ε c (cid:18) I − K k c δ, ∗ B (cid:19) (cid:16) S k c δB (cid:17) − S k m δB , (cid:101) F = (cid:101) F + 1 δε c (cid:18) I − K k c δ, ∗ B (cid:19) (cid:16) S k c δB (cid:17) − [ (cid:101) F ] . (7) Lemma 2.1 (small-volume expansion) . As ωδc − → , A ωδ/cB admits the following asymptotic expan-sion: A ωδ/cB = A B + (cid:0) ωδc − (cid:1) log (cid:0) ωδc − (cid:1) A B, + O (cid:16)(cid:0) ωδc − (cid:1) (cid:17) , d = 2 , A B + (cid:0) ωδc − (cid:1) A B, + O (cid:16)(cid:0) ωδc − (cid:1) (cid:17) , d = 3 , (8) where A B = (cid:18) ε c + 12 ε m (cid:19) I − (cid:18) ε c − ε m (cid:19) K ∗ B , (9) A B, = 1 ε m K (1) B, ( I − P H ∗ ) + (cid:18) I − K ∗ B (cid:19) (cid:101) S − B S (1) B, (cid:18) ε c I − ε m P H ∗ (cid:19) , and A B, = ε m − ε c ε m ε c (cid:18) I − K ∗ B (cid:19) S − B S B, , where the operators P H ∗ , (cid:101) S B , S (1) B, , S B, and K (1) B, are defined in appendix C.1.Proof. See appendix C.2.The operator A ωδ/cB is not self-adjoint in L so it can not be diagonalised directly to solve (6).However, in the static regime, the operator A B can be expressed simply with K ∗ B , which can besymmetrised in the Hilbert space H ∗ ( ∂B ) (see appendix A.2). Lemma 2.2 (spectral decomposition of K ∗ B ) . K ∗ B is self-adjoint with respect to the inner product (cid:104)· , ·(cid:105) H ∗ ( ∂B ) . Moreover, it is compact, so its spectrum is discrete. The spectral theorem yields thedecomposition K ∗ B = + ∞ (cid:88) j =0 λ j (cid:68) · , (cid:101) φ j (cid:69) H ∗ ( ∂B ) (cid:101) φ j , where { λ j } j ∈ N are the eigenvalues of K ∗ B and { (cid:101) φ j } j ∈ N their associated normalised eigenvectors. Proposition 2.2 (spectral decomposition of A B ) . The operator A B has the spectral decomposition A B = + ∞ (cid:88) j =0 τ j (cid:68) · , (cid:101) φ j (cid:69) H ∗ ( ∂B ) (cid:101) φ j , where ( λ j , (cid:101) φ j ) j ∈ N are the eigenvalues and normalised eigenfunctions of K ∗ B in H ∗ ( ∂B ) and τ j = (cid:18) ε c − ε m (cid:19) ( λ ( ω ) − λ j ) . roof. Direct consequence of lemma 2.2 and (9).
Corollary 2.1.
The spectral approximation of the static ( ω = 0 ) solution is given by (cid:101) u ( X ) − (cid:101) u in ( X ) = ∞ (cid:88) j =0 τ j (cid:68) (cid:101) F , (cid:101) φ j (cid:69) H ∗ ( ∂B ) δ S B [ (cid:101) φ j ]( X ) , X ∈ R d \ B, where (cid:101) F is defined in proposition 2.1. In this section we want to apply perturbation theory tools to express the solutions of (6) in terms of theeigenvectors of K ∗ B that appear in the spectral decomposition of the limiting problem in proposition2.2, and to replace τ j by a perturbed value τ j ( ωδc − ) . Classical perturbation theory will give us aTaylor expansion for τ j ( ωδc − ) in ωδc − for any j ∈ N but the remainders and validity range of theseexpansions will depend on the index j of the considered eigenvalue. In order to get a meaningfulexpansion of the scattered field we need to work with a finite number of modes. In practice, there is no need to consider the whole spectral decomposition of the field. It has beenempirically reported that only a few modes actually contribute to the scattered field. The numberof modes to consider increases as the source gets closer to the particle. In this section we give amathematical explanation of this phenomenon : the modes (cid:101) φ j are eigenmodes of a pseudo-differentialoperator of order − , and are oscillating functions. As in classical Fourier analysis, the decay with j ofthe coefficients (cid:104) (cid:101) F , (cid:101) φ j (cid:105) H ∗ ( ∂B ) will be determined by the regularity of the function (cid:101) F and the numberof modes to consider will depend on the spatial variations of (cid:101) F over ∂B . In an homogeneous mediumthe incoming field is smooth and therefore we can expect a fast decay of the coefficients. For B , a strictly convex domain in R with C ∞ -smooth boundary, and (cid:101) F ∈ H J ( ∂B ) for some J ∈ N ∗ we have : (cid:68) (cid:101) F , (cid:101) φ j (cid:69) H ∗ ( ∂B ) = o ( j − J/ ) as j → + ∞ . (10)The proof relies on a theorem from [11] which itself uses the computation of the principal symbolof the Neumann-Poincaré operator done in [27]: Theorem 3.1 (from [11], p. 7) . For B , a strictly convex domain in R with C ∞ -smooth boundary, K ∗ B has a finite number of non-positive eigenvalues. We can modify K ∗ B by adding a finite dimensionalsmoothing operator to have a positive definite elliptic pseudo-differential operator of order -1, whichwe denote by (cid:101) K ∗ B . For each real number s ∈ R there exist constants c s , C s ∈ R + such that c s || (cid:101) φ || H s − / ( ∂B ) ≤ || (cid:101) K ∗ B [ (cid:101) φ ] || H s +1 / ( ∂B ) ≤ C s || (cid:101) φ || H s − / ( ∂B ) (11) for all (cid:101) φ ∈ H s − / ( ∂B ) . Moreover there exists j ∈ N such that (cid:101) K ∗ B [ (cid:101) φ j ] = K ∗ B [ (cid:101) φ j ] and λ j > for all j ≥ j . orollary 3.1. The operator K ∗ B : L ( ∂B ) −→ L ( ∂B ) defined by K ∗ B := ( −S B ) K ∗ B ( −S B ) − isself-adjoint and has the same eigenvalues as K ∗ B . Its eigenvectors are (cid:101) ψ j = ( −S B ) [ (cid:101) φ j ] . It can bemodified by adding a finite dimensional smoothing operator to have a positive definite elliptic pseudo-differential operator of order -1, which we denote by (cid:101) K ∗ B . For each real number s ∈ R there existconstants c s , C s ∈ R + such that c s || (cid:101) φ || H s − / ( ∂B ) ≤ || (cid:101) K ∗ B [ (cid:101) φ ] || H s +1 / ( ∂B ) ≤ C s || (cid:101) φ || H s − / ( ∂B ) (12) for all (cid:101) φ ∈ H s − / ( ∂B ) . Moreover there exists j ∈ N such that (cid:101) K ∗ B [ (cid:101) ψ j ] = K ∗ B [ (cid:101) ψ j ] and λ j > for all j ≥ j . Proof. K ∗ B has the same principal symbol as K ∗ B [28, p. 8].We will also need the decay estimate of the eigenvalues of K ∗ B : Theorem 3.2 (from [28]) . For B , a strictly convex domain in R with C ∞ -smooth boundary theeigenvalues of the Neumann-Poincaré operator satisfy: λ j ∼ C B j − / , with C B a constant depending only on B : C B = (cid:18) W ( ∂B ) − πχ ( ∂B )128 π (cid:19) , where W ( ∂B ) and χ ( ∂B ) denote, respectively, the Willmore energy and the Euler characteristic of theboundary surface ∂B .Proof of proposition 3.1. Consider (cid:101) F ∈ H J ( ∂B ) . Since (cid:101) K ∗ B is a positive definite elliptic self-adjointpseudo-differential operator of order − we can write [17, p. 290]: H s ( ∂B ) = (cid:101) K ∗ B (cid:0) H s − ( ∂B ) (cid:1) ⊕ Ker (cid:16) (cid:101) K ∗ B (cid:17) , where Ker ( (cid:101) K B ) denotes the kernel of (cid:101) K ∗ B . The symbol ⊕ is to be understood in the L scalar productsense. Hence for j ≥ j : (cid:68) (cid:101) F , (cid:101) φ j (cid:69) H ∗ ( ∂B ) = − (cid:68) (cid:101) F , S B [ (cid:101) φ j ] (cid:69) L ( ∂B ) = − (cid:68) (cid:101) F , ( −S B ) [ (cid:101) ψ j ] (cid:69) L ( ∂B ) = − (cid:68) ( −S B ) [ (cid:101) F ] , (cid:101) ψ j (cid:69) L ( ∂B ) . where we used the fact that ( −S B ) is self-adjoint in L ( ∂B ) . Since ( −S B ) [ (cid:101) F ] ∈ H J + ( ∂B ) we have ( −S B ) [ (cid:101) F ] = (cid:101) K ∗ B [ (cid:101) G (1) ] + (cid:101) F (1)ker with (cid:101) G (1) ∈ H J − ( ∂B ) . Then (cid:68) ( −S B ) [ (cid:101) F ] , (cid:101) ψ j (cid:69) L ( ∂B ) = (cid:68) (cid:101) K ∗ B [ (cid:101) G (1) ] + (cid:101) F (1)ker , (cid:101) ψ j (cid:69) L ( ∂B ) = λ j (cid:68) (cid:101) G (1) , (cid:101) ψ j (cid:69) L ( ∂B ) + (cid:68) (cid:101) F (1)ker , (cid:101) ψ j (cid:69) L ( ∂B ) . Since the eigenvectors of (cid:101) K ∗ B are orthogonal in L ( ∂B ) we have: (cid:68) (cid:101) F , (cid:101) φ j (cid:69) H ∗ ( ∂B ) = − λ j (cid:68) (cid:101) G (1) , (cid:101) ψ j (cid:69) L ( ∂B ) .
7e can now write (cid:101) G (1) = (cid:101) K ∗ B [ (cid:101) G (2) ] + (cid:101) F (2)ker with (cid:101) G (2) ∈ H J − ( ∂B ) and we have (cid:68) (cid:101) G (1) , (cid:101) ψ j (cid:69) L ( ∂B ) = λ j (cid:68) (cid:101) G (2) , (cid:101) ψ j (cid:69) L ( ∂B ) . Iterating this procedure J − times yields (cid:68) (cid:101) G (1) , (cid:101) ψ j (cid:69) L ( ∂B ) = λ J − j (cid:68) (cid:101) G ( J ) , (cid:101) ψ j (cid:69) L ( ∂B ) . Hence (cid:68) (cid:101)
F , (cid:101) φ j (cid:69) H ∗ ( ∂B ) = − λ Jj (cid:68) (cid:101) G ( J ) , (cid:101) ψ j (cid:69) L ( ∂B ) . (13)We need to control the L -norm of (cid:101) G ( J ) . We can rewrite the orthogonal decomposition as ( −S B ) [ (cid:101) F ] = (cid:16) (cid:101) K ∗ B (cid:17) J [ (cid:101) G ( J ) ] + (cid:101) F (1)ker . Composing by (cid:101) K ∗ B we get: (cid:101) K ∗ B ◦ ( −S B ) [ (cid:101) F ] = (cid:16) (cid:101) K ∗ B (cid:17) J +1 [ (cid:101) G ( J ) ] . Using the right-hand side of (12) with s = J + we get (cid:13)(cid:13)(cid:13)(cid:13)(cid:16) (cid:101) K ∗ B (cid:17) J +1 [ (cid:101) G ( J ) ] (cid:13)(cid:13)(cid:13)(cid:13) H J +1 ( ∂B ) ≤ C J + (cid:13)(cid:13)(cid:13) ( −S B ) [ (cid:101) F ] (cid:13)(cid:13)(cid:13) H J ( ∂B ) . Using J + 1 times the left hand side of (12) with s − = 0 , , . . . , J yields (cid:13)(cid:13)(cid:13) (cid:101) G ( J ) (cid:13)(cid:13)(cid:13) L ( ∂B ) ≤ (cid:32) J (cid:89) s =0 c s + (cid:33) (cid:13)(cid:13)(cid:13)(cid:13)(cid:16) (cid:101) K ∗ B (cid:17) J +1 [ (cid:101) G ( J ) ] (cid:13)(cid:13)(cid:13)(cid:13) H J +1 ( ∂B ) ≤ C J + (cid:32) J (cid:89) s =0 c s + (cid:33) (cid:13)(cid:13)(cid:13) ( −S B ) [ (cid:101) F ] (cid:13)(cid:13)(cid:13) H J ( ∂B ) . Using the Cauchy-Schwartz inequality in (13) and the fact that (cid:107) (cid:101) ψ j (cid:107) L ( ∂B ) = 1 ( S B is an isometry): (cid:12)(cid:12)(cid:12)(cid:12)(cid:68) (cid:101) F , (cid:101) φ j (cid:69) H ∗ ( ∂B ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ Cλ Jj (cid:107) (cid:101) F (cid:107) H J − ( ∂B ) , where C = C J + (cid:18)(cid:81) Js =0 1 c s + 12 (cid:19) is independent of j . Using theorem 3.2 we can see that for j largeenough since λ j ∼ C B j − / we have: (cid:12)(cid:12)(cid:12)(cid:12)(cid:68) (cid:101) F , (cid:101) φ j (cid:69) H ∗ ( ∂B ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ j − J/ C ( C B ) J (cid:107) (cid:101) F (cid:107) H J − ( ∂B ) , and since j − J/ C ( C B ) J = o (cid:0) j − J/ (cid:1) we get the result. In two dimensions, the picture is a slightly different. Indeed, zero is in the essential spectrum of K ∗ D . The eigenspace associated to zero has infinite dimension and there are infinitely many negativeeigenvalues. As a result, K ∗ D can not be modified into a positive operator by a finite dimensionaloperator. However, for a certain class of domains, it is possible to show that there is a finite numberof plasmonic resonances. For example, it was shown in [6] that an algebraic domain of class Q hasasymptotically a finite number of plasmonic resonances. The asymptotic parameter is the deformationfrom the unit circle. For a larger class of domains the decay of the coefficients (cid:104) F, (cid:101) φ j (cid:105) H ∗ ( ∂D ) can bechecked numerically (see section 6). 8 .2 Modal decomposition Since the incoming wave is solution of the homogeneous Helmholtz equation in the background medium,standard elliptic regularity theory gives us u in ∈ C ∞ ( R d ) . Moreover, the particle B is assumed tobe C ∞ , so the source term in equation (6), i.e. the function (cid:101) F , is smooth on ∂B . Therefore usingproposition 3.1 we have a super-polynomial decay of the coefficients (cid:104) (cid:101) F , (cid:101) φ j (cid:105) H ∗ ( ∂B ) , and we can considerthat only a finite number of modes are excited. The number J of modes to consider depends on theincoming field. Proposition 3.2.
Assume that (cid:101) F = (cid:80) Jj =1 (cid:68) (cid:101) F , (cid:101) φ j (cid:69) H ∗ ( ∂B ) (cid:101) φ j on ∂B for some J ∈ N ∗ . The spectralapproximation of the scattered field as ωδc − → is given by (cid:101) u ( X ) − (cid:101) u in ( X ) = J (cid:88) j =0 τ j ( ω ) (cid:68) (cid:101) F , (cid:101) φ j (cid:69) H ∗ ( ∂B ) δ S k m δB [ (cid:101) φ j ]( X ) , X ∈ R d \ B, where τ j ( ω ) = τ j + (cid:0) ωδc − (cid:1) log (cid:0) ωδc − (cid:1) τ j, + O (cid:16)(cid:0) ωδc − (cid:1) (cid:17) , d = 2 ,τ j + (cid:0) ωδc − (cid:1) τ j, + O (cid:16)(cid:0) ωδc − (cid:1) (cid:17) , d = 3 , with τ j, = (cid:68) A B, (cid:101) φ j , (cid:101) φ j (cid:69) H ∗ ( ∂B ) , τ j, = (cid:68) A B, (cid:101) φ j , (cid:101) φ j (cid:69) H ∗ ( ∂B ) , and (cid:101) F is defined in proposition 2.1.Proof. Note that { (cid:101) φ j } j ∈ N forms an orthonormal basis of H ∗ ( ∂B ) . Writing (cid:16) A B + A ωδ/cB − A B (cid:17) [ (cid:101) Ψ] = (cid:101) F and using the decomposition of (cid:101) Ψ in H ∗ ( ∂B ) , (cid:101) Ψ = (cid:80) + ∞ j =0 (cid:68) (cid:101) Ψ , (cid:101) φ j (cid:69) H ∗ ( ∂B ) (cid:101) φ j , yields the following: (cid:68) (cid:101) Ψ , (cid:101) φ j (cid:69) H ∗ ( ∂B ) = τ j + (cid:68)(cid:16) A ωδ/cB − A B (cid:17) (cid:101) φ j , (cid:101) φ j (cid:69) H ∗ ( ∂B ) (cid:68) (cid:101) F , (cid:101) φ j (cid:69) H ∗ ( ∂B ) j ≤ J j > J. Using (3) and (8) concludes the proof.For each normalised eigenfunction of K ∗ B , we consider the corresponding function on ∂D , φ j ( x ) := (cid:101) φ j (cid:18) x − zδ (cid:19) . Here { φ j } j ∈ N are the rescaled non-normalised eigenfunctions of K ∗ D . Let us introduce ϕ j := φ j || φ j || H ∗ ( ∂D ) . Since || (cid:101) φ j || H ∗ ( ∂B ) = 1 , we have (see appendix B) ϕ j = (cid:40) δ − φ j , d = 2 ,δ − / φ j , d = 3 . Going back to the original unscaled problem: 9 roposition 3.3. As ωδc − (cid:28) , the spectral decomposition of the field is as follows u ( x ) = J (cid:88) j =0 τ j ( ω ) (cid:104) F, ϕ j (cid:105) H ∗ ( ∂D ) S k m D [ ϕ j ]( x ) + u in ( x ) , x ∈ R d \ D, J (cid:88) j =0 τ j ( ω ) (cid:104) F, ϕ j (cid:105) H ∗ ( ∂D ) S k c D [ ϕ j ]( x ) , x ∈ D. (14) Proof.
The scaling lemma B.1 gives S k m δB [ (cid:101) φ j ]( X ) = δ − S k m D [ φ j ]( x ) for d = 2 , . From lemma B.2, wehave (cid:104) (cid:101) F , (cid:101) φ j (cid:105) H ∗ ( ∂B ) = δ − (cid:104) F, φ j (cid:105) H ∗ ( ∂D ) for d = 3 and (cid:104) (cid:101) F , (cid:101) φ j (cid:105) H ∗ ( ∂B ) = δ − (cid:104) F, φ j (cid:105) H ∗ ( ∂D ) for d = 2 . In this section we calculate size and frequency dependent plasmonic resonances. Let j ∈ { , .., J } .Recall that τ j ( ω ) = τ j + (cid:0) ωδc − (cid:1) log (cid:0) ωδc − (cid:1) τ j, + O (cid:16)(cid:0) ωδc − (cid:1) (cid:17) , d = 2 ,τ j + (cid:0) ωδc − (cid:1) τ j, + O (cid:16)(cid:0) ωδc − (cid:1) (cid:17) , d = 3 . Definition 4.1.
We say that ω is a static plasmonic resonance if | τ j | = 0 . Definition 4.2.
We say that ω is first-order corrected plasmonic resonance if (cid:12)(cid:12) τ j + ( ωδc − ) log ( ωδc − ) τ j, (cid:12)(cid:12) = 0 or (cid:12)(cid:12) τ j + ( ωδc − ) τ j, (cid:12)(cid:12) = 0 , with d = 2 or d = 3 , respectively. Remark 4.1.
For j = 0 , we have τ = 1 /ε m , which is of size one by assumption. We exclude j = 0 from the set of resonances. For j ≥ we have P H ∗ [ (cid:101) φ j ] = (cid:101) φ j . Let us define α j := (cid:28)(cid:18) I − K ∗ B (cid:19) (cid:101) S − B S (1) B, [ (cid:101) φ j ] , (cid:101) φ j (cid:29) H ∗ ( ∂B ) , d = 2 , (cid:28)(cid:18) I − K ∗ B (cid:19) S − B S B, [ (cid:101) φ j ] , (cid:101) φ j (cid:29) H ∗ ( ∂B ) , d = 3 . Then, we can calculate τ j ( ω ) = ε m − ε c ε m ε c (cid:16) λ ( ω ) − λ j + (cid:0) ωδc − (cid:1) log (cid:0) ωδc − (cid:1) α j (cid:17) + O (cid:16)(cid:0) ωδc − (cid:1) (cid:17) , d = 2 ,ε m − ε c ε m ε c (cid:16) λ ( ω ) − λ j + (cid:0) ωδc − (cid:1) α j (cid:17) + O (cid:16)(cid:0) ωδc − (cid:1) (cid:17) , d = 3 . Lemma 4.1.
We have α j ∈ R and α j := (cid:18) λ j − (cid:19) (cid:68) S (1) B, [ (cid:101) φ j ] , (cid:101) φ j (cid:69) − / , / , d = 2 , (cid:18) λ j − (cid:19) (cid:68) S B, [ (cid:101) φ j ] , (cid:101) φ j (cid:69) − / , / , d = 3 .
10n what follows we use the lower-case character ω for real frequencies and the upper-case character Ω for complex frequencies. Proposition 4.1.
Using the Drude model (1) , the three-dimensional first-order corrected plasmonicresonances Ω ± j ( δ ) := ± Ω (cid:48) j + i Ω (cid:48)(cid:48) j all lie in the lower part of the complex plane and their modulus isbounded. In the case where we take the medium to be vacuum, i.e., ε m = ε we obtain explicitly for | λ j + 1 / | > − (this occurs, for example, when B is a ball [2]): Ω (cid:48) j = (cid:118)(cid:117)(cid:117)(cid:117)(cid:116) ω p ( λ j + 1 / ω p δc − ) α j − T − (cid:104) ω p δc − ) α j (cid:105) and Ω (cid:48)(cid:48) j = − T − (cid:104) ω p δc − ) α j (cid:105) . Moreover, they are bounded | Ω j | ≤ T − (cid:12)(cid:12)(cid:12) ω p δc − ) α j (cid:12)(cid:12)(cid:12) , ω p (cid:112) λ j + 1 / (cid:114)(cid:12)(cid:12)(cid:12) ω p δc − ) α j (cid:12)(cid:12)(cid:12) . Proof.
We have that τ j (Ω j ) = 0 if and only if Ω j + i Ω j T − ω p − − λ j + Ω j δ c α j = 0 , that is ω p (cid:0) Ω (cid:48) j − Ω (cid:48)(cid:48) j (cid:1) − − λ j + (cid:0) Ω (cid:48) j − Ω (cid:48)(cid:48) j (cid:1) δ c α j − ω p T Ω (cid:48)(cid:48) j = 0 , ω p Ω (cid:48) j Ω (cid:48)(cid:48) j + 2Ω (cid:48) j Ω (cid:48)(cid:48) j δ c α j + 1 ω p T Ω (cid:48) j = 0 . Because δω p c − (cid:28) , we get the desired result. Lagrange improved upper-bound for roots of polyno-mials concludes the proof [24]. Definition 4.3.
In three dimensions, we define the resonance radius as R ( δ ) := max j ∈{ ,..,J } max − (cid:12)(cid:12)(cid:12) ω p δc − ) α j (cid:12)(cid:12)(cid:12) , ω p (cid:112) λ j + 1 / (cid:12)(cid:12)(cid:12) ω p δc − ) α j (cid:12)(cid:12)(cid:12) / . Remark 4.2.
This resonance radius gives our method a range of validity. We compute resonantfrequencies in a perturbative quasistatic regime. So by checking that R ( δ ) δc − < , we ensure that the largest plasmonic frequency lies in a region that is still considered as low-frequency for a particle of size δ . If we pick the size to be too large, namely such that R ( δ ) δc − is bigger thanone, it means that the method is not self-consistent, as the largest resonant frequency might not satisfythe ωδc − < / . Proposition 4.2.
In vacuum, and using the Drude model (1) , the two-dimensional first-order correctedplasmonic resonances are the roots (Ω j ) ≤ j ≤ J ∈ C of the following equation Ω j + i Ω j T − ω p − − λ j + (cid:0) Ω j δc − (cid:1) log (cid:0) Ω j δc − (cid:1) α j = 0 . (15)11 emark 4.3. We can compute an approximation of the roots of (15) by computing in the first placethe static resonances (Ω s,j ) ≤ j ≤ J . Solving τ j = 0 yields Ω ± s,j = ± (cid:115) ω p (cid:18) λ j + 12 (cid:19) − T − i T .
Replacing the dynamic frequency in the logarithm by its static approximation, we transform (15) intothe quadratic equation Ω j + i Ω j T − ω p − − λ j + (cid:0) Ω j δc − (cid:1) log (cid:0) Ω s,j δc − (cid:1) α j = 0 . We get Ω ± j ( δ ) = − i T − ± (cid:114) ω p ( λ j + 1 / (cid:104) α j ( ω p δc − ) log (cid:0) Ω ± s,j δc − (cid:1)(cid:105) (cid:104) α j ( ω p δc − ) log (cid:0) Ω ± s,j δc − (cid:1)(cid:105) . (16) Definition 4.4.
In two dimensions, we define the resonance radius as R ( δ ) := max j ∈{ ,..,J } max − (cid:12)(cid:12)(cid:12) α j ( ω p δc − ) log (cid:0) Ω ± s,j δc − (cid:1)(cid:12)(cid:12)(cid:12) , ω p (cid:112) λ j + 1 / (cid:12)(cid:12)(cid:12) α j ( ω p δc − ) log (cid:0) Ω ± s,j δc − (cid:1)(cid:12)(cid:12)(cid:12) / . Quasi-normal modes are formally defined as solutions of the source-free wave equation [25]. Using therepresentation formula (3), we can now define, as in the physics literature, plasmonic quasi-normalmodes ( e ± j ) j ∈ N that oscillate at complex frequencies Ω ± j ( δ ) : e ± j ( x ) = S Ω ± j ( δ ) c − D [ ϕ j ]( x ) , x ∈ R d \ D, S Ω ± j ( δ ) √ ε c µ m D [ ϕ j ]( x ) , x ∈ D. (17)These ( e ± j ) j ∈ N solve the source-free Helmholtz equation and satisfy the radiation condition at infinity,but they diverge exponentially fast as | x | → ∞ . Remark 4.4.
In the physics literature (see [25, equation (1.1)] for instance) one can often find rep-resentations of the scattered field in the form u ( x, ω ) = (cid:88) j α j ( u in , ω ) e ± j ( x ) , where α j are excitation coefficients depending on the source and independent of the space variable x .These representations are problematic for several reasons. The first one is that any representationof this type is not solution to the Helmholtz equation for ω ∈ C as soon as there are two or moremodes oscillating at different frequencies. The second problem is that in these representations, thescattered wave u − u in is not in L ( R d ) and only compact subspaces of R d can be considered. Then, arenormalisation process is necessary for the eigenmodes since they diverge exponentially. Even thoughthe study of these modes individually can give physical insight to a system (like for example by studyingthe mode volume quantity [14]), they cannot be used in frequency domain representation formulae tosolve the scattering problem. Time domain approximation of the scattered field
In the following section we show that even though they are irrelevant for frequency domain represen-tation, quasi-normal modes can be used to approximate the field in the time domain. The idea isto get around costly time domain computations by pre-computing the modes of the system and thenexpressing the response of the system to any source in terms of the modes. In the physics literature(for example [25, eq. (1.2)]) the field in the time domain is expressed under the form u ( x, t ) = (cid:60) (cid:88) j β j ( t ) e ± j ( x ) . (18)The problem with this type of expansions is that if | x | is big then e ± j ( x ) is exponentially large and thecomputation of u ( x, t ) is not very stable if the modes are pre-computed.We will show in this section that it is possible to express the scattered field in the time domainin a similar expansion, but with non-diverging, pre-computable quantities similar to the quasi-normalmodes. Here we state the main result of the paper, theorem 5.1, and discuss the result.
Let Γ k m ( · , s ) , i.e., the Green’s function for the Helmholtz equation introduced in definition A.1, be theincident wave u in in three dimensions. Given a wideband signal (cid:98) f : t (cid:55)→ (cid:98) f ( t ) ∈ C ∞ ([0 , C ]) , for C > ,we want to express the time domain response of the electric field to an oscillating dipole placed at asource point s . This means that for a fixed δ we can pick an excitation signal such that most of thefrequency content is in the low frequencies but large enough to excite the plasmonic resonances. Wecan pick η (cid:28) and ρ ≥ R ( δ ) such that ˆ R \ [ − ρ,ρ ] | f ( ω ) | d ω ≤ η,ρδc ≤ , where f : ω (cid:55)→ f ( ω ) is the Fourier transform of (cid:98) f . In practice we take ρ = R ( δ ) . The incident fieldhas the following form in the time domain: (cid:98) u in ( x, t ) = ˆ R Γ ωc ( x, s ) f ( ω ) e − iωt d ω = (cid:98) f ( t − | x − s | /c )4 π | x − s | . The goal of this section is to establish a resonance expansion for the low-frequency part of the scatteredelectric field in the time domain. Introduce, for ρ > , the truncated inverse Fourier transform of thescattered field u sca given by P ρ [ u sca ] ( x, t ) := ˆ ρ − ρ u sca ( x, ω ) e − iωt d ω. Recall that z is the centre of the resonator and δ its radius. Let us define t ± ( s, x ) := 1 c ( | s − z | + | x − z | ± δ ) ± C , the time it takes to the wideband signal to reach first the scatterer and then the observation point x .The term ± δ/c accounts for the maximal timespan spent inside the particle.13ecall the spectral decomposition in the frequency domain (proposition 3.3) for x ∈ R \ D : u sca ( x, ω ) = (cid:0) u − u in (cid:1) ( x, ω ) = J (cid:88) j =1 τ j ( ω ) (cid:104) F, ϕ j (cid:105) H ∗ ( ∂D ) S ωc D [ ϕ j ]( x ) . Theorem 5.1.
Let N ∈ N . For J ∈ N large enough, the scattered field has the following form in thetime domain for x ∈ R \ D : P ρ [ u sca ] ( x, t ) = O (cid:0) δ ρ − N (cid:1) , t ≤ t − , πi J (cid:88) j =1 C Ω ± j ( δ ) (cid:104) F, ϕ j (cid:105) H ∗ ( ∂D ) e ± j ( x ) e − i Ω ± j ( δ ) t + O (cid:18) δ t ρ − N (cid:19) , t ≥ t +0 . (19) The complex numbers Ω ± j ( δ ) are the resonant frequencies given by proposition 4.1. The fields e j are the classical quasi-normal modes defined in section 4.2. C Ω ± j ( δ ) is a constant depending only on j ,the size δ and the model for ε c ( ω ) : C Ω ± j ( δ ) := ε (cid:16) Ω ± j ( δ ) + i Ω ± j ( δ ) T − − ω p (cid:17)(cid:16) ω p δc − ) α j (cid:17) (cid:0) Ω ± j ( δ ) − Ω ∓ j ( δ ) (cid:1) . Remark 5.1.
The resonant frequencies (cid:8) Ω ± j ( δ ) (cid:9) ≤ j ≤ J have negative imaginary parts, so theorem 5.1expresses the scattered field as the sum of decaying oscillating fields. The imaginary part of Ω ± j ( δ ) accounts for absorption losses in the particle as well as radiative losses. Remark 5.2 (about the remainder ρ ) . Since for a particle of finite size δ our expansion only holds fora range of frequencies ω such that ωδc − < , we cannot compute the full inverse Fourier transformand we have a remainder that depends on the maximum frequency that we can use. Nevertheless thatmaximum frequency ρ behaves as c/δ and we can see that the remainder gets arbitrarily small for smallparticles. For a completely point-like particle one would get a zero remainder. Remark 5.3.
If we had access to the full inverse Fourier transform of the field, of course, since theinverse Fourier transform of a function which is analytic in the upper-half plane is causal we wouldfind that in the case t ≤ ( | s − z | + | x − z | − δ ) /c , (cid:98) u sca ( x, t ) = 0 . Nevertheless, our method gives theresonant frequencies only in the low-frequency regime. Therefore we only have an approximation for the low-frequency part of the scattered field, which does not have a compact support in time. Nevertheless,as shown in the numerical section 6.4.5, the low-frequency part of the scattered field is actually a goodapproximation for the scattered field. There does not seem to be any resonant frequencies for ω > R ( δ ) .This is highly non-trivial and we do not have a mathematical justification for that. Physically though,it can be explained by looking at the Drude model and noting that when ω → ∞ , ε ( ω ) −→ . The metaldoes not really interact with light at high frequencies. Even though | e ± j ( x ) | −→ ∞ when | x | → ∞ , no terms diverge in (19). Indeed we can rewrite: e ± j ( x ) e − i Ω ± j ( δ ) t = e ± j ( x ) e − i Ω ± j ( δ ) t +0 e − i Ω ± j ( δ )( t − t ) = e ± j ( x ) e − i Ω ± j ( δ )( | s − z | + | x − z | +2 δ ) c − + C e − i Ω ± j ( δ )( t − t +0 ) = C u in ,δ e ± j ( x ) e − i Ω ± j ( δ ) | x − z | c − e − i Ω ± j ( δ )( t − t +0 ) , C u in ,δ depends only on the incoming field and the particle size. We can define the following causal plasmonic quasi-normal modes ( E ± j ) j ∈ N at the complex frequency Ω ± j ( δ ) : E ± j ( x ) = S Ω ± j ( δ ) c − D [ ϕ j ]( x ) e − i Ω ± j ( δ ) | x − z | c − , x ∈ R d \ D, S Ω ± j ( δ ) √ ε c µ m D [ ϕ j ]( x ) , x ∈ D. (20) Remark 5.4.
When referring to E ± j , the term mode is inaccurate, as E ± j does not solve the Helmholtzequation. But since the ( E ± j ) j ∈ N are built from modes with a complex phase correction, we still callthem modes in a loose sense of the term. Theorem 5.1 can be re-stated:
Theorem 5.2 (alternative causal expansion) . P ρ [ u sca ] ( x, t ) = O (cid:0) δ ρ − N (cid:1) , t ≤ t − , πi J (cid:88) j =1 β δj ( u in ) E ± j ( x ) e − i Ω ± j ( δ )( t − t +0 ) + O (cid:18) δ t ρ − N (cid:19) , t ≥ t +0 , (21) where β δj ( u in ) = C Ω ± j ( δ ) (cid:104) F, ϕ j (cid:105) H ∗ ( ∂D ) e − i Ω ± j ( δ )( | s − z | +2 δ ) c − + C . Remark 5.5.
Expansion (21) has exactly the same form as the representation formula found in thephysics literature (like equation (18) ) but without any exponentially diverging quantities. The E ± j canbe computed independently of the source, just like regular quasi-normal modes. Before we can prove theorem 5.1 we need the following lemma:
Lemma 5.1. As ωδc − → , F defined in (7) admits the following asymptotic expansion: F ( x ) = 1 δ (cid:20) δ (cid:18) ε c − ε m (cid:19) ν x · ∇ Γ k m ( z − s ) + O (cid:16)(cid:0) ωδc − (cid:1) (cid:17)(cid:21) , x ∈ ∂D. (22) Proof.
See appendix C.3.
Proof of theorem 5.1.
We start by studying the time domain response of a single mode to a causalexcitation at the source point s . According to proposition 3.3 we need to compute the contribution Ξ j of each mode, that is, ˆ ρ − ρ Ξ j ( x, ω ) e − iωt d ω := ˆ ρ − ρ λ ( ω ) − λ j ( ωδ ) (cid:10) ∇ Γ ωc ( z, s ) · ν ( · ) f ( ω ) , ϕ j (cid:11) H ∗ ( ∂D ) S ωc D [ ϕ j ] e − iωt d ω, where λ j ( ωδ ) := λ j − (cid:0) ωδc − (cid:1) α j + O (cid:16)(cid:0) ωδc − (cid:1) (cid:17) . One can then write: (cid:10) ∇ Γ ωc ( z, s ) · ν ( · ) f ( ω ) , ϕ j (cid:11) H ∗ ( ∂D ) S ωc D [ ϕ j ] = f ( ω ) (cid:18) | z − s | − i ωc (cid:19) (cid:18) λ j − (cid:19) ˆ ∂D × ∂D vϕ j ( v ) ϕ j ( y )16 π | x − y || z − s | e i ωc ( | x − y | + | z − s | ) d σ ( v )d σ ( y ) , where we used (cid:104) ν, ϕ j (cid:105) H ∗ ( ∂D ) = (1 / − λ j ) (cid:104) x, ϕ j (cid:105) , − [5]. Since (cid:104) ν, ϕ (cid:105) H ∗ ( ∂D ) = 0 , the zeroth termvanishes in the summation. 15ow we want to apply the residue theorem to get an asymptotic expansion in the time domain.Note that: ˆ ρ − ρ Ξ j ( x, ω ) e − iωt d ω = ˛ C ± Ξ j ( x, Ω) e − i Ω t dΩ − ˆ C ± ρ Ξ j ( x, Ω) e − i Ω t dΩ , where the integration contour C ± ρ is a semicircular arc of radius ρ in the upper (+) or lower (-) half-plane, and C ± is the closed contour C ± = C ± ρ ∪ [ − ρ, ρ ] . The integral on the closed contour is the maincontribution to the scattered field by the mode and can be computed using the residue theorem to get,for ρ ≥ (cid:60) [Ω ± j ( δ )] , ˛ C + Ξ j ( x, Ω) e − i Ω t dΩ = 0 , ˛ C − Ξ j ( x, Ω) e − i Ω t dΩ = 2 πi Res (cid:0) Ξ j ( x, Ω) e − i Ω t , Ω ± j ( δ ) (cid:1) . Since Ω ± j ( δ ) is a simple pole of ω (cid:55)→ λ ( ω ) − λ j ( ωδ ) we can write: ˛ C − Ξ j ( x, Ω) e − i Ω t dΩ = 2 πi Res (cid:0) Ξ j ( x, Ω) , Ω ± j ( δ ) (cid:1) e − i Ω ± j ( δ ) t . To compute the integrals on the semi-circle, we introduce: B j ( y, v, Ω) = λ j − / λ ( ω ) − λ j ( δ Ω) (cid:18) | z − s | − i ωc (cid:19) ˆ ∂D × ∂D vϕ j ( v ) ϕ j ( y )16 π | x − y || z − s | ( y, v ) ∈ ( ∂D ) . Note that B j ( · , · , Ω) behaves like a polynomial in Ω when | Ω | → ∞ . Given the regularity of theinput signal (cid:98) f ∈ C ∞ ([0 , C ]) , the Paley-Wiener theorem [42, p. 161] ensures decay properties of itsFourier transform at infinity. For all N ∈ N ∗ there exists a positive constant C N such that for all Ω ∈ C | f (Ω) | ≤ C N (1 + | Ω | ) − N e C |(cid:61) (Ω) | . Let T := ( | x − y | + | s − v | ) /c . We now re-write the integrals on the semi-circle ˆ C ± ρ Ξ j ( x, Ω) e − i Ω t dΩ = ˆ C ± ρ f (Ω) ˆ ∂D × ∂D B j ( y, v, Ω) e i Ω( T − t ) d σ ( v )d σ ( y )dΩ . We have that t − + C ≤ T ≤ t +0 − C . Two cases arise. Case 1:
For < t < t − , i.e., when the signal emitted at s has not reached the observation point x ,we choose the upper-half integration contour C + . Transforming into polar coordinates, Ω = ρe iθ for θ ∈ [0 , π ] , we get: (cid:12)(cid:12)(cid:12) e i Ω( T − t ) (cid:12)(cid:12)(cid:12) ≤ e − ( t − − t + C ) (cid:61) (Ω) ∀ ( y, v ) ∈ ( ∂D ) , and (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˆ C + ρ Ξ j ( x, Ω) e − i Ω t dΩ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ ˆ π ρ (cid:12)(cid:12) f (cid:0) ρe iθ (cid:1)(cid:12)(cid:12) e − ρ ( t − − t + C ) sin θ ˆ ∂D × ∂D (cid:12)(cid:12) B j (cid:0) y, v, ρe iθ (cid:1)(cid:12)(cid:12) d σ ( v )d σ ( y )d θ, ≤ ρC N (1 + ρ ) − N δ max θ ∈ [0 ,π ] (cid:13)(cid:13) B j (cid:0) · , · , ρe iθ (cid:1)(cid:13)(cid:13) L ∞ ( ∂D × ∂D ) π − e − ρ ( t − − t ) ρ ( t − − t ) , where we used that for θ ∈ [0 , π/ , we have sin θ ≥ θ/π ≥ and − cos θ ≤ − θ/π . The usualway to go forward from here is to take the limit ρ → ∞ , and get that the limit of the integral on the16emi-circle is zero. However, we work in the quasi-static approximation here, and our modal expansionis not uniformly valid for all frequencies. So we have to work with a fixed maximum frequency ρ ..Since N can be taken arbitrarily large and that B j behaves like a polynomial in ρ whose degree doesnot depend on j , we get that, uniformly for j ∈ [1 , J ] : (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˆ C + ρ Ξ j ( x, Ω) e − i Ω t dΩ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = O (cid:18) δ t − − t ρ N (cid:19) . Of course if one has to consider the full inverse Fourier transform of the scattered electromagneticfield, by causality, one should expect the limit to be zero. However, one would need high-frequencyestimates of the electromagnetic field, as well as a modal decomposition that is uniformly valid for allfrequencies. Since our modal expansion is only valid for a limited range of frequencies we get an errorbound that is arbitrarily small if the particle is arbitrarily small, but not strictly zero.
Case 2:
For t > t +0 , we choose the lower-half integration contour C − . Transforming into polarcoordinates, Ω = ρe iθ for θ ∈ [ π, π ] , we get (cid:12)(cid:12)(cid:12) e i Ω( T − t ) (cid:12)(cid:12)(cid:12) ≤ e ( t − t +0 − C ) (cid:61) (Ω) ∀ ( y, v ) ∈ ( ∂D ) , and (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˆ C − ρ Ξ j ( x, Ω) e − i Ω t dΩ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ ˆ ππ ρ (cid:12)(cid:12) f (cid:0) ρe iθ (cid:1)(cid:12)(cid:12) e ρ ( t − t +0 − C ) sin θ ˆ ∂D × ∂D (cid:12)(cid:12) B j (cid:0) y, v, ρe iθ (cid:1)(cid:12)(cid:12) d σ ( v )d σ ( y )d θ, ≤ ρC N (1 + ρ ) − N δ max θ ∈ [ π, π ] (cid:13)(cid:13) B j (cid:0) · , · , ρe iθ (cid:1)(cid:13)(cid:13) L ∞ ( ∂D × ∂D ) π − e − ρ ( t − t +0 ) ρ ( t − t +0 ) . Exactly as in Case , we cannot take the limit ρ → ∞ . Using the fact that N can be taken arbitrarilylarge and that B j behaves like a polynomial in ρ whose degree does not depend on j , we get that,uniformly for j ∈ [1 , J ] : (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˆ C − ρ Ξ j ( x, Ω) e − i Ω t dΩ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = O (cid:18) δ t ρ − N (cid:19) . The result of theorem 5.1 is obtained by summing the contribution of all the modes considered.
Remark 5.6.
The fact that we work with a finite number of modes is necessary for the perturbationtheory of section 3 but also in this section. Indeed, if we consider all the modes there is an accumulationpoint in the poles of the modal expansion of the field, and therefore we cannot apply the residue theorem.
In two dimensions, the Green’s function does not have an explicit phase term, so we need to introduceanother asymptotic parameter (cid:15) > to be able to use the large argument asymptotics of the Hankelfunction. Our new truncated inverse Fourier transform of the scattered field u sca given by P ρ,(cid:15) [ u sca ] ( x, t ) = ˆ − (cid:15) − ρ u sca ( x, ω ) e − iωt d ω + ˆ ρ(cid:15) u sca ( x, ω ) e − iωt d ω. This allows us to define a notion of far field . A point x is far from D if (cid:15) | x − z | c − (cid:29) . We can nowadd two additional hypotheses:• the source is far away from the particle (or equivalently, the incoming wave is a plane wave)17 the observation point is far away from the particle.The incident field has the following form in the time domain: (cid:98) u in ( x, t ) = (cid:98) f (cid:18) t − d · xc (cid:19) . (23)Besides these two assumptions and a difference in the order of the remainder, the result in twodimensions is essentially the same as in three dimensions. Theorem 5.3.
Let N ∈ N . For J large enough the scattered field has the following form in the timedomain for x far away from D : P ρ,(cid:15) [ u sca ] ( x, t ) = O (cid:0) δρ − N (cid:1) , t ≤ t − , πi J (cid:88) j =1 C Ω ± j ( δ ) (cid:104) F, ϕ j (cid:105) H ∗ ( ∂D ) e ± j ( x ) e − i Ω ± j ( δ ) t + O (cid:18) δt ρ − N (cid:19) , t ≥ t +0 , (24) with Ω ± j ( δ ) being the plasmonic resonant frequencies of the particle given by proposition 4.2. C Ω ± j ( δ ) is a constant depending only on j , the size δ and the model for ε c ( ω ) : C Ω ± j ( δ ) := ε (cid:16) Ω ± j ( δ ) + i Ω ± j ( δ ) T − − ω p (cid:17)(cid:16) ω p δc − ) log (cid:0) Ω ± s,j δc − (cid:1) α j (cid:17) (cid:0) Ω ± j ( δ ) − Ω ∓ j ( δ ) (cid:1) . Proof.
The proof is quite similar to the three-dimensional case. It is included in appendix D for thesake of completeness.
The goal of this section is to illustrate the validity of our approach and to show that the approximationseems to be working with less restrictive hypotheses than the ones in theorem 5.3:• for more general shapes (non-convex or non-algebraic)• closer to the particle (outside of the far field approximation).For these simulations we build upon the codes for the layer potentials developed in [40].
Throughout this section, we consider the three domains sketched on Figure 1 to illustrate our results:
Rounded diamond:
The rounded diamond (a) is defined by the parametric curve ζ ( θ ) = 2 (cid:0) e iθ + 0 . e − iθ (cid:1) ,for θ ∈ [0 , π ] . It is an algebraic domain of class Q from [6]. This shape satisfies condition 1, as wellas the hypotheses of theorem 5.3. Narrow ellipse:
The ellipse (b) semi-axes are on the X - and X - axes and are of length a = 1 and b = 5 , respectively. It is algebraic but not asymptotically a circle in the sense of [6].18 ymbol Magnitude ω p · Hz T − s ε . · − Fm − µ π · − Hm − δ − md (1 / √ , / √ z (0 , Table 1: Physical constants and parameters values.
Five-petal flower:
The flower (c) is defined by (cid:37) = 2 + 0 . θ ) in polar coordinates. It hasCartesian equation . (cid:0) X + X (cid:1) − . X (cid:0) X + X (cid:1) + 6 X (cid:0) X + X (cid:1) − . X (cid:0) X + X (cid:1) − (cid:0) X + X (cid:1) / = 0 in the rescaled ( X , X ) plane. So it is not algebraic (due to the non-integer power of the last term)and not convex. We have no theoretical results on the number of modes that radiate. (b)
F , (cid:101) φ j (cid:105) H ∗ ( ∂B ) decay very rapidly when d = 3 . Ina two-dimensional setting, the theoretical framework is not as clear, but we check numerically thatthe contribution the modes decrease quite fast with j . Recall that the weight of the j th mode isgiven by the scalar product (cid:104) (cid:101) F , (cid:101) φ j (cid:105) H ∗ ( ∂B ) , which, in a low-frequency regime, can be approximatedas (cid:104) ν · ∇ u in , (cid:101) φ j (cid:105) H ∗ ( ∂B ) (see lemma D.2). On panel (a) of Figure 2 we show on all examples that (cid:104) ν · ∇ u in , (cid:101) φ j (cid:105) H ∗ ( ∂B ) decays as j grows. We average over all possible directions d of the incident field.19anel (b) of the same picture shows that the modes themselves, S ω p δc − B [ (cid:101) φ j ]( X ) , decrease as j increases.We average here over all observation positions, X belongs to a circle of radius centred at z = 0 . EllipseFlowerDiamond
EllipseFlowerDiamond ( a )
21y setting ω l such that: − ρδc − = ω (cid:48)− L < ω (cid:48)− L +1 < . . . < ω (cid:48)− = − (cid:15)δc − , (cid:15)δc − = ω (cid:48) < . . . < ω (cid:48) L − < ω (cid:48) L = ρδc − , with ω (cid:48) l +1 − ω (cid:48) l = ( ρ − (cid:15) ) δc − /L for every l ∈ [ − L − , − ∪ [1 , L − . We compute the scattered field inthe frequency domain using the representation formula (3). The single layer potential is approximatedusing N = 2 equally-spaced discretisation points along the boundary ∂B . We define the dimensionlessfrequency ω (cid:48) = ωδc − . The reference solution is computed by taking the truncated inverse Fouriertransform P ρ,(cid:15) [ (cid:98) u sca ] ( X, t ) ≈ ( ρ − (cid:15) ) L L (cid:88) l =1 (cid:16) e − iω (cid:48)− l cδ − t u sca (cid:0) X, ω (cid:48)− l cδ − (cid:1) + e − iω (cid:48) l cδ − t u sca (cid:0) X, ω (cid:48) l cδ − (cid:1)(cid:17) . (25) The expansion is obtained by summing the first J = 30 modes. Using theorem 5.3, the modal approx-imation of order J becomes: U J ( X, t ) =2 πi J (cid:88) j =1 ε (cid:16) Ω + j ( δ ) + i Ω + j ( δ ) T − − ω p (cid:17) (cid:104) (cid:101) F , (cid:101) φ j (cid:105) H ∗ ( ∂B ) (cid:16) ω p δc − ) log (cid:0) Ω + s,j δc − (cid:1) α j (cid:17) (cid:0) Ω + j ( δ ) − Ω − j ( δ ) (cid:1) δ S Ω + j ( δ ) δc − B [ (cid:101) φ j ]( X ) e − i Ω + j ( δ ) t + (cid:16) Ω − j ( δ ) + i Ω − j ( δ ) T − − ω p (cid:17) (cid:104) (cid:101) F , (cid:101) φ j (cid:105) H ∗ ( ∂B ) (cid:16) ω p δc − ) log (cid:0) Ω − s,j δc − (cid:1) α j (cid:17) (cid:0) Ω − j ( δ ) − Ω + j ( δ ) (cid:1) δ S Ω − j ( δ ) δc − B [ (cid:101) φ j ]( X ) e − i Ω − j ( δ ) t . (26)The simulation results are shown in figures 5, 6, 7 and 8. To corroborate our pole expansion, weplot the real part of the reference solution (25) against the real part of the asymptotic one (26) for thedifferent domains and from different observation points. We begin with the diamond, since it is the shape that satisfies the hypotheses of theorem 5.3. Figure 5shows the field scattered by the diamond, measured in the far-field at position X = D . The referencesolution is nicely approximated by the sum of four modes (4, 5, 6 and 7). Reference solutionAsymptotic solution t +0
Reference solutionAsymptotic solution (a) X = A
Reference solutionLarge frequency solution t +0
Acknowledgement
This work was supported by the Swiss National Science Foundation grant number 200021-172483. Theauthors thank Habib Ammari for helpful conversations.
Data availability
The data supporting the findings of this study were generated through matlab and are available fromthe corresponding author on request. 26
Properties of the layer potentials
We briefly recall here some basic properties of layer potential. There is an abundant literature on thesubject. For more details we refer to the books [3, 32, 16, 4].
A.1 Definitions and notations
Definition A.1.
Denote by Γ k the outgoing Green’s function for the homogeneous medium, i.e., theunique solution of the Helmholtz operator: (cid:0) ∆ + k (cid:1) Γ k ( · , y ) = δ y ( · ) in R d satisfying the Sommerfeld radiation condition. In three dimensions, Γ k is given by Γ k ( x, y ) = − e ik | x − y | π | x − y | , x, y ∈ R . In two dimensions, it is given by Γ k ( x, y ) = π log | x − y | , if k = 0 , − i H (1)0 ( k | x − y | ) , if k > , for x, y ∈ R , where H (1)0 is the well-known Hankel function of the first kind and order . Lemma A.1.
The Hessian matrix of the outgoing fundamental solution in three dimensions D x Γ k ( x, z ) =( D ) p,q =1 is with entries D pp = e ik | x − z | π | x − z | (cid:2) | x − z | − x p − z p ) + 3 ik ( x p − z p ) | x − z | + k | x − z | − ik | x − z | (cid:3) ,D pq = e ik | x − z | π | x − z | ( x p − z p )( x q − z q ) (cid:2) − ik | x − z | + k | x − z | (cid:3) , for p (cid:54) = q. Definition A.2.
For a function φ ∈ L ( ∂D ) , we define the single-layer potential by S kD [ φ ]( x ) = ˆ ∂D Γ k ( x, y ) φ ( y )d σ ( y ) , x ∈ R d , and the Neumann-Poincaré operator by K k, ∗ D [ φ ]( x ) = ˆ ∂D ∂ Γ k ( x, y ) ∂ν ( x ) φ ( y )d σ ( y ) , x ∈ ∂D, When k = 0 , we just write S D and K ∗ D for simplicity. A.2 The Calderón identity and symmetrisation of K ∗ D Lemma A.2.
We recall the following classical results [21, 9, 32].1. The following Plemelj’s symmetrisation principle identity (also known as Calderón) holds: K D S D = S D K ∗ D on H − ( ∂D ) . (27)
2. If ∂D ∈ C ,α , for some α > , then K ∗ D is compact. Let ( λ j , ϕ j ) j ∈ N , be the eigenvalues andnormalised eigenfunctions of K ∗ D in H ∗ ( ∂D ) . Then λ j ∈ ] − / , / , λ = 1 / and λ j → as j → ∞ . . The operator K ∗ D is self-adjoint in the Hilbert space H ∗ ( ∂D ) which is H − ( ∂D ) equipped withthe following inner product: (cid:104) u, v (cid:105) H ∗ ( ∂D ) = (cid:40) −(cid:104) u, (cid:101) S D [ v ] (cid:105) − , , d = 2 , − (cid:104) u, S D [ v ] (cid:105) − , , d = 3 , where (cid:101) S D [ v ] = (cid:40) S D [ v ] if (cid:104) v, χ ( ∂D ) (cid:105) − , = 0 , − χ ( ∂D ) if v = ϕ , with ϕ being the unique (in the case of a single particle) eigenfunction of K ∗ D associated witheigenvalue / such that (cid:104) ϕ , χ ( ∂D ) (cid:105) − , = 1 . Also, (cid:104)· , ·(cid:105) − , is the duality pairing between H − ( ∂D ) and H ( ∂D ) .4. From [9], we have the following extension of (27) in two dimensions: K D (cid:101) S D = (cid:101) S D K ∗ D , on H − ( ∂D ) .
5. Since K D [ χ ( ∂D )] = 12 χ ( ∂D ) , it holds that ˆ ∂D ϕ j = 0 , for j (cid:54) = 0 .
6. The following trace formulae hold for φ ∈ H − ( ∂D ) : S kD [ φ ] | + = S kD [ φ ] | − ,∂ S kD [ φ ] ∂ν (cid:12)(cid:12)(cid:12)(cid:12) ± = (cid:18) ± I + K k, ∗ D (cid:19) [ φ ] , where I denotes the identity operator.7. The following representation formula holds: K ∗ D [ φ ] = ∞ (cid:88) j =0 λ j (cid:104) φ, ϕ j (cid:105) H ∗ ( ∂D ) ϕ j , ∀ φ ∈ H ∗ ( ∂D ) . A.3 Invertibility of the boundary operators
Lemma A.3.
For k small enough, the three-dimensional single-layer potential S kD : H − / ( ∂D ) → H / ( ∂D ) is invertible. S D is also invertible. Lemma A.4. S D : H − / ( ∂D ) → H / ( ∂D ) is invertible in three dimensions. In two dimensions, the single-layer potential S D : H − / ( ∂D ) → H / ( ∂D ) is, in general, notinvertible. All the proofs for the following lemmas can be found in [5]. Lemma A.5.
For k small enough, the two-dimensional boundary operator (cid:98) S kD : H ∗ ( ∂D ) → H ∗ ( ∂D ) defined as (cid:98) S kD [ φ ]( x ) = S D [ φ ]( x ) + η k ˆ ∂D φ ( y )d σ ( y ) , (28) is invertible and (cid:16) (cid:98) S kD (cid:17) − = (cid:101) S − D − (cid:68) (cid:101) S − D [ · ] , ϕ (cid:69) H ∗ ( ∂D ) ϕ − U k , (29) where U k = (cid:68) (cid:101) S − D [ · ] , ϕ (cid:69) H ∗ ( ∂D ) S D [ ϕ ] + η k ϕ and η k = (1 / π )(log k + γ − log 2) − i/ , with the constant γ beingthe Euler constant. Note that U k = O (1 / log k ) . emma A.6. For k small enough, the two-dimensional single-layer potential S kD : H ∗ ( ∂D ) → H ∗ ( ∂D ) is invertible. B Scaling properties for a finite volume particle
For each function f defined on ∂D , we define a corresponding function on ∂B by (cid:101) f ( X ) = f ( z + δX ) . Lemma B.1.
It holds that K k, ∗ D [ f ]( x ) = K kδ, ∗ B [ (cid:101) f ]( X ) , S kD [ f ]( x ) = δ S kδB [ (cid:101) f ]( X ) . Lemma B.2.
For f, g defined on ∂D , corresponding to (cid:101) f , (cid:101) g , respectively, we have (cid:104) f, g (cid:105) H ∗ ( ∂D ) = δ (cid:68) (cid:101) f , (cid:101) g (cid:69) H ∗ ( ∂B ) , d = 2 ,δ (cid:68) (cid:101) f , (cid:101) g (cid:69) H ∗ ( ∂B ) , d = 3 , || f || H ∗ ( ∂D ) = (cid:40) δ || (cid:101) f || H ∗ ( ∂B ) , d = 2 ,δ / || (cid:101) f || H ∗ ( ∂B ) , d = 3 . Proof.
In three dimensions, by straightforward calculations we have (cid:104) f, g (cid:105) H ∗ ( ∂D ) = ˆ ∂D f ( x ) ˆ ∂D g ( y )4 π | x − y | d σ ( y )d σ ( x )= δ ˆ ∂B f ( z + δX ) ˆ ∂B g ( z + δY )4 π | X − Y | d σ ( Y )d σ ( X )= δ (cid:68) (cid:101) f , (cid:101) g (cid:69) H ∗ ( ∂B ) . Hence, || f || H ∗ ( ∂D ) = δ / || (cid:101) f || H ∗ ( ∂B ) .In the two-dimensional case we write H ∗ ( ∂D ) = H ∗ ( ∂D ) ⊕ { µϕ , µ ∈ C } and treat both cases: g belongs to either H ∗ ( ∂D ) or { µϕ , µ ∈ C } . In the former case, we have (cid:104) f, g (cid:105) H ∗ ( ∂D ) = − π ˆ ∂D f ( x ) ˆ ∂D g ( y ) log( | x − y | )d σ ( y )d σ ( x )= − δ π ˆ ∂B f ( z + δX ) ˆ ∂B g ( z + δY )(log( δ ) + log( | X − Y | ))d σ ( Y )d σ ( X )= δ (cid:68) (cid:101) f , (cid:101) g (cid:69) H ∗ ( ∂B ) . If g = µϕ , we have (cid:104) f, g (cid:105) H ∗ ( ∂D ) = ˆ ∂D µf ( x )d σ ( x )= δ ˆ ∂B µf ( z + δX )d σ ( X )= δ (cid:68) (cid:101) f , (cid:101) g (cid:69) H ∗ ( ∂B ) , where the last equality follows from the fact that δ (cid:101) ϕ is the (unique) eigenfunction of K ∗ B associatedwith eigenvalue / such that (cid:104) δ (cid:101) ϕ , χ ( ∂B ) (cid:105) − , = 1 . Hence, || f || H ∗ ( ∂D ) = δ || (cid:101) f || H ∗ ( ∂B ) .29 Asymptotic expansions
C.1 Asymptotic expansions of the boundary operators
Lemma C.1.
1. The three-dimensional single-layer potential and its inverse admit the followingexpansions in the quasi-static limit kδ → : S kδB = S B + kδ S B, + ( kδ ) S B, + O (cid:16) ( kδ ) (cid:17) , (cid:0) S kδB (cid:1) − = S − B + kδ B B, + ( kδ ) B B, + O (cid:16) ( kδ ) (cid:17) , where, for φ ∈ H − ( ∂B ) , S B,j [ φ ]( x ) = − i π ˆ ∂B ( i | x − y | ) j − j ! φ ( y )d σ ( y ) , x ∈ R , for j ∈ N . Also B B, = −S − B S B, S − B and B B, = −S − B S B, S − B + S − B S B, S − B S B, S − B .2. The two-dimensional single-layer potential and its inverse admit the following expansions in thequasi-static limit kδ → : S kδB = (cid:98) S kδB + ( kδ ) log( kδ ) S (1) B, + O (cid:16) ( kδ ) (cid:17) , (cid:0) S kδB (cid:1) − = L B + U kδ − ( kδ ) log( kδ ) L B S (1) B, L B + O (cid:16) ( kδ ) (cid:17) , where, for φ ∈ H − ( ∂B ) , S (1) B, [ φ ]( x ) = − π ˆ ∂B | x − y | φ ( y )d σ ( y ) , x ∈ R . Also P H ∗ is the orthogonal projection onto H ∗ , L B = P H ∗ (cid:101) S − B .3. The Neumann-Poincaré operator in three dimensions admits the following expansion in the quasi-static limit K kδ, ∗ B = K ∗ B + ( kδ ) K ∗ B, + O (cid:16) ( kδ ) (cid:17) , where, for φ ∈ H − ( ∂B ) , K ∗ B, [ φ ]( x ) = 18 π ˆ ∂B ( x − y ) · ν ( x ) | x − y | φ ( y )d σ ( y ) , x ∈ ∂B.
4. The Neumann-Poincaré operator in two dimensions admits the following expansion in the quasi-static limit K kδ, ∗ B = K ∗ B + ( kδ ) log( kδ ) K (1) B, + O (cid:16) ( kδ ) (cid:17) , where, for φ ∈ H − ( ∂B ) , K (1) B, [ φ ]( x ) = − π ˆ ∂B ∂ | x − y | ∂ν x φ ( y )d σ ( y ) , x ∈ ∂B. Proof.
The proof can be found in [5]. 30 .2 Proof of lemma 2.1
Proof.
Recall that A ωδ/cB = 1 ε m (cid:18) I + K k m δ, ∗ B (cid:19) + 1 ε c (cid:18) I − K k c δ, ∗ B (cid:19) (cid:16) S k c δB (cid:17) − S k m δB . For the three-dimensional case, using lemma C.1 we have by straightforward calculations A ωδ/cB = 1 ε m (cid:20) I + K ∗ B + ( k m δ ) K B, (cid:21) + 1 ε c (cid:20) I − K ∗ B − ( k c δ ) K B, (cid:21) (cid:2) S − B + ( k c δ ) B B, + ( k c δ ) B B, (cid:3) , (cid:2) S B + ( k m δ ) S B, + ( k m δ ) S B, (cid:3) + O (cid:16)(cid:0) ωδc − (cid:1) (cid:17) , = 12 (cid:18) ε c + 1 ε m (cid:19) I − (cid:18) ε c − ε m (cid:19) K ∗ B + (cid:18) ( k m δ ) ε m − ( k c δ ) ε c (cid:19) K B, + 1 ε c (cid:18) − K ∗ B (cid:19) (cid:18) ( k m δ ) S − B S B, , +( k m δ ) S − B S B, − k c δ S − B S B, − k c k m δ S − B S B, S − B S B, + ( k c δ ) S − B S B, S − B S B, , − ( k c δ ) S − B S B, (cid:19) + O (cid:16)(cid:0) ωδc − (cid:1) (cid:17) , = A B + 1 ε c (cid:18) I − K ∗ B (cid:19) ( k m − k c ) δ S − B S B, + O (cid:16)(cid:0) ωδc − (cid:1) (cid:17) , = A B + (cid:0) ωδc − (cid:1) ε m − ε c ε m ε c (cid:18) I − K ∗ B (cid:19) S − B S B, + O (cid:16)(cid:0) ωδc − (cid:1) (cid:17) , where we used S − B S B = I and (cid:18) I − K ∗ B (cid:19) S − B S B, = 0 . For the two-dimensional case, we have I − K k c δ, ∗ B = 12 I − K ∗ B − ( k c δ ) log( k c δ ) K (1) B, + O (cid:0) ( k c δ ) (cid:1) , (cid:16) S k c δB (cid:17) − = L B + U k c δ − (cid:0) ωδc − (cid:1) log (cid:0) ωδc − (cid:1) ε c ε m L B S (1) B, L B + O (cid:16)(cid:0) ωδc − (cid:1) (cid:17) , S k m δB = (cid:101) S B + Υ k m δ + (cid:0) ωδc − (cid:1) log (cid:0) ωδc − (cid:1) S (1) B, + O (cid:16)(cid:0) ωδc − (cid:1) (cid:17) . Also, L B Υ k m δ = P H ∗ (cid:101) S − B Υ k m δ = 0 , where Υ k m δ [ ψ ] = (cid:104) ψ, (cid:101) φ (cid:105) H ∗ ( S B [ (cid:101) φ ] + χ ( ∂B ) + η k m δ ) . Hence, ( S k c δB ) − S k m δB = P H ∗ + U k c δ (cid:101) S B + U k c δ Υ k m δ + (cid:0) ωδc − (cid:1) log (cid:0) ωδc − (cid:1) L B S (1) D, (cid:18) I − ε c ε m P H ∗ (cid:19) + O (cid:16)(cid:0) ωδc − (cid:1) (cid:17) . We have that (cid:18) I − K ∗ B (cid:19) U k c δ = (cid:68) (cid:101) S − B [ · ] , (cid:101) φ (cid:69) H ∗ ( ∂B ) S B [ (cid:101) φ ] + η k c (cid:18) (cid:101) φ − K ∗ B [ (cid:101) φ ] (cid:19) = 0 . .3 Proof of lemma 5.1 Proof.
Using the Taylor expansion Γ k m ( z + δX, y ) (cid:12)(cid:12) X ∈ ∂B = Γ k m ( z, y ) + δX · ∇ Γ k m ( z, y ) + δ X (cid:62) D X Γ k m ( z, y ) X + ..., we compute, for X ∈ ∂B , (cid:101) F ( X ) = (cid:101) F ( X ) + 1 δε c (cid:18) I − K k c δ, ∗ B (cid:19) (cid:16) S k c δB (cid:17) − [ (cid:101) F ]( X )= − δε m ∂ Γ k m ( z + δX, y ) ∂ν X − δε c (cid:18) I − K ∗ B + O (cid:0) ( k c δ ) (cid:1)(cid:19) (cid:0) S − B + k c δ B B, + O (cid:0) ( k c δ ) (cid:1)(cid:1)(cid:2) Γ k m ( z + δX, y ) (cid:3) = − δε m ν X · ∇ Γ k m ( z + δX, y ) − δε c (cid:18) I − K ∗ B (cid:19) S − B (cid:2) Γ k m ( z + δX, y ) (cid:3) − k c ε c (cid:18) I − K ∗ B (cid:19) B B, (cid:2) Γ k m ( z + δX, y ) (cid:3) + O (cid:0) ω δc − (cid:1) = − ε m ν X · ∇ Γ k m ( z, y ) − Γ k m ( z, y ) δε c (cid:18) I − K ∗ B (cid:19) S − B [ χ ( ∂B )] − ∇ Γ k m ( z, y ) ε c · (cid:18) I − K ∗ B (cid:19) S − B [ X ] + O (cid:0) ω δc − (cid:1) = − ε m ν X · ∇ Γ k m ( z, y ) + 1 ε c ν X · ∇ Γ k m ( z, y ) + O (cid:0) ω δc − (cid:1) , where we used (cid:0) I − K ∗ B (cid:1) S − B [ χ ( ∂B )] = 0 and (cid:0) I − K ∗ B (cid:1) B B, = 0 . It is immediate to see that (cid:0) I − K ∗ B (cid:1) S − B [ X ] = − ν X , indeed assuming there exists φ ∈ H ∗ ( ∂B ) such that S − B [ X ] = φ , then S B [ φ ] = X , and ∂ S B [ φ ] /∂ν X | − = ν X which is equivalent to (cid:0) I − K ∗ B (cid:1) [ φ ] = − ν X using jump condi-tions. C.4 Proof of lemma D.2
Proof.
Using the Taylor expansion e ik m d · ( z + δX ) = e ik m d · z + iωc [ d · δX ] e ik m d · z + O (cid:16) ( k m δ ) (cid:17) , we compute, for X ∈ ∂B , (cid:101) F ( X ) = (cid:101) F ( X ) + 1 δε c (cid:18) I − K k c δ, ∗ B (cid:19) (cid:16) S k c δB (cid:17) − [ (cid:101) F ]( X )= − δε m ∂e ik m d · ( z + δX ) ∂ν X − δε c (cid:18) I − K ∗ B + O (cid:16) ( k c δ ) log( k c δ ) (cid:17)(cid:19) (cid:16) L B + U k c δ + O (cid:16) ( k c δ ) log( k c δ ) (cid:17)(cid:17)(cid:104) e ik m d · ( z + δX ) (cid:105) = − δε m ν X · ∇ e ik m d · ( z + δX ) − δε c (cid:18) I − K ∗ B (cid:19) (cid:101) S − B (cid:104) e ik m d · ( z + δX ) (cid:105) + O (cid:0) ω δc − log (cid:0) ωδc − (cid:1)(cid:1) = − ε m ν X · ∇ e ik m d · z − e ik m d · z δε c (cid:18) I − K ∗ B (cid:19) (cid:101) S − B [ χ ( ∂B )] − ∇ e ik m d · z ε c · (cid:18) I − K ∗ B (cid:19) (cid:101) S − B [ X ] + O (cid:0) ω δc − log (cid:0) ωδc − (cid:1)(cid:1) = − ε m ν X · ∇ e ik m d · z + 1 ε c ν X · ∇ e ik m d · z + O (cid:0) ω δc − log (cid:0) ωδc − (cid:1)(cid:1) , (cid:0) I − K ∗ B (cid:1) (cid:101) S − B [ χ ( ∂B )] = 0 and (cid:0) I − K ∗ B (cid:1) U kδ = 0 . It is immediate to see that (cid:0) I − K ∗ B (cid:1) (cid:101) S − B [ X ] = − ν X , indeed assuming there exists φ ∈ H ∗ ( ∂B ) such that (cid:101) S − B [ X ] = φ , then (cid:101) S B [ φ ] = X so ∂ (cid:101) S B [ φ ] /∂ν X (cid:12)(cid:12)(cid:12) − = ν X and using jump conditions for φ (cid:54) = ϕ we get ∂ (cid:101) S B [ φ ] /∂ν X (cid:12)(cid:12)(cid:12) − = ∂ S B [ φ ] /∂ν X | − = (cid:0) − I + K ∗ B (cid:1) [ φ ] = ν X . D Proof of theorem 5.3
We need the following lemma:
Lemma D.1.
The Hankel function has the following asymptotics as x → + ∞ : H (1)0 ( x ) = (cid:114) πx e i ( x − π/ + O (cid:18) x (cid:19) . (30)For x large and y ∈ ∂D : S ωc D [ ϕ j ] ( x ) = − i ˆ ∂D H (1)0 (cid:0) ωc − | x − y | (cid:1) ϕ j ( y )d σ ( y ) ∼ − i √ √ π ˆ ∂D e i ( ωc − | x − y |− π/ ) (cid:112) ωc − | x − y | ϕ j ( y )d σ ( y ) . Lemma D.2. As ωδc − → , F defined in (7) admits the following asymptotic expansion: F ( x ) = f ( ω ) δ (cid:20) iωδc − e iωc − d · z (cid:18) ε c − ε m (cid:19) d · ν x + O (cid:16)(cid:0) ωδc − (cid:1) log (cid:0) ωδc − (cid:1)(cid:17)(cid:21) , x ∈ ∂D. Proof.
See appendix C.4.
Lemma D.3. As ωδc − → , the scalar field admits the following asymptotic expansion: u sca ( x, ω ) ≈ e − iπ/ √ πc J (cid:88) j =1 (cid:104) d · ν, ϕ j (cid:105) H ∗ ( ∂D ) Ξ j ( x, ω ) , where the modes Ξ j are defined by Ξ j ( x, ω ) := ˆ ∂D ϕ j ( y ) (cid:112) | x − y | f ( ω ) √ ω ( λ ( ω ) − λ j ( ωδ )) e iωc − ( | x − y | + d · z ) d σ ( y ) , and λ j ( ωδ ) := λ j − (cid:0) ωδc − (cid:1) log (cid:0) ωδc − (cid:1) α j + O (cid:16)(cid:0) ωδc − (cid:1) (cid:17) .Proof. Since (cid:104) ν, ϕ (cid:105) H ∗ ( ∂D ) = 0 , the zeroth term vanishes in the summation.The goal of this section is to establish a resonance expansion for the low-frequency part of thescattered field in the time domain. Introduce, for < (cid:15) < ρ , the truncated inverse Fourier transformof the scattered field u sca given by P ρ,(cid:15) [ (cid:98) u sca ] ( x, t ) = ˆ − (cid:15) − ρ u sca ( x, ω ) e − iωt d ω + ˆ ρ(cid:15) u sca ( x, ω ) e − iωt d ω. Recall that z is the centre of the resonator and δ its radius. Let us define t ± ( d, x ) := 1 c ( | x − z | + d · z ± δ ) ± C , the time it takes to the signal to reach first the scatterer and then observation point x . The term ± δ/c accounts for the maximal timespan spent inside the particle.33 roof. We have P ρ,(cid:15) [ (cid:98) u sca ] ( x, t ) ∼ J (cid:88) j =1 e − iπ/ √ πc (cid:104) d · ν, ϕ j (cid:105) H ∗ ( ∂D ) (cid:20) ˆ − (cid:15) − ρ Ξ j ( x, ω ) e − iωt d ω + ˆ ρ(cid:15) Ξ j ( x, ω ) e − iωt d ω (cid:21) . (31)For j ≥ , let us compute the contribution of one mode Ξ j ( x, ω ) . We want to apply the residue theoremto get an asymptotic expansion in the time domain. Note that: ˆ − (cid:15) − ρ Ξ j ( x, ω ) e − iωt d ω + ˆ ρ(cid:15) Ξ j ( x, ω ) e − iωt d ω = ˛ C ± Ξ j ( x, Ω) e − i Ω t dΩ − ˆ C ± ρ Ξ j ( x, Ω) e − i Ω t dΩ − ˆ C ± (cid:15) Ξ j ( x, Ω) e − i Ω t dΩ , where the integration contours C ± ρ and C ± (cid:15) are semi-circular arcs of radius ρ and (cid:15) , respectively, inthe upper (+) or lower (-) half-planes, and C ± is the closed contour defined as C ± := C ± ρ ∪ C ± (cid:15) ∪ [ − ρ, − (cid:15) ] ∪ [ (cid:15), ρ ] . The integral on the closed contour is the main contribution to the scattered fieldby the mode and can be computed using the residue theorem to get, for ρ ≥ max j ∈ N (cid:60) [Ω ± j ( δ )] and < (cid:15) ≤ min j ∈ N (cid:60) [Ω ± j ( δ )] , ˛ C + Ξ j ( x, Ω) e − i Ω t dΩ = 0 , ˛ C − Ξ j ( x, Ω) e − i Ω t dΩ = 2 πi Res (cid:0) Ξ j ( x, Ω) e − i Ω t , Ω ± j ( δ ) (cid:1) . Since Ω ± j ( δ ) is a simple pole of ω (cid:55)→ λ ( ω ) − λ j ( ωδ ) we can write: ˛ C − Ξ j ( x, Ω) e − i Ω t dΩ = 2 πi Res (cid:0) Ξ j ( x, Ω) , Ω ± j ( δ ) (cid:1) e − i Ω ± j ( δ ) t . To compute the integrals on the semi-circle, we introduce: X j ( y, Ω) = ϕ j ( y ) (cid:112) | x − y | √ ω ( λ ( ω ) − λ j ) y ∈ ∂D. Given the regularity of the input signal (cid:98) f ∈ C ∞ ([0 , C ]) , the Paley-Wiener theorem [42, p.161]ensures decay properties of its Fourier transform at infinity. For all N ∈ N ∗ there exists a positiveconstant C N such that for all Ω ∈ C | f (Ω) | ≤ C N (1 + | Ω | ) − N e C |(cid:61) (Ω) | . Let T := ( | x − y | + d · z ) /c . We now rewrite the integrals on the large semi-circle ˆ C ± ρ Ξ j ( x, Ω) e − i Ω t dΩ = ˆ C ± ρ f (Ω) ˆ ∂D X j ( y, Ω) e i Ω( T − t ) d σ ( y )dΩ . We have that t − + C ≤ T ≤ t +0 − C . Two cases arise. Case 1:
For < t < t − , i.e., when the signal emitted at s has not reached the observation point x ,we choose the upper-half integration contour C + . Transforming into polar coordinates, Ω = ρe iθ for θ ∈ [0 , π ] , we get: (cid:12)(cid:12)(cid:12) e i Ω( T − t ) (cid:12)(cid:12)(cid:12) ≤ e − ( t − − t + C ) (cid:61) (Ω) ∀ y ∈ ∂D, (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˆ C + ρ Ξ j ( x, Ω) e − i Ω t dΩ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ ˆ π ρ (cid:12)(cid:12) f (cid:0) ρe iθ (cid:1)(cid:12)(cid:12) e − ρ ( t − − t + C ) sin θ ˆ ∂D | X j ( y, ρe iθ ) | d σ ( y )d θ, ≤ ρC N (1 + ρ ) − N δ max θ ∈ [0 ,π ] (cid:12)(cid:12) X j (cid:0) · , ρe iθ (cid:1)(cid:12)(cid:12) L ∞ ( ∂D ) π − e − ρ ( t − − t ) ρ ( t − − t ) , where we used that for θ ∈ [0 , π/ , we have sin θ ≥ θ/π ≥ and − cos θ ≤ − θ/π . The usualway to go forward from here is to take the limit ρ → ∞ , and get that the limit of the integral on thesemi-circle is zero. As in the three-dimensional case, we work in the quasi-static approximation here,and our modal expansion is not uniformly valid for all frequencies. So we have to work with a fixedmaximum frequency ρ . However, the maximum frequency ρ depends on the size of the particle via thehypothesis ρ ≤ cδ − . Since N can be taken arbitrarily large and that X j behaves like a polynomial in ρ whose degree does not depend on j , we get that, uniformly for j ∈ [1 , J ] : (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˆ C + ρ Ξ j ( x, Ω) e − i Ω t dΩ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = O (cid:18) δt − − t ρ − N (cid:19) . For the upper-half semi-circle of radius (cid:15) , we also transform into polar coordinates with the changeof variable
Ω = (cid:15)e iθ , for θ ∈ [0 , π ] , and get: (cid:12)(cid:12)(cid:12)(cid:12) ˆ C + (cid:15) Ξ j ( x, Ω) e − i Ω t dΩ (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:15)C N (1 + ρ ) − N δ max θ ∈ [0 ,π ] (cid:12)(cid:12) X j (cid:0) · , (cid:15)e iθ (cid:1)(cid:12)(cid:12) L ∞ ( ∂D ) π − e − (cid:15) [( t − − t )] (cid:15) ( t − − t ) , Case 2:
For t > t +0 , we choose the lower-half integration contour C − . Transforming into polarcoordinates, Ω = ρe iθ for θ ∈ [ π, π ] , we get (cid:12)(cid:12)(cid:12) e i Ω( T − t ) (cid:12)(cid:12)(cid:12) ≤ e ( t − t +0 − C ) (cid:61) (Ω) ∀ y ∈ ∂D , and (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˆ C − ρ Ξ j ( x, Ω) e − i Ω t dΩ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ ˆ ππ ρ (cid:12)(cid:12) f (cid:0) ρe iθ (cid:1)(cid:12)(cid:12) e ρ ( t − t +0 ) sin θ ˆ ∂D | X j ( y, ρe iθ ) | d σ ( y )d θ, ≤ ρC N (1 + ρ ) − N δ max θ ∈ [ π, π ] (cid:12)(cid:12) X j (cid:0) · , ρe iθ (cid:1)(cid:12)(cid:12) L ∞ ( ∂D ) π − e − ρ ( t − t +0 ) ρ ( t − t +0 ) , Exactly as in Case , we cannot take the limit ρ → ∞ . However, the maximum frequency ρ dependson the size of the particle via the hypothesis ρ ≤ cδ − . Using the fact that N can be taken arbitrarilylarge and that X j behaves like a polynomial in ρ whose degree does not depend on j , we get that,uniformly for j ∈ [1 , J ] : (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˆ C − ρ Ξ j ( x, Ω) e − i Ω t dΩ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = O (cid:18) δt ρ − N (cid:19) . For the lower-half semi-circle of radius (cid:15) , we also transform into polar coordinates with the changeof variable
Ω = (cid:15)e iθ , for θ ∈ [0 , π ] , and get: (cid:12)(cid:12)(cid:12)(cid:12) ˆ C − (cid:15) Ξ j ( x, Ω) e − i Ω t dΩ (cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:15)C N (1 + ρ ) − N δ max θ ∈ [ π, π ] (cid:12)(cid:12) X j (cid:0) · , (cid:15)e iθ (cid:1)(cid:12)(cid:12) L ∞ ( ∂D ) π − e − (cid:15) ( t − t +0 ) (cid:15) ( t − t +0 ) , The result of theorem 5.3 is obtained by summing the contribution of all the modes.35 eferences [1]
H. Ammari, Y. Deng, and P. Millien , Surface plasmon resonance of nanoparticles and ap-plications in imaging , Archive for Rational Mechanics and Analysis, 220 (2016), pp. 109–153.[2]
H. Ammari, B. Fitzpatrick, H. Kang, M. Ruiz, S. Yu, and H. Zhang , Mathematicaland Computational Methods in Photonics and Phononics , vol. 235, Mathematical Surveys andMonographs, 2018.[3]
H. Ammari, J. Garnier, W. Jing, H. Kang, M. Lim, K. Sølna, and H. Wang , Mathe-matical and statistical methods for multistatic imaging , vol. 2098, Springer, 2013.[4]
H. Ammari and H. Kang , Polarization and Moment Tensors With Applications to InverseProblems and Effective Medium Theory , vol. 162, Springer-Verlag New York, 2007.[5]
H. Ammari, P. Millien, M. Ruiz, and H. Zhang , Mathematical analysis of plasmonicnanoparticles: The scalar case , Arch. Ration. Mech. Anal., 224 (2017), pp. 597–658.[6]
H. Ammari, M. Putinar, M. Ruiz, S. Yu, and H. Zhang , Shape reconstruction of nanoparti-cles from their associated plasmonic resonances , Journal de Mathématiques Pures et Appliquées,122 (2019), pp. 23–48.[7]
H. Ammari, M. Ruiz, S. Yu, and H. Zhang , Mathematical analysis of plasmonic resonancesfor nanoparticles: The full Maxwell equations , Journal of Differential Equations, 261 (2016),pp. 3615–3669.[8] ,
Reconstructing fine details of small objects by using plasmonic spectroscopic data , SIAMJournal on Imaging Sciences, 11 (2018), pp. 1–23.[9]
K. Ando and H. Kang , Analysis of plasmon resonance on smooth domains using spectral prop-erties of the Neumann-Poincaré operator , Journal of Mathematical Analysis and Applications,435 (2016), pp. 162–178.[10]
K. Ando, H. Kang, and Y. Miyanishi , Exponential decay estimates of the eigenvalues forthe neumann-poincare operator on analytic boundaries in two dimensions , J. Integral EquationsApplications, 30 (2018), pp. 473–489.[11]
K. Ando, H. Kang, Y. Miyanishi, and T. Nakazawa , Surface localization of plasmons inthree dimensions and convexity , (2020).[12]
K. Ando, H. Kang, Y. Miyanishi, and M. Putinar , Spectral analysis of Neumann-Poincaréoperator , (2020).[13]
P. Y. Chen, D. J. Bergman, and Y. Sivan , Generalizing normal mode expansion of electro-magnetic green’s tensor to open systems , Phys. Rev. Applied, 11 (2019), p. 044018.[14]
K. Cognée, W. Yan, F. La China, D. Balestri, F. Intonti, M. Gurioli, A. Koenderink,and P. Lalanne , Mapping complex mode volumes with cavity perturbation theory , Optica, 6(2019), pp. 269–273.[15]
R. Colom, R. McPhedran, B. Stout, and N. Bonod , Modal expansion of the scattered field:Causality, nondivergence, and nonresonant contribution , Phys. Rev. B, 98 (2018), p. 085418.[16]
D. Colton and R. Kress , Integral equation methods in scattering theory , vol. 72, SIAM, 2013.[17]
J.-P. Demailly , Complex analytic and differential geometry , Citeseer, 1997.3618]
M. B. Doost, W. Langbein, and E. A. Muljarov , Resonant-state expansion applied tothree-dimensional open optical systems , Phys. Rev. A, 90 (2014), p. 013834.[19]
R.-C. Ge and S. Hughes , Design of an efficient single photon source from a metallic nanoroddimer: a quasi-normal mode finite-difference time-domain approach , Opt. Lett., 39 (2014),pp. 4235–4238.[20]
Y.-G. Ji and H. Kang , A concavity condition for existence of a negative value in Neumann-Poincaré spectrum in three dimensions , Proceedings of the American Mathematical Society, 147(2019), pp. 3431–3438.[21]
D. Khavinson, M. Putinar, and H. S. Shapiro , Poincaré’s variational problem in potentialtheory , Archive for rational mechanics and analysis, 185 (2007), pp. 143–184.[22]
K. D. Kokkotas and B. G. Schmidt , Quasi-normal modes of stars and black holes , LivingReviews in Relativity, 2 (1999), pp. 1433–8351.[23]
P. T. Kristensen, R.-C. Ge, and S. Hughes , Normalization of quasinormal modes in leakyoptical cavities and plasmonic resonators , Phys. Rev. A, 92 (2015), p. 053810.[24]
J. L. Lagrange , Traité de la résolution des équations numériques , Paris, 1798.[25]
P. Lalanne, W. Yan, K. Vynck, C. Sauvan, and J.-P. Hugonin , Light interaction withphotonic and plasmonic resonances , Laser Photonics Reviews, 12 (2018).[26]
P. T. Leung and K. M. Pang , Completeness and time-independent perturbation of morphology-dependent resonances in dielectric spheres , J. Opt. Soc. Am. B, 13 (1996), pp. 805–817.[27]
Y. Miyanishi , Weyl’s law for the eigenvalues of the neumann–poincar \ ’e operators in three di-mensions: Willmore energy and surface geometry , arXiv preprint arXiv:1806.03657, (2018).[28] Y. Miyanishi and G. Rozenblum , Eigenvalues of the neumann–poincaré operator in dimension3: Weyl’s law and geometry , St. Petersburg Mathematical Journal, 31 (2020), pp. 371–386.[29]
A. Moiola and E. A. Spence , Acoustic transmission problems: wavenumber-explicit boundsand resonance-free regions , Mathematical Models and Methods in Applied Sciences, 29 (2019),pp. 317–354.[30]
R. M. More , Theory of decaying states , Phys. Rev. A, 4 (1971), pp. 1782–1790.[31]
E. A. Muljarov, W. Langbein, and R. Zimmermann , Brillouin-wigner perturbation theoryin open electromagnetic systems , EPL (Europhysics Letters), 92 (2010), p. 50010.[32]
J.-C. Nédélec , Acoustic and electromagnetic equations: integral representations for harmonicproblems , vol. 144, Springer Science & Business Media, 2001.[33]
M. A. Ordal, L. L. Long, R. J. Bell, S. E. Bell, R. R. Bell, R. W. Alexander, andC. A. Ward , Optical properties of the metals al, co, cu, au, fe, pb, ni, pd, pt, ag, ti, and w inthe infrared and far infrared , Appl. Opt., 22 (1983), pp. 1099–1119.[34]
A. Pick, B. Zhen, O. D. Miller, C. W. Hsu, F. Hernandez, A. W. Rodriguez, M. Sol-jačić, and S. G. Johnson , General theory of spontaneous emission near exceptional points ,Opt. Express, 25 (2017), pp. 12325–12348.[35]
D. A. Powell , Resonant dynamics of arbitrarily shaped meta-atoms , Phys. Rev. B, 90 (2014),p. 075108.[36]
Y. K. Sirenko, S. Ström, and N. P. Yashina , Modeling and analysis of transient processesin open resonant structures: New methods and techniques , vol. 122, Springer, 2007.3737]
B. Stout, R. Colom, N. Bonod, and R. McPhedran , Eigenstate normalization for openand dispersive systems , arXiv preprint arXiv:1903.07183, (2019).[38]
B. Stout and R. McPhedran , Egocentric physics: Just about Mie , EPL (Europhysics Letters),119 (2017), p. 44002.[39]
B. Vial, F. Zolla, A. Nicolet, and M. Commandré , Quasimodal expansion of electromag-netic fields in open two-dimensional structures , Phys. Rev. A, 89 (2014), p. 023829.[40]
H. Wang , Shape identification in electro-sensing . https://github.com/yanncalec/SIES , 2013.[41] W. Yan, R. Faggiani, and P. Lalanne , Rigorous modal analysis of plasmonic nanoresonators ,Phys. Rev. B, 97 (2018), p. 205422.[42]
K. Yosida , Functional Analysis , Classics in Mathematics, Springer Berlin Heidelberg, 6 ed., 1995.[43]
X. Zambrana-Puyalto and N. Bonod , Purcell factor of spherical mie resonators , Phys. Rev.B, 91 (2015), p. 195422.[44]
S. Zaremba , Les fonctions fondamentales de M. Poincaré et la méthode de Neumann pour unefrontière composée de polygones curvilignes , Journal de Mathématiques Pures et Appliquées, 10(1904), pp. 395–444.[45]