Design of a quasi-2D photonic crystal optomechanical cavity with tunable, large $x^2$-coupling

We present the optical and mechanical design of a mechanically compliant quasi-two-dimensional photonic crystal cavity formed from thin-film silicon in which a pair of linear nanoscale slots are used to create two coupled high-$Q$ optical resonances. The optical cavity supermodes, whose frequencies are designed to lie in the $1500$~nm wavelength band, are shown to interact strongly with mechanical resonances of the structure whose frequencies range from a few MHz to a few GHz. Depending upon the symmetry of the mechanical modes and the symmetry of the slot sizes, we show that the optomechanical coupling between the optical supermodes can be either linear or quadratic in the mechanical displacement amplitude. Tuning of the nanoscale slot size is also shown to adjust the magnitude and sign of the cavity supermode splitting $2J$, enabling near-resonant motional scattering between the two optical supermodes and greatly enhancing the $x^2$-coupling strength. Specifically, for the fundamental flexural mode of the central nanobeam of the structure at $10$~MHz the per-phonon linear cross-mode coupling rate is calculated to be $\tilde{g}_{+-}/2\pi = 1$~MHz, corresponding to a per-phonon $x^2$-coupling rate of $\tilde{g}'/2\pi=1$~kHz for a mode splitting $2J/2\pi = 1$~GHz which is greater than the radiation-limited supermode linewidths.

In a recent experimental work, we realized a silicon photonic crystal optomechanical cavity capable of very large position-squared optomechanical coupling [34].Here, we present details of the cavity design and explore the range of possible optomechanical interactions in such a device.As shown schematically in Fig. 1 the structure consists of a double-waveguide photonic crystal cavity in which two individual waveguide cavity modes are coupled by photon tunneling through a mechanically compliant element.The double-waveguide cavity localizes two optical resonances at telecommunication wavelengths λ =1550 nm with high quality factors Q > 10 5 .The optical modes couple efficiently to mechanical modes with frequencies ω m /2π ranging from 6 MHz up to 1 GHz.We propose a tuning scheme based on electrostatic actuation [22], which allows for control of the optical tunneling rate between the two slotted waveguide modes as well as the supermode frequencies.This provides a direct dynamical control over the linear and quadratic optomechanical coupling strengths.Schematic of a photonic crystal implementation of a similar multimode optomechanical system.The structure consists of a pair of top and bottom photonic crystal slabs which are separated from a central photonic crystal slab by nanoscale air slots.A pair of optical waveguide modes localize around each nanoscale slot, propagating along the axial direction (x) of the structure.By varying the photonic crystal unit cell along the length of the structure, one can form optical cavity modes which are localized to a central "defect" region of the structure.Light in the cavity modes surrounding the top and bottom nanoscale air slots (a1, a2) tunnel across the central photonic crystal slab, forming hybridized supermodes (a+, a−).Mechanical motion of the structure includes in-plane flexural motion of the central nanobeam (b3), the top photonic crystal slab (b1), and the bottom photonic crystal slab (b2).
The paper is organized as follows.We start in Sec.II with a brief introduction of the theory of multimode optomechanics.We consider a generic Hamiltonian describing a system of two optical modes and a series of mechanical modes.Coupling between the optical modes occurs via mechanically independent photon tunneling, and via absorption or emission of mechanical phonons.By considering specific symmetries of the optical and mechanical modes we show how different orders of coupling in mechancial amplitude can occur.In Sec.III we detail the design of the specific double-slotted photonic crystal cavity of interest to our work.We show that by controlling the widths of the waveguide slots of the structure it is possible to adjust the photon tunneling rate of the optical modes, and to even completely suppress this tunneling.In Sec.IV we analyze the mechanical resonances of the photonic crystal cavity and identify several mechanical modes with frequencies ranging from MHz to GHz that have the appropriate symmetry for realizing large optomechanical coupling.Finally, in Sec.V we employ a perturbation theory to calculate the strength of the linear self-mode, linear cross-mode, and position-squared coupling a number of the different mechanical resonances of the photonic crystal structure.

