Modal approximation for plasmonic resonators in the time domain: the scalar case

We study the electromagnetic field scattered by a metallic nanoparticle with dispersive material parameters in a resonant regime. We consider the particle placed in a homogeneous medium in a low-frequency regime. We define modes for the non-Hermitian problem as perturbations of electro-static modes, and obtain a modal approximation of the scattered field in the frequency domain. The poles of the expansion correspond to the eigenvalues of a singular boundary integral operator 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 in the time domain. We present numerical simulations in two dimensions to corroborate our results.


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 dispersive obstacle in the transverse magnetic polarisation (see [32, remark 2.1]). We work in a low-frequency regime corresponding to relevant physical applications, such as the scattering of light in the visible/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 field by 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 this paper, we do not deal with the high frequency part of the field that is usually studied with micro-local analysis tools.

Previous work on plasmonic resonances and layer potentials
It has been shown in [1,6,8] 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 boundary of the particle. In a low-frequency regime, i.e., at frequencies corresponding to wavelengths that are orders of magnitude larger than the particle's size, asymptotic analysis techniques, as in [6], yield a hierarchy of boundary integral equations. The asymptotic small parameter is δωc −1 , where δ is the size of the particle, ω the frequency and c the velocity. At leading order the well-known Neumann-Poincaré operator appears [48]. Using the Plemelj symmetrisation principle and the spectral theory of compact selfadjoint operators, the latter can be diagonalised in the appropriate functional spaces [24,35], which allows the scattered field to be decomposed in a basis of orthogonal modes in the static case [10]. The properties of the eigenvalues of the Neumann-Poincaré operator have been extensively studied in the literature, see the review paper [13] and references therein. For a smooth enough boundary, say C 1,α for some α > 0, the operator is compact and its eigenvalues are real numbers converging to zero. The eigenvalues of the Neumann-Poincaré operator in the two-and three-dimensional cases are 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 rate of the eigenvalues depends strongly on the regularity of the boundary. For an analytic boundary, the eigenvalues have an exponential decay rate [11]. In three dimensions, very few surfaces are known to have negative eigenvalues [23]. For a strictly convex C ∞ domain, there are infinitely many positive eigenvalues and a finite number of negatives ones [12]. The eigenvalues rate of decay is much slower than in two dimensions: λ j = O( j −1/2 ) as j → ∞ [31].

Contributions and organisation of the paper
We begin by describing the problem geometry and we formulate the governing equations in Sect. 2. We introduce the layer potential and boundary integral formulation and recall the modal decomposition of the static (ω = 0) solution. In Sect. 3, we prove that in three dimensions, for a strictly convex particle, the modal expansion can be truncated due to the super-polynomial decay of the expansion's coefficients. With a perturbation argument, we deduce from the static (ω = 0) result a modal approximation in the dynamic case (for a small non-zero frequency). The perturbation analysis yields size dependent dynamic complex resonant frequencies. We show that all the resonant frequencies have a negative imaginary part and lie in a bounded region near the origin. Finally, in Sect. 5, using only elementary complex analysis techniques, we give an approximation for the low-frequency part of the scattered field in the time domain as a finite sum of modes oscillating at complex resonant frequencies. We also show with a simple causality argument that the exponential catastrophe is not problematic in practice. In Sect. 6, we implement this expansion in the two-dimensional setting and illustrate the validity of our approach with numerical simulations.

Problem setting
We are interested in the scattering problem of an incident wave illuminating a plasmonic nanoparticle in R d , d = 2, 3. The homogeneous medium is characterised by electric permittivity ε m and magnetic permeability μ m . Let D be a smooth bounded domain in R d , of class C ∞ , characterised by electric permittivity ε c . We assume the particle to be non-magnetic, i.e., μ c = μ m . Let D = z + δ B where B is the reference domain and contains the origin, and D is located at z ∈ R d and has a characteristic size δ 1. 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 0 the speed of light in vacuum c 0 = 1/ √ ε 0 μ 0 and by c the speed of light in the medium c = 1/ √ ε m μ m .
Hereafter we use the Drude model [36] to express the electric permittivity of the particle: where the positive constants ω p and T −1 are the plasma frequency and the collision frequency or damping 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 the zero level set of a bivariate polynomial. A quadrature domain is the conformal image of the unit disc by a rational function.
Remark 1 Algebraic domains are dense among all planar domains, so every smooth curve can be described as a sequence of algebraic curves [7].

Condition 2
In three dimensions, we assume the domain D to be strictly convex: for any two points in 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.

