Manipulating the Quasi–Normal Modes of Radially Symmetric Resonators

The frequency response of a resonator is governed by the locations of its quasinormal modes in the complex frequency plane. The real part of the QNM determines the resonance frequency and the imaginary part determines the width of the resonance. For applications such as energy harvesting and sensing, the ability to manipulate the frequency, linewidth and multipolar nature of resonances is key. Here, we present a simple analytical tool to control the location and polarity of radially symmetric resonators.


Introduction
Resonances are key to a strong electromagnetic response of photonic components.Over the last two decades resonant interaction has been employed in metamaterials [1][2][3], nano-antennas [4,5] and super-scatterers [6,7] to almost arbitrarily manipulate electromagnetic radiation.Recently there has been much interest in 'Mietronics' [8,9].Work in this field has focused on using higher order multipolar modes to achieve greater control of electromagnetic radiation including directional scattering [10], field confinement [11] and frequency selective imaging [12].
Engineering the resonant frequency of any particular multipolar response is crucial for many applications including efficient absorbers [13,14] and thermal emitters for gas sensing [15].Although the theory of Mie resonances is over 100 years old [16], the problem of designing the spectral response of multipolar resonances remains open.
Many approaches rely on a brute force sweep of the design parameters using memory intensive numerical simulations [17,18] to move resonances, which can be numerically demanding and offers little insight into the final structure.More advanced methods based upon machine learning [19,20] have recently been used to manipulate both electric and magnetic multipole moments, however a large amount of training data must be generated.In addition to this, semi-analytic techniques based upon complex analysis [21] have been developed to enhance the Q-factor of resonances.This exploits the insight provided by the quasi-normal mode framework.
First developed by Gamow to describe alpha decay [22], quasi-normal modes have been employed extensively to model electromagnetic systems in terms of their complex frequency resonances [23][24][25].Crucially, for well separated quasi normal mode resonances, the real part of the quasi-normal mode frequency describes the resonance frequency, whereas the imaginary part encodes the spectral width of the resonance.The ability to manipulate the exact location of quasi-normal modes in the complex plane therefore allows one to control the location and width of the resonance simultaneously.Inverse design techniques to manipulate spatial electromagnetic responses are well developed [26][27][28][29], however few systematic techniques exist to manipulate the spectral response of resonators.
In this work, we develop two techniques to control the spectral location, width, and multipolar nature of the resonances of radially symmetric resonators.Building on previous work that developed a design theory for 1D graded structures [30], we propose a technique to calculate the complex shift one should apply to the permittivity of a resonator to place a particular multipolar resonance at a particular complex frequency.We also employ quasi-normal mode perturbation theory to design graded structures with particular resonances at particular frequencies.Both techniques are simple, easily extensible and numerically efficient.

