Quasinormal mode theory of elastic Purcell factors and Fano resonances of optomechanical beams

We introduce a quasinormal mode theory of mechanical open-cavity modes for optomechanical resonators, and demonstrate the importance of using the generalized (complex) effective mode volume and phase of the quasinormal mode. We first generalize and fix the normal mode theories of the elastic Purcell factor, and then show a striking example of coupled quasinormal modes yielding a pronounced Fano resonance. Our theory is exemplified and confirmed by full three-dimensional calculations on optomechanical beams, but the general findings apply to all mechanical cavity modes. This quasinormal mechanical mode formalism, when also coupled with a quasinormal theory of optical cavities, offers a unified framework for describing a wide range of optomechanical structures where dissipation is an inherent part of the resonator modes.


I. INTRODUCTION
The ability to describe optical cavities in terms of normalized modes and cavity figures of merit has played a significant role in laser optics and cavity quantum electrodynamics (cavity-QED). Mode theories not only quantify the underlying physics, but they simplify the numerical modelling requirements and are essential for developing quantized modes in quantum field theory. A striking example is the Purcell factor [1], which elegantly describes the enhanced emission rate of a quantum dipole emitter: where λ 0 is the free space wavelength, Q is the quality factor, and V eff is the effective mode volume. Purcell's theory was originally derived for closed cavity systems, though loss is partly accounted for in the definition of Q. The above formula assumes perfect spatial and polarization alignment of the emitter, which is typically achieved at a field antinode.
Recently, a corrected form for Purcell's formula has been derived in terms of quasinormal modes (QNMs) [2], which are the modes of an open cavity resonators with complex eigenfrequencies; the QNMs yield a generalized (or complex) mode volume [2,3],Ṽ eff , and only the real part is used in Purcell's formula. This subtle "fix" can have profound consequences, and applies to a wide range of lossy cavity structures, including plasmonic resonators [4,5] and hybrid structures of metals and dielectric parts [6]. Moreover, very recently, it was also shown how to quantize these QNMs [7], where the dissipation becomes an essential component in explaining the breakdown of the Jaynes-Cumming model for several modes, causing intrinsic quantum mechanical coupling between classically orthogonal modes.
There are significant analogies between optics and acoustics/mechanics, where the wave equations for acoustics is given in terms of the pressure and velocity fields [8], instead of the electromagnetic fields for optics. In typical resonator structures, both systems yield open cavity modes with complex eigenfrequencies. Moreover, optomechanical structures can support coupling between mechanical and optical modes [9], offering a wide range of applications in optomechanics [10][11][12][13][14][15]. Despite these similarities, the common use of cavity mode theories in optics is less developed in elastics, and one rarely talks about "mechanical mode" effective mode volumes.
In optomechanical systems, the radiation forces exerted by photons are exploited to induce, control, and/or measure mechanical motion in resonators over a wide range of length scales. The applications of optomechanics vary widely [10], from the transduction [16] and storage [17,18] of quantum information, to ultra-sensitive mass sensing [19,20]. Other applications include ground-state cooling [21] and nonlinear optomechanics [22]. The stereotypical optomechanical system consists of a laser-driven cavity whose electromagnetic fields exert a radiation pressure force on a mechanical resonator, which then acts back on the cavity mode, causing the two modes to interact. Modern optomechanical systems can take the form of ultra-thin membranes, micro-ring resonators, and nanostructures acting a photonic crystal [23,24] and a phononic crystal [25][26][27] simultaneously [28][29][30][31][32], which have been shown to have direct applications in on-chip quantum information processing [33][34][35][36][37][38].
In most optomechanical mode theories to date, the optomechanical coupling rate g 0 is rarely taken from a firstprinciples model, which is in contrast to modal methods in optics where it is more common to adopt an analytical approach based on the optical modes of the structure. Yet there is clearly a need to describe emerging effects such as mode-to-mode transcription and reservoir engineering in terms of the underlying mode properties of the mechanical cavity modes, both in classical and quantum mechanical problems of interest. A recent theory paper introduced the idea of an elastic Purcell factor [39], which like its optical counterpart, describes an enhancement of a dipole emitter (but now a force dipole), in terms of Q and V eff ; similar to earlier works on optical cavities, the authors used a "normal mode theory," which is in general incorrect for open cavity modes [2]; however, for high Q cavities, the normal mode approach can be a very good approximation, but the theory is still ambiguous. Experimental measurements have also recently been performed on the acoustic Purcell factor [40], and thus there is clearly a need to develop mode theories for such geometries and emerging material systems.
To account for mode dissipation in optical resonators, the modelling of the dominant cavity mode is usually calculated by implementing outgoing boundary conditions (otherwise it has an infinite lifetime). Often this problem is treated with closed boundary conditions or as a Hermitian eigenvalue problem, but this is inconsistent with a finite loss. In fact, it is now known that all open cavity modes yield spatially divergent modes, which are the QNMs described earlier. Figure 1 shows an example elastic QNM from an optomechanical photonic crystal beam, and we note that the mode diverges for spatial position far down the open beam. These QNMs have recently proven to be very powerful in photonics design and simulations [3,41,42]. For the purpose of field normalization and developing mode theories, both optical and mechanical modes are usually assumed to be lossless and then dissipation is added later through system-bath coupling theory, or phenomenologically; in contrast, the QNM approach includes losses from the beginning since the eigenfrequencies are complex, unlike normal modes. The QNMs also quantify coupling parameters in a more complete way, an example of which has been shown for two coupled QNMs of dielectric-cavity systems [6], resulting in striking interference effects that demonstrate how the phase of the mode must be maintained.
In this work, we introduce an intuitive and accurate QNM description for mechanical modesq m , which have complex eigenfrequencies,Ω m = Ω m − iγ m , and Q m = Ω m /2Γ m . For single mode resonators, This QNM formalism allows a rigorous definition of the effective mode volume V eff,m [1][2][3] for mechanical modes [33], or, equivalently in mechanics, the effective motional mass m eff,m = ρV eff,m [33,43,44], which is more commonly used in the optomechanics literature [31,32,34,35,45,46]. We first present a generalized elastic effective mode volumeṼ eff,m (r) using a QNM normalization, and show the problems with using a normal mode volume V eff . We then use this complex, position dependantṼ eff,m (r) to carry out an analytical Green function expansion [5], which can be used to quickly solve a wide range of force-displacement problems in an analogous way to how the photon Green function is being used to carry out light-matter investigations in optics [47]. Then, using the case of two coupled modes, we demonstrate the accuracy of the QNM theory in explaining complex Fano resonances, and demonstrate the clear failure of the usual normal mode theory for these acoustic modes.
The rest of our paper is organized as follows: In Sec. II, we present the main theory details and important formulas, introducing the elastic wave equations, QNMs, generalized effective mode volume, and elastic Purcell formula. In Sec. III, we present numerical calculations for a fully 3D optomechanical beam, first for a single QNM design, and then for coupled QNMs that show a striking Fano resonance. In both cases, we highlight the failure of using a normal mode theory, which is shown to be drastic in the case of two coupled modes. Finally, in Sec. IV, we present closing discussions and our main conclusions.