Helmholtz equation for a subwavelength resonator
Given an incident wave u in solution to the Helmholtz equation, the scattering problem in the frequency domain can be modelled by subject to the Sommerfeld radiation condition uniformly in x/|x|, for k m > 0. The transmission conditions are given by Here, ∂ · /∂ν denotes the normal derivative on ∂ D, and the + and − subscripts indicate the limits from outside and inside D, respectively.

Definition 1
We denote the contrast λ by Definition 2 (Resonant frequency, mode) We say ω is a resonant frequency if there is a nontrivial solution to Eq. (2) with u in = 0. We call the solution a mode. A subwavelength resonance occurs when a resonant frequency ω satisfies ωδc −1 < 1.

Remark 2
Although the Helmholtz equation does not model light matter interaction in three dimensions, the spectral properties of the Maxwell volume integral operator are closely linked to those of the Neumann-Poincaré operator [5,19]. Therefore it is insightful to study the three-dimensional case. Moreover, the three-dimensional Helmholtz equation models acoustic waves and after adapting the physical parameters one could use this setting to study bubbles.

Layer potential formulation
Let H 1/2 (∂ D) be the usual Sobolev space and let H −1/2 (∂ D) be its dual space with respect to the duality pairing ·, · 1 2 ,− 1 2 . The field u can be represented using the single-layer potentials S k c D and S k m D , introduced in Definition 8, as follows: where the pair and where K k m , * D is the Neumann-Poincaré operator introduced in Definition 8. The trace relations for the single-layer potential are given in Lemma 6.

Scaling and small-volume approximation
The goal of this section is to establish an equivalent formulation for (4) in the form [Ψ ] = F (Proposition 1), in order to write an asymptotic expansion of the operator A ωδc −1 B (Lemma 1) and a spectral decomposition for the limiting operator A 0 B (Proposition 2). The scaling is new in this context, but the asymptotic expansion and the spectral decomposition were first obtained in [6]. We recall them here for the sake of completeness. The proofs are quite lengthy and technical, so they are included 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 Ξ(X ) := Ξ(z + δ X ), X ∈ ∂ B. The scaling properties of the integral operators are given in Appendix B. The solution u becomes where the single-layer potential S kδ B and Neumann-Poincaré operator K kδ, * B are defined by the fundamental solution Γ kδ . The density pair and Since S k c δ is invertible for k c δ small enough (see Lemmas 7 and 10 ), the following proposition holds. Proposition 1 For d = 2, 3, the following equation holds for Ψ : where Lemma 1 (Small-volume expansion) As ωδc −1 → 0, A ωδc −1 B admits the following asymptotic expansion: and Proof See Appendix C.2.

The operator
is not self-adjoint in L 2 so it can not be diagonalised directly to solve (6). However, in the static regime, the operator A 0 B can be expressed simply with K * B , which can be symmetrised in the Hilbert space H * (∂ B) (see Appendix A.2). Lemma 2 (Spectral decomposition of K * B ) K * B is self-adjoint with respect to the inner product ·, · H * (∂ B) . Moreover, it is compact, so its spectrum is discrete. The spectral theorem yields the decomposition where {λ j } j∈N are the eigenvalues of K * B and { φ j } j∈N their associated normalised eigenvectors.
where (λ j , φ j ) j∈N are the eigenvalues and normalised eigenfunctions of K * B in H * (∂ B) and Proof Direct consequence of Lemma 2 and (9).

Corollary 1
The spectral approximation of the static (ω = 0) solution is given by where F is defined in Proposition 1.

Modal decomposition of the field
In this section, we want to apply perturbation theory tools to express the solution of (6) in terms of the eigenvectors of K * B that appear in the spectral decomposition of the limiting problem in Proposition 2, and to replace τ j by a perturbed value τ j (ωδc −1 ). Classical perturbation theory will give us a Taylor expansion for τ j (ωδc −1 ) in ωδc −1 for any j ∈ N but the remainders and validity range of these expansions will depend on the index j of the considered eigenvalue. In order to get a meaningful expansion of the scattered field we need to work with a finite number of modes.

Modal expansion truncation
In practice, there is no need to consider the whole spectral decomposition of the field. It has been empirically reported that only a few modes actually contribute to the scattered field. The number of modes to consider increases as the source gets closer to the particle. In this section we give a mathematical explanation of this phenomenon : the modes φ j are eigenmodes of a pseudo-differential operator of order −1, and are oscillating functions. As in classical Fourier analysis, the decay with j of the coefficients F, φ j H * (∂ B) will be determined by the regularity of the function F and the number of modes to consider will depend on the spatial variations of F over ∂ B. In an homogeneous medium the incoming field is smooth and therefore we can expect a fast decay of the coefficients.
Remark 3 (On high-order modes) The oscillation rate of high-order modes has been quantified in the case of smooth axisymmetric inclusions [39]. Moreover, it has been also shown in that case that the contribution of high-order modes to the electric field outside the particle decays exponentially as the distance from the inclusion increases. This property has also been shown in [12]. This suggests that not only high-order modes are not really excited by a source bounded away from the particle because of their oscillating properties, but they also do not contribute to the scattered field except in the near field.