Finding the Quasi-Normal Modes of Radially Symmetric Resonators
In this section, we briefly review how one can formulate a finite difference scheme to find the quasi-normal modes of a resonator with a radially graded permittivity.Our key results are how these permittivity profiles can be designed, however it will be necessary to find the quasi-normal modes throughout.
The scalar Helmholtz equation governs a wide variety of physical phenomena.For pressure acoustics and a single polarisation of the electromagnetic field in 2D it is exact, and is commonly used an approximation that neglects polarisation in 3D electromagnetism.For simplicity we will work with the scalar Helmholtz equation throughout.We write the scalar Helmholtz equation as describing a single polarisation of the electromagnetic field () in a material with spatially varying permittivity () and wavenumber .Our aim is to find, for a given permittivity profile, the complex eigen-frequencies supported by the structure.Simply re-arranging the Helmholtz equation (1), this problem can be formulated as an eigenvalue problem For this eigenvalue problem to give the complex quasi-normal mode frequencies, one must impose the out-going wave boundary condition.This makes the Laplace operator depend upon the eigenvalue  [23,30], meaning that the solution of Eqn. ( 2) is no longer a straightforward eigenvalue problem [25].Here, we formulate how this eigenvalue problem can be solved for radially symmetric systems.For a given radially symmetric permittivity profile () = () our aim is to find the locations of the quasi-normal modes supported by the resonator.As we are treating radially symmetric resonators, we write the Laplacian in Eqn.(1) in cylindrical or spherical coordinates depending on the number of dimensions we choose to work in.While the calculations are very similar in 2D and 3D, for clarity we proceed with the 3D example providing the 2D calculation in the supplementary material.Separating the angular and radial variables, we write () = ()  ℓ (, ), where   ℓ (, ) are the spherical harmonics.Substituting this into Eqn.(1) turns the Helmholtz equation into We now remove the first derivative term with the substitution () = ()/, to get This has now brought the equation we must solve into the form of the 1D Helmholtz equation, with an additional term related to the angular momentum ℓ.The boundary condition for an outgoing wave that must be imposed is where ℎ (1)  ℓ () is the spherical Hankel function of the first kind.Asymptotically, as  → ∞ the spherical Hankel function goes as ℎ (1)  ℓ () →  − (+1)   /() [31].This means that the boundary condition upon the derivative of the field can be written as
To impose this boundary condition, we shall need to modify the finite difference matrix that will numerically represent the Laplacian operator.We therefore write the boundary condition given by Eqn.(10) in forwards finite difference form At this point, the boundary condition corresponding to quasi-normal modes can be imposed on the Laplacian.However, to find the modes we must remove the dependence upon  from the finite difference Laplacian matrix.This can be achieved by linearising () around a particular frequency  ★ , giving where  = ( ★ ) −  ★   ( ★ ) and  =   ( ★ ).We can now formulate the eigenvalue problem posed by Eqn. ( 2) such that the eigenvalue  only appears on the right-hand side.Defining we can split the Helmholtz equation into This can then be formed into a quadratic eigenvalue problem [32] 0 1 For radially graded structures, this gives both the radial part of the field of the mode and the eigenfrequency .
To verify this method, we compare it to other methods for finding quasi-normal modes.For isotropic cylinders and spheres, Mie theory can be used to find the complex frequencies analytically.Full-wave solvers such as COMSOL Multiphysics [33] can be employed to find the complex eigenfrequencies of spatially varying structures.Both COMSOL and the method we have outlined here require frequencies to search around, so one must already have an idea of roughly where the resonance is located.Figure 1 shows the comparison of the three methods of finding the quasi-normal mode frequencies of isotropic cylinders and spheres.We have chosen the cylinder and sphere to have radius  = 550 nm, with the cylinder having a permittivity of  = 4 and the sphere  = 2. Details of the Mie theory are given in Supplementary Material.With a method of finding both the complex frequency and mode profile of the quasi-normal modes of cylinders and spheres, we proceed to consider two methods for engineering the permittivity such that resonances are placed at particular complex frequencies.This allows one to achieve a high level of control over the spectral response of a scatterer.