A. Wave equations, normal modes, quasinormal modes and Green functions
Vibrational modes of solids can be calculated using the linear theory of elasticity [48], where one assumes infinitesimal deformations and stress forces that do not result in "yielding" (deformation point of no elastic return).
We first express Newton's second law in terms of the displacement u(r, t): where ρ is the mass per unit volume and f (r, t) is the force vector or force per unit volume. Assuming a harmonic solution of the form, u(r, t) = u(r, Ω)e −iΩt , and in the absence of a force excitation, we obtain the wave equation in an analogous form to the Helmholtz equation in optics: The stress tensor σ is expressed as [49] σ ij (r, Ω) = c ijkl ∂ k u l (r, Ω), where c ijkl is the fourth-order elasticity tensor (where one sums over repeated indices as with Einstein summation convention), and ∂ k ≡ ∂ ∂x k is the partial derivative with respect to the x k direction. Subjecting Eq. (3) to periodic or hard-wall boundary conditions would yield "normal modes" q m from the following eigenvalue problem: where the eigenfrequencies are real and do not account for dissipation.
It is also useful to define the corresponding Green function G, obtained from [49]: (6) where the ith component of a unit force in the n direction at location r is given by δ in δ(r − r ). The displacement can be written in terms of G as [49]: whereŝ j is the normal vector of the surface S. For an elastic body surrounded by empty space, the second term vanishes as there is no stress at the surface. In the case of where f d is a point force at position r 0 . Using the completeness relation for normal modes, where I is a unit dyadic of a 3x3 matrix (and m = 1, 2, · · · ), the Green function can be obtained from an expansion over the normal modes [49] G(r, r ; For problems in cavity physics with a few modes of interest, the above theory is not that useful or practical, since one needs a continuum of modes. Instead, similar to opencavity problems in optics, one desires to describe the main physics in terms of just a few discrete resonator modes. In an elastic resonator medium with open boundary conditions, the resonator modesq m are obtained from: whereΩ m is the complex eigenfrequency andq m (r) are the QNMs.
We can now exploit techniques that have been recently developed for obtaining QNMs and QNM Green functions in optics [3,42,50,51]. First, using the completeness relation for QNMs [50] where now m = ±1, ±2, · · · , and we assume this condition is satisfied for spatial positions within or near the cavity [3,42]. We can thus expand G in terms of QNMs through To have quantities that better relate to a mode volume in optics, which is also a key quantity in Purcell's formula, we next define an alternative QNM mode through so thatQ 2 m has units of inverse volume. There have been various approaches to obtain normalized QNMs in optics [2-4, 42, 52-55], and in this work we use which is analogous to the optical QNM normalization introduced by Lai et al. [56], and used to introduce the idea of a "generalized mode volume" with optical cavities [2]. Compared to its optical counterpart, we note the following substitutions: (i) speed of light, c → shear speed of sound in bulk material v b s (as we are only considering transverse modes); and (ii), refractive index n B → material density ρ. Note that Eq. (15) needs a careful regularization if evaluated over a large simulation volume [52].
In terms of these new elastic QNMs, the Green function expansion can now be written as where we now only need to perform a sum over just a few dominate modes of interest. We stress that not only does this theory give a rigorous definition of Purcell's formula (as we show below), but it fully accounts for phase effects and non-Lorentzian features from the complex QNMs.

B. Complex mode volumes and the elastic Purcell formula in terms of QNMs
The common approach to obtaining mode volumes in the literature [31][32][33]39], is to define the effective mode volume, in the present case for acoustic modes, using a normal mode normalization: whereQ m is actually the QNM spatial profile (solution with outgoing boundary conditions), but without the normalization of Eq. (15). However, Eq (17) is problematic as the value diverges as a function of space for any dissipative modes (unless they are lossless). In essence it uses a theory for the normal modes Q m , which are solutions to a Hermitian eignevalue, and so it is incorrect to then use a QNM as the mode solution.
To fix this problem, the QNM normalization condition in Eq. (15) allows one to obtain a complex (or generalized) effective mode volume for a mechanical mode: whereṼ eff,m = V eff,m + iV Im eff,m , for a specific mode m. Note also that we allow the effective mode volume to be a function of space here as it formally characterises the mode strength squared at those positions, and is not directly related to the integrated mode volume (which diverges). Thus one can refer to this quantity as a "localized mode volume," where the volume is a useful figure of merit for certain applications of interest. As we will show later, the phase of the QNM can have profound consequences, and is essential to the general theory.
As explained in the introduction, this modal Purcell factor theory has been well exploited in optical cavity physics and cavity-QED for decades, and the underlying physical insight in terms of cavity mode properties would be hard to underestimate as a design tool. Thus, here we aim to derive the expression for an elastic Purcell factor, similar to the work of Schmidt et al. [39], but now in terms of the more appropriate QNM Green functions.
The mechanical Purcell factor evaluated of an elastic emitter (coupled to a QNM), f d = f d n, oriented along direction n, at some position r 0 and at a frequency Ω can be written as , (19) where G 0 is the Green function in a homogeneous medium.
One can think of this as a generalized enhancement factor or generalized Purcell factor, as no mode expansions have been performed yet (which is necessary to connect to the usual single mode Purcell formula). The terms P and P 0 are the radiated power from a dipole (at position r 0 ) in an inhomogeneous and homogeneous medium, respectively. Using an analogous approach to Poynting's theorem, P 0 is formally defined from [39,57]: which can be derived analytically for an isotropic medium [39]: where in which v l and v s are the scalar longitudinal and shear speed of sound in the material, respectively. We can therefore write Eq. (19) as where we have introduced a numerically determined constant η to account for elastic anisotropy of the medium (for relatively isotropic crystals η ≈ 1). The general-medium Green function can now be used with our generalized effective mode volumeṼ eff and Eq. (23) for a multi-mode approximation of a cavity decay rate enhancement, and for various other problems in optomechanics. We thus obtain the multi-QNM Green function expansion evaluated at a point r 0 within the resonator: For a slightly more familiar form, we can write this in terms of the (complex) effective mode volume as: where we have assumed that the force direction is along the dominant polarization component of the mode, which is the usual assumption in Purcell's formula. Finally, using Eq. (23), we have the enhanced emission rate for a single QNM: in which we assume that the Green function's response is dominated by a single mode, and the response is onresonance (Ω = Ω m ). Equation (26) is the elastic Purcell factor evaluated at the resonant mode Ω m at the source point, whereas (25) is generalized for various positions and frequencies. Our expression is consistent with the expression recently presented by [39], but with a "corrected" and generalized effective mode volume, consistent for an open cavity mode with complex eigenfrequencies and an unconjugated norm. Note also that our Green function explicitly includes the QNM phase, which we show below is essential for describing the response function of several QNMs, which can yield highly non-Lorentzian lineshapes.