Remark 4
In this work we do not use the exponential decay of the modes far away from the particle to justify the truncation of the quasi-static modal expansion because the exponential decay has only been shown for the electro-static modes, i.e., the quantities S 0 D [φ j ] as j → ∞, while our expansions uses dynamic modes S ω D [φ j ]. Nevertheless we think that the property still holds for ω = 0.

The three-dimensional case
Let H S be the standard Sobolev space of order S. Proposition 3 For B, a strictly convex domain in R 3 with C ∞ -smooth boundary, and F ∈ H S (∂ B) for some S ∈ N * we have : The proof relies on a theorem from [12] which itself uses the computation of the principal symbol of the Neumann-Poincaré operator done in [30]: Theorem 1 (From [12], p. 7) For B, a strictly convex domain in R 3 with C ∞ -smooth boundary, K * B has a finite number of non-positive eigenvalues. We can modify K * B by adding a finite dimensional smoothing operator to have a positive definite elliptic pseudo-differential operator of order − 1, which we denote by K * B . For each real number s ∈ R there exist constants c s , C s ∈ R + such that is self-adjoint and has the same eigenvalues as K * B . Its eigenvectors are . It can be modified by adding a finite dimensional smoothing operator to have a positive definite elliptic pseudo-differential operator of order −1, which we denote by K * B . For each real number s ∈ R there exist constants c s , C s ∈ R + such that for all φ ∈ H s−1/2 (∂ B). Moreover there exists j 0 ∈ N such that Proof K * B has the same principal symbol as K * B [31, p. 8].
We will also need the decay estimate of the eigenvalues of K * B : Theorem 2 (From [31]) For B, a strictly convex domain in R 3 with C ∞ -smooth boundary the eigenvalues of the Neumann-Poincaré operator satisfy: with C B a constant depending only on B: where W (∂ B) and χ(∂ B) denote, respectively, the Willmore energy and the Euler characteristic of the boundary surface ∂ B.
where Ker ( K * B ) denotes the kernel of K * B . The symbol ⊕ is to be understood in the L 2 scalar product sense. Hence for j ≥ j 0 : .
where we used the fact that (−S B ) Then .
Since the eigenvectors of K * B are orthogonal in L 2 (∂ B) we have: .
We can now write .
Iterating this procedure S − 1 times yields . Hence We need to control the L 2 -norm of G (S) . We can rewrite the orthogonal decomposition as ker . Composing by K * B we get: Using the right-hand side of (12) with s = S + 1 2 we get .
Using the Cauchy-Schwartz inequality in (13) and the fact that ψ j L 2 (∂ B) = 1 (S B is an isometry): is independent of j. Using Theorem 2 we can see that for j large enough since λ j ∼ C B j −1/2 we have: ,

The two-dimensional case
In two dimensions, the picture is slightly different. The eigenspace associated to eigenvalue zero can have infinite dimension (it is the case for the disk) and there are infinitely many negative eigenvalues. As a result, K * D can not be modified into a positive operator by adding a finite dimensional operator. However, for a certain class of domains, it is possible to show that there is a finite number of plasmonic resonances. For example, it was shown in [7] that an algebraic domain of class Q has asymptotically a finite number of plasmonic resonances. The asymptotic parameter is the deformation from the unit circle. For a larger class of domains the decay of the coefficients F, φ j H * (∂ D) can be checked numerically (see Sect. 6).

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 to be C ∞ , so the source term in Eq. (6), i.e., the function F, is smooth on ∂ B. Therefore using Proposition 3 we have a super-polynomial decay of the coefficients F, φ j H * (∂ B) , and we can consider that only a finite number of modes are excited. The number J of modes to consider depends on the incoming field.

Proposition 4 Assume that F
The spectral approximation of the scattered field as ωδc −1 → 0 is given by and F is defined in Proposition 1.

Proof Note that { φ j } j∈N forms an orthonormal basis of H * (∂ B). Writing
Using (3) and (8) concludes the proof.
For each normalised eigenfunction of K * B , we consider the corresponding function on ∂ D, Here {φ j } j∈N are the rescaled non-normalised eigenfunctions of K * D . Let us introduce .
Going back to the original unscaled problem:

Proposition 5
As ωδc −1 1, the spectral decomposition of the field is as follows Proof The scaling Lemma 11

Size dependent resonant frequencies
In this section we calculate size dependent plasmonic resonances. Let j ∈ {0, . . . , J }. Recall that

