Efficient anchor loss suppression in coupled near-field optomechanical resonators

Elastic dissipation through radiation towards the substrate is a major loss channel in micro- and nanomechanical resonators. Engineering the coupling of these resonators with optical cavities further complicates and constrains the design of low-loss optomechanical devices. In this work we rely on the coherent cancellation of mechanical radiation to demonstrate material and surface absorption limited silicon near-field optomechanical resonators oscillating at tens of MHz. The effectiveness of our dissipation suppression scheme is investigated at room and cryogenic temperatures. While at room temperature we can reach a maximum quality factor of 7.61k ($fQ$-product of the order of $10^{11}$~Hz), at 22~K the quality factor increases to 37k, resulting in a $fQ$-product of $2\times10^{12}$~Hz.


I. INTRODUCTION
The interaction of optical and mechanical fields in microscale devices enables the manipulation and control of mechanical modes vibrating at radio-frequencies. Some remarkable examples resulting from such an optomechanical interaction include the preparation and measurement of harmonic oscillators' quantum ground states [1][2][3], optically induced synchronization between mechanical oscillators [4], phase noise suppression [5] and highly sensitive sensors [6][7][8][9]. A major limitation in these microdevices is mechanical energy loss that leads to reduced sensitivity [10], lower coherence [11], and increased power consumption [12]; it also remains among the most challenging issues in the design and fabrication processes of micromechanical devices.
While mechanical dissipation is ultimately limited by material absorption, it may also depend on the viscosity of the environment gas, anchor losses and surface scattering or absorption. Although the gas damping can be suppressed in vacuum, the necessary anchoring of micromechanical devices leads to radiation of mechanical waves towards the substrate, which is often the dominant dissipation channel [13]. Fortunately, anchor dissipation can be reduced through phononic shielding [14][15][16], tensile-stressed materials [9,17] , mesa isolation [18], and destructive interference of elastic waves [19][20][21]. Phononic shields are very effective for high frequency mechanical modes [22,23], but their very large footprint at low frequencies [16] is not optimum for photonic integrated circuits. The mesa approach is efficient only for out-of-plane mechanical modes and its deep etch of the substrate may also be incompatible with photonic integration. Therefore, although very high Q mechanical resonators have been reported, they can not always be easily integrated with optical cavities using large-scale photonic integration technology. On the other hand, exploring the destructive interference of mechanical waves is a simple method to obtain both low [24] and high [25] frequency high-Q mechanical resonators, without impacting the device's design or footprint. Yet, the constraints of simultaneously supporting optical and mechanical modes still challenge optomechanical device's design.
Near-field optomechanical (NFO) devices [26][27][28][29] can overcome these challenges as their mechanical and optical resonators are separated structures, which interact through the evanescent optical field. Here we demonstrate a silicon NFO device, fabricated in a commercial foundry, that can reach surface-limited Q-factors by efficient suppression of anchor losses through elastic wave interference.

II. DEVICE PRESENTATION AND ANCHOR LOSS SUPPRESSION
Our NFO device is based on a mechanical resonator composed of coupled paddles that interact with a nearby silicon microdisk optical cavity (Fig. 1a). This design allows us to perfectly balance the mechanical waves radiated to the pedestal by each paddle, without changes to the optical cavity design. Furthermore, it sets an interesting platform to study optomechanical arrays, as it allows coupling several mechanical resonators to a single optical cavity [30][31][32].
The mechanical resonator is composed by two square paddles (2 µm × 2 µm), attached on both sides to suspended beams through 4 nanostrings (200 nm wide) separated by a 200 nm gap. The device is defined on a 220 nm silicon-oninsulator (SOI). Because the paddles' motion couple through the supporting beams, symmetric (S) and anti-symmetric (AS) combinations of individual paddle modes are formed, such as the in-plane modes shown in the finite element  method (FEM) numerical simulation of Fig. 1b. In order to reach a perfect balance between mechanical waves radiating to the supporting beams, the length of the back paddle (L) is offset from the front one by a small length δ. FEM calculations of the mechanical modes show that the resonant frequencies of the coupled mechanical modes display the avoided crossing behavior when δ is varied (Fig. 1c), a signature of coupling of the paddles' motion. The calculations also show that the AS mode has consistently higher anchor loss limited Q's when compared to the S mode (Fig. 1d). This difference appears because the S mode induces a larger displacement on the supporting beams, due to the two paddles' in-phase motion, coupling energy from the paddles to the pedestal and leading to a higher loss rate. On the other hand, the AS mode, due to the anti-phase paddle motion, drastically reduces the displacement at the anchor points, minimizing dissipation to the substrate when the two radiated mechanical waves are balanced. Moreover, we observe a rapid increase of this effect when the paddles are more symmetric (δ ≈ 0), indicating that it is possible to eliminate anchor losses in such NFO device simply by balancing the paddles. Note however that, due to geometry asymmetry of the single-sided pedestal supporting the beams, simulations show that the point of minimum dissipation of the AS mode happens with δ ≈ −15 nm. The devices were fabricated through the EpiXfab initiative at IMEC. They are defined on the 220 nm top silicon layer by deep UV (193 nm) photolithography and plasma etching at the foundry. The 2 µm buried oxide is then partially removed using wet-etch (Buffered Oxide Etch) to define the pedestal and grant the device its mechanical degrees of freedom; this post-process step is performed in-house. Each die has a series of devices where the back paddle has its length varied, while the front paddle is fixed.