A. Mechanical mode effective mode volumes and Purcell factors
We now corroborate our mechanical QNM theory against rigorous numerical calculations obtained for threedimensional mechanical cavities of practical interest in optomechanics.
In this work, we are interested in modelling optomechanical crystal (OMC) cavities on dielectric nanobeams. We consider the impressive nanobeam structure developed by Painter and collaborators [34], consisting of periodic holes with a lattice taper region, in which the hole spacing and size changes. The taper region causes a Fabry-Pérot like effect, resulting in spatially overlapping localized optical and mechanical modes. We adapt the original structure to produce higher loss QNMs by using only a small number of holes to form the OMC cavity. This is so that we may test our QNM normalization at low-Q, where the normal mode approximation breaks down more dramatically.
The nanobeam is modeled as anisotropic silicon in free space with elasticity matrix elements (C 11 , C 12 , C 44 ) = (166, 64, 80) GPa. We consider two OMC cavities, one consisting of 3 holes (3H), and one of 5 holes (5H). Cavity design parameters are specified in Fig. 2 along with the beam dimensions. Each cavity design exhibits single mode behaviour over the frequency ranges shown in Fig. 3. We conduct our numerical investigations using the finite element analysis (FEM) commercial software, COMSOL. The beam is simulated to be infinitely long by implementing perfectly matched layers (PMLs), which simulate outgoing boundary conditions.
Employing the eigenfrequency solver, we obtain the dominant QNMs of interest for each cavity atΩ 3H /2π = 5.160 − i0.016 GHz andΩ 5H /2π = 5.168 − i0.002 GHz. The spatial profile of each mode is shown in Figs. 3a and 3c. We evaluate the effective mode volume of the 3H and 5H QNMs at r 3H 0 = (x, y, z) = (0.0, 244.5, 0.0) nm and r 5H 0 = (0.0, 235.5, 0.0) nm, respectively. Using Equation (18) with the QNM normalization in Eq (15), the real part of the generalized effective mode volumeṼ eff is plotted as a function of simulation size against the more commonly used normal mode V NM eff (Note that the surface term in Eq (15) is only applied outside of the hole region as it assumes a constant outgoing medium). The normal mode normalization in Eq (17) assumes that the mode is localized in space. Consequently, when applied to any cavity mode with finite leakage, V NM eff will diverge exponentially when integrated over all space. For less leaky cavities, this divergence is initially quite slow for a small simulation size, whereas the effect is more dramatic for low-Q modes. The QNM normalization, in contrast, converges as the calculation domain is increased, though may eventually oscillate around the correct value and require regularization [2,52]. Using a sufficient calculation domain size, we obtainṼ 3H eff = 0.06497 + 0.00944i µm 3 and V 5H eff = 0.06143−0.00363i µm 3 at r 3H 0 and r 5H 0 , respectively. We now make use of the generalized effective mode volume to investigate the elastic Purcell effect. The numerically exact Purcell factor is calculated from where P inhomo is the power emission from a point load in the inhomogeneous structure (in this case the beam cavity), and P 0 from a homogeneous sample. We employ the frequency domain solver in COMSOL to obtain full numerical calculations (plotted in Figure 3 b,d). See Appendix A for details on the numerical simulations. Figures 3 b,d show the semi-analytical enhancement rate F P (Ω) calculated using Eq. (25) withṼ 3H eff (r 3H 0 ) and V 5H eff (r 5H 0 ) for the dominant QNM modes of interest,Ω 3H andΩ 5H , respectively. In order to determine η, anisotropic material parameters are applied to a numerical model of a homogeneous sample in which the simulation parameters produce a power spectrum agreeing with the analytical expression in Eq. (21) when isotropic material parameters are used. For our frequency range of interest, we find a relatively constant value of η = P iso 0 (Ω)/P aniso 0 (Ω) ≈ 0.9568, where P iso 0 and P aniso 0 are the radiated power from a point load in an isotropic and anisotropic homogeneous medium, respectively. Also plotted in Fig. 3 b,d is the numerically obtained Purcell factor F exact P calculated from fully threedimensional frequency domain point load simulations. To account for numerical fluctuations (see Appendix A), the average of multiple simulations were conducted in which the mesh parameters are slightly changed. Within this fluctuation, we find excellent agreement between the full numerical simulations and the semi-analytical Green function expansion, with the single-mode approximation being sufficient in describing the dominant resonance response of the system. It is worth mentioning that the modal description provides the enhancement of an emitter at any location and frequency given a sufficient number of QNMs (and often just one QNM) which can be calculated in minutes on a single computer. In contrast, normally one must perform a frequency domain calculation for one frequency at a single point load position, which can take anywhere from hours to days for a sufficient spectrum.