Definition 3
We say that ω is a static plasmonic resonance if τ j = 0.

Definition 4
We say that ω is first-order corrected plasmonic resonance if

Remark 5
For j = 0, we have τ 0 = 1/ε m , which is of size one by assumption. We exclude j = 0 from the set of resonances.
Then, we can calculate

Lemma 3
We have α j ∈ R and In what follows we use the lower-case character ω for real frequencies and the upper-case character Ω for complex frequencies.
Proposition 6 Using the Drude model (1), the three-dimensional first-order corrected plasmonic resonances Ω ± j (δ) := ±Ω j + iΩ j all lie in the lower part of the complex plane and their modulus is bounded. In the case where we take the medium to be vacuum, i.e., ε m = ε 0 we obtain explicitly for |λ j + 1/2| > 10 −2 (this occurs, for example, when B is a ball [2]): Moreover, they are bounded Because δω p c −1 1, we get the desired result. Lagrange improved upper-bound for roots of polynomials concludes the proof [27].

Definition 5
In three dimensions, we define the resonance radius as Remark 6 This resonance radius gives our method a range of validity. We compute resonant frequencies in a perturbative quasi-static regime. So by checking that 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 −1 is bigger than one, it means that the method is not self-consistent, as the largest resonant frequency might not satisfy the ωδc −1 < 1/2.

Proposition 7
In vacuum, and using the Drude model (1), the two-dimensional first-order corrected plasmonic resonances are the roots Ω j 1≤ j≤J ∈ C of the following equation

Remark 7
We can compute an approximation of the roots of (15) by computing in the first place the static resonances Ω s, j 1≤ j≤J . Solving τ j = 0 yields Replacing the dynamic frequency in the logarithm by its static approximation, we transform (15) into the quadratic equation We get

Definition 6
In two dimensions, we define the resonance radius as

Plasmonic quasi-normal modes
Quasi-normal modes are formally defined as solutions of the source-free wave equation [28]. Using the representation formula (3), we can now define, as in the physics literature, plasmonic quasi-normal modes (e ± j ) j∈N that oscillate at complex frequencies Ω ± j (δ): 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 8
In the physics literature (see [28, equation (1.1)] for instance) one can often find representations of the scattered field in the form where α j are excitation coefficients depending on the source and independent of the space variable x. The first problem of these representations is that the eigenmodes diverge exponentially and a renormalisation process is necessary. Numerical instability will appear far away from the particle. Even though the study of these modes individually can give physical insight to a system (like for example by studying the mode volume quantity [16]), they cannot be used in frequency domain representation formulae to solve the scattering problem. Moreover, this type of representation formulae is a generalisation of the well-known eigenmode expansions for waves in bounded domains with classical Dirichlet or Neumann conditions. For these formulae to be valid one needs to address the question of completeness of the eigenbasis. Here we have shown the completeness of the charge distribution modes on the boundary of the particle only. Moreover, our perturbative analysis only holds for a finite number of modes and cannot be used to show completeness of the perturbed modes as soon as ω = 0.

Time domain approximation of the scattered field
In the following section we show that even though they are irrelevant for frequency domain representation, quasi-normal modes can be used to approximate the field in the time domain.
The idea is to get around costly time domain computations by pre-computing the modes of the system and then expressing the response of the system to any source in terms of the modes. In the physics literature (for example [28, eq. (1. 2)]) the field in the time domain is expressed under the form The problem with this type of expansions is that if |x| is big then e ± j (x) is exponentially large and the computation 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 domain in a similar expansion, but with non-diverging, pre-computable quantities similar to the quasi-normal modes.

The three-dimensional case
Here we state the main result of the paper, Theorem 3, and discuss the result.

The modal approximation
Let Γ k m (·, s), i.e., the Green's function for the Helmholtz equation introduced in Definition 7, be the incident wave u in in three dimensions. Given a wideband signal , for C 1 > 0, we want to express the time domain response of the electric field to an oscillating dipole placed at a source point s. This means that for a fixed δ we can pick an excitation signal such that most of the frequency content is in the low frequencies but large enough to excite the plasmonic resonances. We can pick η 1 and ρ ≥ R(δ) such that The incident field has the following form in the time domain: The goal of this section is to establish a resonance expansion for the low-frequency part of the scattered electric field in the time domain. Introduce, for ρ > 0, the truncated inverse Fourier transform of the scattered field u sca given by Recall that z is the centre of the resonator and δ its radius. Let us define the time it takes to the wideband signal to reach first the scatterer and then the observation point x. The term ±2δ/c accounts for the maximal timespan spent inside the particle.
Recall the spectral decomposition in the frequency domain (Proposition 5) for x ∈ R 3 \D: Theorem 3 Let N ∈ N. The scattered field has the following form in the time domain for x ∈ R 3 \D: .

