Eﬀicient and intuitive method for the analysis of light scattering by a resonant nanostructure

: We present a semi-analytical formalism capable of handling the coupling of electromagnetic sources, such as point dipoles or free-propagating fields, with various kinds of dissipative resonances with radiation leakage, Ohmic losses or both. Due to its analyticity, the approach is very intuitive and physically-sound. It is also very economic in computational resources, since once the resonances of a plasmonic or photonic resonator are known, their excitation coefficients are obtained analytically, independently of the polarization, frequency or location of the excitation source. To evidence that the present formalism is very general and versatile, we implement it with the commercial software COMSOL, rather than with our in-house numerical tools.


Introduction
Efficient interconversion of free-propagating and confined fields is a cornerstone of wave engineering for monitoring absorption, radiation, emission of acoustic, quantum, radio or optical waves. The interconversion can be efficiently performed by exciting resonances that are either delocalized or localized in space and by tuning the frequency of the incident wave to match the resonance energy. For electromagnetic waves in closed systems without absorption, the resonances are the eigenstates of Hermitian time-evolution operators. These so-called normal modes have real eigenfrequencies and infinite lifetimes and form a complete orthogonal set [1]. As a consequence, they can be used as a convenient basis to describe the physics of closed and lossless systems. For instance, the density of electromagnetic states reduces to a sum of Dirac distributions each centered at the eigenfrequencies.
However, in most (if not all) cases, the system is dissipative and losses its energy by radiation (open system) and/or by absorption. A non-conservative system is not described by Hermitian operators and the usual normal-mode expansion is therefore not applicable. The physics of dissipative systems has to be described with the eigenstates of the wave equation that obey outgoing-wave boundary conditions at infinity. Consistently with the literature on open systems [2,3], we will refer to these modes as quasi-normal modes (QNMs), to distinguish them from the usual normal modes. QNMs possess a finite lifetime τ and their eigenfrequency ω  is a complex quantity with Im( ω  ) = 1/τ. The exp(iωt) notation is adopted for the time-harmonic fields. The eigenfrequency is a pole of the Green tensor of the system and it admits a classical interpretation in terms of the quality factor the resonance, which can be attached to energy balance arguments on the energy stored and leaked at resonance at least for small Ohmic losses [4]. Note that the theory of nonconservative systems is not completed yet. The orthogonality and the completeness of QNMs has been demonstrated only for very simple open dielectric systems [2,3,5]. Recently, the orthogonality of QNMs has been demonstrated for any three-dimensional resonant system with absorption, provided that spectral dispersion can be neglected [10]. The completeness of QNMs in the more general case is still an open question. Nevertheless, electromagnetic QNMs are at the core of many fascinating optical phenomena and designs: Fig. 1 illustrates our purpose on the textbook case of "Mie" resonance of a gold nanorod, which exhibits a characteristic Lorenzian response due to the complex pole associated to the electric-dipole resonance of the nanorod. Despite the widespread use of optical resonance for various designs, a difficulty arises when handling the QNMs in practice. This difficulty has prevented a widespread use of QNMs for design until now. Actually, because energy is continuously lost, the QNM amplitudes decay in time,