Shifting Resonances By Finding 'Eigen-Permittivities'
Our first strategy for imposing a given quasi-normal mode frequency on a structure is calculating an overall shift to the permittivity to place a quasi-normal mode at a particular frequency.To this end, we split the permittivity into a a spatially varying part () and a background shift   .This allows us to rearrange the Helmholtz equation ( 1) into an eigenvalue problem for the background permittivity − 1 To use this, one can choose a complex  at which one wants the quasi-normal mode to occur, then solve the eigenvalue problem to find the corresponding (complex) permittivity shift   .Unlike previous results [30], for radially symmetric resonators the Laplacian contains information about the angular momentum of the mode that is placed at the desired .One can therefore control both the spectral properties and the multipolar nature of the resonance simultaneously.It is important to note that in order to solve the eigenvalue problem posed by Eqn.(17), the finite difference Laplacian operator must be modified to have the outgoing wave boundary condition at the edge of the resonator, as was discussed in Section 2. An example of this procedure for a cylinder is shown in simulations are employed to calculate the scattered energy from the cylinder.It is well-known that the spectral response of a resonator with several well spaced quasi-normal modes can be approximated as a series of Lorentzians [23], with the peak locations corresponding to the real part of the quasi-normal modes and the widths corresponding to the imaginary parts.Thus, to verify that moving the resonances has had the desired effect on the scattering from the resonator, we plot Lorentzians of the form over the scattered power data.Throughout, the central frequency  0 is taken from the real part of the mode with the width Γ corresponding to the expected imaginary part and  is the amplitude scale.At  = 150 THz there is a peak corresponding to the ℓ = 3 quasi-normal mode, with a width corresponding to the imaginary part of the quasi-normal mode frequency.The behaviour between the main three resonances can also be understood using the quasi-normal mode framework.For example, at  = 120 − 12 THz there is a ℓ = 0 mode, which contributes the broad bump to the spectrum at around 120 THz.Examining the field profile of the resonator, shown inset, we observe that the multipolar nature of the mode is as expected.Other peaks in the spectrum are explained by the presence of the ℓ = 2 and ℓ = 4 modes in the vicinity of the mode we have moved.
We also apply this method to move the ℓ = 1 dipole mode of a sphere of initial permittivity  = 2. Initially, this mode is located at  = 172 + 50 THz.Setting the target frequency to be  = 150 − 5 THz, corresponding to a lower resonance frequency and 10 times smaller linewidth, we solve the eigenvalue problem posed by Eqn.(17) to obtain a permittivity offset of   = 0.96 − 0.88.The modes of the new system, as well as the scattered power are shown in Figure 3.

Shifting Resonances By Radially Grading the Permittivity
The next method we propose employs a radial grading of the permittivity to provide control over the spectral location of the quasi-normal modes of a resonator.For this design strategy we employ perturbation theory, which for radially symmetric resonators is well developed [34,35].
A small change in the permittivity of the the resonator () is connected to a change in the position of the quasi-normal mode by the following expression [36,37] where  is the radius of the resonator and   () is the field profile associated with the mode.Re-writing the normalisation factor on the denominator as ⟨  |  ⟩ and changing the permittivity at a single point   so that () = Δ( −   ), we find that the gradient of the location of the quasi-normal mode with respect to the permittivity is The introduction of boundary terms into the norm ⟨  |  ⟩ is necessary as the self-adjointness of the wave operator now depends upon the boundary conditions [23].Now, say we would like to place a particular mode at a given complex frequency  target .Defining the figure of merit as we find that the gradient of the figure of merit with respect to the permittivity is This result is key: we have an analytic expression for how to change the permittivity at every radial position in order to minimise (or maximise) our chosen figure of merit.Like the adjoint method for designing graded wave-shaping devices [27,38], this provides a numerically efficient tool for designing graded structures.In this case, we can design graded cylinders or spheres with with specified complex quasi-normal mode frequencies a chosen multipolar character.
This method is demonstrated in Figure 4. Starting from an isotropic cylinder of permittivity  = 4, the mode locations   and fields   are found using the matrix formulation presented in Section 2. Choosing to move the ℓ = 2 mode from  ∼ 250 − 45 THz to  target = 100 − 2 THz, the figure of merit we seek to optimise is defined by Eqn.(21).Evaluating the gradient of the figure of merit according to Eqn. (23), the permittivity profile is iteratively updated using gradient descent [39]  +1 () =   () +  F   , (24) where the index  indicates the iteration number.This expression allows for the entire profile to be updated each step.The progress of the optimisation is shown in Figure 4 a), where the mode moves through the complex plane as the structure is updated.The designed permittivity grading is shown in Figure 4 b) and c), with the mode locations of the structure shown in panel d).We note that in order to move the quasi-normal mode towards the real frequency axis, gain is required so that the scattered power is increased.To verify that the mode has been moved, power scattered from the graded structure has been calculated with finite element full-wave simulations using COMSOL Multiphysics [33], with the field distributions associated with each peak shown inset.A Lorentzian is plotted in Figure 4 e) with a width and central frequency corresponding to the desired location of the quasi-normal mode in the complex plane.
Next, we consider grading the permittivity of an initially isotropic sphere of permittivity  = 2 to move the ℓ = 4 resonance to  target = 250 − 5 THz.Using the same framework, the progression of the optimisation is shown in Figure 5 a), with the designed permittivity profile shown in Figure 5 b).The locations of the new modes as well as the scattered power, calculated using full-wave simulations, is shown in Figure 5  can note that most of the grading occurs at distance / > 0.6.This is because we are working with modes that vanish at the center of the resonator and manipulating the permittivity in the region where the field is larger gives greater control over the mode location.