Remark 9
The resonant frequencies Ω ± j (δ) 1≤ j≤J have negative imaginary parts, so Theorem 3 expresses 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 10 (About the remainder ρ) Since for a particle of finite size δ our expansion only holds for a range of frequencies ω such that ωδc −1 < 1, we cannot compute the full inverse Fourier transform and we have a remainder that depends on the maximum frequency that we can use. Nevertheless that maximum frequency ρ behaves as c/δ and we can see that the remainder gets arbitrarily small for small particles. For a completely point-like particle one would get a zero remainder.

Remark 11
If we had access to the full inverse Fourier transform of the field, of course, since the inverse Fourier transform of a function which is analytic in the upper-half plane is causal we would find that in the case t ≤ (|s − z| + |x − z| − 2δ) /c, u sca (x, t) = 0. Nevertheless, our method gives the resonant 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 good approximation 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 ω → ∞, ε(ω) −→ 1. The metal does not really interact with light at high frequencies.

Alternative formulation with non-diverging causal quasi-normal modes
Even though |e ± j (x)| −→ ∞ when |x| → ∞, no terms diverge in (19). Indeed we can rewrite: where 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 (δ): Remark 12 When referring to E ± j , the term mode is inaccurate, as E ± j does not solve the Helmholtz equation. But since the (E ± j ) j∈N are built from modes with a complex phase correction, we still call them modes in a loose sense of the term.
Theorem 3 can be re-stated: (21) has exactly the same form as the representation formula found in the physics literature (like Eq. (18)) but without any exponentially diverging quantities. The E ± j can be computed independently of the source, just like regular quasi-normal modes.

Proof of Theorem 3
Before we can prove Theorem 3 we need the following lemma: (7) admits the following asymptotic expansion: Proof See Appendix C.3.
Proof (Proof of Theorem 3) We start by studying the time domain response of a single mode to a causal excitation at the source point s. According to Proposition 5 we need to compute the contribution Ξ j of each mode, that is, where λ j (ωδ) := λ j − ωδc −1 2 α j + O ωδc −1 3 . One can then write: . Since ν, ϕ 0 H * (∂ D) = 0, the zeroth term vanishes in the summation. Now we want to apply the residue theorem to get an asymptotic expansion in the time domain. Note that: where the integration contour C ± ρ is a semicircular arc of radius ρ in the upper (+) or lower (−) half-plane, and C ± is the closed contour The integral on the closed contour is the main contribution to the scattered field by the mode and can be computed using the residue theorem to get, for ρ ≥ [Ω ± j (δ)], we can write: To compute the integrals on the semi-circle, we introduce: Note that B j (·, ·, Ω) behaves like a polynomial in Ω when |Ω| → ∞. Given the regularity of the input signal f ∈ C ∞ 0 ([0, C 1 ]), the Paley-Wiener theorem [46, p. 161] ensures decay properties of its Fourier transform at infinity. For all N ∈ N * there exists a positive constant Let T := (|x − y| + |s − v|)/c. We now re-write the integrals on the semi-circle: We have that t − 0 + C 1 ≤ T ≤ t + 0 − C 1 . Two cases arise. Case 1: For 0 < t < t − 0 , 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: where we used that for θ ∈ [0, π/2], we have sin θ ≥ 2θ/π ≥ 0 and − cos θ ≤ −1 + 2θ/π. The usual way to go forward from here is to take the limit ρ → ∞, and get that the limit of the integral on the semi-circle is zero. However, 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 fixed maximum frequency ρ. Since N can be taken arbitrarily large and that B j behaves like a polynomial in ρ whose degree does not depend on j, we get that, uniformly for j ∈ [1, J ]: Of course if one has to consider the full inverse Fourier transform of the scattered electromagnetic field, by causality, one should expect the limit to be zero. However, one would need high-frequency estimates of the electromagnetic field, as well as a modal decomposition that is uniformly valid for all frequencies. Since our modal expansion is only valid for a limited range of frequencies we get an error bound that is arbitrarily small if the particle is arbitrarily small, but not strictly zero.
Exactly as in Case 1, we cannot take the limit ρ → ∞. Using the fact that N can be taken arbitrarily large and that B j behaves like a polynomial in ρ whose degree does not depend on j, we get that, uniformly for j ∈ [1, J ]: The result of Theorem 3 is obtained by summing the contribution of all the modes considered.

Remark 14
The fact that we work with a finite number of modes is necessary for the perturbation theory of Sect. 3 but also in this section. Indeed, if we consider all the modes there is an accumulation point in the poles of the modal expansion of the field, and therefore we cannot apply the residue theorem.