II. MULTIMODE HAMILTONIAN
We consider a multimode cavity optomechanical system consisting of two spatially separated optical modes a 1,2 of frequencies ω 1,2 and independent mechanical modes b k of frequencies ω m,k .The individual optical cavity modes are coupled to each other through photon tunneling at a rate J. Conceptually, our structure can be viewed as an on-chip generalization of the membrane-in-middle setups [3,4] in a chip-scale architecture, as it is schematically represented in Fig. 1.The corresponding Hamiltonian can be written as, where g ij,k are the coupling strengths between the optical modes a i , a j and the mechanical mode b k .In the particular case of symmetric optical cavities, a 1 and a 2 are degenerate with ω 1 = ω 2 = ω 0 .Since the individual optical modes are spatially separated, we can to a good approximation neglect the terms proportional to g ij,k with i = j and therefore write g ij,k simply as g i,k .By introducing the supermode basis a ± = (a 1 ± a 2 )/ √ 2, we can diagonalize Ĥopt which yields for the total Hamiltonian, where the frequency difference between the supermodes at zero mechanical displacement is ω + (0) − ω − (0) = 2J.The first term inside the brackets of Eq. ( 2) describes the linear self-mode optomechanical coupling of the a ± supermodes to the mechanical mode of interest.The last term describes the linear cross-mode optomechanical coupling, i.e. the coupling between the a ± supermodes mediated by the mechanical vibrations.In the a ± basis we have for the linear self-mode optomechanical coupling to mechanical mode b k , and the linear cross-mode optomechanical coupling, Following the approach of Ref. [4], we further diagonalize the full Hamiltonian assuming a quasi-static approximation for the mechanical motion.The resulting eigenfrequencies of the a ± ({x k }) supermodes are, where the mechanical displacements x k = x zpf,k ( b † k + bk ) are regarded as a quasi-static variables.Focusing on a single mechanical mode b k and Taylor expanding the optical supermode frequencies as a function of small displacement x k around equilibrium position x k,eq yields, where and Here g ±,k ({x k,eq }) is the linear (self-mode) coupling coefficient and g ±,k ({x k,eq }) is the quadratic coupling coefficient of the a ± ({x k }) supermode to the k th mechanical mode b k .In what follows we will be primarily interested in the linear and quadratic coupling coefficients around the symmetric equilibrium position {x k,eq } = {0}, where the linear coupling is trivially g ±,k ({0}) = g ±,k = (g 1,k + g 2,k )/2 and the quadratic coupling can be related to the linear cross-mode coupling of the a ± supermodes, We return to the more general result for x k,eq = {0} in Sec.V D.
From the above expressions, mechanical modes such that g 1,k = −g 2,k will have vanishing linear self-mode couplings to the a ± supermodes, while mechanical modes such that g 1,k = g 2,k will have vanishing linear cross-mode coupling strength.All intermediate cases such as |g 1,k | = |g 2,k | can of course occur in general.In the following, we introduce a photonic crystal optomechanical resonator supporting multiple mechanical modes of different symmetries and with optimized overlap with the optical modes.We show that by engineering the bandstructure and defect of the photonic crystal the splitting 2J between the optical supermodes can be tuned to arbitrarily small values, which greatly enhances the x 2 -coupling strength.

III. MULTIMODE PHOTONIC CRYSTAL OPTOMECHANICAL CAVITY A. Bandstructure properties
In this paper, we design our structure assuming a silicon thin film device layer of thickness t = 220 nm, Young's modulus Y = 169 GPa, mass density ρ = 2329 kg/m 3 and refractive index n = 3.42.Our initial geometry is a quasi-two-dimensional (quasi-2D) periodic photonic crystal membrane patterned with a triangular lattice of holes (lattice constant a = 480 nm, circular hole radius r = 0.3a) This structure has an in-plane photonic bandgap (a pseudo-bandgap) for guided slab modes of predominantly TE polarization (electric field polarized in the plane of the slab) around a free-space wavelength of λ ≈ 1550 nm.Referring to Fig. 1(b), addition of the slots breaks the translational periodicity in the transverse direction (y-direction), leaving a periodic structure of lattice constant a in the longitudinal (x) direction.Introducing two air slots of width s 1 = s 2 = s ≈ 100 nm into the photonic crystal membrane results in a pair of optical waveguides with guided modes localized to each of the individual air slots due to the pseudo-bandgap.As shown in Fig. 1(b), this splits the triangular lattice into two outer slabs and a central beam.
The optical bandstructure of the guided modes of the double-slotted waveguide structure with a central beam consisting of a single row of air holes is plotted in Fig. 2. The optical bandstructures here and in what follows are computed using the MIT Photonic Bands package [35].To simplify the bandstructure we only plot the optical bands with TE-like polarization (more accurately we plot those modes with even vector symmetry about the mid-plane of the photonic crystal thin-film slab).In addition to the symmetry about the mid-plane of the silicon thin-film slabcorresponding to the σ z mirror operator -the double-slotted waveguide structure also has a symmetry plane about the y = 0 plane as indicated by the green line in Fig. 1(b).The mirror operator corresponding to this symmetry we label σ y , and the modes of the double-slotted photonic crystal waveguide can be categorized by their even and odd parity under σ y .Here we use the labeling convention that the waveguide supermode with even E y electric field profile is called the "even" mode, while the waveguide supermode with odd E y electric field profile is termed the "odd" mode (note that classifying the modes by their vector symmetry would swap the mode labels).The even and odd waveguide bands are shown as solid black and red lines in Fig. 2. The back and red dashed curves correspond to the unguided (in the transverse y-direction) modes of the triangular photonic crystal slabs surrounding each air slot.The shaded grey region corresponds to the region above the light cone of the air cladding surrounding the silicon slab, in which a continuum of radiation modes freely propagate out of the plane of the slab (z direction).The insets to the right of the bandstructure plot in Fig. 2 show the E y field profiles of the even and odd waveguide modes at the edge of the first Brillouin zone (X-point).
In the hopes of achieving large x 2 -coupling strengths, our focus will be to design photonic crystal cavity modes with minimal optical splitting 2J (see Eq. ( 9)).An obvious way to minimize the optical coupling between the slot Γ X waveguide modes would be to increase the separation between the slots.This decreases the photon tunneling rate J, and hence decreases the frequency splitting between the waveguide modes.For instance, we show in Fig. 3 that the frequency splitting between the odd and even waveguide modes at the X-point can be decreased from ∼ 3 THz in the case of waveguide slots separated by a single row of holes, all the way down to ∼ 68 GHz for waveguide slots separated by five rows of holes.In the latter case the bands become nearly degenerate over a significant fraction of the Brillouin zone, whereas in the more strongly-coupled case of a single row of holes the even and odd waveguide bands have different slopes and even cross near the X-point.The mode profiles of the even and odd waveguide supermodes are plotted in Fig. 2 for a separation between slots of a single row of holes, with the odd waveguide modes having a node in the center of the central beam, and thus, more of their energy in the air slots.As described below, for a cavity based upon the coupled waveguides, this feature allows one to control the relative mode frequencies of the even and odd cavity supermodes by changing the size of the air slot gaps in the structure.
In addition to slot width and slot separation, the fact that the even waveguide modes have more of their energy in the central beam may also be exploited to tailor the relative frequencies of the waveguide supermodes.We consider a middle slab consisting of a single row of holes (nanobeam) and analyze the impact of the ellipticity and lattice constant on the frequency of the waveguide modes.We parametrize the ellipticity of the holes by the aspect ratio η = r y /r x , where r y (r x ) are the semi-axis of the ellipse in y (x) direction, so that η = 1 corresponds to circular holes.Further, by setting r y = r 0 / √ η and r x = r 0 √ η where r is circular radius of the unperturbed cells, we keep the air filling fraction invariant between the elliptical holes and the circular holes.Fig. 4(a) shows the frequency shift of the bands at the X-point due to the variation in the aspect ratio of the holes in nanobeam.Only the even band is The three insets illustrate the shape of holes for three different aspect ratios η = 0.55, 1 and 1.7 with (a, r, s1, s2) = (480 nm, 0.3a, 100 nm, 100 nm).(b) Increase of the optical waveguide supermodes' X-point frequency with a decrease of the lattice constant from a = 480 nm to aD = 468 nm = 0.975a calculated with (r, η, s1, s2) = (0.3a, 0.55, 100 nm, 100 nm).
influenced while the odd band does not shift.For completeness we show in Fig. 4(b) the scaling of the even and odd waveguide modes at the X-point versus a scaling of the in-plane lattice constant (i.e., slab thickness constant).As one might expect, the splitting between the waveguide supermodes of different parity are relatively unaffected by the lattice scaling.