III. MEASUREMENT SCHEME AND OPTICAL CHARACTERIZATION
In order to readout the motion of the paddles, the front paddle is designed to be 200 nm away from a 5 µm radius disk optical cavity (Fig. 1a) supporting whispering gallery modes [26] (WGM). Optical readout is possible because the motion of the paddles modulates the frequencies of the WGM's through evanescent field perturbation. The figure of merit of this interaction is the optomechanical coupling rate [11], g 0 = (∂ω/∂x)x zpf , which measures the amount of optical frequency shift caused by a displacement with a quantum-mechanical zero-point fluctuation amplitude (x zpf = /4πm eff f m , where is the Planck's constant, m eff the effective mass and f m the resonance frequency). Although there are usually two main contributions to the optical resonance shift, namely boundary motion [33] and photo-elastic effect [34], we consider only the former because the mechanical modes we are interested induce negligible strain throughout the paddle volume. We estimated g 0 of the in-plane modes by employing perturbation theory [33] to calculate ∂ω/∂x for the various radial orders of the optical transverse electric (TE) modes -also calculated using FEM. The g 0 values -for the perfectly balanced device -ranges from tens to several hundred Hz, as shown in Fig. 1e (see appendix B for more information).
In order to probe the devices we couple light from an external cavity tunable laser into the disk resonator through a tapered fiber (Fig. 2a). The transmitted signal reveals the microdisk optical resonances (Fig. 2b), which usually exhibit frequency splitting due to counter-propagating optical mode coupling, induced by surface roughness and paddle-induced scattering [35] (Fig. 2c). Nevertheless, these modes present reasonably high loaded optical quality factors of about Q opt = 40k, similar to other devices fabricated in the same foundry [22]. The optical mode used to probe the mechanical motion is chosen by monitoring the radio-frequency power spectrum strength at the resonant mechanical frequency, which depends on g 0 , optical linewidth, and taper-cavity loading conditions [36]. We identify which optical mode is excited by comparing the measured optical free spectral range (FSR) with FEM simulations over a broad range (1460 nm to 1610 nm -see appendix A for more information), which indicates that we are coupling light to a higher radial order TE modes; indeed, good agreement between experimental and theoretical values for g 0 also supports this assumption (Fig. 1e). Measurement of g 0 and calibration of the PSD into a displacement noise spectrum density [m/Hz 1/2 ] are performed using a calibrated phase-modulation tone close to the mechanical resonances [36]. The mechanical modes are identified by directly comparing measured frequencies with those from FEM simulations, which agree within a 2% margin. All measurements are performed with the optical cavity undercoupled to the taper. Using a homodyne detection scheme [37] we are able to probe the samples with the laser tuned to the center of the optical resonance, thus avoiding any optomechanical backaction that could affect the mechanical quality factor. With the g 0 measured on the experiments (≈ 2π × 450 Hz) we estimate that even for a 10% of optical linewidth deviation from the resonance we would obtain only a 0.3% error in our mechanical Q's, assuming that all the light is coupled into the cavity (critical coupling with correct polarization). Because we control the polarization of the incident light such that there is very little light inside the cavity, the optomechanical feedback in case of laser deviation from resonance is even more negligible.