B. Coupled mechanical quasinormal modes
Having demonstrated the validity of the generalized effective mode volume with the single mode approximation, we now consider coupled QNMs. The optomechanical cavity used above can support both high-Q and low-Q modes, depending on the design and quantity of holes used. We consider an OMC nanobeam structure with two cavities, one designed for relatively high-Q modes (5H), and one for low-Q modes (3H). Figure 4a shows the nanobeam structure, using the same cavity design parameters outlined in Fig. 2 and a separation of 4 µm between the two cavities. We look at two resonant a b d c QNMs of interest that are close in frequency with spectral overlap. The first of which, QNM 1 with eigenfrequencỹ Ω 1 /2π = 5.172 − i0.012 GHz (Q 1 = 216), is dominated by the 3-hole cavity (mode profile is shown in Figure 4d). The second, QNM 2 withΩ 2 /2π = 5.171 − i0.002 GHz (Q 2 = 1293), is dominated by the 5-hole cavity (Fig. 4e).
We evaluate the the generalized effective mode volume near the antinode of the 3-hole cavity at r c 0 = (2300, 235, 0) nm, and use the multi-QNM Green function expansion to describe the frequency response at this position in Fig. 4b. Here we seem the hybridization of the individual modes studied earlier (Ω 3H andΩ 5H ), where QNM 2 exhibits a Fano-resonance that results in an interference effect in the total decay rate. Note that we have used a 3rd mode in our approximation, QNM 3 (Ω 3 /2π = 5.232 − i0.041 GHz, V eff,3 (r c 0 ) = −0.280 − 0.006 µm 3 ), to allow for a total decay rate that is positive and well behaved in a large frequency range of interest (while the total decay rate is physically meaningful, contributions from individual modes may not always be); note, however, this mode has negligible contribution to the main dominant response (between 5.15 and 5.19 GHz), and the two-QNM description sufficiently describes the system response. Indeed, we once again have excellent agreement with FEM point load simulations (see Fig. 4b).
For a more intuitive understanding of the role of the mode phase φ m , with two QNMs, we can write Eq. (25) in terms of the two dominant QNMs (QNM 1 and QNM 2) as: where we use the normalized Lorentzian function, and we assume the force dipole is projected along the dominant field direction, namely |Q| 2 = |Q ·n f | 2 , though this can easily be generalized. To better clarify the underlying physics of the various phase terms, we also define two other functions, one that neglects the sin contributions: and one that neglects the phase completely: All three agree (Eqs. (28), (30), (31)) only when cos(2φ 1 ) = cos(2φ 2 ) = 1, and sin(2φ 1 ) = sin(2φ 2 ) = 0. Note also that Eq. (31) has the same form as a normal mode solution, though its effective mode volume is different, and the later will in general be overstated (yielding a smaller Purcell factor value). For our numerical example, the QNM phase values at the point of interest are 2φ 1 (r c 0 ) = 0.3204 and 2φ 2 (r c 0 ) = 2.4643 for QNM 1 and QNM 2, respectively. While the QNM 1 phase shift φ 1 is relatively small, the near 180 o phase shift ofQ 2 2 (r c 0 ) results in the negative contribution to the overall enhancement. From these phase values, the cosine values are cos(2φ 1 ) ≈ 0.95 and cos(2φ 1 ) ≈ −0.7, so the latter will contribute as a negative Lorentzian lineshape. Figure 2 shows the Purcell factor predictions from Eqs. (28), (30), and (31), which clearly demonstrates the role of the phase terms. We can see that neglecting the phase of the mode fails entirely in describing the interference effect (Figure 5a). Considering only the cosine terms still provides a negative contribution from QNM 2 (this is equivalent to only using the real part ofṼ eff,1 ), however, the asymmetry of the lineshape requires the full phase. In fact, the sin terms are sin(2φ 1 ) ≈ 0.31 and sin(2φ 1 ) ≈ 0.36, which are significant and certainly cannot be ignored in general. Indeed, the pronounced Fano feature can only be correctly described when the complete phase is used (see Fig. 5b).
Equation (28) provides a clear understanding of phase interactions between coupled modes. However, rather than work with the complex phases, it can be convenient to work in the complex effective mode volume picture, i.e Eq. (25), which is in the spirit of Purcell's formula. A simple way to do this is to write Eq. (28) in terms of the real and imaginary parts of the generalized effective mode volume: where we have effectively replaced the sin and cos terms, as well as |Q m | 2 , with the real and imaginary parts of 1/Ṽ eff,m . It is now easier to see that the normal mode solution using the entirely real V NM eff,m [39] will always be a simple sum of two Lorentzians: as plotted in Fig. 4c, which shows a drastic failure of the normal mode theory. Note that = Eq. (33) for a single mode is equivalent to the one used in the introduction to the elastic Purcell effect [39]. The calculated generalized effective mode volumes of the QNMs of are interest areṼ eff,1 (r c 0 ) = 0.058 − 0.019i µm 3 andṼ eff,2 (r c 0 ) = −0.399−0.321i µm 3 for QNM 1 and QNM 2, respectively. The negative nature ofṼ eff at r c 0 is simply a result of the phase shift caused by the interaction between the two cavity resonances. This is also seen with QNMs in optics [6,41], where the QNM phase causes the Fanolike resonance. Using quantized QNMs in quantum optics gives essentially the same result as classical QNM theory in the bad cavity limit (within numerical precision), but where the interpretation is now through off-diagonal mode coupling [7]; here the dissipation-induced interference cannot be explained by normal mode quantum theories such as the dissipative Jaynes-Cumming model. In this regard, it would be very interesting to develop a quantized QNM theory for mechanical modes as well.