B. Optical cavity
To form an optical cavity from the coupled waveguide system described above we need to find a way of closing the ends of the waveguides.To this end, we introduce a "defect" in the waveguide structure by modifying the geometrical parameters of the waveguide along its propagation (x) axis.This defect region of the waveguide is then embedded between two "mirror" sections of the waveguide as depicted in Fig. 1.Localized cavity modes result for a waveguide modification that pushes the X-point waveguide supermode frequencies in the defect unit cells inside the pseudo-bandgap of the unperturbed mirror unit cells [36,37].
Here the cavity is designed by combining the two defects shown in Fig. 4. The lattice constant is decreased from a = 480 nm to a D = 468 nm quadratically over N D = 7 cells at the center of the structure.Simultaneously with the scaling of the lattice constant, the aspect ratio of the air holes in the central nanobeam are increased from η = 0.55 to circular holes with η = 1 while keeping a constant air filling fraction.Fig. 5(a) shows the modification of the waveguide mode frequencies from the unperturbed mirror cells to the central defect cell.The change in the lattice constant pushes the mode frequencies inside the bandgap and the change of hole ellipticity reduces the frequency splitting between the even and odd supermodes.We simulated the full cavity structure using the COMSOL Multiphysics [38] finiteelement-method (FEM) solver.The length of the entire cavity consists of a total of N x = 42 waveguide unit cells along the x-axis, with N D = 7 central defect unit cells.Each of the outer photonic crystal slabs have N y = 9 rows of air holes in the transverse y-direction.Fig. 5(b) and Fig. 5(c) show the E y field profiles for the odd and the even cavity supermodes, respectively.With the given parameters, the odd mode has a simulated free-space wavelength of λ = 1539 nm and a radiation-limited optical Q-factor of Q = 4.2 × 10 6 .The even mode has a free-space wavelength of λ = 1539.8nm and a substantially lower radiation-limited Q-factor of Q = 4 × 10 5 .This asymmetry in Q-factor results from both the larger coupling to even waveguide modes away from the X-point in the mirror section, a result of the non-monotonic dispersion of the even band, and the reduced vertical radiation loss for a mode of odd in-plane symmetry [39].
C. Dependence of the cavity supermode frequency splitting on the slot size Since the optical modes of interest are mostly localized in the air slots, their frequencies are strongly impacted by the width of the air slots s.Increasing s causes the effective refractive index of the optical waveguide modes to decreases, yielding a blue shift of the resulting optical cavity frequencies.FEM simulations of the waveguide supermodes of the double-slotted photonic crystal structure for slot widths ranging from s = 90 nm to s = 100 nm are plotted in Fig. 6 (a).The frequencies increase approximately linearly with s with slightly different slopes, resulting in a crossing of the cavity supermodes at s ∼ 95 nm.This crossing is made more apparent in Fig. 6(b), where the relative splitting between the even (ω + ) and odd (ω − ) supermode frequencies, ∆ω = ω + − ω − = 2J, is plotted.
As we see, by solely modifying the slots width, it is possible to change the splitting between the cavity supermodes (∆ω = ω + − ω − ) from positive to negative.This is again explained by looking at the even and odd field profiles shown in Fig. 5(b) and 5(c).Since the odd mode has a node in the middle of the central nanobeam, its effective refractive index is more sensitive to a change in the air region than the even cavity supermode.Therefore, for an equal change of the slot widths the frequency of the odd cavity supermode shifts more than the frequency of the even cavity supermode.It is also worth noting that the crossing of the cavity supermodes is only possible because the odd and even waveguide bands cross (see Fig. 5(a)), with the even waveguide band having a higher frequency than the odd waveguide band slightly away from the X-point.
Figure 7(a-c) shows anti-crossing curves obtained where one slot width (s 1 ) is fixed and the other is swept from s 2 = s 1 − 5 nm to s 2 = s 1 + 5 nm.In this case, the cavity mode eigenfrequencies anti-cross with a minimal splitting 2J achieved for the symmetric s 2 = s 1 situation.The odd and even supermode branches of the anti-crossing curves were identified in Fig. 7 by looking at the parity of electric field profile of the simulated eigenmodes at the center of the anti-crossing (here we label the branches by the parity of the supermodes at the center of the anti-crossing curve).Notice that the lower frequency branch for slot width s 1 = 90 nm (Fig. 7(a)) is the odd cavity supermode, whereas the higher frequency branch for slot width s 1 = 100 nm (Fig. 7(c The splitting between the odd (red circles) and even (black squares) cavity supermodes is inverted between (a) and (c).In (b), the splitting is reduced from 2J/2π > 100 GHz to 2J/2π = 17 GHz.(d) Plot of the Q-factor of the even and odd supermode branches versus the second slot width s2 for fixed slot width s1 = 95 nm.The parameters of the structure, save the slot width, are the same as in Fig. 5 for all simulations in (a-d).
s 1 = 95 nm, the splitting at the center of the anti-crossing curve is approximately zero (2J/2π = 17 GHz), as it must be if the parity of the upper and lower frequency supermode branches swap.The possibility to control the supermode frequency splitting in this way provides a unique way to tune the strength of the x 2 -coupling according to Eq. ( 9).
In addition to tuning the cavity supermode splitting, adjusting the air slot sizes may also be used to tune the absolute and relative optical Q-factor of the two cavity supermodes.Figure 7(d) plots the simulated radiation-limited Q-factor for the two cavity supermodes as the slot widths are adjusted and the cavity modes sweep through the anti-crossing point.At the center of the anti-crossing curve the two cavity supermodes are to a good approximation even and odd parity modes around the y = 0 plane, which as we pointed out earlier results in a significant difference in their Q-factor.Motion of the outer photonic crystal slabs or the central nanobeam will also change the relative air slot sizes, thus changing the radiation damping of the optical cavity supermodes.As noted in recent theoretical work [40], this dissipative optomechanical coupling can lead to interference of the quantum noise entering the cavity mode system.This in turn can be used to cool the coupled mechanical resonator to the quantum ground-state even in the unresolved sideband regime (bad-cavity limit).More relevantly, in the case of a predominantly x 2 coupled system, this effect can be used to reduce the parasitic linear back-action and enable continuous x 2 measurements of the mechanical motion.

IV. MECHANICAL RESONANCES
Having considered the optical cavity modes of the double-slotted planar photonic crystal, we now analyze the mechanical modes of the structure.In order to support mechanical resonances the photonic crystal slabs are suspended.The optical modes can interact with both flexural and localized acoustic modes of the central nanobeam.In Sec.IV A, we present the flexural modes of the structure.We show that higher orders flexural modes are found to exist up to 1 GHz with significant optomechanical coupling, making them suitable for operation in the resolved sideband regime, where the optical linewidth κ is much smaller than the mechanical frequency ω m [1].In Sec.IV B, we show that the defect developed to form the optical cavity also gives rise to a localized acoustic resonance of a few GHz frequency.

A. Flexural mechanical resonances
The outer slabs and nanobeam behave as three independent mechanical resonators supporting various in-plane and out-of-plane flexural mechanical resonances.Here, we focus on the in-plane flexural modes that are asymmetric with respect to the y = 0 plane mirror operator (σ y = −1), and symmetric with respect to the x = 0 plane mirror operator (σ x = +1), about the center of the structure.These flexural modes are represented in Fig. 8 with exaggerated deformation profiles.Since our main focus will be on the fundamental resonances, we denote the fundamental in-plane flexural modes of the two outer slabs and nanobeam as b 1 , b 2 and b 3 respectively.Their respective frequencies are denoted ω b1 , ω b2 and ω b3 .
In our design, the outer slabs are suspended by tethers of length l t =2.5 µm and width w t =150 nm, yielding fundamental in-plane flexural resonance of ω m /2π 6 MHz.As shown in Fig. 8 (a), these modes correspond to a uniform displacement of the whole slabs.The displacement of one outer slab causes a uniform change of the width of the adjacent slot, and hence a change of the optical supermode frequencies.The in-plane slab modes provide degrees of freedom for the electromechanical tuning of the slot widths.
In Fig. 8 (b), we plot the displacement profiles of the first three lowest frequency (10.8 MHz, 56 MHz and 130 MHz) nanobeam in-plane flexural modes of symmetry {σ x = +1, σ y = −1}.The y-polarized electric field profiles |E y | 2 are plotted in the inset for both the odd and even optical supermodes a ± .The finite extent of the optical modes along the x-axis of the photonic crystal structure limits the region of the nanobeam that will contribute to the optomechanical interaction.As a result, the nanobeam displacement amplitude x 3 can be approximated by a net effective displacement of the whole nanobeam x3 ≈ x 3 , causing one slot width to change by an amount +x 3 and the other to change by −x 3 .Because of this asymmetric displacement, the optomechanical couplings of these flexural mode to the individual slot modes a 1 and a 2 are expected to be equal and of opposite sign.This favors the quadratic and linear cross-mode interaction terms introduced in Eqs. ( 4) and (9).
Higher order flexural modes of the nanobeam have been identified with frequencies up to 1 GHz and are summarized in table I. Assuming a moderate optical quality factor of Q = 5×10 5 , the resolved-sideband regime condition (κ < ω m ) could be met with a flexural mode of frequency ω m /2π = 400 MHz.

B. Localized phononic crystal resonance
In Sec.III B we described how to localize optical waveguide modes of the double-slotted photonic crystal waveguide propagating by engineering a perturbation to the waveguide unit cell in the propagation direction.In particular, we analyzed the photonic bandstructure of the waveguide unit cell and designed a defect based on a combination of change in the lattice constant and change in the central nanobeam hole aspect ratio.Here we study the phononic bandstructure of the nanobeam unit cell and show that our choice of photonic crystal defect parameters makes the nanobeam compatible with the localization of an GHz-frequency acoustic resonance.
Figure 9(a) and 9(b) show the FEM-simulated acoustic bandstructure of the nanobeam unit cell and the frequency shift of the breathing mode band at the Γ-point as the nanobeam transitions from the mirror unit cell geometry to the defect unit cell geometry.The breathing mode band is shown as a solid red curve.The nanobeam unit cell and the corresponding normalized displacement field profile Q(r) of the breathing mode are depicted in Figs.'9(e-g).The localized breathing mode is drawn from the Γ-point of the bandstructure in order to have a significant optomechanical coupling to the optical mode [41].TABLE I. In-plane flexural modes of the nanobeam.We consider the modes with symmetric displacement with respect to σx.The geometric parameters of the nanobeam are the same as in Fig. 5. 3.1 3.4 3.6 3.9 3.9 3.8 3.8 3.8 3.8 3.6

ωm/2π
Figure 9(c) and 9(d) detail separately the shifts of Γ-point frequency of the breathing mode caused by the perturbation in the lattice constant and hole ellipticity of the nanobeam.Increasing the aspect ratio η will push the breathing mode frequency down into the band gap of the unperturbed mirror cells, while decreasing the lattice constant slightly increases the frequency of the Γ-point mode.As summarized in Fig. 9(a), the same defect used to localize the photonic crystal resonances satisfies the conditions to localize a phononic crystal resonance.Note that, in contrast to the mechanical flexural modes, the displacement of the localized breathing mode is symmetric with respect to σ y so the optomechanical coupling to this mode is expected to be the same for both a 1 and a 2 .Therefore, the breathing mode is expected to have negligible linear cross-mode and quadratic coupling strengths.Nevertheless, the breathing mode presents the advantage of being well in the resolved sideband regime and could potentially be used as an auxiliary mechanical mode in multimode optomechanically induced transparency schemes as proposed in Refs.[9,10].

V. OPTOMECHANICAL COUPLING RELATIONS
With knowledge of the optical and mechanical properties of the double slotted waveguide cavity, we can now turn to the calculation of the different optomechanical coupling factors.We utilize a perturbation theory of Maxwell's equations [42] suitable for dealing with both spatial shifts in the dielectric boundaries of the cavity structure as well as stress-induced modifications of the local dielectric constant of the deformed structure.Applying first-order perturbation theory to the numerically computed unperturbed optical field profiles and mechanical field profiles allows us to evaluate both the linear self-mode coupling g ± and the linear cross-mode coupling g +− for the nanobeam flexural modes [5,7,43].By considering the perturbation theory to second-order yields the strength of the quadratic coupling in our structure.Finally, we consider the modification of the coupling strengths due to deviations of the structure from the symmetric equilibrium position of the central nanobeam and outer slabs.

A. Linear self-mode optomechanical coupling
Maxwell's equations in a source-free, linear dielectric medium, yields the following eigenvalue equation for the electric field, where we have used the Dirac notation |E for the electric field E(r) eigenstate with harmonic time dependence e −iωt .
Here, c is the speed of light in vacuum, and (r) is a dielectric constant which is a function of the spatial coordinate r (and most generally a tensor).We are interested in the change in the modal frequency due to an infinitesimal perturbation δα to the dielectric structure.The first-order correction term to the mode frequency is expressed as where (0) (r) is the unperturbed dielectric constant of the structure, |E (0) and ω (0) are electric field and frequency of the harmonic optical mode of interest, * , and δα depends upon the type of perturbation to the dielectric structure.The change in dielectric constant due to mechanical displacement arises from two main contributions.The first contribution is due to shifting of the interface boundary between two dielectric media.In this case the dielectric function is a high-contrast step function which translates displacements normal to the boundary into local modifications of the dielectric seen by the electric field.As proposed by Johnson, et al., in Ref. [42], the appropriate perturbation theory in this so-called moving boundary (MB) problem is given by, δω where 1(2) is the dielectric constant of medium 1 (2) at any point in the boundary surface A between two media of ⊥ |) is the magnitude of the unperturbed electric (displacement) field polarized in the plane (out of the plane) of the boundary surface A between medium 1 and medium 2, and n(r) is the outward unit vector normal pointing from medium 1 into medium 2 on boundary A.
Here, q(r) is the normalized displacement field of the mechanical mode of interest with maximum displacement equal to unity, max |q(r)| = 1.We can also define an effective mass of the mechanical mode in terms of q, where ρ is the mass density of the dielectric material defining the optomechanical structure.This effective mass is the appropriate motional mass for evaluating the zero-point fluctuation amplitude, x zpf = h/2m eff ω m , of the generalized amplitude coordinate corresponding to the point of maximum amplitude of the mechanical mode.
The second contribution to the linear self-mode coupling is due to the photoelastic effect, resulting from the change of the dielectric constant due to the local strain induced by the mechanical displacement.The first-order perturbation to the dielectric tensor is given by, where (0) is the unperturbed dielectric tensor, 0 is the permittivity of free space, p is the fourth rank photoelastic tensor, and S is the symmetric strain tensor.For an isotropic medium this simplifies to, in index notation.In matrix form, The resulting first-order photoelastic (PE) correction to the optical frequency is In the structures studied here, which are made by etching patterns into a thin-film of silicon, the only two media are silicon and vacuum.As such, for the PE contribution to the linear self-mode coupling we utilize the photoelastic tensor coefficients for silicon in evaluating the integral in the numerator [44]: p 11 = −0.0101,p 12 = 0.009 and p 44 = −0.051.We begin by considering the calculation of g +,b1 and g −,b1 for our double-slotted photonic crystal device, i.e., the linear optomechanical couplings of the optical supermodes a + and a − to the fundamental in-plane flexural mode of either outer slab (we choose b 1 in this case).Table II displays the numerically computed coefficients using the perturbation theory described above in terms of the unperturbed optical and mechanical fields.g +,b1 and g −,b1 can also be approximated by fitting the anti-crossing curves of Fig. 7 using the dispersion relation given in Eq. ( 5).The approximate dispersion relation fit values for the linear couplings are also shown in Tab.II, and compare well to the exact perturbation theory values despite the fact that the couplings derived from the dispersion relations using Eqs.( 3) and ( 5) neglect cross-coupling terms between a 1 and a 2 mediated by the mechanics.
The localized breathing mode of the central nanobeam was found by FEM simulations to be at a mechanical frequency of ω m /2π ≈ 4 GHz, with linear coupling rates of g+ /2π =249 kHz and g− /2π =163 kHz to the a ± supermodes, respectively, where we have used the notation g± = g ± x zpf .
TABLE II.Strength of the linear optomechanical coupling of the optical supermodes to the fundamental in-plane flexural modes of the outer slabs for three different slot widths.The second and third columns display the linear coupling strengths calculated numerically using the perturbation theory.The fourth column gives the values of the optomechanical coupling constant obtained by fitting the anti-crossing curves shown in Fig. 7.The geometric parameters of the nanobeam are the same as in Fig. 8.By analogy with Eq. ( 11), the first order perturbation term for the linear cross-mode coupling g ij , where i = j, can be written as [24] g

First order perturbation theory
In the case of the double-slotted photonic crystal of this work, we have for the shifting boundaries contribution to the cross-mode coupling between the supermodes a + and a − at the symmetric ({x k,eq } = 0) equilibrium position (center of the anti-crossing curve of Fig. 7): Note that for the flexural mechanical modes of the photonic crystal structure (either slab or central nanobeam modes) we expect this to be the dominant contribution to the optomechanical coupling.Expanding √ 2 in terms of the slot modes, and neglecting the cross terms such as E * 1 • E 2 due to the small spatial overlap between the fields of the modes localized in separated slots, we obtain Eq. ( 4) again, g +−,k = (g 1,k − g 2,k )/2.
Consider now the flexural modes of the central nanobeam.At the symmetric ({x k,eq } = 0) equilibrium position, the nanobeam's in-plane flexural modes are such that g 1 = −g 2 .Therefore, at the center of the anti-crossing curve g +− is maximal and equal to the linear coupling of the nanobeam mode to the a 1,2 slot modes.Table III shows the numerically computed linear cross-mode coupling rate gij = g ij x zpf for the fundamental and higher order nanobeam in-plane flexural modes along with their respective frequencies simulated for slot sizes s 1 = s 2 = 95 nm using Eq.(19).Due to the tight localization of optical modes (see Fig. 8) we find there is still significant coupling to higher order flexural modes of frequencies all the way up to 1 GHz.For the numerical simulations of the a ± optical supermodes of the double-slotted photonic crystal structure we also find that the radiation-limited optical quality factor is theoretically equal to 5 × 10 6 and 3 × 10 5 for the odd and even modes respectively.Therefore, as noted earlier, we can expect mechanical modes of frequencies ω m /2π > 300 MHz to be in the resolved-sideband regime.

C. Quadratic optomechanical coupling
By extending the perturbation theory to the second order, it is also possible to calculate the x 2 -coupling strength [45][46][47][48].We obtain, for a given optical mode a i δω In the case of the supermodes a + and a − of the symmetric double-slotted photonic crystal structure ({x k,eq } = 0), the first term vanishes and the only contribution to the x 2 -coupling comes from the second term.For optical splittings such that 2J ω 0 , δω which is what we obtained in Eq. ( 9).In Eq. ( 21) we only consider the contribution from the fundamental optical cavity supermodes because the frequency splitting between them is relatively small.Note that another approach [46] has shown that using a large number of spatially overlapping optical modes rather than decreasing the splitting of just two optical modes (as in our case) can also lead to significant x 2 -coupling strengths.The values of g = g (x zpf ) 2 are summarized in Tab.III for the nanobeam in-plane flexural modes up to 884 MHz.Here we assume an optical a ± supermode frequency splitting of 2J/2π = 1 GHz, which is close to the minimum splitting based on the estimated optical quality factors which allows the optical supermodes to be selectively excited and interrogated.
TABLE III.Linear cross-mode optomechanical (vacuum) coupling rates g+− of the optical supermodes to the nanobeam's inplane fundamental and higher order flexural modes.The x 2 -coupling rate g + is inferred using Eq. ( 9) for a minimum splitting of 2J/2π = 1 GHz.The geometric parameters of the nanobeam are the same as in Fig. 8.In the analysis of the optomechanical coupling coefficients described above in Section V A-V C we considered a symmetric double-slotted structure with equal slot widths s 1 = s 2 (equilibrium position {x k,eq } = 0), and thus the calculations were done for the optical supermodes a ± .In practice, we could find this symmetric condition by tuning one of the slabs to adjust the relative slots sizes until the optical frequency spectrum was at the center of the anticrossing curve as shown in Fig. 7. Far away from the center of the anti-crossing, however, the optical supermodes correspond more closely to the individual slot modes a 1 and a 2 , and we expect different optomechanical coupling strengths.Here we describe how the optomechanical coupling coefficients change upon a large, static displacement of the central nanobeam which takes us far from the symmetric condition near the center of the anti-crossing curve.
From the approximate analytical expression of the supermode dispersion (see Eq. 5), we derive here approximate expressions for g ±,b3 (x 3,eq ), g +−,b3 (x 3,eq ) and g ±,b3 (x 3,eq ) as a function of the static displacement amplitude x 3,eq of the fundamental nanobeam mode: The asterisk ( * ) correspond to numerical FEM simulations of the coupling rates using the perturbation theory with different in-plane static displacement x3 of the nanobeam from its symmetric equilibrium position.The geometrical parameters of the simulated double-slotted structure are: (a, r, s1, s2) = (480 nm, 0.3a, 90 nm+x3, 90 nm−x3).The solid lines are theoretical fits based on Eqs.(22)(23)(24).
As discussed in Sec.IV A, x 3,eq can be approximated by a static displacement x3 of the whole nanobeam.Using the perturbative calculation of the optomechanical coupling coefficients from the numerically simulated optical and mechanical fields one can then obtain g ±,b3 (x 3 ), g +−,b3 (x 3 ) and g ±,b3 (x 3 ) by simulating a structure with the nanobeam displaced from its equilibrium position by x3 (this becomes the new "unperturbed" structure in our perturbative calculations).Here we consider a structure with nominal slot widths of s = 90 nm at equilibrium, and the displacement of the nanobeam from its equilibrium position is swept from x3 = −3 nm to x3 = 3 nm in steps of 0.5 nm.At each position we calculate the coupling coefficients between the optical supermodes and the fundamental in-plane flexural mode of the nanobeam.The results of these simulations and calculations are plotted in Fig. 10.A fit to the numerically calculated coefficients using Eqs.(22)(23)(24) show very good agreement.
These results confirm our previous speculations that far from the symmetric equilibrium position the linear optomechanical coupling is the dominant optomechanical interaction between the optical supermodes and the flexural modes of the central nanobeam, while at symmetric equilibrium position of the beam (x 3 = 0) the optomechanical interaction is predominantly x 2 -coupling.The linear cross-mode coupling between the optical supermodes is also maximal at x3 = 0, however its effects are strongly suppressed in the present case since the splitting between the optical resonances (here 2J/2π = 117 GHz) is much larger than the mechanical frequency of the fundamental in-plane nanobeam mode (10.8 MHz) [29,49].Experimentally, a static displacement of the whole nanobeam by x3 can be mimicked by displacing both outer slabs in the same direction by an amount x3 .Conversely, a change in the equilibrium slot size can be achieved by displacing both outer slabs in opposite directions.In our recent experimental realization of the double-slotted photonic crystal cavity structure we used this tuning degree of freedom by integrating a set of independent capacitors on the outer slabs [34].

VI. SUMMARY
We presented a general formalism for studying linear and quadratic coupling of mechanical motion with optical fields in a multimoded optical cavity.This formalism is used to model and design a double-slotted photonic crystal cavity structure previously studied by us experimentally in Ref. [34].The device supports two high-Q optical resonances at telecommunication wavelengths and mechanical resonances spanning both the unresolved sideband and resolved sideband regimes.We find that depending on the symmetry of the photonic crystal structure, the optical supermodes can be made to interact either linearly or quadratically with mechanical motion.In particular, we show that the splitting between the optical supermodes can be strongly suppressed, which provides a significant enhancement of the quadratic (x 2 ) coupling.The linear and quadratic optomechanical coupling coefficients are calculated using perturbation theory and finite-element-method simulations of the (unperturbed) optical modes and mechanical displacement fields.The simulated optical quality factors of the photonic crystal cavity are calculated to be of order Q ≥ 5 × 10 5 , placing mechanical modes of frequencies ω m /2π ≥ 400 MHz in the resolved sideband regime.With such parameters, it is shown that the designed photonic crystal device can achieve zero-point x 2 -couplings as large as g /2π = 1 kHz for the 10.8 MHz fundamental in-plane flexural mode of the structure, and g /2π = 1 Hz for a higher order flexural mode of frequency 884 MHz in the resolved sideband regime, several orders of magnitude larger than in any other proposed optomechanical system to date.Notably, the splitting of the photonic crystal cavity supermodes may be tuned by adjusting the cavity slot widths to a minimal value approaching that of the cavity mode linewidths (∼ 1 GHz), greatly enhancing the achievable x 2 coupling.Although QND measurements of the stored mechanical energy will require substantially lower optical cavity losses [4,5,34], with this scale of x 2 -coupling it is feasible to consider a number of other interesting experiments with slightly less restrictive parameters.These include the quantum measurement of phonon or photon shot noise [50], or utilizing interference of quantum noise in the bad-cavity limit [40], a continuous position measurement of x 2 .The latter measurement requires additional dissipative optomechanical coupling, which the photonic crystal cavity structure of this work also possesses.

FIG. 1 .
FIG. 1.(a) Schematic of a multimode membrane-in-the-middle optomechanical system consisting of a central movable membrane (b3) and two movable end mirrors (b1, b2).Due to tunneling (rate J) of light through the partially transmitting central membrane, the left and right individual optical cavity modes (a1, a2; not shown) hybridize into the supermodes a+ and a−.(b) Schematic of a photonic crystal implementation of a similar multimode optomechanical system.The structure consists of a pair of top and bottom photonic crystal slabs which are separated from a central photonic crystal slab by nanoscale air slots.A pair of optical waveguide modes localize around each nanoscale slot, propagating along the axial direction (x) of the structure.By varying the photonic crystal unit cell along the length of the structure, one can form optical cavity modes which are localized to a central "defect" region of the structure.Light in the cavity modes surrounding the top and bottom nanoscale air slots (a1, a2) tunnel across the central photonic crystal slab, forming hybridized supermodes (a+, a−).Mechanical motion of the structure includes in-plane flexural motion of the central nanobeam (b3), the top photonic crystal slab (b1), and the bottom photonic crystal slab (b2).

FIG. 2 .
FIG.2.Bandstructure of a triangular lattice of air holes in a silicon slab with two line defects (air slots).We plot the TE-like bands with even (vector) symmetry about the slab mid-plane.The blue line and gray shaded region delimit the light cone of the air cladding surrounding the silicon slab.We focus on the fundamental waveguide modes (solid lines) inside the pseudo-bandgap of the triangular lattice.The y-polarization of the electric field, Ey, of the odd and even modes at the X-point are shown in the insets to the right of the bandstructure.Red and blue correspond to the positive and negative normalized amplitude of the y-polarized electric field.The bandstructure here is computed for a silicon slab with the following set of parameters: refractive index n = 3.42, thickness t = 220 nm, lattice constant a = 480 nm, hole radius r = 0.3a and slot widths s1 = s2 = 100 nm.

FIG. 3 .FIG. 4 .
FIG. 3. Simulated bandstructures of the coupled linear waveguide modes for three different separations of the line defects: (a) one, (b) three and (c) five rows of holes in the central nanobeam.The waveguide modes of interest are shown as solid red (odd modes) and black (even modes) lines.The figure shows the decrease in frequency splitting between the targeted optical modes at the X-point as we increase the number of rows in the central nanobeam.The band diagrams are calculated for the same geometrical parameter as in Fig. 2.

FIG. 5 .
FIG. 5. (a) Left: optical bandstructure of the mirror section of the cavity structure.Right: waveguide supermode frequencies at the X-point as the waveguide unit cell transitions from the outer mirror section to the center of the defect region.As described in the main text, this transition involves a change in the lattice constant (a) and in the nanobeam hole ellipticity (η).In the mirror cell, we use a lattice constant a = 480 nm and elliptical holes with major axis of ry = 194 nm and minor axis of rx = 107 nm (η = 0.55).In the defect cell, we use a = 0.975a and η = 1 (circular holes).The slot widths are s1 = s2=100 nm.The cavity is designed with total of Nx = 42 waveguide unit cells along the x-axis, with a central defect region consisting of ND = 7 defect cells.The outer slabs contain Ny = 9 rows of holes in the transverse y-direction.(b) Plot of the FEM-simulated amplitude of the y-polarization of the electric field Ey(r) of the even cavity mode.(c) Plot of the FEM-simulated odd cavity mode.In (b) and (c), red and blue correspond to positive and negative Ey field amplitudes, respectively.

FIG. 6 .
FIG. 6.(a) Cavity eigenfrequencies for a slot width varying from s = 90 nm to s = 100 nm.The frequency of the odd mode is more influenced by a change of the slot width than the even mode, leading to a change in the spectral ordering of the odd and even eigenmodes.(b) Extracted frequency splitting ∆ω = 2J as function of the slot width.Arbitrarily small splittings are achieved around s = 95 nm.The parameters of the structure are identical to Fig. 5.

FIG. 7 .
FIG.7.Simulated anti-crossing curves obtained by varying s2 while s1 is kept constant at (a) 90 nm (b) 95 nm and (c) 100 nm.The splitting between the odd (red circles) and even (black squares) cavity supermodes is inverted between (a) and (c).In (b), the splitting is reduced from 2J/2π > 100 GHz to 2J/2π = 17 GHz.(d) Plot of the Q-factor of the even and odd supermode branches versus the second slot width s2 for fixed slot width s1 = 95 nm.The parameters of the structure, save the slot width, are the same as in Fig.5for all simulations in (a-d).

3 FIG. 8 .
FIG. 8. Normalized displacement profile of (a) the in-plane slab modes and (b) the nanobeam first and higher order in-plane flexural modes.The inset on top of (b) shows the profiles of |Ey| 2 along the waveguides for both the odd and even symmetry optical supermodes.The deformations are exaggerated for clarity.The photonic crystal parameters are the same as in Fig. 5.The central nanobeam is 731 nm wide and 24 µm long.The outer slabs are suspended by tethers of length lt = 2.5 µm and width wt = 150 nm.

FIG. 9 .
FIG. 9. (a,b) Simulated phononic bandstructure of the nanobeam and defect mode drawn from the Γ-point.The breathing mode band is specified by the solid red line.The even (red lines) and odd (black lines) symmetry acoustic modes are defined with respect to σy mirror operator.(a) Shift of the Γ-point frequency of the breathing mode band for a defect formed by both variations in lattice constant and ellipticity of holes.The defect parameters are the same as in Fig. 5. (c) and (d) show the shift of the breathing mode frequency at the Γ point due to variations of the lattice constant and of the holes aspect ratio, respectively.(e) Mechanical beam and (f) normalized displacement field Q(r) of the localized breathing mode.The color scale indicates the magnitude of Q(r).(g) Exaggerated deformation of the structure due to the breathing mode.All acoustic mode simulations were performed using COMSOL [38].