IV. RESULTS AND DISCUSSION
The room temperature (RT) mechanical quality factors (Q m ) are obtained from the measured power spectral density (PSD) of the two in-plane mechanical modes while the sample is in a vacuum chamber (10 −5 mbar). In Fig. 3(a,b) (orange curves) we show the measured PSD's and their fitted Lorentzians for the device with δ = −50 nm, which has the highest AS mode quality factor, Q RT m,AS = (7.61 ± 0.07)k, and Q RT m,S = (4.53 ± 0.04)k for the S mode. The mechanical resonance frequencies of these coupled modes are around f m = 56 MHz with a frequency splitting of ∆f m = 940 kHz, also in good agreement with the FEM simulations; this corresponds to a f Q-product of 4 × 10 11 Hz. A confirmation that the δ = −50 nm device is indeed the one with best anchor loss suppression can also be drawn from the frequency and optomechanical couplings measured for devices with varying balance between paddles, as shown in Fig. 3(c,d). This device not only has the smallest frequency difference but it also has S/AS modes with almost identical measured optomechanical coupling rates (g 0 ≈ 2π × 450 Hz), which is expected due to its similar masses and frequencies of S/AS modes for the most balanced device.
Comparison of the S and AS modes' quality factors, for devices with different δ (Fig. 3e), shows consistently higher quality factors for the AS modes, suggesting an important contribution from anchor loss. Nevertheless, its modest two-fold improvement compared to the S mode indicates that other loss mechanisms are also playing an important role in the overall dissipation.  Figure 3. a,b) S (a) and AS (b) modes calibrated power spectrum density at room (orange) and cryogenic (green) temperatures; data in colors and fitted Lorentzians in black; fm ≈ 56 MHz. c-f) Asymmetry dependence of frequency difference between S (blue) and AS (red) modes at room temperature (c), g0 measured at room temperature (d), Q-factor at room (e) and cryogenic (f ) temperatures; marks are data while lines are fitted curves of the coupled oscillators model. g) Temperature dependence of quality factor in log × log scale; blue and red lines are the estimated S and AS modes TED limited Q's; the gray area is bounded by the upper and lower AKE limits. In a,b,g): data of device with δ = −50 nm.
In order to assess the role played by the temperature dependent dissipation mechanisms, we insert the sample in a cold-finger cryostat to cool it down to 22 K. At these low temperatures (LT), we observe a high enhancement of 385% for the AS mode quality factor, reaching Q LT m,AS = (37.0 ± 0.6)k (Fig. 3b, green curve), while the S mode increases only by 80%, up to Q LT m,S = (8.2 ± 0.1)k (Fig. 3a, green curve), for the same device with δ = −50 nm. These Q-enhancements are accompanied by a small (0.7%) increase in mechanical frequencies due to an expected material stiffening at LT, resulting in a Qf -product of 2 × 10 12 Hz.
Such a high contrast between the S/AS modes' quality factors at LT clearly indicates the efficiency of the destructive interference scheme. Nonetheless, the calculated anchor loss limited Q m is still much higher than the measured values for the AS mode at LT. This suggests that either destructive interference between the mechanical waves is unbalanced, intrinsic material dissipation, such as thermoelastic damping (TED) [38] or Akhiezer effect (AKE) [39,40], was reached or surface-related dissipation [41,42] is playing an important role at 22 K.
The Q-limit imposed by destructive interference should not depend on temperature and its possible failure is assessed using a simple analytical model of three coupled mass-spring lumped oscillators [20]. The fitted model, which is detailed in the appendix C, predicts not only the frequency and damping dependence on the asymmetry δ, but also takes into account an overall temperature-dependent intrinsic loss (e.g., material and surface losses). The solid lines in Fig. 3(c-f) represent the fitted model prediction, resulting in an fitted anchor loss, for both LT and RT data, that agrees within 1%, while the intrinsic loss varies by a factor 4 (see appendix C for more details). Therefore our model indicates that, despite the similar fitting of LT and RT anchor loss, the observed sharpening in the LT data is reasonably explained by a reduction of the LT intrinsic damping. Also, the slowly varying AS mode Q-enhancement -towards the higher symmetry region (δ ≈ −50 nm) -suggests that the Q is not being limited by variations in the device geometry. Indeed, numerical simulations confirm that small (±100 nm) variations on any of the device transverse dimensions would not quench the AS quality factors to the LT measured levels. Hence, we infer that the AS mode LT Q-factor is not being limited by failure of the destructive interference scheme.
Intrinsic material dissipation, such as the thermoelastic damping (TED, caused by heat flow between mechanically heated/cooled regions) and the Akhiezer effect (AKE, caused by strain-induced perturbation of thermal phonon equilibrium), have a distinct temperature-dependence [38][39][40]43]. Therefore, we investigate their role by increasing the cryostat base temperature from 22 K up to 200 K, while monitoring the mechanical quality factors (Fig. 3g). The measured Q temperature dependence reveals a capped behavior below ∼ 100 K and a power-law reduction up to ∼ 200 K; above this temperature there is an apparent increase in both S/AS quality factors.
The large mismatch between the period of the mechanical oscillation ( 20 ns) and the thermal lifetimes ( µs) in these devices results in very small TED dissipation. From FEM calculations, shown in Fig. 3g (blue and red lines), we expected room-temperature Q limit to be around 10 5 , while the expected low temperature TED-limited Q's are well above the 10 8 level, without ever decreasing below the RT limit (see appendix D for more details). Hence TED cannot explain neither the temperature dependency nor the LT and RT Q limits.
The AKE, while irrelevant at low-temperatures, could be playing an important role at room-temperatures. We verify it by considering the model by Woodruff and Ehrenreich [40] with temperature dependent material properties [44][45][46] (see appendix E for the used values). Among the parameters involved in this model, the Grüneisen parameter (γ) has the largest range of reported values [47]; room-temperature values range from γ RT = 0.17 up to γ RT = 1.5, resulting in an upper and a lower limit for the AKE, respectively. The gray area in Fig. 3g is bounded by these upper and lower limits, and the lower-limit is roughly within a factor 2 above the RT measured values. Despite this mismatch and the large uncertainty in the Grüneisen parameter, Q-factor temperature dependence follows the AKE power-law, which suggests the Akhiezer effect also plays an important role in the overall damping close to room temperature.
Surface-related dissipation, scattering and absorption, are known to be significant in thin silicon devices [41,42], and is found here to play a role across the whole measured temperature range. The first hint of these surfaceeffects at LT arises from the data shown in Fig. 3f, which was actually taken during a second cool-down of the sample. Its maximum Q is 35% smaller than that shown in Fig. 3g (first cool-down), suggesting that some sort of surface modification occurred between measurements. We verified this assumption by measuring the impact of surface treatment of our samples -chemical cleaning [48] (Fig. 4b) and local laser annealing [49] (Fig. 4c) -on the Q m temperature-dependence. Up to 45% improvement over the highest Q m shown in Fig. 3f was observed for the better balanced device (δ = −50 nm) after cleaning. Other tested devices have also shown similar improvement on the quality-factor, as observed on the measurements of the device with δ = −75 nm in Fig. 4c. Although further investigation would be required to further reduce surface losses, the observed Q-enhancement after surface treatment confirm that the LT Q-limit is unlikely due to failure of the destructive interference effect. The second hint suggesting surface losses at higher temperatures arises from the consistent Q-minimum around 220 K, which is coarsely observed in Fig. 3g, but consistent in all liquid-nitrogen cool-downs (Fig. 4). Similar minima have been observed in other structures [41,42,[50][51][52][53] and, although typically associated with dislocation relaxation [50] or surface related effects [41,42], their actual origin still lack proper explanation and its detailed investigation lies beyond the scope of our paper.

V. CONCLUSION
In summary, we have demonstrated the use of destructive interference of elastic waves as an effective approach to obtain high-Q near-field optomechanical resonators using a scalable foundry technology. Our data clearly shows that symmetric modes are highly susceptible to anchor loss, while the antisymmetric modes may be limited by other mechanisms, such as surface and Akhiezer effect. To further increase these devices' performance, dedicated surface treatment steps should be carefully considered during the fabrication process [49,50,53,54], which plays an important role across all temperature ranges. We expect these devices to serve as platforms for studying arrays of optically coupled mechanical resonators, as well as very sensitive force sensors.

A. OPTICAL MODE POLARIZATION AND RADIAL ORDER
We determined the polarization of the optical modes (TE/TM) by comparing finite element method (FEM) numerical solutions of the disk's modes to the measured spectra. We use FEM to solve the modes' azimuthal numbers (m) for various wavelength values (λ 0 ) and calculate their separation (free-spectral-range or FSR). This gives us the dispersion curves shown as solid lines in Fig. 5a, where the radial order increases from bottom up. Then we determine the experimental families of modes from the experimental transmission data, calculate their FSR and plot their dispersion (FSR vs. λ 0 ) over the numerical solution; these are shown as dots in Fig. 5. We marked the optical mode used to probe the mechanical resonators with a green star. We calibrate the wavelength vector of the experimental data using a Mach-Zehnder interferometer and an acetylene cell. This results in very good agreement with the numerically calculated dispersion of a 4.9 µm radius silicon disk, which is only 2% smaller than the nominal disk, below any fabrication precision. We repeated the procedure with TM-like numerical and experimental data (not shown) which also agree very well.
We note that the fact that there are families of modes with smaller FSR than the one used to probe the mechanical resonator (green star) is evidence that we in fact use a high radial order mode to perform the measurements.

B. DEPENDENCE OF G0 ON THE OPTICAL MODE RADIAL ORDER
The moving boundary optomechanical coupling rate (g 0 ) can be estimated through perturbation theory with equation 1 [33]: where x zpf = /(2m ef f Ω m ) is the quantum zero-point fluctuation of a mechanical mode with effective mass m ef f and frequency Ω m , ω 0 is the optical unperturbed resonance frequency, U n is the normalized mechanical displacement perpendicular to the surface S, ∆ = 0 (n 2 in −n 2 out ) is the difference of refractive index inside and outside the material, is the difference of n −2 inside and outside the material, E t is the field tangent to the surface, D n is the electric displacement normal to the surface and E is the total electric field.
We recall that the optical effective mode volume can be defined as [55]: where |E| 2 max is the maximum field intensity and n max is the refractive index of the medium at the point where the field is maximum.
We may then rewrite equation 1 as: where we defined and |Ẽ t | 2 = |E t | 2 /|E| 2 max and |D n | 2 = |D n | 2 /|E| 2 max . In equation 3 only ρ and V ef f depend on the field distribution, hence they must be the terms that determine the dependence of g 0 on the modes' order. Using FEM simulations, solving for the azimuthal number -same as solving for the effective index in a wave-guide -and maintaining the wavelength constant, we calculate V ef f and ρ for modes with radial orders varying from 1 up to 8. Fig. 6a shows the results normalized by the values obtained for the first order mode. While V ef f changes by at most 20%, ρ increases greatly, almost 16 times for the 8 th order mode. Hence the g 0 dependence presented on Fig. 1e of the main text. The behavior of ρ can be explained by the evanescent field of different radial order modes as a function of the distance from the disk border (Fig. 6b). Although the field intensity for each order may vary at the edge of the disk, at the closest paddle's position (dashed line in Fig. 6b) the decay rate difference results in higher field intensities for higher order modes.
We note that the surface integral in ρ should be calculated on all the surfaces of the two paddles but, because of the field decay rate, only the surface closest to the disk contributes significantly to g 0 . This is the reason why the measure of g 0 allows us to determine the better balanced device in our sample. That is because only for this device the amount of amplitude of motion of the front paddle is approximately the same for both S and AS modes.
Nevertheless, one may ask why we measure our devices with the mode indicated in Fig. 5 if there are at least two other modes with higher radial order, hence higher g 0 . The answer is because the signal we obtain doesn't depend only on g 0 but also on the optical line-width, and we observe that the higher the mode order the larger the line-width becomes. Then, as we mentioned on the main text, we chose that particular mode for it was the one that yielded the best signal for our samples. We used a mass-spring lumped coupled oscillators model to explain the tuning fork effect that suppresses anchor loss in our devices. The model consists of two mass-spring oscillators (m 1 and m 2 in Fig. 7), which represent the two paddles in our devices, coupled to a third mass-spring oscillator (m 3 in Fig. 7), which represents the pedestal of our devices. We also include a direct coupling between oscillators 1 and 2, besides the indirect coupling by oscillator 3.
We use the following simplified mathematical model for this system: where we allow the frequencies Ω i = Ω i + iΓ i to be complex, κ 12 = κ 21 is the coupling between resonators 1 and 2 and κ 13 = κ 31 = κ 23 = κ 32 are the couplings of oscillators 1 and 2 with 3, t is the time and x i are the amplitudes of motion of the oscillators. Solving the system of equations 5 we obtain a set of three complex coupled mode eigenfrequencies, two for the symmetric (S) and anti-symmetric (AS) combinations of modes of the oscillators 1 and 2, and a third for the pedestal (not shown). By varying the mass of one of them (e.g. m 2 ) we are able to reproduce the expected avoided crossing of the coupled mode resonant frequencies and the shape of the damping rate observed in our FEM simulations (Fig. 8).
With the eigenvectors of this system we are able to reproduce the δ dependency of g 0 , as presented in figure 3d of the main text. Note that the example of solutions shown in Fig. 8 have minimum AS-mode damping when m 1 = m 2 , hence we add a factor to one of the resonances to take into account the symmetry braking caused by the pedestal in the real devices.
The frequencies of oscillators 1 and 2 (paddles) are estimated as those of doubly clamped silicon beams' first flexural mode: where α 1 = 4.73 is a constant derived from the beam theory [56], Y is the Young's modulus of silicon, t is the string width (200 nm), w is the SOI thickness (220 nm) and L the string length (4 µm). We artificially modify the mass (m) to account for the extra mass the paddles add to the strings, without any changes to the elastic properties. We also included a scaling constant to this frequency to account for minor deviations from the experimental values, which does not affect the general behavior of coupled mode frequencies and damping rates. The frequency of oscillator 3 (pedestal) is estimated to be 76 MHz from FEM simulations of the device without the paddles and nanostrings. This model takes into account losses due to radiation to the substrate (Γ 3 and κ i3 ) and intrinsic channels (Γ 1,2 ), but it does not consider the nature of the intrinsic losses neither does it account for their temperature dependence. We use this model to fit our experimental asymmetry (δ in the main text) dependent data for both room and low temperature, extracting from it the expected material and/or surface limited Q's for those conditions. Table I resumes the resulting parameters that best fit our data, which were used to plot the curves in Fig. 3c,e,f of the main text.

D. THERMOELASTIC DAMPING -TED
To assess the role of thermoelastic damping (TED) in our devices we performed FEM simulations of this effect. We verified the validity of our FEM model by comparing the numerical solution of a doubly clamped beam to the analytical model presented by Roukes and Lifshitz [38]. The result is presented in Fig. 9a, where we artificially varied the material Young's modulus in both models. Once the FEM model was validated we applied it to our devices, using material properties for temperatures varying from 300 K down to 25 K (see section E). The results for TED lmited Q's are in Fig. 9b, where we observe that at room temperature the minimum Q is 80k and it increases to almost 10 10 at 25 K. Anyhow, this effect's model predicts Q limits much higher than the typical measured in our devices, at all temperatures. Fig. 10 summarizes the properties of Silicon used to estimate TED and AKE mechanical Q limits as a function of temperature. For TED we use SOI thermal conductivity from Asheghi et al [57], heat capacity from Desai [44] and coefficient of thermal expansion from Lyon et al [58]. For AKE we use the heat capacity from Desai [44], sound speed and thermal phonon life-time from Lambade et al [45] and the Grüneisen parameter from Philip and Breazeale [46].  Figure 10. Silicon properties versus temperature. a) Heat capacity [44]; b) Grüneisen parameter [46]; c) Coeffcient of thermal expansion [58]; d) Average Debye sound speed [45]; e) Mean thermal phonon life-time [45]; f) Thermal conductivity [57].

E. TEMPERATURE DEPENDENT SILICON PROPERTIES AND THE AKHIEZER MODEL
For the AKE damping we used the expression derived by Woodruff and Ehrenreich [40].
where for sound speed we use the Debye average of the values available for phonons propagating in the [110] direction (v D ). For thermal phonon life-time we use the values for those interacting with acoustic phonons propagating in the [110] direction with polarization in the [110] direction. We do not use the conductivity due to the problems in its relation to the phonon life-time discussed by Ilisavskii [59]. We note that these values are not available for temperatures below 80 K, hence we extrapolate these parameters to obtain estimates below this temperature.