IV. CONCLUSIONS
We have introduced a QNM formalism for mechanical cavity modes and shown that the commonly used normal mode description is problematic for cavity resonances with finite loss. Instead, we have presented and employed a complex, position dependant effective mode volume for mechanical modes using a QNM normalization which can be used to solve a range of force-displacement problems. For validation of the theory, we carried out an analytical Green function expansion using QNMs with the an elastic Purcell factor expression, and found excellent agreement with rigorous numerical simulations for 3D optomechanical beams.
We then demonstrated the accuracy of the QNM theory in explaining interference effects of coupled cavity modes, and pointed out the drastic failure of the usual normal mode theory. Specifically, we explicitly showed the role of the QNM phase in yielding a Fano-like resonance and explained this analytically and numerically from interference effects that are completely absent in a normal mode theory. To the best of our knowledge, this is the first quantitative theoretical description of coupled mechanical modes exhibiting Fano-like features. This QNM approach should serve as a robust and valuable tool in the understanding and and development of emerging optomechanical technologies.

ACKNOWLEDGMENTS
This work was funded by the Natural Sciences and Engineering Research Council of Canada, the Canadian Foundation for Innovation and Queen's University, Canada. We gratefully acknowledge Mohsen Kamandar Dezfouli and Chelsea Carlson for useful discussions. This research was enabled in part by computational support provided by the Centre for Advanced Computing (http://cac.queensu. ca) and Compute Canada (www.computecanada.ca). Figure 6. The COMSOL simulation geometry used for coupled modes consists of 30.5 µm long nanobeam (see Figures 2 and 4a for beam design) with a 2.25 µm PML in every direction, where the radial PML is separated by a 1 µm buffer from the beam center. a Domain highlighted in blue indicates the embedded nanobeam with silicon material parameters. The surrounding domains (grey) use fictitious material parameters approximating a vacuum for inhomogeneous simulations, and use silicon material parameters for the homogeneous simulation. b Regions highlighted in blue show the PML domain, where we have cut out a quadrant of the radial PML to show the interior domains. c Domains highlighted in yellow show the radial PML regions, using a swept mesh with 5 layers. Remaining domains (grey) use a free tetrahedral mesh (see Table I for mesh parameters).
where v is the velocity vector. Power flow calculations in COMSOL were found to be sensitive to mesh geometry, with P inhomo and P 0 fluctuating dramatically with small changes in mesh parameters. However, P ih /P 0 was found to be more convergent provided that the mesh geometry in the simulation of the beam be exactly identical to the simulation mesh of the homogeneous medium. This was achieved by using the same geometry and mesh points for both P inhomo and P 0 simulations, with the material parameters surrounding the beam (see Figure 6a) changed from silicon (homogeneous case) to a fictitious material with elasticity (C 11 , C 12 , C 44 ) = (0, 0, 0) GPa and a density of 0.001 kg/m 3 in order to approximate a vacuum (inhomogeneous case). The fictitious material is necessary as meshes in the COMSOL solid mechanics solver must be assigned a material, and the chosen elasticity and density for the vacuum approximation suffice. In fact, we found that using any density less than 0.01 kg/m 3 has negligible effect on the calculated P inhomo and the calculated eigenfrequencies of the QNMs (which agree with the in-vacuum simulations).   (and computationally feasible) solutions provided that an appropriate mesh element size is used for the point load. For absorbing boundary conditions, PMLs with polynomial coordinate stretching are used with the scaling factor and curvature parameter set to 1. Figures 6b and 6c outline the PML simulation domains and their meshing type (see figure caption text). Figure 7 shows a parameter sweep of the maximum mesh element size assigned to the point load, where we can see convergence with small fluctuation for element sizes larger than 0.1 nm. For calculations in this work, a point load mesh element size of around 0.2 nm is used. Small shifts in this parameter (in the region of convergence) effectively minutely changes the simulation mesh, resulting in the fluctuations around the solution. In order to account for this, we take the average solution of simulations with slightly varied mesh sizes.