The two-dimensional case
In two dimensions, the Green's function does not have an explicit phase term, so we need to introduce another asymptotic parameter > 0 to be able to use the large argument asymptotics of the Hankel function. Our new truncated inverse Fourier transform of the scattered field u sca given by This allows us to define a notion of far field. A point x is far from D if |x − z|c −1 1. We can now add two additional hypotheses: -the source is far away from the particle (or equivalently, the incoming wave is a plane wave) -the observation point is far away from the particle.
The incident field has the following form in the time domain: Besides these two assumptions and a difference in the order of the remainder, the result in two dimensions is essentially the same as in three dimensions.

Theorem 5
Let N ∈ N. The scattered field has the following form in the time domain for x far away from D: with Ω ± j (δ) being the plasmonic resonant frequencies of the particle given by Proposition 7. C Ω ± j (δ) is a constant depending only on j, the size δ and the model for ε c (ω): .

Proof
The proof is quite similar to the three-dimensional case. It is included in Appendix D for the sake of completeness.

Numerical simulations
The goal of this section is to illustrate the validity of our approach and to show that the approximation seems to be working with less restrictive hypotheses than the ones in Theorem 5: -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 layer potentials developed in [44].

Domains and physical parameters
Throughout this section, we consider the three domains sketched on Fig. 1 to illustrate our results: Rounded diamond: The rounded diamond (a) is defined by the parametric curve ζ(θ) = 2 e iθ + 0.066e −3iθ , for θ ∈ [0, 2π]. It is an algebraic domain of class Q from [7]. This shape satisfies condition 1, as well as the hypotheses of Theorem 5. Narrow ellipse: The ellipse (b) semi-axes are on the X 1 -and X 2 -axes and are of length a = 1 and b = 5, respectively. It is algebraic but not asymptotically a circle in the sense of [7].

Five-petal flower:
The flower (c) is defined by = 2 + 0.6 cos(5θ) in polar coordinates. It has Cartesian equation in the rescaled (X 1 , X 2 ) 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. All three domains D = z + δ B are centred at the origin (z = 0) for simplicity. We set the size of the nanoparticle to be δ = 10 −8 m. The numerics are performed on the rescaled domain B and the homogeneous medium is taken to be vacuum (ε m = ε 0 and μ m = μ 0 ). The physical parameter values are summarised in Table 1.

Modes contribution decay
It was shown in Sect. 3.1 that the scalar products F, φ j H * (∂ B) decay very rapidly when d = 3. In a two-dimensional setting, the theoretical framework is not as clear, but we check numerically that the contribution the modes decrease quite fast with j. Recall that the weight of the j th mode is given by the scalar product F, φ j H * (∂ B) , which, in a low-frequency regime, can be approximated as ν · ∇u in , φ j H * (∂ B) (see Lemma 15). On panel (a) of Fig. 2 we show on all examples that ν · ∇u in , φ j H * (∂ B) decays as j grows. We average over all possible directions d of the incident field. Panel (b) of the same picture shows that the modes themselves, S ω p δc −1 B [ φ j ](X ), decrease as j increases. We average here over all observation positions, X belongs to a circle of radius 100 centred at z = 0.

Plasmonic resonances
We plot the first-order corrected plasmonic resonances with positive real parts on Fig. 3. The resonance radius R(δ) from definition 6 is drawn as a red vertical line on the three subplots and is shown to encompass all the low-frequency resonances. Fig. 3 We plot, for the diamond (a), the ellipse (b) and the flower (c), the two-dimensional first-order corrected resonances with positive real parts from (16): 20. These resonances lie in the lower part of the complex plane and their real part is between ω p /4 and ω p (and smaller than R(δ)). Their negative counterparts are symmetric with respect to the imaginary axis We can then verify a posteriori that our choice of size δ is consistent by checking that R(δ) is still in the low-frequency region, see Table 2.

Validation of theorem 5
In this section, we validate the two-dimensional approximation of the scattered wave in the time domain given in Theorem 5 by plotting the asymptotic result against full numerical simulations.
We sketch the simulation setting with the ellipse in Fig. 4a. We define three observation points A, B and C on a circle of radius 150 nm (|X | = 15) and one observation point D on a circle of radius 3000 nm (|X | = 300). They are characterised by their angle with respect to the x-axis: The nanoparticle is illuminated by a plane wave of the form u in (X ) = e ik m d·δ X f (ω) where f is the Fourier transform of a bump function compactly supported in the interval [0, C 1 ], with C 1 = 8 fs. We plot the time domain incoming wave in Fig. 4b. To ease the notations we drop the tilde subscript in the following and write u(X ) instead of u(X ).