Conclusion
We have developed two techniques to solve the problem of designing cylindrical and spherical resonators with multipolar resonances at desired complex frequencies.The first approach involves formulating an eigenvalue problem for a complex permittivity shift of the resonator so that a given mode is placed at a particular location.The second method uses quasi-normal mode perturbation theory to establish a connection between a small change in the permittivity distribution and a small change in the location of a quasi-normal mode.This is then used to form an analytic expression for the gradient of a figure of merit, taken to be the difference between the location of the mode and a desired location.With this, the permittivity distribution of the resonator is iteratively updated until the mode is at the desired spectral location.These methods have the benefit of controlling the resonance frequency of the mode, its linewidth and multipolar nature simultaneously.We expect our methods to find utility in a range of problems from metamaterial design and sensing to optical computing and communication.Future developments of our methods might involve including the effects of material dispersion as well as manipulating multiple modes simultaneously.The ability to arrange the several resonances at desired frequencies would be useful in the design of super-scatterers.In addition, control over the polarisation of the scattered mode might allow the design of sensors that are polarisation dependent or enable metasurface functionality to be multiplexed.

Fig. 1 .
Fig. 1.A comparison of different methods for finding the quasi-normal modes of isotropic cylinders or spheres.a) shows the comparison for a cylinder, and b) for a sphere.Both the cylinder and the sphere have radius  = 550 nm.

Figure 2 .Fig. 2 .
Fig. 2. Using the eigen-permittivity method to place the ℓ = 3 mode of an isotropic cylinder.We seek to move the ℓ = 3 mode to 150-i2 THz, requiring a background permittivity of   = 4.15 + 0.007.a) The locations of the quasi-normal modes.b) Full-wave simulations verifying the scattering behaviour of the cylinder.Peak locations and widths can be related directly to pole locations, and the fields || shown inset verify the multipolar nature of the modes.As the inset plots are of the field norm, the ℓ = 3 lobe should exhibit 6 amplitude peaks.

Fig. 3 .
Fig. 3. Finding an eigen-permittivity to place the dipole ℓ = 1 mode of an isotropic sphere at  = 150 − 5 THz.A background permittivity shift of   = 0.96 − 0.88.Locations of the modes are shown in a).b) Scattered power is then calculated using full-wave simulations, with the mode profile at 150 THz shown inset.

Fig. 4 .
Fig. 4. Designing a graded cylinder with a quadrupole resonance at 100-i2 THz.a) The iterative movement of the ℓ = 2 mode from 152-i15 THz to the desired complex frequency of 100 -i2 THz.Panels b) and c) show the resulting permittivity distribution, which has the pole structure shown in d).Full-wave simulations of the scattered power indicate peaks corresponding to each pole, with inset electric fields show that each mode exhibits the expected multipolar nature.

Fig. 5 .
Fig. 5. Design of a graded sphere with a ℓ = 4 resonance at 250 -i5 THz.a) Progression of the mode over the iterative optimisation, giving b) the final radial permittivity grading.c) Mode locations of the final graded sphere correspond directly to peaks in the d) scattered power.Electric field profiles associated with each mode are shown inset.

Funding.
We acknowledge financial support from the Engineering and Physical Sciences Research Council (EPSRC) of the United Kingdom, via the EPSRC Centre for Doctoral Training in Metamaterials (Grant No. EP/L015331/1).J.R.C. also wishes to acknowledge financial support from Defence Science Technology Laboratory (DSTL).S.A.R.H. acknowledges financial support from the Royal Society (URF\R\211033).