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.
© 2013 Optical Society of America
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 . 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 of the resonance, which can be attached to energy balance arguments on the energy stored and leaked at resonance at least for small Ohmic losses . Note that the theory of non-conservative 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 . 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, . As a consequence, the wavevector is also complex with a positive imaginary part and, since the boundary conditions impose outgoing spherical waves at infinity, the field diverges exponentially as . 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 , see  for photonic-crystal cavities and [8,9] for metallic nanoantennas.
The normalization problem has been solved recently for the general case of three-dimensional absorptive and dispersive systems , 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 closed-form 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  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 . 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 . 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  used the field calculated in the perfectly-matched-layer (PML) region, for PMLs implemented as closed-form constitutive tensors . 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 .
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 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. The derivation is supported by comparative tests provided with COMSOL Multiphysics.
Section 5 concludes the work.
2. 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. The QNMs supported by such a resonant system are electromagnetic field distributions that are solutions of Maxwell’s equations in the absence of source,
In , 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 radiated by the source at the frequency ω verifies,10] is the derivation of a closed-form expression for the complex coupling coefficient αmEquation (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 conditionEquation (4) is not as straightforward as it seems, simply because the QNM diverges exponentially for |r| → ∞. In , 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 , 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 non-hermitian Hamiltonians . 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 . 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.
3. 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 .
3.1 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 , the permittivity is often decomposed as , where represents a background permittivity and is null everywhere except in the resonant structure. Note that 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
3.2 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, , with ωp = 1.26 × 1016 s−1 and γ = 1.41 × 1014 s−1, to define an analytical continuation of the gold permittivity for complex frequencies. The values of ωp and γ are fitted from tabulated data in , and 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/1 ≈921 −48.6i nm (Q = 9.5) estimated from the data obtained with the fine mesh in Fig. 1b. We then use an iterative method to calculate the pole. The method relies on the fact that, as the frequency approaches the pole frequency , one may look for the pole by searching the complex zero of . We use a Padé approximation around the pole , , and iteratively determine a0, a1 and with the method described in the Appendix 2.
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 Jz = 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), |Ez(r1)| 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 × 1015 + i0.1046416051442 × 1015 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.
3.3 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 . 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 is proportional to the QNM . From Eq. (3), the proportionality factor isEq. (7) by the current source J and using r = r0, one gets . By introducing this value into Eq. (7), one obtains an expression for the normalized modeEquation (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.
Figure 2 shows the absolute value of different distributions of the Ez-component. The distributions are invariant by rotation along the z-axis, and only one half of a cross-section plane is shown. They are obtained for the pole estimates using Eq. (8), with (n = 1,…4) taken from Table 1 and for . 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 z-polarized dipole 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 .Eq. (3) and the QNM field is the one calculated with Eq. (8) for. An excellent agreement, comparable with that obtained in  with a normalization technique relying on the numerical calculation of the integral of Eq. (4), is achieved with the black circles that are directly obtained with fully-vectorial computations performed with COMSOL.
4. Absorption and scattering cross sections
In  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 .
4.1 Problem formulation
As shown in Fig. 3(a), we assume that the nanorod is illuminated by an incident field, denoted by , a plane wave for instance. Following Section 3.1, we consider the classical scattering formulation, for which we define the background permittivity and the scatterer permittivity change . From Eqs. (5) and (6), by difference, one obtains the following equations for the scattered fields andEquation (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 , 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 can be well represented by a superposition of QNMsEq. (3) with the new distribution source , one straightforwardly obtains for the complex coupling coefficient
4.2 Absorption cross-section
The absorption cross-section σA of the plasmonic nanorod can be calculated with the following formula Eqs. (11) and (12) for the scattered field. In Eq. (13), S0 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 fully-vectorial calculations (circles) obtained with COMSOL for the fine mesh with the scattered field formulation.
4.3 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  , 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 expressionEquation (14) is not common but is directly found from the classical formula by applying the Lorentz reciprocity theorem. Then σS is obtained by 
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. Away from the resonance, the accuracy reduces; negative values are even predicted for λ > 1.1 µm in Fig. 3©. One should not be surprised by the negative values, simply because σS is expressed as a difference of two terms (see Eq. (15)) that are each approximately predicted. As we checked, the inaccuracy is not due to the single mode approximation (including other higher-order QNMs in the expansion has a negligible impact on the prediction) but to a non-resonant background. Far away from any resonance, the scattering is dominantly due to a continuum of states (of plane waves in a uniform material). This background, which is for instance responsible for the well-known 1/λ4 Rayleigh scattering of small particles, is difficult to represent as a sum of QNMs. Such a difficulty also holds for Purcell-factor calculations. When the dipole source is located at very small distances from a metal particle and experiences a direct electrostatic coupling with the metal (quenching), the QNM expansion becomes inaccurate, see Section 3.3. of the Suppl. Inf. in . Finally, note that the Purcell factor and the extinction cross-section are decomposed over the QNM set as a sum of real numbers that are not necessarily all positive. As evidenced in , 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 ).
4.4 Purcell factor and Green tensor
The scattering formulation can also be used to calculate the field radiated by a point dipole located at . The total field radiated by the dipole at the frequency ω can be written asEq. (5) with a dipole source. The field expression of Eq. (16) is directly related to the Green tensor. It can be used to calculate the total spontaneous decay-rate P(ω) of the dipoleEq. (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  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 .
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 .
For the sake of simplicity, the formalism has been tested on a very simple example, a gold nanorod in a uniform dielectric background. In , 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 complexEquations (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 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 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 around the pole, with a0, a1 and three complex variables that are iteratively estimated. We use Z = calculated at point M located on-axis 30 nm below the nanorod to obtain the results in Table 1. We start by estimating the first complex frequency from Fig. 1 (2πc/ ≈ 921 −48.6i nm), and then calculate Z1, Z2, Z3 for three frequencies ω1=(1−10−5), ω2=,ω3=(1+10−5). By solving a linear system of three equations with three unknowns, we get a new frequency given byTable 1.
The authors thank Anthony Jouanin for fruitful discussions and assistance. Qiang Bai acknowledges a fellowship from the LabEx LAPHIA.
References and links
1. P. M. Morse and H. Feshbach, Methods of Theoretical Physics (McGraw-Hill, 1953).
2. R. K. Chang and A. J. Campillo, Optical Processes in Microcavities, Chap. 1 (World Scientific, 1996).
4. P. Lalanne, C. Sauvan, and J. P. Hugonin, “Photon confinement in photonic crystal nanocavities,” Laser & Photon. Rev. 2(6), 514–526 (2008). [CrossRef]
5. P. T. Leung, K. M. Pang, and K. Young, “Two-component wave formalism in spherical open systems,” J. Phys. Math. Gen. 39(1), 247–267 (2006). [CrossRef]
6. E. M. Purcell, “Spontaneous emission probabilities at radio frequencies,” Phys. Rev. 69, 681 (1946).
9. S. Derom, R. Vincent, A. Bouhelier, and G. Colas des Francs, “Resonance quality, radiative/ohmic losses and modal volume of Mie plasmons,” Europhys. Lett. 98(4), 47008 (2012). [CrossRef]
10. C. Sauvan, J. P. Hugonin, I. S. Maksymov, and P. Lalanne, “Theory of the spontaneous optical emission of nanosize photonic and plasmon resonators,” Phys. Rev. Lett. 110(23), 237401 (2013). [CrossRef]
11. F. L. Teixeira and W. C. Chew, “General closed-form PML constitutive tensors to match arbitrary bianisotropic and dispersive linear media,” IEEE Microw. Guided Wave Lett. 8(6), 223–225 (1998). [CrossRef]
12. J. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. 114(2), 185–200 (1994). [CrossRef]
13. A COMSOL model file associated to the article can be found at http://www.lp2n.institutoptique.fr/Membres-Services/Responsables-d-equipe/LALANNE-Philippe.
15. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Pparticles (Wiley 1983).
16. E. D. Palik, Handbook of Optical Constants of Solids, Part II (Academic Press, 1985).
17. W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes: The Art of Scientific Computing (Cambridge University Press, 1989).