Reference solution
We call reference solution the low-frequency part of the scattered field in the time domain. We first uniformly discretise our frequency domain I ω in L = 10 4 points, with We compute the scattered field in the frequency domain using the representation formula (3). The singlelayer potential is approximated using N = 2 8 equally-spaced discretisation points along the boundary ∂ B. We define the dimensionless frequency ω = ωδc −1 . The reference solution is computed by taking the truncated inverse Fourier transform

Asymptotic solution
The expansion is obtained by summing the first J = 30 modes. Using Theorem 5, the modal approximation of order J becomes: The simulation results are shown in Figs. 5, 6, 7 and 8 . To corroborate our pole expansion, we plot the real part of the reference solution (25) against the real part of the asymptotic one (26) for the different domains and from different observation points.

Comparison in the far-field for the diamond
We begin with the diamond, since it is the shape that satisfies the hypotheses of Theorem 5. Figure 5 shows the field scattered by the diamond, measured in the far-field at position X = D. The reference solution is nicely approximated by the sum of four modes (4, 5, 6 and 7). Figure 6 shows the field scattered by the ellipse, measured at position X = A on panel (a) and X = C on panel (b). In both cases the time domain scattered wave (blue line) is well approximated by the sum of decaying modes (orange symbols). Although we compute the first 30 terms of the modal expansion, the actual number of modes which contribute significantly to approximate the reference solution is much smaller. Indeed, only 1 mode is necessary to reconstruct more than 99% of the signal in Fig. 6.

Extension to a nearer-field for the ellipse and flower
When the observation point is at X = B, we illustrate in Fig. 7 that two modes are needed to match the reference solution for the ellipse. Mode 1, corresponding to a dipole which radiates most of the energy along the x-axis, is associated to the eigenvalue λ 1 = 0.33. Mode 2 corresponds to the dipole which radiates most of the energy along the y-axis and is associated to the eigenvalue λ 1 = −0.33. Mode 1 oscillates slightly faster than mode 2, resulting in the double oscillation visible on the lower plot. These numerical simulations are in line with [9]. Even relatively close to the particle (the observation distance is about a tenth of the wavelength), only two modes radiate in the far-field. Figure 8 shows that even for the non-algebraic flower shape, the scattered wave (blue line) is well approximated by the sum of a small number of decaying modes (orange symbols). As anticipated by Fig. 2, the modes decay being faster for the ellipse than it is for the flower, a larger number of modes is needed for the latter. In Fig. 8, eight modes were needed to reconstruct more than 99% of the reference solution (and five modes sufficed for 95%).

Reference solution Asymptotic solution
Reference solution Asymptotic solution

About the high frequencies
On Fig. 9 we show that the low-frequency part of the time domain solution is actually a good approximation of the full solution, as mentioned in Remark 11. It is completely non-trivial, as we have no information on the localisation of poles for the resolvent in the frequency domain outside the low-frequency range. It seems that there are no more resonances in the high-frequency range due to the dispersive nature of the material. This will be investigated in future work.

About the computational cost
We note that, because a small number of modes usually suffices to approximate the reference solution, the computation cost of the asymptotic solution is relatively cheap. The time needed to compute the reference solution and the asymptotic one are linear in L(= 10 4 ) and J (= 30), respectively. Thus, the time to compute the asymptotic solution is much smaller than the time to compute the reference solution, namely, a hundred times smaller. Moreover, the modes can be pre-calculated and one can compute at a low cost the response of the particle to any given illumination in the time domain.

Concluding remarks
In this paper, we have shown that it is possible to define quasi-normal modes (similar to the ones found in the physics literature) for small plasmonic particles using the spectral decomposition of the Neumann-Poincaré operator and some perturbative spectral analysis.
We have proved that, in a three-dimensional setting, only a few modes are necessary to represent the solutions of the scattering problem by a strictly convex plasmonic particle and that these types of representations can give a very good approximation of the field in the time domain. Our numerical simulations have corroborated the validity of this approach in the two-dimensional case. This theoretical and numerical framework can be adapted to handle more complex systems with multiple particles (see [6]). This work needs to be extended to solutions of Maxwell's equations and to dielectric structures. This will be the subject of forthcoming papers.

A.1 Definitions and notations
Definition 7 Denote by Γ k the outgoing Green's function for the homogeneous medium, i.e., the unique solution of the Helmholtz operator: satisfying the Sommerfeld radiation condition. In three dimensions, Γ k is given by In two dimensions, it is given by 0 is the well-known Hankel function of the first kind and order 0.

Lemma 5 The Hessian matrix of the outgoing fundamental solution in three dimensions
Definition 8 For a function φ ∈ L 2 (∂ D), we define the single-layer potential by and the Neumann-Poincaré operator by When k = 0, we just write S D and K * D for simplicity.