( )
Im 0 m ω >  . As a consequence, the wavevector is also complex with a positive imaginary part and, since the boundary conditions impose outgoing spherical waves exp(-ikr)/r  at infinity, the field diverges exponentially as → ∞ r . Indeed, QNM normalization becomes problematical. In particular, the divergence poses a serious problem when theoretically attempting to derive a self-consistent formalism capable of modeling how incident waves exchange their energy with the QNMs, and vice-versa. The issue has been risen especially in the context of the so-called Purcell effect [6], see [7] for photonic-crystal cavities and [8,9] for metallic nanoantennas.
The normalization problem has been solved recently for the general case of threedimensional absorptive and dispersive systems [10], allowing some of the authors to propose a Purcell formula that is accurate for open systems with strong radiation leakage or absorption and with dispersive materials. The new formalism removes a longstanding theoretical problem associated to the definition of the mode volume for dissipative (non-Hermitian) systems and predicts unexpected behaviors, such as the fact that a detuning between the emitter and cavity frequencies does not necessarily result in a Lorentzian lineshape response, as could be naively anticipated form a pole-type response. It is very general since it encompasses situations for which several QNMs are required to predict the response of the system and it is highly accurate, as was evidenced by comparison with fully-vectorial numerical results for distinct plasmonic and photonic nanocavities of current interest in nanophotonics.
The new formalism relies on the numerical calculation of a few dominant QNMs, but once the QNMs are known, it is fully analytical and thus extremely fast. For instance, the Purcell effect (or the power radiated by the dipole source) is explicitly known with a closedform expression for any frequency detuning, location or orientation of the dipole source, in sharp contrast with full-electromagnetic treatments that require new calculations for each case. Thus the new formalism has the potential to become a valuable tool for design.
The formalism in [10] was developed to describe the coupling between a dipole source and a resonance; the theoretical derivations of the main results were emphasized and the numerical implementation was performed with in-house software with a complete know about. This article complements the previous work with two important new results: • we generalize the approach to encompass more general interconversion processes between free-propagating and localized fields. In particular we derive new formulas for the coupling of QNMs with plane waves (or free-propagative beams in general), and obtain analytical formulas for the extinction and absorption cross sections, two fundamental quantities of scattering processes.
• we also propose a method to normalize the QNMs that is more general than the one in [10]. The net benefit is the derivation of a versatile and practical method for normalizing the QNM, which is suitable for any software and in particular commercial electromagnetic solvers. We illustrate how to calculate normalized QNMs with COMSOL Multiphysics (http://en.wikipedia.org/wiki/COMSOL_Multiphysics). The manuscript is organized as follows. Section 2 is a reminder of important results obtained in [10]. It is presented here for the sake of self-consistency.
In Section 3, we consider the general problem of the calculation and normalization of QNMs. Although the calculation is not an original task in itself, it is often conducted with dedicated softwares by specialists of computational electrodynamics. Here, we show how to calculate them with standard tools, taking the example of COMSOL Multiphysics for which we provide a model sheet. We then tackle the important problem of the normalization. The method developed in [10] used the field calculated in the perfectly-matched-layer (PML) region, for PMLs implemented as closed-form constitutive tensors [11]. Here we provide a simpler and more general method that is independent of the numerical technique used for satisfying the outgoing wave conditions (transparent boundary conditions, PMLs with linear or even non-linear coordinate transforms or wave absorbers). The new method is expected to remain valid even for the split-field PML often used with FDTD [12].
In Section 4, we generalize the formalism that was primarily devoted to the coupling with a dipole source to encompass more general interconversion processes between freepropagating and localized fields. In particular we derive new formulas for the coupling of QNMs with plane waves (or free-propagative beams in general), and obtain analytical formulas for the extinction and absorption cross sections, two fundamental quantities of scattering processes. The derivation is supported by comparative tests provided with COMSOL Multiphysics. Section 5 concludes the work.

Theoretical background
Throughout the manuscript, the exp(iωt) notation is adopted for the harmonic fields (consistently with the COMSOL software). We are interested in the response of a resonant optical system that consists of structured absorbing and dispersive media with open boundaries. The system is characterized by the position-and frequency-dependent permittivity and permeability tensors ε(r,ω) and µ(r,ω). Note that the materials can be anisotropic and magnetic; the sole assumption is that the materials are reciprocal ε Τ = ε and µ Τ = µ, where the superscript denotes matrix transposition.
and that satisfy outgoing wave boundary conditions (the Sommerfeld radiation condition at |r| → ∞). Hereafter, the tilde grapheme is used to quote quantities related to QNMs and we will denote by   1. (a) Excitation of a 30-nm-wide and 100-nm-long gold nanorod by either a plane wave or a dipole located on axis at a 10-nm distance above the rod. The dipole and the plane wave are polarized along the rod axis. The gold nanorod is embedded in a host medium of refractive index n = 1.5. (b) Dipole illumination: Total spontaneous decay-rate of the dipole normalized by the decay rate in a bulk material with n = 1.5. (c) Plane-wave illumination: Extinction cross-section. In (b), blue triangles and black circles are fully-vectorial computational results obtained with COMSOL for a coarse mesh with 41543 elements and a fine mesh with 875575 elements, respectively. The results obtained with the coarse mesh are in good agreement with the ones obtained with the fine mesh. In (c), black circles are fully-vectorial computational results obtained with COMSOL and a fine mesh. In (b) and (c), the red solid curves represent the predictions of the method proposed in the present work with Eq. (9) and (14), assuming that the excitation is due to a single resonance, namely the fundamental electric-dipole resonance of the nanorod. The gold permittivity is given by a Drude model (see text).
In [10], we studied the general problem of the radiation by an electric current source J(r) located in the vicinity of the resonant system. The field ( ) ( ) , , ω ≡ Ψ r E H radiated by the source at the frequency ω verifies, Assuming that

( )
,ω Ψ r can be expanded as a superposition of QNMs, where ( ) m f ω is a non-resonant background that becomes negligible as m ω ω →  and comes from the fact that the QNMs of dispersive media are not orthogonal. In the particular case of non-dispersive media, ( ) m f ω is strictly equal to zero. Equation (3) is obtained by applying the general form of the Lorentz reciprocity theorem (valid for absorbing and dispersive media) to Eqs. (1) and (2), and it is valid provided that the QNMs obey the normalization condition Equation (4) is not as straightforward as it seems, simply because the QNM diverges exponentially for |r| → ∞. In [10], it was shown that the normalization cannot be performed by calculating the integral in the real space. The normalization has to be performed in the complex space by transforming resonance states, which are formally scattering or continuum states, into square-integrable states. The integral can be calculated by using classically available Perfectly-Matched-Layers (PMLs) provided that the integral calculation be performed over all the numerical space, including the PML region. According to [10], it appears that this approach is valid for PMLs implemented as closed-form constitutive tensors. Related complex continuations can be found in quantum mechanics for systems with nonhermitian Hamiltonians [13]. In the next Section, we develop a much-more-general numerical normalization method that does not require the calculation of the integral of Eq. (4). The expansion of a field into a small set of resonances can be problematical, as we do not even know if the QNMs form a complete set. Anyway, the prime motivation is not to derive a novel method to calculate the exact value of the radiated field, but rather to have a simple decomposition of the response into a small set of resonances. Such an expansion is extremely useful to get a deep insight into the problem and, since it can be handled analytically, it enables also ultra-efficient electromagnetic simulations specifically tailored for design and fast optimizations. Hereafter, a single QNM will be considered for simplicity. Thus we will drop the subscript m and denote by ω  and α the unique pole and coupling coefficient, respectively. For multiple-resonance systems, the new expressions that we derive for α hereafter can be safely used. They are strictly valid for non-dispersive materials. For dispersive materials, it can be shown that they are only approximate [10]. However, they remain highly accurate and using the exact formula does not lead an improvement of the accuracy of the prediction, at least in all cases we have tested so far.

Calculation and normalization of the QNMs
In this Section, we propose a simple and versatile technique to calculate and normalize the QNMs. We believe that the approach can be implemented with various (if not all) numerical software implementing the outgoing wave conditions for |r| → ∞. For the sake of illustration, we use the COMSOL Multiphysics software and its RF toolbox for which we provide a model sheet for the nanorod geometry of Fig. 1 [14].

Scattering formulation
A classical technique to calculate a QNM of a given resonant system consists in exciting the system with a source (a plane wave, a dipole source located nearby …) at a frequency ω close to the QNM frequency ω  , to calculate some quantity representative of the response of the system and to iteratively change the excitation frequency ω in the complex plane so that the response diverges. Then the frequency ω will tend towards ω  and the response will tend toward the QNM. By doing so, we can take into account the metal dispersion easily, provided that we use an analytical continuation of the permittivity in the complex frequency plane. In the classical scattering formulation [15], the permittivity ( ) ,ω ε r is often decomposed represents a background permittivity and ( ) ,ω Δε r is null everywhere except in the resonant structure. Note that ( ) , b ω ε r does not necessarily correspond to a homogeneous medium and that the generalization to magnetic materials is straightforward. In the absence of the resonant structure, the background field satisfies the following equations whereas in the presence of the structure the total field (E, H), satisfies ( ) where J(r) is the source term, potentially located at infinity. The field

Pole calculation
To calculate the pole, one needs to iteratively solve Eqs. (5) and (6) for several values of ω that progressively converge towards ω  . In our simulations, we use a Drude model, The values of ω p and γ are fitted from tabulated data in [16], and ( , ) ω μ r is simply equal to the permeability of the vacuum. Actually, COMSOL does not handle complex frequencies and we use a trick relying on effective dielectric and magnetic tensors to artificially mimic complex frequencies with real ones, see details in Appendix 1. The trick is likely to be usable with any software that considers real frequencies and complex permittivity and permeability. In our simulations, we have used in the three-dimensional space either a coarse mesh with 41543 elements or a fine mesh with 875575 elements.
As an initial guess value for the iterative pole search, we take the complex frequency 2πc/  Table 1 illustrates the rapid convergence achieved for the coarse mesh. The CPUtime per iteration is about 30s (it slightly increases as ω ω →  ) on a 4-core desktop computer equipped with a 3Ghz-processor. The data are obtained for a dipole current J z = 10 −11 A⋅m, located on axis at a 10-nm-distance above the nanorod and can be reproduced by strictly following the procedure in Appendix 2. The other numerical parameters are those of the caption of Fig. 1. The iterative procedure is converging fast. In a few iterations (2 on the example of Table 1), |E z (r 1 )| is increased by five orders of magnitude. For the last pole estimate obtained with Eq. (A3) for n = 4, COMSOL cannot compute the radiated field in less than 60 minutes. We believe that this trend that we have systematically observed in all our calculations is due to the numerical algorithm used by COMSOL to invert a matrix that becomes singular as ω ω →  . The resonance frequency of Table 1 is calculated with a coarse mesh; a more accurate value obtained with a thinner mesh is 2.047183121604 × 10 15 + i0.1046416051442 × 10 15 s −1 . This is intentional, since our purpose is to show that the present semi-analytical method is very fast. Actually as shown in Figs. 1(b), 1(c), 2, 3(b) and 3(c) the main physical quantities attached to the scattering by resonance are accurately predicted even if the resonance frequency is only approximately calculated with the coarse discretization. This evidences the soundness of the approach.

Normalization of the QNMs
The QNM normalization can be performed by calculating the volume integral defined in Eq. (4) using PMLs implemented as closed-form constitutive tensors. Using the fact that the integrant in Eq. (4) is an invariant of coordinate transforms in Maxwell's equations and provided that the complex stretching coefficient of the PML is large enough, the integral over the whole computational domain (including the PML region) does not depends on the actual choice for the PML [10]. In this Section, we describe a more general method that can be implemented with any type of PML, transparent boundary conditions or even with methods that do not require numerical boundary condition.
As discussed earlier, for ω ω ≈  , the system response is dominated by the excitation of the QNM. It is thus an excellent approximation to consider that the scattered field Equation (8) is very useful as it allows us to calculate a normalized QNM from the sole knowledge of the field scattered by a resonant structure at any frequency ω close to the QNM frequency ω  . Note that in practice (particularly for nanostructures with a low Q-factor), the frequency ω has to be taken in the complex plane. Fast convergence towards a stable distribution is achieved. To further check the accuracy of the normalization, we calculate the total normalized spontaneous decay-rate P(ω) of a zpolarized dipole ( ) 0 δ − J r r located on axis at a 10-nm-distance above the rod. The red solid curve in Fig. 1(b) is obtained from the semi-analytical expression [10].
where n = 1.5 is the refractive index of the background, c is the speed of light in vacuum, α is given by Eq.  ω ω =  (see Table 1

Absorption and scattering cross sections
In [10] and up to now, we were concerned with the coupling of one or several resonances with a dipole source. In this Section, we generalize the QNM formalism to encompass the important cases of the excitation of resonant systems by free-space beams. The generalization, as we shall see now, allows us to analytically handle the textbook case of extinction and absorption cross-section coefficients [15].

Problem formulation
As shown in Fig. 3(a), we assume that the nanorod is illuminated by an incident field, denoted by ( ) Equation (10) tells us that the field scattered by the resonant structure at frequency ω can be seen as the field radiated in the presence of the structure by a source distribution Δε r E r , a known quantity that is solely depending on the incident illumination. Incidentally, note that Eq. (10) is used by COMSOL with the so-called scattered formulation.
As before, we solve Eq. (10) by assuming that the scattered field Then identifying the source J used to derive the coupling coefficient in Eq. (3) with the new distribution source Δε r E r , one straightforwardly obtains for the complex coupling coefficient where the integral runs over the volume of the scatterer and ( ) g m ω is a non-resonant term.
The coefficient β m is easily calculated as an overlap integral between the normalized QNM and the incident field. Note that the integral has a simpler expression for a homogeneous scatterer where ( ) ,ω Δε r is independent of r and can be put outside the integral.

Absorption cross-section
The absorption cross-section σ A of the plasmonic nanorod can be calculated with the following formula [15] ( ) ( ) ( ) considering the approximate expressions of Eqs. (11) and (12) for the scattered field. In Eq. (13), S 0 is the time-averaged Poynting vector of the incident plane-wave and the integral is performed over the volume V of the resonant structure, where the QNM expansion is a very good approximation of the scattered field. Figure 3(b) shows the absorption cross-section spectra of the gold nanorod for a plane wave incident normally to the rod axis with a polarization parallel to the axis. The solid-red curve represents the predictions obtained with Eq. (13), assuming that only the fundamental electric dipole resonance is excited. Quantitative agreement is achieved with the fullyvectorial calculations (circles) obtained with COMSOL for the fine mesh with the scattered field formulation.

Scattering cross-section
Starting from Eq. (11), there are several ways to calculate the scattering cross-section σ S or the extinction cross-section σ E , and because Eq. (11) is approximate they do not have all the same accuracy. For instance as |r| → ∞, the scattered field vanishes as 1/|r| 2 far away from the scatterer, whereas the finite sum in Eq. (11) exponentially diverges. As we checked, calculating the scattering cross-section with the classical formula [15] #194819  E H n , where A is a closed surface surrounding the scatterer, is inaccurate when A is far away from the scatterer and it slightly depends on of the specific choice of A when A is close to the scatterer. From our tests, we find that best accuracy and consistency are achieved by first calculating the extinction cross-section with the following expression (14) which is valid if the resonant structure is immersed in a transparent medium, ( ) (14) is not common but is directly found from the classical formula ( )  E H E H n by applying the Lorentz reciprocity theorem.
Then σ S is obtained by [15] S E A . The solid-red curve in Fig. 3(c) shows σ S calculated with Eqs. (14) and (15). Again quantitative agreement is achieved with fully-vectorial data obtained with the fine mesh, especially for the peak values for decomposed over the QNM set as a sum of real numbers that are not necessarily all positive. As evidenced in [10], only the total sum is positive but every QNM contribution might be either positive or negative, and this may lead to complex Fano-like responses (see Figs. 3b and S1 in [10]).

Purcell factor and Green tensor
The scattering formulation can also be used to calculate the field radiated by a point dipole We calculated the predictions given by Eq. (17) for a dipole located on axis at a 10-nm distance above the rod, and found that the predictions of Eqs. (9) and (17) are indistinguishable. The reason comes from the fact that the scattered field and the total field are almost identical for dipoles that are located closed to the rod and that efficiently excite the dipole resonance of the rod. However, Eq. (17) that relies on an expansion of the scattered field is more accurate than Eq. (9) that relies on an expansion of the total field when the coupling between the rod resonance and the dipole source is small. The inaccuracy reported in Fig. S3 of the Suppl. Inf. in [10] with Eq. (9) when the direct coupling to free-space becomes dominant as the nanorod-source separation-distance increases is removed with the present scattering formulation. The main advantage of Eq. (9) is that it leads to a fully analytical result to define a mode volume for absorptive and dispersive resonators [10].

Conclusion
Dissipative resonances are essential ingredients for many physical effects and devices. We have presented a numerical formalism that is able to analytically handle the coupling between optical resonances and various kinds of sources, such as Dirac electrical dipoles − Eqs. (3) and (9) − or free-propagating fields − Eqs. (13) and (14). The formalism encompasses radiation leakage and Ohmic losses. It has been implemented with the COMSOL software, but the implementation is very versatile and can be made with many other in-house or commercial software. We emphasize that the new normalization procedure of Section 3.3 is much simpler and more general than the one initially presented in [10]. For the sake of simplicity, the formalism has been tested on a very simple example, a gold nanorod in a uniform dielectric background. In [10], it was shown that photonic-crystal cavities and metallic nanoantennas can be analyzed with a high accuracy and efficiencies. These are just examples and we believe that much more complex resonances, such as those produced by disordered systems, may be analyzed as well.
The present formalism is very economic in computational resources. Perhaps it is worth emphasizing that the calculation of all the solid curves in Figs. 1(b), 1(c), 3(b) and 3(c) is performed in 4 minutes without using symmetry. As shown in Section 3, the computation can be performed with efficient iterative algorithms that converge in a few iterations. Once the quasi-normal modes are known, then their excitation coefficients, see Eqs. (3) and (12), are obtained analytically for any source frequency, polarization, location … For instance, if one needs to study the extinction cross-section of a nanoparticle for various frequencies, angles of incidence or polarizations, one generally needs to repeat many independent fully-vectorial computations for each parameters, in sharp contrast with the present formalism.

Appendix 1: Implementation of complex frequencies with COMSOL
Let us consider the following set of equations, where ω is potentially complex Equations (18) and (19) are fully equivalent and admit the same solution. The trick is likely to be usable with any software that considers real frequencies and complex permittivity and permeability. The real frequency L ω may take any arbitrary value. In our calculations, we simply take . Note that, during the iterative calculation discussed in Section 3 and Appendix 2, the real frequency L ω at which the calculation is performed with COMSOL is fixed; only the effective permittivity and permeability are iteratively changed.

Appendix 2: Iterative approach for the pole calculation
The pole computation relies on an iterative procedure. Assuming that the inverse of one of the electromagnetic field component denoted Z(ω) can be approximated by  Table 1. We start by estimating the first complex frequency 1 ω  from Fig. 1  ω ω ω ω , we define the new pole estimate 2 ω  as the frequency that corresponds to the smaller value of |Z|. Note that 2 ω  might be equal to 1 ω  . From the quadruplet, we further eliminate the frequency with the largest |Z|. By iterating the process, we generate a frequency series, 1 ω  , 2 ω  , 3 4 , ...
ω ω   , which rapidly converges toward the actual pole value ω  in general, as illustrated in Table 1.