The operator K *
D is self-adjoint in the Hilbert space H * (∂ D) which is H − 1 2 (∂ D) equipped with the following inner product: with ϕ 0 being the unique (in the case of a single particle) eigenfunction of K * D associated with eigenvalue 1/2 such that ϕ 0 , χ(∂ D) − 6. The following trace formulae hold for φ ∈ H − 1 2 (∂ D): where I denotes the identity operator. 7. The following representation formula holds:

A.3 Invertibility of the boundary operators
Lemma 7 For k small enough, the three-dimensional single-layer potential In two dimensions, the single-layer potential S D : is, in general, not invertible. All the proofs for the following lemmas can be found in [6].

Lemma 9
For k small enough, the two-dimensional boundary operator S k D : is invertible and

B Scaling properties for a finite volume particle
For each function f defined on ∂ D, we define a corresponding function on ∂ B by f (X ) := f (z + δ X ).

Lemma 11 It holds that
Lemma 12 For f , g defined on ∂ D, corresponding to f , g, respectively, we have Proof In three dimensions, by straightforward calculations we have .

C.1 Asymptotic expansions of the boundary operators
Lemma 13 1. The three-dimensional single-layer potential and its inverse admit the following expansions in the quasi-static limit kδ → 0: The two-dimensional single-layer potential and its inverse admit the following expansions in the quasi-static limit kδ → 0: The Neumann-Poincaré operator in three dimensions admits the following expansion in the quasi-static limit where, for φ ∈ H − 1 2 (∂ B),

The Neumann-Poincaré operator in two dimensions admits the following expansion in the quasi-static limit
Proof The proof can be found in [6].

C.2 Proof of lemma 1
Proof Recall that For the three-dimensional case, using Lemma 13 we have by straightforward calculations For the two-dimensional case, we have Hence, We have that

C.3 Proof of Lemma 4
Proof Using the Taylor expansion we compute, for X ∈ ∂ B, where we used 1

C.4 Proof of Lemma 15
Proof Using the Taylor expansion we compute, for X ∈ ∂ B, where we used 1 2 It is immediate to see that 1 2 /∂ν X − = ν X and using jump conditions for

D Proof of Theorem 5
We need the following lemma:
The goal of this section is to establish a resonance expansion for the low-frequency part of the scattered field in the time domain. Introduce, for 0 < < ρ, the truncated inverse Fourier transform of the scattered field u sca given by Recall that z is the centre of the resonator and δ its radius. Let us define the time it takes to the signal to reach first the scatterer and then observation point x. The term ±2δ/c accounts for the maximal timespan spent inside the particle.
Proof We have For j ≥ 1, let us compute the contribution of one mode Ξ j (x, ω). We want to apply the residue theorem to get an asymptotic expansion in the time domain. Note that: where the integration contours C ± ρ and C ± are semi-circular arcs of radius ρ and , respectively, in the upper (+) or lower (-) half-planes, and C ± is the closed contour defined as The integral on the closed contour is the main contribution to the scattered field by the mode and can be computed using the residue theorem to get, for ρ ≥ max j∈{1,...,J } [Ω ± j (δ)] and 0 < ≤ min j∈{1,...,J } [Ω ± j (δ)], Since Ω ± j (δ) is a simple pole of ω → 1 λ(ω) − λ j (ωδ) we can write: To compute the integrals on the semi-circle, we introduce: Given the regularity of the input signal f ∈ C ∞ 0 ([0, We have that t − 0 + C 1 ≤ T ≤ t + 0 − C 1 . Two cases arise. Case 1: For 0 < t < t − 0 , 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: |X j (y, ρe iθ )|dσ (y)dθ, where we used that for θ ∈ [0, π/2], we have sin θ ≥ 2θ/π ≥ 0 and − cos θ ≤ −1 + 2θ/π. The usual way to go forward from here is to take the limit ρ → ∞, and get that the limit of the integral on the semi-circle is zero. As in the three-dimensional case, we work in the quasistatic approximation here, and our modal expansion is not uniformly valid for all frequencies.
So we have to work with a fixed maximum frequency ρ. However, the maximum frequency ρ depends on the size of the particle via the hypothesis ρ ≤ cδ −1 . 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 ]: For the upper-half semi-circle of radius , we also transform into polar coordinates with the change of variable Ω = e iθ , for θ ∈ [0, π], and get: .
Exactly as in Case 1, we cannot take the limit ρ → ∞. However, the maximum frequency ρ depends on the size of the particle via the hypothesis ρ ≤ cδ −1 . Using the fact that 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 ]: For the lower-half semi-circle of radius , we also transform into polar coordinates with the change of variable Ω = e iθ , for θ ∈ [0, π], and get: . The result of Theorem 5 is obtained by summing the contribution of all the modes.